Nonlinear Regression

Nonlinear Regression

Nonlinear regression fits a response $y$ to predictors $\mathbf{x}$ through a known function $f(\mathbf{x},\mathbf{\theta })$ that depends on parameters $\mathbf{\theta }$ in a non-linear way:
$y_i=f({\mathbf{x}}_i,\mathbf{\theta })+{\epsilon }_i,i=1,\dots ,n,$
where $\mathbf{\theta }\in {\mathbb{R}}^q$ is the parameter vector and ${\epsilon }_i$ is a noise term, typically assumed $\mathrm{WN}(0,{\sigma }^2)$. The standard estimator is nonlinear least squares (NLS), which minimizes the residual sum of squares over the parameters:
$\hat{\mathbf{\theta }}=\arg {\min }_{\mathbf{\theta }}{\sum }_{i=1}^n[y_i-f({\mathbf{x}}_i,\mathbf{\theta }){]}^2.$
Unlike linear regression, this minimization has no closed-form solution in general and must be solved iteratively. The standard algorithms are Gauss-Newton, which linearizes $f$ around the current parameter estimate and solves the resulting linear regression at each step, and Levenberg-Marquardt, a damped variant that interpolates between Gauss-Newton and gradient descent.

An important distinction is between linearizable and intrinsically nonlinear models. A model like $y={\beta }_0{e}^{{\beta }_{1}t}\cdot \epsilon$ with multiplicative log-normal noise becomes $\log y=\log {\beta }_0+{\beta }_{1}t+\log \epsilon$, which is linear in the parameters and can be fit with OLS after the log transform. Such linearizable models inherit the exact finite-sample theory of linear regression. Intrinsically nonlinear models — those where no transformation reduces $f$ to linear in the parameters — require iterative NLS and only have asymptotic (large-$n$) inferential theory.

Local Minima $Q(\mathbf{\theta })$ is non-convex in general, so iterative methods can converge to local minima rather than the global one. The choice of starting values is therefore critical. Standard practice is to try multiple starting points, examine the profile likelihood to understand the surface, or use global search methods (simulated annealing, genetic algorithms, basin-hopping) for problems where convergence is unreliable.

The NLS objective

Intuition

Nonlinear regression is needed whenever the underlying mechanism produces a response that is not a linear combination of inputs. Population growth that saturates follows a logistic curve; chemical reactions decay exponentially; physical oscillations follow sinusoids whose frequency itself is unknown. In each case, the relationship between predictors and response is governed by a parametric shape whose parameters enter the model non-linearly, so the closed-form OLS formula does not apply.

The Gauss-Newton algorithm is the workhorse for these problems. The trick is to recognize that near the current parameter estimate ${\mathbf{\theta }}^{(m)}$, the function $f$ can be approximated by its first-order Taylor expansion, which is linear in the parameter increment $\mathbf{\theta }-{\mathbf{\theta }}^{(m)}$. Solving this linearized problem by ordinary least squares produces the next iterate ${\mathbf{\theta }}^{(m+1)}$, and the process repeats until convergence. So nonlinear regression is essentially a sequence of linear regressions on the linearized residuals, each one bringing the estimate closer to the NLS minimum.

The non-convexity of $Q$ is what makes NLS qualitatively harder than OLS. The OLS objective is a convex quadratic with a unique minimum that any decent optimizer finds; the NLS objective can have many local minima, and which one Gauss-Newton converges to depends entirely on where it starts. This is why parameter initialization is so important in nonlinear regression and why production code typically tries several starts and reports the best.

Properties

(Gauss-Newton Algorithm) Linearize, solve a linear regression, iterate.

At iteration $m$, expand $f({\mathbf{x}}_i,\mathbf{\theta })\approx f({\mathbf{x}}_i,{\mathbf{\theta }}^{(m)})+{\mathbf{F}}_i^{(m)}(\mathbf{\theta }-{\mathbf{\theta }}^{(m)})$ where ${\mathbf{F}}_i^{(m)}=\partial f/\partial {\mathbf{\theta }}^{\top }{|}_{{\mathbf{\theta }}^{(m)}}$ is the gradient row vector. Substituting into the NLS objective produces a linear regression of the current residuals ${\mathbf{r}}^{(m)}=\mathbf{y}-\mathbf{f}({\mathbf{\theta }}^{(m)})$ on the Jacobian ${\mathbf{F}}^{(m)}$, with OLS solution that gives the parameter update
${\mathbf{\theta }}^{(m+1)}={\mathbf{\theta }}^{(m)}+({\mathbf{F}}^{(m)\top }{\mathbf{F}}^{(m)}{)}^{-1}{\mathbf{F}}^{(m)\top }{\mathbf{r}}^{(m)}.$
This is identical to the OLS update for a linear regression with design matrix ${\mathbf{F}}^{(m)}$ and response ${\mathbf{r}}^{(m)}$. Iterating produces a sequence of parameter estimates that converge (locally) to the NLS minimum.

(Levenberg-Marquardt Damping) Interpolate between Gauss-Newton and gradient descent.

Gauss-Newton converges fast near the minimum but can overshoot or fail to converge when started far away. The Levenberg-Marquardt modification adds a damping factor $\lambda$ to the Hessian approximation:
${\mathbf{\theta }}^{(m+1)}={\mathbf{\theta }}^{(m)}+({\mathbf{F}}^{(m)\top }{\mathbf{F}}^{(m)}+\lambda \mathbf{I}{)}^{-1}{\mathbf{F}}^{(m)\top }{\mathbf{r}}^{(m)}.$
With $\lambda =0$ this is pure Gauss-Newton; as $\lambda \to \infty$ it becomes gradient descent with a small step. The implementation adapts $\lambda$ at every step: shrink it after a successful iteration, grow it after a failed one. This is the algorithm scipy’s optimize.curve_fit uses by default and is generally more robust than vanilla Gauss-Newton.

(Asymptotic Distribution) Large-$n$ inference uses the Jacobian instead of the design matrix.

Under regularity conditions, as $n\to \infty$,
$\hat{\mathbf{\theta }}\stackrel{\text{approx}}{\sim }\mathcal{N}(\mathbf{\theta },{\sigma }^2({\mathbf{F}}^{\top }\mathbf{F}{)}^{-1}),$
where $\mathbf{F}$ is the $n\times q$ Jacobian matrix of partial derivatives evaluated at $\hat{\mathbf{\theta }}$. The structure mirrors the linear-regression result $\hat{\mathbf{\beta }}\sim \mathcal{N}(\mathbf{\beta },{\sigma }^2({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1})$ with the design matrix $\mathbf{X}$ replaced by the Jacobian $\mathbf{F}$. So standard errors, confidence intervals, and hypothesis tests are computed exactly as in linear regression, but using the Jacobian. The approximation is asymptotic, so for small samples the reported standard errors should be treated with caution.

(Residual Variance) ${\hat{\sigma }}^2=\Vert \mathbf{y}-\mathbf{f}(\hat{\mathbf{\theta }}){\Vert }^2/(n-q)$ adjusts for fitted parameters.

Same form as in linear regression: residual sum of squares divided by degrees of freedom. The denominator $n-q$ subtracts the number of fitted parameters, which is the count of independent constraints the optimization has imposed on the residuals.

(Linearizable vs. Intrinsically Nonlinear) Some non-linear models become linear after a transformation.

Examples of linearizable models: $y={\beta }_0{e}^{{\beta }_{1}x}$ with multiplicative log-normal noise becomes linear in $\log y$; $y={\beta }_0+{\beta }_1/x$ becomes linear in the new predictor $1/x$; $y=({\beta }_0+{\beta }_{1}x{)}^{-1}$ becomes linear in $1/y$. The advantage of recognizing the linearization is that OLS gives an exact closed-form fit with exact finite-sample inference. The limit is that the noise structure changes under the transformation: log-transforming Gaussian noise gives non-Gaussian noise on the log scale, and vice versa, so the right noise assumption depends on which scale is “natural” for the application. Intrinsically nonlinear models (no transformation linearizes them) include the logistic curve $y=K/(1+e^{-\beta (t-t_0)})$ and the sinusoid with unknown frequency $y=A\cos (2\pi ft+\varphi )$.

(Local Minima) The objective surface is non-convex; multiple starts are mandatory.

The NLS sum of squares $Q(\mathbf{\theta })={\sum }_i(y_i-f({\mathbf{x}}_i,\mathbf{\theta }){)}^2$ has the same kind of multimodality that any non-convex function can have: local minima where the gradient vanishes but the value is higher than the global minimum. Iterative methods like Gauss-Newton converge to whichever minimum is closest in the basin of attraction of the starting point. Standard mitigations: try several starting values, use a global search like simulated annealing for the initial pass, examine the loss surface visually for low-dimensional problems, or check the residual pattern at convergence (correctly fit residuals look like white noise; locally minimized residuals retain visible structure).

(Model Selection) Use information criteria to compare nested NLS fits.

The same AIC and BIC apply to NLS:
$\mathrm{AIC}=n\log ({\hat{\sigma }}^2)+2q,\mathrm{BIC}=n\log ({\hat{\sigma }}^2)+q\log n,$
where $q$ is the parameter count. BIC penalizes complexity more heavily and tends to select simpler models. Both criteria assume Gaussian residuals and a correctly specified family; for non-Gaussian noise or a misspecified family, cross-validation is the safer choice.

Examples

Logistic growth model

Cumulative adoption of a technology is modeled as $y_t=K/(1+{\theta }_2{e}^{-{\theta }_{3}t})+{\epsilon }_t$, where $K$ is the carrying capacity, ${\theta }_2$ controls the initial level, and ${\theta }_3$ is the growth rate.

Why this is nonlinear and how to fit it $(K,{\theta }_2,{\theta }_3)$ — no transformation reduces it to linear-in-parameters form. Fitting requires NLS, typically Gauss-Newton with smart starting values: $\hat{K}$ initialized from the observed maximum, ${\hat{\theta }}_3$ from the slope of the log-rate at the inflection point, and ${\hat{\theta }}_2$ from the early-data exponential growth rate. Asymptotic standard errors come from ${\hat{\sigma }}^2({\mathbf{F}}^{\top }\mathbf{F}{)}^{-1}$ at the converged estimate, and residual diagnostics (autocorrelation, normality) check whether the fit is honest.

The model is intrinsically nonlinear in the parameter vector

Linearizing exponential growth

Bacterial population data follows $N(t)=N_0{e}^{rt}\cdot {\epsilon }_t$ with multiplicative noise. Estimate $N_0$ and $r$.

Solution $\log N(t)=\log N_0+rt+\log {\epsilon }_t$. Now the model is linear in the parameters $(\log N_0,r)$, and ordinary OLS fits both with closed-form estimators. The fitted intercept is $\widehat{\log N_0}$ and the slope is $\hat{r}$. Recover ${\hat{N}}_0=e^{\widehat{\log N_0}}$. The standard errors come from the OLS theory on the log-transformed problem, and the residuals — which should be checked for normality on the log scale, not the original scale — provide diagnostics. This is much cleaner than iterative NLS on the multiplicative version, which would have to handle the heteroskedasticity of the original noise.

Take logs:

Why local minima matter

You fit a sum-of-Gaussians model $y=A_1{e}^{-(x-{\mu }_1{)}^2/2{\sigma }_1^2}+A_2{e}^{-(x-{\mu }_2{)}^2/2{\sigma }_2^2}+\epsilon$ to spectroscopic data with two visible peaks at $\mu \approx 3$ and $\mu \approx 7$.

What can go wrong with random initialization ${\hat{\mu }}_1=0,{\hat{\mu }}_2=5$ with default amplitudes, Gauss-Newton might converge to a local minimum where both Gaussians collapse onto one peak, or where one Gaussian has zero amplitude and the other tries to fit both peaks. The global minimum (one Gaussian per peak) is found only by initializing ${\hat{\mu }}_1\approx 3$ and ${\hat{\mu }}_2\approx 7$. Standard practice in spectroscopy is therefore to estimate peak locations from a peak-finding pre-pass (or by visual inspection) and pass those as starting values for NLS, never to rely on default initialization.

Starting from

See Also

• 1.2. Linear Regression: the linear case, which Gauss-Newton reduces to at each iteration.
• 1.2a. Multiple Linear Regression: matrix-form OLS that the linearization solves at every step.
• 1.5. Regularization: ridge/lasso penalties, which apply to NLS just as they do to OLS.
• 1.6. Cross-Validation: model selection when the noise distribution is uncertain.
• 2.4 Time Series Regression: the sinusoidal NLS model and Nyquist-limited frequency estimation as a time-series application.