{
 "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": []
  }
 ]
}