2.3 Modèles non-linéaires: Méthode de Levenberg-Marquardt.

Minimisation d'une forme quadratique et d'une forme quelconque

Considérons tout d'abord une fonction scalaire $f$ quelconque, dépendante de plusieurs variables réelles, c'est-à-dire définie sur un espace affine $\vec{x}$, où les coordonnées sont exprimées dans un système d'origine $\vec{x_{0}}$, et écrivons le début du développement de Taylor de $f$ au voisinage de ce point origine:

$$f(\vec{x}) = f(\vec{x_0}) +\sum_i\frac{\partial f}{\partial x_i} (\vec{x_0}) x... ...i,j}\frac{\partial^2 f}{\partial x_i \partial x_j} (\vec{x_0}) x_i x_j + \ldots$$

$$\simeq c - \vec{b}\cdot \vec{x} + \frac{1}{2} \vec{x}{\bf H}\vec{x}$$ (2.27)

avec $c \equiv f(\vec{x_0})$, $\vec{b}\equiv -\vec{\nabla} f(\vec{x_0})$, et ${\rm H}_{ij}\equiv \frac{\partial^2 f}{\partial x_i \partial x_j}(\vec{x_0})$ est appelée la matrice hessienne ou hessien de $f$ en $\vec{x_0}$.

Lorsque l'approximation (2.27) est exacte, on dit que $f$ est une forme quadratique; son gradient s'en déduit immédiatement:

$$\vec{\nabla}f = {\bf H} \vec{x} - \vec{b}$$ (2.28)

et donc le gradient s'annulera -c'est-à-dire la forme atteindra un extremum- pour une valeur de $\vec{x}$ solution de : ${\bf H}\vec{x} = \vec{b}$.

Supposons maintenant que (2.27) soit seulement une (bonne) approximation de la forme $f$, on déduit facilement de (2.28) une formule directe pour converger d'un point initial $\vec{x_{0}}$ vers un minimum de $f$:

$$\vec{x}_{min}= \vec{x_{0}} - {\bf H}^{-1}{\vec{\nabla}f(\vec{x_{0}})}$$ (2.29)

En revanche, il se peut que (2.27) soit une mauvaise approximation de la forme que nous essayons de minimiser. Dans ce cas, tout ce qu'on peut tenter est de se déplacer d'un pas de longueur fixée arbitrairement dans la direction opposée au gradient (méthode des plus fortes pentes ou ``steepest descent method''); autrement dit:

$$\vec{x}_{1}= \vec{x_{0}} - {\rm Cte} \, {\vec{\nabla}f(\vec{x_{0}})}$$ (2.30)

Gradient et Hessien du $\chi ^{2}$

Si un modèle à ajuster est donné, pour un jeu de paramètres $(a_{1},a_{2},\ldots,a_{M})=\vec{a}$ choisis, par $y=y(x,\vec{a})$, la valeur de son $\chi ^{2}$ par rapport à $N$ mesures $y_{i}(x_{i})$ est une forme (non forcément quadratique) définie sur l'espace vectoriel de dimension $M$ des paramètres à ajuster par:

$$\chi^{2}(\vec{a})=\sum_{i=1}^{N}\left[ \frac{y_{i}-y(x_{i},\vec{a})} {\sigma_{i}}\right]^{2}$$ (2.31)

Les composantes du gradient du $\chi^2$ sont donc:

$$\frac{\partial \chi^{2}}{\partial a_{k}}=-2\sum_{i=1}^{N} \... ...igma_{i}^{2}} \frac{\partial y(x_{i},\vec{a})}{\partial a_{k}}$$ (2.32)

et le hessien s'en déduit par une dérivation partielle de plus, soit (en négligeant les termes de second ordre dans le modèle $y$):

$$\frac{\partial^{2} \chi^{2}}{\partial a_{k}\partial a_{l}}=2... ... a_{k}}\frac{\partial y(x_{i},\vec{a})}{\partial a_{l}}\right]$$ (2.33)

Pour simplifier les notations, il est utile de poser : $\beta_{k}\equiv -\frac{1}{2} \partial\chi^{2}/\partial a_{k}$ et $\alpha_{kl}\equiv \frac{1}{2} \partial^{2}\chi^{2}/\partial a_{k}\partial a_{l}$, ou vectoriellement $\vec{\beta}=-\frac{1}{2}\vec{\nabla}\chi^2$ et $[\alpha]=\frac{1}{2}{\bf H}$. Dans le contexte de l'ajustement aux moindres carrés, cette matrice $[\alpha]$ est appelée matrice de courbure.

Avec ces notations, et en posant $\delta \vec{a}=\vec{a}_{min}-\vec{a}_{0}$, autrement dit $\delta \vec{a}$ est l'incrémentation à réaliser sur le jeu de paramètres initial $\vec{a_{0}}$ dans le schéma de convergence (2.29), la minimisation d'un $\chi ^{2}$ supposé quadratique revient à résoudre le système linéaire:

$$[\alpha]\delta \vec{a}=\vec{\beta}$$ (2.34)

En revanche, lorsque le $\chi ^{2}$ est «loin» d'être une forme quadratique, on utilisera la méthode des plus fortes pentes (2.30) qui s'écrit simplement:

$$\delta \vec{a} = {\rm Cte} \, \vec{\beta}$$ (2.35)

Quand changer de méthode pour minimiser le $\chi ^{2}$? La stratégie de Levenberg-Marquardt

Levenberg et Marquardt ont proposé une méthode efficace et astucieuse pour passer continûment du schéma d'inversion du hessien à celui des plus fortes pentes. Ce dernier sera utilisé loin du minimum et on tend à lui substituer le schéma d'inversion du hessien au fur et à mesure que l'on approche du minimum. Cette méthode de Levenberg-Marquardt a fait ses preuves et fonctionne remarquablement bien pour des modèles et domaines de la physique fort variés, si bien qu'elle constitue désormais le standard pour résoudre les problèmes d'ajustement aux moindres carrés de modèles non-linéaires.

On peut brièvement décrire la méthode de Levenberg-Marquardt comme une stratégie de recherche du $\chi^2$ minimum utilisant au mieux les schémas (2.34) et (2.35), et cela grâce à deux idées déterminantes. La première idée aboutit à modifier le schéma des plus fortes pentes (2.35) en remplaçant la constante (le pas) par un vecteur dont on choisit judicieusement les composantes. On peut interpréter ce choix comme une "mise à l'échelle", pour chacun des paramètres, du pas que l'on va effectuer dans la direction du minimum du $\chi^2$. On réalise ce choix en remarquant que cette constante de proportion entre une dérivée par rapport à $a_{k}$ et une différence finie en $a_{k}$ a naturellement la dimension de $a_{k}^{2}$. Par ailleurs, on postule qu'un ordre de grandeur de cette constante peut être donné par une composante de la matrice de courbure $[\alpha]$; or la seule composante de $[\alpha]$ dépendante de $a_{k}$ qui ait la dimension requise est $1/\alpha_{kk}$, et le schéma (2.35) doit donc être modifié pour s'écrire en composantes:

$$\delta a_{l}= \frac{1}{\lambda \alpha_{ll}} \beta_{l}$$ (2.36)

où $\lambda$ est un facteur $>1$ permettant de réduire globalement (et non composante par composante) le pas si celui-ci s'avérait trop grand (comme cela se fait dans la méthode des plus fortes pentes).

La deuxième idée consiste alors à poser $[\alpha]'= [\alpha] + \lambda {\bf I_{d}}$, où ${\bf I_{d}}$ est la matrice identité $M\times M$. Les deux schémas (2.34) et (2.36) sont avantageusement remplacés par l'unique formulation: trouver l'incrémentation $\delta \vec{a}$ solution du système

$$[\alpha]' \delta \vec{a}=\vec{\beta}$$ (2.37)

Lorsque $\lambda$ est grand, la matrice $[\alpha]'$ est à diagonale dominante et le système précédent équivaut à (2.36); lorsqu'au contraire $\lambda$ tend vers 0, ce système équivaut à (2.34). ^2.6

Incertitudes sur les paramètres ajustés et matrice de covariance

Lorsqu'on a trouvé un minimum acceptable (voir calcul de confiance section 2.1) du chi-carré pour un jeu de $M$ paramètres $\vec{a}_{min}$, la variation de $\chi^2$ autour de ce minimum $\chi^2_{min}$ pour une variation $\delta \vec{a}$ des paramètres ajustés est donnée par: (en appliquant l'équation (2.27) au $\chi^2$ et puisque $\vec{\nabla}\chi^2(\vec{a}_{min})=0$)

$$\chi^2= \chi^{2}_{min}+\delta\vec{a}\,[\alpha]\,\delta\vec{a}$$ (2.38)

On va s'intéresser en particulier à la variation du $\chi^2$ lorsqu'on fait arbitrairement varier un seul paramètre $a_1$ , les autres paramètres restant fixés à leurs valeurs ajustées de $\vec{a}_{min}$. Notons $\chi^2_{M-1}$ le moindre chi-carré à $M-1$ degrés de liberté obtenu en fixant le paramètre $a_1$ à sa valeur arbitraire et soit $\vec{a}$ le nouveau jeu de paramètres qui minimise ce chi-carré. Posons $\Delta \chi^2_{1}\equiv \chi^2_{M-1}-\chi^2_{min}$ et $\delta \vec{a}=\vec{a}-\vec{a}_{min}$ (remarquons qu'aucune des composantes de $\delta \vec{a}$ n'est nulle a priori). On montre que ce $\Delta \chi^2_{1}(a_{1})$ est distribué comme le carré d'une variable aléatoire à distribution normale^2.7. Autrement dit, on aura formellement $\Delta \chi^2_{1} < 1$ pour $\delta a_{1} < 1 \sigma$(68.3% des cas), $\Delta \chi^2_{1} < 4$ pour $\delta a_{1} < 2 \sigma$(95.4% des cas), $\Delta \chi^2_{1} < 9$ pour $\delta a_{1} < 3 \sigma$(99.73% des cas), etc ...

On peut par ailleurs relier l'incertitude $\delta a_{1}$ sur le paramètre $a_{1}$ à $\Delta \chi^2_{1}$ en remarquant que $\vec{\nabla}\chi^2(\vec{a})=0$ sur toutes ses composantes sauf la première, et comme, d'après (2.34) $[\alpha] \delta \vec{a}=-\frac{1}{2}\vec{\nabla}\chi^2$, on aura, en posant comme en section 2.2.2 la matrice de covariance ${\bf C}=[\alpha]^{-1}$,

$$\delta a_1=-\frac{1}{2}\frac{\partial \chi^2}{\partial a_1}{\rm C_{11}}$$ (2.39)

On déduit de (2.39) et (2.38):

$$\Delta \chi^2_1= \delta\vec{a}\,[\alpha]\,\delta\vec{a}=(\delta a_{1})^{2}/{\rm C_{11}}$$ (2.40)

soit

$$\delta a_{1} = \pm \sqrt{\Delta \chi^2_1}\sqrt{{\rm C_{11}}}$$ (2.41)

et, donc, si l'on définit l'incertitude sur le paramètre $a_{l}$ par $\delta_{l}\equiv 1\sigma$, soit $\Delta \chi^2_{l} \equiv 1$, on aura finalement (comme dans le cas linéaire) :

$$\delta a_{l} = \sqrt{{\rm C_{ll}}}$$ (2.42)