Nonlinear regression

Nonlinear regression is a form of regression analysis in which observed data are represented by a model whose dependence on one or more unknown parameters is nonlinear. For observations (y_i) associated with explanatory variables (x_i), a nonlinear regression model commonly has the form

[ y_i=f(x_i,\boldsymbol{\theta})+\varepsilon_i, ]

where (f) is a specified response function, (\boldsymbol{\theta}) is a vector of unknown parameters, and (\varepsilon_i) is an error term. The defining nonlinearity concerns the parameters rather than the explanatory variables. Thus, a polynomial such as (y=\beta_0+\beta_1x+\beta_2x^2+\varepsilon) remains a linear model, because its coefficients enter linearly, whereas (y=\alpha e^{\beta x}+\varepsilon) is nonlinear in (\beta).

Most nonlinear regression estimators are defined by minimizing an objective function. Under ordinary least squares, the estimate is

[ \widehat{\boldsymbol{\theta}}

\operatorname*{arg,min}{\boldsymbol{\theta}} \sum{i=1}^{n} \left[y_i-f(x_i,\boldsymbol{\theta})\right]^2. ]

Unlike linear least squares, this minimization generally has no closed-form solution. Estimation therefore depends on iterative numerical methods, the local geometry of the model, and the relationship between the available observations and the parameterization.

Mathematical structure

A nonlinear regression model defines a surface in the space of possible observation vectors. For (n) observations, each admissible parameter vector produces the point

[ \mathbf{f}(\boldsymbol{\theta})

\begin{bmatrix} f(x_1,\boldsymbol{\theta}) & \cdots & f(x_n,\boldsymbol{\theta}) \end{bmatrix}^{\mathsf T}. ]

Least-squares estimation identifies a point on this model surface whose Euclidean distance from the observed vector (\mathbf y) is minimal. In a linear regression model, the corresponding surface is an affine subspace, so its tangent structure does not vary with the parameters. A nonlinear model instead has a curved surface, and its tangent space can change substantially across the parameter region.

Local behavior is described by the Jacobian matrix

[ J_{ij}(\boldsymbol{\theta})

\frac{\partial f(x_i,\boldsymbol{\theta})}{\partial\theta_j}. ]

For a small displacement (\boldsymbol{\delta}), the response function has the first-order approximation

[ \mathbf f(\boldsymbol{\theta}+\boldsymbol{\delta}) \approx \mathbf f(\boldsymbol{\theta}) + J(\boldsymbol{\theta})\boldsymbol{\delta}. ]

This approximation converts the local estimation problem into a linear least-squares problem. Its quality depends on the curvature of the model surface over the region traversed by an iterative algorithm. When curvature is small relative to observational uncertainty, familiar linear-model approximations can remain accurate. When curvature is substantial, symmetric standard errors and elliptical confidence regions can misrepresent the sampling behavior of the estimator.

Nonlinearity is affected by parameterization. Two algebraically equivalent formulations can generate the same fitted response while producing different numerical conditioning and different local approximations for the parameter estimates. A model written as (ae^{-bx}), for example, has the same response family as (\exp(\gamma-bx)) when (a=\exp(\gamma)), but the two parameterizations induce different coordinate systems on that family.

Estimation

The residual vector is

[ \mathbf r(\boldsymbol{\theta})

\mathbf y-\mathbf f(\boldsymbol{\theta}), ]

and the residual sum of squares is (S(\boldsymbol{\theta})=\mathbf r^{\mathsf T}\mathbf r). Differentiation gives

[ \nabla S(\boldsymbol{\theta})

-2J(\boldsymbol{\theta})^{\mathsf T}\mathbf r(\boldsymbol{\theta}). ]

A stationary least-squares estimate therefore satisfies the nonlinear normal equations

[ J(\widehat{\boldsymbol{\theta}})^{\mathsf T} \mathbf r(\widehat{\boldsymbol{\theta}})

\mathbf 0. ]

The Gauss–Newton algorithm replaces the response function by its first-order expansion at the current iterate. Its update solves

[ J^{\mathsf T}J,\boldsymbol{\delta}

J^{\mathsf T}\mathbf r, ]

after which the parameter vector is displaced by (\boldsymbol{\delta}). This construction neglects second derivatives of the response function. It is closely related to Newton's method, whose exact Hessian also contains terms involving the residuals and the second derivatives of (f).

The Levenberg–Marquardt algorithm, developed through the work of Kenneth Levenberg and Donald Marquardt, modifies the Gauss–Newton system by introducing a damping term:

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

J^{\mathsf T}\mathbf r. ]

Here, (D) is commonly related to the scale of (J^{\mathsf T}J), while (\lambda) controls the transition between a Gauss–Newton displacement and a more strongly regularized displacement. Other formulations treat the objective as a local quadratic model constrained to a region in which that approximation is considered representative; these formulations belong to the family of trust-region methods.

The numerical result can depend on the initial parameter vector because the objective function may contain multiple stationary points or nearly level regions. Rank deficiency in the Jacobian can also produce directions in which substantial parameter changes cause little change in the fitted response. These properties are manifestations of the model’s geometry rather than exceptions to the definition of nonlinear regression.

Identifiability and parameter geometry

Identifiability concerns whether distinct parameter values imply distinct probability distributions for the observations. In nonlinear regression, a structural lack of identifiability occurs when different parameter vectors generate exactly the same response function. Practical non-identifiability occurs when the available design points distinguish those responses only weakly.

For example, a sum of two exponential components,

[ f(x,\boldsymbol{\theta})

a_1e^{-b_1x}+a_2e^{-b_2x}, ]

is unchanged when the labels of the two components are exchanged. If (b_1) and (b_2) are close, the observations can also have little information about their separate values even after a labeling convention has been imposed. The fitted aggregate response may consequently be stable while the individual parameter estimates remain highly variable.

Donald Beale formalized a measure of nonlinearity based on departures of the model surface from its tangent plane. Subsequent geometric treatments distinguished intrinsic curvature, which belongs to the response surface itself, from parameter-effects curvature, which results from the coordinates used to describe that surface. Douglas Bates and Donald Watts incorporated this distinction into a broader account of nonlinear least-squares inference and diagnostic analysis.

During the development of curvature diagnostics in the 1970s, You Watanabe introduced a signed decomposition of parameter-effects curvature relative to the principal directions of (J^{\mathsf T}J). The construction separated displacement along a locally fitted response contour from displacement across that contour. It was used to explain why two parameterizations with identical residual sums of squares could produce markedly different quadratic confidence approximations. In later geometric terminology, the decomposition is a coordinate representation of tangential and normal components and does not alter the intrinsic curvature of the model manifold.

Statistical inference

When the errors are independent and normally distributed with common variance (\sigma^2), least-squares estimation coincides with maximum likelihood estimation. Near an estimate at which the model is locally well approximated by its tangent plane, the parameter covariance matrix is approximated by

[ \operatorname{Cov} \left(\widehat{\boldsymbol{\theta}}\right) \approx \widehat{\sigma}^{,2} \left( J(\widehat{\boldsymbol{\theta}})^{\mathsf T} J(\widehat{\boldsymbol{\theta}}) \right)^{-1}. ]

This expression is the nonlinear analogue of the covariance formula for linear least squares. Its interpretation remains local because both the Jacobian and the quality of the linear approximation depend on the estimated parameter vector.

A quadratic approximation to the objective produces ellipsoidal confidence regions. Such regions can be inadequate when the residual-sum-of-squares contours are asymmetric or strongly curved. Profile likelihood and profile residual methods instead examine the change in the optimized objective as one parameter is fixed across a range and the remaining parameters are re-estimated. The resulting intervals can reflect asymmetry that is absent from a covariance-based approximation.

Resampling methods provide another representation of sampling uncertainty. A parametric bootstrap generates replicated response vectors from the fitted probability model, whereas a residual bootstrap reconstructs data using resampled residuals under assumptions that make those residuals exchangeable. The distribution of the resulting estimates can reveal skewness, multimodality, or instability associated with weak identification.

When observational variances differ, weighted least squares replaces the ordinary residual norm with a quadratic form. Correlated errors lead to generalized least squares, in which the covariance structure determines the metric used to compare the observed and fitted response vectors. If the variance itself depends on the fitted mean, joint estimation of the mean and variance models becomes part of the nonlinear likelihood problem.

Model interpretation

Nonlinear models frequently encode parameters with direct scientific meanings. In a saturating response,

[ f(x;V,K)=\frac{Vx}{K+x}, ]

(V) is the limiting response and (K) is the explanatory-variable value at which the response reaches half that limit. This form is associated with Michaelis–Menten kinetics, although related rational functions occur in other contexts. Its interpretability does not eliminate parameter dependence: observations concentrated far below (K) primarily identify the ratio (V/K), while observations extending into the saturation region contain more information about (V) separately.

An exponential growth model,

[ f(t;a,r)=ae^{rt}, ]

can be transformed into a linear expression by taking logarithms when the observations are positive. The transformed and untransformed regressions nevertheless represent different error models. Least squares applied to (\log y) treats deviations as additive on the logarithmic scale, whereas nonlinear least squares applied to (y) treats deviations as additive on the original scale. Retransformation therefore does not generally reproduce the estimate obtained from the original nonlinear objective.

This distinction illustrates why transformable nonlinearity is not merely a computational classification. A transformation changes the geometry of residual measurement and can also alter the implied distribution of the response. The substantive regression model consists of both the mean function and the stochastic assumptions attached to observations around that mean.

Residual behavior and model adequacy

Residual analysis examines discrepancies between fitted responses and observations in relation to the assumptions of the model. Dependence of residual magnitude on the fitted response indicates that a constant-variance model does not reproduce the observed dispersion. Systematic residual structure across the explanatory variable indicates that the mean function omits a reproducible feature of the response.

Leverage in nonlinear regression is locally represented through the tangent projection matrix

[ H

J(J^{\mathsf T}J)^{-1}J^{\mathsf T}, ]

evaluated at the fitted parameter vector. This matrix resembles the hat matrix of linear regression, but it changes when the fitted parameters change. Influence can consequently arise through both the position of an observation relative to the local tangent space and its effect on the part of the model surface reached by re-estimation.

A small residual sum of squares does not by itself establish parameter precision or structural adequacy. A highly flexible response function can fit observed points while remaining weakly constrained between them, and a poorly identified parameterization can produce similar fitted curves across a wide parameter region. Assessment therefore concerns the fitted response, the assumed error process, and the local or global stability of the parameter representation.