Gauss-Newton

Algorithme de Gauss-Newton

L'algorithme de Gauss-Newton est une méthode de résolution des problèmes de moindres carrés non-linéaires. Elle peut être vue comme une modification de la méthode de Newton dans le cas multidimensionnel afin de trouver le minimum d'une fonction (à plusieurs variables). Mais l'algorithme de Gauss-Newton est totalement spécifique à la minimisation d'une somme de fonctions au carré et présente le grand avantage de ne pas nécessiter les dérivées secondes, parfois complexes à calculer.

Les problèmes de moindres carrés non-linéaires surviennent par exemple dans les problèmes de régressions non-linéaires, où des paramètres du modèle sont recherchés afin de correspondre au mieux aux observations disponibles.

Cette méthode est due au mathématicien renommé Carl Friedrich Gauss.

Sommaire

1 Algorithme
2 Remarques
3 Exemple
4 Propriété de convergence
5 Dérivation à partir de la méthode de Newton
6 Versions améliorées
7 Algorithmes associés
8 Notes et références
9 Voir aussi
- 9.1 Articles connexes

Algorithme

Soit m fonctions $r i$ ( $i=1,\ldots,m$ ) de n variables ${\boldsymbol \beta}=(\beta_1, \beta_2, \dots, \beta_n),$ avec m≥n, l'algorithme de Gauss–Newton doit trouver le minimum de la somme des carrés

$S(\boldsymbol \beta)= \sum_{i=1}^m r_i^2(\boldsymbol \beta).$ ^[1]

En débutant avec l'approximation initiale $\boldsymbol \beta^0$ du minimum, la méthode procède par itérations:

$\boldsymbol \beta^{s+1} = \boldsymbol \beta^s+\delta\boldsymbol\beta,$

où l'incrément $\delta\boldsymbol\beta$ vérifie les équations normales

$\left(\mathbf{J_r^T J_r} \right)\delta\boldsymbol\beta= - \mathbf{ J_r^T r}.$

Ici, on note par r le vecteur des fonctions r_i, et par J_r la matrice jacobienne m×n de r par rapport à β, tous les deux évalués en β^s. La matrice transposée est notée à l'aide de l'exposant T.

Dans les problèmes d'ajustement des données, où le but est de trouver les paramètres $\boldsymbol \beta$ d'un certain modèle $y=f(x, \boldsymbol \beta)$ permettant le meilleur ajustement aux observations $(x i, y i),$ , les fonctions r_i sont les résidus

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

Alors, l'incrément $\delta\boldsymbol\beta$ peut s'exprimer en fonction de la jacobienne de la fonction f:

$\left( \mathbf{ J_f^T J_f} \right)\delta\boldsymbol\beta= {\mathbf{J_f^T r}}.$

Dans tous les cas, une fois connu l'estimation à l'étape s, les équations normales permettent de trouver l'estimation à l'étape suivante; pour résumer, on a:

$\boldsymbol \beta^{s+1} = \boldsymbol \beta^s - \left(\mathbf{J_r^T J_r} \right)^{-1} \mathbf{ J_r^T r}$

L'ensemble du terme de droite est calculable car ne dépend que de $\boldsymbol\beta^s$ et permet de trouver l'estimation suivante.

Remarques

L'hypothèse m≥n est nécessaire, car dans le cas contraire la matrice $\mathbf{J_r^T J_r}$ serait non inversible et les équations normales ne pourraient être résolues.

L'algorithme de Gauss–Newton peut être dérivé par approximation linéaire du vecteur de fonctions $r i .$ En utilisant le Théorème de Taylor, on peut écrire qu'à chaque itération

$\mathbf{r}(\boldsymbol \beta_0)\approx \mathbf{r}(\boldsymbol \beta^s)+\mathbf{J_r}(\boldsymbol \beta^s)\delta\boldsymbol\beta$

avec $\delta\boldsymbol \beta=\boldsymbol \beta_0-\boldsymbol \beta^s$ ; notons que $\boldsymbol\beta_0$ représente la vraie valeur des paramètres pour laquelle les résidus $\mathbf{r}(\boldsymbol \beta_0)$ s'annulent. Trouver l'incrément $\delta\boldsymbol \beta$ revient à résoudre

$-\mathbf{r}(\boldsymbol \beta^s) \approx \mathbf{J_r}(\boldsymbol \beta^s)\delta\boldsymbol\beta$

ce qui peut se faire par la technique classique de régression linéaire et qui fournit exactement les équations normales.

Les équations normales sont un système de m équations linéaires d'inconnu $\delta \boldsymbol\beta$ . Ce système peut se résoudre en une étape, en utilisant la factorisation de Cholesky ou, encore mieux, la décomposition QR de J_r. Pour de grands systèmes, une méthode itérative telle que la méthode du gradient conjugué peut être plus efficace. S'il existe une dépendance linéaire entre les colonnes J_r, la méthode échouera car $\mathbf{J_r^T J_r}$ deviendra singulier.

Exemple

Courbe calculée avec $\hat\beta_1=0.362$ et $\hat\beta_2=0.556$ (en bleu) contre les données observées (en rouge).

Dans cet exemple, l'algorithme de Gauss–Newton est utilisé pour ajuster un modèle en minimisant la somme des carrés entre les observations et les prévisions du modèle.

Dans une expérience de biologie, on étudie la relation entre la concentration du substrat [S] et la vitesse de réaction dans une réaction enzymatique à partir de données reportées dans le tableau suivant.

i	1	2	3	4	5	6	7
[S]	0.038	0.194	0.425	0.626	1.253	2.500	3.740
rate	0.050	0.127	0.094	0.2122	0.2729	0.2665	0.3317

On souhaite ajuster les données à la courbe de la forme:

$\text{rate}=\frac{V_{max}[S]}{K_M+[S]}$

L'estimation par moindres carrés porte sur les paramètres $V m a x$ et $K M$ .

On note $x i$ et $y i$ les valeurs de $[S]$ et la vitesse de réaction, pour $i=1, \dots, 7.$ On pose $β 1 = V m a x$ et $β 2 = K M .$ Nous allons chercher $β 1$ et $β 2$ pour minimiser la somme des carrés des résidus

$r_i = y_i - \frac{\beta_1x_i}{\beta_2+x_i}$ ( $i=1,\dots, 7$ )

La jacobienne $\mathbf{J_r}$ du vecteur des résidus $r i$ par rapport aux inconnus $β j$ est une matrice $7\times 2$ dont la ligne n° i est

$\frac{\partial r_i}{\partial \beta_1}= -\frac{x_i}{\beta_2+x_i},\ \frac{\partial r_i}{\partial \beta_2}= \frac{\beta_1x_i}{\left(\beta_2+x_i\right)^2}.$

Commençant avec l'estimation initiale $β 1$ =0,9 et $β 2$ =0,2, il suffit de 5 itérations de l'algorithme de Gauss–Newton pour atteindre les estimations optimales $\hat\beta_1=0,362$ et $\hat\beta_2=0,556$ . Le critère de la somme des carrés des résidus chute de 1,202 à 0,0886 en 5 itérations. Le tableau suivant détaille les cinq itérations:

Itération	Estimation	Somme des carrés des résidus
1	[0,9;0,2]	1,4455000
2	[0,33266;0,26017]	0,0150721
3	[0,34281;0,42608]	0,0084583
4	[0,35778;0,52951]	0,0078643
5	[0,36141;0,55366]	0,0078442
6	[ 0,3618;0,55607]	0,0078440

La figure ci-contre permet de juger de l'adéquation du modèle aux données en comparant la courbe ajustée (bleue) aux observations (rouge).

Propriété de convergence

On peut démontrer que l'incrément $δβ$ est une direction de descente pour S ^[2], et que si l'algorithme converge, alors la limite est un point stationnaire pour la somme des carrés S. Toutefois, la convergence n'est pas garantie, pas plus qu'une convergence locale contrairement à la méthode de Newton.

La vitesse de convergence de l'algorithme de Gauss–Newton peut approcher la vitesse quadratique.^[3] L'algorithme peut converger lentement voire ne pas converger du tout si le point de départ de l'algorithme est trop loin du minimum ou si la matrice $\mathbf{J_r^T J_r}$ est mal conditionnée.

L'algorithme peut donc échouer à converger. Par exemple, le problème avec $m = 2$ équations et $n = 1$ variable, donné par

$\begin{align} r_1(\beta) &= \beta + 1 \\ r_2(\beta) &= \lambda \beta^2 + \beta - 1. \end{align}$

L'optimum se situe en $β = 0$ . Si $λ = 0$ alors le problème est en fait linéaire et la méthode trouve la solution en une seule itération. Si |λ| < 1, alors la méthode converge linéairement et les erreurs décroissent avec un facteur |λ| à chaque itération. Cependant, si |λ| > 1, alors la méthode ne converge même pas localement.^[4]

Dérivation à partir de la méthode de Newton

Dans ce qui suit, l'algorithme de Gauss–Newton sera tiré de l'algorithme d'optimisation de Newton; par conséquent, la vitesse de convergence sera au plus quadratique.

La relation de récurrence de la méthode de Newton pour minimiser une fonction S de paramètres $\boldsymbol\beta$ , est

$\boldsymbol\beta^{s+1} = \boldsymbol\beta^s - \mathbf H^{-1} \mathbf g \,$

où g représente le gradient de S et H sa matrice hessienne. Puisque $S = \sum_{i=1}^m r_i^2$ , le gradient est

$g_j=2\sum_{i=1}^m r_i\frac{\partial r_i}{\partial \beta_j}.$

Les éléments de la Hessienne sont calculés en dérivant les éléments du gradient, $g j$ , par rapport à $β k$

$H_{jk}=2\sum_{i=1}^m \left(\frac{\partial r_i}{\partial \beta_j}\frac{\partial r_i}{\partial \beta_k}+r_i\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k} \right).$

La méthode de Gauss–Newton est obtenue en ignorant les dérivées d'ordre supérieur à deux. La Hessienne est approchée par

$H_{jk}\approx 2\sum_{i=1}^m J_{ij}J_{ik}$

où $J_{ij}=\frac{\partial r_i}{\partial \beta_j}$ est l'élément $(i, j)$ de la jacobienne $\mathbf{J_r}$ . Le gradient et la hessienne approchée sont alors

$\mathbf g=2\mathbf{J_r^Tr, \quad H \approx 2J_r^TJ_r}.\,$

Ces expressions sont replacées dans la relation de récurrence initiale afin d'obtenir la relation récursive

$\boldsymbol{\beta}^{s+1} = \boldsymbol\beta^s+\delta\boldsymbol\beta;\ \delta\boldsymbol\beta= - \mathbf{\left( J_r^T J_r \right)^{-1} J_r^T r}.$

La convergence de la méthode n'est pas toujours garantie. L'approximation

$\left|r_i\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k}\right| \ll \left|\frac{\partial r_i}{\partial \beta_j}\frac{\partial r_i}{\partial \beta_k}\right|$

doit être vraie pour pouvoir ignorer les dérivées du second ordre. Cette approximation peut être valide dans deux cas, pour lesquels on peut s'attendre à obtenir la convergence:^[5]

Les valeurs de la fonction $r i$ sont petites en magnitude, au moins près du minimum;
Les fonctions sont seulement faiblement non-linéaires, si bien que $\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k}$ est relativement petit en magnitude.

Versions améliorées

Avec la méthode de Gauss–Newton, la somme des carrés S peut ne pas décroître à chaque itération. Toutefois, puisque $\delta\boldsymbol\beta$ est une direction de descente, à moins que $S(\boldsymbol \beta^s)$ soit un point stationnaire, il se trouve que $S(\boldsymbol \beta^s+\alpha\ \delta\boldsymbol\beta) < S(\boldsymbol \beta^s)$ pour tout $α > 0$ suffisamment petit. Ainsi, en cas de divergence, une solution est d'employer une fraction, $α$ , de l'incrément, $\delta\boldsymbol\beta$ dans la formule de mise à jour

$\boldsymbol \beta^{s+1} = \boldsymbol \beta^s+\alpha\ \delta\boldsymbol\beta$ .

En d'autres termes, le vecteur d'incrément est trop long, mais pointe bien vers le bas, si bien que parcourir une fraction du chemin fait décroître la fonction objectif S. Une valeur optimale pour $α$ peut être trouvée en utilisant un algorithme de Recherche linéaire: la magnitude de $α$ est déterminée en trouvant le minimum de S en faisant varier α sur une grille de l'intervalle $0 < α < 1$ .

Dans le cas où la direction de l'incrément est telle que α est proche de zéro, une méthode alternative pour éviter la divergence est l'algorithme de Levenberg-Marquardt. Les équations normales sont modifiées de telle sorte que l'incrément est décalé en direction de la descente la plus forte

$\left(\mathbf{J^TJ+\lambda D}\right)\delta\boldsymbol\beta=\mathbf{J}^T \mathbf{r}$ ,

où D est une matrice diagonale positive. Remarquons que lorsque D est la matrice identité et que $\lambda\to+\infty$ , alors $\delta\boldsymbol\beta/\lambda\to \mathbf{J}^T \mathbf{r}$ , par conséquent la direction de $\delta\boldsymbol\beta$ s'approche de la direction du gradient $\mathbf{J}^T \mathbf{r}$ .

Le paramètre de Marquardt, $λ$ , peut aussi être optimisé par une méthode de recherche linéaire, mais ceci rend la méthode fort inefficace dans la mesure où le vecteur d'incrément doit être re-calculé à chaque fois que $λ$ change.

Algorithmes associés

Dans une méthode quasi-Newton, comme celle due à Davidon, Fletcher et Powell, une estimation de la matrice Hessienne, $\frac{\partial^2 S}{\partial \beta_j \partial\beta_k}$ , est approchée numériquement en utilisant les premières dérivées $\frac{\partial r_i}{\partial\beta_j}$ .

Une autre méthode pour résoudre les problèmes de moindres carrés en utilisant seulement les dérivées premières est la descente de gradient. Toutefois, cette méthode ne prend pas en compte les dérivées secondes, même sommairement. Par conséquent, cette méthode s'avère particulièrement inefficace pour beaucoup de fonctions.

Notes et références

↑ (en) A. Björck, Numerical methods for least squares problems, SIAM, Philadelphia (ISBN 0-89871-360-9)
↑ Björck, p260
↑ Björck, p341, 342
↑ « ».
↑ (en) Jorge Nocedal, Numerical optimization, New York: Springer (ISBN 0387987932)

Voir aussi

Regardez d'autres dictionnaires:

Gauss–Newton algorithm — The Gauss–Newton algorithm is a method used to solve non linear least squares problems. It can be seen as a modification of Newton s method for finding a minimum of a function. Unlike Newton s method, the Gauss–Newton algorithm can only be used… … Wikipedia
Gauss-Newton-Verfahren — Das Gauß Newton Verfahren (nach Carl Friedrich Gauß und Isaac Newton) ist ein numerisches Verfahren zur Lösung nichtlinearer Minimierungsprobleme, die durch Anwendung der Methode der kleinsten Quadrate auf nichtlineare Ausgleichsprobleme… … Deutsch Wikipedia
Algorithme de Gauss-Newton — En mathématiques, l algorithme de Gauss Newton est une méthode de résolution des problèmes de moindres carrés non linéaires. Elle peut être vue comme une modification de la méthode de Newton dans le cas multidimensionnel afin de trouver le… … Wikipédia en Français
Algoritmo de Gauss-Newton — En matemáticas, el algoritmo de Gauss Newton se utiliza para resolver problemas no lineales de mínimos cuadrados. Es una modificación del método de optimización de Newton que no usa segundas derivadas y se debe a Carl Friedrich Gauss. El problema … Wikipedia Español
Generalized Gauss–Newton method — The generalized Gauss–Newton method is a generalization of the least squares method originally described by Carl Friedrich Gauss and of Newton s method due to Isaac Newton to the case of constrained nonlinear least squares problems … Wikipedia
Newton's method in optimization — A comparison of gradient descent (green) and Newton s method (red) for minimizing a function (with small step sizes). Newton s method uses curvature information to take a more direct route. In mathematics, Newton s method is an iterative method… … Wikipedia
Gauss (homonymie) — Cette page d’homonymie répertorie les différents sujets et articles partageant un même nom. Carl Friedrich Gauss (1777 1855), mathématicien, astronome et physicien allemand. Le gauss, une unité de mesure du champ magnétique, noté G. GAUSS, un… … Wikipédia en Français
Newton's method — In numerical analysis, Newton s method (also known as the Newton–Raphson method), named after Isaac Newton and Joseph Raphson, is a method for finding successively better approximations to the roots (or zeroes) of a real valued function. The… … Wikipedia
Gauss' law for gravity — In physics, Gauss law for gravity, also known as Gauss flux theorem for gravity, is a law of physics which is essentially equivalent to Newton s law of universal gravitation. Its form is mathematically similar to Gauss law for electricity; in… … Wikipedia
Newton polynomial — In the mathematical field of numerical analysis, a Newton polynomial, named after its inventor Isaac Newton, is the interpolation polynomial for a given set of data points in the Newton form. The Newton polynomial is sometimes called Newton s… … Wikipedia

Dictionnaires et Encyclopédies sur 'Academic'

Gauss-Newton

Algorithme de Gauss-Newton

Sommaire

Algorithme

Remarques

Exemple

Propriété de convergence

Dérivation à partir de la méthode de Newton

Versions améliorées

Algorithmes associés

Notes et références

Voir aussi

Articles connexes

Regardez d'autres dictionnaires:

Share the article and excerpts

Dictionnaires et Encyclopédies sur 'Academic'

Wikipédia en Français

Gauss-Newton

Algorithme de Gauss-Newton

Sommaire

Algorithme

Remarques

Exemple

Propriété de convergence

Dérivation à partir de la méthode de Newton

Versions améliorées

Algorithmes associés

Notes et références

Voir aussi

Articles connexes

Regardez d'autres dictionnaires:

Share the article and excerpts

Direct link