强化学习的学习之路（四十五）2021-02-14共轭梯度法实现策略上的单调提升(Monotonic Improvement with Dual gradient descent)

作为一个新手，写这个强化学习-基础知识专栏是想和大家分享一下自己学习强化学习的学习历程，希望对大家能有所帮助。这个系列后面会不断更新，希望自己在2021年能保证平均每日一更的更新速度，主要是介绍强化学习的基础知识，后面也会更新强化学习的论文阅读专栏。本来是想每一篇多更新一点内容的，后面发现大家上CSDN主要是来提问的，就把很多拆分开来了（而且这样每天任务量也小一点哈哈哈哈偷懒大法）。但是我还是希望知识点能成系统，所以我在目录里面都好按章节系统地写的，而且在github上写成了书籍的形式，如果大家觉得有帮助，希望从头看的话欢迎关注我的github啊，谢谢大家！另外我还会分享深度学习-基础知识专栏以及深度学习-论文阅读专栏，很早以前就和小伙伴们花了很多精力写的，如果有对深度学习感兴趣的小伙伴也欢迎大家关注啊。大家一起互相学习啊！可能会有很多错漏，希望大家批评指正！不要高估一年的努力，也不要低估十年的积累，与君共勉！

共轭梯度法实现策略上的单调提升(Monotonic Improvement with Dual gradient descent)

我们先介绍一下共轭梯度法，再介绍一下共轭梯度法在策略的单调提升上的应用。

首先，开门见山地介绍一些共轭梯度法：共轭梯度法的目的：让每一次优化变得“极致”。

前面我们一般都用策略梯度法去解决优化问题。虽然梯度下降法的每一步都朝着局部最优的方向前进，但它在不同的迭代轮数会选择非常相近的方向，这说明当某一次选择了一个更新方向和步长后，这个方向并没有被更新完，未来还会存在这个方向的残差。如果把参数更新的轨迹显示出来，我们可以看到有时轨迹会走成 下面这种zig-zag 的形状：

于是，我们对优化提出了更高的要求：每朝一个方向走，就把这个方向优化到极致，使得之后的更新再也不需要朝这个方向进行。我们引入一个变量一一误差，假设最优的参数为$x^{\ast },$ 当前第 $t$ 轮的参数为 $x_t,$ 误差可以定义为
$$e_t=x^{\ast }-x_t$$
这个误差表示了参数的最优点和当前点之间的距离。那么目标就更明确了，我们希望每一步优化后，当前的误差和刚才优化的方向正交。现在令 $r_t$ 表示第 $t$ 轮的更新方向，就可以得到下面的公式：

​ $r_t^{\mathrm{T}}{e}_{t+1}=0$
假设每一轮的优化量都和误差正交，那么如果我们的优化空间有 $d$ 维，理论上最多只需要迭代 $d$ 轮就可以求解出来，这样在优化时间上就有了保证。如果我们直接使用这个公式，就会发现一个问题: 公式中需要知道误差。换句话说，我们需要知道最优点在哪儿。如果我们知道最优点在哪儿，就不用优化了。 可是不知道最优点儿，这个方法又无法直接使用，仿佛陷入了一个死循环，那么接下来我们就利用数学工具来试着推导解决这个问题。

前面我们提到优化方向和误差正交，如果使用了共轭这个工具，现在两者的关系将变为共轭正交，也就是存在一个矩阵 $A(A$就是轭 )，使得优化方向和误差满足正交的性质:
$$r_t^{\mathrm{T}}A{e}_{t+1}=0$$
等等，这不就说明轭的存在使原本正交的两项变得不正交了吗? 但其实如果我们在上面的原始公式中间加一个特殊的矩阵，例如单位阵：
$$r_t^{\mathrm{T}}\mathbf{I}{\mathbf{e}}_{t+1}=0$$
加入单位阵并不会改变结果，所以原本正交的二者依然正交。所以换个角度理解，我们以前简单的正交都是在单位阵这个简单且性质优良的“轭”的下实现的。所以回到正交中，如果使用单位阵作为轭, 那么绑定在一起的就是常规意义上正交的一对向量; 如果使用其他矩阵，那么共轭正交的向量肯定也会满足其他的性质。

明确了共轭梯度法的目标和特点，我们就要开始推导算法公式了。共轭梯度法属于线搜索的一种，因此和梯度下降法类似，我们的总体思路不变，优化过程分如下两步:

• 确定优化方向

• 确定优化步长。

我们先来介绍优化步长的计算方法。假设当前的参数为 $X_t,$ 我们已经得到了优化方向 $r_t,$ 下面要确定的就是步长 ${\alpha }_t,$ 根据前面提过的共轭正交公式
$${\mathbf{r}}_t^{\mathrm{T}}\mathbf{A}{\mathbf{e}}_{t+1}=0$$
我们可以开始推导:
$$\begin{array}{rl}{r}_t^{\mathrm{T}}\mathbf{A}{\mathbf{e}}_{t+1}& ={\mathbf{r}}_t^{\mathrm{T}}\mathbf{A}\left[{\mathbf{e}}_t+{\mathbf{X}}_t-{\mathbf{X}}_{t+1}\right]\\ & ={\mathbf{r}}_t^{\mathrm{T}}\mathbf{A}\left[{\mathbf{e}}_t+{\alpha }_t{\mathbf{r}}_t\right]\\ & ={\mathbf{r}}_t^{\mathrm{T}}\mathbf{A}{\mathbf{e}}_t+{\alpha }_t{\mathbf{r}}_t^{\mathrm{T}}\mathbf{A}{\mathbf{r}}_t=0\end{array}$$
于是可以整理得到
$$\begin{array}{l}{\alpha }_t=-\frac{r_t^{\mathrm{T}}A{e}_t}{r_t^{\mathrm{T}}A{r}_t}\\ {\alpha }_t=-\frac{r_t^{\mathrm{T}}A\left({\mathbf{X}}^{\ast }-{\mathbf{X}}_t\right)}{{\mathbf{r}}_t^{\mathrm{T}}\mathbf{A}{r}_t}\end{array}$$
我们知道 $A{X}^{\ast }=b,$ 第 $t$ 轮的梯度 $g_t=A{X}_t-b,$ 于是公式最终变为
$${\alpha }_t=\frac{{\mathbf{r}}_t^{\mathrm{T}}{\mathbf{g}}_t}{{\mathbf{r}}_t^{\mathrm{T}}\mathbf{A}{r}_t}$$
到这里，我们利用矩阵 $A$ 成功地把公式中的误差 $e$ 抵消，于是步长变得可解。

完成了步长的求解，接下来就要回到第一步看看优化方向的计算方法。我们要解决的主要问题是如何让优化方向和误差正交。由于每一次的优化后，剩下的误差和本次的优化正交（共轭正交 ) , 所以可以看出每一个优化方向彼此间都是正交的。那么，我们接下来就用Gram-Schmidt 方法来构建这些彼此正交的向量。

在线性代数课上, 我们曾经学过一个向量正交化的方法一 Gram-Schmidt 方法。这个算法的输入是 $N$ 维空间中 $N$ 个线性无关的向量，由于向量间线性无关，任何一个向量都无法通过其他向量表达出来。算法的输出是 $N$ 个相互正交的向量，也就是我们最终想要的向量组合。它的具体算法如下：

令输入向量为 $u_1,u_2,\cdots ,u_N,$ 输出向量为 $d_1,d_2,\cdots ,d_N,$ 那 么有:

• 对于第一个向量，我们保持它不变: $u_{1}=d_{1} $

• 对于第二个向量，我们去掉其中和第一个向量共线的部分，令去掉的比例为 ${\beta }_i,$ 所以第二个向量 $d_2$ 等于 $u_2+{\beta }_1{d}_1$

• 对于第三个向量，我们去掉其中和第一、第二个向量共线的部分: $d_3={\mathbf{u}}_3+$ ${\sum }_{i=1}^2{\beta }_{3,i}{d}_i$

• 对于第 $N$ 个向量，我们去掉其中和第一、第二、第 $N-1$ 个向量共线的部分 :${\mathbf{d}}_N={\mathbf{u}}_N+{\sum }_{i=1}^{N-1}{\beta }_{N,i}{\mathbf{d}}_i$

  那么我们怎么求这些比例项 $\beta$ 呢? 我们利用前面提到的性质，向量之间正交（这里还是共轭正交 )，于是有

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

​ 进一步展开，可以得到：
$${\mathbf{d}}_l^{\mathrm{T}}\mathbf{A}{\mathbf{d}}_t={\mathbf{d}}_l^{\mathrm{T}}\mathbf{A}\left({\mathbf{u}}_t+\sum _{i=1}^{t-1}{\beta }_{t,i}{\mathbf{d}}_i\right)={\mathbf{d}}_l^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_t+{\mathbf{d}}_l^{\mathrm{T}}\mathbf{A}\sum _{i=1}^{t-1}{\beta }_{t,i}{\mathbf{d}}_i$$

​ 利用正交的性质，可以将公式化简为
$${\mathbf{d}}_l^{\mathrm{T}}\mathbf{A}{u}_t+{\mathbf{d}}_l^{\mathrm{T}}\mathbf{A}{\beta }_{t,l}{\mathbf{d}}_l=0$$
​ 所以最终我们得到
$${\beta }_{t,l}=-\frac{d_l^{\mathrm{T}}A{u}_l}{d_l^{\mathrm{T}}\mathbf{A}{d}_l}$$
求得了这个变量, 就可以依此求出当前向量与其他向量的共线部分, 从而确保所有向量相互正交的性质。当然，我们也看到，为了实现正交的效果，算法的复杂度为 $O\left(N^2\right),$ 在实现速度上不是很快。

接下来我们需要解决如下两个问题：

• 我们要用什么向量构建这些正交向量?
• 随着向量数量的增加，我们需要计算的参数越来越多，我们一共要计算 $O\left(N^2\right)$ 数量级的 $\beta ,$ 这个计算的数目还是有点多，能不能再减少呢?

要解决上面两个问题我们需要先证明下面三个推论：

1. 第 $t$ 步计算的梯度 $g_t$ 和前 $t-1$ 步的更新 { d j } j = 1 t − 1 \left\{d_{j}\right\}_{j=1}^{t-1} {dj​}j=1t−1​ 正交
2. 第 $t$ 步计算的梯度 $g_t$ 和前 $t-1$ 步的梯度 { g t } j = 1 t − 1 \left\{\boldsymbol{g}_{t}\right\}_{j=1}^{t-1} {gt​}j=1t−1​ 正交
3. 第 $t$ 步计算的梯度 $g_t$ 和前 $t-2$ 步的更新 { d j } j = 1 t − 2 \left\{\boldsymbol{d}_{j}\right\}_{j=1}^{t-2} {dj​}j=1t−2​ 共轭正交。

先来证明第一个推论：前面提到共轭梯度法的基本思想：每一轮把某一方向优化彻底，保证后面的优化不再对这个方向做任何操作。假设算法一共进行了 $T$ 轮迭代，我们就可以用这 $T$ 轮求出的优化方向组合成最终的误差。假设我们的初始参数值为 $x_1,$ 那么它到最优点 $x^{\ast }$ 的距离为 $e_1,$ 这也是当前的误差。第 $t$ 步求出的优化方向为 $d_t,$ 每一轮的优化步长为 ${\gamma }_i,$ 根据上面的定义，可以得到如下公式：
$$e_1=\sum _{i=1}^T{\gamma }_i{d}_i$$
这样我们就用更新方向表示了误差。对于不同迭代轮数的误差，我们也可以用类似的公式表示。上面的公式利用所有的更新组合成误差，我们反过来思考也可以理解，通过每一步的更新，算法最终收敛，误差变为 $0。$ 所以，这里的 $\gamma$ 和前面提到的步长 $\alpha$之间是等价的关系。
下面就来证明 :

$$\text{ 当 }i<j\text{ 时， }$$
$$\begin{array}{l}{\mathbf{d}}_i^{\mathrm{T}}{\mathbf{g}}_j=0\end{array}$$
证明过程如下：

​ $${\mathbf{d}}_i^{\mathrm{T}}{\mathbf{g}}_j={\mathbf{d}}_i^{\mathrm{T}}\left(\mathbf{A}{\mathbf{X}}_j\right)$$

​ $$\begin{array}{l}=d_i^{\mathrm{T}}\left(\mathbf{A}{\mathbf{X}}_j-\mathbf{A}{\mathbf{X}}^{\ast }\right)\\ =d_i^{\mathrm{T}}\mathbf{A}{\mathbf{e}}_j\\ =d_i^{\mathrm{T}}\mathbf{A}\left({\sum }_{t=j}^T{\gamma }_t{\mathbf{d}}_t\right)\end{array}$$

因为前面已经说明任意两个优化方向是相互正交的 ( Gram-Schmidt 方法 ), 所以这个式子的结果为 0 。利用共轭正交的性质，我们得到了更好的正交性质，这就是共轭梯度法最精髓的地方下面证明第二个推论。我们直接给出推导过程：
$${\mathbf{g}}_i^{\mathrm{T}}{\mathbf{g}}_j={\left({\mathbf{d}}_i+\sum _{j=1}^{i-1}{\beta }_j{\mathbf{d}}_j\right)}^{\mathrm{T}}{\mathbf{g}}_j=0$$
此时我们发现，由于每一步的梯度与此前的梯度正交，因此我们可以使用它作为线性无关的向量组，从而组成共轭正交的向量组。第三个推论将以第二个推论为基础，也就是
$${\mathbf{g}}_i^{\mathrm{T}}{\mathbf{g}}_j=0(i<j)$$
将公式进行变换，可以得到
$$\begin{array}{rl}{\mathbf{g}}_{j+1}^{\mathrm{T}}{\mathbf{g}}_i& ={\left(\mathbf{A}{\mathbf{X}}_{j+1}\right)}^{\mathrm{T}}{\mathbf{g}}_i\\ & ={\left(\mathbf{A}\left({\mathbf{X}}_j+{\alpha }_j{\mathbf{d}}_j\right)\right)}^{\mathrm{T}}{\mathbf{g}}_i\\ & ={\left(\mathbf{A}{\mathbf{X}}_j\right)}^{\mathrm{T}}{\mathbf{g}}_i+{\left(\mathbf{A}{\alpha }_j{\mathbf{d}}_j\right)}^{\mathrm{T}}{\mathbf{g}}_i\end{array}$$
​ $=g_j^{\mathrm{T}}{\mathbf{g}}_i+{\alpha }_j{\mathbf{d}}_j^{\mathrm{T}}{\mathbf{A}}^{\mathrm{T}}{\mathbf{g}}_i(\mathbf{A}$ 是对称矩阵 $)$
$$=g_j^{\mathrm{T}}{g}_i+{\alpha }_j{d}_j^{\mathrm{T}}A{g}_i$$
对公式进行变换，可以得到
$${\mathbf{d}}_j^{\mathrm{T}}\mathbf{A}{\mathbf{g}}_i=\frac{1}{{\alpha }_j}\left[{\mathbf{g}}_{j+1}^{\mathrm{T}}{\mathbf{g}}_i-{\mathbf{g}}_j^{\mathrm{T}}{\mathbf{g}}_i\right]$$
此时我们发现公式的左边就是求解 $\beta$ 的分子，那么公式右边等于什么呢? 我们令 $j<i,$于是有:

• 当 $j=i-1$ 时，等式右边为 $\frac{1}{{\alpha }_j}\left[{\mathbf{g}}_i^{\mathrm{T}}{\mathbf{g}}_i-{\mathbf{g}}_{i-1}^{\mathrm{T}}{\mathbf{g}}_i\right],$ 可以看出当 $g_i$ 不等于 0 时，这部分公式不等于 0
• 当 $j<i-1$ 时，等式右边为 $\frac{1}{{\alpha }_j}\left[{\mathbf{g}}_{j+1}^{\mathrm{T}}{\mathbf{g}}_i-{\mathbf{g}}_j^{\mathrm{T}}{\mathbf{g}}_i\right],$ 由于 $j+1<i,j<i,$ 根据前面的推断，这部分公式为 $0_{}$
  所以可以得到，使用 Gram-Schmidt 方法可以保证更新方向共轭正交。对于第 $t$ 轮优化，我们只计算 ${\beta }_{t-1}$ 即可，其他的 $\beta$ 值均为 $0,$ 不用计算，这样计算量得到了极大的降低。

到这里我们完成了共轭梯度法的理论推导，由于算法比较复杂，我们需要对算法进行回顾：

• 最初，我们希望我们的更新方向与误差正交，从而得出了 $d_i^{\mathrm{T}}{e}_i=0$ 。
• 由于误差无法提前知道，直接计算正交比较困难，所以我们将正交改为共轭正交，于是公式变为 $d_i^{\mathrm{T}}A{e}_i=0_{\circ }$
• 我们希望使用函数的梯度作为优化方向的基础，每一次求出当前的梯度，就用梯度减去之前更新方向的成分，以保证每一次的优化都保持共轭正交。
• 我们采用 Gram-Schmidt 方法确保更新方向能够维持共轭正交的性质，根据证明，我们发现在使用梯度时，每一轮只需要计算一个 $\beta$ 值就能实现共轭正交的效果，计算量大大降低。 $\beta$ 的求解方法为 ${\beta }_{t,l}=-\frac{d_1^{\mathrm{T}}A{u}_l}{d_l^{T}A{d}_l}$ 。
• 完成了更新方向的计算，我们根据前面的推论，计算每一步的步长: ${\alpha }_t=$ $\frac{r_{\mathrm{t}}^{\mathrm{T}}{g}_t}{r_t^{\mathrm{T}}A{r}_t}$

然后我们在这里放一下共轭梯度的实现代码，代码来自于Baselines项目：

import numpy as np


def conjugate_gradient(f_ax, b_vec, cg_iters=10, callback=None, verbose=False, residual_tol=1e-10):
    """
    conjugate gradient calculation (Ax = b), bases on
    https://epubs.siam.org/doi/book/10.1137/1.9781611971446 Demmel p 312

    :param f_ax: (function) The function describing the Matrix A dot the vector x
                 (x being the input parameter of the function)
    :param b_vec: (numpy float) vector b, where Ax = b
    :param cg_iters: (int) the maximum number of iterations for converging
    :param callback: (function) callback the values of x while converging
    :param verbose: (bool) print extra information
    :param residual_tol: (float) the break point if the residual is below this value
    :return: (numpy float) vector x, where Ax = b
    """
    first_basis_vect = b_vec.copy()  # the first basis vector
    residual = b_vec.copy()  # the residual
    x_var = np.zeros_like(b_vec)  # vector x, where Ax = b
    residual_dot_residual = residual.dot(residual)  # L2 norm of the residual

    fmt_str = "%10i %10.3g %10.3g"
    title_str = "%10s %10s %10s"
    if verbose:
        print(title_str % ("iter", "residual norm", "soln norm"))

    for i in range(cg_iters):
        if callback is not None:
            callback(x_var)
        if verbose:
            print(fmt_str % (i, residual_dot_residual, np.linalg.norm(x_var)))
        z_var = f_ax(first_basis_vect)
        v_var = residual_dot_residual / first_basis_vect.dot(z_var)
        x_var += v_var * first_basis_vect
        residual -= v_var * z_var
        new_residual_dot_residual = residual.dot(residual)
        mu_val = new_residual_dot_residual / residual_dot_residual
        first_basis_vect = residual + mu_val * first_basis_vect

        residual_dot_residual = new_residual_dot_residual
        if residual_dot_residual < residual_tol:
            break

    if callback is not None:
        callback(x_var)
    if verbose:
        print(fmt_str % (i + 1, residual_dot_residual, np.linalg.norm(x_var)))
    return x_var

关于共轭梯度法的介绍到这就告一段落了，接下来我们回顾一下我们前面策略提升的问题：

针对这个问题，我们利用共轭梯度法将将优化目标中的约束条件转化为优化目标的一部分，从而将对目标函数在约束下进行优化的过程等价为直接优化对偶目标函数的行为。原始的优化目标在引入拉格朗日乘子后可以转化为如下的形式：

​ $$\mathcal{L}\left({\theta }^{\prime },\lambda \right)=\sum _t{E}_{s_t\sim p_o\left(s_t\right)}\left[E_{a_t\sim {\pi }_{\theta }\left(a_t,\mid s_t\right)}\left[\frac{{\pi }_{{\theta }^{\prime }}\left({\mathbf{a}}_t\mid {\mathbf{s}}_t\right)}{{\pi }_{\theta }\left({\mathbf{a}}_t\mid {\mathbf{s}}_t\right)}{\gamma }^t{A}^{{\pi }_{\theta }}\left({\mathbf{s}}_t,{\mathbf{a}}_t\right)\right]\right]-\lambda \left(D_{\mathrm{K}\mathrm{L}}\left({\pi }_{{\theta }^{\prime }}\left({\mathbf{a}}_t\mid {\mathbf{s}}_t\right)\Vert {\pi }_{\theta }\left({\mathbf{a}}_t\mid {\mathbf{s}}_t\right)\right)-ϵ\right)$$

进而可以进行迭代求解：

从直观上理解，将约束引入目标函数后, 乘子 $\lambda$ 也是一个需要优化的参数。当约束违反程度大的时候, 后面拉格朗日项就变成特别大的负项，如果需要最大化这个目标，则需要将乘子变大，反之依然。通过对这个乘子的调节，进而修正约束部分的重要性，从而达到自适应优化的目标。

上一篇：强化学习的学习之路（四十四）2021-02-13 Monotonic Improvement with KL Divergence
下一篇：强化学习的学习之路（四十六）2021-02-15自然梯度法实现策略上的单调提升(Monotonic Improvement with Natural gradient descent)