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.