Moindres carrés non linéaires

Les moindres carrés non linéaires est une forme des moindres carrés spécialisée dans l'estimation d'un modèle non linéaire en n paramètres à partir de m observations (m > n). Une façon d'estimer ce genre de problème est de considérer des itérations successives se basant sur une version linéarisée du modèle initial.

La théorie

Considérons un jeu de m couples d'observations, $(x_1, y_1), (x_2, y_2),\dots,(x_m, y_m),$ et une fonction de régression du type $y=f(x, \boldsymbol \beta),$. Cette fonction dépend des explicatives x mais aussi du vecteur des n paramètres $\boldsymbol \beta = (\beta_1, \beta_2, \dots, \beta_n),$ avec $m\ge n.$ On souhaite trouver le vecteur de paramètres $\boldsymbol \beta$ qui ajuste au mieux les données, au sens des moindres carrés:

$S=\sum_{i=1}^{m}r_i^2$

est minimisée en $\boldsymbol \beta$, où les résidus r_i sont donnés par

$r_i= y_i - f(x_i, \boldsymbol \beta)$

pour $i=1, 2,\dots, m.$

Le minimum de la somme des carrés des résidus S est atteint lorsque le Gradient s'annule (condition nécessaire). Puisque le problème est formulé avec n paramètres, il y a donc n équations normales:

$\frac{\partial S}{\partial \beta_j}=2\sum_i r_i\frac{\partial r_i}{\partial \beta_j}=0 \ (j=1,\ldots,n).$

Dans un système non linéaire, les dérivées $\frac{\partial r_i}{\partial \beta_j}$ dépendent aussi bien des variables explicatives que des paramètres: il faut donc renoncer à résoudre les équations normales aussi simplement que dans le cas linéaire. On a alors recours à une résolution numérique, à l'aide d'un procédé itératif

$\boldsymbol{\beta}^{k+1} = \boldsymbol{\beta}^k + \Delta \boldsymbol{\beta}.$

qui fournit des approximations successives $\boldsymbol{\beta}^k$ de plus en plus proches de la vraie valeur (inconnue) des paramètres, $\boldsymbol{\beta}_0$.

À chaque itération, le modèle initial est linéarisé par un développement de Taylor autour de $\boldsymbol{\beta}^k$ comme suit:

$f(x_i,\boldsymbol \beta_0) \approx f(x_i,\boldsymbol \beta^k) +\sum_j \frac{\partial f(x_i,\boldsymbol \beta^k)}{\partial \beta_{0,j}} \left(\beta_{0,j} -\beta^{k}_j \right) \approx f(x_i,\boldsymbol \beta^k) +\sum_j J_{ij} \Delta\beta_j.$

La Matrice jacobienne, J, dépend des données et de l'approximation en cours, aussi change-t-elle d'une itération à l'autre. Ainsi, en terme du modèle linéarisé, $\frac{\partial r_i}{\partial \beta_j}=-J_{ij}$ et les résidus sont donnés par

$r_i=\Delta y_i- \sum_{j=1}^{n} J_{ij}\Delta\beta_j; \ \Delta y_i=y_i- f(x_i,\boldsymbol \beta^k).$

Les équations normales deviennent

$-2\sum_{i=1}^{m}J_{ij} \left( \Delta y_i-\sum_{s=1}^{n} J_{is}\Delta \beta_s \right)=0$

ou encore

$\sum_{i=1}^{m}\sum_{s=1}^{n} J_{ij}J_{is}\Delta \beta_s=\sum_{i=1}^{m} J_{ij}\Delta y_i \; (j=1,n).\,$

Matriciellement, on en arrive à

$\mathbf{\left(J^TJ\right)\Delta \boldsymbol \beta=J^T\Delta} y.$

La linéarisation permet donc d'écrire:

$\boldsymbol{\beta}^{k+1} = \boldsymbol{\beta}^k + \left(\mathbf{J^TJ}\right)^{-1} \mathbf{J^T\Delta} y.$

Il faut remarquer que l'ensemble du terme de droite dépend seulement de l'itération en cours, à savoir $\boldsymbol{\beta}^k$, et permet donc de trouver la prochaine itération $\boldsymbol{\beta}^{k+1}$.

On peut aisément généraliser l'approche précédente en considérant une somme pondérée des carrés des résidus

$S=\sum_{i=1}^{m}W_{ii}r_i^2.$

Idéalement, chaque élément de la matrice diagonale de pondération W devrait être égal à l'inverse de la variance de l'observation ^[1] Les équations normales deviennent alors

$\mathbf{\left(J^TWJ\right)\Delta \boldsymbol \beta=J^TW\Delta y}$

ce qui procure la base de l'algorithme d'optimisation de Gauss-Newton.

Différences entre les moindres carrés linéaires et non-linéaires

Il y a de nombreuses divergences entre les moindres carrés linéaires (MCL) et non-linéaires (MCN):

• Les MCN est un procédé itératif, qui nécessite donc un point de départ et des critères d'arrêt. Les MCL sont directs (algèbre linéaire);
• Les MCN nécessitent de calculer la matrice jacobienne (dérivées premières). Une expression analytique peut être compliquée à obtenir: si c'est le cas, une différentiation numérique s'impose;
• La divergence est un problème courant des MCN: en effet, il n'est pas rare de voir augmenter la fonction objectif (somme des carrés des résidus) d'une itération à l'autre. Cela peut être dû au manque de précision de l'approximation linéaire par le développement de Taylor;
• Pour les MCL, la solution est unique mais pas pour les MCN: plusieurs minima (locaux) peuvent exister.

Interprétation géométrique

Dans le cas des moindres carrés linéaires, la fonction objectif S est une fonction quadratique des paramètres

$S=\sum_i W_{ii} \left(y_i-\sum_jX_{ij}\beta_j \right)^2$

Lorsqu'il y a un seul paramètre à estimer β, la fonction S est une parabole en β. Pour deux paramètres ou plus, le contour de S est constitué d'ellipses concentriques, à condition que la matrice $\mathbf{X^TWX}$ est définie positive. Le minimum, atteint pour la valeur optimale des paramètres, est le centre de ces ellipses concentriques.

Dans le cas non linéaire, le contour en ellipses concentriques n'est vrai qu'au voisinage du minimum, puisque dans ce cas l'approximation linéaire de Taylor s'avère être une bonne approximation de la fonction objectif.

$S \approx\sum_i W_{ii} \left(y_i-\sum_j J_{ij}\beta_j \right)^2$

Plus les paramètres s'éloignent de leur valeur optimale, plus le contour dévie de sa forme ellipsoïdale. Ceci signifie qu'il est essentiel de choisir l'approximation initiale $\boldsymbol{\beta}^0$ du procédé itératif proche des valeurs optimales, qui sont par définition inconnues.

Astuces calculatoires

Choix du vecteur d'incrément

Dans le procédé itératif

$\boldsymbol{\beta}^{k+1} = \boldsymbol{\beta}^k +\Delta \boldsymbol{\beta}.$

il est impératif de se prémunir contre la divergence. Plusieurs astuces concernant le choix de $\Delta \boldsymbol{\beta}$ peuvent être considérées:

• changer sa norme sans pour autant changer sa direction. On peut alors introduire une proportion f (entre 0 et 1) et amender le procédé itératif en $\boldsymbol\beta^{k+1}=\boldsymbol\beta^k+f\Delta \boldsymbol\beta$. Par exemple, on peut diviser par deux la valeur de f jusqu'à observer une réelle diminution de la fonction objectif. On peut aussi rechercher la valeur de f par Recherche linéaire^[2]: pour une grille de valeurs de f, on cherche la valeur de f provoquant la diminution la plus importante de S. Mais une telle méthode est couteuse en calculs, car il faut à chaque fois re-calculer la fonction objectif;
• si la direction de l'incrément est trop éloigné de sa direction "optimale" et que la méthode précédente échoue, il faudra peut-être changer légèrement la direction du vecteur d'incrément $\Delta \boldsymbol{\beta}$. Pour cela, les équations normales sont transformées en $\mathbf{\left( J^TWJ +\lambda I \right)\Delta \boldsymbol \beta=\left( J^TW \right) \Delta y}$, où λ est le paramètre de Marquardt^[3] et I la matrice identité. On consultera avec profit la page consacrée à la méthode de Levenberg-Marquardt.

Décompositions

QR

Le minimum de la somme des carrés des résidus peut se trouver sans former les équations normales. Les résidus du modèle linéarisé s'écrivent

$\mathbf{r=\Delta y-J\Delta\boldsymbol\beta}$

La jacobienne fait l'objet d'une décomposition orthogonale, comme par exemple la Décomposition QR:

$\mathbf{J=QR}$

où Q est une matrice orthogonale $m \times m$ et où R est une matrice $m \times n$, partitionnée en un bloc $\mathbf\R_n$, de dimension $n \times n$, et en un bloc nul, de dimension $m-n \times n$ zero block. De plus, $\mathbf\R_n$ est triangulaire supérieure.

$\mathbf{R}= \begin{bmatrix} \mathbf{R}_n \\ \mathbf{0}\end{bmatrix}$

Le vecteur de résidu est pré-multiplié par $\mathbf Q^T$.

$\mathbf{Q^Tr=Q^T \Delta y -R \Delta\boldsymbol\beta}= \begin{bmatrix} \mathbf{\left(Q^T \Delta y -R \Delta\boldsymbol\beta \right)}_n \\ \mathbf{\left(Q^T \Delta y \right)}_{m-n}\end{bmatrix}$

La somme des carrés reste inchangée puisque $S=\mathbf{r^T Q Q^Tr = r^Tr}$. Le minimum de S est atteint lorsque le bloc supérieur est nul. Par conséquent, le vecteur d'incrément recherché $\Delta\boldsymbol\beta$ se trouve en résolvant

$\mathbf{R_n \Delta\boldsymbol\beta =\left(Q^T \Delta y \right)_n}$

La résolution est facile d'accès car la matrice R a été prise triangulaire supérieure.

Décomposition en valeurs singulières

Une variante de la méthode précédente fait intervenir la Décomposition en valeurs singulières, dans laquelle R est diagonalisée:

$\mathbf{J=U \boldsymbol\Sigma V^T}$

où $\mathbf U$ est orthogonal, $\boldsymbol\Sigma$ est une matrice diagonale de valeurs singulières et $\mathbf V^T$ est la matrice orthogonale des vecteurs propres de $\mathbf\ {R_n}$. Dans cette configuration, l'incrément est donné par

$\mathbf{\boldsymbol\Delta\beta=V \boldsymbol\Sigma^{-1}\left( U^T \boldsymbol\Delta y \right)}_n$

La relative simplicité de cette expression est très utile dans l'analyse théorique des moindres carrés. Cette méthode est largement détaillée dans Lawson et Hanson^[4].

Critères de convergence

Le critère de convergence le plus immédiat est que la somme des carrés ne doit pas décroître d'une itération à l'autre. Ce critère est souvent difficile à implémenter. Aussi, on préfère le critère de convergence suivant

$\left|\frac{S^k-S^{k+1}}{S^k}\right|<0,0001$

Évidemment, la valeur 0,0001 est arbitraire et doit être changée. En particulier, il faudra augmenter cette valeur lorsque les erreurs expérimentales sont importantes. Un critère alternatif est

$\left|\frac{\Delta \beta_j}{\beta_j}\right|<0,001, \, j=1,n$

Là encore, la valeur numérique est arbitraire.

Minimums multiples

On peut trouver plusieurs minimums de S dans certaines circonstances:

• Un paramètre intervenant à une certaine puissance. Par exemple, pour ajuster des données à une courbe de Lorentz:

$f(x_i, \boldsymbol \beta)=\frac{\alpha}{1+\left(\frac{\gamma-x_i}{\beta} \right)^2}$

Il y a alors deux solution pour le paramètre β: $\hat \beta$ ainsi que $-\hat \beta$ donnent la même valeur optimale de la fonction critère;

• Deux paramètres peuvent s'interchanger sans changer le critère. Par exemple, cela survient lorsqu'on rencontre le produit de deux paramètres: ainsi αβ donnera la même chose que βα;
• Un paramètre intervient dans une fonction périodique, comme par exemple dans $\sin \beta\,$. Dans ce cas, $\hat \beta +2n \pi$ donne la même valeur critère.

Tous les minimums multiples ne donnent pas la même valeur de la fonction objectif. Un autre problème concerne les minimums locaux. Par exemple, dans le modèle

$f(x_i, \beta)=\left(1-3\beta +\beta^3 \right)x_i$

il y a un minimum local en $\beta\,$ = 1 et un minimum global en $\hat \beta\,$ = -3^[5].

Pour être certain d'avoir obtenu un minimum global, il est souhaitable de recommencer la procédure de minimisation en changeant le point de départ. Quand on obtient le même résultat quel que soit le point de départ, on peut alors penser obtenir un minimum global.

L'existence de minimums multiples a une conséquence importante: la fonction objectif admet une valeur maximum quelque part entre les deux minimums. Les équations normales en ce maximum fait intervenir des matrices non définies positives. Une telle situation est à proscrire, en particulier comme initialisation du procédé itératif. Par exemple, pour le cas de l'ajustement du modèle de Lorentz, le cas β = 0 est à éviter.

Autres méthodes

Linéarisation

Un modèle non linéaire peut parfois se transformer en modèle linéaire. Par exemple, lorsque le modèle est une simple fonction exponentielle,

$f(x_i,\boldsymbol \beta)= \alpha e^{\beta x_i}$

on peut obtenir un modèle linéaire par transformation logarithmique.

$\log f(x_i,\boldsymbol \beta)=\log \alpha + \beta x_i$

La somme des carrés devient

$S=\sum_i (\log y_i-\log \alpha - \beta x_i)^2\!$

Toutefois, si on n'a aucun renseignement sur la structure des aléas, cette transformation peut être problématique: de toute évidence, les erreurs expérimentales sur y ne sont pas les mêmes que sur log y. Estimer le modèle initial et celui linéarisé donnera des estimations différentes et des variances estimées. En pratique, le modèle exponentiel s'estime dans une procédure à part.

Quelques ajustements

On peut considérer quelques ajustements

• Calcul de la matrice jacobienne par approximation numérique. Dans certains modèles, obtenir des dérivées analytiques peut s'avérer délicat. On doit avoir recours à une approximation numérique de la jacobienne; par exemple, l'approximation de son entrée (i,j) est:

$\mathbf{J}_{i,j} = \frac{\partial f(x_i, \boldsymbol \beta)}{\partial \beta_j} \approx \frac{\delta f(x_i, \boldsymbol \beta)}{\delta \beta_j}$

Cette approximation s'obtient par le calcul de $f(x_i, \boldsymbol \beta)\,$ pour $\beta_j\,$ et $\beta_j+\delta \beta_j\,$. L'incrément,$\delta \beta_j\,$, doit être choisi ni trop grand ni trop petit (afin d'éviter les erreurs d'arrondi);

• L'inclusion des dérivées secondes dans le développement de Taylor. On obtient la méthode classique de Newton.

$f(x_i, \boldsymbol \beta)=f^k(x_i, \boldsymbol \beta) +\sum_j J_{ij} \Delta \beta_j + \frac{1}{2}\sum_j\sum_k \Delta\beta_j \Delta\beta_k H_{jk_{(i)}},\ H_{jk_{(i)}}=\frac{\partial^2 f(x_i, \boldsymbol \beta)}{\partial \beta_j \partial \beta_k }$

La matrice H est la Matrice hessienne. Bien que présentant de bien meilleures propriétés de convergence près de l'optimum, ce modèle se comporte mal quand les paramètres sont loin de leur valeur optimale. Le calcul de la Hessienne ajoute à la complexité de l'algorithme. Cette méthode n'est pas utilisée en général;

• On peut remplacer l'algorithme de Newton par un algorithme de pseudo-Newton, où on calcule numériquement la hessienne par approximations successives. C'est l'algorithme de Davidon-Fletcher-Powell (DFP);

Voir aussi

• Optimisation non linéaire
• Optimisation (mathématiques)
• Algorithme de Levenberg-Marquardt

Références

• C. T. Kelley, Iterative Methods for Optimization, SIAM Frontiers in Applied Mathematics, no 18, 1999, ISBN: 0-89871-433-8. Online copy

1. ↑ Ceci implique que les observations ne soient pas corrélées. Dans le cas contraire, on utilise plutôt la quantité

   $S=\sum_k \sum_j r_k W_{kj} r_j\,$

   Dans ce cas, la matrice idéale de pondération devrait être égale à l'inverse de la Matrice de variance-covariance des observations.
2. ↑ M.J. Box, D. Davies and W.H. Swann, Non-Linear optimisation Techniques, Oliver & Boyd, 1969
3. ↑ Cette technique a été proposée indépendamment par Levenberg (1944), Girard (1958), Wynne (1959), Morrison (1960) et Marquardt (1963). C'est le nom de ce dernier que la littérature scientifique a généralement conservé.
4. ↑ C.L. Lawson and R.J. Hanson, Solving Least Squares Problems, Prentice-Hall,1974
5. ↑ P. Gans, Data Fitting in the Chemical Sciences, Wiley, 1992, p71.