The Levenberg-Marquardt algorithm for nonlinear least ...

1 The Levenberg-Marquardt algorithm fornonlinear least squares curve-fitting problems Henri P. GavinDepartment of Civil and Environmental EngineeringDuke UniversitySeptember 18, 2020 AbstractThe Levenberg-Marquardt algorithm was developed in the early 1960 s to solvenonlinear least squares problems. least squares problems arise in the context of fit-ting a parameterized mathematical model to a set of data points by minimizing anobjective expressed as the sum of the squares of the errors between the model functionand a set of data points. If a model is linear in its parameters, the least squares ob-jective is quadratic in the parameters. This objective may be minimized with respectto the parameters in one step via the solution to a linear matrix equation.

2 If the fitfunction is not linear in its parameters, the least squares problem requires an itera-tive solution algorithm . Such algorithms reduce the sum of the squares of the errorsbetween the model function and the data points through a sequence of well-chosenupdates to values of the model parameters. The Levenberg-Marquardt algorithm com-bines two numerical minimization algorithms: the gradient descent method and theGauss- newton method. In the gradient descent method, the sum of the squared er-rors is reduced by updating the parameters in the steepest-descent direction. In theGauss- newton method, the sum of the squared errors is reduced by assuming the leastsquares function is locally quadratic in the parameters, and finding the minimum ofthis quadratic.

3 The Levenberg-Marquardt method acts more like a gradient-descentmethod when the parameters are far from their optimal value, and acts more like theGauss- newton method when the parameters are close to their optimal value. This doc-ument describes these methods and illustrates the use of software to solve nonlinearleast squares curve-fitting fitting a model function y(t;p) of an independent variabletand a vector ofnparameterspto a set ofmdata points (ti,yi), it is customary and convenient to minimizethe sum of the weighted squares of the errors (or weighted residuals) between the datayiand the curve-fit function y(t;p). 2(p) =m i=1[y(ti) y(ti;p) yi]2(1)= (y y(p))TW(y y(p))(2)=yTW y 2yTW y+ yTW y(3)where yiis the measurement error for datumy(ti).

4 Typically the weighting matrixWisdiagonal withWii= 1/ 2yi. More formally,Wcan be set to the inverse of the measurement2 The Levenberg-Marquardt algorithm for nonlinear least squareserror covariance matrix, in the unusual case that it is known. More generally, the weightsWii, can be set to pursue other curve-fitting goals. This scalar-valued goodness-of-fit measureis called thechi-squared error criterionbecause the sum of squares of normally-distributedvariables is distributed as the chi-squared the function y(t;p) is nonlinear in the model parametersp, then the minimization of 2with respect to the parameters must be carried out iteratively. The goal of each iterationis to find a perturbationhto the parameterspthat reduces Gradient Descent MethodThe steepest descent method is a general minimization method which updates parame-ter values in the downhill direction: the direction opposite to the gradient of the objectivefunction.

5 The gradient descent method converges well for problems with simple objectivefunctions [6, 7]. For problems with thousands of parameters, gradient descent methods aresometimes the only viable gradient of the chi-squared objective function with respect to the parameters is p 2= 2(y y(p))TW p(y y(p))(4)= 2(y y(p))TW[ y(p) p](5)= 2(y y)TW J(6)where them nJacobian matrix [ y/ p] represents the local sensitivity of the function yto variation in the parametersp. For notational simplicity the variableJwill be usedfor [ y/ p]. Note that in models that are linear in the parameters, y=Xp, the Jacobian[ y/ p] is the matrix of model basis vectorsX. The parameter updatehthat moves theparameters in the direction of steepest descent is given byhgd= JTW(y y),(7)where the positive scalar determines the length of the step in the steepest-descent gauss - newton MethodThe gauss - newton method is a method for minimizing a sum-of-squares objective func-tion.

6 It presumes that the objective function is approximately quadratic in the parametersnear the optimal solution [2]. For moderately-sized problems the gauss - newton methodtypically converges much faster than gradient-descent methods [8].The function evaluated with perturbed model parameters may be locally approximatedthrough a first-order Taylor series expansion. y(p+h) y(p) +[ y p]h= y+Jh,(8) Gavin September 18, Gavin3 Substituting the approximation y(p+h) y(p) +Jhinto equation (3) for 2(p+h), 2(p+h) yTW y+ yTW y 2yTW y 2(y y)TW Jh+hTJTW Jh.(9)The first-order Taylor approximation (8) results in an approximation for 2that is quadraticin the perturbationh. The Hessian of the chi-squared fit criterion is approximatelyJTW parameter updatehthat minimizes 2is found from 2/ h= 0: h 2(p+h) 2(y y)TW J+ 2hTJTW J,(10)and the resulting normal equations for the gauss - newton update are[JTW J]hgn=JTW(y y).

7 (11)Note that the right hand side vectors in normal equations for the gradient descent method(7) and the gauss - newton method (11) are Levenberg-Marquardt MethodThe Levenberg-Marquardt algorithm adaptively varies the parameter updates betweenthe gradient descent update and the gauss - newton update,[JTW J+ I]hlm=JTW(y y),(12)where small values of thedamping parameter result in a gauss - newton update and largevalues of result in a gradient descent update. The damping parameter is initializedto be large so that first updates are small steps in the steepest-descent direction. If anyiteration happens to result in a worse approximation ( 2(p+hlm)> 2(p)), then isincreased. Otherwise, as the solution improves, is decreased, the levenberg -Marquardtmethod approaches the gauss - newton method, and the solution typically accelerates to thelocal minimum [6, 7, 8].

8 In marquardt s update relationship [8], the values of are normalized to the values ofJTW J.[JTW J+ diag(JTW J)]hlm=JTW(y y),(13) ImplementationMany variations of the Levenberg-Marquardt have been published in papers and incode, , [4, 6, 10, 11]. This document borrows from some of these. In iterationi, thestephis evaluated by comparing 2(p) to 2(p+h). The step is accepted if the metric i[9] is greater than a user-specified threshold, 4>0. This metric is a measure of theactual improvement in 2as compared to the improvement of an LM update assuming theapproximation (8) were exact. i(hlm) = 2(p) 2(p+hlm)(y y)TW(y y) (y y Jhlm)TW(y y Jhlm)(14)= 2(p) 2(p+hlm)hTlm( ihlm+JTW(y y(p)))if using eq n (12) forhlm(15)= 2(p) 2(p+hlm)hTlm( idiag(JTWJ)hlm+JTW(y y(p)))if using eq n (13) forhlm(16) Gavin September 18, 20204 The Levenberg-Marquardt algorithm for nonlinear least squaresIf in an iteration i(h)> 4thenp+his sufficiently better thanp,pis replaced byp+h, and is reduced by a factor.

9 Otherwise is increased by a factor, and the algorithm proceedsto the next Initialization and update of the L-M parameter, , and the may select one of three methods for initializing and updating 0= o; ois user-specified [8].use eq n (13) forhlmand eq n (16) for if i(h)> 4:p p+h; i+1= max[ i/L ,10 7];otherwise: i+1= min [ iL ,107];2. 0= omax[diag[JTW J]]; ois eq n (12) forhlmand eq n (15) for =((JTW(y y(p)))Th)/(( 2(p+h) 2(p))/2 + 2(JTW(y y(p)))Th);if i( h)> 4:p p+ h; i+1= max [ i/(1 + ),10 7];otherwise: i+1= i+| 2(p+ h) 2(p)|/(2 );3. 0= omax[diag[JTW J]]; ois user-specified [9].use eq n (12) forhlmand eq n (15) for if i(h)> 4:p p+h; i+1= imax [1/3,1 (2 i 1)3] ; i= 2;otherwise: i+1= i i; i+1= 2 i;For the examples in section , method 1 [8] withL 11 andL 9 exhibits goodconvergence Computation and rank-1 update of the Jacobian,[ y/ p]In the first iteration, in every 2niterations, and in iterations where 2(p+h)> 2(p),the Jacobian (J Rm n) is numerically approximated using forward differences,Jij= yi pj= y(ti;p+ pj) y(ti;p)|| pj||,(17)or central differences (default)Jij= yi pj= y(ti;p+ pj) y(ti;p pj)2|| pj||,(18)where thej-th element of pjis the only non-zero element and is set to j(1 +|pj|).

10 Forproblems with many parameters, a finite differences Jacobian is computationally the Jacobian is re-computed using finite differences only occasionally, convergence can beachieved with fewer function evaluations. So in all other iterations, the Jacobian is updatedusing the Broyden rank-1 update formula,J=J+(( y(p+h) y(p) Jh)hT)/(hTh).(19)This rank-1 Jacobian update equation (19) requires no additional function Gavin September 18, Convergence criteriaConvergence is achieved whenoneof the following three criteria is satisfied, Convergence in the gradient: max JTW(y y) < 1; Convergence in parameters: max|hi/pi|< 2; or Convergence in 2: uses the value of thereduced 2, 2 = 2/(m n+ 1)< , iterations terminate when the iteration count exceeds a pre-specified AnalysisOnce the optimal curve-fit parameterspfitare determined, parameter statistics arecomputed for the converged solution.