Gauss–Newton algorithm

From Wikipedia, the free encyclopedia

The Gauss–Newton algorithm is a method used to solve non-linear least squares problems. It can be seen as 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 due to the renowned mathematician Carl Friedrich Gauss.

1 The algorithm
2 Notes
3 Example
4 Convergence properties
5 Derivation from Newton's method
6 Improved versions
7 Related algorithms
8 References and notes

[edit] The algorithm

Given m functions $r i$ ( $i=1,\ldots,m$ ) of n variables ${\boldsymbol \beta}=(\beta_1, \beta_2, \dots, \beta_n),$ with m≥n, the Gauss–Newton algorithm finds the minimum of the sum of squares

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

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

$\boldsymbol \beta^{s+1} = \boldsymbol \beta^s+\delta\boldsymbol\beta,$

with the increment $\delta\boldsymbol\beta$ satisfying the normal equations

$\left(\mathbf{J_r^T J_r} \right)\delta\boldsymbol\beta= - \mathbf{ J_r^T r}.$

Here, r is the vector of functions r_i, and J_r is the m×n Jacobian matrix of r with respect to β, both evaluated at β^s. The superscript T denotes the matrix transpose.

In data fitting, where the goal is to find the parameters $\boldsymbol \beta$ such that a given model function $y=f(x, \boldsymbol \beta)$ fits best 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 increment $\delta\boldsymbol\beta$ can be expressed in terms of the Jacobian of the function f, as

$\left( \mathbf{ J_f^T J_f} \right)\delta\boldsymbol\beta= {\mathbf{J_f^T r}}.$

[edit] Notes

The assumption m≥n in the algorithm statement is necessary, as otherwise the matrix $\mathbf{J_r^T J_r}$ is not invertible and the normal equations cannot be solved.

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\boldsymbol\beta$

with $\delta\boldsymbol \beta=\boldsymbol \beta-\boldsymbol \beta^s.$ The task of finding $\delta\boldsymbol \beta$ minimizing the sum of squares of the right-hand side 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, $\delta \boldsymbol\beta$ . They may be solved in one step, using Choleski factorization, 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 $\mathbf{J_r^T J_r}$ becomes singular.

[edit] 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_{max}[S]}{K_M+[S]}$

that fits best the data in the least squares sense, with the parameters $V m a x$ 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 $β 1 = V m a x$ and $β 2 = K M .$ We will find $β 1$ and $β 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 $β 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 $β 1$ =0.9 and $β 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.202 to 0.0886 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.

[edit] Convergence properties

It can be shown that the increment $δβ$ is a descent direction for S ^[2], 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^T J_r}$ is ill-conditioned.

The Gauss–Newton algorithm may fail to converge. 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 $β = 0$ . If $λ = 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]

[edit] 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 is at most quadratic.

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 $β 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 $\mathbf{J_r}$ . The gradient and the approximate Hessian can be written in matrix notation as

$\mathbf g=2\mathbf{J_r^Tr, \quad H \approx 2J_r^TJ_r}.\,$

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

$\boldsymbol{\beta}^{s+1} = \boldsymbol\beta^s+\delta\boldsymbol\beta;\ \delta\boldsymbol\beta= - \mathbf{\left( J_r^T J_r \right)^{-1} J_r^T 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.^[5]

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.

[edit] Improved versions

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

$\boldsymbol \beta^{s+1} = \boldsymbol \beta^s+\alpha\ \delta\boldsymbol\beta$ .

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 $α$ can be found by using a line search algorithm, that is, the magnitude of $α$ is determined by finding the the value that minimizes S, usually using a direct search method in the interval $0 < α < 1$ .

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, 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\boldsymbol\beta=\mathbf{J}^T \mathbf{r}$ ,

where D is a positive diagonal matrix. Note that when D is the identity matrix and $\lambda\to+\infty$ , then $\delta\boldsymbol\beta/\lambda\to \mathbf{J}^T \mathbf{r}$ , therefore the direction of $\delta\boldsymbol\beta$ approaches the direction of the gradient $\mathbf{J}^T \mathbf{r}$ .

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

[edit] Related algorithms

In a quasi-Newton method, such as that due to Davidon, Fletcher and Powell 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.

Another method for solving least squares 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.

[edit] References and notes

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

v • d • e Least squares and regression analysis

Least squares	Linear least squares - Non-linear least squares - Partial least squares -Total least squares - Gauss–Newton algorithm - Levenberg–Marquardt algorithm

Regression analysis	Linear regression - Nonlinear regression - Linear model - Generalized linear model - Robust regression - Least-squares estimation of linear regression coefficients- Mean and predicted response - Poisson regression - Logistic regression - Isotonic regression - Ridge regression - Segmented regression - Nonparametric regression - Regression discontinuity

Statistics	Gauss–Markov theorem - Errors and residuals in statistics - Goodness of fit - Studentized residual - Mean squared error - R-factor (crystallography) - Mean squared prediction error - Minimum mean-square error - Root mean square deviation - Squared deviations - M-estimator

Applications	Curve fitting - Calibration curve - Numerical smoothing and differentiation - Least mean squares filter - Recursive least squares filter - Moving least squares - BHHH algorithm

Categories: Optimization algorithms

Gauss–Newton algorithm

From Wikipedia, the free encyclopedia

Contents

[edit] The algorithm

[edit] Notes

[edit] Example

[edit] Convergence properties

[edit] Derivation from Newton's method

[edit] Improved versions

[edit] Related algorithms

[edit] References and notes

Views

Navigation

Interaction

Search

Languages