吴恩达机器学习笔记（4）—多变量线性回归：正规方程（附代码）

到目前为止，我们都在使用梯度下降算法来求得参数 $\theta$ 的最优解，但是对于某些线性回归问题，最小二乘法的矩阵解法——正规方程则是更好的解决方案，不用再通过迭代算法，而是可以直接一次性求解 $\theta$ 的最优解。

一、正规方程

利用多元微分学知识，我们知道对于代价函数：
$$J(\theta )=J({\theta }_0,{\theta }_1,...,{\theta }_n)$$

如果它是连续的，则要求出它的最小值，只需要令各偏导为零：
$$\frac{\partial J}{\partial {\theta }_j}=0,j=0,1,2,...,n$$

或写作向量形式：
$$\frac{\partial J}{\partial \theta }=\vec{0}$$

就能解出令 $J(\theta )$ 最小化的 $\theta$ 值。如果真的做完微积分并且求解出参数 ${\theta }_0$ 到 ${\theta }_n$ ，这个偏微分最终可能很复杂，并且遍历所有的偏微分也很耗费时间。

由此，我们将代价函数转化为有确定解的代数方程组（其方程式数目正好等于未知数的个数），这个方程组就是正规方程。正规方程如下：
$$\theta =(X^{T}X{)}^{-1}{X}^{T}y$$

求解例子：

运用正规方程方法求解 $\theta$ 值：

二、数学推导1

下面我们就对多元线性回归的代价函数进行求解：
$$J(\theta )=\frac{1}{2m}\sum _{i=1}^m{\left({\theta }^T{x}^{(i)}-y^{(i)}\right)}^2$$

于是其偏导函数为：
$$\frac{\partial J}{\partial \theta }=\frac{1}{m}\sum _{i=1}^m\left({\theta }^T{x}^{(i)}-y^{(i)}\right)x^{(i)}$$

要使之为零向量，只能是：
$${\theta }^T{x}^{(i)}=y^{(i)},i=1,2,...,m$$

恒成立。写作矩阵为：
$${\theta }^T{X}^T=y^T或X\theta =y$$

其中，
$$X_{m\times (n+1)}=\left[\begin{array}{cccc}{x}_0^{(1)}& x_1^{(1)}& \cdots & x_n^{(1)}\\ x_0^{(2)}& x_1^{(2)}& \cdots & x_n^{(2)}\\ ⋮& ⋮& \ddots & ⋮\\ x_0^{(m)}& x_1^{(m)}& \cdots & x_n^{(m)}\end{array}\right]=\left[\begin{array}{c}{x}^{(1{)}^T}\\ x^{(2{)}^T}\\ ⋮\\ x^{(m{)}^T}\end{array}\right],y=\left[\begin{array}{c}{y}^{(1)}\\ y^{(2)}\\ ⋮\\ y^{(m)}\end{array}\right]$$

两边同时乘以 $X^T$，假设 $X^{T}X$ 可逆，解得：
$$\theta =(X^{T}X{)}^{-1}{X}^{T}y$$

三、数学推导2

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

首先将代价函数表示为：
$$\begin{array}{rl}J(\theta )& =\frac{1}{2m}(X\theta -y{)}^T(X\theta -y)\\ & =\frac{1}{2}({\theta }^T{X}^T-y^T)(X\theta -y)\\ & =\frac{1}{2}({\theta }^T{X}^{T}X\theta -{\theta }^T{X}^{T}y-y^{T}X\theta +y^{T}y)\end{array}$$

接下来对 $J(\theta )$ 求偏导，需要用到矩阵的求导法则（证明过程略去不表）：

1. 当 $f(x)=Ax$ 时，
   $$\frac{\partial f(x)}{\partial x^T}=\frac{\partial (Ax)}{\partial x^T}=A$$
2. 当 $f(x)=x^{T}Ax$ 时，
   $$\frac{\partial f(x)}{\partial x}=\frac{\partial (x^{T}Ax)}{\partial x}=Ax+A^{T}x$$
3. 当 $f(x)=a^{T}x$ 时，
   $$\frac{\partial a^{T}x}{\partial x}=\frac{\partial x^{T}a}{\partial x}=a$$
4. 当 $f(x)=x^{T}Ay$ 时，
   $$\frac{\partial x^{T}Ay}{\partial x}=Ay$$

分别用法则 2、4、3 求导，得到：
$$\begin{array}{rl}\frac{\partial J(\theta )}{\partial \theta }& =\frac{1}{2m}(2X^{T}X\theta -X^{T}y-(y^{T}X{)}^T+0)\\ & =\frac{1}{2m}(2X^{T}X\theta -X^{T}y-X^{T}y+0)\\ & =\frac{1}{m}(X^{T}X\theta -X^{T}y)\end{array}$$

令偏导为零，解得：
$$\theta =(X^{T}X{)}^{-1}{X}^{T}y$$

四、梯度下降 vs. 正规方程

观察到在正规方程的结果中，$X^{T}X$ 是一个 $(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^{T}X$ 可逆的假设，如若不可逆，我们只需将代码中的 inv() 换成 pinv() 求出伪逆矩阵即可。在 python 中有两个函数可以求解矩阵的逆，一个被称为pinv()，另一个是inv()，这两者之间的差异是技术性的，一个是所谓的伪逆，另一个被称为逆。使用pinv() 函数可以计算出 $\theta$ 的值，即便矩阵是不可逆的。

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

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

解决方法对应为：

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