最小二乘估计(Least Squares Estimator)的公式的推导

Tags: , ,

最近在学习ML(Machine Learning),注意到了一个有趣的东西:Least Squares Estimator

先从简单说起吧。看下面的式子:

\[ y = ax + e \]

这是一个非常简单的直线方程。如果赋予y、a、x、b具体的意义,这个式子就有意思了:

  1. 假设x是一个统计变量(预先就知道的),譬如,x代表人的年龄。

  2. 假设y是关于x的一个label量(预先就知道的),譬如,y代表的是年龄为x时的人的智商。

  1. 假设y和x存在线性关系,那么可以有 y = ax。这个式子表明年龄为x时,智商为ax。

  2. 当x、y的取值只有一对时,a = y/x,但当x、y不只一对时,y = ax可能会无解(因为求解的是方程组 \( y_{i} = ax_{i} \) 了)

  3. 为了使方程组 \( y_{i} = ax_{i} \) 可以求解,需要把方程组扩展成 \( y_{i} = ax_{i} + e_{i} \) 。

  4. \( y_{i} = ax_{i} + e_{i} \)使得我们有机会求出a,但同时也产生了很多个\( e_{i} \)。每对都有它自己的error系数的话,这个a的意义就减弱了。

  5. 为了使得a变得更有意义,我们希望每个error系数尽可能地小(无限逼近0最好了),同时又能求出唯一的a。

  6. 又因为现实生活中,智商肯定不只跟年龄x有关系,还和其他参数有关系,那么可以再把公式扩展成:

\[ y_{i} = a_{1}x_{i1} + a_{2}x_{i2} + \cdots + a_{k}x_{ik} + e_{i} , 1\le i\le n, k\ge 1 \]

现在,把上式写成矩阵形式:

\[ \vec y = X\vec a + \vec e \]

\[ \left[ \begin{matrix} y_{1}\\ y_{2}\\ \vdots \\ y_{n}\\ \end{matrix} \right] = \left[ \begin{matrix} x_{11}& x_{12}& \cdots & x_{1k}\\ x_{21}& x_{22}& \cdots & x_{2k}\\ \vdots & \vdots & \ddots & \vdots \\ x_{n1}& x_{n2}& \cdots & x_{nk}\\ \end{matrix} \right] \left[ \begin{matrix} a_{1}\\ a_{2}\\ \vdots \\ a_{k}\\ \end{matrix} \right] + \left[ \begin{matrix} e_{1}\\ e_{2}\\ \vdots \\ e_{n}\\ \end{matrix} \right] \]

再回到上面的第7点:为了使得\(\vec a\)变得更有意义,我们希望\(\vec e\)的每个分量尽可能地小。明确这一点非常重要。

那么,这个目标完成情况应该如何衡量?其实很简单,既然\(\vec e\)是一个向量(n维空间),那么\(\vec e\)的长度就是我们需要的指标:

\[ |\vec e| = \sqrt { \sum ^{n}_{i=1}e_{i}^{2} } \]

开根号是不必要的,我们可以换成下面这个指标:

\[ |\vec e|^{2} = \sum ^{n}_{i=1}e_{i}^{2} = \vec e\vec e = \vec e^{T}\vec e \]

小结一下:当\( \vec e^{T}\vec e \)取得最小值时,\(\vec a\)能取得最优解。

继续推导。

由上文可知:

\( \vec e = \vec y - X\vec a \)

\( \vec e^{T} = (\vec y - X\vec a)^{T} \)

\( \vec e^{T}\vec e = (\vec y - X\vec a)^{T}(\vec y - X\vec a) \)

\( = (\vec y^{T} - \vec a^{T}X^{T})(\vec y - X\vec a) \)

\( = \vec y^{T}\vec y - \vec a^{T}X^{T}\vec y - \vec y^{T}X\vec a + \vec a^{T}X^{T}X\vec a \)

注意,中间的2个子项是可以合并的。首先,仔细观察\( \vec a^{T}X^{T}\vec y \)这个子项,发现它是一个,那么就有:

\( \vec a^{T}X^{T}\vec y = (\vec a^{T}X^{T}\vec y)^{T} \)

(一个数值可认为是一个1维的方阵,1维方阵的转置矩阵是它本身)

而又有:

\( (\vec a^{T}X^{T}\vec y)^{T} = \vec y^{T}(\vec a^{T}X^{T})^{T} \)

\( = \vec y^{T}(X\vec a) = \vec y^{T}X\vec a \)

得:

\( \vec a^{T}X^{T}\vec y = \vec y^{T}X\vec a \)

所以上面的方程可变为:

\[ \vec e^{T}\vec e = \vec y^{T}\vec y - 2\vec y^{T}X\vec a + \vec a^{T}X^{T}X\vec a \]

如何让\( \vec e^{T}\vec e \)取得最小值?此时需要使用新的招数:矩阵微分

矩阵微分

矩阵微分公式:

设:

\[ \vec y = A\vec x \]

y是一个\(m \times 1\)的矩阵,A是一个\(m \times n\)的矩阵,x是一个\(n \times 1\)的矩阵。

则有:

\[ \frac {\partial \vec y}{\partial \vec x} = A 【公式1】 \]

这是如何得到的呢?实际上超级简单,上面这个式子指的是,\(\vec y \)的每一个分量对\(\vec x \)的每一个分量的微分,结果显然就是一个\(m \times n\)矩阵。

扩展公式:

设:

\[ \alpha = \vec y^{T}A\vec x \]

则有:

\[ \frac {\partial \alpha }{\partial \vec x} = \vec y^{T}A 【公式2】 \]

\[ \frac {\partial \alpha }{\partial \vec y} = \vec x^{T}A^{T} 【公式3】 \]

设:

\[ \alpha = \vec x^{T}A\vec x \]

且A是对称矩阵,

则有:

\[ \frac {\partial \alpha }{\partial \vec x} = 2\vec x^{T}A 【公式4】 \]

应用矩阵微分公式

再来看下刚才的\( \vec e^{T}\vec e \)方程:

\[ \vec e^{T}\vec e = \vec y^{T}\vec y - 2\vec y^{T}X\vec a + \vec a^{T}X^{T}X\vec a \]

对等号右边的式子求关于\(\vec a\)的微分,得到:

\( \frac {\partial \vec y^{T}\vec y}{\partial \vec a} - 2\frac {\partial \vec y^{T}X\vec a}{\partial \vec a} + \frac {\partial \vec a^{T}X^{T}X\vec a}{\partial \vec a} \)

当这个式子(导数)等于0时, 就得到了\( \vec e^{T}\vec e \)的最小值。

显然,第一个子项为0,所以可把它去掉,得到:

\( - 2\frac {\partial \vec y^{T}X\vec a}{\partial \vec a} + \frac {\partial \vec a^{T}X^{T}X\vec a}{\partial \vec a} = 0\)

\( 2\frac {\vec y^{T}X\vec a}{\partial \vec a} = \frac {\vec a^{T}X^{T}X\vec a}{\partial \vec a} \)

观察左边的式子,和上文的【公式2】是一样的,所以有:

\( 2\frac {\vec y^{T}X\vec a}{\partial \vec a} = 2\vec y^{T}X \)

观察右边的式子,符合上文的【公式4】,所以有:

\( \frac {\vec a^{T}X^{T}X\vec a}{\partial \vec a} = 2\vec a^{T}X^{T}X \)

综上,得:

\( 2\vec y^{T}X = 2\vec a^{T}X^{T}X \)

\( \vec y^{T}X = \vec a^{T}X^{T}X \)

\( (\vec y^{T}X)^{T} = (\vec a^{T}X^{T}X)^{T} \)

\( X^{T}\vec y = X^{T}X\vec a \)

\[ \vec a = (X^{T}X)^{-1}X^{T}\vec y \]

这个东西就是所谓的最小二乘估计(Least Squares Estimator)了。

特殊情况下的最小二乘估计

上文的讨论中没有考虑到一种特殊情况:X是一个方阵。当X是方阵时,上面的公式可进一步简化:

\[ \vec a = (X^{T}X)^{-1}X^{T}\vec y \]

\[ \vec a = X^{-1}(X^{T})^{-1}X^{T}\vec y \]

\[ \vec a = X^{-1}\vec y \]

还是搞不明白X不是方阵时,公式为什么要那么复杂?答案就在于,X不是方阵时,\( X^{-1} \)和\( (X^{T})^{-1} \)都不成立,导致最小二乘估计公式不能简化。

关于这个问题的进一步讨论,请阅读我的另一篇文章线性代数之伪逆矩阵

实例

使用ML中常用的Iris数据集http://archive.ics.uci.edu/ml/datasets/Iris来验证下上面的公式是否可用。

Iris鸢尾花卉数据集,是一类多重变量分析的数据集。通过花萼长度,花萼宽度,花瓣长度,花瓣宽度4个属性预测鸢尾花卉属于(Setosa,Versicolour,Virginica)三个种类中的哪一类。

整个数据集可在这里浏览。

首先是为3个类别各赋予1个标签值:

  • Iris-setosa = 1
  • Iris-versicolor = 2
  • Iris-virginica = 3

然后从整个数据集中挑选训练数据集(Train Dataset),譬如从3个类别中各取出前10个数据项。

5.1,  3.5,  1.4,  0.2,  1
4.9,  3.0,  1.4,  0.2,  1
4.7,  3.2,  1.3,  0.2,  1
4.6,  3.1,  1.5,  0.2,  1
5.0,  3.6,  1.4,  0.2,  1
5.4,  3.9,  1.7,  0.4,  1
4.6,  3.4,  1.4,  0.3,  1
5.0,  3.4,  1.5,  0.2,  1
4.4,  2.9,  1.4,  0.2,  1
4.9,  3.1,  1.5,  0.1,  1
7.0,  3.2,  4.7,  1.4,  2
6.4,  3.2,  4.5,  1.5,  2
6.9,  3.1,  4.9,  1.5,  2
5.5,  2.3,  4.0,  1.3,  2
6.5,  2.8,  4.6,  1.5,  2
5.7,  2.8,  4.5,  1.3,  2
6.3,  3.3,  4.7,  1.6,  2
4.9,  2.4,  3.3,  1.0,  2
6.6,  2.9,  4.6,  1.3,  2
5.2,  2.7,  3.9,  1.4,  2
6.3,  3.3,  6.0,  2.5,  3
5.8,  2.7,  5.1,  1.9,  3
7.1,  3.0,  5.9,  2.1,  3
6.3,  2.9,  5.6,  1.8,  3
6.5,  3.0,  5.8,  2.2,  3
7.6,  3.0,  6.6,  2.1,  3
4.9,  2.5,  4.5,  1.7,  3
7.3,  2.9,  6.3,  1.8,  3
6.7,  2.5,  5.8,  1.8,  3
7.2,  3.6,  6.1,  2.5,  3

此时,已经得到了\( X \)、\( \vec y \)的值了,\( X \)是上面这个表格的前4列(30x4的矩阵),\( \vec y \)是第5列(30x1的矩阵)。

我们的目标是求出系数\( \vec a \),它是一个4x1的矩阵。

根据前文推导出来的公式:

\( \vec a = (X^{T}X)^{-1}X^{T}\vec y \)

因为矩阵比较庞大的关系,只能直接给出\( (X^{T}X)^{-1}X^{T} \)的结果了,读者们也可以自己做下计算:

X` = 
5.1  4.9  4.7  4.6  5  5.4  4.6  5  4.4  4.9  7  6.4  6.9  5.5  6.5  5.7  6.3  4.9  6.6  5.2  6.3  5.8  7.1  6.3  6.5  7.6  4.9  7.3  6.7  7.2
3.5  3  3.2  3.1  3.6  3.9  3.4  3.4  2.9  3.1  3.2  3.2  3.1  2.3  2.8  2.8  3.3  2.4  2.9  2.7  3.3  2.7  3  2.9  3  3  2.5  2.9  2.5  3.6
1.4  1.4  1.3  1.5  1.4  1.7  1.4  1.5  1.4  1.5  4.7  4.5  4.9  4  4.6  4.5  4.7  3.3  4.6  3.9  6  5.1  5.9  5.6  5.8  6.6  4.5  6.3  5.8  6.1
0.2  0.2  0.2  0.2  0.2  0.4  0.3  0.2  0.2  0.1  1.4  1.5  1.5  1.3  1.5  1.3  1.6  1  1.3  1.4  2.5  1.9  2.1  1.8  2.2  2.1  1.7  1.8  1.8  2.5

X`X = 
1051.290    532.720     723.350     230.480   
532.720     281.280     345.450     108.190   
723.350     345.450     550.410     182.580   
230.480     108.190     182.580     62.220


(X`X)^-1 = 
0.541034    -0.573805   -0.641245   0.875298
-0.573805   0.633743    0.631976    -0.830927
-0.641245   0.631976    0.920228    -1.42389
0.875298    -0.830927   -1.42389    2.396868


((X`X)^-1)X` = 
0.028273    0.206968    0.048125    -0.076847   -0.083211   -0.056253   -0.097334   -0.032575   -0.006168   -0.002067   0.162628    0.053786    0.125186    0.228843    0.273287    -0.270475   -0.098417   0.033124    0.094950    -0.011335   -0.144267   -0.018560   0.174707    -0.270956   0.001741    -0.003648   -0.181042   -0.178793   0.046731    0.106397
0.010276    -0.191835   -0.013523   0.106879    0.131031    0.115039    0.150711    0.067480    0.031694    0.017830    -0.181668   -0.046873   -0.144359   -0.250620   -0.294553   0.267479    0.117184    -0.036068   -0.122374   0.028729    0.190919    0.027353    -0.189075   0.26628 0.008918    -0.033594   0.204029    0.134858    -0.090344   -0.072185
-0.054892   -0.242631   -0.080010   0.104963    0.072430    -0.003185   0.060144    0.038057    0.014794    0.054978    -0.134766   -0.076454   -0.092183   -0.243448   -0.301346   0.404405    0.092525    -0.012496   -0.017495   -0.032696   0.007320    -0.025114   -0.217735   0.383162    -0.067400   0.105802    0.158253    0.386076    0.057919    -0.288185
0.041703    0.282107    0.083251    -0.205964   -0.12892    0.024128    -0.073167   -0.105123   -0.072449   -0.183062   0.131452    0.130738    0.081924    0.323375    0.408249    -0.628974   -0.084976   -0.007234   -0.066687   0.110491    0.221148    0.125436    0.354307    -0.554733   0.211204    -0.204767   -0.121187   -0.676157   -0.157020   0.617249

a = ((X`X)^-1)X`y = 
-0.291   
0.441   
0.637   
-0.094  

a的值为:

\[ \vec a = \left[ \begin{matrix} -0.291\\ 0.441\\ 0.637\\ -0.094\\ \end{matrix} \right] \]

然后就是测试a的可靠度如何。方法就是把a用到剩余的其他数据项里,算出predict值,并和实际的值做比较,看预测正确的有多少个,公式为:

\[ \vec y_{predict} = X_{test}\vec a \]

结果是:

0.997   
1.103   
0.809   
0.763   
0.822   
1.200   
0.939   
0.923   
1.072   
1.119   
0.992   
1.066   
0.867   
1.007   
1.294   
0.868   
1.026   
0.967   
0.859   
1.044   
0.971   
0.846   
1.241   
1.125   
0.887   
0.702   
0.752   
0.887   
0.852   
0.952   
0.888   
0.505   
0.940   
1.051   
1.364   
0.790   
1.192   
0.946   
1.026   
0.873   
1.563   
2.140   
1.678   
2.366   
1.820   
2.089   
2.419   
2.021   
1.892   
1.854 
2.583   
1.886   
2.250   
2.341   
2.033   
2.074   
2.182   
2.398   
2.258   
1.623   
1.775   
1.721   
1.874   
2.543   
2.477   
2.470   
2.270   
1.862   
2.183   
1.928   
2.236   
2.346   
1.894   
1.567   
2.114   
2.227   
2.173   
2.092   
1.426   
2.066   
2.580   
2.526   
2.650   
2.441   
2.570   
2.709   
2.766   
3.496   
3.085   
2.268   
2.818   
2.539   
3.074   
2.310   
2.939   
2.969   
2.319   
2.500   
2.742  
2.772   
2.788   
3.266   
2.733   
2.509   
2.807   
2.752   
3.008   
2.839   
2.465   
2.602   
2.759   
2.392   
2.573   
2.975   
2.902   
2.470   
2.276   
2.556   
2.919   
2.686 

再四舍五入一下,得到:

1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
2
3
2
2
2
2
2
2
2
2
2
2
2
2
3
2
2
2
2
2
2
2
2
2
2
2
2
2
2
1
2
3
3
3
2
3
3
3
3
3
2
3
3
3
2
3
3
2
3
3
3
3
3
3
3
3
3
3
3
2
3
3
2
3
3
3
2
2
3
3
3

此时就可以算准确率了,经过比较,上面的predict值正确的有109个,总共的测试项有120,准确率高达90.83%哦。

参考资料

https://economictheoryblog.com/2015/02/19/ols_estimator/

http://www.atmos.washington.edu/~dennis/MatrixCalculus.pdf

(未经授权禁止转载)
Written on May 6, 2016