凸优化

尽管优化方法可以最小化深度学习中的损失函数值，但本质上优化方法达到的目标与深度学习的目标并不相同。

• 优化方法目标：训练集损失函数值
• 深度学习目标：测试集损失函数值（泛化性）

优化在深度学习中的挑战

• 1. 局部最小值

eg:以该函数为例,局部最优值

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

d2l.set_figsize((4.5, 2.5))
x = np.arange(-1.0, 2.0, 0.1)
fig,  = d2l.plt.plot(x, f(x))
fig.axes.annotate('local minimum', xy=(-0.3, -0.25), xytext=(-0.77, -1.0),
                  arrowprops=dict(arrowstyle='->'))
fig.axes.annotate('global minimum', xy=(1.1, -0.95), xytext=(0.6, 0.8),
                  arrowprops=dict(arrowstyle='->'))
d2l.plt.xlabel('x')
d2l.plt.ylabel('f(x)');

• 2. 鞍点

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

• 3. 梯度消失

x = np.arange(-2.0, 5.0, 0.01)
fig, = d2l.plt.plot(x, np.tanh(x))
d2l.plt.xlabel('x')
d2l.plt.ylabel('f(x)')
fig.axes.annotate('vanishing gradient', (4, 1), (2, 0.0) ,arrowprops=dict(arrowstyle='->'))

Jensen 不等式

性质

1. 无局部极小值
   证明：假设存在 是局部最小值，则存在全局最小值 , 使得 , 则对 :

2. 与凸集的关系
   对于凸函数 ，定义集合 ，则集合 为凸集
   证明：对于点 , 有 , 故
3. 二阶条件
   是凸函数

类似于运筹学的优化问题

梯度下降

从一维说起：
证明：沿梯度反方向移动自变量可以减小函数值

泰勒展开：

代入沿梯度方向的移动量 ：

def f(x):
# 原函数
    return x**2  # Objective function

def gradf(x):
# 导函数
    return 2 * x  # Its derivative

def gd(eta):
# 输入为学习率
    x = 10
    results = [x]
    for i in range(10):
        x -= eta * gradf(x)
        results.append(x)
    print('epoch 10, x:', x)
    return results
def show_trace(res):
# 展示下降的梯度
    n = max(abs(min(res)), abs(max(res)))
    f_line = np.arange(-n, n, 0.01)
    d2l.set_figsize((3.5, 2.5))
    d2l.plt.plot(f_line, [f(x) for x in f_line],'-')
    d2l.plt.plot(res, [f(x) for x in res],'-o')
    d2l.plt.xlabel('x')
    d2l.plt.ylabel('f(x)')
show_trace(res)
# 通过更改学习率，输出展示梯度下降的图像
show_trace(gd(theta))

theta=0.2

theta=0.05

theta=1.1

多维情况

2维的情况

def train_2d(trainer, steps=20):
    x1, x2 = -5, -2
    results = [(x1, x2)]
    for i in range(steps):
        x1, x2 = trainer(x1, x2)
        results.append((x1, x2))
    print('epoch %d, x1 %f, x2 %f' % (i + 1, x1, x2))
    return results
def show_trace_2d(f, results): 
    d2l.plt.plot(*zip(*results), '-o', color='#ff7f0e')
    x1, x2 = np.meshgrid(np.arange(-5.5, 1.0, 0.1), np.arange(-3.0, 1.0, 0.1))
    d2l.plt.contour(x1, x2, f(x1, x2), colors='#1f77b4')
    d2l.plt.xlabel('x1')
    d2l.plt.ylabel('x2')
eta = 0.1
def f_2d(x1, x2):  # 目标函数
    return x1 ** 2 + 2 * x2 ** 2

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

show_trace_2d(f_2d, train_2d(gd_2d))

eta = 0.1

自适应方法——牛顿法

在 处泰勒展开：

最小值点处满足: , 即我们希望 , 对上式关于 求导，忽略高阶无穷小，有：

应用：

c = 0.15 * np.pi

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

def gradf(x):
    return np.cos(c * x) - c * x * np.sin(c * x)

def hessf(x):
    return - 2 * c * np.sin(c * x) - x * c**2 * np.cos(c * x)
def newton(eta=1):
    x = 10
    results = [x]
    for i in range(10):
        x -= eta * gradf(x) / hessf(x)
        results.append(x)
    print('epoch 10, x:', x)
    return results
show_trace(newton())

预处理 （Heissan阵辅助梯度下降）

随机梯度下降(SGD)

参数更新
对于有 个样本对训练数据集，设 是第 个样本的损失函数, 则目标函数为:

其梯度为:

使用该梯度的一次更新的时间复杂度为

随机梯度下降更新公式 :

且有：

def f(x1, x2):
    return x1 ** 2 + 2 * x2 ** 2  # Objective

def gradf(x1, x2):
    return (2 * x1, 4 * x2)  # Gradient

def sgd(x1, x2):  # Simulate noisy gradient
    global lr  # Learning rate scheduler
    (g1, g2) = gradf(x1, x2)  # Compute gradient
    (g1, g2) = (g1 + np.random.normal(0.1), g2 + np.random.normal(0.1))
    eta_t = eta * lr()  # Learning rate at time t
    return (x1 - eta_t * g1, x2 - eta_t * g2)  # Update variables

eta = 0.1
lr = (lambda: 1)  # Constant learning rate
show_trace_2d(f, train_2d(sgd, steps=50))

另外，我们还可将学习率设置为动态的

def exponential():
    global ctr
    ctr += 1
    return math.exp(-0.1 * ctr)

ctr = 1
lr = exponential  # Set up learning rate
show_trace_2d(f, train_2d(sgd, steps=1000))

exponential

def polynomial():
    global ctr
    ctr += 1
    return (1 + 0.1 * ctr)**(-0.5)

ctr = 1
lr = polynomial  # Set up learning rate
show_trace_2d(f, train_2d(sgd, steps=50))

polynomial

简洁实现

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):
            # 除以2是为了和train_ch7保持一致, 因为squared_loss中除了2
            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))
    d2l.set_figsize()
    d2l.plt.plot(np.linspace(0, num_epochs, len(ls)), ls)
    d2l.plt.xlabel('epoch')
    d2l.plt.ylabel('loss')
train_pytorch_ch7(optim.SGD, {"lr": 0.05}, features, labels, 10)

image.png

小批量梯度下降法

小批量梯度下降，是对批量梯度下降以及随机梯度下降的一个折中办法。其思想是：每次迭代 使用 batch_size 个样本来对参数进行更新。
  这里我们假设 batchsize=10 ，样本数 m=1000 。
  伪代码形式为：

1

优点：
  （1）通过矩阵运算，每次在一个batch上优化神经网络参数并不会比单个数据慢太多。
  （2）每次使用一个batch可以大大减小收敛所需要的迭代次数，同时可以使收敛到的结果更加接近梯度下降的效果。(比如上例中的30W，设置batch_size=100时，需要迭代3000次，远小于SGD的30W次)
  （3）可实现并行化。
  缺点：
  （1）batch_size的不当选择可能会带来一些问题。

指数加权平均

平均值

指数加权平均本质上就是一种近似求平均的方法。

指数加权平均的优势

我们可以看到指数加权平均的求解过程实际上是一个递推的过程，那么这样就会有一个非常大的好处，每当我要求从0到某一时刻（n）的平均值的时候，我并不需要像普通求解平均值的作为，保留所有的时刻值，类和然后除以n。
而是只需要保留0-(n-1)时刻的平均值和n时刻的温度值即可。也就是每次只需要保留常数值，然后进行运算即可，这对于深度学习中的海量数据来说，是一个很好的减少内存和空间的做法。
梯度下降优化法经历了SGD→SGDM→NAG→AdaGrad→AdaDelta→Adam→Nadam这样的发展历程。之所以会不断地提出更加优化的方法，究其原因，是引入了动量（Momentum）这个概念。

数学推导：

我们可以对 展开：

Supp 上界

Approximate Average of Steps

令 ，那么 。因为

所以当 时，，如 。如果把 当作一个比较小的数，我们可以在近似中忽略所有含 和比 更高阶的系数的项。例如，当 时，

因此，在实际中，我们常常将 看作是对最近 个时间步的 值的加权平均。例如，当 时， 可以被看作对最近20个时间步的 值的加权平均；当 时， 可以看作是对最近10个时间步的 值的加权平均。而且，离当前时间步 越近的 值获得的权重越大（越接近1）。

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

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

Another version:

由指数加权移动平均的形式可得，速度变量 实际上对序列 做了指数加权移动平均。换句话说，相比于小批量随机梯度下降，动量法在每个时间步的自变量更新量近似于将前者对应的最近 个时间步的更新量做了指数加权移动平均后再除以 。所以，在动量法中，自变量在各个方向上的移动幅度不仅取决当前梯度，还取决于过去的各个梯度在各个方向上是否一致。在本节之前示例的优化问题中，所有梯度在水平方向上为正（向右），而在竖直方向上时正（向上）时负（向下）。这样，我们就可以使用较大的学习率，从而使自变量向最优解更快移动。

PS:

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

Code:

def get_data_ch7():  
    data = np.genfromtxt('/home/kesci/input/airfoil4755/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()

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
d2l.train_ch7(sgd_momentum, init_momentum_states(),
              {'lr': 0.02, 'momentum': 0.5}, features, labels)

优化算法的框架

优化算法伪代码

• 1. 最速梯度下降法
     回归函数及目标函数

以均方误差作为目标函数（损失函数），目的是使其值最小化，用于优化上式。

对目标函数求导

沿导数相反方向移动theta

• 2.随机梯度下降法
  SGD是最速梯度下降法的变种。
  使用最速梯度下降法，将进行N次迭代，直到目标函数收敛，或者到达某个既定的收敛界限。每次迭代都将对m个样本进行计算，计算量大。
  SGD每次迭代仅对一个样本计算梯度，直到收敛。

• 3.动量（Momentum）梯度下降法
  SGD方法的一个缺点是，其更新方向完全依赖于当前的batch，因而其更新十分不稳定。解决这一问题的一个简单的做法便是引入momentum。momentum即动量，它模拟的是物体运动时的惯性，即更新的时候在一定程度上保留之前更新的方向，同时利用当前batch的梯度微调最终的更新方向。这样一来，可以在一定程度上增加稳定性，从而学习地更快，并且还有一定摆脱局部最优的能力：

红色箭头显示了一个带动量的小批量梯度下降的方向。蓝色的点显示了每个步骤的梯度（关于当前的小批）的方向，很明显的就是每一次更新的时候动量的存在矫正了“偏”的参数更新方向

与标准的梯度下降法的关系：

1、当β=0时，就变成了标准的梯度下降法

2、在我们进行动量梯度下降算法的时候，由于使用了指数加权平均的方法，原来在纵轴方向上的上下波动，经过平均以后，接近于0，纵轴上的波动变得非常的小（因为其是取前1/(1−β) 次的梯度）；但在横轴方向上，所有的微分都指向横轴方向，因此其平均值仍然很大。最终实现红色线所示的梯度下降曲线。

3、一般情况下β的取值越大越好，因为越大可以取到更多次的过去的值，可以让下降更加光滑，但是β的取值不能够太，一般取到0.9
4、动量梯度下降法既可以与 batchgradient descent, mini-batch gradient descent or stochastic gradient descent联合使用
总结：动量梯度下降法是在梯度上面进行优化

• 4.AdaGrad

在之前介绍过的优化算法中，目标函数自变量的每一个元素在相同时间步都使用同一个学习率来自我迭代。举个例子，假设目标函数为，自变量为一个二维向量，该向量中每一个元素在迭代时都使用相同的学习率。例如，在学习率为的梯度下降中，元素和都使用相同的学习率来自我迭代：

在[“动量法”]一节里我们看到当和的梯度值有较大差别时，需要选择足够小的学习率使得自变量在梯度值较大的维度上不发散。但这样会导致自变量在梯度值较小的维度上迭代过慢。动量法依赖指数加权移动平均使得自变量的更新方向更加一致，从而降低发散的可能。它根据自变量在每个维度的梯度值的大小来调整各个维度上的学习率，从而避免统一的学习率难以适应所有维度的问题 。

Algorithm

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

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

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

Feature

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

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)

RMSProp

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

Algorithm

我们在“动量法”一节里介绍过指数加权移动平均。不同于AdaGrad算法里状态变量是截至时间步所有小批量随机梯度按元素平方和，RMSProp算法将这些梯度按元素平方做指数加权移动平均。具体来说，给定超参数计算

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

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

照例，让我们先观察RMSProp算法对目标函数中自变量的迭代轨迹。回忆在“AdaGrad算法”一节使用的学习率为0.4的AdaGrad算法，自变量在迭代后期的移动幅度较小。但在同样的学习率下，RMSProp算法可以更快逼近最优解。

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

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

alpha, beta = 0.4, 0.9
d2l.show_trace_2d(f_2d, d2l.train_2d(rmsprop_2d))

# 初始化算法
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['beta'], 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)

AdaDelta

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

Algorithm

AdaDelta算法也像RMSProp算法一样，使用了小批量随机梯度按元素平方的指数加权移动平均变量。在时间步0，它的所有元素被初始化为0。给定超参数，同RMSProp算法一样计算

与RMSProp算法不同的是，AdaDelta算法还维护一个额外的状态变量，其元素同样在时间步0时被初始化为0。我们使用来计算自变量的变化量：

其中是为了维持数值稳定性而添加的常数，如。接着更新自变量：

最后，我们使用来记录自变量变化量按元素平方的指数加权移动平均：

可以看到，如不考虑的影响，AdaDelta算法与RMSProp算法的不同之处在于使用来替代超参数。

PS:

AdaDelta算法需要对每个自变量维护两个状态变量，即和。按AdaDelta算法中的公式实现该算法。

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[:] = 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

Adam

Adam算法在RMSProp算法基础上对小批量随机梯度也做了指数加权移动平均 [1]。下面我们来介绍这个算法。

Algorithm

Adam算法使用了动量变量和RMSProp算法中小批量随机梯度按元素平方的指数加权移动平均变量，并在时间步0将它们中每个元素初始化为0。给定超参数（算法作者建议设为0.9），时间步的动量变量即小批量随机梯度的指数加权移动平均：

和RMSProp算法中一样，给定超参数（算法作者建议设为0.999），
将小批量随机梯度按元素平方后的项做指数加权移动平均得到：

由于我们将和中的元素都初始化为0，
在时间步我们得到。将过去各时间步小批量随机梯度的权值相加，得到 。需要注意的是，当较小时，过去各时间步小批量随机梯度权值之和会较小。例如，当时，。为了消除这样的影响，对于任意时间步，我们可以将再除以，从而使过去各时间步小批量随机梯度权值之和为1。这也叫作偏差修正。在Adam算法中，我们对变量和均作偏差修正：

接下来，Adam算法使用以上偏差修正后的变量和，将模型参数中每个元素的学习率通过按元素运算重新调整：

其中是学习率，是为了维持数值稳定性而添加的常数，如。和AdaGrad算法、RMSProp算法以及AdaDelta算法一样，目标函数自变量中每个元素都分别拥有自己的学习率。最后，使用迭代自变量：

PS:

我们按照Adam算法中的公式实现该算法。其中时间步通过hyperparams参数传入adam函数。

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.999, 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