Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Least Squares & Projection

Ordinary least squares chooses the linear combination of design columns closest to the observed response vector. Its geometry explains the normal equations, and its sampling properties explain both its strengths and the limits of its optimality.

The Linear Mean Model

Condition on a fixed design X∈Rn×pX\in\mathbb R^{n\times p} and suppose

Y=Xβ+ε,E[ε∣X]=0,Cov⁡(ε∣X)=σ2In.Y=X\beta+\varepsilon,\qquad E[\varepsilon\mid X]=0,\qquad \operatorname{Cov}(\varepsilon\mid X)=\sigma^2I_n.

The rows are xi⊤x_i^\top, so E[Yi∣X]=xi⊤βE[Y_i\mid X]=x_i^\top\beta. An intercept is represented by a column of ones. Assume full column rank for the inverse formulas below.

These mean and covariance assumptions do not require Gaussian errors. Normality will be added only when stating an exact normal sampling distribution.

Deriving the Normal Equations

Minimize

S(β)=∥y−Xβ∥2=y⊤y−2β⊤X⊤y+β⊤X⊤Xβ.S(\beta)=\lVert y-X\beta\rVert^2 =y^\top y-2\beta^\top X^\top y+\beta^\top X^\top X\beta.

Differentiate:

∇S(β)=−2X⊤(y−Xβ),∇2S(β)=2X⊤X.\nabla S(\beta)=-2X^\top(y-X\beta),\qquad \nabla^2S(\beta)=2X^\top X.

The stationary condition is

X⊤Xβ^=X⊤y.X^\top X\widehat\beta=X^\top y.

If rank⁡(X)=p\operatorname{rank}(X)=p, the Hessian is positive definite and the unique minimizer is

β^=(X⊤X)−1X⊤y.\boxed{\widehat\beta=(X^\top X)^{-1}X^\top y.}

Computation: use QR, SVD, or a least-squares solver rather than explicitly calculating the inverse. Forming X⊤XX^\top X squares the condition number and can amplify numerical difficulties.

Geometry: Orthogonal Projection

The fitted vector and residual are

y^=Hy,H=X(X⊤X)−1X⊤,e=(I−H)y.\widehat y=Hy,\qquad H=X(X^\top X)^{-1}X^\top,\qquad e=(I-H)y.
PropertyMeaning
y^∈col⁡(X)\widehat y\in\operatorname{col}(X)The fit is a linear combination of the design columns.
X⊤e=0X^\top e=0The residual is orthogonal to every design column.
H⊤=HH^\top=HThe hat matrix is symmetric.
H2=HH^2=HProjecting an already fitted vector does nothing.
He=0He=0The projection removes the residual component.
Why This Is the Closest Fitted Vector

For any candidate XbXb,

y−Xb=e+(y^−Xb).y-Xb=e+(\widehat y-Xb).

The two terms are orthogonal because e⊥col⁡(X)e\perp\operatorname{col}(X). Hence

∥y−Xb∥2=∥e∥2+∥y^−Xb∥2≥∥e∥2.\lVert y-Xb\rVert^2=\lVert e\rVert^2+\lVert\widehat y-Xb\rVert^2\geq\lVert e\rVert^2.

Thus y^\widehat y is the unique closest vector in the column space. Coefficients are unique only when the design has full column rank.

If an intercept is included, 1⊤e=0\mathbf1^\top e=0, so residuals sum to zero and the fitted responses have mean yˉ\bar y. These properties need not hold for a fit without an intercept or for a penalized intercept.

Rank-Deficient Designs

Least-squares solutions always exist. If XX is rank deficient, there are infinitely many coefficient solutions differing by vectors in the null space of XX. They all produce the same fitted vector:

y^=XX+y,\widehat y=XX^+y,

where X+X^+ is the Moore-Penrose pseudoinverse. The minimum-norm coefficient solution is β^=X+y\widehat\beta=X^+y. The fitted subspace has dimension r=rank⁡(X)r=\operatorname{rank}(X), so tr⁡(XX+)=r\operatorname{tr}(XX^+)=r.

Sampling Distribution

Write A=(X⊤X)−1X⊤A=(X^\top X)^{-1}X^\top. Since β^=AY\widehat\beta=AY and AX=IAX=I,

E[β^∣X]=β,Cov⁡(β^∣X)=A(σ2I)A⊤=σ2(X⊤X)−1.E[\widehat\beta\mid X]=\beta,\qquad \operatorname{Cov}(\widehat\beta\mid X)=A(\sigma^2I)A^\top =\sigma^2(X^\top X)^{-1}.

With Gaussian errors,

β^∣X∼N ⁣(β,σ2(X⊤X)−1).\boxed{\widehat\beta\mid X\sim N\!\left(\beta,\sigma^2(X^\top X)^{-1}\right).}

For any fixed vector aa,

E[a⊤β^∣X]=a⊤β,Var⁡(a⊤β^∣X)=σ2a⊤(X⊤X)−1a.E[a^\top\widehat\beta\mid X]=a^\top\beta,\qquad \operatorname{Var}(a^\top\widehat\beta\mid X)=\sigma^2a^\top(X^\top X)^{-1}a.

At a new design row x0⊤x_0^\top, this gives the variance of the estimated mean x0⊤β^x_0^\top\widehat\beta. Predicting an independent future response adds its noise variance σ2\sigma^2 under the same model.

For the fitted vector,

Cov⁡(Y^∣X)=σ2H,\operatorname{Cov}(\widehat Y\mid X)=\sigma^2H,

and, when n>pn>p, σ^2=RSS⁡/(n−p)\widehat\sigma^2=\operatorname{RSS}/(n-p) is unbiased under the stated mean and covariance assumptions.

Gauss-Markov: Best Linear Unbiased Estimation

What the theorem guarantees

Smallest covariance among estimators that are linear in the response and unbiased for every coefficient vector, under the stated assumptions.

What it leaves open

A biased estimator can have lower mean squared error. Ridge is outside the class compared by Gauss-Markov.

Among estimators β~=CY\widetilde\beta=CY that are linear in YY and unbiased for every β\beta, OLS has the smallest covariance matrix in the positive-semidefinite ordering:

Cov⁡(β~∣X)−Cov⁡(β^∣X)⪰0.\operatorname{Cov}(\widetilde\beta\mid X)-\operatorname{Cov}(\widehat\beta\mid X)\succeq0.

This means every fixed linear combination a⊤β^a^\top\widehat\beta has variance no larger than its counterpart a⊤β~a^\top\widetilde\beta.

Proof of Gauss-Markov

Write C=A+DC=A+D, where A=(X⊤X)−1X⊤A=(X^\top X)^{-1}X^\top. Unbiasedness requires CX=ICX=I, so DX=0DX=0.

Then AD⊤=(X⊤X)−1(DX)⊤=0AD^\top=(X^\top X)^{-1}(DX)^\top=0, and similarly DA⊤=0DA^\top=0. Therefore

Cov⁡(CY∣X)=σ2(A+D)(A+D)⊤=σ2(X⊤X)−1+σ2DD⊤.\begin{aligned} \operatorname{Cov}(CY\mid X) &=\sigma^2(A+D)(A+D)^\top\\ &=\sigma^2(X^\top X)^{-1}+\sigma^2DD^\top. \end{aligned}

For any aa,

a⊤σ2DD⊤a=σ2∥D⊤a∥2≥0.a^\top\sigma^2DD^\top a=\sigma^2\lVert D^\top a\rVert^2\geq0.

The difference is positive semidefinite. Gaussianity is not used.

Scope: “best” means smallest variance among linear unbiased estimators under the given covariance assumptions. It does not compare OLS with biased ridge estimators by mean squared error. It also does not establish the same result under heteroscedastic or correlated errors.

Worked Example

Fit an intercept and slope to x=(−1,0,1)x=(-1,0,1) and y=(1,0,1)y=(1,0,1). Then

X⊤X=(3002),X⊤y=(20),β^=(2/30).X^\top X=\begin{pmatrix}3&0\\0&2\end{pmatrix},\qquad X^\top y=\begin{pmatrix}2\\0\end{pmatrix},\qquad \widehat\beta=\begin{pmatrix}2/3\\0\end{pmatrix}.

The fitted vector is (2/3,2/3,2/3)⊤(2/3,2/3,2/3)^\top and e=(1/3,−2/3,1/3)⊤e=(1/3,-2/3,1/3)^\top. Both residual sums ∑iei\sum_i e_i and ∑ixiei\sum_i x_i e_i equal zero. RSS is 2/32/3, so the residual variance estimate is 2/32/3 with one residual degree of freedom.