Menu
Home
People
Places
Arts
History
Plants & Animals
Science
Life & Culture
Technology
Reference.org
Levenberg–Marquardt algorithm
open-in-new
<h2 id="the-problem">The problem</h2> <p>The primary application of the Levenberg–Marquardt algorithm is in the least-squares curve fitting problem: given a set of m {\displaystyle m} empirical pairs ( x i , y i ) {\displaystyle \left(x_{i},y_{i}\right)} of independent and dependent variables, find the parameters β {\displaystyle {\boldsymbol {\beta }}} of the model curve f ( x , β ) {\displaystyle f\left(x,{\boldsymbol {\beta }}\right)} so that the sum of the squares of the deviations S ( β ) {\displaystyle S\left({\boldsymbol {\beta }}\right)} is minimized: </p> β ^ ∈ argmin β S ( β ) ≡ argmin β ∑ i = 1 m [ y i − f ( x i , β ) ] 2 , {\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},} which is assumed to be non-empty. <h2 id="the-solution">The solution</h2> <p>Like other numeric minimization algorithms, the Levenberg–Marquardt algorithm is an <a href="/facts/Iteration/AMc1p7BL">iterative</a> procedure. To start a minimization, the user has to provide an initial guess for the parameter vector β {\displaystyle {\boldsymbol {\beta }}} . In cases with only one minimum, an uninformed standard guess like β T = ( 1 , 1 , … , 1 ) {\displaystyle {\boldsymbol {\beta }}^{\text{T}}={\begin{pmatrix}1,\ 1,\ \dots ,\ 1\end{pmatrix}}} will work fine; in cases with <a href="/facts/Local_minimum/ADlgnyzV">multiple minima</a>, the algorithm converges to the global minimum only if the initial guess is already somewhat close to the final solution. </p><p>In each iteration step, the parameter vector β {\displaystyle {\boldsymbol {\beta }}} is replaced by a new estimate β + δ {\displaystyle {\boldsymbol {\beta }}+{\boldsymbol {\delta }}} . To determine δ {\displaystyle {\boldsymbol {\delta }}} , the function f ( x i , β + δ ) {\displaystyle f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)} is approximated by its <a href="/facts/Gradient/TsR7uukj">linearization</a>: </p> f ( x i , β + δ ) ≈ f ( x i , β ) + J i δ , {\displaystyle f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)\approx f\left(x_{i},{\boldsymbol {\beta }}\right)+\mathbf {J} _{i}{\boldsymbol {\delta }},} <p>where </p> J i = ∂ f ( x i , β ) ∂ β {\displaystyle \mathbf {J} _{i}={\frac {\partial f\left(x_{i},{\boldsymbol {\beta }}\right)}{\partial {\boldsymbol {\beta }}}}} <p>is the <a href="/facts/Gradient/TsR7uukj">gradient</a> (row-vector in this case) of f {\displaystyle f} with respect to β {\displaystyle {\boldsymbol {\beta }}} . </p><p>The sum S ( β ) {\displaystyle S\left({\boldsymbol {\beta }}\right)} of square deviations has its minimum at a zero <a href="/facts/Gradient/TsR7uukj">gradient</a> with respect to β {\displaystyle {\boldsymbol {\beta }}} . The above first-order approximation of f ( x i , β + δ ) {\displaystyle f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)} gives </p> S ( β + δ ) ≈ ∑ i = 1 m [ y i − f ( x i , β ) − J i δ ] 2 , {\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},} <p>or in vector notation, </p> S ( β + δ ) ≈ ‖ y − f ( β ) − J δ ‖ 2 = [ y − f ( β ) − J δ ] T [ y − f ( β ) − J δ ] = [ y − f ( β ) ] T [ y − f ( β ) ] − [ y − f ( β ) ] T J δ − ( J δ ) T [ y − f ( β ) ] + δ T J T J δ = [ y − f ( β ) ] T [ y − f ( β ) ] − 2 [ y − f ( β ) ] T J δ + δ T J T J δ . {\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}}} <p>Taking the derivative of this approximation of S ( β + δ ) {\displaystyle S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)} with respect to δ {\displaystyle {\boldsymbol {\delta }}} and setting the result to zero gives </p> ( J T J ) δ = J T [ y − f ( β ) ] , {\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],} <p>where J {\displaystyle \mathbf {J} } is the <a href="/facts/Jacobian_matrix_and_determinant/0tNv6o8j">Jacobian matrix</a>, whose i {\displaystyle i} -th row equals J i {\displaystyle \mathbf {J} _{i}} , and where f ( β ) {\displaystyle \mathbf {f} \left({\boldsymbol {\beta }}\right)} and y {\displaystyle \mathbf {y} } are vectors with i {\displaystyle i} -th component f ( x i , β ) {\displaystyle f\left(x_{i},{\boldsymbol {\beta }}\right)} and y i {\displaystyle y_{i}} respectively. The above expression obtained for β {\displaystyle {\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 m × n {\displaystyle m\times n} , where n {\displaystyle n} is the number of parameters (size of the vector β {\displaystyle {\boldsymbol {\beta }}} ). The matrix multiplication ( J T J ) {\displaystyle \left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} \right)} yields the required n × n {\displaystyle n\times n} square matrix and the matrix-vector product on the right hand side yields a vector of size n {\displaystyle n} . The result is a set of n {\displaystyle n} linear equations, which can be solved for δ {\displaystyle {\boldsymbol {\delta }}} . </p><p>Levenberg's contribution is to replace this equation by a "damped version": </p> ( J T J + λ I ) δ = J T [ y − f ( β ) ] , {\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],} <p>where I {\displaystyle \mathbf {I} } is the identity matrix, giving as the increment δ {\displaystyle {\boldsymbol {\delta }}} to the estimated parameter vector β {\displaystyle {\boldsymbol {\beta }}} . </p><p>The (non-negative) damping factor λ {\displaystyle \lambda } is adjusted at each iteration. If reduction of S {\displaystyle S} is rapid, a smaller value can be used, bringing the algorithm closer to the <a href="/facts/Gauss%E2%80%93Newton_algorithm/7c5en5cU">Gauss–Newton algorithm</a>, whereas if an iteration gives insufficient reduction in the residual, λ {\displaystyle \lambda } can be increased, giving a step closer to the gradient-descent direction. Note that the <a href="/facts/Gradient/TsR7uukj">gradient</a> of S {\displaystyle S} with respect to β {\displaystyle {\boldsymbol {\beta }}} equals − 2 ( J T [ y − f ( β ) ] ) T {\displaystyle -2\left(\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]\right)^{\mathrm {T} }} . Therefore, for large values of λ {\displaystyle \lambda } , the step will be taken approximately in the direction opposite to the gradient. If either the length of the calculated step δ {\displaystyle {\boldsymbol {\delta }}} or the reduction of sum of squares from the latest parameter vector β + δ {\displaystyle {\boldsymbol {\beta }}+{\boldsymbol {\delta }}} fall below predefined limits, iteration stops, and the last parameter vector β {\displaystyle {\boldsymbol {\beta }}} is considered to be the solution. </p><p>When the damping factor λ {\displaystyle \lambda } is large relative to ‖ J T J ‖ {\displaystyle \|\mathbf {J} ^{\mathrm {T} }\mathbf {J} \|} , inverting J T J + λ I {\displaystyle \mathbf {J} ^{\mathrm {T} }\mathbf {J} +\lambda \mathbf {I} } is not necessary, as the update is well-approximated by the small gradient step λ − 1 J T [ y − f ( β ) ] {\displaystyle \lambda ^{-1}\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]} . </p><p>To make the solution scale invariant Marquardt's algorithm solved a modified problem with each component of the gradient scaled according to the curvature. This provides larger movement along the directions where the gradient is smaller, which avoids slow convergence in the direction of small gradient. Fletcher in his 1971 paper <i>A modified Marquardt subroutine for non-linear least squares</i> simplified the form, replacing the identity matrix I {\displaystyle \mathbf {I} } with the diagonal matrix consisting of the diagonal elements of J T J {\displaystyle \mathbf {J} ^{\text{T}}\mathbf {J} } : </p> [ J T J + λ diag ( J T J ) ] δ = J T [ y − f ( β ) ] . {\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].} <p>A similar damping factor appears in <a href="/facts/Tikhonov_regularization/eq6kH8PX">Tikhonov regularization</a>, which is used to solve linear <a href="/facts/Ill-posed_problems/U7J4Vde8">ill-posed problems</a>, as well as in <a href="/facts/Ridge_regression/Le6KJgJ1">ridge regression</a>, an <a href="/facts/Estimation_theory/4nbgFDec">estimation</a> technique in <a href="/facts/Statistics/DNPsKGYU">statistics</a>. </p> <h3>Choice of damping parameter</h3> <p>Various more or less heuristic arguments have been put forward for the best choice for the damping parameter λ {\displaystyle \lambda } . 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 <a href="/facts/Gradient_descent/pFFrek0F">steepest descent</a>, in particular, very slow convergence close to the optimum. </p><p>The absolute values of any choice depend on how well-scaled the initial problem is. Marquardt recommended starting with a value λ 0 {\displaystyle \lambda _{0}} and a factor ν > 1 {\displaystyle \nu >1} . Initially setting λ = λ 0 {\displaystyle \lambda =\lambda _{0}} and computing the residual sum of squares S ( β ) {\displaystyle S\left({\boldsymbol {\beta }}\right)} after one step from the starting point with the damping factor of λ = λ 0 {\displaystyle \lambda =\lambda _{0}} and secondly with λ 0 / ν {\displaystyle \lambda _{0}/\nu } . If both of these are worse than the initial point, then the damping is increased by successive multiplication by ν {\displaystyle \nu } until a better point is found with a new damping factor of λ 0 ν k {\displaystyle \lambda _{0}\nu ^{k}} for some k {\displaystyle k} . </p><p>If use of the damping factor λ / ν {\displaystyle \lambda /\nu } results in a reduction in squared residual, then this is taken as the new value of λ {\displaystyle \lambda } (and the new optimum location is taken as that obtained with this damping factor) and the process continues; if using λ / ν {\displaystyle \lambda /\nu } resulted in a worse residual, but using λ {\displaystyle \lambda } resulted in a better residual, then λ {\displaystyle \lambda } is left unchanged and the new optimum is taken as the value obtained with λ {\displaystyle \lambda } as damping factor. </p><p>An effective strategy for the control of the damping parameter, called <i>delayed gratification</i>, 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.<a class="footnote-ref" id="fnref:7" href="#fn:7"><sup>7</sup></a> 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.<a class="footnote-ref" id="fnref:8" href="#fn:8"><sup>8</sup></a> </p> <h3>Geodesic acceleration</h3> <p>When interpreting the Levenberg–Marquardt step as the velocity v k {\displaystyle {\boldsymbol {v}}_{k}} along a <a href="/facts/Geodesic/fbWR6l80">geodesic</a> path in the parameter space, it is possible to improve the method by adding a second order term that accounts for the acceleration a k {\displaystyle {\boldsymbol {a}}_{k}} along the geodesic </p> v k + 1 2 a k {\displaystyle {\boldsymbol {v}}_{k}+{\frac {1}{2}}{\boldsymbol {a}}_{k}} <p>where a k {\displaystyle {\boldsymbol {a}}_{k}} is the solution of </p> J k a k = − f v v . {\displaystyle {\boldsymbol {J}}_{k}{\boldsymbol {a}}_{k}=-f_{vv}.} <p>Since this geodesic acceleration term depends only on the <a href="/facts/Directional_derivative/Fjc3JSVn">directional derivative</a> f v v = ∑ μ ν v μ v ν ∂ μ ∂ ν f ( x ) {\displaystyle f_{vv}=\sum _{\mu \nu }v_{\mu }v_{\nu }\partial _{\mu }\partial _{\nu }f({\boldsymbol {x}})} along the direction of the velocity v {\displaystyle {\boldsymbol {v}}} , it does not require computing the full second order derivative matrix, requiring only a small overhead in terms of computing cost.<a class="footnote-ref" id="fnref:9" href="#fn:9"><sup>9</sup></a> Since the second order derivative can be a fairly complex expression, it can be convenient to replace it with a <a href="/facts/Finite_difference/aY1txFhj">finite difference</a> approximation </p> f v v i ≈ f i ( x + h δ ) − 2 f i ( x ) + f i ( x − h δ ) h 2 = 2 h ( f i ( x + h δ ) − f i ( x ) h − J i δ ) {\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}}} <p>where f ( x ) {\displaystyle f({\boldsymbol {x}})} and J {\displaystyle {\boldsymbol {J}}} have already been computed by the algorithm, therefore requiring only one additional function evaluation to compute f ( x + h δ ) {\displaystyle f({\boldsymbol {x}}+h{\boldsymbol {\delta }})} . The choice of the finite difference step h {\displaystyle h} can affect the stability of the algorithm, and a value of around 0.1 is usually reasonable in general.<a class="footnote-ref" id="fnref:10" href="#fn:10"><sup>10</sup></a> </p><p>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 </p> 2 ‖ a k ‖ ‖ v k ‖ ≤ α {\displaystyle {\frac {2\left\|{\boldsymbol {a}}_{k}\right\|}{\left\|{\boldsymbol {v}}_{k}\right\|}}\leq \alpha } <p>where α {\displaystyle \alpha } is usually fixed to a value lesser than 1, with smaller values for harder problems.<a class="footnote-ref" id="fnref:11" href="#fn:11"><sup>11</sup></a> </p><p>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.<a class="footnote-ref" id="fnref:12" href="#fn:12"><sup>12</sup></a> </p> <h2 id="example">Example</h2> <p>In this example we try to fit the function y = a cos ( b X ) + b sin ( a X ) {\displaystyle y=a\cos \left(bX\right)+b\sin \left(aX\right)} using the Levenberg–Marquardt algorithm implemented in <a href="/facts/GNU_Octave/vVVq5evM">GNU Octave</a> as the <i>leasqr</i> function. The graphs show progressively better fitting for the parameters a = 100 {\displaystyle a=100} , b = 102 {\displaystyle 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 ( β x ) {\displaystyle \cos \left(\beta x\right)} has minima at parameter value β ^ {\displaystyle {\hat {\beta }}} and β ^ + 2 n π {\displaystyle {\hat {\beta }}+2n\pi } . </p> <h2 id="see-also">See also</h2> <ul><li><a href="/facts/Trust_region/9xmnBqxY">Trust region</a></li> <li><a href="/facts/Nelder%E2%80%93Mead_method/6aXswEax">Nelder–Mead method</a></li> <li>Variants of the Levenberg–Marquardt algorithm have also been used for solving nonlinear systems of equations.<a class="footnote-ref" id="fnref:13" href="#fn:13"><sup>13</sup></a></li></ul> <h2 id="further-reading">Further reading</h2> <ul><li>Moré, Jorge J.; Sorensen, Daniel C. (1983). <a href="https://digital.library.unt.edu/ark:/67531/metadc283525/m2/1/high_res_d/metadc283525.pdf">"Computing a Trust-Region Step"</a> (PDF). <i>SIAM J. Sci. Stat. Comput</i>. 4 (3): 553–572. <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1137%2F0904038">10.1137/0904038</a>.</li> <li>Gill, Philip E.; Murray, Walter (1978). "Algorithms for the solution of the nonlinear least-squares problem". <i><a href="/facts/SIAM_Journal_on_Numerical_Analysis/uUEWPahC">SIAM Journal on Numerical Analysis</a></i>. 15 (5): 977–992. <a href="/facts/Bibcode_(identifier)/9HtdQSGB">Bibcode</a>:<a href="https://ui.adsabs.harvard.edu/abs/1978SJNA...15..977G">1978SJNA...15..977G</a>. <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1137%2F0715063">10.1137/0715063</a>.</li> <li>Pujol, Jose (2007). "The solution of nonlinear inverse problems and the Levenberg-Marquardt method". <i>Geophysics</i>. 72 (4). SEG: W1 – W16. <a href="/facts/Bibcode_(identifier)/9HtdQSGB">Bibcode</a>:<a href="https://ui.adsabs.harvard.edu/abs/2007Geop...72W...1P">2007Geop...72W...1P</a>. <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1190%2F1.2732552">10.1190/1.2732552</a>.</li> <li>Nocedal, Jorge; Wright, Stephen J. (2006). <i>Numerical Optimization</i> (2nd ed.). Springer. <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 978-0-387-30303-1.</li></ul> <h2 id="external-links">External links</h2> <ul><li>Detailed description of the algorithm can be found in <a href="https://numerical.recipes/book.html">Numerical Recipes in C, Chapter 15.5: Nonlinear models</a></li> <li>C. T. Kelley, <i>Iterative Methods for Optimization</i>, SIAM Frontiers in Applied Mathematics, no 18, 1999, <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 0-89871-433-8. <a href="http://www.siam.org/books/textbooks/fr18_book.pdf">Online copy</a></li> <li><a href="https://web.archive.org/web/20140301154319/http://www3.villanova.edu/maple/misc/mtc1093.html">History of the algorithm in SIAM news</a></li> <li><a href="http://ananth.in/docs/lmtut.pdf">A tutorial by Ananth Ranganathan</a></li> <li>K. Madsen, H. B. Nielsen, O. Tingleff, <i><a href="http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf">Methods for Non-Linear Least Squares Problems</a></i> (nonlinear least-squares tutorial; L-M code: <a href="https://archive.today/20180516200006/http://www2.imm.dtu.dk/projects/hbn_software/marquardt.m">analytic Jacobian</a> <a href="https://archive.today/20180516200045/http://www2.imm.dtu.dk/projects/hbn_software/SMarquardt.m">secant</a>)</li> <li>T. Strutz: <i>Data Fitting and Uncertainty (A practical introduction to weighted least squares and beyond).</i> 2nd edition, Springer Vieweg, 2016, <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 978-3-658-11455-8.</li> <li>H. P. Gavin, <a href="http://people.duke.edu/~hpgavin/ce281/lm.pdf"><i>The Levenberg-Marquardt method for nonlinear least-squares curve-fitting problems</i></a> (<a href="/facts/MATLAB/qPjLISCk">MATLAB</a> implementation included)</li></ul> <h2 id="references">References</h2> <ol> <li id="fn:1"><p>Levenberg, Kenneth (1944). "A Method for the Solution of Certain Non-Linear Problems in Least Squares". Quarterly of Applied Mathematics. 2 (2): 164–168. doi:10.1090/qam/10666. <a href="/wiki/Kenneth_Levenberg" target="_blank">/wiki/Kenneth_Levenberg</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></p></li> <li id="fn:2"><p>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. hdl:10338.dmlcz/104299. <a href="/wiki/Donald_Marquardt" target="_blank">/wiki/Donald_Marquardt</a> <a href="#fnref:2" class="footnote-back-ref">↩</a></p></li> <li id="fn:3"><p>Girard, André (1958). "Excerpt from Revue d'optique théorique et instrumentale". Rev. Opt. 37: 225–241, 397–424. <a href="#fnref:3" class="footnote-back-ref">↩</a></p></li> <li id="fn:4"><p>Wynne, C. G. (1959). "Lens Designing by Electronic Digital Computer: I". Proc. Phys. Soc. Lond. 73 (5): 777–787. Bibcode:1959PPS....73..777W. doi:10.1088/0370-1328/73/5/310. <a href="/wiki/Bibcode_(identifier)" target="_blank">/wiki/Bibcode_(identifier)</a> <a href="#fnref:4" class="footnote-back-ref">↩</a></p></li> <li id="fn:5"><p>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. <a href="#fnref:5" class="footnote-back-ref">↩</a></p></li> <li id="fn:6"><p>Wiliamowski, Bogdan; Yu, Hao (June 2010). "Improved Computation for Levenberg–Marquardt Training" (PDF). IEEE Transactions on Neural Networks and Learning Systems. 21 (6). <a href="https://www.eng.auburn.edu/~wilambm/pap/2010/Improved%20Computation%20for%20LM%20Training.pdf" target="_blank">https://www.eng.auburn.edu/~wilambm/pap/2010/Improved%20Computation%20for%20LM%20Training.pdf</a> <a href="#fnref:6" class="footnote-back-ref">↩</a></p></li> <li id="fn:7"><p>Transtrum, Mark K; Machta, Benjamin B; Sethna, James P (2011). "Geometry of nonlinear least squares with applications to sloppy models and optimization". Physical Review E. 83 (3). APS: 036701. arXiv:1010.1449. Bibcode:2011PhRvE..83c6701T. doi:10.1103/PhysRevE.83.036701. PMID 21517619. S2CID 15361707. <a href="/wiki/ArXiv_(identifier)" target="_blank">/wiki/ArXiv_(identifier)</a> <a href="#fnref:7" class="footnote-back-ref">↩</a></p></li> <li id="fn:8"><p>Transtrum, Mark K; Sethna, James P (2012). "Improvements to the Levenberg-Marquardt algorithm for nonlinear least-squares minimization". arXiv:1201.5885 [physics.data-an]. <a href="/wiki/ArXiv_(identifier)" target="_blank">/wiki/ArXiv_(identifier)</a> <a href="#fnref:8" class="footnote-back-ref">↩</a></p></li> <li id="fn:9"><p>"Nonlinear Least-Squares Fitting". GNU Scientific Library. Archived from the original on 2020-04-14. <a href="https://web.archive.org/web/20200414204913/https://www.gnu.org/software/gsl/doc/html/nls.html" target="_blank">https://web.archive.org/web/20200414204913/https://www.gnu.org/software/gsl/doc/html/nls.html</a> <a href="#fnref:9" class="footnote-back-ref">↩</a></p></li> <li id="fn:10"><p>Transtrum, Mark K; Sethna, James P (2012). "Improvements to the Levenberg-Marquardt algorithm for nonlinear least-squares minimization". arXiv:1201.5885 [physics.data-an]. <a href="/wiki/ArXiv_(identifier)" target="_blank">/wiki/ArXiv_(identifier)</a> <a href="#fnref:10" class="footnote-back-ref">↩</a></p></li> <li id="fn:11"><p>Transtrum, Mark K; Sethna, James P (2012). "Improvements to the Levenberg-Marquardt algorithm for nonlinear least-squares minimization". arXiv:1201.5885 [physics.data-an]. <a href="/wiki/ArXiv_(identifier)" target="_blank">/wiki/ArXiv_(identifier)</a> <a href="#fnref:11" class="footnote-back-ref">↩</a></p></li> <li id="fn:12"><p>Transtrum, Mark K; Sethna, James P (2012). "Improvements to the Levenberg-Marquardt algorithm for nonlinear least-squares minimization". arXiv:1201.5885 [physics.data-an]. <a href="/wiki/ArXiv_(identifier)" target="_blank">/wiki/ArXiv_(identifier)</a> <a href="#fnref:12" class="footnote-back-ref">↩</a></p></li> <li id="fn:13"><p>Kanzow, Christian; Yamashita, Nobuo; Fukushima, Masao (2004). "Levenberg–Marquardt methods with strong local convergence properties for solving nonlinear equations with convex constraints". Journal of Computational and Applied Mathematics. 172 (2): 375–397. Bibcode:2004JCoAM.172..375K. doi:10.1016/j.cam.2004.02.013. <a href="https://doi.org/10.1016%2Fj.cam.2004.02.013" target="_blank">https://doi.org/10.1016%2Fj.cam.2004.02.013</a> <a href="#fnref:13" class="footnote-back-ref">↩</a></p></li> </ol>