20 KiB
线性回归问题总结
目录
- 摘要
- 线性回归概述
- 假设函数与平方损失函数
- 为什么是最小二乘法
- 多项式回归
- 梯度下降法
- 一般描述
- 学习率的选择
- 在多变量下存在的问题
- 特征缩放
- 均值归一化
- 局部最优
- 平方损失函数是凸函数
- 规范方程法
- 直觉理解
- 严格数学推导
- 不可逆
- 过拟合
- 过拟合概述与产生原因，解决方案
- 正则化的损失函数
- 规范方程法：一定是可逆阵
- 梯度下降法：仍然是凸函数
摘要
这篇文章主要介绍线性回归问题。任何一种统计学习方法都由三个要素组成，即模型，策略以及算法。所谓模型，就是本文中的假设函数(hypothesis function)；策略或者评价准则，即本文中的损失函数；而算法，则是选取最优模型的方法。本文将从这三个方面对线性回归问题进行总结。
第一部分阐述了线性回归问题的基本概念，给出了它的假设函数以及损失函数的数学表达式，并且指出平方损失函数在本质上与最大似然估计无异。接下来的两个部分，分别阐述了求解该模型的两种算法，即梯度下降法与规范方程法，给出了两种方法的数学推导与数学表达式，并就实际应用中可能存在的问题进行讨论。最后，介绍了在统计学习中常见的一个问题——过拟合，并且就解决过拟合问题的常见方法——正则化(regularization)进行了深入的探讨。
这篇文章的重点在于：平方损失函数与最大似然估计的联系，平方损失函数是凸函数的证明，规范方程法的数学推导，正则化后的平方损失函数仍然是凸函数的证明以及正则化的规范方程一定可逆的证明。
线性回归概述
假设函数与平方损失函数。
线性回归问题是这样一个问题，对于输入$x$（$x$是一个$n$维向量，表示$n$个特征），假设输出$y$是$x$的线性函数。不妨设假设函数为
\hat{y} = h_\theta(x) = \theta_0 + \theta_1x_1 + \cdots + \theta_nx_n
为了简单起见，不妨引入一个$x_0 = 1$，这样
\hat{y} = \Sigma_{i = 0}^n \theta_ix_i = \theta^Tx
这里$\hat{y}$表示对$y$的估计值。因此，线性回归的任务，就是找到一个合适的列向量$\theta$，使得假设函数尽可能反映输入样本的特征。为了衡量假设函数与输入样本的符合程度，需要定义一个新的函数，这里称为损失函数或者代价函数(cost function)。
在线性回归问题中，常用的损失函数为平方损失函数(square cost function)，定义为J(\theta)
J(\theta) = \frac{1}{2m}\Sigma_{i = 1}^m (\hat{y}^{(i)} - y^{(i)})^2 = \frac{1}{2m}\Sigma_{i = 1}^m (h_\theta(x^{(i)}) - y^{(i)})^2
容易注意到，平方损失函数的主体是用最小二乘法给出的，它反映了样本预测值与实际值之间的偏离程度。假设函数越符合样本，则损失函数的值就会越小；反之，假设函数越偏离样本，损失函数的值就会越大。因此可以用损失函数的取值来定量衡量参数$\theta$的好坏，最优的一组$\theta$值，应该恰好使损失函数取到最小值。这样，就将参数优化问题，转化为了求解最值问题。
为什么使用平方损失函数？
使用平方损失函数并非偶然，而是由最大似然估计推导而出的。下面给出推导过程。
在线性回归问题中，模型假设实际上是输出$y$乃以$h_\theta(x)$为均值的正态分布，即
y|x = h_\theta(x) + \epsilon,\ \epsilon \sim N(0, \sigma^2)
所以有
y|x \sim N(h_\theta(x), \sigma^2)\\\\
f(y | x) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y - \theta^Tx)^2}{2\sigma^2}}
因此可以写出似然函数L(\theta)
L(\theta) = \Pi_{i = 1}^m f(y^{(i)}|x^{(i)}) = (\frac{1}{\sqrt{2\pi}\sigma})^m\cdot e^{-\frac{1}{2\sigma^2}\Sigma_{i = 1}^m (y^{(i)} - \theta^Tx^{(i)})^2}\\\\
lnL(\theta) = -mln(\sqrt{2\pi}\sigma) - \frac{1}{2\sigma^2}\Sigma_{i = 1}^m(y^{(i)} - \theta^Tx^{(i)})^2
这样，对似然函数求最大值，就等价于求
\Sigma_{i = 1}^m(y^{(i)} - \theta^Tx^{(i)})^2
的最小值，此时上面的平方损失函数也取到最小值。因此可以看到，平方损失函数实际上与最大似然估计是一脉相承的。
多项式回归。
实际上，线性回归还可以应用在非线性的场合，比方说有一个非线性的假设函数
h_\theta(x) = \theta_0 + \theta_1x + \theta_2x^2
不妨设$x_1 = x, x_2 = x^2$，这样就可以将上面的假设函数写作
h_\theta(x_1, x_2) = \theta_0 + \theta_1x_1 + \theta_2x_2
这样，就把非线性的假设函数转化为了线性回归问题，从而可以利用线性回归的方法来解决。该问题被称为多项式回归问题(polynomial regression)。
理论上来讲，一般的可导函数总是可以通过泰勒展开，写成多项式的形式，因此多项式回归理论上可以解决很多复杂的问题。然而，在实际应用中，当选择的多项式形式过于复杂时，可能会产生过拟合(overfitting)的问题，将会在后面阐述。
梯度下降法
梯度下降法的基本表述。
梯度下降法(gradient descent)是求解最优化问题的常见数值方法。它的基本思想是在每次迭代过程中，使得被优化的参数沿着函数下降速度最快的方向（即梯度的方向）移动一个微小的步长，随着迭代的进行，被优化函数的取值会越来越小，直至收敛至极小值点。
对于任意被优化函数$J(\theta)$，梯度下降法的参数更新公式为
\theta_j := \theta_j - \alpha\frac{\partial}{\partial \theta_j}J(\theta)
这里的$\alpha$称为学习率(learning rate)，它其实就是梯度下降法采用的步长。采用梯度下降法来求解平方损失函数的最小值，则有
\theta_j := \theta_j - \frac{\alpha}{m}\Sigma_{i = 1}^m(h_\theta(x^{(i)}) - y^{(i)})x_j^{(i)}
通过选择合适的迭代次数和学习率，就可以得到使$J(\theta)$最小的$\theta$了。
梯度下降法存在的问题若干
学习率$\alpha$的选择。
学习率的选择在梯度下降法中至关重要。如果学习率太小，则需要更多次迭代才能找到最优解，需要较长的学习时间；而如果学习率太大，也有可能导致收敛速度慢，甚至可能会发散。这是因为学习率太大时，不再满足全微分方程中微小变化量的条件，因此每次迭代后被优化函数都未必可以取到比上一次迭代更小的值。学习率的选择对$J(\theta)$收敛速度的影响如下图所示：
可以看到，当$\alpha$很小时（如图中红色曲线所示），损失函数$J(\theta)$收敛速度较慢。随着$\alpha$逐渐增大，曲线逐渐变得陡峭，收敛速度越来越快。可是当$\alpha$很大时（$\alpha = 1.3$，图中黑色曲线）收敛速度反而相对此前变慢了；$\alpha$继续增大将导致代价函数发散（图中洋红色曲线）。
因此，理想的方案是选择这样一个$\alpha$，它应该处于收敛与发散的临界点，一方面保证被优化函数最终收敛，另一方面还要有较快的收敛速度。
局部最优解。
从上面的讨论中可以看出，梯度下降法最终将收敛到某个局部最优点，而无法保证收敛到全局最优点。实际上，当选择不同的初始值时，使用梯度下降法往往会收敛到不同的极值点，如下图所示：
因此，使用梯度下降法无法保证得到的结果一定是全局最优解。然而，需要指出，对于线性回归与平方损失函数，则不存在这个问题，因为平方损失函数是凸函数。证明如下：
平方损失函数$J(\theta)$满足
J(\theta) = \frac{1}{2m}\Sigma_{i = 1}^m (h_\theta(x^{(i)}) - y^{(i)})^2
所以
\frac{\partial}{\partial \theta_j}J(\theta) = \frac{1}{m}\Sigma_{i = 1}^m(h_\theta(x^{(i)}) - y^{(i)})x_j^{(i)}\\\\
\frac{\partial^2}{\partial \theta_j\partial \theta_k}J(\theta) = \frac{1}{m}\Sigma_{i = 1}^mx_j^{(i)}x_k^{(i)}
因此，可以写出Hessian矩阵H:
H = \frac{1}{m}\left[
\begin{matrix}
\Sigma_{i = 1}^mx_0^{(i)^2} & \Sigma_{i = 1}^mx_0^{(i)}x_1^{(i)} & \cdots & \Sigma_{i = 1}^mx_0^{(i)}x_n^{(i)}\\\\
\Sigma_{i = 1}^mx_1^{(i)}x_0^{(i)} & \Sigma_{i = 1}^mx_1^{(i)^2} & \cdots & \Sigma_{i = 1}^mx_1^{(i)}x_n^{(i)}\\\\
\vdots & \vdots & & \vdots&\\\\
\Sigma_{i = 1}^mx_n^{(i)}x_0^{(i)} & \Sigma_{i = 1}^mx_n^{(i)}x_1^{(i)} & \cdots & \Sigma_{i = 1}^mx_n^{(i)^2}\\\\
\end{matrix}
\right]
容易看出，
H = \frac{1}{m}X^TX
因此$H$显然是一个半正定矩阵，根据凸函数判别的充分条件，可以得出平方损失函数为凸函数。
凸函数具有一系列的优良性质，比如说只有一个极值点，即全局最小值点。因此对凸函数使用梯度下降法总是可以收敛至全局最优点。
多特征的线性回归
在实际应用中，一个问题往往会有成百上千个特征，甚至更多。这些不同特征的有截然不同的取值范围，如果不同特征取值范围差别很大，则会导致梯度下降法收敛速度慢的问题。
考虑存在两个特征$x_1, x_2$，它们各自的取值范围分别是$0 < x_1 < 1$，$0 < x_2 < 1000$，两者相差了1000倍。在使用梯度下降法时，参数$\theta_2$的变化量也将大约是$\theta_1$的1000倍，这将导致损失函数等高线图呈现扁平的椭圆状，梯度下降路径相对取值范围更大的特征对应的参数来回波动，从而导致收敛速度缓慢。
为了解决上面的问题，需要对取值相差比较大的特征进行特征缩放(feature scaling)与均值规范(mean normalization)，具体的操作是求出每一个特征的均值$\mu$以及范围（最大值减去最小值），随后对每一个特征点进行以下更新
x_j^{(i)} := \frac{x_j^{(i)} - \mu^{(i)}}{range^{(i)}}
在实际应用中，也可以用标准差(standard deviation)来代替范围进行更新。这样更新后的特征取值范围都在$[-1, 1]$之间，绘出的等高线图接近于（高维）球形，收敛速度就要快多了。
规范方程
规范方程的直觉理解
前面提到假设函数的表达式为
h_\theta(x) = \theta^Tx = x^T\theta
设样本容量为$m$，则输入$X$以及对应的输出$Y$都可以写成矩阵的形式
X = \left[
\begin{matrix}
x^{(1)^T}\\\\
x^{(2)^T}\\\\
\vdots\\\\
x^{(m)^T}
\end{matrix}
\right],\ \
Y = \left[
\begin{matrix}
y^{(1)}\\\\
y^{(2)}\\\\
\vdots\\\\
y^{(m)}
\end{matrix}
\right],\ \
\theta = \left[
\begin{matrix}
\theta_0\\\\
\theta_1\\\\
\vdots\\\\
\theta_n
\end{matrix}
\right]
这样，就可以把假设函数也写成矩阵形式
Y = X\theta
为了得到$\theta$，只需求解该矩阵方程，为此将等式的两边同时乘以X^T
X^TY = X^TX\theta
所以
\theta = (X^TX)^{-1}X^TY
这就是规范方程法的数学表述。以下就这个结论给出严格推导。
规范方程法的推导
平方损失函数(cost function)$J(\theta)$可以写作
J(\theta) = \frac{1}{2m}\Sigma_{i = 1}^m[h_\theta(x^{(i)}) - y^{(i)}]^2
参数$\theta$应该是使$J(\theta)$最小的一组列向量。由多元函数极值可知，$\theta$应该满足：
\frac{\partial}{\partial \theta_j}J(\theta) = \frac{1}{m}\Sigma_{i = 1}^m[h_\theta(x^{(i)}) - y^{(i)}]x_j^{(i)} = 0
解上述方程即可以得到$\theta$的解析解。
为了简单起见，下面用矩阵形式来对$J(\theta)$进行表述。所以$J(\theta)$可以写作
J(\theta) = \frac{1}{2m}(X\theta - Y)^T(X\theta - Y) = \frac{1}{2m}(\theta^TX^TX\theta - 2Y^TX\theta + Y^TY)
对$\theta$求导，所以
\frac{d}{d\theta}J(\theta) = \frac{1}{2m}(2X^TX\theta - 2X^TY) = 0
从而可得
\theta = (X^TX)^{-1}X^TY
该结论与上面的直觉理解完全一致。
规范方程法存在的问题
需要注意的是，在实际应用中，$X^TX$往往是不可逆的，比如当特征的数量多于样本容量时。因此，往往是通过对$X^TX$求伪逆，来代替$(X^TX)^{-1}$。
梯度下降法与规范方程法的比较
在通常情况下，当样本的特征比较少时，比如特征数少于10000时，选择规范方程法会有比较好的效果；而当样本数量太多，则应该选择梯度下降法来对模型进行求解。
这是因为，在规范方程法中，涉及到矩阵的求逆运算，其时间复杂度是$O(n^3)$，因此当模型的特征太多时，该逆矩阵的计算会消耗大量的时间。反之，即使特征数量非常多，梯度下降法也仍然具有较好的性能。
过拟合
过拟合概述
过拟合是统计学习中非常常见的一个问题，它是指模型的选择过于复杂，因而可以在训练集上得到很好的训练效果，使得损失函数接近于零，可是却难以将该模型泛化(generalize)到其他应用场合，比如说在一个新的测试集上预测准确率就会非常低下。欠拟合，正常拟合，以及过拟合的关系如下图所示
可以看到，第一张图中的模型选择过于简单，用线性函数去拟合一个多项式回归问题，发生了欠拟合；第三张图的模型则过于复杂，用了一个非常高阶的多项式来拟合一个低阶问题，发生了过拟合。
关于过拟合，我想到这样一个直观的理解：训练模型的过程就好比一个人复习考试的过程，一种复习方案是深入理解题目背后的原理和知识点，另一种复习方案则是把每个题目都记忆下来。显然后者在训练集上会有出色的表现，如果记忆力足够好，甚至可以达到100%的正确率，可是在实际考试过程中，遇到的都是他从未记忆过的题目，此时就束手无策了。第二种复习方案就是这里的过拟合。
过拟合会发生，主要是因为模型选择过于复杂，选择的特征太多。因此，一种解决过拟合的方案，就是手动剔除一些不重要的特征，简化模型，或者可以使用模型选择算法(model selection algorithm)，后者将在后面的文章中提到。另一种更加常见的方案，就是下面要重点阐述的正则化。
正则化的损失函数
考虑这样一种情况，当前模型的假设函数为
h_\theta(x) = \theta_0 + \theta_1x + \theta_2x^2 + \theta_3x^3 + \theta_4x^4
为了避免发生过拟合，我可以限制$\theta_3, \theta_4$的大小，使得高阶项具有相对更小的权重。具体的方案是在损失函数中对$\theta_3, \theta_4$添加罚项(penalty)，比如
J(\theta) = \frac{1}{2m}\Sigma_{i = 1}^m[h_\theta(x^{(i)}) - y^{(i)}]^2 + 1000\theta_3^2 + 1000\theta_4^2
这样，为了使$J(\theta)$取到最小值，$\theta_3, \theta_4$的取值必然相对较小，因此在一定程度上简化了模型。
在实际应用中，很难知道哪一些特征会导致模型过拟合，因此常见的方案是对每一个$\theta_j$，都添加一个上述的罚项，也称为正则化项(regularized term)。这样，正则化的损失函数就有以下的形式：
J(\theta) = \frac{1}{2m}[\Sigma_{i = 1}^m(h_\theta(x^{(i)}) - y^{(i)})^2 + \lambda\Sigma_{i = 1}^n\theta_i^2]
需要注意的是，这里并没有对$\theta_0$引入正则化项，老师说这是一种约定俗成，并没有特别的原因。使用梯度下降法或者规范方程法，就可以求得在正则化下的最优参数$\theta$了。
引入正则化项的规范方程推导
为了令引入正则化项的$J(\theta)$取到最小值，同样可以利用多元函数极值，令
\frac{\partial}{\partial \theta_j}J(\theta) = 0
逐个求解出$\theta_j$，就可以得到最优的一组$\theta$，在此不再赘述。
为了简化问题，仍然用矩阵形式表示$J(\theta)$，所以
\begin{aligned}
J(\theta) = &\frac{1}{2m}[(X\theta - Y)^T(X\theta - Y) + \lambda(I\theta)^TI\theta]\\\\
= &\frac{1}{2m}(\theta^TX^TX\theta - 2Y^TX\theta + Y^TY + \lambda(I\theta)^TI\theta]
\end{aligned}
其中
I = \left[
\begin{matrix}
0 & 0 & \cdots & 0\\\\
0 & 1 & \cdots & 0\\\\
\vdots & \vdots & & \vdots\\\\
0 & 0 & \cdots & 1
\end{matrix}
\right].
所以
\frac{\partial}{\partial \theta}J(\theta) = \frac{1}{2m}(2X^TX\theta - 2X^TY + 2\lambda I\theta) = 0
从而可以解出
\theta = (X^TX + \lambda I)^{-1}X^TY = (X^TX + \lambda \left[
\begin{matrix}
0 & 0 & \cdots & 0\\\\
0 & 1 & \cdots & 0\\\\
\vdots & \vdots & & \vdots\\\\
0 & 0 & \cdots & 1
\end{matrix}
\right])^{-1}X^TY
需要说明的是，正则化项的引入不仅解决了过拟合的问题，还使得新的矩阵
X^TX + \lambda I
一定是可逆矩阵，而此前的
X^TX
则不能保证它的可逆性。下面给出证明。
前面已经指出
X = \left[
\begin{matrix}
x^{(1)^T}\\\\
x^{(2)^T}\\\\
\vdots\\\\
x^{(m)^T}
\end{matrix}
\right] = \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 & & \vdots\\\\
x_0^{(m)} & x_1^{(m)} & \cdots & x_n^{(m)}\\\\
\end{matrix}
\right]
所以
X^T = \left[
\begin{matrix}
x_0^{(1)} & x_0^{(2)} & \cdots & x_0^{(m)}\\\\
x_1^{(1)} & x_1^{(2)} & \cdots & x_1^{(m)}\\\\
\vdots & \vdots & & \vdots\\\\
x_n^{(1)} & x_n^{(2)} & \cdots & x_n^{(m)}\\\\
\end{matrix}
\right]
容易注意到，$X^TX$的对角线元素一定大于或者等于零。且对角线上第一个元素为
a_{11} = \Sigma_{i = 1}^mx_0^{(i)}\cdot x_0^{(i)} = m
因此$X^T + \lambda I$的对角线元素$a_{ii}$一定大于$0$，对于任意的$\lambda > 0$成立。所以
\left| X^T + \lambda I \right| = \Pi_{i = 1}^{n + 1}a_{ii} > 0
因此$X^T + \lambda I$一定是一个可逆矩阵。实际上，还可以得到更进一步的结论，即$X^T + \lambda I$一定是一个正定矩阵。这是因为$X^TX$和$I$都是半正定矩阵，即对于任意的非零向量$\eta$，都有
\eta^TX^TX\eta \ge 0,\ \eta^TI\eta \ge 0
因此两者之和也必定是一个半正定矩阵。由于此前已经得到$X^T + \lambda I$可逆，故不存在为零的特征值，所以所有特征值都必定大于零，因此$X^T + \lambda I$是一个正定矩阵。
引入正则化项的梯度下降法
关于梯度下降法自然不必多说，这里主要证明，引入正则化项后的平方损失函数仍然是凸函数，因此使用梯度下降法仍然可以收敛到全局最优点。
再次给出$J(\theta)$的表达式
J(\theta) = \frac{1}{2m}[\Sigma_{i = 1}^m(h_\theta(x^{(i)}) - y^{(i)})^2 + \lambda\Sigma_{i = 1}^n\theta_i^2]
所以有
\frac{\partial}{\partial \theta_0}J(\theta) = \frac{1}{m}\Sigma_{i = 1}^m(h_\theta(x^{(i)}) - y^{(i)})x_0^{(i)}\\\\
\frac{\partial}{\partial \theta_j}J(\theta) = \frac{1}{m}\Sigma_{i = 1}^m(h_\theta(x^{(i)}) - y^{(i)})x_j^{(i)} + \frac{1}{m}\theta_j\ (j \ne 0)
\frac{\partial^2}{\partial \theta_0^2}J(\theta) = \frac{1}{m}\Sigma_{i = 1}^mx_0^{(i)}x_0^{(i)}\\\\
\frac{\partial^2}{\partial \theta_j^2}J(\theta) = \frac{1}{m}\Sigma_{i = 1}^mx_j^{(i)}x_j^{(i)} + \frac{\lambda}{m}\ (j \ne 0)\\\\
\frac{\partial^2}{\partial \theta_j\partial \theta_k}J(\theta) = \frac{1}{m}\Sigma_{i = 1}^mx_j^{(i)}x_k^{(i)} \ (j \ne k)
因此，可以写出Hessian矩阵
H = \frac{1}{m}\left[
\begin{matrix}
\Sigma_{i = 1}^mx_0^{(i)^2} & \Sigma_{i = 1}^mx_0^{(i)}x_1^{(i)} & \cdots & \Sigma_{i = 1}^mx_0^{(i)}x_n^{(i)}\\\\
\Sigma_{i = 1}^mx_1^{(i)}x_0^{(i)} & \Sigma_{i = 1}^mx_1^{(i)^2} + \lambda & \cdots & \Sigma_{i = 1}^mx_1^{(i)}x_n^{(i)} \\\\
\vdots & \vdots & & \vdots&\\\\
\Sigma_{i = 1}^mx_n^{(i)}x_0^{(i)} & \Sigma_{i = 1}^mx_n^{(i)}x_1^{(i)} & \cdots & \Sigma_{i = 1}^mx_n^{(i)^2} + \lambda\\\\
\end{matrix}
\right]
容易看出，
H = \frac{1}{m}(X^TX + \lambda I)
上面已经证明了，该矩阵是一个正定矩阵。根据凸函数的充分条件，可以得到正则化后的平方损失函数仍然是一个凸函数。