12. 利用梯度下降法解线性回归

来源：原始 Jupyter Notebook。内容已转换为网页阅读格式；下载原文件。

一元线性回归的梯度下降

下面使用梯度下降算法来求解前面的线性回归问题，我们知道，对一元线性回归，损失函数为：

loss = J(a, b) = \frac{1}{2m} \sum_{i=1}^{m}(y_i - (ax_i + b))^2

这是一个二元二次方程，我们对该方程求梯度，也就是分别对参数 a 和 b 求偏导，得到参数的更新公式：

\begin{align} a &:= a - \eta \frac{\partial J(a, b)}{\partial a} \\ b &:= b - \eta \frac{\partial J(a, b)}{\partial b} \\ \end{align}

将损失函数带入上式求导，得到：

\begin{align} a &:= a - \eta \frac{1}{m} \sum_{i=1}^{m}(y_i - (ax_i + b))x_i \\ b &:= b - \eta \frac{1}{m} \sum_{i=1}^{m}(y_i - (ax_i + b)) \end{align}

我们使用代码来实现这个过程，先来看一个最简单的例子，手工造几个散点，然后使用梯度下降法来求解线性回归。如下所示：

import matplotlib.pyplot as plt
import numpy as np

X = np.array([0,1,2])
Y = np.array([0.9, 3.1, 5.1])

plt.plot(X, Y, 'k.')
plt.show()

运行结果

Notebook 运行结果

很显然这是一个一元线性回归，根据上面的介绍，我们可以定义出损失函数和参数 a, b 的更新公式：

# 损失函数
def loss(X, Y, a, b):
    s = 0
    for i in range(X.size):
        s += ((a*X[i] + b)-Y[i])**2
    return s * (1./2*Y.size)

# 参数 a 的更新公式
def update_a(X, Y, a, b, eta):
    s = 0
    for i in range(X.size):
        s += ((a*X[i] + b)-Y[i])*X[i]
    return a - eta * s * (1./Y.size)

# 参数 b 的更新公式
def update_b(X, Y, a, b, eta):
    s = 0
    for i in range(X.size):
        s += ((a*X[i] + b)-Y[i])
    return b - eta * s * (1./Y.size)

我们选取参数的初始值 $(a, b) = (0, 0)$ ，并令学习率 $\eta = 0.01$ ，然后进行第一次迭代，并计算损失函数的值和上次损失函数的值之间的差值（这个差值可以用于判断梯度下降过程是否已经收敛，如果两次损失函数的差值足够小，就认为收敛）：

a = 0
b = 0
eta = 0.01

a_new = update_a(X, Y, a, b, eta)
b_new = update_b(X, Y, a, b, eta)
loss_new = loss(X, Y, a_new, b_new)
diff = abs(loss(X, Y, a, b) - loss(X, Y, a_new, b_new))
print("a = {0}, b = {1}, loss = {2}, diff = {3}".format(a_new, b_new, loss_new, diff))

运行结果

a = 0.04433333333333332, b = 0.03033333333333333, loss = 52.07898433333334, diff = 2.566015666666658

我们不妨画出迭代一次后的图像，很显然，这个差值还比较大，我们需要继续这个过程：

plt.plot(X, Y, 'k.')

x = np.linspace(0, 10, 10)
y = a_new*x+b_new

plt.plot(x, y, 'g-')
plt.show()

运行结果

Notebook 运行结果

我们把这个过程改为迭代 10 次：

for i in range(10):
    a_new = update_a(X, Y, a, b, eta)
    b_new = update_b(X, Y, a, b, eta)
    loss_new = loss(X, Y, a_new, b_new)
    diff = abs(loss(X, Y, a, b) - loss(X, Y, a_new, b_new))
    print("a = {0}, b = {1}, loss = {2}, diff = {3}".format(a_new, b_new, loss_new, diff))
    a = a_new
    b = b_new

运行结果

a = 0.04433333333333332, b = 0.03033333333333333, loss = 52.07898433333334, diff = 2.566015666666658
a = 0.08762444444444442, b = 0.059919999999999994, loss = 49.63396493034815, diff = 2.4450194029851886
a = 0.12989817037037035, b = 0.08877788888888888, loss = 47.30423372457149, diff = 2.3297312057766604
a = 0.17117875530864196, b = 0.11692446162962962, loss = 45.08435194574915, diff = 2.219881778822341
a = 0.21148986477053494, b = 0.14437676279358025, loss = 42.969137414855254, diff = 2.115214530893894
a = 0.25085459939642357, b = 0.17115142985127244, loss = 40.95365243850431, diff = 2.0154849763509475
a = 0.2892955084413038, b = 0.1972647028921288, loss = 39.03319227448755, diff = 1.9204601640167596
a = 0.32683460293836075, b = 0.22273243411212781, loss = 37.20327414148798, diff = 1.8299181329995662
a = 0.3634933685482668, b = 0.24757009707495625, loss = 35.45962674730114, diff = 1.7436473941868442
a = 0.3992927781017128, b = 0.27179279575205734, loss = 33.79818031109894, diff = 1.6614464362021977

可以看到 loss 和 diff 都变得越来越小，说明损失函数在不断收敛，画出现在的图像：

plt.plot(X, Y, 'k.')

x = np.linspace(0, 10, 10)
y = a_new*x+b_new

plt.plot(x, y, 'g-')
plt.show()

运行结果

Notebook 运行结果

比迭代一次靠谱多了，我们再继续这个过程，迭代 1000 次：

for i in range(1000):
    a_new = update_a(X, Y, a, b, eta)
    b_new = update_b(X, Y, a, b, eta)
    loss_new = loss(X, Y, a_new, b_new)
    diff = abs(loss(X, Y, a, b) - loss(X, Y, a_new, b_new))
    a = a_new
    b = b_new
    
print("a = {0}, b = {1}, loss = {2}, diff = {3}".format(a_new, b_new, loss_new, diff))

运行结果

a = 2.083663961231524, b = 0.9559983763104101, loss = 0.010980851812161636, diff = 5.500906925493335e-06

diff 的值已经趋近于 0，说明已经达到收敛，这个时候的图形，就是我们根据梯度下降法求得的结果：

plt.plot(X, Y, 'k.')

x = np.linspace(0, 10, 10)
y = a_new*x+b_new

plt.plot(x, y, 'g-')
plt.show()

运行结果

Notebook 运行结果

多元线性回归的梯度下降

上面是对一元线性回归的梯度下降法的介绍，回到线性回归的一般形式：

\hat{y} = \theta_0 + \theta_1x_1 + \theta_2x_2 + ... + \theta_nx_n = \theta^Tx

其对应的损失函数为：

loss = J(\theta_0, \theta_1, \dots, \theta_n) = J(\theta) = \frac{1}{2m} \sum_{i=1}^{m}(y^{(i)} - \theta^Tx^{(i)})^2

对 $J(\theta)$ 进行求导，得到参数 $\theta$ 的更新公式：

\left\{ \begin{align} \theta_0 &:= \theta_0 - \eta \frac{1}{m} \sum_{i=1}^{m}(y^{(i)} - \theta^Tx^{(i)})x_0^{(i)} \\ \theta_1 &:= \theta_1 - \eta \frac{1}{m} \sum_{i=1}^{m}(y^{(i)} - \theta^Tx^{(i)})x_1^{(i)} \\ \theta_2 &:= \theta_2 - \eta \frac{1}{m} \sum_{i=1}^{m}(y^{(i)} - \theta^Tx^{(i)})x_2^{(i)} \\ \vdots \\ \theta_n &:= \theta_n - \eta \frac{1}{m} \sum_{i=1}^{m}(y^{(i)} - \theta^Tx^{(i)})x_n^{(i)} \\ \end{align} \right.

其中， $m$ 表示样本个数， $n$ 表示特征维度， $x^{(i)}$ 表示第 i 个样本， $x^{(i)}_n$ 表示第 i 个样本的第 n 个特征。为了方便表达和计算，损失函数也可以写成下面的矩阵形式：

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

根据前面的矩阵求导，得到梯度下降的更新公式为：

\begin{aligned} \theta &:= \theta - \eta \frac{\partial}{\partial\theta}J(\theta) \\ &= \theta - \eta \frac{1}{m} X^T(\rm{X}\theta-\bf{y}) \end{aligned}

TODO：使用矩阵的迹（Trace）求矩阵导数。

我们把上面的代码改成矩阵计算的形式：

import numpy as np
import matplotlib.pyplot as plt

# 计算损失函数 J 的梯度
def gradient_function(theta, X, y):
    diff = np.dot(X, theta) - y
    return (1./y.size) * np.dot(np.transpose(X), diff)

# 梯度下降
def gradient_descent(X, y, eta, theta_init, max_iter):
    theta = theta_init
    i = 0
    gradient = gradient_function(theta, X, y)
    while not np.all(np.absolute(gradient) <= 1e-10):
        theta = theta - eta * gradient
#         print('theta:', theta[0,0], theta[1,0])
        i += 1
        if i > max_iter:
            break
        gradient = gradient_function(theta, X, y)
#         print('gradient:', gradient[0,0], gradient[1,0])
    return theta

# 损失函数
def loss_function(theta, X, y):
    diff = np.dot(X, theta) - y
    return (1./2*y.size) * np.dot(np.transpose(diff), diff)

# 数据集
X0 = np.ones((3, 1))
X1 = np.array([0,1,2]).reshape(-1, 1)
X = np.hstack((X0, X1))
y = np.array([0.9, 3.1, 5.1]).reshape(-1, 1)

# 学习率
eta = 0.01

# 最大迭代次数
max_iter = 10000

# 初始值
theta_init = np.array([0, 0]).reshape(-1, 1)

# 使用梯度下降求解线性回归
theta = gradient_descent(X, y, eta, theta_init, max_iter)
print('theta:', theta[0,0], theta[1,0])
print('loss:', loss_function(theta, X, y))

# 绘制图形
plt.plot(X1, y, 'k.')
x = np.linspace(0, 10, 10)
y = theta[0,0]*x+theta[1,0]
plt.plot(x, y, 'g-')
plt.show()

运行结果

theta: 0.9333333336913473 2.099999999741958
loss: [[0.01]]

运行结果

Notebook 运行结果

疑惑一

为什么上面两种方法求出来的解不一样？

疑惑二

在前面介绍正规方程法解一元线性回归时，使用了下面的数据集，这个数据集不知道为什么用梯度下降法无法收敛？

import numpy as np
import matplotlib.pyplot as plt

# 计算损失函数 J 的梯度
def gradient_function(theta, X, y):
    diff = np.dot(X, theta) - y
    return (1./y.size) * np.dot(np.transpose(X), diff)

# 梯度下降
def gradient_descent(X, y, eta, theta_init, max_iter):
    theta = theta_init
    i = 0
    gradient = gradient_function(theta, X, y)
    while not np.all(np.absolute(gradient) <= 1e-10):
        theta = theta - eta * gradient
        i += 1
        if i > max_iter:
            break
        gradient = gradient_function(theta, X, y)
    return theta

# 损失函数
def loss_function(theta, X, y):
    diff = np.dot(X, theta) - y
    return (1./2*y.size) * np.dot(np.transpose(diff), diff)

# 数据集
X0 = np.ones((6, 1))
X1 = np.array([39.93, 42.05, 43.18, 44.68, 49.87, 53.57]).reshape(-1, 1)
X = np.hstack((X0, X1))
y = np.array([199,   290,   298,   310,   399,   420]).reshape(-1, 1)

# 学习率
eta = 0.00001

# 最大迭代次数
max_iter = 10000

# 初始值
theta_init = np.array([1, 1]).reshape(-1, 1)

# 使用梯度下降求解线性回归
theta = gradient_descent(X, y, eta, theta_init, max_iter)
print('theta:', theta[0,0], theta[1,0])
print('loss:', loss_function(theta, X, y))

# 绘制图形
plt.plot(X1, y, 'k.')
x = np.linspace(30, 60, 100)
y = theta[0,0]*x+theta[1,0]
plt.plot(x, y, 'g-')
plt.show()

运行结果

theta: 0.7489393245009788 7.078701795630663
loss: [[32892.8042478]]

运行结果

Notebook 运行结果

一元线性回归的梯度下降#

多元线性回归的梯度下降#

疑惑一#

疑惑二#

一元线性回归的梯度下降

多元线性回归的梯度下降

疑惑一

疑惑二