8.5. Matrices

Dans ce paragraphe, les matrices seront représentées par des listes de listes. Par exemple, la matrice \(\begin{pmatrix}1&2&3\\4&5&6\end{pmatrix}\) sera représentée par la liste de listes [[1, 2, 3], [4, 5, 6]] [1].

8.5.1. Produit matriciel

In [1]: def produit(A, B):
   ...:     return [[sum(L[k] * B[k][j] for k in range(len(L))) for j in range(len(B[0]))] for L in A]
   ...: 

In [2]: A = [[1, 2, 3], [4, 5, 6]]

In [3]: B = [[1, 2], [3, 4], [5,6]]

In [4]: produit(A, B)
Out[4]: [[22, 28], [49, 64]]

In [5]: produit(B, A)
Out[5]: [[9, 12, 15], [19, 26, 33], [29, 40, 51]]

8.5.2. Opérations élémentaires

On définit plusieurs opérations élémentaires sur les lignes d’une matrice.

• l’échange de lignes \(L_i\leftrightarrow L_j\)

  In [6]: def echange_lignes(M, i, j):
     ...:     M[i], M[j] = M[j], M[i]
     ...:     return M
     ...: 
  

• la transvection \(L_i\leftarrow L_i+\lambda L_j\)

  In [7]: def transvection_ligne(M, i, j, l):
     ...:     M[i] = [M[i][k] + l * M[j][k] for k in range(len(M[i]))]
     ...:     return M
     ...: 
  

• la dilatation \(L_i\leftarrow\lambda L_i\)

  In [8]: def dilatation_ligne(M, i, l):
     ...:     M[i] = [coeff * l for coeff in M[i]]
     ...:     return M
     ...: 
  

Avertissement

Les fonctions précédentes, modifient la matrice donnée en argument puisqu’une liste est un objet mutable.

8.5.3. Algorithme du pivot de Gauss

A l’aide des opérations élémentaires précédemment définies, on peut alors définir une fonction appliquant l’algorithme du pivot de Gauss à une matrice pour la mettre sous forme échelonnée.

Pour des raisons de stabilité numérique, on recherche le pivot de valeur absolue maximale.

In [9]: def recherche_pivot_lignes(M, i):
   ...:     m = abs(M[i][i])
   ...:     j = i
   ...:     for k in range(i + 1, len(M)):
   ...:         if abs(M[i][j]) > m:
   ...:             j = k
   ...:     return j
   ...: 

In [10]: def pivot_lignes(M):
   ....:     for i in range(len(M)):
   ....:         j = recherche_pivot_lignes(M, i)
   ....:         if j != i:
   ....:             echange_lignes(M, i, j)
   ....:         if M[i][i] != 0:
   ....:             for j in range(i + 1, len(M)):
   ....:                 transvection_ligne(M, j, i, -M[j][i] / M[i][i])
   ....:     return M
   ....: 

Note

Le test if M[i][i] != 0:, s’il est correct en théorie, est en fait ridicule en pratique. Puisque l’on ne travaille qu’avec des valeurs approchées, un pivot nul en théorie (si l’on effectuait des calculs exacts) ne sera jamais nul en pratique.

In [11]: M = [[1, 2, 3, 4], [5, 6, 7, 8], [6, 8, 10, 12], [4, 4, 4, 4]]

In [12]: pivot_lignes(M)
Out[12]: 
[[1, 2, 3, 4],
 [0.0, -4.0, -8.0, -12.0],
 [0.0, 0.0, 0.0, 0.0],
 [0.0, 0.0, 0.0, 0.0]]

Note

On pourrait alors utiliser la forme échelonnée pour calculer le rang d’une matrice : il suffirait alors de compter le nombre de lignes non nulles. Mais à nouveau, il n’est pas évident de savoir en pratique si une ligne est réellement nulle puisqu’on a accès qu’à des valeurs approchées de ses coefficients.

8.5.4. Résolution de systèmes linéaires

On considère un système de Cramer sous forme matricielle \(AX=B\) où \(A\) est une matrice inversible, \(B\) une matrice colonne donnée et \(X\) une matrice colonne inconnue. Pour résoudre ce système, il suffit dans un premier temps de mettre la matrice \(\begin{pmatrix}A\mid B\end{pmatrix}\) sous forme échelonnée. On peut utiliser la fonction pivot_lignes précédemment définie mais on aura également besoin d’une fonction permettant de concaténer une matrice carrée (sous forme d’une liste de listes) et une matrice colonne (sous forme d’une liste).

In [13]: def concatenation_vecteur(A, B):
   ....:     return [A[i] + [B[i]] for i in range(len(A))]
   ....: 

Une fois que le pivot de Gauss a été effectué sur la matrice \(\begin{pmatrix}A\mid B\end{pmatrix}\), il faut effectuer un pivot « à rebours » pour déterminer la solution du système \(AX=B\).

In [14]: def pivot_lignes_rebours(M):
   ....:     for i in reversed(range(len(M))):
   ....:         dilatation_ligne(M, i, 1 / M[i][i])
   ....:         for j in range(i):
   ....:             transvection_ligne(M, j, i, -M[j][i])
   ....:     return M
   ....: 

La matrice colonne solution est alors la dernière colonne de la matrice obtenue, qu’il faut donc extraire.

In [15]: def extract_vecteur(M):
   ....:     return [L[-1] for L in M]
   ....: 

On peut alors définir une fonction d’arguments une matrice inversible \(A\) et une matrice colonne \(B\) renvoyant l’unique solution du système \(AX=B\).

In [16]: def resolution(A, B):
   ....:     M = concatenation_vecteur(A, B)
   ....:     pivot_lignes(M)
   ....:     pivot_lignes_rebours(M)
   ....:     return extract_vecteur(M)
   ....: 

In [17]: A = [[1, -1, 2], [3, 2, 1], [2, -3, -2]]

In [18]: B = [5, 10, -10]

In [19]: resolution(A, B)
Out[19]: [1.0, 2.0, 3.0]

8.5.5. Inversion d’une matrice

On peut également utiliser l’algorithme du pivot de Gauss pour inverser une matrice : on transforme une matrice inversible en la matrice identité en effectuant l’algorithme du pivot de Gauss puis l’algorithme du pivot de Gauss « à rebours ». On récpercute les opérations effectuées sur une matrice identité de même taille que \(A\), qui est alors transformée en l’inverse de la matrice initiale. Pour effectuer aissément les mêmes opérations sur les lignes d’une matrice \(A\) et la matrice identité \(I\), on forme la matrice \(\begin{pmatrix}A\mid I\end{pmatrix}\).

In [20]: def concat_identite(A):
   ....:     return [A[i] + [1 if j== i else 0 for j in range(len(A))] for i in range(len(A))]
   ....: 

Après les pivots, il reste à extraire la matrice inverse.

In [21]: def extract_inverse(M):
   ....:     return [L[len(M):] for L in M]
   ....: 

On peut alors proposer la fonction suivante.

In [22]: def inverse(A):
   ....:     M = concat_identite(A)
   ....:     pivot_lignes(M)
   ....:     pivot_lignes_rebours(M)
   ....:     return extract_inverse(M)
   ....: 

In [23]: A = [[1, 5, 6], [2, 11, 19], [3, 19, 47]]

In [24]: B = inverse(A)

In [25]: B
Out[25]: [[156.0, -121.0, 29.0], [-37.0, 29.0, -7.0], [5.0, -4.0, 1.0]]

In [26]: produit(A, B)
Out[26]: [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]

In [27]: produit(B, A)
Out[27]: [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]

8.5.6. Calcul du déterminant

On peut également se servir du pivot de Gauss pour calculer le déterminant d’une matrice carrée. En effet, le déterminant est invariant par transvection et échange de lignes et le déterminant d’une matrice triangulaire est le produit de ses coefficients diagonaux [2].

In [28]: def determinant(M):
   ....:     pivot_lignes(M)
   ....:     p = 1
   ....:     for i in range(len(M)):
   ....:         p *= M[i][i]
   ....:     return p
   ....: 

In [29]: M = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]

In [30]: determinant(M)
Out[30]: -0.0

[1]

Le module numpy possède un type matrix permettant de simplifier grandement les fonctions suivantes. Il possède d’ailleurs également un sous module numpy.linalg regroupant de nombreuses fonctions ayant trait à l’algèbre linéaire sur les matrices.

À faire

mettre un lien vers chapitre numpy ?

[2]

On pourrait penser à calculer le déterminant via la formule qui l’exprime en fonction des coefficients de la matrice ou à l’aide d’un développement par rapport à une ligne ou une colonne mais on verra dans le chapitre ??? que c’est nettemment moins efficace que le pivot de Gauss

À faire

Mettre référence vers le chapitre complexité