Moindres carrés non linéaires

Les moindres carrés non linéaires est une forme des moindres carrés adaptée pour 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.

Considérons un jeu de m couples d'observations, (x₁, y₁), (x₂, y₂),...,(x_m, y_m), et une fonction de régression du type y = f (x, β). Cette fonction dépend des explicatives x mais aussi du vecteur des n paramètres β = (β₁, β₂, ..., β_n) avec m ≥ n. On souhaite trouver le vecteur de paramètres β qui ajuste au mieux les données, au sens des moindres carrés :

${\displaystyle S=\sum _{i=1}^{m}r_{i}^{2}}$

est minimisée en β, où les résidus r_i sont donnés par

${\displaystyle r_{i}=y_{i}-f(x_{i},{\boldsymbol {\beta }}),\ 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 :

${\displaystyle {\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 ${\displaystyle {\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

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

qui fournit des approximations successives β^k de plus en plus proches de la vraie valeur (inconnue) des paramètres, β₀.

À chaque itération, le modèle initial est linéarisé par un développement de Taylor autour de β^k comme suit :

${\displaystyle 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 _{j}}}\left(\beta _{j}-\beta _{j}^{k}\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 termes du modèle linéarisé, ${\displaystyle {\frac {\partial r_{i}}{\partial \beta _{j}}}=-J_{ij}}$ et les résidus sont donnés par

${\displaystyle 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

${\displaystyle -2\sum _{i=1}^{m}J_{ij}\left(\Delta y_{i}-\sum _{s=1}^{n}J_{is}\Delta \beta _{s}\right)=0}$

ou encore

${\displaystyle \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 à

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

La linéarisation permet donc d'écrire :

${\displaystyle {\boldsymbol {\beta }}^{k+1}={\boldsymbol {\beta }}^{k}+\left(\mathbf {J^{T}J} \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 β^k, et permet donc de trouver la prochaine itération β^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

${\displaystyle 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

${\displaystyle \mathbf {\left(J^{T}WJ\right)\Delta {\boldsymbol {\beta }}=J^{T}W\Delta y} }$

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

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

• les MCN sont un procédé itératif, qui nécessitent 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.

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

${\displaystyle S=\sum _{i}W_{ii}\left(y_{i}-\sum _{j}X_{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 X^TWX soit 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.

${\displaystyle 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 β⁰ du procédé itératif proche des valeurs optimales, qui sont par définition inconnues.

Dans le procédé itératif

${\displaystyle {\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 Δβ peuvent être considérées :

• changer sa norme sans pour autant changer sa direction. On peut alors introduire une proportion p (entre 0 et 1) et amender le procédé itératif en ${\displaystyle {\boldsymbol {\beta }}^{k+1}={\boldsymbol {\beta }}^{k}+p\Delta {\boldsymbol {\beta }}}$. Par exemple, on peut diviser par deux la valeur de p jusqu'à observer une réelle diminution de la fonction objectif. On peut aussi rechercher la valeur de p par recherche linéaire^[2] : pour une grille de valeurs de p, on cherche la valeur de p 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 Δβ. Pour cela, les équations normales sont transformées en ${\displaystyle \mathbf {\left(J^{T}WJ+\lambda I\right)\Delta {\boldsymbol {\beta }}=\left(J^{T}W\right)\Delta y} }$, où ${\displaystyle \lambda }$ 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.

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

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

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

${\displaystyle \mathbf {J} =\mathbf {QR} }$

où Q est une matrice orthogonale m × m et où R est une matrice m × n, partitionnée en un bloc R_n, de dimension n × n, et en un bloc nul, de dimension (m–n) × n. De plus, R_n est triangulaire supérieure.

${\displaystyle \mathbf {R} ={\begin{bmatrix}\mathbf {R} _{n}\\\mathbf {0} \end{bmatrix}}}$

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

${\displaystyle \mathbf {Q^{T}r=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 =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é ${\displaystyle \Delta {\boldsymbol {\beta }}}$ se trouve en résolvant

${\displaystyle \mathbf {R} _{n}\Delta {\boldsymbol {\beta }}=\left(\mathbf {Q} ^{T}\Delta \mathbf {y} \right)_{n}}$

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

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

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

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

${\displaystyle \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].

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

${\displaystyle \left|{\frac {S^{k}-S^{k+1}}{S^{k}}}\right|<\varepsilon }$

La valeur du seuil est arbitraire et doit être adaptée en fonction du problème. En particulier, il faudra augmenter cette valeur lorsque les erreurs expérimentales sont importantes. Un critère alternatif est

${\displaystyle \left|{\frac {\Delta \beta _{j}}{\beta _{j}}}\right|<\varepsilon ,\,j=1,n.}$

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 :

${\displaystyle 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 β : ${\displaystyle {\hat {\beta }}}$ ainsi que ${\displaystyle -{\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 ${\displaystyle \alpha \beta }$ donnera la même chose que ${\displaystyle \beta \alpha }$;
• Un paramètre intervient dans une fonction périodique, comme dans ${\displaystyle \sin \beta \,}$. Dans ce cas, ${\displaystyle {\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

${\displaystyle f(x_{i},\beta )=\left(1-3\beta +\beta ^{3}\right)x_{i}}$

il y a un minimum local en β = 1 et un minimum global en ${\displaystyle {\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. On utilise typiquement des points de départ aléatoires, par exemple déterminés par l'algorithme de Metropolis-Hastings ; c'est la méthode du recuit simulé.

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.

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,

${\displaystyle f(x_{i},{\boldsymbol {\beta }})=\alpha e^{\beta x_{i}}}$

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

${\displaystyle \log f(x_{i},{\boldsymbol {\beta }})=\log \alpha +\beta x_{i}}$

La somme des carrés devient

${\displaystyle 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.

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:

${\displaystyle \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 ${\displaystyle f(x_{i},{\boldsymbol {\beta }})\,}$ pour ${\displaystyle \beta _{j}\,}$ et ${\displaystyle \beta _{j}+\delta \beta _{j}\,}$. L'incrément, ${\displaystyle \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.

${\displaystyle 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);

La bibliothèque logicielle en C++ Ceres Solver a été spécifiquement conçue pour résoudre le problème des moindres carrés non linéaires^[6]. Elle est utilisée dans divers produits Google comme Street View^[7].

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

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

• Portail des probabilités et de la statistique