线性回归 - 宝可梦案例

线性回归 - 宝可梦案例

一、案例分析

假设有10个已知的 x_data 和 y_data ，x 和 y之间的关系是
y_data = b + w * x_data ，需要通过线性回归学习得到w与b。
其中已知：
进化前cp值，x_data为[338., 333., 328., 207., 226., 25., 179., 60., 208., 606.]；
进化后cp值，y_data为[640., 633., 619., 393., 428., 27., 193., 66., 226., 1591.]

模型使用：linear model
工具使用：jupyter

二、案例实现

（1）准备工作

需要用到的库有：numpy、 matplotlib.pyplot

• 导入库函数

import numpy as np
import matplotlib.pyplot as plt

• Matplotlib没有中文字体等的动态解决方法

plt.rcParams['font.sans-serif'] = ['Simhei']  # 显示中文
from pylab import mpl
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题

（2）进行数据的导入

# 进化前的cp值
x_data = [338., 333., 328., 207., 226., 25., 179., 60., 208., 606.]
# 进化后的cp值
y_data = [640., 633., 619., 393., 428., 27., 193., 66., 226., 1591.]

（3）绘制损失函数的热点图

什么是损失函数热点图？？
即通过将所有的可能范围内的w和b的组合情况下的损失结果直接在图中通过热力图的形式显示出来，并作为之后的训练过程的背景。
达到的效果：
线性回归中常用的损失函数 Loss function 是：

L(w,b) = ${\sum }_{n=1}^{10}({\check{y}}^n-(b+w\cdot x_{cp}){)}^2$

x = np.arange(-200, -100, 1)   # x轴坐标-200到200，间距为1
y = np.arange(-5, 5, 0.1)      # y轴坐标-5到5，间距为0.1
#按照x，y轴的长度生成一个二维矩阵，元素全0填充
L = np.zeros((len(x), len(y)))
# 损失值Loss
for i in range(len(x)):
    for j in range(len(y)):
        b = x[i]  # b表示横轴
        w = y[j]  # w表示纵轴
        # 每一个w和b组合 => 一个Loss Function
        for n in range(len(x_data)):
            L[j][i] = L[j][i] + (y_data[n] - b - w * x_data[n] ) ** 2 # 损失累加
    # 计算这种w和b的组合下，损失的平均值
        L[j][i]  /=  len(x_data)

(4)线性回归梯度下降计算w和b

需要计算 Loss 关于 w 和 b 的偏微分，模拟梯度下降过程。
由Loss function推导公式如下：

$\frac{\partial L}{\partial w}={\sum }_1^{n}2\cdot [y\_true-(w\cdot x^i+b)]\cdot (-x^i)$

$\frac{\partial L}{\partial b}={\sum }_1^{n}2\cdot [y\_true-(w\cdot x^i+b)]\cdot (-1)$

# linear regression
#1、给参数赋任意初始值
b = -120    #bias
w = -4      #表示x特征的权重
lr = 0.0000001 #学习率
iteration = 100000 #迭代次数
b_history = [b]   #保存训练过程的所有b
w_history = [w]  #保存训练过程的所有w
#2、模拟梯度下降过程
for i in range(iteration):
    grad_b = 0.0 # b的偏微分
    grad_w = 0.0 # w的偏微分
    #求w和b的偏导
    for j in range(len(x_data)):
        grad_b += 2 * (y_data[j] - w * x_data[j] -b) * (-1)
        grad_w += 2 * (y_data[j] - w * x_data[j] -b) * (-x_data[j])
    # update w、b
    b -= lr * grad_b
    w -= lr * grad_w
#保存参数w、b的历史值
    b_history.append(b)
    w_history.append(w)

（5）绘制结果图

利用 matplotlib 里的函数即可

# 填充背景，按照损失的大小绘制等高线
plt.contourf(x, y, L, 50, alpha=0.5, cmap=plt.get_cmap('jet'))
# 绘制出最佳结果点，用x标记
plt.plot([-188.4], [2.67], 'x', ms=12, mew=3, color='orange')
#绘制出w和b历史轨迹
plt.plot(b_history, w_history, 'o-', ms=3, lw=1.5, color='black')
# 图像的基本参数设置
plt.xlim(-200, -100)
plt.ylim(-5, 5)
plt.xlabel('b')
plt.ylabel('w')
plt.title('线性回归')
plt.show()

三、结果显示及优化

（1）第一次结果

这张图是 学习率 Ir 为 0.0000001时 损失函数随w、b变化图，黑色线区域为计算的 w、b的历史轨迹，黄色叉叉处是我们需要得到的最优解。

结果可知，经过100000次的update，我们最终的参数离最佳解仍然非常的遥远。
结论：学习率learning rate不够大，将其增大到Ir = 0.000001（初始的10倍）

（2）优化后第二次结果

结果可知，经过100000次的update，参数离最佳解更近了，但出现了剧烈的震荡，且仍未达到最佳解。

结论：学习率learning rate不够大，将其增大到Ir = 0.00001（初始的100倍）

（3）优化后第三次结果

结果可知，我们的损失函数变化过程中，并没有和理想中的一样 靠近最佳解，反而出现了很大的震荡，说明由于学习速率过大，导致更新参数的时候，一次性步长过大，就很难收敛到我们设定好的理想损失值处，从而呈现出了发散状态。

结论：学习速率过大了，因此上一次直接增大lr为初始值的100倍是不合理的。

解决方案：

• 很明显，若直接将学习率增大到100倍，会导致学习率过大，但是增大到10倍则学习率过小，可以尝试继续采用更新学习率的方式，一步一步的调整学习率，缩小范围，从而得到最佳解。
• 直接采用 Adagrad 算法更新参数 ，怎么实现呢？接下来，我们来学习尝试一下。

四、采用Adagrad 算法更新参数

（1）什么是Adagrad 算法？

直观看公式：
如果常规的梯度下降表示方法是：

$w^{t+1}=w^t-{\eta }^t\cdot g^t$

那么Adagrad算法的梯度下降就是：

$w^{t+1}=w^t-\frac{{\eta }^t}{{\sigma }^t}\cdot g^t$

通俗来讲：Adagrad算法就是加入了 σ ，使得原本不变的学习率开始时刻随着 w 与 b 的梯度变化着。

所以公式中的 ${\sigma }^t$的含义就是：之前参数的所有微分的均方根

（2）Adagrad算法举例及化简

$w^1⟵w^0⟵\frac{{\eta }^0}{{\sigma }^0}\cdot g^0{\sigma }^0=\sqrt{(g^0{)}^2}$
$w^2⟵w^1⟵\frac{{\eta }^1}{{\sigma }^1}\cdot g^1{\sigma }^1=\sqrt{\frac{1}{2}[(g^0{)}^2+(g^1{)}^2]}$
$w^3⟵w^2⟵\frac{{\eta }^2}{{\sigma }^2}\cdot g^2{\sigma }^2=\sqrt{\frac{1}{3}[(g^0{)}^2+(g^1{)}^2+(g^2{)}^2]}$
$\downarrow \downarrow \downarrow \downarrow$
$w^{t+1}⟵w^t⟵\frac{{\eta }^t}{{\sigma }^t}\cdot g^t{\sigma }^t=\sqrt{\frac{1}{t+1}{\sum }_{i=0}^t(g^i{)}^2}$

注意：根号中累加后需要除以计数，然后才开根号（即求均方根）

最后公式可以化简得到：$w^{t+1}=w^t-\frac{\eta }{\sqrt{{\sum }_{i=0}^t(g^i{)}^2}}\cdot g^t$
(化简步骤在本文后的总结疑问中会提到)

（3）通过Adagrad算法实现参数更新

# 1.给参数附任意初始值
b = -120
w = -4
lr = 1  # 学习率
iteration = 100000  # 迭代次数

b_history = [b]  # 保存训练过程的所有b
w_history = [w]  # 保存训练过程的所有w

lr_w = 0
lr_b = 0
# 2.模拟回归过程
for i in range(iteration):
    grad_w = 0.0  # w的偏微分
    grad_b = 0.0  # b的偏微分
    for j in range(len(x_data)):
        grad_w += 2 * (y_data[j] - w * x_data[j] - b) * (-x_data[j])
        grad_b += 2 * (y_data[j] - w * x_data[j] - b) * (-1)

    # 采用化简后的Adagrad算法公式更新参数
    lr_w = lr_w + grad_w ** 2
    lr_b = lr_b + grad_b ** 2
    w = w - lr / np.sqrt(lr_w) * grad_w
    b = b - lr / np.sqrt(lr_b) * grad_b

    # 保存历史参数
    w_history.append(w)
    b_history.append(b)

（4）使用算法后的结果展示

该图是采用Adagrad算法更新参数函数的损失变化图，如图可以直观看到，最终我们的参数达到了我们期望的最优解，即完成目标。

五、案例总结

（1）个人感受

本次案例 通过少量的数据 重现了线性回归的模型

涉及的知识点：

• 损失函数Loss
• 梯度下降（Gradient Descent）
• • Loss关于w、b的偏微分公式
• Numpy、matplotlib函数的使用
• Adagrad算法

通过这次的案例，让我体会了线性回归的整个过程，更深刻得理解了其中w、b的update过程（梯度下降），以及利用梯度下降和Adagrad算法找到最优解的过程，很奇妙。

（2）过程遇到的疑惑及解答

1、为什么Adagrad算法可以达到目的，算法如何化简得到公式：$w^{t+1}=w^t-\frac{\eta }{\sqrt{{\sum }_{i=0}^t(g^i{)}^2}}\cdot g^t$，化简后的公式有什么用？
答：
普通的梯度下降需要人工手动调整学习率，很明显这是一个麻烦的方法，因为学习率过大或者过小都是不行的，很难手动找到那个合适的值。
而Adagrad算法能够让学习率随着迭代次数增加而变小，初始迭代时，使用较大的学习率加速下降；迭代几次后，减小学习率防止振荡和越过；且对每个参数都使用了不同的却合适的学习率。
公式化简得到 $w^{t+1}=w^t-\frac{\eta }{\sqrt{{\sum }_{i=0}^t(g^i{)}^2}}\cdot g^t$：
初始公式：$w^{t+1}=w^t-\frac{{\eta }^t}{{\sigma }^t}\cdot g^t$

而 ${\eta }^t=\frac{\eta }{\sqrt{t+1}}{\sigma }^t=\sqrt{\frac{1}{t+1}\cdot {\sum }_{i=0}^t(g^i{)}^2}$

将它们代入初始公式中，化简相约，即可得到化简公式：

$w^{t+1}=w^t-\frac{\eta }{\sqrt{{\sum }_{i=0}^t(g^i{)}^2}}\cdot g^t$
而我们就是利用了化简后的公式去计算更新参数的：

# 采用化简后的Adagrad算法公式更新参数
lr_w = lr_w + grad_w ** 2
lr_b = lr_b + grad_b ** 2
w = w - lr / np.sqrt(lr_w) * grad_w
b = b - lr / np.sqrt(lr_b) * grad_b

2、初始加入数据时，数字后 有点 和 没点 有什么区别吗？

答：
其实很简单，有点表示float类型，没点表示int类型，然而在我们的计算中，要保证精确，自然是使用float类型，因此，加上点更合适。
3、除了用常规公式求偏导，其实还有另一种（下图）求偏导的方式，是怎么做到的呢？

答：
常规求偏导是：

对比两种方法，一种用了循环逐一求，一种直接使用了函数对矩阵进行计算，说明一下w的求导变换：
变换后的是：grad_w = -2.0 * np.dot(y_d - y_hat, x_d)
其中：y_hat = w * x_data[j] -b，表示的是预测的y值，则y_d - y_hat就是实际值与预测值间的差值，y_d - y_hat = y_data[j] - w * x_data[j] -b，np.dot()将y_d - y_hat与x_d每一项分别相乘再相加，即达到了(y_data[j] - w * x_data[j] -b) * x_data[j]的效果，最后乘以-2 即 得到2 * (y_data[j] - w * x_data[j] -b) * (-x_data[j])。
b的求导变换也类似。
所以两种方法的式子最后结果都一样，只是运用的方法不同，常规求导更好理解，但另一种方法显得更高级，各有千秋。

六、完整代码

import numpy as np
import matplotlib.pyplot as plt
# matplotlib没有中文字体，动态解决
plt.rcParams['font.sans-serif'] = ['Simhei']  # 显示中文
from pylab import mpl
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
# 进化前的cp值
x_data = [338., 333., 328., 207., 226., 25., 179., 60., 208., 606.]
# 进化后的cp值
y_data = [640., 633., 619., 393., 428., 27., 193., 66., 226., 1591.]
x = np.arange(-200, -100, 1)   # x轴坐标-200到200，间距为1
y = np.arange(-5, 5, 0.1)      # y轴坐标-5到5，间距为0.1
#按照x，y轴的长度生成一个二维矩阵，元素全0填充
L = np.zeros((len(x), len(y)))
# 损失值Loss
for i in range(len(x)):
    for j in range(len(y)):
        b = x[i]  # b表示横轴
        w = y[j]  # w表示纵轴
        # 每一个w和b组合 => 一个Loss Function
        for n in range(len(x_data)):
            L[j][i] = L[j][i] + (y_data[n] - b - w * x_data[n] ) ** 2 # 损失累加
    # 计算这种w和b的组合下，损失的平均值
        L[j][i]  /=  len(x_data)
# 1.给参数附任意初始值
b = -120
w = -4
lr = 1  # 学习率
iteration = 100000  # 迭代次数

b_history = [b]  # 保存训练过程的所有b
w_history = [w]  # 保存训练过程的所有w

lr_w = 0
lr_b = 0
# 2.模拟回归过程
for i in range(iteration):
    grad_w = 0.0  # w的偏微分
    grad_b = 0.0  # b的偏微分
    for j in range(len(x_data)):
        grad_w += 2 * (y_data[j] - w * x_data[j] - b) * (-x_data[j])
        grad_b += 2 * (y_data[j] - w * x_data[j] - b) * (-1)

    # 采用化简后的Adagrad算法公式更新参数
    lr_w = lr_w + grad_w ** 2
    lr_b = lr_b + grad_b ** 2
    w = w - lr / np.sqrt(lr_w) * grad_w
    b = b - lr / np.sqrt(lr_b) * grad_b

    # 保存历史参数
    w_history.append(w)
    b_history.append(b)
# 绘制变化图
# 填充等高线
plt.contourf(x, y, L, 50, alpha=0.5, cmap=plt.get_cmap('jet'))
plt.plot([-188.4], [2.67], 'x', ms=12, mew=3, color='orange')
plt.plot(b_history, w_history, 'o-', ms=3, lw=1.5, color='black')
plt.xlim(-200, -100)
plt.ylim(-5, 5)
plt.xlabel('b')
plt.ylabel('w')
plt.title('线性回归')
plt.show()