Inverse et déterminant d'une matrice de flottants
Bonjour!
Je suis pas trop sûr du forum où poster cette question, j'espère que je ne me suis pas trompé.
Je rencontre un petit problème avec un algorithme que je suis en train de coder.
En gros j'ai une matrice Q de flottants (des arrondis) :
- Si la matrice est inversible (donc son déterminant != 0), je l'inverse et je m'en sers pour résoudre un système.
- Si la matrice n'est pas inversible, j'applique une autre méthode un peu moins efficace.
A priori rien de compliqué, mais je viens de réaliser qu'à cause des arrondis, j'imagine que mon déterminant n'est jamais égal à zéro. Du coup j'ai commencé à écrire un "if (det(Q) < epsilon)" mais est-ce que j'ai le droit en fait?
Dans le cas de flottants, est-ce qu'une matrice avec un déterminant de 0,0001 est inversible ou pas?
J'ai un peu peur que les erreurs que je rencontre viennent de là, et j'imagine que c'est un problème courant mais je n'ai pas réussi à trouver d'infos là dessus (je dois pas taper les bons mots-clés).
Merci d'avance!
Smelk
Inverse et déterminant d'une matrice de flottants
Bonjour, :D
Tu as pressenti le coeur de la difficulté:
Citation:
Envoyé par
Smelk
... je viens de réaliser qu'à cause des arrondis, j'imagine que mon déterminant n'est jamais égal à zéro. Du coup j'ai commencé à écrire un "if (det(Q) < epsilon)" mais est-ce que j'ai le droit en fait? ...
mais il faut cependant sortir d'une l'impasse, celle de la tolérance à accorder à la vérification de l'égalité T = 0 qui, en l'absence de toute autre donnée, n'a aucun sens. Ainsi lorsque tu demandes:
Citation:
Envoyé par
Smelk
... est-ce qu'une matrice avec un déterminant de 0,0001 est inversible ou pas? ...
la réponse peut être affirmative ou non, selon l'ordre de grandeur des éléments non-nuls présents.
1°) Le cas le plus simple est celui de l'égalité entre deux réels positifs: A = B , réductible à la précédente: (A - B) = 0 .
En convenant de noter (EA, EB) les incertitudes affectant chacun des termes, il est conforme à l'intuition que dans le cas de valeurs (A et B) suffisamment proches, les domaines correspondants: [A - EA ; A + EA] , [B - EB ; B + EB]
se recouvrent au moins partiellement (intersection non vide), d'où par exemple (si A > B) la condition: A - EA < B + EB ;
l'écart entre les deux termes admet ainsi une limite supérieure: Abs(A - B) < EA + EB qu'il convient d'adapter numériquement au problème à résoudre.
a) dans le cas de données expérimentales, les incertitudes (EA, EB) découlent de l'estimation de la précision des mesures;
b) s'il s'agit de calculs théoriques, la précision est dans le meilleur des cas limitée par un seuil irréductible lié au codage des nombres flottants: c'est le "epsilon machine" caractéristique du logiciel, que l'on peut tester en cherchant le plus petit nombre vérifiant: 1 + (Emach) > 1 ;
Pour le format Extended, il vaut 2-63 = 1E-19 dans le cas du Pascal et apparentés; il est plus élevé dans la plupart des autres langages.
2°) Soit maintenant le déterminant d'ordre deux: D = (<W , X> , <Y , Z>) = WZ - XY , dont on envisage l'annulation: D = 0 équivalant à WZ = XY .
Si l'on admet pour les 4 termes la même incertitude relative e = EW/Abs(W) = EX/Abs(X) = EY/Abs(Y) = EZ/Abs(Z) , il vient:
EWZ = Abs(W)*EZ + Abs(Z)*EW = Abs(W)*Abs(Z)*(e + e) = 2*e*Abs(WZ)
EXY = Abs(X)*EY + Abs(Y)*EX = Abs(X)*Abs(Y)*(e + e) = 2*e*Abs(XY) ,
et l'équation (D = 0) implique, par analogie avec ce qui a été donné au (1°), et en supposant toujours (WZ > XY):
WZ - EWZ < XY + EXY soit encore: WZ - XY < EWZ + EXY = 2*e*(Abs(WZ) + Abs(XY) ;
il vient donc plus généralement: Abs(D) < 2*e*(Abs(WZ) + Abs(XY)) < 4*e*Max(Abs(WZ), Abs(XY)) .
L'incertitude s'évalue en fonction du produit des termes diagonaux.
3°) Dans le cas d'un déterminant d'ordre (n) plus élevé, on ne peut raisonnablement que conjecturer une estimation de l'incertitude, compte tenu de la complexité des résultats - à moins que tes matrices soient d'un type très particulier.
Il faut supposer tous les éléments indépendants et affectés de la même incertitude relative (e).
Ainsi un déterminant d'ordre (3) comporte la somme algébrique de 3! = 6 termes, dont chacun est le produit de 3 éléments de la matrice, et se trouve donc affecté d'une incertitude relative égale à (3*e); d'où une incertitude absolue ED ~ (3*6*e)*M
dans laquelle intervient la "valeur moyenne" d'un produit M = M1 * M2 * M3
estimable à partir des moyennes (Mj) des valeurs absolues des éléments d'une même colonne.
Cette dernière évaluation me très paraît incertaine, et très excessive: cela ne relève plus d'un simple calcul d'incertitudes, lorsque se superposent un grand nombre de termes indépendants;
Une autre solution consisterait à entreprendre une série de calculs de déterminants (D), pour une matrice A dont les éléments sont faiblement et aléatoirement décalés par rapport à ceux d'une matrice de référence (A°), de déterminant connu (D°); cela devrait conduire à un échantillon de valeur centrale (Dm ~ D°), d'écart-type (SD) et donc à l'intervalle de confiance: [Dm - (2*SD) ; Dm) + (2*SD)].
De tels essais permettraient de remplacer la condition (D = 0) en pratique inutilisable par: Abs(D) < ED = 2 * SD .
Inverse et déterminant d'une matrice de flottants (suite)
Exemple: soit la matrice (A°) dont le déterminant dépend d'un paramètre (x) et s'annule pour x = 5 (D = 60 - 12*x):
Code:
1
2
3
4
|
<< 1 , 4 , 7 >
A° = < 2 , x , 8 >
< 3 , 6 , 9 >> |
puis une seconde matrice (B) dont tous les éléments résultent d'un processus pseudo-aléatoire, mais restent proches de ceux de la précédente, et ne s'en écartent pas de plus de 1.7 % (en plus ou en moins).
Il interviendra par conséquent une double boucle d'instructions du type:
Code:
1
2
3
4
5
6
|
FOR i:= 1 TO 3 DO
FOR j:= 1 TO 3 DO BEGIN
p:= Random; q:= (2 * p) - 1;
p:= 1 + (0.017 * q); B[i, j]:= p * A0[i, j]
END; |
On envisage, pour quelques valeurs de (x), une liste de 11 déterminants aléatoires, de laquelle on déduira la moyenne (Dm), l'écart-type (SD) et l'incertitude définie par la relation: (ED = 2 * SD).
Code:
1
2
3
4
5
6
7
8
9
10
11
12
|
x D° Dm Ed
4.0 12.0 11.95 1.31
4.5 6.0 6.49 1.04
4.8 2.4 2.30 1.49
4.9 1.2 1.04 1.65
5.0 0.0 -0.01 1.55
5.1 -1.2 -1.57 1.77
5.2 -2.4 -2.38 1.91
5.5 -6.0 -5.98 1.90
6.0 -12.0 -12.46 1.88 |
On fait ainsi apparaître les cas pour lesquels le déterminant peut-être considéré comme nul (x = 4.9 , 5.0 et 5.1) ainsi que les cas discutables (x= 4.8 ou 5.2).
Inverse et déterminant d'une matrice de flottants
Une variante permet de s'affranchir de l'hypothèse d'une incertitude relative constante, en consignant:
a) les données expérimentales dans une première matrice (A),
b) les incertitudes absolues correspondantes dans une seconde (B);
puis en calculant les déterminants d'une suite de matrices aléatoires (C) construites terme à terme à partir des précédentes:
Code:
1
2
3
4
5
6
|
FOR i:= 1 TO 3 DO
FOR j:= 1 TO 3 DO BEGIN
p:= Random; q:= (2 * p) - 1;
p:= q * B[i, j]; C[i, j]:= A[i, j] + p
END; |
De l'ensemble des valeurs obtenues, on déduit les résultats déjà mentionnés (Dm, SD et ED).
# Exemple 1:
Code:
1
2
3
4
5
6
7
8
|
A = ( < 800 , 700 , 950 > ,
< 999 , 550 , 875 > ,
< 950 , 822 , 948 > ) ; D<sub>A</sub> = 44 402 700 ;
B = ( < 8 , 7 , 10 > ,
< 10 , 6 , 9 > ,
< 10 , 8 , 9 > ) ; |
On a obtenu pour N = 20 déterminants les valeurs suivantes:
Dm = 45 198 013 ; SD = 3 146 340 et ED = 6 292 680 .
Dans ce cas le déterminant (DA) dépasse très largement l'incertitude absolue (ED), la matrice correspondante (A) est inversible.
# Exemple 2:
Code:
1
2
3
4
5
6
7
8
|
A = ( < 800 , 700 , 950 > ,
< 835 , 550 , 875 > ,
< 671 , 822 , 948 > ) ; D<sub>A</sub> = 55 500 ;
B = ( < 8 , 7 , 10 > ,
< 8 , 6 , 9 > ,
< 7 , 8 , 9 > ) ; |
On trouve cette fois:
Dm = 474 328 ; SD = 2 583 424 et ED = 5 166 848 .
Le déterminant, très inférieur en valeur absolue à l'incertitude, n'a plus de signification et doit être considéré comme nul.
■ Le test numérique proposé ici a le mérite de montrer indirectement ce qu'il peut advenir des valeurs calculées du déterminant; il me paraît cependant assez lourd, donc difficilement applicable au-delà de l'ordre 4 (tu n'as rien dit sur ce point, concernant tes systèmes d'équations).
On peut cependant envisager une méthode générale, en partant de l'indépendance des (n2) données contenues dans la première matrice (A), et du fait que la variation de l'une d'entre elles: bij = a'ij - aij
entraîne une variation du déterminant qui lui est proportionnelle: D'A - DA = bij * Cofij(A) d'après la formule de Laplace.
Cela conduit tout naturellement au calcul de la comatrice correspondante, matrice des cofacteurs de (A): Com(A) = [Cofij(A)] ,
puis à la matrice (C) résultant du produit terme à terme (produit matriciel d'Hadamard) de Com(A) avec la matrice des incertitudes (B):
C = (Com(A) ║ B) .
Il ne reste plus, à ce stade, qu'à calculer la norme de Frobenius de la matrice obtenue:
NF(C) = (Si=1nSj=1n(cij)2)1/2 , susceptible de constituer un excellent jalon.
Reprenons pour illustration numérique les deux exemples précédents:
# Ex 1: D1 = Det(A1) = 44 402 700 ; N1 = NF(C1) = 4 943 371 :
déterminant incontestablement non-nul, puisque très supérieur en valeur absolue au second terme.
# Ex 2: D2 = Det(A2) = 55 500 ; N2 = NF(C2) = 4 640 530 ;
déterminant quasi-nul, puisque inférieur à la limite associée.
Remarquer que les normes obtenues (~ 5E6) présentent le même ordre de grandeur que les incertitudes dont il était initialement question (5 à 6E6).
On dispose désormais, pour toute matrice, d'un test de quasi-nullité du déterminant: Abs(EA) < NF(Com(A) ║ B)
dispensant de toute étude statistique sur une suite matricielle pseudo-aléatoire.
Ma réponse était longue et hésitante, :D parce que ce domaine ne m'est pas familier, et que l'un de mes programmes se plantait sur un problème stupide d'arithmétique modulaire, exigeant de savoir compter jusqu'à 5 ... :ptdr:
# J'ai cherché du côté des matrices aléatoires (random matrices): beaucoup de liens, pointant sur des articles à priori intéressants, mais de niveau inabordable (au moins pour moi) - rien d'utile malheureusement, et le plus souvent incompréhensible :calim2: .
Inverse et déterminant d'une matrice de flottants
J'ai repris pour vérification la première matrice donnée en exemple au message (#3):
Code:
1
2
3
4
|
<< 1 , 4 , 7 > << 9*x - 48 , 6 , 12 - 3*x >
A° = < 2 , x , 8 > ; Com(A°) = < 6 , -12 , 6 > ; B = 0.017 * A° ;
< 3 , 6 , 9 >> < 32 - 7*x , 6 , x - 8 >> |
On obtient alors: NF(Com(A°) ║ B) = 0.017 * NF(Com(A°) ║ A°) = .017 * (1188*x2 - 9720*x + 28080)(1/2)
et pour les valeurs précédemment envisagées du paramètre (x):
Code:
1
2
3
4
5
6
7
8
9
10
11
12
|
x D° NF
4.0 12.0 1.54
4.5 6.0 1.56
4.8 2.4 1.59
4.9 1.2 1.61
5.0 0.0 1.63
5.1 -1.2 1.65
5.2 -2.4 1.67
5.5 -6.0 1.75
6.0 -12.0 1.90 |
Mêmes observations, mêmes conclusions.
La variation observée pour la norme (NF) est par ailleurs beaucoup plus régulière que celle de la dispersion des échantillons pseudo-aléatoires successifs (SD = ED / 2).
Inverse et déterminant d'une matrice de flottants
Bonsoir, :D
Ce problème intéressant mais inédit m'a paru apparenté à la loi de probabilité de Gauss et à un mouvement brownien, sans que je puisse donner, pour sa solution, une formulation claire et facilement compréhensible.
J'ai donc indiqué une série d'opérations intuitives, en m'assurant du terme exact désignant chacune d'entre elles, afin que chacun se procure éventuellement la documentation nécessaire.
Mais à vouloir trop bien faire, peut-être que la présentation apparaissait résolument glaçante pour un lecteur débutant :aie:, je reviendrai donc sur la réponse apportée, afin que Smelk ne reste pas sur son impression de découragement.
Il y a eu des commentaires inattendus.
Citation:
Envoyé par
souviron34
... Les flottants en question étant des arrondis (on pourrait le faire sans cela, mais cela nécessiterait un peu plus de calculs), on connait la "taille" de l'arrondi : nombre de chiffres après la virgule = d .. (qui doit être 2 ou 3 je suppose)
Il suffit de multiplier tous les flottants par 10d , et résoudre le calcul du déterminant.. avec des nombres entiers[/SUP] ...
Tous les résultats (DA, EA, NF) seront effectivement multipliés par le même facteur: 10n*D dans le cas de matrices d'ordre (n); je me suis assuré de leur homogénéité. Et la question de l'incertitude ne disparaîtra pas pour autant: la donnée u = 0.523 ±0.007 équivaut à u' = 523 ±7 .
Mais n'est-ce pas une fuite en avant, avec d'éventuelles complications ? Des données comportant 4 chiffres significatifs conduiront à des résultats à 12 chiffres pour des matrices (3x3); que se passera-t-il à des ordres plus élevés ? Devra-t-on passer en précision absolue ?
D'autant plus que si chaque terme n'est connu qu'à (10-4) près, le résultat sera affecté d'une incertitude relative (n) fois plus grande (estimation très grossière - voir les calculs précédents) ... Quel est l'intérêt de travailler sur de grands nombres entiers ?
Citation:
Envoyé par
tbc92
... Mais il peut y avoir des effets de bord plus vicieux. Si je me souviens bien, le signe du déterminant, ça permet de savoir si le repère formé par nos 3 vecteurs est direct (pouce/index/majeur de la main droite), ou bien indirect (pouce/index/majeur de la main gauche). Et pour un simple problème d'arrondi, un repère peut paraître direct, alors qu'il est indirect.
Que le déterminant de 3 vecteurs orthogonaux puisse changer de signe, cela suppose une incertitude absolue énorme sur les données !
Voilà un exemple simple:
Code:
1
2
3
4
|
<< 1 , 0 , 0 > << e , e , e > << 1 , 0 , 0 >
A = < 0 , 1 , 0 > ; Da = 1 ; B = < e , e , e > ; Com(A) = < 0 , 1 , 0 >
< 0 , 0 , 1 >> < e , e , e >> < 0 , 0 , 1 >> |
Il vient dans ces conditions: C = Com(A)║B = e * A et NF(C) = (3 * e2)1/2 = 31/2 * e ;
il faudrait donc pour une éventuelle inversion de signe: e>~ DA / 31/2 ~ 0.6 , soit 60 % .
Inverse et déterminant d'une matrice de flottants
@ Smelk et sans doute quelques autres, découragés comme lui par l'aspect rébarbatif de la solution proposée (#5).
# Une incertitude absolue sur une grandeur (X) est la limite supérieure des valeurs absolues des écarts observables sur la grandeur considérée:
EX = Max(Abs(X' - X)) .
Elle correspond, dans le cas des données expérimentales, à la précision des mesures effectuées; dans le cas par exemple de X = 1.545 , on aura en l'absence de toute autre indication: EX = 0.0005 (selon la norme AFNOR).
Elle admet pour limite inférieure la précision du stockage en mémoire des nombres flottants: si (Nm) bits sont réservés à leur mantisse, elle doit valoir
EX = X*2-Nm (si je me trompe :D, souviron34 apportera les corrections nécessaires).
Voir Virgule flottante, IEEE 754
# Le déterminant (D) d'une matrice, calculé à partir des (n2) éléments de l'objet, est lui aussi entaché d'une incertitude (ED) qu'il importe d'évaluer, puisqu'on peut le considérer comme quasiment nul à la condition: Abs(D) <~ ED . Là réside l'originalité et l'intérêt de la question posée.
Il ne s'agit pas d'un simple calcul de petites variations: la réponse eût été donnée en deux lignes; il faut évaluer l'erreur moyenne (ou maximale, peu importe) résultant de la superposition de nombreuses petites variations, dont on ne connaît que la limite supérieure (en valeur absolue).
1°) Notion de comatrice (matrice des cofacteurs)
a) C'est la notion au départ la plus difficile, en raison du caractère abstrait de sa notation; et pourtant le calcul de ses éléments s'apparente étroitement à ce que tu connais déjà. Cela vaudrait le coup, pour toi, de t'investir sur ce point.
Par exemple le calcul d'un déterminant (3x3) fait intervenir trois cofacteurs (sans qu'ils soient expressément nommés); ainsi pour la matrice (A) provisoirement notée sans indices
Code:
1
2
3
4
|
<< r , u , x >
A = < s , v , y >
< t , w , z >> |
le développement du déterminant par rapport aux éléments de la première colonne conduit à l'expression:
D = r*(v*z - w*y) + s*(w*x - u*z) + t*(u*y - v*x)
à rapprocher de l'expression théorique:
D = a11*Cof11(A) + a21*Cof21(A) + a31*Cof31(A) .
Le développement par rapport aux éléments de la 2de ou 3me colonne conduirait semblablement aux six autres cofacteurs.
b) Soit maintenant une autre matrice A', ne différant de la précédente que par l'un de ses éléments (celui en ligne 2, colonne 1, par exemple):
Code:
1
2
3
4
|
<< r , u , x >
A' = < s' , v , y >
< t , w , z >> |
l"écart entre la nouvelle valeur du déterminant D' = a11*Cof11(A) + a'21*[Cof21(A) + a31*Cof31(A)
et l'ancienne découle immédiatement de la définition donnée plus haut:
D' - D = (a'21 - a21)*[Cof21(A) .
2°) Le déterminant comme fonction. Chaque cofacteur apparaît comme la dérivée partielle du déterminant par rapport à l'élément correspondant: Cofij(A) = (dD/daij) (désolé pour la notation !).
On pourrait donc envisager l'expression approchée de la variation du déterminant en fonction des variations élémentaires de chacun des éléments de la matrice (hij = a'ij - aij) par la combinaison linéaire des (n2) produits:
D' - D = Si=1nSj=1n(Comij(A) * hij) ,
qui correspond à la formule citée sur ce lien: Det(A + H) = Det(A) + Tr(tCom(A).H) + o(║H║) . (1)
Cependant le problème est d'établir une limite supérieure pour Abs(D' - D), connaissant celles de Abs(hij); d'où la nécessité de consigner:
a) dans une matrice (B) les incertitudes absolues affectant les élément de (A): bij = E(aij) = Max(Abs(hij)) ;
b) dans une autre matrice (C) les produits correspondants cij = Cofij(A) * bij .
3°) On pourrait à ce stade en déduire l'incertitude absolue sur le déterminant: ED = Max(Abs(E' - E)) = Si=1nSj=1n(Abs(cij)) ;
un tel résultat résulte néanmoins d'une majoration outrancière, parce que le fait que toutes les variations (hij) soient simultanément extrêmales et de même signe a très peu de chances de se produire.
C'est pourquoi l'expression finalement retenue: ED = (Si=1nSj=1n(Abs(cij)2))1/2 = (Si=1nSj=1n(cij2))1/2 = NF(C) (1)
constitue une bien meilleure estimation: à précision fixe, (ED) régresse proportionnellement de (n2) à (n).
C'est à cette étape que doit intervenir la loi de probabilité de Gauss. Mais cela nous mènerait vraiment trop loin :D .
(1): Tr(M) désigne la trace d'une matrice, tM sa transposée et (X.Y) le produit matriciel standard.
Tr(tX.Y) correspond à la somme de tous les produits (xij*yij) - à vérifier stylo en main ! Il s'agit d'une extension du produit scalaire.
Au sujet de la norme de Frobenius, j'aurais dû parler de "norme euclidienne", expression beaucoup plus parlante:
NF(X) = ║X║ = Tr(tX.X)1/2 .
Voir Norme (mathématiques)
Norme matricielle
Normes
Le site de Gilles Dubois devrait être consulté beaucoup plus souvent par tous ceux qui, comme moi :D , ont quelques lacunes en maths.
