Levenberg–Marquardt algorithm

In mathematics and computing, the Levenberg–Marquardt algorithm (LMA or just LM), also known as the damped least-squares (DLS) method, is used to solve non-linear least squares problems. These minimization problems arise especially in least squares curve fitting.

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 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,^[1] while working at the Frankford Army Arsenal. It was rediscovered in 1963 by Donald Marquardt,^[2] who worked as a statistician at DuPont, and independently by Girard^[3] Wynne^[4] and Morrison.^[5]

The problem

The primary application of the Levenberg–Marquardt algorithm is in the least-squares curve fitting problem: given a set of $m$ empirical datum pairs (x_i, y_i) of independent and dependent variables, find the parameters β of the model curve f(x, β) so that the sum of the squares of the deviations S(β) is minimized:

{\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},

which is assumed to be non-empty.

The solution

Like other numeric minimization algorithms, the Levenberg–Marquardt algorithm is an iterative procedure. To start a minimization, the user has to provide an initial guess for the parameter vector β. In cases with only one minimum, a uninformed standard guess like β^T = (1, 1, ..., 1) will work fine; in cases with multiple minima, the algorithm converges to the global minimum only if the initial guess is already somewhat close to the final solution.

In each iteration step, the parameter vector β is replaced by a new estimate β + δ. To determine δ, the function $f(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }})$ is approximated by its linearization:

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

where

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

is the gradient (row-vector in this case) of f with respect to β.

The sum $S({\boldsymbol {\beta }})$ of square deviations has at its minimum a zero gradient with respect to β. The above first-order approximation of $f(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }})$ gives

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

or in vector notation,

{\begin{aligned}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{aligned}}

Taking the derivative of $S({\boldsymbol {\beta }}+{\boldsymbol {\delta }})$ with respect to δ and setting the result to zero gives

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

where $\mathbf {J}$ is the Jacobian matrix, whose i-th row equals $\mathbf {J} _{i}$ , and where $\mathbf f$ and $\mathbf {y}$ are vectors with i-th component $f(x_{i},{\boldsymbol {\beta }})$ and $y_{i}$ respectively. This is a set of linear equations, which can be solved for δ.

Levenberg's contribution is to replace this equation by a "damped version":

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

where I is the identity matrix, giving as the increment δ to the estimated parameter vector β.

The (non-negative) damping factor λ is adjusted at each iteration. If reduction of 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, λ can be increased, giving a step closer to the gradient-descent direction. Note that the gradient of S with respect to β equals $-2{\big (}\mathbf {J} ^{\mathrm {T} }[\mathbf {y} -\mathbf {f} ({\boldsymbol {\beta }})]{\big )}^{\mathrm {T} }$ . Therefore, for large values of λ, the step will be taken approximately in the direction of the gradient. If either the length of the calculated step δ or the reduction of sum of squares from the latest parameter vector β + δ fall below predefined limits, iteration stops, and the last parameter vector β is considered to be the solution.

Levenberg's algorithm has the disadvantage that if the value of damping factor λ is large, inverting J^TJ + λI is not used at all. Marquardt 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, Marquardt replaced the identity matrix I with the diagonal matrix consisting of the diagonal elements of J^TJ, resulting in the Levenberg–Marquardt algorithm:

[\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 }})].

A similar damping factor appears in Tikhonov regularization, which is used to solve linear ill-posed problems, as well as in ridge regression, an estimation technique in statistics.

Choice of damping parameter

Various more or less heuristic arguments have been put forward for the best choice for the damping parameter λ. Theoretical arguments exist showing why some of these choices guarantee local convergence of the algorithm; however, these choices can make the global convergence of the algorithm suffer from the undesirable properties of steepest descent, in particular, very slow convergence close to the optimum.

The absolute values of any choice depends on how well-scaled the initial problem is. Marquardt recommended starting with a value λ₀ and a factor ν > 1. Initially setting λ = λ₀ and computing the residual sum of squares S(β) after one step from the starting point with the damping factor of λ = λ₀ and secondly with λ₀/ν. If both of these are worse than the initial point, then the damping is increased by successive multiplication by ν until a better point is found with a new damping factor of λ₀ν^k for some k.

If use of the damping factor λ/ν results in a reduction in squared residual, then this is taken as the new value of λ (and the new optimum location is taken as that obtained with this damping factor) and the process continues; if using λ/ν resulted in a worse residual, but using λ resulted in a better residual, then λ is left unchanged and the new optimum is taken as the value obtained with λ as damping factor.

Example

Poor fit

Better fit

Best fit

In this example we try to fit the function $y=a\cos(bX)+b\sin(aX)$ using the Levenberg–Marquardt algorithm implemented in GNU Octave as the leasqr function. The graphs show progressively better fitting for the parameters a = 100, b = 102 used in the initial curve. Only when the parameters in the last graph are chosen closest to the original, are the curves fitting exactly. This equation is an example of very sensitive initial conditions for the Levenberg–Marquardt algorithm. One reason for this sensitivity is the existence of multiple minima — the function $\cos(\beta x)$ has minima at parameter value ${\hat {\beta }}$ and ${\hat {\beta }}+2n\pi$ .

References

↑ Levenberg, Kenneth (1944). "A Method for the Solution of Certain Non-Linear Problems in Least Squares". Quarterly of Applied Mathematics. 2: 164–168.
↑ Marquardt, Donald (1963). "An Algorithm for Least-Squares Estimation of Nonlinear Parameters". SIAM Journal on Applied Mathematics. 11 (2): 431–441. doi:10.1137/0111030.
↑ Girard, André (1958). "Excerpt from Revue d'optique théorique et instrumentale". Rev. Opt. 37: 225–241, 397–424.
↑ Wynne, C. G. (1959). "Lens Designing by Electronic Digital Computer: I". Proc. Phys. Soc. Lond. 73 (5): 777. doi:10.1088/0370-1328/73/5/310.
↑ Morrison, David D. (1960). "Methods for nonlinear least squares problems and convergence proofs". Proceedings of the Jet Propulsion Laboratory Seminar on Tracking Programs and Orbit Determination: 1–9.
↑ 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.

External links

Descriptions

Detailed description of the algorithm can be found in Numerical Recipes in C, Chapter 15.5: Nonlinear models
C. T. Kelley, Iterative Methods for Optimization, SIAM Frontiers in Applied Mathematics, no 18, 1999, ISBN 0-89871-433-8. Online copy
History of the algorithm in SIAM news
A tutorial by Ananth Ranganathan
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, The Levenberg-Marquardt method for nonlinear least squares curve-fitting problems (MATLAB implementation included)

Implementations

Levenberg-Marquardt is a built-in algorithm in SciPy, GNU Octave, Scilab, Mathematica , Matlab, NeuroSolutions, Origin, Fityk, IGOR Pro, LabVIEW and SAS numerical computing environments. There also exist numerous software libraries which allow to use LM algorithm in standalone applications. Some of them support only basic unconstrained optimization, whilst other ones support different combinations of box and linear constraints.

Box and linearly constrained implementations:

ALGLIB has box and linearly constrained implementation of improved LM in C# and C++. Improved algorithm takes less time to converge and can use either Jacobian or exact Hessian.
Artelys Knitro is a non-linear solver with an implementation of the box-constrained Levenberg–Marquardt algorithm. It is written in C and has interfaces to C++/C#/Java/Python/MATLAB/R.
ceres is a non-linear minimization library with an implementation of the box-constrained Levenberg–Marquardt algorithm. It is written in C++ and uses eigen.
IDL, add-on MPFIT supports box constraints.
levmar is an implementation in C/C++ with support for box and general linear constraints, distributed under the GNU General Public License.
- levmar includes a MEX file interface for MATLAB.
- Perl (PDL), python, Haskell and .NET interfaces to levmar are available: see PDL::Fit::Levmar or PDL::Fit::LM, levmar (for python), HackageDB levmar and LevmarSharp.
R (programming language) has the minpack.lm package (box constrained).

Unconstrained implementations:

The oldest implementation still in use is lmdif, from MINPACK, in Fortran, in the public domain. See also:
- lmfit, a self-contained C implementation of the MINPACK algorithm, with an easy-to-use wrapper for curve fitting, liberal licence (freeBSD).
- eigen, a C++ linear-algebra library, includes an adaptation of the minpack algorithm in the "NonLinearOptimization" module.
- The GNU Scientific Library has a C interface to MINPACK.
- C/C++ Minpack includes the Levenberg–Marquardt algorithm.
- Python library scipy, module scipy.optimize.leastsq provides wrapper for the MINPACK routines.
sparseLM is a C implementation aimed at minimizing functions with large, arbitrarily sparse Jacobians. Includes a MATLAB MEX interface. Unconstrained only.
SMarquardt.m is a stand-alone routine for Matlab or Octave.
InMin library contains a C++ implementation of the algorithm based on the eigen C++ linear-algebra library. It has a pure C-language API, as well as a Python binding.
NMath has an implementation for the .NET Framework.
gnuplot uses its own implementation gnuplot.info.
Java programming language implementations: 1) Javanumerics, 2) Apache Commons Math, 3) finmath lib.
OOoConv implements the L–M algorithm as an OpenOffice.org Calc spreadsheet.
LMAM/OLMAM Matlab toolbox implements Levenberg–Marquardt with adaptive momentum for training feedforward neural networks.
GADfit is a Fortran implementation of global fitting based on a modified Levenberg–Marquardt. Uses automatic differentiation. Allows fitting functions of arbitrary complexity, including integrals.

This article is issued from Wikipedia. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.

[Levenberg-1] Levenberg, Kenneth (1944). "A Method for the Solution of Certain Non-Linear Problems in Least Squares". Quarterly of Applied Mathematics. 2: 164–168.

[Marquardt-2] Marquardt, Donald (1963). "An Algorithm for Least-Squares Estimation of Nonlinear Parameters". SIAM Journal on Applied Mathematics. 11 (2): 431–441. doi:10.1137/0111030.

[Girard-3] Girard, André (1958). "Excerpt from Revue d'optique théorique et instrumentale". Rev. Opt. 37: 225–241, 397–424.

[Wynne-4] Wynne, C. G. (1959). "Lens Designing by Electronic Digital Computer: I". Proc. Phys. Soc. Lond. 73 (5): 777. doi:10.1088/0370-1328/73/5/310.

[Morrison-5] Morrison, David D. (1960). "Methods for nonlinear least squares problems and convergence proofs". Proceedings of the Jet Propulsion Laboratory Seminar on Tracking Programs and Orbit Determination: 1–9.

[6] 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.