Discussion :

[Dam2227]

Bonjour.

Je fais de la spectroscopie et pour analyser des spectres, je dois construire un hamiltonien qui reproduit les niveaux énergétiques de ma molécules et ainsi calculer les constantes fondamentales de cette molécules.

Le Hamiltonien est ici une matrice d'une certaine taille et j'utilise la méthode de Gauss-Newton pour fitter mes données.

Le programme marche a priori bien. Je l'ai testé en simulant des transitions (des raies d'un spectres) d'une molécule dont les constantes fondamentales sont déjà connues et il me donne les résultats attendus.

Malheureusement, j'ai un petit problème. En effet, lorsque je fitte mes données, il est instable. Je m'explique. Le jeu de constantes converge vers la solution en quelques itérations et juste avant d'arriver au jeu de valeurs finales (pour lesquelles ma RMSE est minimale), une des constante devient instable et 'saute' vers une autre valeur qui est très différente de la valeur attendue. Voici un exemple :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
 
rmse=0.063452	1
rmse=0.039652	2
rmse=0.010345	3
rmse=0.003369	4
rmse=0.002870	5
rmse=0.002838	6
rmse=0.032519	7
rmse=0.008480	8
rmse=0.003243	9
rmse=0.002861	10

Vous pouvez voir qu'il y a un saut entre les itérations 6 et 7. Et ce saut se répète continuellement : le programme reconverge toujours vers la même rmse (vers les 0.002834) puis diverge d'un coup etc...

Maintenant, comme je l'ai dit précédemment, je sais exactement qu'elle constante me pose problème. Malheureusement, je ne peux pas m'en passer, car elle a un sens physique important. Au final, le programme a été testé et devrait bien se comporter. Je me demandais donc s'il était possible que ce problème soit purement numérique et si oui s'il y a des solutions qui peuvent être envisagées pour y remédier.

Merci d'avance!

[souviron34]

je dirais qu'il y a plusieurs origines possibles :

• il peut s'agir de l'algo d'ajustement de la valeur : une fois l'écart calculé, comment est ajustée la valeur ?
  
  
• il peut s'agir de l'algo d'ajustement de l'ensemble des valeurs : comment est dispatché l'écart, ou comment est réparti la correction à appliquer ?
  
  
• il peut s'agir de la signifcation physique (ou de l'équation la décrivant) ?

J'aurais tendance à dire qu'il y a un bug dans l'ajustement de la valeur...


Mais sans en savoir plus sur l'équation de départ, j ne peux en dire plus..

[Nebulix]

C'est quoi la methode de Gauss Newton ?
Si elle utilise des interpolations, il est possible qu'à petite échelle le bruit sur les données produise une divergence.
J'utilise le "simplexe" de Nelder-Mead, je te le conseille.

[Dam2227]

Alors, la méthode de Gauss Newton est très bien expliquée à cette page : http://fr.wikipedia.org/wiki/Algorithme_de_Gauss-Newton.

Je ne suis pas concerné par le bruit dans mon problème. En effet, je récupère l'énergie de la transition (ce que je vois sur le spectre) en allant chercher la valeur à la main. Cette transition est le résultat d'une transition entre deux niveaux énergétiques. Le but du programme est de calculer ces niveaux énergétiques en fonction des constantes fondamentales de la molécule.

Je calcule dans un premier temps le jacobien du système (dérivée de l'énergie par rapport aux différents paramètres : J). Ensuite, je calcule l'ajustement en utilisant une fonction de résolution de problème linéaire (j'utilise Matlab) :

Code :

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

1
2
 
Delta = J\r;

où r est le résidu (transitions observées - transitions calculées).
Toutes mes constantes sont dans le vecteur appelé a et je l'incrémente simplement en ajoutant Delta :

Code :

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

1
2
 
a=a+Delta;

Je ne pense pas que ce soit la signification physique qui bug ou l'équation le décrivant. Le programme a été tester avant et de plus, l'équation physique le décrivant provient d'un article publié il y a des décennies, utilisé par pas mal de monde.

L'opérateur \ de matlab qui permet de faire des résolution de problème linéaire est supposé robuste et je n'ai pas de souvenirs d'avoir eu des bugs. De plus, j'ai aussi utilisé cette ligne de code :

Code :

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

1
2
 
Delta = inv(((J.')*w*J))*(J.')*w*r;

qui fait la même chose que \, mais en code. Et j'ai le même résultat. (w est le poids associé à chaque mesure. Dans mon cas, toutes les mesures ont le même poids donc utiliser \ ou cette formule doit revenir au même, ce que je constate bien).

Maintenant, je pourrais toujours programmer de la base un algo résolvant des équations linéaires. Mais est-ce que ça en vaut vraiment la peine?

[Jerome Briot]

Juste une remarque de terminologie :

Code :

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

1
2
 
Delta = inv(((J.')*w*J))*(J.')*w*r;

qui fait la même chose que \, mais en code.

Pour être exact, l'opérateur backslash \ (qui est en gros un des fondements de MATLAB) détermine l'algorithme optimale en fonction du type de système linéaire à résoudre (analyse de la matrice "à inverser) et choisi une méthode (LU, QR, Cholesky,...)

=> How does the backslash operator work when A is full?

On ne peut donc pas réellement dire que \ soit l'équivalent de inv (qui calcul l'inverse de la matrice du système.

[Nebulix]

Si tu utilises une "boîte noire" et que ça ne marche pas, il n'y a pas grand chose à faire...

Maintenant, je pourrais toujours programmer de la base un algo résolvant des équations linéaires. Mais est-ce que ça en vaut vraiment la peine?

Je me suis donné cette peine très relative (recopier quelques dizaines de lignes de code) il y a 25 ans.
Depuis je n'ai plus de problème de ce type...
à toi de voir

[Dam2227]

Pardon Dut. Je voulais juste pointer le fait que les deux méthodes me donnaient le même résultat.

Je crois qu'au final je vais essayer de coder moi-même quelque chose et voir ce que ça me donne parce que j'ai toujours l'impression que c'est une erreur purement numérique (d'après tous les tests que j'ai bien pu faire).

Merci en tout cas pour vos suggestions. Si vous en avez d'autres, d'hésitez pas^^.

[Dam2227]

Bonjour.

Alors voilà, j'ai des nouvelles concernant cette instabilité dans le fit. Apparemment, elle ne viendrait pas du programme lui-même, mais bien de la physique du système. Après avoir examiné mon jacobien, j'ai remarqué que la dérivée liée à la constante instable produit des valeurs très faibles au fur et à mesure que le fit converge. En fait, ça revient à dire que le fond du puit de potentiel (par rapport à cette constante) est très plat. Une fois arrivé en bas, la dérivée ne sait plus trop vers où pointer et donc diverge (fait de grands sauts). Une solution consisterait à calculer les dérivées secondes (bien que dans mon cas ce soit très très pénibles car j'ai de nombreux paramètres), ou bien je peux fixer la constante sur une valeur qui minimise le fit et lui associer une barre d'erreur liée à la rmse finale (par exemple, quand le 4ème chiffre significatif de la rmse change).

Voilà voilà. Le point positif dans tout ça, c'est que j'ai essayé les méthodes de SVD, LU et QR. Je n'ai pas perdu mon temps .

Merci pour vos répondes!

[souviron34]

ben sans calculer une dérivée seconde mais l'approximer tu peux calculer une pente, et détecter quand elle tend vers 0..

[Dam2227]

Bonjour souviron34. Je sais quand la pente tend vers 0 car il suffit que je regarde dans le Jacobien. Je pensais à la dérivée seconde dans le but de pouvoir empêcher les valeurs de mon jacobien (colonne correspondant à la dérivée par rapport à la constante problématique) de tendre vers 0. Au final, je pense que le plus simple sera de fixer la constante et essayer de trouver son incertitude par une méthode directe quelconque.

En fait, même quand je détecte que la pente tend vers 0, arrêter le programme ne va pas me donner la valeur minimale de la rmse. Tout simplement parce que le fit diverge sur les bords du fond plat et donc ne va jamais au milieu de ce fond (qui correspond à la valeur de ma constante qui minimise mon système). Il faut donc que je dise au programme, à la main, quelle doit être la valeur de ma constante (à la barre d'erreur près).

Voilà voilà, je ne sais pas si je suis très clair dans mes explications.

++

[souviron34]

Nul n'est besoin d'atteindre un hypothétique "milieu" d'un fond plat..

Un simple seuil suffit...

Code :

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

1
2
3
4
5
pente = Delta (valeurT2 - ValeurT1)
 
si pente < seuil
  on a atteint le minimum
fin si

Je disais ça pour éviter un calcul de maths qui peut être complexe.

Mon approche est physique.


Ce qui revient à ce que tu dis, mais sans le formaliser comme un Jacobien ou autre chose..

Comme dans toute modélisation physique, les équations ne font que décrire la théorie..

La pratique est légèrement différente.


En ce sens, toute méthode d'approximation est bonne, ce qui est couvert par le terme "fit", ou plutôt "best fit".


A ce moment-là, soit la méthode ce calcul est stable, soit elle est instable.

Il semble que pour ton problème la méthode utilisée soit instable.

Par conséquent il faut en changer.


Que ce soit via un moindre carrés, une recherche dichotomique, une convergence, etc, ..

De toutes façons, tu auras forcément, vu que tu travailles avec de la réalité et non une théorie, quelque part une valeur "en dur" (éventuellement relative) pour arrêter le calcul..

[Dam2227]

Salut souviron34. Oui, je comprend bien ce que tu veux dire et je suis d'accord avec toi. Dans tous les cas, le but, est de modifier le programme pour pouvoir prendre en compte ce problème. Maintenant, il reste à savoir quelle méthode est la meilleure. Il va falloir que j'en discute avec mon superviseur et voir ce qu'il en pense. Quand je parlais de fixer la constante, c'est juste que lorsque je la fixe, je peux vraiment trouver la rmse minimale. L'inconvénient, c'est que je suis obligé d'intervenir. Dans la solution que tu proposes, je n'atteindrais jamais cette rmse minimale (parce que ça diverge avant), mais là, le programme se débrouille tout seul. Reste à savoir ce qui est le plus convenable de faire dans de telles circonstances et là, je m'en remet au patron .

Merci pour tout en tout cas!