用梯度下降法简单实现一元线性拟合

核心概念

• 梯度：函数在某一点的最大方向导数，表示函数在该点上升最快的方向

• 梯度向量：函数的某一点处沿着一个向量的方向导数最大（即取得梯度），则该向量为梯度向量

• 训练轮数：迭代次数，完整的使用训练集中的所有样本进行一次训练称之为一轮训练

• 批次大小：每次迭代使用的样本数

• 学习率：每次迭代的步长

• 损失函数：用于评估模型的好坏的某一种误差度量函数

• 梯度下降法：沿着梯度向量的反方向移动，直至找到函数的最小值。在训练模型中，这里的梯度指损失函数的梯度

数学原理

已知数据集：

$$D=\{(x_1,y_1),(x_2,y_2),\cdots ,(x_n,y_n)\}$$

假如我们需要拟合一条直线，那模型的方程应该是：

$$f(x)=kx+b$$

我们可以通过最小化残差平方和来确定$k$和$b$的值，这样一来损失函数：

$$g(k,b)=\sum _{i=1}^n(f(x_i)-y_i{)}^2=\sum _{i=1}^n(k{x}_i+b-y_i{)}^2$$

这里需要知道的是原本模型的参数就作为损失函数的自变量。

求损失函数的偏导数和梯度向量：

$$\begin{array}{rl}\frac{\partial g(k,b)}{\partial k}& =\sum _{i=1}^{n}2(f(x_i)-y_i)\cdot x_i\\ \frac{\partial g(k,b)}{\partial b}& =\sum _{i=1}^{n}2(f(x_i)-y_i)\\ \nabla g(k,b)& =\left[\begin{array}{c}\frac{\partial g(k,b)}{\partial k}\\ \frac{\partial g(k,b)}{\partial b}\end{array}\right]\end{array}$$

假如训练轮数是$E$，学习率是$\eta$，那么梯度下降法的迭代公式是：

$$\left[\begin{array}{c}{k}_i\\ b_i\end{array}\right]=\left[\begin{array}{c}{k}_{i-1}\\ b_{i-1}\end{array}\right]-\eta \nabla g(k_{i-1},b_{i-1})$$

在上述步骤中需要初始化$k=k_0$和$b=b_0$，然后不断迭代$E$次，时间复杂度是$O(n\cdot E)$。

代码实现

Python代码如下：

import matplotlib.pyplot as plt
import random

# 模拟梯度下降法
example_func = lambda x: 25 * x - 19                  # 模板函数，样本点添加了噪声后可能会有所偏差

points = [(x, example_func(x) + random.gauss(0, 1)) for x in (random.uniform(-10, 10) for _ in range(100))]  # 生成数据集
k, b = 5, 3                                           # f(x) = k * x + b，初始化时k和b的值随意设置
learning_rate = 5e-6                                  # 学习率
epochs = 200000                                       # 训练轮数

# 损失函数: g(k, b) = \sum_{i=1}^{n} (f(x_i) - y_i)^2   # 残差平方和

for epoch in range(epochs):                           # O(N * E)
    k_grad, b_grad = 0, 0
    for x, y in points:                               # 把整个数据集看作一个batch
        k_grad += 2 * (k * x + b - y) * x             # d(g(k, b)) / d(k)
        b_grad += 2 * (k * x + b - y)                 # d(g(k, b)) / d(b)
    # 沿着梯度向量的反方向移动，梯度向量：[[k_grad], [b_grad]]
    k -= learning_rate * k_grad
    b -= learning_rate * b_grad
    if epoch % 20000 == 0:
        print(f"Epoch {epoch}: k = {k}, b = {b}")
print(f"Final: k = {k}, b = {b}")

# 最小二乘法对照
x_avg = sum(x for x, _ in points) / len(points)
y_avg = sum(y for _, y in points) / len(points)
k_best = sum((x - x_avg) * (y - y_avg) for x, y in points) / sum((x - x_avg) ** 2 for x, _ in points)
b_best = y_avg - k_best * x_avg
print(f"Best: k = {k_best}, b = {b_best}")

plt.plot([x for x, _ in points], [y for _, y in points], 'ro')           # 散点
plt.plot([x for x, _ in points], [k * x + b for x, _ in points], '-')    # 直线
plt.show()

需要注意的是，这里没有划分训练集和验证集，因为可以用最小二乘法来验证模型的好坏。
代码中的example_func是一个模板函数，由于后续添加了噪声，所以拟合的直线不会完全与模板函数重合。
代码中的learning_rate和epochs是需要调整的超参数。如果学习率过高，可能会导致梯度爆炸。
代码中的points是一个数据集，这里的数据集是随机生成的，实际应用中应该是真实的数据集。
代码处理可视化的部分使用了matplotlib库，需要提前安装。其余部分均是Python内置的库。
代码中的计算过程均在CPU完成，仅供学习参考。

控制台输出：

Epoch 0: k = 5.617168567733822, b = 2.9736399575387327
Epoch 20000: k = 25.004648079906683, b = -19.018192641629117
Epoch 40000: k = 25.00464807956941, b = -19.018192687654867
Epoch 60000: k = 25.00464807956941, b = -19.018192687654867
Epoch 80000: k = 25.00464807956941, b = -19.018192687654867
Epoch 100000: k = 25.00464807956941, b = -19.018192687654867
Epoch 120000: k = 25.00464807956941, b = -19.018192687654867
Epoch 140000: k = 25.00464807956941, b = -19.018192687654867
Epoch 160000: k = 25.00464807956941, b = -19.018192687654867
Epoch 180000: k = 25.00464807956941, b = -19.018192687654867
Final: k = 25.00464807956941, b = -19.018192687654867
Best: k = 25.004648079569343, b = -19.018192687656658

绘制的图像如下：