Algorithme de Gauss-Jordan de Wikipédia : un truc qui cloche ?

Le problème exposé dans ce sujet a été résolu.

Edit important : l’implémentation fonctionne

Les erreurs corrigées sont :

  1. Oubli de stocker des coefficients provenant de la matrice d’entrée originale (je ne parle donc pas de la matrice d’entrée identité) dans une variable temporaire, ce qui fait qu’à la fin d’une opération (e.g. une division) effectuée sur ses lignes, ces coefficients voyaient leur valeur initiale remplacée par le résultat de cette opération, et n’étaient donc pas utilisable pour modifier la matrice d’entrée identité.

  2. La soustraction sur la matrice d’entrée identité était en partie incorrecte.

Voici le code exempt d’erreurs :

 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
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
    /**
      * Returns the identity matrix of the specified dimension
      * @param size the number of columns (i.e. the number of rows) of the desired identity matrix
      * @return the identity matrix of the specified dimension
      */
    def getIdentityMatrix(size : Int): scala.collection.mutable.Seq[scala.collection.mutable.Seq[Double]] = {
      scala.collection.mutable.Seq.tabulate(size)(r => scala.collection.mutable.Seq.tabulate(size)(c => if(r == c) 1.0 else 0.0))
    }

      /**
        * This algorithm processes column by column.
        * STEP 1. It finds the greatest coefficient for the current column (called 'a') and, if it equals 0, returns NULL (since the matrix
        * can't be inverted) ; otherwise (STEP 2.), it swaps the pivot's line with this new line and the pivot becomes the adequate coefficient
        * of this new line
        * STEP 3. It divides the pivot's line by the pivot
        * STEP 4. It sets each of the current column's coefficient to 0 by subtracting the corresponding lines by the pivot's line
        * @return
        */
      def getGaussJordanInvertedMatrix: (Matrix, Matrix) = {

        // We get first the matrix to be inverted, second the identity one
        val mutable_being_inversed_matrix : collection.mutable.Seq[collection.mutable.Seq[Double]] = scala.collection.mutable.Seq(content.map(ms => scala.collection.mutable.Seq(ms:_*)):_*)
        val identity_matrix : collection.mutable.Seq[collection.mutable.Seq[Double]] = getIdentityMatrix(content.length)  // We get the identity matrix. It will be modified
                                                                          // as the original matrix will.

        var id_last_pivot : Int = 0  // ID of the last pivot, i.e. ID of the current column
        content.indices.foreach(general_id_column => {
          println("Current column : " + general_id_column)

          //  STEP 1.
          val id_line_with_max_coefficient_in_this_column = (id_last_pivot until content.length).maxBy(id_line_in_this_column => Math.abs(mutable_being_inversed_matrix(id_line_in_this_column)(general_id_column)))

          if(mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column)(general_id_column) == 0) {
            println("The Gauss-Jordan elimination's algorithm returns an error : indeed, the matrix can't be inverted")

          } else {

            //  STEP 2.
            val tmp_line : scala.collection.mutable.Seq[Double] = mutable_being_inversed_matrix(id_last_pivot)
            mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column)
            mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column) = tmp_line

            val identity_tmp_line : scala.collection.mutable.Seq[Double] = identity_matrix(id_last_pivot)
            identity_matrix(id_last_pivot) = identity_matrix(id_line_with_max_coefficient_in_this_column)
            identity_matrix(id_line_with_max_coefficient_in_this_column) = identity_tmp_line
            println("\nSWAP DONE")
            println(Console.BLUE + "Original matrix :\n" + Console.RESET + mutable_being_inversed_matrix.mkString("\n"))
            println(Console.RED + "Identity matrix :\n" + Console.RESET + identity_matrix.mkString("\n"))

            //  STEP 3.
            val tmp = mutable_being_inversed_matrix(id_last_pivot)(general_id_column)
            mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_last_pivot).map(coefficient => coefficient / tmp)
            identity_matrix(id_last_pivot) = identity_matrix(id_last_pivot).map(coefficient => coefficient / tmp)

            println("\nDIVISION DONE")
            println(Console.BLUE + "Original matrix :\n" + Console.RESET + mutable_being_inversed_matrix.mkString("\n"))
            println(Console.RED + "Identity matrix :\n" + Console.RESET + identity_matrix.mkString("\n"))

            //  STEP 4.
            content.indices.foreach(id_line => {
              val tmp = mutable_being_inversed_matrix(id_line)(general_id_column)

              if(id_line != id_last_pivot) {
                content.indices.foreach(id_column => {
                  mutable_being_inversed_matrix(id_line)(id_column) -= mutable_being_inversed_matrix(id_last_pivot)(id_column) * tmp
                  identity_matrix(id_line)(id_column) -= identity_matrix(id_last_pivot)(id_column) * tmp
                })
              }

            })

            println("\nSUBTRACTION & MULTIPLICATION DONE")
            println(Console.BLUE + "Original matrix :\n" + Console.RESET + mutable_being_inversed_matrix.mkString("\n"))
            println(Console.RED + "Identity matrix :\n" + Console.RESET + identity_matrix.mkString("\n"))
            println()

            id_last_pivot += 1

          }

        })

        (new Matrix(identity_matrix), new Matrix(mutable_being_inversed_matrix))
      }

Hello !

J’ai essayé d’implémenter l’algorithme qui permet d’inverser une matrice, cf. Wikipédia : https://fr.wikipedia.org/wiki/Élimination_de_Gauss-Jordan .

Il consiste à passer sur chaque colonne de la matrice. On trouve la ligne contenant le pivot (on s’arrange pour que ce pivot soit le coefficient maximal de la colonne en cours à partir de la ligne n°<ID de la ligne contenant le dernier pivot trouvé +1 i.e. ID de la colonne en cours +1, c’est la même chose> via un échange de lignes), on divise cette même ligne par le pivot (comme ça le pivot de vient égal à 1), puis on soustrait les lignes notées A à la ligne du pivot multipliée par le coefficient de la colonne en cours à la ligne a€A (de sorte à avoir in fine une colonne à 0).

MAIS. Que se passe-t-il si les conditions suivantes sont réunies ?

  1. On s’apprête à travailler sur la dernière colonne, c’est-à-dire que le dernier pivot trouvé se situe sur la dernière ligne et dernière colonne (i.e. : dernière case de la matrice)

  2. Le travail ou les travaux de la ou des ligne(s) précédente(s) ont impliqué que ce pivot est à 0.

  3. Bien évidemment, grâce au max que j’ai expliqué dans la première phrase de mon post, aucun échange de ligne ne pourra avoir lieu : l’algorithme se trouve alors incapable de trouver un pivot != 0, et on échoue alors l’inversion de la matrice.

Le point n°2 peut avoir lieu dans la mesure où la soustraction de l’avant-dernière ligne, quand on se situe sur l’avant-dernière colonne, donne un 0 dans la dernière case de la matrice. En fait, ça peut aussi être une soustraction qui précède celle-ci qui toutes deux donnent encore 0.

Exemple

En appliquant cet algorithme sur :

List(1.2, 4.3, 8.1, 0.0)

List(0.0, 0.0, 9.1, 0.0)

List(0.0, -5.4, 0.3, 0.0)

List(0.0, 0.3, 0.8, 1.0)

, on obtient :

ArrayBuffer(1.0, 0.0, 0.0, 0.0)

ArrayBuffer(0.0, 1.0, 0.0, 3.3333333333333335)

ArrayBuffer(0.0, 0.0, 1.0, 0.0)

ArrayBuffer(0.0, 0.0, 0.0, 0.0)

Mon implémentation

 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
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
  /**
    * This algorithm processes column by column.
    * STEP 1. It finds the greatest coefficient for the current column (called 'a') and, if it equals 0, returns NULL (since the matrix
    * can't be inverted) ; otherwise (STEP 2.), it swaps the pivot's line with this new line and the pivot becomes the adequate coefficient
    * of this new line
    * STEP 3. It divides the pivot's line by the pivot
    * STEP 4. It sets each of the current column's coefficient to 0 by subtracting the corresponding lines by the pivot's line
    * @return
    */
  def getGaussJordanInvertedMatrix: (Matrix, Matrix) = {

    // We get first the matrix to be inverted, second the identity one
    val mutable_being_inversed_matrix : collection.mutable.Seq[collection.mutable.Seq[Double]] = scala.collection.mutable.Seq(content.map(ms => scala.collection.mutable.Seq(ms:_*)):_*)
    val identity_matrix : collection.mutable.Seq[Seq[Double]] = getIdentityMatrix(content.length).to[collection.mutable.Seq]  // We get the identity matrix. It will be modified
                                                                      // as the original matrix will.

    var id_last_pivot : Int = 0  // ID of the last pivot, i.e. ID of the current column
    content.indices.foreach(general_id_column => {

      //  STEP 1.
      val id_line_with_max_coefficient_in_this_column = (id_last_pivot until content.length).maxBy(id_line_in_this_column => mutable_being_inversed_matrix(id_line_in_this_column)(general_id_column))

      if(mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column)(general_id_column) == 0) {
        null

      } else {

        //  STEP 2.
        val tmp_line : scala.collection.mutable.Seq[Double] = mutable_being_inversed_matrix(id_last_pivot)
        mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column)
        mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column) = tmp_line

        //  STEP 3.
        mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_last_pivot).map(coefficient => coefficient / mutable_being_inversed_matrix(id_last_pivot)(general_id_column))

        //  STEP 4.
        content.indices.foreach(id_line => {
          if(id_line != id_last_pivot) {
            content.indices.foreach(id_column => {
              mutable_being_inversed_matrix(id_line)(id_column) = mutable_being_inversed_matrix(id_line)(id_column) - mutable_being_inversed_matrix(id_line)(general_id_column) * mutable_being_inversed_matrix(id_last_pivot)(id_column)
            })
          }
        })

        id_last_pivot += 1

      }

    })

    (new Matrix(identity_matrix), new Matrix(mutable_being_inversed_matrix))
  }

Exécution de mon implémentation : https://jsfiddle.net/ueuvf3ra/1/

+0 -0

Là j’avoue que je sèche, car je pense avoir pas trop mal compris l’algo. J’ai écrit cette implémentation, commentée :

 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
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
  /**
    * This algorithm processes column by column.
    * STEP 1. It finds the greatest coefficient for the current column (called 'a') and, if it equals 0, returns NULL (since the matrix
    * can't be inverted) ; otherwise (STEP 2.), it swaps the pivot's line with this new line and the pivot becomes the adequate coefficient
    * of this new line
    * STEP 3. It divides the pivot's line by the pivot
    * STEP 4. It sets each of the current column's coefficient to 0 by subtracting the corresponding lines by the pivot's line
    * @return
    */
  def getGaussJordanInvertedMatrix: (Matrix, Matrix) = {

    // We get first the matrix to be inverted, second the identity one
    val mutable_being_inversed_matrix : collection.mutable.Seq[collection.mutable.Seq[Double]] = scala.collection.mutable.Seq(content.map(ms => scala.collection.mutable.Seq(ms:_*)):_*)
    val identity_matrix : collection.mutable.Seq[Seq[Double]] = getIdentityMatrix(content.length).to[collection.mutable.Seq]  // We get the identity matrix. It will be modified
                                                                      // as the original matrix will.

    var id_last_pivot : Int = 0  // ID of the last pivot, i.e. ID of the current column
    content.indices.foreach(general_id_column => {

      //  STEP 1.
      val id_line_with_max_coefficient_in_this_column = (id_last_pivot until content.length).maxBy(id_line_in_this_column => mutable_being_inversed_matrix(id_line_in_this_column)(general_id_column))

      if(mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column)(general_id_column) == 0) {
        null

      } else {

        //  STEP 2.
        val tmp_line : scala.collection.mutable.Seq[Double] = mutable_being_inversed_matrix(id_last_pivot)
        mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column)
        mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column) = tmp_line

        //  STEP 3.
        mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_last_pivot).map(coefficient => coefficient / mutable_being_inversed_matrix(id_last_pivot)(general_id_column))

        //  STEP 4.
        content.indices.foreach(id_line => {
          if(id_line != id_last_pivot) {
            content.indices.foreach(id_column => {
              mutable_being_inversed_matrix(id_line)(id_column) = mutable_being_inversed_matrix(id_line)(id_column) - mutable_being_inversed_matrix(id_line)(general_id_column) * mutable_being_inversed_matrix(id_last_pivot)(id_column)
            })
          }
        })

        id_last_pivot += 1

      }

    })

    (new Matrix(identity_matrix), new Matrix(mutable_being_inversed_matrix))
  }
+0 -0

Logiquement, si ça arrive, c’est que tu as une ligne qui n’est pas linéairement indépendante des autres. Donc la matrice n’est pas inversible. Ça vient forcément de l’algo.

Après pour ton code, il est assez long et je ne maîtrise pas du tout scala.

Sans gestion des zeros, ça donne ça :

 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
a = [[1, 2, 3],
    [4, 7, 3],
    [1, 5, 7]]

def subRow( a, b, c ):
   for i in range(len(a)):
       a[i] = a[i] - c*b[i]

   return a

def mult(a, c):
   return list(map(lambda x: c*x, a))

def normalize( i , a ):
   c = a[i]

   for j in range(len(a)):
       a[j] = a[j]/c
   return a

for i in range(len(a)):
   a[i] = normalize(i, a[i])

   for j in range(len(a)):
       if j != i :
           a[j] = subRow(a[j], a[i], a[j][i])

La gestion des zeros étant plus chiante que difficile

 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
28
29
30
a = [[1.2, 4.3, 8.1,0],
    [0, 0, 9.1, 0],
    [0, -5.4, 0.3, 0],
    [0, 0.3, 0.8, 0.1]]

def subRow( a, b, c ):
   for i in range(len(a)):
       a[i] = a[i] - c*b[i]

   return a

def mult(a, c):
   return list(map(lambda x: c*x, a))

for i in range(len(a)):
   if a[i][i] == 0:
       j=i+1
       while j <= len(a) and a[j][i] == 0:
           j+=1
       if j > len(a) :
           print("Inversion impossible")
           break
       a[i],a[j]=a[j],a[i]

   a[i] = mult(a[i], 1/a[i][i])

   for j in range(len(a)):
       if j != i :
           a[j] = subRow(a[j], a[i], a[j][i])
print(a)

Version gérant les zéros je crois

+0 -0

@ache : tu ne calcules pas l’inverse (du coup il est difficile de tester ton code…) et ta version gérant les 0 renvoie un IndexError car ta condition en 18 est trop large (et de fait ta condition en 20 est trop restreinte).

@The-Aloha-Protocol : ton code est par super lisible comme les lignes font 3 km de long et je ne connais pas scala non plus… Une façon de déboguer ton code serait de sortir ta matrice à chaque étape pour être sûr qu’il n’y en a pas une qui foire et de vérifier que le pivot choisi est le bon.

+0 -0

Oui effectivement pour la condition. Pour l’inverse, peu importe c’est super simple à rajouter :

 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
28
29
30
31
32
33
34
35
36
37
a = [[1.2, 4.3, 8.1,0],
    [0, 0, 9.1, 0],
    [0, -5.4, 0.3, 0],
    [0, 0.3, 0.8, 0.1]]

def subRow( a, b, c ):
   for i in range(len(a)):
       a[i] = a[i] - c*b[i]

   return a

def mult(a, c):
   return list(map(lambda x: c*x, a))


ID=[]
for j in range(len(a)):
   ID.append(list([0]*len(a)))
   ID[j][j] = 1
   a[j] += ID[j]

for i in range(len(a)):
   if a[i][i] == 0:
       j=i+1
       while j < len(a) and a[j][i] == 0:
           j+=1
       if j > len(a) :
           print("Inversion impossible")
           break
       a[i],a[j]=a[j],a[i]

   a[i] = mult(a[i], 1/a[i][i])

   for j in range(len(a)):
       if j != i :
           a[j] = subRow(a[j], a[i], a[j][i])
print(a)

Si tu pouvais sortir la matrice à chaque étape ou comparer avec la version python, ça serrait cool

+0 -0

Ton problème est lorsque tu fais les soustractions des lignes entre elles. Typiquement, quand tu es sûr la deuxième ligne :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
ArrayBuffer(ArrayBuffer(1.0, 3.5833333333333335, 6.75, 0.0),
            ArrayBuffer(0.0, 1.0, 2.666666666666667, 3.3333333333333335),
            ArrayBuffer(0.0, -5.4, 0.3, 0.0),
            ArrayBuffer(0.0, 0.0, 9.1, 0.0))

SUBTRACTION & MULTIPLICATION DONE
ArrayBuffer(ArrayBuffer(1.0, 0.0, 6.75, 0.0),
            ArrayBuffer(0.0, 1.0, 2.666666666666667, 3.3333333333333335),
            ArrayBuffer(0.0, 0.0, 0.3, 0.0),
            ArrayBuffer(0.0, 0.0, 9.1, 0.0))

Sur toutes les lignes, seulement la deuxième colonne est affectée et pas le reste…

+0 -0 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
DIVISION DONE
ArrayBuffer(1.0, 3.5833333333333335, 6.75, 0.0)
ArrayBuffer(0.0, 0.0, 9.1, 0.0)
ArrayBuffer(0.0, -5.4, 0.3, 0.0)
ArrayBuffer(0.0, 0.3, 0.8, 1.0)

SUBTRACTION & MULTIPLICATION DONE
ArrayBuffer(1.0, 3.5833333333333335, 6.75, 0.0)
ArrayBuffer(0.0, 0.0, -52.324999999999996, 0.0)
ArrayBuffer(0.0, 13.950000000000001, -1.7249999999999999, 0.0)
ArrayBuffer(0.0, -0.7749999999999999, -4.6000000000000005, 1.0)

Ça me parait foireux, la première colonne ne contient déjà que des 0 sur les 3 lignes donc ces 3 lignes ne devrait pas changer…

Et j’ai l’impression que tu sors pas l’inverse à la fin, ça serait un moyen simple pour toi de contrôler si l’algo marche ou pas.

+0 -0

O_o

Ton code agis bizarement. Pour chaque ligne tu fais la première -1 à chaque terme multiplié terme à terme avec les autres. Ça n’a pas de sens >_<

Tu sembles confondre les opérations.

Ce que tu fais :

  • Ajouter une constante à chaque terme d’une liste (en gros, tu ajoutes la ligne $\alpha (1,1,1,1)$ qui n’existe probablement pas encore
  • Multiplié deux listes liste terme à terme.

Ce que tu as le droit de faire:

  • Multiplié une ligne par une constante différente de 0
  • Ajouter/soustraire une liste à une autre

Ce que tu dois faire c’est normaliser la ligne courrante (disons j). Ça c’est bon. Ensuite il va falloir parcourir les lignes afin de d’obtenir que des 0 (sauf sur la ligne j) dans la colonne j, pour cela, tu soustraits la ligne normalisée (j donc) multipliée par le coefficient à annuler.

+0 -0

Edit: j’ai vraiment l’impression que le problème vient de :

multipliée par A[i,j]

, qui se situe ici : https://fr.wikipedia.org/wiki/Élimination_de_Gauss-Jordan . Vous confirmez que j correspond bien au numéro de la colonne en cours de traitement, et que A[i,j] utilisée sur chaque ligne i ne dépend, parmi toutes les colonnes, QUE de la colonne j pour chaque ligne ?

======

@ache, je ne comprends pas ce que tu essaies de me faire remarquer. Je m’explique.

Le bout de code auquel tu fais référence me semble être la partie soustraction/multiplication. C’est-à-dire :

1
2
3
4
5
6
7
content.indices.foreach(id_line => {
        if(id_line != id_last_pivot) {
          content.indices.foreach(id_column => {
            mutable_being_inversed_matrix(id_line)(id_column) -= mutable_being_inversed_matrix(id_line)(general_id_column) * mutable_being_inversed_matrix(id_last_pivot)(id_column)
          })
        }
      })

Pour chaque ligne tu fais la première -1 à chaque terme multiplié terme à terme avec les autres

Non ; je fais : chaque ligne MOINS la case de cette ligne à la colonne courante FOIS la ligne du pivot. Dans l’algo de Wiki, c’est ce qui correspond à :

1
 |   |   |   |   Soustraire à la ligne i la ligne r multipliée par A[i,j] (de façon à annuler A[i,j])

Ajouter une constante à chaque terme d’une liste (en gros, tu ajoutes la ligne α(1,1,1,1) qui n’existe probablement pas encore

Je n’ajoute jamais de ligne dans mon code. Sauf si par "ajout" tu entends "soustraction" ? Mais ça c’est normal.

Multiplié deux listes liste terme à terme.

Non ; je soustrais deux listes terme à terme (c’est la définition que j’ai donnée à "Soustraire la ligne i à la ligne r […]").

Ce que tu as le droit de faire:

Multiplié une ligne par une constante différente de 0 Ajouter/soustraire une liste à une autre

C’est ce que je fais.

Ensuite il va falloir parcourir les lignes afin de d’obtenir que des 0 (sauf sur la ligne j) dans la colonne j, pour cela, tu soustraits la ligne normalisée (j donc) multipliée par le coefficient à annuler.

C’est bien ce que je fais.

+0 -0
Connectez-vous pour pouvoir poster un message.
Connexion

Pas encore membre ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte