path: root/PVCM/cama/fr/ma1 np06 Linalg pour le calcul matriciel.ipynb

{
 "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": 1,
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "# Linalg (linear algebra)\n",
    "\n",
    "Numpy intègre le calcul matriciel (ou l'algèbre linéaire) dans sa sous-bibliothèque [numpy.linalg](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html). Pour être efficace il est\n",
    "important que Numpy soit relié aux bibliothèques [Lapack](https://www.netlib.org/lapack/) et [BLAS](https://www.netlib.org/blas/) (la version d'Intel est [MKL](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html)). Ces bibliothèques sont imbatables, Numpy relié à ces bibliothèques sera largement plus rapide qu'un programme dans n'importe quel autre langage qui ne les utilise pas (vous relevez le défi ?).\n",
    "\n",
    "La bibliothèque [Scipy](https://scipy.org/) a aussi une sous-bibliothèque [linalg](https://docs.scipy.org/doc/scipy/reference/linalg.html) qui est très proches. Si vous ne trouvez pas\n",
    "ce que vous cherchez dans la version de Numpy, il peut être intéressant de regarder si Scipy l'a."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "blas_mkl_info:\n",
      "  NOT AVAILABLE\n",
      "blis_info:\n",
      "  NOT AVAILABLE\n",
      "openblas_info:\n",
      "    libraries = ['openblas', 'openblas']\n",
      "    library_dirs = ['/usr/local/lib']\n",
      "    language = c\n",
      "    define_macros = [('HAVE_CBLAS', None)]\n",
      "    runtime_library_dirs = ['/usr/local/lib']\n",
      "blas_opt_info:\n",
      "    libraries = ['openblas', 'openblas']\n",
      "    library_dirs = ['/usr/local/lib']\n",
      "    language = c\n",
      "    define_macros = [('HAVE_CBLAS', None)]\n",
      "    runtime_library_dirs = ['/usr/local/lib']\n",
      "lapack_mkl_info:\n",
      "  NOT AVAILABLE\n",
      "openblas_lapack_info:\n",
      "    libraries = ['openblas', 'openblas']\n",
      "    library_dirs = ['/usr/local/lib']\n",
      "    language = c\n",
      "    define_macros = [('HAVE_CBLAS', None)]\n",
      "    runtime_library_dirs = ['/usr/local/lib']\n",
      "lapack_opt_info:\n",
      "    libraries = ['openblas', 'openblas']\n",
      "    library_dirs = ['/usr/local/lib']\n",
      "    language = c\n",
      "    define_macros = [('HAVE_CBLAS', None)]\n",
      "    runtime_library_dirs = ['/usr/local/lib']\n",
      "Supported SIMD extensions in this NumPy install:\n",
      "    baseline = SSE,SSE2,SSE3\n",
      "    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2\n",
      "    not found = AVX512F,AVX512CD,AVX512_KNL,AVX512_KNM,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as lin\n",
    "\n",
    "np.set_printoptions(precision=3)\n",
    "np.show_config()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "lang": "fr"
   },
   "source": [
    "## Opérations de base\n",
    "\n",
    "On a vu que les opérateurs +, -, \\* et /  sont appliqués terme à terme ce qui est juste pour + et - dans le cadre du calcul matriciel mais pas pour \\* et /.\n",
    "\n",
    "* Le __produit scalaire__ utilise la méthode `dot` ou l'opérateur `@`\n",
    "* La division que l'on peut imaginer comme \n",
    "   * la __résolution__ d'un système matriciel utilise la fonction `solve`\n",
    "   * le calcul de l'__inverse__ de la matrice utilise la fonction `inv`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "multiplication terme à terme : \n",
      " [[ 1  4]\n",
      " [ 9 16]]\n",
      "produit matriciel : \n",
      " [[ 7 10]\n",
      " [15 22]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1,2],[3,4]])\n",
    "print(\"multiplication terme à terme : \\n\",A * A) # tous les opérateurs sont appliqués terme à terme\n",
    "print(\"produit matriciel : \\n\", A.dot(A))        # A @ A can be used too"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Résolution de système matriciel :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x =  [3. 7.]\n",
      "verification :  [17. 37.]\n"
     ]
    }
   ],
   "source": [
    "b = np.array([17,37])\n",
    "x = lin.solve(A, b) # bien mieux que de calculer la matrice inverse (plus rapide et plus stable)\n",
    "print(\"x = \", x)\n",
    "print(\"verification : \", A @ x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Si on le désire vraiment, on peut calculer la matrice inverse (mais c'est plus long) :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A⁻¹ :\n",
      " [[-2.   1. ]\n",
      " [ 1.5 -0.5]]\n",
      "x =  [3. 7.]\n"
     ]
    }
   ],
   "source": [
    "print(\"A⁻¹ :\\n\", lin.inv(A))       # la matrice inverse\n",
    "print(\"x = \", lin.inv(A).dot(b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Enfin __la transposée__ est simplement `T` :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 3],\n",
       "       [2, 4]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Extractions\n",
    "\n",
    "On peut récupérer\n",
    "\n",
    "* la diagonale d'une matrice avec la fonction `diag` (attention le résultat est un vecteur si l'argument est une matrice et une matrice si l'argument est un vecteur)\n",
    "* la matrice triangulaire inférieur avec la fonction `tril` (l pour lower) et supérieure avec `triu` (u pour upper). Il est possible de décaler la diagonale du triangle, voir la doc\n",
    "\n",
    "On peut aussi construire une matrice triangulaire avec des 1 et 0 avec `tri`et donc aussi une matrice triangulaire quelconque :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4., 0., 0.],\n",
       "       [8., 9., 0.],\n",
       "       [5., 5., 4.]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.tri(3,3) * np.random.randint(1, 10, size=(3,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Et voici comme extraire une matrice tridiagonale :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3, 7, 0, 0, 0],\n",
       "       [6, 1, 8, 0, 0],\n",
       "       [0, 3, 5, 9, 0],\n",
       "       [0, 0, 4, 4, 7],\n",
       "       [0, 0, 0, 6, 3]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.random.randint(1, 10, size=(5,5))\n",
    "np.tril(np.triu(A, k=-1), k=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Opérations sur la matrice\n",
    "\n",
    "La bibliothèque offre des fonctions de\n",
    "\n",
    "* décomposition (LU, Choleski, QR, SVD...)\n",
    "* valeurs et vecteurs propres\n",
    "* norme, déterminant, conditionnement et rang"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.349  0.585  0.132  0.359  0.624]\n",
      " [-0.697 -0.557  0.221 -0.298  0.257]\n",
      " [-0.465 -0.061 -0.052  0.668 -0.576]\n",
      " [-0.349  0.585  0.132 -0.564 -0.447]\n",
      " [-0.232  0.036 -0.956 -0.134  0.115]] \n",
      "\n",
      " [[ -8.602  -7.44  -10.927 -13.252 -10.23 ]\n",
      " [  0.      7.527  -0.039   3.907   7.559]\n",
      " [  0.      0.      1.61   -3.513  -0.301]\n",
      " [  0.      0.      0.      4.333   0.656]\n",
      " [  0.      0.      0.      0.      1.302]]\n",
      "\n",
      "Vérification :\n",
      " [[3. 7. 4. 8. 9.]\n",
      " [6. 1. 8. 5. 3.]\n",
      " [4. 3. 5. 9. 4.]\n",
      " [3. 7. 4. 4. 7.]\n",
      " [2. 2. 1. 6. 3.]]\n"
     ]
    }
   ],
   "source": [
    "Q,R = lin.qr(A)   # Q est orthogonale, R est triangulaire supérieur\n",
    "print(Q, '\\n\\n', R)\n",
    "print(\"\\nVérification :\\n\", Q.dot(R))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "#### Valeurs propres et vecteurs propres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([23.099+0.j   , -3.442+1.822j, -3.442-1.822j, -1.408+0.j   ,\n",
       "         1.192+0.j   ]),\n",
       " array([[-0.547+0.j   , -0.191+0.047j, -0.191-0.047j, -0.51 +0.j   ,\n",
       "          0.487+0.j   ],\n",
       "        [-0.457+0.j   , -0.035-0.466j, -0.035+0.466j, -0.462+0.j   ,\n",
       "         -0.423+0.j   ],\n",
       "        [-0.476+0.j   ,  0.441+0.128j,  0.441-0.128j,  0.41 +0.j   ,\n",
       "         -0.519+0.j   ],\n",
       "        [-0.447+0.j   , -0.536+0.j   , -0.536-0.j   , -0.167+0.j   ,\n",
       "         -0.101+0.j   ],\n",
       "        [-0.257+0.j   ,  0.435+0.233j,  0.435-0.233j,  0.575+0.j   ,\n",
       "          0.552+0.j   ]]))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lin.eig(A) # donne les valeurs propres et les vecteurs propres de A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "#### Déterminant, norme etc. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Déterminant : -587.9999999999999\n",
      "Norme 2 : 26.343879744638983 \n",
      "Norme 1 : 32.0\n",
      "Conditionnement : 35.929347867977604\n"
     ]
    }
   ],
   "source": [
    "print(\"Déterminant :\", lin.det(A))\n",
    "print(\"Norme 2 :\", lin.norm(A), \"\\nNorme 1 :\", lin.norm(A, 1) )\n",
    "print(\"Conditionnement :\", lin.cond(A,2))  # I choose norm 2 to compute the condition number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ]
}