共轭梯度法(Conjugate Gradient Method)

共轭梯度法是求解大型线性方程组的一种很有效的方法,也可以通过改进用来求解非线性优化问题。下面先介绍线性线性共轭梯度法,再介绍非线性共轭梯度法。

线性共轭梯度法

共轭梯度法是一个迭代求解的算法。假如我们要求解线性方程组 $Ax=b$,其中 $A$ 是 $n\times n$ 的对称正定矩阵。这实际上等价于一个最小化问题,即 $$Ax=b \Leftrightarrow \min \phi(x)=\cfrac{1}{2}x^TAx-b^Tx$$ 它们有相同的解

由于 $\phi(x)$ 是凸函数,所以它的局部最优点也是全局最优点,记为 $x^{\star}$
则有 $\phi’(x^{\star}) = 0$
又因为 $\phi’(x) = Ax - b$,所以优化问题和线性方程组有相同的解。

共轭方向法(Conjugate Direction Method)

共轭梯度法其实是特殊情况下的共轭方向法,所以我们先介绍共轭方向法。
如果一个非零向量的集合 $\{p_0,p_1,\cdots,p_l\}$ 满足以下条件,我们就说它关于对称正定矩阵 $A$ 是共轭的 $$p_i^TAp_j=0 \quad \text{for all}\quad i\neq j$$很容易验证这样的一组向量也是线性独立的。
下面是共轭方向法的过程:

  • 给定一个初始解 $x_0\in\mathbb{R}^n$ 和一组共轭方向 $\{p_0,p_1,\cdots,p_{n-1}\}$
  • 通过式 $x_{k+1}=x_k+\alpha_kp_k$ 生成一个解序列 $\{x_k\}$
  • 其中 $\alpha_k$ 是二次方程 $\phi(\cdot)$ 沿 $x_k+\alpha_kp_k$ 方向的极小值点,$\alpha_k$ 可以这样计算 $$\alpha_k=-\cfrac{r_k^Tp_k}{p_k^TAp_k}$$

记 $f(\alpha)=\phi(x_k+\alpha p_k)$,要求解使 $\phi(\cdot)$ 最小的 $\alpha_k$,可以令 $f’(\alpha)=0$ 来求解
$f(\alpha)=\cfrac{1}{2}(x_k+\alpha p_k)^TA(x_k+\alpha p_k)-b^T(x_k+\alpha p_k)$
对于第一项,根据复合求导法则,外部求导结果为 $A(x_k+\alpha p_k)$,内部求导结果为 $p_k^T$。所以,$f’(\alpha)=p_k^TA(x_k+\alpha p_k) -b^Tp_k$,令 $f’(\alpha)=0$,则可以得到上面的解,其中 $r_k=\phi’(x_k)=Ax_k-b^Tx_k$

二次终止定理

定理: 对于任一 $x_0\in\mathbb{R}^n$,用共轭方向法生成的序列 $\{x_k\}$ 最多用 $n$ 步就可以收敛到线性方程组 $Ax=b$ 的解 $x^{\star}$
证明: 因为共轭方向 $\{p_i\}$ 是线性独立的,它们一定可以张成整个 $n$ 维空间 $\mathbb{R}^n$,则 $$x^{\star}-x_0=\sigma_0p_0+\sigma_1p_1+\cdots+\sigma_{n-1}p_{n-1}$$ 给上式左乘 $p_k^TA$,并利用共轭的性质,可以得到 $$\sigma_k=\cfrac{p_k^TA(x^{\star}-x_0)}{p_k^TAp_k}$$ 如果 $x_k$ 是由共轭方向法生成的,则有 $$x_k=x_0+\alpha_0p_0+\alpha_1p_1+\cdots+\alpha_{n-1}p_{n-1}$$
给上式左乘 $p_k^TA$,并利用共轭的性质,则有 $$p_k^TA(x_k-x_0)=0$$ 因此有 $$p_k^TA(x^{\star}-x_0)=p_k^TA(x^{\star}-x_k)=p_k^TA(b-Ax_k)=-p_k^Tr_k$$ 则有 $$\sigma_k=-\cfrac{p_k^Tr_k}{p_k^TAp_k}=\alpha_k$$ 即证明了通过共轭方向法可以在最多 $n$ 步时得到解 $x^{\star}$

共轭方向法的实质

如果矩阵 $A$ 是对角矩阵,则函数 $\phi(\cdot)$ 的等值线是轴与坐标轴平行的椭圆,如下图:


这时我们可以沿每个坐标轴方向进行一维最小化,则可以在 $n$ 步之内找到最优点,每一步都在一个坐标轴方向上达到了最优。

如果矩阵 $A$ 不是对角矩阵,即椭圆的轴与坐标轴不平行,则可以通过修正使得新的矩阵变成对角矩阵
定义新的变量 $\hat{x}=S^{-1}x$,其中 $S=[p_0\quad p_1\quad \cdots\quad p_{n-1}]$ 是由 $A$ 的一组共轭向量 $\{p_0,p_1,\cdots,p_{n-1}\}$ 组成的矩阵,定义 $$\hat{\phi}(\hat{x})=\phi(S\hat{x})=\cfrac{1}{2}\hat{x}^T(S^TAS)\hat{x}-(S^Tb)^T\hat{x}$$ 由共轭的性质可知,矩阵 $S^TAS$ 是一个对角矩阵,那么就可以在最多 $n$ 步一维最小化后找到最优解。
反过来思考共轭方向法,它可以看作是沿椭圆的每个轴进行一维最小化,在 $n$ 步内达到最优点(个人感觉,有待进一步研究)。

扩张子空间定理

定理:令 $x_0\in\mathbb{R}^n$ 为任意一个初始解,设序列 $\{x_k\}$ 是由共轭方向法产生的,则有$$r_k^Tp_i=0, \quad\forall i=0,\cdots,k-1$$ 并且 $x_k$ 是 $\phi(x)$ 在集合 $\{x|x=x_0+span\{p_0,p_1,\cdots,p_{k-1}\}\}$ 上的极小值点。

证明略。扩张子空间告诉我们每一次迭代后,残差 $r_k$ 与所有迭代过的方向都是正交的,且 $x_k$ 是在当前子空间中的最优点。这说明每一次迭代都在它的方向上达到了最优,故而最多 $n$ 次迭代后就能达到全局最优。

共轭梯度法

共轭梯度法实际上是共轭方向法,只不过共轭方向的选取方法比较特殊:只依赖前一次的方向 $p_{k-1}$ 来计算出下一次的方向 $p_k$
每一步的方向都是最速下降方向 $-r_k$ 和 $p_{k-1}$ 的线性组合:$$p_k=-r_k+\beta_kp_{k-1}$$ 其中 $\beta_k$ 可以由 $p_k$ 和 $p_{k-1}$ 关于 $A$ 共轭这个性质计算出来,即 $$\beta_k=\cfrac{r_k^TAp_{k-1}}{p_{k-1}^TAp_{k-1}}$$

$p_k$ 和 $p_{k-1}$ 关于 $A$ 共轭,则有 $p_kAp_{k-1}=0$
即 $(-r_k+\beta_kp_{k-1})Ap_{k-1}=0$
就可以解出来上面的结果

我们选择 $p_0=-r_0$
共轭梯度算法的过程如下所示:

给定初始解 $x_0$
令 $r_0\leftarrow Ax_0-b,\quad p_0\leftarrow -r_0,\quad k\leftarrow 0$
$\mathbf{while} \quad r_k\neq 0$
$$\begin{array}{r,l}\alpha_k&\leftarrow&-\cfrac{r_k^Tp_k}{p_k^TAp_k} \\ x_{k+1}&\leftarrow&x_k+\alpha_kp_k \\ r_{k+1}&\leftarrow&Ax_{k+1}-b \\ \beta_{k+1}&\leftarrow&\cfrac{r_{k+1}^TAp_k}{p_k^TAp_k} \\ p_{k+1} &=& -r_{k+1}+\beta_{k+1}p_k \\ k&\leftarrow&k+1 \end{array}$$
$\mathbf{end(while)}$

在实际中,常常用下面两个公式更新 $\alpha$ 和 $\beta$:$$\alpha_k\leftarrow-\cfrac{r_k^Tr_k}{p_k^TAp_k}\qquad \beta_{k+1}\leftarrow\cfrac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k}$$

小于 $n$ 次的情况

上面说的都是最多 $n$ 次就能找到最优解,那么迭代次数与什么有关呢?
我们可以看出,共轭方向法每次都在一个方向上达到最优,但如果椭球在这个方向上没有分量呢?(相当于一个高维空间的低维物体,高维分量为0)那么无论高维取什么值,函数都在扩充空间上都是最优的(因为在低维时已达到最优,而高维分量为0,对结果不影响)。
有定理表明,当 $A$ 有 $r$ 个非零奇异值时,算法最多进行 $r$ 次迭代就能达到最优点。
为什么还有“最多”呢?因为有可能我们选取的解恰好有可能在某个方向上已经是最优的了。设想最优解是 $(0,0,1)^T$,我们一开始选择了 $(0,0,0)^T$,那么算法只需要一次就可以完成。

非线性共轭梯度法

上面的算法可以用来对凸二次函数进行最小化,对它进行修改可以使其适用于一般的非线性函数 $f$。

The Fletcher-Reeves Method

对上面的算法主要做了两处改动:

  • 步长 $\alpha_k$ 的选择:进行线搜索,得到一个 $f$ 在 $p_k$ 方向上的近似极小点
  • 残差$r$:改为非线性函数 $f$ 的梯度 $\nabla f$

当线搜索是精确线搜索时,则方向 $p_k$ 可以保证是下降方向;而当线搜索是非精确线搜索时,方向 $p_k$ 就不一定是下降方向了,此时可以限制 $\alpha_k$ 满足强Wolfe条件,则得到的方向 $p_k$ 一定是下降方向。

关于强Wolfe条件和线搜索的内容可以参考这篇文章

The Polak-Ribiere Method

PR-CG算法和FR-CG算法的不同在于选取 $\beta_k$ 的不同:$$\beta_{k+1}^{\text{FR}}=\cfrac{\nabla f_{k+1}^T\nabla f_{k+1}}{\nabla f_k^T\nabla f_k}\qquad \beta_{k+1}^{\text{PR}}=\cfrac{\nabla f_{k+1}^T(\nabla f_{k+1}-\nabla f_k)}{\nabla f_k^T\nabla f_k}$$ 当 $f$ 是严格凸的二次函数,且使用精确线搜索时,两种算法是相同的。
数值结果表明,PR-CG算法更加地高效,强健性也更好。

重启

在非线性共轭梯度法中,常常每过 $n$ 步就进行重启,重新设置 $\beta_k=0$,即选取最速下降方向。
一个常用的重启判定策略是:当连续两次的梯度值离正交很远时(即夹角偏离90度很大),进行一次重启。

参考资料

  1. NJU《最优化理论》课程讲义