Linalg (linear algebra)¶

Numpy intègre le calcul matriciel (ou l'algèbre linéaire) dans sa sous-bibliothèque numpy.linalg. Pour être efficace il est important que Numpy soit relié aux bibliothèques Lapack et BLAS (la version d'Intel est MKL). 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 ?).

La bibliothèque Scipy a aussi une sous-bibliothèque linalg qui est très proches. Si vous ne trouvez pas ce que vous cherchez dans la version de Numpy, il peut être intéressant de regarder si Scipy l'a.

In [1]:
import numpy as np
import numpy.linalg as lin

np.set_printoptions(precision=3)
np.show_config()

blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
blas_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
lapack_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
    not found = AVX512F,AVX512CD,AVX512_KNL,AVX512_KNM,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL

Opérations de base¶

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 /.

  • Le produit scalaire utilise la méthode dot ou l'opérateur @
  • La division que l'on peut imaginer comme
    • la résolution d'un système matriciel utilise la fonction solve
    • le calcul de l'inverse de la matrice utilise la fonction inv
In [2]:
A = np.array([[1,2],[3,4]])
print("multiplication terme à terme : \n",A * A) # tous les opérateurs sont appliqués terme à terme
print("produit matriciel : \n", A.dot(A))        # A @ A can be used too

multiplication terme à terme : 
 [[ 1  4]
 [ 9 16]]
produit matriciel : 
 [[ 7 10]
 [15 22]]

Résolution de système matriciel :

In [3]:
b = np.array([17,37])
x = lin.solve(A, b) # bien mieux que de calculer la matrice inverse (plus rapide et plus stable)
print("x = ", x)
print("verification : ", A @ x)

x =  [3. 7.]
verification :  [17. 37.]

Si on le désire vraiment, on peut calculer la matrice inverse (mais c'est plus long) :

In [4]:
print("A⁻¹ :\n", lin.inv(A))       # la matrice inverse
print("x = ", lin.inv(A).dot(b))

A⁻¹ :
 [[-2.   1. ]
 [ 1.5 -0.5]]
x =  [3. 7.]

Enfin la transposée est simplement T :

In [5]:
A.T
Out[5]:
array([[1, 3],
       [2, 4]])

Extractions¶

On peut récupérer

  • 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)
  • 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

On peut aussi construire une matrice triangulaire avec des 1 et 0 avec triet donc aussi une matrice triangulaire quelconque :

In [6]:
np.tri(3,3) * np.random.randint(1, 10, size=(3,3))
Out[6]:
array([[4., 0., 0.],
       [8., 9., 0.],
       [5., 5., 4.]])

Et voici comme extraire une matrice tridiagonale :

In [7]:
A = np.random.randint(1, 10, size=(5,5))
np.tril(np.triu(A, k=-1), k=1)
Out[7]:
array([[3, 7, 0, 0, 0],
       [6, 1, 8, 0, 0],
       [0, 3, 5, 9, 0],
       [0, 0, 4, 4, 7],
       [0, 0, 0, 6, 3]])

Opérations sur la matrice¶

La bibliothèque offre des fonctions de

  • décomposition (LU, Choleski, QR, SVD...)
  • valeurs et vecteurs propres
  • norme, déterminant, conditionnement et rang
In [8]:
Q,R = lin.qr(A)   # Q est orthogonale, R est triangulaire supérieur
print(Q, '\n\n', R)
print("\nVérification :\n", Q.dot(R))

[[-0.349  0.585  0.132  0.359  0.624]
 [-0.697 -0.557  0.221 -0.298  0.257]
 [-0.465 -0.061 -0.052  0.668 -0.576]
 [-0.349  0.585  0.132 -0.564 -0.447]
 [-0.232  0.036 -0.956 -0.134  0.115]] 

 [[ -8.602  -7.44  -10.927 -13.252 -10.23 ]
 [  0.      7.527  -0.039   3.907   7.559]
 [  0.      0.      1.61   -3.513  -0.301]
 [  0.      0.      0.      4.333   0.656]
 [  0.      0.      0.      0.      1.302]]

Vérification :
 [[3. 7. 4. 8. 9.]
 [6. 1. 8. 5. 3.]
 [4. 3. 5. 9. 4.]
 [3. 7. 4. 4. 7.]
 [2. 2. 1. 6. 3.]]

Valeurs propres et vecteurs propres¶

In [9]:
lin.eig(A) # donne les valeurs propres et les vecteurs propres de A
Out[9]:
(array([23.099+0.j   , -3.442+1.822j, -3.442-1.822j, -1.408+0.j   ,
         1.192+0.j   ]),
 array([[-0.547+0.j   , -0.191+0.047j, -0.191-0.047j, -0.51 +0.j   ,
          0.487+0.j   ],
        [-0.457+0.j   , -0.035-0.466j, -0.035+0.466j, -0.462+0.j   ,
         -0.423+0.j   ],
        [-0.476+0.j   ,  0.441+0.128j,  0.441-0.128j,  0.41 +0.j   ,
         -0.519+0.j   ],
        [-0.447+0.j   , -0.536+0.j   , -0.536-0.j   , -0.167+0.j   ,
         -0.101+0.j   ],
        [-0.257+0.j   ,  0.435+0.233j,  0.435-0.233j,  0.575+0.j   ,
          0.552+0.j   ]]))

Déterminant, norme etc.¶

In [10]:
print("Déterminant :", lin.det(A))
print("Norme 2 :", lin.norm(A), "\nNorme 1 :", lin.norm(A, 1) )
print("Conditionnement :", lin.cond(A,2))  # I choose norm 2 to compute the condition number

Déterminant : -587.9999999999999
Norme 2 : 26.343879744638983 
Norme 1 : 32.0
Conditionnement : 35.929347867977604
In [ ]: