Content deleted Content added
No edit summary Tags: Mobile edit Mobile web edit |
→Further reading: rm deadlink that just goes to the same place as the doi |
||
| (34 intermediate revisions by 27 users not shown) | |||
Line 1:
{{short description|Algorithm used to solve non-linear least squares problems}}
▲The LMA is used in many software applications for solving generic curve-fitting problems. However, as with many fitting algorithms, the LMA finds only a [[local minimum]], which is not necessarily the [[global minimum]]. The LMA interpolates between the [[Gauss–Newton algorithm]] (GNA) and the method of [[gradient descent]]. The LMA is more [[Robustness (computer science)|robust]] than the GNA, which means that in many cases it finds a solution even if it starts very far off the final minimum. For well-behaved functions and reasonable starting parameters, the LMA tends to be a bit slower than the GNA. LMA can also be viewed as [[Gauss–Newton]] using a [[trust region]] approach.
The algorithm was first published in 1944 by [[Kenneth Levenberg]],<ref name="Levenberg"/> while working at the [[Frankford Arsenal|Frankford Army Arsenal]]. It was rediscovered in 1963 by [[Donald Marquardt]],<ref name="Marquardt"/> who worked as a [[statistician]] at [[DuPont]], and independently by Girard,<ref name="Girard"/> Wynne<ref name="Wynne"/> and Morrison.<ref name="Morrison"/>
The LMA is used in many software applications for solving generic curve-fitting problems. By using the Gauss–Newton algorithm it often converges faster than first-order methods.<ref>{{cite journal|title=Improved Computation for Levenberg–Marquardt Training|last1=Wiliamowski|first1=Bogdan|last2=Yu|first2=Hao|journal=IEEE Transactions on Neural Networks and Learning Systems|volume=21|issue=6|date=June 2010|url=https://www.eng.auburn.edu/~wilambm/pap/2010/Improved%20Computation%20for%20LM%20Training.pdf}}</ref> However, like other iterative optimization algorithms, the LMA finds only a [[local minimum]], which is not necessarily the [[global minimum]].
== The problem ==
Line 30 ⟶ 31:
&= \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{align}</math>
Taking the derivative of this approximation of <math>S\left (\boldsymbol\beta + \boldsymbol\delta\right )</math> with respect to {{tmath|\boldsymbol\delta}} and setting the result to zero gives
:<math>\left (\mathbf J^{\mathrm T} \mathbf J\right )\boldsymbol\delta = \mathbf J^{\mathrm T}\left [\mathbf y - \mathbf f\left (\boldsymbol\beta\right )\right ],</math>
where <math>\mathbf J</math> is the [[Jacobian matrix and determinant|Jacobian matrix]], whose {{tmath|i}}-th row equals <math>\mathbf J_i</math>, and where <math>\mathbf f\left (\boldsymbol\beta\right )</math> and <math>\mathbf y</math> are vectors with {{tmath|i}}-th component
<math>f\left (x_i, \boldsymbol\beta\right )</math> and <math>y_i</math> respectively. The above expression obtained for {{tmath|\boldsymbol\beta}} comes under the Gauss–Newton method. The Jacobian matrix as defined above is not (in general) a square matrix, but a rectangular matrix of size <math>m \times n</math>, where <math>n</math> is the number of parameters (size of the vector <math>\boldsymbol\beta</math>). The matrix multiplication <math>\left (\mathbf J^{\mathrm T} \mathbf J\right)</math> yields the required <math>n \times n</math> square matrix and the matrix-vector product on the right hand side yields a vector of size <math>n</math>. The result is a set of <math>n</math> linear equations, which can be solved for {{tmath|\boldsymbol\delta}}.
Levenberg's contribution is to replace this equation by a "damped version":
Line 44:
where {{tmath|\mathbf I}} is the identity matrix, giving as the increment {{tmath|\boldsymbol\delta}} to the estimated parameter vector {{tmath|\boldsymbol\beta}}.
The (non-negative) damping factor {{tmath|\lambda}} is adjusted at each iteration. If reduction of {{tmath|S}} is rapid, a smaller value can be used, bringing the algorithm closer to the [[Gauss–Newton algorithm]], whereas if an iteration gives insufficient reduction in the residual, {{tmath|\lambda}} can be increased, giving a step closer to the gradient-descent direction. Note that the [[gradient]] of {{tmath|S}} with respect to {{tmath|\boldsymbol\
When the damping factor {{tmath|\lambda}} is large relative to <math> \| \mathbf J^{\mathrm T} \mathbf J \| </math>, inverting <math> \mathbf J^{\mathrm T} \mathbf J + \lambda \mathbf I </math> is not necessary, as the update is well-approximated by the small gradient step <math> \lambda^{-1} \mathbf J^{\mathrm T}\left [\mathbf y - \mathbf f\left (\boldsymbol\beta\right )\right ]</math>.
Levenberg's algorithm has the disadvantage that if the value of damping factor {{tmath|\lambda}} is large, inverting {{tmath|\mathbf J^\text{T}\mathbf J + \lambda\mathbf I}} is not used at all. R. Fletcher provided the insight that we can scale each component of the gradient according to the curvature, so that there is larger movement along the directions where the gradient is smaller. This avoids slow convergence in the direction of small gradient. Therefore, Fletcher in his 1971 paper ''A modified Marquardt subroutine for non-linear least squares'' replaced the identity matrix {{tmath|\mathbf I}} with the diagonal matrix consisting of the diagonal elements of {{tmath|\mathbf J^\text{T}\mathbf J}}, thus making the solution scale invariant:▼
▲
:<math>\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 ].</math>
Line 58 ⟶ 60:
If use of the damping factor {{tmath|\lambda / \nu}} results in a reduction in squared residual, then this is taken as the new value of {{tmath|\lambda}} (and the new optimum ___location is taken as that obtained with this damping factor) and the process continues; if using {{tmath|\lambda / \nu}} resulted in a worse residual, but using {{tmath|\lambda}} resulted in a better residual, then {{tmath|\lambda}} is left unchanged and the new optimum is taken as the value obtained with {{tmath|\lambda}} as damping factor.
An effective strategy for the control of the damping parameter, called ''delayed gratification'', consists of increasing the parameter by a small amount for each uphill step, and decreasing by a large amount for each downhill step. The idea behind this strategy is to avoid moving downhill too fast in the beginning of optimization, therefore restricting the steps available in future iterations and therefore slowing down convergence.<ref name="Transtrum2011"/> An increase by a factor of 2 and a decrease by a factor of 3 has been shown to be effective in most cases, while for large problems more extreme values can work better, with an increase by a factor of 1.5 and a decrease by a factor of 5.<ref name="Transtrum2012"/>
=== Geodesic acceleration ===
When interpreting the Levenberg–Marquardt step as the velocity <math>\boldsymbol{v}_k</math> along a [[geodesic]] path in the parameter space, it is possible to improve the method by adding a second order term that accounts for the acceleration <math>\boldsymbol{a}_k</math> along the geodesic
:<math>
\boldsymbol{v}_k + \frac{1}{2} \boldsymbol{a}_k
</math>
where <math>\boldsymbol{a}_k</math> is the solution of
:<math>
\boldsymbol{J}_k \boldsymbol{a}_k = -f_{vv} .
</math>
Since this geodesic acceleration term depends only on the [[directional derivative]] <math>f_{vv} = \sum_{\mu\nu} v_{\mu} v_{\nu} \partial_{\mu} \partial_{\nu} f (\boldsymbol{x})</math> along the direction of the velocity <math>\boldsymbol{v}</math>, it does not require computing the full second order derivative matrix, requiring only a small overhead in terms of computing cost.<ref>{{cite web|url=https://www.gnu.org/software/gsl/doc/html/nls.html|title=Nonlinear Least-Squares Fitting|publisher=GNU Scientific Library|archive-url=https://web.archive.org/web/20200414204913/https://www.gnu.org/software/gsl/doc/html/nls.html|archive-date=2020-04-14}}</ref> Since the second order derivative can be a fairly complex expression, it can be convenient to replace it with a [[finite difference]] approximation
:<math>
\begin{align}
f_{vv}^i &\approx \frac{f_i(\boldsymbol{x} + h \boldsymbol{\delta}) - 2 f_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{align}
</math>
where <math>f(\boldsymbol{x})</math> and <math>\boldsymbol{J}</math> have already been computed by the algorithm, therefore requiring only one additional function evaluation to compute <math>f(\boldsymbol{x} + h \boldsymbol{\delta})</math>. The choice of the finite difference step <math>h</math> can affect the stability of the algorithm, and a value of around 0.1 is usually reasonable in general.<ref name="Transtrum2012"/>
Since the acceleration may point in opposite direction to the velocity, to prevent it to stall the method in case the damping is too small, an additional criterion on the acceleration is added in order to accept a step, requiring that
:<math>
\frac{2 \left\| \boldsymbol{a}_k \right\|}{\left\| \boldsymbol{v}_k \right\|} \le \alpha
</math>
where <math>\alpha</math> is usually fixed to a value lesser than 1, with smaller values for harder problems.<ref name="Transtrum2012"/>
The addition of a geodesic acceleration term can allow significant increase in convergence speed and it is especially useful when the algorithm is moving through narrow canyons in the landscape of the objective function, where the allowed steps are smaller and the higher accuracy due to the second order term gives significant improvements.<ref name="Transtrum2012"/>
==Example==
Line 71 ⟶ 110:
== See also ==
* [[Trust region]]
* [[Nelder–Mead method]]
* Variants of the Levenberg–Marquardt algorithm have also been used for solving nonlinear systems of equations.<ref>{{cite journal |doi=10.1016/j.cam.2004.02.013|title=Levenberg–Marquardt methods with strong local convergence properties for solving nonlinear equations with convex constraints|journal=Journal of Computational and Applied Mathematics|volume=172|issue=2|pages=375–397|year=2004|last1=Kanzow|first1=Christian|last2=Yamashita|first2=Nobuo|last3=Fukushima|first3=Masao|bibcode=2004JCoAM.172..375K|doi-access=free}}</ref>
==References==
{{Reflist|refs=
<ref name="Levenberg">{{cite journal
| last=Levenberg |first=Kenneth |
| year = 1944
| title = A Method for the Solution of Certain Non-Linear Problems in Least Squares
Line 83 ⟶ 122:
| volume = 2
|issue=2 | pages = 164–168
|doi=10.1090/qam/10666 | doi-access=free
}}</ref> <ref name="Girard">{{cite journal
| last=Girard |first=André
Line 113 ⟶ 153:
}}</ref>
<ref name="Marquardt">{{cite journal
| last=Marquardt |first=Donald |
| year = 1963
| title = An Algorithm for Least-Squares Estimation of Nonlinear Parameters
Line 121 ⟶ 161:
| issue = 2
| pages = 431–441
|hdl=10338.dmlcz/104299 | hdl-access= free
}}</ref> <ref name="Transtrum2011">{{cite journal|title=Geometry of nonlinear least squares with applications to sloppy models and optimization|year=2011|journal=Physical Review E|volume=83|pages=036701|issue=3|publisher=APS|last1=Transtrum|first1=Mark K|last2=Machta|first2=Benjamin B|last3=Sethna|first3=James P|doi=10.1103/PhysRevE.83.036701|pmid=21517619|arxiv=1010.1449|bibcode=2011PhRvE..83c6701T|s2cid=15361707}}</ref>
<ref name="Transtrum2012">{{cite arXiv|title=Improvements to the Levenberg-Marquardt algorithm for nonlinear least-squares minimization|year=2012|eprint=1201.5885|last1=Transtrum|first1=Mark K|last2=Sethna|first2=James P|class=physics.data-an}}</ref>
}}
Line 140 ⟶ 183:
* {{cite journal
| last1=Gill |first1= Philip E.
| last2=Murray |first2=Walter |
| year = 1978
| title = Algorithms for the solution of the nonlinear least-squares problem
Line 151 ⟶ 194:
}}
* {{cite journal
|last = Pujol
|year = 2007
|publisher = SEG
| ▲ |bibcode= 2007Geop...72W...1P
▲ }}
* {{cite book
| last1 = Nocedal | first1 = Jorge
Line 175 ⟶ 218:
== External links ==
* Detailed description of the algorithm can be found in [
* 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]
* K. Madsen, H. B. Nielsen, O. Tingleff, ''[http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf Methods for Non-Linear Least Squares Problems]''
* 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
{{Optimization algorithms|unconstrained}}
{{DEFAULTSORT:Levenberg-Marquardt algorithm}}
| |||