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 ∈ R n × p X\in\mathbb R^{n\times p} X ∈ R n × p and suppose
Y = X β + ε , E [ ε ∣ X ] = 0 , Cov ( ε ∣ X ) = σ 2 I n . Y=X\beta+\varepsilon,\qquad E[\varepsilon\mid X]=0,\qquad \operatorname{Cov}(\varepsilon\mid X)=\sigma^2I_n. Y = Xβ + ε , E [ ε ∣ X ] = 0 , Cov ( ε ∣ X ) = σ 2 I n . The rows are x i ⊤ x_i^\top x i ⊤ , so E [ Y i ∣ X ] = x i ⊤ β E[Y_i\mid X]=x_i^\top\beta E [ Y i ∣ X ] = x i ⊤ β . 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. S ( β ) = ∥ y − Xβ ∥ 2 = y ⊤ y − 2 β ⊤ X ⊤ y + β ⊤ X ⊤ Xβ . Differentiate:
∇ S ( β ) = − 2 X ⊤ ( y − X β ) , ∇ 2 S ( β ) = 2 X ⊤ X . \nabla S(\beta)=-2X^\top(y-X\beta),\qquad \nabla^2S(\beta)=2X^\top X. ∇ S ( β ) = − 2 X ⊤ ( y − Xβ ) , ∇ 2 S ( β ) = 2 X ⊤ X . The stationary condition is
X ⊤ X β ^ = X ⊤ y . X^\top X\widehat\beta=X^\top y. X ⊤ X β = X ⊤ y . If rank ( X ) = p \operatorname{rank}(X)=p rank ( X ) = p , the Hessian is positive definite and the unique minimizer is
β ^ = ( X ⊤ X ) − 1 X ⊤ y . \boxed{\widehat\beta=(X^\top X)^{-1}X^\top y.} β = ( X ⊤ X ) − 1 X ⊤ y . Computation: use QR, SVD, or a least-squares solver rather than explicitly calculating the inverse. Forming X ⊤ X X^\top X X ⊤ X squares the condition number and can amplify numerical difficulties.
Under Y ∣ X ∼ N ( X β , σ 2 I ) Y\mid X\sim N(X\beta,\sigma^2I) Y ∣ X ∼ N ( Xβ , σ 2 I ) ,
ℓ ( β , σ 2 ) = − n 2 log ( 2 π σ 2 ) − ∥ y − X β ∥ 2 2 σ 2 . \ell(\beta,\sigma^2)=-\frac n2\log(2\pi\sigma^2)
-\frac{\lVert y-X\beta\rVert^2}{2\sigma^2}. ℓ ( β , σ 2 ) = − 2 n log ( 2 π σ 2 ) − 2 σ 2 ∥ y − Xβ ∥ 2 . For fixed positive variance, maximizing over β \beta β minimizes squared error. When variance is also estimated and RSS is positive, the Gaussian MLE is σ ^ M L E 2 = RSS / n \widehat\sigma^2_{\mathrm{MLE}}=\operatorname{RSS}/n σ MLE 2 = RSS / n .
The mean model also implies the moment condition
E [ X ⊤ ( Y − X β ) ∣ X ] = 0. E[X^\top(Y-X\beta)\mid X]=0. E [ X ⊤ ( Y − Xβ ) ∣ X ] = 0. Setting its observed counterpart to zero gives the same normal equations. This moment argument only needs the mean restriction; Gaussianity is unnecessary.
Geometry: Orthogonal Projection ¶ The fitted vector and residual are
y ^ = H y , H = X ( X ⊤ X ) − 1 X ⊤ , e = ( I − H ) y . \widehat y=Hy,\qquad H=X(X^\top X)^{-1}X^\top,\qquad e=(I-H)y. y = Hy , H = X ( X ⊤ X ) − 1 X ⊤ , e = ( I − H ) y . Property Meaning y ^ ∈ col ( X ) \widehat y\in\operatorname{col}(X) y ∈ col ( X ) The fit is a linear combination of the design columns. X ⊤ e = 0 X^\top e=0 X ⊤ e = 0 The residual is orthogonal to every design column. H ⊤ = H H^\top=H H ⊤ = H The hat matrix is symmetric. H 2 = H H^2=H H 2 = H Projecting an already fitted vector does nothing. H e = 0 He=0 He = 0 The projection removes the residual component.
For any candidate X b Xb X b ,
y − X b = e + ( y ^ − X b ) . y-Xb=e+(\widehat y-Xb). y − X b = e + ( y − X b ) . The two terms are orthogonal because e ⊥ col ( X ) e\perp\operatorname{col}(X) e ⊥ col ( X ) . Hence
∥ y − X b ∥ 2 = ∥ e ∥ 2 + ∥ y ^ − X b ∥ 2 ≥ ∥ e ∥ 2 . \lVert y-Xb\rVert^2=\lVert e\rVert^2+\lVert\widehat y-Xb\rVert^2\geq\lVert e\rVert^2. ∥ y − X b ∥ 2 = ∥ e ∥ 2 + ∥ y − X b ∥ 2 ≥ ∥ e ∥ 2 . Thus y ^ \widehat y 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 1 ⊤ e = 0 , so residuals sum to zero and the fitted responses have mean y ˉ \bar y 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 X X X is rank deficient, there are infinitely many coefficient solutions differing by vectors in the null space of X X X . They all produce the same fitted vector:
y ^ = X X + y , \widehat y=XX^+y, y = X X + y , where X + X^+ X + is the Moore-Penrose pseudoinverse. The minimum-norm coefficient solution is β ^ = X + y \widehat\beta=X^+y β = X + y . The fitted subspace has dimension r = rank ( X ) r=\operatorname{rank}(X) r = rank ( X ) , so tr ( X X + ) = r \operatorname{tr}(XX^+)=r tr ( X X + ) = r .
Sampling Distribution ¶ Write A = ( X ⊤ X ) − 1 X ⊤ A=(X^\top X)^{-1}X^\top A = ( X ⊤ X ) − 1 X ⊤ . Since β ^ = A Y \widehat\beta=AY β = A Y and A X = I AX=I A X = I ,
E [ β ^ ∣ X ] = β , Cov ( β ^ ∣ X ) = A ( σ 2 I ) 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}. E [ β ∣ X ] = β , Cov ( β ∣ X ) = A ( σ 2 I ) A ⊤ = σ 2 ( X ⊤ 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).} β ∣ X ∼ N ( β , σ 2 ( X ⊤ X ) − 1 ) . For any fixed vector a a a ,
E [ a ⊤ β ^ ∣ X ] = a ⊤ β , Var ( a ⊤ β ^ ∣ X ) = σ 2 a ⊤ ( X ⊤ X ) − 1 a . 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. E [ a ⊤ β ∣ X ] = a ⊤ β , Var ( a ⊤ β ∣ X ) = σ 2 a ⊤ ( X ⊤ X ) − 1 a . At a new design row x 0 ⊤ x_0^\top x 0 ⊤ , this gives the variance of the estimated mean x 0 ⊤ β ^ x_0^\top\widehat\beta x 0 ⊤ β . Predicting an independent future response adds its noise variance σ 2 \sigma^2 σ 2 under the same model.
For the fitted vector,
Cov ( Y ^ ∣ X ) = σ 2 H , \operatorname{Cov}(\widehat Y\mid X)=\sigma^2H, Cov ( Y ∣ X ) = σ 2 H , and, when n > p n>p n > p , σ ^ 2 = RSS / ( n − p ) \widehat\sigma^2=\operatorname{RSS}/(n-p) σ 2 = 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 β ~ = C Y \widetilde\beta=CY β = C Y that are linear in Y Y Y 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. Cov ( β ∣ X ) − Cov ( β ∣ X ) ⪰ 0. This means every fixed linear combination a ⊤ β ^ a^\top\widehat\beta a ⊤ β has variance no larger than its counterpart a ⊤ β ~ a^\top\widetilde\beta a ⊤ β .
Write C = A + D C=A+D C = A + D , where A = ( X ⊤ X ) − 1 X ⊤ A=(X^\top X)^{-1}X^\top A = ( X ⊤ X ) − 1 X ⊤ . Unbiasedness requires C X = I CX=I CX = I , so D X = 0 DX=0 D X = 0 .
Then A D ⊤ = ( X ⊤ X ) − 1 ( D X ) ⊤ = 0 AD^\top=(X^\top X)^{-1}(DX)^\top=0 A D ⊤ = ( X ⊤ X ) − 1 ( D X ) ⊤ = 0 , and similarly D A ⊤ = 0 DA^\top=0 D A ⊤ = 0 . Therefore
Cov ( C Y ∣ X ) = σ 2 ( A + D ) ( A + D ) ⊤ = σ 2 ( X ⊤ X ) − 1 + σ 2 D D ⊤ . \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} Cov ( C Y ∣ X ) = σ 2 ( A + D ) ( A + D ) ⊤ = σ 2 ( X ⊤ X ) − 1 + σ 2 D D ⊤ . For any a a a ,
a ⊤ σ 2 D D ⊤ a = σ 2 ∥ D ⊤ a ∥ 2 ≥ 0. a^\top\sigma^2DD^\top a=\sigma^2\lVert D^\top a\rVert^2\geq0. a ⊤ σ 2 D D ⊤ a = σ 2 ∥ D ⊤ a ∥ 2 ≥ 0. 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) x = ( − 1 , 0 , 1 ) and y = ( 1 , 0 , 1 ) y=(1,0,1) y = ( 1 , 0 , 1 ) . Then
X ⊤ X = ( 3 0 0 2 ) , X ⊤ y = ( 2 0 ) , β ^ = ( 2 / 3 0 ) . 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}. X ⊤ X = ( 3 0 0 2 ) , X ⊤ y = ( 2 0 ) , β = ( 2/3 0 ) . The fitted vector is ( 2 / 3 , 2 / 3 , 2 / 3 ) ⊤ (2/3,2/3,2/3)^\top ( 2/3 , 2/3 , 2/3 ) ⊤ and e = ( 1 / 3 , − 2 / 3 , 1 / 3 ) ⊤ e=(1/3,-2/3,1/3)^\top e = ( 1/3 , − 2/3 , 1/3 ) ⊤ . Both residual sums ∑ i e i \sum_i e_i ∑ i e i and ∑ i x i e i \sum_i x_i e_i ∑ i x i e i equal zero. RSS is 2 / 3 2/3 2/3 , so the residual variance estimate is 2 / 3 2/3 2/3 with one residual degree of freedom.