深度学习中的优化器比较

一. 几个数学概念

1) 梯度（一阶导数）

考虑一座在 (x1, x2) 点高度是 f(x1, x2) 的山。那么，某一点的梯度方向是在该点坡度最陡的方向，而梯度的大小告诉我们坡度到底有多陡。

2) Hesse 矩阵（二阶导数）

Hesse 矩阵常被应用于牛顿法解决的大规模优化问题(后面会介绍)，主要形式如下：

当 f(x) 为二次函数时，梯度以及 Hesse 矩阵很容易求得。二次函数可以写成下列形式：

其中 x为列向量，A 是 n 阶对称矩阵，b 是 n 维列向量， c 是常数。f(x) 梯度是 Ax+b, Hesse 矩阵等于 A

3) Jacobi 矩阵

Jacobi 矩阵实际上是向量值函数的梯度矩阵，假设F:Rn→Rm 是一个从n维欧氏空间转换到m维欧氏空间的函数。这个函数由m个实函数组成:。这些函数的偏导数(如果存在)可以组成一个m行n列的矩阵(m by n)，这就是所谓的雅可比矩阵。

总结一下

a) 如果 f(x) 是一个标量函数，那么雅克比矩阵是一个向量，等于 f(x) 的梯度， Hesse 矩阵是一个二维矩阵。如果 f(x) 是一个向量值函数，那么Jacobi 矩阵是一个二维矩阵，Hesse 矩阵是一个三维矩阵。

b) 梯度是 Jacobian 矩阵的特例，梯度的 jacobian 矩阵就是 Hesse 矩阵（一阶偏导与二阶偏导的关系）。

二.优化算法介绍

补充在前：实际上在我使用LSTM为流量基线建模时候，发现有效的激活函数是elu、relu、linear、prelu、leaky_relu、softplus，对应的梯度算法是adam、mom、rmsprop、sgd，效果最好的组合是：prelu+rmsprop。

1.批量梯度下降（Batch gradient descent，BGD）

θ=θ−η⋅∇θJ(θ) 

每迭代一步，都要用到训练集的所有数据，每次计算出来的梯度求平均 

η代表学习率LR

2.随机梯度下降（Stochastic Gradient Descent，SGD）

θ=θ−η⋅∇θJ(θ;x(i);y(i)) 

通过每个样本来迭代更新一次，以损失很小的一部分精确度和增加一定数量的迭代次数为代价，换取了总体的优化效率的提升。增加的迭代次数远远小于样本的数量。

缺点：

• 对于参数比较敏感，需要注意参数的初始化 
• 容易陷入局部极小值 
• 当数据较多时，训练时间长 
• 每迭代一步，都要用到训练集所有的数据。

import tensorflow as tf
tf.compat.v1.train.GradientDescentOptimizer()

3. 小批量梯度下降（Mini Batch Gradient Descent，MBGD）

θ=θ−η⋅∇θJ(θ;x(i:i+n);y(i:i+n)) 

为了避免SGD和标准梯度下降中存在的问题，对每个批次中的n个训练样本，这种方法只执行一次更新。【每次更新全部梯度的平均值】

在深度学习中为了解决学习率的灵活性，防止迭代到一定次数时，因学习率过大使得函数值在某个值的附近来回波动而不能继续先最优解的方向迭代，可以对学习率进行指数退化。

import tensorflow as tf
'''
tf.train.exponential_decay
#该函数返回以下结果
decayed_learning_rate = learning_rate *
         decay_rate ^ (global_step / decay_steps)
##例： 以0.96为基数，每100000 步进行一次学习率的衰退
'''
global_step = tf.Variable(0, trainable=False)
starter_learning_rate = 0.1
learning_rate = tf.train.exponential_decay(starter_learning_rate, global_step,
                                           100000, 0.96, staircase=True)
# Passing global_step to minimize() will increment it at each step.
learning_step = (
    tf.train.GradientDescentOptimizer(learning_rate)
    .minimize(...my loss..., global_step=global_step)
)

4.近端梯度下降法 (proximal gradient descent)

近端梯度下降方法，这个是tf取的名字，其实在文章中，作者简称该梯度下降算法为folos，之所以叫接近的梯度下降算法，是因为作者在原梯度下降算法的基础上加入了修正，进而解决了一些不可微的问题，

近端梯度下降法是众多梯度下降 (gradient descent) 方法中的一种，其英文名称为proximal gradident descent，其中，术语中的proximal一词比较耐人寻味，将proximal翻译成“近端”主要想表达"（物理上的）接近"。与经典的梯度下降法和随机梯度下降法相比，近端梯度下降法的适用范围相对狭窄。对于凸优化问题，当其目标函数存在不可微部分（例如目标函数中有  -范数或迹范数）时，近端梯度下降法才会派上用场。

特定的凸优化问题

一般而言，近端梯度下降法常用于解决以下这类优化问题：

假设目标函数  是由  和  叠加而成，其中，限定  是可微的凸函数、  是不可微 (或局部不可微) 的凸函数。

以线性回归问题为例

给定  ，若线性回归表达式为  ，其中，变量  表示系数向量。优化的目标函数可以写成如下形式：

其中，  表示  -范数；  表示  -范数；正则项中的  -范数是为了保证变量 是稀疏的。在这个目标函数中，正则项中的  -范数便是一个不可微的凸函数。由此便引申出采用近端梯度下降法对这类优化问题进行求解。

近端算子 (proximal operator)

一般来说，对于任意向量  ，如果不可微函数为，则近端算子 (proximal operator) 为

其中，  表示关于变量  和函数  的近端算子； 表示关于变量 软阈值函数（下面将着重介绍），选用符号  是因为软阈值soft-thresholding的英文首字母为s。

这条公式想表达的意思是：对于任意给定的  ，我们希望找到使得  最小化的解  ，其中，  是一个新增参数，表示近端梯度下降的步长 (step size)。实际上，近端算子只和不可微的凸函数  有关 (如公式(2)中的)。

当不可微的凸函数的形式为  时，则对应的软阈值函数为

如公式(2)，当  （  -范数作为正则项，其正则项参数为  ）时，软阈值函数则为，其定义为

公式(4)给出的软阈值函数恰好与公式(2)中的目标函数相对应。至于为什么软阈值函数的形式是这样，其推导则需自己百度查找。

近端梯度下降 (proximal gradient descent)

对于优化问题  ，变量  的迭代递推公式为

其中，变量上标的  表示当前迭代次数。

迭代递推公式证明过程（涉及知识：泰勒展开）:

注意：由于公式第三行中的  和第四行中的  均与决策变量  无关，因此公式第三行等于公式第四行。

线性回归问题的迭代过程

回到公式(1)，我们定义了函数  的表达式为

由于  ，根据近端梯度下降法的迭代递推公式可以推导出变量的近端梯度更新公式为

这条更新公式通常也被称为迭代软阈值算法 (iterative soft-thresholding algorithm, ISTA)。

众所周知，如果将这里的  换成  ，则整个线性回归的优化问题就会变成一个最小二乘问题。获取变量  的最优解也会变得非常简单，只需要令  即可，此时的最优解为

其中，  表示正则项参数。

import tensorflow as tf
tf.compat.v1.train.ProximalGradientDescentOptimizer()

当然，这里还有一个很值得思考的问题：

• 对于上述线性回归问题，为何选用近端梯度下降法，而非次梯度 (subgradient descent) 下降法呢？

实际上，次梯度下降法也能求解这种线性回归问题。一般而言，对于损失函数

我们可以得到正则项的梯度为  然而不幸的是，当  时，正则项不可微。为了消除这个影响，我们可以采用次梯度  ，当再遇到  时，次梯度就会等于  .

相应地，损失函数的次梯度为

与梯度下降法类似，我们可以用次梯度下降法使得损失函数最小化。但仔细看一下损失函数次梯度的表达式  ，我们会发现：次梯度法可能不会像我们所期许的那样，得到真正的稀疏解 ，次梯度法只会使解的部分数值 (即  的部分元素) 尽可能小 (尽可能接近  )。因此，只有选用近端梯度下降这类方法，我们才能获得稀疏解，并同时提升迭代过程的收敛性。

5.指数加权平均的概念

指数加权平均(exponentially weighted averges)，也叫指数加权移动平均，是一种常用的序列数据处理方式。

它的计算公式如下：

其中，

• θ_t：为第 t 天的实际观察值，
• V_t: 是要代替 θ_t 的估计值，也就是第 t 天的指数加权平均值，
• β： 为 V_{t-1} 的权重，是可调节的超参。( 0 < β < 1 )

直接看上面的数据图会发现噪音很多，

这时，我们可以用 指数加权平均 来提取这组数据的趋势。

 

可以看出，红色的数据比蓝色的原数据更加平滑，少了很多噪音，并且刻画了原数据的趋势。

指数加权平均，作为原数据的估计值，不仅可以 1. 抚平短期波动，起到了平滑的作用，2. 还能够将长线趋势或周期趋势显现出来。

所以应用比较广泛，在处理统计数据时，在股价等时间序列数据中，CTR 预估中，美团外卖的收入监控报警系统中的 hot-winter 异常点平滑，深度学习的优化算法中都有应用。

import matplotlib.pyplot as plt
import numpy as np
# from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits import mplot3d
np.random.seed(28)
#生成每个月的平均气温数据
average_temperature = [13,10,13,18,23,28,33,36,33,28,]
temperature1 = np.random.normal(loc=13,scale=3,size=(31,))
temperature2 = np.random.normal(loc=10,scale=3,size=(28,))
temperature3 = np.random.normal(loc=13,scale=3,size=(31,))
temperature4 = np.random.normal(loc=18,scale=3,size=(30,))
temperature5 = np.random.normal(loc=23,scale=3,size=(31,))
temperature6 = np.random.normal(loc=28,scale=3,size=(30,))
temperature7 = np.random.normal(loc=33,scale=3,size=(31,))
temperature8 = np.random.normal(loc=36,scale=3,size=(31,))
temperature9 = np.random.normal(loc=33,scale=3,size=(30,))
temperature10 = np.random.normal(loc=28,scale=3,size=(31,))
temperature11 = np.random.normal(loc=23,scale=3,size=(30,))
temperature12 = np.random.normal(loc=18,scale=3,size=(31,))
temperature_data = np.hstack((temperature1,temperature2,temperature3,temperature4,temperature5,temperature6,temperature7,
                  temperature8,temperature9,temperature10,temperature11,temperature12))
days = np.arange(1,len(temperature_data)+1)

#求指数加权平均数
def exp_mov_avg(data,beta,amend=False):
    '''
    :param data: 一维数据
    :param beta: 加权系数
    :param amend: 是否修正
    :return:
    '''
    v=[]
    v_t = 0
    if amend:
        for i,e in enumerate(data):
            print(i)
            v_t = beta * v_t + (1 - beta) * e
            #修正数值
            amend_v = v_t/(1-beta**(i+1))
            v.append(amend_v)
    else:
        for e in data:
            v_t = beta*v_t+(1-beta)*e
            v.append(v_t)
    return v

mov_avg = exp_mov_avg(temperature_data,0.5)
mov_avg2 = exp_mov_avg(temperature_data,0.9)
mov_avg3 = exp_mov_avg(temperature_data,0.98)
mov_avg4 = exp_mov_avg(temperature_data,0.9,True)
plt.scatter(days,temperature_data)
plt.plot(days,mov_avg,'y-',label='beta=0.5')
plt.plot(days,mov_avg2,'r-',label='beta=0.9')
plt.plot(days,mov_avg3,'g-',label='beta=0.98')
plt.plot(days,mov_avg4,color='purple',ls='-',label='beta=0.9 and amend')
plt.legend()
plt.show()

'''使用tensorflow的api'''
# import tensorflow as tf
# tf.compat.v1.enable_eager_execution()
# var = tf.Variable(0,dtype=tf.float32)
# ema=tf.compat.v1.train.ExponentialMovingAverage(0.9)
# mov_avg = []
# for e in temperature_data:
#     new_var = tf.constant(e,dtype=tf.float32)
#     var.assign(new_var)
#     ema.apply([var])
#     mov_avg.append(ema.average(var).numpy())
# 
# plt.scatter(days,temperature_data)
# plt.plot(days,mov_avg,'r-',label='beta=0.9')
# plt.legend()
# plt.show()

β 如何选择

根据前面的计算式子：

将 V_{100} 展开得到：

 
从这里我们就可已看出指数加权平均的名称由来，第100个数据其实是前99个数据加权和，而前面每一个数的权重呈现指数衰减，即越靠前的数据对当前结果的影响较小 

这里可以看出，V_t 是对每天温度的加权平均，之所以称之为指数加权，是因为加权系数是随着时间以指数形式递减的，时间越靠近，权重越大，越靠前，权重越小。

再来看下面三种情况：

当 β = 0.9 时，指数加权平均最后的结果如图红色线所示，代表的是最近 10 天的平均温度值；

当 β = 0.98 时，指结果如图绿色线所示，代表的是最近 50 天的平均温度值；

当 β = 0.5 时，结果如下图黄色线所示，代表的是最近 2 天的平均温度值；

β 越小，噪音越多，虽然能够很快的适应温度的变化，但是更容易出现奇异值。

β 越大，得到的曲线越平坦，因为多平均了几天的温度，这个曲线的波动更小。但有个缺点是，因为只有 0.02 的权重给了当天的值，而之前的数值权重占了 0.98 ，曲线进一步右移，在温度变化时就会适应地更缓慢一些，会出现一定延迟。

通过上面的内容可知，β 也是一个很重要的超参数，不同的值有不同的效果，需要调节来达到最佳效果，一般 0.9 的效果就很好。

缺点：存在开始数据的过低问题，可以通过偏差修正，但是在深度学习的优化算法中一般会忽略这个问题 
 
当t不断增大时，分母逐渐接近1，影响就会逐渐减小了

优点：【相较于滑动窗口平均】 

• 1.占用内存小，每次覆盖即可 
• 2.运算简单

6.Momentum（动量梯度下降法）

上面提到了一些 指数加权平均 的应用，这里我们着重看一下在优化算法中的作用。

动量优化算法在梯度下降法的基础上进行改变，具有加速梯度下降的作用。类别：（1）标准动量优化方法Momentum,(2)NAG动量优化方法（NAG在Tensorflow中与Momentum合并在同一函数tf.train.MomentumOptimizer中，可以通过参数配置启用。）

以 Momentum 梯度下降法为例，

Momentum 梯度下降法，就是计算了梯度的指数加权平均数，并以此来更新权重，它的运行速度几乎总是快于标准的梯度下降算法。

这是为什么呢？

让我们来看一下这个图，

例如这就是我们要优化的成本函数的形状，图中红点就代表我们要达到的最小值的位置，假设我们从左下角这里出发开始用梯度下降法，那么蓝色曲线就是一步一步迭代，一步一步向最小值靠近的轨迹。可以看出这种上下波动，减慢了梯度下降法的速度，而且无法使用更大的学习率，因为如果用较大的学习率，可能会偏离函数的范围。如果有一种方法，可以使得在纵轴上，学习得慢一点，减少这些摆动，但是在横轴上，学习得快一些，快速地从左向右移移向红点最小值，那么训练的速度就可以加快很多。

这个方法就是动量 Momentum 梯度下降法，它在每次计算梯度的迭代中，对 dw 和 db 使用了指数加权平均法的思想。

momentum是模拟物理里动量的概念，积累之前的动量来替代真正的梯度。引入一个累计历史梯度信息动量加速SGD。公式如下： 
 

然而网上更多的是另外一种版本，即去掉（1-β） 

 

相当于上一版本上本次梯度的影响权值*1/(1-β) 

两者效果相当，只不过会影响一些最优学习率的选取 

优点

• 下降初期时，使用上一次参数更新，下降方向一致，乘上较大的μ能够进行很好的加速
• 下降中后期时，在局部最小值来回震荡的时候，gradient→0，β得更新幅度增大，跳出陷阱
• 在梯度改变方向的时候，μ能够减少更新

即在正确梯度方向上加速，并且抑制波动方向张的波动大小，在后期本次计算出来的梯度会很小，以至于无法跳出局部极值，Momentum方法也可以帮助跳出局部极值参数设置 。加上动量项就像从山顶滚下一个球，求往下滚的时候累积了前面的动量(动量不断增加)，因此速度变得越来越快，直到到达终点。同理，在更新模型参数时，对于那些当前的梯度方向与上一次梯度方向相同的参数，那么进行加强，即这些方向上更快了；对于那些当前的梯度方向与上一次梯度方向不同的参数，那么进行削减，即这些方向上减慢了。因此可以获得更快的收敛速度与减少振荡。

β的常用值为0.9，即可以一定意义上理解为平均了前10/9次的梯度。 至于LR学习率的设置，后面所有方法一起总结吧。

from matplotlib import pyplot as plt
import numpy as np

fig = plt.figure()
x = np.arange(-0.8, 1.2, 0.025)
plt.plot(x,-x**3-x**2+2*x**4)
plt.title("y = 2*x^4-x^3-x^2")
def f(x):
    return 2*x**4-x**3-x**2
def h(x):
    return 8*x**3 - 3*x**2 - 2*x
η = 0.05
α = 0.9
v = 0
x = -0.8
iters = 0
X = []
Y = []
while iters<12:
    iters+=1
    X.append(x)
    Y.append(f(x))
    v = α*v - η*h(x)
    x = x + v
    print(iters,x)
plt.plot(X,Y)
plt.show()

import tensorflow as tf
tf.compat.v1.train.MomentumOptimizer(self, learning_rate, momentum,use_locking=False, name="Momentum", use_nesterov=False)

7.Nesterov accelerated gradient (NAG)

由上面Momentum公式，可以画出Momentum第一和第二次参数更新的方向：

Momentum的想法很简单，就是多更新一部分上一次迭代的更新量，来平滑这一次迭代的梯度。从物理的角度上解释，就像是一个小球滚落的时候会受到自身历史动量的影响，所以才叫动量（Momentum）算法。这样做直接的效果就是使得梯度下降的的时候转弯掉头的幅度不那么大了，于是就能够更加平稳、快速地冲向局部最小点。

然后NAG就对Momentum说：“既然我都知道我这一次一定会走αv_t-1的量，那么我何必还用现在这个位置的梯度呢？我直接先走到αv_t-1之后的地方，然后再根据那里的梯度再前进一下，岂不美哉？”所以就有了下面的公式：

跟上面Momentum公式的唯一区别在于，梯度不是根据当前参数位置W_t，而是根据先走了本来计划要走的一步后，达到的参数位置W_t - αv_t-1计算出来的。

 

动量解决SGD的两个问题：（1）SGD引入的噪声（2）Hessian矩阵病态（SGD收敛过程的梯度相比正常来回摆动幅度较大，NAG是Momentum变种，Nesterov accelerated gradient（NAG，涅斯捷罗夫梯度加速）不仅增加了动量项，并且在计算参数的梯度时，在损失函数中减去了动量项，即计算∇θJ(θ−γνt−1)，这种方式预估了下一次参数所在的位置。 NAG的计算在模型参数施加当前速度之后，可以理解为在Momentum 中引入了一个校正因子。在Momentum中，小球会盲目的跟从下坡的梯度，易发生错误。因此，需要提前知道下降的方向，同时，在快到目标点时速度会有所下降，以不至于超出。    可以表示小球下一个大概的位置，从而知道下一个位置的梯度，然后使用当前位置来更新参数。

优点： 

这种基于预测的更新方法，使我们避免过快地前进，并提高了算法地响应能力，NAG 可以使 RNN 在很多任务上有更好的表现。大大改进了 RNN 在一些任务上的表现【为什么对RNN好呢，不懂啊】 

没有对比就没有伤害，NAG方法收敛速度明显加快。波动也小了很多。实际上NAG方法用到了二阶信息，所以才会有这么好的结果。先按照原来的梯度走一步的时候已经求了一次梯度，后面再修正的时候又求了一次梯度，所以是二阶信息。 

为了从另一个角度更加深入地理解这个算法，我们可以对NAG原来的更新公式进行变换，得到这样的等效形式:

NAG的原始形式到等效形式的推导。由

其中，和分别是这一次和上一次的更新方向，表示目标函数在处的梯度，超参数是对上一次更新方向的衰减权重，所以一般是0到1之间，是学习率。

可得

记:

上式代入上上式，就得到了NAG等效形式的第二个式子：

对展开可得:

于是我们可以写出的形式，然后用减去消去后面的无穷多项，就得到了NAG等效形式的第一个式子：

最终我们就得到了NAG的等效形式：

这个NAG的等效形式与Momentum的区别在于，本次更新方向多加了一个，它的直观含义就很明显了：如果这次的梯度比上次的梯度变大了，那么有理由相信它会继续变大下去，那我就把预计要增大的部分提前加进来；如果相比上次变小了，也是类似的情况。这样的解释听起来好像和原本的解释一样玄，但是读者可能已经发现了，这个多加上去的项不就是在近似目标函数的二阶导嘛！所以NAG本质上是多考虑了目标函数的二阶导信息，怪不得可以加速收敛了！其实所谓“往前看”的说法，在牛顿法这样的二阶方法中也是经常提到的，比喻起来是说“往前看”，数学本质上则是利用了目标函数的二阶导信息。

结论：

在原始形式中，Nesterov Accelerated Gradient（NAG）算法相对于Momentum的改进在于，以“向前看”看到的梯度而不是当前位置梯度去更新。经过变换之后的等效形式中，NAG算法相对于Momentum多了一个本次梯度相对上次梯度的变化量，这个变化量本质上是对目标函数二阶导的近似。由于利用了二阶导的信息，NAG算法才会比Momentum具有更快的收敛速度。

参数设置： 
同Momentum，并把use_nesterov设为True

import tensorflow as tf
tf.compat.v1.train.MomentumOptimizer(learning_rate,momentum,use_nesterov=True)

momentum项和nesterov项都是为了使梯度更新更加灵活，对不同情况有针对性。但是，人工设置一些学习率总还是有些生硬，接下来介绍几种自适应学习率的方法（AdaGrad，RMSProp，Adam，AdaDelta等。优点：无需人为调节学习率，可以自动调节。缺点；随着迭代次数增多，学习率越来越小，最终趋于0。）

8.Adagrad

前面的一系列优化算法有一个共同的特点，就是对于每一个参数都用相同的学习率进行更新。但是在实际应用中各个参数的重要性肯定是不一样的，所以我们对于不同的参数要动态的采取不同的学习率，让目标函数更快的收敛。 

adagrad方法是将每一个参数的每一次迭代的梯度取平方累加再开方，用基础学习率除以这个数，来做学习率的动态更新。然后给不同的参数不同的学习率【这样每一个参数的学习率就与他们的梯度有关系了，那么每一个参数的学习率就不一样了！也就是所谓的自适应学习率】 

假设N元函数f(x)，针对一个自变量研究Adagrad梯度下降的迭代过程，

 

可以看出，Adagrad算法中有自适应调整梯度的意味（adaptive gradient），学习率需要除以一个东西，这个东西就是前n次迭代过程中偏导数的平方和再加一个常量最后开根号。

举例：

使用Adagrad算法求y = x²的最小值点

导函数为g(x) = 2x

初始化x⁽⁰⁾ = 4，学习率η=0.25，ε=0.1

from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
x = np.arange(-4, 4, 0.025)
plt.plot(x,x**2)
plt.title("y = x^2")
def f(x):
    return x**2
def h(x):
    return 2*x
η = 0.25
ε = 0.1
x = 4
iters = 0
sum_square_grad = 0
X = []
Y = []
while iters<12:
    iters+=1
    X.append(x)
    Y.append(f(x))
    sum_square_grad += h(x)**2
    x = x - η/np.sqrt(sum_square_grad+ε)*h(x)
    print(iters,x)
plt.plot(X,Y,"ro")
ax = plt.subplot()
for i in range(len(X)):
    ax.text(X[i], (X[i])**2, "({:.3f},{:.3f})".format(X[i], (X[i])**2), color='red')
plt.show()

公式的另一种写法：

 其中，⊙是按元素相乘，r为梯度累积变量，r的初始值为0。ε为全局学习率，需要自己设置。δ为小常数，为了防止分母为零。

参数设置： 

• 只需要设置初始学习率，后面学习率会自我调整，越来越小

优点：

• 前期Gt较小的时候， regularizer较大，能够放大梯度
• 后期Gt较大的时候，regularizer较小，能够约束梯度
• 学习率将随着梯度的倒数增长，也就是说较大梯度具有较小的学习率，而较小的梯度具有较大的学习率，可以解决普通的sgd方法中学习率一直不变的问题
• 适合处理稀疏梯度:相当于为每一维参数设定了不同的学习率：压制常常变化的参数，突出稀缺的更新。能够更有效地利用少量有意义样本。在数据分布稀疏的场景上，能更好利用稀疏梯度的信息，比标准的SGD算法更有效地收敛。

缺点： 

• 主要缺陷来自分母项的对梯度平方不断累积，随之时间步地增加，分母项越来越大，最终导致学习率收缩到太小无法进行有效更新。如果初始梯度很大的话，会导致整个训练过程的学习率一直很小，从而导致学习时间变长。

import tensorflow as tf
tf.compat.v1.train.AdagradOptimizer()
tf.compat.v1.train.ProximalAdagradOptimizer() #用来优化目标函数中存在有不可部分的或是加了L1或L2正则的。
tf.compat.v1.train.AdagradDAOptimizer() 

注意：AdagradDA全称是Adagrad双重平均算法（Adagrad Dual Averaging algorithm），这个算法我实在是看不懂了，光是数学原理就要好多页，这里贴出论文的链接http://www.jmlr.org/papers/volume12/duchi11a/duchi11a.pdf

Adagrad的一大优势时可以避免手动调节学习率，比如设置初始的缺省学习率为0.01，然后就不管它，另其在学习的过程中自己变化。当然它也有缺点，就是它计算时要在分母上计算梯度平方的和，由于所有的梯度的平方必为正数，这样就造成在训练的过程中，分母累积的和会越来越大。这样学习到后来的阶段，网络的更新能力会越来越弱，能学到的更多知识的能力也越来越弱，因为学习率会变得极其小【就会提前停止学习】，为了解决这样的问题又提出了Adadelta算法。

9.RMSprop

由于AdaGrad算法的机制，导致每个元素的学习率在迭代过程中只能降低或者不变，因此很可能出现早期迭代到不好的极值点之后，由于学习率太小而无法冲出这个极值点导致最后收敛到的解不优，为了解决这一问题，RMSProp是基于AdaGrad算法做了一点小修改，其更新公式为：

 

其中，ϵ是防止分母爆0的常数。

这样，就有了一个改进版的AdaGrad。

该方法即Tieleman&Hinton的RMSProp，由于RMSProp和AdaDelta是同年出现的，

Matthew D. Zeiler并不知道这种改进的AdaGrad被祖师爷命名了。

RMSProp利用了二阶信息做了Gradient优化，在BatchNorm之后，对其需求不是很大。

但是没有根本实现自适应的学习率，依然需要线性搜索初始学习率，然后对其逐数量级下降。

另外，RMSProp的学习率数值与MomentumSGD差别甚大，需要重新线性搜索初始值。

注意：建议取值为1。
特点：

• 其实RMSprop依然依赖于全局学习率
• RMSprop算是Adagrad的一种发展，和Adadelta的变体，效果趋于二者之间
• 适合处理非平稳目标（也就是与时间有关的）
• 对于RNN效果很好，因为RMSprop的更新只依赖于上一时刻的更新，所以适合。

from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
x = np.arange(-4, 4, 0.025)
plt.plot(x,x**2)
plt.title("y = x^2")
def f(x):
    return x**2
def h(x):
    return 2*x
g = 1
x = 4
ρ = 0.9
η = 0.01
ε = 10e-10
iters = 0
X = []
Y = []
while iters<12:
    iters+=1
    X.append(x)
    Y.append(f(x))
    g = ρ*g+(1-ρ)*h(x)**2
    x = x - η/np.sqrt(g+ε)*h(x)
    print(iters,x)
plt.plot(X,Y,"ro")
ax = plt.subplot()
for i in range(len(X)):
    ax.text(X[i], (X[i])**2, "({:.3f},{:.3f})".format(X[i], (X[i])**2), color='red')
plt.show()

import tensorflow as tf
tf.compat.v1.train.RMSPropOptimizer(learning_rate)

10.Adadelta

Adagrad会累加之前所有的梯度平方，而Adadelta只累加固定大小的项【其实就是相当于指数滑动平均，只用了前多少步的梯度平方平均值】，并且也不直接存储这些项，仅仅是近似计算对应的平均值【这也就是指数滑动平均的优点】 

整理得更新公式为：

优点： 

• 不用依赖于全局学习率了 
• 训练初中期，加速效果不错，很快 
• 避免参数更新时两边单位不统一的问题 

缺点： 

• 训练后期，反复在局部最小值附近抖动

import tensorflow as tf
tf.compat.v1.train.AdadeltaOptimizer()

11.Adam

Adam = Adaptive + Momentum，顾名思义Adam集成了SGD的一阶动量和RMSProp的二阶动量。 

 

 就像Adadelta和RMSprop一样Adam会存储之前衰减的平方梯度，同时它也会保存之前衰减的梯度。经过一些处理之后再使用类似Adadelta和RMSprop的方式更新参数。

特点：

• 结合了Adagrad善于处理稀疏梯度和RMSprop善于处理非平稳目标的优点
• 对内存需求较小
• 为不同的参数计算不同的自适应学习率
• 也适用于大多非凸优化
• 适用于大数据集和高维空间

from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
x = np.arange(-4, 4, 0.025)
plt.plot(x,x**2)
plt.title("y = x^2")
def f(x):
    return x**2
def h(x):
    return 2*x
x = 4
m = 0
v = 0
β1 = 0.9
β2 = 0.999
η = 0.001
ε = 10e-8
iters = 0
X = []
Y = []
while iters<12:
    iters+=1
    X.append(x)
    Y.append(f(x))
    m = β1*m + (1-β1)*h(x)
    v = β2*v + (1-β2)*h(x)**2
    m_het = m/(1-β1**iters)
    v_het = v/(1-β2**iters)
    x = x - η/np.sqrt(v_het+ε)*m_het
    print(iters,x)
plt.plot(X,Y,"ro")
ax = plt.subplot()
for i in range(len(X)):
    ax.text(X[i], (X[i])**2, "({:.3f},{:.3f})".format(X[i], (X[i])**2), color='red')
plt.show()

import tensorflow as tf
tf.compat.v1.train.AdamOptimizer()

12.Adamax

13.Nadam

tensorflow的优化器及Keras中的默认参数

import tensorflow as tf
tf.compat.v1.train.GradientDescentOptimizer()
tf.compat.v1.train.AdadeltaOptimizer()
tf.compat.v1.train.AdagradOptimizer()
tf.compat.v1.train.AdagradDAOptimizer()
tf.compat.v1.train.MomentumOptimizer()
tf.compat.v1.train.AdamOptimizer()
tf.compat.v1.train.FtrlOptimizer()
tf.compat.v1.train.ProximalGradientDescentOptimizer()
tf.compat.v1.train.ProximalAdagradOptimizer()
tf.compat.v1.train.RMSPropOptimizer()

tf.keras.optimizers.SGD(lr=0.001,momentum=0.9)
tf.keras.optimizers.Adagrad(lr=0.01,epsilon=1e-8)
tf.keras.optimizers.Adadelta(lr=0.01,rho=0.95,epsilon=1e-8)
tf.keras.optimizers.RMSprop(lr=0.001,rho=0.9,epsilon=1e-8)
tf.keras.optimizers.Adam(lr=0.001,beta_1=0.9,beta_2=0.999,epsilon=1e-8)

tensorflow中分布式Tensorflow中同步梯度更新（tf.compat.v1.train.SyncReplicasOptimizer）

在tensorflow的训练中，分布式可以大大的加快模型训练速度，但是分布式怎么分配和参数设定，都和SyncReplicasOptimizer这个函数有很大关系。

首先这个SyncReplicasOptimizer函数是专门出来分布式深度学习中的同步梯度下降的(要用异步梯度的可以直接无视)

def __init__(self,
opt: Any,
replicas_to_aggregate: Any,
total_num_replicas: Any = None,
variable_averages: Any = None,
variables_to_average: Any = None,
use_locking: bool = False,
name: str = "sync_replicas") -> None
Construct a sync_replicas optimizer.
Params:
opt – The actual optimizer that will be used to compute and apply the gradients. Must be one of the Optimizer classes.
replicas_to_aggregate – number of replicas to aggregate for each variable update.
total_num_replicas – Total number of tasks/workers/replicas, could be different from replicas_to_aggregate. If total_num_replicas > replicas_to_aggregate: it is backup_replicas + replicas_to_aggregate. If total_num_replicas replicas_to_aggregate: Replicas compute multiple batches per update to variables.
variable_averages – Optional `ExponentialMovingAverage` object, used to maintain moving averages for the variables passed in `variables_to_average`.
variables_to_average – a list of variables that need to be averaged. Only needed if variable_averages is passed in.
use_locking – If True use locks for update operation.
name – string. Optional name of the returned operation.

 opt：将用于计算和应用梯度的实际优化器。必须是Optimizer类之一。
 variable_averages：可选的`ExponentialMovingAverage`对象，用于保持传入的变量的移动平均值
variables_to_average：需要平均的变量列表。只要如果传入variable_averages则需要。
use_locking：如果True使用锁定进行更新操作。
name：string。返回操作的可选名称。

除开opt之外，其他可以先不指定，我们这里重点看看replicas_to_aggregate和total_num_replicas这两个参数，这是非常重要的。

replicas_to_aggregate：是指我们所需要通过分布式的worker计算完成之后所需要的副本数，这里可以认为就是梯度的个数

total_num_replicas：是指我们指定出入计算样本的数量（几个机器就分配几个batch_size） 

#同步模式计算更新梯度
rep_op = tf.train.SyncReplicasOptimizer(optimizer,
                                        replicas_to_aggregate=len(worker_hosts),
                                        total_num_replicas=len(worker_hosts),
                                        use_locking=True)

我们可以看出，replicas_to_aggregate和 total_num_replicas都是等于集群中worker的数量的，但是源码中对于这两个参数又添加了额外的解释，(replicas_to_aggregate以下简称副本数，total_num_replicas以下简称计算数)

replicas_to_aggregate: number of replicas to aggregate for each variable
update.
total_num_replicas: Total number of tasks/workers/replicas, could be
different from replicas_to_aggregate.
If total_num_replicas > replicas_to_aggregate: it is backup_replicas +
replicas_to_aggregate.
If total_num_replicas < replicas_to_aggregate: Replicas compute
multiple batches per update to variables.

 什么意思呢？意思就是说，副本数可以和计算数可以是不相等的。就好比有三个人交了三张试卷一样，三个人就是计算数，三张试卷就是副本数。

 源码中说，If total_num_replicas > replicas_to_aggregate: it is backup_replicas + replicas_to_aggregate.

这是指如果计算数大于副本数，那么多的那个就是backup worker，即备份工作节点。

比如：四个人交了三张试卷，没交试卷的我们就认为是备份工作节点   

那这又是什么玩意儿呢，图啥呢？

在同步梯度下降中，所有worker都要等大家一起计算完成之后才能进行计算，因此会出现等待机制。就像四个人做试卷，但是每个人的速度不一样，但是要同时交卷，所以我们要等待交卷最慢的那个人算完之后才能算完，这样很浪费算的很快的那个人的能力，所以就有了backup的提出，4个人交卷我们只要三个人的试卷，最慢的那个我们不要了，直接清除掉，开始下一轮的考试~

因此，加入backup worker的方法可以加速分布式同步梯度下法的训练模型，那么backup取几个呢？Revisiting Distributed Synchronous SGD这篇论文中认为是20：1的比例。

还有一种情况，那就是计算数大于副本数。当计算数小于副本数，那么他会多计算其他worker的样本，就好比三个人(给三人一人发一个试卷，分配的是没多余的)，但是要交四份试卷，那就先各算各的，先计算完成的再从已知的三份试卷中拿一个计算，然后再等着这三个人一共上交了四份试卷之后，这样就ok了。有什么用呢？我也不清楚，我们看下一步

当我们说这计算数和副本数的时候总是和机器的物理数量有关系（gpu数量），那么其实这两个都可以不等于机器的物理数量。例如：

gpu数量：3           那么len（worker）就是3

计算数：10           那么这10就指的是一次给分10个batch_size，然后让3个worker去算

副本数：8             那么这是个batch_size算完之后我只要8个梯度，然后平均，参数服务器的什么后续操作

import numpy as np
import tensorflow as tf

F = tf.flags.FLAGS
# Define parameters
FLAGS = tf.flags.FLAGS
tf.flags.DEFINE_float('learning_rate', 0.00003, 'Initial learning rate.')
tf.flags.DEFINE_integer('steps_to_validate', 1000,
                        'Steps to validate and print loss')

# For distributed
tf.flags.DEFINE_string("ps_hosts", "",
                       "Comma-separated list of hostname:port pairs")
tf.flags.DEFINE_string("worker_hosts", "",
                       "Comma-separated list of hostname:port pairs")
tf.flags.DEFINE_string("job_name", "", "One of 'ps', 'worker'")
tf.flags.DEFINE_integer("task_index", 0, "Index of task within the job")
tf.flags.DEFINE_integer("issync", 0, "是否采用分布式的同步模式，1表示同步模式，0表示异步模式")

# Hyperparameters
learning_rate = FLAGS.learning_rate
steps_to_validate = FLAGS.steps_to_validate

def main(_):
    ps_hosts = FLAGS.ps_hosts.split(",")
    worker_hosts = FLAGS.worker_hosts.split(",")
    cluster = tf.train.ClusterSpec({"ps": ps_hosts, "worker": worker_hosts})
    server = tf.train.Server(cluster, job_name=FLAGS.job_name, task_index=FLAGS.task_index)

    issync = FLAGS.issync
    if FLAGS.job_name == "ps":
        server.join()
    elif FLAGS.job_name == "worker":
        with tf.device(tf.train.replica_device_setter(
                worker_device="/job:worker/task:%d" % FLAGS.task_index,
                cluster=cluster)):
            global_step = tf.Variable(0, name='global_step', trainable=False)

            input = tf.placeholder("float")
            label = tf.placeholder("float")

            weight = tf.get_variable("weight", [1], tf.float32, initializer=tf.random_normal_initializer())
            biase = tf.get_variable("biase", [1], tf.float32, initializer=tf.random_normal_initializer())
            pred = tf.matmul(input, weight) + biase

            loss_value = loss(label, pred)
            optimizer = tf.train.GradientDescentOptimizer(learning_rate)

            grads_and_vars = optimizer.compute_gradients(loss_value)
            if issync == 1:
                # 同步模式计算更新梯度
                rep_op = tf.train.SyncReplicasOptimizer(optimizer,
                                                        replicas_to_aggregate=len(
                                                            worker_hosts),
                                                        replica_id=FLAGS.task_index,
                                                        total_num_replicas=len(
                                                            worker_hosts),
                                                        use_locking=True)
                train_op = rep_op.apply_gradients(grads_and_vars,
                                                  global_step=global_step)
                init_token_op = rep_op.get_init_tokens_op()
                chief_queue_runner = rep_op.get_chief_queue_runner()
            else:
                # 异步模式计算更新梯度
                train_op = optimizer.apply_gradients(grads_and_vars,
                                                     global_step=global_step)

            init_op = tf.initialize_all_variables()

            saver = tf.compat.v1.train.Saver()
            tf.summary.scalar('cost', loss_value)
            summary_op = tf.summary.merge_all()

        sv = tf.train.Supervisor(is_chief=(FLAGS.task_index == 0),
                                 logdir="./checkpoint/",
                                 init_op=init_op,
                                 summary_op=None,
                                 saver=saver,
                                 global_step=global_step,
                                 save_model_secs=60)

        with sv.prepare_or_wait_for_session(server.target) as sess:
            # 如果是同步模式
            if FLAGS.task_index == 0 and issync == 1:
                sv.start_queue_runners(sess, [chief_queue_runner])
                sess.run(init_token_op)
            step = 0
            while step < 1000000:
                train_x = np.random.randn(1)
                train_y = 2 * train_x + np.random.randn(1) * 0.33 + 10
                _, loss_v, step = sess.run([train_op, loss_value, global_step],
                                           feed_dict={input: train_x, label: train_y})
                if step % steps_to_validate == 0:
                    w, b = sess.run([weight, biase])
                    print("step: %d, weight: %f, biase: %f, loss: %f" % (step, w, b, loss_v))

        sv.stop()


def loss(label, pred):
    return tf.square(label - pred)


if __name__ == "__main__":
    tf.app.run()

14.牛顿法——二阶优化方法

1）牛顿法（Newton's method）

牛顿法是一种在实数域和复数域上近似求解方程的方法。方法使用函数f (x)的泰勒级数的前面几项来寻找方程f (x) = 0的根。牛顿法最大的特点就在于它的收敛速度很快。

具体步骤：

首先，选择一个接近函数 f (x)零点的 x₀，计算相应的 f (x₀) 和切线斜率f  ' (x₀)（这里f ' 表示函数 f  的导数）。然后我们计算穿过点(x_0,  f  (x₀)) 并且斜率为f '(x₀)的直线和 x 轴的交点的x坐标，也就是求如下方程的解：

我们将新求得的点的 x 坐标命名为x₁，通常x₁会比x₀更接近方程f  (x) = 0的解。因此我们现在可以利用x₁开始下一轮迭代。迭代公式可化简为如下所示：

已经证明，如果f  ' 是连续的，并且待求的零点x是孤立的，那么在零点x周围存在一个区域，只要初始值x₀位于这个邻近区域内，那么牛顿法必定收敛。 并且，如果f  ' (x)不为0, 那么牛顿法将具有平方收敛的性能. 粗略的说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍。下图为一个牛顿法执行过程的例子。

　　由于牛顿法是基于当前位置的切线来确定下一次的位置，所以牛顿法又被很形象地称为是"切线法"。牛顿法的搜索路径（二维情况）如下图所示：

　　牛顿法搜索动态示例图：

 

在最速下降法（梯度下降）中，我们看到，该方法主要利用的是目标函数的局部性质，具有一定的“盲目性”。牛顿法则是利用局部的一阶和二阶偏导信息，推测整个目标函数的形状，进而可以求得出近似函数的全局最小值，然后将当前的最小值设定近似函数的最小值。相比最速下降法，牛顿法带有一定对全局的预测性，收敛性质也更优良。牛顿法的主要推导过程如下：

第一步，利用 Taylor 级数求得原目标函数的二阶近似：

 

第二步，把 x 看做自变量， 所有带有 x^k 的项看做常量，令一阶导数为 0 ，即可求近似函数的最小值：

 

即：

 

第三步，将当前的最小值设定近似函数的最小值（或者乘以步长)。

牛顿法是二阶算法，因为该算法使用了海森矩阵（Hessian matrix）求权重的二阶偏导数。牛顿法的目标就是采用损失函数的二阶偏导数寻找更好的训练方向。

因此牛顿下降法是用二次曲面去拟合当前的局部曲面,现在采用如下表示：f(wi) = fi、ᐁf(wi) = gi 和 Hf(wi) = Hi。在w0点使用泰勒级数展开式二次逼近函数f.

H0为函数f在点w0的海森矩阵，对w求导，得到g,通过将g设定为0，我们就可以找到f(w)的最小值，也就得到了以下方程式。

因此，从参数向量w0开始，牛顿法服从以下方式进行迭代：

向量 Hi-1·gi也就是所说的牛顿下降步（Newton's step）。注意，参数的这些变化将朝着极大值而不是极小值逼近，出现这样的情况是因为海森矩阵非正定。

因此在不能保证矩阵正定的情况下，损失函数并不能保证在每一次迭代中都是减少的。为了防止上述问题，牛顿法的方程式通常可以修改为：

牛顿法主要存在的问题是：

• Hesse 矩阵不可逆时无法计算
• 矩阵的逆计算复杂为 n 的立方，当问题规模比较大时，计算量很大，解决的办法是采用拟牛顿法如 BFGS, L-BFGS, DFP, Broyden’s Algorithm 进行近似。
• 如果初始值离局部极小值太远，Taylor 展开并不能对原函数进行良好的近似

import numpy as np
def Jacobian(x):
    return np.array([x[0], 0.4 * x[1], 1.2 * x[2]])

def Hessian(x):
    return np.array([[1, 0, 0], [0, 0.4, 0], [0, 0, 1.2]])

def Newton(x0):
    i = 0
    iMax = 10
    x = x0
    Delta = 1
    alpha = 1

    while i < iMax and Delta > 10 ** (-5):
        p = np.dot(np.linalg.inv(Hessian(x)), Jacobian(x))
        xOld = x
        x = x + alpha * p
        Delta = sum((x - xOld) ** 2)
        i += 1
    print(x)
x0 = np.array([-2, 2, -2])
Newton(x0)

关于牛顿法和梯度下降法的效率对比：

从本质上去看，牛顿法是二阶收敛，梯度下降是一阶收敛，所以牛顿法就更快。如果更通俗地说的话，比如你想找一条最短的路径走到一个盆地的最底部，梯度下降法每次只从你当前所处位置选一个坡度最大的方向走一步，牛顿法在选择方向时，不仅会考虑坡度是否够大，还会考虑你走了一步之后，坡度是否会变得更大。所以，可以说牛顿法比梯度下降法看得更远一点，能更快地走到最底部。（牛顿法目光更加长远，所以少走弯路；相对而言，梯度下降法只考虑了局部的最优，没有全局思想。）

根据wiki上的解释，从几何上说，牛顿法就是用一个二次曲面去拟合你当前所处位置的局部曲面，而梯度下降法是用一个平面去拟合当前的局部曲面，通常情况下，二次曲面的拟合会比平面更好，所以牛顿法选择的下降路径会更符合真实的最优下降路径。

 

注：红色的牛顿法的迭代路径，绿色的是梯度下降法的迭代路径。

牛顿法的优缺点总结：

• 优点：二阶收敛，收敛速度快（迭代次数少）；
• 缺点：牛顿法是一种迭代算法，每一步都需要求解目标函数的Hessian矩阵的逆矩阵，计算比较复杂。

15.拟牛顿法（Quasi-Newton method）

拟牛顿法是求解非线性优化问题最有效的方法之一，于20世纪50年代由美国Argonne国家实验室的物理学家W.C.Davidon所提出来。Davidon设计的这种算法在当时看来是非线性优化领域最具创造性的发明之一。不久R. Fletcher和M. J. D. Powell证实了这种新的算法远比其他方法快速和可靠，使得非线性优化这门学科在一夜之间突飞猛进。

拟牛顿法的本质思想是改善牛顿法每次需要求解复杂的Hessian矩阵的逆矩阵的缺陷，它使用正定矩阵来近似Hessian矩阵的逆，从而简化了运算的复杂度。拟牛顿法和最速下降法一样只要求每一步迭代时知道目标函数的梯度。通过测量梯度的变化，构造一个目标函数的模型使之足以产生超线性收敛性。这类方法大大优于最速下降法，尤其对于困难的问题。另外，因为拟牛顿法不需要二阶导数的信息，所以有时比牛顿法更为有效。如今，优化软件中包含了大量的拟牛顿算法用来解决无约束，约束，和大规模的优化问题。

具体步骤：

拟牛顿法的基本思想如下。首先构造目标函数在当前迭代x_k的二次模型：

这里B _k是一个对称正定矩阵，于是我们取这个二次模型的最优解作为搜索方向，并且得到新的迭代点：

其中我们要求步长a_k 满足Wolfe条件。这样的迭代与牛顿法类似，区别就在于用近似的Hesse矩阵B_k   代替真实的Hesse矩阵。所以拟牛顿法最关键的地方就是每一步迭代中矩阵B_k的更新。现在假设得到一个新的迭代x_k+1，并得到一个新的二次模型：

我们尽可能地利用上一步的信息来选取B _k。具体地，我们要求  

从而得到

这个公式被称为割线方程。 常用的拟牛顿法有DFP算法和BFGS算法。

因为需要很多的操作求解海森矩阵的值还有计算矩阵的逆，应用牛顿法的计算量是十分巨大的。因此有一种称之为拟牛顿法（quasi-Newton）或

变量矩阵法来解决这样的缺点。这些方法并不是直接计算海森矩阵然后求其矩阵的逆，拟牛顿法是在每次迭代的时候计算一个矩阵，其逼近海森矩阵的逆。

最重要的是，该逼近值只是使用损失函数的一阶偏导来计算。

海森矩阵由损失函数的二阶偏导组成，拟牛顿法背后的思想主要是仅使用损失函数的一阶偏导数，通过另一矩阵G逼近海森矩阵的逆。拟牛顿法的公式可以表示为：

学习速率η可以设定为固定常数，也可以通过单变量函数优化得到。其中矩阵G逼近海森矩阵的逆，且有不同的方式进行逼近。通常，最常用的两种方式是Davidon–Fletcher–Powell formula (DFP) 和 the Broyden–Fletcher–Goldfarb–Shanno formula (BFGS)。

拟牛顿法的训练过程流程图如下所示，从图中我们可以看出来模型是通过第一次计算拟牛顿训练方向而优化参数的，然后再寻找适当的学习速率。
拟牛顿法适用于绝大多数案例中：它比梯度下降和共轭梯度法收敛更快，并且也不需要确切地计算海森矩阵及其逆矩阵。

16.Levenberg–Marquardt Algorithm

Levenberg–Marquardt algorithm 能结合以上两种优化方法的优点，并对两者的不足做出改进。与 line search 的方法不同，LMA 属于一种“信赖域法”(trust region)，牛顿法实际上也可以看做一种信赖域法，即利用局部信息对函数进行建模近似，求取局部最小值。所谓的信赖域法，就是从初始点开始，先假设一个可以信赖的最大位移 s（牛顿法里面 s 为无穷大），然后在以当前点为中心，以 s 为半径的区域内，通过寻找目标函数的一个近似函数（二次的）的最优点，来求解得到真正的位移。在得到了位移之后，再计算目标函数值，如果其使目标函数值的下降满足了一定条件，那么就说明这个位移是可靠的，则继续按此规则迭代计算下去；如果其不能使目标函数值的下降满足一定的条件，则应减小信赖域的范围，再重新求解。LMA 最早提出是用来解决最小二乘法曲线拟合的优化问题的。

Levenberg-Marquardt 算法，也称为衰减最小二乘法（damped least-squares method），该算法的损失函数采用平方误差和的形式。该算法的执行也不需要

计算具体的海森矩阵，它仅仅只是使用梯度向量和雅克比矩阵（Jacobian matrix）。

这里解释一下雅克比矩阵（Jacobian matrix）：Jacobi 矩阵实际上是向量值函数的梯度矩阵，假设F:Rn→Rm 是一个从n维欧氏空间转换到m维欧氏空间的函数。

这个函数由m个实函数组成:

这些函数的偏导数(如果存在)可以组成一个m行n列的矩阵(m by n)，这就是所谓的雅可比矩阵：

该算法的损失函数如下方程式所示为平方误差和的形式：

在上式中，m是数据集样本的数量。
我们可以定义损失函数的雅可比矩阵以误差对参数的偏导数为元素，公式如下：

其中m是数据集样本的数量，n是神经网络的参数向量。那么雅可比矩阵是m*n阶矩阵，损失函数的梯度向量就可以按如下计算出来：

E在这里是所有误差项的向量。
最终，我们可以利用一下表达式逼近海森矩阵：

其中λ为衰减因子，它确保了海森矩阵的证定性（positiveness），I是单位矩阵。下面的表达式定义了Levenberg-Marquardt 算法中参数的更新和优化过程：

当衰减参数λ为0时，Levenberg-Marquardt 算法就是使用海塞矩阵逼近值的牛顿法。而当λ很大时，该算法就近似于采用很小学习速率的梯度下降法。如果进行迭代导致了

损失函数上升，衰减因子λ就会下降，从而Levenberg-Marquardt 算法更接近于牛顿法。该过程经常用于加速收敛到极小值点。

第一步就是计算损失函数、梯度和海森矩阵逼近值，随后再在每次迭代降低损失中调整衰减参数。

Levenberg-Marquardt 算法是为平方误差和函数所定制的。这就让使用这种误差度量的神经网络训练地十分迅速。但是其缺点是：

• 1. 它不能用于平方根误差或交叉熵误差（cross entropy error）等函数。
• 2. 该算法与正则项不兼容。
• 3. 对于大型数据集或神经网络，雅可比矩阵变得十分巨大，因此也需要大量的内存。
• 在大型数据集或神经网络中并不推荐采用Levenberg-Marquardt 算法。

import numpy as np
def function(p):
    r = np.array([10 * (p[1] - p[0] ** 2), (1 - p[0])])
    fp = np.dot(np.transpose(r), r)  # = 100*(p[1]-p[0]**2)**2 + (1-p[0])**2
    J = (np.array([[-20 * p[0], 10], [-1, 0]]))
    grad = np.dot(np.transpose(J), np.transpose(r))
    return fp, r, grad, J


def lm(p0, tol=10 ** (-5), maxits=100):
    nvars = np.shape(p0)[0]
    nu = 0.01
    p = p0
    fp, r, grad, J = function(p)
    e = sum(np.dot(np.transpose(r), r))
    nits = 0
    while nits < maxits and np.linalg.norm(grad) > tol:
        nits += 1
        fp, r, grad, J = function(p)
        H = np.dot(np.transpose(J), J) + nu * np.eye(nvars)

        pnew = np.zeros(np.shape(p))
        nits2 = 0
        while (p != pnew).all() and nits2 < maxits:
            nits2 += 1
            dp, resid, rank, s = np.linalg.lstsq(H, grad)
            pnew = p - dp
            fpnew, rnew, gradnew, Jnew = function(pnew)
            enew = np.sum(np.dot(np.transpose(rnew), rnew))
            rho = np.linalg.norm(np.dot(np.transpose(r), r) - np.dot(np.transpose(rnew), rnew))
            rho /= np.linalg.norm(np.dot(np.transpose(grad), pnew - p))

            if rho > 0:
                update = 1
                p = pnew
                e = enew
                if rho > 0.25:
                    nu = nu / 10
            else:
                nu = nu * 10
                update = 0
        print(fp, p, e, np.linalg.norm(grad), nu)
        
p0 = np.array([-1.92, 2])
lm(p0)

17.共轭梯度法（Conjugate gradient）

共轭方向法是介于最速下降法和Newton法之间的一种方法。克服了最速下降法的锯齿现象，从而提高了收敛速度；同时，共轭方向法的迭代公式比较简单，不必求目标函数的Hesse矩阵，比Newton法减少了计算量和存储量。是一种比较实用而且有效的方法。

共轭梯度法是共轭方向法的一种，意思是搜索方向都互相共轭

共轭向量及其性质

定义1：(共轭方向) Q是n×n 对称正定矩阵，若 nn 维向量空间中的非零向量 p0，p1，⋅⋅⋅，pm−1p0，p1，···，pm−1 满足

则称p0，p1，⋅⋅⋅，pm−1p0，p1，···，pm−1 是Q共轭向量或称向量p0，p1，⋅⋅⋅，pm−1p0，p1，···，pm−1是Q共轭的(简称共轭)，称p0，p1，⋅⋅⋅，pm−1p0，p1，···，pm−1的方向是Q共轭方向。
当Q=I(单位矩阵)Q=I(单位矩阵)时，公式(1)变为

即向量p0，p1，⋅⋅⋅，pm−1互相正交。由此可见，正交是共轭的一种特殊情况，共轭是正交的推广。
定理1： 若非零向量p0，p1，⋅⋅⋅，pm−1是Q共轭的，则他们是线性无关的。
推论1： 在n维的向量空间中，非零的共轭向量个数不超过n。
定义2： 设p0，p1，⋅⋅⋅，pm−1是Rn中线性无关向量，x0∈Rⁿ，Pⁱ∈Rⁿ那么由形式为

z的向量构成的集合为由点x0和向量p0，p1，⋅⋅⋅，pm−1所生成的线性流行。记为L[x0;p0,p1,···,pm−1]。

基本思想

在考虑普通函数之前，我们首先用2元正定二次函数进行讲解。首先考虑如下的正定二次函数

要求的目标函数f(x)的最优值，根据最速下降法的思想，我们首先选定一个初始点x0x0，然后沿着该点的最速下降方向p0=−∇f(x0)做直线搜索，得到点x1，由最速下降法的性质(可以想一下最速下降法中的锯齿现象)可知

在第二次迭代过程中我们不适用−∇f(x1)作为这次迭代的搜索方向，我们想直在第二次迭代之后能直接到达最优点x^*（x^*表示是当前迭代找到的最优点），那么这次的迭代方向p1应该满足什么条件呢？
首先根据最速下降法的迭代公式我们有

其中t₁是最优步长因子,显然在未到达最优点x^*之前，t₁是不等于0的。对目标函数求梯度，有梯度函数

对于最优点x^*，我们有

即

在上式两边同乘以P₀^T，由于公式(4)，并且t₁≠0我们可以得到

 

由公式(7)我们知道，p0p0和p1p1为QQ的共轭向量。 
现在我们假设:

 

 

 

综上所述，对于n元正定二次函数，我们可以从任意点出发，然后沿着这n个共轭方向最多做n次直线搜索(原因将在下一小节进行讲解)，就可以求的目标函数的最优点。

优化模型建立

虽然梯度下降法的每一步都是朝着局部最优的方向前进的，但是它在不同的迭代轮数中会选择非常近似的方向，说明这个方向的误差并没通过一次更新方向和步长更新完，在这个方向上还存在误差，因此参数更新的轨迹是锯齿状。共轭梯度法的思想是，选择一个优化方向后，本次选择的步长能够将这个方向的误差更新完，在以后的优化更新过程中不再需要朝这个方向更新了。由于每次将一个方向优化到了极小，后面的优化过程将不再影响之前优化方向上的极小值，所以理论上对N维问题求极小只用对N个方向都求出极小就行了。为了不影响之前优化方向上的更新量，需要每次优化方向共轭正交。假定每一步的优化方向用p_k表示，可得共轭正交

由此可得，每一步优化后，当前的误差和刚才的优化方向共轭正交。

若为N维空间优化问题，则每次优化方向可以组成这个空间中的一组基底。P={p₁,p2,…,p_n} 。

算法只需要解决两个问题：

1. 优化方向
2. 优化步长

优化方向确定

优化步长的选取

三个推论

最终简化

算法在三中基本推导完毕，但是在工程应用中如果每次进行p_k的正交化需要对之前所有的优化向量求解β，现简化如下：

最终的推导结束。整理为如下的伪代码（用于正定二次函数的共轭梯度法）

对于正定二次函数(3)，从任意点出发，顺次沿nn个QQ共轭方向做直线搜索，最多经过nn次就可以到达极小点。 
下面是对于正定二次函数的共轭方向法的算法描述。

import numpy as np
np.random.seed(12)
A = np.random.randint(8,18,size=(2,2))
A = A+A.T #把A变成对称矩阵
B = np.random.randint(-5,8,size=(2,1))

def func(X):
    return 0.5*np.dot(X.T,np.dot(A,X)) - np.dot(B.T,X)

def d_fun(X):
    return np.dot(A,X) - B

#最速下降法
def gd(x,func,d_f,lr=0.005,e=0.0001):
    init_x = x
    init_y = func(init_x)
    ls = 1
    x_point = init_x
    y_point = init_y

    while ls > e:
        init_x -= lr*d_f(init_x)
        new_y = func(init_x)
        ls = np.abs(new_y-init_y)[0,0]
        init_y = new_y
        x_point = np.hstack((x_point,init_x))
        y_point = np.hstack((y_point,init_y))
    return x_point.T,y_point.T

import copy
#共轭下降法
def cad(init_x,d_f):
    print(init_x)
    x_0 = copy.deepcopy(init_x)
    r_0 = -d_f(x_0)
    p_0 = r_0
    r_k = r_0
    p_k = p_0
    x_k = x_0
    k = 1
    n = len(x_0)
    print(n)
    x_point = x_k
    while True:
        alpha_k = np.dot(r_k.T,r_k)[0,0]/np.dot(p_k.T,np.dot(A,p_k))[0,0]
        x_k += alpha_k*p_k
        r_k2 = r_k - alpha_k*np.dot(A,p_k)
        x_point = np.hstack((x_point,x_k))
        if k >= n:
            break
        beta_k = np.dot(r_k2.T,r_k2)[0,0]/np.dot(r_k.T,r_k)[0,0]
        p_k = r_k2 + beta_k*p_k
        r_k = r_k2
        k += 1
    x_point[:,0] = init_x.reshape(-1)
    return x_point.T

x_p = cad(np.array([[0],[0]],dtype=np.float),d_fun)
x_list,y_list = gd(np.array([[0],[0]],dtype=np.float),func,d_fun)
x1_max, x2_max = np.max(x_list,axis=0)
x1_min, x2_min = np.min(x_list,axis=0)
x1 = np.linspace(x1_min-1,x1_max+1,100)
x2 = np.linspace(x2_min-1, x2_max+1,100)
X1,X2 = np.meshgrid(x1,x2)
X = np.array(list(zip(X1.flatten(),X2.flatten())))
Y = np.array(list(map(lambda x:func(x.reshape(2,1)),X))).reshape(*X1.shape)
contour = plt.contour(X1,X2,Y,levels=list(np.linspace(np.min(y_list)+0.2,np.max(y_list),5)),colors='k')
plt.clabel(contour)
plt.scatter(x_p[:,0],x_p[:,1],c='g')
plt.plot(x_p[:,0],x_p[:,1],c='g',label='cad')
plt.plot(x_list[:,0],x_list[:,1],'r-',label='gd')
plt.legend()
plt.show()

在算法描述的第4步中，只是说要求得共轭方向，并没有说明如何求共轭方向，不同求共轭方向的方法对应了不同的共轭方向法。所以可以按照上述方法，依次构造出共轭向量，最多经过n次迭代就能找到最优点x^*。 
共轭梯度法的算法描述如下。

用于非二次函数的共轭梯度法

用于一般函数的Fletcher-Reeves共轭梯度法描述如下。

总结

共轭梯度法是一种典型的共轭方向法，它的搜索方向是负梯度方向和上一次搜索方向的一个组合

关于β_k-1，有两种常用的计算公式

（PRP公式）

（PRP公式）

预处理共轭法改善了G的条件数，使算法的收敛速度加快

预处理的方法是寻找一个非奇异矩阵C，使得C^-TGC^-1的条件数小于G的条件数

再开始共轭梯度法是满足某一条件后重新使用最速下降方向作为搜索方向，这个条件包括迭代n步，或共轭梯度方向是上升方向

FR共轭梯度法：

# coding=utf-8
# FR-CG共轭梯度法
from linear_search.Function import *
from linear_search.wolfe import *
import numpy as np


def conjugate(_f, _x):
    fun = Function(_f)
    x = array(_x)
    d = -fun.grad(x)
    while fun.norm(x) > 0.01:
        alpha = wolfe(_f, x, d)
        g = np.mat(fun.grad(x))
        beta = 1 / np.dot(g, g.T)
        x = x + alpha * d
        g = np.mat(fun.grad(x))
        beta = beta * np.dot(g, g.T)
        d = array(-g + beta * d)[0]
    return x

 再开始共轭梯度法：

# coding=utf-8
# 再开始共轭梯度法
from linear_search.Function import *
from linear_search.wolfe import *
import numpy as np


def restart_conjugate(_f, _x, n):
    fun = Function(_f)
    x = array(_x)
    while True:
        d = -fun.grad(x)
        k=0
        if(np.linalg.norm(d)<0.01):
            break
        while fun.norm(x) > 0.01:
            g = np.mat(fun.grad(x))

            alpha = wolfe(_f, x, d)
            x = x + alpha * d
            k=k+1

            g1= np.mat(fun.grad(x))

            if np.dot(g1,g.T)/np.dot(g1,g1.T)>0.1 or k>=n:
                if np.linalg.norm(g1)<0.01:
                    return x
                break
            else:
                beta = np.dot(g1, g1.T) / np.dot(g, g.T)
                d = array(-g + beta * d)[0]
                if np.dot(mat(d),g1.T)>0:
                    break
    return x

其他的文档介绍

共轭梯度法是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点，该算法希望能加速梯度下降的收敛速度，同时避免使用海森矩阵进行求值、存储和求逆获得必要的优化信息。共轭梯度法不仅是解决大型线性方程组最有用的方法之一，也是解大型非线性最优化最有效的算法之一。 在各种优化算法中，共轭梯度法是非常重要的一种。其优点是所需存储量小，具有步收敛性，稳定性高，而且不需要任何外来参数。

在共轭梯度训练算法中，因为是沿着共轭方向（conjugate directions）执行搜索的，所以通常该算法要比沿着梯度下降方向优化收敛得更迅速。共轭梯度法的训练方向

是海森矩阵共轭的。我们用d表示训练方向向量，然后从初始参数向量w0和初始训练方向向量d0=-g0开始，共轭梯度法所构建的训练方向序列为：

 

上式中，Y称之为共轭参数，并且有一些方法计算这个参数。两种常用的方法源自Fletcher、Reeves 和 Polak 、Ribiere。对于所有的共轭梯度算法，训练方向

周期性地重置为负梯度方向。参数通过下面的表达式得以更新和优化，通常学习速率η可使用单变量函数优化方法求得。

共轭梯度法的训练过程流程图如下表示，从图中我们可以看出来模型是通过第一次计算共轭梯度训练方向而优化参数的，然后再寻找适当的学习速率。

共轭梯度法已经证实其在神经网络中要比梯度下降法有效得多。并且由于共轭梯度法并没有要求使用海森矩阵，所以在大规模神经网络中其还可以做到很好的性能。

还是用最开始的目标函数      来编写共轭梯度法的优化代码：

import numpy as np

def Jacobian(x):
    # return array([.4*x[0],2*x[1]])
    return np.array([x[0], 0.4 * x[1], 1.2 * x[2]])

def Hessian(x):
    # return array([[.2,0],[0,1]])
    return np.array([[1, 0, 0], [0, 0.4, 0], [0, 0, 1.2]])


def CG(x0):
    i = 0
    k = 0

    r = -Jacobian(x0)
    p = r

    betaTop = np.dot(r.transpose(), r)
    beta0 = betaTop

    iMax = 3
    epsilon = 10 ** (-2)
    jMax = 5

    # Restart every nDim iterations
    nRestart = np.shape(x0)[0]
    x = x0

    while i < iMax and betaTop > epsilon ** 2 * beta0:
        j = 0
        dp = np.dot(p.transpose(), p)
        alpha = (epsilon + 1) ** 2
        # Newton-Raphson iteration
        while j < jMax and alpha ** 2 * dp > epsilon ** 2:
            # Line search
            alpha = np.dot(Jacobian(x).transpose(), p) / (np.dot(p.transpose(), np.dot(Hessian(x), p)))
            print("N-R", x, alpha, p)

            x = x + alpha * p
            j += 1
        print(x)
        # Now construct beta
        r = -Jacobian(x)
        print("r: ", r)

        betaBottom = betaTop
        betaTop = np.dot(r.transpose(), r)
        beta = betaTop / betaBottom
        print("Beta: ", beta)

        # Update the estimate
        p = r + beta * p
        print("p: ", p)

        print("----")
        k += 1

        if k == nRestart or np.dot(r.transpose(), p) <= 0:
            p = r
            k = 0
            print("Restarting")
            
        i += 1
    print(x)
    
x0 = np.array([-2, 2, -2])
CG(x0)

总结

提速可以归纳为以下几个方面： 

• 使用momentum来保持前进方向(velocity)； 
• 为每一维参数设定不同的学习率：在梯度连续性强的方向上加速前进； 
• 用历史迭代的平均值归一化学习率：突出稀有的梯度；

综上分析，可以得出如下几个结论：

• AdaGrad、RMSProp、AdaDelta和Adam几个优化算法，目标函数自变量中每个元素都分别拥有自己的学习率；
• AdaGrad目标函数自变量中各个元素的学习率只能保持下降或者不变，因此当学习率在迭代早期降得较快且当前解依然不佳时，由于后期学习率过小，可能较难找到一个有用的解；
• RMSProp和AdaDelta算法都是解决AdaGrad上述缺点的改进版本，本质思想都是利用最近的时间步的小批量随机梯度平方项的加权平均来降低学习率，从而使得学习率不是单调递减的（当最近梯度都较小的时候能够变大）。不同的是，RMSProp算法还是保留了传统的学习率超参数，可以显式指定。而AdaDelta算法没有显式的学习率超参数，而是通过ΔxΔx做运算来间接代替学习率；
• Adam算法可以看成是RMSProp算法和动量法的结合。

 

如图，收敛速度最慢的是梯度下降算法，但该算法同时也只要求最少内存。相反，Levenberg-Marquardt 算法可能是收敛速度最快的，但其同时也要求最多的内存。

比较折中的方法是拟牛顿法。
总而言之，如果我们的神经网络有数万参数，为了节约内存，我们可以使用梯度下降或共轭梯度法。如果我们需要训练多个神经网络，并且每个神经网络都只有数百参数、数千样本，那么我们可以考虑 Levenberg-Marquardt 算法。而其余的情况，拟牛顿法都能很好地应对。

二.在线算法（OGD =》 FOBOS =》 RDA =》 FTRL）

在线学习算法相关背景

在线学习算法强调的是训练的实时性，面向流式数据，每次训练不使用全量数据，而是以之前训练好的参数为基础，每次利用一个样本更新一次模型，是属于增量学习的一部分，从而快速更新模型，提高模型的时效性。

在线学习算法与离线学习关注点的不同

统计学习的先验假设是数据存在自己一定的分布，我们的目的是寻找与实际分布距离最小的策略来泛化未知的结果。数据由真实模型产生，如果能有无限数据、并在包含有真实模型的空间里求解，也许我们能算出真是模 型。但实际上我们只有有限的有噪音的数据，这又限制我们只能使用相对简单的模型。所以，理想的算法是能够用不多的数据来得到一个不错的模型。

离线学习算法强调的是数据中不断训练已有的模型，不断迭代达到真实的数据分布，在训练过程中，所有数据都是可见的，目标是一定的（就是最终的那个真实分布），其中可以采用不同的更新策略和采样策略，所以有了批量梯度下降和随机梯度下降等算法。

但是，在线学习算法的关注点是不同的，在线学习算法的限制条件是只能看到过去的数据和当前的数据，看不到未来的数据，所以我们训练的策略可能是任意的数据分布，通过不断的训练可能颠覆之前的训练结果，所以在线学习甚至是为了破坏我们之前的策略精心设计的。

在线学习关注点在于，追求对所知道的所有知识所能设计最优的策略，那么同这个最优的策略的差距成为后悔（regret）:后悔没有从一开始就选择这个策略，当然，我们希望的是，随着时间的增加，这个差异会不断的变小。因为我们不对数据进行任何假设，所以策略是否完美并不是我们关心的（比如回答所有问题），我们追求的是，没有后悔（no-regret）

batch模式和delta模式

梯度下降可以分成两种模式，batch模式和delta模式。batch模式的时效性比delta模式要低一些。分析一下batch模式，比如昨天及昨天的数据训练成了模型M，那么今天的每一条训练数据在训练过程中都会更新一次模型M，从而生成今天的模型M1。

batch学习可以认为是离线学习算法，强调的是每次训练都需要使用全量的样本，因而可能会面临数据量过大的问题。一般进行多轮迭代来向最优解靠近。online learning没有多轮的概念，如果数据量不够或训练数据不够充分，通过copy多份同样的训练数据来模拟batch learning的多轮训练也是有效的方法。

delta模式可以认为是在线学习算法，没有多轮的概念，如果数据量不够或训练数据不够充分，通过copy多份同样的训练数据来模拟batch learning的多轮训练也是有效的方法。所以，OGD和SGD都属于在线学习算法，因为每次更新模型只用一个样本。SGD则每次只针对一个观测到的样本进行更新。通常情况下，SGD能够比GD“更快”地令 逼近最优值。当样本数 特别大的时候，SGD的优势更加明显，并且由于SGD针对观测到的“一条”样本更新 ，很适合进行增量计算，实现梯度下降的Online模式（OGD, OnlineGradient Descent）。

在线凸优化

现在做在线学习和CTR常常会用到逻辑回归（ Logistic Regression），而传统的批量（batch）算法无法有效地处理超大规模的数据集和在线数据流，google先后三年时间（2010年-2013年）从理论研究到实际工程化实现的FTRL（Follow-the-regularized-Leader）算法，在处理诸如逻辑回归之类的带非光滑正则化项（例如1范数，做模型复杂度控制和稀疏化）的凸优化问题上性能非常出色，据闻国内各大互联网公司都第一时间应用到了实际产品中，我们的系统也使用了该算法。这里对FTRL相关发展背景和工程实现的一些指导点做一些介绍，凸优化的理论细节不做详细介绍，感兴趣可以去查阅相应paper，相关paper列表会在文后附上。机器学习并非本人在校时的专业方向，不过在校期间积累的基础不算太差，而且很多东西也是相通的，钻研一下基本意思都还能搞明白。当然，有不准确的地方欢迎大家讨论指正。

本文主要会分三个部分介绍，如果对理论产生背景不感兴趣的话，可以直接看第3部分的工程实现（这一部分google13年那篇工程化的paper介绍得很详细）：

1. 相关背景：包括通用性的问题描述、批量算法、传统在线学习算法等
2. 简单介绍与FTRL关系比较密切的Truncated Gradient、FOBOS以及RDA（Regularized Dual Averaging）等算法
3. FTRL理论公式以及工程实现（对前因后果和理论方面不感兴趣的可以直接看这一小节的工程实现部分）

本章只是简单介绍，如果想要深入了解在线凸优化（OCO），强烈推荐阅读Elad Hazand的著作和Zinkevich的Paper

WM的理论证明可以参考Littlestone 94,Freund 99，虽然在上个世纪已经完成，但是将其理论拓展到一般的凸的函数还是在03年由Zinkevich完成的。

对于loss函数+正则化的结构风险最小化的优化问题（逻辑回归也是这种形式）有两种等价的描述形式，以1范数为例，分别是：

a、无约束优化形式的soft regularization formulation：

     

 

b、带约束项的凸优化问题convex constraint formulation：

批量（batch）算法

 批量算法中每次迭代对全体训练数据集进行计算（例如计算全局梯度），优点是精度和收敛还可以，缺点是无法有效处理大数据集（此时全局梯度计算代价太大），且没法应用于数据流做在线学习。这里分无约束优化形式和约束优化（与上面问题描述可以对应起来）两方面简单介绍一下一些传统批量算法。

• a、无约束优化形式：1、全局梯度下降，很常用的算法，就不细说了，每一步求一个目标函数的全局梯度，用非增学习率进行迭代；2、牛顿法（切线近似）、LBFGS（割线拟牛顿，用之前迭代结果近似Hessian黑塞矩阵的逆矩阵，BFGS似乎是几个人名的首字母的简称）等方法。牛顿和拟牛顿等方法一般对于光滑的正则约束项（例如2范数）效果很好，据说是求解2范数约束的逻辑回归类问题最好的方法，应用也比较广，但是当目标函数带L1非光滑、带不可微点的约束项后，牛顿类方法比较无力，理论上需要做修改。感兴趣的可以去查查无约束优化的相关数值计算的书，我也没有更深入研究相关细节，这里不做重点关注。
• b、不等式约束凸优化形式：1、传统的不等式约束优化算法内点法等；2、投影梯度下降（约束优化表示下），gt是subgradient，直观含义是每步迭代后，迭代结果可能位于约束集合之外，然后取该迭代结果在约束凸集合上的投影作为新的迭代结果（第二个公式中那个符号标识向X的投影）：

在线算法

如上所述，批量算法有自身的局限性，而在线学习算法的特点是：每来一个训练样本，就用该样本产生的loss和梯度对模型迭代一次，一个一个数据地进行训练，因此可以处理大数据量训练和在线训练。常用的有在线梯度下降（OGD）和随机梯度下降（SGD）等，本质思想是对上面【问题描述】中的未加和的单个数据的loss函数 L（w，zi)做梯度下降，因为每一步的方向并不是全局最优的，所以整体呈现出来的会是一个看似随机的下降路线。典型迭代公式如下：

这里使用混合正则化项：，例如可能是1范数与2范数强凸项的混合（后面会看到其实很多都是这种混合正则化的格式，而且是有一定直观含义的）。迭代公式中：gt是loss函数（单点的loss，未加和）的subgradient，与gt相加的那一项是混合正则化项中的第二项的梯度，投影集合C是约束空间（例如可能是1范数的约束空间），跟上面介绍的投影梯度下降类似的做法。

梯度下降类的方法的优点是精度确实不错，但是不足相关paper主要提到两点：

• 1、简单的在线梯度下降很难产生真正稀疏的解，稀疏性在机器学习中是很看重的事情，尤其我们做工程应用，稀疏的特征会大大减少predict时的内存和复杂度。这一点其实很容易理解，说白了，即便加入L1范数（L1范数能引入稀疏解的简单示例可以产看PRML那本书的第二章，我前面一篇blog的ppt里也大概提了），因为是浮点运算，训练出的w向量也很难出现绝对的零。到这里，大家可能会想说，那还不容易，当计算出的w对应维度的值很小时，我们就强制置为零不就稀疏了么。对的，其实不少人就是这么做的，后面的Truncated Gradient和FOBOS都是类似思想的应用；
• 2、对于不可微点的迭代会存在一些问题，具体有什么问题，有一篇paper是这么说的：the iterates of the subgradient method are very rarely at the points of non-differentiability。我前后看了半天也没看明白，有熟悉的同学可以指导一下。

三种做稀疏解的途径

1）、简单加入L1范数

• 局限如上面所提，a+b两个float数很难绝对等于零，无法产生真正稀疏的特征权重

2）、在1范数的基础上做截断，最直观没技术含量的思路，那就设定一个阈值，做截断来保证稀疏，可以结合L1范数

• 简单截断方法，每online训练K个数据截断一次，对OGD的迭代结果，每K步做一次截断置零：

但是简单截断方法有问题：权重小，可能是确实是无用特征，还或者可能是该特征才刚被更新一次（例如训练刚开始的阶段、或者训练数据中包含该特征的样本数本来就很少），另外，简单rounding技术太aggressive了，可能会破坏在线训练算法的理论完备性。

• 简单截断基础上，不太aggressive的Truncated gradient（TG） （09年的工作），其实后面的FOBOS也可以归为这一类：

参数λ 和g_i控制模型稀疏，值越大越稀疏。 TG和简单截断法的区别在于 

3）、Black-box wrapper approaches：

• 黑盒的方法去除一些特征，然后重新训练的看被消去的特征是否有效。
• 需要在数据集上对算法跑多次，所以不太实用

下面会提一下FOBOS（Forward-Backward Splitting method ，其实应该叫FOBAS的，历史原因）以及RDA，因为后面的FTRL其实相当于综合了这两种算法的优点。

原始的OGD

优化函数：

其中x表示模型的参数，xi表示模型的第i个参数。

L(x)是逻辑回归的似然函数的负对数。

最右边那一项是L1正则项。

可以考虑加上L2正则来防止过拟合。

假设模型参数为x，第t个样本的特征表示为vt

则样本label为1的概率为：

该样本对应的对数损失函数为：

令

    

则lt关于xt的求导过程如下：

上式即为第t个样本的梯度g_t。

原始的OGD使用的迭代公式为：

其中et表示非增的学习率，也叫做步长，

OGD的算法比较简单，不够好，即使是加了L1正则化也可能不产生稀疏解。

FOBOS

FOBOS (Forward-Backward Splitting)是由John Duchi和Yoram Singer提出的^[11]。从全称上来看，该方法应该叫FOBAS，但是由于一开始作者管这种方法叫FOLOS（Forward Looking Subgradients），为了减少读者的困扰，作者干脆只修改一个字母，叫FOBOS。

该算法又称为前后项算法，其权重的更新方式分为两个部分 ,可以看作truncated gradient的一种特殊形式

基本思想：跟projected subgradient方法类似，不过将每一个数据的迭代过程，分解成一个经验损失梯度下降迭代和一个最优化问题。分解出的第二个最优化问题，有两项：第一项2范数那一项表示不能离第一步loss损失迭代结果太远，第二项是正则化项，用来限定模型复杂度抑制过拟合和做稀疏化等。这个最优化问题有一些特殊的性质，从而保证了最终结果的稀疏性和理论上的完备，具体细节感兴趣的可以查看对应paper。我这里更多关注直观含义和工程实现，忽略理论方面的内容。

 

第一步为标准的梯度下降，然后将结果进行正则化的微调，满足最速下降切保证了一定的稀疏性。其中ψ(W)ψ(W)是正则项可以对权重进行约束，常用的选择为L1，则对应的方法为L1-FOBOS。

AOGD

此算法在网上极少相关文档，还不能详细介绍。

其迭代伪代码为：

RDA

RDA（Regularized dual averaging），微软10年的工作

• 非梯度下降类方法，属于更加通用的一个primal-dual algorithmic schema的一个应用
• 克服了SGD类方法所欠缺的exploiting problem structure，especially for problems with explicit regularization。
• 能够更好地在精度和稀疏性之间做trade-off

RDA（正则对偶平均）是微软的研究成果，其权重更新策略为

权重更新包括三个部分第一部分线性加权，即对历史梯度进行加权平均；第二部分正则项部分对特征进行稀疏化；第三部分严格递增序列相当于额外的正则项。
实验证明该方法能够产生较好的稀疏和精度。

FTRL

主要用于CTR预测的在线训练，成千上万维度导致大量稀疏特征。一般希望模型参数更加稀疏，但是简单的L1正则无法真正做到稀疏，一些梯度截断方法（TG）的提出就是为了解决这个问题，在这其中FTRL是兼备精度和稀疏性的在线学习方法。FTRL的基本思想是将接近于0的梯度直接置零，计算时直接跳过以减少计算量。

在上文中我们从原理上定性比较了L1-FOBOS和L1-RDA在稀疏性上的表现。有实验证明，L1-FOBOS这一类基于梯度下降的方法有比较高的精度，但是L1-RDA却能在损失一定精度的情况下产生更好的稀疏性。那么这两者的优点能不能在一个算法上体现出来？这就是FTRL要解决的问题。

FTRL（Follow the Regularized Leader）是由Google的H. Brendan McMahan在2010年提出的[1]，后来在2011年发表了一篇关于FTRL和AOGD、FOBOS、RDA比较的论文[2]，2013年又和Gary Holt, D. Sculley, Michael Young等人发表了一篇关于FTRL工程化实现的论文[3]。

本节我们首先从形式上将L1-FOBOS和L1-RDA进行统一，然后介绍从这种统一形式产生FTRL算法，以及介绍FTRL算法工程化实现的算法逻辑。

1.  L1-FOBOS和L1-RDA在形式上的统一

L1-FOBOS在形式上，每次迭代都可以表示为（这里我们令，是一个随变化的非增正序列）：

 

 

把这两个公式合并到一起，有：

 

 

通过这个公式很难直接求出的解析解，但是我们可以按维度将其分解为个独立的最优化步骤：

由于与变量无关，因此上式可以等价于：

 

 

再将这N个独立最优化子步骤合并，那么L1-FOBOS可以写作：

 

 

而对于L1-RDA，我们可以写作：

 

 

这里；如果令，，上面两个式子可以写作：

 

   公式(1)     公式(2)

 

需要注意，与论文[2]中的Table 1不同，我们并没有将L1-FOBOS也写成累加梯度的形式。

比较公式(1)和公式(2)，可以看出L1-FOBOS和L1-RDA的区别在于：(1) 前者对计算的是累加梯度以及L1正则项只考虑当前模的贡献，而后者采用了累加的处理方式；(2) 前者的第三项限制的变化不能离已迭代过的解太远，而后者则限制不能离点太远。

2. FTRL算法原理

FTRL综合考虑了FOBOS和RDA对于正则项和W限制的区别，其特征权重的更新公式为：

 

   公式(3)

 

注意，公式(3)中出现了L2正则项，在论文[2]的公式中并没有这一项，但是在其2013年发表的FTRL工程化实现的论文[3]却使用到了L2正则项。事实上该项的引入并不影响FRTL的稀疏性，后面的推导过程会显示这一点。L2正则项的引入仅仅相当于对最优化过程多了一个约束，使得结果求解结果更加“平滑”。

公式(3)看上去很复杂，更新特征权重貌似非常困难的样子。不妨将其进行改写，将最后一项展开，等价于求下面这样一个最优化问题：

 

 

上式中最后一项相对于来说是一个常数项，并且令，上式等价于：

 

 

针对特征权重的各个维度将其拆解成N个独立的标量最小化问题：

 

 

到这里，我们遇到了与上一篇RDA中类似的优化问题，用相同的分析方法可以得到：

 

   公式(4)

 

从公式(4)可以看出，引入L2正则化并没有对FTRL结果的稀疏性产生任何影响。

3.权重更新策略

工程优化：

1. 预测的内存方面：L1范式加策略，训练结果w很稀疏，在用w做predict的时候节省了内存，很直观

2. 训练的内存方面：

1）在线丢弃训练数据中很少出现的特征，即稀疏特征处理：

• Possion Inclusion(通过一定概率扔掉特征，对某一维度特征所来的训练样本，以p的概率接受并更新模型)
• BloomFilter Iclusion(通过碰撞优化特征，用bloom filter从概率上做某一特征出现k次才更新)

2）浮点数重新编码

• 特征权重不需要用32bit或64bit的浮点数存储，存储浪费空间
• 16bit encoding，但是要注意处理rounding技术对regret带来的影响

3）训练若干相似model

• 对同一份数据序列，同时训练多个相似的model
• 这些model有各自独享的一些feature，也有一些共享的feature
• 出发点：有的特征维度可以是各个模型独享的，而有的各个模型共享的特征，可以用同样的数据训练

4）Single Value Structure（据说有公司已经在实际中这么搞，大数据量下也能够保证不错的auc）

• 多个model公用一个feature存储（例如放到cbase或redis中），各个model都更新这个共有的feature结构
• 对于某一个model，对于他所训练的特征向量的某一维，直接计算一个迭代结果并与旧值做一个平均

5） 使用正负样本的数目来计算梯度的和（所有的model具有同样的N和P） 

６) Subsampling Training Data

• 在实际中，CTR远小于50%，所以正样本更加有价值。通过对训练数据集进行subsampling，可以大大减小训练数据集的大小
• 正样本全部采（至少有一个广告被点击的query数据），负样本使用一个比例r采样（完全没有广告被点击的query数据）。但是直接在这种采样上进行训练，会导致比较大的biased prediction
• 解决办法：训练的时候，对样本再乘一个权重。权重直接乘到loss上面，从而梯度也会乘以这个权重，假设采样概率为r，此时相当于每个样本都有了权重w_t（正样本权重为1，负样本权重为1/r），同时相当于梯度也会乘上该系数，不会改变最终的损失。
• El(wt)=stwtl(wt)+(1−st)0=l(wt)，其中st=1/wt为该样本被选中的概率
• 先采样减少负样本数目，在训练的时候再用权重弥补负样本，非常不错的想法。

4. Per-Coordinate Learning Rates

前面介绍了FTRL的基本推导，但是这里还有一个问题是一直没有被讨论到的：关于学习率的选择和计算。事实上在FTRL中，每个维度上的学习率都是单独考虑的（Per-Coordinate Learning Rates）。

在一个标准的OGD里面使用的是一个全局的学习率策略，这个策略保证了学习率是一个正的非增长序列，对于每一个特征维度都是一样的。

考虑特征维度的变化率：如果特征1比特征2的变化更快，那么在维度1上的学习率应该下降得更快。我们很容易就可以想到可以用某个维度上梯度分量来反映这种变化率。在FTRL中，维度i上的学习率是这样计算的：

 

 公式(5)

 

由于，所以公式(4)中有，这里的和是需要输入的参数，公式(4)中学习率写成累加的形式，是为了方便理解后面FTRL的迭代计算逻辑。

5. FTRL算法逻辑与工程实现

大家对上面那一大坨前因后果和公式都不感兴趣，ok，没关系，google非常贴心地在13年给出了一篇工程性很强的paper，其实大部分公司使用FTRL的，根本不会关心上面那一大段东西，直接按着伪代码写，调调参，看结果很不错就可以了。我们公司开始就是这么搞的，哈哈，不过人总是要有点儿好奇心的不是，深究一下前因后果和基本的理论公式感觉还是挺不同的。

逻辑回归下的per-coordinate FTRL_Proximal的伪代码如下，在公式表达的基础上做了一些变换和实现上的trick，细节paper里有，大家在自己做实现的时候，可以在实际数据集上再并行加加速：

算法流程：

根据前面的学习，可以使用LR作为基本学习器，使用FTRL作为在线最优化的方法来获取LR的权重系数，从而达到在不损失精度的前提下获得稀疏解的目标。

工程化实现的几个核心点是：

• 梯度计算代码
• 损失函数计算代码
• 权重更新计算代码

算法代码实现如下所示，这里只给出了核心部分代码，主要做了以下优化:

• 修复原代码更新梯度时候的逻辑错误。
• 支持pypy3加速，5.6G的训练数据，使用Mac单机可以9分钟跑完一个模型，一个epoch之后的logloss结果为0.3916；这对于刷竞赛或者实际工程中模型训练已经是比较理想的性能了。
• 使用16位浮点数保存权重数据，降低模型文件大小。
• 使用json保存模型非0参数w,z,n，进一步压缩线上使用的时候占据的内存空间，可以进一步考虑压缩模型文件大小，比如使用bson编码+gzip压缩后保存等。

from datetime import datetime
from csv import DictReader
from math import exp, log, sqrt
import gzip
import random
import json
import argparse


class FTRLProximal(object):
    """  FTRL Proximal engineer project with logistic regression  Reference:  https://static.googleusercontent.com/media/research.google.com/zh-CN//pubs/archive/41159.pdf   """

    def __init__(self, alpha, beta, L1, L2, D,
                 interaction=False, dropout=1.0,
                 dayfeature=True,
                 device_counters=False):

        # parameters
        self.alpha = alpha
        self.beta = beta
        self.L1 = L1
        self.L2 = L2
        self.dayfeature = dayfeature
        self.device_counters = device_counters

        # feature related parameters
        self.D = D
        self.interaction = interaction
        self.dropout = dropout

        # model
        self.n = [0.] * D
        self.z = [0.] * D
        self.w = [0.] * D

    def _indices(self, x):
        '''  A helper generator that yields the indices in x  The purpose of this generator is to make the following  code a bit cleaner when doing feature interaction.  '''

        for i in x:
            yield i

        if self.interaction:
            D = self.D
            L = len(x)
            for i in range(1, L):  # skip bias term, so we start at 1
                for j in range(i + 1, L):
                    # one-hot encode interactions with hash trick
                    yield abs(hash(str(x[i]) + '_' + str(x[j]))) % D

    def predict(self, x, dropped=None):
        """  use x and computed weight to get predict  :param x:  :param dropped:  :return:  """
        # wTx is the inner product of w and x
        wTx = 0.
        for j, i in enumerate(self._indices(x)):

            if dropped is not None and dropped[j]:
                continue

            wTx += self.w[i]

        if dropped is not None:
            wTx /= self.dropout

        # bounded sigmoid function, this is the probability estimation
        return 1. / (1. + exp(-max(min(wTx, 35.), -35.)))

    def update(self, x, y):
        """  update weight and coordinate learning rate based on x and y  :param x:  :param y:  :return:  """

        ind = [i for i in self._indices(x)]

        if self.dropout == 1:
            dropped = None
        else:
            dropped = [random.random() > self.dropout for i in range(0, len(ind))]

        p = self.predict(x, dropped)

        # gradient under logloss
        g = p - y

        # update z and n
        for j, i in enumerate(ind):

            # implement dropout as overfitting prevention
            if dropped is not None and dropped[j]:
                continue

            g_i = g * i
            sigma = (sqrt(self.n[i] + g_i * g_i) - sqrt(self.n[i])) / self.alpha
            self.z[i] += g_i - sigma * self.w[i]
            self.n[i] += g_i * g_i

            sign = -1. if self.z[i] < 0 else 1.  # get sign of z[i]

            # build w on the fly using z and n, hence the name - lazy weights -
            if sign * self.z[i] <= self.L1:
                # w[i] vanishes due to L1 regularization
                self.w[i] = 0.
            else:
                # apply prediction time L1, L2 regularization to z and get
                self.w[i] = (sign * self.L1 - self.z[i]) \
                            / ((self.beta + sqrt(self.n[i])) / self.alpha + self.L2)

    def save_model(self, save_file):
        """  保存weight数据到本地  :param save_file:  :return:  """
        with open(save_file, "w") as f:
            w = {k: v for k, v in enumerate(self.w) if v != 0}
            z = {k: v for k, v in enumerate(self.z) if v != 0}
            n = {k: v for k, v in enumerate(self.n) if v != 0}
            data = {
                'w': w,
                'z': z,
                'n': n
            }
            json.dump(data, f)

    def load_weight(self, model_file, D):
        """  loada weight data  :param model_file:  :return:  """
        with open(model_file, "r") as f:
            data = json.load(f)
            self.w = data.get('w', [0.] * D)
            self.z = data.get('z', [0.] * D)
            self.n = data.get('n', [0.] * D)

    @staticmethod
    def loss(y, y_pred):
        """  log loss for LR model  :param y:  :param y_pred:  :return:  """
        p = max(min(y_pred, 1. - 10e-15), 10e-15)
        return -log(p) if y == 1. else -log(1. - p)


def data(f_train, D, dayfilter=None, dayfeature=True, counters=False):
    ''' GENERATOR: Apply hash-trick to the original csv row  and for simplicity, we one-hot-encode everything   INPUT:  path: path to training or testing file  D: the max index that we can hash to   YIELDS:  ID: id of the instance, mainly useless  x: a list of hashed and one-hot-encoded 'indices'  we only need the index since all values are either 0 or 1  y: y = 1 if we have a click, else we have y = 0  '''

    device_ip_counter = {}
    device_id_counter = {}

    for t, row in enumerate(DictReader(f_train)):
        # process id
        ID = row['id']
        del row['id']

        # process clicks
        y = 0.
        if 'click' in row:
            if row['click'] == '1':
                y = 1.
            del row['click']

        # turn hour really into hour, it was originally YYMMDDHH

        date = row['hour'][0:6]
        row['hour'] = row['hour'][6:]

        if dayfilter != None and not date in dayfilter:
            continue

        if dayfeature:
            # extract date
            row['wd'] = str(int(date) % 7)
            row['wd_hour'] = "%s_%s" % (row['wd'], row['hour'])

        if counters:
            d_ip = row['device_ip']
            d_id = row["device_id"]
            try:
                device_ip_counter[d_ip] += 1
                device_id_counter[d_id] += 1
            except KeyError:
                device_ip_counter[d_ip] = 1
                device_id_counter[d_id] = 1
            row["ipc"] = str(min(device_ip_counter[d_ip], 8))
            row["idc"] = str(min(device_id_counter[d_id], 8))

        # build x
        x = [0]  # 0 is the index of the bias term
        for key in row:
            value = row[key]
            # one-hot encode everything with hash trick
            index = abs(hash(key + '_' + value)) % D
            x.append(index)
        yield t, ID, x, y

import tensorflow as tf
tf.compat.v1.train.FtrlOptimizer(learning_rate)

总结

实际的互联网广告应用需要的是快速地进行model的更新。为了保证快速的更新，训练样本是一条一条地过来的，每来一个样本，model的参数对这个样本进行一次迭代，从而保证了model的及时更新，这种方法叫做OGD（Online gradient descent）。

传统Batch算法优点是精度和收敛还可以，缺点是无法有效处理大数据集（此时全局梯度计算代价太大），且没法应用于数据流做在线学习

SGD存在的问题上面主要列了1）精度低；2）收敛慢；3）几乎得不到稀疏解。其中对online learning最重要的问题是SGD很难得到需要的正则化设计的解，特别是几乎得不到稀疏解

RDA，2010微软提出，特点：相对FOBOS，在精度与稀疏性之间做平衡，其中实验表明，在L1正则下，RDA比FOBOS可以更加有效地得到稀疏解。

 

与其他在线算法的迭代公式的对比（其实OGD如何一步步到类似形式的迭代公式的过程，限于时间，这里就不细说了，最后我会附一篇自己做分享会时做的ppt，里面有，感兴趣的可以下载看看），不同的方法在这种统一的描述形式下，区别点仅在第二项和第三项的处理方式：

不同的方法按照统一的描述形式，区别点主要在第二项和第三项：

• 第一项：梯度或累积梯度；

• 第二项：L1正则化项的处理。

 

• 第三项：这个累积加和限定了新的迭代结果x不要离已迭代过的解太远（也即FTRL-Proximal中proximal的含义），或者离0太远（central），这一项其实也是low regret的需求

这小结做了在线最优化算法的整理和总结，沿着稀疏性的主线，先后介绍了简单截断法、TG、FOBOS、RDA以及FTRL。从类型上来看，简单截断法、TG、FOBOS属于同一类，都是梯度下降类的算法,并且TG在特定条件可以转换成简单截断法和FOBOS；RDA属于简单对偶平均的扩展应用；FTRL可以视作RDA和FOBOS的结合，同时具备二者的优点。目前来看，RDA和FTRL是最好的稀疏模型Online Training的算法。

对于Online模式的并行化计算，一方面可以参考ParallelSGD的思路，另一方面也可以借鉴batch模式下对高维向量点乘以及梯度分量并行计算的思路。总之，在理解算法原理的基础上将计算步骤进行拆解，使得各节点能独自无关地完成计算最后汇总结果即可。

三.相关注意问题

1.关于批量梯度下降的batch_size选择问题

训练集较小【<2000】：直接使用batch梯度下降，每次用全部的样本进行梯度更新 
训练集较大：batch_size一般设定为[64,512]之间，设置为2的n次方更符合电脑内存设置，代码会运行快一些 
此外还要考虑GPU和CPU的存储空间和训练过程的波动问题 
 
batch_size越小，梯度的波动越大，正则化的效果也越强，自然训练速度也会变慢，实验时应该多选择几个batch_size进行实验，以挑选出最优的模型。

2.关于批量梯度下降的weight_decay【似乎与L2正则有关系，待补充】

3.关于优化算法选择的经验之谈

1. Adam在实际应用中效果良好，超过了其他的自适应技术。
2. 如果输入数据集比较稀疏，SGD、NAG和动量项等方法可能效果不好。因此对于稀疏数据集，应该使用某种自适应学习率的方法，且另一好处为不需要人为调整学习率，使用默认参数就可能获得最优值。
3. 如果想使训练深层网络模型快速收敛或所构建的神经网络较为复杂，则应该使用Adam或其他自适应学习速率的方法，因为这些方法的实际效果更优。
4. SGD通常训练时间更长，但是在好的初始化和学习率调度方案的情况下，结果更可靠。
5. Adadelta，RMSprop，Adam是比较相近的算法，在相似的情况下表现差不多。在想使用带动量的RMSprop，或者Adam的地方，大多可以使用Nadam取得更好的效果。

4.训练优化器的目的

加速收敛 2. 防止过拟合 3. 防止局部最优

5.选用优化器的目的

在构建神经网络模型时，选择出最佳的优化器，以便快速收敛并正确学习，同时调整内部参数，最大程度地最小化损失函数。

6.为什么神经网络的训练不采用二阶优化方法 (如Newton, Quasi Newton)？

7. 优化SGD的其手段

 
 

import tensorflow as tf
import os
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
from matplotlib.animation import FuncAnimation
%matplotlib tk
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'
plt.style.use('seaborn-white')

def func(x,y):
    return x**2-y**2

#制作画图用的网格数据
x1 = np.linspace(-1,1,100)
x2 = np.linspace(-1,2,100)
X1,X2 = np.meshgrid(x1,x2)
X = np.array(list(zip(X1.flatten(),X2.flatten())))
Y = np.array(list(map(lambda x:func(x[0],x[1]),X))).reshape(*X1.shape)
x = tf.Variable(initial_value=-0.5,dtype=tf.float32)
y = tf.Variable(initial_value=0.005,dtype=tf.float32)
z = tf.subtract(tf.square(x),tf.square(y))

def op_compare():
    sp = tf.Variable(1,dtype=tf.int64,trainable=False)
    sp_add = sp.assign(tf.add(sp,1))

    op_list = [
        tf.compat.v1.train.AdagradDAOptimizer(learning_rate=0.25,global_step=sp), #这个优化器一定要放在第一，否则会报错，也搞不懂为什么。
        tf.compat.v1.train.GradientDescentOptimizer(learning_rate=0.25),
        tf.compat.v1.train.MomentumOptimizer(learning_rate=0.25,momentum=0.9),
        tf.compat.v1.train.MomentumOptimizer(learning_rate=0.25,momentum=0.9,use_nesterov=True,name='NAG'),
        tf.compat.v1.train.AdadeltaOptimizer(0.25),
        tf.compat.v1.train.RMSPropOptimizer(0.25),
        tf.compat.v1.train.AdamOptimizer(0.25)

    ]
    
    data_list = []
    with tf.compat.v1.Session() as sess:
        for i,init_op in enumerate(op_list):
            print('第{}个优化器'.format(i+1))
            op = init_op.minimize(z)
            sess.run(tf.compat.v1.global_variables_initializer())
            current_z = sess.run(z)
            diff = 2
            stop_e = 0.0001
            x_list = []
            print(sess.run(sp))
            res_x, res_y, res_z = sess.run([x,y,z])
            print(res_x,res_y)
            x_list.append([res_x,res_y,res_z])
            print(op.name)
            while diff>stop_e:
                sess.run([op,sp_add])
                res_x, res_y, res_z = sess.run([x,y,z])
                if res_z==np.inf or res_z==-np.inf or res_z<np.min(Y)-5:
                    break
                x_list.append([res_x,res_y,res_z])
                diff = abs(res_z-current_z)
                current_z = res_z
            data = np.array(x_list)
            data_list.append((op.name,data))
        return data_list

data_list = op_compare()

fig = plt.figure(figsize=(18,16))
ax = Axes3D(fig)
ax.plot_surface(X1,X2,Y,cmap=plt.cm.plasma)
count = 0
colors = ['r','g','b','y','purple','brown','pink']

line_list = [
    ax.plot([], [], [], color='{}'.format(c),marker='o', ls='-', markevery=[-1],markersize=5, animated=False,label='{}'.format(dt[0])) 
    for dt,c in zip(data_list,colors)
]   #重新设置曲线的值

def init():
    ax.set_ylim([-1, 2])
    ax.set_xlim([-1, 1])
    ax.set_zlim([np.min(Y),np.max(Y)])
    return line_list               #返回曲线

def update(n,data,lines):
#     print(n)
    for dt,ln in zip(data,lines):
        show_data = dt[1]
        id = n+1
        if n >= len(show_data)-1:
            id = -1
        xdata = show_data[:id,0]
        ydata = show_data[:id,1]
        zdata = show_data[:id,2]
        ln[0].set_data([xdata,ydata])
        ln[0].set_3d_properties(zdata)
    
    return lines
plt.legend()

#获取数据的最大长度
record_len = []
for d in data_list:
    record_len.append(len(d[1]))
max_len = max(record_len)
ani = FuncAnimation(fig, update, frames=max_len, fargs=(data_list,line_list),
                    init_func=init,interval=100,)#这里的frames在调用update函数是会将frames作为实参传递给“n”

plt.show()

附录