Levenberg-Marquardt Method

The Levenberg-Marquardt (LM) method is a robust optimization algorithm for solving nonlinear least squares problems. It combines the advantages of gradient descent and the Gauss-Newton method, making it one of the most widely used algorithms for parameter estimation in curve fitting, computer vision, and machine learning.

Problem Formulation

Nonlinear Least Squares Problem

Given a nonlinear model function $f(x;\theta )$ with parameters $\theta \in {\mathbb{R}}^n$ and observed data points $(x_i,y_i)$ for $i=1,2,...,m$, the residuals are:

$r_i(\theta )=y_i-f(x_i;\theta )$

The nonlinear least squares problem aims to find the parameters ${\theta }^{\ast }$ that minimize the sum of squared residuals:

${\min }_{\theta }{\sum }_{i=1}^m{r}_i(\theta {)}^2={\min }_{\theta }|\mathbf{r}(\theta ){|}^2$

where $\mathbf{r}(\theta )=[r_1(\theta ),r_2(\theta ),...,r_m(\theta ){]}^T$ is the residual vector.

The Algorithm

Motivation and Intuition

The Levenberg-Marquardt method addresses the limitations of both:

1. Gradient Descent: Robust far from the solution but slow near the solution
2. Gauss-Newton Method: Fast near the solution but can be unstable far from it

LM combines these approaches by adaptively switching between them based on the progress of the optimization.

Update Equation

The core of the LM algorithm is the damped update equation:

$({\mathbf{J}}^T\mathbf{J}+\lambda \mathbf{D})\Delta \theta =-{\mathbf{J}}^T\mathbf{r}$

where:

• $\mathbf{J}$ is the Jacobian matrix of the residuals
• $\mathbf{r}$ is the residual vector
• $\lambda >0$ is the damping parameter
• $\mathbf{D}$ is a scaling matrix, often the identity matrix $\mathbf{I}$ or the diagonal of ${\mathbf{J}}^T\mathbf{J}$
• $\Delta \theta$ is the parameter update

The new parameter estimate becomes:

${\theta }_{k+1}={\theta }_k+\Delta \theta$

$J=\left[\begin{array}{cccc}\frac{\partial f_1}{\partial x_1}& \frac{\partial f_1}{\partial x_2}& \cdots & \frac{\partial f_1}{\partial x_n}\\ \frac{\partial f_2}{\partial x_1}& \frac{\partial f_2}{\partial x_2}& \cdots & \frac{\partial f_2}{\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial f_m}{\partial x_1}& \frac{\partial f_m}{\partial x_2}& \cdots & \frac{\partial f_m}{\partial x_n}\end{array}\right]$

Extra

$\text{Question 7 [15 marks]}$

Consider using the Gauss-Newton method for nonlinear least squares problems.

(a) Derive the Gauss-Newton algorithm for solving nonlinear least squares problems of the form: ${\min }_{\theta }f(\theta )=\frac{1}{2}{\sum }_{i=1}^m{r}_i(\theta {)}^2$ where $r_i(\theta )$ are the residual functions. Show how this leads to the iterative update formula: ${\theta }_{k+1}={\theta }_k-(J_k^T{J}_k{)}^{-1}{J}_k^T{r}_k$ where $J_k$ is the Jacobian of the residuals evaluated at ${\theta }_k$ and $r_k$ is the vector of residuals. [5]

(b) Consider the nonlinear curve fitting problem: $y={\theta }_1{e}^{{\theta }_{2}x}+{\theta }_3$ Given the following data points: $(x_1,y_1)=(0,5)$, $(x_2,y_2)=(1,3)$, $(x_3,y_3)=(2,2)$, $(x_4,y_4)=(3,1.5)$

Starting with initial parameters ${\theta }_0=[4,-0.5,1{]}^T$, perform one iteration of the Gauss-Newton method. Calculate the residuals, Jacobian, and the parameter update. [7]

(c) Compare the Gauss-Newton method with Newton’s method for optimization. Discuss the advantages and limitations of the Gauss-Newton method, particularly when the residuals are significant at the solution. [3]

$\text{Solution to Question 7 - Gauss-Newton Method}$

(a) For a nonlinear least squares problem: ${\min }_{\theta }f(\theta )=\frac{1}{2}{\sum }_{i=1}^m{r}_i(\theta {)}^2=\frac{1}{2}||r(\theta )|{|}_2^2$

The gradient of $f$ is: $\nabla f(\theta )={\sum }_{i=1}^m{r}_i(\theta )\nabla r_i(\theta )=J(\theta {)}^{T}r(\theta )$

where $J(\theta )$ is the Jacobian matrix with elements $J_{ij}=\frac{\partial r_i}{\partial {\theta }_j}$.

The Hessian of $f$ is: ${\nabla }^{2}f(\theta )=J(\theta {)}^{T}J(\theta )+{\sum }_{i=1}^m{r}_i(\theta ){\nabla }^2{r}_i(\theta )$

The Gauss-Newton method approximates the Hessian by ignoring the second-term (which involves second derivatives of residuals): ${\nabla }^{2}f(\theta )\approx J(\theta {)}^{T}J(\theta )$

Using Newton’s method with this approximation: ${\theta }_{k+1}={\theta }_k-[{\nabla }^{2}f({\theta }_k){]}^{-1}\nabla f({\theta }_k)$ ${\theta }_{k+1}={\theta }_k-[J_k^T{J}_k{]}^{-1}{J}_k^T{r}_k$

This is the Gauss-Newton update formula.

(b) For the model $y={\theta }_1{e}^{{\theta }_{2}x}+{\theta }_3$, the residuals are: $r_i(\theta )=y_i-({\theta }_1{e}^{{\theta }_2{x}_i}+{\theta }_3)$

With initial parameters ${\theta }_0=[4,-0.5,1{]}^T$, let’s calculate the residuals: $r_1=5-(4e^{-0.5\cdot 0}+1)=5-(4+1)=0$ $r_2=3-(4e^{-0.5\cdot 1}+1)=3-(4\cdot 0.607+1)=3-3.428=-0.428$ $r_3=2-(4e^{-0.5\cdot 2}+1)=2-(4\cdot 0.368+1)=2-2.472=-0.472$ $r_4=1.5-(4e^{-0.5\cdot 3}+1)=1.5-(4\cdot 0.223+1)=1.5-1.892=-0.392$

So $r_0=[0,-0.428,-0.472,-0.392{]}^T$

Now, calculate the Jacobian. The partial derivatives are: $\frac{\partial r_i}{\partial {\theta }_1}=-e^{{\theta }_2{x}_i}$ $\frac{\partial r_i}{\partial {\theta }_2}=-{\theta }_1{x}_i{e}^{{\theta }_2{x}_i}$ $\frac{\partial r_i}{\partial {\theta }_3}=-1$

Evaluating at ${\theta }_0=[4,-0.5,1{]}^T$: $J_0=\left[\begin{array}{ccccccccc}-1& 0& -1-0.607& -0.607& -1-0.368& -0.736& -1-0.223& -0.669& -1\end{array}\right]$

Computing $J_0^T{J}_0$: $J_0^T{J}_0=\left[\begin{array}{ccccccc}1.5426& 1.5426& 2.1981.5426& 2.0461& 2.0122.198& 2.012& 4\end{array}\right]$

Computing $J_0^T{r}_0$: $J_0^T{r}_0=\left[\begin{array}{c}0.55070.96241.292\end{array}\right]$

Solving the system $(J_0^T{J}_0)\Delta \theta =J_0^T{r}_0$: $\Delta \theta \approx \left[\begin{array}{c}0.38350.21740.1692\end{array}\right]$

Therefore: ${\theta }_1={\theta }_0-\Delta \theta =\left[\begin{array}{c}4-0.51\end{array}\right]-\left[\begin{array}{c}0.38350.21740.1692\end{array}\right]=\left[\begin{array}{c}3.6165-0.71740.8308\end{array}\right]$

(c) Comparison of Gauss-Newton and Newton’s Method:

Advantages of Gauss-Newton:

• Only requires first derivatives (Jacobian) while Newton’s method needs second derivatives (Hessian)
• More computationally efficient for least squares problems
• Works well when residuals are small at the solution

Limitations:

• Can converge slowly or diverge when residuals are large at the solution
• The approximation of the Hessian (ignoring second-order terms) becomes poor when residuals are significant
• Doesn’t handle rank deficiency in the Jacobian well

When residuals are significant at the solution, the ignored second-order term in the Hessian approximation becomes important. In such cases, modified approaches like the Levenberg-Marquardt algorithm (which adds damping to Gauss-Newton) perform better by improving stability and convergence.