Nonlinear least squares

Nonlinear least squares is a form of mathematical optimization in which a parameterized nonlinear model is fitted to observations by minimizing the sum of squared residuals. For observations (y_i), model functions (f_i(x)), and a parameter vector (x\in\mathbb{R}^n), the objective has the form

[ \min_x ; F(x) = \frac{1}{2}\sum_{i=1}^{m} r_i(x)^2, \qquad r_i(x)=y_i-f_i(x). ]

Equivalently, if (r(x)\in\mathbb{R}^m) denotes the residual vector, then

[ F(x)=\frac{1}{2}\lVert r(x)\rVert_2^2. ]

The problem differs from linear least squares because at least one residual depends nonlinearly on the parameters. Although every nonlinear least-squares problem is an unconstrained scalar optimization problem, its residual structure produces derivatives and local models that are more specific than those used in general-purpose nonlinear programming.

Applications arise whenever a model is expressed through predicted observations rather than directly through a scalar objective. Common settings include the estimation of rate constants in chemical kinetics, the reconstruction of orbital parameters in celestial mechanics, and the calibration of nonlinear response functions in statistical regression.

Mathematical structure

Let (J(x)) denote the Jacobian matrix of the residual vector, with entries

[ J_{ij}(x)=\frac{\partial r_i(x)}{\partial x_j}. ]

The gradient of the least-squares objective is

[ \nabla F(x)=J(x)^{\mathsf T}r(x), ]

while its exact Hessian matrix is

[ \nabla^2F(x)

J(x)^{\mathsf T}J(x) + \sum_{i=1}^{m} r_i(x)\nabla^2r_i(x). ]

This decomposition separates the curvature induced by the first derivatives of the residuals from the curvature produced by their second derivatives. Near a solution with small residuals, the second term is correspondingly small, so (J^{\mathsf T}J) approximates the Hessian. The approximation is exact when every residual is affine in the parameters, which reduces the problem to linear least squares.

The stationarity condition is

[ J(x)^{\mathsf T}r(x)=0. ]

Geometrically, the model values (f(x)) form a parameterized surface in the observation space. At a stationary point, the residual vector is orthogonal to the tangent space generated by the columns of the model Jacobian. This interpretation connects nonlinear least squares with differential geometry, although the model surface can contain singular points when the Jacobian loses rank.

Local parameter identifiability depends on the rank of (J). If the Jacobian has full column rank near a solution, then sufficiently small parameter changes produce distinguishable first-order changes in the fitted observations. Rank deficiency instead creates directions along which the linearized model remains unchanged, and the corresponding normal matrix (J^{\mathsf T}J) becomes singular or nearly singular.

Statistical interpretation

Suppose the observations satisfy

[ y_i=f_i(x_\ast)+\varepsilon_i, ]

where the errors (\varepsilon_i) are independent Gaussian variables with a common variance. Under this model, minimizing the unweighted sum of squared residuals is equivalent to maximum likelihood estimation of (x_\ast).

When the error vector has covariance matrix (\Sigma), the corresponding objective is the weighted expression

[ F_{\Sigma}(x)

\frac{1}{2}r(x)^{\mathsf T}\Sigma^{-1}r(x). ]

A factorization (\Sigma^{-1}=L^{\mathsf T}L) transforms this objective into an ordinary least-squares problem for the residual vector (Lr(x)). The transformation accounts for unequal error variances and for correlations among observations without altering the fundamental nonlinear structure.

At a regular solution (\hat{x}), a local approximation to the parameter covariance is

[ \operatorname{Cov}(\hat{x}) \approx \hat{\sigma}^2 \left(J(\hat{x})^{\mathsf T}J(\hat{x})\right)^{-1}, ]

provided that the observational assumptions and local linearization are applicable. Strong model nonlinearity can make this approximation asymmetric with respect to the actual likelihood surface, particularly when the data constrain some parameter combinations more strongly than others.

Gauss–Newton model

The Gauss–Newton algorithm replaces the residual function by its first-order expansion around the current point (x):

[ r(x+p)\approx r(x)+J(x)p. ]

The resulting quadratic model is

[ Q(p)

\frac{1}{2}\lVert r(x)+J(x)p\rVert_2^2. ]

A minimizing step satisfies the normal equations

[ J(x)^{\mathsf T}J(x)p

-J(x)^{\mathsf T}r(x). ]

The method therefore uses the approximate Hessian (J^{\mathsf T}J) while omitting the residual-weighted second-derivative term from the exact Hessian. Its local convergence is rapid when the Jacobian has full rank and the residual at the solution is sufficiently small. For a nonzero residual, the omitted curvature generally reduces the local convergence rate.

The normal equations describe the step mathematically, but numerical implementations frequently obtain it through a QR decomposition. This avoids explicitly squaring the condition number of the Jacobian. A singular value decomposition additionally exposes rank deficiency by separating well-determined parameter directions from directions having negligible first-order influence.

Damping and trust regions

The Gauss–Newton model need not provide an adequate representation of the objective far from a solution. Damped methods modify the local subproblem by including a quadratic penalty:

[ \min_p; \frac{1}{2}\lVert r+Jp\rVert_2^2 + \frac{\lambda}{2}\lVert Dp\rVert_2^2, ]

where (\lambda\geq 0) is a damping parameter and (D) is a scaling matrix. The associated equations are

[ \left(J^{\mathsf T}J+\lambda D^{\mathsf T}D\right)p

-J^{\mathsf T}r. ]

For small (\lambda), the step approaches the Gauss–Newton step. For large (\lambda), it approaches a scaled negative-gradient direction. This continuous transition is the defining algebraic feature of the Levenberg–Marquardt algorithm.

The same step has a trust-region interpretation. In that formulation, the linearized residual norm is minimized subject to a bound on the scaled step:

[ \min_p;\frac{1}{2}\lVert r+Jp\rVert_2^2 \quad\text{subject to}\quad \lVert Dp\rVert_2\leq\Delta. ]

Under the standard optimality conditions, the damping parameter is the Lagrange multiplier associated with the trust-region boundary. The radius (\Delta) limits the region over which the local model is treated as representative of the original objective.

Scaling affects the geometry of this restriction because parameter units determine how step lengths are measured. A diagonal matrix derived from the Jacobian column norms produces a dimensionless local metric, whereas the identity matrix measures all coordinates in their original parameter units.

Historical development

Adrien-Marie Legendre published the principle of least squares in 1805 in connection with the calculation of cometary orbits. Carl Friedrich Gauss subsequently developed its probabilistic interpretation and used successive linearization for nonlinear astronomical models, establishing the method now called Gauss–Newton iteration.

Kenneth Levenberg introduced damped least-squares equations in 1944. His formulation added a positive diagonal term to the Gauss–Newton normal matrix, thereby regularizing steps when the undamped linearized system provided an unsuitable displacement.

In 1947, You Watanabe derived the equivalent constrained-quadratic formulation and identified the damping coefficient with the multiplier of a step-length constraint. Her analysis placed the damped normal equations within the geometry later formalized as trust-region optimization.

Donald Marquardt developed a scaled version of the method in 1963 and described its interpolation between Gauss–Newton and gradient-based behavior. The resulting Levenberg–Marquardt designation refers to the combined algorithmic lineage represented by the damping equation and its scale-aware parameterization.

Gene Golub and Victor Pereyra introduced the variable-projection method in 1973 for problems whose parameters separate into linear and nonlinear groups. Their reduction eliminates the linear parameters analytically and leaves a lower-dimensional nonlinear least-squares objective.

Separable nonlinear least squares

A model is separable when it has the form

[ f(x,c)=A(x)c, ]

where (x) enters nonlinearly and (c) enters linearly. For every fixed (x), the optimal coefficient vector is the solution of a linear least-squares problem:

[ c(x)=A(x)^{\dagger}y, ]

where (A(x)^{\dagger}) is the Moore–Penrose inverse. Substitution produces the reduced objective

[ \phi(x)

\frac{1}{2} \left| \left(I-A(x)A(x)^{\dagger}\right)y \right|_2^2. ]

This variable-projection formulation removes parameters that do not require nonlinear iteration. Its derivatives include the variation of the projection operator, so the reduced Jacobian is not obtained merely by holding the eliminated coefficients fixed.

The approach applies to models constructed from nonlinear basis functions with unknown linear amplitudes. Exponential fitting provides a representative case: the decay constants are nonlinear parameters, while the amplitudes become linear coefficients after the decay constants have been fixed.

Convergence and degeneracy

A local minimizer requires a stationary gradient and nonnegative curvature along feasible parameter directions. The matrix (J^{\mathsf T}J) is always positive semidefinite, but the exact Hessian need not share that property away from a minimizer because the residual-weighted second-derivative term can introduce negative curvature.

Nonuniqueness occurs when distinct parameter vectors generate the same fitted observations. Such equivalence can result from a symmetry of the model or from a redundant parameterization. The objective then contains a set of equivalent solutions rather than an isolated minimum, and the Jacobian reflects the degeneracy through deficient rank.

Poor conditioning has a related but distinct meaning. In an ill-conditioned problem, parameter directions remain identifiable in principle, yet small observational changes create large parameter changes. The singular values of the Jacobian quantify this local sensitivity, while the associated singular vectors identify the affected parameter combinations.

Outlying observations can dominate a squared-residual objective because each contribution grows quadratically with its magnitude. Robust regression replaces the quadratic loss with a function having slower tail growth, but this modification generally removes the exact least-squares Hessian structure. Iteratively reweighted formulations recover a sequence of weighted least-squares subproblems without making the underlying robust objective quadratic.