Algoritmo de Gauss-Newton

En matemáticas , el algoritmo de Gauss-Newton es un método para resolver problemas de mínimos cuadrados no lineales . Puede verse como una modificación del método de Newton en el caso multidimensional para encontrar el mínimo de una función (con varias variables). Pero el algoritmo de Gauss-Newton es totalmente específico para la minimización de una suma de funciones cuadradas y tiene la gran ventaja de no requerir las segundas derivadas , que a veces son complejas de calcular.

Los problemas de mínimos cuadrados no lineales surgen, por ejemplo, en los problemas de regresión no lineal, donde se buscan los parámetros del modelo para que coincidan mejor con las observaciones disponibles.

Este método se debe a Carl Friedrich Gauss .

Algoritmo

Sean m funciones ( ) de n variables con m ≥ n , el algoritmo de Gauss-Newton debe encontrar el mínimo de la suma de cuadrados: $Rhode Island$ $i = 1, \ ldots, m$ ${{\ boldsymbol \ beta}} = (\ beta _ {1}, \ beta _ {2}, \ dots, \ beta _ {n}),$

S ({\ boldsymbol \ beta}) = \ sum _ {{i = 1}} ^ {m} r_ {i} ^ {2} ({\ boldsymbol \ beta}).

Suponiendo un valor inicial del mínimo, el método procede por iteraciones : ${\ boldsymbol \ beta} ^ {0}$

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

donde el incremento satisface las ecuaciones normales $\ delta {\ boldsymbol \ beta}$

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

Aquí, denotamos por r el vector de funciones r i , y por J r la matriz jacobiana m × n de r con respecto a β , ambos evaluados en β s . La matriz transpuesta se denota con el exponente T .

En los problemas de ajuste de datos, donde el objetivo es encontrar los parámetros de un determinado modelo que permitan el mejor ajuste a las observaciones , las funciones r i son los residuos ${\ boldsymbol \ beta}$ $y = f (x, {\ boldsymbol \ beta})$ $(x_ {i}, y_ {i}),$

r_ {i} ({\ boldsymbol \ beta}) = y_ {i} -f (x_ {i}, {\ boldsymbol \ beta}).

Entonces, el incremento se puede expresar en función del jacobiano de la función f : $\ delta {\ boldsymbol \ beta}$

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

En todos los casos, una vez que se conoce la estimación en el paso s , las ecuaciones normales permiten encontrar la estimación en el siguiente paso; para resumir, tenemos:

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

El término completo de la derecha es computable porque depende solo y permite encontrar la siguiente estimación. ${\ boldsymbol \ beta} ^ {s}$

Observaciones

El supuesto m ≥ n es necesario, porque de lo contrario la matriz sería no invertible y las ecuaciones normales no podrían resolverse. ${\ mathbf {J_ {r} ^ {T} J_ {r}}}$

El algoritmo de Gauss-Newton se puede derivar por aproximación lineal del vector de funciones. Usando el teorema de Taylor , podemos escribir eso en cada iteración $r_ {i}.$

{\ mathbf {r}} ({\ boldsymbol \ beta} _ {0}) \ approx {\ mathbf {r}} ({\ boldsymbol \ beta} ^ {s}) + {\ mathbf {J_ {r}} } ({\ boldsymbol \ beta} ^ {s}) \ delta {\ boldsymbol \ beta}

con ; tenga en cuenta que representa el valor real de los parámetros para los cuales los residuos se cancelan entre sí. Encontrar el incremento es como resolver $\ delta {\ boldsymbol \ beta} = {\ boldsymbol \ beta} _ {0} - {\ boldsymbol \ beta} ^ {s}$ ${\ boldsymbol \ beta} _ {0}$ ${\ mathbf {r}} ({\ boldsymbol \ beta} _ {0})$ $\ delta {\ boldsymbol \ beta}$

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

que se puede hacer mediante la técnica clásica de regresión lineal y que proporciona exactamente las ecuaciones normales.

Las ecuaciones normales son un sistema de m ecuaciones lineales de desconocido . Este sistema se puede resolver en un solo paso, utilizando la factorización de Cholesky o, mejor aún, la descomposición QR de J r . Para sistemas grandes, un método iterativo como el método de gradiente conjugado puede ser más eficiente. Si hay una dependencia lineal entre las columnas J r , el método fallará porque se volverá singular. $\ delta {\ boldsymbol \ beta}$ ${\ mathbf {J_ {r} ^ {T} J_ {r}}}$

Finalmente, observemos que el método de Gauss-Newton es efectivo cuando el error cuadrático final es bajo, o cuando la no linealidad es “no muy pronunciada”. El método es particularmente sensible a la presencia de puntos "aberrantes" (es decir, ubicados lejos de la curva del modelo).

Ejemplo

En este ejemplo, el algoritmo de Gauss-Newton se utiliza para ajustar un modelo minimizando la suma de cuadrados entre las observaciones y las predicciones del modelo.

En un experimento de biología, estudiamos la relación entre la concentración del sustrato [ S ] y la velocidad de reacción en una reacción enzimática a partir de los datos presentados en la siguiente tabla.

I	1	2	3	4	5	6	7
[S]	0,038	0,194	0,425	0,626	1,253	2500	3.740
omitido	0,050	0,127	0.094	0.2122	0.2729	0.2665	0.3317

Queremos ajustar los datos a la curva de la forma:

{\ text {rate}} = {\ frac {V _ {{\ max}} [S]} {K_ {M} + [S]}}

La estimación de mínimos cuadrados cubre los parámetros y . $V _ {{\ max}}$ $K_ {M}$

Denotamos y los valores de y la velocidad de reacción, por lo que establecemos y vamos a encontrar y minimizar la suma de los cuadrados de los residuales $x_ {i}$ $y_i$ $[S]$ $i = 1, \ puntos, 7.$ $\ beta _ {1} = V _ {{\ max}}$ $\ beta _ {2} = K_ {M}.$ $\ beta_1$ $\ beta _ {2}$

r_ {i} = y_ {i} - {\ frac {\ beta _ {1} x_ {i}} {\ beta _ {2} + x_ {i}}}

( )

i = 1, \ puntos, 7

El jacobiano del vector de los residuales con respecto a las incógnitas es una matriz cuya línea n ° i es ${\ mathbf {J_ {r}}}$ $Rhode Island$ $\ beta _ {j}$ $7 \ por 2$

{\ frac {\ parcial r_ {i}} {\ parcial \ beta _ {1}}} = - {\ frac {x_ {i}} {\ beta _ {2} + x_ {i}}}, \ { \ frac {\ parcial r_ {i}} {\ parcial \ beta _ {2}}} = {\ frac {\ beta _ {1} x_ {i}} {\ left (\ beta _ {2} + x_ { i} \ derecha) ^ {2}}}.

Comenzando con la estimación inicial = 0,9 y = 0,2, se necesitan 5 iteraciones del algoritmo de Gauss-Newton para alcanzar las estimaciones óptimas y . El criterio de la suma de cuadrados de los residuos cae de 1,202 a 0,0886 en 5 iteraciones. La siguiente tabla detalla las cinco iteraciones: $\ beta_1$ $\ beta _ {2}$ ${\ hat \ beta} _ {1} = 0.362$ ${\ hat \ beta} _ {2} = 0.556$

Iteración	Estimar	Suma de cuadrados de residuos
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 figura opuesta permite juzgar la idoneidad del modelo a los datos comparando la curva ajustada (azul) con las observaciones (rojo).

En un segundo ejemplo, tenemos una medida de espectrometría o difractometría que exhibe un pico asimétrico, con ruido. Queremos ajustar un modelo de pico asimétrico compuesto por dos medios gaussianos que tienen la misma posición (expectativa) y la misma altura, pero diferentes "anchos" (desviaciones estándar): la función del modelo es de la forma

{\ Displaystyle \ left \ {{\ begin {alineado} f (x) = \ & \ mathrm {A} (2) \ times \ exp \ left ({\ frac {(x- \ mathrm {A} (1)) ) ^ {2}} {\ mathrm {A} (3)}} \ right) {\ text {si}} x \ leqslant \ mathrm {A} (1) {\ text {; }} \\ f (x) = \ & \ mathrm {A} (2) \ times \ exp \ left ({\ frac {(x- \ mathrm {A} (1)) ^ {2}} {\ mathrm {A} (4)}} \ right) {\ text {si}} x \ geqslant \ mathrm {A} (1) {\ text {.}} \ End {alineado}} \ right.}

La regresión consiste en ajustar los parámetros A (1), A (2), A (3) y A (4).

Primero, determinamos numéricamente la derivada y la segunda derivada para obtener los parámetros iniciales del pico, con el algoritmo de Savitzky-Golay :

cima de la curva estimada por el mínimo de la segunda derivada, A (1) = -0.03;
valor de la curva suavizada en esta ubicación: A (2) = 9,7;
medio ancho a media altura a la izquierda igual ah g = (0,84 - 0,03) (desde la posición del punto de inflexión);
medio ancho a media altura a la derecha igual ah d = (0,45 + 0,03).

Para un gaussiano puro, la mitad del ancho a la mitad de la altura h está relacionada con la desviación estándar σ por:

2 h = 2,35 × σ

y tenemos

A (3, 4) = 2σ 2 ,

sea A (3, 4) = 8 ( h g, d / 2.35) 2 .

El algoritmo converge en 5 pasos (variación de la desviación cuadrática normalizada menor que 10 −7 ) con los siguientes resultados:

A (1) = 0,00404, para un valor teórico de 0;
A (2) = 9,83, "10;
A (3) = 1,02, "1;
A (4) = 0,313, "0,3.

En la figura opuesta, los puntos experimentales (200 en número) forman la curva azul, el modelo ajustado está representado por la curva roja.

Si partimos del conjunto de parámetros arbitrarios iniciales [1, 1, 1, 1], el algoritmo converge en diez pasos con la condición de utilizar un factor de amortiguación α ajustado automáticamente en cada paso (ver más abajo).

Propiedad de convergencia

Se puede demostrar que el incremento es una dirección de descenso de S , y si los converge el algoritmo, entonces el límite es un punto estacionario para la suma de cuadrados S . Sin embargo, la convergencia no está garantizada, no más que una convergencia local contraria al método de Newton. $\ delta \ beta$

La velocidad de convergencia del algoritmo de Gauss-Newton puede aproximarse a la velocidad cuadrática. El algoritmo puede converger lentamente o incluso no converger en absoluto si el punto de partida del algoritmo está demasiado lejos del mínimo o si la matriz está mal acondicionada . ${\ mathbf {J_ {r} ^ {T} J_ {r}}}$

Por tanto, es posible que el algoritmo no converja. Por ejemplo, el problema con ecuaciones y variables, dado por $m = 2$ $n = 1$

{\ begin {alineado} r_ {1} (\ beta) & = \ beta +1 \\ r_ {2} (\ beta) & = \ lambda \ beta ^ {2} + \ beta -1. \ end {alineado }}

Lo óptimo está en . Si entonces el problema es lineal y el método encuentra la solución en una sola iteración. Si | λ | <1, entonces el método converge linealmente y los errores disminuyen con un factor | λ | en cada iteración. Sin embargo, si | λ | > 1, entonces el método ni siquiera converge localmente. $\ beta = 0$ $\ lambda = 0$

Derivación del método de Newton

En lo que sigue, el algoritmo de Gauss-Newton se tomará del algoritmo de optimización de Newton; por lo tanto, la velocidad de convergencia será como máximo cuadrática.

La relación de recurrencia del método de Newton para minimizar una función S de parámetros es ${\ boldsymbol \ beta}$

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

donde g representa el gradiente de S y H su matriz de Hesse . Dado que , el gradiente es $S = \ sum _ {{i = 1}} ^ {m} r_ {i} ^ {2}$

g_ {j} = 2 \ sum _ {{i = 1}} ^ {m} r_ {i} {\ frac {\ parcial r_ {i}} {\ parcial \ beta _ {j}}}.

Los elementos del hessiano se calculan derivando los elementos del gradiente, con respecto a $g_ {j}$ $\ beta _ {k}$

H _ {{jk}} = 2 \ sum _ {{i = 1}} ^ {m} \ izquierda ({\ frac {\ parcial r_ {i}} {\ parcial \ beta _ {j}}} {\ frac {\ parcial r_ {i}} {\ parcial \ beta _ {k}}} + r_ {i} {\ frac {\ parcial ^ {2} r_ {i}} {\ parcial \ beta _ {j} \ parcial \ beta _ {k}}} \ derecha).

El método de Gauss-Newton se obtiene ignorando las derivadas de orden mayor que dos. El Hessian es abordado por

H _ {{jk}} \ approx 2 \ sum _ {{i = 1}} ^ {m} J _ {{ij}} J _ {{ik}}

donde está el elemento del jacobiano . El gradiente y el hessiano aproximado son entonces $J _ {{ij}} = {\ frac {\ parcial r_ {i}} {\ parcial \ beta _ {j}}}$ $(yo, j)$ ${\ mathbf {J_ {r}}}$

{\ mathbf g} = 2 {\ mathbf {J_ {r} ^ {T} r, \ quad H \ approx 2J_ {r} ^ {T} J_ {r}}}. \,

Estas expresiones se reemplazan en la relación de recurrencia inicial para obtener la relación recursiva

{\ 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 convergencia del método no siempre está garantizada. La aproximación

\ izquierda | r_ {i} {\ frac {\ parcial ^ {2} r_ {i}} {\ parcial \ beta _ {j} \ parcial \ beta _ {k}}} \ derecha | \ ll \ izquierda | { \ frac {\ parcial r_ {i}} {\ parcial \ beta _ {j}}} {\ frac {\ parcial r_ {i}} {\ parcial \ beta _ {k}}} \ derecha |

debe ser verdadera para ignorar las derivadas de segundo orden. Esta aproximación puede ser válida en dos casos, para los cuales se puede esperar obtener convergencia:

los valores de la función son pequeños en magnitud, al menos cerca del mínimo; $Rhode Island$
las funciones son solo débilmente no lineales, por lo que su magnitud es relativamente pequeña. ${\ frac {\ parcial ^ {2} r_ {i}} {\ parcial \ beta _ {j} \ parcial \ beta _ {k}}}$

Versiones mejoradas

Con el método de Gauss-Newton, es posible que la suma de cuadrados S no disminuya con cada iteración. Sin embargo, dado que es una dirección de descenso, a menos que sea un punto estacionario, resulta que para todo lo suficientemente pequeño, $\ delta {\ boldsymbol \ beta}$ $S ({\ boldsymbol \ beta} ^ {s})$ $\ alpha> 0$

S ({\ boldsymbol \ beta} ^ {s} + \ alpha \ \ delta {\ boldsymbol \ beta}) <S ({\ boldsymbol \ beta} ^ {s})

Por tanto, en caso de discrepancia, una solución es utilizar un factor de amortiguación del incremento en la fórmula de actualización: $\alfa$ $\ delta {\ boldsymbol \ beta}$

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

En otras palabras, el vector de la subasta es demasiado largo, pero está apuntando hacia abajo, por lo que navegar por una fracción de la forma en que disminuye la función objetivo S . Se puede encontrar un valor óptimo para usando un algoritmo de búsqueda lineal : la magnitud de se determina encontrando el mínimo de S variando α en una cuadrícula del intervalo . También podemos reevaluar α de forma sencilla en cada paso, por ejemplo mediante dicotomía : $\alfa$ $\alfa$ $0 <\ alpha <1$

reduciéndolo (por ejemplo, aplicando un factor de 0,5) hasta que tengamos S (β s + αδβ s ) <S (β s );
incrementándolo la próxima vez que se cumpla la condición (por ejemplo, acercándolo a 1 tomando (1 + α) / 2).

Este enfoque no es óptimo, pero reduce el esfuerzo requerido para determinar α.

En el caso de que la dirección del incremento sea tal que α esté cerca de cero, un método alternativo para evitar la divergencia es el algoritmo de Levenberg-Marquardt . Las ecuaciones normales se cambian para que el incremento se cambie hacia el descenso más empinado.

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

donde D es una matriz diagonal positiva. Observe que cuando D es la matriz identidad y que , entonces , la dirección de se acerca a la dirección del gradiente . $\ lambda \ a + \ infty$ $\ delta {\ boldsymbol \ beta} / \ lambda \ to {\ mathbf {J}} ^ {T} {\ mathbf {r}}$ $\ delta {\ boldsymbol \ beta}$ ${\ mathbf {J}} ^ {T} {\ mathbf {r}}$

El parámetro de Marquardt ,, también se puede optimizar mediante un método de búsqueda lineal, pero esto hace que el método sea muy ineficaz, ya que el vector de incremento debe recalcularse cada vez que cambia. $\ lambda$ $\ lambda$

Algoritmos asociados

En un método cuasi-Newton, como el de Davidon, Fletcher y Powell, una estimación de la matriz hessiana, se aproxima numéricamente usando las primeras derivadas . ${\ frac {\ parcial ^ {2} S} {\ parcial \ beta _ {j} \ parcial \ beta _ {k}}}$ ${\ frac {\ parcial r_ {i}} {\ parcial \ beta _ {j}}}$

Otro método para resolver problemas de mínimos cuadrados utilizando solo las primeras derivadas es el algoritmo de gradiente . Sin embargo, este método no tiene en cuenta las segundas derivadas, ni siquiera de forma resumida. Por lo tanto, este método es particularmente ineficaz para muchas funciones.

Notas y referencias

(fr) Este artículo está tomado parcial o totalmente del artículo de Wikipedia en inglés titulado “ Algoritmo de Gauss-Newton ” ( ver la lista de autores ) .

(en) Åke Björck , Métodos numéricos para problemas de mínimos cuadrados , Filadelfia, SIAM ,1996, 408 p. ( ISBN 978-0-89871-360-2 , LCCN 96003908 , presentación en línea )
R. Fletcher , Métodos prácticos de optimización , Nueva York, John iley and Sons,1987
Björck, pág. 260
Björck, pág. 341-342
(in) Roger Fletcher , Métodos prácticos de optimización , Nueva York, 2 ,1987, 436 p. ( ISBN 978-0-471-91547-8 ) , pág. 113
(in) Jorge Nocedal y Stephen J. Wright , Optimización numérica , Nueva York, Springer ,1999, 2 nd ed. , 636 p. ( ISBN 978-0-387-98793-4 , presentación en línea )

enlaces externos

Calculando incertidumbres Un libro de 246 páginas que estudia la regresión lineal y no lineal. El algoritmo se detalla y se aplica al experimento de biología tratado como ejemplo en este artículo (página 79 con las incertidumbres sobre los valores estimados).