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

np.set_printoptions(precision=3, linewidth=150, suppress=True)

Systèmes matriciels¶

Un système de plusieurs équations à autant d'inconnues peut s'écrire sous forme d'un système matriciel de la forme $A {\bf x} = {\bf b}$ c.a.d.

$$ \begin{bmatrix} a_{11} a_{12} \ldots a_{1n} \\ a_{21} a_{22} \ldots a_{2n} \\ \vdots \\ a_{n1} a_{n2} \ldots a_{nn} \\ \end{bmatrix} \; \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \\ \end{bmatrix} = \begin{bmatrix} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \\ \end{bmatrix} $$

Ainsi si votre épicier ne vous donne pas le prix des pommes, poires et bananes mais que vous en achetez 3 fois des quantités différentes, vous pouvez retrouver le prix de chaque fruit. Il suffit de résoudre les 3 équations qui correspondent aux 3 achats ou le système matriciel 3x3 correpondant. Pour résourdre ce dernier la facon intuitive est d'inverser A et de calculer $A^{-1}\, {\bf b}$ :

In [2]:
# A est la quantité de chaque fruit achetée
# x est le prix de chaque fruit
# b est la somme qu'on a payé pour chaque course

A = np.array([[6,5,4], [5,3,2], [7,3,2]])
b = np.array([11.7, 7.9, 9.5])

x = lin.inv(A) @ b
x
Out[2]:
array([0.8, 0.9, 0.6])

Les systèmes matriciels sont probablement l'application la plus importante des matrices car quoi de plus important que de connaitre le prix des bananes ?

Résolution d'un système matriciel¶

On regardera les méthodes suivantes :

  • Méthode du pivot de Gauss (rendre la matrice triangulaire)
  • Décomposition LU (effet de bord de la résolution par le pivot de Gauss)
  • Méthode de Gauss-Jordan (rendre la matrice diagonale)

Note : Si vous trouvez la suite difficile, faites une pause et regardez cette explication sur une matrice 3x3. Elle commence par triangulariser la matrice comme le fait la méthode du pivot de Gauss, puis elle continue pour diagonaliser la matrice comme le fait Gauss-Jordan et finalement elle arrive à l'inverse de la matrice. C'est très bien pour comprendre ou faire les calculs à la main sur une petite matrice, c'est horrible à programmer. Ce que je vous montre c'est la transformation de l'intuition qu'offre cette vidéo en un calcul matriciel rapide et simple à programmer.

Méthode du pivot de Gauss¶

On transforme A en une matrice triangulaire qui permet ensuite de résoudre le système en O(n²) opérations.

Pour cela on commence à mettre des 0 sur la première colonne en dessous de la diagonale. Pour cela il suffit de multiplier A par la matrice $E_1$ suivante :

$$ E_1 = \begin{bmatrix} \;1 \quad 0\; 0 \ldots 0 \\ \frac{-a_{21}}{a_{11}} \, 1\; 0 \ldots 0 \\ \frac{-a_{31}}{a_{11}} \, 0\; 1 \ldots 0 \\ \vdots \\ \frac{-a_{n1}}{a_{11}}\; 0\; 0 \ldots 1 \\ \end{bmatrix} $$

$E_2$ sera la matrice équivalente avec des termes $\frac{-a_{k2}}{a_{22}}$ sous la diagonales afin des 0 dans la 2e colonne de A sous la diagonale, etc.

Si vous ne connaissez pas cette méthode, je vous conseille de vérifier à la main sur une matrice 4x4 que $E_1 \, A$ met bien des zéros là où il faut. Ensuite faites le travail sur la 2e colonne.

Bien sur si on multiplie A par $E_1$, il faut aussi multiplier ${\bf b}$ par $E_1$ pour que le système matriciel soit équivalent. Cette multiplication $E_1 \, \bf b$ peut se calculer nettement plus rapidement qu'avec un produit matrice vecteur. Donnez la bonne formule et vérifiez dans le code ci-dessous.

Système matriciel avec matrice triangulaire¶

Pour finir il faut résoudre U x = c avec U une matrice triangulaire supérieur. Cela se fait en partant de la dernière ligne qui donne directement x[-1] = c[-1] / U[-1,-1]. Une fois x[-1] connu, on en déduit la valeur de x[-2] avec l'avant dernière ligne, etc.

On faisant les calculs on constate qu'il faut environ n²/2 additions et multiplications.

Là encore, si vous ne vous rappelez pas comment on fait, faites le à la main avec une matrice 3x3. Vous pouvez aussi ajouter des print dans le code qui suit pour voir les matrices intermédiaires.

In [3]:
# Cf analyse du code ci-dessous

def solve_gauss(A, b):
    A = A.astype(np.float64)   # si A a des entiers, on va avoir des erreurs de calculs
    for i in range(len(A)-1):
        coefs = - A[i+1:,i] / A[i,i]  # Normalement il faut tester que A[i,i] != 0
        A[i+1:, i:] += np.outer(coefs, A[i, i:]) # ou np.einsum('i,j -> ij', coefs, A[i, i:])
        b[i+1:] += coefs * b[i]
    # A est maintenant triangulaire surpérieur
    res = np.zeros(len(b), dtype=b.dtype)
    res[-1] = b[-1] / A[-1,-1]
    for i in range(len(A)-1)[::-1]:
        res[i] = (b[i] - A[i,i+1:] @ res[i+1:]) / A[i,i]
    return res

Analyse du code¶

Comme E = Id + Coefs, on a E A = A + Coefs A et donc A += Coef A. Ici Coefs est la matrice E ci-dessus sans les 1 sur la diagonale.

Le problème est que Coef @ A fait n³ opérations ce qui est beaucoup trop et surtout inutile sachant que Coef est un vecteur dans une matrice de zéros. Si ce vecteur est dans la colonne $i$ alors seuls les termes de la ligne $i$ de A serviront pour générer ce produit matriciel (à faire sur un papier). De plus on a res[j,k] = Coef[j,i] * A[i,k] ce qui s'écrit en Numpy res = np.outer(coefs, A[i,:]) si coefs est la colonne $i$ de Coefs.

Comme la matrice Coefs n'a des valeurs non nulles qu'en dessous de la diagonale, cela veut dire que le résultat res de la multiplication n'a que des zéros au dessus de la ligne $i$. Donc autant ne faire que le calcul pour ce qui est en dessus de la ligne $i$ : A[i+1:, i:] += np.outer(coefs, A[i, i:]) avec coefs un vecteur de taille $n-i$.

In [4]:
A = 10 * np.random.random(size=(5,5))
b = A.sum(axis=1)
print(f"A\n {A} \nb\n {b}")
solve_gauss(A, b)

A
 [[2.878 0.086 7.04  4.388 7.731]
 [1.997 0.405 9.546 0.681 3.954]
 [9.664 5.849 3.799 5.096 8.049]
 [7.136 3.182 6.452 1.741 3.594]
 [1.064 3.951 8.733 7.093 9.787]] 
b
 [22.122 16.584 32.457 22.105 30.627]
Out[4]:
array([1., 1., 1., 1., 1.])

Complexité du pivot de Gauss¶

Exercice 10.1¶

Calculer la complexité de solve_gauss. Vous devez trouver n³/3 multiplications et autant d'additions plus des termes de degrés inférieurs. Pour compter les opérations je vous suggère de faire un dessin et d'imaginer votre dessin en 3 dimensions.

Décomposition LU (Lower, Upper)¶

Si on a besoin de résoudre plusieurs systèmes matriciels avec A, alors il est préférable de décomposer A en un produit d'une matriciel triangulaire inférieur et d'une matrice triangulaire supérieure.

$A = L\, U$

Pour cela il suffit de relancer le pivot de Gauss mais au lieu de modifier b à chaque itération, on fabrique à la volée l'inverse de la matrices $E_n \ldots E_2\, E_1$ :

$$ E_n \ldots E_2\, E_1 \, A\, x = E_n \ldots E_2\, E_1 \, b \quad \textrm{donc} \\ (E_n \ldots E_2\, E_1)^{-1} \, E_n \ldots E_2\, E_1 \, A\, x = b \\ E_1^{-1} \, E_2^{-1} \ldots E_n^{-1} \; E_n \ldots E_2\, E_1 \, A\, x = b \\ $$

Ce calcul est très simple car ces matrices inverses ont la même forme que les $E_k$ avec simplement le valeurs opposées sous la diagonale :

$$ E_1^{-1} = \begin{bmatrix} \;1 \quad 0\; 0 \ldots 0 \\ \frac{a_{21}}{a_{11}} \, 1\; 0 \ldots 0 \\ \frac{a_{31}}{a_{11}} \, 0\; 1 \ldots 0 \\ \vdots \\ \frac{a_{n1}}{a_{11}}\; 0\; 0 \ldots 1 \\ \end{bmatrix} $$

et le produit $E^{-1} = E_1^{-1} \,E_1^{-1} \,E_2^{-1} \,E_3^{-1} \, \ldots \,E_n^{-1}$ est simplement la concaténation des colonnes (je vous laisse vérifier) :

$$ E^{-1} = \begin{bmatrix} \;1 \quad 0\; \; 0 \ldots 0 \\ \frac{a_{21}}{a_{11}} \; \; 1\; \; 0 \ldots 0 \\ \frac{a_{31}}{a_{11}} \, \frac{a_{32}}{a_{22}} \; 1 \ldots 0 \\ \vdots \\ \frac{a_{n1}}{a_{11}}\; \frac{a_{n2}}{a_{22}}\; \frac{a_{n3}}{a_{33}} \ldots 1 \\ \end{bmatrix} $$

et ainsi on a $L = E^{-1}$ et $U = E_n \ldots E_2\, E_1 \, A$ puisque tout le travail fait sur A est justement pour mettre des 0 sous la diagonale et donc en faire une matrice triangulaire supérieure.

Une fois la décomposition LU faites, résoudre $L\, U \, {\bf x} = {\bf b}$ se fait en 2 étapes :

  • on résoud $L\, {\bf y} = {\bf b}$ ce qui donne ${\bf y}$
  • puis on résoud $U\, {\bf x} = {\bf y}$ et on a notre solution ${\bf x}$

en O(n²) opérations à chaque fois puisqu'il s'agit de matrices triangulaires.

In [5]:
def LU(A):
    L = np.diag([1.,] * len(A))
    for i in range(len(A)-1):
        E = np.diag([1.,] * len(A))
        E[i+1:,i] = - A[i+1:,i] / A[i,i]
        L[i+1:,i] = -E[i+1:,i]
        A[i:, i:] = E[i:,i:] @ A[i:,i:]
    return L, A
In [6]:
A = 10 * np.random.random(size=(5,5))
print("A\n",A)
L,U = LU(A.copy())  # Attention, notre fonction modifie A donc si on veut le réutiliser il faut une copie
print(f"L\n{L}\nU\n{U}")
A - (L @ U)

A
 [[7.512 7.513 5.909 4.879 0.641]
 [8.468 6.44  4.16  1.435 4.431]
 [6.661 4.823 2.936 3.073 4.768]
 [5.927 6.654 1.168 9.273 0.706]
 [0.631 1.075 2.517 6.229 7.547]]
L
[[  1.      0.      0.      0.      0.   ]
 [  1.127   1.      0.      0.      0.   ]
 [  0.887   0.906   1.      0.      0.   ]
 [  0.789  -0.358 122.023   1.      0.   ]
 [  0.084  -0.219 -40.943  -0.357   1.   ]]
U
[[   7.512    7.513    5.909    4.879    0.641]
 [   0.      -2.029   -2.502   -4.066    3.708]
 [   0.       0.      -0.036    2.432    0.838]
 [   0.       0.       0.    -292.814 -100.761]
 [   0.       0.       0.       0.       6.662]]
Out[6]:
array([[ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0., -0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0., -0., -0.],
       [ 0.,  0.,  0., -0., -0.]])

Gauss Jordan¶

L'idée est de diagonaliser A par des multiplications matricielles similaires à celles de Gauss mais qui annulent aussi les termes au dessus de la diagonale :

$$ E_3 = \begin{bmatrix} 1 \; 0\; \frac{-a_{13}}{a_{33}} \; 0 \ldots 0 \\ 0 \; 1\; \frac{-a_{23}}{a_{33}} \; 0 \ldots 0 \\ 0 \; 0\quad 1 \quad 0 \ldots 0 \\ 0 \; 0\; \frac{-a_{43}}{a_{33}} 1 \ldots 0 \\ \vdots \\ 0 \; 0\; \frac{-a_{n3}}{a_{33}} 0 \ldots 1 \\ \end{bmatrix} $$

Cette méthode est plus compacte à écrire que Gauss car il n'y a pas remonté finale, mais elle est plus lente car on modifie tout A et b à chaque itération.

In [7]:
def solve_gauss_jordan(A,b):
    for i in range(len(A)):
        d1 = np.diag([1.,] * len(A))
        d1[:,i] = - A[:,i] / A[i,i]
        A = d1 @ A
        b = d1 @ b
    return b / np.diag(A)
In [8]:
A = 10 * np.random.random(size=(5,5))
b = A.sum(axis=1)
solve_gauss_jordan(A, b)
Out[8]:
array([1., 1., 1., 1., 1.])

Comparaison de la vitesse de méthodes¶

In [9]:
n = 500
A = 10 * np.random.random(size=(n,n))
b = A.sum(axis=1)
In [10]:
%timeit solve_gauss_jordan(A, b)

2.63 s ± 218 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [11]:
%timeit solve_gauss(A, b)   # 2.5 fois plus rapide que Gauss-Jordan

101 ms ± 3.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [12]:
%timeit B = lin.inv(A)  # 6 fois plus rapide que mon pivot de Gauss

16.4 ms ± 227 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [13]:
%timeit lin.solve(A, b) # un peu plus rapide

14.2 ms ± 277 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Pour n=5000, lin.solve est 3 fois plus rapide que lin.inv sur ma machine. Et plus on augmente $n$ et plus l'écart se creuse

Note: Les calculs utilisent tous les processeurs de ma machine, Numpy sait paralléliser et utilise les bonnes bibliothèques (BLAS et Lapack) s'il est bien installé.

Morale¶

  • Une bonne méthode est mieux qu'une mauvaise (et oui !).
  • Une bonne bibliothèque fait toute la différence (vous ne ferez pas mieux que BLAS et Lapack).

Erreurs d'arrondi¶

In [14]:
np.set_printoptions(precision=20, linewidth=150, suppress=True)

np.finfo(np.float32)
Out[14]:
finfo(resolution=1e-06, min=-3.4028235e+38, max=3.4028235e+38, dtype=float32)
In [15]:
e = 1E-6
A = np.array([[e, 1], [1, 2]], dtype='float32')
b = np.array([1., 3.], dtype='float32')
print(f"A\n {A} \nb\n {b}")

A
 [[0.000001 1.      ]
 [1.       2.      ]] 
b
 [1. 3.]
In [16]:
x = solve_gauss(A.copy(),b)
print("x =", x)
A @ x                # on ne retrouve pas exactement b

x = [1.013279 0.999999]
Out[16]:
array([1.      , 3.013277], dtype=float32)

Il y a un problème d'arrondi évident sachant que la solution est [1 + 2*e, 1 - e] :

In [17]:
x = np.array([1 + 2*e, 1 - e], dtype='float32')
print("La solution est :", x)
A @ x             # on retrouve b

La solution est : [1.000002 0.999999]
Out[17]:
array([1., 3.], dtype=float32)

Le problème est qu'on divise les valeurs de la première colonne de A par le premier pivot : 1E-6.

Aussi on multiplie par 1E6 ce qui va faire des grandes valeurs or la précision est au niveau du 6e chiffre, on pert les suivants comme le montre la simple addition qui suit :

In [18]:
x = np.array([1E6], dtype='float32')   # 1 million
y = np.array([1E-2], dtype='float32')  # 0.01
x+y                                    # résultat : 1 million sans rien après la virgule
Out[18]:
array([1000000.], dtype=float32)

À mélanger des grands nombres avec des petits, on fait des erreurs d'arrondi.

C'est vrai pour les calculs matriciels qu'on fait, c'est toujours vrai et c'est pour cela qu'on aime normaliser les données sur lesquelles on travaille lorsqu'on le peut.

Solution au problème d'arrondi dans le cas du pivot de Gauss¶

Il suffit d'échanger les lignes (ou colonnes mais c'est plus compliqué) pour avoir la plus grande valeur, en valeur absolue, comme pivot. C'est la méthode du pivot partiel.

La méthode du pivot total consiste à prendre comme pivot la plus grande valeur parmis la ligne et la colonne concernées.

Exercice pour le TP¶

Programmer la méthode du pivot partiel et vérifier ce que cela donne sur l'exemple ci-dessus qui pose des problèmes d'arrondi.

In [3]:
import numpy as np
A = np.array([1,2,3])
In [4]:
np.isin?
In [6]:
np.isin(5,A)
Out[6]:
array(False)
In [7]:
B = np.array([1,5,6])
In [12]:
np.isin(B,A)
Out[12]:
array([ True, False, False])
In [ ]: