机器学习：第三篇使用随机梯度下降法(SGD)求解方程

使用随机梯度下降法(SGD)求解一元三次方程的根。在知乎上看的。

求解分析

原题目是求解: $a^3+a^2$ 这里为了方便于我们理解，我把方程简化成2元方程： $a^2 = 49$ 。使用随机梯度下降法求解方程。

我们可以把我们平时遇到的神经网络看成是一个函数 $y=f_\theta(x)$ ，其中y是神经网络的输出， $\theta$ 是神经网络的参数，而x是神经网络的输入。

我们就可以把题目看成是有一个神经网络，假定函数为 $y=f_\theta(x)=\theta^2-49$ ，要求的这样一个最佳的 $\theta$ 使得无论输入x的任意值，使得y都逼近于0。

我们首先需要定义一个损失函数来对参数进行评估： $cost(x)=f_\theta^2(x)$ ，当cost(x)的值越小，那么 $\theta$ 的表征预测值和准确值就越接近。

神经网络的目标就是求得一个 $\theta$ 使得上面的损失函数最小。

使用随机梯度下降法迭代求取 $\theta_n$ ，直到我们的损失函数取得范围内的最小值。

$\theta_n=\theta_{n-1}-lr*grad(cost(x))$
其中lr表示学习率，就是梯度下降的速率，而grad表示cost对 $\theta$ 的梯度。

分析总结

$y=f_\theta(x)=\theta^2-49$

$cost(x)=f_\theta^2(x)$

$grad(cost(x)) = 2 * f_\theta(x) * grad(f_\theta(x)) = 2 * f_\theta(x) * 2\theta$

$\theta_n = \theta_{n-1} - lr * 2 * f_{\theta_{n-1}}(x) * 2 \theta_{n-1}$

算法实现

import random

def f(a):
    return a**2 - 49

def cost(a):
    return f(a)**2

def grad(a):
    return 2*f(a)*2*a

def solve(eps, lr, max_iterations=1000):
    a = random.random()  #产生随机值
    i = 1
    print('Initial value:', a)
    curr_cost = cost(a)  #计算损失值
    d_cost = abs(curr_cost)  #对损失值取绝对值
    while d_cost > eps and max_iterations > 0:
        a -= grad(a) * lr
        print('Iteration {}: {}'.format(i, a))
        max_iterations -= 1
        prev_cost, curr_cost = curr_cost, cost(a)
        d_cost = abs(curr_cost - prev_cost)
        i += 1
    return a if max_iterations > 0 else None

eps = 0.001
lr = 0.00001
print('Solution: ', solve(eps, lr))

eps表示表征预测值与准确值之间的可接受误差。
max_iterations表示最大训练次数

运行结果：

Iteration 995: 3.9199363201818906
Iteration 996: 3.9252100612710863
Iteration 997: 3.9304844015103737
Iteration 998: 3.9357593147627243
Iteration 999: 3.9410347748443146
Iteration 1000: 3.9463107555250176
Solution:  None

很明显，学习率太低，迭代了1000次也没用迭代到我们想要的结果。调大学习率lr=0.0001，然后

Iteration 303: 6.991042411421022
Iteration 304: 6.991392875178312
Iteration 305: 6.991729652432557
Iteration 306: 6.992053275734823
Iteration 307: 6.99236425706317
Solution:  6.99236425706317

这次迭代了307次，计算出了结果是6.99236425706317，误差在 $\pm0.001$ 内。

那么问题来了，既然调大学习率可以减少迭代次数，那么我们将学习率继续调大会怎么样？比如将学习率设定为0.01，那么我们看下运行结果

Iteration 995: 9.571797108320075
Iteration 996: -6.745934546769428
Iteration 997: -7.688305739056321
Iteration 998: -4.579141017036237
Iteration 999: -9.713542729358231
Iteration 1000: 7.907955353113504
Solution:  None

可以看到，预测值在反复横跳，就是不收敛，最终结果发散，无法获得准确的值。

学习率自适应

所以学习率在神经网络中是一个非常重要的参数，那么既然学习率过大会不收敛，学习率过高会增加训练次数，那么自适应的学习率就显得很重要了，一般简单的做法是，先设定一个较大值，训练n次之后，学习率减半，如此往复。或者是在n次cost的值都大于上一次之后，学习率减半。

代码如下：

# 如果连续patience次迭代cost都大于上一次迭代，那么学习率自动变小到decay倍
def solve(eps, lr, max_iterations=1000, patience=1, decay=0.5):
    a = random.random()
    i = 1
    print('Initial value:', a)
    print('Initial lr:', lr)
    curr_cost = cost(a)
    d_cost = abs(curr_cost)
    prev_costs = []
    while d_cost > eps and max_iterations > 0:
        a -= grad(a) * lr
        print('Iteration {}: {}'.format(i, a))
        max_iterations -= 1

        if len(prev_costs) > 0 and curr_cost < prev_costs[-1]:
            prev_costs.clear()
        elif len(prev_costs) == patience and curr_cost > prev_costs[-1]:
            lr *= decay
            print('New lr:', lr)
            prev_costs.clear()

        prev_cost, curr_cost = curr_cost, cost(a)
        d_cost = abs(curr_cost - prev_cost)
        prev_costs.insert(0, prev_cost)
        i += 1
    return a if max_iterations > 0 else None

结果就不放了，这个函数比较简单，也没有太好的优化效果。

Pytorch实现

然后是pytorch的实现方法，运用pytorch中的优化器和学习率适应算法，例如Adam等，还有自动求导等等，可以让代码更加简单。

import torch

class Model(torch.nn.Module):
    def __init__(self):
        super(Model, self).__init__()
        self.a = torch.nn.Parameter(torch.rand(1, requires_grad=True))

    def forward(self, x):
        return self.a**2 - 49


model = Model()
criterion = torch.nn.MSELoss(reduction='sum')
optimizer = torch.optim.SGD(model.parameters(), lr=5e-5)
scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, factor=0.5, patience=1)

for i in range(300):
    y = model(torch.rand(1))
    loss = criterion(y, torch.zeros(1))
    optimizer.zero_grad() #将导数置0
    loss.backward()  #反向传播求梯度
    optimizer.step()  #更新优化器参数
    print('Iteration {}: lr={}, a={}'.format(i+1, optimizer.param_groups[0]['lr'], model.state_dict()['a'].item()))
    if loss.item()<0.0001:  #当预测值在误差内提前打断学习
        break
    scheduler.step(loss)  #更新学习率参数

首先是torch.nn.Parameter(torch.rand(1, requires_grad=True))，可以把这个函数理解为类型转换函数，将一个不可训练的类型Tensor转换成可以训练的类型parameter并将这个parameter绑定到这个module里面(net.parameter()中就有这个绑定的parameter，所以在参数优化的时候可以进行优化的)，所以经过类型转换这个self.a变成了模型的一部分，成为了模型中根据训练可以改动的参数了。使用这个函数的目的也是想让某些变量在学习的过程中不断的修改其值以达到最优化。
torch.nn.MSELoss是均方损失类，它的公式是 $loss(x_i,y_i)=(x_i-y_i)^2$ ，它有三个类型的参数，分别是size_average,reduce以及reduction，以后再介绍，其中reduction="sum"表示计算输入值的加权和，并且返回一个标量。loss = criterion(y, torch.zeros(1))就是计算预测值的输出y与我们的期望值0之间的均方误差。

优化器

torch.optim是一个优化算法的集合包，可以很容易的优化我们的神经网络。要构造一个Optimizer，你必须给它一个包含参数（必须都是Variable对象）进行优化。然后，您可以指定optimizer的参数选项，比如学习率，权重衰减等。具体参考torch.optim中文文档。

Stochastic Gradient Descent (SGD)

SGD是最基础的优化方法，普通的训练方法, 需要重复不断的把整套数据放入神经网络NN中训练, 这样消耗的计算资源会很大。当我们使用SGD会把数据拆分后再分批不断放入 NN 中计算。每次使用批数据, 虽然不能反映整体数据的情况, 不过却很大程度上加速了 NN 的训练过程, 而且也不会丢失太多准确率。

Momentum

Momentum 传统的参数 W 的更新是把原始的 W 累加上一个负的学习率(learning rate) 乘以校正值 (dx)。此方法比较曲折。

我们把这个人从平地上放到了一个斜坡上, 只要他往下坡的方向走一点点, 由于向下的惯性, 他不自觉地就一直往下走, 走的弯路也变少了。这就是 Momentum 参数更新。

AdaGrad 优化学习率

AdaGrad 优化学习率，使得每一个参数更新都会有自己与众不同的学习率。与momentum类似，不过不是给喝醉酒的人安排另一个下坡, 而是给他一双不好走路的鞋子, 使得他一摇晃着走路就脚疼, 鞋子成为了走弯路的阻力, 逼着他往前直着走。

RMSProp

RMSProp 有了 momentum 的惯性原则 , 加上 adagrad 的对错误方向的阻力, 我们就能合并成这样. 让 RMSProp同时具备他们两种方法的优势. 不过细心的同学们肯定看出来了, 似乎在 RMSProp 中少了些什么. 原来是我们还没把 Momentum合并完全, RMSProp 还缺少了 momentum 中的这一部分。所以, 我们在 Adam 方法中补上了这种想法。

Adam

Adam 计算m 时有 momentum 下坡的属性, 计算 v 时有 adagrad 阻力的属性, 然后再更新参数时把 m 和 V 都考虑进去. 实验证明, 大多数时候, 使用 adam 都能又快又好的达到目标, 迅速收敛. 所以说, 在加速神经网络训练的时候, 一个下坡, 一双破鞋子, 功不可没。

代码如下，可以把上面的SGD优化器替换看神经网络运行的效果(LR和其他参数根据实际情况调节)：

# SGD 就是随机梯度下降
optimizer = torch.optim.SGD(net_SGD.parameters(), lr=LR)
# momentum 动量加速,在SGD函数里指定momentum的值即可
optimizer = torch.optim.SGD(net_Momentum.parameters(), lr=LR, momentum=0.8)
# RMSprop 指定参数alpha
optimizer = torch.optim.RMSprop(net_RMSprop.parameters(), lr=LR, alpha=0.8)
# Adam 参数betas=(0.9, 0.99)
optimizer = torch.optim.Adam(net_Adam.parameters(), lr=LR, betas=(0.8, 0.99))