Algoritmo de Levenberg-Marquardt

Algoritmo utilizado para resolver problemas de mínimos cuadrados no lineales

En matemáticas e informática, el algoritmo de Levenberg-Marquardt ( LMA o simplemente LM ), también conocido como el método de mínimos cuadrados amortiguados ( DLS ), se utiliza para resolver problemas de mínimos cuadrados no lineales . Estos problemas de minimización surgen especialmente en el ajuste de curvas de mínimos cuadrados . El LMA interpola entre el algoritmo de Gauss-Newton (GNA) y el método de descenso de gradiente . El LMA es más robusto que el GNA, lo que significa que en muchos casos encuentra una solución incluso si comienza muy lejos del mínimo final. Para funciones de buen comportamiento y parámetros iniciales razonables, el LMA tiende a ser más lento que el GNA. El LMA también puede verse como Gauss-Newton utilizando un enfoque de región de confianza .

El algoritmo fue publicado por primera vez en 1944 por Kenneth Levenberg ^[1] mientras trabajaba en el Arsenal del Ejército de Frankford . Fue redescubierto en 1963 por Donald Marquardt ^[2] , que trabajaba como estadístico en DuPont , y de forma independiente por Girard, ^[3] Wynne ^[4] y Morrison ^{[5] .}

El LMA se utiliza en muchas aplicaciones de software para resolver problemas genéricos de ajuste de curvas. Al utilizar el algoritmo de Gauss-Newton, a menudo converge más rápido que los métodos de primer orden. ^[6] Sin embargo, al igual que otros algoritmos de optimización iterativos, el LMA solo encuentra un mínimo local , que no es necesariamente el mínimo global .

El problema

La aplicación principal del algoritmo de Levenberg-Marquardt es el problema de ajuste de curvas de mínimos cuadrados: dado un conjunto de pares empíricos de variables independientes y dependientes, encuentre los parámetros de la curva del modelo de modo que la suma de los cuadrados de las desviaciones se minimice: ${\displaystyle m}$ ${\displaystyle \left(x_{i},y_{i}\right)}$ ${\displaystyle {\boldsymbol {\beta }}}$ ${\displaystyle f\left(x,{\boldsymbol {\beta }}\right)}$ ${\displaystyle S\left({\boldsymbol {\beta }}\right)}$

${\displaystyle {\hat {\boldsymbol {\beta }}}\in \operatorname {argmin} \limits _{\boldsymbol {\beta }}S\left({\boldsymbol {\beta }}\right)\equiv \operatorname {argmin} \limits _{\boldsymbol {\beta }}\sum _{i=1}^{m}\left[y_{i}-f\left(x_{i},{\boldsymbol {\beta }}\right)\right]^{2},}$que se supone que no está vacío.

La solución

Al igual que otros algoritmos de minimización numérica, el algoritmo de Levenberg-Marquardt es un procedimiento iterativo . Para iniciar una minimización, el usuario debe proporcionar una estimación inicial para el vector de parámetros ⁠ ⁠ ${\displaystyle {\boldsymbol {\beta }}}$ . En los casos con un solo mínimo, una estimación estándar no informada como funcionará bien; en los casos con múltiples mínimos , el algoritmo converge al mínimo global solo si la estimación inicial ya está algo cerca de la solución final. ${\displaystyle {\boldsymbol {\beta }}^{\text{T}}={\begin{pmatrix}1,\ 1,\ \dots ,\ 1\end{pmatrix}}}$

En cada paso de iteración, el vector de parámetros ⁠ ⁠ ${\displaystyle {\boldsymbol {\beta }}}$ se reemplaza por una nueva estimación ⁠ ⁠ ${\displaystyle {\boldsymbol {\beta }}+{\boldsymbol {\delta }}}$ . Para determinar ⁠ ⁠ ${\displaystyle {\boldsymbol {\delta }}}$ , la función se aproxima mediante su linealización : ${\displaystyle f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)}$

${\displaystyle f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)\approx f\left(x_{i},{\boldsymbol {\beta }}\right)+\mathbf {J} _{i}{\boldsymbol {\delta }},}$

dónde

${\displaystyle \mathbf {J} _{i}={\frac {\partial f\left(x_{i},{\boldsymbol {\beta }}\right)}{\partial {\boldsymbol {\beta }}}}}$

es el gradiente (vector fila en este caso) de ⁠ ⁠ ${\displaystyle f}$ con respecto a ⁠ ⁠ ${\displaystyle {\boldsymbol {\beta }}}$ .

La suma de las desviaciones cuadradas tiene su mínimo en un gradiente cero con respecto a ⁠ ⁠ . La aproximación de primer orden anterior de da ${\displaystyle S\left({\boldsymbol {\beta }}\right)}$ ${\displaystyle {\boldsymbol {\beta }}}$ ${\displaystyle f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)}$

${\displaystyle S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)\approx \sum _{i=1}^{m}\left[y_{i}-f\left(x_{i},{\boldsymbol {\beta }}\right)-\mathbf {J} _{i}{\boldsymbol {\delta }}\right]^{2},}$

o en notación vectorial,

${\displaystyle {\begin{aligned}S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)&\approx \left\|\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)-\mathbf {J} {\boldsymbol {\delta }}\right\|^{2}\\&=\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)-\mathbf {J} {\boldsymbol {\delta }}\right]^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)-\mathbf {J} {\boldsymbol {\delta }}\right]\\&=\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]-\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}-\left(\mathbf {J} {\boldsymbol {\delta }}\right)^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]+{\boldsymbol {\delta }}^{\mathrm {T} }\mathbf {J} ^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}\\&=\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]-2\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}+{\boldsymbol {\delta }}^{\mathrm {T} }\mathbf {J} ^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}.\end{aligned}}}$

Tomando la derivada de esta aproximación de con respecto a ⁠ ⁠ y fijando el resultado en cero se obtiene ${\displaystyle S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)}$ ${\displaystyle {\boldsymbol {\delta }}}$

${\displaystyle \left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} \right){\boldsymbol {\delta }}=\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right],}$

donde es la matriz jacobiana , cuya fila ⁠ ⁠ -ésima es igual a , y donde y son vectores con el componente ⁠ ⁠ -ésimo y respectivamente. La expresión anterior obtenida para ⁠ ⁠ se aplica al método de Gauss-Newton. La matriz jacobiana definida anteriormente no es (en general) una matriz cuadrada, sino una matriz rectangular de tamaño , donde es el número de parámetros (tamaño del vector ). La multiplicación de matrices produce la matriz cuadrada requerida y el producto matriz-vector en el lado derecho produce un vector de tamaño . El resultado es un conjunto de ecuaciones lineales, que se pueden resolver para ⁠ ⁠ . ${\displaystyle \mathbf {J} }$ ${\displaystyle i}$ ${\displaystyle \mathbf {J} _{i}}$ ${\displaystyle \mathbf {f} \left({\boldsymbol {\beta }}\right)}$ ${\displaystyle \mathbf {y} }$ ${\displaystyle i}$ ${\displaystyle f\left(x_{i},{\boldsymbol {\beta }}\right)}$ ${\displaystyle y_{i}}$ ${\displaystyle {\boldsymbol {\beta }}}$ ${\displaystyle m\times n}$ ${\displaystyle n}$ ${\displaystyle {\boldsymbol {\beta }}}$ ${\displaystyle \left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} \right)}$ ${\displaystyle n\times n}$ ${\displaystyle n}$ ${\displaystyle n}$ ${\displaystyle {\boldsymbol {\delta }}}$

La contribución de Levenberg es reemplazar esta ecuación por una "versión amortiguada":

${\displaystyle \left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} +\lambda \mathbf {I} \right){\boldsymbol {\delta }}=\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right],}$

donde ⁠ ⁠ ${\displaystyle \mathbf {I} }$ es la matriz identidad, que da como incremento ⁠ ⁠ ${\displaystyle {\boldsymbol {\delta }}}$ al vector de parámetros estimado ⁠ ⁠ ${\displaystyle {\boldsymbol {\beta }}}$ .

El factor de amortiguamiento (no negativo) ⁠ ⁠ ${\displaystyle \lambda }$ se ajusta en cada iteración. Si la reducción de ⁠ ⁠ ${\displaystyle S}$ es rápida, se puede usar un valor más pequeño, acercando el algoritmo al algoritmo de Gauss-Newton , mientras que si una iteración da una reducción insuficiente en el residuo, ⁠ ⁠ ${\displaystyle \lambda }$ se puede aumentar, dando un paso más cercano a la dirección de descenso del gradiente. Tenga en cuenta que el gradiente de ⁠ ⁠ ${\displaystyle S}$ con respecto a ⁠ ⁠ ${\displaystyle {\boldsymbol {\beta }}}$ es igual a . Por lo tanto, para valores grandes de ⁠ ⁠ , el paso se dará aproximadamente en la dirección opuesta al gradiente. Si la longitud del paso calculado ⁠ ⁠ o la reducción de la suma de cuadrados del último vector de parámetros ⁠ ⁠ caen por debajo de los límites predefinidos, la iteración se detiene y el último vector de parámetros ⁠ ⁠ se considera la solución. ${\displaystyle -2\left(\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]\right)^{\mathrm {T} }}$ ${\displaystyle \lambda }$ ${\displaystyle {\boldsymbol {\delta }}}$ ${\displaystyle {\boldsymbol {\beta }}+{\boldsymbol {\delta }}}$ ${\displaystyle {\boldsymbol {\beta }}}$

Cuando el factor de amortiguamiento ⁠ ⁠ ${\displaystyle \lambda }$ es grande en relación con , no es necesario invertir, ya que la actualización se aproxima bien mediante el pequeño paso de gradiente . ${\displaystyle \|\mathbf {J} ^{\mathrm {T} }\mathbf {J} \|}$ ${\displaystyle \mathbf {J} ^{\mathrm {T} }\mathbf {J} +\lambda \mathbf {I} }$ ${\displaystyle \lambda ^{-1}\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]}$

Para que la solución fuera invariante en escala, el algoritmo de Marquardt resolvió un problema modificado con cada componente del gradiente escalado de acuerdo con la curvatura. Esto proporciona un mayor movimiento a lo largo de las direcciones donde el gradiente es menor, lo que evita una convergencia lenta en la dirección del gradiente pequeño. Fletcher, en su artículo de 1971 A modified Marquardt subrutine for non-linear minimum squares (Una subrutina Marquardt modificada para mínimos cuadrados no lineales), simplificó la forma, reemplazando la matriz identidad ⁠ ⁠ ${\displaystyle \mathbf {I} }$ por la matriz diagonal que consiste en los elementos diagonales de ⁠ ⁠ ${\displaystyle \mathbf {J} ^{\text{T}}\mathbf {J} }$ :

${\displaystyle \left[\mathbf {J} ^{\mathrm {T} }\mathbf {J} +\lambda \operatorname {diag} \left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} \right)\right]{\boldsymbol {\delta }}=\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right].}$

Un factor de amortiguamiento similar aparece en la regularización de Tikhonov , que se utiliza para resolver problemas lineales mal planteados , así como en la regresión de cresta , una técnica de estimación en estadística .

Elección del parámetro de amortiguación

Se han propuesto varios argumentos más o menos heurísticos para la mejor elección del parámetro de amortiguamiento . Existen argumentos teóricos que muestran por qué algunas de estas elecciones garantizan la convergencia local del algoritmo; sin embargo, estas elecciones pueden ${\displaystyle \lambda }$ hacer que la convergencia global del algoritmo sufra las propiedades indeseables del descenso más pronunciado , en particular, una convergencia muy lenta cerca del óptimo.

Los valores absolutos de cualquier elección dependen de qué tan bien escalado esté el problema inicial. Marquardt recomendó comenzar con un valor ⁠ ⁠ ${\displaystyle \lambda _{0}}$ y un factor ⁠ ⁠ ${\displaystyle \nu >1}$ . Inicialmente, se establece y calcula la suma residual de cuadrados después de un paso desde el punto de partida con el factor de amortiguamiento de y, en segundo lugar, con ⁠ ⁠ . Si ambos son peores que el punto inicial, entonces el amortiguamiento se incrementa mediante multiplicaciones sucesivas por ⁠ ⁠ hasta que se encuentra un punto mejor con un nuevo factor de amortiguamiento de ⁠ ⁠ para algún ⁠ ⁠ . ${\displaystyle \lambda =\lambda _{0}}$ ${\displaystyle S\left({\boldsymbol {\beta }}\right)}$ ${\displaystyle \lambda =\lambda _{0}}$ ${\displaystyle \lambda _{0}/\nu }$ ${\displaystyle \nu }$ ${\displaystyle \lambda _{0}\nu ^{k}}$ ${\displaystyle k}$

Si el uso del factor de amortiguamiento ⁠ ⁠ ${\displaystyle \lambda /\nu }$ da como resultado una reducción del cuadrado del residuo, entonces este se toma como el nuevo valor de ⁠ ⁠ ${\displaystyle \lambda }$ (y la nueva ubicación óptima se toma como la obtenida con este factor de amortiguamiento) y el proceso continúa; si el uso de ⁠ ⁠ ${\displaystyle \lambda /\nu }$ dio como resultado un residuo peor, pero el uso de ⁠ ⁠ ${\displaystyle \lambda }$ dio como resultado un residuo mejor, entonces ⁠ ⁠ ${\displaystyle \lambda }$ se deja sin cambios y el nuevo óptimo se toma como el valor obtenido con ⁠ ⁠ ${\displaystyle \lambda }$ como factor de amortiguamiento.

Una estrategia eficaz para el control del parámetro de amortiguamiento, llamada gratificación retrasada , consiste en aumentar el parámetro en una pequeña cantidad para cada paso ascendente y disminuirlo en una gran cantidad para cada paso descendente. La idea detrás de esta estrategia es evitar moverse cuesta abajo demasiado rápido al comienzo de la optimización, restringiendo así los pasos disponibles en futuras iteraciones y, por lo tanto, ralentizando la convergencia. ^[7] Se ha demostrado que un aumento por un factor de 2 y una disminución por un factor de 3 son eficaces en la mayoría de los casos, mientras que para problemas grandes pueden funcionar mejor valores más extremos, con un aumento por un factor de 1,5 y una disminución por un factor de 5. ^[8]

Aceleración geodésica

Al interpretar el paso de Levenberg-Marquardt como la velocidad a lo largo de una trayectoria geodésica en el espacio de parámetros, es posible mejorar el método agregando un término de segundo orden que tenga en cuenta la aceleración a lo largo de la trayectoria geodésica. ${\displaystyle {\boldsymbol {v}}_{k}}$ ${\displaystyle {\boldsymbol {a}}_{k}}$

${\displaystyle {\boldsymbol {v}}_{k}+{\frac {1}{2}}{\boldsymbol {a}}_{k}}$

¿Dónde está la solución de? ${\displaystyle {\boldsymbol {a}}_{k}}$

${\displaystyle {\boldsymbol {J}}_{k}{\boldsymbol {a}}_{k}=-f_{vv}.}$

Dado que este término de aceleración geodésica depende únicamente de la derivada direccional a lo largo de la dirección de la velocidad , no requiere calcular la matriz de derivadas de segundo orden completa, lo que requiere solo una pequeña sobrecarga en términos de costo computacional. ^[9] Dado que la derivada de segundo orden puede ser una expresión bastante compleja, puede ser conveniente reemplazarla con una aproximación de diferencias finitas. ${\displaystyle f_{vv}=\sum _{\mu \nu }v_{\mu }v_{\nu }\partial _{\mu }\partial _{\nu }f({\boldsymbol {x}})}$ ${\displaystyle {\boldsymbol {v}}}$

${\displaystyle {\begin{aligned}f_{vv}^{i}&\approx {\frac {f_{i}({\boldsymbol {x}}+h{\boldsymbol {\delta }})-2f_{i}({\boldsymbol {x}})+f_{i}({\boldsymbol {x}}-h{\boldsymbol {\delta }})}{h^{2}}}\\&={\frac {2}{h}}\left({\frac {f_{i}({\boldsymbol {x}}+h{\boldsymbol {\delta }})-f_{i}({\boldsymbol {x}})}{h}}-{\boldsymbol {J}}_{i}{\boldsymbol {\delta }}\right)\end{aligned}}}$

donde y ya han sido calculados por el algoritmo, por lo que se requiere solo una evaluación de función adicional para calcular . La elección del paso de diferencia finita puede afectar la estabilidad del algoritmo, y un valor de alrededor de 0,1 suele ser razonable en general. ^[8] ${\displaystyle f({\boldsymbol {x}})}$ ${\displaystyle {\boldsymbol {J}}}$ ${\displaystyle f({\boldsymbol {x}}+h{\boldsymbol {\delta }})}$ ${\displaystyle h}$

Dado que la aceleración puede apuntar en dirección opuesta a la velocidad, para evitar que se detenga el método en caso de que la amortiguación sea demasiado pequeña, se agrega un criterio adicional sobre la aceleración para aceptar un paso, que requiere que

${\displaystyle {\frac {2\left\|{\boldsymbol {a}}_{k}\right\|}{\left\|{\boldsymbol {v}}_{k}\right\|}}\leq \alpha }$

donde generalmente se fija en un valor menor que 1, con valores más pequeños para problemas más difíciles. ^[8] ${\displaystyle \alpha }$

La adición de un término de aceleración geodésica puede permitir un aumento significativo en la velocidad de convergencia y es especialmente útil cuando el algoritmo se mueve a través de cañones estrechos en el paisaje de la función objetivo, donde los pasos permitidos son más pequeños y la mayor precisión debido al término de segundo orden brinda mejoras significativas. ^[8]

Ejemplo

Mal ajuste

Mejor ajuste

Mejor ajuste

En este ejemplo, intentamos ajustar la función utilizando el algoritmo Levenberg–Marquardt implementado en GNU Octave como la función leasqr . Los gráficos muestran un ajuste progresivamente mejor para los parámetros , utilizados en la curva inicial. Solo cuando los parámetros en el último gráfico se eligen más cercanos al original, las curvas se ajustan exactamente. Esta ecuación es un ejemplo de condiciones iniciales muy sensibles para el algoritmo Levenberg–Marquardt. Una razón para esta sensibilidad es la existencia de múltiples mínimos: la función tiene mínimos en el valor del parámetro y . ${\displaystyle y=a\cos \left(bX\right)+b\sin \left(aX\right)}$ ${\displaystyle a=100}$ ${\displaystyle b=102}$ ${\displaystyle \cos \left(\beta x\right)}$ ${\displaystyle {\hat {\beta }}}$ ${\displaystyle {\hat {\beta }}+2n\pi }$

Véase también

• Región de confianza
• Método Nelder-Mead
• También se han utilizado variantes del algoritmo de Levenberg-Marquardt para resolver sistemas de ecuaciones no lineales. ^[10]

Referencias

1. Levenberg, Kenneth (1944). "Un método para la solución de ciertos problemas no lineales mediante mínimos cuadrados". Quarterly of Applied Mathematics . 2 (2): 164–168. doi : 10.1090/qam/10666 .
2. Marquardt, Donald (1963). "Un algoritmo para la estimación de parámetros no lineales por mínimos cuadrados". Revista SIAM de Matemáticas Aplicadas . 11 (2): 431–441. doi :10.1137/0111030. hdl : 10338.dmlcz/104299 .
3. Girard, André (1958). "Extracto de Revue d'optique théorique et instrumentale ". Rev. Optar . 37 : 225–241, 397–424.
4. Wynne, CG (1959). "Diseño de lentes mediante computadora digital electrónica: I". Proc. Phys. Soc. Lond . 73 (5): 777–787. Bibcode :1959PPS....73..777W. doi :10.1088/0370-1328/73/5/310.
5. Morrison, David D. (1960). "Métodos para problemas de mínimos cuadrados no lineales y pruebas de convergencia". Actas del Seminario del Laboratorio de Propulsión a Chorro sobre Programas de Seguimiento y Determinación de Órbitas : 1–9.
6. Wiliamowski, Bogdan; Yu, Hao (junio de 2010). "Computación mejorada para el entrenamiento de Levenberg–Marquardt" (PDF) . Transacciones IEEE sobre redes neuronales y sistemas de aprendizaje . 21 (6).
7. Transtrum, Mark K; Machta, Benjamin B; Sethna, James P (2011). "Geometría de mínimos cuadrados no lineales con aplicaciones a modelos descuidados y optimización". Physical Review E . 83 (3). APS: 036701. arXiv : 1010.1449 . Bibcode :2011PhRvE..83c6701T. doi :10.1103/PhysRevE.83.036701. PMID  21517619. S2CID  15361707.
8. Transtrum, Mark K; Sethna, James P (2012). "Mejoras en el algoritmo Levenberg-Marquardt para la minimización de mínimos cuadrados no lineales". arXiv : 1201.5885 [physics.data-an].
9. "Ajuste de mínimos cuadrados no lineales". Biblioteca científica GNU. Archivado desde el original el 14 de abril de 2020.
10. Kanzow, Christian; Yamashita, Nobuo; Fukushima, Masao (2004). "Métodos de Levenberg-Marquardt con fuertes propiedades de convergencia local para resolver ecuaciones no lineales con restricciones convexas". Revista de Matemática Computacional y Aplicada . 172 (2): 375–397. Bibcode :2004JCoAM.172..375K. doi : 10.1016/j.cam.2004.02.013 .

Lectura adicional

• Moré, Jorge J.; Sorensen, Daniel C. (1983). "Computing a Trust-Region Step" (PDF) . SIAM J. Sci. Stat. Comput . 4 (3): 553–572. doi :10.1137/0904038.
• Gill, Philip E.; Murray, Walter (1978). "Algoritmos para la solución del problema de mínimos cuadrados no lineales". Revista SIAM de Análisis Numérico . 15 (5): 977–992. Bibcode :1978SJNA...15..977G. doi :10.1137/0715063.
• Pujol, Jose (2007). "La solución de problemas inversos no lineales y el método de Levenberg-Marquardt". Geofísica . 72 (4). SEG: W1–W16. Bibcode :2007Geop...72W...1P. doi :10.1190/1.2732552.
• Nocedal, Jorge; Wright, Stephen J. (2006). Optimización numérica (2.ª ed.). Springer. ISBN 978-0-387-30303-1.

Enlaces externos

• Se puede encontrar una descripción detallada del algoritmo en Recetas numéricas en C, Capítulo 15.5: Modelos no lineales.
• CT Kelley, Métodos iterativos para optimización , SIAM Frontiers in Applied Mathematics, n.º 18, 1999, ISBN 0-89871-433-8 . Copia en línea 
• Historia del algoritmo en las noticias de SIAM
• Un tutorial de Ananth Ranganathan
• K. Madsen, HB Nielsen, O. Tingleff, Métodos para problemas de mínimos cuadrados no lineales (tutorial de mínimos cuadrados no lineales; código LM: secante jacobiana analítica)
• T. Strutz: Ajuste de datos e incertidumbre (Una introducción práctica a los mínimos cuadrados ponderados y más allá). 2.ª edición, Springer Vieweg, 2016, ISBN 978-3-658-11455-8 . 
• HP Gavin, El método Levenberg-Marquardt para problemas de ajuste de curvas por mínimos cuadrados no lineales ( implementación en MATLAB incluida)