8.2 The Levenberg-Marquardt Method

One of the methods that is currently implemented in FEBio's optimization method is the constrained Levenberg-Marquardt method via the levmar library. (see http://users.ics.forth.gr/~lourakis/levmar/ for more information on this library.)

The Levenberg-Marquardt method is a numerical algorithm that minimizes a function that is defined as a sum of squares of nonlinear functions, i.e. the objective function as defined above. It combines the steepest-descent method with a Gauss-Newton method to find the parameters that minimize the objective function.

The LM method requires a set of measured values \((x_{i},y_{i})\) and an initial guess for the a vector. It then tries to find a better estimate for a by replacing it with \(\mathbf{a}+\delta\). The function \(f(x_{i};\mathbf{a}+\delta)\) is linearly approximated.

\[ f\left(x_{i};\mathbf{a}+\delta\right)\approx f\left(x_{i};\mathbf{a}\right)+\mathbf{J}_{i}\delta \]

where \(\mathbf{J}_{i}\)is the Jacobian of \(f\) with respect to \(\delta\). Substituting this in the objective function and minimizing with respect to \(\delta\) leads to,

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

where y is the vector of \(y_{i},\)and f is the vector of \(f(x_{i};\boldsymbol{a})\).

The main idea of the LM method is to replace this linear equation with the following.

\[ (\mathbf{J^{\mathrm{T}}J\mathrm{+\mu(\mathbf{J}^{\mathrm{T}}\mathbf{J})_{\mathit{ii}})\delta=}J^{\mathrm{T}}\mathrm{(\mathbf{y}-\mathbf{f}(\boldsymbol{a}))}} \]

Here, \(\mu\)is a damping parameter that is controlled by the algorithm. When \(\mu\)is small, the method approximates Gauss-Newton, when \(\mu\)is large it is closer to a steepdest-descent method. The algorithm will modify try to modify \(\mu\)such that an improvement to the parameter vector a can be found in each iteration. The method will terminate when the value of the objective function falls below a user-specified tolerance (the obj_tol parameter in FEBio).

The evaluation of the Jacobian requires evaluating the derivatives of \(f\) with respect to a. These derivatives are approximated via forward difference formulas. For example, the \(k\)-th component of the gradient is approximated as follows.

\[ \frac{\partial f}{\partial a_{k}}\approx\frac{1}{\delta a_{k}}\left[f\left(a_{1},\cdots,a_{k}+\delta a_{k},\cdots,a_{m}\right)-f\left(a_{1},\cdots,a_{k},\cdots,a_{m}\right)\right] \]

The value for \(\delta a_{k}\)is determined from the following formula.

\[ \delta a_{k}=\varepsilon\left(1+a_{k}\right) \]

where, \(\varepsilon\) is the forward difference scale factor (the fdiff_scale _option in FEBio). In FEBio, the initial value for the damping parameter \(\mu\)can be set with the tau _parameter.