Projection Matrices and Least Squares

You have 100 data points and want to fit a straight line. A line has 2 degrees of freedom; 100 points almost never lie exactly on one. There is no exact solution — no line passes through all of them. So what does "the best line" even mean, and how do you compute it without guessing and checking? The answer turns out to be geometric: you're not solving an equation, you're casting a shadow.

The analogy: a lamp shining straight down onto a table

Imagine a point floating in the air above a table, lit from directly overhead. Its shadow on the table is the closest point on the table to the floating point — no other point on the table is nearer. That's projection: given a point outside a flat surface, its projection onto that surface is the nearest point that actually lies on it.

Least squares is exactly this, with the "table" being every combination of your predictor variables can possibly produce, and the "floating point" being your actual observed outcomes, which almost certainly can't be produced by any combination of predictors. The best fit is the shadow of your data on the table of achievable fits — and "best" here means closest in the sum-of-squared-differences sense, which is precisely what a perpendicular shadow minimises.

Writing it down

Say you have $n$ observations and $p$ predictor variables, stacked into a matrix $X$ (n rows, p columns) and a target vector $y$ ($n$ numbers). You want coefficients $\beta$ so that $X\beta$ — the predictions — are as close to $y$ as possible. "As close as possible" in the squared sense is exactly the geometric picture: $X\beta$ must be the point in the column space of $X$ (every vector reachable by combining $X$'s columns) that is nearest to $y$. Nearest, for a flat space, means the leftover error $y - X\beta$ must be perpendicular to every column of $X$:

$$X^\top (y - X\hat\beta) = 0 \quad\Longrightarrow\quad \hat\beta = (X^\top X)^{-1}X^\top y .$$

In words: $X^\top$ means "take the dot product of the error with each column of $X$" — setting all of those to zero is the algebraic statement of "the error is perpendicular to the whole fitted space." Solving that system for $\hat\beta$ gives the famous normal-equations formula.

The projection matrix $P$ packages this into one operator that goes straight from $y$ to the fitted values without stopping at $\beta$:

$$P = X(X^\top X)^{-1}X^\top, \qquad \hat y = Py .$$

In words: $P$ is a matrix that, applied to any vector, returns its shadow on the column space of $X$. Two properties make it a projection and nothing else: $P^2 = P$ (projecting an already-projected point does nothing — the shadow of a shadow is itself) and $P^\top = P$ (it's symmetric, which is the algebraic fingerprint of "perpendicular drop," not some slanted cast). The leftover $I - P$ is also a projection — onto everything perpendicular to the column space — and it produces the residuals: $e = (I-P)y$.

This explorer shows a general linear map stretching a circle into an ellipse. A projection is the special, degenerate case: it flattens the circle onto a line or plane entirely — one direction survives untouched (anything already in the column space), every perpendicular direction collapses to zero. Try pushing the transform toward a flat ellipse to see the collapsing in action.

Worked example 1: fitting a line through three points, by hand

Fit $y = \beta_0 + \beta_1 x$ through $(1,1)$, $(2,2)$, $(3,4)$. Here $X = \begin{pmatrix}1&1\\1&2\\1&3\end{pmatrix}$, $y=\begin{pmatrix}1\\2\\4\end{pmatrix}$.

Step 1. $X^\top X = \begin{pmatrix}3 & 6\\6 & 14\end{pmatrix}$ — top-left is "count of points" (3), the rest come from summing $x$ and $x^2$ ($1+2+3=6$, $1+4+9=14$).

Step 2. $X^\top y = \begin{pmatrix}7\\17\end{pmatrix}$ — sum of $y$ is $1+2+4=7$; sum of $x\cdot y$ is $1+4+12=17$.

Step 3. Invert the 2x2: $\det = 3(14)-6(6)=42-36=6$, so $(X^\top X)^{-1} = \tfrac{1}{6}\begin{pmatrix}14 & -6\\-6 & 3\end{pmatrix}$.

Step 4. $\hat\beta = \tfrac16\begin{pmatrix}14 & -6\\-6&3\end{pmatrix}\begin{pmatrix}7\\17\end{pmatrix} = \tfrac16\begin{pmatrix}14(7)-6(17)\\-6(7)+3(17)\end{pmatrix}=\tfrac16\begin{pmatrix}98-102\\-42+51\end{pmatrix}=\tfrac16\begin{pmatrix}-4\\9\end{pmatrix}=\begin{pmatrix}-0.667\\1.5\end{pmatrix}$.

So the fitted line is $\hat y = -0.667 + 1.5x$. Check: at $x=1$, $\hat y = 0.833$ (actual 1, residual $0.167$); at $x=2$, $\hat y=2.333$ (actual 2, residual $-0.333$); at $x=3$, $\hat y=3.833$ (actual 4, residual $0.167$). Multiply residuals by the $x$-column: $0.167(1) - 0.333(1) + 0.167(1) = 0.001 \approx 0$ (rounding); multiply by $x$ itself: $0.167(1)-0.333(2)+0.167(3) = 0.167-0.667+0.5=0$. The perpendicularity condition holds exactly — that's the geometry confirming the algebra.

Worked example 2: the projection matrix directly

Same $X$. Compute $P = X(X^\top X)^{-1}X^\top$. Using $(X^\top X)^{-1}=\tfrac16\begin{pmatrix}14&-6\\-6&3\end{pmatrix}$ from above, first form $X(X^\top X)^{-1}$, a 3x2 matrix, then multiply by $X^\top$ (2x3) to get a 3x3 matrix $P$. Skipping the arithmetic (it's mechanical matrix multiplication), the diagonal of $P$ comes out to $(0.833, 0.333, 0.833)$ — these are the leverage values, and they sum to exactly 2, the number of parameters ($p=2$), which is a general fact: $\text{trace}(P) = p$ always. Applying $P$ to $y=(1,2,4)$ reproduces the fitted values $(0.833, 2.333, 3.833)$ found above — same answer, computed via the operator instead of via $\beta$.

What this means in practice

• Every OLS regression is a projection. The fitted values are literally $Py$; the residuals are $(I-P)y$; and those two pieces are geometrically perpendicular, which is why $R^2$ decomposes cleanly into "explained" and "unexplained" sums of squares.
• The trace of $P$ counts effective parameters. This generalises beyond plain regression — ridge regression and smoothing splines have a "projection-like" matrix whose trace gives an effective degrees of freedom, even when it isn't exactly $p$.
• Numerically, nobody inverts $X^\top X$ directly. It's condition-number-sensitive (see the matrix condition number); production code uses QR or Cholesky factorisations of $X^\top X$ instead, which are algebraically equivalent but numerically stable.