TheGauss–Newton algorithm is used to solvenon-linear least squares problems, which is equivalent to minimizing a sum of squared function values. It is an extension ofNewton's method for finding aminimum of a non-linearfunction. Since a sum of squares must be nonnegative, the algorithm can be viewed as using Newton's method to iteratively approximatezeroes of the components of the sum, and thus minimizing the sum. In this sense, the algorithm is also an effective method forsolving overdetermined systems of equations. It has the advantage that second derivatives, which can be challenging to compute, are not required.^[1]

Non-linear least squares problems arise, for instance, innon-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 mathematiciansCarl Friedrich Gauss andIsaac Newton, and first appeared in Gauss's 1809 workTheoria motus corporum coelestium in sectionibus conicis solem ambientum.^[2]

Given${\displaystyle m}$ functions${\displaystyle \mathbf{\text{r}}=(r_1,\dots ,r_m)}$ (often called residuals) of${\displaystyle n}$ variables${\displaystyle \mathit{\beta }=({\beta }_1,\dots {\beta }_n),}$ with${\displaystyle m\ge n,}$ the Gauss–Newton algorithmiteratively finds the value of${\displaystyle \beta }$ that minimize the sum of squares^[3]$${\displaystyle S(\mathit{\beta })=\sum _{i=1}^m{r}_i(\mathit{\beta }{)}^2.}$$

Starting with an initial guess${\displaystyle {\mathit{\beta }}^{(0)}}$ for the minimum, the method proceeds by the iterations$${\displaystyle {\mathit{\beta }}^{(s+1)}={\mathit{\beta }}^{(s)}-{\left({{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{r}}\right)}^{-1}{{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}\mathbf{r}\left({\mathit{\beta }}^{(s)}\right),}$$

where, ifr andβ arecolumn vectors, the entries of theJacobian matrix are$${\displaystyle {\left({\mathbf{J}}_{\mathbf{r}}\right)}_{ij}=\frac{\mathrm{\partial }{r}_i\left({\mathit{\beta }}^{(s)}\right)}{\mathrm{\partial }{\beta }_j},}$$

and the symbol${\displaystyle {}^{\mathrm{T}}}$ denotes thematrix transpose.

At each iteration, the update${\displaystyle \mathrm{\Delta }={\mathit{\beta }}^{(s+1)}-{\mathit{\beta }}^{(s)}}$ can be found by rearranging the previous equation in the following two steps:

• ${\displaystyle \mathrm{\Delta }=-{\left({{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{r}}\right)}^{-1}{{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}\mathbf{r}\left({\mathit{\beta }}^{(s)}\right)}$
• ${\displaystyle {{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{r}}\mathrm{\Delta }=-{{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}\mathbf{r}\left({\mathit{\beta }}^{(s)}\right)}$

With substitutions$A={{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{r}}$,${\displaystyle \mathbf{b}=-{{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}\mathbf{r}\left({\mathit{\beta }}^{(s)}\right)}$, and${\displaystyle \mathbf{x}=\mathrm{\Delta }}$, this turns into the conventional matrix equation of form${\displaystyle A\mathbf{x}=\mathbf{b}}$, which can then be solved in a variety of methods (seeNotes).

Ifm =n, the iteration simplifies to

$${\displaystyle {\mathit{\beta }}^{(s+1)}={\mathit{\beta }}^{(s)}-{\left({\mathbf{J}}_{\mathbf{r}}\right)}^{-1}\mathbf{r}\left({\mathit{\beta }}^{(s)}\right),}$$

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

In data fitting, where the goal is to find the parameters${\displaystyle \mathit{\beta }}$ such that a given model function${\displaystyle \mathbf{f}(\mathbf{x},\mathit{\beta })}$ best fits some data points${\displaystyle (x_i,y_i)}$, the functions${\displaystyle r_i}$are theresiduals:$${\displaystyle r_i(\mathit{\beta })=y_i-f\left(x_i,\mathit{\beta }\right).}$$

Then, the Gauss–Newton method can be expressed in terms of the Jacobian${\displaystyle {\mathbf{J}}_{\mathbf{f}}=-{\mathbf{J}}_{\mathbf{r}}}$ of the function${\displaystyle \mathbf{f}}$ as$${\displaystyle {\mathit{\beta }}^{(s+1)}={\mathit{\beta }}^{(s)}+{\left({{\mathbf{J}}_{\mathbf{f}}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{f}}\right)}^{-1}{{\mathbf{J}}_{\mathbf{f}}}^{\mathrm{T}}\mathbf{r}\left({\mathit{\beta }}^{(s)}\right).}$$

Note that${\displaystyle {\left({{\mathbf{J}}_{\mathbf{f}}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{f}}\right)}^{-1}{{\mathbf{J}}_{\mathbf{f}}}^{\mathrm{T}}}$ is the leftpseudoinverse of${\displaystyle {\mathbf{J}}_{\mathbf{f}}}$.

The assumptionm ≥n in the algorithm statement is necessary, as otherwise the matrix${\displaystyle {{\mathbf{J}}_{\mathbf{r}}}^T{\mathbf{J}}_{\mathbf{r}}}$ is not invertible and the normal equations cannot be solved (at least uniquely).

The Gauss–Newton algorithm can be derived bylinearly approximating the vector of functionsr_i. UsingTaylor's theorem, we can write at every iteration:$${\displaystyle \mathbf{r}(\mathit{\beta })\approx \mathbf{r}\left({\mathit{\beta }}^{(s)}\right)+{\mathbf{J}}_{\mathbf{r}}\left({\mathit{\beta }}^{(s)}\right)\mathrm{\Delta }}$$

with${\displaystyle \mathrm{\Delta }=\mathit{\beta }-{\mathit{\beta }}^{(s)}}$. The task of finding${\displaystyle \mathrm{\Delta }}$ minimizing the sum of squares of the right-hand side; i.e.,$${\displaystyle min{\Vert \mathbf{r}\left({\mathit{\beta }}^{(s)}\right)+{\mathbf{J}}_{\mathbf{r}}\left({\mathit{\beta }}^{(s)}\right)\mathrm{\Delta }\Vert }_2^2,}$$

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

The normal equations aren simultaneous linear equations in the unknown increments${\displaystyle \mathrm{\Delta }}$. They may be solved in one step, usingCholesky decomposition, or, better, theQR factorization of${\displaystyle {\mathbf{J}}_{\mathbf{r}}}$. For large systems, aniterative method, such as theconjugate gradient method, may be more efficient. If there is a linear dependence between columns ofJ_r, the iterations will fail, as${\displaystyle {{\mathbf{J}}_{\mathbf{r}}}^T{\mathbf{J}}_{\mathbf{r}}}$ becomes singular.

When${\displaystyle \mathbf{r}}$ is complex${\displaystyle \mathbf{r}:{\mathbb{C}}^n\to \mathbb{C}}$ the conjugate form should be used:${\displaystyle {\left({\overline{{\mathbf{J}}_{\mathbf{r}}}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{r}}\right)}^{-1}{\overline{{\mathbf{J}}_{\mathbf{r}}}}^{\mathrm{T}}}$.

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.

It is desired to find a curve (model function) of the form$${\displaystyle \text{rate}=\frac{V_{\text{max}}\cdot [S]}{K_M+[S]}}$$

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

Denote by${\displaystyle x_i}$ and${\displaystyle y_i}$ the values of[S] andrate respectively, with${\displaystyle i=1,\dots ,7}$. Let${\displaystyle {\beta }_1=V_{\text{max}}}$ and${\displaystyle {\beta }_2=K_M}$. We will find${\displaystyle {\beta }_1}$ and${\displaystyle {\beta }_2}$ such that the sum of squares of the residuals$${\displaystyle r_i=y_i-\frac{{\beta }_1{x}_i}{{\beta }_2+x_i},(i=1,\dots ,7)}$$

is minimized.

The Jacobian${\displaystyle {\mathbf{J}}_{\mathbf{r}}}$ of the vector of residuals${\displaystyle r_i}$ with respect to the unknowns${\displaystyle {\beta }_j}$ is a${\displaystyle 7\times 2}$ matrix with the${\displaystyle i}$-th row having the entries$${\displaystyle \frac{\mathrm{\partial }{r}_i}{\mathrm{\partial }{\beta }_1}=-\frac{x_i}{{\beta }_2+x_i};\frac{\mathrm{\partial }{r}_i}{\mathrm{\partial }{\beta }_2}=\frac{{\beta }_1\cdot x_i}{{\left({\beta }_2+x_i\right)}^2}.}$$

Starting with the initial estimates of${\displaystyle {\beta }_1=0.9}$ and${\displaystyle {\beta }_2=0.2}$, after five iterations of the Gauss–Newton algorithm, the optimal values${\displaystyle {\hat{\beta }}_1=0.362}$ and${\displaystyle {\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 with the observed data.

The Gauss-Newton iteration is guaranteed to converge toward a local minimum point${\displaystyle \hat{\beta }}$ under 4 conditions:^[4] The functions${\displaystyle r_1,\dots ,r_m}$ are twice continuously differentiable in an open convex set${\displaystyle D\ni \hat{\beta }}$, the Jacobian${\displaystyle {\mathbf{J}}_{\mathbf{r}}(\hat{\beta })}$ is of full column rank, the initial iterate${\displaystyle {\beta }^{(0)}}$ is near${\displaystyle \hat{\beta }}$, and the local minimum value${\displaystyle |S(\hat{\beta })|}$ is small. The convergence is quadratic if${\displaystyle |S(\hat{\beta })|=0}$.

It can be shown^[5] that the increment Δ is adescent direction forS, and, if the algorithm converges, then the limit is astationary point ofS. For large minimum value${\displaystyle |S(\hat{\beta })|}$, however, convergence is not guaranteed, not evenlocal convergence as inNewton's method, or convergence under the usual Wolfe conditions.^[6]

The rate of convergence of the Gauss–Newton algorithm can approachquadratic.^[7] The algorithm may converge slowly or not at all if the initial guess is far from the minimum or the matrix${\displaystyle {\mathbf{J}}_{\mathbf{r}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{r}}}$ isill-conditioned. For example, consider the problem with${\displaystyle m=2}$ equations and${\displaystyle n=1}$ variable, given by

For,${\displaystyle \beta =0}$ is a local optimum. If${\displaystyle \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.^[8]

The Gauss-Newton iteration$${\displaystyle {\mathbf{x}}^{(k+1)}={\mathbf{x}}^{(k)}-J({\mathbf{x}}^{(k)}{)}^{†}\mathbf{f}({\mathbf{x}}^{(k)}),k=0,1,\dots }$$is an effective method for solvingoverdetermined systems of equations in the form of${\displaystyle \mathbf{f}(\mathbf{x})=\mathbf{0}}$ with$${\displaystyle \mathbf{f}(\mathbf{x})=\left[\begin{array}{c}{f}_1(x_1,\dots ,x_n)\\ ⋮\\ f_m(x_1,\dots ,x_n)\end{array}\right]}$$and${\displaystyle m>n}$ where${\displaystyle J(\mathbf{x}{)}^{†}}$ is theMoore-Penrose inverse (also known aspseudoinverse) of theJacobian matrix${\displaystyle J(\mathbf{x})}$ of${\displaystyle \mathbf{f}(\mathbf{x})}$. It can be considered an extension ofNewton's method and enjoys the same local quadratic convergence^[4] toward isolated regular solutions.

If the solution doesn't exist but the initial iterate${\displaystyle {\mathbf{x}}^{(0)}}$ is near a point${\displaystyle \hat{\mathbf{x}}=({\hat{x}}_1,\dots ,{\hat{x}}_n)}$ at which the sum of squares$\sum _{i=1}^m|f_i(x_1,\dots ,x_n){|}^2\equiv \Vert \mathbf{f}(\mathbf{x}){\Vert }_2^2$ reaches a small local minimum, the Gauss-Newton iteration linearly converges to${\displaystyle \hat{\mathbf{x}}}$. The point${\displaystyle \hat{\mathbf{x}}}$ is often called aleast squares solution of the overdetermined system.

In what follows, the Gauss–Newton algorithm will be derived fromNewton'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.^[9]

The recurrence relation for Newton's method for minimizing a functionS of parameters${\displaystyle \mathit{\beta }}$ is$${\displaystyle {\mathit{\beta }}^{(s+1)}={\mathit{\beta }}^{(s)}-{\mathbf{H}}^{-1}\mathbf{g},}$$

whereg denotes thegradient vector ofS, andH denotes theHessian matrix ofS.

Since$S=\sum _{i=1}^m{r}_i^2$, the gradient is given by$${\displaystyle g_j=2\sum _{i=1}^m{r}_i\frac{\mathrm{\partial }{r}_i}{\mathrm{\partial }{\beta }_j}.}$$

Elements of the Hessian are calculated by differentiating the gradient elements,${\displaystyle g_j}$, with respect to${\displaystyle {\beta }_k}$:$${\displaystyle H_{jk}=2\sum _{i=1}^m\left(\frac{\mathrm{\partial }{r}_i}{\mathrm{\partial }{\beta }_j}\frac{\mathrm{\partial }{r}_i}{\mathrm{\partial }{\beta }_k}+r_i\frac{{\mathrm{\partial }}^2{r}_i}{\mathrm{\partial }{\beta }_j\mathrm{\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$${\displaystyle H_{jk}\approx 2\sum _{i=1}^m{J}_{ij}{J}_{ik},}$$

where$J_{ij}=\mathrm{\partial }{r}_i/\mathrm{\partial }{\beta }_j$ are entries of the JacobianJ_r. Note that when the exact hessian is evaluated near an exact fit we have near-zero${\displaystyle r_i}$, so the second term becomes near-zero as well, which justifies the approximation. The gradient and the approximate Hessian can be written in matrix notation as$${\displaystyle \mathbf{g}=2{{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}\mathbf{r},\mathbf{H}\approx 2{{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{r}}.}$$

These expressions are substituted into the recurrence relation above to obtain the operational equations$${\displaystyle {\mathit{\beta }}^{(s+1)}={\mathit{\beta }}^{(s)}+\mathrm{\Delta };\mathrm{\Delta }=-{\left({{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{r}}\right)}^{-1}{{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}\mathbf{r}.}$$

Convergence of the Gauss–Newton method is not guaranteed in all instances. The approximation$${\displaystyle \left|r_i\frac{{\mathrm{\partial }}^2{r}_i}{\mathrm{\partial }{\beta }_j\mathrm{\partial }{\beta }_k}\right|\ll \left|\frac{\mathrm{\partial }{r}_i}{\mathrm{\partial }{\beta }_j}\frac{\mathrm{\partial }{r}_i}{\mathrm{\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:^[10]

1. The function values${\displaystyle r_i}$ are small in magnitude, at least around the minimum.
2. The functions are only "mildly" nonlinear, so that$\frac{{\mathrm{\partial }}^2{r}_i}{\mathrm{\partial }{\beta }_j\mathrm{\partial }{\beta }_k}$ is relatively small in magnitude.

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

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 functionS. An optimal value for${\displaystyle \alpha }$ can be found by using aline search algorithm, that is, the magnitude of${\displaystyle \alpha }$ is determined by finding the value that minimizesS, usually using adirect search method in the interval or abacktracking line search such asArmijo-line search. Typically,${\displaystyle \alpha }$ should be chosen such that it satisfies theWolfe conditions or theGoldstein conditions.^[11]

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 theLevenberg–Marquardt algorithm, atrust region method.^[3] The normal equations are modified in such a way that the increment vector is rotated towards the direction ofsteepest descent,$${\displaystyle \left({\mathbf{J}}^{\mathrm{T}}\mathbf{J}\mathbf{+}\lambda \mathbf{D}\right)\mathrm{\Delta }=-{\mathbf{J}}^{\mathrm{T}}\mathbf{r},}$$

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

The so-called Marquardt parameter${\displaystyle \lambda }$ may also be optimized by a line search, but this is inefficient, as the shift vector must be recalculated every time${\displaystyle \lambda }$ is changed. A more efficient strategy is this: When divergence occurs, increase the Marquardt parameter until there is a decrease inS. 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 ofS then becomes a standard Gauss–Newton minimization.

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

In order to make this kind of approach work, one needs at least an efficient method for computing the product$${\displaystyle {{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{r}}\mathbf{p}}$$

for some vectorp. Withsparse matrix storage, it is in general practical to store the rows of${\displaystyle {\mathbf{J}}_{\mathbf{r}}}$ 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 definesc_i as rowi of the matrix${\displaystyle {\mathbf{J}}_{\mathbf{r}}}$, the following simple relation holds:$${\displaystyle {{\mathbf{J}}_{\mathbf{r}}}^{\mathrm{T}}{\mathbf{J}}_{\mathbf{r}}\mathbf{p}=\sum _i{\mathbf{c}}_i\left({\mathbf{c}}_i\cdot \mathbf{p}\right),}$$

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 forparallel computations. Note that every rowc_i is the gradient of the corresponding residualr_i; with this in mind, the formula above emphasizes the fact that residuals contribute to the problem independently of each other.

In aquasi-Newton method, such as that due toDavidon, Fletcher and Powell or Broyden–Fletcher–Goldfarb–Shanno (BFGS method) an estimate of the full Hessian$\frac{{\mathrm{\partial }}^{2}S}{\mathrm{\partial }{\beta }_j\mathrm{\partial }{\beta }_k}$ is built up numerically using first derivatives$\frac{\mathrm{\partial }{r}_i}{\mathrm{\partial }{\beta }_j}$ only so that aftern 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 isgradient 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.

The following implementation inJulia provides one method which uses a provided Jacobian and another computing withautomatic differentiation.

"""    gaussnewton(r, J, β₀, maxiter, tol)Perform Gauss–Newton optimization to minimize the residual function `r` with Jacobian `J` starting from `β₀`. The algorithm terminates when the norm of the step is less than `tol` or after `maxiter` iterations."""functiongaussnewton(r,J,β₀,maxiter,tol)β=copy(β₀)for_in1:maxiterJβ=J(β);Δ=-(Jβ'*Jβ)\(Jβ'*r(β))β+=Δifsqrt(sum(abs2,Δ))<tolbreakendendreturnβendimportAbstractDifferentiationasAD,Zygotebackend=AD.ZygoteBackend()# other backends are available"""    gaussnewton(r, β₀, maxiter, tol)Perform Gauss–Newton optimization to minimize the residual function `r` starting from `β₀`. The relevant Jacobian is calculated using automatic differentiation. The algorithm terminates when the norm of the step is less than `tol` or after `maxiter` iterations."""functiongaussnewton(r,β₀,maxiter,tol)β=copy(β₀)for_in1:maxiterrβ,Jβ=AD.value_and_jacobian(backend,r,β)Δ=-(Jβ[1]'*Jβ[1])\(Jβ[1]'*rβ)β+=Δifsqrt(sum(abs2,Δ))<tolbreakendendreturnβend

• 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.{{cite book}}: CS1 maint: publisher location (link)

• 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).

• 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.

• Optimization algorithms and methods
• Least squares
• Statistical algorithms

• CS1 maint: archived copy as title
• Articles with short description
• Short description is different from Wikidata
• CS1 maint: publisher location