深度学习中的优化算法

优化与深度学习

在深度学习问题中，我们通常会预先定义一个损失函数，有了损失函数，我们就可以使用优化算法试图将其最小化。在优化中，这样的损失函数通常被称作优化问题的目标函数。

虽然优化为深度学习提供了最小化损失函数的方法，但本质上，优化与深度学习的目标是有区别的，由于优化算法的目标函数通常是一个基于训练数据集的损失函数，优化的目标在于降低训练误差。而深度学习的目标在于降低泛化误差。为了降低泛化误差，除了使用优化算法降低训练误差以外，还需注意应对过拟合。

优化在深度学习中的挑战

%matplotlib inline
import sys
sys.path.append(r'.')
import d2lzh
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import torch

• 局部最小值与全局最小值

给定函数：

[f(x)=x*cos(pi x) ]

def f(x):
    return x*np.cos(np.pi*x)

d2lzh.set_figsize((4.5,2.5))
x=np.arange(-1,2,0.1)
fig=plt.plot(x,f(x))
fig[0].axes.annotate('local minimum',xy=(-0.3,-0.25),xytext=(-0.77,-1.0),arrowprops=dict(arrowstyle='->'))
fig[0].axes.annotate('global ',xy=(1.1,-0.95),xytext=(0.6,0.8),arrowprops=dict(arrowstyle='->'))

Text(0.6, 0.8, 'global ')

鞍点

梯度接近或变成0可能是由于在局部最优解附件造成的。事实上，另一种可能性是当前解在鞍点附件。

x=np.arange(-2,2,0.1)
fig=plt.plot(x,x**3)
fig[0].axes.annotate('saddle point',xy=(0,-0.1),xytext=(0,-5.5),arrowprops=dict(arrowstyle='->'))

Text(0, -5.5, 'saddle point')

x,y=np.mgrid[-1:1:31j,-1:1:31j]
z=x**2-y**2
ax=plt.figure().add_subplot(111,projection='3d')
ax.plot_wireframe(x,y,z,**{'rstride':2,'cstride':2})
ax.plot([0],[0],[0],'rx')
ticks=[-1,0,1]
plt.xticks(ticks)
plt.yticks(ticks)
ax.set_zticks(ticks)
plt.xlabel('x')
plt.ylabel('y')

Text(0.5, 0, 'y')

梯度下降和随机梯度下降

一维梯度下降

一维函数的梯度是一个标量，也称为导数，即自变量增长的速率，导数的相反数即自变量减小的速率，因此，顺着导数相反数的方向，函数值就会减小。

假设连续可导的函数f，输入，输出都是标量，给定绝对值足够小的数(epsilon)，根据泰勒展开公式，得到如下近似：
(f(x+epsilon) approx f(x)+epsilon f'(x)).

接着，找到一个常数(eta>0)，使得(|eta f'(x)|)足够小，那么可以将(epsilon)替换为(-eta f'(x))并得到：

(f(x-eta f'(x)) approx f(x)-eta f'(x)^2)，如果 f'(x)不为0，那么(eta f'(x)^2>0)，所以：

(f(x-eta f'(x))< f(x))，这意味着，如果通过(x leftarrow x-eta f'(x))来迭代x,f(x)的值可能会降低。

以$$f(x)=x^2$$为例，看一看梯度下降是如何工作的,使用x=10为初始值，$$eta =0.2$$，使用梯度下降迭代10次。

def gd(eta):
    x=10
    results=[x]
    for i in range(10):
        x-=eta*2*x
        results.append(x)
    print('epoch 10, x:',x)
    return results

res=gd(0.2);res

epoch 10, x: 0.06046617599999997





[10,
 6.0,
 3.5999999999999996,
 2.1599999999999997,
 1.2959999999999998,
 0.7775999999999998,
 0.46655999999999986,
 0.2799359999999999,
 0.16796159999999993,
 0.10077695999999996,
 0.06046617599999997]

def show_trace(res):
    n=max(abs(min(res)),abs(max(res)),10)
    f_line=np.arange(-n,n,0.1)
    d2lzh.set_figsize()
    plt.plot(f_line,[x*x for x in f_line])
    plt.plot(res,[x*x for x in res],'-o')
    plt.xlabel('x')
    plt.ylabel('y')

show_trace(res)

学习率

上述梯度下降算法中的正数$$eta$$通常叫做学习率，这是一个超参数，需要人工设定。如果使用过小的学习率，会导致x更新缓慢从而需要更多的迭代才能得到较好的解。

show_trace(gd(0.05))

epoch 10, x: 3.4867844009999995

可以发现使用过小的学习率，最终的值依然与最优解存在较大的偏差。

show_trace(gd(1.1))

epoch 10, x: 61.917364224000096

如果使用过大的学习率，(|eta f'(x)|)可能会过大从而使得前面的一阶泰勒展开公式不再成立：这时我们无法保证迭代x会降低f(x)的值。

多维梯度下降

可以通过$$x leftarrow x-eta abla f(x)$$来迭代，不断降低目标函数的f值。

下面我们构造一个输入为二维向量(x=[x_1,x_2]^T)和输出为标量的目标函数(f(x)=x_1^2+x_2^2)。那么梯度( abla f(x) =[2x_1,4x_2]^T)，观察梯度下降从初始位置[-5,-2]开始对自变量x的迭代轨迹。先定义两个辅助函数，第一个函数使用给定的自变量更新函数，从初始位置[-5,2]开始迭代自变量20次，第二个函数对自变量x的迭代轨迹进行可视化。

def train_2d(trainer):
    """
    trainer:梯度下降，迭代x,y
    """
    x1,x2,s1,s2=-5,-2,0,0
    results=[(x1,x2)]
    for i in range(20):
        x1,x2,s1,s2=trainer(x1,x2,s1,s2)
        results.append((x1,x2))
    print('epoch %d, x1 %f, x2 %f'%(i+1,x1,x2))
    return results

def show_trace_2d(f,results):
    plt.plot(*zip(*results),'-o',color='orange')
    x1,x2=np.meshgrid(np.arange(-5.5,1.0,0.1),np.arange(-3,1,0.1))
    plt.contour(x1,x2,f(x1,x2))
    plt.xlabel('x1')
    plt.xlabel('y1')

eta=0.1
def f_2d(x1,x2):
    return x1**2+2*x2**2
def gd_2d(x1,x2,s1,s2):
    return (x1-eta*2*x1,x2-eta*4*x2,0,0)

show_trace_2d(f_2d,train_2d(gd_2d))

epoch 20, x1 -0.057646, x2 -0.000073

随机梯度下降

在深度学习中，目标函数通常是训练数据集中有关各个样本的损失函数的平均，设(f_i(x))是有关索引i的训练数据样本的损失函数，n是训练数据样本数，x是模型的参数向量，那么目标函数的定义为：

(f(x)=frac{1}{n}Sigma_{i=1}^n f_i(x))

目标函数在x处的梯度为：

( abla f(x)=frac {1}{n} Sigma_{i=1}^{n} abla f_i(x))

如果使用梯度下降，每次自变量迭代的计算开销是O(n)，它随着n线性增长。因此，当训练数据样本数很大时，梯度下降每次迭代的计算开销很高。

随机梯度下降（stochastic gradient descent, SGD)减少了每次迭代的计算开销。在随机梯度下降的每次迭代中，我们随机均匀采样的一个样本索引(i in {1,2,....n})，并计算梯度( abla f_i(x))来迭代x:

(x leftarrow x-eta abla f_i(x))

可以看到每次迭代的计算开销从梯度下降的O(n)降到了O(1)。值得强调的是，随机梯度( abla f(x))是对梯度( abla f(x))的无偏估计：

(E_i abla f_i(x)=frac {1} {n} Sigma_{i=1}^n abla f_i(x)= abla f(x))

下面，通过在梯度中添加均值为0的随机噪声来模拟随机梯度下降，依此来比较它与梯度下降的区别：

def sgd_2d(x1,x2,s1,s2):
    return (x1-eta*(2*x1+np.random.normal(0.1)),
           x2-eta*(4*x2+np.random.normal(0.1)),0,0)

show_trace_2d(f_2d,train_2d(sgd_2d))

epoch 20, x1 -0.105398, x2 0.045164

可以看到，随机梯度下降中自变量的迭代轨迹相对于梯度下降中的来说更曲折。这是由于试验所添加的噪声使得模拟的随机梯度的准确度下降。在实际中，这些噪声通常指训练数据中的无意义的干扰。

小批量随机梯度下降

在每一次迭代中，梯度下降使用整个训练数据集来计算梯度，因此它有时也被称为批量梯度下降，而随机梯度下降在每次迭代中只随机采样一个样本计算梯度。我们还可以在每轮迭代中随机均匀采样多个样本来组成一个小批量，然后用这个小批量来计算梯度。下面来描述小批量随机梯度下降：

设目标函数f(x):(R^d ightarrow R)。在迭代开始前的时间步设为0.该时间步的自变量记为(x_0 in R^d)，通常由随机初始化得到。在接下来的每一个时间步t>0中，小批量随机梯度下降随机均匀采用一个由训练数据样本索引组成的小批量(B_t)。可以通过重复采样或者不重复采样得到一个小批量中的各个样本。前者允许同一小批量中出现重复的样本，后者则不允许如此，且更常见。对于这两种，都可以使用：

(g_t leftarrow abla f_{B_t}(x_{t-1})=frac {1}{|B|} Sigma_{i in B_t} abla f_i(x_{t-1}))

来计算时间步t的小批量(B_t)上目标函数位于(x_{t-1})处的梯度(g_t)，这里(|B|)表示批量大小，即小批量中样本的个数，是一个超参数。同随机梯度一样，重复采样所得的小批量随机梯度(g_t)也是对梯度( abla f(x_{t-1})的无偏估计。给定学习率(eta_t)（正数），小批量随机梯度下降对自变量的迭代如下：

(x_t leftarrow x_{t-1}-eta_t g_t)

基于随机采样得到的梯度的方差在迭代过程中无法减小，因此在实际中，小批量随机梯度下降的学习率可以在迭代过程中自我衰减，例如(eta_t=eta t^{alpha})（通常(alpha)取-1，或-0.5），(eta_t=eta alpha^t)((alpha)取0.95）或者每次迭代若干次后将学习率衰减一次。如此一来，学习率和小批量随机梯度乘积的方差会减小。而梯度下降在迭代过程中一直使用目标函数的真实梯度，无需自我衰减学习率。

小批量随机梯度下降中每次迭代的计算开销为O(B)。当批量大小为1时，该算法即为随机梯度下降，当批量大小等于训练数据样本时候，该算法即为梯度下降。当批量较小时，每次迭代中使用的样本少，这会导致并行处理和内存使用效率低。这使得在计算同样数目样本的情况下比使用更大批量时所花的时间更多。当批量大时，每个小批量梯度可能含有更多冗杂信息。为了得到较好的解，批量交大时比批量较小时所需计算的所需计算的样本数目可能更多，例如增大迭代周期数。

读取数据

def get_data_ch7():
    data=np.genfromtxt('Datasets/airfoil_self_noise/airfoil_self_noise.dat',delimiter='	')
    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

torch.Size([1500, 5])

labels.shape

torch.Size([1500])

def sgd(params,states,hyperparams):
    for p in params:
        p.data-=hyperparams['lr']*p.grad.data

import time
def train_ch7(optimizer_fn,states,hyperparams,features,labels,batch_size=10,num_epochs=2):
    net,loss=d2lzh.linreg,d2lzh.squared_loss
    w=torch.nn.Parameter(torch.tensor(np.random.normal(0,0.01,size=(features.shape[1],1)),dtype=torch.float32),requires_grad=True)
    b=torch.nn.Parameter(torch.zeros(1,dtype=torch.float32),requires_grad=True)

    def eval_loss():
        return loss(net(features,w,b),labels).mean().item()
    ls=[eval_loss()] #记录当下loss值
    data_iter=torch.utils.data.DataLoader(torch.utils.data.TensorDataset(features,labels),batch_size,shuffle=True) #生成data_iterator
    for _ in range(num_epochs):
        start=time.time()
        for batch_i,(x,y) in enumerate(data_iter):
            l=loss(net(x,w,b),y).mean()
            if w.grad is not None:
                w.grad.data.zero_()
                b.grad.data.zero_()
            l.backward()
            optimizer_fn([w,b],states,hyperparams)
            if (batch_i+1)*batch_size % 100 ==0:
                ls.append(eval_loss())
    print('loss:%f, %f sec per epoch'% (ls[-1],time.time()-start))
    d2lzh.set_figsize()
    plt.plot(np.linspace(0,num_epochs,len(ls)),ls)
    plt.xlabel('epoch')
    plt.ylabel('loss')

def train_sgd(lr,batch_size,num_epochs=2):
    train_ch7(sgd,None,{'lr':lr},features,labels,batch_size,num_epochs)

train_sgd(1,1500,6)

loss:0.247343, 0.014991 sec per epoch

train_sgd(0.005,1)

loss:0.242825, 0.289010 sec per epoch

train_sgd(0.05,10)

loss:0.243564, 0.040000 sec per epoch

简洁实现

在Pytorch里可以通过创建optimizer实例来调用优化算法，这让实现更简洁。下面实现一个通用的训练函数，它通过优化算法的函数optimizer_fn和超参数optimizer_hyperparams来创建optimizer实例。

from torch import nn
def train_pytorch_ch7(optimizer_fn,optimizer_hyperparams,features,labels,batch_size=10,num_epochs=2):
    net=nn.Sequential(nn.Linear(features.shape[-1],1))
    loss=nn.MSELoss()
    optimizer=optimizer_fn(net.parameters(),**optimizer_hyperparams)

    def eval_loss():
        return loss(net(features).view(-1),labels).item()/2
    
    ls=[eval_loss()]
    data_iter=torch.utils.data.DataLoader(torch.utils.data.TensorDataset(features,labels),batch_size,shuffle=True)
    for _ in range(num_epochs):
        start=time.time()
        for batch_i,(x,y) in enumerate(data_iter):
            l=loss(net(x).view(-1),y)/2 

            optimizer.zero_grad()
            l.backward()
            optimizer.step()

            if (batch_i+1)*batch_size%100==0:
                ls.append(eval_loss())
    print('loss:%f,%f sec per epoch'%(ls[-1],time.time()-start))
    d2lzh.set_figsize()
    plt.plot(np.linspace(0,num_epochs,len(ls)),ls)
    plt.xlabel('epoch')
    plt.ylabel('loss')

from torch import optim

train_pytorch_ch7(optim.SGD,{'lr':0.05},features,labels,10)

loss:0.245635,0.041000 sec per epoch

动量法

目标函数有关自变量的梯度代表了目标函数在自变量当前位置下降最快的方向，因此，梯度下降又称为最陡下降。在每次迭代中，梯度下降根据自变量当前位置，沿着当前位置的梯度更新自变量。然而，如果自变量的迭代方向仅仅取决于自变量的当前位置，这可能会带来一些问题。

现在考虑一个输入和输出分别为二维向量(x=[x_1,x_2]^T)和标量的目标函数(f(x)=0.1x_1^2+2x_2^2)。下面实现基于这个目标函数的梯度下降，并演示使用学习率为0.4时自变量的迭代轨迹。

eta=0.4
def f_2d(x1,x2):
    return 0.1*x1**2+2*x2**2

def gd_2d(x1,x2,s1,s2):
    return (x1-eta*2*x1,x2-eta*4*x2,0,0)

show_trace_2d(f_2d,train_2d(gd_2d))

epoch 20, x1 -0.000000, x2 -0.000073

可以看到，同一个位置，目标函数在竖直方向（x2)比水平方向（x1)的斜率的绝对值更大。因此，给定学习率，梯度下降迭代自变量时会使得自变量在竖直方向比在水平方向移动的幅度更大。因此，我们需要一个较小的学习率从而避免自变量在竖直方向上越过目标函数的最优解。然而，这会造成自变量在水平方向朝最优解移动变慢。

eta=0.6
show_trace_2d(f_2d,train_2d(gd_2d))

epoch 20, x1 -0.000000, x2 -1673.365109

加大学习率，可以发现，此时自变量在竖直方向上越过最优解，并逐渐发散。

动量法

动量法的提出就是为了解决梯度下降的上述问题。设时间步t的自变量为(x_t)，学习率为(eta_t)，在时间步0，动量法创建速度变量(v_0)，并将其元素初始化为0，在时间步t>0，动量法对每次迭代的步骤修改如下：

(v_t leftarrow gamma v_{t-1}+eta_t g_t)

(x_t leftarrow x_{t-1}-v_t)

动量参数0<=(gamma)<1,当为0时，动量法等价于小批量随机梯度下降。

def momentum_2d(x1,x2,v1,v2):
    v1=gamma*v1+eta*0.2*x1
    v2=gamma*v2+eta*4*x2
    return x1-v1,x2-v2,v1,v2

eta,gamma=0.4,0.5
show_trace_2d(f_2d,train_2d(momentum_2d))

epoch 20, x1 -0.062843, x2 0.001202

eta=0.6
show_trace_2d(f_2d,train_2d(momentum_2d))

epoch 20, x1 0.007188, x2 0.002553

可以看到用较大的学习率0.6，也不再发散了。

指数加权移动平均

为了从数学上理解动量法，先学习下指数加权移动平均。给定超参数(0<=gamma<1)，当前时间步t的变量(y_t)是上一时间步t-1的变量和当前时间步另一变量(x_t)的线性组合：

(y_t=gamma y_{t-1}+(1-gamma)x_t)

对(y_t)展开：

(y_t=(1-gamma)x_t+gamma y_{t-1})

(=(1-gamma)x_t+(1-gamma)gamma x_{t-1}+(1-gamma) gamma^2 x_{t-1}+gamma^3 y_{t-3} ...)

令(n=1/(1-gamma))，那么有((1-1/n)^n=gamma^{1/(1-gamma)}),考虑到
(lim_{n ightarrow oo}(1-1/n)^n=exp(-1) approx 0.3679)

所以当(gamma ightarrow 1)，(gamma^{1/(1-gamma)}=exp(-1))，如 当(gamma =0.95)，带入右式，得(0.95^20 approx exp(-1))，如果我们认为exp(-1)是一个比较小的数，那么我们就可以忽略所有含有(gamma^{1/(1-gamma)})或比它更高阶的项，即：

(y_t approx 0.05 Sigma_{i=1}^{19} 0.95^i x_{t-i})

因此，在实际中，我们常常将(y_t)看作是对最近(1/(1-gamma))个时间步的(x_t)值的加权平均。例如当(gamma)为0.95时，就是对最近20个时间步的(x_t)的值加权平均；当为0.9，就是对最近的10个进行加权平均。

x=np.random.randn(20)
def gety(yt_1,gamma,xt):
    y_t=gamma*yt_1+(1-gamma)*xt
    return y_t
y=[0]
for i in x:
    y_i=gety(y[-1],0.9,i)
    y.append(y_i)
plt.plot(range(20),x)
plt.plot(range(21),y,'-o')

[<matplotlib.lines.Line2D at 0x1a23de5c1c8>]

由指数加权移动平均理解动量法

现在我们对动量法的速度分量做变形：

(v_t leftarrow gamma v_{t-1}+(1-gamma)(frac {eta_t} {1-gamma} g_t))

由指数加权移动平均的型式可得，速度变量实际是对序列({ eta_{t-i}g_{t-i}/(1-gamma):i=0,...1/(1-gamma)-1 })做了指数加权移动平均。也就是说，相比于小批量随机梯度下降，动量法在每个时间步的自变量更新量近似于将最近(1/(1-gamma))个时间步的普通更新量（学习率与梯度的乘积）做了指数加权移动平均后除以(1-gamma)。所以，在动量法中，自变量在各个方向上的移动幅度不仅取决于当前梯度，还取决于过去各个梯度在各个方向上是否一致。在之前示例的优化问题中，所有梯度在水平方向上为正（向右），而在竖直方向上时正时负，因此，我们就可以使用较大的学习率，从而使自变量向最优解更快移动。

从零开始实现

相对于小批量随机梯度下降，动量法需要对每一个自变量维护一个同它一样形状的速度变量，且超参数里多了动量超参数。实现中，我们将速度变量用更广义的状态变量states表示。

features,labels=d2lzh.get_data_ch7()

def init_momentum_states():
    v_w=torch.zeros((features.shape[-1],1),dtype=torch.float32)
    v_b=torch.zeros(1,dtype=torch.float32)
    return v_w,v_b

def sgd_momentum(params,states,hyperparams):
    for p,v in zip(params,states):
        v.data=hyperparams['momentum']*v.data+hyperparams['lr']*p.grad.data
        p.data-=v.data

d2lzh.train_ch7(sgd_momentum,init_momentum_states(),{'lr':0.02,'momentum':0.5},features,labels)

loss:0.242594, 0.042000 sec per epoch

d2lzh.train_ch7(sgd_momentum,init_momentum_states(),{'lr':0.02,'momentum':0.9},features,labels)

loss:0.253085, 0.046000 sec per epoch

目标函数值在后期迭代过程中不够平滑，直觉上，10倍小批量梯度比2倍小批量梯度打了5倍，可试着将学习率减小到原来的1/5.

d2lzh.train_ch7(sgd_momentum,init_momentum_states(),{'lr':0.004,'momentum':0.9},features,labels)

loss:0.243076, 0.045998 sec per epoch

简洁实现

d2lzh.train_pytorch_ch7(torch.optim.SGD,{'lr':0.004,'momentum':0.9},features,labels)

loss:0.245053,0.043995 sec per epoch

ADAGRAD算法

在之前的优化算法中，目标函数自变量的每一个元素在相同时间步都使用同一个学习率来自我迭代。举个例子，假设目标函数为f,自变量为一个二维向量([x_1,x_2]^T)，该向量中每一个元素在迭代时都使用相同的学习率。例如，在学习率为(eta)的梯度下降中，元素(x_1)，(x_2)都使用相同的学习率来自我迭代：

(x_1 leftarrow x_1-eta frac {partial f} {partial x_1}, x_2 leftarrow x_2-eta frac {partial f} {partial x_2})

动量法中，当(x_1 x_2)的梯度值由较大差别，需要选择足够小的学习率使得自变量在梯度值较大的维度上不发散。但这样会导致自变量在梯度值较小的维度上迭代过慢。动量法依赖于指数加权移动平均使得自变量的更新方向更加一致，从而降低发散的可能。而 AdaGrad算法，根据自变量在每个维度的梯度值的大小来调整各个维度上的学习率，从而避免统一的学习率难以适应所有维度的问题. 相比较于动量法，AdaGrad的学习率显然对所有时间步的梯度进行了考虑（平方累加），而动量法仅仅考虑了(1/(1-gamma))个时间步的梯度的加权平均

算法

AdaGrad算法会使用一个小批量随机梯度(g_t)按元素平方的累加变量(s_t)。在时间步0，AdaGrad将(s_0)中每个元素初始化为0，在时间步t，首先将小批量随机梯度(g_t)按元素平方后累加到变量(s_t)：

(s_t leftarrow s_{t-1}+g_t odot g_t)

其中，(odot)是按元素相乘。接着，将目标函数自变量中每个元素的学习率通过按元素运算重新调整：

(x_t leftarrow x_{t-1}- frac {eta} {sqrt {s_t +epsilon}} odot g_t)

(eta)：学习率。(epsilon):为维持数值稳定性而添加的常数。这里开方，除法，乘法的运算都是按元素运算的。这些按元素运算使得目标函数自变量中每个元素都分别拥有自己的学习率。

特点

小批量随机梯度按元素平方的累加变量(s_t)出现在学习率的分母项中。因此，如果目标函数有关自变量中的某个元素的偏导数一直比较大，那么该元素的学习率将下降较快,因为学习率不断在除以较大的值；反之，如果目标函数有关自变量中某个元素的偏导数一直较小，那么该元素的学习率将下降较慢，因为学习率不断除以较小值。然后由于(s_t)一直在累加按元素平方的梯度，自变量中每个元素的学习率在迭代过程中一直在降低（或不变）。所以，当学习率在迭代早期降得较快且当前解依然不佳时，AdaGrad算法在迭代后期由于学习率过小，可能难以找到一个有用的解。

仍然以目标函数(f(x)=0.1x_1^2+2x_2^2)为例子，观察AdaGrad算法对自变量的迭代轨迹。学习率依然取0.4

import math

def adagrad_2d(x1,x2,s1,s2):
    g1,g2,eps=0.2*x1,4*x2,1e-6
#     print(g1,g2)
    s1+=g1**2
    s2+=g2**2
    x1-=eta/math.sqrt(s1+eps)*g1
    x2-=eta/math.sqrt(s2+eps)*g2
    return x1,x2,s1,s2

def f_2d(x1,x2):
    return 0.1*x1**2+2*x2**2

eta=0.4
d2lzh.show_trace_2d(f_2d,d2lzh.train_2d(adagrad_2d))

epoch 20, x1 -2.382563, x2 -0.158591

eta=2
d2lzh.show_trace_2d(f_2d,d2lzh.train_2d(adagrad_2d))

epoch 20, x1 -0.002295, x2 -0.000000

将学习率增大到2，可以发现，自变量较为迅速的逼近了最优解。

从零开始实现

features,labels=d2lzh.get_data_ch7()

def init_adagrad_states():
    s_w=torch.zeros((features.shape[-1],1),dtype=torch.float32)
    s_b=torch.zeros(1,dtype=torch.float32)
    return s_w,s_b

def adagrad(params,states,hyperparams):
    eps=1e-6
    for p,s in zip(params,states):
        s.data+=(p.grad.data**2)
        p.data-=hyperparams['lr']*p.grad.data/torch.sqrt(s+eps)

d2lzh.train_ch7(adagrad,init_adagrad_states(),{'lr':0.1},features,labels)

loss:0.243944, 0.048000 sec per epoch

简洁实现

d2lzh.train_pytorch_ch7(torch.optim.Adagrad,{'lr':0.1},features,labels)

loss:0.242530,0.047001 sec per epoch

• 小结

（1）AdaGrad算法在迭代过程中不断调整学习率，并让目标函数自变量中每个元素都分别拥有自己的学习率。

（2）使用AdaGrad算法时，自变量中每个元素的学习率在迭代过程中一直在降低（或不变）

RMSPROP算法

AdaGrad算法，因为调整学习率时，分母上的变量(s_t)一直在累加按元素平方的小批量随机梯度，所以目标函数自变量每个元素的学习率在迭代过程中一直在降低（或不变）。因此，当学习率在迭代早期降得快且解依然不佳时，AdaGrad算法在迭代后期由于学习率过小，可能较难找到一个有用的解，为解决这个问题，RMSPROP算法对AdaGrad算法做了一点小小的修改。

算法

不同于AdaGrad算法里状态状态变量(s_t)是截至时间步t所有小批量随机梯度(g_t)按元素平方和，RMSPROP算法将这些梯度按元素平方做指数加权移动平均。具体来讲，给定超参数(0<=gamma<1)，RMSPROP算法在时间步t>0计算：

(s_t leftarrow gamma s_{t-1}+(1-gamma)g_t odot g_t)

和AdaGrad算法一样，RMSPRPOP算法将目标函数自变量中每个元素的学习率通过按元素运算重新调整，然后更新自变量：

(x_t leftarrow x_{t-1}- frac {eta} { sqrt {s_t+epsilon}} odot g_j)

其中，(eta)是学习率，(epsilon)是为了维持数值稳定性而添加的常数，如(10^{-6})。因为RMSPROP算法是状态变量对平方项的指数加权移动平均，所以可以看作是(1/(1-gamma）)个时间步的小批量随机梯度平方项的加权平均。如此一来，自变量每个元素的学习率在迭代过程中就不再一直降低（或不变）

先观察(f(x)=0.1x_1^2+2x_2^2)中自变量的迭代轨迹。在相同的学习率下，RMSPROP可以更快逼近最优解。

def rmsprop_2d(x1,x2,s1,s2):
    g1,g2,eps=0.2*x1,4*x2,1e-6
    s1=gamma*s1+(1-gamma)*g1**2
    s2=gamma*s2+(1-gamma)*g2**2
    x1-=eta/math.sqrt(s1+eps)*g1
    x2-=eta/math.sqrt(s2+eps)*g2
    return x1,x2,s1,s2

def f_2d(x1,x2):
    return 0.1*x1**2+2*x2**2

eta,gamma=0.4,0.9
d2lzh.show_trace_2d(f_2d,d2lzh.train_2d(rmsprop_2d))

epoch 20, x1 -0.010599, x2 0.000000

从零开始实现

features,labels=d2lzh.get_data_ch7()
def init_rmsprop_states():
    s_w=torch.zeros((features.shape[-1],1),dtype=torch.float32)
    s_b=torch.zeros(1,dtype=torch.float32)
    return s_w,s_b

def rmsprop(params,states,hyperparams):
    gamma,eps=hyperparams['gamma'],1e-6
    for p,s in zip(params,states):
        s.data=gamma*s.data+(1-gamma)*(p.grad.data)**2
        p.data-=hyperparams['lr']*p.grad.data/torch.sqrt(s+eps)

d2lzh.train_ch7(rmsprop,init_rmsprop_states(),{'lr':0.01,'gamma':0.9},features,labels)

loss:0.246476, 0.051000 sec per epoch

简洁实现

d2lzh.train_pytorch_ch7(torch.optim.RMSprop,{'lr':0.01,'alpha':0.9},features,labels)

loss:0.243686,0.049001 sec per epoch

ADADELTA 算法

除了RMSProp算法外，另一个常用优化算法AdaDelta算法也针对AdaGrad算法在迭代后期可能较难找到有用解的问题做了改进。AdaDeleta算法没有学习率这一个超参数。

算法

AdaDelta 算法也像RMSProp算法一样，使用小批量随机梯度(g_t)按元素平方的指数加权移动平均变量(s_t)。在时间步0，它的所有元素为0，给定超参数(0 le ho < 1)，在时间步t>0：

(s_t leftarrow ho s_{t-1}+(1- ho)g_t odot g_t)

与RMSPRrop不同的是，AdaDelta算法还维护一个额外的状态变量(Delta x_t)，起元素同样在时间步0时被初始化为0.我们使用(Delta x_{t-1})来计算自变量的变化量：

(g_t' leftarrow sqrt {frac { Delta x_{t-1}+epsilon} {s_t +epsilon} } odot g_t)

接着更新自变量：

(x_t leftarrow x_{t-1}-g_t').

最后，使用(Delta x_t)来记录自变量(g_t')按元素平方的指数加权移动平均：

(Delta x_t leftarrow ho Delta x_{t-1} +(1- ho)g_t' odot g_t')

可以看到，如不考虑(epsilon)的影响，AdaDelta算法与RMSProp算法的不同之处在于用(sqrt {Delta x_{t-1} })来代替学习率(eta)。

从零开始

features,labels=d2lzh.get_data_ch7()

def init_adadelta_states():
    s_w,s_b=torch.zeros((features.shape[1],1),dtype=torch.float32),torch.zeros(1,dtype=torch.float32)
    delta_w,delta_b=torch.zeros(features.shape[1],1,dtype=torch.float32),torch.zeros(1,dtype=torch.float32)
    return ((s_w,delta_w),(s_b,delta_b))

def adadelta(params,states,hyperparams):
    rho,eps=hyperparams['rho'],1e-5
    for p,(s,delta) in zip(params,states):
        # s,delta 分别为(s_w,s_b)，（delta_w,delta_b)
        s=rho*s+(1-rho)*(p.grad.data**2)
        g=p.grad.data*torch.sqrt((delta+eps)/(s+eps))
        p.data-=g
        delta=rho*delta+(1-rho)*g*g

d2lzh.train_ch7(adadelta,init_adadelta_states(),{'rho':0.9},features,labels)

loss:0.244290, 0.053967 sec per epoch

简洁实现

d2lzh.train_pytorch_ch7(torch.optim.Adadelta,{'rho':0.9},features,labels)

loss:0.284449,0.062043 sec per epoch

Adam算法

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

算法

Adam 算法使用了动量变量(v_t)和RMSProp算法中的(s_t)，并在时间步0将它们初始化为0.给定超参数(0 le eta_1 <1)（建议为0.9），时间步t的动量变量(v_t)即小批量随机梯度(g_t)的指数加权移动平均：

(v_t leftarrow eta_1 v_{t-1} +(1-eta_1)g_t)

(v_t=eta_1 v_{t-1} +(1-eta_1)g_t)

(=eta_1^2 v_{t-2}+eta_1 (1-eta_1)g_{t-1} +(1-eta_1)g_t)

(=eta_1^3 v_{t-3}+eta_1^2(1-eta_1)g_{t-2} +eta_1 (1-eta_1)g_{t-1} +(1-eta_1)g_t)

(...)

如此得到：(v_t=v_0+(1-eta_1)Sigma_{t=1}^t eta_1^{t-i}g_i)

和RMSProp算法一样，给定超参数(0 le eta_2 s_{t-1}+(1-eta_2) g_t odot g_t).

(s_t leftarrow eta_2 s_{t-1} + (1-eta_2)g_t odot g_t)

由于将(v_0)，(s_0)中的元素均初始化为0,故有(v_t=(1-eta_1)Sigma_{t=1}^t eta_1^{t-i}g_i)

将过去各个时间步小批量随机梯度权值相加，等比数列求和，得到((1-eta_1)Sigma_{i=1}^{t}eta_1^{t-i}=eta_1^{t-1}*(1-eta_1^{-t})/(1-eta_1^{-1})*(1-eta_1))

(=(eta_1^{t-1}-eta_1^{-1})/(eta_1-1)*eta_1*(1-eta_1)=1-eta_1^t)

需要注意的是，当t较小时，过去各时间不小批量随机梯度权值之和会较小，比如，当(eta_1)为0.9时，(v_1=0.1g_1)，为消除这样的影响，对任意时间步t，可将(v_t)再除以一个(1-eta_1^t)，从而使得过去各个时间步小批量随机梯度权重之和为1.这叫做偏差修正。

(hat{v_t} leftarrow frac {v_t} {1-eta_1^t})

(hat{s_t} leftarrow frac {s_t} {1-eta_2^t})

接下来，Adam算法使用偏差修正后的变量(hat{v_t},hat{s_t})，将模型参数中的每个元素的学习率通过按元素运算重新调整：

(g_t' leftarrow frac {eta hat {v_t}} {sqrt {hat {s_t}}+epsilon})

(x_t leftarrow x_{t-1}-g_t')

从零开始

features,labels=d2lzh.get_data_ch7()

def init_adam_states():
    v_w,v_b=torch.zeros((features.shape[-1],1),dtype=torch.float32),torch.zeros(1,dtype=torch.float32)
    s_w,s_b=torch.zeros((features.shape[-1],1),dtype=torch.float32),torch.zeros(1,dtype=torch.float32)
    return ((v_w,s_w),(v_b,s_b))

def adam(params,states,hyperparams):
    beta1,beta2,eps=0.9,0.99,1e-6
    for p,(v,s) in zip(params,states):
        v=beta1*v+(1-beta1)*p.grad.data
        s=beta2*s+(1-beta2)*p.grad.data**2
        v_bias_corr=v/(1-beta1**hyperparams['t'])
        s_bias_corr=s/(1-beta2**hyperparams['t'])
        p.data-=hyperparams['lr']*v_bias_corr/(torch.sqrt(s_bias_corr)+eps)
        hyperparams['t']+=1

d2lzh.train_ch7(adam,init_adam_states(),{'lr':0.01,'t':1},features,labels)

loss:0.249869, 0.059039 sec per epoch

简洁实现

d2lzh.train_pytorch_ch7(torch.optim.Adam,{'lr':0.01},features,labels)

loss:0.246421,0.057026 sec per epoch

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