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

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

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

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

{∂∂a0F(a0,a1,a2)=2a0+2a1x1ˉ+2a2x2ˉ−2yˉ=0∂∂a1F(a0,a1,a2)=2a1x12ˉ+2a2x1x2ˉ+2a0x1ˉ−2x1yˉ=0∂∂a2F(a0,a1,a2)=2a2x22ˉ+2a1x1x2ˉ+2a0x2ˉ−2x2yˉ=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.

求得方程组的解是：

{a0=yˉ−a1x1ˉ−a2x2ˉa1=(x1yˉ−x1ˉyˉ)(x22ˉ−(x2ˉ)2)−(x2yˉ−x2ˉyˉ)(x1x2ˉ−x1ˉx2ˉ)(x12ˉ−(x1ˉ)2)(x22ˉ−(x2ˉ)2)−(x1x2ˉ−x1ˉx2ˉ)2a2=(x1yˉ−x1ˉyˉ)(x1x2ˉ−x1ˉx2ˉ)−(x2yˉ−x2ˉyˉ)(x12ˉ−(x1ˉ)2)(x1x2ˉ−x1ˉx2ˉ)2−(x22ˉ−(x2ˉ)2)(x12ˉ−(x1ˉ)2) \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.

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

{a0+a1x1ˉ+a2x2ˉ=yˉa1x12ˉ+a2x1x2ˉ+a0x1ˉ=x1yˉa2x22ˉ+a1x1x2ˉ+a0x2ˉ=x2yˉ \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{aligned}{ \left[ \begin{array}{ccc} 1 & \bar{x_1} & \bar{x_2} \ \bar{x_1} & \bar{x_1^2} & \bar{x_1x_2} \ \bar{x_2} & \bar{x_1x_2} & \bar{x_2^2} \end{array} \right ] \left[ \begin{array}{ccc} a_0 \ a_1 \ a_2 \end{array} \right ]

\left[ \begin{array}{ccc} \bar{y} \ \bar{x_1y} \ \bar{x_2y} \end{array} \right ] }\end{aligned} $$

所以：

$$ \begin{aligned}{ \left[ \begin{array}{ccc} a_0 \ a_1 \ a_2 \end{array} \right ]

\left[ \begin{array}{ccc} 1 & \bar{x_1} & \bar{x_2} \ \bar{x_1} & \bar{x_1^2} & \bar{x_1x_2} \ \bar{x_2} & \bar{x_1x_2} & \bar{x_2^2} \end{array} \right ]^{-1} \left[ \begin{array}{ccc} \bar{y} \ \bar{x_1y} \ \bar{x_2y} \end{array} \right ] }\end{aligned} $$

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 就完成了多元线性方程组的求解。

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

y^=a0+a1x1+a2x2+...+anxn \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} \]

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

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

$$ \begin{aligned}{ \left[ \begin{array}{ccc} \hat{y_1} \ \hat{y_2} \ … \ \hat{y_k} \end{array} \right ]

\left[ \begin{array}{ccc} 1 & x_{11} & … & x_{1n} \ 1 & x_{21} & … & x_{2n} \ … & … & … & … \ 1 & x_{k1} & … & x_{kn} \end{array} \right ] \left[ \begin{array}{ccc} a_0 \ a_1 \ … \ a_n \end{array} \right ] }\end{aligned} $$

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

y^=Xa \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=AT+BT (A+B)^T = A^T+B^T 和 (AB)T=BTAT (AB)^T = B^TA^T 可得：

loss=(y−Xa)T(y−Xa)=(yT−aTXT)(y−Xa)=yTy−yTXa−aTXTy+aTXTXa \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}

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

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

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

XTy=XTXa⇒a=(XTX)−1XTy \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}

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

a=X†y \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]]