最小二乘估计矩阵形式的推导

最近写文章有用到一些算法，自己推一下，顺便mark下来。

这么久我才发现csdn居然都能写Tex了(666)。

考虑一般线性回归模型(OLR)

考虑只含有一个指标的一般线性回归模型(ordinary linear regression model)有如下形式：

yi=β0+β1xi1+ϵ，i=1,2,…,ny_i=\beta_0+\beta_1x_{i1}+\epsilon，i=1,2,\dots,n yi​=β0​+β1​xi1​+ϵ，i=1,2,…,n

显然这是基于nnn个观测数据或者叫样本的模型形式。其中β0\beta_0β0​称为截距项系数，β1\beta_1β1​称为x1x_1x1​的回归系数，它们都是未知的常值参数。ϵ\epsilonϵ是不能被观测到的随机误差项，并且满足E(ϵ)=0E(\epsilon)=0E(ϵ)=0,Var(ϵ)=σ2>0\mathrm{Var(\epsilon)}=\sigma^2>0Var(ϵ)=σ2>0。其实是有x0x_0x0​的,只是通常认为x0=1x_0=1x0​=1。还有一个关键的假设就是xxx不是随机变量(xxx要都随机了,这模型就没法玩了)。

实际上我们所研究的问题往往包含多个指标。那么这些指标(x1,x2,...,xp)(x_1,x_2,...,x_p)(x1​,x2​,...,xp​)就对对应(β0,β1,...,βp)(\beta_0,\beta_1,...,\beta_p)(β0​,β1​,...,βp​)个回归系数，这个时候模型的形式就变成了多元线性回归模型：

yi=β0+β1xi1+β2xi2+⋯+βpxip+ϵi，i=1,2,…,ny_i=\beta_0+\beta_1x_{i1}+\beta_2x_{i2}+\cdots+\beta_px_{ip}+\epsilon_i ， i=1,2,\dots,n yi​=β0​+β1​xi1​+β2​xi2​+⋯+βp​xip​+ϵi​，i=1,2,…,n

所以为了简化计算和书写方便，我们可以把它写成矩阵的形式：

Y=Xβ+ϵY=X\boldsymbol{\beta}+\boldsymbol{\epsilon} Y=Xβ+ϵ

Y=[y1y2⋮yn]X=[1x11⋯x1p1x21⋯x2p⋮⋮⋮⋮1xn1⋯xnp]β=[β0β1⋮βp]ε=[ε1ε2⋮εn]Y=\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{bmatrix} X=\begin{bmatrix} 1 & x_{11} & \cdots & x_{1p} \\ 1 & x_{21} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots &\vdots \\ 1 & x_{n1} & \cdots & x_{np} \\ \end{bmatrix} \boldsymbol{\beta}=\begin{bmatrix} \beta_0\\ \beta_1\\ \vdots\\ \beta_p\\ \end{bmatrix} \boldsymbol{\varepsilon}=\begin{bmatrix} \varepsilon_1\\ \varepsilon_2\\ \vdots\\ \varepsilon_n\\ \end{bmatrix} Y=⎣⎢⎢⎢⎡​y1​y2​⋮yn​​⎦⎥⎥⎥⎤​X=⎣⎢⎢⎢⎡​11⋮1​x11​x21​⋮xn1​​⋯⋯⋮⋯​x1p​x2p​⋮xnp​​⎦⎥⎥⎥⎤​β=⎣⎢⎢⎢⎡​β0​β1​⋮βp​​⎦⎥⎥⎥⎤​ε=⎣⎢⎢⎢⎡​ε1​ε2​⋮εn​​⎦⎥⎥⎥⎤​

其中XXX称为设计矩阵(只是习惯叫法)，YYY就不多说了。同样也有一些前提：XXX必须是列满秩；随机误差向量ε\boldsymbol{\varepsilon}ε要满足高斯-马尔科夫条件(1829年，高斯提供了最小二乘法的优化效果强于其他方法的证明，因此被称为高斯-马尔可夫定理)：

(i) E(ε)=0E(\boldsymbol{\varepsilon})=0E(ε)=0

(ii)Var(ε)=σ2I\mathrm{Var(\boldsymbol{\varepsilon)}}=\sigma^2\boldsymbol{I}Var(ε)=σ2I

最小二乘估计

最小二乘估计法(LSE)(LSE)(LSE)，它和机器学习领域的梯度下降法还是有一定的区别的(后者没有这么多假设，实用性更广泛)，准确的来讲LESLESLES只是一种算法，因为随机误差向量ϵ\boldsymbol{\epsilon}ϵ并不能被观测，所以回归方程不存在解，我们只能尽可能的去接近真实值从而解出全局最优解，即确定一个β^\hat{\boldsymbol{\beta}}β^​使得ε=Y−Xβ\boldsymbol{\varepsilon}=Y-X\boldsymbol{\beta}ε=Y−Xβ各元素的平方和达到最小，可以记为：

Q(β)=∑i=1nεi2=εTε=(Y−Xβ)T(Y−Xβ)=(YTY−2βTXTY+βTXTXβ)\begin{aligned} Q(\boldsymbol{\beta}) &=\sum_{i=1}^n\varepsilon_i^2\\ &=\boldsymbol{\varepsilon}^T\boldsymbol{\varepsilon}\\ &=(Y-X\boldsymbol{\beta})^T(Y-X\boldsymbol{\beta})\\ &=(Y^TY-2\boldsymbol{\beta}^TX^TY+\boldsymbol{\beta}^TX^TX\boldsymbol{\beta}) \end{aligned} Q(β)​=i=1∑n​εi2​=εTε=(Y−Xβ)T(Y−Xβ)=(YTY−2βTXTY+βTXTXβ)​

令：

∂Q(β)∂β=−2XTY+2XTXβ=0\frac{\partial{Q(\boldsymbol{\beta})}}{\partial\beta}=-2X^TY+2X^TX\boldsymbol{\beta}=0 ∂β∂Q(β)​=−2XTY+2XTXβ=0

这里需要一些矩阵求导的概念，接下来我们就可以得到一个叫做正规方程的东西：

XTXβ=XTYX^TX\boldsymbol{\beta}=X^TY XTXβ=XTY

由rank(XTX)=rank(X)=p+1\mathrm{rank}(X^TX)=\mathrm{rank}(X)=p+1rank(XTX)=rank(X)=p+1知XTXX^TXXTX是正定矩阵，所以XXX^XXX存在逆矩阵，那么正规方法就有唯一解了：

β^=(β^0,β^1,⋯,β^p)T=(XTX)−1XTY\hat{\boldsymbol{\beta}}=(\hat{\beta}_0,\hat{\beta}_1,\cdots,\hat{\beta}_p)^T=(X^TX)^{-1}X^TY β^​=(β^​0​,β^​1​,⋯,β^​p​)T=(XTX)−1XTY

此时β\boldsymbol{\beta}β的估计就得到了，如果再把它带回到模型中去就有：

Y^=Xβ^=X(XTX)−1XTY=SY\hat{Y}=X\hat{\boldsymbol{\beta}}=X(X^TX)^{-1}X^TY=SYY^=Xβ^​=X(XTX)−1XTY=SY

一般统计学上称SSS是YYY的帽子矩阵，这个称呼是因为有SSS的存在使YYY带上了帽子（总感觉怪怪的？）接下来看残差：

ε^=Y−Y^=(I−H)Y\hat{\boldsymbol{\varepsilon}}=Y-\hat{Y}=(I-H)Y ε^=Y−Y^=(I−H)Y

III是nnn阶的单位矩阵，显然残差的总和为0，是因为Q(β)Q(\boldsymbol{\beta})Q(β)对截距项求偏导数等于0时：

−2∑i=1n[yi−(β0+∑i=1pβixi)]=0-2\sum_{i=1}^n[y_i-(\beta_0+\sum_{i=1}^p\beta_ix_i)]=0 −2i=1∑n​[yi​−(β0​+i=1∑p​βi​xi​)]=0

这个式子很明显表达了当存在截距项时，残差和必然为0，这也是为什么200年前拉普拉斯放弃了最小一乘法。也可以证明最小二乘法得到的估计和最大似然估计的结果是相同的，都是无偏估计。关于最小二乘法的BLUE性质不是本文的重点不再赘述。

补充几个推导过程中用到的矩阵求偏导法则

∂xTa∂x=∂aTx∂x=a\frac{\partial x^Ta}{\partial x}=\frac{\partial a^Tx}{\partial x}=a ∂x∂xTa​=∂x∂aTx​=a

∂xTAx∂x=Ax+ATx\frac{\partial x^TAx}{\partial x}=Ax+A^Tx ∂x∂xTAx​=Ax+ATx

如果AAA是对称的:Ax+ATx=2AxAx+A^Tx=2AxAx+ATx=2Ax.

至此推导过程完毕。

参考文献：梅长林,王宁《近代回归分析方法》[M],科学出版社，.

如果觉得《最小二乘估计矩阵形式的推导》对你有帮助，请点赞、收藏，并留下你的观点哦！