梯度下降和正规方程

监督学习：给定一组样本（训练集），交给机器训练，让机器通过训练集来预测数据。

符号说明：

• $m$ 训练样本的数量
• $x$ 输入变量/特征
• $y$ 输出变量/预测的目标变量
• $(x,y)$ 一组训练样本
• $(x^i,y^i)$ 第$i$组训练样本

---

线性回归模型

监督学习的目的：根据训练数据得到一个假设函数 $h$ ，这个假设函数可以实现：输入一个变量/特征，输出预测的目标变量。

在一开始，认为假设函数 $h(x) = \theta_0+\theta_1 x$ ，也就是预测的目标变量是关于输入变量的线性函数。我们先拟合线性函数，在后续再对模型加以改进，引入更加复杂的模型。

这种线性模型被称作线性回归，上述的模型只有一个变量 $x$ ，也被称为一元线性回归（单变量线性回归）。

---

代价函数

为了让我们的线性回归模型 $h(x) = \theta_0+\theta_1 x$ 更加拟合训练数据，我们需要合理选择两个参数 $\theta_0,\theta_1$ ，使得这条直线上面的点更加接近数据集，也就是我们要最小化如下的式子：

$\min \sum_{i=1}^m(h(x^i)-y^i)^2$

这表示我们对于相同的数据，线性回归预测的数据和真实值的误差要尽可能地小。由于有 $m$ 个训练数据，所以要对 $m$ 求和，平方是为了处理负数的情况。

我们也可以通过误差的平均值来刻画拟合程度，也就是把上式除以 $m$ ，但为了后续计算方便（主要是为了方便求导），我们除以 $2m$ ，也就是下述式子：

$\min \frac{1}{2m}\sum_{i=1}^m(h(x^i)-y^i)^2$

我们称如下关于 $\theta_0,\theta_1$ 的二元函数 $J(\theta_0,\theta_1)$ 为代价函数，我们的目的是为了最小化这个函数：

$J(\theta_0,\theta_1) = \frac{1}{2m}\sum_{i=1}^m(h(x^i)-y^i)^2$

这个代价函数对于大多数线性回归模型都是十分合理的。

---

梯度下降法

梯度下降法是一种最小化函数的方法，它不仅仅实用于上面提到的代价函数 $J(\theta_0,\theta_1)$ ，也适用于更一般的多元函数 $J(\theta_0,\theta_1,\dots,\theta_n)$ 。下面就以 $J(\theta_0,\theta_1)$ 为例来介绍梯度下降法。

梯度下降法的步骤如下：

• 选定 $\theta_0,\theta_1$ 的初值，通常情况下令 $\theta_0=0,\theta_1=0$ ；
• 不断改变 $\theta_0,\theta_1$ 的值，使得函数 $J(\theta_0,\theta_1)$ 向着减少的方向移动，直到找到最小值为止（再改变 $\theta_0,\theta_1$ 的值也不能使函数值减少）。

其对应的数学逻辑，是下面的这个公式：

$\theta _j : = \theta _j -\alpha \frac{\partial }{\partial \theta_j}J(\theta_0,\theta_1)\quad \text{(for j=0 and j=1)}$

这个公式表示每次操作都将 $\theta_0,\theta_1$ 按照后续设定的方向进行更新。其中 $\alpha$ 是一个正参数，称作学习率，它决定梯度下降的速率，当 $\alpha$ 越大，表示梯度下降越迅速。后面的一项是一个导数项，我们在后续会进行推导以及确定 $\alpha$ 的值。

需要注意的是，这个公式对于 $\theta_0,\theta_1$ 的更新是同时的。换句话说，下面的这个更新步骤是不正确的：

$\begin{aligned} \theta _0 &: = \theta _0 -\alpha \frac{\partial }{\partial \theta_j}J(\theta_0,\theta_1) \\ \theta _1 &: = \theta _1 -\alpha \frac{\partial }{\partial \theta_j}J(\theta_0,\theta_1) \end{aligned}$

而应该是：

$\begin{aligned} temp_0 &: = \theta _0 -\alpha \frac{\partial }{\partial \theta_j}J(\theta_0,\theta_1) \\ temp_1 &: = \theta _1 -\alpha \frac{\partial }{\partial \theta_j}J(\theta_0,\theta_1)\\ \theta_0 &: = temp_0 \\ \theta_1 &: = temp_1 \end{aligned}$

下面的临时变量 $temp_0,temp_1$ 确保了 $\theta_0,\theta_1$ 的更新是同时的，而不像第一个步骤一样，更新 $\theta_0$ 后会使用新的 $\theta_0$ 的值来更新 $\theta_1$。

为了理解后面的导数项，我们从更简单的一元函数来讨论，那么上面的梯度下降公式就变成：

$\theta _1 : = \theta _1 -\alpha \frac{\mathrm{d} }{\mathrm{d} \theta_1}J(\theta_1)$

由于函数变为一元的，我们把偏导数符号 $\partial$ 变为 $\mathrm{d}$。

对于函数上的一点 $(\theta,f(\theta))$ ，导数项 $k=\frac{\mathrm{d} }{\mathrm{d} \theta}J(\theta)$ 反映的是函数在这一点上切线的斜率。显然，可以分为以下三种情况：

• $k > 0$ 则说明在该点的很小的去心邻域 $(\theta-\epsilon,\theta+\epsilon)$ 是局部递增的，那么 $\theta$ 应该向减小的方向移动；
• $k < 0$ 则说明在该点的很小的去心邻域 $(\theta-\epsilon,\theta+\epsilon)$ 是局部递减的，那么 $\theta$ 应该向增大的方向移动；
• $k = 0$ 则说明该函数在 $\theta$ 处取得极小值，此时不需要进行移动。

关于学习率 $\alpha$ 的选择，假如 $\alpha$ 过小，则梯度下降的速率过慢，无法快速收敛到极小值；假如 $\alpha$ 过大，则可能越过最小值点，导致离极小值越来越远，无法收敛到极小值。

注意：此处使用“极小值”的叙述而不是“最小值”。根据我们的梯度下降算法的步骤，从一个起点开始梯度下降，最终只能走到一个“山谷”，这是一个局部性质。为了得到整个函数的“最小值”，我们需要取多个不同的起点进行梯度下降，将最后得到的“极小值”进行比较，得到“最小值”。

对于多元函数，我们采用“隔离法”的策略：以其中一个变量为主体求偏导，就可以来以此更新这个变量的值了。

---

梯度下降法的三种变体

批量梯度下降法（Batch Gradient Descent, BGD）

由于梯度下降法每次寻找的只是极小值，为了寻找到最小值，我们需要选取多个不同的起点进行梯度下降。最稳妥的方法，每次迭代都选择所有样本来计算梯度，这样最小的极小值一定是最小值。这就是 批量梯度下降法。

这个做法哪怕不是凸函数，只要收敛，一定能收敛到全局最优。代价是计算效率低：每次更新都需要计算所有样本的梯度

随机梯度下降法（Stochastic Gradient Descent, SGD）

我们当然可以每次只使用一个训练样本来更新参数，这就是 随机梯度下降法。每次只需要计算一个样本的梯度，可以大大的节省时间和空间，但是收敛不稳定甚至可能不收敛。对于凸函数，随机梯度下降法可以达到全局最优。

小批量梯度下降法（Mini-batch Gradient Descent）

为了在稳定性和时间复杂度之间取得平衡，实际工程中采用的比较多的方法是 小批量梯度下降法，将样本随机分成许多大小较小的 批量。每次迭代时，选取一个批量来计算函数梯度，以此估计用全样本计算的结果。假如批量的大小为 $B\ll N$，可以看出，当 $B=1$ 时就退化成SGD算法，当 $B=N$ 时就退化成GD算法。对于MBGD，反复随机抽样进行迭代大概率可以在平均意义上消除梯度估计的偏差，并且还可以通过调整 $B$ 来控制随机性。

由于每一个批量都是独立的，在计算梯度的时候可以进行并行化操作，进一步加快梯度下降的速度。

---

一元线性回归的梯度下降

我们现在就可以使用梯度下降法来求出我们的线性回归模型中的参数 $\theta_0,\theta_1$ 了。

首先来处理后面的偏导数，这需要用到一点微积分的知识：

$\begin{aligned} \frac{\partial }{\partial \theta_j}J(\theta_0,\theta_1) &= \frac{\partial }{\partial \theta_j} \frac{1}{2m}\sum_{i=1}^m(h(x^i)-y^i)^2\\ &=\frac{\partial }{\partial \theta_j} \frac{1}{2m}\sum_{i=1}^m(\theta_0+\theta_1 x^i-y^i)^2\\ \frac{\partial }{\partial \theta_0}J(\theta_0,\theta_1) &= \frac{1}{m}\sum_{i=1}^m(h(x^i)-y^i)\\ \frac{\partial }{\partial \theta_1}J(\theta_0,\theta_1) &= \frac{1}{m}\sum_{i=1}^m(h(x^i)-y^i)x^i\\ \end{aligned}$

由此得到我们的线性回归算法，我们只需要反复执行下述操作，直到我们找到极小值：

$\begin{aligned} temp_0 &: = \theta _0 -\alpha \frac{1}{m}\sum_{i=1}^m(h(x^i)-y^i) \\ temp_1 &: = \theta _1 -\alpha \frac{1}{m}\sum_{i=1}^m(h(x^i)-y^i)x^i\\ \theta_0 &: = temp_0 \\ \theta_1 &: = temp_1 \end{aligned}$

在应用梯度下降法求解线性回归问题的时候，我们不需要取多个不同的起点进行梯度下降。这是因为线性回归问题的代价函数是一个“凸函数”，“极小值”和“最小值”是相同的。

---

多元线性回归

在更多的情况下，我们会有不止一种的变量特征，对此，我们引入多变量的回归分析。

符号说明

• $n$ 变量/特征的个数；
• $x^i$ 一个训练样本的一组特征向量，其维数应该为 $n$ ；
• $x_j$ 一个特征在不同训练样本中的特征量组成的向量，其维数应该为 $m$ ；
• $x^i_j$ 表示在第 $i$ 组特征向量中的第 $j$ 个特征量。

既然如此，我们的线性回归方程也要做出相应的改变，从一元线性函数变成 $n$ 元线性函数，如下：

$h(x) = \theta_0 +\theta_1 x_1+\theta_2 x_2+\dots +\theta_n x_n$

这个函数输入一个特征列向量 $x$ ，表示一个训练样本的特征值。为了方便起见，令 $x_0 = 1$，称其为偏置项：

$x = \begin{pmatrix} x_0 \\ x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix}$

为了简化表示，我们定义如下的向量：

$\theta = \begin{pmatrix} \theta_0 \\ \theta_1 \\ \theta_2 \\ \vdots \\ \theta_n \end{pmatrix}$

那么就有：

$h(x)=\theta^Tx$

---

多元线性回归的梯度下降

我们有公式：

$h(x) = \theta_0 +\theta_1 x_1+\theta_2 x_2+\dots +\theta_n x_n$

其中 $x$ 是一个列向量。对于 $m$ 个数据，用这个公式计算其预测值。故其代价函数为：

$J(\theta) = \frac{1}{2m} \sum_{i=1}^{m}(h(x^i)-y^i)^2$

每个 $x^i$ 都是一个 $n+1$ 维的列向量，表示其特征（包含偏置项）。

梯度下降公式也差不多：

$\begin{aligned} \theta _j &: = \theta _j -\alpha \frac{\partial }{\partial \theta_j}J(\theta)\quad \text{(for j=0,\dots ,n)} \\ &:=\theta_j - \alpha \frac{1}{m} \sum_{i=1}^m(h(x^i)-y^i)x_j^i \quad \text{(for j=0,\dots ,n)} \end{aligned}$

---

特征缩放

当我们的不同的特征变量的取值范围差别很大的时候，梯度下降法就会变得更加不容易收敛，如下面的特征变量及其取值范围。

特征	取值范围
房屋面积	30 ~ 300 平米
卧室数量	1 ~ 5 个

这两个特征放在一起线性回归，会使得面积特征的影响更大，而卧室数量特征值几乎起不到作用。

我们可以把每个特征直接压缩到 $[-1,1]$ 的范围中，也就是进行下面的变换：

$x^{\prime} = \frac{x}{max \left | x\right |}$

再用这些特征进行线性回归，就能较快的收敛了。

还有一种特征缩放的方式，我们称其为归一化，它的变换如下：

$x^{\prime} = \frac{x-\mu}{s}$

其中 $\mu$ 表示 $x$ 的平均值，$s$ 表示 $x$ 的最大值和最小值之差，也就是极差。这种特征缩放的好处是：可以让取值范围中心对称。

---

学习率

这里重点探讨学习率 $\alpha$ 的选取，它直接影响了梯度下降法的效果。

遗憾的是，我们没有一种普适的算法来直接求出最适合一个问题的学习率，学习率的选取通常是进行不断的尝试。

在梯度下降法中，我们需要保证每次迭代后，代价函数 $J(\theta)$ 的值要减小。换句话说，假如建立一个 $J(\theta)$ 的值关于迭代次数的函数关系，这个函数应当严格单调递减，并且最终会收敛到一个最小值。

假如这个函数没有单调递减，我们可以认定学习率 $\alpha$ 的选取太大了，应该选取更小的值。原因是小的 $\alpha$ 只会影响梯度下降的速率，并不会影响梯度下降的最终结果。在实际应用中，可以设定不同的迭代次数，画出迭代次数和 $J(\theta)$ 的函数图像，从而更加直观地判断学习率 $\alpha$ 的选取是否合理。

在实际应用中，我们可以尝试选取几个 $\alpha$ 的值，判断其是否合理。经验地，我们可以按数量级取，如$0.001,0.01,0.1,1,\dots$

---

特征和多项式回归

有时候，对相同的问题进行线性回归，特征的选取对问题的复杂程度有着很大的影响。

举个简单的例子，现在需要预测房屋的售价，可以以房屋的临街距离和房屋的宽度为特征来建立二元线性回归。但事实上，临街距离和宽度这两个特征综合起来体现的就是房屋的面积，于是我们可以只以房屋的面积为特征来建立一元线性回归，这样就简化了问题。

因此，我们可以寻找特征之间的关系，将多个特征综合成一个特征，从而减少问题的维度。

这么做的一个直接的好处就是方便我们接下来的多项式回归：在现实生活中，绝大部分的问题并没有线性回归那么简单，而是呈现出非线性结构。我们可以考虑使用多项式回归，如二次函数、三次函数等。

假如特征 $x$ 和目标变量 $y$ 之间的关系明显为非线性关系，我们考虑使用三次函数来拟合：

$h(x) = \theta_0+\theta_1 x+\theta_2 x^2+\theta_3 x^3$

我们只需要做如下换元，即可将其变为三元线性回归：

$\begin{cases} x_1 = x \\ x_2 = x^2 \\ x_3 = x^3 \end{cases}$

这就变为：

$h(x) = \theta_0+\theta_1 x_1+\theta_2 x_2+\theta_3 x_3$

然后使用多元线性回归的方法来解决就可以了。

需要注意的是，这么换元就体现了特征缩放的重要性了。举个例子，假如 $x$ 的范围是 $[1,1000]$ ，那么 $x^2 \in [1,1000000]$ ，$x^3 \in [1,1000000000]$ ，是非常恐怖的。这时可以运用特征缩放的技巧，将 $x$ 的范围控制在 $[0,1]$ 就不会出现指数造成的范围扩大了。

---

正规方程

相比于梯度下降法，有一种更加直接的方法来求出代价函数取最小值时参数向量 $\theta$ 的值，这就是下面所讲的正规方程。

假如 $\theta \in \mathbb{R}^{n+1}$，那么如下方法可以求出当 $J(\theta_0,\theta_1,\dots,\theta_n)$ 取最小值时，各参数的值：只需对每个参数求偏导，再令导函数为0，解这 $n+1$ 个方程，就能得到这 $n+1$ 个参数。

$\frac{\partial}{\partial \theta_j}J(\theta) = \dots = 0 \quad \text{for every j}$

有些地方也称这种方法为拉格朗日乘数法

下面运用高等数学和线性代数的知识来尝试使用正规方程法来求出参数向量 $\theta$ 的值：

首先，目标是最小化下面的函数：

$J(\theta) = \frac{1}{2m} \sum_{i=1}^{m}(h(x^i)-y^i)^2$

写成矩阵形式：

$J(\theta) = \frac{1}{2m}(X\theta-y)^T(X\theta-y)$

其中矩阵 $X$ 定义如下：

$X = \begin{pmatrix} (x^1)^T \\ (x^2)^T \\ \vdots \\ (x^m)^T \\ \end{pmatrix}$

每个 $x^i$ 都是一个 $n+1$ 维的列向量，表示其特征（包含偏置项）。$y$ 是一个列向量，表示训练数据的预测值。

对这个矩阵形式的函数来求导：

$\nabla J(\theta) = \frac{1}{m}X^T(X\theta-y)$

把导函数置为0：

$\begin{aligned} \frac{1}{m}X^T(X\theta-y) = 0 \\ X^TX\theta = X^Ty \\ \theta = (X^TX)^{-1}X^Ty \end{aligned}$

在笔者写这篇文章的时候，虽然有简单的线性代数基础，却没有多元函数微积分学的高等数学基础，因此这里以结论的形式呈现。需要了解矩阵求导相关知识的参见。

在 MATLAB/Octave 等语言中，可以很方便地对矩阵进行计算，因此正规方程的实现是很简单的。

---

梯度下降法和正规方程的比较

梯度下降法和正规方程都能求出代价函数取最小值时参数 $\theta$ 的值，什么情况下选择什么方法是需要我们考虑的。

一般情况下，我们会倾向于使用正规方程法，它不需要选择学习率 $\alpha$ ，也不需要进行多次迭代，并且实现相对简单。但是矩阵求逆、矩阵乘法的时间复杂度是很大的，为 $O(n^3)$ ，在训练数据集较为庞大时，梯度下降法会显得更加迅速。

---

正规方程在矩阵不可逆情况下的解决方法

在正规方程 $\theta = (X^TX)^{-1}X^Ty$ 中，很自然而然地会引出一个问题：假如我的 $X^TX$ 不可逆怎么办？

一般来说，矩阵 $X^TX$ 不可逆一般都是由以下情况造成的：

• 出现了线性相关的特征。比如有两个特征，$x_1$ 以米为单位，$x_2$ 以厘米为单位，显然 $x_2 = 100x_1$，它们是线性相关的，所以矩阵 $X$ 不是满秩，故矩阵 $X^TX$ 不可逆；
• 特征太多了。例如当 $m \le n$ 时，会限制 $r(X) \le m \le n+1$ ，这就有 $r(X^TX) \le r(X) \le n+1$，直接导致矩阵 $X^TX$ 不是满秩，也就是不可逆。

当我们发现 $X^TX$ 不可逆时，我们需要检查是否存在线性相关的特征；如果确实不存在线性相关的特征，就尝试删除某些无用的数据，或使用后面会介绍的一种叫做正则化的方法来处理数据集。

在 MATLAB/Octave 编程语言中，正规方程可以使用 pinv(X'*X)*X'*y 来进行，这里哪怕 $X^TX$ 不可逆，也不会影响最终结果的计算。因为 pinv 函数计算的是所谓“伪逆”，这里不再具体叙述。