The nearest point of a flat thing

Six measurements, two unknowns. Each measurement says the two unknowns satisfy one equation, and six equations in two unknowns have no solution unless the measurements were exact — which they never are. The system is not hard to solve. It is inconsistent, and asking a solver for the answer is asking for something that does not exist.

The repair is to change the question, and the change is geometric.

Written as $Ax = b$, the left-hand side ranges over every combination of the columns of $A$ as $x$ varies — what the matrix can reach — which is a flat subspace through the origin. Solving the system means asking whether $b$ lies in that subspace. If it does not — and with more rows than columns it almost never does — the question worth asking instead is which point of the subspace is nearest, and that is a question with an answer.

The residual has to be perpendicular #

Call the nearest point $p$ and the segment from it to $b$ the residual, $r = b - p$. The claim is that $r$ meets the subspace at a right angle, and the argument is one sentence.

Suppose it did not. Then $r$ has a non-zero shadow on some direction $u$ lying in the subspace, and sliding $p$ a little way along $u$ moves it closer to $b$ — the same computation the shadow figure makes, with the foot of the perpendicular being the nearest point of a line. So $p$ was not the nearest point after all. The nearest point is the one from which no direction in the subspace offers an improvement, and that is exactly the point whose residual casts no shadow on any of them.

The generator does not take that on trust. It computes the projection, checks that the residual is perpendicular to each column, and then prices six neighbouring points of the plane and requires every one of them to be further from the target — because “nearest” is a claim about every other point and an assertion about one point is not the same statement.

Why perpendicular means nearest, with the arithmetic #

The geometric argument above is the honest one and there is a two-line version that is worth having because it says how much worse every other point is.

Let $q$ be any other point of the subspace. Then $b - q = (b - p) + (p - q)$, and the second bracket lies in the subspace while the first is perpendicular to it. So the two are perpendicular to each other, and

$$|b - q|^2 = |b - p|^2 + |p - q|^2 .$$

That is Pythagoras, and it settles everything at once. The distance from $b$ to any other point of the subspace exceeds the distance to $p$ by exactly $|p-q|^2$, so $p$ is nearest, it is strictly nearest, and it is the only nearest point. Existence, optimality and uniqueness fall out of one right angle.

The same identity read the other way says $|b|^2 = |p|^2 + |r|^2$: what the model explains and what it does not are perpendicular, and their squared sizes add. That decomposition is the whole of what an analysis of variance is doing, and the reason the pieces add without a cross-term is the right angle and nothing else.

Where the normal equations come from #

Perpendicular to the subspace means perpendicular to each of the columns, and that is a set of equations. Writing $p = Ax$, the condition is that every column $a_j$ satisfies $a_j\cdot(b - Ax) = 0$, which collects into

$$A^{\mathsf T}A\,x = A^{\mathsf T}b .$$

These are the normal equations, and “normal” here means perpendicular rather than usual. They are square — as many equations as unknowns — and they have a solution whenever the columns are independent, so a problem with no solution has been replaced by one that always has exactly one.

Nothing was minimised to get there. No derivative was taken, no calculus appeared, and the condition arrived from the geometry of a right angle. That is worth registering, because the calculus derivation is the one usually given: write the sum of squares as a function of the coefficients, differentiate, set to zero. It reaches the same equations, it takes longer, and it explains nothing about why the answer is a projection.

The line everybody draws #

The most familiar instance has two columns — a column of ones and a column of the $x$ values — so the subspace is the set of all straight lines and the nearest point is the line of best fit.

The vertical segments are the residuals, and they are the components of a single vector in six-dimensional space — the same $r$ as before, with the plane now living inside a space with one dimension per data point. Nobody draws that space, and the picture everybody does draw shows the residuals separately rather than as one arrow. The two pictures are of the same object.

Nothing about the construction restricts it to lines. Add a column of $x^2$ and the subspace becomes the set of all parabolas; the geometry does not change at all, because the argument never mentioned what the columns were.

A fit gets better as columns are added, always, and that is a defect rather than a feature. Enlarging the subspace can only bring the nearest point nearer, so a model with more parameters fits better whatever the parameters mean, and a smaller residual is therefore not evidence of anything on its own. With six points and six columns the residual is zero and the fit is exact and worthless.

Why squares, and not something else #

The choice of penalty is not forced by the data and every alternative gives a different answer, so it is worth being explicit about what squaring buys.

It is what makes the answer a projection. Sum of squares is a squared distance, and a squared distance comes from an inner product. Every consequence above — that the answer exists, that it is unique, that the residual is perpendicular, that the pieces add by Pythagoras — is a consequence of the penalty being a squared length. Choose the sum of absolute residuals instead and the answer is still defined, still sensible, and none of that machinery applies: there is no perpendicularity, the minimiser may not be unique, and finding it is a linear program rather than a linear solve.

It is differentiable. The absolute value has a corner at zero, and a penalty with corners is minimised by an argument about which corner rather than by setting a derivative to nought. A square has no corner anywhere.

It is the answer to a statistical question, if a particular assumption holds. If the errors are independent and normally distributed with equal spread, the coefficients making the observed data most likely are exactly the least-squares ones. That is a genuine justification and it is conditional on an assumption about the errors that is very often unchecked — and the assumption is not idle, since the bell shape arises precisely when an error is a sum of many small independent contributions, which is a description of some measurements and not of others.

And the cost, which is the reason the choice is sometimes wrong: squaring makes a single bad point expensive. A residual of ten contributes a hundred, so one measurement wrong by ten times the typical error dominates the total and drags the line towards itself. Minimising absolute residuals gives a fit that ignores such a point almost entirely. Which is right depends on whether the outlying point is a mistake or a measurement, and no amount of arithmetic decides that.

Nobody solves the normal equations #

The normal equations are how least squares is derived and are not how it is computed, and the reason is one number.

Forming $A^{\mathsf T}A$ squares the ratio between the largest and smallest amounts the matrix stretches by. That ratio governs how much an error in the input is magnified in the answer, so a problem where the columns are somewhat alike — which describes almost every fit with more than three or four columns — becomes a problem where they are alike squared, and digits are lost accordingly. A matrix that would have cost six digits of the answer costs twelve.

The fix is the previous rung. Orthogonalise the columns first, so that the subspace has a perpendicular basis; the projection is then a list of shadows, computed one at a time and with nothing to invert, and recovering the coefficients from them is a back-substitution through the triangular table. Same answer, half the digits lost, and the extra cost is the factor of about three the previous rung priced.

That is a general shape worth carrying: the equations that define a quantity are frequently not the equations anybody should solve for it. The normal equations are correct, they are the shortest route from the geometry to an algorithm, and they are the wrong algorithm.

Where it needs a condition #

The columns must be independent. If they are not, the subspace is smaller than the number of columns suggests, and while the nearest point is still unique, the coefficients producing it are not — a whole family of $x$ gives the same $Ax$. The normal equations then have a singular matrix and no unique solution, which is the arithmetic reporting the situation correctly.

The standard repair is to ask for the shortest such $x$ among all of them, which is a second least-squares problem sitting inside the first and picks out a unique answer. That answer is what the pseudo-inverse computes.

Nearly dependent columns are worse than dependent ones, and this is the situation practice actually meets. The coefficients are then unique and enormously sensitive: two columns that nearly coincide can be traded off against each other almost freely, so a tiny change in the data produces a large change in the coefficients while barely moving the fitted values at all. The fit is stable and the explanation is not, and a reader shown only the fitted line has no way to tell.

And the errors are assumed to be in $b$ alone. The residuals drawn are vertical, which asserts that the $x$ values are known exactly and only the $y$ values are uncertain. If both are measured, the perpendicular distance to the line is the honest penalty, the problem is no longer linear, and the answer is different — that variant has its own name and its own literature, and it is the situation in most laboratories.

What the residual is worth reading #

The residual is what the model could not account for, and the right angle says something sharp about it that is easy to state and easy to forget.

The residual carries no information any combination of the columns could have used. That is what perpendicularity to every column means: whatever is left over is invisible to the subspace, so no adjustment of the coefficients can reach any part of it. A fit does not leave behind “the small errors”; it leaves behind exactly the component of the data that the chosen columns are blind to.

So a pattern in the residuals is evidence about the columns rather than about the data. If the residuals of a straight-line fit rise, then fall, then rise again, the residual vector is perpendicular to the two columns that were supplied and is decidedly not perpendicular to a column of $x^2$ that was not — and the visible pattern is the shadow the missing column would have had. Reading residuals is therefore a search for a direction the subspace was missing, and the search is conducted by eye because nobody knows what to look for in advance.

The same reading explains why the residuals of a good fit look like noise. Noise is the case where the leftover has no consistent shadow on anything at all, which is the honest statement that the subspace was large enough. It is not a claim that the model is true; it is a claim that no direction anybody thought to draw would improve it.

There is a trap in the arithmetic here worth naming. The residuals of a fit always sum to zero when a column of ones is present — that is one of the perpendicularity conditions, not a property of the data — and they are always uncorrelated with every column, for the same reason. So checking either of those is checking that the solve was done correctly, and no more. A diagnostic that the fitting procedure guarantees is not a diagnostic, and both of these are routinely reported as though they said something about the model.

The same right angle, in a space of random quantities #

The most surprising place this picture reappears has no vectors in it.

Take the “vectors” to be random quantities with finite variance, and the inner product of two of them to be the average of their product. It satisfies the three conditions, so there is a geometry: a length is a root-mean-square, and two quantities are perpendicular when their product averages to nought.

Now ask for the best prediction of one quantity $Y$ from another $X$ — the function of $X$ minimising the average squared error. Every function of $X$ is a “vector” in this space, and the functions of $X$ form a subspace. So the best prediction is the nearest point of that subspace, and the error is perpendicular to it.

That nearest point is the conditional expectation, and every property of it that is otherwise a calculation becomes a right angle. The tower property — that averaging a conditional average gives the plain average — is the statement that projecting twice is the same as projecting once. That the error is uncorrelated with any function of $X$ is perpendicularity. That variance decomposes into explained and unexplained parts is the same Pythagoras as three sections ago.

A theorem about prediction and a theorem about the nearest point of a plane are the same theorem, and the only thing shared between the two settings is that both have an inner product. That is what the abstraction is for.

Legendre published it and Gauss said he had it first #

Legendre published the method in 1805, in an appendix on determining the orbits of comets, and the presentation is entirely modern: the criterion, the normal equations, a worked example.

Gauss published in 1809 and remarked that he had been using the method since 1795. Legendre objected, at length and in writing, and the dispute was never settled amicably. There is reasonable evidence Gauss was telling the truth — the method is visible in his unpublished work on the orbit of Ceres in 1801, which is the calculation that made his reputation, since he predicted where the newly-lost asteroid would reappear and it did.

What Gauss added in 1809 is the part that was genuinely his: the statistical justification. Legendre offered least squares as a sensible criterion and defended it as reasonable. Gauss asked which error distribution would make it the most likely answer, found that the answer is the normal one, and thereby turned a convention into a consequence — a manoeuvre that also, incidentally, is where the normal distribution’s central position in statistics begins.

The priority dispute is the famous part and the interesting part is the other one. Two people had the same method; one of them asked what would have to be true for it to be the right method, and that question was worth more than the method.

What the picture cannot show #

The plane in the hero figure is two-dimensional and the target is in three dimensions, which is the only case that can be drawn. The situation the method exists for has one dimension per measurement — six for the fitted line, thousands in practice — and the subspace has one per column. Nothing about the argument changes, and nothing about it can be seen.

The scatter figure hides the same thing in the opposite way. It draws six residuals as six separate segments, which invites reading them as six small independent errors; they are the six coordinates of one vector, and the statement being made about them is a single statement about that vector’s direction. The picture that would show this is the hero, and the picture that shows the data is the scatter, and no single drawing is both.

And no figure here shows what a fit is worth. The residual can be made as small as wanted by adding columns, so its size is not evidence; whether the model means anything is a question about where the columns came from, and that question is not visible in any picture of the answer.

The ladder from here #

Rungs above: the adjoint, which is what $A^{\mathsf T}$ is doing in the normal equations once the transpose is stripped of its coordinates — defined by moving a map across an inner product, and the reason a symmetric matrix has perpendicular eigenvectors. The pseudo-inverse and the singular value decomposition, which handle the dependent case and say exactly how sensitive the coefficients are. Weighted and generalised least squares, where measurements of different quality are given different weights and the inner product changes to match. Bessel’s inequality and Parseval’s identity, which are this projection with infinitely many perpendicular directions. And ridge regression, which deliberately declines to find the nearest point, on the grounds that the nearest point to noisy data is not where the answer is.

Changing the question until it has an answer #

The move at the centre of this essay is worth stating on its own, because it is one of the most productive in the subject and it is almost never presented as a move at all.

The equations had no solution. The response was not to look harder, and not to declare the problem ill-posed; it was to ask a different question — nearest rather than equal — chosen so that the original question’s answer, when it exists, is also the new question’s answer. A consistent system has residual zero, so nothing was lost, and an inconsistent one now has an answer where it had none.

The same manoeuvre is everywhere once it is recognised. A polynomial with no rational root gets a larger field to have a root in. A sequence with no limit gets a limsup. A function with no derivative gets a weak one, and a matrix with no inverse gets a pseudo-inverse. In every case the enlargement is designed so that the old answer survives where it existed, and the test of whether the enlargement is honest is exactly that: the new question must agree with the old one wherever the old one had something to say.