Discussion :

[cyber_N]

HELLO LE MONDE!!

Je souhaiterais savoir comment calculer le déterminant d'une matrice 4x4
(ET NON D'UNE MATRICE 3X3 ou 2x2)

Avez des explications et/ou des liens intéressants ??

Merci encore.
Séb.

[Wavyx]

tu peux commencer par là:
http://mathworld.wolfram.com/Determinant.html

mais il existe pe des methodes spécifiques pour les 4x4 ??

sinon ben ya la méthode générale

[alveric]

Tu as besoin de performances optimisees, ou tu veux juste un truc qui marche ? C'est une matrice entiere, reelle, complexe ?
Si tu veux un truc qui ne marche que pour les 4x4, mais qui le fait bien... bah tape la formule complete !
Ou alors tu prends la formule 7 du lien que Wavyz t'a donne, et tu retombes sur 4 determinants de degre 3, pour lesquels la formule est donnee quelques lignes plus bas.
Ou alors tu fais une triangulation de la matrice, et tu fais le produit de la diagonale grace a la formule 24.

Il y a le choix... Mais je ne me souviens pas de methode particuliere pour les matrice d'ordre 4...

[Matthieu Brucher]

La triangularisation à l'aidede matrices de Houseolder est sans doute une des meilleures méthodes que je connaisse pour calculer un déterminant - complexité en O(n^3) -

[Peltchag]

Bonjour, je fais remonter ce topic, car moi aussi j'aurais besoin d'un algo permettant de calculer le déterminant d'une matrice. Mais il me faudrait une formule générale, car l'ordre de mes matrices va varier de 2 à 6 en général.
J'ai bien entendu regarder le lien de Wavyx, mais il n'est pas évident d'en déduire la forumle générale.
J'ai bien compris qu'il fallait faire :
Somme(i=1 à k) des a(i,j)C(i,j) où C(i,j) = (-1)^(1+j)M(i,j)
mais il faut aussi calculer M(i,j)

Bref, ma date de fin de projet approche à grands pas, et j'ai pas énormément de temps (hélas) à passer à pondre de nouveaux algos, je voudrais donc savoir s'il n'y a pas un algo existant quelque part, voir même une implémentation dans un language connu (comme du C par exemple) qui me permettrait de déduire l'algo

Merci beaucoup

[Wavyx]

en gros il suffit de faire un appel récursif jusqu'à une dimension 3 par exemple.
sinon tu peux regarder dans matlab comment ils ont fait ou encore en cherchant un peu sur le web et sur source forge:
http://sourceforge.net/projects/jlinalg/
http://sourceforge.net/projects/slate/
http://sourceforge.net/projects/ezmatrixcalc/

mais je m'y connais pas en librairie pour les matrices ou l'algèbre linéaire donc pê que qqun aura de meilleurs liens/sources

[mathieu_t]

Franchement la méthode donnée par Wavx est de mon point de vue la plus simple à implémenter, d'autant que ça prend la forme d'un algo récursif (on décompose le déterminant d'une matrice n*n par des déterminants de matrice (n-1)*(n-1) etc ...)

Sinon j'ai trouvé en tapant sous Google "class Matrix C++" le lien suivant : http://www.techsoftpl.com/matrix/download.htm mais il y en a bien d'autres...

Cela calcule le déterminant apparemment...

A+

[Peltchag]

en effet, j'avais pas pense a regarder du cote de sourceforge

je vais essayer d'eplucher ca rapidement

mais c'est vrai qu'au final, je sens qu'il sera plus rapide de chercher l'algo puis de l'implementer, je n'ai pas non plus envie de passer par une librairie externe (je sais, je suis ch... difficile)

merci pour votre aide

[alveric]

Le probleme avec la formule utilisant les mineurs d'ordre n-1, c'est qu'il faut appliquer l'algo recursivement sur les mineurs. Or, ces mineurs sont des portions pour la plupart non-contigues de la matrice d'origine, ce qui signifie qu'il faut reconstruire, a chaque niveau d'appel, une matrice n-1 en recopiant les valeurs voulues...
Du coup, on se retrouve a recopier (n-1)^2 valeurs dans la matrice qui sert de mineur, on applique l'algo sur la nouvelle matrice pour connaitre le determinant, on multiplie si necessaire par -1, et on fait ca n fois !
Autrement dit, C(n) = n * C(n-1) + n*((n-1)^2 copies de valeur) + (n-1 additions) + (n/2 multiplications par -1)
A moins qu'il y ait une optimisation possible pour les copies de valeurs, ou un moyen efficace de travailler directement sur la matrice d'origine, ce dont je doute (mais je peux me tromper), ca devient vite lourd.
De plus, ca signifie qu'il faut avoir en memoire une matrice d'ordre n, une d'ordre n-1... jusqu'a 1.
Bon, pour n <= 6, ca peut etre acceptable, apres tout je n'ai pas fait les tests, mais la triangularisation serait plus efficace a priori (enfin, a mon priori).

[Peltchag]

je suis bien conscient du probleme en effet, mais je ne connais pas d'autre méthode

peux tu m'en dire plus sur la triangularisation ?

[alveric]

Cette methode utilise la propriete de multilinearite du determinant (16 du premier lien de Waxyz), et le fait que le determinant d'une matrice triangulaire est le produit des termes de la diagonale (24).
Le but est d'annuler toute la partie sous la diagonale, par des combinaison lineaire de lignes.
ex:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
|a11 a12 a13|
|a21 a22 a23|
|a31 a32 a33|

A la premiere iteration, tu annules a21 en ajoutant -(a21/a11) fois la premiere ligne a la deuxieme. Tu obtiens

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
|a11 a12 a13|
| 0  b22 b23|
|a31 a32 a33|

Puis tu recommences pour annuler a31, puis a32. Tu obtiens alors

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
|a11 a12 a13|
| 0  b22 b23|
| 0  0   c33|

Et le determinant vaut

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

det(A) = a11 * b22 * c33

Evidemment, si a11 est nul, il faut permuter des lignes... (vu que le determinant est aussi invariant par permutation de lignes)
Je ne sais plus la complexite de l'algo, mais ca doit se retrouver.

[Peltchag]

ok, c'est peut etre en effet plus rapide

aurais tu un bout d'algo ou de code de cette methode ?

[alveric]

Je te l'ai deja presque donne...
Tu commences par la premiere ligne. Si a11 est nul, tu cherches sur la premiere colonne un terme non nul. Si tu en trouves un, tu permutes les deux lignes et passes a la suite. Sinon le determinant est nul et tu peux t'arreter.
Une fois que tu as un a11 non nul, tu cherches a annuler a21, puis a31...an1. (un truc: calcule 1/a11 une seule fois dans une variable temporaire. Les divisions sont toujours des calculs longs pour le CPU, alors autant en faire le moins possible)
Et tu recommences avec b22...
La condition d'arret, c'est soit arriver a un akk nul et tous les ajk , j>k nuls aussi, soit avoir triangule totalement la matrice.

[Peltchag]

oui, je sais, je suis deja parti sur l'algo, mais ca aurait bien pratique si quelqu'un l'avait deja fait, car bcp plus rapide

je suis desole, je suis un peu pris par le temps

[mathieu_t]

mais la triangularisation serait plus efficace a priori

Oui en effet, j'ai regardé rapido un algo que j'ai, ça passe par la triangularisation...

Je ne peux pas poster la classe Matrix que j'ai (code trop long), mais je peux l'envoyer par mails à ceux qui le souhaitent...

En revanche dans le lien que j'ai donné au dessus, ils le font par le pivot (méthode récursive donnée au dessus)...

[mathieu_t]

Je peux quand même donner les bouts intéressants :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// determinant(u, w, v) -- Returns the determinant of a square matrix,
// which is the product of the diagonal weights in W, O(n).
// Determinant of a rectangular matrix is not well defined, 0 is returned.
// Elimination is done on the matrix and singularies are marked away
// if W is empty, otherwise the decomposition cached in (U, W, V) is used.
 
double Matrix::determinant(Matrix& u, Matrix& w, Matrix& v) const {
  if (this->num_rows != this->num_cols) {
    cerr << "Return 0 for determinant of a rectangular matrix." << endl;
    return 0;
  }
  if (w.data == NULL)
    this->singularities(u, w, v);               // cache SVD in U,W,V.
  double product = 1;
  for (unsigned long k = 0; k < w.num_cols; k++)		// product of diagonals in W.
    product *= w(k, k);				// small weights should be zeroed
  return product;				// out, so that determinant = 0.
}

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
// singularities -- Returns the number of singular weights which make
// the well-condition number less than given lower bound threshold. O(n).
// Default lower bound on the well-condition number is 1.0e-6.
// These singular weights are marked by setting them to zero.
// The inverse of these zero weights will be set to zero instead of
// infinite, which corresponds to finding least-norm solution.
 
int Matrix::singularities(Matrix& u, Matrix& w, Matrix& v, double lbound) const {
  if (w.data == NULL)
    this->singular_value_decomposition(u, w, v); // cache SVD in U,W,V
  double max_w = 0;
  for (unsigned long k=0; k < w.num_cols; k++) {// find max absolute weight
    double weight = (double) fabs(w(k, k));		// in diagonal
    if (max_w < weight) max_w = weight;
  }
  double lbound_w = (double) (max_w * lbound);	// double must support "* double"
  unsigned long count = 0;
  for (unsigned long l=0; l < w.num_cols; l++) {// find near-singular weights
    double weight = (double) fabs(w(l, l));		// and zero them out
    if (weight < lbound_w) {
      w(l, l) = 0;
      count++;					// count the number of singular
    }						// diagonals.
  }
  return count;
}

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
// svd(u, w, v) -- Computes the singular value decomposition of THIS matrix, and
// returns the 3 matrices U, W, V, such that: *this = U * W * V.transpose(). O(n^3).
// Small diagonals in W should be zeroed out, to clearly mark singularities
// in the matrix. The range space are columns of U such that the
// corresponding diagonals in W are non zero. The null space are columns
// of V such that the corresponding diagonals in W are zero.
// The vectors in U and V are orthonormal. If the matrix is symmetric, then
// U and V are equal and they contain the eigenvectors,
// W contains the eigenvalues.
// The eigenvalues are positive if the matrix is positive definite.
// The matrices U, W, V are returned for the user to mark out singularities
// based on some lower bound on the well-condition number, such as 1.0e-6.
// These matrices are cached so that subsequent computations such
// as rank, determinant, condition number, inverse, back substitution, will
// not need to compute the expensive SVD again.
// See Numerical Recipees for reference.
 
bool Matrix::singular_value_decomposition(Matrix& u, Matrix& w, Matrix& v) const {
    unsigned long n = this->num_cols;			 // dimension of vectors in V.
    unsigned long m = (n>this->num_rows)? n : this->num_rows; // dimension of vectors in U.
    u.resize(m, n);				 // fill *this into u.
    u.update(*this);
    for (unsigned long i = this->num_rows; i < m; i++)	// clear the extra rows
      for (unsigned long j = 0; j < n; j++)
	u.data[i][j] = 0;
    w.resize(n, n);
    w.fill(0);					// clear w
    v.resize(n, n);
    if (Matrix::svdcmp(u.data, m, n,	// NR does all work.
				 w.data[0], // and stores results
				 v.data)) { // in place
      for (unsigned long k = 1; k < n; k++) {
	w(k, k) = w(0, k);			// fill diagonal of w
	w(0, k) = 0;
      }
      return true;
    } else
      return false;
}

Voilà, en espérant que ça aidera...

A+

[Peltchag]

très intéressant merci

je t'ai envoyé mon mail en MP, pour avoir la suite si possible

[Peltchag]

voila ce que j'ai fait a present :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
    //----- parcours de la matrice pour verrifier que son determinant n'est pas nul -----
    Indice = 1
     Répéter
        //on cherche d'abord sur la ligne
        LigneNulle = Vrai
        //on se place sur la ligne voulu
        Ligne = Indice
        //et on parcourt chaque colonne
        Colonne = 1
         Répéter
            //des qu'on tombe sur un coefficient different de zero
            Si (TabCoeffs(Ligne, Colonne) <> 0) Alors
                //on dit que la ligne n'est pas nulle
                LigneNulle = Faux
             FinSi
            //on passe a la colonne suivante
            Colonne = Colonne + 1
        //on sort s'il n'y a plus de colonnes
        //ou si on sait que la ligne n'est pas nulle
         Jusqu 'à ((Colonne > Dimension) Ou (LigneNulle = Faux))
 
        //si la ligne est bonne, on regarde la colonne
        //par le meme principe
        Si (LigneNulle = Faux) Alors
            ColonneNulle = Vrai
            Colonne = Indice
            Ligne = 1
             Répéter
                Si (TabCoeffs(Ligne, Colonne) <> 0) Alors
                    ColonneNulle = Faux
                 FinSi
                Ligne = Ligne + 1
             Jusqu 'à ((Ligne > Dimension) Ou (ColonneNulle = Faux))
         FinSi
 
        //une fois qu'on a regarde ligne et colonne,
        //on passe a l'indice suivant
        Indice = Indice + 1
    //on s'arrete si on a parcouru toute la matrice
    //ou si on a trouve une ligne ou une colonne nulle
     Jusqu 'à ((Indice > Dimension) Ou (LigneNulle = Vrai) Ou (ColonneNulle = Vrai))
 
    //si une ligne ou une colonne est nulle, le determinant est nul
    Si (LigneNulle = Vrai) Ou (ColonneNulle = Vrai) Alors
        DeterminantCalcule = 0
    //sinon, on se lance dans le calcul
    Else
        //methode par la triangularisation (pivot de gauss)
        //le determinant d'une matrice triangulaire est le produit de sa diagonale
        //on parcourt la matrice ligne par ligne
        Pour LignePivot <- 1 à Dimension
            //si le coefficient de la diagonale est nul
            Si (TabCoeffs(LignePivot, LignePivot) = 0) Alors
                //on cherche a permutter la ligne avec une autre, dont le coeff n'est pas nul
                LignePermutation = LignePivot + 1
                //on parcourt les lignes en dessous de la ligne actuelle
                Tant que ((LignePivot + 1 <= Dimension) et (PermutationEffectuee = Faux))
                    //des qu'un coeff n'est pas nul
                    Si (TabCoeffs(LignePermutation, LignePivot) <> 0) Alors
                        //on inverse les deux lignes
                         PermutterLigne(LignePivot, LignePermutation)
                        PermutationEffectuee = Vrai
                     FinSi
                //on s'arrete si on a parcouru toutes les lignes
                //ou si on a effectue une permutation
                 FinTq
             FinSi
 
            //on met les coefficients dans la colonne du pivot a zero
            Pour Ligne <- LignePivot + 1 à Dimension - 1
                 //Lj <- Lj - (aji/aii) Li
                 OperationSurLigne(Ligne, LignePivot, -TabCoeffs(LignePivot, Ligne) / TabCoeffs(Ligne, Ligne))
             FinPour
         FinPour
 
        //on calcule le determinant en faisant le produit de la diagonale
        DeterminantCalcule = 1
        Pour Ligne <- 1 à Dimension
            DeterminantCalcule = DeterminantCalcule * TabCoeffs(Ligne, Ligne)
        FinPour
     FinSi
 
    Retourner(DeterminantCalcule)

il n'y a pas les declarations de variables ni les fonctions annexes, mais c'est deja assez long comme ca

j'ai pas encore pu vérifier si ca marchait bien, je vais m'y pencher a present

[alveric]

voila ce que j'ai fait a present :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
            //des qu'on tombe sur un coefficient different de zero
            Si (TabCoeffs(Ligne, Colonne) <> 0) Alors

Avec des flottants, il ne faut jamais comparer directement avec zero, mais tester si la valeur absolue est inferieure a une constante (le langage peut en fournir une adaptee a la precision du systeme).

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
       Pour LignePivot <- 1 à Dimension
            Si (TabCoeffs(LignePivot, LignePivot) = 0) Alors
                LignePermutation = LignePivot + 1
                Tant que ((LignePivot + 1 <= Dimension) et (PermutationEffectuee = Faux))
                    Si (TabCoeffs(LignePermutation, LignePivot) <> 0) Alors
                         PermutterLigne(LignePivot, LignePermutation)
                        PermutationEffectuee = Vrai
                     FinSi
                 FinTq
             FinSi

Meme remarque pour le test de nullite d'un flottant.
Apres cette boucle, il faut verifier si PermutationEffectuee est vraie ou non. Si elle est fausse, alors c'est un cas d'arret et le determinant est nul. (cf. ce que j'ai dit plus haut)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
            //on met les coefficients dans la colonne du pivot a zero
            Pour Ligne <- LignePivot + 1 à Dimension - 1
                 //Lj <- Lj - (aji/aii) Li

Le 1/aii devrait etre calcule avant (ok, pour une matrice 4x4, le gain est faible. Mais ici l'algo est ecrit pour etre independant de l'ordre de la matrice, donc ca compte).

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

                OperationSurLigne(Ligne, LignePivot, -TabCoeffs(LignePivot, Ligne) / TabCoeffs(Ligne, Ligne))

"OperationSurLigne" ? Pas tres explicite, comme nom...

[Peltchag]

merci pour les remarques

mais je ne vais pas travailler avec des flottants moi, juste des decimaux (2 ou 5 chiffres apres la virgule, jamais plus). de plus, mon systeme d'equations concerne une production/consommation, et il y aura forcement des valeurs a zero dans la matrice, c'est pourquoi je suis parti sur de tels tests en fait.

pour la procedure, je sais que le nom est un peu sordide, mais ne me souvenant plus du nom mathematique de l'operation, j'ai mis un nom vite fait (flemme oblige), et ai ajoute un commentaire juste au dessus pour dire ce que cetait

je vais tenir compte de tout ce que tu m'as dit, et y apporte une petite retouche, merci bien !