04. 多元线性回归的求解

来源：原始 Jupyter Notebook。内容已转换为网页阅读格式；下载原文件。

从一元线性回归和二元线性回归的解法中可以看到，线性回归的一般思路是这样的：

定义损失函数： $loss = (y - \hat{y})^2$
将线性回归函数 $\hat{y} = a_0 + a_1x_1 + a_2x_2 + ... + a_nx_n$ 代入损失函数，并求解系数 $a_0$ ， $a_1$ ，…， $a_n$ 使得损失函数最小
损失函数的最小值可以通过对损失函数求导来实现，分别对每个系数求偏导，并令其等于 0
得到一个关于各个系数的线性方程组，求解这个线性方程组得到各个系数的值

在求解二元线性回归时，线性方程组是这样的：

\left\{ \begin{aligned} \frac{\partial}{\partial a_0} F(a_0, a_1, a_2) &=& 2a_0 + 2a_1\bar{x_1} + 2a_2\bar{x_2} - 2\bar{y} &=& 0 \\ \frac{\partial}{\partial a_1} F(a_0, a_1, a_2) &=& 2a_1\bar{x_1^2} + 2a_2\bar{x_1x_2} + 2a_0\bar{x_1} - 2\bar{x_1y} &=& 0 \\ \frac{\partial}{\partial a_2} F(a_0, a_1, a_2) &=& 2a_2\bar{x_2^2} + 2a_1\bar{x_1x_2} + 2a_0\bar{x_2} - 2\bar{x_2y} &=& 0 \\ \end{aligned} \right.

求得方程组的解是：

\left\{ \begin{aligned} a_0 &=& \bar{y} - a_1\bar{x_1} - a_2\bar{x_2} \\ a_1 &=& \frac{(\bar{x_1y}-\bar{x_1}\bar{y})(\bar{x_2^2}-(\bar{x_2})^2) - (\bar{x_2y}-\bar{x_2}\bar{y})(\bar{x_1x_2}-\bar{x_1}\bar{x_2})}{(\bar{x_1^2}-(\bar{x_1})^2)(\bar{x_2^2}-(\bar{x_2})^2) - (\bar{x_1x_2}-\bar{x_1}\bar{x_2})^2} \\ a_2 &=& \frac{(\bar{x_1y}-\bar{x_1}\bar{y})(\bar{x_1x_2}-\bar{x_1}\bar{x_2}) - (\bar{x_2y}-\bar{x_2}\bar{y})(\bar{x_1^2} - (\bar{x_1})^2)}{(\bar{x_1x_2}-\bar{x_1}\bar{x_2})^2 - (\bar{x_2^2}-(\bar{x_2})^2)(\bar{x_1^2}-(\bar{x_1})^2)} \end{aligned} \right.

可以看到，使用这种方法来求解线性回归非常麻烦。当我们扩展到多元线性回归时，这个计算量是极大的。那么有没有一种通用的方法来求解线性回归问题呢？我们不妨把上面的方程组化简：

\left\{ \begin{aligned} a_0 + a_1\bar{x_1} + a_2\bar{x_2} &=& \bar{y} \\ a_1\bar{x_1^2} + a_2\bar{x_1x_2} + a_0\bar{x_1} &=& \bar{x_1y} \\ a_2\bar{x_2^2} + a_1\bar{x_1x_2} + a_0\bar{x_2} &=& \bar{x_2y} \\ \end{aligned} \right.

转换为矩阵形式：

\begin{bmatrix} 1 & \bar{x}_1 & \bar{x}_2 \\ \bar{x}_1 & \overline{x_1^2} & \overline{x_1x_2} \\ \bar{x}_2 & \overline{x_1x_2} & \overline{x_2^2} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \end{bmatrix} = \begin{bmatrix} \bar{y} \\ \overline{x_1y} \\ \overline{x_2y} \end{bmatrix}

所以：

\begin{bmatrix} a_0 \\ a_1 \\ a_2 \end{bmatrix} = \begin{bmatrix} 1 & \bar{x}_1 & \bar{x}_2 \\ \bar{x}_1 & \overline{x_1^2} & \overline{x_1x_2} \\ \bar{x}_2 & \overline{x_1x_2} & \overline{x_2^2} \end{bmatrix}^{-1} \begin{bmatrix} \bar{y} \\ \overline{x_1y} \\ \overline{x_2y} \end{bmatrix}

import numpy as np

X1 = np.array([6000, 6000, 8000, 8000, 10000])
X2 = np.array([60, 80, 70, 100, 90])
Y = np.array([340000, 350000, 400000, 450000, 500000])

x1 = np.sum(X1) / X1.size
x2 = np.sum(X2) / X2.size
y = np.sum(Y) / Y.size
print("bar_x1 = {0}".format(x1))
print("bar_x2 = {0}".format(x2))
print("bar_y = {0}".format(y))

x1y = np.sum(np.multiply(X1, Y)) / Y.size
x2y = np.sum(np.multiply(X2, Y)) / Y.size
x1x1 = np.sum(np.multiply(X1, X1)) / Y.size
x1x2 = np.sum(np.multiply(X1, X2)) / Y.size
x2x2 = np.sum(np.multiply(X2, X2)) / Y.size
print("bar_x1y = {0}".format(x1y))
print("bar_x2y = {0}".format(x2y))
print("bar_x1x1 = {0}".format(x1x1))
print("bar_x1x2 = {0}".format(x1x2))
print("bar_x2x2 = {0}".format(x2x2))

# 不使用科学计数法
np.set_printoptions(suppress=True)

mx = np.matrix([1, x1, x2, x1, x1x1, x1x2, x2, x1x2, x2x2]).reshape(3, 3)
my = np.matrix([y, x1y, x2y]).reshape(3, 1)
#a = np.linalg.inv(mx).dot(my)
a = mx.I * my
print("a = {0}".format(a))

运行结果

bar_x1 = 7600.0
bar_x2 = 80.0
bar_y = 408000.0
bar_x1y = 3188000000.0
bar_x2y = 33280000.0
bar_x1x1 = 60000000.0
bar_x1x2 = 620000.0
bar_x2x2 = 6600.0
a = [[62105.2631579 ]
 [   32.10526316]
 [ 1273.68421053]]

这里使用了矩阵求逆和乘法运算，在计算量上比之前的方法并没有减少，但是从编程的角度来说，确实简化了很多人工计算，直接使用 a = mx.I * my 就完成了多元线性方程组的求解。

我们再回过头来看下线性回归函数：

\hat{y} = a_0 + a_1x_1 + a_2x_2 + ... + a_nx_n

既然线性方程组可以转换为矩阵表示，这里的线性函数当然也可以转换为矩阵表示：

\[ \hat{y} = \left[ \begin{array}{ccc} 1 & x_1 & x_2 & ... & x_n \end{array} \right ] \left[ \begin{array}{ccc} a_0 \\ a_1 \\ a_2 \\ ... \\ a_n \end{array} \right ] = \bf{x}^{T}\bf{a} \]

其中， $\bf{x}$ 和 $\bf{a}$ 都是一个 $n \times 1$ 的矩阵，也就是一个列向量，所以用小写加粗字母表示（一般情况，小写字母 $x$ 表示标量，粗体小写字母 $\bf{x}$ 表示（列）向量，大写字母 $X$ 表示矩阵）。

上面只是对某个样本数据的表示，假设有 k 个样本，那么可以写成矩阵形式：

\begin{bmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_k \end{bmatrix} = \begin{bmatrix} 1 & x_{11} & \cdots & x_{1n} \\ 1 & x_{21} & \cdots & x_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{k1} & \cdots & x_{kn} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{bmatrix}

分别用 $\bf{\hat{y}}$ ， $\bf{a}$ 、 $X$ 代表上面的矩阵，上式可以简写成：

\bf{\hat{y}} = \rm{X}\bf{a}

损失函数仍然和之前一样，使用平方损失：

\[ \begin{aligned} loss &= (y - \hat{y})^2 \\ &= \sum_{i=0}^{k}(y_i - \hat{y_i})^2 \\ &= \left[ \begin{array}{ccc} y_1-\hat{y_1} & y_2-\hat{y_2} & ... & y_k-\hat{y_k} \end{array} \right ] \left[ \begin{array}{ccc} y_1-\hat{y_1} \\ y_2-\hat{y_2} \\ ... \\ y_k-\hat{y_k} \end{array} \right ] \\ &= \left[ \begin{array}{ccc} y_1-\hat{y_1} \\ y_2-\hat{y_2} \\ ... \\ y_k-\hat{y_k} \end{array} \right ]^T \left[ \begin{array}{ccc} y_1-\hat{y_1} \\ y_2-\hat{y_2} \\ ... \\ y_k-\hat{y_k} \end{array} \right ] \\ &= (\bf{y}-\bf{\hat{y}})^T (\bf{y}-\bf{\hat{y}}) \\ &= (\bf{y}-\rm{X}\bf{a})^T (\bf{y}-\rm{X}\bf{a}) \end{aligned} \]

根据矩阵的基本计算法则： $(A+B)^T = A^T+B^T$ 和 $(AB)^T = B^TA^T$ 可得：

\begin{aligned} loss &= (\bf{y}-\rm{X}\bf{a})^T (\bf{y}-\rm{X}\bf{a}) \\ &= (\bf{y}^T-\bf{a}^T\rm{X}^T) (\bf{y}-\rm{X}\bf{a}) \\ &= \bf{y}^T\bf{y} - \bf{y}^T\rm{X}\bf{a} - \bf{a}^T\rm{X}^T\bf{y} + \bf{a}^T\rm{X}^T\rm{X}\bf{a} \end{aligned}

求解线性回归，也就是求解系数 $\bf{a}$ 使得 loss 函数取得最小值，为了求函数的最小值，这里也要用到之前说的求导数的方法，不过这里我们是对矩阵求导数。

矩阵求导（Matrix Derivative） 也叫 矩阵微分（Matrix Differential），在机器学习、图像处理、最优化等领域的公式推导中经常用到。关于矩阵求导的概念和原理这里不做过多的介绍，可以参考这里。我们要求损失函数关于向量 $\bf{a}$ 的导数：

$$ \begin{aligned} \frac{\partial loss}{\partial \bf{a}} &= \frac{\partial \bf{y}^T\bf{y}}{\partial \bf{a}}

\frac{\partial \bf{y}^T\rm{X}\bf{a}}{\partial \bf{a}}
\frac{\partial \bf{a}^T\rm{X}^T\bf{y}}{\partial \bf{a}}

\frac{\partial \bf{a}^T\rm{X}^T\rm{X}\bf{a}}{\partial \bf{a}} \ &= 0 - \rm{X}^T\bf{y} - \rm{X}^T\bf{y} + 2\rm{X}^T\rm{X}\bf{a} \ &= 2\rm{X}^TX\bf{a} - 2\rm{X}^T\bf{y} \end{aligned} $$

令 $\frac{\partial loss}{\partial \bf{a}} = 0$ 可得：

\rm{X}^T\bf{y} = \rm{X}^T\rm{X}\bf{a} \Rightarrow \bf{a} = (\rm{X}^T\rm{X})^{-1}\rm{X}^T\bf{y}

其中， $\rm{X}^T\rm{X}$ 必须是可逆的，也就是说必须是满秩矩阵（又叫 非奇异矩阵）或正定矩阵。可以看出，通过一些简单的矩阵运算（转置、求逆、乘法）就可以求出损失函数的最小值，上面这个公式在统计学中被称为 正规方程（normal equation）。其中 $(\rm{X}^T\rm{X})^{-1}\rm{X}^T$ 称为矩阵 $\rm{X}$ 的 伪逆矩阵（pseudoinverse）或 广义逆（Generalized inverse），通常记为 $\rm{X}^{\dagger}$ 。

\bf{a} = \rm{X}^{\dagger}\bf{y}

通过这种方法求解线性回归，代码实现非常简便：

import numpy as np

X = np.matrix([
    [1, 1, 1, 1, 1], 
    [6000, 6000, 8000, 8000, 10000], 
    [60, 80, 70, 100, 90]]).T
Y = np.matrix([340000, 350000, 400000, 450000, 500000]).T

a = (X.T * X).I * X.T * Y
print("a = {0}".format(a))

运行结果

a = [[62105.26315789]
 [   32.10526316]
 [ 1273.68421053]]

或者直接使用 numpy 自带的 pinv() 函数求伪逆矩阵：

import numpy as np

X = np.matrix([
    [1, 1, 1, 1, 1], 
    [6000, 6000, 8000, 8000, 10000], 
    [60, 80, 70, 100, 90]]).T
Y = np.matrix([340000, 350000, 400000, 450000, 500000]).T

a = np.linalg.pinv(X).dot(Y)
print("a = {0}".format(a))

运行结果

a = [[62105.26315789]
 [   32.10526316]
 [ 1273.68421053]]