Gauss-Newton Method

The Gauss-Newton method is a specialized optimization algorithm primarily used for solving nonlinear least squares problems. It’s particularly effective for parameter estimation in curve fitting and model fitting applications, where we need to minimize the sum of squared residuals between observed data and model predictions.

Problem Formulation

Nonlinear Least Squares

Consider a nonlinear function $f(x;\theta )$ with parameters $\theta \in {\mathbb{R}}^n$, and a set of $m$ observed data points $(x_i,y_i)$ for $i=1,2,\dots ,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 ),\dots ,r_m(\theta ){]}^T$ is the residual vector.

Algorithm Description

The Gauss-Newton method is an iterative procedure that approximates the nonlinear least squares problem with a linear least squares problem at each iteration.

Linearization of the Residual

At each iteration $k$, we linearize the residual function around the current parameter estimate ${\theta }_k$:

$r_i(\theta )\approx r_i({\theta }_k)+\nabla r_i({\theta }_k{)}^T(\theta -{\theta }_k)$

where $\nabla r_i({\theta }_k)$ is the gradient of the $i$-th residual with respect to $\theta$, evaluated at ${\theta }_k$.

Jacobian Matrix

The gradients of all residuals form the Jacobian matrix $\mathbf{J}(\theta )\in {\mathbb{R}}^{m\times n}$:

$\mathbf{J}(\theta )=\left[\begin{array}{ccccccccccccc}\frac{\partial r_1}{\partial {\theta }_1}& \frac{\partial r_1}{\partial {\theta }_2}& \cdots & \frac{\partial r_1}{\partial {\theta }_n}\frac{\partial r_2}{\partial {\theta }_1}& \frac{\partial r_2}{\partial {\theta }_2}& \cdots & \frac{\partial r_2}{\partial {\theta }_n}⋮& ⋮& \ddots & ⋮\frac{\partial r_m}{\partial {\theta }_1}& \frac{\partial r_m}{\partial {\theta }_2}& \cdots & \frac{\partial r_m}{\partial {\theta }_n}\end{array}\right]$

For models of the form $f(x;\theta )$, the Jacobian components are:

$J_{ij}(\theta )=\frac{\partial r_i}{\partial {\theta }_j}=-\frac{\partial f(x_i;\theta )}{\partial {\theta }_j}$

Update Rule

At each iteration, the Gauss-Newton update is:

$\Delta \theta =-({\mathbf{J}}^T\mathbf{J}{)}^{-1}{\mathbf{J}}^T\mathbf{r}$

where both $\mathbf{J}$ and $\mathbf{r}$ are evaluated at the current parameter estimate ${\theta }_k$.

The new parameter estimate becomes:

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

Algorithm Steps

1. Choose an initial parameter estimate ${\theta }_0$
2. For iteration $k=0,1,2,\dots$ until convergence:
   • Compute the residual vector $\mathbf{r}({\theta }_k)$
   • Compute the Jacobian matrix $\mathbf{J}({\theta }_k)$
   • Solve the normal equations: $({\mathbf{J}}^T\mathbf{J})\Delta \theta =-{\mathbf{J}}^T\mathbf{r}$
   • Update: ${\theta }_{k+1}={\theta }_k+\Delta \theta$
   • Check convergence criteria
3. Return the final parameter estimate ${\theta }^{\ast }$

Theoretical Justification

Relationship to Newton’s Method

The Gauss-Newton method can be derived as an approximation to Newton’s method for minimizing the least squares objective function:

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

The gradient of $L$ is:

$\nabla L(\theta )=\mathbf{J}(\theta {)}^T\mathbf{r}(\theta )$

The Hessian of $L$ is:

${\nabla }^{2}L(\theta )=\mathbf{J}(\theta {)}^T\mathbf{J}(\theta )+{\sum }_{i=1}^m{r}_i(\theta ){\nabla }^2{r}_i(\theta )$

The Gauss-Newton method approximates the Hessian by dropping the second term:

${\nabla }^{2}L(\theta )\approx \mathbf{J}(\theta {)}^T\mathbf{J}(\theta )$

This approximation is accurate when:

• The residuals $r_i(\theta )$ are small, or
• The functions $r_i(\theta )$ are nearly linear, making ${\nabla }^2{r}_i(\theta )$ small

Interpretation as Linear Least Squares

Each Gauss-Newton iteration solves a linear least squares problem:

${\min }_{\Delta \theta }|\mathbf{J}({\theta }_k)\Delta \theta +\mathbf{r}({\theta }_k){|}^2$

This is obtained by linearizing the residuals and substituting into the original objective function.

Convergence Properties

Convergence Rate

• When the starting point is sufficiently close to the solution and the residuals are small or zero at the solution, Gauss-Newton exhibits quadratic convergence (similar to Newton’s method)
• With larger residuals, convergence is linear
• May not converge if the starting point is far from the solution or if the problem is ill-conditioned

Comparison with Other Methods

Method	Convergence Rate	Residual Size	Second Derivatives	Matrix Inversion
Gauss-Newton	Quadratic/Linear	Small/Large	Not needed	Required
Newton	Quadratic	Any	Required	Required
Gradient Descent	Linear	Any	Not needed	Not required
Levenberg-Marquardt	Varies	Any	Not needed	Required

Limitations and Challenges

Singular or Ill-Conditioned Jacobian

If ${\mathbf{J}}^T\mathbf{J}$ is singular or nearly singular, the update equation becomes numerically unstable. This can occur due to:

• Overparameterization
• Highly correlated parameters
• Lack of sensitivity to certain parameters

Solutions include:

• Regularization (adding a small value to the diagonal of ${\mathbf{J}}^T\mathbf{J}$)
• Using pseudoinverse instead of inverse
• Singular value decomposition (SVD)

Divergence

Gauss-Newton may diverge when:

• The initial guess is far from the solution
• The objective function has significant nonlinearity
• The model is poorly specified relative to the data

Large Residuals

When residuals are large at the solution, the Hessian approximation becomes poor, leading to:

• Slower convergence
• Potential oscillation or divergence
• Suboptimal solutions

Enhancements and Variants

Line Search

Incorporating a line search to determine the step size $\alpha$:

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

where $\alpha$ is chosen to ensure sufficient decrease in the objective function.

Trust Region Methods

Restrict the step size to a region where the linear approximation is trusted:

${\min }_{\Delta \theta }|\mathbf{J}({\theta }_k)\Delta \theta +\mathbf{r}({\theta }_k){|}^2\text{subject to}|\Delta \theta |\le {\Delta }_k$

where ${\Delta }_k$ is the trust region radius adjusted based on model accuracy.

Levenberg-Marquardt Algorithm

A hybrid of Gauss-Newton and gradient descent:

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

where $\lambda$ is a damping parameter that:

• When large, makes the algorithm behave like gradient descent (robust but slow)
• When small, makes the algorithm behave like Gauss-Newton (fast but potentially unstable)

Practical Implementation

Computing the Jacobian

Three main approaches:

1. Analytical Jacobian: Derive and implement explicit formulas for derivatives
   • Most accurate
   • Potentially complex for complicated models
   • Error-prone for large models
2. Numerical Approximation: Use finite differences
   • Forward difference: $\frac{\partial r_i}{\partial {\theta }_j}\approx \frac{r_i(\theta +h{e}_j)-r_i(\theta )}{h}$
   • Central difference: $\frac{\partial r_i}{\partial {\theta }_j}\approx \frac{r_i(\theta +h{e}_j)-r_i(\theta -h{e}_j)}{2h}$
   • Simple to implement but less accurate and more computationally expensive
3. Automatic Differentiation: Leverage specialized software
   • Combines accuracy of analytical with convenience of numerical methods
   • Modern libraries like TensorFlow, PyTorch, JAX support this

Solving the Normal Equations

Methods to solve $({\mathbf{J}}^T\mathbf{J})\Delta \theta =-{\mathbf{J}}^T\mathbf{r}$:

1. Cholesky Decomposition: If ${\mathbf{J}}^T\mathbf{J}$ is positive definite
   • Decompose ${\mathbf{J}}^T\mathbf{J}=\mathbf{L}{\mathbf{L}}^T$
   • Solve $\mathbf{Ly}=-{\mathbf{J}}^T\mathbf{r}$ for $\mathbf{y}$
   • Solve ${\mathbf{L}}^T\Delta \theta =\mathbf{y}$ for $\Delta \theta$
2. QR Decomposition: More stable but more expensive
   • Decompose $\mathbf{J}=\mathbf{QR}$
   • Solve $\mathbf{R}\Delta \theta =-{\mathbf{Q}}^T\mathbf{r}$
3. SVD (Singular Value Decomposition): Most stable but most expensive
   • Decompose $\mathbf{J}=\mathbf{U\Sigma }{\mathbf{V}}^T$
   • Compute $\Delta \theta =-\mathbf{V}{\mathbf{\Sigma }}^{+}{\mathbf{U}}^T\mathbf{r}$

Convergence Criteria

Common stopping conditions:

1. Small parameter change: $|\Delta \theta |<{\epsilon }_1$
2. Small residual change: $|L({\theta }_{k+1})-L({\theta }_k)|<{\epsilon }_2$
3. Small gradient: $|{\mathbf{J}}^T\mathbf{r}|<{\epsilon }_3$
4. Maximum iterations reached: $k>k_{\max }$

Applications

Curve Fitting

Fitting parameterized curves to data points:

• Exponential decay: $f(t;A,\lambda )=A{e}^{-\lambda t}$
• Power laws: $f(x;a,b)=a{x}^b$
• Gaussian peaks: $f(x;\mu ,\sigma ,A)=A{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}$

Model Calibration

Adjusting model parameters to match experimental data:

• Chemical reaction kinetics
• Mechanical system parameters
• Electrical circuit models

System Identification

Determining mathematical models of dynamic systems:

• Transfer function identification
• State-space model parameter estimation
• Time-series modeling

Implementation Example

Python implementation of the Gauss-Newton method for curve fitting:

import numpy as np
from scipy import optimize
 
def gauss_newton(func, x_data, y_data, theta0, jac=None, tol=1e-6, max_iter=100):
    """
    Gauss-Newton algorithm for nonlinear least squares.
    
    Parameters:
    - func: Function that computes the model prediction
    - x_data: Independent variable data
    - y_data: Dependent variable data (observations)
    - theta0: Initial parameter estimate
    - jac: Function to compute the Jacobian (if None, use finite differences)
    - tol: Convergence tolerance
    - max_iter: Maximum number of iterations
    
    Returns:
    - theta: Optimal parameter estimate
    - info: Information about the optimization
    """
    theta = np.asarray(theta0)
    n_params = len(theta)
    n_samples = len(x_data)
    
    # If no Jacobian function is provided, create a numerical approximation
    if jac is None:
        def numerical_jacobian(theta):
            J = np.zeros((n_samples, n_params))
            epsilon = 1e-8  # Small step for finite differences
            
            for j in range(n_params):
                theta_plus = theta.copy()
                theta_plus[j] += epsilon
                
                # Forward difference approximation
                J[:, j] = -(func(x_data, theta_plus) - func(x_data, theta)) / epsilon
            
            return J
        
        jacobian_func = numerical_jacobian
    else:
        jacobian_func = jac
    
    # Gauss-Newton iterations
    for iteration in range(max_iter):
        # Compute residuals
        y_pred = func(x_data, theta)
        residuals = y_data - y_pred
        
        # Compute objective function value
        obj_value = 0.5 * np.sum(residuals**2)
        
        # Compute Jacobian
        J = jacobian_func(theta)
        
        # Compute gradient of objective function
        gradient = -J.T @ residuals
        
        # Compute Gauss-Newton step
        try:
            # Using normal equations: (J^T J) delta_theta = -J^T residuals
            JTJ = J.T @ J
            delta_theta = np.linalg.solve(JTJ, -J.T @ residuals)
        except np.linalg.LinAlgError:
            # If matrix is singular, use pseudoinverse
            delta_theta = -np.linalg.pinv(J) @ residuals
        
        # Update parameters
        theta_new = theta + delta_theta
        
        # Check convergence
        if np.linalg.norm(delta_theta) < tol:
            theta = theta_new
            break
        
        theta = theta_new
    
    # Compute final residuals and objective value
    y_pred = func(x_data, theta)
    residuals = y_data - y_pred
    final_obj_value = 0.5 * np.sum(residuals**2)
    
    info = {
        'iterations': iteration + 1,
        'obj_value': final_obj_value,
        'residuals': residuals,
        'jacobian': jacobian_func(theta)
    }
    
    return theta, info

Related Pages

• Least Squares
• Model Fitting and Regression
• Levenberg-Marquardt Method
• Trust Region Methods
• Newton-Raphson Method
• Numerical Methods