Discussion :

[Gm7468]

Bonjour,

Je suis sur un projet assez conséquent pour réaliser un programme qui fait diverses actions et calculs, dont des manipulations de matrices.

J'utilise donc dans mon programme des matrices, avec des calculs de base (transposée, produit) et des inversions. J'avais choisis BOOST pour faire tout ce dont j'ai besoin à ce niveau, mais je suis maintenant bloqué sur un problème d'inversion. Toute la documentation que j'ai pu trouver aussi bien sur les forums et sur BOOST concernant les inversions de matrices m'a envoyé ici, j'utilise donc InvertMatrix pour mes inversions.

je ne comprends pas la ligne:

// REMEMBER to update "lu.hpp" header includes from boost-CVS

Pour utiliser cette fonction InvertMatrix, il faudrait que je bidouille le lu de BOOST?

Cette fonction me renvoie de bonnes choses à certains endroits, et fait complètement planter l'application à d'autres... Mais je ne comprends pas pourquoi.

Par exemple, pour la matrice:

0.127213 -0.00600185 -0.000149287 0.587534
-0.00600185 0.00553831 -0.00593579 -0.0335789
-0.000149287 -0.00593579 0.0298501 0
0.587534 -0.0335789 0 3

La fonction me renvoie:

82.8377 -11.3315 -1.83901 -16.3502
-11.3315 252.673 50.1882 5.04737
-1.83901 50.1882 43.4716 0.921916
-16.3502 5.04737 0.921916 3.59192

Par contre, Excel me donne:

82.8339694 -11.33084355 -1.83890416 -16.34941688
-11.33084355 252.6730142 50.18818708 5.047245904
-1.83890416 50.18818708 43.47161365 0.921894277
-16.34941688 5.047245904 0.921894277 3.591773088

Ce qui fait une différence de:

0.003730597 -0.000656448 -0.00010584 -0.000783119
-0.000656448 -1.41645E-05 1.29217E-05 0.000124096
-0.00010584 1.29217E-05 -1.36459E-05 2.17228E-05
-0.000783119 0.000124096 2.17228E-05 0.000146912

Bon, au moins l'inversion fonctionne.
mais d'où peut venir la différence? Mes matrices sont en double, ce qui devrait être suffisant pour conserver la précision... non? Faudrait-il que je passe en long double? Y a-t-il plus grand encore? ^^

Cependant, pour la matrice:

29.1694 145.914 2.9161 0.0291432
145.914 729.909 14.5872 0.145783
2.9161 14.5872 0.291527 0.00291348
0.0291432 0.145783 0.00291348 2.91E-05

Le programme bug et renvoie:
Check failed in file c:\qtsdk\mingw\bin\../lib/gcc/mingw32/4.4.0/../../../../include/boost/numeric/ublas/lu.hpp at line 299:
detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cm2)
terminate called after throwing an instance of 'boost::numeric::ublas::internal_logic'
what(): internal logic
...c'est-à-dire au niveau du lu_substitute de InvertMatrix...

Une idée sur le problème?

Excel me renvoie:

29583.05966 -1036.301438 -110716.8234 -13342628.25
-1036.301438 263.7012475 -687.1073365 -214312.0502
-110716.8234 -687.1073363 1034457.92 10747589.39
-13342628.25 -214312.0502 10747589.4 13352224310

La matrice est donc bien inversible, par contre, on voit que l'étendue est très importante, +/- 10^10... Mais ce qui devrait largement tenir dans du double...

Encore une idée?


Merci!

[davcha]

Utilise plutôt ça : http://arma.sourceforge.net/

Sinon, l'explication vient sûrement de la différence entre les algos pour approximer la matrice inverse.

[Gm7468]

Merci pour ta contribution.

Si je comprends bien, ce que tu me proposes là, c'est de changer complètement tout mon code pour passer de boost à armadillo?

Bien que cela représenterait pas mal de travail, je n'ai rien contre. Cependant, est-ce qu'un peu de justifications supplémentaires pourraient être apportées?

Boost a l'air d'être éprouvé, alors pourquoi le balayer si vite?

Je n'ai trouvé aucune info supplémentaire sur cette bibliothèque lorsque je cherchais les inversions de matrice, est-elle peu connue? récente? moins efficace/fiable/... que Boost?

[Invité]

J'aimerais également te proposer gsl, que j'ai trouvé assez facile d'emploi, et suffisamment complet à mon gout.

Quant à l'erreur, comme dit par davcha les algo mis en jeu sont probablement différents.
Cela dit, le déterminant de

Code :

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

1
2
3
4
29.1694 145.914 2.9161 0.0291432;
145.914 729.909 14.5872 0.145783;
2.9161 14.5872 0.291527 0.00291348;
0.0291432 0.145783 0.00291348 2.91E-05

est -1.8620e-15 qui vaut quasiment 0, donc tu y vas un peu au bonheur la chance.

Un exemple simple:

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
A=[
29.1694 145.914 2.9161 0.0291432;
145.914 729.909 14.5872 0.145783;
2.9161 14.5872 0.291527 0.00291348;
0.0291432 0.145783 0.00291348 2.91E-05
]
invA=inv(A)
invInvA=inv(inv(A))
B=[
 29583.05966 -1036.301438 -110716.8234 -13342628.25;
-1036.301438 263.7012475 -687.1073365 -214312.0502;
-110716.8234 -687.1073363 1034457.92 10747589.39;
-13342628.25 -214312.0502 10747589.4 13352224310 
];
invB=inv(B)

et la sortie

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
A =
 
   2.9169e+01   1.4591e+02   2.9161e+00   2.9143e-02
   1.4591e+02   7.2991e+02   1.4587e+01   1.4578e-01
   2.9161e+00   1.4587e+01   2.9153e-01   2.9135e-03
   2.9143e-02   1.4578e-01   2.9135e-03   2.9100e-05
 
invA =
 
   1.6191e+04  -1.2514e+03  -9.9930e+04   5.8695e+04
  -1.2514e+03   2.6025e+02  -5.1384e+02   9.4276e+02
  -9.9930e+04  -5.1384e+02   1.0258e+06  -4.7279e+04
   5.8695e+04   9.4276e+02  -4.7279e+04  -5.8737e+07
 
invInvA =
 
   2.9169e+01   1.4591e+02   2.9161e+00   2.9143e-02
   1.4591e+02   7.2991e+02   1.4587e+01   1.4578e-01
   2.9161e+00   1.4587e+01   2.9153e-01   2.9135e-03
   2.9143e-02   1.4578e-01   2.9135e-03   2.9100e-05
 
invB =
 
   2.9167e+01   1.4590e+02   2.9158e+00   2.9140e-02
   1.4590e+02   7.2984e+02   1.4586e+01   1.4577e-01
   2.9158e+00   1.4586e+01   2.9150e-01   2.9132e-03
   2.9140e-02   1.4577e-01   2.9132e-03   2.9114e-05

on constate que l'inversion de B donne bien A, l'inversion de l'inversion de A donne bien A, mais l'inversion de A (invA) est différente de B...
bref, le comportement sur l'inverse est erratique...

[Gm7468]

Merci également galerien69 pour ta participation.

Je ne pourrai pas utiliser gsl car elle est GPL, j'avais oublié de le mentionner, mais je ne peux pas placer mon application en licence GPL, donc c'est exclu.
Armadillo est LGPL, ce qui serait acceptable.

Je suis d'accord au niveau de la précision, cela ne pose pour le moment pas le plus gros problème, et je verrai à terme comment je veux traiter ces différences, qui, si elles restent minimes, n'ont pas d'importance.
Par contre, le bug du programme est lui totalement bloquant, et il est donc primordial de le supprimer, en corrigeant Boost, ou en passant à autre chose. Mais le choix doit être justifié, et vous êtes là pour m’aiguiller et me conseiller

Merci!

La question reste ouvert...

[Joel F]

l'inversion de matrice c'est un truc instable numeriquement donc cooli sur "on change de lib" :o

ensuite, y a besoin de cette inverse ? en genral on resout Ax=b et on ne calcule *jamais* A-1

[Gm7468]

Il y a effectivement besoin de cette matrice inverse, à mon avis.
C'est pour une compensation par les moindres carrés:
J'ai une matrice A, et une matrice K.
La solution recherchée est:
dX = (At * A)^-1 * (At * K)

Soit N = At * A et nk = At * K,
on a donc dX = N^-1 * nk.

La matrice N^-1 sert donc pour ce calcul de dX, mais également pour d'autres calculs qui nécessitent la matrice Q ( = N^-1).

Si tu vois une solution pour ne pas calculer N^-1, je suis preneur, mais comme il me la faut à plusieurs endroits, elle sera nécessaire, calculée par l'inversion ou par un autre moyen...

[Invité]

La matrice N^-1 sert donc pour ce calcul de dX, mais également pour d'autres calculs qui nécessitent la matrice Q ( = N^-1).

Si tu vois une solution pour ne pas calculer N^-1, je suis preneur, mais comme il me la faut à plusieurs endroits, elle sera nécessaire, calculée par l'inversion ou par un autre moyen...

c'est rare qu'on ait à calculer explicitement une matrice inverse.

Ici, le cas trivial de dX, c'est que
dX = (At * A)^-1 * (At * K) = A^-1 At^-1 At K = A^-1 K
AdX=k
et tu cherches dX
Tu peux y aller à coup de LU, ou simplement à coup de gradient pour trouver dX.

Si tous tes systems sont linéaires, tu peux trouver ta solution X (ici dX) en une vingtaine d'itérations avec un gradient par exemple. Après il reste peut être à conditionner l'initialisation, mais globalement, tu es gagnant à ne pas calculer l'inverse explicitement.

Quels sont tes autres systèmes?

edit: je viens de comprendre que A est probablement non carrée, et que tu prends A^tA pour avoir une matrice carrée. Si c'est le cas, c'est évident que A^tA ne sera pas inversible (je me rappèle plus de la démo...). Tu peux regarder du côté des valeurs singulières.