优化算法集锦

• 梯度下降法
• 原理
• 可视化过程
• 系统实现
• 随机梯度下降法
• 原理
• 可视化过程
• 系统实现
• 小批量梯度下降法
• 原理
• 可视化过程
• 系统实现
• 动量法
• 提出原因
• 原理
• 指数加权平均
• 动量法原理
• 可视化过程
• 系统实现
• AdaGrad 法
• 提出原因
• 原理
• 可视化过程
• 系统实现
• RMSprop法
• 提出原因
• 原理
• 可视化过程
• 系统实现
• Adadealt 法
• 提出原因
• 原理
• 可视化过程
• 系统实现
• Adam 法
• 提出原因
• 原理
• 可视化过程
• 系统实现
• 小结

本文是对于深度学习中优化算法的解释与实现

首先需要明确：优化算法是为了尽可能低降低训练误差(Loss function)而不是泛化误差。

梯度下降法

原理

梯度下降法是深度学习值最常使用到的优化算法，也是后续算法的基石。下面以多维梯度下降法为例简述其原理：

假设对于损失函数 f ( x ) f(x) f(x),其输入 x = [ x 1 , x 2 , x 3 , . . . . . . , x n ] T x=[x_{1},x_{2},x_{3},......,x_{n} ]^{T} x=[x1​,x2​,x3​,......,xn​]T是一个n维向量，其输入为标量（即一个实数）。

则损失函数的梯度可表示为： ▽ f ( x ) = [ ∂ f ( x ) ∂ x 1 , ∂ f ( x ) ∂ x 2 , ∂ f ( x ) ∂ x 3 , . . . . . . , ∂ f ( x ) ∂ x n ] T \bigtriangledown f(x)=[\frac{\partial f(x)}{\partial x_{1}},\frac{\partial f(x)}{\partial x_{2}},\frac{\partial f(x)}{\partial x_{3}},......,\frac{\partial f(x)}{\partial x_{n}}]^{T} ▽f(x)=[∂x1​∂f(x)​,∂x2​∂f(x)​,∂x3​∂f(x)​,......,∂xn​∂f(x)​]T

其中 ∂ f ( x ) ∂ x i \frac{\partial f(x)}{\partial x_{i}} ∂xi​∂f(x)​表示损失函数关于样本 x x x中特征 x i x_{i} xi​的变化率。

依据方向导数性质， f ( x ) f(x) f(x)在某一方向 u u u上的变化率 D u f ( x ) D_{u}f(x) Du​f(x)为：
D u f ( x ) = ▽ f ( x ) ⋅ u = ▽ f ( x ) c o s ( θ ) ∣ ∣ u ∣ ∣ D_{u}f(x)=\bigtriangledown f(x)·u=\bigtriangledown f(x)cos(\theta)||u|| Du​f(x)=▽f(x)⋅u=▽f(x)cos(θ)∣∣u∣∣

为了尽快降低损失，我们需要该方向导数尽可能的为足够小的负数，从而当损失函数沿着该方向移动时，能够保证其下降最快。根据上式明显可以看出：当 θ = π \theta=\pi θ=π时，方向导数值最小，下降最快。故有：
x = x − α ▽ f ( x ) x=x-\alpha \bigtriangledown f(x) x=x−α▽f(x)。

这就是批量梯度下降法的公式。其中 α \alpha α称为学习率，它控制着每一步步长的相对大小。

当学习率过小时，更新速度较慢，此时需要更多迭代周期和更长时间。
当学习率过大时，步长较长，在各方向梯度数值差异不大且损失函数处于更新前期时，往往效果较好。但当各方向梯度数值差异较大时，容易造成某一方向上的振荡，从而不能达到最优值或足够好的次优值。

可视化过程

上面的解释可能较为晦涩，下面我们通过可视化的方式来更直观地感受这一点。
以 f ( x ) = x 1 2 + x 2 2 f(x)=x_{1}^{2}+x_{2} ^{2} f(x)=x12​+x22​为例，假设其初始位置为 x 1 = − 5 , x 2 = − 2 x_{1}=-5,x_{2}=-2 x1​=−5,x2​=−2。
可视化梯度下降代码如下：

import torch
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
def gradient_descent(num_epoch,lr):
    x1,x2=-5,-2
    result=[(x1,x2)]
    for i in range(num_epoch):
        g1=2*x1
        g2=2*x2
        x1-=lr*g1
        x2-=lr*g2
        result.append((x1,x2))
    fig,ax=plt.subplots(figsize=(6,6))
    x1,x2=np.meshgrid(np.arange(-5,2,0.1),np.arange(-2,2,0.1))
    ax.contour(x1,x2,x1**2+x2**2)
    ax.plot(*zip(*result),'-o')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
        
gradient_descent(10,0.1)
#迭代十轮，学习率设为0.1
        

结果如下：

可以看到，该损失函数从(-5,-2)出发，采用0.1的学习率，经过十次梯度下降后，已经接近了最小值点（0,0)。

将迭代次数设为20：

更加接近最优点，再将学习率设为0.2：

可以看到，此时几乎已经达到了最小值。

可得出结论：迭代次数与学习率都影响最后的结果，迭代次数的增加可以使得学习率不变的情况下，损失函数执行更多次梯度下降，从而更接近最优值。而学习率则控制每次下降的步长大小，可以使得在迭代次数不变的情况下，每次迭代梯度下降更多，这也就是为什么当我们调整学习率后，上图中的两次落脚点间距变大。

系统实现

这里我们使用NASA的数据集（共有1500余样本，每个样本有5个特征，以及一个飞机机翼噪声的标签）来实现梯度下降法。
首先处理数据：

def get_data_ch7():
    data=np.genfromtxt('airfoil_self_noise.dat',delimiter='\t')
    data=(data-data.mean(axis=0))/data.std(axis=0)
    return torch.tensor(data[:1500,:-1],dtype=torch.float32),torch.tensor(data[:1500,-1],dtype=torch.float32)
features,labels=get_data_ch7()
features.shape

然后定义模型并训练：

net=torch.nn.Sequential(
       torch.nn.Linear(features.shape[1],1)
    )
def train(net,num_epochs,features,labels,device,optimizer,loss):
    net=net.to(device)
    features=features.to(device)
    labels=labels.to(device)
    def get_loss():
        return loss(net(features).view(-1,1),labels).cpu().item()
    result=[get_loss()]
    for epoch in range(num_epochs):
        y_hat=net(features)
        l=loss(y_hat.view(-1,1),labels)
        optimizer.zero_grad()
        l.backward()
        optimizer.step()
        result.append(get_loss())
    fig,ax=plt.subplots(figsize=(6,6))
    ax.plot(np.linspace(0,num_epochs,len(result)),result)
    ax.set_xlabel('epoch')
    ax.set_ylabel('loss')
    print('Loss:',result[-1])
device=torch.device('cuda' if torch.cuda.is_available() else 'cpu')
optimizer=torch.optim.SGD(net.parameters(),lr=0.1)
loss=torch.nn.MSELoss()
train(net,20,features,labels,device,optimizer,loss) 

此处采用学习率为0.1，迭代周期为20，结果如下：

可以看到，在第15轮迭代时，损失函数基本趋于平缓，最终损失降到了1.0左右。虽然该方法可以有效地降低损失函数值，但其效率较低。通过上述代码可以观察到，梯度下降法在每个迭代周期只进行一次梯度下降，也即每遍历1500个样本进行一次梯度下降，这导致了采用梯度下降法所需要的时间较长。除此之外，实际中样本量可能是百万或千万级，如果采用该方法，则每次迭代自变量的代价过高。针对这个问题，随机梯度下降法给予了一定改进。

随机梯度下降法

原理

随机梯度下降法则相对梯度下降法走向了另一个极端。**随机梯度下降法采用随机均匀采样的一个样本的梯度来代替所有样本的梯度从而更新参数。这样做的原因是每一个样本的梯度其实是总体梯度的无偏估计。**所以相对拥有1500个样本的上例来说，梯度下降法每一个世代进行一次梯度下降，而随机梯度下降每一个世代可以进行1500次梯度下降，效率大大提高。但是，也正因为每次选取的是一个样本的梯度，所以难免会在过程中左右振荡（某些样本具有干扰性），并不会一直指向最优值方向。

可视化过程

代码如下：

def gradient_descent(num_epoch,lr):
    x1,x2=-5,-2
    result=[(x1,x2)]
    for i in range(num_epoch):
        g1=2*x1
        g2=2*x2
        x1-=lr*(g1+np.random.normal(0,0.3))
        x2-=lr*(g2+np.random.normal(0,0.3))
        result.append((x1,x2))
    fig,ax=plt.subplots(figsize=(6,6))
    x1,x2=np.meshgrid(np.arange(-5,2,0.1),np.arange(-2,2,0.1))
    ax.contour(x1,x2,x1**2+x2**2)
    ax.plot(*zip(*result),'-o')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    print('Final (x1,x2):',result[-1])
        
gradient_descent(20,0.2)

可以很明显看到，每次梯度更新方向并不一直指向最小值点，但总体仍然会指向最优点。当到达最优点附近时，并不会停靠在某一位置，而是会在某一小范围内振荡，这导致我们往往无法达到最优点。

系统实现

代码如下：

net=torch.nn.Sequential(
       torch.nn.Linear(features.shape[1],1)
    )
def train(net,num_epochs,features,labels,device,optimizer,loss,batch_size=10):
    net=net.to(device)
    features=features.to(device)
    labels=labels.to(device)
    data_iter=torch.utils.data.DataLoader(torch.utils.data.TensorDataset(features,labels),batch_size,shuffle=True)
    def get_loss():
        return loss(net(features).view(-1,1),labels).cpu().item()
    result=[get_loss()]
    for epoch in range(num_epochs):
        for batch_count,(X,y) in enumerate(data_iter):
            y_hat=net(X)
            l=loss(y_hat.view(-1,1),y)
            optimizer.zero_grad()
            l.backward()
            optimizer.step()
            if (batch_count+1)*batch_size%100==0:
                result.append(get_loss())
    fig,ax=plt.subplots(figsize=(6,6))
    ax.plot(np.linspace(0,num_epochs,len(result)),result)
    ax.set_xlabel('epoch')
    ax.set_ylabel('loss')
    print('Loss:',result[-1])
device=torch.device('cuda' if torch.cuda.is_available() else 'cpu')
optimizer=torch.optim.SGD(net.parameters(),lr=0.0001)
loss=torch.nn.MSELoss()
train(net,2,features,labels,device,optimizer,loss,1)  

这里采用的是2个世代，0.0001的学习率，结果如下：
可以看到，相对于梯度下降法，其只需到第2个世代就可以将损失降到与梯度下降法经过15个世代相同的水平。但是，也正是因为每一个世代都有1500次梯度下降，所以随机梯度下降法往往配合学习率衰减或采用极小的学习率（否则后期会振荡剧烈）。同时，由于其对于世代中每一个样本都需要执行梯度下降，故耗时远大于梯度下降法。

基于效果和效率的考虑，衍生出了梯度下降法和随机梯度下降法的折中方案：小批量梯度下降法。

小批量梯度下降法

原理

该方法为梯度下降法和随机梯度下降法的折中，事实上深度学习领域很多算法都采用了折中思想。

小批量梯度下降法每次从样本中采样一个小批量（包含多个样本），利用小批量中样本的平均梯度去代替全局梯度更新参数，降低损失函数。小批量中样本的平均梯度也是全局梯度的无偏估计。

可视化过程

代码如下：

def gradient_descent(num_epoch,lr):
    x1,x2=-5,-2
    result=[(x1,x2)]
    for i in range(num_epoch):
        g1=2*x1
        g2=2*x2
        x1-=lr*(g1+np.random.normal(0,0.1))
        x2-=lr*(g2+np.random.normal(0,0.1))
        result.append((x1,x2))
    fig,ax=plt.subplots(figsize=(6,6))
    x1,x2=np.meshgrid(np.arange(-5,2,0.1),np.arange(-2,2,0.1))
    ax.contour(x1,x2,x1**2+x2**2)
    ax.plot(*zip(*result),'-o')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    print('Final (x1,x2):',result[-1])
        
gradient_descent(20,0.2)

结果如下：

小批量梯度下降的指向稳定性介于梯度下降和随机梯度下降间。

系统实现

代码如下：

net=torch.nn.Sequential(
       torch.nn.Linear(features.shape[1],1)
    )
def train(net,num_epochs,features,labels,device,optimizer,loss,batch_size=10):
    net=net.to(device)
    features=features.to(device)
    labels=labels.to(device)
    data_iter=torch.utils.data.DataLoader(torch.utils.data.TensorDataset(features,labels),batch_size,shuffle=True)
    def get_loss():
        return loss(net(features).view(-1,1),labels).cpu().item()
    result=[get_loss()]
    for epoch in range(num_epochs):
        for batch_count,(X,y) in enumerate(data_iter):
            y_hat=net(X)
            l=loss(y_hat.view(-1,1),y)
            optimizer.zero_grad()
            l.backward()
            optimizer.step()
            if (batch_count+1)*batch_size%100==0:
                result.append(get_loss())
    fig,ax=plt.subplots(figsize=(6,6))
    ax.plot(np.linspace(0,num_epochs,len(result)),result)
    ax.set_xlabel('epoch')
    ax.set_ylabel('loss')
    print('Loss:',result[-1])
device=torch.device('cuda' if torch.cuda.is_available() else 'cpu')
optimizer=torch.optim.SGD(net.parameters(),lr=0.01)
loss=torch.nn.MSELoss()
train(net,5,features,labels,device,optimizer,loss,10)  

这里我们采用5个世代，学习率为0.01，每个批量含10个样本，结果如下：

可以看到，小批量梯度下降的速度介于梯度下降和随机梯度下降之间。使用的学习率也应介于梯度下降和随机梯度下降所使用的学习率间。

动量法

提出原因

上述三种梯度下降方法理论上可以解决任何的问题，但有时往往效率较低。例如当多方向的梯度差异较大时，优化轨迹往往有振荡现象。

下面以优化 f ( x ) = 0.1 x 1 2 + 2 x 2 2 f(x)=0.1x_{1}^ {2}+2x_{2}^{2} f(x)=0.1x12​+2x22​为例来说明这一点，假设其初始状态位于 x 1 = − 5 , x 2 2 = − 2 x_{1}=-5,x2_{2}=-2 x1​=−5,x22​=−2。

代码如下：

def gradient_descent(num_epoch,lr):
    x1,x2=-5,-2
    result=[(x1,x2)]
    for i in range(num_epoch):
        g1=0.2*x1
        g2=4*x2
        x1-=lr*(g1+np.random.normal(0,0.1))
        x2-=lr*(g2+np.random.normal(0,0.1))
        result.append((x1,x2))
    fig,ax=plt.subplots(figsize=(6,6))
    x1,x2=np.meshgrid(np.arange(-5,2,0.1),np.arange(-2,2,0.1))
    ax.contour(x1,x2,0.1*x1**2+2*x2**2)
    ax.plot(*zip(*result),'-o')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    print('Final (x1,x2):',result[-1])
        
gradient_descent(20,0.4)

结果如下：

可以看到，前期梯度振荡剧烈，这是因为在 x 2 x_{2} x2​方向上的梯度远大于在 x 1 x_{1} x1​方向上的梯度。不过上图看起来并不会影响最终结果，下面我们将学习率调制0.8，得到如下结果：

可以看到，此时 x 2 x_{2} x2​方向上的最终值完全偏离了最优值。

虽然最终情况可以通过调小学习率改善（以下为学习率为0.2的情况）：

但这需要迭代世代的大量增长，代价过高。

为此，我们希望可以找到一种方法，使得梯度下降时 x 1 x_{1} x1​方向上的梯度更大（更快接近最优值）而 x 2 x_{2} x2​方向上的梯度更小（减小振荡），从而允许我们使用较大的学习率。为此，动量法被提出。

原理

指数加权平均

动量法以及后续介绍的算法大多与指数加权平均紧密相连，为此，我们需要先了解指数加权平均。

给定一组样本 ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . . . . , ( x 200 , y 200 ) (x_{1},y_{1}),(x_{2},y_{2}),......,(x_{200},y_{200}) (x1​,y1​),(x2​,y2​),......,(x200​,y200​)，并将其可视化。
代码如下：

x=np.arange(-100,100,1)
y=np.zeros(x.shape)
for i in range(len(x)):
    y[i]=x[i]**2+np.random.normal(10,1)**3.5


fig,ax=plt.subplots(figsize=(6,6))
ax.scatter(x,y)

结果如下：

下面定义指数移动加权平均的式子：

v t = β ∗ v t − 1 + ( 1 − β ) ∗ y t v_{t}=\beta*v_{t-1}+(1-\beta)*y_{t} vt​=β∗vt−1​+(1−β)∗yt​
v 0 = 0 v_{0}=0 v0​=0

我们先将其中的 β \beta β设为0.9，并画出其与x的曲线图，代码如下：

def expon_average(lambd,x):
    v=0
    result=[]
    for i in range(len(x)):
        v=lambd*v+(1-lambd)*y[i]
        result.append(v)
    ax.plot(x,result,color='r')
expon_average(0.5,x)  

结果如下：

可以看到，此时 v t v_{t} vt​粗略地拟合了 y y y。下面以 β = 0.9 \beta=0.9 β=0.9为例解释指数移动加权平均的作用。

根据公式，我们有：

v 1 = β ∗ v 0 + ( 1 − β ) ∗ y 1 v_{1}=\beta*v_{0}+(1-\beta)*y_{1} v1​=β∗v0​+(1−β)∗y1​
v 2 = β ∗ v 1 + ( 1 − β ) ∗ y 2 v_{2}=\beta*v_{1}+(1-\beta)*y_{2} v2​=β∗v1​+(1−β)∗y2​
v 3 = β ∗ v 2 + ( 1 − β ) ∗ y 3 v_{3}=\beta*v_{2}+(1-\beta)*y_{3} v3​=β∗v2​+(1−β)∗y3​
.
.
.
v 198 = β ∗ v 197 + ( 1 − β ) ∗ y 198 v_{198}=\beta*v_{197}+(1-\beta)*y_{198} v198​=β∗v197​+(1−β)∗y198​
v 199 = β ∗ v 198 + ( 1 − β ) ∗ y 199 v_{199}=\beta*v_{198}+(1-\beta)*y_{199} v199​=β∗v198​+(1−β)∗y199​
v 200 = β ∗ v 199 + ( 1 − β ) ∗ y 200 v_{200}=\beta*v_{199}+(1-\beta)*y_{200} v200​=β∗v199​+(1−β)∗y200​

也即：
v 200 = β ∗ v 199 + ( 1 − β ) ∗ y 200 v_{200}=\beta*v_{199}+(1-\beta)*y_{200} v200​=β∗v199​+(1−β)∗y200​
= y 200 ∗ 0.1 + 0.9 ∗ ( 0.9 ∗ v 198 + 0.1 ∗ y 199 ) =y_{200}*0.1+0.9*(0.9*v_{198}+0.1*y_{199}) =y200​∗0.1+0.9∗(0.9∗v198​+0.1∗y199​)
= y 200 ∗ 0.1 + 0.9 ∗ ( 0.9 ∗ ( 0.9 ∗ v 196 + 0.1 ∗ y 197 ) + 0.1 ∗ y 199 ) =y_{200}*0.1+0.9*(0.9*(0.9*v_{196}+0.1*y_{197})+0.1*y_{199}) =y200​∗0.1+0.9∗(0.9∗(0.9∗v196​+0.1∗y197​)+0.1∗y199​)
= 0.1 ∗ y 200 + 0. 9 2 ∗ 0.1 ∗ y 199 + 0. 9 3 ∗ 0.1 ∗ y 198 + . . . . . . 0. 9 n ∗ 0.1 ∗ y 200 − n + 1 =0.1*y_{200}+0.9^{2}*0.1*y_{199}+0.9^{3}*0.1*y_{198}+......0.9^{n}*0.1*y_{200-n+1} =0.1∗y200​+0.92∗0.1∗y199​+0.93∗0.1∗y198​+......0.9n∗0.1∗y200−n+1​

可以看到，对于每一个样本 0.1 ∗ y i 0.1*y_{i} 0.1∗yi​，其权重逐渐减小，为 β 201 − i \beta^{201-i} β201−i。

根据极限 l i m n → ∞ ( 1 − 1 n ) n = 1 e lim_{n\rightarrow \infty}(1-\frac{1}{n})^{n}=\frac{1}{e} limn→∞​(1−n1​)n=e1​，可知当项数超过 1 1 − β \frac{1}{1-\beta} 1−β1​时，其权重衰减至大约 1 3 \frac{1}{3} 31​,此时我们认为，权重值低于此的便不再考虑（因为其分量太小），而余下的项，权重加起来接近于1，且各项权重差别不大，所以该算法其实是对近 1 1 − β \frac{1}{1-\beta} 1−β1​项的近似平均。而由于 v t v{t} vt中t的不断增大，其平均的 1 1 − β \frac{1}{1-\beta} 1−β1​项也随后移动，故该算法被称为指数移动加权平均。

而 β \beta β则控制了平均项数的多少。

当其较大时，平均的项数较少，绘制的指数加权平均曲线变化更为剧烈（或说更为贴近原数据点变化趋势）。
当其较小时，平均项数较多，绘制的指数加权平均曲线变换更为缓慢。

下面我们将 β \beta β分别设为0.95和0.5,观察曲线形状，代码同上，修改参数即可，结果如下：

β = 0.95 \beta=0.95 β=0.95时：

β = 0.5 \beta=0.5 β=0.5时：

可以看到，两者的平缓程度大不相同。

注意到当 β \beta β较大时，在曲线的起始阶段往往拟合值较真实值较小。产生这样现象的原因其实是因为 β \beta β较大时，每一次拟合都大部分依赖于上一次的 v v v值，而 v v v值初始被设置为0，且 y y y的权重较小，故前期当 v v v尚且处于增长阶段时，所拟合出的曲线相对真实位置偏低。

要解决这个问题只需要做偏差修正即可，即将拟合值从 v v v变为 v 1 − β t \frac{v}{1-\beta^{t}} 1−βtv​。当t较小时，该值相对于 v v v值会变大，而当t较大时，该值近似于 v v v值，所以偏差修正几乎只在前期产生作用。

下面使用偏差修正来改善 β = 0.95 \beta=0.95 β=0.95时的情况，代码如下：

def expon_average(lambd,x):
    v=0
    result=[]
    for i in range(len(x)):
        v=lambd*v+(1-lambd)*y[i]
        result.append(v/(1-lambd**(i+1)))
    ax.plot(x,result,color='r')
expon_average(0.95,x)  

结果如下：

可以看到，情况得到了明显改善。

动量法原理

根据上述指数加权平均的思想，为了改善梯度下降中的振荡问题，我们选择对每一个小批量中的梯度做指数加权平均，这样更新时所用到的梯度就不仅仅取决于一个批量，而是多个批量梯度的矢量平均。而上述例子中 x 2 x_{2} x2​的梯度有向上和向下两个方向，多批量的梯度平均使得 x 2 x_{2} x2​方向上的梯度被抵消了一部分，从而减缓了该方向上更新时的振荡。其次， x 1 x_{1} x1​方向上的梯度大多向最优值方向，移动指数加权平均使得该方向上梯度值更加精确地指向最优值方向。

其算法公式如下：

v t = β ∗ v t − 1 + ( 1 − β ) ∗ g t v_{t}=\beta*v_{t-1}+(1-\beta)*g_{t} vt​=β∗vt−1​+(1−β)∗gt​
v 0 = 0 v_{0}=0 v0​=0
p a r a m s = p a r a m s − α ∗ v t params=params-\alpha*v_{t} params=params−α∗vt​

其中 g t g_{t} gt​为当期时间步 t t t的小批量梯度。

可视化过程

仍然采用与上述梯度下降中相同的函数和学习率进行测试，这里 β = 0.9 \beta=0.9 β=0.9。
代码如下：

def gradient_descent(lambd,num_epoch,lr):
    x1,x2=-5,-2
    result=[(x1,x2)]
    g1,g2=0.0,0.0
    for i in range(num_epoch):
        g1=lambd*g1+(1-lambd)*0.2*x1
        g2=lambd*g2+(1-lambd)*4*x2
        x1-=lr*g1
        x2-=lr*g2
        result.append((x1,x2))
    fig,ax=plt.subplots(figsize=(6,6))
    x1,x2=np.meshgrid(np.arange(-5,2,0.1),np.arange(-2,2,0.1))
    ax.contour(x1,x2,0.1*x1**2+2*x2**2)
    ax.plot(*zip(*result),'-o')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    print('Final (x1,x2):',result[-1])
        
gradient_descent(0.9,30,0.8)

结果如下：

可以看到，振荡幅度明显减小，最终值更加接近最优值，所以动量法允许我们使用更大的学习率加快更新速度。

系统实现

代码如下：

net=torch.nn.Sequential(
       torch.nn.Linear(features.shape[1],1)
    )
def train(net,num_epochs,features,labels,device,optimizer,loss,batch_size=10):
    net=net.to(device)
    features=features.to(device)
    labels=labels.to(device)
    data_iter=torch.utils.data.DataLoader(torch.utils.data.TensorDataset(features,labels),batch_size,shuffle=True)
    def get_loss():
        return loss(net(features).view(-1,1),labels).cpu().item()
    result=[get_loss()]
    for epoch in range(num_epochs):
        for batch_count,(X,y) in enumerate(data_iter):
            y_hat=net(X)
            l=loss(y_hat.view(-1,1),y)
            optimizer.zero_grad()
            l.backward()
            optimizer.step()
            if (batch_count+1)*batch_size%100==0:
                result.append(get_loss())
    fig,ax=plt.subplots(figsize=(6,6))
    ax.plot(np.linspace(0,num_epochs,len(result)),result)
    ax.set_xlabel('epoch')
    ax.set_ylabel('loss')
    print('Loss:',result[-1])
device=torch.device('cuda' if torch.cuda.is_available() else 'cpu')
optimizer=torch.optim.SGD(net.parameters(),lr=0.0001,momentum=0.9)
loss=torch.nn.MSELoss()
train(net,10,features,labels,device,optimizer,loss,10) 

这里将 β \beta β设为0.9，批量尺寸设为10，结果如下：

同样也可以得到十分不错的结果。

AdaGrad 法

提出原因

上文我们看到，动量法通过对梯度进行指数加权平均成功缓解了更新时的振荡问题，从而允许我们使用更大的学习率。

本次要介绍的AdaGrad法，通过对于不同梯度使用不同的学习率从而达到减缓振荡的目的，其最终效果与动量法类似。

原理

下面给出公式：

s t = s t − 1 + g t 2 s_{t}=s_{t-1}+g_{t}^{2} st​=st−1​+gt2​
s 0 = 0 s_{0}=0 s0​=0
p a r a m s = p a r a m s − α ∗ g t s t + ϵ params=params-\frac{\alpha*g_{t}}{\sqrt{s_{t}+\epsilon}} params=params−st​+ϵ ​α∗gt​​

其中关于梯度的运算均为按元素运算。 ϵ \epsilon ϵ是一个常量，主要作用是为了保证数值稳定，除此之外没有较大的影响，一般取 1 0 − 6 − 1 0 − 8 10^{-6}-10 ^{-8} 10−6−10−8。

可以看到， s t s_{t} st​每次叠加一个批量的梯度平方，每一维特征都有自己的梯度，将这些梯度平方后加到 s t s_{t} st​上，再将学习率设置为 α s t + ϵ \frac{\alpha}{\sqrt{s_{t}+\epsilon}} st​+ϵ ​α​。这使得各维特征都具有自己的学习率而不必与其它特征保持一致。

一方面，这样做使得梯度较大的特征得到了较小的学习率，梯度较小的学习率得到了较大的学习率，从而均衡了步长，一定程度上缓解了振荡。

另一方面，不论是哪一维特征的梯度，总是伴随着学习率的衰减（因为 s t s_{t} st​每次都会加上一个批量梯度的平方）。如果学习率衰减到较小时仍未达到较好解，那么之后将很难达到最优解。（因为学习率一直在衰减，需要更多的迭代周期）。

可视化过程

代码如下：

def gradient_descent(lambd,num_epoch,lr):
    x1,x2=-5,-2
    result=[(x1,x2)]
    g1,g2=0.0,0.0
    s1,s2=0.0,0.0
    eps=1e-8
    for i in range(num_epoch):
        g1=lambd*g1+(1-lambd)*0.2*x1
        g2=lambd*g2+(1-lambd)*4*x2
        s1+=g1**2
        s2+=g2**2
        x1-=lr*g1/(np.sqrt(s1+eps))
        x2-=lr*g2/(np.sqrt(s2+eps))
        result.append((x1,x2))
    fig,ax=plt.subplots(figsize=(6,6))
    x1,x2=np.meshgrid(np.arange(-5,2,0.1),np.arange(-2,2,0.1))
    ax.contour(x1,x2,0.1*x1**2+2*x2**2)
    ax.plot(*zip(*result),'-o')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    print('Final (x1,x2):',result[-1])
        
gradient_descent(0.9,50,0.3)

结果如下：

可以看到，后期的步幅非常小，这是因为学习率一直在衰减导致的。由此，我们不得不使用更多的时代以达到满意值（这里使用了50代，而动量法只使用了30代），不过这种持续的衰减允许我们使用较大的初始学习率。

系统实现

代码如下：

    labels=labels.to(device)
    data_iter=torch.utils.data.DataLoader(torch.utils.data.TensorDataset(features,labels),batch_size,shuffle=True)
    def get_loss():
        return loss(net(features).view(-1,1),labels).cpu().item()
    result=[get_loss()]
    for epoch in range(num_epochs):
        for batch_count,(X,y) in enumerate(data_iter):
            y_hat=net(X)
            l=loss(y_hat.view(-1,1),y)
            optimizer.zero_grad()
            l.backward()
            optimizer.step()
            if (batch_count+1)*batch_size%100==0:
                result.append(get_loss())
    fig,ax=plt.subplots(figsize=(6,6))
    ax.plot(np.linspace(0,num_epochs,len(result)),result)
    ax.set_xlabel('epoch')
    ax.set_ylabel('loss')
    print('Loss:',result[-1])
device=torch.device('cuda' if torch.cuda.is_available() else 'cpu')
net=net.to(device)
optimizer=torch.optim.Adagrad(net.parameters(),lr=0.01)
loss=torch.nn.MSELoss()
train(net,10,features,labels,device,optimizer,loss,10)      

结果如下：

效果也很不错。

RMSprop法

提出原因

对于Adagrad算法，其通过自动的学习率衰减使得各方向上梯度更为均衡。但往往会因为后期学习率过小从而难以达到最优值或足够优的次优值。

RMSprop算法对于Adagrad算法做了改进，其采用的 s t s_{t} st​不再是梯度平方的叠加，而是梯度的平方的指数加权平均，这使得 s t s_{t} st​的变换更为迟缓。在保持梯度均衡的同时减慢了学习率的衰减，使得梯度下降更容易到达最优值或足够优的此优值。

原理

算法公式如下：
s t = β ∗ s t − 1 + ( 1 − β ) ∗ g t 2 s_{t}=\beta*s_{t-1}+(1-\beta)*g_{t}^{2} st​=β∗st−1​+(1−β)∗gt2​
s 0 = 0 s_{0}=0 s0​=0
p a r a m s = p a r a m s − α ∗ g t s t + ϵ params=params-\frac{\alpha*g_{t}}{\sqrt{s_{t}+\epsilon}} params=params−st​+ϵ ​α∗gt​​

可以看到，此时的 s t s_{t} st​为梯度平方的指数加权平均，使得学习率不再一昧衰减。而是随着梯度的大小变化，当梯度较大时，学习率较小，当梯度较小时，学习率较大。这样灵活的调整使得在避免振荡的同时没有减缓优化速度。

可视化过程

代码如下：

def gradient_descent(lambd,num_epoch,lr):
    x1,x2=-5,-2
    result=[(x1,x2)]
    g1,g2=0.0,0.0
    s1,s2=0.0,0.0
    eps=1e-8
    for i in range(num_epoch):
        g1=0.2*x1
        g2=4*x2
        s1=lambd*s1+(1-lambd)*g1**2
        s2=lambd*s2+(1-lambd)*g2**2
        x1-=lr*g1/(np.sqrt(s1+eps))
        x2-=lr*g2/(np.sqrt(s2+eps))
        result.append((x1,x2))
    fig,ax=plt.subplots(figsize=(6,6))
    x1,x2=np.meshgrid(np.arange(-5,2,0.1),np.arange(-2,2,0.1))
    ax.contour(x1,x2,0.1*x1**2+2*x2**2)
    ax.plot(*zip(*result),'-o')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    print('Final (x1,x2):',result[-1])
        
gradient_descent(0.9,20,0.3)

结果如下：

可以看到，使用相同的学习率，RMSprop只需要20世代便可以达到足够优秀的值，而Adaprop则需要50代。

系统实现

代码如下：

net=torch.nn.Sequential(
       torch.nn.Linear(features.shape[1],1)
    )
def train(net,num_epochs,features,labels,device,optimizer,loss,batch_size=10):
    net=net.to(device)
    features=features.to(device)
    labels=labels.to(device)
    data_iter=torch.utils.data.DataLoader(torch.utils.data.TensorDataset(features,labels),batch_size,shuffle=True)
    def get_loss():
        return loss(net(features).view(-1,1),labels).cpu().item()
    result=[get_loss()]
    for epoch in range(num_epochs):
        for batch_count,(X,y) in enumerate(data_iter):
            y_hat=net(X)
            l=loss(y_hat.view(-1,1),y)
            optimizer.zero_grad()
            l.backward()
            optimizer.step()
            if (batch_count+1)*batch_size%100==0:
                result.append(get_loss())
    fig,ax=plt.subplots(figsize=(6,6))
    ax.plot(np.linspace(0,num_epochs,len(result)),result)
    ax.set_xlabel('epoch')
    ax.set_ylabel('loss')
    print('Loss:',result[-1])
device=torch.device('cuda' if torch.cuda.is_available() else 'cpu')
net=net.to(device)
optimizer=torch.optim.RMSprop(net.parameters(),lr=0.002,alpha=0.9)
loss=torch.nn.MSELoss()
train(net,10,features,labels,device,optimizer,loss,10)   

结果如下：

效果不错，优化速度也很快。

Adadealt 法

提出原因

针对Adagrad算法在后期较难找到优解的问题，Adadealt算法也做了一些改进。比较特殊的是，该算法并没有学习率。

原理

公式如下：

s t = β ∗ s t − 1 + ( 1 − β ) ∗ g t 2 s_{t}=\beta*s_{t-1}+(1-\beta)*g_{t}^{2} st​=β∗st−1​+(1−β)∗gt2​
s 0 = 0 s_{0}=0 s0​=0
d = Δ x t − 1 + ϵ s t + ϵ ∗ g t d=\frac{\sqrt{\Delta x_{t-1}+\epsilon}}{\sqrt{s_{t}+\epsilon}}*g_{t} d=st​+ϵ ​Δxt−1​+ϵ ​​∗gt​
Δ x t = β ∗ Δ x t − 1 + ( 1 − β ) ∗ d 2 \Delta x_{t}=\beta*\Delta x_{t-1}+(1-\beta)*d^{2} Δxt​=β∗Δxt−1​+(1−β)∗d2
p a r a m s = p a r a m s − d params=params-d params=params−d

该算法使用 Δ x t − 1 + ϵ \sqrt{\Delta x_{t-1}+\epsilon} Δxt−1​+ϵ ​来代替学习率。这一项是关于参数变化量的函数。会随着参数变化量自动调整。

可视化过程

代码如下：

def gradient_descent(lambd,num_epoch,lr):
    x1,x2=-5,-2
    result=[(x1,x2)]
    g1,g2=0.0,0.0
    s1,s2=0.0,0.0
    d1,d2=0.0,0.0
    g_1,g_2=0.0,0.0
    eps=1e-8
    for i in range(num_epoch):
        g1=0.2*x1
        g2=4*x2
        s1=lambd*s1+(1-lambd)*g1**2
        s2=lambd*s2+(1-lambd)*g2**2
        g_1=np.sqrt((d1+eps)/(s1+eps))*g1
        g_2=np.sqrt((d2+eps)/(s2+eps))*g2
        d1=lambd*d1+(1-lambd)*g_1**2
        d2=lambd*d2+(1-lambd)*g_2**2
        x1=x1-g_1
        x2=x2-g_2
        result.append((x1,x2))
    fig,ax=plt.subplots(figsize=(6,6))
    x1,x2=np.meshgrid(np.arange(-5,2,0.1),np.arange(-2,2,0.1))
    ax.contour(x1,x2,0.1*x1**2+2*x2**2)
    ax.plot(*zip(*result),'-o')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    print('Final (x1,x2):',result[-1])
        
gradient_descent(0.9,5000,0.1)

结果如下：

效果不错。

系统实现

代码如下：

net=torch.nn.Sequential(
       torch.nn.Linear(features.shape[1],1)
    )
def train(net,num_epochs,features,labels,device,optimizer,loss,batch_size=10):
    net=net.to(device)
    features=features.to(device)
    labels=labels.to(device)
    data_iter=torch.utils.data.DataLoader(torch.utils.data.TensorDataset(features,labels),batch_size,shuffle=True)
    def get_loss():
        return loss(net(features).view(-1,1),labels).cpu().item()
    result=[get_loss()]
    for epoch in range(num_epochs):
        for batch_count,(X,y) in enumerate(data_iter):
            y_hat=net(X)
            l=loss(y_hat.view(-1,1),y)
            optimizer.zero_grad()
            l.backward()
            optimizer.step()
            if (batch_count+1)*batch_size%100==0:
                result.append(get_loss())
    fig,ax=plt.subplots(figsize=(6,6))
    ax.plot(np.linspace(0,num_epochs,len(result)),result)
    ax.set_xlabel('epoch')
    ax.set_ylabel('loss')
    print('Loss:',result[-1])
device=torch.device('cuda' if torch.cuda.is_available() else 'cpu')
net=net.to(device)
optimizer=torch.optim.Adadelta(net.parameters(),rho=0.9)
loss=torch.nn.MSELoss()
train(net,10,features,labels,device,optimizer,loss,10)      

结果如下：

同样能取得不错的效果。

Adam 法

提出原因

Adam算法是动量法和RMSprop算法的结合，是目前使用最广泛的优化算法。

原理

公式如下：
s t = β 1 ∗ s t − 1 + ( 1 − β 1 ) ∗ g t 2 s_{t}=\beta_{1}*s_{t-1}+(1-\beta_{1})*g_{t}^{2} st​=β1​∗st−1​+(1−β1​)∗gt2​
v t = β 2 ∗ v t − 1 + ( 1 − β 2 ) ∗ g t v_{t}=\beta_{2}*v_{t-1}+(1-\beta_{2})*g_{t} vt​=β2​∗vt−1​+(1−β2​)∗gt​
s t h a t = s t 1 − β 1 t s^ {hat}_{t}=\frac{s_{t}}{1-\beta_{1}^{t}} sthat​=1−β1t​st​​
v t h a t = v t 1 − β 2 t v^ {hat}_{t}=\frac{v_{t}}{1-\beta_{2}^{t}} vthat​=1−β2t​vt​​
p a r a m s = p a r a m s − α ∗ v t h a t s t h a t + ϵ params=params-\frac{\alpha*v^{hat}_{t}}{\sqrt{s^{hat}_{t}+\epsilon}} params=params−sthat​+ϵ ​α∗vthat​​

该算法不仅对小批量梯度做了指数加权平均（动量法），还对小批量梯度平方做了指数加权平均（RMSprop)。除此之外，引入了偏差修正，整体效果上来说是诸多算法中最优的一个。

可视化过程

代码如下：

def gradient_descent(lambd1,lambd2,num_epoch,lr):
    x1,x2=-5,-2
    result=[(x1,x2)]
    g1,g2=0.0,0.0
    s1,s2=0.0,0.0
    v1,v2=0.0,0.0
    eps=1e-8
    for i in range(num_epoch):
        g1=0.2*x1
        g2=4*x2
        s1=lambd1*s1+(1-lambd1)*(g1**2)
        s2=lambd1*s2+(1-lambd1)*(g2**2)
        v1=lambd2*v1+(1-lambd2)*g1
        v2=lambd2*v2+(1-lambd2)*g2
        v1=v1/(1-lambd2**(i+1))
        v2=v2/(1-lambd2**(i+1))
        s1=s1/(1-lambd1**(i+1))
        s2=s2/(1-lambd1**(i+1))
        x1-=lr*v1/np.sqrt(s1+eps)
        x2-=lr*v2/np.sqrt(s2+eps)
        result.append((x1,x2))
    fig,ax=plt.subplots(figsize=(6,6))
    x1,x2=np.meshgrid(np.arange(-5,2,0.1),np.arange(-2,2,0.1))
    ax.contour(x1,x2,0.1*x1**2+2*x2**2)
    ax.plot(*zip(*result),'-o')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    print('Final (x1,x2):',result[-1])
        
gradient_descent(0.999,0.9,30,2)

结果如下：

效果不错。

系统实现

代码如下：

net=torch.nn.Sequential(
       torch.nn.Linear(features.shape[1],1)
    )
def train(net,num_epochs,features,labels,device,optimizer,loss,batch_size=10):
    net=net.to(device)
    features=features.to(device)
    labels=labels.to(device)
    data_iter=torch.utils.data.DataLoader(torch.utils.data.TensorDataset(features,labels),batch_size,shuffle=True)
    def get_loss():
        return loss(net(features).view(-1,1),labels).cpu().item()
    result=[get_loss()]
    for epoch in range(num_epochs):
        for batch_count,(X,y) in enumerate(data_iter):
            y_hat=net(X)
            l=loss(y_hat.view(-1,1),y)
            optimizer.zero_grad()
            l.backward()
            optimizer.step()
            if (batch_count+1)*batch_size%100==0:
                result.append(get_loss())
    fig,ax=plt.subplots(figsize=(6,6))
    ax.plot(np.linspace(0,num_epochs,len(result)),result)
    ax.set_xlabel('epoch')
    ax.set_ylabel('loss')
    print('Loss:',result[-1])
device=torch.device('cuda' if torch.cuda.is_available() else 'cpu')
net=net.to(device)
optimizer=torch.optim.Adam(net.parameters(),lr=0.01)
loss=torch.nn.MSELoss()
train(net,10,features,labels,device,optimizer,loss,10)

结果如下：

效果不错。

小结

1.梯度下降法每世代更新一次，随机梯度下降法每世代更新样本数次，小批量梯度下降法介于两者之间。
2.动量法的提出缓解了梯度下降时因各方向梯度差异较大导致的振荡问题。
3.Adagrad算法则通过为不同特征分配不同的学习率来缓解振荡问题。
4.RMSprop算法采用了梯度平方指数加权平均缓解了Adagrad算法难以达到最优解的问题。
5.Adadealt算法则通过另一个方向缓解了该问题，此算法中以参数变化量平方的指数移动加权平均的二次根式来代替学习率。
6.Adam算法综合了动量法与RMSprop算法，是目前应用最广泛的优化算法。