Suppose you weigh yourself on the same bathroom scale three mornings in a row and read 61.2, 60.8, and 61.0 kg. What do you weigh? There is no single number that makes all three readings correct at once — the scale wobbled, your feet shifted, the numbers simply disagree. And yet you know, instinctively, what to do: pick the number that sits closest to all three at the same time. You've just invented least squares.

This is the everyday situation where you have too much information that doesn't quite fit together: more measurements than unknowns, all slightly at odds, no way to please every one exactly. Least squares is the beautiful, once-and-for-all answer to "so what's the best I can do?" — and, remarkably, it's built entirely from projection and orthogonality, the two ideas you already have in hand.

When a system has no solution

A linear system A\vec{x} = \vec{b} is overdetermined when A has more rows than columns — more equations than unknowns. Geometrically, the products A\vec{x}, as \vec{x} ranges over everything, sweep out only the column space of A — the set of all combinations of A's columns. If the target \vec{b} happens to lie outside that column space, then no \vec{x} hits it, and the system has no exact solution at all.

That's the normal state of affairs for real data. Three noisy weighings, a thousand (x,y) points that "should" lie on a line but don't quite — you have far more equations than unknowns, and asking for a perfect fit is asking for the impossible. So we change the question. Instead of demanding A\vec{x} = \vec{b} exactly, we hunt for the \hat{\vec{x}} that makes A\hat{\vec{x}} land as close to \vec{b} as anything possibly can:

The quantity \vec{r} = \vec{b} - A\vec{x} is the residual — the leftover error — and \lVert A\vec{x} - \vec{b}\rVert^2 is the sum of the squared residuals. Least squares is the name of the game: make that squared error as small as it will go.

The picture: project \vec{b} onto the column space

Here is the whole idea in one sentence. Of all the points you can reach — every A\vec{x} in the column space — the one closest to \vec{b} is the orthogonal projection of \vec{b} onto that space. Closest point on a subspace to an outside point? That is exactly what projection computes.

Think of the column space as a flat plane and \vec{b} as a point floating above it. The nearest spot on the plane sits directly beneath \vec{b}, and the little rod connecting them stands straight up — perpendicular to the plane. Any other point on the plane is reached by a slanted, and therefore longer, rod. So the best A\hat{\vec{x}} is the foot of that perpendicular, and the residual \vec{b} - A\hat{\vec{x}} is orthogonal to the entire column space — it points straight off the plane.

Being orthogonal to the whole column space means being orthogonal to each column of A individually. And "orthogonal to every column of A" is precisely the statement that A^{\mathsf{T}} (whose rows are those columns) dotted into the residual gives all zeros:

The normal equations

Take that one perpendicularity condition and multiply it out. Distributing A^{\mathsf{T}} gives A^{\mathsf{T}}\vec{b} - A^{\mathsf{T}}A\hat{\vec{x}} = \vec{0}, and moving one term across leaves the single most important equation in all of applied linear algebra:

These are the normal equations (named for the normal, i.e. perpendicular, residual). Notice the magic: A might be a tall 1000 \times 2 slab with no inverse in sight, but A^{\mathsf{T}}A is a tidy little square 2 \times 2 matrix. As long as the columns of A are independent, A^{\mathsf{T}}A is invertible, and you just solve a small square system for \hat{\vec{x}}. Written out: \hat{\vec{x}} = (A^{\mathsf{T}}A)^{-1}A^{\mathsf{T}}\vec{b} — though in practice you solve the system rather than form that inverse.

Worked example 1 — a line of best fit through three points

Fit a straight line y = mx + c through the three points (1,1), (2,2), (3,2). They aren't collinear, so no line passes through all three — this is a job for least squares. Each point demands m x_i + c = y_i, which stacks into A\vec{x} = \vec{b} with unknowns \vec{x} = \begin{bmatrix} m \\ c \end{bmatrix}:

The first column holds the x-values (the coefficients of m); the second is all 1s (the coefficients of c). Now build the two small pieces:

Solve the 2\times2 normal equations \begin{bmatrix} 14 & 6 \\ 6 & 3 \end{bmatrix}\begin{bmatrix} m \\ c \end{bmatrix} = \begin{bmatrix} 11 \\ 5 \end{bmatrix}. The determinant is 14\cdot 3 - 6\cdot 6 = 6, so

The best-fit line is y = \tfrac12 x + \tfrac23. No line hits all three points, but this one gets provably closer — smaller total squared vertical miss — than any other line you could draw.

Worked example 2 — a system with no exact solution

Least squares isn't only for fitting lines; it rescues any overdetermined system. Take three equations in the two unknowns x and y:

The first two force x = 1 and y = 1, but then the third wants x + y = 3 while it's actually 2. Contradiction — there is no exact solution. In matrix form,

Solving \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}\hat{\vec{x}} = \begin{bmatrix} 4 \\ 4 \end{bmatrix} (determinant 3) gives, by the symmetry of the problem, \hat{x} = \hat{y} = \tfrac43. This isn't a "real" solution — it satisfies none of the three equations exactly — but it's the compromise that's least wrong overall. To confirm the geometry, form the residual and dot it against the columns:

The residual is dead orthogonal to both columns of A — the signature of a genuine least-squares solution.

Watch the line pivot to the best fit

Five data points, and the blue line is the least-squares line of best fit computed live from the normal equations. The short dashed segments are the residuals — the vertical gaps least squares is trying to shrink. Drag the slider to lift or drop the middle point and watch the line swing and the "sum of squared residuals" readout respond. Notice that the line is pulled hard toward a point that strays far off, because a big residual, once squared, dominates the total.

The residuals are measured vertically — straight up and down in the y (the \vec{b}) direction — not as the perpendicular distance from the point to the line. That's the whole reason the algebra came out as projection onto the column space: the error we minimise lives in \vec{b}'s world.

This is linear regression

Everything above is the machinery behind the "line of best fit" — what statisticians call linear regression. Given a cloud of data points (x_i, y_i), you stack the x_is and a column of 1s into A, stack the y_is into \vec{b}, and solve the normal equations for the slope and intercept. The very same procedure, with more columns in A (one per feature), fits a plane, a hyperplane, or a whole polynomial — it's the workhorse of curve-fitting, forecasting, and a great deal of machine learning.

For big problems, forming and solving A^{\mathsf{T}}A directly can be numerically shaky (it squares the "spread" of the data and amplifies rounding error). The standard professional cure is to factor A = QR into an orthogonal Q and a triangular R first, then solve — the subject of QR factorization.