Why Fitting k Parameters Leaves n − k Directions

August 2, 2026 statistics, estimation, math, linear-algebra

Ordinary least squares (OLS) fits a linear model by minimizing squared residuals. Geometrically, fitting k coefficients projects the observations onto a k-dimensional model space, leaving n − k perpendicular directions in which noise can remain. This is the same fact behind n − 1 in sample variance.

The Puzzle

In the previous post, we estimated the mean by choosing the value $\hat m$ that minimized the squared deviations:

\hat m=\arg\min_m\frac1n\sum_{i=1}^n(y_i-m)^2.

The underlying model was

y_i=\mu+\epsilon_i, \qquad \mathbb E[\epsilon_i]=0, \qquad \operatorname{Var}(\epsilon_i)=\sigma^2.

This fit has one parameter, $m$ . After fitting it, the squared residuals have expectation $(n-1)\sigma^2$ , which is why the unbiased variance estimator divides their sum by $n-1$ .

Now suppose the mean is not constant. Each observation is a pair $(x_i,y_i)$ : $x_i$ is an observed predictor value and $y_i$ is the corresponding noisy response. Throughout this post, we treat $x_1,\ldots,x_n$ as known, fixed numbers. The randomness comes from the noise in the responses.

The population mean of the response at predictor value $x_i$ is $\alpha+\beta x_i$ :

y_i=\alpha+\beta x_i+\epsilon_i, \qquad \mathbb E[\epsilon_i]=0, \qquad \operatorname{Var}(\epsilon_i)=\sigma^2.

We do exactly what we did before, except that we now fit a candidate line $a+bx$ instead of a candidate constant $m$ :

(\hat a,\hat b) = \arg\min_{a,b} \frac1n\sum_{i=1}^n(y_i-a-bx_i)^2.

This procedure is ordinary least squares (OLS). The previous post was already its one-parameter version: a constant is just a line with no slope. Here we fit two parameters, the intercept and the slope. Does estimating one additional parameter change $n-1$ to $n-2$ ?

Provided $n>2$ and the predictor is not constant, the unbiased residual-variance estimator is no longer divided by $n-1$ , but by $n-2$ :

Residual Variance in Simple Linear Regression

s^2=\frac{1}{n-2}\sum_{i=1}^n(y_i-\hat a-\hat b x_i)^2.

The optimization has barely changed. We replaced the candidate constant $m$ with the candidate line $a+bx$ . Yet one fitted coefficient became two, and the divisor lost one more. In this post, we will see a geometric proof of why a full-rank linear fit with $k$ coefficients leaves $n-k$ residual degrees of freedom.

The Model Is a Subspace

Write the $n$ observed responses as one vector $\mathbf y\in\mathbb R^n$ . In a linear model, $\mathbf X\in\mathbb R^{n\times k}$ is the design matrix and $\theta\in\mathbb R^k$ is the unknown coefficient vector.

The column space $\operatorname{col}(\mathbf X)$ is the set of all linear combinations of the columns of $\mathbf X$ . As $\vartheta\in\mathbb R^k$ varies, $\mathbf X\vartheta$ runs through exactly this set, so the column space is the set of every response vector the model is capable of fitting. If the $k$ columns are linearly independent, it is a $k$ -dimensional flat subspace of $\mathbb R^n$ .

For a line with an intercept and slope, the two columns are $\mathbf1$ and $\mathbf x$ . Choosing the intercept and slope chooses a linear combination of those two vectors:

\begin{gathered} \mathbf y=\mathbf X\theta+\boldsymbol\epsilon,\\[0.4em] \mathbf X= \begin{pmatrix} 1&x_1\\ 1&x_2\\ \vdots&\vdots\\ 1&x_n \end{pmatrix}, \qquad \theta= \begin{pmatrix} \alpha\\ \beta \end{pmatrix}. \end{gathered}

Throughout the derivation we treat $\mathbf X$ as fixed, assume full column rank with $n>k$ , and assume the mean function is correctly specified. Every observation has the same noise variance, and distinct noise terms are uncorrelated:

\mathbb E[\boldsymbol\epsilon\mid\mathbf X]=\mathbf0, \qquad \mathbb E[\boldsymbol\epsilon\boldsymbol\epsilon^\top\mid\mathbf X] =\sigma^2\mathbf I_n.

OLS Creates the Right Angle

OLS chooses the model vector closest to the observed response $\mathbf y$ . The shortest path from a point to a flat subspace meets that subspace at a right angle. We can verify that claim directly.

Let $\hat{\mathbf y}$ be the fitted point and $\mathbf r$ its residual:

\begin{gathered} \hat\theta=\arg\min_\vartheta\|\mathbf y-\mathbf X\vartheta\|^2,\\ \hat{\mathbf y}=\mathbf X\hat\theta,\\ \mathbf r=\mathbf y-\hat{\mathbf y}. \end{gathered}

The fit only has the freedom to choose a vector from $\operatorname{col}(\mathbf X)$ . The residual $\mathbf r$ is what remains of $\mathbf y$ after making the best possible choice within that space: the part the model cannot cover using its available columns.

Now take any nonzero direction $\mathbf v$ in the model space. Moving the fitted point by $t\mathbf v$ gives another allowed fitted point, $\hat{\mathbf y}+t\mathbf v$ . The vector from this new fitted point to the observed response is the new residual:

\mathbf y-(\hat{\mathbf y}+t\mathbf v)=\mathbf r-t\mathbf v.

The squared length of the new residual is the new squared distance to $\mathbf y$ . Expanding it gives

\|\mathbf r-t\mathbf v\|^2 =\|\mathbf r\|^2 -2t\,\mathbf r^\top\mathbf v +t^2\|\mathbf v\|^2.

Suppose $\mathbf r^\top\mathbf v\neq0$ . Choose $t$ to remove the component of $\mathbf r$ pointing along $\mathbf v$ . Substituting that value into the expansion gives

\begin{gathered} t=\frac{\mathbf r^\top\mathbf v}{\|\mathbf v\|^2},\\[0.4em] \|\mathbf r-t\mathbf v\|^2 =\|\mathbf r\|^2 -\frac{(\mathbf r^\top\mathbf v)^2}{\|\mathbf v\|^2} <\|\mathbf r\|^2. \end{gathered}

The left-hand side is the squared distance from $\mathbf y$ to the new fitted point $\hat{\mathbf y}+t\mathbf v$ ; the right-hand side is the squared distance from $\mathbf y$ to the OLS fit $\hat{\mathbf y}$ . The inequality says the new point would be closer, contradicting the fact that OLS chose the closest point. Therefore $\mathbf r^\top\mathbf v=0$ for every direction $\mathbf v$ in the model space: the residual is perpendicular to that space. The normal equations state the same conclusion in calculus notation:

\begin{gathered} \mathbf r\perp\operatorname{col}(\mathbf X),\\ \mathbf X^\top(\mathbf y-\mathbf X\hat\theta)=\mathbf0,\\ \mathbf X^\top\mathbf r=\mathbf0. \end{gathered}

Each entry of $\mathbf X^\top\mathbf r$ is the inner product of the residual with one column of $\mathbf X$ . All are zero: the residual is perpendicular to every fitted direction.

Recall the variance question we started with: the residual sum of squares is $\mathrm{RSS}=\|\mathbf r\|^2$ , and we want to know how much noise remains in it after fitting the model. The right-angle argument identifies the part left over, but it does not yet explain why its expected squared length is exactly $(n-k)\sigma^2$ . To obtain that factor, we need an operator that isolates the residual part and lets us count its surviving noise directions.

The Hat Matrix Performs the Projection

The right-angle argument tells us what OLS does. Solving the normal equations gives the operator that actually does it:

\begin{aligned} \mathbf X^\top\mathbf X\hat\theta&=\mathbf X^\top\mathbf y,\\ \hat\theta&=(\mathbf X^\top\mathbf X)^{-1}\mathbf X^\top\mathbf y,\\ \hat{\mathbf y} &=\underbrace{\mathbf X(\mathbf X^\top\mathbf X)^{-1}\mathbf X^\top}_{\mathbf H}\mathbf y =\mathbf H\mathbf y. \end{aligned}

$\mathbf H$ is the hat matrix, named because it puts the hat on $\mathbf y$ . It projects any response vector onto $\operatorname{col}(\mathbf X)$ . Its companion extracts the perpendicular part:

\hat{\mathbf y}=\mathbf H\mathbf y, \qquad \mathbf r=(\mathbf I-\mathbf H)\mathbf y.

The hat matrix did not create the geometry. Least squares did. What $\mathbf H$ adds is a concrete projection operator that we can calculate with.

Why Only Noise Remains in the Residual

Recall that the true model decomposes the response as $\mathbf y=\mathbf X\theta+\boldsymbol\epsilon$ . The vector $\mathbf X\theta$ contains the noise-free population mean at each observed predictor value; this is the signal. The vector $\boldsymbol\epsilon$ is the noise.

Substitute this decomposition into the residual operator:

\begin{aligned} \mathbf r &=(\mathbf I-\mathbf H)(\mathbf X\theta+\boldsymbol\epsilon)\\ &=\underbrace{(\mathbf I-\mathbf H)\mathbf X\theta}_{\mathbf0} +(\mathbf I-\mathbf H)\boldsymbol\epsilon\\ &=(\mathbf I-\mathbf H)\boldsymbol\epsilon. \end{aligned}

The mean-response term $\mathbf X\theta$ disappears because it already lies in the fitted space, where $\mathbf H$ reproduces it exactly. The residual is therefore not “the part of the data we threw away.” It is precisely the component of the noise that points outside the model space.

Before counting those directions, it is worth making this cancellation visible. With three observations, an intercept-and-slope model spans a literal plane in $\mathbb R^3$ . Holding the noise fixed, changing the true intercept ( $\alpha$ ) or slope ( $\beta$ ) moves both the generated response and its fitted component within that plane, leaving the residual unchanged. The $\delta$ slider adds a possible noise component perpendicular to the plane; because the model has no direction with which to cover it, the entire component passes into the residual.

intercept α = 0.60slope β = 0.50outside component δ = 0.80

response space for three observations

the same residual as three coordinates

‖r‖² = Σrᵢ² = 0.640

changing α or β moves the fit within the plane; only δ changes r

With three observations, the intercept and slope span the plane col(X) inside ℝ³. Changing α or β moves the fitted vector within that plane and leaves the residual unchanged. The added vector δq is perpendicular to both columns of X, so the model cannot cover it: it passes entirely into r.

This is the first half of the story: fitting removes the part of the response that lies in the model space, so the residual contains only the noise component perpendicular to that space. What remains is to count how much noise energy is expected to live there.

Counting the Noise Directions

This is the math-heavy section of the post, but the takeaway is simple: after the projection, the residual keeps exactly $n-k$ orthogonal noise coordinates, each contributing $\sigma^2$ in expectation. Everything below is just making that sentence precise.

We could continue in the original observation coordinates and expand $\mathbf r=(\mathbf I-\mathbf H)\boldsymbol\epsilon$ there. But each residual coordinate would mix several noise coordinates, and squaring them would produce many cross-terms. The geometry would be hidden inside the algebra.

To see the difficulty, write $\mathbf M=\mathbf I-\mathbf H$ . Then

\begin{aligned} r_i &=\sum_{\ell=1}^n M_{i\ell}\epsilon_\ell,\\ \|\mathbf r\|^2 &=\sum_{i=1}^n \left(\sum_{\ell=1}^n M_{i\ell}\epsilon_\ell\right)^2\\ &=\sum_{i=1}^n\sum_{\ell=1}^n\sum_{m=1}^n M_{i\ell}M_{im}\epsilon_\ell\epsilon_m,\\[0.4em] \mathbb E[\|\mathbf r\|^2\mid\mathbf X] &=\sigma^2\sum_{i=1}^n\sum_{\ell=1}^n M_{i\ell}^2. \end{aligned}

The cross-terms disappear under white noise, but the last double sum still does not visibly say $(n-k)\sigma^2$ . We could finish by uncovering the projector identities hidden inside $\mathbf M$ ; the trace proof in Appendix A takes exactly that route. But there is a more revealing option: instead of simplifying the projector in awkward coordinates, choose coordinates in which the projector is already simple.

Choose an orthonormal basis $\mathcal B=(\mathbf u_1,\ldots,\mathbf u_n)$ aligned with the two perpendicular spaces: the first $k$ basis vectors span $\operatorname{col}(\mathbf X)$ , and the remaining $n-k$ span its orthogonal complement. Because the basis is orthonormal, this changes only the coordinates used to describe the residual vector, not its squared length $\|\mathbf r\|^2$ . In these coordinates, the projection becomes simple bookkeeping.

Because $\mathbf H$ projects onto $\operatorname{col}(\mathbf X)$ , it leaves each of the first $k$ basis vectors unchanged. Each of the remaining basis vectors is perpendicular to that space, so $\mathbf H$ sends it to zero:

\mathbf H\mathbf u_j= \begin{cases} \mathbf u_j,&j\le k,\\ \mathbf0,&j>k. \end{cases}

The notation $[\mathbf H]_{\mathcal B}$ means “the matrix of the transformation $\mathbf H$ when vectors are described using coordinates in the basis $\mathcal B$ .” The action above means that, in these coordinates, $\mathbf H$ returns the first $k$ coordinates unchanged and replaces the remaining $n-k$ coordinates with zero. Its residual counterpart $\mathbf I-\mathbf H$ does the reverse:

[\mathbf H]_{\mathcal B}= \begin{pmatrix} \mathbf I_k&\mathbf0\\ \mathbf0&\mathbf0 \end{pmatrix}, \qquad [\mathbf I-\mathbf H]_{\mathcal B}= \begin{pmatrix} \mathbf0&\mathbf0\\ \mathbf0&\mathbf I_{n-k} \end{pmatrix}.

Thus $\mathbf H$ keeps the first $k$ coordinates and deletes the remaining $n-k$ , while $\mathbf I-\mathbf H$ does the reverse.

Now expand the noise in this basis:

\boldsymbol\epsilon =\sum_{j=1}^n z_j\mathbf u_j, \qquad z_j=\mathbf u_j^\top\boldsymbol\epsilon =\sum_{i=1}^n u_{ij}\epsilon_i,

where $u_{ij}$ is the $i$ th entry of $\mathbf u_j$ . Each $z_j$ is therefore a weighted sum of the original noise coordinates. White noise means that the $\epsilon_i$ are uncorrelated, each has mean zero, and each has variance $\sigma^2$ .

Key Step

$z_j$ is the amount of noise pointing along the unit direction $\mathbf u_j$ . The next calculation shows that every unit direction contains expected noise energy $\sigma^2$ , even though changing the direction changes the weighted combination $z_j=\sum_i u_{ij}\epsilon_i$ . This is what will make the rest of the proof a dimension count: every residual direction that survives contributes one $\sigma^2$ .

Now compute its expected squared magnitude:

\begin{aligned} \mathbb E[z_j^2\mid\mathbf X] &=\sum_{i=1}^n u_{ij}^2\,\mathbb E[\epsilon_i^2\mid\mathbf X] +\underbrace{\sum_{i\ne\ell}u_{ij}u_{\ell j}\, \mathbb E[\epsilon_i\epsilon_\ell\mid\mathbf X]}_{0}\\ &=\sigma^2\sum_{i=1}^n u_{ij}^2 =\sigma^2\|\mathbf u_j\|^2 =\sigma^2. \end{aligned}

Three facts work together in this calculation. First, the original noise coordinates are uncorrelated, so the cross-terms vanish. Second, each coordinate has variance $\sigma^2$ , leaving $\sigma^2\sum_i u_{ij}^2$ . Third, $\mathbf u_j$ is a unit vector, so $\sum_i u_{ij}^2=1$ . Since $z_j$ also has mean zero, its expected square equals its variance.

Finally, $\mathbf I-\mathbf H$ deletes the first $k$ coordinates and keeps the remaining $n-k$ :

\begin{aligned} \mathbf r &=\sum_{j=k+1}^n z_j\mathbf u_j,\\ \|\mathbf r\|^2 &=\sum_{j=k+1}^n z_j^2,\\ \mathbb E[\|\mathbf r\|^2\mid\mathbf X] &=\underbrace{\sigma^2+\cdots+\sigma^2}_{n-k\text{ residual directions}} =(n-k)\sigma^2. \end{aligned}

That is the result:

Variance is noise energy per direction. OLS uses $k$ directions for the fit and leaves $n-k$ directions in the residuals.

The divisor counts residual directions, not discarded observations. Therefore, for $n>k$ , the unbiased residual-variance estimator is

Residual Variance in Full-Rank OLS

s^2=\frac{1}{n-k}\sum_{i=1}^n r_i^2, \qquad \mathbb E[s^2\mid\mathbf X]=\sigma^2.

Why n − 1 Is the Same Story

The previous post is now just the case $k=1$ . A mean-only model has $\mathbf X=\mathbf1$ , so its fitted space is the single direction containing all constant vectors. Its hat matrix is the averaging operator:

\begin{aligned} \mathbf H &=\mathbf1(\mathbf1^\top\mathbf1)^{-1}\mathbf1^\top =\frac1n\mathbf1\mathbf1^\top,\\ \mathbf H\mathbf y&=(\hat m,\ldots,\hat m)^\top,\\ \mathbf1^\top\mathbf r&=\sum_i r_i=0. \end{aligned}

Fitting the mean uses one direction. The residual must lie in the perpendicular hyperplane whose coordinates sum to zero, leaving $n-1$ directions and therefore

\mathbb E\!\left[\sum_i(y_i-\hat m)^2\right]=(n-1)\sigma^2.

For a line with an intercept and slope, the residuals instead satisfy two orthogonality constraints,

\sum_i r_i=0, \qquad \sum_i x_i r_i=0,

leaving $n-2$ directions. The two cases are the same projection at different dimensions.

Case	Fitted space	Fitted dimension	Residual dimension	Divisor
Mean only	$\operatorname{span}(\mathbf1)$	$1$	$n-1$	$n-1$
Intercept + slope	$\operatorname{span}(\mathbf1,\mathbf x)$	$2$	$n-2$	$n-2$
Full-rank OLS	$\operatorname{col}(\mathbf X)$	$k$	$n-k$	$n-k$

What the Design Changes, and What It Does Not

The result makes a stronger claim than “divide by $n-k$ .” The predictor values determine which $k$ -dimensional space OLS fits, so different designs can produce very different fitted lines and residual patterns. But under white noise, the orientation of that space does not affect the expected total residual energy. Every full-rank $k$ -dimensional fit leaves $n-k$ noise directions.

The animation below checks that design-invariance numerically. Each trial draws fresh white noise for two different designs. Their individual residual sums of squares fluctuate, but both running averages converge to the same target, $(n-k)\sigma^2$ .

draws: 1

Design A: spread out

x's evenly spaced

this draw: Σr² = 7.4avg = 7.39

Design B: clustered + one far point

7 bunched x's, one far out

this draw: Σr² = 10.0avg = 10.03

Two designs, same noise (σ² = 1), n = 8, k = 2. Press Run: each draw adds fresh noise to the true line (faint dashed) and refits (solid). The residual patterns differ: in Design B the far-out predictor value strongly anchors the line, so its residual stays small while the cluster absorbs the slack. Yet the running average of Σr² converges to the same (n − k)σ² = 6 in both. The design decides where the residuals go; only k decides how much of the noise they keep.

This is why the divisor depends on the dimension $k$ , not on the particular values inside $\mathbf X$ .

The Takeaway

If you retain one picture, let it be this: OLS can move only within the $k$ -dimensional space $\operatorname{col}(\mathbf X)$ . It chooses the closest point there, leaving a perpendicular residual in the remaining $n-k$ directions. White noise supplies expected energy $\sigma^2$ to each direction, so

\mathbb E[\mathrm{RSS}\mid\mathbf X]=(n-k)\sigma^2, \qquad s^2=\frac{\mathrm{RSS}}{n-k}.

The $1$ in $n-1$ and the $k$ in $n-k$ count fitted directions, not discarded observations.

Appendix B gives a less intuitive but direct algebraic derivation by extending the exact optimization identity from the previous post to $k$ dimensions.

Appendix A: The Same Count Using the Trace

The geometric proof is complete. This appendix records the same direction count algebraically.

$\mathbf H$ is an orthogonal projector, not an orthogonal matrix. Its symmetry guarantees an orthonormal eigenbasis. Along the $k$ fitted directions, $\mathbf H$ has eigenvalue $1$ ; along the $n-k$ residual directions, it has eigenvalue $0$ . The residual projector $\mathbf I-\mathbf H$ reverses those eigenvalues, so its trace is $n-k$ .

The defining properties of the projector and cyclicity of the trace give the fitted count directly:

\begin{aligned} \mathbf H^\top&=\mathbf H, &\mathbf H^2&=\mathbf H,\\[0.3em] \operatorname{tr}(\mathbf H)&=\operatorname{rank}(\mathbf H)=k,\\[0.3em] \operatorname{tr}(\mathbf H) &=\operatorname{tr}\!\left( \mathbf X(\mathbf X^\top\mathbf X)^{-1}\mathbf X^\top \right)\\ &=\operatorname{tr}\!\left( (\mathbf X^\top\mathbf X)^{-1}\mathbf X^\top\mathbf X \right) =\operatorname{tr}(\mathbf I_k)=k. \end{aligned}

Now write the residual sum of squares as a quadratic form:

\begin{aligned} \mathrm{RSS} &=\mathbf r^\top\mathbf r\\ &=\boldsymbol\epsilon^\top(\mathbf I-\mathbf H)^\top (\mathbf I-\mathbf H)\boldsymbol\epsilon\\ &=\boldsymbol\epsilon^\top(\mathbf I-\mathbf H)\boldsymbol\epsilon,\\[0.5em] \mathbb E[\mathrm{RSS}\mid\mathbf X] &=\sigma^2\operatorname{tr}(\mathbf I-\mathbf H)\\ &=\sigma^2\bigl(n-\operatorname{tr}(\mathbf H)\bigr) =(n-k)\sigma^2. \end{aligned}

Appendix B: The Previous Proof in k Dimensions

The previous post used an exact scalar identity. The same identity exists for a $k$ -parameter least-squares fit because $L$ is quadratic and its gradient vanishes at $\hat\theta$ . Here

L(\vartheta)=\frac1n\|\mathbf y-\mathbf X\vartheta\|^2,

the normalized least-squares loss. This has the same minimizer as the unnormalized squared distance above, but the factor $\frac1n$ matters for the Hessian scale. We also use the standard conditional unbiasedness fact $\mathbb E[\hat\theta\mid\mathbf X]=\theta$ under the assumptions stated at the start of the post. For OLS, the Hessian of the loss and the covariance of the estimator contain inverse factors that cancel when the expected quadratic gap is written as a trace:

\begin{aligned} V(\mu)&=V(\hat m)+(\hat m-\mu)^2,\\[0.4em] L(\theta) &=L(\hat\theta) +\tfrac12(\theta-\hat\theta)^\top\nabla^2L\,(\theta-\hat\theta),\\[0.4em] \nabla^2L&=\frac2n\mathbf X^\top\mathbf X,\\ \operatorname{Cov}(\hat\theta\mid\mathbf X) &=\sigma^2(\mathbf X^\top\mathbf X)^{-1},\\[0.4em] \mathbb E[L(\theta)-L(\hat\theta)\mid\mathbf X] &=\tfrac12\operatorname{tr}\!\left( \nabla^2L\,\operatorname{Cov}(\hat\theta\mid\mathbf X) \right)\\ &=\frac{\sigma^2}{n}\operatorname{tr}(\mathbf I_k) =\frac{k\sigma^2}{n}. \end{aligned}

This is the previous proof lifted from one dimension to $k$ . It computes the expected reduction in normalized squared error obtained by fitting the $k$ coefficients. Since

L(\theta)=\frac1n\|\boldsymbol\epsilon\|^2, \qquad L(\hat\theta)=\frac{\mathrm{RSS}}{n}, \qquad \mathbb E[L(\theta)\mid\mathbf X]=\sigma^2,

the result above immediately gives

\mathbb E[L(\hat\theta)\mid\mathbf X] =\sigma^2-\frac{k\sigma^2}{n} =\frac{n-k}{n}\sigma^2,

or, after multiplying by $n$ , $\mathbb E[\mathrm{RSS}\mid\mathbf X]=(n-k)\sigma^2$ , exactly the result in Appendix A.

The two appendices therefore count complementary pieces of the same noise. Appendix B counts the expected $k\sigma^2$ absorbed by the fitted projector $\mathbf H$ ; Appendix A counts the expected $(n-k)\sigma^2$ left by the residual projector $\mathbf I-\mathbf H$ . Their traces add to $n$ :

\operatorname{tr}(\mathbf H)+\operatorname{tr}(\mathbf I-\mathbf H) =k+(n-k)=n.

References

Hoaglin, D. C., & Welsch, R. E. (1978). The hat matrix in regression and ANOVA. The American Statistician, 32(1), 17–22.
Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning (2nd ed.), §7.5–7.6. Springer.
Efron, B., & Hastie, T. (2016). Computer Age Statistical Inference. Cambridge University Press.

The Puzzle

The Model Is a Subspace

OLS Creates the Right Angle

The Hat Matrix Performs the Projection

Why Only Noise Remains in the Residual

Counting the Noise Directions

Why n − 1 Is the Same Story

What the Design Changes, and What It Does Not

The Takeaway

Appendix A: The Same Count Using the Trace

Appendix B: The Previous Proof in k Dimensions

References

Discussion