The Levenberg-Marquardt algorithm

Non linear Least Squares

Nonlinear least squares regression extends linear least squares regression for use with a much larger and more general class of functions.

A nonlinear model is any model of the form $y=f(x;\beta )+ϵ$ our y in this case could take the form $y=\frac{{\beta }_{0}X}{{\beta }_1+X}$

Consider a set of m possibly non linear equations in n variables $x=(x_1,\dots x_n)$ which can be written as $f_i(x)=0,i=1,..,m$. for any x we call $f_i(x)$ which is the ith residual.

we can seek an approximate solution by finding x that minimizes the sum of squares of the residuals

simply put this translates to

the optimal condition to be met is the partial derivatives of $f(x)$ with respect to each $x_1,\dots .,x_n$ must vanish at $\hat{x}$

Gauss–Newton Method

The Gauss-Newton method is a modification of Newton’s method (a root-finding algorithm). Instead of using the second derivative (Hessian matrix), it uses a line search for the search direction.

In Newton’s method, the step size update requires calculating the full Hessian matrix $H$, which involves second-order derivatives of the total error function. This can be complex. The Gauss-Newton method simplifies this by using the Jacobian matrix to approximate the Hessian.

Mathematical Description

Sum of Squared Errors: $E(x,w)=\frac{1}{2}\sum _{m=1}^M{e}_m^2$ Where:

• ( x ) is the input vector
• ( w ) is the weight vector
• ( e ) is the error for the m-th data point defined as: $e_m=y_m-f(x_m,w)$

Gradient of the Error Function: The gradient ( g ) is defined as: $g=\frac{\partial E(x,w)}{\partial w}$

Integrating the error function and the gradient function gives us: $g_i=\frac{\partial E}{\partial w_i}=\frac{\partial \left(\frac{1}{2}\sum _{m=1}^M{e}_m^2\right)}{\partial w_i}$

The relationship between the gradient vector and the Jacobian matrix is: $g=Je$

Hessian Matrix: The Hessian matrix ( H ) is given by: $h_{i,j}=\frac{{\partial }^{2}E}{\partial w_i\partial w_j}=\frac{{\partial }^2\left(\frac{1}{2}\sum _{m=1}^M{e}_m^2\right)}{\partial w_i\partial w_j}$

The relationship between the Hessian and the Jacobian matrix is approximated as: $H\approx J^{T}J$

Gauss-Newton Algorithm: The update rule for the Gauss-Newton algorithm is: $w_{k+1}=w_k-(J_k^T{J}_k{)}^{-1}{J}_k{e}_k$

Pseudo code

1. Initialize weights $w_0$.
2. Compute the Jacobian matrix $J$ and the error vector $e$.
3. Update the weights using: $w_{k+1}=w_k-(J^{T}J{)}^{-1}{J}^{T}e$.
4. Repeat until convergence.

Implementation in F#

For this as well we will do the implementation in F#.

1
2
3
4
5
6
7
8
9
10
11
let gaussNewton (r: Vector<float> -> Vector<float>) (J: Vector<float> -> Matrix<float>) (beta0: Vector<float>) (maxIter: int) (tol: float) =    
    let rec iterate beta iter =
        if iter > maxIter then beta 
        else
            let Jbeta = J beta
            let rBeta = r beta
            let delta = -(Jbeta.Transpose() * Jbeta).Inverse() * (Jbeta.Transpose() * rBeta)
            let newBeta = beta + delta
            if delta.L2Norm() < tol then newBeta 
            else iterate newBeta (iter + 1) 
    iterate beta0 0

The Levenberge-Marquardt Algorithm.

If you take a closer look at the Levenberge-Marquardt method, it combines the advantages of gradient decent (which is stability) and gauss-Newton methods(which is speed). The Levenberge-Marquardt instead of using line search like the Gauss-Newton methods it instead uses a trust region strategy.

Without going further into the details using trust region avoids one of the weaknesses of Gauss–Newton, which is its behavior when the Jacobian $J(x)$ is rank-deficient.How the Levenberg-Marquardt algorithm avoids this is in how the hessian matrix $J^{T}J$ is approximated,it introduces a different approximation unlike the gauss newton method. $\begin{array}{rl}H& \approx J^{T}J+\mu I\\ \text{Where}:\\ \mu & \text{ is positive, combination coefficient}\\ I& \text{ is the Identity Matrix}\end{array}$

this makes the update rule for the Levenberg-Marquardt to be defined as :

Pseudo Code

1. Initialize weights $w_0$, damping parameter $\mu >0$, and tolerance $ϵ$.
2. Repeat until convergence:
   • Step 1: Compute the Jacobian matrix $J_k$ and the error vector $e_k$ at $w_k$.
   • Step 2: Compute the approximate Hessian matrix:
     $H=J_k^T{J}_k+\mu I$
   • Step 3: Compute the update direction:
     $\mathrm{\Delta }{w}_k=-H^{-1}{J}_k^T{e}_k$
   • Step 4: Update the weights:
     $w_{k+1}=w_k+\mathrm{\Delta }{w}_k$
   • Step 5: Evaluate the error function at the new weights $w_{k+1}$:
     • If the error decreases:
       • Accept the new weights $w_{k+1}$.
       • Decrease $\mu$ (e.g., $\mu =\mu /\beta ),where(\beta >1$.
     • Else:
       • Reject $w_{k+1}$ and increase $\mu$ (e.g.,$\mu =\mu \cdot \beta$.
3. Check for convergence:
   • Stop if $|\mathrm{\Delta }{w}_k|<ϵ$ or the error does not improve significantly.

Implementation in F#

We can proceed to do a basic implementation with LU decomposition.For the purpose of understanding.I’d recommend using the Fsharp.stats library for non linear fitting given its robustness.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
let levenbergMarquardt jacobian  func (xData: float[][]) (yData: float[]) (initialParameters: float[]) (maxIter: int) (lambda: float) =
    let mutable parameters = initialParameters
    let mutable lambda = lambda
    let nParams = initialParameters.Length

    for iter in 1 .. maxIter do
        let J = Array.init xData.Length (fun i -> jacobian xData.[i] parameters)
        let r = Array.init xData.Length (fun i -> yData.[i] - func xData.[i] parameters)

        let JT = Array.transpose J
        let JTJ = Array.init nParams (fun i -> 
            Array.init nParams (fun j -> Array.sumBy (fun k -> JT.[i].[k] * J.[k].[j]) [| 0 .. J.Length - 1 |]))
        let JTr = Array.init nParams (fun i -> Array.sumBy (fun k -> JT.[i].[k] * r.[k]) [| 0 .. J.Length - 1 |])

        //we Modify the hessian with levenberg regularization
        let H = Array.init nParams (fun i -> 
            Array.init nParams (fun j -> if i = j then JTJ.[i].[j] + lambda else JTJ.[i].[j]))

        //solve ln eq with LU decomposition
        let deltaW = solveWithLU H JTr

        parameters <- Array.map2 (+) parameters deltaW

        //lambda tuning
        let newError = Array.sumBy (fun i -> (yData.[i] - func xData.[i] parameters) ** 2.0) [| 0 .. xData.Length - 1 |]
        let oldError = Array.sumBy (fun i -> (yData.[i] - func xData.[i] (Array.map2 (-) parameters deltaW)) ** 2.0) [| 0 .. xData.Length - 1 |]

        if newError < oldError then
            lambda <- lambda / 10.0
        else
            lambda <- lambda * 10.0

    parameters

To verify that our implementation works, we fit a quadratic function to sample data points.

1
2
3
4
5
6
7
8
9
10
11
12
13
let modelObjfunc (x: float[]) (parameters: float[]) : float =   
    parameters.[0] * x.[0] ** 2.0 + parameters.[1] * x.[0] + parameters.[2]

let jacobian (x: float[]) (parameters: float[]) =
    [| x.[0] ** 2.0; x.[0]; 1.0 |]  

let xData = [| [| 1.0 |]; [| 2.0 |]; [| 3.0 |]; [| 4.0 |] |]
let yData = [| 2.2; 2.8; 3.6; 4.5 |]
let initialParameters = [| 0.5; 0.5; 0.5 |]
let maxIterations = 500
let lambda = 0.01

let optimizedParams = levenbergMarquardt jacobian  modelObjfunc xData yData initialParameters maxIterations lambda

The results from the above are: ` [|0.075; 0.395; 1.725|]` After running the optimization, we can plot the fitted curve using Plotly.NET

Conclusions.

The Levenberg-Marquardt algorithm is best applied in curve fitting. It blends the best of gradient descent (which is great for slow, steady learning) and Gauss-Newton (which is faster but less stable). By adapting its approach based on how close it is to the best solution, it improves accuracy and speed. However, like many optimization techniques, it usually finds the closest good solution (local minimum) rather than the absolute best one (global minimum).Another challenge is when the residuals are large in converges slower or fails to converge at all.

References

The Levenberg-Marquardt Method and its Implementation in Python