Gauss–Newton algorithm

The Gauss–Newton algorithm is used to solve non-linear least squares problems. It is a modification of Newton's method for finding a minimum of a function. Unlike Newton's method, the Gauss–Newton algorithm can only be used to minimize a sum of squared function values, but 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.

Description

Given m functions r = (r₁, …, r_m) (often called residuals) of n variables β = (β₁, …, β_n), with m ≥ n, the Gauss–Newton algorithm iteratively finds the value of the variables which minimizes the sum of squares^[1]

S(\boldsymbol \beta)= \sum_{i=1}^m r_i^2(\boldsymbol \beta).

Starting with an initial guess $\boldsymbol \beta^{(0)}$ for the minimum, the method proceeds by the iterations

\boldsymbol \beta^{(s+1)} = \boldsymbol \beta^{(s)} - \left(\mathbf{J_r}^\mathsf{T} \mathbf{J_r} \right)^{-1} \mathbf{ J_r} ^\mathsf{T} \mathbf{r}(\boldsymbol \beta^{(s)})

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

(\mathbf{J_r})_{ij} = \frac{\partial r_i (\boldsymbol \beta^{(s)})}{\partial \beta_j},

and the symbol $^\mathsf{T}$ denotes the matrix transpose.

If m = n, the iteration simplifies to

\boldsymbol \beta^{(s+1)} = \boldsymbol \beta^{(s)} - \left( \mathbf{J_r} \right)^{-1} \mathbf{r}(\boldsymbol \beta^{(s)})

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

In data fitting, where the goal is to find the parameters β such that a given model function y = f(x, β) best fits some data points (x_i, y_i), the functions r_i are the residuals

r_i(\boldsymbol \beta)= y_i - f(x_i, \boldsymbol \beta).

Then, the Gauss-Newton method can be expressed in terms of the Jacobian J_f of the function f as

\boldsymbol \beta^{(s+1)} = \boldsymbol \beta^{(s)} + \left(\mathbf{J_f}^\mathsf{T} \mathbf{J_f} \right)^{-1} \mathbf{ J_f} ^\mathsf{T}\mathbf{r}(\boldsymbol \beta^{(s)}).

Notes

The assumption m ≥ n in the algorithm statement is necessary, as otherwise the matrix J_r^TJ_r 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 r_i. Using Taylor's theorem, we can write at every iteration:

\mathbf{r}(\boldsymbol \beta)\approx \mathbf{r}(\boldsymbol \beta^s)+\mathbf{J_r}(\boldsymbol \beta^s)\Delta

with $\Delta=\boldsymbol \beta-\boldsymbol \beta^s.$ The task of finding Δ minimizing the sum of squares of the right-hand side, i.e.,

\mathbf{min}\|\mathbf{r}(\boldsymbol \beta^s)+\mathbf{J_r}(\boldsymbol \beta^s)\Delta\|_2^2

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

The normal equations are m linear simultaneous equations in the unknown increments, Δ. They may be solved in one step, using Cholesky decomposition, or, better, the QR factorization of J_r. 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 J_r, the iterations will fail as J_r^TJ_r becomes singular.

Example

Calculated curve obtained with

\hat\beta_1=0.362

and

\hat\beta_2=0.556

(in blue) versus the observed data (in red).

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 $[S]$ and reaction rate in an enzyme-mediated reaction, the data in the following table were obtained.

$i$	1	2	3	4	5	6	7
$[S]$	0.038	0.194	0.425	0.626	1.253	2.500	3.740
rate	0.050	0.127	0.094	0.2122	0.2729	0.2665	0.3317

It is desired to find a curve (model function) of the form

\text{rate}=\frac{V_\text{max}[S]}{K_M+[S]}

that fits best the data in the least squares sense, with the parameters $V_\text{max}$ and $K_M$ to be determined.

Denote by $x_i$ and $y_i$ the value of $[S]$ and the rate from the table, $i=1, \dots, 7.$ Let $\beta_1=V_\text{max}$ and $\beta_2=K_M.$ We will find $\beta_1$ and $\beta_2$ such that the sum of squares of the residuals

r_i = y_i - \frac{\beta_1x_i}{\beta_2+x_i} \; (i=1,\dots, 7)

is minimized.

The Jacobian $\mathbf{J_r}$ of the vector of residuals $r_i$ in respect to the unknowns $\beta_j$ is an $7\times 2$ matrix with the $i$ -th row having the entries

\frac{\partial r_i}{\partial \beta_1}= \frac{x_i}{\beta_2+x_i},\ \frac{\partial r_i}{\partial \beta_2}= -\frac{\beta_1x_i}{\left(\beta_2+x_i\right)^2}.

Starting with the initial estimates of $\beta_1=0.9$ and $\beta_2=0.2$ , after five iterations of the Gauss–Newton algorithm the optimal values $\hat\beta_1=0.362$ and $\hat\beta_2=0.556$ 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 versus the observed data.

Convergence properties

It can be shown^[2] that the increment Δ is a descent direction for S, and, if the algorithm converges, then the limit is a stationary point of S. However, convergence is not guaranteed, not even local convergence as in Newton's method.

The rate of convergence of the Gauss–Newton algorithm can approach quadratic.^[3] The algorithm may converge slowly or not at all if the initial guess is far from the minimum or the matrix $\mathbf{J_r^\mathsf{T} J_r}$ is ill-conditioned. For example, consider the problem with $m=2$ equations and $n=1$ variable, given by

\begin{align} r_1(\beta) &= \beta + 1 \\ r_2(\beta) &= \lambda \beta^2 + \beta - 1. \end{align}

The optimum is at $\beta = 0$ . (Actually the optimum is at $\beta = -1$ for $\lambda = 2$ , because $S(0) = 1^2 + (-1)^2 = 2$ , but $S(-1) = 0$ .) If $\lambda = 0$ 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.^[4]

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.^[5]

The recurrence relation for Newton's method for minimizing a function S of parameters, $\boldsymbol\beta$ , is

\boldsymbol\beta^{(s+1)} = \boldsymbol\beta^{(s)} - \mathbf H^{-1} \mathbf g \,

where g denotes the gradient vector of S and H denotes the Hessian matrix of S. Since $S = \sum_{i=1}^m r_i^2$ , the gradient is given by

g_j=2\sum_{i=1}^m r_i\frac{\partial r_i}{\partial \beta_j}.

Elements of the Hessian are calculated by differentiating the gradient elements, $g_j$ , with respect to $\beta_k$

H_{jk}=2\sum_{i=1}^m \left(\frac{\partial r_i}{\partial \beta_j}\frac{\partial r_i}{\partial \beta_k}+r_i\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k} \right).

The Gauss–Newton method is obtained by ignoring the second-order derivative terms (the second term in this expression). That is, the Hessian is approximated by

H_{jk}\approx 2\sum_{i=1}^m J_{ij}J_{ik}

where $J_{ij}=\frac{\partial r_i}{\partial \beta_j}$ are entries of the Jacobian J_r. The gradient and the approximate Hessian can be written in matrix notation as

\mathbf g=2\mathbf{J}_\mathbf{r}^\mathsf{T}\mathbf{r}, \quad \mathbf{H} \approx 2 \mathbf{J}_\mathbf{r}^\mathsf{T}\mathbf{J_r}.\,

These expressions are substituted into the recurrence relation above to obtain the operational equations

\boldsymbol{\beta}^{(s+1)} = \boldsymbol\beta^{(s)}+\Delta;\quad \Delta = -\left( \mathbf{J_r}^\mathsf{T}\mathbf{J_r} \right)^{-1} \mathbf{J_r}^\mathsf{T}\mathbf{r}.

Convergence of the Gauss–Newton method is not guaranteed in all instances. The approximation

\left|r_i\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k}\right| \ll \left|\frac{\partial r_i}{\partial \beta_j}\frac{\partial r_i}{\partial \beta_k}\right|

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.^[6]

The function values $r_i$ are small in magnitude, at least around the minimum.
The functions are only "mildly" non linear, so that $\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k}$ is relatively small in magnitude.

Improved versions

With the Gauss–Newton method the sum of squares of the residuum S may not decrease at every iteration. However, since Δ is a descent direction, unless $S(\boldsymbol \beta^s)$ is a stationary point, it holds that $S(\boldsymbol \beta^s+\alpha\Delta) < S(\boldsymbol \beta^s)$ for all sufficiently small $\alpha>0$ . Thus, if divergence occurs, one solution is to employ a fraction, $\alpha$ , of the increment vector, Δ in the updating formula

\boldsymbol \beta^{s+1} = \boldsymbol \beta^s+\alpha\ \Delta

In other words, the increment vector is too long, but it points in "downhill", so going just a part of the way will decrease the objective function S. An optimal value for $\alpha$ can be found by using a line search algorithm, that is, the magnitude of $\alpha$ is determined by finding the value that minimizes S, usually using a direct search method in the interval $0<\alpha<1$ .

In cases where the direction of the shift vector is such that the optimal fraction, $\alpha$ , is close to zero, an alternative method for handling divergence is the use of the Levenberg–Marquardt algorithm, also known as the "trust region method".^[1] The normal equations are modified in such a way that the increment vector is rotated towards the direction of steepest descent,

\left(\mathbf{J^TJ+\lambda D}\right)\Delta=-\mathbf{J}^T \mathbf{r}

where D is a positive diagonal matrix. Note that when D is the identity matrix I and $\lambda\to+\infty$ , then $\lambda\Delta=\lambda\left(\mathbf{J^TJ}+\lambda \mathbf{I}\right)^{-1}\left(-\mathbf{J}^T \mathbf{r}\right)=\left(\mathbf{I}-\mathbf{J^TJ}/ \lambda+\cdots\right)\left(-\mathbf{J}^T \mathbf{r}\right)\to -\mathbf{J}^T \mathbf{r}$ , therefore the direction of Δ approaches the direction of the negative gradient $-\mathbf{J}^T \mathbf{r}$ .

The so-called Marquardt parameter, $\lambda$ , may also be optimized by a line search, but this is inefficient as the shift vector must be re-calculated every time $\lambda$ 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

For large scale optimization, the Gauss-Newton method is of special interest because it is often (though certainly not always) true that the matrix $\mathbf{J}_\mathbf{r}$ is more sparse than the approximate Hessian $\mathbf{J}_\mathbf{r}^\mathsf{T}\mathbf{J_r}$ . 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.

In order to make this kind of approach work, one needs at minimum an efficient method for computing the product

\mathbf{J}_\mathbf{r}^\mathsf{T}\mathbf{J_r} \mathbf p

for some vector p. With sparse matrix storage, it is in general practical to store the rows of $\mathbf{J}_\mathbf{r}$ in a compressed form (eg, without zero entries), making a direct computation of the above product tricky due to the transposition. However, if one defines c_i as row i of the matrix $\mathbf{J}_\mathbf{r}$ , the following simple relation holds

\mathbf{J}_\mathbf{r}^\mathsf{T}\mathbf{J_r} \mathbf p = \sum_i \mathbf c_i (\mathbf c_i \cdot \mathbf p)

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 c_i is the gradient of the corresponding residual r_i; with this in mind, the formula above emphasizes the fact that residuals contribute to the problem independently of each other.

Other applications

The Gauss–Newton algorithm is a popular method for solving nonlinear inverse problems. A particular application is generating computational models of oil and gas reservoirs for consistency with observed production data.^[7]

Related algorithms

In a quasi-Newton method, such as that due to Davidon, Fletcher and Powell or Broyden–Fletcher–Goldfarb–Shanno (BFGS method) an estimate of the full Hessian, $\frac{\partial^2 S}{\partial \beta_j \partial\beta_k}$ , is built up numerically using first derivatives $\frac{\partial r_i}{\partial\beta_j}$ 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.

Notes

1 2 Björck (1996)
↑ Björck (1996) p260
↑ Björck (1996) p341, 342
↑ Fletcher (1987) p.113
↑ http://www.henley.ac.uk/web/FILES/maths/09-04.pdf
↑ Nocedal (1999) p259
↑ "History matching production data and uncertainty assessment with an efficient TSVD parameterization algorithm". Journal of Petroleum Science and Engineering 113: 54–71. doi:10.1016/j.petrol.2013.11.025.

References

Björck, A. (1996). Numerical methods for least squares problems. SIAM, Philadelphia. ISBN 0-89871-360-9.
Fletcher, Roger (1987). Practical methods of optimization (2nd ed.). New York: John Wiley & Sons. ISBN 978-0-471-91547-8. .
Nocedal, Jorge; Wright, Stephen (1999). Numerical optimization. New York: Springer. ISBN 0-387-98793-2.

Isaac Newton

Publications	De analysi per aequationes numero terminorum infinitas (1669, published 1711) Method of Fluxions (1671) De motu corporum in gyrum (1684) Philosophiæ Naturalis Principia Mathematica (1687) General Scholium (1713) Opticks (1704) The Queries (1704) Arithmetica Universalis (1707)

Other writings	Notes on the Jewish Temple Quaestiones quaedam philosophicae The Chronology of Ancient Kingdoms Amended (1728) An Historical Account of Two Notable Corruptions of Scripture (1754)

Newtonianism	Bucket argument Newton's inequalities Newton's law of cooling Newton's law of universal gravitation Post-Newtonian expansion Parameterized post-Newtonian formalism Newton–Cartan theory Schrödinger–Newton equation Gravitational constant Newton's laws of motion Newtonian dynamics Newton's method in optimization Gauss–Newton algorithm Truncated Newton method Newton's rings Newton's theorem about ovals Newton–Pepys problem Newtonian potential Newtonian fluid Classical mechanics Newtonian fluid Corpuscular theory of light Leibniz–Newton calculus controversy Rotating spheres Newton's cannonball Newton–Cotes formulas Newton's method Newton fractal Generalized Gauss–Newton method Newton's identities Newton polynomial Newton's theorem of revolving orbits Newton–Euler equations Newton number Kissing number problem Power number Solar mass Dynamics Absolute time and space Finite difference Table of Newtonian series Impact depth Structural coloration Inertia

Life	Cranbury Park Woolsthorpe Manor Early life of Isaac Newton Later life of Isaac Newton Religious views of Isaac Newton Isaac Newton's occult studies The Mysteryes of Nature and Art Scientific revolution Copernican Revolution

Friends and family	Catherine Barton John Conduitt William Clarke Benjamin Pulleyn William Stukeley William Jones Isaac Barrow Abraham de Moivre John Keill

Discoveries and inventions	Calculus Newton disc Newton polygon Newton–Okounkov body Newton's reflector Newtonian telescope Newton scale Newton's metal Newton's cradle Sextant

Phrases	Hypotheses non fingo Standing on the shoulders of giants

Theory expansions	Kepler's laws of planetary motion Problem of Apollonius

Related	Writing of Principia Mathematica Newton (Blake) List of things named after Isaac Newton Isaac Newton in popular culture Elements of the Philosophy of Newton Isaac Newton S/O Philipose Newton (unit)

Least squares and regression analysis

Computational statistics

Correlation and dependence

Regression analysis

Regression as a
statistical model

Linear regression	Simple linear regression Ordinary least squares Generalized least squares Weighted least squares General linear model

Predictor structure	Polynomial regression Growth curve (statistics) Segmented regression Local regression

Non-standard	Nonlinear regression Nonparametric Semiparametric Robust Quantile Isotonic

Non-normal errors	Generalized linear model Binomial Poisson Logistic

Decomposition of variance

Model exploration

Background

Design of experiments

Numerical approximation

Applications

Regression analysis category
Statistics category
Statistics portal
Statistics outline
Statistics topics

Optimization: Algorithms, methods, and heuristics

Unconstrained nonlinear: Methods calling …

… functions

… and gradients

Convergence	Trust region Wolfe conditions

Quasi–Newton	BFGS and L-BFGS DFP Symmetric rank-one (SR1)

Other methods	Gauss–Newton Gradient Levenberg–Marquardt Conjugate gradient Truncated Newton

… and Hessians

Newton's method

The graph of a strictly concave quadratic function is shown in blue, with its unique maximum shown as a red dot. Below the graph appears the contours of the function: The level sets are nested ellipses.

Constrained nonlinear

General	Barrier methods Penalty methods

Differentiable	Augmented Lagrangian methods Sequential quadratic programming Successive linear programming

Convex optimization

Convex
minimization

Linear and
quadratic

Interior point	Affine scaling Ellipsoid algorithm of Khachiyan Projective algorithm of Karmarkar

Basis-Exchange	Simplex algorithm of Dantzig Revised simplex algorithm Criss-cross algorithm Principal pivoting algorithm of Lemke

Combinatorial

Paradigms

Graph
algorithms

Minimum spanning tree	Bellman–Ford Borůvka Dijkstra Floyd–Warshall Johnson Kruskal

Network flows

Metaheuristics

Evolutionary algorithm Hill climbing Local search Simulated annealing Tabu search

Categories
- Algorithms and methods
- Heuristics
Software

This article is issued from Wikipedia - version of the Friday, February 12, 2016. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.