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

前言

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

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

准备 training set

首先需要生成 training set。

$$y=wx+b+ϵ$$

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

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

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

梯度下降

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

$$J(w_1,b_1)=\frac{1}{2m}\sum _{i=1}^m(w_1{x}^{(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 _{i=1}^m(w_1{x}^{(i)}+b_1-y^{(i)})x^{(i)}$$

$$\frac{\partial J}{\partial b_1}=\frac{1}{m}\sum _{i=1}^m(w_1{x}^{(i)}+b_1-y^{(i)})$$

化简得到：

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

$$\frac{\partial J}{\partial b_1}=\frac{1}{m}(w_1\sum _{i=1}^m{x}^{(i)}+m{b}_1-\sum _{i=1}^m{y}^{(i)})$$

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

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

$$\{\begin{array}{l}{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{array}$$

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

教程中并没有给出判断收敛的定量指标，这里我采取的策略是检测 $|\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$ 来说，搜索的范围应当是 $(-\mathrm{\infty },\mathrm{\infty })$，为了获得尽可能精确的结果，步进长度应当 $\to 0$，这在程序上根本无法实现……

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

简直——太神奇了!

2024-09-12 UPD

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