Algorithme d'inversion d'une matrice carrée par opérations élémentaires sur les lignes

Bonjour à tous,

Je souhaite écrire une implémentation de l’algorithme qui inverse une matrice en procédant à une succession d’opérations sur les lignes. C’est quelque chose que j’ai déjà vu une fois à la fac, mais la prof ne nous avait pas demandé d’en faire un programme et on a vu ça vite fait.

Pour ce faire, je me suis basé sur ce cours : https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/matrix-inverse .

Résumons grossièrement

On dispose de trois opérations de base

1. Échanger deux lignes

2. Multiplier chaque coefficient d’une ligne par un nombre (différent de 0)

3. Remplacer une ligne par la somme d’elle-même et le multiple d’une autre ligne

Idée générale

L’idée est d’appliquer une succession de ces trois opérations de base sur notre matrice N°1 ; ces applications, comment dire, "personnalisées" pour cette matrice. En parallèle, on applique exactement ces mêmes opérations sur une matrice identité qu’on s’est préalablement donnée (matrice N°2). A la toute fin du processus, la matrice N°1 doit être devenue une matrice identité, et notre matrice d’identité (N°2) doit être devenue une matrice non-identité et il s’agit de la matrice inverse de la matrice N°1.

Le processus

Reste à spécifier comment sont utilisées ces opérations de base (j’ai déjà précisé que cela dépendait de la matrice (N°1) traitée, mais ce n’est bien sûr pas une indication suffisante pour savoir comment inverser une matrice). Le processus se déroule en trois étapes.

Les explications ci-dessous proviennent du cours mentionné ci-dessus, et nécessitent de connaître la notion de pivot de matrice (définie dans ce même cours).

Etape N°1

Opération de base utilisée : échanges de deux lignes successifs,

The first thing we will do is to examine the value for the pivot coefficient for the current column (the column being processed). In order for our technique to work, this value has to be different than zero. If it is different than zero then we can directly move to step two. If it is equal to zero, we will need to find another row in the matrix for this column which value is different than zero, and swap these two rows (operation 1). However more than one row may have values different than zero ? So which one should we select ? We will pickup the row which absolute value is the highest.

If we can’t find another value for the pivot coefficient (that is if all the other values in the column are zero) then the matrix can’t be inverted and is singular

Le code source est fourni sur le cours.

Etape 2 : j’ai besoin d’aide !

L’idée est apparemment, grâce à des opérations de base, de mettre à 0 toutes les valeurs de chaque colonne (excepté leur pivot) des matrices N°1 et identité (N°2).

SAUF QUE j’ai beau lire et relire le cours, je ne comprends absolument pas ce qu’ils expliquent. Je trouve ça très flou, et l’exemple qu’ils donnent semble être partiel donc ça n’aide pas…

If the current column we are processing has a valid pivot coefficient, then we can proceed to the next step. We will reduce each coefficient of the column (apart from the pivot coefficient which we won’t touch for now) to zero. In English that won’t make much sense but here it is. Lets call ’m’, the current cofficient (for the current column/row) that we want to reduce to zero. We will subtract m muliplied by each coefficient of the row which index is the same as the index of the current column, from each coefficients of the current row (operation 2 and 3). The technique will work if m itself is divided by the pivot coefficient value (line 4). This is illustrated for one coefficient in the following example. The second column is the current column. We want to reduce the first coefficient (row 0) from the current column (4.3) to zero. For each coefficient of the first row (say C0j), we take the corresponding coefficient (C1j) from the second row (because we are processing the second column), multiply this coefficient by the coefficient we want to reduce by zero (C02) divided by the pivot coefficient from the second column (C12) and subtract the resulting number from C0j.

Bon et ils donnent le code source mais voilà quoi… Il n’a pas l’air très compliqué mais… J’aimerais plutôt comprendre sans devoir le lire…

Notez que je ne demande pas d’explications du genre "Mais pourquoi cette succession précise d’opérations de base permet-elle de mettre toutes ces valeurs à 0 ? Ce n’est pas magique je suppose ?!". Cela fera effectivement l’objet d’un second topic (le cours ne l’explique en effet pas).

Etape N°3

Convertir tous les pivots à 1. Bon là je pense comprendre. Je n’ai pas cherché à lire plus que ça cette étape du cours, car j’aimerais déjà comprendre l’étape 2 mais ça semble parler de diviser la valeur de chaque colonne par le pivot de cette dernière, et sûrement d’arrondir donc ça va, ça me semble compréhensible et aisément réalisable.

Ce que j’ai actuellement écrit

Etape 1 qui fonctionne

Le code qui suit est une simple implémentation de ce que j’ai cité et expliqué précédemment (cf. partie "Etape 1").

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

  def getGaussJordanInvertedMatrix: (Matrix, Matrix) = {
    var identity_matrix : mutable.Seq[Seq[Double]] = getIdentityMatrix(content.length)  // We get the identity matrix. It will be modified
                                                                      // as the original matrix will.

    var inverted_matrix_content : scala.collection.mutable.Seq[Seq[Double]] = scala.collection.mutable.Seq[Seq[Double]](content: _*)

    content.zipWithIndex.foreach({ case (sub_sequence, first_index) =>
      sub_sequence.indices.foreach(last_index => {  // These first three lines simply indicates if either or not the
                                                    // current case is a pivot (i.e. if it belongs to the diagonal)

        if(getDiagonalValuesCoordinates(content.length).contains((first_index, last_index))) {

          if(inverted_matrix_content(first_index)(last_index) == 0) {  // Swapping (on both original and identity
                                                                       // matrices) is done only if the pivot = 0

            val next_row_index = (first_index until content.length).maxBy(content(_)(last_index).abs)
            inverted_matrix_content = inverted_matrix_content.updated(next_row_index, inverted_matrix_content(first_index)).updated(first_index, inverted_matrix_content(next_row_index))
            identity_matrix = identity_matrix.updated(next_row_index, identity_matrix(first_index)).updated(first_index, identity_matrix(next_row_index))

          }

        }
      })
    })

    (new Matrix(identity_matrix), new Matrix(inverted_matrix_content))
  }

Etape 2 : pouvez-vous me réexpliquer cela s’il vous plaît ?

Comme je l’ai déjà dit, le cours est pour moi assez incompréhensible au niveau de cette étape, qu’il présente en outre comme un truc magique et compliqué… Y a moyen de reformuler ça svp ?

J’en ferai une implémentation selon VOS explications (limitez-vous à l’étape 2 svp, car le cours me paraît compréhensible concernant l’étape 3) et je vous la montrerai !

Merci d’avance !

Salut,

J’ai rien panné au texte, mais à tous les coups c’est la méthode de Gauss en fait. Prenons une matrice $3\times3$ pour l’exemple dont je note les éléments $a_{ij}$ et les différents lignes $L_i$. Il faut tout d’abord échanger les lignes de sorte à ce qu’il n’y ait plus de $0$ sur la diagonale (si ce n’est pas possible, c’est que la matrice n’est pas inversible).

Une fois que c’est fait, faire l’opération suivante : $L_2\leftarrow L_2 - \dfrac{a_{21}}{a_{11}}L_1$. Une fois réalisée, tu auras $a_{21}=0$. Ensuite, faire l’opération $L_3\leftarrow L_3-\dfrac{a_{31}}{a_{11}}L_1-\dfrac{a_{32}}{a_{22}}L_2$. Tu auras alors $a_{31}=a_{32}=0$ (sur une matrice plus grande, il faut répéter ces opérations sur les lignes suivantes jusqu’à obtenir une matrice diagonale supérieure).

Enfin, il suffit de faire la même chose que précédemment mais sur la partie supérieure de la matrice : $L_2\leftarrow L_2-\dfrac{a_{23}}{a_{33}}L_3$ (ce qui conduit à $a_{23}=0$) puis $L_1\leftarrow L_1-\dfrac{a_{12}}{a_{22}}L_2-\dfrac{a_{13}}{a_{33}}L_3$.

À la fin du processus, ta matrice est diagonale.

Reagarde sur Wikipedia : Élimination de Gauss-Jordan. Il y a des explications et même du pseudo-code.
Cette méthode est quand même moins pénible que la méthode de Cramer.

Non, juste de la matrice N°1.

Bon, oublie temporairement tout ce que tu as lu (j’insiste, il y a plusieurs manières de faire et ça serait pas une bonne idée de vouloir à tout prix refaire ce qu’ils expliquent). Je te donne 3 lignes de 3 nombres, tu as le droit aux opérations suivantes :

(i) Échanger deux lignes.
(ii) Choisir un nombre $x$ et une ligne, et multiplier tous les nombres de la ligne par $x$.
(iii) Choisir une ligne source, une liste destination et un nombre $x$, et ajouter à chaque nombre de la ligne destination $x$ fois le nombre sur la même colonne et la ligne source.

À partir de $\begin{bmatrix}7&2&3\\4&5&6\\2&2&2\end{bmatrix}$, tu as le droit d’appliquer un nombre quelconque de ces opérations. Ton objectif est d’obtenir la matrice $\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix}$. Comment tu procèdes (en pseudo-code ou en français) ?

Merci pour vos messages et excusez-moi de répondre avec du retard, je manque un peu de temps libre en ce moment.

@Lucas-84 je pense que je peux diviser la ligne n°1 par 7, la ligne n°2 par 5, la ligne n°3 par 2. Une fois ces divisions faites, j’aurais des 1 en diagonale, comme sur l’identité. Puis je soustrais la ligne n°2 ainsi obtenue par 4/5 x la ligne n°3, de sorte à avoir un 0 en (2, 1) (2è ligne, 1ère colonne). Puis je m’occupe de la même façon des autres cases où il faut mettre un 0 à la place du coefficient qu’elles contiennent, en ayant en tête que la modification d’une ligne pour mettre un 0 aura une répercussion sur les autres opérations à effectuer sur cette même ligne pour mettre les 0 sur les cases != 0 restantes (e.g. ayant mis un 0 en (2,1) grâce à une soustraction avec 4/5 x ligne N°3, je devrais tenir compte de cette opération pour mettre un 0 en (2, 3) ).

Pseudo-code inspiré de celui présenté dans Wikipédia

1. Dans un premier temps, je m’assure que la diagonale ne contient que des coefficients != 0. Si je rencontre un coefficient == 0 en diagonale, j’échange la ligne qui le contient avec la ligne contenant un coefficient != 0 et qui soit LE coefficient le plus grand de toutes les lignes.

2. A cette étape, ma diagonale contient bien des pivots != 0.

3. Je divise chaque ligne par son pivot, de sorte que la diagonale ne contiennent que des 1.

4. Je me débrouille pour obtenir des 0 partout, ailleurs bien entendu que sur la diagonale. Volontairement, je ne détaille pas trop ce point (j’y reviendrai plus tard).

Est-ce que ça vous paraît bon ? Pourquoi doit-on choisir le coefficient maximal de toutes les lignes, si un pivot == 0 (étape 1) ??? Tant que le pivot est != 0, ce serait bon non ?

Edit important

En fait je me suis trompé : les divisions et échanges de lignes pour obtenir les 1 en diagonale doivent se faire après les opérations apportant des 0, car à l’inverse on virerait les 1 …

Pourquoi doit-on choisir le coefficient maximal de toutes les lignes, si un pivot == 0 (étape 1) ??? Tant que le pivot est != 0, ce serait bon non ?

Ce sont de simple questions d’algèbre en virgule flottante. Si on travaille avec des nombres réels tels que définis mathématiquement, ça n’a effectivement aucune importance. le problème, c’est qu’un ordinateur travaille avec des nombres qui ont une précision limitée, si tu utilises un petit nombre comme pivot, tu vas te retrouver avec des nombres gigantesques sur ta ligne après normalisation par le pivot, et tu vas te prendre des erreurs d’arrondis dans les dents lors des soustractions des lignes les unes avec les autres. De manière générale, des problèmes similaires peuvent se poser si tu as de forts contrastes dans les valeurs présentes dans ta matrice, la notion technique derrière est ce qu’on appelle le conditionnement de la matrice que tu cherches à inverser.

D’ailleurs, placer les nombres maximums sur la diagonale (et en général de réorganiser la matrice pour minimiser les problèmes potentiels d’arrondis) est fait même si il n’y a pas déjà de zéro sur la diagonale par n’importe quel solveur direct digne de ce nom. Dans ton cas où c’est juste un code pour jouer un peu, ce n’est évidemment clairement pas une priorité.

@Lucas-84 je pense que je peux diviser la ligne n°1 par 7, la ligne n°2 par 5, la ligne n°3 par 2. Une fois ces divisions faites, j’aurais des 1 en diagonale, comme sur l’identité. Puis je soustrais la ligne n°2 ainsi obtenue par 4/5 x la ligne n°3,

4/5ème de la ligne 1, tu veux dire ?

de sorte à avoir un 0 en (2, 1) (2è ligne, 1ère colonne). Puis je m’occupe de la même façon des autres cases où il faut mettre un 0 à la place du coefficient qu’elles contiennent, en ayant en tête que la modification d’une ligne pour mettre un 0 aura une répercussion sur les autres opérations à effectuer sur cette même ligne pour mettre les 0 sur les cases != 0 restantes (e.g. ayant mis un 0 en (2,1) grâce à une soustraction avec 4/5 x ligne N°3, je devrais tenir compte de cette opération pour mettre un 0 en (2, 3) ).

Tu es conscient qu’après avoir fait une opération comme "soustraire à la ligne 2 la ligne 1 multipliée par 4/5", ton invariant "je n’ai que des 1 sur la diagonale" n’est plus vérifié ? Du coup je vois pas trop l’intérêt de ton étape de normalisation a priori (c’est le même problème dans ton "pseudo-code inspiré de Wikipedia").
Nvm, c’est ce que tu as remarqué dans ton EDIT

Du coup, l’étape où tu mets tous les coefficients à zéro, c’est justement le coeur de l’algo. Tu peux pas trop le passer sous silence.

Dans un premier temps, je m’assure que la diagonale ne contient que des coefficients != 0. Si je rencontre un coefficient == 0 en diagonale, j’échange la ligne qui le contient avec la ligne contenant un coefficient != 0 et qui soit LE coefficient le plus grand de toutes les lignes.

Toutes lignes confondues ? En faisant ça tu peux recréer des zéros plus haut sur la diagonale (imagine quelque chose comme $\begin{bmatrix}0&1&1\\0&0&1\\10&10&1\end{bmatrix}$, si je te suis tu vas faire swap(1, 3) puis swap(1, 2), et te retrouver avec un 0 en (0,0)).

J’ajouterais que c’est assez dérangeant que dans certains cours sur Gauss on mélange sans sourciller les principes mathématiques avec des considérations numériques, je pense que ça embrouille plus qu’autre chose.

@Lucas-84 : je suis d’accord avec tout ton message (pour le "Toutes lignes confondues ?" je m’étais mal exprimé effectivement, c’est "Toutes lignes confondues en partant de la i-ième du haut", et "i" est le numéro de ligne du dernier pivot trouvé).

Du coup j’ai bien compris l’algo de wikipédia. Le souci que j’ai depuis le début, c’est que j’aimerais l’écrire en fonctionnel… uniquement avec des map et reduce, ainsi qu’un max au besoin. Sauf que comme les différentes étapes (mettre des 0, normaliser, soustraire…) sont imbriquées, ça me semble hyper chaud d’utiliser des map et des reduce. Vous avez une idée ?

Ah super ! Merci beaucoup !