ML学习笔记 #03 正规方程：多元线性回归

前文提到线性回归问题时，我们说在数学上也可以用最小二乘法（Least Square Method）来解决，实际上其思路也是最小化平方误差代价函数（也称残差函数）。

到目前为止，我们都在使用梯度下降算法，但是对于某些线性回归问题，最小二乘法的矩阵解法——正规方程则是更好的解决方案。

正规方程 | Normal Equation

利用多元微分学知识，我们知道对于代价函数：

\[ J(\theta)=J(\theta_0, \theta_1,\cdots, \theta_n) \] 如果它是连续的，则要求出它的最小值，只需要令各偏导为零： \[ \frac{\partial J}{\partial \theta_j}=0,\quad j=0,1,\cdots,n \] 或写作向量形式： \[ \frac{\partial J}{\partial \theta}=\vec 0 \]

就能解出令 \(J(\theta)\) 最小化的 \(\theta\) 值。

由此，我们将代价函数转化为有确定解的代数方程组（其方程式数目正好等于未知数的个数），这个方程组就是正规方程（Normal Equation）。

数学推导 1

下面我们就对多元线性回归的代价函数进行求解：

\[ J(\theta)=\frac{1}{2m}\sum_{i=1}^m\left(\theta^Tx^{(i)}-y^{(i)}\right)^2 \] 于是其偏导函数为： \[ \frac{\partial J}{\partial \theta}=\frac{1}{m}\sum_{i=1}^m\left(\theta^Tx^{(i)}-y^{(i)}\right)x^{(i)} \] 要使之为零向量，只能是： \[ \theta^Tx^{(i)}=y^{(i)},\quad i=1,2,\cdots,m \] 恒成立。写作矩阵为： \[ \theta ^TX^T=y^T\text{ ，或 }X\theta =y \] 其中， \[ X_{(n+1)\times m}=\left[ \begin{matrix} x_{0}^{(1)}& x_{1}^{(1)}& \cdots& x_{n}^{(1)}\\ x_{0}^{(2)}& x_{1}^{(2)}& \cdots& x_{n}^{(2)}\\ \vdots& \vdots& \ddots& \vdots\\ x_{0}^{(m)}& x_{1}^{(m)}& \cdots& x_{n}^{(m)}\\ \end{matrix} \right] =\left[ \begin{array}{c} {x^{(1)}}^T\\ {x^{(2)}}^T\\ \vdots\\ {x^{(m)}}^T\\ \end{array} \right] ,\quad y=\left[ \begin{array}{c} y^{(1)}\\ y^{(2)}\\ \vdots\\ y^{(m)}\\ \end{array} \right] \] 两边同时乘以 \(X^T\)，假设 \(X^TX\) 可逆，解得： \[ \theta=(X^TX)^{-1}X^Ty \] ### 数学推导 2

前面的推导中，在向量形式的偏导函数中发现了简化条件，将零向量提出来单独求解。下面介绍另一个纯矩阵形式的解法。

首先将代价函数表示为： \[ \begin{aligned} J(\theta )&=\frac{1}{2m}\left( X\theta -y \right) ^T\left( X\theta -y \right)\\ &=\frac{1}{2}\left( \theta ^TX^T-y^T \right) \left( X\theta -y \right)\\ &=\frac{1}{2}\left( \theta ^TX^TX\theta -\theta ^TX^Ty-y^TX\theta +y^Ty \right)\\ \end{aligned} \] 接下来对 \(J(\theta )\) 求偏导，需要用到矩阵的求导法则（证明过程略去不表）：

1. 当 \(f(x) = Ax\) 时， \[ \frac{\partial f (x)}{\partial x^T} = \frac{\partial (Ax)}{\partial x^T} =A \]

2. 当 \(f(x) = x^TAx\) 时， \[ \frac{\partial f (x)}{\partial x} = \frac{\partial (x^TAx)}{\partial x} =Ax+A^Tx \]

3. 当 \(f(x) = a^Tx\) 时， \[ \frac{\partial a^Tx}{\partial x} = \frac{\partial x^Ta}{\partial x} =a \]

4. 当 \(f(x) = x^TAy\) 时， \[ \frac{\partial x^TAy}{\partial x} = Ay \]

分别用法则 2、4、3 求导，得到： \[ \begin{aligned} \frac{\partial J\left( \theta \right)}{\partial \theta}&=\frac{1}{2m}\left( 2X^TX\theta -X^Ty-(y^TX)^T+0 \right)\\ &=\frac{1}{2m}\left( 2X^TX\theta -X^Ty-X^Ty+0 \right)\\ &=\frac{1}{m}\left(X^TX\theta -X^Ty\right)\\ \end{aligned} \] 令偏导为零，解得： \[ \theta =\left( X^TX \right) ^{-1}X^Ty \]

梯度下降 vs. 正规方程

观察到在正规方程的结果中，\(X^TX\) 是一个 \((n+1)\times(n+1)\) 的矩阵，因此直接取逆计算 \(\theta\) 的复杂度是 \(O(n^3)\) 。如果 \(n\) 不是很大，这是有效的，但是如果 \(n\) 达到了 \(10^4,10^5\) 或更高，就需要使用梯度下降了。

下面从其他方面对两种算法进行比较：

区别	梯度下降	正规方程
学习率 \(\alpha\)	需要选择	不需要
迭代	需要多次迭代	一次运算得出
\(n\) 的取值	当 \(n\) 大时也能较好适用	当 \(n\) 小于时 \(10^4\) 还是可以接受的
特征缩放	特征取值范围相差大时需要	不需要缩放
适用情形	适用于各种类型的模型	只适用于线性模型

这里提及适用情形，是因为随着问题的深入，算法将越发复杂。例如在分类算法中的逻辑回归等模型，就无法使用正规方程求解。

代码实现

下面仍以 Coursera 上的多元线性回归数据集 ex1data2.txt 为例实现，代码非常简洁：

import numpy as np

# load data, data.shape = (47, 3)
data = np.genfromtxt("ex1data2.txt", delimiter=',')
(m, n) = data.shape
X = np.c_[np.ones(m), data[:, :-1]]
y = data[:, -1]

# Normal Equation
theta = np.linalg.inv(X.T @ X) @ X.T @ y
print(theta)

# predict
predict = np.array([1, 1650, 3])
print(predict @ theta)

很快就计算完了，预测在 $( x_1=1650,x_2=3 ) $ 时的房价为 293081.46433489426。大约相当于 3000 次梯度下降迭代的精度。

不可逆情形

前一节的推导基于 \(X^TX\) 可逆（Invertible）的假设，如若不可逆（Non-invertible，也称 Singular），我们只需将代码中的 inv() 换成 pinv() 求出伪逆矩阵即可。

通常导致矩阵不可逆的原因可能有：

• 存在冗余特征（特征之间不相互独立）；
• 特征数 \(n\) 远大于样本数 \(m\)（样本数不足以刻画这么多特征）。

解决方法对应为：

• 删除冗余特征（线性相关特征只保留其一）；
• 削减非必要的特征，或正则化方法（Regularization），后文将介绍。