Nonlinear least squares

The formula

The objective, same as every least-squares problem:

S(θ) = Σ rᵢ² ,   rᵢ = yᵢ − f(xᵢ; θ)

Read: pick the parameter vector θ that minimizes the summed squared gaps between what you measured and what the model f predicts. The difference from linear regression is that θ sits inside f nonlinearly — a time constant in an exponent, a center frequency in a denominator — so no single matrix solve produces the answer.

The Jacobian, the fit's local sensitivity matrix:

Jᵢⱼ = ∂f(xᵢ; θ) / ∂θⱼ

Read: row i says how data point i responds to a small nudge of each parameter. Everything the iteration knows about the model, it knows through J.

The Gauss–Newton step ^[1]:

(Jᵀ·J)·Δθ = Jᵀ·r ,   then   θ ← θ + Δθ

Read: linearize the model at the current guess, solve the resulting linear least-squares problem for a correction, apply it, repeat. Nonlinear least squares is ordinary least squares run in a loop.

The Levenberg–Marquardt step ^[2]^[3]^[4]:

(Jᵀ·J + λ·diag(Jᵀ·J))·Δθ = Jᵀ·r

Read: add an adjustable damping λ to the diagonal. With λ near zero this is Gauss–Newton; with λ large it becomes a short, cautious step down the gradient. Shrink λ after a step that reduced S, grow it after a step that didn't. This is the algorithm inside nearly every "Fit" button you have ever pressed ^[7].

Parameter uncertainty at the solution:

Cov(θ̂) ≈ s²·(Jᵀ·J)⁻¹ ,   s² = S(θ̂) / (n − p)

Read: the same covariance formula as linear regression, evaluated on the linearization at the converged answer — which means the quoted standard errors are only as good as that local linear approximation.

Where you meet it

Time-constant extraction on the bench. A thermocouple pulled from a bath decays toward ambient as T(t) = A·e^(−t/τ) + T_amb. Fit 25 samples over 60 seconds (synthetic data, truth A = 80 °C, τ = 12 s, ambient 25 °C, noise 0.5 °C) and LM returns A = 80.35 ± 0.35, τ = 11.80 ± 0.12 s, T_amb = 25.20 ± 0.17 with residual scatter s ≈ 0.46 °C. The ± 0.12 s is what goes in the sensor spec, and it came from the covariance formula above, not from squinting at the plot.
Q from a VNA trace. A resonator's |S21| versus frequency is a Lorentzian; center frequency, peak level, and loaded Q enter it nonlinearly. The marker-based 3 dB method uses three points of the trace — a least-squares fit uses all 1601, which is why the fit's Q repeats run-to-run and the marker Q doesn't. Equivalent-circuit extraction from a full S-parameter sweep is the same job with more parameters.
Acceleration models at the reliability review. Arrhenius rate versus 1/T with unknown activation energy, Weibull life parameters, creep laws — the review board wants the parameter and its confidence interval, and both come out of the same (JᵀJ)⁻¹.
FRF fits on the shaker table. Pulling natural frequency and damping ratio out of a measured frequency response — resonance frequency in the denominator, damping in a squared term — is nonlinear least squares whether the software calls it circle fitting, rational-fraction fitting, or curve fitting.

How it works

The core move is linearize-and-iterate. Gauss–Newton is Newton's optimization method with the expensive second-derivative term of the Hessian thrown away — legitimate exactly when residuals are small, because that term is proportional to the residuals ^[1]. Consequence worth knowing: on a model that genuinely fits, convergence near the answer is nearly quadratic (my thermocouple example converges in 6–8 iterations from a rough guess); on a model that fits poorly, the discarded term isn't small and even the final iterations crawl.

Starting values are not optional. The iteration converges to the valley you started in, and nonlinear objectives can have several ^[7]. The exponential example is instructive: started at (A, τ, T_amb) = (50, 3, 50) — every value wrong by 2× or worse — Gauss–Newton lands on the right answer in 6 iterations. Started with τ guessed too long (200 s instead of 12), it diverges outright, and even LM slides onto a plateau where τ runs toward infinity: for huge τ the exponential looks like a straight line, the data can't tell τ = 10⁴ from τ = 10⁹, and JᵀJ goes numerically singular (condition number ~10¹⁸ on that path, versus ~10² near the true answer). Get starting values from physics: eyeball the asymptote for the offset, the 63% point for τ, the peak location for a center frequency. Software converging to a wrong-but-converged answer is a quieter failure than divergence, and the fitted-curve overlay plot is the only reliable detector.

What λ actually buys you. Marquardt's observation was geometric: the Gauss–Newton step and the steepest-descent step typically point nearly 90° apart, and each fails in the direction the other handles well ^[3]^[4]. The damping term steers continuously between them — and because it also fattens the diagonal, it regularizes the near-singular JᵀJ that ill-conditioned problems produce. LM buys robustness against a bad step, not against a bad basin. It will not rescue a start in the wrong valley.

The log-linearization shortcut is quietly wrong. The classic dodge for y = A·e^(−t/τ) is to take logs and run a straight-line fit on ln y versus t. That fit minimizes squared errors in log space, which weights the noisy tail of the decay as heavily as the clean early samples: with constant 2-unit noise, a point at y = 100 has log-scale scatter of 0.02 while a point at y = 4 has 0.5. Over 400 synthetic datasets (truth τ = 10), the log-line route averaged τ = 10.55 with scatter ±1.06; true nonlinear least squares averaged τ = 9.997 ± 0.22 — the shortcut is biased and nearly five times noisier for these conditions. Log-linearizing is fine for generating a starting value. It is not the fit.

Some parameter sets can't be separated no matter the algorithm. A sum of two exponentials with similar time constants, or amplitude-plus-offset when the record never reaches the asymptote, gives near-parallel Jacobian columns — the multi-parameter version of collinearity. Symptoms: standard errors bigger than the parameters, estimates that swing with one added point, λ that never shrinks. That is the data's confession that it doesn't contain the answer; the fix is a longer record or a reparameterized model, not more iterations.

The error bars are linearized error bars. s²(JᵀJ)⁻¹ trusts the tangent plane at the solution. For strongly curved problems or short records the true confidence region is banana-shaped and asymmetric, and the symmetric ± understates one side. When a parameter drives a safety margin, profile the fit — re-run with that parameter pinned across a range and watch S — instead of trusting the quoted ± blindly. And nonlinear fits inherit least squares' soft spot for outliers with fewer diagnostic tools to catch them ^[7]: plot the residuals every time.

History

The linearize-and-iterate idea is as old as least squares itself: Gauss's 1809 Theoria Motus, the book that put orbit determination on a least-squares footing, already contains the method now called Gauss–Newton ^[1]^[8]. What Gauss did not have was a machine that could iterate fast enough to expose the method's failure modes.

That machine arrived in the 1940s and 1950s, and the fix came from industry, twice. Kenneth Levenberg, working at the Army's Frankford Arsenal, published the damping idea in 1944 in the Quarterly of Applied Mathematics — add λ to the diagonal so a wild linearized step becomes a cautious gradient step ^[2]^[5]. The paper sat mostly unread. Through the late 1950s at least three others independently reinvented essentially the same trick ^[2]^[5].

The version everyone runs today belongs to Donald Marquardt, a statistician in DuPont's Engineering Department in Wilmington. By his own account he had been fighting nonlinear fits since the early 1950s on some of the first industrial digital computers, watching the textbook Taylor-series and gradient methods fail "for different reasons" — and noticing, first by plotting, that the two methods' correction vectors sat nearly at right angles to each other ^[4]. His 1963 SIAM paper framed the damped step as an optimum interpolation between them and added the scaling that makes the step insensitive to parameter units ^[2]^[3]^[4]. The problems that drove him were bench problems — his prior paper had fitted electron paramagnetic resonance spectra — and the method's practical calling card was exactly what a test engineer wants: it converges from rougher starting guesses than its predecessors ^[4]. By 1977 the paper had passed 995 citations ^[4]; it remains one of the most-used algorithms in engineering software, sitting behind scipy.optimize.curve_fit, MATLAB's lsqcurvefit, and the fit routine of nearly every instrument vendor ^[6].

Related tools

/tools/rc-charge-time — the single-exponential model whose τ this method extracts from real data
/tools/half-life-decay — same exponential in decay clothing; fitting logged counts is the log-linearization trap above
/tools/q-factor-bandwidth — the 3-point marker version of the Q a Lorentzian fit gets from the whole trace
/tools/vswr-return-loss — the S-parameter quantities that equivalent-circuit fits are run against
/tools/arrhenius-rate — activation-energy models the reliability board wants fitted with confidence intervals
/tools/vibration-natural-freq — the natural frequency an FRF curve fit estimates from shaker data

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →