[AI] 基于梯度下降法的一次函数拟合

前言

分析了下自己的技能树，打算点一些 AI 方向的 perk。久闻吴恩达机器学习的大名，因此先从他的教程看起。

目前学到梯度下降，但还没有敲过任何代码——正好今天时间比较充裕，遂决定进行一些小小的实践。

准备 training set

首先需要生成 training set。

\[y=wx+b+\epsilon \]

如上。我的想法是，先随机生成 \(w\) 和 \(b\)，对于在 \([xMin, xMax]\) 等间距取值的 \(x\)，计算 \(y=wx+b\)，这样 \(y\) 和 \(x\) 就成一次函数关系了。

不过实际情况一般会更复杂一些——数据并非严格地遵循一次函数的规律，而是近似地符合该规律，因此我在 \(y\) 上面叠加了一个扰动量 \(\epsilon\)，并引入了扰动因子 \(disturbanceFac\) 来控制扰动的大小。

当扰动因子为 \(1\) 时，生成的数据如下：

梯度下降

我们的模型可以表示为 \(f(x)=w_1x+b_1\)，使用平方误差代价函数，得到：

\[J(w_1,b_1)=\frac{1}{2m}\sum\limits_{i=1}^{m}(w_1x^{(i)}+b_1-y^{(i)})^2 \]

分别计算出 \(J(w_1,b_1)\) 对 \(w_1\) 和 \(b_1\) 的偏导数：

\[\frac{\partial{J}}{\partial{w_1}}=\frac{1}{m}\sum\limits_{i=1}^{m}(w_1x^{(i)}+b_1-y^{(i)})x^{(i)} \]

\[\frac{\partial{J}}{\partial{b_1}}=\frac{1}{m}\sum\limits_{i=1}^{m}(w_1x^{(i)}+b_1-y^{(i)}) \]

化简得到：

\[\frac{\partial{J}}{\partial{w_1}}=\frac{1}{m}(w_1\sum\limits_{i=1}^{m}{x^{(i)}}^2+b_1\sum\limits_{i=1}^{m}x^{(i)}-\sum\limits_{i=1}^{m}x^{(i)}y^{(i)}) \]

\[\frac{\partial{J}}{\partial{b_1}}=\frac{1}{m}(w_1\sum\limits_{i=1}^{m}x^{(i)}+mb_1-\sum\limits_{i=1}^{m}y^{(i)}) \]

其中的 \(\sum\limits_{i=1}^{m}x^{(i)}\)、\(\sum\limits_{i=1}^{m}y^{(i)}\)、\(\sum\limits_{i=1}^{m}{x^{(i)}}^2\)、\(\sum\limits_{i=1}^{m}x^{(i)}y^{(i)}\) 可以通过预计算得到，不需要重复计算。

每次 epoch，我们遵循如下规则对 \(w_1\) 和 \(b_1\) 进行更新：

\[\begin{cases} w_1\leftarrow{}w_1-\alpha{}\frac{\partial{J}}{\partial{w_1}}\\ b_1\leftarrow{}b_1-\alpha{}\frac{\partial{J}}{\partial{b_1}} \end{cases} \]

重复上述过程，直至收敛。

教程中并没有给出判断收敛的定量指标，这里我采取的策略是检测 \(|\frac{\partial{J}}{\partial{w_1}}|\) 和 \(|\frac{\partial{J}}{\partial{b_1}}|\) 是否均小于预先设定的阈值 \(convThreshold\)，如果是则认为收敛，否则继续迭代。

完整代码

import matplotlib.pyplot as plt
import numpy as np
import random

# 0. initial params
wMin, wMax = 0, 10.0
bMin, bMax = 0, 10.0
xMin, xMax = 1, 10
m = 16
alpha = 0.001
convThreshold = 0.001
disturbanceFac = 1

# 1. generate training set
# y = wx + b + eps
# eps: disturbance, range [-epsK, epsK], epsK = disturbanceFac * (xMax - xMin) * w/ (2 * (m - 1))
w = random.uniform(wMin, wMax)
b = random.uniform(bMin, bMax)
epsK = disturbanceFac * (xMax - xMin) * w / (2 * (m - 1))
eps = np.array([random.uniform(-epsK, epsK) for _ in range(m)])

print("(w, b) = ({}, {})".format(w, b))

x = np.linspace(xMin, xMax, m)
y = w * x + b + eps

plt.scatter(x, y, marker="x", c="r")
plt.grid()

print("Press enter to start training..")
input()

# 2. perform gradient descent algorithm to predict w, b
# w1, b1 are estimated values of w, b
w1, b1 = 0.0, 0.0

# 2.1 precompute some reusable terms
xSum = 0
ySum = 0
xSqSum = 0
xySum = 0
for i in range(m):
    xSum += x[i]
    ySum += y[i]
    xSqSum += x[i] * x[i]
    xySum += x[i] * y[i]

# 2.2 perform algorithm
epoch = 0
partialW1 = 100
partialB1 = 100
while abs(partialW1) > convThreshold or abs(partialB1) > convThreshold:
    print("epoch #{}".format(epoch))
    print("{} {}".format(w1, b1))

    partialW1 = 1 / m * (w1 * xSqSum + b1 * xSum - xySum)
    partialB1 = 1 / m * (w1 * xSum + m * b1 - ySum)

    w1, b1 = w1 - alpha * partialW1, b1 - alpha * partialB1

    epoch += 1

# 3. display result
x1 = np.linspace(xMin, xMax, m)
y1 = w1 * x1 + b1
plt.plot(x1, y1)
plt.show()

测试结果

\(disturbanceFac=1\)

\(disturbanceFac=3\)

\(disturbanceFac=5\)

一点感悟

话说 C 语言的 HelloWorld 就是打印 HelloWorld，OpenGL 的 HelloWorld 是绘制三角形，单片机的 HelloWorld 是点灯……本文大概就是 AI 的 HelloWorld了。

放在以前，要想求解这样的问题，我只能用暴力搜索解决，然而一旦试图去使用暴搜，就会陷入这样的困境：搜索的范围该设多大？搜索的精度（步进长度）又该是多少？对于任意取值的 \(w,b\) 来说，搜索的范围应当是 \((-\infty{},\infty{})\)，为了获得尽可能精确的结果，步进长度应当 \(\rightarrow0\)，这在程序上根本无法实现……

但是，用梯度下降法 “一瞬间” 就能得到精度极高的结果。

简直——太神奇了!

2024-09-12 UPD

突然想到，目前唯一的问题在于它是一种启发式算法——启发式算法都很玄学，因为其运行效率依赖于数据的特征以及选取的参数，不具备时间复杂度这样能相对确切地描述算法运行效率的指标。