Pytorch学习(六）—— 优化算法介绍

优化与深度学习

深度学习中的优化算法通常只考虑最小化目标函数，优化为深度学习提供了最小化目标函数的方法。由于优化算法的的目标函数通常是一个基于训练数据集的损失函数，优化的目标在于降低训练误差。而深度学习的目标在于降低泛化误差，为了降低泛化误差，除了使用优化算法降低训练误差外，还需要注意应对过拟合。

深度学习中绝大多数目标函数都很复杂，因此，很多优化问题并不存在解析解，而需要使用基于数值方法的优化算法找到近似解，即数值解。以下优化算法都是基于数值方法的算法，为了求得最小化目标函数的数值解，通过优化算法的有限次迭代模型参数来尽可能的降低损失函数的值。

深度学习模型的目标函数可能有若干局部最优值，当一个优化问题的数值解在局部最优解附近时，由于目标函数有关解的梯度接近或变成零，最终迭代求得得数值解可能只令目标函数局部最小化而非全局最小化。除此之外，另一种可能得情况是当前解在鞍点附近。在下图鞍点位置，目标函数在 x 轴方向上是局部最小值，但在 y 轴方向上是局部最大值。

梯度下降和随机梯度下降

假设连续可导的函数 $f$ 输入和输出都是标量，给定绝对值足够小的数 $\xi$，根据泰勒展开公式，我们得到以下近似: $f(x+\xi )\approx f(x)+\xi f^{{}^{\prime }}(x)$ ，这里 $f^{{}^{\prime }}(x)$ 是函数$f$ 在 $x$ 处的梯度，即导数。接下来，找到一个常数 $\eta >0$ ，使得 $|\eta f^{{}^{\prime }}(x)|$ 足够小，那么可以将 $\xi$ 替换为 $-\eta f^{{}^{\prime }}(x)$ 并得到:
                   $f(x-\eta f^{{}^{\prime }}(x))\approx f(x)-\eta f^{{}^{\prime }}(x{)}^2$。

如果导数 $f^{\prime }(x)\ne 0$，那么 $\eta f^{{}^{\prime }}(x{)}^2>0$，即:

                   $f(x-\eta f^{{}^{\prime }}(x))\lesssim f(x)$。

这意味着，如果通过$x=x-\eta f^{{}^{\prime }}(x)$ 来迭代 $x$，函数 $f(x)$ 的值可能会降低。因此，在梯度下降中，选取一个初始值 $x$ 和 $\eta >0$，然后通过上式不断迭代 $x$, 直到达到停止条件。

梯度下降算法中的 $\eta$ 通常叫做学习率，是一个人工设定的超参。学习率过大或过小都会有问题，一个合适的学习率通常是需要多次实验才能找到。如果使用过小的学习率，会导致 $x$ 更新缓慢从而需要更多的迭代才能得到较好的解。如果使用过大的学习率，$|\eta f^{{}^{\prime }}(x)|$ 可能会过大从而使泰勒展开式不再成立，这时我们无法保证迭代 $x$ 会降低 $f(x)$ 的值。

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

                   $f(x)=\frac{1}{n}{\sum }_{i=1}^n{f}_i(x)$，

记$\nabla f(x)$为目标函数 $f(x)$ 有关 $x$ 的梯度偏导数，则目标函数 $f(x)$ 在 $x$ 处的梯度计算为:

                   $\nabla f(x)=\frac{1}{n}{\sum }_{i=1}^n\nabla f_i(x)$。

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

随机梯度下降减少了每次迭代的计算开销。在随机梯度下降的每次迭代中，我们随机采样一个样本 $i\in (1,...,n)$，并计算梯度$\nabla f_i(x)$来迭代$x$: $x=x-\eta \nabla f_i(x)$。这样每次迭代的计算开销就从$O(n)$降到了常数$O(1)$。

小批量随机梯度下降

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

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

                   $g_t=\nabla f_{B_t}(x_{t-1})=\frac{1}{|B|}{\sum }_{i\in B_t}\nabla f_i(x_{t-1})$

来计算时间步 $t$ 的小批量 $B_t$ 上目标函数位于$x_{t-1}$处的梯度$g_t$。这里的 $|B|$ 代表批量大小，即小批量中样本的个数，也是一个超参数。给定学习率${\eta }_t$，小批量随机梯度下降对自变量的迭代为:

                   $x_t=x_{t-1}-{\eta }_t{g}_t$。

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

对于批量梯度下降算法来说，迭代过程中使用的是目标函数的真实梯度，可理解为每次迭代，梯度下降的方向都会使得目标函数更接近最优解。然而对于小批量随机梯度下降算法，每次迭代随机选取$|B|$个样本的梯度平均值，并不代表真实的梯度下降方向，而是在保证整体下降的趋势情况下略微波动。如果学习率大小不变的话，最终会导致目标函数在接近最优解的附近来回波动(蓝色曲线所示)。如果让学习率随着损失函数接近最优解或者迭代次数增加而慢慢减小时，就能使目标函数在更接近最优解的范围内波动。(绿色曲线所示)

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

def sgd(params, lr=0.05):
    for p in params:
        p.data -= lr * p.grad.data

动量法

目标函数有关自变量的梯度代表了目标函数在自变量当前位置下降最快的方向。因此梯度下降也叫做最陡下降。在每次迭代中，梯度下降根据自变量当前位置，沿着当前位置的梯度更新自变量。

然而，当输入为多维向量，且在同一位置，目标函数不同维度梯度值相差较大时，对于给定学习率，梯度下降迭代自变量时，会使自变量在某一维度上移动幅度过大(过小)。

为了解决上述问题，引入了动量法，借助指数加权平均使得自变量更新方向更加一致，从而降低函数发散的可能。设时间步$t$的自变量为$x_t$, 学习率为${\eta }_t$。在时间步$0$，动量法创建速度变量$v_0$，并将其元素初始化成$0$。在时间步 $t>0$，动量法对每次迭代的步骤做如下修改:

                   $v_t=\gamma v_{t-1}+{\eta }_t{g}_t$，
                   $x_t=x_{t-1}-v_t$

其中，动量超参数 $\gamma$ 满足 $0\le \gamma <1$。当$\gamma =0$ 时，动量法等价于小批量随机梯度下降。

为了从数学上理解动量法，我们先解释一下指数加权平均: 给定超参数$0\le \gamma <1$，当前时间步$t$ 的变量$y_t$是上一时间步$t-1$的变量$y_{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 )\ast \gamma x_{t-1}+{\gamma }^2{y}_{t-2}$
                     …

令$n=\frac{1}{(1-\gamma )}$，那么$(1-\frac{1}{n}{)}^n={\gamma }^{\frac{1}{(1-\gamma )}}$。因为：

                   $\underset{n\to \infty }{\lim }(1-\frac{1}{n}{)}^n=\frac{1}{e}$

所以当$\gamma \to 1$时，$\frac{1}{(1-\gamma )}=\frac{1}{e}$。如果把$\frac{1}{e}$当作一个比较小的数，我们可以在近似中忽略所有含有$\frac{1}{(1-\gamma )}$和比$\frac{1}{(1-\gamma )}$更高阶的系数的项。例如，当 $\gamma =0.95$ 时:

                   $y_t\approx 0.05{\sum }_{i=0}^{19}0.95^i{x}_{t-i}$。

因此，在实际中，我们通常将$y_t$看作是对最近$\frac{1}{(1-\gamma )}$个时间步的$x_t$值的加权平均。例如，当$\gamma =0.95$时，$y_t$可以看作对最近20个时间步的$x_t$值的加权平均。

对动量法的速度变量做如下变形：

                   $v_t=\gamma v_{t-1}+(1-\gamma )(\frac{{\eta }_t{g}_t}{1-\gamma })$

由指数加权平均的形式可得，速度变量$v_t$实际上对序列{$\frac{{\eta }_{t-i}{g}_{t-i}}{1-\gamma }:i=0,1,...,\frac{1}{1-\gamma }-1$}做了指数加权平均。换句话说，相比小批量随机梯度下降，动量法在每个时间步的自变量更新量近似于将最近$\frac{1}{1-\gamma }$个时间步的普通更新量做了指数加权移动平均后再除以$1-\gamma$。所以，在动量法中，自变量在各个方向上的移动幅度不仅取决于当前梯度，还取决于过去的各个梯度在各个方向上是否一致。

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

此处实际实现为加权和

AdaGrad算法

动量法依赖指数加权平均使得自变量的更新方向更加一致，从而降低发散的可能。AdaGrad算法则根据自变量在每个维度的梯度值大小来调整各个维度上的学习率，从而避免统一的学习率难以适应所有维度的问题。

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

                   $s_t=s_{t-1}+g_t⨀g_t$

其中$⨀$是按照元素相乘。接着，将目标函数中每个元素的学习率通过按元素运算重新调整一下:
                   $x_t=x_{t-1}-\frac{\eta }{\sqrt{s_t+ϵ}}⨀g_t$

其中 $\eta$ 是学习率，$ϵ$ 是为了维持数值稳定性而添加的常数，如$10^{-6}$。这里开方、除法和乘法的运算都是按元素运算的。

需要注意的是，如果目标函数有关自变量中某个元素的偏导数一直都较大，那么该元素的学习率将会下降较快；反之，如果目标函数有关自变量中某个元素的偏导数一直都很小，那么该元素的学习率将下降较慢。然而，由于$s_t$一直在累加按元素平方的梯度，自变量中每个元素的学习率在迭代过程中一直降低。所以，当学习率在迭代早期降得较快且当前解依然不佳时，在迭代后期由于学习率过小，可能很难找到一个有用的解。

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算法存在的后期学习率过低问题，RMSProp算法对其进行修改。不同于AdaGrad算法里状态变量$s_t$是截至时间步$t$所有小批量随机梯度$g_t$按元素平方和，RMSProp算法将这些梯度按元素平方做指数加权移动平均。具体来说，给定超参数$0\le \gamma <1$，RMSProp算法在时间步$t>0$计算:

                   $s_t=\gamma s_{t-1}+(1-\gamma )g_t⨀g_t$

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

                   $x_t=x_{t-1}-\frac{\eta }{\sqrt{s_t+ϵ}}⨀g_t$

因为RMSProp算法的状态变量$s_t$是对平方项$g_t⨀g_t$的指数加权平均，所以可以看作是最近$\frac{1}{(1-\gamma )}$个时间步的小批量随机梯度平方项的加权平均。在同样的学习率下，RMSProp算法可以更快逼近最优解，并且自变量每个元素的学习率在迭代过程中就不再一直降低。

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)

AdaDelta算法

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

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

                   $s_t=\rho s_{t-1}+(1-\rho )g_t⨀g_t$

与RMSProp算法不同的是，AdaDelta还维护了一个额外的状态变量$\Delta x_t$，其元素同样在时间步0时被初始化为0。我们使用$\Delta x_{t-1}$ 来计算自变量的变化量:

                   $g_t^{{}^{\prime }}=\sqrt{\frac{\Delta x_{t-1}+ϵ}{s_t+ϵ}}⨀g_t$

接着更新自变量:

                   $x_t=x_{t-1}-g_t^{{}^{\prime }}$

最后，使用$\Delta x_t$来记录自变量变化量$g_t^{{}^{\prime }}$按元素平方的指数加权移动平均:

                   $\Delta x_t=\rho \Delta x_{t-1}+(1-\rho )g_t^{{}^{\prime }}⨀g_t^{{}^{\prime }}$

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算法的结合，使用了动量变量$v_t$和RMSProp算法中小批量随机梯度按元素平方的指数加权移动平均变量$s_t$，并在时间步0中将它们中每个元素初始化为0。给定超参数$0\le {\beta }_1<1$(算法作者建议设置为0.9)，时间步$t$的动量变量$v_t$即小批量随机梯度$g_t$的指数加权移动平均:

                   $v_t={\beta }_1{v}_{t-1}+(1-{\beta }_1)g_t$

和RMSProp一样，给定超参数$0\le {\beta }_2<1$(算法作者建议设为0.999)，将小批量随机梯度按元素平方后的项$g_t⨀g_t$做指数加权移动平均得到$s_t$:

                   $s_t={\beta }_2{s}_{t-1}+(1-{\beta }_2)g_t⨀g_t$

由于我们将$v_0$和$s_0$中的元素都初始化为0，在时间步$t$我们得到$v_t=(1-{\beta }_1){\sum }_{i=1}^t{\beta }_1^{t-i}{g}_i$。将过去各个时间步小批量随机梯度的权值相加，得到$(1-{\beta }_1){\sum }_{i=1}^t{\beta }_1^{t-1}=1-{\beta }_1^t$。需要注意的是当$t$较小时，过去各个时间步小批量随机梯度权值之和会比较小。例如，当${\beta }_1=0.9$时，$v1=0.1g_1$。为了消除这样的影响，对于任意时间步$t$，我们可以将$v_t$再除以$1-{\beta }_1^t$，从而使过去各时间步小批量随机梯度权值和为1。这也叫作偏差修正。在Adam算法中，我们对变量$v_t$和$s_t$均做偏差修正:

                   ${\hat{v}}_t=\frac{v_t}{1-{\beta }_1^t}$ ，

                   ${\hat{s}}_t=\frac{s_t}{1-{\beta }_2^t}$

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

                   $g_t^{{}^{\prime }}=\frac{\eta {\hat{v}}_t}{\sqrt{{\hat{s}}_t+ϵ}}$

和AdaGrad算法、RMSProp算法以及AdaDelta算法一样，目标函数自变量中每个元素都分别拥有自己的学习率。最后使用$g_t^{{}^{\prime }}$迭代自变量:

                   $x_t=x_{t-1}-g_t^{{}^{\prime }}$

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

Pytorch实现

Pytorch中定义了以上所有优化函数的优化器实例，只需设置相关超参即可。