共轭梯度法（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《最优化理论》课程讲义