The Gauss–Newton algorithm is used to solve non-linear least squares problems, which is equivalent to minimizing a sum of squared function values. It is an extension of Newton's method for finding a minimum of a non-linear function. Since a sum of squares must be nonnegative, the algorithm can be viewed as using Newton's method to iteratively approximate zeroes of the components of the sum, and thus minimizing the sum. In this sense, the algorithm is also an effective method for solving overdetermined systems of equations. It has the advantage that second derivatives, which can be challenging to compute, are not required.

Non-linear least squares problems arise, for instance, in non-linear regression, where parameters in a model are sought such that the model is in good agreement with available observations.

The method is named after the mathematicians Carl Friedrich Gauss and Isaac Newton, and first appeared in Gauss's 1809 work Theoria motus corporum coelestium in sectionibus conicis solem ambientum.

Description

the Gauss–Newton algorithm iteratively that minimize the sum of squares

where, if r and β are column vectors, the entries of the Jacobian matrix

denotes the matrix transpose.

can be found by rearranging the previous equation in the following two steps:

, which can then be solved in a variety of methods (see Notes).

If , the iteration simplifies to

which is a direct generalization of Newton's method in one dimension.

are the residuals

is the left pseudoinverse.

Notes

The assumption is not invertible and the normal equations cannot be solved (at least uniquely).

The Gauss–Newton algorithm can be derived by linearly approximating the vector of functions ri. Using Taylor's theorem

is a linear least-squares problem, which can be solved explicitly, yielding the normal equations in the algorithm.

The normal equations are n. They may be solved in one step, using Cholesky decomposition, or, better, the QR factorization. For large systems, an iterative method, such as the conjugate gradient method, may be more efficient. If there is a linear dependence between columns of Jr becomes singular.

.

Example

In this example, the Gauss–Newton algorithm will be used to fit a model to some data by minimizing the sum of squares of errors between the data and model's predictions.

In a biology experiment studying the relation between substrate concentration and reaction rate in an enzyme-mediated reaction, the data in the following table were obtained.

to be determined.

the values of and rate

is minimized.

are obtained. The sum of squares of residuals decreased from the initial value of 1.445 to 0.00784 after the fifth iteration. The plot in the figure on the right shows the curve determined by the model for the optimal parameters with the observed data.

Convergence properties

under 4 conditions:.

It can be shown that the increment Δ is a descent direction for , and, if the algorithm converges, then the limit is a stationary point of , however, convergence is not guaranteed, not even local convergence as in Newton's method, or convergence under the usual Wolfe conditions.

The rate of convergence of the Gauss–Newton algorithm can approach quadratic. is ill-conditioned

, then the problem is in fact linear and the method finds the optimum in one iteration. If |λ| < 1, then the method converges linearly and the error decreases asymptotically with a factor |λ| at every iteration. However, if |λ| > 1, then the method does not even converge locally.

Solving overdetermined systems of equations

is an effective method for solving overdetermined systems is the Moore-Penrose inverse (also known as pseudoinverse) of the Jacobian matrix. It can be considered an extension of Newton's method and enjoys the same local quadratic convergence toward isolated regular solutions.

is often called a least squares solution of the overdetermined system.

Derivation from Newton's method

In what follows, the Gauss–Newton algorithm will be derived from Newton's method for function optimization via an approximation. As a consequence, the rate of convergence of the Gauss–Newton algorithm can be quadratic under certain regularity conditions. In general (under weaker conditions), the convergence rate is linear.

The recurrence relation for Newton's method for minimizing a function S

where g denotes the gradient vector of S, and H denotes the Hessian matrix of S.

are entries of the Jacobian Jr

that needs to hold to be able to ignore the second-order derivative terms may be valid in two cases, for which convergence is to be expected:

1. are small in magnitude, at least around the minimum.

2. is relatively small in magnitude.

Improved versions

With the Gauss–Newton method the sum of squares of the residuals S

In other words, the increment vector is too long, but it still points "downhill", so going just a part of the way will decrease the objective function S can be found by using a line search is determined by finding the value that minimizes S, usually using a direct search method or a backtracking line search such as Armijo-line search should be chosen such that it satisfies the Wolfe conditions or the Goldstein conditions.

In cases where the direction of the shift vector is such that the optimal fraction α is close to zero, an alternative method for handling divergence is the use of the Levenberg–Marquardt algorithm, a trust region method. The normal equations are modified in such a way that the increment vector is rotated towards the direction of steepest descent

where D is a positive diagonal matrix. Note that when D is the identity matrix I, therefore the direction.

is changed. A more efficient strategy is this: When divergence occurs, increase the Marquardt parameter until there is a decrease in S. Then retain the value from one iteration to the next, but decrease it if possible until a cut-off value is reached, when the Marquardt parameter can be set to zero; the minimization of S then becomes a standard Gauss–Newton minimization.

Large-scale optimization

is more sparse. In such cases, the step calculation itself will typically need to be done with an approximate iterative method appropriate for large and sparse problems, such as the conjugate gradient method.

for some vector p. With sparse matrix in a compressed form (e.g., without zero entries), making a direct computation of the above product tricky due to the transposition. However, if one defines ci as row i

so that every row contributes additively and independently to the product. In addition to respecting a practical sparse storage structure, this expression is well suited for parallel computations. Note that every row ci is the gradient of the corresponding residual ri; with this in mind, the formula above emphasizes the fact that residuals contribute to the problem independently of each other.

Related algorithms

In a quasi-Newton method, such as that due to Davidon, Fletcher and Powell or Broyden–Fletcher–Goldfarb–Shanno (BFGS method only so that after n refinement cycles the method closely approximates to Newton's method in performance. Note that quasi-Newton methods can minimize general real-valued functions, whereas Gauss–Newton, Levenberg–Marquardt, etc. fits only to nonlinear least-squares problems.

Another method for solving minimization problems using only first derivatives is gradient descent. However, this method does not take into account the second derivatives even approximately. Consequently, it is highly inefficient for many functions, especially if the parameters have strong interactions.

Example implementations

Julia

The following implementation in Julia provides one method which uses a provided Jacobian and another computing with automatic differentiation.

Notes

References

• .

External links

• Probability, Statistics and Estimation The algorithm is detailed and applied to the biology experiment discussed as an example in this article (page 84 with the uncertainties on the estimated values).

Implementations

• Artelys Knitro is a non-linear solver with an implementation of the Gauss–Newton method. It is written in C and has interfaces to C++/C#/Java/Python/MATLAB/R.