DL_Notebook 8.2_Momentum/AdaGrad

深度学习最全优化方法总结比较（SGD，Adagrad，Adadelta，Adam，Adamax，Nadam）
Notebook 8.1中MBGD，依然存在问题：

1. learning rate的选择：
2. 容易收敛到局部最优

so各种优化方法出现。

ill-conditioned problem

矩阵的条件数（condition number）度量了 问题的解对于输入数据中的扰动以及解决问题过程的舍入误差roundoff errors的敏感度。
条件数会告诉你解决一个线性系统会在多大程度上放大数据中的噪声，Think of it as a measure of amplification, a gain.
如果cond = 1e8, 那么当解决一个线性系统时，这个解会把噪声放大约1e8倍。e.g.: 数据是双精度的，噪声是最低有效位那种量级（too小小小小），那条件数影响不大； however，来自现实世界的数据，误差可能在1%左右，当噪声被放大1e8倍，噪声将主导最后的解。

Condition number of Hessian Matrix：
$$con{d}_H=\frac{{\lambda }_{max}}{{\lambda }_{min}}$$

最大学习率

maximum learning rate

• 对于$f(x)$，根据凸优化结论，学习率$\eta <\frac{1}{L}$，取=时，有最快收敛速度，其中$L=ma{x}_x{\nabla }^{2}f(x)$

  比如对一个凸优化表面，大概路线（快+收敛）：
  

• 如果要收敛，需要满足 ：$\eta <\frac{2}{L}$，大于$\frac{2}{L}$会发散

• $\frac{1}{L}<\eta <\frac{2}{L}$时，大概路线如下（不那么快+收敛）：
  

ill-conditioned例子，看看会对优化带来什么影响：
$$\begin{gathered} f(\mathbf{x})=0.1x_1^2+2x_2^2 \\ con{d}_H=\frac{4}{0.2}=20\to \text{ill-conditioned} \end{gathered}$$

%matplotlib inline
import sys
sys.path.append("/home/kesci/input") 
import d2lzh1981 as d2l
import torch

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 * 0.2 * x1, x2 - eta * 4 * x2, 0, 0)

d2l.show_trace_2d(f_2d, d2l.train_2d(gd_2d))

代码中，$\eta =0.4$，同时满足在$x_1,x_2$两个方向上的收敛条件。

• $x_2$方向：$f$对$x_2$求二阶导得$L=4$；$0.25=\frac{1}{L}<\eta <\frac{2}{L}=0.5$。所以可以看到，梯度下降路线在$x_2$方向上，在最优解的两侧来回走Z字路线（不那么快+收敛）。
• $x_1$方向：$f$对$x_1$求二阶导得$L=0.2$；$\eta <\frac{1}{L}=5<\frac{2}{L}=10$，学习率0.4相对于5、10取得比较小，所以优化速度较慢。
• 如果想要快速下降，根据$x_1$的计算结果，$\eta$应该取得接近5？但显然为了满足$x_2$方向的收敛条件0.5，不可以取这么大，快速和收敛难以两全，凸显出ill-conditioned的矛盾。

下面将学习率调得稍大一点，$\eta =0.6$，此时自变量在竖直方向不断越过最优解并逐渐发散。

Supp: Preconditioning

在二阶优化中，我们使用Hessian matrix的逆矩阵(或者pseudo inverse)来左乘梯度向量$i.e.{\Delta }_x=H^{-1}\mathbf{g}$，这样的做法称为precondition
相当于将$H$映射为一个单位矩阵，拥有分布均匀的Spectrum（${\lambda }_{max}$与${\lambda }_{min}$非常接近），也即我们去优化的等价标函数的Hessian matrix为良好的单位矩阵（identity matrix）。
e.g.： 原本的优化问题是扁平的优化表面，preconditioning之后变成了接近圆的优化表面，这样在两个方向（例子中的$x_1,x_2$）都可以选择较大的学习率，收敛速度会更快，不会被某一个方向限制。

Solution to ill-condition

• Preconditioning gradient vector: 被应用于 Adam, RMSProp, AdaGrad, Adelta, KFC, Natural gradient 以及其他二阶优化算法，包括牛顿法。
• Averaging history gradient: like momentum, which allows larger learning rates to accelerate convergence; applied in Adam, RMSProp, SGD momentum.

先讲第二种解决办法，momentum，再讲preconditioning。

Momentum

模拟物理里动量的概念，累计之前的动量来代替真正的梯度。更新公式：
Version 1：
$$\begin{array}{rl}{\mathbf{m}}_t& \leftarrow \beta {\mathbf{m}}_{t-1}+{\eta }_t{\mathbf{g}}_t,\\ {\mathbf{x}}_t& \leftarrow {\mathbf{x}}_{t-1}-{\mathbf{m}}_t,\end{array}$$

Another version （2）：
$$\begin{array}{rl}{\mathbf{m}}_t& \leftarrow \beta {\mathbf{m}}_{t-1}+(1-\beta ){\mathbf{g}}_t,\\ {\mathbf{x}}_t& \leftarrow {\mathbf{x}}_{t-1}-{\alpha }_t{\mathbf{m}}_t,\\ {\alpha }_t& =\frac{{\eta }_t}{1-\beta }\end{array}$$
其中，动量超参数$\beta =0$时，momentum算法相当于小批量随机梯度下降（mini-batch gradient descent）。

对version1，在实践中，通常$\beta =0.5,0.9,0.99$ ，分别对应最大 2 倍、10 倍、100 倍的步长（for more details，搜索Exponential Moving Average（EMA）指数加权移动平均）
和自适应学习率一样，$\beta$ 也可以使用某种策略在训练时进行自适应调整；一般初始值是一个较小的值，随后会慢慢变大。

momentum轨迹，例子f(x)未变：

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

eta, beta = 0.4, 0.5
d2l.show_trace_2d(f_2d, d2l.train_2d(momentum_2d))

momentum（橙） v.s SGD（红）

可以看到使用较小的学习率$\eta =0.4，\beta =0.5$时，动量法在竖直方向上的移动更加平滑，且在水平方向上更快逼近最优解。
$\eta =0.6$时，也不像之前GD那样会发散：

SGD 的更新步长只是梯度乘以学习率；而momentum，步长还取决于历史梯度序列的大小和排列；当许多连续的梯度指向相同的方向时，步长会被不断增大。

在pytorch中，torch.optim.SGD可以实现momentum，指出momentum超参数即可（即前面的$\beta$）

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

AdaGrad

由于ill-conditioned problem，需要选择足够小的学习率使得自变量在梯度值较大的维度上不发散，但这样会导致自变量在梯度值较小的维度上迭代过慢（SGD）。
ill-condition解决办法之一 ，momentum依赖EMA使得自变量的更新方向更加一致，从而降低发散的可能。
AdaGrad算法，属于另一类解决办法preconditiong，它根据自变量在每个维度的梯度值的大小来调整各个维度上的学习率，从而避免统一的学习率难以适应所有维度的问题 。

Algorithm

${\mathbf{s}}_t$：square gradient ，截至时间步$t$所有小批量随机梯度${\mathbf{g}}_t$按元素平方和。$t=0$时，${\mathbf{s}}_t$中每个元素初始化为0。
$$\begin{array}{rl}{\mathbf{s}}_t& \leftarrow {\mathbf{s}}_{t-1}+{\mathbf{g}}_t\odot {\mathbf{g}}_t\\ {\mathbf{x}}_t& \leftarrow {\mathbf{x}}_{t-1}-\frac{\eta }{\sqrt{{\mathbf{s}}_t+ϵ}}\odot {\mathbf{g}}_t\end{array}$$
其中，$\eta$人工设置；常数$ϵ$保证分母非0，维持数值稳定性，如$10^{-6}$。
这里开方、除法和乘法的运算都是按元素运算的。这些按元素运算使得目标函数自变量中每个元素都分别拥有自己的学习率。（$\frac{\eta }{\sqrt{{\mathbf{s}}_t+ϵ}}$看作是学习率）

Feature

• 自变量中每个元素的学习率在迭代过程中一直在降低（或不变）。
• 约束项regularizer$-\frac{1}{\sqrt{{\mathbf{s}}_t+ϵ}}$的存在，使得：如果目标函数有关自变量中某个元素的偏导数一直都较大，那么该元素的学习率将下降较快；反之，如果目标函数有关自变量中某个元素的偏导数一直都较小，该元素的学习率将下降较慢。
• 适合处理稀疏梯度

缺点：

• 由公式可以看出，仍依赖于人工设置一个全局学习率 $\eta$，如果过大，会使regularizer过于敏感，对梯度的调节太大
• 学习率是单调递减的，训练中后期，分母越来越大，可能导致训练困难，甚至提前结束

例子，还是那个f(x)：

%matplotlib inline
import math
import torch
import sys
sys.path.append("/home/kesci/input") 
import d2lzh1981 as d2l

def adagrad_2d(x1, x2, s1, s2):
    g1, g2, eps = 0.2 * x1, 4 * x2, 1e-6  # 前两项为自变量梯度
    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
d2l.show_trace_2d(f_2d, d2l.train_2d(adagrad_2d))

epoch 20, x1 -2.382563, x2 -0.158591

$\eta =0.4$，轨迹平滑，学习率递减，迭代后期自变量移动幅度较小：

$\eta =2$，相对于momentum和SGD来说是一个非常大的学习率，可以看到以很快的速度到达最优点，且比较平衡：

Implement

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_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)

训练：

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

或者使用pytorch提供的AdaGrad算法：

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