几种常见梯度优化方法

优化算法是机器学习领域的重要内容，本文介绍几种常见的无约束的优化算法。关于无约束问题优化方法的一般讨论请参考此文。

• 梯度下降法 （Gradient Descent）
• 动量法 （Momentum）
• 共轭梯度法（Conjugate Gradient）
• 自然梯度法（Natural Gradient）

梯度下降法 （Gradient Descent）

凡是接触过机器学习、强化学习、凸优化等领域的读者都会对梯度下降法非常熟悉。鉴于此原因，这里也不做详细介绍，仅给出基本方法。关于梯度下降在凸优化中的应用请参考此文。

学过数学的人都知道，函数的梯度方向表示了函数增长速度最快的方向，那么和它相反的方向就可以看作是函数值减少速度最快的方向。就机器学习模型优化问题而言，当目标被设定为求解目标函数的最小值的时候，只要朝着梯度下降的方向前进，就能不断接近最优值。

这里给出一种最简单的方法——固定步长方法：

def gd(x_start, step, g):
    x = x_start\
    for i in range(20):
        grad = g(x)
        x -= grad * step
        print '[ Epoch {0} ] grad = {1}, x = {2}'.format(i,grad,x)
        if abs(grad) < 1e-6:
            break
    return x

这里 x 保存当前优化过程中的参数值，g 是待优化函数 f(x) 的导数，step 是步长。步长的选择十分重要，步长太小导致收敛速度太慢；而步长太大则导致不收敛问题。关于调整步长目前已经有成熟的方法，这里不作细致讨论。

动量法 （Momentum）

动量代表了优化过程中累计的能量，它将在后面的优化过程中持续作用，推动目标值前进。动量会以一定的形式衰减。基于上一节的梯度下降方法，这里给出动量法的简单实现：

def momentum(x_start, step, g,discount=0.7):
    x = np.array(x_start,dtype = 'floast64')
    pre_grad = np.zeros_like(x)
    for i in range(50):
        grad = g(x)
        pre_grad = pre_grad * discount + grad * step
        x -= pre_grad
        print '[ Epoch {0} ] grad = {1}, x = {2}'.format(i,grad,x)
        if abs(grad) < 1e-6:
            break
    return x

可以看出，代码中多了 pre_grad。这个变量用于存储历史动量的积累，每一轮都会乘以一个打折量 discount 。动量法可以帮助目标值穿越“狭窄山谷”形状的优化曲面，从而达到最终的最优点。很多时候，目标函数在不同维度上的“陡峭”程度差别很大，固定步长的方法可在每一个维度上以相同的步长进行迭代，这就导致某些维度上更新乏力，而某些维度上大量进行无用功。

使用了动量后，历史的更新会以衰减的形式不断作用在这些方向上，反复做的无用功可以相互抵消，而更新不足的方向则在加强。然而，动量优化也存在一些问题，前几轮迭代累计不足，目标值在反复做无用功的方向震荡更大，带来更强烈的抖动。又有科研人员发明了解决方法——Nesterov算法。我们对上述代码稍加改造：

import numpy as np
def nesterov(x_start, step, g, discount = 0.7):
    x = np.array(x_start,dtype = 'floast64')
    pre_grad = np.zeros_like(x)
    for i in range(50):
        x_future = x - step * discount * pre_grad
        grad = g(x_future)
        pre_grad = pre_grad * discount + grad
        x -= pre_grad * step
        print '[ Epoch {0} ] grad = {1}, x = {2}'.format(i,grad,x)
        if abs(grad) < 1e-6:
            break
    return x

可以看出，这里用 x_future 存储下一步的变量，然后计算更新后的优化点的梯度。这是与动力法的主要区别。想象一个场景，优化点累积到了某个抖动方向后的梯度后，对动量算法来说，虽然梯度指向累计梯度相反的方向，但是刚开始动量不够大，所以优化点还是会在累计的方向上前进一小段。对 Nesterov 方法来说，如果在该方向前进一小段，则梯度中指向累积梯度相反的方向的量会大很多，最终两个方向梯度抵消，反而使抖动方向的量迅速减少。

共轭梯度法（Conjugate Gradient）

共轭梯度法是用来求解正定二次规划的一种优化方法，它具有二次终止性。前面的讨论中，我们并没有给出优化问题的具体形式，本节我们首先将待优化问题形式定义为：

minx12XTAX−bTX$$\underset{x}{min}\frac{1}{2}{X}^{T}AX-b^{T}X$$

其中 A$A$ 为半正定矩阵，b$b$ 为已知量，对公式进行求导并令导数为 0，得：

AX∗=b$$A{X}^{\ast }=b$$

其中 X∗$X^{\ast }$ 为最优值。我们虽然可以直接求解上述方程，但是对于高纬度问题计算量较大，所以我们希望通过优化的方式找到最优解。

约束与共轭

梯度下降法的每一步都朝着局部最优的方向前进，但是几轮迭代它前进的方向可能非常相似，因为这个方向没有更新完，未来还会存在这个方向的残差。于是，我们对优化提出了更高要求——能不能再选中优化方向和步长上更新智能，每朝一个方向走，就把这个方向走到极致（极致就是再往前走就会远离目标）。我们引入一个误差变量 et$e_t$

et=x∗−xt$$e_t=x^{\ast }-x_t$$

误差反映了参数的最优点和当前点之间的距离（向量），那么我们的目标可以定义为：每一步优化后，当前的误差和刚才优化的方向正交。优化的方向不一定是梯度方向。

我们令 rt$r_t$ 表示第 t$t$ 轮的更新方向，可以得到：

rTtet+1=0$$r_t^T{e}_{t+1}=0$$

有了这个约束，如果我们的优化空间有 d$d$ 维，理论上最多只需要迭代 d$d$ 轮就可以求解出来。然而公式中的误差项我们并不知道。未解决此问题，我们引入共轭这个概念。在上述正交约束中加入一个矩阵 A$A$ 使得：

rTtAet+1=0$$r_t^{T}A{e}_{t+1}=0$$

现状两者关系变成了共轭正交。之前所说的正交可以理解为 A=I$A=I$ 时的共轭正交。共轭正交向量还具有其它性质。我们可以将共轭正交的关系看作是 rt$r_t$ 与经过线性变换后的 et+1$e_{t+1}$ 构成的正交关系。这个新的正交实际上只是多绕了一个弯，与原来的正交差距不大，有个这个设定，下面的内容就比较好理解了。

最优步长

共轭梯度法属于线搜索的一种，类似于梯度下降，优化过程分两步：

1. 确定优化方向；
2. 确定优化步长。

这里我们先来介绍如何确定优化步长，优化方向的确定留到后面。假设当前参数为 Xt$X_t$：

rTtAet+1=rTtA[et+Xt−Xt+1]=rTtA[et+αtrt]=rTtAet+αtrTtArt=0$$\begin{array}{rl}{r}_t^{T}A{e}_{t+1}& =r_t^{T}A[e_t+X_t-X_{t+1}]\\ & =r_t^{T}A[e_t+{\alpha }_t{r}_t]\\ & =r_t^{T}A{e}_t+{\alpha }_t{r}_t^{T}A{r}_t\\ & =0\end{array}$$

αt=−rTtAetrTtArt=−rTtA(X∗−Xt)rTtArt$${\alpha }_t=-\frac{r_t^{T}A{e}_t}{r_t^{T}A{r}_t}=-\frac{r_t^{T}A(X^{\ast }-X_t)}{r_t^{T}A{r}_t}$$

回想我们对优化问题的定义可知 AX∗=b$A{X}^{\ast }=b$，第 t$t$ 轮的梯度 gt=AXt−b$g_t=A{X}_t-b$ ，于是：

αt=rTtgtrTtArt$${\alpha }_t=\frac{r_t^T{g}_t}{r_t^{T}A{r}_t}$$

到这里，公式中的未知量 et$e_t$ 已被抵消，公式变得可解。

Gram-Schmidt方法

下面我们讨论如何确定优化方向。我们要解决的主要问题是如何让优化方向和误差正交。由于每一次优化后，剩下的误差和本次的优化方向正交（共轭正交），所以可以看出每一个优化方向都是彼此正交的。现在问题变成了如何构造这些正交的方向。

Gram-Schmidt 方法是线性代数中常见的一种向量正交化方法。这个算法的大意是在 N$N$ 维空间中输入 N$N$ 个线性无关的向量，输出是 N$N$ 个相互正交的向量，也就是我们最终想要的向量组合。假设输入的 N$N$ 个向量为 u1,u2,…,uN$u_1,u_2,\dots ,u_N$ 输出的 N$N$ 个向量为 d1,d2,…,dN$d_1,d_2,\dots ,d_N$， 它的具体做法如下：

1. 对于第一个向量，保持不变：d1=u1$d_1=u_1$；
2. 对于第二个向量，去掉其与第一个向量共线的部分：d2=u2+β1d1$d_2=u_2+{\beta }_1{d}_1$；
3. 对于第三个向量，去掉其与第一、二个向量共线的部分：d3=u3+∑2i=1β3,idi$d_3=u_3+\sum _{i=1}^2{\beta }_{3,i}{d}_i$；
4. 对于第N$N$个向量，去掉其与前 N−1$N-1$ 个向量共线的部分：dN=uN−∑N−1i=1βN,idi$d_N=u_N-\sum _{i=1}^{N-1}{\beta }_{N,i}{d}_i$

如何确定这些 β$\beta$ 呢？利用共轭正交的性质有：

dTlAdt=0,(l=1,2,…,t−1)$$d_l^{T}A{d}_t=0,(l=1,2,\dots ,t-1)$$

进一步展开：

dTlAdt=dTlA(ut+∑i=1t−1βt,idi)=dTlAut+dTlA∑i=1t−1βt,idi$$d_l^{T}A{d}_t=d_l^{T}A(u_t+\sum _{i=1}^{t-1}{\beta }_{t,i}{d}_i)=d_l^{T}A{u}_t+d_l^{T}A\sum _{i=1}^{t-1}{\beta }_{t,i}{d}_i$$

利用正交的性质（任意两个互异输出均相互共轭正交），公式可进一步简化为：

dTlAut+dTlAβt,ldl=0$$d_l^{T}A{u}_t+d_l^{T}A{\beta }_{t,l}{d}_l=0$$

所以我们最终得到：

βt,l=−dTlAutdTlAdl$${\beta }_{t,l}=-\frac{d_l^{T}A{u}_t}{d_l^{T}A{d}_l}$$

共轭梯度

到这里，还有两个问题待解决：

1. 要用什么向量构造这些正交向量？
2. 随着向量数量的增加，我们需要计算的参数越来越多。这个计算量能不能减少呢？

解决上述两个问题的方法就是——用每一步的梯度构建这些正交向量，同时利用梯度的特点减少计算量。这便是共轭梯度法名称的由来。

最后我们来解决第二个问题。要解决此问题，我需要证明三个推论，这里我们依次证明。

推论一：第 t$t$ 步计算的梯度 gt$g_t$ 和前 t−1$t-1$ 步的更新 {dj}t−1j=1$\{d_j{\}}_{j=1}^{t-1}$ 正交。
证明：最优点 x∗$x^{\ast }$ 到初始参数值 x1$x_1$ 的距离为 e1$e_1$。

e1=∑i=1Tαidi$$e_1=\sum _{i=1}^T{\alpha }_i{d}_i$$

这样我们就用更新方向表示了误差，那么：当 i<j$i<j$ 时

dTigj=dTi(AXj−b)=dTi(AXj−AX∗)=dTiAej=dTiA∑i=1Tαidi$$d_i^T{g}_j=d_i^T(A{X}_j-b)=d_i^T(A{X}_j-A{X}^{\ast })=d_i^{T}A{e}_j=d_i^{T}A\sum _{i=1}^T{\alpha }_i{d}_i$$

由于利用 Gram-Schmidt 方法构造的更新方向均是互相共轭正交的，所以

dTigj=0,i<j$$d_i^T{g}_j=0,i<j$$

证毕。

推论二：第 t$t$ 步计算的梯度 gt$g_t$ 和前 t−1$t-1$ 步的梯度{gj}t−1j=1$\{g_j{\}}_{j=1}^{t-1}$ 正交。
证明：我们可以直接利用推论一证明推论二。注意这里的 gi$g_i$ 相当于Gram-Schmidt方法中提到的输入 ui$u_i$。

gTigj=(di−∑k=1i−1βkdk)Tgj=0$$g_i^T{g}_j=(d_i-\sum _{k=1}^{i-1}{\beta }_k{d}_k{)}^T{g}_j=0$$

证毕。
此时我们发现，每一步的梯度与此前的梯度正交，因为我们能用梯度作为线性无关的向量组，从而组成共轭正交的向量组。

推论三：第 t$t$ 步计算的梯度 gt$g_t$ 和前 t−2$t-2$ 步的更新 {dj}t−2j=1$\{d_j{\}}_{j=1}^{t-2}$ 正交。
证明：

Xi+1=Xi+αidi$$X_{i+1}=X_i+{\alpha }_i{d}_i$$

dTiAgjgTjAdi=(gTjAdi)T( A 是对称矩阵 )=gTjA1αi(Xi+1−Xi)=1αigTj(AXi+1−AXi)=1αigTj(gi+1−gi)=0,i+1<j (推论二) $$\begin{array}{rl}{d}_i^{T}A{g}_j& =(g_j^{T}A{d}_i{)}^T\text{( }A\text{ 是对称矩阵 )}\\ g_j^{T}A{d}_i& =g_j^{T}A\frac{1}{{\alpha }_i}(X_{i+1}-X_i)\\ & =\frac{1}{{\alpha }_i}{g}_j^T(A{X}_{i+1}-A{X}_i)\\ & =\frac{1}{{\alpha }_i}{g}_j^T(g_{i+1}-g_i)\\ & =0,i+1<j\text{ (推论二) }\end{array}$$

证毕。

所以可以看出，使用 Gram-Schmidt 方法可以保证更新方向共轭正交，对于第 t$t$ 轮优化，我们只计算 βt−1${\beta }_{t-1}$ 即可，其它 β$\beta$ 均是0。

下面来简单看下这个算法的实现：

import numpy as np
def cg(f_Ax,b,cg_iters=10,callback = None,verbose=False,residual_tol=1e-10):
    p = b.copy()
    r = b.copy()
    x = np.zeros_like(b)
    rdotr = r.dot(r)
    for i in range(cg_iters):
        z = f_Ax(p)
        v = rdotr / p.dot(z)
        x += v*p
        r -= v*z
        newrdotr = r.dot(r)
        mu = newrdotr/rdotr
        p = r + mu*p
        rdotr = newrdotr
        if rdotr < residual_tol:
            break
    return x

代码中用到了残差 r$r$ ，真正的优化方向是 p$p$，优化步长是 v$v$。

自然梯度法（Natural Gradient）

当优化问题的两个坐标轴的尺度差异比较大时，使用一个统一的学习率会产生问题：某一个坐标轴可能会发散。尺度不同在优化问题中很难避免，虽然每一轮迭代对参数的更新量相差不多，但他们对模型的影响完全不同。梯度下降法针对的目标是参数，我们更新的目标也是参数，而优化的真正目标是找到最优的模型。每一轮迭代中，对参数的改进和对模型的改进是不同的。梯度下降法只考虑了在梯度方向上对参数进行更新，并没有考虑模型层面的更新，因此就会出现模型更新不均匀的现象。

自然梯度法不简单地使用学习率对参数更新进行量化，而是对模型效果进行量化。假设我们的待优化模型为 f(w)$f(w)$，其中参数为 w$w$。按照梯度下降法的思想，我们定义每一轮的迭代都要解决这样一个子问题：

minΔwf(w+Δw)s.t. |Δw|<ϵ$$\begin{array}{rl}& \underset{\mathrm{\Delta }w}{min}f(w+\mathrm{\Delta }w)\\ & \text{s.t. }|\mathrm{\Delta }w|<ϵ\end{array}$$

将原函数进行一阶Taylor展开，问题变为：

minΔwf(w)+∇wf(w)Δws.t. |Δw|<ϵ$$\begin{array}{rl}& \underset{\mathrm{\Delta }w}{min}f(w)+{\mathrm{\nabla }}_{w}f(w)\mathrm{\Delta }w\\ & \text{s.t. }|\mathrm{\Delta }w|<ϵ\end{array}$$

如果我们改为对模型的距离进行约束，就可以得到：

minΔwf(w)+∇wf(w)Δws.t. KL(f(w),f(w+Δw))<ϵ$$\begin{array}{rl}& \underset{\mathrm{\Delta }w}{min}f(w)+{\mathrm{\nabla }}_{w}f(w)\mathrm{\Delta }w\\ & \text{s.t. }\mathrm{K}\mathrm{L}(f(w),f(w+\mathrm{\Delta }w))<ϵ\end{array}$$

这里模型的变化用KL散度进行衡量（注意KL散度没有对称性）：

KL(f(x,w1)∥f(x,w2))=∫xf(x,w1)logf(x,w1)f(x,w2)dx$$\mathrm{K}\mathrm{L}(f(x,w_1)\Vert f(x,w_2))={\int }_{x}f(x,w_1)\log \frac{f(x,w_1)}{f(x,w_2)}\mathrm{d}x$$

它可以比较两个概率分布之间的差异。这边是自然梯度法的优化形式。可以看出，有了模型层面的约束，每一轮迭代无论参数发生多大变化，模型的变化都会限制在一定的范围内。然而KL散度的计算并不容易。这个问题如何解决呢？

Fisher信息矩阵

KL散度可以变成熵和交叉熵的运算：

KL(f(w)∥f(w+Δw))=Efw[logf(w)f(w+Δw)]=Efw[logf(w)]−Efw[logf(w+Δw)]$$KL(f(w)\Vert f(w+\mathrm{\Delta }w))=E_{f_w}\left[\log \frac{f(w)}{f(w+\mathrm{\Delta }w)}\right]=E_{f_w}[\log f(w)]-E_{f_w}[\log f(w+\mathrm{\Delta }w)]$$

通过对上式右边第二项进行Taylor展开，我们最终可以得到：

KL(f(w)∥f(w+Δw))≈Efw[logf(w)]−Efw[logf(w)+∇wlogf(w)Δw+12ΔwT∇2wlogf(w)Δw)]=−Efw[∇wlogf(w)Δw]−Efw[12ΔwT∇2wlogf(w)Δw)]=−[∫xf(w)1f(w)∇wf(w)dx]Δw−12ΔwT[∫xf(w)∇2wlogf(w)dx]Δw$$\begin{array}{rl}KL(f(w)\Vert f(w+\mathrm{\Delta }w))& \approx E_{f_w}[\log f(w)]-E_{f_w}[\log f(w)+{\mathrm{\nabla }}_w\log f(w)\mathrm{\Delta }w+\frac{1}{2}\mathrm{\Delta }{w}^T{\mathrm{\nabla }}_w^2\log f(w)\mathrm{\Delta }w)]\\ & =-E_{f_w}[{\mathrm{\nabla }}_w\log f(w)\mathrm{\Delta }w]-E_{f_w}[\frac{1}{2}\mathrm{\Delta }{w}^T{\mathrm{\nabla }}_w^2\log f(w)\mathrm{\Delta }w)]\\ & =-[{\int }_{x}f(w)\frac{1}{f(w)}{\mathrm{\nabla }}_{w}f(w)\mathrm{d}x]\mathrm{\Delta }w-\frac{1}{2}\mathrm{\Delta }{w}^T[{\int }_{x}f(w){\mathrm{\nabla }}_w^2\log f(w)\mathrm{d}x]\mathrm{\Delta }w\end{array}$$

我们对 f(w)$f(w)$ 的定义一般都是一个连续、可微、有界、性质优良的函数，所以第一项积分和微分可以互换，上式最终变为：

KL(f(w)∥f(w+Δw))=−12ΔwTEfw[∇2wlogf(w)]Δw$$KL(f(w)\Vert f(w+\mathrm{\Delta }w))=-\frac{1}{2}\mathrm{\Delta }{w}^T{E}_{f_w}[{\mathrm{\nabla }}_w^2\log f(w)]\mathrm{\Delta }w$$

这里包含一个二阶导的期望值，仍然比较复杂。

为了进一步简化，我们引入Fisher信息矩阵。我们首先定义 Score函数 为对数似然函数的一阶导数：

lfw=∇wf(x,w)$$l_{f_w}={\mathrm{\nabla }}_{w}f(x,w)$$

那么可以看出，Score函数的期望值是0，即：

Efw[lfw]=∫xf(w)∇wlogf(w)dx=0$$E_{f_w}[l_{f_w}]={\int }_{x}f(w){\mathrm{\nabla }}_w\log f(w)\mathrm{d}x=0$$

Fisher信息矩阵可以通过Score函数定义：

Ifw=Efw[lfwlTfw]$$I_{f_w}=E_{f_w}[l_{f_w}{l}_{f_w}^T]$$

可以证明：在概率分布函数具备良好性质时，Fisher信息矩阵和KL散度的二阶导数的相反数相等。证明过程非常直接，此处从略。

自然梯度法的目标公式

通过前面的推演，我们的目标函数变为：

minΔwf(w)+∇wf(w)Δws.t. 12ΔwTIfwΔw<ϵ$$\begin{array}{rl}& \underset{\mathrm{\Delta }w}{min}f(w)+{\mathrm{\nabla }}_{w}f(w)\mathrm{\Delta }w\\ & \text{s.t. }\frac{1}{2}\mathrm{\Delta }{w}^T{I}_{f_w}\mathrm{\Delta }w<ϵ\end{array}$$

该问题可以通过拉格朗日乘子法表示为：

minΔwf(w)+∇wf(w)Δw+λ[12ΔwTIfwΔw−ϵ]$$\underset{\mathrm{\Delta }w}{min}f(w)+{\mathrm{\nabla }}_{w}f(w)\mathrm{\Delta }w+\lambda [\frac{1}{2}\mathrm{\Delta }{w}^T{I}_{f_w}\mathrm{\Delta }w-ϵ]$$

求导取极限点，可以得到：

∇wf(w)+λIfwΔw=0Δw=−1λI−1fw∇wf(w)$$\begin{array}{rl}& {\mathrm{\nabla }}_{w}f(w)+\lambda I_{f_w}\mathrm{\Delta }w=0\\ & \mathrm{\Delta }w=-\frac{1}{\lambda }{I}_{f_w}^{-1}{\mathrm{\nabla }}_{w}f(w)\end{array}$$

公式中 1λ$\frac{1}{\lambda }$ 可以作为梯度下降法的学习率类似的分量，那么自然梯度法的优化方向就可以看作 I−1fw∇wf(w)$I_{f_w}^{-1}{\mathrm{\nabla }}_{w}f(w)$，与梯度下降法不同，它需要额外求解 Fisher 信息矩阵的逆。

感谢冯超——《强化学习精要：核心算法与TensorFlow实现》（电子工业出版社）