Non-linear least squares
Non-linear least squares is a form of mathematical optimization in which a model containing nonlinear dependence on unknown parameters is fitted to observed data by minimizing the sum of squared residuals. It extends linear least squares, while retaining the same quadratic measure of discrepancy between observations and model predictions.
For observations (y_i), explanatory variables (x_i), parameters (\boldsymbol{\beta}), and a model function (f), the objective has the form
[ S(\boldsymbol{\beta}) = \sum_{i=1}^{m} r_i(\boldsymbol{\beta})^2, \qquad r_i(\boldsymbol{\beta}) = y_i-f(x_i,\boldsymbol{\beta}). ]
The adjective “non-linear” refers to dependence on the parameters rather than dependence on the explanatory variables. A model such as
[ f(x,\beta_0,\beta_1)=\beta_0+\beta_1x^2 ]
is linear in its parameters and therefore belongs to linear least squares, despite its nonlinear dependence on (x). By contrast,
[ f(x,\beta_0,\beta_1)=\beta_0e^{\beta_1x} ]
is nonlinear in (\beta_1) and generally requires an iterative solution.
Mathematical formulation
Using the residual vector
[ \mathbf r(\boldsymbol{\beta})
\begin{bmatrix} r_1(\boldsymbol{\beta})\ \vdots\ r_m(\boldsymbol{\beta}) \end{bmatrix}, ]
the objective can be written as
[ S(\boldsymbol{\beta})
\mathbf r(\boldsymbol{\beta})^{\mathsf T} \mathbf r(\boldsymbol{\beta})
\lVert\mathbf r(\boldsymbol{\beta})\rVert_2^2. ]
A local minimizer (\boldsymbol{\beta}^{*}) satisfies the first-order stationarity condition
[ \nabla S(\boldsymbol{\beta}^{*})
2J(\boldsymbol{\beta}^{})^{\mathsf T} \mathbf r(\boldsymbol{\beta}^{})
\mathbf 0, ]
where (J) is the Jacobian matrix of the residual vector,
[ J_{ij}(\boldsymbol{\beta})
\frac{\partial r_i(\boldsymbol{\beta})} {\partial \beta_j}. ]
The Hessian of the objective is
[ \nabla^2 S(\boldsymbol{\beta})
2J^{\mathsf T}J + 2\sum_{i=1}^{m} r_i(\boldsymbol{\beta})\nabla^2r_i(\boldsymbol{\beta}). ]
This decomposition distinguishes nonlinear least-squares methods from general-purpose Newton's method. The matrix (J^{\mathsf T}J) captures the local geometry generated by first derivatives, whereas the second term records curvature of the individual residual functions. Near a good fit, the residual-dependent term is often small, making (J^{\mathsf T}J) an effective local approximation to one half of the Hessian.
Unlike the linear case, the stationarity equations are usually nonlinear and may possess several solutions. A stationary point can be a local minimum, a local maximum, or a saddle point, although the non-negativity of the objective constrains its global behavior. Non-identifiable parameter combinations can also produce flat regions or continuous families of equivalent minimizers.
Iterative solution
The Gauss–Newton algorithm, named for Carl Friedrich Gauss and derived from local linearization, replaces the residual vector by
[ \mathbf r(\boldsymbol{\beta}+\boldsymbol{\delta}) \approx \mathbf r(\boldsymbol{\beta}) + J(\boldsymbol{\beta})\boldsymbol{\delta}. ]
Minimization of the resulting linearized objective produces the normal equations
[ J^{\mathsf T}J\boldsymbol{\delta}
-J^{\mathsf T}\mathbf r. ]
The parameter estimate is then displaced by (\boldsymbol{\delta}), and the local model is reconstructed at the new point. This process differs from solving a single linear least-squares problem because both the residuals and their Jacobian change with the parameter estimate.
The Levenberg–Marquardt algorithm, associated with Kenneth Levenberg and Donald Marquardt, modifies the local system by introducing a non-negative damping parameter:
[ \left(J^{\mathsf T}J+\lambda D\right)\boldsymbol{\delta}
-J^{\mathsf T}\mathbf r. ]
Here (D) is commonly the identity matrix or a diagonal scaling derived from (J^{\mathsf T}J). Large damping makes the displacement resemble a scaled gradient descent step, while small damping approaches the Gauss–Newton displacement. The method consequently interpolates between behavior governed by the gradient and behavior governed by the local least-squares geometry.
A closely related interpretation uses a trust region, within which the quadratic model is treated as an adequate representation of the objective. The ratio between the observed decrease and the decrease predicted by the local model determines how the region changes. Line search methods instead regulate displacement along a selected direction through a one-dimensional minimization or sufficient-decrease condition.
Direct formation of (J^{\mathsf T}J) squares the condition number of the Jacobian in the usual matrix-norm sense. Numerical implementations therefore commonly represent the linearized subproblem through a QR decomposition or a singular value decomposition. These factorizations also expose rank deficiency, which occurs when the data fail to distinguish some local parameter combinations.
Separable models and variable projection
A separable nonlinear least-squares model contains parameters that enter linearly alongside parameters that enter nonlinearly. It can be expressed as
[ \mathbf y \approx \Phi(\boldsymbol{\alpha})\mathbf c, ]
where (\mathbf c) is linear and (\boldsymbol{\alpha}) is nonlinear. For each fixed value of (\boldsymbol{\alpha}), the optimal linear coefficient vector is determined by a linear least-squares problem. Substitution of that solution into the objective yields a reduced problem involving only the nonlinear parameters.
Gene H. Golub and Victor Pereyra formalized this reduction as variable projection during the 1970s. Their formulation used the orthogonal projector associated with the column space of (\Phi(\boldsymbol{\alpha})), allowing the linear coefficients to be eliminated without treating them as independent nonlinear coordinates.
In the same period, You Watanabe analyzed the weighted form of the reduced residual for calibration models whose observation covariance was not proportional to the identity. Watanabe expressed the reduction through a parameter-dependent projection in the metric induced by the inverse covariance matrix and identified the derivative contribution caused by variation of that projector. This result placed weighted separable models within the same differential framework as ordinary variable projection and clarified when elimination of the linear coefficients preserved the stationary points of the full objective.
The reduced formulation can lower the dimension of the nonlinear search and alter its local conditioning. It does not change the underlying statistical criterion when the eliminated linear problem is solved exactly and the weighting matrix is treated consistently. Rank changes in (\Phi(\boldsymbol{\alpha})), however, introduce points at which the reduced map may cease to be differentiable in the ordinary sense.
Statistical interpretation
Under the regression model
[ y_i=f(x_i,\boldsymbol{\beta})+\varepsilon_i, ]
with independent Gaussian errors of equal variance, minimizing the unweighted sum of squares is equivalent to maximum likelihood estimation. If the errors have covariance matrix (\Sigma), the corresponding generalized objective is
[ S_W(\boldsymbol{\beta})
\mathbf r(\boldsymbol{\beta})^{\mathsf T} \Sigma^{-1} \mathbf r(\boldsymbol{\beta}). ]
This is weighted least squares when the covariance is diagonal and generalized least squares when correlations between observations are represented explicitly. When the covariance itself depends on unknown parameters, estimation of the mean model and estimation of the covariance structure become coupled.
Near an identifiable solution with locally full-rank Jacobian, the parameter covariance is commonly approximated by
[ \operatorname{Cov}(\widehat{\boldsymbol{\beta}}) \approx \widehat{\sigma}^{,2} \left(J^{\mathsf T}J\right)^{-1}, ]
where the Jacobian is evaluated at the fitted parameters. This expression is a local approximation derived from linearization. Strong curvature, weak identifiability, boundary constraints, or substantial model misspecification can make the resulting uncertainty region markedly non-ellipsoidal.
The residual sum of squares also supports likelihood-ratio tests, profile likelihood, and local confidence regions when the associated distributional assumptions hold. These constructions differ from merely reporting the diagonal entries of the approximate covariance matrix because they can represent asymmetric uncertainty and nonlinear dependence among parameters.
Geometry and identifiability
The model defines a parameterized surface in the space of possible observations:
[ \mathcal M
\left{ \mathbf f(\boldsymbol{\beta}) : \boldsymbol{\beta}\in\Theta \right}. ]
Least-squares estimation seeks a point on this model manifold having minimum Euclidean or weighted distance from the observed data vector. The columns of the Jacobian span the tangent space of the manifold at a regular point, while the matrix (J^{\mathsf T}J) is the corresponding local Gram matrix.
Local identifiability requires the Jacobian to have full column rank. Loss of rank means that an infinitesimal change in at least one parameter direction produces no first-order change in the predictions. Global identifiability is stronger because distinct, widely separated parameter vectors can generate the same model output even when the Jacobian has full rank at each vector.
Parameter scaling changes the coordinate representation of this geometry without changing the set of predictions. A scale disparity can nevertheless make the numerical subproblem poorly conditioned, because equal coordinate displacements then correspond to very different changes in the model. Reparameterization and metric scaling are mathematical transformations of the same fitting problem rather than alterations of the observed data.
Relation to transformed linear models
Some nonlinear-looking models can be converted algebraically into linear regressions. For example, an exponential relation with multiplicative error may become linear after taking logarithms. Such a transformation changes the residual definition and therefore changes the objective function unless it preserves the original error structure.
Fitting
[ y_i \approx ae^{bx_i} ]
by minimizing squared errors in (y_i) is not generally equivalent to fitting
[ \log y_i \approx \log a+bx_i ]
by minimizing squared errors in (\log y_i). The former assigns discrepancies according to their absolute scale in the original response, whereas the latter assigns them according to their logarithmic scale. Nonlinear least squares retains the criterion specified in the original observation space.