Gauss–Newton algorithm

The Gauss–Newton algorithm is a method used to solve non-linear least squares problems. Non-linear least squares is the form of Least squares analysis which is used to fit a set of m observations with a model that is non-linear in n unknown It can be seen as a modification of Newton's method for finding a minimum of a function. In Mathematics, Newton's method is a well-known Algorithm for finding roots of Equations in one or more Dimensions It can also be used In Mathematics, maxima and minima, known collectively as extrema, are the largest value (maximum or smallest value (minimum that The Mathematical concept of a function expresses dependence between two quantities one of which is given (the independent variable, argument of the 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. In statistics nonlinear regression is a form of Regression analysis in which observational data are modeled by a function which is a nonlinear combination of the model parameters

The method is due to the renowned mathematician Carl Friedrich Gauss. Johann Carl Friedrich Gauss (ˈɡaʊs, Gauß Carolus Fridericus Gauss ( 30 April 1777 – 23 February 1855) was a German

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

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. In Computational mathematics, an iterative method attempts to solve a problem (for example an equation or system of equations by finding successive Approximations In Vector calculus, the Jacobian is shorthand for either the Jacobian matrix or its Determinant, the Jacobian determinant. The superscript T denotes the matrix transpose. This article is about the Matrix Transpose operator For other uses see Transposition In Linear algebra, the transpose of a

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}}.$

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. In Statistics and optimization, the concepts of statistical error and residual are easily confused with each other

The Gauss–Newton algorithm can be derived by linearly approximating the vector of functions $r i . In Mathematics, a linear approximation is an approximation of a general function using a Linear function (more precisely an Affine function$ 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. In Calculus, Taylor's theorem gives a sequence of approximations of a Differentiable function around a given point by Polynomials (the Taylor Linear least squares is an important Computational problem that arises primarily in applications when it is desired to fit a linear Mathematical model to

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. Linear least squares is an important Computational problem that arises primarily in applications when it is desired to fit a linear Mathematical model to In Linear algebra, the QR decomposition (also called the QR factorization) of a matrix is a decomposition of the matrix into an orthogonal For large systems, an iterative method, such as the conjugate gradient method, may be more efficient. In Computational mathematics, an iterative method attempts to solve a problem (for example an equation or system of equations by finding successive Approximations In Mathematics, the conjugate gradient method is an Algorithm for the Numerical solution of particular systems of linear equations, namely those If there is a linear dependence between columns of J_r, the iterations will fail as $\mathbf{J_r^T J_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_{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.

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. In optimization, a descent direction is a vector \mathbf{p}\in\mathbb R^n that in the sense below moves us closer towards a local minimum \mathbf{x}^* In Mathematics, particularly in Calculus, a stationary point is an input to a function where the Derivative is zero (equivalently the However, convergence is not guaranteed, not even local convergence as in Newton's method. In Numerical analysis, an Iterative method is called locally convergent if the successive Approximations produced by the method are guaranteed to converge In Mathematics, Newton's method is a well-known Algorithm for finding roots of Equations in one or more Dimensions It can also be used

The rate of convergence of the Gauss–Newton algorithm can approach quadratic. In Numerical analysis (a branch of Mathematics) the speed at which a convergent Sequence approaches its limit is called the rate of convergence ^[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. In Numerical analysis, the condition number associated with a problem is a measure of that problem's amenability to digital computation that is hownumerically well-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]

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. In Mathematics, Newton's method is a well-known Algorithm for finding roots of Equations in one or more Dimensions It can also be used 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. In Vector calculus, the gradient of a Scalar field is a Vector field which points in the direction of the greatest rate of increase of the scalar In Mathematics, the Hessian matrix is the Square matrix of second-order Partial derivatives of a function. 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.

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 (unconstrained optimization, the line search strategy is one of two basic iterative approaches to finding a local minimum \mathbf{x}^* In (unconstrained optimization, the line search strategy is one of two basic iterative approaches to finding a local minimum \mathbf{x}^*

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". In Mathematics and computing the Levenberg–Marquardt algorithm (or LMA) provides a numerical solution to the problem of minimizing a function generally ^[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. For the analytical method called "steepest descent" see Method of steepest descent. 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}$ . Direction is the information contained in the relative position of one point with respect to another point without the Distance information

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.

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. In optimization, quasi-Newton methods (also known as variable metric methods are well-known algorithms for finding local Maxima and minima of functions The Davidon-Fletcher-Powell formula ( DFP) finds the solution to the secant equation that is closest to the current estimate and satisfies the curvature condition (see below

Another method for solving least squares problems using only first derivatives is gradient descent. For the analytical method called "steepest descent" see Method of steepest descent. However, this method does not take into account the second derivatives even approximately. Consequently, it is highly inefficient for many functions.

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. John Wiley & Sons Inc, also referred to as Wiley, is a global Publishing company that markets its products to professionals and consumers students and instructors 113, ISBN 978-0-471-91547-8 .
^ Nocedal, Jorge; Wright, Stephen (1999). Numerical optimization. New York: Springer. ISBN 0387987932.

© 2009 citizendia.org; parts available under the terms of GNU Free Documentation License, from http://en.wikipedia.org
network: | |

Contents