Discussion :

[steph42]

Bonjour à tous,

je cherche des méthodes (algorithmes) de résolutions de systèmes (très largement) surdéterminés (de l'ordre de 300 000 équations pour 4 inconnues, chacune ayant son importance!!).

Je me sers de Matlab pour résoudre, et j'ai testé:
le \ (qui revient à résoudre par la minimisation de l'erreur des moindres carrés),
la svd, en fait spm_svd puisque svd seul ne fonctionne pas (out of memory...) (svd = single decomposition value), qui permet de passer de la matrice A (de taille (300000, 4) du système de départ à : A = U*S*V' où U(300000, p) V(300000, 4) et S(p, p) qui est une matrice diagonale avec les p valeurs propres sur la diagonale. J'ai pris pour seuil 1 (les valeurs inférieures à 1 sont mises à 0 lors du calcul des 3 matrices).

J'ai entendu parler de la méthode de régularisation de Tikhonov, qui dérive des moindres carrés, mais comment choisir l'opérateur de régularisation?

D'autre part, existe -t'il d'autres méthodes de résolutions de systèmes surdéterminés?

En effet, les 2 méthodes testées jusqu'à maintenant ne donnent que des valeurs bien évidemment "approchées", et je cherche la "meilleure", d'un point de vue numérique, mais aussi d'un point de vue "vitesse de calcul".

Merci à tous.
P.S: j'espère avoir été assez clair...

[souviron34]

a mon humble avis il y en a 2 grandes classes :

• les Simplex
• les Monte-Carlo

[Nemerle]

bonjour, je ne connais pas du tout ce type de sujet, mais cela m'intéresse. Quelles forment ont tes équations?

[steph42]

Re-bonjour,

merci pour ta réponse souviron. Je t'avoue que je ne vois pas trop comment m'en sortir avec de tels algorithmes... et surtout comment les mettre en oeuvre pour 300 000 équations...

Précision du problème:
les équations sont linéaires. On a tout bêtement un système matriciel de la forme:
Ax=b où je cherche x. si x=[x1 x2 x3 x4]
on a 300 000 fois (i désignant un numéro de ligne):
ai1*x1+ai2*x2+ai3*x3+ai4*x4 = bi

(Dans certains cas, je peux avoir jusqu'à 1 680 000 équations pour 4 inconnues)

est ce que d'un point de vue théorique on peut (étant donné que l'on se trouve dans le cas d'un système surdéterminé) utiliser la méthode de Gauss (substitution)? même question pour la factorisation LU (méthode de Crout), pour la méthode de Cholesky, pour la méthode de Householder?

Et si oui, connaissez vous l'algorithme de ces méthodes?
Merci beaucoup!

[Nemerle]

je vais poursuivre par une bétise: pourquoi ne suffit-il pas de résoudre les 4 premières équations (le système 4x4 donnant les xi de façon uniques), puis de checker les équations restantes?????

[steph42]

La méthode que tu proposes donne des résultats totalement faux.

En effet, chaque équation contient une information de localisation spatiale. Donc ne prendre que les 4 premières équations, ou les 4 dernières reviendrait à choisir un a priori sur la localisation des points qui suivront l'équation linéaire.

Je m'explique: ton modèle fonctionne parfaitement pour les 4 premiers points (ou 4 premières équations), et il ne fonctionne pas pour les 299 996 autres! et je t'assure que c'est le cas. Alors que faire? prendre au "pif" 4 équations et résoudre? le problème est le même: il y aura 299 996 points pour lesquels le système ne fonctionnera pas (certainement moins, mais il y aura énormément de points pour lesquels ça ne collera pas)

Il faut donc trouver le modèle le moins mauvais possible, et qui donne des résultats pas trop mauvais pour CHACUN des points...

J'espère que ce que j'ai dit est compréhensible...

[pseudocode]

je cherche des méthodes (algorithmes) de résolutions de systèmes (très largement) surdéterminés (de l'ordre de 300 000 équations pour 4 inconnues, chacune ayant son importance!!).

Et l'approche bete et mechante "Moindre Carrré avec Jacobienne du systeme" ? Ca ne fait qu'une matrice 4*4 a inverser....

http://www.developpez.net/forums/sho...d.php?t=373845

[steph42]

J'ai déjà fait par la minimisation des moindres carré, et par svd. Ce que je cherche, ce sont d'autres méthodes (algorithmes, cadre d'utilisation, et si possible code ou fonctions déjà existant).

En effet, les moindres carrés tournent très rapidement (je suis sous matlab), et les résultats sont meilleurs qu'avec la svd (fonction spm_svd).

Mais la question bête et méchante que l'on m'a posé: "pourquoi n'avez vous testé que ces 2 méthodes si les résultats ne sont pas suffisamment bons? "

Comme je suis bête et discipliné, je cherche donc d'autres méthodes.

[pseudocode]

Mais la question bête et méchante que l'on m'a posé: "pourquoi n'avez vous testé que ces 2 méthodes si les résultats ne sont pas suffisamment bons? "

Pas "suffisamment bons" par rapport a quoi ?

je suppose que les "equations" sont les individus d'une population de résultats d'une d'experience. Ou, en plus clair, des valeurs de y avec "y = f(x1,x2,x3,x4)".

Si les mesures "y=f(X)" sont exactes alors, comme le dit Nemerle, il suffit d'en prendre 4 au pif et de résoudre le systeme d'equation

Si les mesures "y=f(X)" ont une erreur, alors il faut choisir une strategie d'élimination ou de réduction de l'erreur. Les "moindres carrés" en sont une, mais si vous estimez que ce n'est pas "suffisamment bon" alors il vous faut modeliser ce qu'est une "bonne" solution.

[souviron34]

je dirais que la methode standard pour ce genre de probleme est le Simplex..

En gros, si l'on prend un exemple d'une courbe f(x,y) (donc representation dans un espace de 2 dims), cela revient a prendre 3 points "aleatoires" dans cet espace et a iterer (avec methode ) pour que les 3 points convergent vers le minimum.

Donc en generalisant a un espace de dimension N, on prend N+1 points...

+ "Simplex"

[Nemerle]

ok, je comprends (un peu mieux)... ton problème n'est pas un simple système à 4 inconnues...

Il me semble que tu cherches les "meilleurs" candidats qui satisfassent au mieux des 300 000 contraintes. Mais la question rigoureuse est: comment définis-tu la notion de satisfaction? Comme pseudocode, par une notion de distance? ou par une notion de diffusion de distribution? par une percolation? Par un faisceau?

"tout" est possible...

[steph42]

Bonjour,

merci à tous pour vos réponses.

En fait, pour faire simple: j'ai 4 vecteurs de 300 000 points. chacun étant "prédictif" (selon la littérature) d'un vecteur b. J'ai donc eu l'idée d'essayer d'estimer au mieux le vecteur b avec une combinaison linéaire des 4 vecteurs de départ. Je créé la matrice A = [V1 V2 V3 V4].

Les moindres carrés donnent une solution. Cependant, en calculant le b théorique à partir des 4 coefficients trouvés, on constate de grosses différences de valeurs avec le b mesuré.
Les résultats ne sont donc pas suffisamment bons par rapport aux "vraies" valeurs.

Il est très facile de trouver des coefficients pour mon GLM qui foncionnent pour un "certain" nombre de points, mais dans ce cas, l'erreur est telle sur les autres points, que la solution n'est pas acceptable au final.

Il faut donc trouver une méthode qui minimise l'erreur pour chacun des points, et dont (en plus) la résultante (la somme?) de ces erreurs soit aussi la plus minime possible. Les moindres carrés sont une solution, qui ne donne pas entière satisfaction. Je cherche donc des solutions mathématiques meilleures, s'il en existe.

Je vais me pencher sur l'algorithme du simplex, mais j'ai peur que l'erreur sur certains points fasse "exploser" l'erreur globale. Mais peut être que je n'ai pas compris la "philosophie" de la méthode.

[FR119492]

Salut !

Pour répondre tout d'abord à ton message de 17h48, les méthodes de Gauss, Crout et Cholesky sont inutilisables dans ton cas, car, par leur principe même, elles ne s'appliquent qu'à des matrices carrées. La méthode de Householder pourrait probablement aller, mais je ne l'aime pas (tu me diras que ce n'est pas un argument!).

Ma préférence va donc à la méthode SVD qui trouve, pour autant que je sache, la "solution" x telle que la norme euclidienne de A*x-B soit minimale. Avec des systèmes plus petits que le tien, ça marche très bien.

L'explication de ton message "out of memory" pourrait être la suivante: si tu as fait [U,S,V] = svd(A) pour une matrice A(300000,4), tu as tenté de construire une matrice U(300000,300000) qui ne tient évidemment pas dans la mémoire de ton ordinateur. Tu devrais essayer la version "économique, soit [U,S,V]=svd(A,0), qui te donnera une matrice U de la même taille que A.

Si ça ne résout pas ton problème, fais-le moi savoir, et on creusera un peu plus.

Bonne chance
Jean-Marc Blanc

[souviron34]

l'utilisation du simplex est justement exactement pour ton probleme :

En fait son fonctionnement est (en gros) le suivant :

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
 
Pour iter = 0  a L iter
 
    estimation des coeffs et des "ecarts"
 
    Pour i = 0 a N equations
 
        Pour j = 0 a N equations         
            resolution du systeme de m equations a m inconnues
            calcul de l'erreur pour les (N-m) equations
        Fin pour
 
        calcul de l'erreur globale pour cette equation i
 
    Fin pour
 
 
    Calcul de l'erreur globale de cette iteration
 
    Delta valeurs = Delta valeurs iter precedente ponderes par erreur - erreur precedente
 
Fin pour

Bon c'est du lourd, mais c'est tres utilise en physique (en particulier interactions a N corps) pour justement resoudre un grand nombre d'equations dites "avec contraintes" (un exemple dans mon domaine d'origine : une analyse spectrale d'un amas de 100 000 etoiles. On a par exemple 10 inconnues (10 longueurs d'ondes). On a 100 000 spectres, donc 100 000 equations, a 10 inconnues...)

[steph42]

Merci FR119492.

En fait, j'ai utilisé spm_svd, qui permet en outre de fixer un seuil pour la mise à 0 des éléments, lors de la diagonalisation. Ca me donne satisfaction. Je vais cependant tester ce que tu m'as conseillé, car svd et spm_svd ne sont peut être pas codés de la même façon.

Merci souviron34.

Je ne pense pas que les métodes de Monte Carlo soient les plus adaptées. Qu'en penses tu?

Pour ce qui est de l'algo du simplex, en effet, nos problèmatiques sont très proches. Mais je ne sais pas si j'ai bien compris ce que tu m'as dit:

On parcourt la matrice (de taille (300 000, 4)) et on résoud tous les systèmes de p(p=4 pour moi) équations à p inconnues. on calcule ce que donne le modèle théorique et on "note" l'erreur avec les valeurs mesurées. Ensuite, on choisit les coefficients correspondant à l'erreur la plus "petite".

Si c'est bien cela, j'ai une question. Lorsque l'on calcule les valeurs théoriques, on le fait pour les 4 équations en question, ou pour les 300 000?

Mon but en fait, et de trouver des coefficients qui permettent de bien approximer les valeurs mesurées, mais qui permettent aussi d'approximer les valeurs mesurées sur un autre nuage de points (donc sur 300 000 nouvelles équations).

J'espère être clair...

[souviron34]

Je ne pense pas que les métodes de Monte Carlo soient les plus adaptées. Qu'en penses tu?

En effet je pense qu'elles sont moins bien adaptees que le Simplex.. Mais mon utilisation du Simplex remonte a loin... En fait au debut de ma carriere

Pour ce qui est de l'algo du simplex, en effet, nos problèmatiques sont très proches. Mais je ne sais pas si j'ai bien compris ce que tu m'as dit:

On parcourt la matrice (de taille (300 000, 4)) et on résoud tous les systèmes de p(p=4 pour moi) équations à p inconnues. on calcule ce que donne le modèle théorique et on "note" l'erreur avec les valeurs mesurées. Ensuite, on choisit les coefficients correspondant à l'erreur la plus "petite".

Si c'est bien cela, j'ai une question. Lorsque l'on calcule les valeurs théoriques, on le fait pour les 4 équations en question, ou pour les 300 000?

Mon but en fait, et de trouver des coefficients qui permettent de bien approximer les valeurs mesurées, mais qui permettent aussi d'approximer les valeurs mesurées sur un autre nuage de points (donc sur 300 000 nouvelles équations).

J'espère être clair...

En fait, c'est a peu pres ca.. Ta question disparait a cause de l'algo en lui-meme :

et on résoud tous les systèmes de p(p=4 pour moi) équations à p inconnues. on calcule ce que donne le modèle théorique et on "note" l'erreur avec les valeurs mesurées. Ensuite, on choisit les coefficients correspondant à l'erreur la plus "petite".

on resoud toutes les combinaisons de p equations a p inconnues, en prenant comme premiere equation d'abord la premiere de toutes, puis ainsi de suite..

L'algo calcule alors le determinant (et donc la solution exacte) de ces p equations. Puis il applique les coefficients a toutes les equations restantes, et en tire un facteur d'erreur par equation. Une certaine formule donne la combinaison de ces facteurs d'erreurs comme le facteur d'erreur associe a l'equation prise comme premiere equation.

Puis on itere en changeant de premiere equation.

A la suite de ces N iterations, on a un tableau d'erreurs de N valeurs, permettant de determiner quelle est (dans cette iteration) la meilleure combinaison d'equations et donc les meilleures valeurs des coefficients.


La prise en compte de ce facteur d'erreur (qui en fait est un vecteur de p valeurs) et sa comparaison avec l'etape d'iteration precedente donne un vecteur "deplacement" de p valeurs, que l'on applique aux coefficients.

Et on re-itere...

[Nemerle]

a-y est, z'ai compris.... ça n'a effectivement rien à voir avec la résolution de systèmes linéaires.

Je vais réfléchir au problème mais quelques rqs:

1. Il est clair que puisque tu cherches une approximation, tu dois choisir ce qu'est ta fonction de performance.
2. Tu peux choisir effectivement les moindres carrés en brut, mais tu peux aussi à partir des moindres carrés choisir quelque-chose de plus précis: par exemple tu regardes la distribution de tes erreurs pour tes 300 000 equations, et tu demandes à ce que celle-ci soit la plus "large" possible ou encore la plus pointue. La méthode Simplex n'est qu'une méthode parmis d'autres: elle cherche à minimiser la norme quadratique par quadruplets de points de la distribution.
3. Tu peux changer de norme: plutôt que de prendre L2 (i.e. les moindres carrées), tu prends Lp avec p<>2, ou encore la norme Somme|x-xi|... encore une fois cela dépend de ce que tu veux.
4. Tu peux aussi faire une régression linéaire par bascule:

- tu prends les eq. simplifiées ai1x1=bi, et tu fais une régression sur x1 qui te donne une valeur val_x1
- tu prends les eq. simplifiées ai2x2=bi-val_x1, et tu fais une régression sur x2 qui te donne une valeur val_x2
- etc...

5. Il y a d'autres méthodes possibles pour ta fonction de performance, je crois me souvenir de trucs avec des espaces de kirnov (??)....


Mais finalement, tout revient encore une fois à choisir ta fonction de performance. Moi personnellement, je chercherais plus à avoir une forme de distribtion satisfaisante plutôt que de minimiser ma somme de distances au points...

PS: une fois la fct de perf choisie, tu peux t'amuser aussi à implémenter un algorithme génétique, cela marche très bien dans ce genre de problème...

[steph42]

Merci souviron34.

Alors je sais que je suis lourd, mais, encore une petite question...

Si j'ai 300 000 équations... Je fixe la première, puis je cherche toutes les combinaisons possibles de 3 équations parmi les 299 999 restantes? Ca me fait donc un 300 000 * 299 999 * 299 998 /6 possibilités... donc autant d'erreurs globales à calculer??
puisque l'on incrémente juste sur la première équation, ce nombre de possibilités (de choix d'équations parmi les restantes) va aller en diminuant, mais ça fait quand meme un nombre énorme d'erreurs à stocher, pour chaque iteration...?

Si c'est bien ça, au final je prends les 4 coefficients (p=4 pour moi!) qui me donnent les meilleurs résultats au final?

Merci nemerle.
Je suis d'accord avec toutes tes remarques... J'ajouterai qu'il y a encore la méthode des équations normales: on multiplie la matrice du système par sa transposée (et on fait pareil de l'autre côté), pour se ramener à la résolution d'un système carré.

J'ai trouvé un peu de doc sur la régularisation de Tikhonov. Mais je ne vois pas trop comment l'implémanter. Tout ce que j'ai trouvé est assez "théorique"...

Pour ce qui est de la régression à bascule, je n'ai pas trop compris en fait. J'ai
A=
a11 a12 ... a1n
a21 a22 ... a2n
. ... .
am1 am2 ... amn

x =
x1
x2
...
xn

b=
b1
b2
...
bm

Je résouds: d'abord:
a11x1 = b1
a21x1 = b2
...
am1x1 = bm
régression et je trouve une valeur "moyenne" de x1 en fait.
après je résouds:
a12x2 = b1-val_x1
a22x2 = b2-valx1
...
am2x2 = bm-val x1
régression et je trouve une valeur "moyenne" de x2?

après:
a13x3 = b1-valx1-valx2? ou b1-valx2?
...
am3x3 = b1-valx1-valx2 ? ou b1-valx2?

Merci!

Je tacherai après avoir tout programmé et testé, de vous dire ce qui marche le mieux, et ce qui est le plus rapide.

[Matthieu Brucher]

En fait, tu devrais plutôt écrire tout ça sous la forme d'une fonction de coût (ce que tu fait pour la résolution classique en minimisant un terme quadratique de l'erreur de reconstruction).
Tu peux prendre ensuite un critère semi-quadratique, une gaussienne généralisée, ou une norme 1 carrément, sachant que ce ne sont des fonctions que de 4 variables et qu'elles sont linéaires, c'est très rapide à optimiser.

[souviron34]

Si c'est bien ça, au final je prends les 4 coefficients (p=4 pour moi!) qui me donnent les meilleurs résultats au final?
.

voui en gros c'est ca.. Et l'algo (il y en a de deja tout codes) le fait pour toi..

Je me souviens qu'il y en a (en Fortran) dans la biblio. de maths du CERN (j'ai plus son nom, la, mais c'est quelque chose comme lib..ml ) mais on en trouve en C je crois.. regarde Google....

c'est vrai que ca prend du temps, mais on va quand meme plus vite qu'avant : quand j'avais tente ca, en 1982, sur un HP1000, 100 000 equations avec 10 inconnues, il fallait 1 an de CPU, et 1h de Cray.. Je pense qu'aujourdhui sur un PC moyen, cela devrait prendre de 3 a 5 heures...