[Matrices] Comment calculer le Déterminant d'une matrice 4x4

HELLO LE MONDE!! :D

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.
:wink: :wink:

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 :roll:

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 ! :wink:
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...

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) -

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 :)

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 ;)

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+

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) :P

merci pour votre aide :)

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).

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 ?

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:

---

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:

---

1 2 3

|a11 a12 a13| | 0 b22 b23| |a31 a32 a33|

---

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

Code:

---

1 2 3

|a11 a12 a13| | 0 b22 b23| | 0 0 c33|

---

Et le determinant vaut

Code:

---

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.

ok, c'est peut etre en effet plus rapide :)

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

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.

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 :P

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

Citation:

---

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)...

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

Code:

---

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:

---

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:

---

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+

très intéressant merci :)

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

voila ce que j'ai fait a present :

Code:

---

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 :)

Citation:

---

Envoyé par Peltchag

voila ce que j'ai fait a present :

Code:

---

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).

Citation:

---

Envoyé par Peltchag

Code:

---

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)

Citation:

---

Envoyé par Peltchag

Code:

---

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).

Citation:

---

Envoyé par Peltchag

Code:

---

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

---

---

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

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 !