Discussion :

[Vaffreingue]

Bonjour, je développe un programme de calcul matriciel.
Je m'aide de maple pour vérifier mes resultats.
Cependant, pour certains calculs j'obtiens des résultats pour le moins étranges.

Petit exemple :

A=
0 1 1 1
0 -2 -1 1
0 -1 1 0
-1 1 0 -1

Inv(A) =
-0.4 -0.6 -0.2 -1
0.2 -0.2 -0.4 0
0.2 -0.2 0.6 0
0.6 0.4 -0.2 0

Si je fais les opérations à la main je trouve bien :

A.Inv(A) =
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1

par exemple :

(1;1) => 1 = 0*(-0.4) + 1*(0.2) + 1*(0.2) + 1*(0.6)
(3;2) => 0 = 0*(-0.6) + (-1)*(-0.2) + 1*(-0.2) + 0*(0.4)

Pour l'instant pas de soucis
Sauf que quand c'est mon programme qui fait le calcul j'obtient :

A.Inv(A) =
1 0 5.55112e-17 0
0 1 -5.55112e-17 0
0 0 1 0
0 5.55112e-17 0 1

et sur maple :


Donc, première conclusion : Vive le calcul à la main
Et deuxième... Qu'est ce qu'il se passe ici ?

Je dois avouer que je ne sais pas d'ou viens l'erreur, si nécessaire je peux montrer mon code, mais je ne suis pas sur que l'erreur vienne de celui-ci ...

Enfin si quelqu'un a une idée, je suis preneur !

Merci d'avance

PS: je viens de voir que pour des matrices plus petites ça ne le fait jamais et la probabilité d'erreur augmente avec la taille de la matrice ...

[gbdivers]

Bonjour

Il n'y a pas d'erreur (a priori), c'est juste un problème de représentation des nombres réels en informatique. Lorsque tu as 5.55112e-17, c'est à dire 0,000000...000555 (avec 17 zéros), on peut considérer que ça fait... 0
Donc A.inv(A) est bien :

1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1

Le mieux est probablement d'arrondir les chiffres lors de l'affichage à quelques décimales (1 ou 2)

[Rewpparo]

Il n'y a pas d'erreur, mais les types a virgule flottante (float, double) ne sont pas précis. Tu auras régulièrement des erreurs dans le genre, il faut savoir les anticiper, elles sont inévitables.
Mais franchement, a 10^-17, ca devrait le faire.

[Vaffreingue]

D'accord, dans ce cas, dans ma classe, il est préférable que je fasse l'arrondi de quelle manière ?

Edit: En plus, je ne vois pas pourquoi maple fait aussi l'erreur, alors que ce logiciel fait des calculs très précis je trouve aussi des approximations de ce genre ...

[Ekleog]

Plusieurs solutions :

Soit un cast vers l'entier le plus proche, si il est vraiment proche :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
double my_round(double value) {
   static const double epsilon = 1e-16; // A definir comme tu preferes
   const int rounded = static_cast<int>(value);
   if (v - rounded > 0.5) ++rounded;
   // rounded contient l'entier le plus proche. Il est aussi possible d'utiliser std::round de <cmath> pour faire ça, pour ce message je n'ai simplement pas pris le temps de vérifier le prototype exact.
   double error = rounded - value;
   if (-epsilon < error && error < epsilon) // Si on est proche de l'entier le plus proche
      return rounded; // On arrondit vers lui
   else // Sinon
      return value; // On ne change rien
}

Soit on fait du calcul formel, qui a l'avantage d'être infiniment précis (autant que le calcul à la main sans arrondi, quoi) -- mais là, bon courage !

[gbdivers]

Lors de l'affichage avec cout, tu peux préciser la précision souhaitée :

Code :

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

1
2
3
4
double f = 3.14159;
cout.unsetf(ios::floatfield);            // floatfield not set
cout.precision(5);
cout << f << endl;

voir http://www.cplusplus.com/reference/i...ase/precision/

Plusieurs solutions :

Soit un cast vers l'entier le plus proche, si il est vraiment proche :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
double my_round(double value) {
   static const double epsilon = 1e-16; // A definir comme tu preferes
   const int rounded = static_cast<int>(value);
   if (v - rounded > 0.5) ++rounded;
   // rounded contient l'entier le plus proche. Il est aussi possible d'utiliser std::round de <cmath> pour faire ça, pour ce message je n'ai simplement pas pris le temps de vérifier le prototype exact.
   double error = rounded - value;
   if (-epsilon < error && error < epsilon) // Si on est proche de l'entier le plus proche
      return rounded; // On arrondit vers lui
   else // Sinon
      return value; // On ne change rien
}

Soit on fait du calcul formel, qui a l'avantage d'être infiniment précis (autant que le calcul à la main sans arrondi, quoi) -- mais là, bon courage !

Il vaut mieux utiliser les fonctions d'arrondis, non ?
floor et ceil

[Vaffreingue]

Ekleog : peux-tu détailler un peu cette fonction s'il te plait ? Parce que j'aimerai n'utiliser que des fonctions "handmade"

gbdivers : en fait c'est surtout la valeur qui m'interesse, l'affichage après si la valeur est "exacte" (enfin si ça existe ) devrait aller

[Flob90]

Bonjour,

1/ Les arrondis ca se fait uniquement sur l'affichage. Les incertitudes de mesure/calcul ca se détermine/propage via des forumles, si tu commences à arrondir (*) de partout ca n'est plus possible.

2/ A l'affichage tu peux décider d'arrondir à une précision donnée, mais utiliser un affichage du type A +/- delta, est tout aussi correct.

3/ Si tu veux faire quelque chose d'efficase tu n'auras pas le choix de passer par les expressions templates, et avec ca ton problème va se simplifier.

Avec les expressions templates tu vas obtenir un arbre représentant ton expression, arbre que tu évalueras avec () (en général), il est (probablement) possible lors de l'évaluation de vérifier si il y a une expression du type mult<T,inv<U> >, et si c'est le cas de comparer (à epsilon près avec un epsilon bien déterminé : propagation des erreurs) l'évaluation de l'expression T et de l'expression U, et si elles sont "égales" de retourner l'identité.

Je vais essayer de poster un code "from scratch" qui montre l'idée, mais tu peux déjà aller regarder du coté de boost::proto, il donne des exemples liés au calcul algébrique.

(*) Pour être exact, la valeur retenue doit être cohérante avec l'incertitude, ca peut amener certains arrondies/troncatures, mais on arrondie pas dans le but de compenser l'incertitude.

PS: Personnelement, si tu débutes en C++ et que tu n'ambionnes pas quelque chose de très efficase, je me contenterais de garder dans ma classe matrice quelque chose représentant l'incertitude, et ferais une fonction d'affichage à une précision donnée ou alors un affichage de type A +/- delta.

[Vaffreingue]

Flob90 : c'est vrai, vu que je débute je n'ai pas compris grand chose aux templates
je pense que je ne vais pas me compliquer trop la vie, un e-17 ça reste potable et au moment de mon affichage je me débrouillerai

cela dit, par rapport à la propagation d'incertitude je viens de voir que pour une matrice 10*10 j'avais des "0" estimés qui valaient après calcul à 0,7 ou un chiffre beaucoup plus grand qu'avant... ce qui est très facheux!

[Flob90]

C'est étrange, les incertitudes grandisent en se propagant, mais pour passer de 10^-17 à 0.1, il faut énormément d'opérations. Mapple aussi te donne 0.7 ? Dans ce cas c'est peut-être toi qui a calculé un 0 alors que ca n'en était pas un.

[Vaffreingue]

Oui dans l'ensemble Maple me donne a peu près les mêmes valeurs que mon programme, je suis en train de le retaper voir si je n'ai pas laissé une instructions qui merdouille quelque part.

En fait comme le calcul matriciel ça peut se faire en récursif pour certaines opérations il arrive de trouver des calculs assez longs en effet

[Flob90]

Oui, sauf que pour passer de 10^-17 à 0.1, il faut, environ, 10^32 opérations, je doute que tu fasses autant d'opérations pour calculer un coefficient de la matrice.

[Vaffreingue]

Pas faux, donc c'est bien une erreur ailleurs ! Je vais revérifier tout ça, merci de me l'avoir montré !

edit : j'ai fait un ptit test des valeurs erronées :
-matrice 3*3 environ 1 fois sur 6 j'ai un e-15
-pour 4*4 quasiment tout le temps du e-13
-pour 5*5 tout le temps aussi du e-12
-pour 6*6 du e-10
-pour 7*7 du e-7
-idem pour 8*8
-pour 9*9 du e-6
-pour 10*10 du e-5

-pour 11*11 ma première ligne est :

-1.42729e+10
0.000221252
-4.24228e+10
1.42124e+11
-2.16188e+10
-6.07565e+09
4.87028e+10
-5.93413e+09
-1.61922e+09
-2.31562e+10
5.9881e+09

le premier terme est bon, mais les autres devraient etre nuls, et là carrément non

[Jean-Marc.Bourguet]

Il y a des algo qui se comportent mieux. Cherche un cours d'analyse numérique (qui est la branche qui étudie les algos de calculs numériques, en particulier l'effet des imprécisions).

[Invité]

-pour 11*11 ma première ligne est :

(...)

le premier terme est bon, mais les autres devraient etre nuls, et là carrément non

Que calcules tu exactement?

Si tu fais des sommes et des produits matriciels, sur des matrices à coefficients bornés, les résultats que tu observes doivent tenir à un bug. L'erreur d'arrondi n'augmente pas linéairement avec le nombre de calculs, mais elle est toujours nettement inférieure à la valeur absolue des plus grands coefficients de tes matrices.

Les choses peuvent très mal tourner quand tu fais des calculs plus complexes (tout ce qui contient un mélange d'additions et de divisions : recherches de valeurs propres, calcul de pivot, ou inversion). Dans ce cas, les erreurs peuvent s'amplifier et devenir énormes.

En général, cela se produit parce que tes matrices sont mal conditionnées (c'est à dire que le rapport entre leur plus petite et leur plus grande valeur propre est grand). Ca devient intenable quand on a des matrices non inversibles (ou quasiment...). Dans ce cas, le conditionnement est infini.

Le problème est malheureusement très courant avec des matrices "de la vraie vie", par exemple des variances-covariances statistiques, qui sont toujours mal conditionnées (sinon, cela voudrait dire que tous les axes de la matrice sont indépendants, et celle ci n'aurait aucun intérêt). Si tu développes un programme de calcul matriciel, tu vas le rencontrer assez vite.

Une "astuce" permettant de stabiliser les calculs consiste à ajouter à ta matrice de départ une diagonale dont les termes sont très petits (1e-6 par exemple). Celle ci ne perturbera pas les résultats, mais stabilisera le calcul. Sur de petites matrices (quelques dizaines de lignes et de colonnes), c'est généralement suffisant.

Francois

[Vaffreingue]

Bonjour !

Pour l'instant, je vais m'en tenir à cette solution : après chaque opération je vais tronquer la valeur à e-6 et puis normalement ça devrait suffir pour l'instant

Merci à vous pour vos réponses !

[Jean-Marc.Bourguet]

Bonjour !

Pour l'instant, je vais m'en tenir à cette solution : après chaque opération je vais tronquer la valeur à e-6 et puis normalement ça devrait suffir pour l'instant

Merci à vous pour vos réponses !

Pas une bonne idee du tout ca va vraisemblablement te causer d'autres problemes (p.e. tu ne vas pas etre capable de traiter une matrice dont tous les termes sont petits par rapport a 1e-6 et pourtant suffisemment grand pour ne jamais poser de probleme si on traite les choses comme il faut).

[Aleph69]

Bonsoir,

il faut que tu choisisses si tu veux faire du calcul formel (précis mais lent) ou du calcul numérique (rapide mais imprécis). Dans le second cas, tu ne dois jamais arrondir toi-même les calculs. Si tu considères que la précision de tes calculs n'est pas suffisante, deux solutions s'offrent à toi :
1. modifier l'algorithme pour améliorer sa stabilité,
2. changer de format de flottant
Deux remarques par ailleurs :
i. on n'inverse jamais une matrice : on la factorise puis on inverse les facteurs,
ii. le résultat obtenu avec la matrice 11x11 est suspect. Es-tu sûr que ta matrice est bien inversible?
N'hésite pas à fournir ta matrice et le code source si besoin.