Legendre与Gauss在19世纪初提出了最小二乘的思想，1900年Markov证明了最小二乘估计的性质良好，在此之后最小二乘就开始广泛应用于线性模型的估计了。对于线性模型 y = X β + ϵ , E ϵ = 0 , C o v ( ϵ ) = σ 2 I y=X\beta + \epsilon,E\epsilon=0,Cov(\epsilon)=\sigma^2I y=Xβ+ϵ,Eϵ=0,Cov(ϵ)=σ2I
其中 y , ϵ y,\epsilon y,ϵ为 n × 1 n\times 1 n×1的向量， X X X为 n × p n \times p n×p的Design Matrix，如果 r a n k ( X ) ≥ p rank(X)\ge p rank(X)≥p，称这个线性模型为满秩的；否则称之为降秩的。这部分我们将介绍普通最小二乘法（OLS）、带约束的最小二乘法、广义最小二乘法（GLS）、稳健性、两步法、最小二乘法的几何解释以及常用数值算法，这一篇介绍OLS。
OLS的思路是 min β Q = ∥ e ∥ 2 = ( y − X β ) ′ ( y − X β ) = y ′ y − 2 y ′ X β + β ′ X ′ X β \min_{\beta}\ \ Q = \left\| e \right\|^2 = (y-X\beta)'(y-X\beta)=y'y-2y'X\beta+\beta'X'X\beta βmin Q=∥e∥2=(y−Xβ)′(y−Xβ)=y′y−2y′Xβ+β′X′Xβ
计算 Q Q Q关于 β \beta β的梯度 ∇ β Q = − 2 X ′ y + 2 X ′ X β = 0 ⇒ X ′ X β = X ′ y \nabla_{\beta} Q=-2X'y+2X'X\beta=0 \Rightarrow X'X\beta = X'y ∇βQ=−2X′y+2X′Xβ=0⇒X′Xβ=X′y
这个方程叫做OLS的正则方程，求解这个方程可以得到系数的OLS估计，并且基于这个方程还可以获得残差的性质。 X ′ y X'y X′y在 X ′ X' X′的列空间中，因此这个方程是相容的，可以用系数矩阵的广义逆表示解： β ^ = ( X ′ X ) − X ′ y \hat{\beta} = (X'X)^{-}X'y β^=(X′X)−X′y
计算 Q Q Q关于 β \beta β的Hessian矩阵， H β Q = 2 X ′ X ≥ 0 H_{\beta}Q = 2X'X\ge 0 HβQ=2X′X≥0
因此 β ^ \hat{\beta} β^使 Q Q Q取最小值，并且最小值点唯一。
下面考虑广义逆的确定。假设 r a n k ( X ) ≥ p rank(X)\ge p rank(X)≥p，则 X ′ X X'X X′X是满秩的方阵， β ^ = ( X ′ X ) − 1 X ′ y \hat{\beta} = (X'X)^{-1}X'y β^=(X′X)−1X′y
假设 r a n k ( X ) < p rank(X)<p rank(X)<p，则 X ′ X X'X X′X降秩，它的逆不存在，此时不存在 β \beta β的线性无偏估计。 证明 假设 A y Ay Ay是线性无偏估计，则 E ( A y ) = A X β = β ⇒ A X = I p ⇒ r a n k ( A X ) = p E(Ay) = AX\beta = \beta \Rightarrow AX = I_p \Rightarrow rank(AX)=p E(Ay)=AXβ=β⇒AX=Ip⇒rank(AX)=p，然而 r a n k ( A X ) ≤ r a n k ( X ) < p rank(AX)\le rank(X)<p rank(AX)≤rank(X)<p，这就出现了矛盾。
如果我们就停留在上一节的讨论，那么OLS的局限性是很强的，因为我们有时并不一定是要关注 β \beta β，也不一定总是需要线性无偏估计。
在 r a n k ( X ) < p rank(X)<p rank(X)<p这种情况中，称 β \beta β是不可估计，但我们可以考察 β \beta β的某种线性组合（参考RCD的线性组合与contract理论） c ′ β c'\beta c′β，如果 ∃ a \exists a ∃a， n × 1 n \times 1 n×1的列向量，使得 E [ a ′ y ] = c ′ β E[a'y]=c'\beta E[a′y]=c′β，就称 c ′ β c'\beta c′β是可估函数。显然 c c c属于 X ′ X' X′的列空间。如果 c 1 ′ β c_1'\beta c1′β与 c 2 ′ β c_2'\beta c2′β中的 c 1 , c 2 c_1,c_2 c1,c2线性无关，则称 c 1 ′ β c_1'\beta c1′β与 c 2 ′ β c_2'\beta c2′β线性无关。因为 c c c属于 X ′ X' X′的列空间，所以一组线性无关的可估函数最多有 r a n k ( X ) rank(X) rank(X)个，并且据此可以得出 ∃ a \exists a ∃a, c = X ′ a c=X'a c=X′a，进而 c ′ β ^ = c ′ ( X ′ X ) − X ′ y = a ′ X ( X ′ X ) − X ′ y c'\hat{\beta} = c'(X'X)^{-}X'y=a'X(X'X)^{-}X'y c′β^=c′(X′X)−X′y=a′X(X′X)−X′y
参考矩阵分析与多元统计那个系列，包含广义逆的 X ( X ′ X ) − X ′ X(X'X)^{-}X' X(X′X)−X′一项与广义逆的选取无关，这个性质保证 c ′ β ^ c'\hat{\beta} c′β^是一个良定义。另外， E ( c ′ β ^ ) = a ′ X ( X ′ X ) − X ′ E y = a ′ X ( X ′ X ) − X ′ X β = a ′ X β = c ′ β E(c'\hat{\beta}) = a'X(X'X)^{-}X'Ey=a'X(X'X)^{-}X'X\beta=a'X\beta=c'\beta E(c′β^)=a′X(X′X)−X′Ey=a′X(X′X)−X′Xβ=a′Xβ=c′β
说明 c ′ β ^ c'\hat{\beta} c′β^是可估函数 c β c\beta cβ的无偏估计。基于这些分析，我们可以自信地定义 c ′ β ^ c'\hat{\beta} c′β^为 c ′ β c'\beta c′β的OLS估计。OLS是具有唯一性的，但在回归那个系列我们讨论过，线性无偏估计不具有唯一性，但Gauss-Markov定理指出OLS估计是最优线性无偏估计（Best Linear Unbiased Estimator，BLUE）：
Gauss-Markov定理 c ′ β ^ c'\hat{\beta} c′β^是 c ′ β c'\beta c′β的BLUE。 证明 无偏性说明过了，下面讨论最优性（所有线性无偏估计中方差最小）。计算 c ′ β ^ c'\hat{\beta} c′β^的方差： V a r ( c ′ β ^ ) = V a r ( a ′ X ( X ′ X ) − X ′ y ) = σ 2 a ′ X ( X ′ X ) − X ′ X ( X ′ X ) − X ′ a = σ 2 a ′ X ( X ′ X ) − X ′ a = σ 2 c ′ ( X ′ X ) − c Var(c'\hat{\beta}) = Var(a'X(X'X)^{-}X'y) =\sigma^2a'X(X'X)^{-}X'X(X'X)^{-}X'a \\= \sigma^2a'X(X'X)^{-}X'a = \sigma^2 c'(X'X)^{-}c Var(c′β^)=Var(a′X(X′X)−X′y)=σ2a′X(X′X)−X′X(X′X)−X′a=σ2a′X(X′X)−X′a=σ2c′(X′X)−c
假设 b ′ y b'y b′y是 c ′ β c'\beta c′β的另一个线性无偏估计，则 E ( b ′ y ) = b ′ X β = c ′ β E(b'y)=b'X\beta=c'\beta E(b′y)=b′Xβ=c′β，也就是 c = X ′ b c=X'b c=X′b，计算 V a r ( b ′ y ) = σ 2 b ′ b Var(b'y) = \sigma^2b'b Var(b′y)=σ2b′b
考虑两个方差的差， V a r ( b ′ y ) − V a r ( c ′ β ^ ) = σ 2 [ b ′ b − c ′ ( X ′ X ) − c ] = σ 2 [ b ′ − c ′ ( X ′ X ) − ] [ b − ( X ′ X ) − c ] = σ 2 ∥ b − ( X ′ X ) − c ∥ ≥ 0 Var(b'y)-Var(c'\hat{\beta}) =\sigma^2[b'b-c'(X'X)^{-}c] \\ = \sigma^2[b'-c'(X'X)^{-}][b-(X'X)^{-}c]=\sigma^2 \left\| b-(X'X)^{-}c\right\| \ge 0 Var(b′y)−Var(c′β^)=σ2[b′b−c′(X′X)−c]=σ2[b′−c′(X′X)−][b−(X′X)−c]=σ2∥∥b−(X′X)−c∥∥≥0
所以OLS是BLUE。
模型的残差为 e = y − X β ^ = ( I − P X ) y e=y-X\hat{\beta}=(I-P_X)y e=y−Xβ^=(I−PX)y，其中 P X = X ( X ′ X ) − X ′ P_X=X(X'X)^{-}X' PX=X(X′X)−X′是到 X X X列空间中的投影矩阵，可以计算 E e ^ = ( I − P X ) X β = X β − X ( X ′ X ) − X ′ X β = 0 C o v ( e ^ ) = σ 2 ( I − P X ) ( I − P X ) ′ = σ 2 ( I − P X ) E\hat{e}=(I-P_X)X\beta=X\beta - X(X'X)^{-}X'X\beta=0 \\ Cov(\hat{e})=\sigma^2(I-P_X)(I-P_X)'=\sigma^2(I-P_X) Ee^=(I−PX)Xβ=Xβ−X(X′X)−X′Xβ=0Cov(e^)=σ2(I−PX)(I−PX)′=σ2(I−PX)
基于这个结果可以构造 σ 2 \sigma^2 σ2的无偏估计为 σ ^ 2 = e ^ ′ e ^ n − r a n k ( X ) \hat\sigma^2=\frac{\hat{e}'\hat{e}}{n-rank(X)} σ^2=n−rank(X)e^′e^
证明 首先 e ^ ′ e ^ = y ′ ( I − P X ) y \hat{e}'\hat{e}=y'(I-P_X)y e^′e^=y′(I−PX)y，下面计算 E [ e ^ ′ e ^ ] = ( X β ) ′ ( I − P X ) + t r [ ( I − P X ) C o v ( y ) ] = σ 2 t r ( I − P X ) = σ 2 ( n − r a n k ( X ) ) E[\hat{e}'\hat{e}]=(X\beta)'(I-P_X)+tr[(I-P_X)Cov(y)] \\=\sigma^2tr(I-P_X)=\sigma^2(n-rank(X)) E[e^′e^]=(Xβ)′(I−PX)+tr[(I−PX)Cov(y)]=σ2tr(I−PX)=σ2(n−rank(X))
第一行到第二行应用的第一个结论是 ( I − P X ) X = 0 (I-P_X)X=0 (I−PX)X=0，用到的第二个结论是如果 E X = μ , C o v ( X ) = Σ EX=\mu,Cov(X)=\Sigma EX=μ,Cov(X)=Σ，则 E [ X ′ A X ] = μ ′ A μ + t r ( A Σ ) E[X'AX]=\mu'A\mu+tr(A\Sigma) E[X′AX]=μ′Aμ+tr(AΣ)
下面证明这个恒等式。计算 X ′ A X = ( X − μ + μ ) ′ A ( X − μ + μ ) = ( X − μ ) ′ A ( X − μ ) + μ ′ A ( X − μ ) + ( X − μ ) ′ A μ + μ ′ A μ X'AX=(X-\mu+\mu)'A(X-\mu+\mu) \\ = (X-\mu)'A(X-\mu)+\mu'A(X-\mu)+(X-\mu)'A\mu+\mu'A\mu X′AX=(X−μ+μ)′A(X−μ+μ)=(X−μ)′A(X−μ)+μ′A(X−μ)+(X−μ)′Aμ+μ′Aμ
接下来求期望， E [ μ ′ A ( X − μ ) ] = E [ μ ′ A X ] − μ ′ A μ = 0 E[\mu'A(X-\mu)]=E[\mu'AX]-\mu'A\mu=0 E[μ′A(X−μ)]=E[μ′AX]−μ′Aμ=0
类似的，第三项的期望也为0，计算第一项的期望， E [ ( X − μ ) ′ A ( X − μ ) ] = E t r [ ( X − μ ) ′ A ( X − μ ) ] = E t r [ A ( X − μ ) ( X − μ ) ′ ] = t r A E [ ( X − μ ) ( X − μ ) ′ ] = t r ( A Σ ) E[(X-\mu)'A(X-\mu)]=Etr[(X-\mu)'A(X-\mu)] \\ = Etr[A(X-\mu)(X-\mu)'] = trAE[(X-\mu)(X-\mu)'] = tr(A\Sigma) E[(X−μ)′A(X−μ)]=Etr[(X−μ)′A(X−μ)]=Etr[A(X−μ)(X−μ)′]=trAE[(X−μ)(X−μ)′]=tr(AΣ)
这就完成了整个证明。在回归那个系列给出过另一种证明，但思路类似，都是把数量用trace表示，在利用trace中矩阵乘法满足交换律的技巧。
需要注意的是OLS对随机误差的分布形式是没有要求的，只有当我们试图对OLS估计量做统计推断的时候，我们才需要考虑随机误差的分布形式。这是OLS与MLE一个很大的区别，MLE为了获得估计量就需要在一开始引入某种特定的分布，再去最大化给定数据在这种分布下的似然。
现在假设随机误差服从正态分布，则OLS有一些优越的性质，这些性质在MATH 564、566、571A的博客中都有过证明了，所以这里就简单归纳一下：
OLS也是MLE，且 c ′ β ^ ∼ N ( c ′ β , σ 2 c ′ ( X ′ X ) − c ) c'\hat{\beta} \sim N(c'\beta,\sigma^2c'(X'X)^{-}c) c′β^∼N(c′β,σ2c′(X′X)−c) σ 2 \sigma^2 σ2的MLE是 e ^ ′ e ^ / n \hat{e}'\hat{e}/n e^′e^/n，且 ( n − r a n k ( X ) ) σ ^ 2 ∼ σ 2 χ n − r a n k ( X ) 2 (n-rank(X))\hat{\sigma}^2\sim \sigma^2 \chi^2_{n-rank(X)} (n−rank(X))σ^2∼σ2χn−rank(X)2 c ′ β ^ c'\hat{\beta} c′β^与 σ ^ 2 \hat{\sigma}^2 σ^2互相独立 c ′ β ^ c'\hat{\beta} c′β^是唯一的UMVUE