Linear least squares

The formula

The setup: you have n data points and a model that is linear in its unknown coefficients. Stack the data into a matrix equation:

y ≈ X·β

Read: the vector of measurements y should be approximately the model matrix X times the coefficient vector β. Each row of X is one measurement's known values (for a straight line, a row is [1, x_i]); each entry of β is a coefficient you want (intercept, slope, and so on).

The criterion:

S = Σ (y_i − ŷ_i)²  →  minimize over β

Read: add up the squared vertical misses between each measurement and the model's prediction, and choose the coefficients that make that total smallest.

The closed-form answer — the normal equations:

Xᵀ·X·β̂ = Xᵀ·y

Read: the best-fit coefficients solve a small square linear system built from your data. For k coefficients this is a k × k solve, no matter how many thousands of points you took.

For the workhorse case — a straight line y = m·x + b — the normal equations collapse to two formulas worth memorizing:

m = Σ (x_i − x̄)·(y_i − ȳ) / Σ (x_i − x̄)²
b = ȳ − m·x̄

Read: the slope is the covariance of x and y divided by the variance of x, and the fitted line always passes through the point (x̄, ȳ) — the centroid of your data.

Where you meet it

Transducer calibration on the bench. You put five known pressures on a 4–20 mA transmitter and log the current: at 0, 5, 10, 15, 20 psi you read 4.02, 8.05, 11.96, 16.03, 20.01 mA. The least-squares line through those points is I = 0.7992·P + 4.022, and that slope and offset go straight into the data-acquisition scaling — not the two-point endpoints, because the fit averages down the noise in all five readings.
Pulling Young's modulus off a tensile test. The modulus is the slope of the elastic region of the stress–strain record. Nobody eyeballs it; the test software runs a least-squares fit over the linear segment, and the argument at the review board is about which segment, not about the fit.
Drift characterization on a test stand. A load cell or thermocouple channel logged over an eight-hour soak gives you drift in units per hour as a fitted slope, with the residuals telling you whether the drift is actually linear or something worse.
Arrhenius life projections in a reliability review. Plot log of failure rate against 1/T from accelerated-life ovens, fit a line, and the slope hands you activation energy. Every acceleration factor briefed to a customer rode in on a least-squares slope.

How it works

Squaring the residuals does two jobs at once. It makes the objective smooth and bowl-shaped, so calculus gives a unique closed-form minimum instead of an iterative search. And it makes the algebra linear: set the derivative of S with respect to each coefficient to zero and the normal equations fall out. A useful byproduct — the residuals of a fit with an intercept always sum to exactly zero, which is a free sanity check on any implementation.

"Linear" means linear in the coefficients, not in x. Fitting y = a + b·x + c·x² is still linear least squares — the columns of X are just 1, x, x². So is fitting a sum of known basis functions like a·sin(ωt) + b·cos(ωt) when ω is known. What is not linear least squares is fitting y = a·e^(b·x) directly, because b sits inside the exponential. The common workaround — take the log and fit a line — works, but understand what you traded: minimizing squared error in log y weights the small-y points far more heavily than minimizing squared error in y itself. Sometimes that is what you want. Know which one you asked for.

Squared error means outliers vote with squared strength. A point 10 units off the trend pulls on the fit 100 times harder than a point 1 unit off. One dropped sample, one bad channel, one operator bumping the stand, and your slope moves. Plot the residuals every time. A residual plot that looks like structureless noise means the model form is adequate; a bow, a trend, or one point sitting far from the rest means it is not, and no summary statistic will tell you that. A high R² will happily coexist with a systematically wrong model.

The normal equations are the textbook route, not the numerical one. Forming Xᵀ·X squares the condition number of the problem — a fitting matrix with condition number 3.2 × 10³ becomes a normal-equations matrix at about 1.05 × 10⁷, and you have burned half your available digits before solving anything. Every serious library (LAPACK, NumPy's lstsq, MATLAB's backslash) solves the least-squares problem through QR factorization or the SVD instead. High-order polynomial fits on raw x values are the classic way to hit this wall; centering and scaling x before fitting buys back most of the damage.

A fit without uncertainty numbers is half an answer. The scatter of the residuals estimates the measurement noise ^[5]:

s = sqrt( Σ r_i² / (n − 2) )          (straight-line case; n − k in general)
sigma_m = s / sqrt( Σ (x_i − x̄)² )    (standard error of the slope)

Read: the residual standard deviation s is the typical miss, and the slope's uncertainty shrinks as you spread your x values wider and take more of them. For the transmitter example above, s ≈ 0.038 mA and the slope is 0.7992 ± 0.0024 mA/psi — which is the number a review board should actually see, because it says whether the sensitivity meets its tolerance, not just what it measured. Note the Σ (x_i − x̄)² in the denominator: bunching all your calibration points in the middle of the range is the fastest way to buy an expensive, uncertain slope.

Know the assumptions you are leaning on. Ordinary least squares treats x as exact and puts all the error in y. If your independent variable is also a noisy measurement — common when both axes come off instrumentation — the fitted slope is biased low, and you are in errors-in-variables territory (total least squares). It also assumes the errors have roughly constant variance; if scatter grows with signal level, weighted least squares (divide each residual by its expected sigma) is the honest version. And a fit is an interpolation license, not an extrapolation license: the polynomial that hugs your data inside the tested range can do anything it pleases outside it.

History

The method appeared in print in 1805, in a nine-page appendix Legendre attached to a book on computing comet orbits — four pages of method, five pages applying it to the shape of the Earth and the definition of the metre ^[1]^[4]^[6]. Legendre stated the idea plainly: among all possible coefficient choices, take the one minimizing the sum of squares of the errors. The technique spread through European astronomy and geodesy within a decade ^[1].

Then Gauss published his own treatment in 1809, in Theoria Motus, his treatise on celestial mechanics — and remarked that he had been using the principle since 1795 ^[1]^[3]. Legendre took it about as well as you would expect a man watching his best-known result get back-dated by a rival, and spent years defending his priority ^[2]. The dispute never fully resolved, but the accounting most historians settle on is: Legendre published first and clearly; Gauss almost certainly did have it earlier and, more importantly, built the theory under it — connecting least squares to probability and the normal error curve, which Legendre never did ^[1]^[2].

Gauss also had the better demo. When the newly discovered asteroid Ceres vanished behind the Sun in 1801, the 24-year-old Gauss computed where it would reappear. His prediction disagreed with everyone else's; on December 7, 1801, von Zach recovered Ceres almost exactly where Gauss said to look ^[1]^[3]. By the early 1820s Gauss had proven the result now taught as the Gauss–Markov theorem: for a linear model with zero-mean, uncorrelated, equal-variance errors, least squares gives the minimum-variance linear unbiased estimate ^[1]^[7]. Two centuries on, it is still the default first tool applied to nearly every calibration curve and test log in the business.

Related tools

/tools/loop-4-20ma — the slope-and-offset scaling that a transmitter calibration fit feeds
/tools/stress-strain — Young's modulus is a least-squares slope through the elastic region
/tools/strain-gauge-bridge — bridge output vs. load is calibrated with a linear fit
/tools/arrhenius-rate — activation energy comes from a line fit on log-rate vs. 1/T
/tools/snr-enob — ADC sine-fit testing per IEEE 1057 is a least-squares problem under the hood

Sources

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