Voir la source de Moindres carrés non linéaires

Les '''moindres carrés non linéaires''' est une forme des [[Méthode des moindres carrés|moindres carrés]] adaptée pour l'estimation d'un modèle non linéaire en ''n'' paramètres à partir de ''m'' observations  (''m''&nbsp;>&nbsp;''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 ==

{{article détaillé|Méthode des moindres carrés}}
Considérons un jeu de {{mvar|m}} couples d'observations, {{math|(''x''{{ind|1}}, ''y''{{ind|1}}), (''x''{{ind|2}}, ''y''{{ind|2}}),...,(''x{{ind|m}}'', ''y{{ind|m}}'')}}, et une fonction de régression du type {{math|''y'' {{=}} ''f'' (''x'', '''''β''''')}}. Cette fonction dépend des ''explicatives'' {{mvar|x}} mais aussi du vecteur des ''n'' paramètres {{math|'''''β''''' {{=}} (''β''{{ind|1}}, ''β''{{ind|2}}, ..., ''β{{ind|n}}'')}} avec {{math|''m'' ≥ ''n''}}. On souhaite trouver le vecteur de paramètres {{mvar|'''β'''}} qui ajuste au mieux les données, au sens des moindres carrés :
:<math>S=\sum_{i=1}^{m}r_i^2</math>
est minimisée en {{mvar|'''β'''}}, où les résidus {{mvar|r<sub>i</sub>}} sont donnés par
:<math>r_i= y_i - f(x_i, \boldsymbol \beta), \ i=1, 2,\dots, m.</math>

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'' :
:<math>\frac{\partial S}{\partial \beta_j}=2\sum_i r_i\frac{\partial r_i}{\partial \beta_j}=0 \ (j=1,\ldots,n).</math>

Dans un système non linéaire, les dérivées <math>\frac{\partial r_i}{\partial \beta_j}</math> 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
:<math>\boldsymbol{\beta}^{k+1} = \boldsymbol{\beta}^k + \Delta \boldsymbol{\beta}.</math>
qui fournit des approximations successives {{mvar|'''β'''{{exp|k}}}} de plus en plus proches de la vraie valeur (inconnue) des paramètres, {{math|'''''β'''''{{ind|0}}}}.

À chaque [[itération]], le modèle initial est linéarisé par un [[Série de Taylor|développement de Taylor]] autour de {{mvar|'''β'''{{exp|k}}}} comme suit :
:<math>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^{k}_j \right) \approx f(x_i,\boldsymbol \beta^k) +\sum_j J_{ij} \Delta\beta_j. </math>
La [[matrice jacobienne]] {{math|'''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é, <math>\frac{\partial r_i}{\partial \beta_j}=-J_{ij}</math> et les résidus sont donnés par
:<math>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).</math>
Les équations normales deviennent
:<math>-2\sum_{i=1}^{m}J_{ij} \left( \Delta y_i-\sum_{s=1}^{n} J_{is}\Delta \beta_s \right)=0</math>
ou encore
:<math>\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).\,</math>
Matriciellement, on en arrive à
:<math>\mathbf{\left(J^TJ\right)\Delta \boldsymbol \beta=J^T\Delta} y.</math>
La linéarisation permet donc d'écrire :
:<math>\boldsymbol{\beta}^{k+1} = \boldsymbol{\beta}^k + \left(\mathbf{J^TJ}\right)^{-1} \mathbf{J^T\Delta} y.</math>
Il faut remarquer que l'ensemble du terme de droite dépend seulement de l'itération en cours, à savoir {{mvar|'''β'''{{exp|k}}}}, et permet donc de trouver la prochaine itération {{math|'''β'''{{exp|''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
:<math>S=\sum_{i=1}^{m}W_{ii}r_i^2.</math>
Idéalement, chaque élément de la [[matrice diagonale]] de pondération {{math|'''W'''}} devrait être égal à l'inverse de la  [[Variance (statistiques et probabilités)|variance]] de l'observation <ref>Ceci implique que les observations ne soient pas corrélées. Dans le cas contraire, on utilise plutôt la quantité
:<math>S=\sum_k \sum_j r_k W_{kj} r_j\,</math>
Dans ce cas, la matrice idéale de pondération devrait être égale à l'inverse de la [[matrice de variance-covariance]] des observations.</ref>  
Les équations normales deviennent alors
:<math>\mathbf{\left(J^TWJ\right)\Delta \boldsymbol \beta=J^TW\Delta y}</math>
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 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.

=== 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
:<math>S=\sum_i W_{ii} \left(y_i-\sum_jX_{ij}\beta_j \right)^2</math>
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'[[Ellipse (mathématiques)|ellipses]] concentriques, à condition que la matrice {{math|'''X'''{{exp|''T''}}'''WX'''}} 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. 
:<math>S \approx\sum_i W_{ii} \left(y_i-\sum_j J_{ij}\beta_j \right)^2</math>
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 {{math|'''''β'''''{{exp|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
:<math>\boldsymbol{\beta}^{k+1} = \boldsymbol{\beta}^k +\Delta \boldsymbol{\beta}.</math>
il est impératif de se prémunir contre la divergence. Plusieurs astuces concernant le choix de {{math|Δ'''''β'''''}} 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 <math>\boldsymbol\beta^{k+1}=\boldsymbol\beta^k+p\Delta \boldsymbol\beta</math>. 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 (optimisation)|recherche linéaire]]<ref name=BDS>M.J. Box, D. Davies and W.H. Swann, Non-Linear optimisation Techniques, Oliver & Boyd, 1969</ref> : 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 {{math|Δ'''''β'''''}}. Pour cela, les équations normales sont transformées en <math>\mathbf{\left( J^TWJ +\lambda I \right)\Delta \boldsymbol \beta=\left( J^TW \right) \Delta y}</math>, où <math>\lambda</math> est le paramètre de Marquardt<ref>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é.</ref> et '''I''' la [[matrice identité]]. On consultera avec profit la page consacrée à la méthode de [[Algorithme de Levenberg-Marquardt|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
:<math>\mathbf{r}=\Delta \mathbf{y}-\mathbf{J}\Delta\boldsymbol\beta</math>
La jacobienne fait l'objet d'une décomposition orthogonale, comme la [[décomposition QR]] :
:<math>\mathbf{J}=\mathbf{QR}</math>
où {{math|'''Q'''}} est une matrice orthogonale {{math|''m'' × ''m''}} et où {{math|'''R'''}} est une matrice {{math|''m'' × ''n''}}, partitionnée en un bloc {{math|'''R'''{{ind|''n''}}}}, de dimension {{math|''n'' × ''n''}}, et en un bloc nul, de dimension {{math|(''m''–''n'') × ''n''}}. De plus, {{math|'''R'''{{ind|''n''}}}} est triangulaire supérieure.
:<math>\mathbf{R}= \begin{bmatrix}
\mathbf{R}_n \\
\mathbf{0}\end{bmatrix}</math>

Le vecteur de résidu est pré-multiplié par {{math|'''Q'''{{exp|''T''}}}}. 
:<math>\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}</math>

La somme des carrés reste inchangée puisque {{math|''S'' {{=}}'''r'''{{exp|''T''}} '''Q Q'''{{exp|''T''}}'''r''' {{=}} '''r'''{{exp|''T''}}'''r'''}}. Le minimum de ''S'' est atteint lorsque le bloc supérieur est nul. Par conséquent, le vecteur d'incrément recherché <math>\Delta\boldsymbol\beta</math> se trouve en résolvant
:<math>\mathbf{R}_n \Delta\boldsymbol\beta =\left(\mathbf{Q}^T \Delta \mathbf{y} \right)_n</math>
La résolution est facile d'accès car la matrice {{math|'''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 {{math|'''R'''}} est diagonalisée :
:<math>\mathbf{J=U \boldsymbol\Sigma V^T}</math>
où {{math|'''U'''}} est orthogonal, {{math|'''Σ'''}} est une matrice diagonale de valeurs singulières et {{math|'''V'''{{exp|''T''}}}} est la [[matrice orthogonale]] des vecteurs propres de {{math|'''R'''{{ind|''n''}}}}. Dans cette configuration, l'incrément est donné par
:<math>\mathbf{\boldsymbol\Delta\beta=V \boldsymbol\Sigma^{-1}\left( U^T  \boldsymbol\Delta y \right)}_n</math>
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<ref name=LH>C.L. Lawson and R.J. Hanson, Solving Least Squares Problems, Prentice-Hall, 1974</ref>.

=== 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
:<math>\left|\frac{S^k-S^{k+1}}{S^k}\right|<\varepsilon</math>
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
:<math>\left|\frac{\Delta \beta_j}{\beta_j}\right|<\varepsilon, \, j=1,n.</math>

== 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 :
:<math>f(x_i, \boldsymbol \beta)=\frac{\alpha}{1+\left(\frac{\gamma-x_i}{\beta} \right)^2}</math>
Il y a alors deux solution pour le paramètre β : <math>\hat \beta</math> ainsi que <math>-\hat \beta</math> 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 <math>\alpha \beta</math> donnera la même chose que <math>\beta \alpha</math>;
* Un paramètre intervient dans une [[fonction périodique]], comme dans <math>\sin \beta\,</math>. Dans ce cas, <math>\hat \beta +2n \pi</math> 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
:<math>f(x_i, \beta)=\left(1-3\beta +\beta^3 \right)x_i</math> 
il y a un minimum local en {{math|''β'' {{=}} 1}} et un minimum global en <math>\hat \beta\,</math> = -3<ref>P. Gans, Data Fitting in the Chemical Sciences, Wiley, 1992, {{p.|71}}.</ref>.

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 {{math|''β'' {{=}} 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]],
:<math>f(x_i,\boldsymbol \beta)= \alpha e^{\beta x_i}</math>
on peut obtenir un modèle linéaire par transformation logarithmique.
:<math>\log f(x_i,\boldsymbol \beta)=\log \alpha + \beta x_i</math>
La somme des carrés devient
:<math>S=\sum_i (\log y_i-\log \alpha - \beta x_i)^2\!</math>

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:
:<math>\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}</math>
Cette approximation s'obtient par le calcul de <math>f(x_i, \boldsymbol \beta)\,</math> pour <math>\beta_j\,</math> et <math>\beta_j+\delta \beta_j\,</math>. L'incrément, <math>\delta \beta_j\,</math>, doit être choisi ni trop grand ni trop petit (afin d'éviter les [[Erreur d'arrondi|erreurs d'arrondi]]);
* L'inclusion des dérivées secondes dans le développement de Taylor. On obtient la méthode classique de Newton.
:<math>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 } </math>
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);


== Implémentations de résolution ==
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<ref>Sameer Agarwal, Keir Mierle et al., [http://ceres-solver.org Ceres Solver].</ref>. Elle est utilisée dans divers produits [[Google]] comme Street View<ref>{{Lien web
 | titre = Users
 | url = http://ceres-solver.org/users.html
 | site = Ceres Solver
 | consulté le = 15 octobre 2016
}}.</ref>.

== Voir aussi ==

=== Articles connexes ===
* [[Optimisation non linéaire]]
* [[Optimisation (mathématiques)]]
* [[Algorithme de Levenberg-Marquardt]]

=== Bibliographie ===
* C. T. Kelley, ''Iterative Methods for Optimization'', SIAM Frontiers in Applied Mathematics, no 18, 1999, {{ISBN|0-89871-433-8}}. [http://www.siam.org/books/textbooks/fr18_book.pdf  Online copy]

=== Notes et références ===
<references/>

{{Portail|probabilités et statistiques}}

[[Catégorie:Algorithme d'optimisation]]
[[Catégorie:Régression]]