numpy

La bibliothèque numpy

Une des meilleures bibliothèques python pour le calcul numérique (donc l'algèbre linéaire) est numpy, que l'on peut installer ainsi sous une debian :

  $ aptitude install python-numpy

dans un shell, en étant administrateur.

Création d'une matrice

Une fois python lancé, voici comment importer la bibliothèque numpy, et créer une première matrice :

  >>> from numpy import *
  >>> A = matrix((1,2,3,4,5,6,7,8,9)).reshape((3,3))
  >>> A
  matrix([ [1, 2, 3],
          [4, 5, 6],
          [7, 8, 9] ])

Dans ce qui précède, on crée une matrice ligne de 9 éléments, que l'on redimensionne en $3 \times 3$ (par la méthode reshape). Rappelons qu'il existe d'autres manières de créer une matrice.

Enfin, pour accéder à un élément, on procède ainsi :

  >>> A[2,2] = 1
  >>> A
  matrix([ [1, 2, 3],
           [4, 5, 6],
           [7, 8, 1] ])

(Attention : l'élément supérieur gauche est en position (0,0)).

Matrices particulières

Dans ce qui suit, on crée une matrice nulle, une matrice identité, une autre constituée de uns, et une dernière remplie aléatoirement. Les arguments correspondent à la taille de la matrice obtenue.

  >>> zeros((3,3))
  array([ [ 0.,  0.,  0.],
          [ 0.,  0.,  0.],
          [ 0.,  0.,  0.] ])
  >>> eye(3)
  array([ [ 1.,  0.,  0.],
          [ 0.,  1.,  0.],
          [ 0.,  0.,  1.] ])
  >>> ones((2,4))
  array([ [ 1.,  1.,  1.,  1.],
          [ 1.,  1.,  1.,  1.] ])
  >>> random.random((4,4))
  array([ [ 0.47238805,  0.81230427,  0.30001673,  0.66223623],
          [ 0.86574785,  0.61391398,  0.67799879,  0.05343759],
          [ 0.56630056,  0.3685825 ,  0.84244261,  0.99121238],
          [ 0.85312013,  0.41931625,  0.10471793,  0.41117828] ])

Opérations matricielles

Considérons les deux matrices suivantes :

  >>> A = matrix((1,2,3,4,5,6,7,8,9)).reshape((3,3))
  >>> B = matrix((1,2,1,2,1,2,1,1,1)).reshape((3,3))
  >>> A
  matrix([ [1, 2, 3],
           [4, 5, 6],
           [7, 8, 9] ])

On a accès aux opérations matricielles élémentaires : somme et multiplication matricielle, multiplication par un scalaire...

  >>> A+B
  matrix([ [ 2,  4,  4],
           [ 6,  6,  8],
           [ 8,  9, 10] ])
  >>> 2*A
  matrix([ [ 2,  4,  6],
           [ 8, 10, 12],
           [14, 16, 18] ])
  >>> A*B
  matrix([ [ 8,  7,  8],
           [20, 19, 20],
           [32, 31, 32] ])

On a aussi la possibilité de transposer les matrices (symétrie par rapport à la diagonale descendante), ou d'en calculer sa trace (somme des éléments diagonaux) :

  >>> A.transpose()
  matrix([ [1, 4, 7],
          [2, 5, 8],
          [3, 6, 9] ])
  >>> A.trace()
  matrix([ [15] ])

Les fonctions importantes

Les calculs d'inverses, de déterminants, etc. sont accessibles avec le paquetage numpy.linalg :

  >>> from numpy.linalg import *
  >>> A = matrix((1,2,3,4,5,6,7,8,1)).reshape((3,3))
  >>> print det(A)
  24.0
  >>> inv(A)
  matrix([ [-1.79166667,  0.91666667, -0.125     ],
           [ 1.58333333, -0.83333333,  0.25      ],
           [-0.125     ,  0.25      , -0.125     ] ])
  >>> A*inv(A)
  matrix([ [  1.00000000e+00,   0.00000000e+00,   0.00000000e+00],
           [ -1.11022302e-16,   1.00000000e+00,   8.32667268e-17],
           [  3.88578059e-16,  -4.44089210e-16,   1.00000000e+00] ])

Dans ce qui précède, on a importé le paquetage linalg, affiché le déterminant et l'inverse de A, et calculé $ A \times A^{-1}$. Cependant, on aimerait obtenir un résultat plus propre, plus ressemblant à l'identité matricielle. Pour ce faire, on change les options d'affichage :

  >>> set_printoptions(suppress = True)
  >>> A*inv(A)
  matrix([ [ 1.,  0.,  0.],
           [-0.,  1.,  0.],
           [ 0., -0.,  1.] ])
  >>> print inv(A)
  [ [-1.79166667  0.91666667 -0.125     ]
    [ 1.58333333 -0.83333333  0.25      ]
    [-0.125       0.25       -0.125     ] ]

Produit Matriciel

Le produit classique

Présentation

Le produit de deux matrices ne peut se définir que si le nombre de colonnes de la première matrice est le même que le nombre de lignes de la deuxième matrice.

Si $A = (a_{ij})$ est une matrice de type $(m,n)$ et $B = (b_{ij})$ est une matrice de type $(n,p)$, alors leur produit, noté $AB = (c_{ij})$ est la matrice de type $(m,p)$ donnée par :

$c_{ij}= \sum_{k=1}^n a_{ik}\cdot b_{kj} = a_{i1} \cdot b_{1j} + a_{i2} \cdot b_{2j}+ \ldots + a_{in} \cdot b_{nj}$

pour chaque $1 \le i \le m$ et $1 \le j \le p$.

La figure suivante (extraite de http://fr.wikipedia.org/wiki/Produit_matriciel) montre comment calculer les coefficients $c_{12}$ et $c_{33}$ de la matrice produit $AB$ si $A$ est une matrice de type (4,2), et $B$ est une matrice de type (2,3).

Par exemple :

• ${c_{12}} = \sum_{r=1}^2 a_{1r}b_{r2} = a_{11}b_{12}+a_{12}b_{22}$,
• ${c_{33}} = \sum_{r=1}^2 a_{3r}b_{r3} = a_{31}b_{13}+a_{32}b_{23}$.

Un mot de complexité

La complexité de la multiplication naïve est en $O(n^3)$ :

• Pour chaque $c_{i,j}$, on effectue $n$ multiplications et $n-1$ additions.
• On a donc $2n-1$ opérations par coefficient, et $n^2$ coefficients.
• On trouve au final $2 n^3 - n^2$ opérations pour la multiplication, qui est du même ordre de grandeur que $n^3$ : on note cela $O(n^3)$.

Notons, pour finir cette section, que les additions sont généralement ignorées dans le calcul de la complexité : elles sont beaucoup plus rapides que la multiplication, en particulier si la taille des entrées est supérieure à la taille du mot machine.

Travaux pratiques

1. Importer le module python numpy.
2. (Re)découvrir la programmation orientée objets en python : ici.
3. Programmer une classe 'Matrice' munie d'une méthode 'produit'. Celle-ci implante la multiplication entre la matrice courante et la matrice donnée en argument de la méthode, si les dimensions coïncident.
4. Appliquer la méthode à des exemples de taille croissante que vous laisserez dans votre code.

L'algorithme de Strassen

Présentation

Cet algorithme, dû à Volker Strassen (1969), permet d'obtenir le produit de deux matrices de taille $n$ en $O(n^{2,807})$. Ce résultat prouve que la multiplication naïve, et donc le pivot de Gauss, ne sont pas optimaux.

Description de l'algorithme

Soient $A$ et $B$ deux matrices carrées de taille $n$, dont on souhaite calculer le produit $C = A \times B$. On suppose, pour se simplifier la vie, que $n$ est de la forme une puissance de 2 : si tel n'est pas le cas, il suffit de compléter les matrices avec des 0.

Les trois matrices A, B et C sont divisées en matrices par blocs de taille égale :

$ \mathbf{A} = \left(\begin{array}{cc} \mathbf{A}_{1,1} & \mathbf{A}_{1,2} \\ \mathbf{A}_{2,1} & \mathbf{A}_{2,2} \end{array}\right) \mbox { , } \mathbf{B} = \left(\begin{array}{cc} \mathbf{B}_{1,1} & \mathbf{B}_{1,2} \\ \mathbf{B}_{2,1} & \mathbf{B}_{2,2} \end{array}\right) \mbox { , } \mathbf{C} = \left(\begin{array}{cc} \mathbf{C}_{1,1} & \mathbf{C}_{1,2} \\ \mathbf{C}_{2,1} & \mathbf{C}_{2,2} \end{array}\right)$

où les $\mathbf{A}_{i,j}, \mathbf{B}_{i,j}$, et $\mathbf{C}_{i,j}$ ont une taille $\frac{n}{2}$. On a alors, en identifiant $A\times B$ et $C$ :

• $\mathbf{C}_{1,1} = \mathbf{A}_{1,1} \times \mathbf{B}_{1,1} + \mathbf{A}_{1,2} \times \mathbf{B}_{2,1}$
• $\mathbf{C}_{1,2} = \mathbf{A}_{1,1} \times \mathbf{B}_{1,2} + \mathbf{A}_{1,2} \times \mathbf{B}_{2,2}$
• $\mathbf{C}_{2,1} = \mathbf{A}_{2,1} \times \mathbf{B}_{1,1} + \mathbf{A}_{2,2} \times \mathbf{B}_{2,1}$
• $\mathbf{C}_{2,2} = \mathbf{A}_{2,1} \times \mathbf{B}_{1,2} + \mathbf{A}_{2,2} \times \mathbf{B}_{2,2}$

On constate que cette méthode nécessite 8 multiplications de matrices pour calculer les $C_{i,j}$, et ce comme dans le produit classique.

L'idée de Strassen est de réussir à obtenir ces $C_{i,j}$ avec seulement 7 multiplications, en introduisant les matrices suivantes :

• $\mathbf{M}_{1} = (\mathbf{A}_{1,1} + \mathbf{A}_{2,2}) \times (\mathbf{B}_{1,1} + \mathbf{B}_{2,2})$
• $\mathbf{M}_{2} = (\mathbf{A}_{2,1} + \mathbf{A}_{2,2}) \times \mathbf{B}_{1,1}$
• $\mathbf{M}_{3} = \mathbf{A}_{1,1} \times (\mathbf{B}_{1,2} - \mathbf{B}_{2,2})$
• $\mathbf{M}_{4} = \mathbf{A}_{2,2} \times (\mathbf{B}_{2,1} - \mathbf{B}_{1,1})$
• $\mathbf{M}_{5} = (\mathbf{A}_{1,1} + \mathbf{A}_{1,2})\times \mathbf{B}_{2,2}$
• $\mathbf{M}_{6} = (\mathbf{A}_{2,1} - \mathbf{A}_{1,1})\times (\mathbf{B}_{1,1} + \mathbf{B}_{1,2})$
• $\mathbf{M}_{7} = (\mathbf{A}_{1,2} - \mathbf{A}_{2,2})\times (\mathbf{B}_{2,1} + \mathbf{B}_{2,2})$

Les $C_{i,j}$ sont alors exprimées comme

• $\mathbf{C}_{1,1} = \mathbf{M}_{1} + \mathbf{M}_{4} - \mathbf{M}_{5} + \mathbf{M}_{7}$
• $\mathbf{C}_{1,2} = \mathbf{M}_{3} + \mathbf{M}_{5}$
• $\mathbf{C}_{2,1} = \mathbf{M}_{2} + \mathbf{M}_{4}$
• $\mathbf{C}_{2,2} = \mathbf{M}_{1} - \mathbf{M}_{2} + \mathbf{M}_{3} + \mathbf{M}_{6}$

On fait plus d'additions et de soustractions matricielles (qui ont une complexité linéaire) que dans la méthode naïve, mais une multiplication de moins : le gain est énorme.

Le procédé est répété jusqu'à ce que les matrices $A$ et $B$ soient tailles raisonnables.

Complexité

Le nombre de multiplications $T(n)$ est ici de l'ordre de $n^{\log_{2}7}\approx n^{2.807}$, solution de l'équation par récurrence $T(n) = 7T\left(\frac{n}{2}\right) + O(n^2)$ : pour faire une multiplication de taille $n$, on fait 7 multiplications de taille $\frac{n}{2}$ et quelques opérations annexes.

Discussion

La constante dans la complexité est proche de 50 : l'algorithme de Strassen n'est efficace que pour des matrices de grandes tailles.

De plus, cet algorithme est peu efficace sur les matrices creuses (qui ont peu de coefficients non-nuls), courantes en analyse numérique.

Enfin, l'algorithme n'est pas très stable numériquement.

Importance du résultat

Le travail de Strassen a permis de développer un nouveau champ de la recherche en informatique théorique : le calcul rapide du produit matriciel.

Plusieurs résultats importants suivront : Samuel Winograd publie un algorithme en 1980 qui utilise moins d'additions que celui de Strassen. Et, surtout, Don Coppersmith et Winograd publient en 1987 l'algorithme du même nom, améliorant grandement la complexité de la multiplication (en $O(n^{2,376})$).

Ce dernier algorithme est le plus efficace en théorie ; cependant, rien n'indique qu'il est optimal, en terme de complexité : on pense pouvoir atteindre l'exposant 2.

D'un point de vue théorique, donc, cet algorithme est d'importance, mais il n'a aucun intérêt pratique : la constante dans le grand $O$ est beaucoup trop grande.

Travaux pratiques

1. Mettre en place l'algorithme de Strassen, pour améliorer votre méthode de multiplication matricielle.
2. Comparer les temps de calcul pour des matrices d'une certaine taille. On pourra utiliser PyX pour illustrer cela.

Résolution de systèmes linéaires

Résolution avec numpy

On montre, pour commencer, comment faire pour résoudre numériquement un système linéaire, en utilisant numpy.

On commence par définir la matrice $A$ et le vecteur $b$ du système $AX=b$ qui nous intéresse :

  >>> from numpy import *
  >>> A = matrix([[1,2,3],[4,5,6],[7,8,9]])
  >>> A
  matrix([[1, 2, 3],
          [4, 5, 6],
          [7, 8, 9]])
  >>> b=matrix((1,2,3)).reshape((3,1))
  >>> b
  matrix([[1],
          [2],
          [3]])

On importe ensuite le paquetage numpy.linalg, qui contient la fonction solve permettant cette résolution, et on calcule $X$ :

  >>> from numpy.linalg import *
  >>> X=solve(A,b)
  >>> X
  matrix([[-0.26193963],
          [ 0.52387926],
          [ 0.0713937 ]])

Pour finir, on vérifie que le produit $AX$ redonne bien $b$...

  >>> A*X
  matrix([[ 1.],
          [ 2.],
          [ 3.]])

Méthode du pivot de Gauss

Présentation

Pour résoudre le système $Ax = b$, où $A$ est une matrice quelconque, on cherche à se ramener à un système plus simple que l'on sait résoudre : un système triangulaire (méthode du pivot de Gauss) ou un système diagonal (méthode de Gauss-Jordan) Pour cela on transforme un système en un système équivalent obtenu par :

• échange de deux lignes,
• multiplication d'une ligne par un nombre non nul,
• addition d'un multiple d'une ligne à une autre ligne.

Exemple

Soit à résoudre le système suivant de trois équations à trois inconnues

$ \left\{ \begin{array}{rrrrrrrr} x & +& 2y& + & 2z& = & 2 &L_1\\ x & +& 3y& - & 2 z & = & -1&L_2\\ 3x & +& 5y & + & 8z & = & 8 &L_3\end{array} \right.$

Le premier pivot est $a_{11} = 1 $. On supprime l'inconnue $x$ de $L_2$ et $L_3$ en utilisant la transformation :

• $L'_2 \leftarrow L_2 - \frac{a_{21}}{a_{11}} \cdot L_1$
• $L'_3 \leftarrow L_3 - \frac{a_{31}}{a_{11}} \cdot L_1$

On a :

$ \left\{ \begin{array}{rrrrrrrrr} x & +& 2y& + & 2z &= & 2 & & &L'_1\\ & & y& - & 4z &= &-3 & & & L'_2 \leftarrow L_2 - \frac{a_{21}}{a_{11}} \cdot L_1\\ & & -y& + & 2z &= &2 & & & L'_3 \leftarrow L_3 - \frac{a_{31}}{a_{11}} \cdot L_1\end{array} \right.$

Le second pivot est $a'_{22} = 1 $.

On supprime l'inconnue $y$ de $L_3$ en utilisant la transformation :

• $L"_3 = L'_3 - \frac{a'_{32}}{a'_{22}} \cdot L'_2$

Ce qui donne, cette fois-ci...

$ \left\{ \begin{array}{rrrrrrrrrr} x & +& 2y& + & 2z &= & 2 & & & L"_1\\ & & y& - & 4z &= &-3 & & & L"_2 \\ & & & & -2z &= &-1 & & & L"_3 = L'_3 - \frac{a'_{32}}{a'_{22}} \cdot L'_2 \end{array} \right.$

Une fois la valeur de $z$ obtenue, il reste à remonter l'information pour obtenir les valeurs des variables $y$ puis $x$.

Formalisation

La méthode de Gauss est itérative : on considère les suites $(A^n)$ et $(b^n)$ de matrices telle que $A^0 = A$ et $b^0 = b $ (attention à la notation utilisée : $X^n$ n'est pas une n^ième puissance, mais le n^ième terme d'une suite).

On construit successivement la séquence de matrices à $n$ ligne et $n+1$ colonnes $(A^0, b^0)$, $(A^1, b^1)$, ..., $(A^n, b^n)$ telle qu'à l'itération $i$, la matrice $(A^i,b^i)$ a la forme

$\left(\begin{array}{rrrr|rrrr|r} a_{1,1}^i & a_{1,2}^i & \ldots & a_{1,i}^i & a_{1,i+1}^i & \ldots & a_{1,n-1}^i & a_{1,n}^i & b_1^i\\ 0 & a_{2,2}^i & \ldots & a_{2,i}^i & a_{2,i+1}^i & \ldots & a_{2,n-1}^i & a_{2,n}^i &b_2^i \\ \vdots & \ddots & \ddots & \vdots & & \vdots & \vdots & \vdots & \vdots \\ 0 & \ldots & 0 & a_{i,i}^i & a_{i,i+1}^i & \ldots & a_{i,n-1}^i & a_{i,n}^i & b_n^i \\0 & 0 & \ldots & O & a_{i+1,i+1}^i& \ldots & a_{i+1,n-1}^i & a_{i+1,n}^i & b_{n+1}^i\\ \vdots & \vdots & & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \ldots & O & a_{n-1,i+1}^i& \ldots & a_{n-1,n-1}^i & a_{n-1,n}^i & b_{n-1}^i\\ 0 & 0 & \ldots & O & a_{n,i+1}^i & \ldots & a_{n,n-1}^i & a_{n,n}^i & b_{n}^i\\ \end{array}\right)$

On obtient la matrice $A^i$ en faisant apparaître des 0 à la colonne $i+1$. Pour cela chaque ligne $L_j$, $j > i+1$, est remplacée par une combinaison linéaire d'elle même avec la ligne $L_{i+1}$ comme suit :

$L_j = L_j - \frac{a_{i+1,j}^i}{a_{i+1,i+1}^i}\cdot L_{i+1}$ si $a_{i+1,i+1}^i$ est non nul.

Travaux pratiques

1. Programmez l'algorithme du pivot de Gauss.
2. Retrouvez les résultats de l'exemple ci-dessus.

Méthode de Gauss-Jordan

Présentation

C'est une variante de la méthode du pivot de Gauss : à la l'étape $i$, on combine toutes les lignes (sauf la ligne i) avec la ligne $i$ (au lieu de ne le faire que pour les lignes d'indice supérieurs à $i$). L'intéret est que cela fait apparaître des 0 sur toute la colonne sauf au niveau du pivot $a_{ii}^i$

Exemple

Soit à résoudre le système de trois équations à trois inconnues suivant :

$ \left\{ \begin{array}{rrrrrrrl} 2x & +& y& - & 4z& = & 8 &L_1\\ 3x & +& 3y& - & 5z & = & 14&L_2\\ 4x & +& 5y & - & 2z & = & 16 &L_3\end{array} \right.$

1. on choisit comme premier pivot le nombre $a_{11} = 2 $
2. on normalise la ligne $L_1$ en la divisant par 2 et on obtient $L'_1$
3. on applique les transformations suivantes pour les autres lignes :
   • $L'_2 \leftarrow L_2 -3 \cdot L'_1$
   • $L'_3 \leftarrow L_3 -4 \cdot L'_1$

On trouve alors...

$ \left\{ \begin{array}{rrrrrrrl} x & +& 1/2y& - & 2z &= & 4 & L'_1 \leftarrow L_1/2\\ & & 3/2 \cdot y& + & z &= &2 & L'_2 \leftarrow L_2 -3 \cdot L'_1\\ & & 3y& + & 6z &= &0 & L'_3 \leftarrow L_3 -4 \cdot L'_1 \end{array}\right.$

1. on choisit comme second pivot le nombre $a_{22} = 3/2 $
2. on normalise la ligne $L'_2$ en la multipliant par 2/3 et on obtient $L"_2$
3. on applique les transformations suivantes pour les autres lignes :
   • $L"_1 \leftarrow L'_1 -1/2 \cdot L"_2$
   • $L"_3 \leftarrow L'_3 -3 \cdot L"_2$

Cela donne :

$ \left\{ \begin{array}{rrrrrrrl} x & & & & -7/3 \cdot z &= & 10/3 & L"_1 \leftarrow L'_1 -1/2 \cdot L"_2 \\ & & y& + & 4/3 \cdot z &= &4/3 & L"_2 \leftarrow 2/3 \cdot L'_2 \\ & & & & 4z &= &-4 & L"_3 \leftarrow L'_3 -3 \cdot L"_2 \end{array}\right.$

1. on choisit enfin comme troisième pivot le nombre $a_{33} = 4 $
2. on normalise la ligne $L"_3$ en la divisant par 4 et on obtient $L_3^{3}$
3. on applique les transformations suivantes pour les autres lignes :
   • $L_1^{3} \leftarrow L"_1 + 7/3 \cdot L_2^{3}$
   • $L_2^{3} \leftarrow L"_2 - 2/3 \cdot L_2^{3}$

D'où, finalement,

$ \left\{ \begin{array}{rrrrrrrl} x & & & & & = & 1 & L_1^{3} \leftarrow L"_1 + 7/3 \cdot L_2^{3}\\ & & y& & &= &2 & L_2^{3} \leftarrow L"_2 - 2/3 \cdot L_2^{3} \\ & & & & z &= &-1 & L_3^{3} \leftarrow L"_3/4 \cdot L"_2 \end{array}\right.$

et on trouve directement les racines sur la 4^ème colonne.

Algorithme

On a le pseudo algorithme suivant :

  pour i = 1 à n 
    recherche du pivot (non nul) sur la ligne p, p >= i 
    échange éventuel des lignes i et p -> le pivot est Aii
    division de la ligne Li par Aii pour obtenir L'i
    pour k = 1 à n sauf i
      L'k = L_k - Aki .L'i

Les solutions sont dans la colonne $n+1$.

Travaux Pratiques.

• Implanter l'algorithme de Gauss-Jordan et s'assurer de son fonctionnement sur les exemples du TP.
• Résoudre avec l'algorithme implanté les deux systèmes équivalents

$\left(\begin{array}{cc}10^{-12} & 1 \\ 1 & 1 \end{array}\right) \cdot \left(\begin{array}{c}x \\y \end{array}\right)=\left(\begin{array}{c} 1 \\2 \end{array}\right)$ et $\left(\begin{array}{cc}1 & 1 \\ 10^{-12} & 1 \end{array}\right) \cdot \left(\begin{array}{c}x \\y \end{array}\right)=\left(\begin{array}{c} 1 \\2 \end{array}\right)$.

• On cherche à résoudre le système $Ax = b$ où $A$ est donnée par

$ \left(\begin{array}{cccc}10 & 7 & 8 & 7 \\7 & 5 & 6 & 5 \\9 & 6 & 10& 9 \\ 7 & 5 & 9 & 10 \end{array} \right)$

$\left( \begin{array}{c} 32 \\ 23 \\ 33 \\ 31 \end{array} \right)$ $\left( \begin{array}{c} 32,1 \\ 22,9 \\ 33,1 \\ 30,9 \end{array} \right)$.

Calcul de déterminants

Présentation

D'après wikipedia...

En mathématiques, initialement introduit en algèbre pour déterminer le nombre de solutions d'un système d'équations linéaires, le déterminant se révèle un outil très puissant dans de nombreux domaines (étude d'endomorphisme, recherche de valeurs propres, calcul différentiel).

C'est ainsi qu'on définit le déterminant d'un système d'équations, le déterminant d'un endomorphisme ou le déterminant d'un système de vecteurs.

Méthode de calcul

On a déjà vu quelle était la fonction numpy qui permettait le calcul du déterminant.

En pratique, on peut procéder à un algorithme du type pivot, pour obtenir une forme triangulaire de la matrice $A$ qui nous intéresse. Le déterminant de $A$ est alors égal au produit des éléments diagonaux. Cette méthode est a peu près la seule à n'être pas prohibitif en temps de calcul.

Le seul problème arrive si la matrice n'est constituée que d'entiers : le pivot introduit artificiellement des fractions, ce qui rallonge les calculs et peut être source d'erreurs. Cela est d'autant plus paradoxal qu'il existe une formule, inutile en pratique, affirmant que le déterminant s'exprime comme sommes et produits des coefficients de la matrice : tout devrait se passer chez les entiers.

Il existe cependant une variante à la méthode du pivot pour le calcul du déterminant, utile dans le cas des matrices à coefficients entiers : la méthode de Jordan-Bareiss. Dans ce cas, on ne travaille qu'avec des entiers :

• on est ainsi à l'abri des erreurs d'arrondis,
• l'algorithme devrait être plus rapide.

En python

La bibliothèque sympy a implanté l'algorithme de Bareis :

  >>> from sympy import *
  >>> A = Matrix((-5,2,4,-4,
                  10,-3,-6,8,
                  11,-4,-6,8,22,
                  -8,-14,17)).reshape(4,4)
  >>> A
  [-5,  2,   4, -4]
  [10, -3,  -6,  8]
  [11, -4,  -6,  8]
  [22, -8, -14, 17]
  >>> A.det_bareis()
  -2
  >>> help(A.det_bareis)

Algorithme de Jordan-Bareiss

L'algorithme

Soit $M$ une matrice carrée de taille $n$, à coefficients entiers, dont on souhaite calculer le déterminant.

On suppose que les mineurs principaux $x_{k,k}^k$ de $M$ sont tous non nuls (le mineur principal d'ordre k de A est le déterminant obtenu à partir des k premières lignes et des k premières colonnes de A).

  $M_{0,0} \leftarrow 1$
  pour k de 0 à n-1 faire
      pour i de k+1 à n faire
          pour j de k+1 à n faire
              $M_{i,j} = \frac{M_{i,j}M_{k,k}-M_{i,k}M_{k,j}}{M_{k-1,k-1}}$

La division ci-dessus est exacte : $M_{i,j}$ sera entier. Le déterminant est égal à $M_{n-1,n-1}$.

Remarquons, pour bien comprendre l'algorithme, et bien pouvoir le mettre en oeuvre, que les coefficients de $M$ vont de $(1,1)$ à $(n,n)$. Donc...

• Le $M_{0,0}$ est donc en-dehors de la matrice, et ne sert qu'à éviter un problème pour la dernière ligne de l'algorithme,
• On prendra garde à ce que, en python :
  • le range(x,y) va de x à y-1,
  • les coefficients des matrices ne vont pas de 1 à $n$, mais de 0 à $n-1$ : l'algorithme est donc à adapter légèrement.

Travaux pratiques

1. Mettre en place un algorithme de calcul de déterminant pour une matrice quelconque.
2. Implantez l'algorithme de Jordan-Bareiss.
3. Vérifiez vos résultats en utilisant la fonction det de numpy.