init: initial commit

1 files changed, 589 insertions, 0 deletions

diff --git a/PVCM/cama/en/ma32 Conditionnement d'une matrice.ipynb b/PVCM/cama/en/ma32 Conditionnement d'une matrice.ipynb
new file mode 100644
index 0000000..9ff805a
--- /dev/null
+++ b/PVCM/cama/en/ma32 Conditionnement d'une matrice.ipynb
@@ -0,0 +1,589 @@
+{
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2,
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import numpy.linalg as lin\n",
+    "\n",
+    "np.set_printoptions(precision=3, linewidth=150, suppress=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Conditionnement d'une matrice\n",
+    "\n",
+    "Soit la matrice A suivante :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[10,  7,  8,  7],\n",
+       "       [ 7,  5,  6,  5],\n",
+       "       [ 8,  6, 10,  9],\n",
+       "       [ 7,  5,  9, 10]])"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A = np.array([[10, 7, 8, 7], [7, 5, 6, 5], [8, 6, 10, 9], [7, 5, 9, 10]])\n",
+    "A"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Une matrice symétrique qui n'a rien de méchant a priori. Son déterminant est 1."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.9999999999999869"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "lin.det(A)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Construisons **b** de telle sorte que la solution du système matriciel A **x** = **b** soit [1,1,1,1] :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[32 23 33 31]\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "array([1., 1., 1., 1.])"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "b = A.sum(axis=1)\n",
+    "print(b)\n",
+    "x = lin.solve(A, b)\n",
+    "x"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Perturbons légèrement **b**. Dans le cas d'une expérience cela s'appelle une erreur de mesure. En informatique\n",
+    "cela peut être le résultat d'erreurs d'arrondi."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.0033319453118976702"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "bp = [32.1, 22.9, 33.1, 30.9]\n",
+    "eb = lin.norm(b - bp) / lin.norm(b) # une erreur se mesure par rapport à la valeur de la donnée\n",
+    "eb"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "On a une erreur sur **b** de l'ordre de 0,3 %. On la note $ ||{\\bf \\delta b}|| \\, / \\,||{\\bf b}||$.\n",
+    "\n",
+    "On pourrait espérer une erreur sur le résultat du même ordre de\n",
+    "grandeur.\n",
+    "Regardons la solution **x** de notre système matriciel perturbé :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([  9.2, -12.6,   4.5,  -1.1])"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "xp = lin.solve(A, bp)\n",
+    "xp"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Cette solution n'a rien à voir avec [1, 1, 1, 1]."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "8.19847546803699"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ex = lin.norm(x - xp) / lin.norm(x)\n",
+    "ex"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "L'erreur est de l'ordre de 8 pour 1.\n",
+    "\n",
+    "L'erreur est 2460 fois plus grande que l'erreur sur **b**."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2460.567236431514"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ex / eb"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Pourquoi ?\n",
+    "\n",
+    "On a \n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "& A ({\\bf x} + {\\bf \\delta x}) = {\\bf b} + {\\bf \\delta b} \\quad \\textrm{et donc} \\\\\n",
+    "& A \\, {\\bf \\delta x} = {\\bf \\delta b} \\; \\textrm{ puisque } A {\\bf x} = {\\bf b} \\quad \\textrm{et finalement}\\\\\n",
+    "& {\\bf \\delta x} = A^{-1} \\, {\\bf \\delta b}\n",
+    "\\end{align}\n",
+    "$$\n",
+    "\n",
+    "Comme A et son inverse sont des applications linéaires on a\n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "& ||{\\bf b}|| \\le ||A|| \\, ||{\\bf x}||\n",
+    "\\quad \\textrm{et} \\quad ||{\\bf \\delta x}|| \\le ||A^{-1}|| \\, ||{\\bf \\delta b}||\n",
+    "\\end{align}\n",
+    "$$\n",
+    "\n",
+    "donc \n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "\\frac{||{\\bf \\delta x}||}{||{\\bf x}||}  \\le ||A^{-1}|| \\, \\frac{||{\\bf \\delta b}||}{||{\\bf x}||}\n",
+    "\\le ||A^{-1}|| \\, ||A|| \\, \\frac{||{\\bf \\delta b}||}{||{\\bf b}||}\n",
+    "\\end{align}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "3009.5787080586942"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "lin.norm(lin.inv(A)) * lin.norm(A)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Le problème est là.\n",
+    "\n",
+    "On appelle cela le conditionnement de A :\n",
+    "    \n",
+    "cond(A) = $||A^{-1}|| \\, ||A||$\n",
+    "\n",
+    "**Une matrice mal conditionnée va générer des erreurs de calcul lors de la résolution du système matriciel.**"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2984.0927016757555"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.linalg.cond(A) # scipy n'a pas le conditionnement mais numpy l'a. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Bizarre cette différence avec le calul ci-dessus qui a donné 3009."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Perturbons la matrice\n",
+    "\n",
+    "Que ce passe-t-il si on perturbe A et non b ?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[ 0.098,  0.43 ,  0.206,  0.09 ],\n",
+       "       [-0.153,  0.292, -0.125,  0.784],\n",
+       "       [ 0.927, -0.233,  0.583,  0.058],\n",
+       "       [ 0.136,  0.851, -0.858, -0.826]])"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Erreur relative sur A : 0.06868857112100454\n"
+     ]
+    }
+   ],
+   "source": [
+    "np.random.seed(0)\n",
+    "\n",
+    "dA = 2 * np.random.random(size = A.shape) - 1\n",
+    "display(dA)\n",
+    "ea = lin.norm(dA) / lin.norm(A)\n",
+    "print('Erreur relative sur A :', ea)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[10.098,  7.43 ,  8.206,  7.09 ],\n",
+       "       [ 6.847,  5.292,  5.875,  5.784],\n",
+       "       [ 8.927,  5.767, 10.583,  9.058],\n",
+       "       [ 7.136,  5.851,  8.142,  9.174]])"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "Ap = A + dA\n",
+    "Ap"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([-12.365,  15.574,  10.146,  -5.94 ])"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "xp = lin.solve(Ap, b)\n",
+    "xp"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "11.432687335993894"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ex = lin.norm(xp - x) / lin.norm(x)\n",
+    "ex"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "166.44235204505293"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ex / ea"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "On note que l'erreur est nettement moins grande. La raison est qu'on n'a pas trouvé l'erreur sur A qui\n",
+    "perturbera le plus possible le résultat. En fait ce n'est pas que le conditionnement de A qui compte,\n",
+    "l'erreur est aussi importante. Deux erreurs de même norme pertuberont différemment le résultat.\n",
+    "\n",
+    "\n",
+    "Notons aussi que dans ce cas les maths sont un peu différente mais on retrouve le conditionnement de A :\n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "& (A + \\Delta A) \\, ({\\bf x} + {\\bf \\delta x}) = {\\bf b} \\quad \\textrm{et donc} \\\\\n",
+    "& A \\, {\\bf \\delta x} + \\Delta A \\, ({\\bf x} + {\\bf \\delta x}) = 0 \\; \\textrm{ puisque } A {\\bf x} = {\\bf b} \\quad \\textrm{et finalement}\\\\\n",
+    "& {\\bf \\delta x} = -A^{-1} \\,\\Delta A \\, ({\\bf x} + {\\bf \\delta x}) \\quad \\textrm{et} \\\\\n",
+    "& ||{\\bf \\delta x}|| \\le ||A^{-1}|| \\, ||\\Delta A|| \\, ||{\\bf x} + {\\bf \\delta x}||\n",
+    "\\end{align}\n",
+    "$$\n",
+    "\n",
+    "Ainsi\n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "\\frac{||{\\bf \\delta x}||}{||{\\bf x} + {\\bf \\delta x}||}  \n",
+    "\\le ||A^{-1}|| \\, ||\\Delta A|| =  ||A^{-1}|| \\, ||A|| \\, \\frac{||\\Delta A||}{||A||}\n",
+    "\\end{align}\n",
+    "$$\n",
+    "\n",
+    "ou\n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "\\frac{||{\\bf \\delta x}||}{||{\\bf x} + {\\bf \\delta x}||}  \n",
+    "\\le cond(A) \\, \\frac{||\\Delta A||}{||A||}\n",
+    "\\end{align}\n",
+    "$$\n",
+    "\n",
+    "L'erreur ne se mesure pas par rapport à **x** mais par rapport à ${\\bf x} + {\\bf \\delta x}$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Propriétés\n",
+    "\n",
+    "* $cond(A) \\ge 1$ car $Id = A\\, A^{-1}$ et donc $1 \\le ||A||\\, ||A^{-1}||$ (pour les normes subordonnées car $||Id||_F = \\sqrt{n}$)\n",
+    "* $cond(A) = cond(A^{-1})$ par définition du conditionnement\n",
+    "* $\\displaystyle cond_2(A) = \\frac{\\max_i |\\lambda_i|}{\\min_i |\\lambda_i|}$ si la matrice est normale \n",
+    "  où le 2 indique qu'on utilise la norme 2 et $\\lambda_i$ sont les valeurs propres de A"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(2984.092701676269+0j)"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "vp = lin.eigvals(A)\n",
+    "vp.max() / vp.min()  # et voici la différence avec le calcul ||A|| ||inv(A)|| ci-dessus\n",
+    "                     # ca sent l'erreur numérique"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "* si A est unitaire ou orthogonale alors $cond_2(A) = 1$ (une rotation ou symétrie ne va pas agrandir l'erreur)\n",
+    "* le conditionnement n'est pas modifié par transformation unitaire "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Préconditionnement\n",
+    "\n",
+    "Si le conditionnement n'est pas modifié par une transformation unitaire, il l'est par d'autres transformations.\n",
+    "Ainsi\n",
+    "\n",
+    "$$\n",
+    "\\forall A, \\exists B \\; \\textrm{appelée matrice de préconditionnement t.q.} \\quad cond(B\\, A) \\le cond(A)\n",
+    "$$\n",
+    "\n",
+    "aussi on lieu de résoudre $A {\\bf x} = {\\bf b}$ on résoud $B\\, A {\\bf x} = B\\, {\\bf b}$\n",
+    "\n",
+    "Toute la difficulté consiste à trouver une matrice de préconditionnement B qui soit simple à calculer.\n",
+    "\n",
+    "#### Exercice  11.1\n",
+    "\n",
+    "Quelle est la meilleure matrice de préconditionnement de A ? Pourquoi ?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ]
+}
\ No newline at end of file