Ver el código de «Algoritmo de Levenberg-Marquardt»

En [[matemáticas]] y computación, 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 se usa en muchas aplicaciones de software para resolver problemas genéricos de ajuste de curvas. Sin embargo, como ocurre con muchos algoritmos de ajuste, el LMA solo encuentra un [[Extremos de una función|mínimo local]], que no es necesariamente el [[Extremos de una función|mínimo global]]. El LMA interpola entre el [[algoritmo de Gauss-Newton]] (GNA) y el método de  [[Gradient descent|descenso de]] [[gradiente]]. El LMA es más [[Robustness (computer science)|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 de inicio razonables, el LMA tiende a ser un poco más lento que el GNA. El LMA también se puede ver como [[Algoritmo de Gauss-Newton|Gauss-Newton]] utilizando un enfoque de [[región de confianza]].

El algoritmo fue publicado por primera vez en 1944 por [[Kenneth Levenberg]],<ref>{{Cita publicación|url=https://www.ams.org/qam/1944-02-02/S0033-569X-1944-10666-0/|título=A method for the solution of certain non-linear problems in least squares|apellidos=Levenberg|nombre=Kenneth|fecha=1944|publicación=Quarterly of Applied Mathematics|volumen=2|número=2|páginas=164–168|idioma=en|issn=0033-569X|doi=10.1090/qam/10666|pmid=}}</ref> mientras trabajaba en el Arsenal del Ejército de Frankford. Fue redescubierto en 1963 por [[Donald Marquardt]],<ref>{{Cita publicación|url=https://epubs.siam.org/doi/10.1137/0111030|título=An Algorithm for Least-Squares Estimation of Nonlinear Parameters|apellidos=Marquardt|nombre=Donald W.|fecha=1 de junio de 1963|publicación=Journal of the Society for Industrial and Applied Mathematics|volumen=11|número=2|páginas=431–441|issn=0368-4245|doi=10.1137/0111030|pmid=}}</ref> quien trabajó como [[estadístico]] en [[DuPont Corporation|DuPont]], e independientemente por Girard<ref>{{Cita publicación|url=|título=Excerpt from Revue d'optique théorique et instrumentale|apellidos=Girard|nombre=André|fecha=1958|publicación=Revue d'optique théorique et instrumentale|volumen=37|páginas=225–241, 397–424|fechaacceso=|doi=|pmid=}}</ref> Wynne<ref>{{Cita publicación|url=https://doi.org/10.1088/0370-1328/73/5/310|título=Lens Designing by Electronic Digital Computer: I|apellidos=Wynne|nombre=C G|fecha=1 de mayo de 1959|publicación=Proceedings of the Physical Society|volumen=73|número=5|páginas=777–787|idioma=en|issn=0370-1328|doi=10.1088/0370-1328/73/5/310|pmid=}}</ref> y Morrison.<ref>{{Cita publicación|url=|título=Methods for nonlinear least squares problems and convergence proofs|apellidos=Morrison|nombre=David|fecha=1960|publicación=Proceedings of the Jet Propulsion Laboratory Seminar on Tracking Programs and Orbit Determination|páginas=1-9|fechaacceso=|doi=|pmid=}}</ref>

== El problema ==
La aplicación principal del algoritmo de Levenberg-Marquardt se encuentra en el problema de ajuste de curvas de mínimos cuadrados: dado un conjunto de <math>m</math>pares de datos empíricos ( ''x <sub>i</sub>'' , ''y <sub>i</sub>'' ) de variables independientes y dependientes, encuentre los parámetros '''''β''''' de la curva del modelo ''f'' ( ''x'' , '''''β''''' ) de modo que la suma de los cuadrados de las desviaciones ''S'' ( '''''β''''' ) se minimice:

: <math>\hat{\boldsymbol\beta} \in \operatorname{argmin}\limits_{\boldsymbol\beta} S(\boldsymbol\beta) \equiv \operatorname{argmin}\limits_{\boldsymbol\beta} \sum_{i=1}^m [y_i - f(x_i, \boldsymbol\beta)]^2,</math> 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 [[Iteración|iterativo]]. Para iniciar una minimización, el usuario debe proporcionar una estimación inicial del vector de parámetros '''''β'''''. En los casos con solo un mínimo, una conjetura estándar no informada como '''''β''''' <sup>T</sup> = (1, 1, ..., 1) funcionará bien; en los casos con [[Extremos de una función|múltiples mínimos]], el algoritmo converge al mínimo global solo si la estimación inicial ya es algo cercana a la solución final.

En cada paso de iteración, el vector parámetro '''''β''''' se reemplaza por una nueva estimación '''''β''''' + '''''δ'''''. Para determinar '''''δ''''' , la función <math>f(x_i, \boldsymbol\beta + \boldsymbol\delta)</math>se aproxima por su linealización:

: <math>f(x_i, \boldsymbol\beta + \boldsymbol\delta) \approx f(x _i, \boldsymbol\beta) + \mathbf J_i \boldsymbol\delta,</math>

donde

: <math>\mathbf J_i = \frac{\partial f(x_i, \boldsymbol\beta)}{\partial \boldsymbol\beta}</math>

es el [[gradiente]] (fila-vector en este caso) de ''f'' con respecto a '''''β''''' .

La suma <math>S(\boldsymbol\beta)</math>de desviaciones cuadradas tiene su mínimo en un [[gradiente]] de cero con respecto a '''''β''''' . La anterior aproximación de primer orden de <math>f(x_i, \boldsymbol\beta + \boldsymbol\delta)</math>da

: <math>S(\boldsymbol\beta + \boldsymbol\delta) \approx \sum_{i=1}^m [y_i - f(x_i, \boldsymbol\beta) - \mathbf J_i \boldsymbol\delta]^2,</math>

o en notación vectorial,

: <math>\begin{align}
 S(\boldsymbol\beta + \boldsymbol\delta) &\approx \|\mathbf y - \mathbf f(\boldsymbol\beta) - \mathbf J\boldsymbol\delta\|^2\\
  &= [\mathbf y - \mathbf f(\boldsymbol\beta) - \mathbf J\boldsymbol\delta]^{\mathrm T} [\mathbf y - \mathbf f(\boldsymbol\beta) - \mathbf J\boldsymbol\delta]\\
  &= [\mathbf y - \mathbf f(\boldsymbol\beta)]^{\mathrm T} [\mathbf y - \mathbf f(\boldsymbol\beta)] - [\mathbf y - \mathbf f(\boldsymbol\beta)]^{\mathrm T} \mathbf J \boldsymbol\delta - (\mathbf J \boldsymbol\delta)^{\mathrm T} [\mathbf y - \mathbf f(\boldsymbol\beta)] + \boldsymbol\delta^{\mathrm T} \mathbf J^{\mathrm T} \mathbf J \boldsymbol\delta\\
  &= [\mathbf y - \mathbf f(\boldsymbol\beta)]^{\mathrm T} [\mathbf y - \mathbf f(\boldsymbol\beta)] - 2[\mathbf y - \mathbf f(\boldsymbol\beta)]^{\mathrm T} \mathbf J \boldsymbol\delta + \boldsymbol\delta^{\mathrm T} \mathbf J^{\mathrm T} \mathbf J\boldsymbol\delta.
\end{align}</math>

Tomando la derivada de <math>S(\boldsymbol\beta + \boldsymbol\delta)</math>con respecto a '''δ''' y establecer el resultado a cero da

: <math>(\mathbf J^{\mathrm T} \mathbf J)\boldsymbol\delta = \mathbf J^{\mathrm T} [\mathbf y - \mathbf f(\boldsymbol\beta)],</math>

donde <math>\mathbf J</math> es la [[Jacobiano|matriz jacobiana]], cuya fila ''i'' es igual a <math>\mathbf J_i</math>, y donde <math>\mathbf f</math>y <math>\mathbf y</math>son vectores con ''i''-componente <math>f(x_i, \boldsymbol\beta)</math>y <math>y_i</math>respectivamente. Este es un conjunto de ecuaciones lineales, que se pueden resolver para '''''δ''''' .

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

: <math>(\mathbf J^{\mathrm T} \mathbf J + \lambda \mathbf I) \boldsymbol\delta = \mathbf J^{\mathrm T} [\mathbf y - \mathbf f(\boldsymbol\beta)],</math>

donde '''I''' es la matriz de identidad, dando como incremento '''''δ''''' al vector de parámetros estimado '''''β'''''.

El factor de amortiguamiento (no negativo) λ se ajusta en cada iteración. Si la reducción de ''S'' es rápida, se puede usar un valor más pequeño, acercando el algoritmo al algoritmo de [[Algoritmo de Gauss-Newton|Gauss-Newton]], mientras que si una iteración produce una reducción insuficiente del residual, λ puede aumentarse, dando un paso más cerca del descenso del gradiente dirección. Tenga en cuenta que el [[gradiente]] de ''S'' con respecto a '''''β''''' es igual a <math>-2\big(\mathbf J^{\mathrm T} [\mathbf y - \mathbf f(\boldsymbol\beta)]\big)^{\mathrm T}</math>. Por lo tanto, para valores grandes de ''λ'', el paso se tomará aproximadamente en la dirección del gradiente. Si la longitud del paso calculado '''''δ''''' o la reducción de la suma de cuadrados del último vector de parámetros '''''β''''' + '''''δ''''' cae por debajo de los límites predefinidos, la iteración se detiene, y el último vector de parámetros '''''β''''' se considera la solución.

Si el valor del factor de amortiguamiento λ es grande en relación con la norma de '''J'''<sup>T</sup>'''J''', no es necesario resolver el sistema con matriz '''J'''<sup>T</sup>'''J''' + λ'''I''' , dado que esta matriz se puede aproximar mediante λ'''I'''. En este caso el algoritmo realiza una actualización en la dirección del gradiente, con un paso pequeño:  <math>\delta = \frac{1}{\lambda} \mathbf J^{\mathrm T} [\mathbf y - \mathbf f(\boldsymbol\beta)].</math>

R. Fletcher proporcionó la idea de que podemos escalar cada componente del gradiente según la curvatura, de modo que haya un movimiento más grande a lo largo de las direcciones donde el gradiente es más pequeño. Esto evita la convergencia lenta en la dirección del gradiente pequeño. Por lo tanto, Fletcher en su artículo de 1971 "''Una subrutina de Marquardt modificada para mínimos cuadrados no lineales''", reemplazó la matriz de identidad '''I''' con la matriz diagonal que consiste en los elementos diagonales de '''J'''<sup>T</sup>'''J''', lo que hace que la escala de la solución sea invariante:

: <math>[\mathbf J^{\mathrm T} \mathbf J + \lambda \operatorname{diag}(\mathbf J^{\mathrm T} \mathbf J)] \boldsymbol\delta = \mathbf J^{\mathrm T} [\mathbf y - \mathbf f(\boldsymbol\beta)].</math>

Un factor de amortiguación similar aparece en la [[Regularización de Tíjonov|regularización de Tikhonov]], que se utiliza para resolver problemas lineales [[Problema bien definido|mal planteados]], así como en la [[Regularización de Tíjonov|regresión de crestas]], una técnica de [[Estimación estadística|estimación]] en [[estadística]].

=== Elección del parámetro de amortiguación ===
Se han presentado varios argumentos más o menos heurísticos para la mejor opción para el parámetro de amortiguación λ. Existen argumentos teóricos que muestran por qué algunas de estas opciones garantizan la convergencia local del algoritmo; sin embargo, estas elecciones pueden hacer que la convergencia global del algoritmo sufra de las propiedades indeseables del [[Gradient descent|descenso más pronunciado]], en particular, una convergencia muy lenta cerca del óptimo.

Los valores absolutos de cualquier elección dependen de qué tan escalada esté el problema inicial. Marquardt recomendó comenzar con un valor λ <sub>0</sub> y un factor ν>1. Configurando inicialmente λ = λ<sub>0</sub> y calculando la suma residual de cuadrados ''S'' ('''''β''''') después de un paso desde el punto de inicio con el factor de amortiguamiento de λ = λ<sub>0</sub> y en segundo lugar con λ<sub>0</sub>/ν. Si ambos de estos son peores que el punto inicial, entonces la amortiguación se incrementa por multiplicación sucesiva por ν hasta que se encuentra un punto mejor con un nuevo factor de amortiguación de λ<sub>0</sub>ν<sup>''k''</sup> para algunos ''k'' .

Si el uso del factor de amortiguación λ / ν da como resultado una reducción en el residuo al cuadrado, esto se toma como el nuevo valor de λ (y la nueva ubicación óptima se toma como la obtenida con este factor de amortiguación) y el proceso continúa; si el uso de λ/ν resultó en un peor residuo, pero el uso de λ resultó en un mejor residual, entonces λ se mantiene sin cambios y el nuevo óptimo se toma como el valor obtenido con λ como factor de amortiguamiento.

== Ejemplo ==
[[Archivo:Lev-Mar-poor-fit.png|miniaturadeimagen| Mal encajado ]]
[[Archivo:Lev-Mar-better-fit.png|miniaturadeimagen| Mejor ajuste ]]
[[Archivo:Lev-Mar-best-fit.png|miniaturadeimagen| Mejor ajuste ]]
En este ejemplo tratamos de ajustar la función <math>y = a \cos(bX) + b \sin(aX)</math>utilizando el algoritmo de Levenberg-Marquardt implementado en [[GNU Octave|Octave GNU]] como la función ''leasqr''. Las gráficas muestran una mejor adaptación progresiva para los parámetros ''a'' = 100, ''b'' = 102 utilizados en la curva inicial. Solo cuando los parámetros en el último gráfico son elegidos más cercanos al original, las curvas se ajustan exactamente. Esta ecuación es un ejemplo de condiciones iniciales muy sensibles para el algoritmo de Levenberg-Marquardt. Una de las razones de esta sensibilidad es la existencia de múltiples mínimos: la función <math>\cos(\beta x)</math>tiene mínimos en el valor del parámetro <math>\hat\beta</math>y <math>\hat\beta + 2n\pi</math>.

== Véase también ==

* [[Trust region|Región de confianza]]
* [[Método Nelder-Mead|Método Nelder – Mead]] (también conocido como simplex)
* Las variantes del algoritmo de Levenberg-Marquardt también se han utilizado para resolver sistemas de ecuaciones no lineales.<ref>Christian Kanzow, Nobuo Yamashita, Masao Fukushima "Levenberg-Marquardt methods with strong local convergence properties for solving nonlinear equations with convex constraints," JCAM, 172(2):375-97, Dec. 2004, {{Doi|10.1016/j.cam.2004.02.013}}.</ref>

== Referencias ==
<references responsive=""></references>

== Otras lecturas ==
{{Refbegin}}
* {{cite journal|last1=Moré|first1=Jorge J.|last2=Sorensen|first2=Daniel C.|year=1983|title=Computing a Trust-Region Step|journal=SIAM J. Sci. Stat. Comput.|pages=553–572|issue=4}}
* {{cite journal|last1=Gill|first1=Philip E.|last2=Murray|first2=Walter|authorlink2=Walter Murray (mathematician)|year=1978|title=Algorithms for the solution of the nonlinear least-squares problem|journal=[[SIAM Journal on Numerical Analysis]]|doi=10.1137/0715063|volume=15|issue=5|pages=977–992}}
* {{Enlace roto|1={{cite journal|last=Pujol|first=Jose|year=2007|title=The solution of nonlinear inverse problems and the Levenberg-Marquardt method|journal=Geophysics|publisher=SEG|doi=10.1190/1.2732552|volume=72|number=4|pages=W1–W16|url=http://link.aip.org/link/?GPY/72/W1/1}} |2=http://link.aip.org/link/?GPY%2F72%2FW1%2F1 |bot=InternetArchiveBot }}
* {{cite book|last1=Nocedal|first1=Jorge|last2=Wright|first2=Stephen J.|year=2006|title=Numerical Optimization|edition=2nd|publisher=Springer|isbn=0-387-30303-0}}
{{Refend}}

== Enlaces externos ==

=== Descripciones ===

* Detailed description of the algorithm can be found in [http://www.nrbook.com/a/bookcpdf.php Numerical Recipes in C, Chapter 15.5: Nonlinear models] {{Wayback|url=http://www.nrbook.com/a/bookcpdf.php |date=20120304160342 }}
* C. T. Kelley, ''Iterative Methods for Optimization'', SIAM Frontiers in Applied Mathematics, no 18, 1999, {{ISBN|0-89871-433-8}}. [http://www.siam.org/books/textbooks/fr18_book.pdf Online copy]
* [https://web.archive.org/web/20140301154319/http://www3.villanova.edu/maple/misc/mtc1093.html History of the algorithm in SIAM news]
* [http://ananth.in/docs/lmtut.pdf A tutorial by Ananth Ranganathan]
* [http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf Methods for Non-Linear Least Squares Problems] by K. Madsen, H. B. Nielsen, O. Tingleff is a tutorial discussing nonlinear least-squares in general and the Levenberg–Marquardt method in particular
* T. Strutz: ''Data Fitting and Uncertainty (A practical introduction to weighted least squares and beyond).'' 2nd edition, Springer Vieweg, 2016, {{ISBN|978-3-658-11455-8}}.
* H. P. Gavin, [http://people.duke.edu/~hpgavin/ce281/lm.pdf ''The Levenberg-Marquardt method for nonlinear least squares curve-fitting problems''] ([[MATLAB]] implementation included)

=== Implementaciones ===
Levenberg-Marquardt es un algoritmo incorporado en entornos de computación numérica [[SciPy]], [[GNU Octave|Octave GNU]], [[Scilab]], [[Mathematica]], [[MATLAB|Matlab]], [[NeuroSolutions]], [[Origin (data analysis software)|Origin]], [[Fityk]], [[IGOR Pro]], [[LabVIEW]] y [[SAS (software)|SAS]]. También existen numerosas bibliotecas de software que permiten utilizar el algoritmo LM en aplicaciones independientes. Algunos de ellos admiten solo la optimización básica sin restricciones, mientras que otros admiten diferentes combinaciones de restricciones de caja y lineales.

Caja e implementaciones linealmente restringidas:

* ALGLIB tiene una [http://www.alglib.net/optimization/levenbergmarquardt.php implementación] limitada y lineal de LM mejorada en C # y C++. Algoritmo mejorado toma menos tiempo para converger y puede usar Jacobian o [[Matriz hessiana|Hessian]] exacto.
* Artelys Knitro es un solucionador no lineal con una implementación del algoritmo Levenberg-Marquardt restringido por caja. Está escrito en C y tiene interfaces con C ++ / C # / Java / Python / MATLAB / R.
* [http://ceres-solver.org/ ceres] es una biblioteca de minimización no lineal con una implementación del algoritmo Levenberg-Marquardt restringido por cuadros. Está escrito en C ++ y usa [http://eigen.tuxfamily.org/index.php?title=Main_Page eigen] .
* [[IDL (programming language)|IDL]] , el complemento [http://cow.physics.wisc.edu/~craigm/idl/fitting.html MPFIT] admite restricciones de cuadro.
* [http://www.ics.forth.gr/%7elourakis/levmar/ levmar] es una implementación en [[C (lenguaje de programación)|C]]/[[C++]] con soporte para restricciones de caja y lineales generales, distribuida bajo la [[GNU General Public License|Licencia Pública General de GNU]] .
** levmar incluye una interfaz de archivo MEX para [[MATLAB]] .
** [[Perl]] (PDL), [[python]] , [[Haskell]] y [[Microsoft .NET|.]] [[Microsoft .NET|Las]] interfaces [[Microsoft .NET|NET]] para levmar están disponibles: vea [http://www.johnlapeyre.com/pdl/index.html PDL :: Fit :: Levmar] o [https://metacpan.org/module/PDL::Fit::LM PDL :: Fit :: LM] , [https://github.com/bjodah/levmar levmar (para python)] , [http://hackage.haskell.org/package/levmar HackageDB levmar] y [https://github.com/AvengerDr/LevmarSharp LevmarSharp] .
* [[R (lenguaje de programación)]] tiene el paquete [https://cran.r-project.org/web/packages/minpack.lm/index.html minpack.lm] (cuadro restringido).

Implementaciones sin restricciones:

* La implementación más antigua aún en uso es [http://www.netlib.org/minpack/ lmdif] , de [[MINPACK]], en [[Fortran]] , en el [[dominio público]] . Ver también:
** [https://archive.today/20130714113627/http://apps.jcns.fz-juelich.de/lmfit lmfit] , una implementación en [[C (lenguaje de programación)|C]] autocontenida del algoritmo MINPACK, con una envoltura fácil de usar para el ajuste de curvas, licencia liberal ([[freeBSD]]).
** [http://eigen.tuxfamily.org/index.php?title=Main_Page eigen] , una biblioteca de [[álgebra lineal]] de C ++, incluye una adaptación del algoritmo minpack en el módulo "NonLinearOptimization".
** La [[GNU Scientific Library|Biblioteca Científica de GNU]] tiene una interfaz en C para MINPACK.
** [http://devernay.free.fr/hacks/cminpack.html C / C ++ Minpack] incluye el algoritmo de Levenberg-Marquardt.
** Python library [[SciPy|scipy]] , módulo <code>scipy.optimize.leastsq</code> proporciona un contenedor para las rutinas de MINPACK .
* [http://www.ics.forth.gr/%7elourakis/sparseLM/ sparseLM] es una implementación de [[C (lenguaje de programación)|C]] destinada a minimizar funciones con grandes jacobianos arbitrariamente [[Matriz dispersa|dispersos]]. Incluye una interfaz de MATLAB MEX. Sólo sin restricciones.
* [https://web.archive.org/web/20130722142233/http://www2.imm.dtu.dk/~hbni/Software/SMarquardt.m SMarquardt.m] es una rutina independiente para Matlab u Octave.
* La biblioteca [http://www.bnikolic.co.uk/inmin/inmin-library.html InMin] contiene una implementación en C ++ del algoritmo basada en la biblioteca de álgebra lineal [http://eigen.tuxfamily.org/index.php?title=Main_Page eigen] C ++. Tiene una API de lenguaje C puro, así como un enlace de Python.
* NMath tiene una implementación para [[Microsoft .NET|NET Framework]] .
* [[gnuplot]] utiliza su propia implementación [http://www.gnuplot.info/ gnuplot.info] .
* Implementaciones de [[Java (lenguaje de programación)|lenguaje de programación Java]]: 1) [http://scribblethink.org/Computer/Javanumeric/index.html Javanumerics] , 2) [https://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math4/fitting/leastsquares/LevenbergMarquardtOptimizer.html Apache Commons Math] {{Wayback|url=https://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math4/fitting/leastsquares/LevenbergMarquardtOptimizer.html |date=20190223040656 }} , 3) [http://finmath.net/finmath-lib/ finmath lib] .
* [http://oooconv.free.fr/fitoo/fitoo_en.html OOoConv] implementa el algoritmo L – M como una hoja de cálculo de OpenOffice.org Calc.
* [https://github.com/namp/lmam-olmam-matlab-toolbox La caja de herramientas LMAM/OLMAM Matlab] implementa Levenberg-Marquardt con un impulso adaptativo para el entrenamiento de redes neuronales avanzadas.
* [https://raullaasner.github.io/gadfit GADfit] es una implementación de adaptación global de Fortran basada en un Levenberg-Marquardt modificado. Utiliza la diferenciación automática. Permite ajustar funciones de complejidad arbitraria, incluyendo integrales.

{{Control de autoridades}}
[[Categoría:Algoritmos estadísticos]]