Algoritmo de Gauss-Newton

El algoritmo de Gauss-Newton se utiliza para resolver problemas de mínimos cuadrados no lineales , lo que equivale a minimizar la suma de los cuadrados de una función. Es una extensión del método de Newton para encontrar el mínimo de una función no lineal . Dado que la suma de cuadrados debe ser no negativa, el algoritmo puede considerarse como una aplicación del método de Newton para aproximar iterativamente los ceros de las componentes de la suma, minimizando así dicha 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 requiere el cálculo de segundas derivadas, que pueden ser 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 manera 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 ]}

Dado${\displaystyle m}$funciones${\displaystyle {\textbf {r}}=(r_{1},\ldots ,r_{m})}$(a menudo llamados residuos) de${\displaystyle n}$variables${\displaystyle {\boldsymbol {\beta }}=(\beta _{1},\ldots \beta _{n}),}$con${\displaystyle m\geq n,}$El algoritmo de Gauss-Newton encuentra iterativamente el valor de${\displaystyle \beta }$que minimizan la suma de cuadrados ^{[ 3 ]}$${\displaystyle S({\boldsymbol {\beta }})=\sum _{i=1}^{m}r_{i}({\boldsymbol {\beta }})^{2}.}$$

Comenzando con una suposición inicial${\displaystyle {\boldsymbol {\beta }}^{(0)}}$Para el mínimo, el método procede mediante iteraciones $${\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${\displaystyle ^{\operatorname {T} }}$denota la transpuesta de la matriz .

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

• ${\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${\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)}$, y${\displaystyle \mathbf {x} =\Delta }$, esto se convierte en la ecuación matricial convencional de forma${\displaystyle A\mathbf {x} =\mathbf {b} }$, que luego se puede resolver mediante diversos métodos (véase Notas ).

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${\displaystyle {\boldsymbol {\beta }}}$de tal manera que una función modelo dada${\displaystyle \mathbf {f} (\mathbf {x} ,{\boldsymbol {\beta }})}$Se ajusta mejor a algunos puntos de datos${\displaystyle (x_{i},y_{i})}$, las funciones${\displaystyle r_{i}}$son los residuos : $${\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.${\displaystyle \mathbf {J_{f}} =-\mathbf {J_{r}} }$de la función${\displaystyle \mathbf {f} }$como $${\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${\displaystyle \left(\mathbf {J_{f}} ^{\operatorname {T} }\mathbf {J_{f}} \right)^{-1}\mathbf {J_{f}} ^{\operatorname {T} }}$es la pseudoinversa izquierda de${\displaystyle \mathbf {J_{f}} }$.

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

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${\displaystyle \Delta ={\boldsymbol {\beta }}-{\boldsymbol {\beta }}^{(s)}}$La tarea de encontrar${\displaystyle \Delta }$minimizando la suma de los cuadrados del lado derecho; es decir, $${\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 de mínimos cuadrados lineales , que se puede resolver explícitamente, obteniendo así las ecuaciones normales en el algoritmo.

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

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

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

En un experimento biológico que estudiaba la relación entre la concentración del 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 ajusta a los datos en el sentido de mínimos cuadrados, con los parámetros${\displaystyle V_{\text{max}}}$y${\displaystyle K_{M}}$Por determinar.

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

se minimiza.

El jacobino${\displaystyle \mathbf {J_{r}} }$del vector de residuos${\displaystyle r_{i}}$con respecto a las incógnitas${\displaystyle \beta _{j}}$es un${\displaystyle 7\times 2}$matriz con el${\displaystyle i}$-fila que tiene las entradas $${\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}}}.}$$

Comenzando con las estimaciones iniciales de${\displaystyle \beta _{1}=0.9}$y${\displaystyle \beta _{2}=0.2}$Después de cinco iteraciones del algoritmo de Gauss-Newton, los valores óptimos${\displaystyle {\hat {\beta }}_{1}=0.362}$y${\displaystyle {\hat {\beta }}_{2}=0.556}$Se obtienen los siguientes resultados. La suma de los cuadrados de los residuos disminuyó del valor inicial de 1,445 a 0,00784 tras 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.

Se garantiza que la iteración de Gauss-Newton convergerá hacia un punto mínimo local.${\displaystyle {\hat {\beta }}}$bajo 4 condiciones: ^{[ 4 ]} Las funciones${\displaystyle r_{1},\ldots ,r_{m}}$son dos veces continuamente diferenciables en un conjunto convexo abierto${\displaystyle D\ni {\hat {\beta }}}$, el jacobino${\displaystyle \mathbf {J} _{\mathbf {r} }({\hat {\beta }})}$es de rango de columna completo, la iteración inicial${\displaystyle \beta ^{(0)}}$está cerca${\displaystyle {\hat {\beta }}}$y el valor mínimo local${\displaystyle |S({\hat {\beta }})|}$es pequeño. La convergencia es cuadrática si${\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. Para un valor mínimo grande${\displaystyle |S({\hat {\beta }})|}$Sin embargo, la convergencia no está garantizada, ni siquiera la convergencia local como en el método de Newton , ni la convergencia bajo las condiciones habituales de Wolfe. ^{[ 6 ]}

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

Para${\displaystyle \lambda <1}$,${\displaystyle \beta =0}$es un óptimo local. Si${\displaystyle \lambda =0}$En ese caso, el problema es lineal y el método encuentra el óptimo en una iteración. Si |λ| < 1, el método converge linealmente y el error disminuye asintóticamente con un factor |λ| en cada iteración. Sin embargo, si |λ| > 1, el método ni siquiera converge localmente. ^{[ 8 ]}

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

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

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. En consecuencia, la tasa de convergencia del algoritmo de Gauss-Newton puede ser cuadrática bajo ciertas condiciones de regularidad. En general (bajo condiciones más débiles), la tasa de convergencia es lineal. ^{[ 9 ]}

Relación de recurrencia para el método de Newton para minimizar una función S de parámetros${\displaystyle {\boldsymbol {\beta }}}$es $${\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 hessiana de S.

Desde${\textstyle S=\sum _{i=1}^{m}r_{i}^{2}}$, el gradiente viene dado por $${\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,${\displaystyle g_{j}}$, con respecto a${\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 de derivada de segundo orden (el segundo término en esta expresión). Es decir, el hessiano se aproxima mediante $${\displaystyle H_{jk}\approx 2\sum _{i=1}^{m}J_{ij}J_{ik},}$$

dónde${\textstyle J_{ij}={\partial r_{i}}/{\partial \beta _{j}}}$son entradas del jacobiano J _r . Nótese que cuando el hessiano exacto se evalúa cerca de un ajuste exacto tenemos casi cero${\displaystyle r_{i}}$, por lo que el segundo término también se vuelve casi cero, lo que justifica la aproximación. El gradiente y el hessiano aproximado se pueden escribir en notación matricial como $${\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 operacionales. $${\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 de 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 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${\displaystyle r_{i}}$son de pequeña magnitud, al menos en torno al mínimo.
2. Las funciones son solo "ligeramente" no lineales, de modo que${\textstyle {\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k}}}}$es relativamente pequeño en magnitud.

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${\displaystyle S\left({\boldsymbol {\beta }}^{s}\right)}$es un punto estacionario, sostiene que${\displaystyle S\left({\boldsymbol {\beta }}^{s}+\alpha \Delta \right)<S\left({\boldsymbol {\beta }}^{s}\right)}$para todos suficientemente pequeños${\displaystyle \alpha >0}$Por lo tanto, si se produce una divergencia, una solución es emplear una fracción.${\displaystyle \alpha }$del vector de incremento Δ en la fórmula de actualización: $${\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. Un valor óptimo para${\displaystyle \alpha }$se puede encontrar utilizando un algoritmo de búsqueda lineal , es decir, la magnitud de${\displaystyle \alpha }$se determina encontrando el valor que minimiza S , generalmente utilizando un método de búsqueda directa en el intervalo${\displaystyle 0<\alpha <1}$o una búsqueda lineal con retroceso como la búsqueda lineal de Armijo . Normalmente,${\displaystyle \alpha }$debe elegirse de manera que satisfaga las condiciones de Wolfe o las condiciones de Goldstein . ^{[ 11 ]}

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 rota 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. Nótese que cuando D es la matriz identidad I y${\displaystyle \lambda \to +\infty }$, entonces${\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} }$Por lo tanto, la dirección de Δ se aproxima a la dirección del gradiente negativo.${\displaystyle -\mathbf {J} ^{\operatorname {T} }\mathbf {r} }$.

El llamado parámetro de Marquardt${\displaystyle \lambda }$También se puede optimizar mediante una búsqueda lineal, pero esto es ineficiente, ya que el vector de desplazamiento debe recalcularse cada vez.${\displaystyle \lambda }$se modifica. Una estrategia más eficiente es la siguiente: cuando se produce una divergencia, aumente el parámetro de Marquardt hasta que S disminuya . Luego, mantenga el valor de una iteración a la siguiente, pero redúzcalo si es posible hasta alcanzar un valor límite, momento en el que el parámetro de Marquardt puede establecerse en cero; la minimización de S se convierte entonces en una minimización estándar de Gauss-Newton.

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${\displaystyle \mathbf {J} _{\mathbf {r} }}$es más disperso que el Hessiano aproximado${\displaystyle \mathbf {J} _{\mathbf {r} }^{\operatorname {T} }\mathbf {J_{r}} }$En tales casos, el cálculo del paso en sí normalmente deberá realizarse con un método iterativo aproximado apropiado para problemas grandes y dispersos, como el método del gradiente conjugado .

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 el almacenamiento de matrices dispersas , en general es práctico almacenar las filas de${\displaystyle \mathbf {J} _{\mathbf {r} }}$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${\displaystyle \mathbf {J} _{\mathbf {r} }}$Se cumple la siguiente relación simple: $${\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 esta forma, cada fila contribuye de manera aditiva e independiente al producto. Además de respetar una estructura de almacenamiento dispersa práctica, esta expresión es idónea para cálculos paralelos . Cabe destacar que cada fila c _i es el gradiente del residuo r _i correspondiente ; teniendo esto en cuenta, la fórmula anterior subraya que los residuos contribuyen al problema de forma independiente entre sí.

En un método cuasi-Newton , como el de Davidon, Fletcher y Powell o el de Broyden-Fletcher-Goldfarb-Shanno ( método BFGS ), se obtiene una estimación de la matriz hessiana completa.${\textstyle {\frac {\partial ^{2}S}{\partial \beta _{j}\partial \beta _{k}}}}$se construye numéricamente utilizando primeras derivadas${\textstyle {\frac {\partial r_{i}}{\partial \beta _{j}}}}$De esta forma, tras n ciclos de refinamiento, el método se aproxima lo máximo posible al método de Newton en cuanto a rendimiento. Cabe destacar que los métodos cuasi-Newton pueden minimizar funciones reales generales, mientras que los métodos de Gauss-Newton, Levenberg-Marquardt, etc., solo son adecuados para problemas de mínimos cuadrados no lineales.

Otro método para resolver problemas de minimización utilizando únicamente derivadas de primer orden es el descenso de gradiente . Sin embargo, este método no tiene en cuenta las derivadas de segundo orden, ni siquiera de forma aproximada. Por consiguiente, resulta muy ineficiente para muchas funciones, especialmente si los parámetros presentan fuertes interacciones.

• 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.{{cite book}}: CS1 mantenimiento: ubicación del editor ( enlace )

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