Least squares: the line that argues with every point

Give someone a scatter of dots and a ruler and ask them to draw "the line through the data." They'll manage it — a person can eyeball a decent fit in about a second, and two people will draw nearly the same line. That's the strange part. Nobody agreed on a definition of best, yet everyone's hand converges on roughly the same answer. Least squares is what happens when you insist on writing that shared intuition down as an equation, and then follow the equation exactly where it leads. It leads somewhere lovely.

The method is old and slightly contested — Legendre published it in 1805, Gauss claimed he'd been using it since 1795, and the resulting priority spat was as bitter as mathematicians get. It has aged better than almost anything else in this series: every linear regression, every trend line in a spreadsheet, every "line of best fit" in a physics lab is this one idea. Let's earn it.

What "best" has to mean

We have points (x_1, y_1), \ldots, (x_n, y_n) and we want a line \hat y = ax + b. For each point the line makes a prediction \hat y_i = a x_i + b, and misses by a residual

r_i = y_i - (a x_i + b).

A good line makes the residuals small. But "small" for a whole collection of numbers needs a single scalar to minimize, and the choice of scalar is the whole game. Minimize the sum of absolute residuals \sum |r_i| and you get a different, perfectly respectable line (it's called least absolute deviations, and it shrugs off outliers). Minimize the largest residual \max_i |r_i| and you get yet another. Least squares picks

L(a, b) = \sum_{i=1}^n r_i^2 = \sum_{i=1}^n \big(y_i - a x_i - b\big)^2,

the sum of squared residuals. Why squares? Three answers, in increasing order of honesty.

The cheap answer: squaring kills the sign, so positive and negative misses can't cancel, and it's smooth — no ugly corner like |r| has at zero, which means we can differentiate and set the derivative to zero. Convenient.

The better answer is a picture we'll get to: squared error is Euclidean distance, and minimizing it turns the whole problem into geometry.

The honest answer is a statistical one. Suppose each y_i is a true linear response plus independent noise \varepsilon_i \sim \mathcal{N}(0, \sigma^2). The probability of seeing the data you saw, as a function of the parameters, is the product of Gaussian densities:

\mathcal{L}(a,b) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi}\,\sigma} \exp\!\left(-\frac{(y_i - a x_i - b)^2}{2\sigma^2}\right).

Take the logarithm — it's monotone, so it doesn't move the maximum — and everything but one term is a constant:

\log \mathcal{L}(a,b) = -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - a x_i - b)^2 + \text{const}.

Maximizing the likelihood is exactly minimizing the sum of squared residuals. So least squares isn't an arbitrary aesthetic preference. If you believe your noise is Gaussian — and thanks to the central limit theorem, "a sum of many small independent nudges" is Gaussian more often than you'd think — then the least-squares line is the maximum-likelihood line. Squaring is what a bell curve looks like once you take its log.

The arena: matrices

Two parameters is a warm-up. The real version has many features, so we stack everything. Write the parameters as a vector \beta, put the data in a matrix X (one row per point, one column per feature, with a column of ones absorbing the intercept), and the whole model is one line:

\hat y = X\beta, \qquad L(\beta) = \lVert X\beta - y \rVert^2.

That \lVert \cdot \rVert^2 is the squared Euclidean length of the residual vector r = X\beta - y — the sum of its squared components, which is exactly \sum r_i^2. Everything from here is about one question: which \beta makes that length smallest? We answer it twice.

Derivation one: set the gradient to zero

Expand the squared norm using \lVert v \rVert^2 = v^\top v:

L(\beta) = (X\beta - y)^\top (X\beta - y) = \beta^\top X^\top X \beta - 2\beta^\top X^\top y + y^\top y.

This is a quadratic in \beta — a bowl. Differentiate with respect to \beta (the vector calculus rules \nabla_\beta\, \beta^\top A \beta = 2A\beta for symmetric A, and \nabla_\beta\, \beta^\top c = c, are all you need) and set it to zero:

\nabla_\beta L = 2 X^\top X \beta - 2 X^\top y = 0.

Rearrange, and out fall the normal equations:

\boxed{\,X^\top X\, \beta = X^\top y\,}

When X^\top X is invertible, the solution is \hat\beta = (X^\top X)^{-1} X^\top y. And because the Hessian is 2X^\top X, which is positive semidefinite (for any v, v^\top X^\top X v = \lVert Xv\rVert^2 \ge 0), the bowl only ever curves upward — the stationary point is a genuine minimum, not a saddle or a max. One equation, guaranteed to be the best. That's the calculus story, and it's complete. But it doesn't tell you why the answer looks the way it does. For that, drop the algebra and look.

Derivation two: the answer is a shadow

Here's the reframing that makes least squares feel inevitable. The vector X\beta — as \beta ranges over all possible parameter vectors — traces out every linear combination of the columns of X. That set is a flat subspace: the column space of X, written \mathrm{col}(X). Our target y almost certainly does not lie in that subspace (if it did, the fit would be perfect and there'd be nothing to minimize). So the question becomes purely geometric:

Of all the points in the subspace \mathrm{col}(X), which one is closest to y?

And you already know the answer, because you know it about shadows. The closest point in a plane to a point floating above it is the foot of the perpendicular — the spot directly below, where a dropped plumb line meets the floor. The closest point is the orthogonal projection of y onto the subspace.

Call that closest point \hat y = X\hat\beta. "Closest" means the residual y - \hat y is perpendicular to the subspace — perpendicular to every column of X at once. Perpendicular means dot products vanish, and "dot product with every column" is exactly what X^\top computes:

X^\top (y - X\hat\beta) = 0 \quad\Longrightarrow\quad X^\top X \hat\beta = X^\top y.

The same normal equations, derived without a single derivative — just the geometric fact that the shortest way to a plane is straight down. The word "normal" in "normal equations" is not a synonym for ordinary; it means perpendicular. The equations are the statement that the residual is normal to the column space.

Multiply through and you can write the fit as \hat y = X(X^\top X)^{-1} X^\top y = P y, where P = X(X^\top X)^{-1} X^\top is the projection matrix (the "hat matrix" — it puts the hat on y). Project twice and nothing changes, because you're already in the subspace: P^2 = P. That single algebraic identity is the whole geometry, compressed.

Where it hurts: outliers

Squaring the residuals had a cost, and it's worth staring at before you trust the method blindly. A point twice as far from the line contributes four times as much to the loss. A point ten times as far contributes a hundred times as much. Least squares is therefore easy to bully: one wild outlier, sitting far off the trend, can grab the line and drag it noticeably off course, because bending toward that one point buys a huge reduction in its enormous squared residual — at the mere cost of small increases elsewhere.

This is the flip side of the Gaussian story. A bell curve has thin tails; it thinks large deviations are nearly impossible. Feed it a large deviation anyway and it contorts the fit rather than believe the point is an anomaly. If your noise really is heavy-tailed, least absolute deviations or a robust loss will serve you better. Knowing why least squares is fragile is what tells you when to reach for something else.

The best way to feel this is to do it. The least-squares playground lets you drag the data points around and watch the best-fit line recompute in real time, with the residual sticks and the loss surface redrawing as you go. Grab one point and haul it away from the crowd — watch the whole line swing to chase it, and watch the bottom of the loss bowl slide across the parameter plane. It's the fastest way I know to build a gut sense for both the elegance and the fragility of squaring your errors.

Where this goes

Least squares is the seed. Grow the story into a full statistical model — with assumptions, uncertainty on the coefficients, and the Gauss–Markov theorem that says this estimator is the best of its kind — and you have linear regression. Notice that the normal equations need X^\top X to be invertible and involve a matrix solve that's awkward when X is enormous, and you'll want to find the minimum of that bowl by rolling downhill instead of solving in one shot — which is gradient descent, and the same bowl-shaped loss surface you can already see moving in the playground. Same minimum, two roads to it.