path: root/PVCM/cama/en/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": "en"
   },
   "source": [
    "# Linalg (linear algebra)\n",
    "\n",
    "Numpy integrates matrix calculus (or linear algebra) in its sub-library [numpy.linalg](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html). To be effective it is\n",
    "important that Numpy is linked to the [Lapack](https://www.netlib.org/lapack/) and [BLAS](https://www.netlib.org/blas/) libraries (Intel's version is [ MKL](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html)). These libraries are unbeatable, Numpy linked to these libraries will be much faster than a program in any other language that does not use them (you take up the challenge?).\n",
    "\n",
    "The [Scipy](https://scipy.org/) library also has a [linalg](https://docs.scipy.org/doc/scipy/reference/linalg.html) sub-library which is very similar. If you can't find\n",
    "what you're looking for in Numpy's version, it might be worth looking to see if Scipy has it."
   ]
  },
  {
   "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": {
    "lang": "en"
   },
   "source": [
    "## Basic operations\n",
    "\n",
    "We have seen that the operators +, -, \\*and / are applied term by term which is correct for + and - in the context of matrix calculation but not for \\* and /.\n",
    "\n",
    "* The __dot product__ uses the `dot` method or the `@` operator\n",
    "* The division that we can imagine as\n",
    "   * the __resolution__ of a matrix system uses the `solve` function\n",
    "   * the calculation of the __inverse__ of the matrix uses the `inv` function"
   ]
  },
  {
   "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": "en"
   },
   "source": [
    "Matrix system resolution:"
   ]
  },
  {
   "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": "en"
   },
   "source": [
    "If you really want to, you can calculate the inverse matrix (but it takes longer):"
   ]
  },
  {
   "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": "en"
   },
   "source": [
    "Finally __the transpose__ is simply `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": "en"
   },
   "source": [
    "## Extraction\n",
    "\n",
    "We can extract\n",
    "\n",
    "* the diagonal of a matrix with the `diag` function (note the result is a vector if the argument is a matrix and a matrix if the argument is a vector)\n",
    "* the lower triangular matrix with the function `tril` (l for lower) and upper with `triu` (u for upper). It is possible to shift the diagonal of the triangle, see the doc\n",
    "\n",
    "We can also construct a triangular matrix with 1s and 0s with `tri` and therefore also any triangular matrix:"
   ]
  },
  {
   "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": "en"
   },
   "source": [
    "And here is how to extract a tridiagonal matrix:"
   ]
  },
  {
   "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": "en"
   },
   "source": [
    "## Matrix operations\n",
    "\n",
    "The library provides functions for\n",
    "\n",
    "* decomposition (LU, Choleski, QR, SVD...)\n",
    "* eigenvalues and eigenvectors\n",
    "* norm, determinant, conditioning and rank"
   ]
  },
  {
   "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": "en"
   },
   "source": [
    "#### Eigenvalues ​​and eigenvectors"
   ]
  },
  {
   "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": "en"
   },
   "source": [
    "#### Determinant, norm 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": []
  }
 ]
}