从线性回归说起

机器学习要解决的问题,可以分为回归(Regression)和分类(Classification)两种类型,为了解决这些问题,数据科学领域的专家们研究出了很多算法,线性回归很可能是其中最简单的一种,而工作、生活中许多简单问题往往用不上多么复杂的模型,让我们从最简单的线性回归说起。

线性回归:快速过一遍

现在我们有一堆数据 \((X, y)\)\(X\) 是研究对象的特征,\(y\) 是想要知道的结果,比如我们想要根据身高(\(X\))来预测体重(\(y\)),假设他们之间有较强的线性关系,那么我们可以用一条直线来拟合这些数据点:

\[ h_\theta=\theta_0+\theta_1x_1=X^T\theta \]

其中 \(\theta\)\(2\times1\) 阶矢量

\[ \theta=\left[ \begin{matrix} \theta_0\\ \theta_1 \end{matrix}\right] \]

\(X^T\) 表示\(X\)的转置矩阵(\(m\times2\), \(m\) 是训练的数据量,\(x_i^j\) 表示的是第 \(j\) 个数据的第 \(i\) 个特征),其中\[X^T\] 第一列的所有数据都是1(因为 \(\theta_0\) 的系数是1):

\[ X^T=\left[\begin{matrix} x_0^1, x_1^1\\ x_0^2, x_1^2\\ ...\\ x_0^m, x_1^m\\ \end{matrix}\right] \]

在已知 \(\theta\) 时,我们可以根据下式预测某个数据 \(x^{(i)}\) 对应的 \(y^{(i)}\)

\[ y^{(i)}=h_\theta(x^{(i)})=\theta_0+\theta_1x_1^{(i)}=\theta^Tx^{(i)} \]

既然是拟合的直线,其结果和实际的数据就会存在偏差,我们用损失函数(Cost Function)来表示预测值(\(h_\theta\))与实际值(\(y\))的偏离程度:

\[ J(\theta)=\frac{1}{2m}\sum_{i=1}^{m} (h_\theta(x^{(i)})-y^{(i)})^2 \]

你可能会感到奇怪,为什么分母是 \(2m\) 而不是 \(m\) ? 这并没什么特别的理由,只是在求损失函数的偏微分时可以和指数2结合,使最后的结果更加简洁(往下看就知道了,理论上分母是 \(m, 2m, 20m\) 都可以)。

为了使预测值尽可能接近实际值,我们要让 \(J(\theta)\) 达到最小值,数学上可以证明如果 \(J(\theta)\) 存在极值,那么一定在\(\frac{\partial J(\theta)}{\partial \theta}=0\) 时得到,\(\frac{\partial J(\theta)}{\partial \theta}\) 其实就是 \(J(\theta)\) 函数的梯度,或者更通俗一点,斜率。在机器学习里,使用梯度下降算法更新参数 \(\theta\) 是很常见的:

\[ \theta_j = \theta_j-\alpha\frac{\partial J}{\partial \theta_j}=\theta_j-\frac{\alpha}{m}\sum_{i=1}^{m}(h_\theta(x^{(i)})-y^{(i)})x_j \]

上式的 \(\alpha\) 取正值,给定了下降的速度,你也可以理解为机器学习的速度,显然 \(\alpha\) 越大,梯度下降的越快,但这通常并不是一件好事,过大的下降速度往往会导致迭代算法直接略过最优点从而无法收敛。如果你不断用合理的 \(\alpha\) 上式进行迭代,损失函数的值最后会无限接近于极值点。

我们当然可以用一个循环遍历 \(\theta_{j=1:n}\) (本例中 \(n=2\) ),当 \(n\) 很大时,使用循环来遍历很花时间,实际应用都会用矩阵计算来提升效率:

\[ \theta=\theta-\frac{\alpha}{m}[\sum_{i=1}^{m}(X^T\theta-y)]X \]

上式中的 \(\sum\) 表示将 \(X^T\theta-y\) 的所有元素相加,这种一次同时更新所有参数的形式叫做批量梯度下降法(Batch gradient descent),不断重复直至 \(\theta\) 收敛。(建议自己推导一下「=」右边的矩阵计算出来的结果是什么)

我用很短的篇幅简单介绍了线性回归模型实现的过程,你可能注意到有很多「坑」还没填,比如损失函数是怎么来的?变量不只一个怎么办?还有其他方法可以求解线性回归的模型吗?那就请接着往下看吧。

多变量的线性回归

我们之前说的根据身高(\(X\))来预测体重(\(y\))是一个单变量的线性回归问题,我想用身高+性别来预测体重怎么办?其实非常简单,引入新的特征 \(x_2\) 表示性别(1表示男性,-1表示女性)即可,此时对应的「直线」是:

\[ h_\theta=\theta_0+\theta_1x_1+\theta_2x_2=X^T\theta \]

其中 \(\theta\)\(3 \times1\) 阶矢量

\[ \theta=\left[ \begin{matrix} \theta_0\\ \theta_1\\ \theta_2 \end{matrix}\right] \]

\(X^T\)\(m\times3\)):

\[ X^T=\left[\begin{matrix} x_0^1, x_1^1,x_2^1\\ x_0^2, x_1^2,x_2^2\\ ...\\ x_0^m, x_1^m,x_2^m\\ \end{matrix}\right] \]

将公式中的 \(\theta\)\(X\) 用新的矩阵代替即可,当变量继续增加时,我们只要相应地修改矩阵即可,非常的简单。

损失函数是怎么来的

我们先定义模型预测的误差 \(\epsilon^{(i)}=y^{(i)}-h_{\theta}(x^{(i)})\),对于误差 $^{(i)} $,假设其服从正态分布: \[ \epsilon^{(i)} \sim\mathcal{N}(0,\sigma^2) \] 根据正态分布的性质,我们知道 \(y^{(i)} \sim\mathcal{N}(h_{\theta}(x^{(i)}),\sigma^2)\),据此写出\(y^{(i)}\) 的概率密度函数: \[ p(y^{(i)}|x^{(i)};\theta)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y^{(i)}-h_{\theta}(x^{(i)})^2}{2\sigma^2}} \] 让我们来捋一捋,这里 \(p\) 描述的是以 \(\theta\) 为参数,在 \(x^{(i)}\) 已知的情况下,预测结果等于 \(y^{(i)}\) 的概率,注意只针对 \((x^{(i)}, y^{(i)})\) 这一个数据,我们关心的是对于所有的训练数据,预测结果全中的概率,考虑到训练数据之间相互独立,所有训练数据预测结果都对的概率为: \[ p=\prod_{i=1}^{m}p(y^{(i)}|x^{(i)};\theta) \] 接下来要做的就是让这个概率尽可能的大。熟悉的味道,已知结果 \(y^{(i)}=h_{\theta}(x^{(i)})\),参数 \(\theta\) 未定,此处需要进行最大似然估计,如果你看过我之前的文章似然与极大似然估计,下面这一步是无比顺溜的: \[ \begin{split} \mathcal{L}&=\prod_{i=1}^{m}p(y^{(i)}|x^{(i)};\theta)\\ &=\prod_{i=1}^{m}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y^{(i)}-h_{\theta}(x^{(i)})^2}{2\sigma^2}}\\ \end{split} \]\(\mathcal{L}\) 进行对数化: \[ \begin{split} \log{\mathcal{L}}&=\sum_{i=1}^{m}\log{\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y^{(i)}-h_{\theta}(x^{(i)})^2}{2\sigma^2}}}\\ &=\sum_{i=1}^{m}\left(\log{\frac{1}{\sqrt{2\pi}\sigma}+\log{e^{-\frac{(y^{(i)}-h_{\theta}(x^{(i)})^2}{2\sigma^2}}}}\right)\\ &=m\log{\frac{1}{\sqrt{2\pi}\sigma}}-\sum_{i=1}^{m}\frac{(y^{(i)}-h_{\theta}(x^{(i)})^2}{2\sigma^2}\\ &=const-\frac{1}{\sigma^2}\cdot\frac{1}{2}\sum_{i=1}^{m}(y^{(i)}-\theta^Tx^{(i)})^2 \end{split} \] 最大化似然函数 \(\mathcal{L}\) 等价于最小化 \(\frac{1}{2}\sum_{i=1}^{m}(y^{(i)}-\theta^Tx^{(i)})^2\) ,也等价于最小化 \(\frac{1}{2m}\sum_{i=1}^{m}(y^{(i)}-\theta^Tx^{(i)})^2\),正是我们在上一节所说的损失函数。

通过这一部分的推导,我们发现所谓的梯度下降使损失函数最小化的过程实质上是最大化似然函数的过程。

为什么梯度下降法能减小损失函数

考虑 \(\theta_j\) 的迭代过程(变量只有一个),损失函数是一个存在极小值点的二次函数,用初中生都能听懂的话来说就是一条开口向上的抛物线:

  1. 当起点在右侧时,\(\frac{\partial J}{\partial \theta_j}>0\), 要损失函数降低就必须让 \(\theta_j\)减小以接近极值点,此时\(\theta_j^{new} = \theta_j-\alpha\frac{\partial J}{\partial \theta_j}<\theta_j\)
  2. 当起点在左侧时,\(\frac{\partial J}{\partial \theta_j}<0\), 要损失函数降低就必须让 \(\theta_j\)增大以接近极值点,此时\(\theta_j^{new} = \theta_j-\alpha\frac{\partial J}{\partial \theta_j}>\theta_j\)

综上,我们发现在梯度下降算法里,梯度的正负可以告知我们现在的点位于损失函数的哪一边,为我们确定下一步迭代的方向,增大还是减小;而迭代的幅度其实是由 \(\alpha\) 控制的,\(\alpha\) 的值越小,算法收敛的可能性越高,而且结果越接近理论值,但是迭代的幅度相应也越小,高精度是以更多的迭代次数为代价的

正规方程 Normal Equation

除了用梯度下降法求解极值点以外,对于简单的现象问题还可以直接用线性方程求解,它的好处是无需迭代、不用担心 \(\alpha\) 的大小问题,让我们回到损失函数先: \[ J(\theta)=\frac{1}{2m}\sum_{i=1}^{m} (h_\theta(x^{(i)})-y^{(i)})^2 \]

将上式写成矩阵的形式: ​ \[ \begin{split} J&=\frac{1}{2m}\sum_{i=1}^{m} (h_\theta(x^{(i)})-y^{(i)}) \cdot(h_\theta(x^{(i)})-y^{(i)})\\ &=\frac{1}{2m}(X\theta-y)^T(X\theta-y) \end{split} \] \(J\) 最后的结果是一个数,但也可以理解为\(1\times 1\) 的矩阵,我们知道要求得损失函数的极小值需要令 \(J\) 的梯度为0,由于 \(m\) 并不影响 \(J\) 的极值点,我们对 \(\frac{1}{2}(X\theta-y)^T(X\theta-y)\) 求极值即可: \[ \begin{split} \nabla_\theta J&=\frac{1}{2}\nabla_\theta(X\theta-y)^T(X\theta-y)\\ &=\frac{1}{2}\nabla_\theta(\theta^TX^T-y^T)(X\theta-y)\\ &=\frac{1}{2}\nabla_\theta(\theta^TX^TX\theta-y^TX\theta-\theta^TX^Ty+y^Ty)\\ &=\frac{1}{2}\nabla_\theta(\theta^TX^TX\theta-y^TX\theta-\theta^TX^Ty)\\ &=\frac{1}{2}\nabla_\theta \text{tr}(\theta^TX^TX\theta-y^TX\theta-\theta^TX^Ty)\\ &=\frac{1}{2}\nabla_\theta (\text{tr} \theta^TX^TX\theta-2\text{tr} y^TX\theta)\\ &=\frac{1}{2}(2X^TX\theta-2X^Ty)\\ &=X^TX\theta-X^Ty\\ \end{split} \] ps. 上式的推导过程需要用到矩阵的迹(trace)相关性质以及矩阵函数求导的知识,限于本文的篇幅这里不再介绍,有需要的话可以参阅 Andrew Ng 的讲义

\(\nabla_\theta J=0\),可以得到下式: \[ \begin{split} &X^TX\theta-X^Ty=0\\ \Leftrightarrow &X^TX\theta=X^Ty\\ \Leftrightarrow&\theta=(X^TX)^{-1}X^Ty\\ \end{split} \] 正规方程求得的是解析解,无需迭代,也不需要进行 feature scaling,看上去简直完美,相比之下梯度下降算法毫无可取之处?理想丰满抵不过现实残酷,求解正规方程需要求矩阵的逆,学过线性代数的同学都知道不是所有的矩阵都能求逆,而且求逆也是一个很花时间的过程(\((X^TX)^{-1}\) 的求解复杂度为 \(o(n^3)\)),在矩阵很大时,求逆的计算量往往大于梯度下降算法,因此实际问题还是以梯度下降这种迭代方法为主,当研究对象的特征较小时可以考虑使用正规方程求解。