The Least-Squares Solution

You point a telescope at the sky on five different nights and jot down where a comet appears. You only need three numbers to pin down its path — but you have five noisy measurements, and no path on Earth threads all five dots exactly. The measurements disagree, because every one carries a little error. So which orbit do you report?

This is the everyday shape of an over-determined inverse problem: more equations than unknowns, and no model that satisfies them all. Writing the measurements as d = Gm, there is simply no m with Gm exactly equal to d. The honest goal is then the model that comes closest: minimise the residual r = d - Gm in the least-squares sense,

\hat m = \arg\min_m \|d - Gm\|^2.

"Least squares" is exactly that phrase read backwards: make the sum of the squared misfits — the squared length of the residual vector — as small as it can be. Big misses are punished hard (squaring a big number makes it huge), so the fit refuses to let any single point stray too far.

The picture: projection onto the column space

Here is the geometry that makes the whole method click. As m ranges over every possible model, Gm traces out the column space of G — a flat subspace (a plane, say) living inside the tall measurement space. Your data d is a point that, because of noise, floats off that plane. It is unreachable: no model lands on it.

The closest reachable point is the foot of the perpendicular from d down onto the plane — the orthogonal projection of d onto the column space of G. Drop a plumb-line from the data to the subspace; where it lands is G\hat m, and the residual d - G\hat m is that plumb-line itself — perpendicular to the subspace, hence perpendicular to every column of G. That perpendicularity is not a side-effect; it is the solution.

The normal equations

Turn the picture into algebra. "Residual ⟂ every column of G" means each column, dotted with the residual, is zero — and stacking those dot products is exactly G^{\mathsf T}(d - Gm) = 0. (The transpose is what turns columns into the rows that do the dotting.) Rearranged, that is the normal equations:

G^{\mathsf T}G\,\hat m = G^{\mathsf T}d \quad\Longrightarrow\quad \hat m = (G^{\mathsf T}G)^{-1}G^{\mathsf T}d.

Notice the shapes fix themselves: if G is n\times p (tall), then G^{\mathsf T}G is a small square p\times p matrix — the over-determined problem has been squeezed down to an ordinary square system you can actually solve. This is the same closed form as the normal equation in regression: fitting a line, or training linear regression by minimising a squared-error cost, is a small over-determined inverse problem wearing different clothes. It works cleanly when G has full column rank; when it does not — or is badly conditioned — we will need the generalized inverse and regularization.

Worked example: fit a horizontal line to three readings

Simplest over-determined problem there is: estimate one constant c from three noisy readings 3, 5, 4. The model says every reading should equal c, so

G = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix},\quad d = \begin{bmatrix} 3 \\ 5 \\ 4 \end{bmatrix},\quad Gc = d \text{ is over-determined.}

The normal equations give G^{\mathsf T}G = 3 and G^{\mathsf T}d = 3+5+4 = 12, so 3\,\hat c = 12 and \hat c = 4. The least-squares estimate of a single constant is just the mean — a reassuring sanity check that the method does the obvious thing on the simplest case. The residuals (-1, 1, 0) sum to zero: exactly the "residual ⟂ the column of ones" condition.

Worked example: fit a line through four points

Now fit y = a + bx to the four points (0,1),(1,1),(2,3),(3,3). The model is

G = \begin{bmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{bmatrix},\quad m = \begin{bmatrix} a \\ b \end{bmatrix},\quad d = \begin{bmatrix} 1 \\ 1 \\ 3 \\ 3 \end{bmatrix}.

Build the two small matrices. With \sum x = 6, \sum x^2 = 14, \sum y = 8, \sum xy = 0+1+6+9 = 16,

G^{\mathsf T}G = \begin{bmatrix} 4 & 6 \\ 6 & 14 \end{bmatrix},\quad G^{\mathsf T}d = \begin{bmatrix} 8 \\ 16 \end{bmatrix}.

Solve the 2\times2 system: 4a + 6b = 8 and 6a + 14b = 16. Eliminating gives b = 0.8 and a = 0.8, so the best-fit line is y = 0.8 + 0.8x. No line hits all four points, but this one balances the four misfits so their squares add to the least possible total.

Minimise the squared residuals

Five data points and a line you control (slope and intercept). The vertical segments are the residuals; the readout is their sum of squares \|r\|^2. Hunt for the smallest value — the faint line is the true least-squares optimum (\|r\|^2 = 3.6), and you will find you cannot beat it. Watch how the total rises steeply the moment any single point is left far from the line: that steepness is the squaring doing its work.

The formula \hat m = (G^{\mathsf T}G)^{-1}G^{\mathsf T}d is beautiful on paper and treacherous in a computer. The trouble is the product G^{\mathsf T}G. A matrix's sensitivity to noise is measured by its condition number, and forming G^{\mathsf T}G roughly squares it:

\kappa(G^{\mathsf T}G) \approx \kappa(G)^2.

If G is already ill-conditioned — say \kappa(G)=10^6, not unusual — then G^{\mathsf T}G has condition number around 10^{12}, and in ordinary double precision (about 16 digits) you have thrown away most of your accuracy before you even solve. Tiny measurement noise gets blown up into a wildly wrong \hat m. This is why practitioners almost never form G^{\mathsf T}G in real code. They solve the least-squares problem instead through a QR decomposition or the SVD, both of which reach the same \hat m while keeping the condition number at \kappa(G) rather than its square. Same answer in exact arithmetic — vastly more trustworthy answer in a real machine.

On 1 January 1801 an Italian astronomer spotted a new object — the dwarf planet Ceres — then lost it in the Sun's glare after tracking it for only 40 days across a sliver of sky. Would it ever be found again? A 24-year-old Carl Friedrich Gauss took the handful of noisy telescope observations and, using his method of least squares, computed where Ceres would reappear. Months later astronomers pointed their telescopes exactly where he said — and there it was. The feat made Gauss a European celebrity almost overnight.

There was a quarrel, too: the Frenchman Adrien-Marie Legendre published least squares first, in 1805, and was furious when Gauss claimed to have used it since 1795. Historians now credit both with independent discovery. Either way, that 200-year-old idea — minimise the sum of squares — is still the single most-used data-fitting method in all of science and engineering, quietly running inside everything from spreadsheet trendlines to spacecraft navigation.

See it explained