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

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

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

\[ y = ax + e \]

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

假设x是一个统计变量（预先就知道的），譬如，x代表人的年龄。
假设y是关于x的一个label量（预先就知道的），譬如，y代表的是年龄为x时的人的智商。

假设y和x存在线性关系，那么可以有 y = ax。这个式子表明年龄为x时，智商为ax。
当x、y的取值只有一对时，a = y/x，但当x、y不只一对时，y = ax可能会无解（因为求解的是方程组 \( y_{i} = ax_{i} \) 了）
为了使方程组 \( y_{i} = ax_{i} \) 可以求解，需要把方程组扩展成 \( y_{i} = ax_{i} + e_{i} \) 。
\( y_{i} = ax_{i} + e_{i} \)使得我们有机会求出a，但同时也产生了很多个\( e_{i} \)。每对都有它自己的error系数的话，这个a的意义就减弱了。
为了使得a变得更有意义，我们希望每个error系数尽可能地小（无限逼近0最好了），同时又能求出唯一的a。
又因为现实生活中，智商肯定不只跟年龄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】 \]

（MatrixCalculus.pdf的Proposition 5）

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

扩展公式：

设：

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

其中的\( \vec y\)是m x 1，\( \vec x\)是n x 1, \( A \)是m x n, 因此\( \alpha \) 是一个标量(scalar)。

可以得到：

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

（MatrixCalculus.pdf的Proposition 7）

把\( \vec y^{T}A \)看成一个A'，就是公式1了，很好理解。

然后还有：

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

（MatrixCalculus.pdf的Proposition 7）

这个是基于\( \alpha \) 是一个标量(scalar)的事实，标量转置后不变：

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

再应用公式1就得到了公式3。

再设：

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

且A是对称矩阵，

则有：

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

（MatrixCalculus.pdf的Proposition 8和9）

Proposition 9的证明基于Proposition 8，Proposition 8的证明比较复杂建议看原文。Proposition 9只是基于A是对称矩阵A的转置等于A这个事实。

公式2、3、4都可以在文末的MatrixCalculus.pdf链接里找到推导过程。下文将会用到这几条公式，敬请注意。

应用矩阵微分公式

再来看下刚才的\( \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 \]

结果是：

再四舍五入一下，得到：