Least squares is a projection

Too many equations, not enough unknowns

Fit a line to seven points and you have written down seven equations in two unknowns, the slope and the intercept. Unless the points happen to be exactly collinear, no line passes through all of them at once, so the system is inconsistent: b is not in the column space of A. The same thing happens in a camera pose problem with hundreds of measurements and six unknowns, or in a regression with a million rows and a handful of features. Inconsistent is the normal case, not the exception.

The move that rescues it is to stop insisting on equality and start minimising a distance. Among all the vectors Ax the matrix can produce, find the one nearest to b. That single reframing is called least squares, and it is the workhorse behind fitting, calibration and the loss functions of machine learning.

It helps to name the two regimes. A system is underdetermined when you have fewer equations than unknowns; then there are usually infinitely many solutions and the interesting question is which one is smallest — the least-norm problem of Part 18. A system is overdetermined when you have more equations than unknowns, which is the case here. Then the columns of A are too few to span the whole of the measurement space, so a generic right-hand side b lands outside the column space and there is no solution at all. Overdetermined systems are what data produces; the whole art is deciding what to do when equality is impossible.

Least squares is the standard answer because the squared distance is the one error measure that is both geometric and convenient. Geometrically, it is ordinary Euclidean distance in the measurement space, so “closest” has its everyday meaning. Algebraically, squaring makes the error a smooth quadratic function of the unknowns, which is what lets a single linear system produce the optimum instead of a search. Minimising the sum of absolute gaps would give a perfectly sensible answer too, but that function has corners and no closed form; the sum of squares has a bowl, and a bowl has a bottom you can compute.

One more framing before the pictures. Think of A as a machine that consumes a small vector x and produces a big vector Ax. The reachable outputs — the whole image of the machine — form a subspace of the measurement space. The data b is a point in that same space, usually outside the reachable set. Least squares asks for the nearest reachable point, and the answer will turn out to be the perpendicular shadow of b onto the reachable set.

The name comes from the method of least squares that Legendre and Gauss published around the turn of the nineteenth century for fitting orbits to imperfect astronomical observations. Their problem had exactly this shape: many noisy angle measurements, a handful of orbital parameters, and no way to satisfy all the measurements at once. The solution they popularised — minimise the sum of squared errors — is still the first tool reached for whenever a model has more data than freedom. Only later did it get the clean subspace interpretation this part is built around.

Drag the data, watch the fit follow

Start with the picture everyone meets first. Each data point is a pair (xᵢ, yᵢ), and a candidate line y = mx + c predicts mxᵢ + c at the point's x. The prediction is usually wrong; the gap between the measured yᵢ and the predicted value is the residual, drawn below as a dashed vertical segment. The sum of squared residuals is exactly that list of gaps, each squared and added up.

For a line fit the matrix has two columns — one holding the x-values and one holding all ones for the intercept — and the unknown vector is the pair (m, c):

Read that matrix by columns and the geometry clicks. The first column is the vector of x-coordinates; the second is the vector of all ones. A candidate line's prediction vector is m times the first column plus c times the second — a linear combination of just two directions. So the predictions reachable by some line form a two-dimensional plane sitting inside the n-dimensional measurement space. The data vector b is a single point in that space, and it is almost never on the plane. Least squares is the search for the point of the plane nearest to it.

The fitted intercept has a meaning worth stating: the least-squares line always passes through the centre of mass (x̄, ȳ) of the data. That is not a coincidence but a consequence of the second column being all ones, and it is why restricting the slope slider in the last demo to lines through the centroid loses no least-squares solution.

The least-squares solve lstsq(A, b) hands back the (m, c) that makes the total squared error as small as possible. Drag any handle and watch the line, the residual segments and the running sum all update together. Because the handles are seeded, a reload shows the same starting cloud.

Notice what the fit does not do: it does not pass through the points, and it cannot, because no line can. It chooses the compromise that leaves the smallest possible total of squared vertical gaps. The next two sections explain why that compromise is a projection, and why the algebra that computes it always looks like AᵀA x = Aᵀb.

One experiment makes the minimisation tangible before the algebra arrives. Drag a single point far from the others and watch two things happen at once: the fitted line tilts toward the outlier, and the sum of squares jumps because that point's residual is squared. A point twice as far contributes four times as much, which is why a lone outlier can dominate a least-squares fit — a real weakness of the method, and the reason robust regression replaces the square with a gentler penalty when outliers are expected.

Perpendicular to everything A can make

Here is the key geometric claim. The set of vectors Ax, as x ranges over every possible choice, is the column space of A — the span of its columns. We want the point in that space closest to b. From Part 9, the closest point in a subspace is the orthogonal projection, and the error vector connecting b to its projection is perpendicular to the subspace. So the best Ax satisfies a clean condition: the residual b − Ax is orthogonal to every vector in the column space, which means orthogonal to every column of A.

That is the whole derivation. The right-hand equation is the normal equation: a small, square, n×n system for the unknown coefficients, where n is the number of columns of A (two, for a line). The left-hand condition is the geometric content — every column dotted with the residual gives zero. When the columns are independent, AᵀA is invertible and the solution is unique.

The projection itself can be written as a matrix once you solve that system. Substituting x = (AᵀA)⁻¹Aᵀb into Ax gives

The matrix P sends any vector straight onto the column space and leaves anything already in the column space untouched. That second property is the algebraic fingerprint of a projection, P² = P: projecting a vector that is already in the subspace does nothing. It is worth checking on the demo above with a single column, where P collapses to the rank-one matrix aaᵀ / a·a.

Seen through the four subspaces of Part 7, the normal equation says something tidy about each side. The vector Aᵀb lives in the row space, and the vector Aᵀ(Ax) lives in the row space too; the equation sets them equal. Meanwhile Aᵀ(b − Ax) = 0 says the residual lies in the left null space, the orthogonal complement of the column space. The residual and the fitted values are therefore perpendicular by construction, and the measurement vector splits into a part the model can explain plus a part it cannot.

That split obeys Pythagoras. Because Ax and b − Ax are perpendicular, their squared lengths add:

Minimising the residual is the same as maximising the length of the projection, since the left-hand side is fixed by the data. The fitted vector is the largest piece of b the columns can reproduce, and the residual is the irreducibly orthogonal remainder. Replacing the flat line with a richer subspace can only make that remainder shorter, which is the linear-algebra reason a bigger model fits training data at least as well — and the reason it can also overfit.

This is also why the numerically careful way to solve it is by QR rather than by forming AᵀA explicitly: squaring the matrix squares its condition number, a theme that returns in Part 19. Geometrically, though, both routes compute the same projection.

One column, one perpendicular residual

The cleanest place to see the projection is in two dimensions, with a single column. Let a be the one column of A. Its column space is the line through the origin in the direction of a, shaded below. The vector b points somewhere off that line; the closest point on the line is the shadow of b cast perpendicular to it, and that shadow is Ax for a single number x.

Projecting onto a line is one dot product away. The shadow is (a·b / a·a) a, so the coefficient is x = a·b / a·a. Written in matrix form, A is the 2×1 column matrix, AᵀA is the 1×1 matrix holding a·a, and Aᵀb is 1×1 holding a·b — the normal equation exactly. Drag a and b and watch the residual stay perpendicular to the line.

The readout checks the orthogonality numerically: a·(b − Ax) stays at zero to within rounding, no matter where you drag. That is the normal equation in action. The projection Ax is the best approximation to b that this column space can offer, and the leftover residual is pure error the model cannot explain.

Two edge cases are worth dragging into view. If you pull b until it lies exactly on the shaded line, the projection equals b and the residual shrinks to zero — the system suddenly has an exact solution, and least squares agrees with ordinary solving. If you pull a — the column — toward the origin, the line's direction becomes unstable and a small wobble in a swings the shadow a long way; that sensitivity is the conditioning story of Part 19 in miniature.

The picture generalises wordlessly to higher dimensions. Replace the line with a plane (two independent columns), or a three-dimensional subspace (three columns), and the same statement holds: the best Ax is the perpendicular foot of b on the subspace, and the residual is perpendicular to the whole subspace, not merely to one column. Orthogonality to every column is the same as orthogonality to every linear combination of them, which is why testing against the columns is enough.

There is a second reading of the same drawing that returns to the data. Each of the n measurement rows is a point, and the line we fit is a candidate through them. The residual panel above and this projection panel are the same computation seen at two scales: the vertical gaps are the components of the residual vector, and the projection is the reason those particular gaps, and not some other set, are the smallest possible in total.

The residual sum is a bowl with one low point

Projection language is elegant, but it is worth seeing the minimisation directly. Restrict attention to lines that pass through the data's centre of mass, so the intercept is pinned and the only free choice is the slope s. Then the sum of squared residuals becomes a function of that single number:

Each term is a square, so E(s) is a parabola — a smooth bowl with exactly one lowest point. The slope at that bottom is the least-squares slope, and it is precisely what lstsq returns. Drag the slider to move the trial slope left and right; the ball on the curve shows the current error, and the marked point is the normal-equation solution. Every other slope sits higher, which is the concrete meaning of “least” in least squares.

Two facts fall out of the bowl picture. First, the minimum is unique whenever there is genuine spread in the x-values; if every x is identical, the curve flattens and the slope is undetermined, the singularity that Part 18 resolves with the pseudoinverse. Second, the curve is steep when the model is badly wrong and shallow near the bottom, which is why optimisation methods in Part 21 slow down as they approach a good fit.

The calculus version of the same story sets the derivative of E to zero and lands back on the normal equation. Differentiating a sum of squares brings down a factor of two and leaves a sum of residuals times the direction each residual moves when the slope changes; setting that to zero says the residual is orthogonal to the column of slope-influences, which is exactly one row of Aᵀ(b − Ax) = 0. The geometry and the derivative are not two methods but one statement wearing different notation.

This bowl is also the simplest member of the family of loss surfaces that modern optimisation walks down. A linear model has a convex quadratic loss with a single global minimum and a closed-form solution; a neural network has a non-convex surface with many directions and no formula, so it uses gradients and steps instead. Yet the local behaviour near an optimum is still bowl-shaped — a quadratic approximation — and that is why the normal equation keeps reappearing inside methods like Gauss–Newton. Understand the exact bowl here and you understand the shape every solver is trying to imitate.

Do not read the residual as failure. In a measurement setting the residual is where the noise lives: the part of the data the model's two parameters cannot reproduce. A residual that looks like structure — a curve, a cluster, a trend — is a signal that the model is too simple, and you should add a column for the missing effect. A residual that looks like scattered fuzz is the model doing its job. This habit of reading the leftover is one of the most useful things least squares gives you, and it is why calibration and regression reports always plot residuals rather than only quoting a single error number.

One projection, two worlds

The same three ingredients — a design matrix, a noisy target, and a squared-error minimisation — appear under many names. In robotics the matrix encodes a measurement model and the unknown is a physical quantity to be recovered. In machine learning the matrix is a layer's features and the unknown is a set of weights. In both cases the normal equation is the ideal, and the practical question is only how to solve it at the scale and conditioning of the problem at hand.

These sources cover the same material from three directions: the geometric intuition, the classroom derivation, and the numerical practice. Read the first for the picture, the second for the algebra, and the third to see why practitioners rarely form AᵀA by hand.

• Grant Sanderson, "Essence of Linear Algebra", 3Blue1Brown — the dot-product and projection chapters that make the least-squares picture geometric.
• Gilbert Strang, 18.06 Linear Algebra, MIT OpenCourseWare — Lectures 16 and 17 derive the normal equation and the projection matrix.
• Immersive Math, Chapter 4: Solving systems of linear equations — the inconsistent case in a draggable textbook.
• Sheldon Axler, Linear Algebra Done Right, the chapter on orthogonal projection — the subspace-distance viewpoint done carefully.