线性回归的模型描述

在监督学习模型里，我们通常有一个已知的数据集，称为训练集(Training Set)。而监督学习要做的就是让计算机根据训练集数据实现输出。

在继续进行之前，我们先做如下记号规定

$m$ ：训练样本的数量
$\boldsymbol{x}$ ：输入变量
$\boldsymbol{y}$ ：输出变量
$(x,y)$ ：一个训练样本
$(x^{(i)},y^{(i)})$ ：训练集中索引为 $i$ 的训练样本
$h$ ：模型使用的假设函数(Hypothesis Function)，构造从 $\boldsymbol{x}$ 到 $\boldsymbol{y}$ 的映射

假设函数可以有多种形式，当 $h$ 为线性函数形式时，亦即

$h(x)=\theta_0+\theta_1x$

此时对应的模型称为线性回归模型(Linear Regression)，其对应的假设函数也可记为 $h_{\theta}(x)$ 。在本小节中，我们讨论的都是只包含一个输入变量的线性回归问题，因此也称为单变量线性回归模型。

[{"url":"/img/mlearning/m2f1.png","alt":"线性回归模型","title":""}]

代价函数

代价函数的数学定义

在线性回归模型的假设函数中，参数 $\theta_0,\theta_1$ 是可以任意选取的，因此如何有效选取这些参数是个重要的问题。

为了选取适当的 $\theta_0,\theta_1$ ，我们需要让 $h(x)$ 对于训练集中的的每一个 $x^{(i)}$ ，输出一个值 $y$ 且使得 $y$ 与 $y^{(i)}$ 的差距尽可能小。这个问题也可以被叙述成最小化问题，我们需要找到一对值 $(\theta_0,\theta_1)$ ，使得各个训练数据的输出误差都尽可能小。用数学语言表述如下

$\text{Find a pair of value }(\theta_0,\theta_1)\quad s.t.\min \frac{1}{2m}\sum_{i=1}^{m}\Big(h_{\theta}(x^{(i)})-y^{(i)}\Big)^2$

在上式子中， $\theta_0,\theta_1$ 是影响后一部分函数取值的变量，因此我们定义代价函数(Cost Function)

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

代价函数也称作平方误差函数，因为求和部分是误差项 $h_{\theta}(x^{(i)})-y^{(i)}$ 的平方，使得求和内容均非负。平方误差函数解决大多数线性回归问题的较优选择。故问题转化为

$\text{Find a pair of value }(\theta_0,\theta_1)\quad s.t.\min J(\theta_0,\theta_1)$

代价函数的讨论

我们先考虑 $\theta_0=0$ 的情况，也就是只保留线性部分。此时代价函数是关于 $\theta_1$ 的函数

$\begin{aligned}J(\theta_1)&=\frac{1}{2m}\sum_{i=1}^{m}\Big(h_{\theta}(x^{(i)})-y^{(i)}\Big)^2\\ &=\frac{1}{2m}\sum_{i=1}^{m}\Big(\theta_1x^{(i)}-y^{(i)}\Big)^2\end{aligned}$

通过描点绘图等方法不难看出，代价函数 $J$ 是关于 $\theta_1$ 的二次函数，且二次项系数为正。由二次函数的图像特点可以知道，此时一定存在一个 $\theta_1$ 使得函数取得最小值，这一点是比较重要的。

[{"url":"/img/mlearning/m2f2.png","alt":"简化版代价函数","title":""}]

在简化版代价函数的基础上，引入 $\theta_0$ 导致代价函数 $J$ 变成了一个多元函数。为了直观地解决最小化问题，我们可以通过绘制等值线图如下。

[{"url":"/img/mlearning/m2f3.png","alt":"代价函数的等值线图","title":""}]

在每一条等值线上，不同的 $\theta_0,\theta_1$ 对应的代价函数值是一样的。如果能够借助计算机工具绘制出这些等值线图，我们便可以直观找出满足条件的参数 $\theta_0,\theta_1$ 。然而，可视化通常是困难的。

梯度下降

梯度下降的数学定义

现在来介绍一种数学方法，能够实现代价函数的最小化。我们不妨先看看更一般的情况

代价函数 $J$ 是关于 $n$ 个参数 $\theta_0,\theta_1,\cdots,\theta_{n-1}$ 的函数
需要找出一对值 $(\theta_0,\theta_1,\cdots,\theta_{n-1})$ ，使得 $J$ 最小。

梯度下降(Gradient Descent)法的流程如下

任意选取一对值 $(\theta_0,\theta_1,\cdots,\theta_{n-1})$
通过某种变换改变 $(\theta_0,\theta_1,\cdots,\theta_{n-1})$ 的值
重复上述过程直到我们找到 $J$ 的最小值或局部最小值

现在回到线性回归模型上，我们规定梯度下降法中的某种变换如下

$\theta_{j}:=\theta_j-\alpha\frac{\partial}{\partial \theta_j}J(\theta_0,\theta_1),\quad(j=0,1)$

在上式中，符号:=表示右式赋值给左式， $\partial$ 是偏导符号。我们要求 $\theta_0$ 和 $\theta_1$ 是一起变换的，因为 $J$ 是同时关于 $\theta_0,\theta_1$ 的函数。

该算法由梯度(Gradient)的数学性质保证。一个多元函数 $f$ 在某点 $P$ 的梯度 $\nabla f|_{P}$ 对应着其增长率最大的方向，其负方向 $-\nabla f|_{P}$ 对应着函数降低最快的方向。

学习率

在变换中， $\alpha$ 表示学习率(Learning Rate)，其不同取值对算法有不同影响。

如果 $\alpha$ 很小，那么梯度下降的过程会很“缓慢”
如果 $\alpha$ 很大，那么梯度下降的过程会很“迅速”，但可能不收敛

[{"url":"/img/mlearning/m2f4.png","alt":"学习率对梯度下降的影响","title":""}]

如何选取学习率 $\alpha$ 的取值也是一个问题，本节不涉及。

梯度下降法的局限性

梯度下降法只能保证找到局部最小值，不能保证找到最小值。当前状态已经处于局部最小值时，根据多元函数的性质，函数在该点对各个变元的偏导为0，因此算法只能收敛于此局部最小点。
固定 $\alpha$ 时，梯度下降法的迭代过程会越来越缓慢。此时可以适当增大学习率 $\alpha$ 。

线性回归的梯度下降法

数学模型

为了编写线性回归的梯度下降算法，我们必须要搞清楚代价函数 $J$ 对 $\theta_0,\theta_1$ 的偏导是什么。为此，我们需要进行适当的数学推导。展开求和项后求偏导，并代入算法，结果如下

$\begin{cases}\theta_0:=\theta_0-\alpha\cdot\dfrac{1}{m}\displaystyle\sum_{i=1}^{m}\Big(h_{\theta}(x^{(i)})-y^{(i)}\Big)\\\\ \theta_1:=\theta_1-\alpha\cdot\dfrac{1}{m}\displaystyle\sum_{i=1}^{m}\Big(h_{\theta}(x^{(i)})-y^{(i)}\Big)x^{(i)}\end{cases}$

线性回归模型的代价函数总是一个凸函数(Convex Function)，直观上看是一个开口向上的碗型形状，且只有一个全局最小值点。这也就是说，线性回归模型是可以利用梯度下降法找到全局最优解的。

[{"url":"/img/mlearning/m2f5.png","alt":"线性回归模型的代价函数","title":""}]

梯度下降算法是一个迭代算法。当其每一次迭代都遍历所有的训练集数据时，称该梯度算法为Batch梯度下降；除了梯度下降法，还有正规方程组法等可以求解代价函数的最小值。

代码实现

本代码基于Jupyter Notebook书写，点此处查看相关下载、使用说明等。语言使用Python。

import numpy as np 
import matplotlib.pyplot as plt 

# sample length
m = 18

# generate matrix X
X0 = np.ones((m, 1))
X1 = np.arange(1, m+1).reshape(m, 1)
X = np.hstack((X0, X1))

# matrix y
y = np.array([2,3,3,5,8,10,10,13,15,15,16,19,19,20,22,22,25,28]).reshape(m,1)

# learning rate alpha
alpha = 0.01

def cost_function(theta, X, y):
    diff = np.dot(X, theta) - y
    return (1./2*m) * np.dot(np.transpose(diff), diff)

def gradient_function(theta, X, y):
    diff = np.dot(X, theta) - y
    return (1./m) * np.dot(np.transpose(X), diff)

def gradient_descent(X, y, alpha):
    theta = np.array([1, 1]).reshape(2, 1)
    gradient = gradient_function(theta, X, y)
    while not np.all(np.absolute(gradient) <= 1e-5):
        theta = theta - alpha * gradient
        gradient = gradient_function(theta, X, y)
    return theta

[theta0, theta1] = gradient_descent(X, y, alpha)
plt.figure()
plt.scatter(X1,y)
plt.plot(X1, theta0 + theta1*X1, color='r')
plt.title('Linear Regression with Gradient Descent')
plt.grid(True)
plt.show()