Algoritmo de Gauss-Newton

Algoritmo matemático

Ajuste de una curva ruidosa mediante un modelo de pico asimétrico con parámetros minimizando la suma de residuos cuadrados en los puntos de la cuadrícula , utilizando el algoritmo de Gauss-Newton. Arriba: datos sin procesar y modelo. Abajo: Evolución de la suma normalizada de los cuadrados de los errores. ${\displaystyle f_{\beta }(x)}$ ${\displaystyle \beta }$ ${\displaystyle r_{i}(\beta )=y_{i}-f_{\beta }(x_{i})}$ ${\displaystyle x_{i}}$

El algoritmo de Gauss-Newton se utiliza para resolver problemas de mínimos cuadrados no lineales , lo que equivale a minimizar una suma de valores de función al cuadrado. Es una extensión del método de Newton para encontrar el mínimo de una función no lineal . Dado que una suma de cuadrados no debe ser negativa, se puede considerar que el algoritmo utiliza el método de Newton para aproximar iterativamente los ceros de los componentes de la suma y, por lo tanto, minimizar la suma. En este sentido, el algoritmo también es un método eficaz para resolver sistemas de ecuaciones sobredeterminados. Tiene la ventaja de que no se requieren segundas derivadas, que pueden resultar difíciles de calcular. ^[1]

Los problemas de mínimos cuadrados no lineales surgen, por ejemplo, en la regresión no lineal , donde se buscan parámetros en un modelo de modo que el modelo concuerde bien con las observaciones disponibles.

El método lleva el nombre de los matemáticos Carl Friedrich Gauss e Isaac Newton , y apareció por primera vez en la obra de Gauss de 1809 Theoria motus corporum coelestium in sectionibus conicis solem ambientum . ^[2]

Descripción

Funciones dadas (a menudo llamadas residuales) de variables con el algoritmo de Gauss-Newton encuentran iterativamente el valor de que minimiza la suma de cuadrados ^[3] ${\displaystyle m}$ ${\displaystyle {\textbf {r}}=(r_{1},\ldots ,r_{m})}$ ${\displaystyle n}$ ${\displaystyle {\boldsymbol {\beta }}=(\beta _{1},\ldots \beta _{n}),}$ ${\displaystyle m\geq n,}$ ${\displaystyle \beta }$ $${\displaystyle S({\boldsymbol {\beta }})=\sum _{i=1}^{m}r_{i}({\boldsymbol {\beta }})^{2}.}$$

Comenzando con una estimación inicial del mínimo, el método avanza mediante iteraciones ${\displaystyle {\boldsymbol {\beta }}^{(0)}}$ $${\displaystyle {\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}-\left(\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {J_{r}} \right)^{-1}\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right),}$$

donde, si r y β son vectores columna , las entradas de la matriz jacobiana son $${\displaystyle \left(\mathbf {J_{r}} \right)_{ij}={\frac {\partial r_{i}\left({\boldsymbol {\beta }}^{(s)}\right)}{\partial \beta _{j}}},}$$

y el símbolo denota la transposición de la matriz . ${\displaystyle ^{\operatorname {T} }}$

En cada iteración, la actualización se puede encontrar reorganizando la ecuación anterior en los dos pasos siguientes: ${\displaystyle \Delta ={\boldsymbol {\beta }}^{(s+1)}-{\boldsymbol {\beta }}^{(s)}}$

• ${\displaystyle \Delta =-\left(\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {J_{r}} \right)^{-1}\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)}$
• ${\displaystyle \mathbf {J_{r}} ^{\operatorname {T} }\mathbf {J_{r}} \Delta =-\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)}$

Con sustituciones , y , esto se convierte en la ecuación matricial convencional de la forma , que luego se puede resolver mediante una variedad de métodos (ver Notas). ${\textstyle A=\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {J_{r}} }$ ${\displaystyle \mathbf {b} =-\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)}$ ${\displaystyle \mathbf {x} =\Delta }$ ${\displaystyle A\mathbf {x} =\mathbf {b} }$

Si m = n , la iteración se simplifica a

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

que es una generalización directa del método de Newton en una dimensión.

En el ajuste de datos, donde el objetivo es encontrar los parámetros tales que una función de modelo determinada se ajuste mejor a algunos puntos de datos , las funciones son los residuales : ${\displaystyle {\boldsymbol {\beta }}}$ ${\displaystyle \mathbf {f} (\mathbf {x} ,{\boldsymbol {\beta }})}$ ${\displaystyle (x_{i},y_{i})}$ ${\displaystyle r_{i}}$ $${\displaystyle r_{i}({\boldsymbol {\beta }})=y_{i}-f\left(x_{i},{\boldsymbol {\beta }}\right).}$$

Entonces, el método de Gauss-Newton se puede expresar en términos del jacobiano de la función como ${\displaystyle \mathbf {J_{f}} =-\mathbf {J_{r}} }$ ${\displaystyle \mathbf {f} }$ $${\displaystyle {\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}+\left(\mathbf {J_{f}} ^{\operatorname {T} }\mathbf {J_{f}} \right)^{-1}\mathbf {J_{f}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right).}$$

Tenga en cuenta que es la pseudoinversa izquierda de . ${\displaystyle \left(\mathbf {J_{f}} ^{\operatorname {T} }\mathbf {J_{f}} \right)^{-1}\mathbf {J_{f}} ^{\operatorname {T} }}$ ${\displaystyle \mathbf {J_{f}} }$

Notas

La suposición m ≥ n en el enunciado del algoritmo es necesaria, ya que de lo contrario la matriz no es invertible y las ecuaciones normales no se pueden resolver (al menos de forma única). ${\displaystyle \mathbf {J_{r}} ^{T}\mathbf {J_{r}} }$

El algoritmo de Gauss-Newton se puede derivar aproximando linealmente el vector de funciones r _i . Usando el teorema de Taylor , podemos escribir en cada iteración: $${\displaystyle \mathbf {r} ({\boldsymbol {\beta }})\approx \mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)+\mathbf {J_{r}} \left({\boldsymbol {\beta }}^{(s)}\right)\Delta }$$

con . La tarea de encontrar minimizando la suma de cuadrados del lado derecho; es decir, ${\displaystyle \Delta ={\boldsymbol {\beta }}-{\boldsymbol {\beta }}^{(s)}}$ ${\displaystyle \Delta }$ $${\displaystyle \min \left\|\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)+\mathbf {J_{r}} \left({\boldsymbol {\beta }}^{(s)}\right)\Delta \right\|_{2}^{2},}$$

es un problema lineal de mínimos cuadrados , que se puede resolver explícitamente, produciendo las ecuaciones normales en el algoritmo.

Las ecuaciones normales son n ecuaciones lineales simultáneas en incrementos desconocidos . Se pueden resolver en un solo paso, usando la descomposición de Cholesky o, mejor, la factorización QR de . Para sistemas grandes, un método iterativo , como el método del gradiente conjugado , puede ser más eficiente. Si hay una dependencia lineal entre columnas de J _r , las iteraciones fallarán, ya que se vuelve singular. ${\displaystyle \Delta }$ ${\displaystyle \mathbf {J_{r}} }$ ${\displaystyle \mathbf {J_{r}} ^{T}\mathbf {J_{r}} }$

Cuando es complejo se debe utilizar la forma conjugada: . ${\displaystyle \mathbf {r} }$ ${\displaystyle \mathbf {r} :\mathbb {C} ^{n}\to \mathbb {C} }$ ${\displaystyle \left({\overline {\mathbf {J_{r}} }}^{\operatorname {T} }\mathbf {J_{r}} \right)^{-1}{\overline {\mathbf {J_{r}} }}^{\operatorname {T} }}$

Ejemplo

Curva calculada obtenida con y (en azul) versus los datos observados (en rojo) ${\displaystyle {\hat {\beta }}_{1}=0.362}$ ${\displaystyle {\hat {\beta }}_{2}=0.556}$

En este ejemplo, se utilizará el algoritmo de Gauss-Newton para ajustar un modelo a algunos datos minimizando la suma de cuadrados de los errores entre los datos y las predicciones del modelo.

En un experimento de biología que estudia la relación entre la concentración de sustrato [ S ] y la velocidad de reacción en una reacción mediada por enzimas, se obtuvieron los datos de la siguiente tabla.

Se desea encontrar una curva (función modelo) de la forma $${\displaystyle {\text{rate}}={\frac {V_{\text{max}}\cdot [S]}{K_{M}+[S]}}}$$

que mejor se ajuste a los datos en el sentido de mínimos cuadrados, con los parámetros y por determinar. ${\displaystyle V_{\text{max}}}$ ${\displaystyle K_{M}}$

Denota por y los valores de [ S ] y tasa respectivamente, con . Deja y . Encontraremos y tal que la suma de los cuadrados de los residuos ${\displaystyle x_{i}}$ ${\displaystyle y_{i}}$ ${\displaystyle i=1,\dots ,7}$ ${\displaystyle \beta _{1}=V_{\text{max}}}$ ${\displaystyle \beta _{2}=K_{M}}$ ${\displaystyle \beta _{1}}$ ${\displaystyle \beta _{2}}$ $${\displaystyle r_{i}=y_{i}-{\frac {\beta _{1}x_{i}}{\beta _{2}+x_{i}}},\quad (i=1,\dots ,7)}$$

se minimiza.

El jacobiano del vector de residuos con respecto a las incógnitas es una matriz cuya fila -ésima tiene las entradas ${\displaystyle \mathbf {J_{r}} }$ ${\displaystyle r_{i}}$ ${\displaystyle \beta _{j}}$ ${\displaystyle 7\times 2}$ ${\displaystyle i}$ $${\displaystyle {\frac {\partial r_{i}}{\partial \beta _{1}}}=-{\frac {x_{i}}{\beta _{2}+x_{i}}};\quad {\frac {\partial r_{i}}{\partial \beta _{2}}}={\frac {\beta _{1}\cdot x_{i}}{\left(\beta _{2}+x_{i}\right)^{2}}}.}$$

A partir de las estimaciones iniciales de y , después de cinco iteraciones del algoritmo de Gauss-Newton, se obtienen los valores óptimos de y . La suma de los cuadrados de los residuos disminuyó del valor inicial de 1,445 a 0,00784 después de la quinta iteración. El gráfico de la figura de la derecha muestra la curva determinada por el modelo para los parámetros óptimos con los datos observados. ${\displaystyle \beta _{1}=0.9}$ ${\displaystyle \beta _{2}=0.2}$ ${\displaystyle {\hat {\beta }}_{1}=0.362}$ ${\displaystyle {\hat {\beta }}_{2}=0.556}$

Propiedades de convergencia

Se garantiza que la iteración de Gauss-Newton convergerá hacia un punto mínimo local bajo 4 condiciones: ^[4] Las funciones son dos veces diferenciables continuamente en un conjunto convexo abierto , el jacobiano tiene rango de columna completo, la iteración inicial es cercana y el local El valor mínimo es pequeño. La convergencia es cuadrática si . ${\displaystyle {\hat {\beta }}}$ ${\displaystyle r_{1},\ldots ,r_{m}}$ ${\displaystyle D\ni {\hat {\beta }}}$ ${\displaystyle \mathbf {J} _{\mathbf {r} }({\hat {\beta }})}$ ${\displaystyle \beta ^{(0)}}$ ${\displaystyle {\hat {\beta }}}$ ${\displaystyle |S({\hat {\beta }})|}$ ${\displaystyle |S({\hat {\beta }})|=0}$

Se puede demostrar ^[5] que el incremento Δ es una dirección de descenso para S y , si el algoritmo converge, entonces el límite es un punto estacionario de S. Sin embargo, para valores mínimos grandes , la convergencia no está garantizada, ni siquiera la convergencia local como en el método de Newton , o la convergencia en las condiciones habituales de Wolfe. ^[6] ${\displaystyle |S({\hat {\beta }})|}$

La tasa de convergencia del algoritmo de Gauss-Newton puede acercarse a la cuadrática . ^[7] El algoritmo puede converger lentamente o no converger en absoluto si la estimación inicial está lejos del mínimo o la matriz está mal condicionada . Por ejemplo, considere el problema con ecuaciones y variables, dado por ${\displaystyle \mathbf {J_{r}^{\operatorname {T} }J_{r}} }$ ${\displaystyle m=2}$ ${\displaystyle n=1}$ $${\displaystyle {\begin{aligned}r_{1}(\beta )&=\beta +1,\\r_{2}(\beta )&=\lambda \beta ^{2}+\beta -1.\end{aligned}}}$$

El óptimo está en . (En realidad, el óptimo es para , porque , pero .) Si , entonces el problema es lineal y el método encuentra el óptimo en una iteración. Si |λ| < 1, entonces el método converge linealmente y el error disminuye asintóticamente con un factor |λ| en cada iteración. Sin embargo, si |λ| > 1, entonces el método ni siquiera converge localmente. ^[8] ${\displaystyle \beta =0}$ ${\displaystyle \beta =-1}$ ${\displaystyle \lambda =2}$ ${\displaystyle S(0)=1^{2}+(-1)^{2}=2}$ ${\displaystyle S(-1)=0}$ ${\displaystyle \lambda =0}$

Resolver sistemas de ecuaciones sobredeterminados.

La iteración de Gauss-Newton es un método eficaz para resolver sistemas de ecuaciones sobredeterminados en forma de con y dónde está la inversa de Moore-Penrose (también conocida como pseudoinversa ) de la matriz jacobiana de . Puede considerarse una extensión del método de Newton y disfruta de la misma convergencia cuadrática local ^[4] hacia soluciones regulares aisladas. $${\displaystyle \mathbf {x} ^{(k+1)}=\mathbf {x} ^{(k)}-J(\mathbf {x} ^{(k)})^{\dagger }\mathbf {f} (\mathbf {x} ^{(k)})\,,\quad k=0,1,\ldots }$$ ${\displaystyle \mathbf {f} (\mathbf {x} )=\mathbf {0} }$ $${\displaystyle \mathbf {f} (\mathbf {x} )={\begin{bmatrix}f_{1}(x_{1},\ldots ,x_{n})\\\vdots \\f_{m}(x_{1},\ldots ,x_{n})\end{bmatrix}}}$$ ${\displaystyle m>n}$ ${\displaystyle J(\mathbf {x} )^{\dagger }}$ ${\displaystyle J(\mathbf {x} )}$ ${\displaystyle \mathbf {f} (\mathbf {x} )}$

Si la solución no existe pero la iteración inicial está cerca de un punto en el que la suma de cuadrados alcanza un mínimo local pequeño, la iteración de Gauss-Newton converge linealmente a . El punto a menudo se denomina solución de mínimos cuadrados del sistema sobredeterminado. ${\displaystyle \mathbf {x} ^{(0)}}$ ${\displaystyle {\hat {\mathbf {x} }}=({\hat {x}}_{1},\ldots ,{\hat {x}}_{n})}$ ${\textstyle \sum _{i=1}^{m}|f_{i}(x_{1},\ldots ,x_{n})|^{2}\equiv \|\mathbf {f} (\mathbf {x} )\|_{2}^{2}}$ ${\displaystyle {\hat {\mathbf {x} }}}$ ${\displaystyle {\hat {\mathbf {x} }}}$

Derivación del método de Newton

A continuación, el algoritmo de Gauss-Newton se derivará del método de Newton para la optimización de funciones mediante una aproximación. Como consecuencia, la tasa de convergencia del algoritmo de Gauss-Newton puede ser cuadrática bajo ciertas condiciones de regularidad. En general (en condiciones más débiles), la tasa de convergencia es lineal. ^[9]

La relación de recurrencia del método de Newton para minimizar una función S de parámetros es ${\displaystyle {\boldsymbol {\beta }}}$ $${\displaystyle {\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}-\mathbf {H} ^{-1}\mathbf {g} ,}$$

donde g denota el vector gradiente de S y H denota la matriz de Hesse de S.

Dado que , el gradiente está dado por ${\textstyle S=\sum _{i=1}^{m}r_{i}^{2}}$ $${\displaystyle g_{j}=2\sum _{i=1}^{m}r_{i}{\frac {\partial r_{i}}{\partial \beta _{j}}}.}$$

Los elementos del Hessiano se calculan diferenciando los elementos del gradiente, con respecto a : ${\displaystyle g_{j}}$ ${\displaystyle \beta _{k}}$ $${\displaystyle 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).}$$

El método de Gauss-Newton se obtiene ignorando los términos derivados de segundo orden (el segundo término de esta expresión). Es decir, el hessiano se aproxima por $${\displaystyle H_{jk}\approx 2\sum _{i=1}^{m}J_{ij}J_{ik},}$$

¿Dónde están las entradas del jacobiano J _r ? Tenga en cuenta que cuando la arpillera exacta se evalúa cerca de un ajuste exacto tenemos cerca de cero , por lo que el segundo término también se vuelve cercano a cero, lo que justifica la aproximación. El gradiente y el hessiano aproximado se pueden escribir en notación matricial como ${\textstyle J_{ij}={\partial r_{i}}/{\partial \beta _{j}}}$ ${\displaystyle r_{i}}$ $${\displaystyle \mathbf {g} =2{\mathbf {J} _{\mathbf {r} }}^{\operatorname {T} }\mathbf {r} ,\quad \mathbf {H} \approx 2{\mathbf {J} _{\mathbf {r} }}^{\operatorname {T} }\mathbf {J_{r}} .}$$

Estas expresiones se sustituyen en la relación de recurrencia anterior para obtener las ecuaciones operativas. $${\displaystyle {\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}+\Delta ;\quad \Delta =-\left(\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {J_{r}} \right)^{-1}\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {r} .}$$

La convergencia del método Gauss-Newton no está garantizada en todos los casos. la aproximación $${\displaystyle \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|}$$

que debe cumplirse para poder ignorar los términos de la derivada de segundo orden puede ser válido en dos casos, para los cuales se espera convergencia: ^[10]

1. Los valores de la función son de magnitud pequeña, al menos alrededor del mínimo. ${\displaystyle r_{i}}$
2. Las funciones son sólo "levemente" no lineales, por lo que su magnitud es relativamente pequeña. ${\textstyle {\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k}}}}$

Versiones mejoradas

Con el método de Gauss-Newton, la suma de los cuadrados de los residuos S puede no disminuir en cada iteración. Sin embargo, dado que Δ es una dirección de descenso, a menos que sea un punto estacionario, se cumple que para todos es suficientemente pequeño . Por tanto, si se produce divergencia, una solución es emplear una fracción del vector de incremento Δ en la fórmula de actualización: ${\displaystyle S\left({\boldsymbol {\beta }}^{s}\right)}$ ${\displaystyle S\left({\boldsymbol {\beta }}^{s}+\alpha \Delta \right)<S\left({\boldsymbol {\beta }}^{s}\right)}$ ${\displaystyle \alpha >0}$ ${\displaystyle \alpha }$ $${\displaystyle {\boldsymbol {\beta }}^{s+1}={\boldsymbol {\beta }}^{s}+\alpha \Delta .}$$

En otras palabras, el vector de incremento es demasiado largo, pero aún apunta "cuesta abajo " , por lo que recorrer solo una parte del camino disminuirá la función objetivo S. Se puede encontrar un valor óptimo de usando un algoritmo de búsqueda de líneas , es decir, la magnitud de se determina encontrando el valor que minimiza S , generalmente usando un método de búsqueda directa en el intervalo o una búsqueda de líneas de retroceso como la búsqueda de líneas de Armijo. . Normalmente, debe elegirse de modo que satisfaga las condiciones de Wolfe o las condiciones de Goldstein . ^[11] ${\displaystyle \alpha }$ ${\displaystyle \alpha }$ ${\displaystyle 0<\alpha <1}$ ${\displaystyle \alpha }$

En los casos en que la dirección del vector de desplazamiento es tal que la fracción óptima α es cercana a cero, un método alternativo para manejar la divergencia es el uso del algoritmo de Levenberg-Marquardt , un método de región de confianza . ^[3] Las ecuaciones normales se modifican de tal manera que el vector de incremento se gira hacia la dirección de descenso más pronunciado , $${\displaystyle \left(\mathbf {J^{\operatorname {T} }J+\lambda D} \right)\Delta =-\mathbf {J} ^{\operatorname {T} }\mathbf {r} ,}$$

donde D es una matriz diagonal positiva. Tenga en cuenta que cuando D es la matriz identidad I y , entonces , la dirección de Δ se aproxima a la dirección del gradiente negativo . ${\displaystyle \lambda \to +\infty }$ ${\displaystyle \lambda \Delta =\lambda \left(\mathbf {J^{\operatorname {T} }J} +\lambda \mathbf {I} \right)^{-1}\left(-\mathbf {J} ^{\operatorname {T} }\mathbf {r} \right)=\left(\mathbf {I} -\mathbf {J^{\operatorname {T} }J} /\lambda +\cdots \right)\left(-\mathbf {J} ^{\operatorname {T} }\mathbf {r} \right)\to -\mathbf {J} ^{\operatorname {T} }\mathbf {r} }$ ${\displaystyle -\mathbf {J} ^{\operatorname {T} }\mathbf {r} }$

El llamado parámetro de Marquardt también se puede optimizar mediante una búsqueda de líneas, pero esto es ineficiente, ya que el vector de desplazamiento debe recalcularse cada vez que se cambia. Una estrategia más eficiente es esta: cuando ocurre divergencia, aumenta el parámetro de Marquardt hasta que haya una disminución en S. Luego conserve el valor de una iteración a la siguiente, pero redúzcalo si es posible hasta alcanzar un valor de corte, cuando el parámetro Marquardt se pueda establecer en cero; la minimización de S se convierte entonces en una minimización estándar de Gauss-Newton. ${\displaystyle \lambda }$ ${\displaystyle \lambda }$

Optimización a gran escala

Para la optimización a gran escala, el método de Gauss-Newton es de especial interés porque a menudo (aunque ciertamente no siempre) es cierto que la matriz es más dispersa que la hessiana aproximada . En tales casos, el cálculo del paso en sí generalmente deberá realizarse con un método iterativo aproximado apropiado para problemas grandes y dispersos, como el método del gradiente conjugado . ${\displaystyle \mathbf {J} _{\mathbf {r} }}$ ${\displaystyle \mathbf {J} _{\mathbf {r} }^{\operatorname {T} }\mathbf {J_{r}} }$

Para que este tipo de enfoque funcione, se necesita al menos un método eficiente para calcular el producto. $${\displaystyle {\mathbf {J} _{\mathbf {r} }}^{\operatorname {T} }\mathbf {J_{r}} \mathbf {p} }$$

para algún vector p . Con un almacenamiento de matriz escaso , en general es práctico almacenar las filas en forma comprimida (por ejemplo, sin entradas cero), lo que dificulta el cálculo directo del producto anterior debido a la transposición. Sin embargo, si se define c _i como la fila i de la matriz , se cumple la siguiente relación simple: ${\displaystyle \mathbf {J} _{\mathbf {r} }}$ ${\displaystyle \mathbf {J} _{\mathbf {r} }}$ $${\displaystyle {\mathbf {J} _{\mathbf {r} }}^{\operatorname {T} }\mathbf {J_{r}} \mathbf {p} =\sum _{i}\mathbf {c} _{i}\left(\mathbf {c} _{i}\cdot \mathbf {p} \right),}$$

de modo que cada fila contribuya de forma aditiva e independiente al producto. Además de respetar una estructura práctica de almacenamiento disperso, esta expresión es muy adecuada para cálculos paralelos . Tenga en cuenta que cada fila ci _es el gradiente del correspondiente r _i residual ; Teniendo esto en cuenta, la fórmula anterior enfatiza el hecho de que los residuos contribuyen al problema independientemente unos de otros.

Algoritmos relacionados

En un método cuasi-Newton , como el de Davidon, Fletcher y Powell o Broyden-Fletcher-Goldfarb-Shanno ( método BFGS ) , se construye numéricamente una estimación del hessiano completo utilizando primeras derivadas solo de modo que después de n ciclos de refinamiento el El método se aproxima mucho al método de Newton en cuanto a rendimiento. Tenga en cuenta que los métodos cuasi-Newton pueden minimizar funciones generales de valores reales, mientras que Gauss-Newton, Levenberg-Marquardt, etc. se ajustan sólo a problemas de mínimos cuadrados no lineales. ${\textstyle {\frac {\partial ^{2}S}{\partial \beta _{j}\partial \beta _{k}}}}$ ${\textstyle {\frac {\partial r_{i}}{\partial \beta _{j}}}}$

Otro método para resolver problemas de minimización utilizando sólo primeras derivadas es el descenso de gradiente . Sin embargo, este método no tiene en cuenta las segundas derivadas ni siquiera de forma aproximada. En consecuencia, es muy ineficiente para muchas funciones, especialmente si los parámetros tienen interacciones fuertes.

Implementaciones de ejemplo

Julia

La siguiente implementación en Julia proporciona un método que utiliza un jacobiano proporcionado y otro cálculo con diferenciación automática .

"""  gaussnewton(r,J,β₀,maxiter,tol)Realice la optimización de Gauss-Newton para minimizar la función residual `r` con `J` jacobiano a partir de `β₀`. El algoritmo termina cuando la norma del paso es menor que "tol" o después de iteraciones "maxiter". """ función gaussnewton ( r , J , β₀ , maxiter , tol ) β = copiar ( β₀ ) para _ en 1 : maxiter Jβ = J ( β ); Δ = - ( Jβ '* Jβ ) \ ( Jβ '* r ( β )) β += Δ si sqrt ( suma ( abs2 , Δ )) < tol romper fin terminar regresar β fin                             importe AbstractDifferentiation como AD , Zygote backend = AD . ZygoteBackend () # otros backends están disponibles       """  gaussnewton(r,β₀,maxiter,tol)Realice la optimización de Gauss-Newton para minimizar la función residual `r` a partir de `β₀`. El jacobiano relevante se calcula mediante diferenciación automática. El algoritmo termina cuando la norma del paso es menor que "tol" o después de iteraciones "maxiter". """ función gaussnewton ( r , β₀ , maxiter , tol ) β = copiar ( β₀ ) para _ en 1 : maxiter rβ , Jβ = AD . value_and_jacobian ( backend , r , β ) Δ = - ( Jβ [ 1 ] '* Jβ [ 1 ]) \ ( Jβ [ 1 ] '* rβ ) β += Δ if sqrt ( suma ( abs2 , Δ )) < tol romper fin fin regresar β fin                              

Notas

1. Mittelhammer, Ron C.; Molinero, Douglas J.; Juez, George G. (2000). Fundamentos econométricos. Cambridge: Prensa de la Universidad de Cambridge. págs. 197-198. ISBN 0-521-62394-4.
2. Floudas, Christodoulos A .; Pardalos, Panos M. (2008). Enciclopedia de optimización . Saltador. pag. 1130.ISBN​ 9780387747583.
3. Björck (1996)
4. JE Dennis, Jr. y RB Schnabel (1983). Métodos numéricos para optimización sin restricciones y ecuaciones no lineales . Reproducción SIAM 1996 de la edición Prentice-Hall 1983. pag. 222.
5. Björck (1996), pág. 260.
6. Mascarenhas (2013), "La divergencia de los métodos BFGS y Gauss Newton", Programación matemática , 147 (1): 253–276, arXiv : 1309.7922 , doi : 10.1007/s10107-013-0720-6, S2CID  14700106
7. Björck (1996), pág. 341, 342.
8. Fletcher (1987), pág. 113.
9. "Copia archivada" (PDF) . Archivado desde el original (PDF) el 4 de agosto de 2016 . Consultado el 25 de abril de 2014 .{{cite web}}: CS1 maint: archived copy as title (link)
10. Nocedal (1999), pág. 259.
11. Nocedal, Jorge. (1999). Optimización numérica . Wright, Stephen J., 1960-. Nueva York: Springer. ISBN 0387227423. OCLC  54849297.

Referencias

• Björck, A. (1996). Métodos numéricos para problemas de mínimos cuadrados . SIAM, Filadelfia. ISBN 0-89871-360-9.
• Fletcher, Roger (1987). Métodos prácticos de optimización (2ª ed.). Nueva York: John Wiley & Sons . ISBN 978-0-471-91547-8..
• Nocedal, Jorge; Wright, Stephen (1999). Optimización numérica . Nueva York: Springer. ISBN 0-387-98793-2.

enlaces externos

• Probabilidad, Estadística y Estimación El algoritmo se detalla y se aplica al experimento de biología discutido como ejemplo en este artículo (página 84 con las incertidumbres sobre los valores estimados).

Implementaciones

• Artelys Knitro es un solucionador no lineal con implementación del método de Gauss-Newton. Está escrito en C y tiene interfaces para C++/C#/Java/Python/MATLAB/R.