神经网络基础系列之《优化器——原理及应用》

神经网络在每次迭代中，梯度下降根据⾃变量当前位置，沿着当前位置的梯度更新⾃变量。然而，如果⾃变量的 迭代⽅向仅仅取决于⾃变量当前位置，这可能会带来⼀些问题。比如下图：

fig1

给定学习率，梯度下降迭代⾃变量时会使⾃变量在竖直⽅向⽐在⽔平⽅向移动幅度更⼤。那么，我们需要⼀个较小的学习率从而避免⾃变量在竖直⽅向上越过⽬标函数最优解。然而，这会造成⾃变量在⽔平⽅向上朝最优解移动变慢。

为了解决梯度更新的问题，提出了很多优化器算法，下面一一介绍。

给定学习率，梯度下降迭代⾃变量时会使⾃变量在竖直⽅向⽐在⽔平⽅向移动幅度更⼤。那么，我们需要⼀个较小的学习率从而避免⾃变量在竖直⽅向上越过⽬标函数最优解。然而，这会造成⾃变量在⽔平⽅向上朝最优解移动变慢。

为了解决梯度更新的问题，提出了很多优化器算法，下面一一介绍。

SGD随机梯度下降

就是对minibtch做梯度下降。其优缺点如下：

  • 大部分时候你向着全局最小值靠近，有时候你会远离最小值，因为那个样本恰好给你指的方向不对，因此随机梯度下降法是有很多噪声的，平均来看，它最终会靠近最小值，不过有时候也会方向错误，因为随机梯度下降法永远不会收敛，而是会一直在最小值附近波动。一次性只处理了一个训练样本，这样效率过于低下。

  • 实践中最好选择不大不小的 mini-batch，得到了大量向量化，效率高，收敛快。

下面以LR算法为例，说明SGD的计算过程。

  • 模拟数据并归一化

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    %matplotlib inline
    x = [30,35,37,59,70,76,88,100]
    y = [1100,1423,1377,1800,2304,2588,3495,4839]

    x_max = max(x)
    x_min = min(x)
    y_max = max(y)
    y_min = min(y)

    for i in range(0,len(x)):
        x[i] = (x[i] - x_min)/(x_max - x_min)
        y[i] = (y[i] - y_min)/(y_max - y_min)
        
    print (x)
    print (y)

  • 数据等高线

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    import numpy as np

    def calc_loss(a,b,x,y):
        tmp = y - (a * x + b)
        tmp = tmp ** 2  # 对矩阵内的每一个元素平方
        SSE = sum(tmp) / (2 * len(x))
        return SSE

    def draw_hill(x,y):
        a = np.linspace(-20,20,100)
        b = np.linspace(-20,20,100)
        x = np.array(x)
        y = np.array(y)

        allSSE = np.zeros(shape=(len(a),len(b)))
        for ai in range(0,len(a)):
            for bi in range(0,len(b)):
                a0 = a[ai]
                b0 = b[bi]
                SSE = calc_loss(a=a0,b=b0,x=x,y=y)
                allSSE[ai][bi] = SSE

        a,b = np.meshgrid(a, b)
        return [a,b,allSSE]

    [ha,hb,hallSSE] = draw_hill(x,y)
    hallSSE = hallSSE.T# 重要，将所有的losses做一个转置。原因是矩阵是以左上角至右下角顺序排列元素，而绘图是以左下角为原点。

    print (ha.shape)
    print (hb.shape)
    print (hallSSE.shape)

  • 绘制学习率的曲线

    1
    2
    3
    4
    import matplotlib.pyplot as plt
    rate = 0.1 # learning rate
    fig = plt.figure(1, figsize=(12, 8))
    fig.suptitle('learning rate: %.2f method:momentum'%(rate), fontsize=15)

  • 绘制曲面

    1
    2
    3
    4
    from mpl_toolkits.mplot3d import Axes3D
    ax = fig.add_subplot(2, 2, 1, projection='3d')
    ax.set_top_view()
    ax.plot_surface(ha, hb, hallSSE, rstride=2, cstride=2, cmap='rainbow')

  • 绘制等高线图

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    plt.subplot(2,2,2)
    ta = np.linspace(-20, 20, 100)
    tb = np.linspace(-20, 20, 100)
    plt.contourf(ha,hb,hallSSE,15,alpha=0.5,cmap=plt.cm.hot)
    C = plt.contour(ha,hb,hallSSE,15,colors='black')
    plt.clabel(C,inline=True)
    plt.xlabel('a')
    plt.ylabel('b')

    plt.ion() # iteration on

  • 求loss的梯度

    1
    2
    3
    4
    5
    def da(y,y_p,x):
        return (y-y_p)*(-x)

    def db(y,y_p):
        return (y-y_p)*(-1)

  • 进行迭代训练

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    import random

    def shuffle_data(x,y):
        # 随机打乱x，y的数据，并且保持x和y一一对应
        seed = random.random()
        random.seed(seed)
        random.shuffle(x)
        random.seed(seed)
        random.shuffle(y)

    def get_batch_data(x,y,batch=3):
        shuffle_data(x,y)
        x_new = x[0:batch]
        y_new = y[0:batch]
        return [x_new,y_new]

    a = 10.0
    b = -20.0
    rate = 0.2

    all_loss = []
    all_step = []
    last_a = a
    last_b = b
    for step in range(1,200):
        loss = 0
        all_da = 0
        all_db = 0
        shuffle_data(x,y)
        [x_new,y_new] = get_batch_data(x,y,batch=4)
        for i in range(0,len(x_new)):
            y_p = a*x_new[i] + b
            loss = loss + (y_new[i] - y_p)*(y_new[i] - y_p)/2
            all_da = all_da + da(y_new[i],y_p,x_new[i])
            all_db = all_db + db(y_new[i],y_p)
        #loss_ = calc_loss(a = a,b=b,x=np.array(x),y=np.array(y))
        loss = loss/len(x_new)

        # 绘制图1中的loss点
        ax.scatter(a, b, loss, color='black')
        # 绘制图2中的loss点
        plt.subplot(2,2,2)
        plt.scatter(a,b,s=5,color='blue')
        plt.plot([last_a,a],[last_b,b],color='aqua')
        # 绘制图3中的回归直线
        plt.subplot(2, 2, 3)
        plt.plot(x, y)
        plt.plot(x, y, 'o')
        x_ = np.linspace(0, 1, 2)
        y_draw = a * x_ + b
        plt.plot(x_, y_draw)
        # 绘制图4的loss更新曲线
        all_loss.append(loss)
        all_step.append(step)
        plt.subplot(2,2,4)
        plt.plot(all_step,all_loss,color='orange')
        plt.xlabel("step")
        plt.ylabel("loss")


        # print('a = %.3f,b = %.3f' % (a,b))
        last_a = a
        last_b = b
        a = a - rate*all_da
        b = b - rate*all_db

        if step%1 == 0:
            print("step: ", step, " loss: ", loss)
            plt.show()
            plt.pause(0.01)
    plt.show()

Momentum 动量算法

如果把梯度下降法想象成一个小球从山坡到山谷的过程，那么前面几篇文章的小球是这样移动的：从A点开始，计算当前A点的坡度，沿着坡度最大的方向走一段路，停下到B。在B点再看一看周围坡度最大的地方，沿着这个坡度方向走一段路，再停下。确切的来说，这并不像一个球，更像是一个正在下山的盲人，每走一步都要停下来，用拐杖来来探探四周的路，再走一步停下来，周而复始，直到走到山谷。而一个真正的小球要比这聪明多了，从A点滚动到B点的时候，小球带有一定的初速度，在当前初速度下继续加速下降，小球会越滚越快，更快的奔向谷底。momentum 动量法就是模拟这一过程来加速神经网络的优化的。可以用下图说明：

fig2

根据图示，其计算过程如下：

  • A为起始点，首先计算A点的梯度∇a，然后下降到B点：$\Theta_{new}=\Theta−\alpha∇a$，其中$\Theta$为参数，$\alpha$为学习率。

  • 到了B点需要加上A点的梯度，这里梯度需要有一个衰减值$\gamma$，推荐取0.9。这样的做法可以让早期的梯度对当前梯度的影响越来越小，如果没有衰减值，模型往往会震荡难以收敛，甚至发散。所以B点的参数更新公式是这样的：$v_t=\gamma v_{t-1} + \alpha∇b$，$\Theta_{new}=\Theta - v_t$。其中$v_{t-1}$表示之前所有步骤所累积的动量和。

Momentum算法的特点如下：

  • 下降初期时，使用上一次参数更新，下降方向一致，乘上较大的$\gamma$能够进行很好的加速

  • 下降中后期时，在局部最小值来回震荡的时候，$∇b->0$，$\gamma$使得更新幅度增大，跳出陷阱

  • 在梯度改变方向的时候，$\gamma$能够减少更新 总而言之，momentum项能够在相关方向加速SGD，抑制振荡，从而加快收敛

下面以LR算法为例，说明动量的计算过程：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
def da(y,y_p,x):
    return (y-y_p)*(-x)

def db(y,y_p):
    return (y-y_p)*(-1)

a = 10.0
b = -20.0
all_loss = []
all_step = []
last_a = a
last_b = b
va = 0
vb = 0
gamma = 0.9
for step in range(1,100):
    loss = 0
    all_da = 0
    all_db = 0
    for i in range(0,len(x)):
        y_p = a*x[i] + b
        loss = loss + (y[i] - y_p)*(y[i] - y_p)/2  
        # loss = (y-ax[i]-b)(y-ax[i]-b)/2=(y^2 + (ax[i])^2 + b^2 - 2yax[i] - 2yb + 2ax[i]b/2
        # 对loss中的参数a求偏导：(2ax[i]^2 - 2yx[i] + 2x[i]b)/2 = x[i](ax[i]+b-y) = x[i](y[i]-y) = -x[i](y - y[i]) 
        all_da = all_da + da(y[i],y_p,x[i])
        # 对loss中的参数b求偏导：(2b - 2y + 2ax[i])/2 = (ax[i] + b - y) = (y - y[i]) * (-1)
        all_db = all_db + db(y[i],y_p)
    #loss_ = calc_loss(a = a,b=b,x=np.array(x),y=np.array(y))
    loss = loss/len(x)

    # 绘制图1中的loss点
    ax.scatter(a, b, loss, color='black')
    # 绘制图2中的loss点
    plt.subplot(2,2,2)
    plt.scatter(a,b,s=5,color='blue')
    plt.plot([last_a,a],[last_b,b],color='aqua')
    # 绘制图3中的回归直线
    plt.subplot(2, 2, 3)
    plt.plot(x, y)
    plt.plot(x, y, 'o')
    x_ = np.linspace(0, 1, 2)
    y_draw = a * x_ + b
    plt.plot(x_, y_draw)
    # 绘制图4的loss更新曲线
    all_loss.append(loss)
    all_step.append(step)
    plt.subplot(2,2,4)
    plt.plot(all_step,all_loss,color='orange')
    plt.xlabel("step")
    plt.ylabel("loss")


    # print('a = %.3f,b = %.3f' % (a,b))
    last_a = a
    last_b = b
    va = gamma * va+ rate*all_da
    vb = gamma * vb+ rate*all_db
    a = a - va
    b = b - vb

    if step%1 == 0:
        print("step: ", step, " loss: ", loss)
        plt.show()
        plt.pause(0.01)
plt.show()
# plt.pause(9999)

总结：优化算法中，⽬标函数⾃变量的每⼀个元素在相同时间步都使⽤同⼀个学习率来⾃我迭代。在“动量法”⾥我们看到当x1和x2的梯度值有较⼤差别时，需要选择⾜够小的学习率使得⾃变量在梯度值较⼤的维度上不发散。但这样会导致⾃变量在梯度值较小的维度上迭代过慢。动量法依赖指数加权移动平均使得⾃变量的更新⽅向更加⼀致，从而降低发散的可能。

AdaGrad算法

AdaGrad思路基本是借鉴L2 Regularizer，不过此时调节的不是𝑊，而是𝐺𝑟𝑎𝑑𝑖𝑒𝑛𝑡。AdaGrad算法根据⾃变量在每个维度的梯度值的⼤小来调整各个维度上的学习率，从而避免统⼀的学习率难以适应所有维度的问题。

这个算法的优点是可以对低频的参数做较大的更新，对高频的做较小的更新，也因此，对于稀疏的数据它的表现很好，很好地提高了 SGD 的鲁棒性，例如识别 Youtube 视频里面的猫，训练 GloVe word embeddings，因为它们都是需要在低频的特征上有更大的更新。

它的缺点是当学习率在迭代早期降得较快且当前解依然不佳时，AdaGrad算法在迭代后期由于学习率过小，可能较难找到⼀个有⽤的解。具体计算过程如下：

  • $s_t = s_{t-1} + g_t \odot g_t$

  • $x_t=x_{t-1}-\frac{\eta }{s_t + \varepsilon } \odot g_t$

在时间步0，AdaGrad将$s_0$中每个元素初始化为0。在时间步t，⾸先将小批量随机梯度$g_t$按元素平⽅后累加到变量$s_t$，接着，我们将⽬标函数⾃变量中每个元素的学习率通过按元素运算重新调整⼀下。其中$\eta$是学习率，$\varepsilon$是为了维持数值稳定性而添加的常数，如10的-6次方。这⾥开⽅、除法和乘法的运算都是按元素运算的。这些按元素运算使得⽬标函数⾃变量中每个元素都分别拥有⾃⼰的学习率。一般$\varepsilon$选取0.01

需要强调的是，小批量随机梯度按元素平⽅的累加变量$s_t$出现在学习率的分⺟项中。因此：

  • 如果⽬标函数有关⾃变量中某个元素的偏导数⼀直都较⼤，那么该元素的学习率将下降较快；也就是说，在训练前期，梯度较小，使得Regularizer项很大，放大梯度。[激励阶段]

  • 反之，如果⽬标函数有关⾃变量中某个元素的偏导数⼀直都较小，那么该元素的学习率将下降较慢；也就是说，训练后期，梯度较大，使得Regularizer项很小，缩小梯度。[惩罚阶段]

可以这样理解：AdaGrad过程，是一个递推过程，次从𝜏=1，推到𝜏=𝑡，把沿路的𝐺𝑟𝑎𝑑𝑖𝑒𝑛𝑡的平方根，作为Regularizer。由于Regularizer是专门针对Gradient的，所以有利于解决Gradient Vanish/Expoloding问题。

下图为AdaGrad优化的一个过程：

fig3

AdaGrad的缺点是：

  • 由公式可以看出，仍依赖于人工设置一个全局学习率。$\eta$设置过大的话，会使regularizer过于敏感，对梯度的调节太大

  • 中后期，分母上梯度平方的累加将会越来越大，使$g_t->0$，使得训练提前结束

下面继续以LR为例，说明其计算过程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
a = 10.0
b = -20.0
all_loss = []
all_step = []
last_a = a
last_b = b
n = np.array([0,0])
epsilon = 1e-8
for step in range(1,500):
    loss = 0
    all_da = 0
    all_db = 0
    for i in range(0,len(x)):
        y_p = a*x[i] + b
        loss = loss + (y[i] - y_p)*(y[i] - y_p)/2
        all_da = all_da + da(y[i],y_p,x[i])
        all_db = all_db + db(y[i],y_p)
    #loss_ = calc_loss(a = a,b=b,x=np.array(x),y=np.array(y))
    loss = loss/len(x)

    # 绘制图1中的loss点
    ax.scatter(a, b, loss, color='black')
    # 绘制图2中的loss点
    plt.subplot(2,2,2)
    plt.scatter(a,b,s=5,color='blue')
    plt.plot([last_a,a],[last_b,b],color='aqua')
    # 绘制图3中的回归直线
    plt.subplot(2, 2, 3)
    plt.plot(x, y)
    plt.plot(x, y, 'o')
    x_ = np.linspace(0, 1, 2)
    y_draw = a * x_ + b
    plt.plot(x_, y_draw)
    # 绘制图4的loss更新曲线
    all_loss.append(loss)
    all_step.append(step)
    plt.subplot(2,2,4)
    plt.plot(all_step,all_loss,color='orange')
    plt.xlabel("step")
    plt.ylabel("loss")

    # print('a = %.3f,b = %.3f' % (a,b))
    last_a = a
    last_b = b

    #-- 参数更新
    n[0] = n[0]+np.square(all_da)
    n[1] = n[1]+np.square(all_db)
    rate_new = rate/(np.sqrt(n + epsilon))
    print('rate_new a:',rate_new[0],' b:',rate_new[1])
    a = a - (rate/(np.sqrt(n[0] + epsilon)))*all_da
    b = b - (rate/(np.sqrt(n[1] + epsilon)))*all_db

    if step%1 == 0:
        print("step: ", step, " loss: ", loss)
        plt.show()
        plt.pause(0.01)
plt.show()

RMSProp算法

当学习率在迭代早期降得较快且当前解依然不佳时，AdaGrad算法在迭代后期由于学习率过小，可能较难找到⼀个有⽤的解。为了解决这⼀问题，RMSProp算法对AdaGrad算法做了⼀点小小的修改：

  • $s_t = \lambda s_{t-1} + (1-\lambda)g_t \odot g_t$

  • $x_t=x_{t-1}-\frac{\eta }{s_t + \varepsilon } \odot g_t$

因为RMSProp算法的状态变量st是对平⽅项$g_t \odot g_t$的指数加权移动平均，所以可以看作是最近$\frac{1}{1-\lambda}$个时间步的小批量随机梯度平⽅项的加权平均。如此⼀来，⾃变量每个元素的学习率在迭代过程中就不再⼀直降低（或不变）。

下图为RMSProp优化的一个过程：

fig4

RMSProp有如下特点：

  • 训练初中期，加速效果不错，很快

  • 训练后期，反复在局部最小值附近抖动（下面实验中，如果step设置为500，可能出现loss非常大的情况）

  • 其实RMSprop依然依赖于全局学习率

  • RMSprop算是Adagrad的一种发展，和Adadelta的变体，效果趋于二者之间

  • 适合处理非平稳目标 - 对于RNN效果很好

  • RMSProp利用了二阶信息做了Gradient优化，在BatchNorm之后，对其需求不是很大。

下面继续以LR为例，说明其计算过程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
a = 10.0
b = -20.0
all_loss = []
all_step = []
last_a = a
last_b = b
lambdas = 0.9
n = np.array([0,0])
epsilon = 1e-8
for step in range(1,200):
    loss = 0
    all_da = 0
    all_db = 0
    for i in range(0,len(x)):
        y_p = a*x[i] + b
        loss = loss + (y[i] - y_p)*(y[i] - y_p)/2
        all_da = all_da + da(y[i],y_p,x[i])
        all_db = all_db + db(y[i],y_p)
    #loss_ = calc_loss(a = a,b=b,x=np.array(x),y=np.array(y))
    loss = loss/len(x)

    # 绘制图1中的loss点
    ax.scatter(a, b, loss, color='black')
    # 绘制图2中的loss点
    plt.subplot(2,2,2)
    plt.scatter(a,b,s=5,color='blue')
    plt.plot([last_a,a],[last_b,b],color='aqua')
    # 绘制图3中的回归直线
    plt.subplot(2, 2, 3)
    plt.plot(x, y)
    plt.plot(x, y, 'o')
    x_ = np.linspace(0, 1, 2)
    y_draw = a * x_ + b
    plt.plot(x_, y_draw)
    # 绘制图4的loss更新曲线
    all_loss.append(loss)
    all_step.append(step)
    plt.subplot(2,2,4)
    plt.plot(all_step,all_loss,color='orange')
    plt.xlabel("step")
    plt.ylabel("loss")

    # print('a = %.3f,b = %.3f' % (a,b))
    last_a = a
    last_b = b

    #-- 参数更新
    n[0] = lambdas * n[0] + (1 - lambdas) * np.square(all_da)
    n[1] = lambdas * n[1] + (1 - lambdas) * np.square(all_db)
    rate_new = rate/(np.sqrt(n + epsilon))
    print('rate_new a:',rate_new[0],' b:',rate_new[1])
    a = a - (rate/(np.sqrt(n[0] + epsilon)))*all_da
    b = b - (rate/(np.sqrt(n[1] + epsilon)))*all_db

    if step%1 == 0:
        print("step: ", step, " loss: ", loss)
        plt.show()
        plt.pause(0.01)
plt.show()

AdaDelta算法

AdaDelta基本思想是用一阶的方法，近似模拟二阶牛顿法。除了RMSProp算法以外，另⼀个常⽤优化算法AdaDelta算法也针对AdaGrad算法在迭代后期可能较难找到有⽤解的问题做了改进。有意思的是，AdaDelta算法没有学习率这⼀超参数。其计算过程如下：

  • $s_t = \lambda s_{t-1} + (1-\lambda)g_t \odot g_t$

  • $g’_t=\sqrt{\frac{\Delta x_t + \varepsilon }{s_t + \varepsilon}} \odot g_t$

  • $\Delta x_t=px_{t-1}-(1-p)g’_t \odot g’_t$

  • $x_t = x_{t-1} - g’_t$

AdaDelta算法的第一步和RMSProp算法⼀样，与RMSProp算法不同的是，AdaDelta算法还维护⼀个额外的状态变量$\Delta x_t$来计算自变量的变化量（按$g’_t$元素平⽅的指数加权移动平均），其元素同样在时间步0时被初始化为0。

可以这样理解：相比于同样是基于Gradient的Regularizer，不过只取最近的w个状态，这样不会让梯度被惩罚至0。

与RMSProp相比，AdaDelta算法没有学习率超参数，它通过使用有关自变量更新量平方的指数加权移动平均的项来替代RMSProp算法中的学习率。

AdaDelta的特点：

  • 从多个数据集情况来看，AdaDelta在训练初期和中期，具有非常不错的加速效果。

  • 但是到训练后期，进入局部最小值雷区之后，AdaDelta就会反复在局部最小值附近抖动。主要体现在验证集错误率上，脱离不了局部最小值吸引盆。这时候，切换成动量SGD，如果把学习率降低一个量级，就会发现验证集正确率有2%~5%的提升。[注]：使用Batch Norm之后，这样从AdaDelta切到SGD会导致数值体系崩溃，原因未知。

下面继续以LR为例，说明其计算过程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
a = 10.0
b = -20.0
all_loss = []
all_step = []
last_a = a
last_b = b
n = np.array([0,0])
theta = np.array([0,0]).astype(np.float32) # 每一次a,b迭代的更新值

apple = np.array([0,0]).astype(np.float32)
pear = np.array([0,0]).astype(np.float32)
# 迭代
for step in range(1,500):
    loss = 0
    all_da = 0
    all_db = 0
    all_d = np.array([0,0]).astype(np.float32)
    for i in range(0,len(x)):
        y_p = a*x[i] + b
        loss = loss + (y[i] - y_p)*(y[i] - y_p)/2
        all_da = all_da + da(y[i],y_p,x[i])
        all_db = all_db + db(y[i],y_p)
    #loss_ = calc_loss(a = a,b=b,x=np.array(x),y=np.array(y))
    all_d = np.array([all_da,all_db])
    loss = loss/len(x)

    # 绘制图1中的loss点
    ax.scatter(a, b, loss, color='black')
    # 绘制图2中的loss点
    plt.subplot(2,2,2)
    plt.scatter(a,b,s=5,color='blue')
    plt.plot([last_a,a],[last_b,b],color='aqua')
    # 绘制图3中的回归直线
    plt.subplot(2, 2, 3)
    plt.plot(x, y)
    plt.plot(x, y, 'o')
    x_ = np.linspace(0, 1, 2)
    y_draw = a * x_ + b
    plt.plot(x_, y_draw)
    # 绘制图4的loss更新曲线
    all_loss.append(loss)
    all_step.append(step)
    plt.subplot(2,2,4)
    plt.plot(all_step,all_loss,color='orange')
    plt.xlabel("step")
    plt.ylabel("loss")

    # print('a = %.3f,b = %.3f' % (a,b))
    last_a = a
    last_b = b

    #-- 参数更新
    apple = gamma*apple + (1-gamma)*(all_d**2) # apple with all_d of this step
#     rms_apple = np.sqrt(apple + epsilon)
    rms_apple = apple + epsilon

    pear = gamma*pear + (1-gamma)*(theta**2) # pear with theta of last step
#     rms_pear = np.sqrt(pear + epsilon)

#     theta = -(rms_pear/rms_apple)*all_d
    theta = -np.sqrt(rms_pear/rms_apple) * all_d
    [a,b] = [a,b] + theta

    if step%1 == 0:
        print("step: ", step, " loss: ", loss,"rms_pear: ",rms_pear," rms_apple",rms_apple)
        plt.show()
        plt.pause(0.01)
plt.show()

Adam算法

Adam算法在RMSProp算法基础上对小批量随机梯度也做了指数加权移动平均。

  • $v_t=\beta_1v_{t-1} + (1-\beta_1)g_t$

    • Adam算法使⽤了动量变量$v_t$和RMSProp算法中小批量随机梯度按元素平⽅的指数加权移动平均变量$s_t$，并在时间步0将它们中每个元素初始化为0。给定超参数0 ≤ $\beta_1$ < 1（算法作者建议设为0.9），时间步t的动量变量$v_t$即小批量随机梯度$g_t$的指数加权移动平均

  • $s_t=\beta_2s_{t-1}+(1-\beta_2)g_t\odot g_t$

    • 和RMSProp算法中⼀样，给定超参数0 ≤ $\beta_2$ < 1（算法作者建议设为0.999），将小批量随机梯度按元素平⽅后的项$g_t\odot g_t$做指数加权移动平均得到$s_t$

  • $v’_t=\frac{v_t}{1-\beta^t_1}$、$s’_t=\frac{s_t}{1-\beta^t_2}$

    • 因为当t较小时，过去各时间步小批量随机梯度权值之和会较小。例如，当$\beta_1=0.9$时，$v_1=0.1g_1$。为了消除这样的影响，对于任意时间步t，我们可以将$v_t$再除以${1-\beta^t_1}$，从而使过去各时间步小批量随机梯度权值之和为1。这也叫作偏差修正。在Adam算法中，我们对变量$v_t$和$s_t$均作偏差修正得到

  • $g’_t=\frac{\eta v’_t}{\sqrt{s’_t}+\varepsilon }$

    • 接下来，Adam算法使⽤以上偏差修正后的变量$v_t$和$s_t$，将模型参数中每个元素的学习率通过按元素运算重新调整得到（5）。其中$\eta$是学习率，$\varepsilon$是为了维持数值稳定性而添加的常数，如10的-8次方。和AdaGrad算法、RMSProp算法以及AdaDelta算法⼀样，⽬标函数⾃变量中每个元素都分别拥有⾃⼰的学习率

  • $x_t=x_{t-1}-g’_t$

    • 使用$g’_t$迭代自变量

Adam算法有如下特点：

  • 除了像Adadelta和RMSprop一样存储了过去梯度的平方$s_t$的指数衰减平均值 ，也像momentum一样保持了过去梯度$v_t$的指数衰减平均值
  • 如果$v_t$和$s_t$被初始化为0向量，那它们就会向0偏置，所以做了偏差校正，通过计算偏差校正后的$v_t$和$s_t$来抵消这些偏差

下面继续以LR为例，说明其计算过程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
a = 10.0
b = -20.0
all_loss = []
all_step = []
last_a = a
last_b = b
m = 0.0
v = 0.0
theta = np.array([0,0]).astype(np.float32)

rate = 0.001
beta1 = 0.9
beta2 = 0.999
epsilon = 1e-8

for step in range(1,500):
    loss = 0
    all_da = 0
    all_db = 0
    for i in range(0,len(x)):
        y_p = a*x[i] + b
        loss = loss + (y[i] - y_p)*(y[i] - y_p)/2
        all_da = all_da + da(y[i],y_p,x[i])
        all_db = all_db + db(y[i],y_p)
    loss = loss/len(x)
    all_d = np.array([all_da,all_db]).astype(np.float32)
    # 绘制图1中的loss点
    ax.scatter(a, b, loss, color='black')
    # 绘制图2中的loss点
    plt.subplot(2,2,2)
    plt.scatter(a,b,s=5,color='blue')
    plt.plot([last_a,a],[last_b,b],color='aqua')
    # 绘制图3中的回归直线
    plt.subplot(2, 2, 3)
    plt.plot(x, y)
    plt.plot(x, y, 'o')
    x_ = np.linspace(0, 1, 2)
    y_draw = a * x_ + b
    plt.plot(x_, y_draw)
    # 绘制图4的loss更新曲线
    all_loss.append(loss)
    all_step.append(step)
    plt.subplot(2,2,4)
    plt.plot(all_step,all_loss,color='orange')
    plt.xlabel("step")
    plt.ylabel("loss")


    # print('a = %.3f,b = %.3f' % (a,b))
    last_a = a
    last_b = b

    m = beta1*m + (1-beta1)*all_d
    v = beta2*v + (1-beta2)*(all_d**2)

    m_ = m/(1 - beta1)
    v_ = v/(1 - beta2)


    theta = -(rate*m_/(np.sqrt(v_) + epsilon))


    [a,b] = [a,b] + theta


    if step%1 == 0:
        print("step: ", step, " loss: ", loss)
        plt.show()
        plt.pause(0.01)
plt.show()

如何选择优化算法

  • 如果数据是稀疏的，就用自适用方法，即 Adagrad, Adadelta, RMSprop, Adam。

  • RMSprop, Adadelta, Adam 在很多情况下的效果是相似的。

  • Adam 就是在 RMSprop 的基础上加了 bias-correction 和 momentum，

  • 随着梯度变的稀疏，Adam 比 RMSprop 效果会好。

  • 整体来讲，Adam 是最好的选择。

  • 很多论文里都会用 SGD，没有 momentum 等。SGD 虽然能达到极小值，但是比其它算法用的时间长，而且可能会被困在鞍点。

  • 如果需要更快的收敛，或者是训练更深更复杂的神经网络，需要用一种自适应的算法。

keras中优化器的用法

以全连接网络做分类任务为例子，说明keras中optimizer的用法

生成数据

1
2
3
4
5
6
X = np.linspace(-1, 1, 200) #在返回（-1, 1）范围内的等差序列
np.random.shuffle(X)    # 打乱顺序
Y = 0.5 * X + 2 + np.random.normal(0, 0.05, (200, )) #生成Y并添加噪声

X_train, Y_train = X[:160], Y[:160]     # 前160组数据为训练数据集
X_test, Y_test = X[160:], Y[160:]      #后40组数据为测试数据集

构建模型

1
2
3
4
5
6
from keras import optimizers
from keras.layers import Dense, Activation
from keras.models import Sequential

model = Sequential()
model.add(Dense(input_dim=1, units=1))

使用SGD

1
2
3
4
5
6
# clipvalue=0.5表示保留(-0.5, 0.5)之间的梯度，其他的需要做梯度裁剪
# momentum=0.9使用了动量平滑，可以设置为0.0表示纯SGD，nesterov=True表示使用Nesterov动量
sgd = optimizers.SGD(lr=0.01, decay=1e-6, momentum=0.9, nesterov=True, clipvalue=0.5)
model.compile(loss='mean_squared_error', optimizer=sgd)

model.summary()

训练并测试

1
2
3
4
5
6
7
8
9
10
11
print('Training -----------')
for step in range(501):
    cost = model.train_on_batch(X_train, Y_train)
    if step % 50 == 0:
        print("After %d trainings, the cost: %f" % (step, cost))
        
print('\nTesting ------------')
cost = model.evaluate(X_test, Y_test, batch_size=40)
print('test cost:', cost)
W, b = model.layers[0].get_weights()
print('Weights=', W, '\nbiases=', b)

使用RMSprop

1
rmsprop = optimizers.RMSprop(lr=0.001, rho=0.9, epsilon=1e-06)

使用Adagrad

1
adagrad = optimizers.Adagrad(lr=0.01, epsilon=1e-06)

使用Adadelta

1
adadelta = optimizers.Adadelta(lr=1.0, rho=0.95, epsilon=1e-06)

使用Adam

1
adam = optimizers.Adam(lr=0.001, beta_1=0.9, beta_2=0.999, epsilon=1e-08)

BERT中使用的优化器解析

bert使用的optimizer叫做AdamWeightDecayOptimizer，先放代码为敬：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
import tensorflow as tf

def create_optimizer(loss, init_lr, num_train_steps, num_warmup_steps, use_tpu):
  """Creates an optimizer training op."""
  global_step = tf.train.get_or_create_global_step()

  learning_rate = tf.constant(value=init_lr, shape=[], dtype=tf.float32)

  # Implements linear decay of the learning rate.
  learning_rate = tf.train.polynomial_decay(
      learning_rate,
      global_step,
      num_train_steps,
      end_learning_rate=0.0,
      power=1.0,
      cycle=False)

  # Implements linear warmup. I.e., if global_step < num_warmup_steps, the
  # learning rate will be `global_step/num_warmup_steps * init_lr`.
  if num_warmup_steps:
    global_steps_int = tf.cast(global_step, tf.int32)
    warmup_steps_int = tf.constant(num_warmup_steps, dtype=tf.int32)

    global_steps_float = tf.cast(global_steps_int, tf.float32)
    warmup_steps_float = tf.cast(warmup_steps_int, tf.float32)

    warmup_percent_done = global_steps_float / warmup_steps_float
    warmup_learning_rate = init_lr * warmup_percent_done

    is_warmup = tf.cast(global_steps_int < warmup_steps_int, tf.float32)
    learning_rate = (
        (1.0 - is_warmup) * learning_rate + is_warmup * warmup_learning_rate)

  # It is recommended that you use this optimizer for fine tuning, since this
  # is how the model was trained (note that the Adam m/v variables are NOT
  # loaded from init_checkpoint.)
  optimizer = AdamWeightDecayOptimizer(
      learning_rate=learning_rate,
      weight_decay_rate=0.01,
      beta_1=0.9,
      beta_2=0.999,
      epsilon=1e-6,
      exclude_from_weight_decay=["LayerNorm", "layer_norm", "bias"])

  if use_tpu:
    optimizer = tf.contrib.tpu.CrossShardOptimizer(optimizer)

  tvars = tf.trainable_variables()
  grads = tf.gradients(loss, tvars)

  # This is how the model was pre-trained.
  (grads, _) = tf.clip_by_global_norm(grads, clip_norm=1.0)

  train_op = optimizer.apply_gradients(
      zip(grads, tvars), global_step=global_step)

  # Normally the global step update is done inside of `apply_gradients`.
  # However, `AdamWeightDecayOptimizer` doesn't do this. But if you use
  # a different optimizer, you should probably take this line out.
  new_global_step = global_step + 1
  train_op = tf.group(train_op, [global_step.assign(new_global_step)])
  return train_op

class AdamWeightDecayOptimizer(tf.train.Optimizer):
  """A basic Adam optimizer that includes "correct" L2 weight decay."""

  def __init__(self,
               learning_rate,
               weight_decay_rate=0.0,
               beta_1=0.9,
               beta_2=0.999,
               epsilon=1e-6,
               exclude_from_weight_decay=None,
               name="AdamWeightDecayOptimizer"):
    """Constructs a AdamWeightDecayOptimizer."""
    super(AdamWeightDecayOptimizer, self).__init__(False, name)

    self.learning_rate = learning_rate
    self.weight_decay_rate = weight_decay_rate
    self.beta_1 = beta_1
    self.beta_2 = beta_2
    self.epsilon = epsilon
    self.exclude_from_weight_decay = exclude_from_weight_decay

  def apply_gradients(self, grads_and_vars, global_step=None, name=None):
    """See base class."""
    assignments = []
    for (grad, param) in grads_and_vars:
      if grad is None or param is None:
        continue

      param_name = self._get_variable_name(param.name)

      m = tf.get_variable(
          name=param_name + "/adam_m",
          shape=param.shape.as_list(),
          dtype=tf.float32,
          trainable=False,
          initializer=tf.zeros_initializer())
      v = tf.get_variable(
          name=param_name + "/adam_v",
          shape=param.shape.as_list(),
          dtype=tf.float32,
          trainable=False,
          initializer=tf.zeros_initializer())

      # Standard Adam update.
      next_m = (
          tf.multiply(self.beta_1, m) + tf.multiply(1.0 - self.beta_1, grad))
      next_v = (
          tf.multiply(self.beta_2, v) + tf.multiply(1.0 - self.beta_2,
                                                    tf.square(grad)))

      update = next_m / (tf.sqrt(next_v) + self.epsilon)

      # Just adding the square of the weights to the loss function is *not*
      # the correct way of using L2 regularization/weight decay with Adam,
      # since that will interact with the m and v parameters in strange ways.
      #
      # Instead we want ot decay the weights in a manner that doesn't interact
      # with the m/v parameters. This is equivalent to adding the square
      # of the weights to the loss with plain (non-momentum) SGD.
      if self._do_use_weight_decay(param_name):
        update += self.weight_decay_rate * param

      update_with_lr = self.learning_rate * update

      next_param = param - update_with_lr

      assignments.extend(
          [param.assign(next_param),
           m.assign(next_m),
           v.assign(next_v)])
    return tf.group(*assignments, name=name)

  def _do_use_weight_decay(self, param_name):
    """Whether to use L2 weight decay for `param_name`."""
    if not self.weight_decay_rate:
      return False
    if self.exclude_from_weight_decay:
      for r in self.exclude_from_weight_decay:
        if re.search(r, param_name) is not None:
          return False
    return True

  def _get_variable_name(self, param_name):
    """Get the variable name from the tensor name."""
    m = re.match("^(.*):\\d+$", param_name)
    if m is not None:
      param_name = m.group(1)
    return param_name

可以看到里面使用了warmup+Adam+weight_dacay，现在一一进行说明。

（1）warmup机制在初期使用warmup_step阶段使用较小的学习率（随着warmup_step增大逐步增大到init_lr）

（2）Weight decay是在每次更新的梯度基础上减去一个梯度（ $\Theta$为模型参数向量，$\bigtriangledown f_t(\theta_t) $为t时刻loss函数的梯度，$\alpha$为学习率）:$\theta_t=(1-\lambda ) \theta_t - \alpha \bigtriangledown f_t(\theta_t)$

（3）L2 regularization是给参数加上一个L2惩罚($f_t(\theta )$为loss函数)：$f^{reg}_t(\theta )=f_t(\theta )+\frac{\lambda ‘}{2}||\theta ||^2$。(当$\lambda ‘=\frac{\lambda }{\alpha }$时，与weight decay等价，仅在使用标准SGD优化时成立)

（4）Adam自动调整学习率，大幅提高了训练速度，也很少需要调整学习率，但是有相当多的资料报告Adam优化的最终精度略低于SGD。问题出在哪呢，其实Adam本身没有问题，问题在于目前大多数DL框架的L2 regularization实现用的是weight decay的方式，而weight decay在与Adam共同使用的时候有互相耦合。理由如下图，其中红色是传统的Adam+L2 regularization的方式，绿色是bert使用的接入weight decay的方式，能够完成梯度下降与weight decay的解耦：

fig5

fig6

大部分的模型都会有L2 regularization约束项，因此很有可能出现adam的最终效果没有sgd的好。如果在tf里面要对不同区域的tensor做不同的L2 regularization调整的化可以参考https://zhuanlan.zhihu.com/p/40814046 先做adam后在手工更新L2 regularization梯度的方法。

推荐：一个框架看懂优化算法之异同 SGD/AdaGrad/Adam：https://zhuanlan.zhihu.com/p/32230623