共轭梯度法的推导

在数值线性代数中，共轭梯度法是一种求解对称正定线性方程组

$\boldsymbol{Ax}=\boldsymbol{b}$

的迭代方法。共轭梯度法可以从不同的角度推导而得，包括作为求解最优化问题的共轭方向法的特例，以及作为求解特征值问题的 Arnoldi/Lanczos 迭代的变种。

本条目记述这些推导方法中的重要步骤。

从 Arnoldi/Lanczos 迭代推导

共轭梯度法可以看作 Arnoldi/Lanczos 迭代应用于求解线性方程组时的一个变种。

一般 Arnoldi 方法

Arnoldi 迭代从一个向量 $\boldsymbol{r}_0$ 开始，通过定义 $\boldsymbol{v}_i=\boldsymbol{w}_i/\lVert\boldsymbol{w}_i\rVert_2$，其中

$\boldsymbol{w}_i=\begin{cases} \boldsymbol{r}_0 & \text{if }i=1\text{,}\\ \boldsymbol{Av}_{i-1}-\sum_{j=1}^{i-1}(\boldsymbol{v}_j^\mathrm{T}\boldsymbol{Av}_{i-1})\boldsymbol{v}_j & \text{if }i>1\text{,} \end{cases}$

逐步建立 Krylov 子空间

$\mathcal{K}(\boldsymbol{A},\boldsymbol{r}_0)=\{\boldsymbol{r}_0,\boldsymbol{Ar}_0,\boldsymbol{A}^2\boldsymbol{r}_0,\ldots\}$

的一组标准正交基 $\{\boldsymbol{v}_1,\boldsymbol{v}_2,\boldsymbol{v}_3,\ldots\}$。

换言之，对于 $i>1$，$\boldsymbol{v}_i$ 由将 $\boldsymbol{Av}_{i-1}$ 相对于 $\{\boldsymbol{v}_1,\boldsymbol{v}_2,\ldots,\boldsymbol{v}_{i-1}\}$ 进行 Gram-Schmidt 正交化然后归一化得到。

写成矩阵形式，迭代过程可以表示为

$\boldsymbol{AV}_i=\boldsymbol{V}_{i+1}\boldsymbol{\tilde{H}}_i\text{,}$

其中

$\begin{align} \boldsymbol{V}_i&=\begin{bmatrix} \boldsymbol{v}_1 & \boldsymbol{v}_2 & \cdots & \boldsymbol{v}_i \end{bmatrix}\text{,}\\ \boldsymbol{\tilde{H}}_i&=\begin{bmatrix} h_{11} & h_{12} & h_{13} & \cdots & h_{1,i}\\ h_{21} & h_{22} & h_{23} & \cdots & h_{2,i}\\ & h_{32} & h_{33} & \cdots & h_{3,i}\\ & & \ddots & \ddots & \vdots\\ & & & h_{i,i-1} & h_{i,i}\\ & & & & h_{i+1,i} \end{bmatrix}=\begin{bmatrix} \boldsymbol{H}_i\\ h_{i+1,i}\boldsymbol{e}_i^\mathrm{T} \end{bmatrix}\text{,}\\ h_{ji}&=\begin{cases} \boldsymbol{v}_j^\mathrm{T}\boldsymbol{Av}_i & \text{if }j\leq i\text{,}\\ \lVert\boldsymbol{w}_{i+1}\rVert_2 & \text{if }j=i+1\text{,}\\ 0 & \text{if }j>i+1\text{.} \end{cases} \end{align}$

当应用于求解线性方程组时，Arnoldi 迭代对应于初始解 $\boldsymbol{x}_0$ 的残量 $\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}_0$ 开始迭代，在每一步迭代之后计算 $\boldsymbol{y}_i=\boldsymbol{H}_i^{-1}(\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{e}_1)$ 和新的近似解 $\boldsymbol{x}_i=\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{y}_i$.

直接 Lanzcos 方法

在余下的讨论中，我们假定 $\boldsymbol{A}$ 是对称正定矩阵。由于 $\boldsymbol{A}$ 的对称性, 上 Hessenberg 矩阵 $\boldsymbol{H}_i=\boldsymbol{V}_i^\mathrm{T}\boldsymbol{AV}_i$ 变成对阵三对角矩阵。于是它可以被更明确地表示为

$\boldsymbol{H}_i=\begin{bmatrix} a_1 & b_2\\ b_2 & a_2 & b_3\\ & \ddots & \ddots & \ddots\\ & & b_{i-1} & a_{i-1} & b_i\\ & & & b_i & a_i \end{bmatrix}\text{.}$

这使得迭代中的 $\boldsymbol{v}_i$ 有一个简短的三项递推式。Arnoldi 迭代从而简化为 Lanczos 迭代。

由于 $\boldsymbol{A}$ 对称正定，$\boldsymbol{H}_i$ 同样也对称正定。因此，$\boldsymbol{H}_i$ 可以通过不选主元的 LU 分解分解为

$\boldsymbol{H}_i=\boldsymbol{L}_i\boldsymbol{U}_i=\begin{bmatrix} 1\\ c_2 & 1\\ & \ddots & \ddots\\ & & c_{i-1} & 1\\ & & & c_i & 1 \end{bmatrix}\begin{bmatrix} d_1 & b_2\\ & d_2 & b_3\\ & & \ddots & \ddots\\ & & & d_{i-1} & b_i\\ & & & & d_i \end{bmatrix}\text{,}$

其中 $c_i$ 和 $d_i$ 有简单的递推式：

$\begin{align} c_i&=b_i/d_{i-1}\text{,}\\ d_i&=\begin{cases} a_1 & \text{if }i=1\text{,}\\ a_i-c_ib_i & \text{if }i>1\text{.} \end{cases} \end{align}$

改写 $\boldsymbol{x}_i=\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{y}_i$ 为

$\begin{align} \boldsymbol{x}_i&=\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{H}_i^{-1}(\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{e}_1)\\ &=\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{U}_i^{-1}\boldsymbol{L}_i^{-1}(\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{e}_1)\\ &=\boldsymbol{x}_0+\boldsymbol{P}_i\boldsymbol{z}_i\text{,} \end{align}$

其中

$\begin{align} \boldsymbol{P}_i&=\boldsymbol{V}_{i}\boldsymbol{U}_i^{-1}\text{,}\\ \boldsymbol{z}_i&=\boldsymbol{L}_i^{-1}(\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{e}_1)\text{.} \end{align}$

此时需要观察到

$\begin{align} \boldsymbol{P}_i&=\begin{bmatrix} \boldsymbol{P}_{i-1} & \boldsymbol{p}_i \end{bmatrix}\text{,}\\ \boldsymbol{z}_i&=\begin{bmatrix} \boldsymbol{z}_{i-1}\\ \zeta_i \end{bmatrix}\text{.} \end{align}$

实际上，$\boldsymbol{p}_i$ 和 $\zeta_i$ 同样有简短的递推式：

$\begin{align} \boldsymbol{p}_i&=\frac{1}{d_i}(\boldsymbol{v}_i-b_i\boldsymbol{p}_{i-1})\text{,}\\ \alpha_i&=-c_i\zeta_{i-1}\text{.} \end{align}$

通过这个形式，我们得到 $\boldsymbol{x}_i$ 的一个简单的递推式：

$\begin{align} \boldsymbol{x}_i&=\boldsymbol{x}_0+\boldsymbol{P}_i\boldsymbol{z}_i\\ &=\boldsymbol{x}_0+\boldsymbol{P}_{i-1}\boldsymbol{z}_{i-1}+\zeta_i\boldsymbol{p}_i\\ &=\boldsymbol{x}_{i-1}+\zeta_i\boldsymbol{p}_i\text{.} \end{align}$

以上的递推关系立即导出比共轭梯度法稍微更复杂的直接 Lanczos 方法。

从正交性和共轭性导出共轭梯度法

如果我们允许缩放 $\boldsymbol{p}_i$ 并通过常数因子补偿缩放的系数，我们可能可以的到以下形式的更简单的递推式：

$\begin{align} \boldsymbol{x}_i&=\boldsymbol{x}_{i-1}+\alpha_{i-1}\boldsymbol{p}_{i-1}\text{,}\\ \boldsymbol{r}_i&=\boldsymbol{r}_{i-1}-\alpha_{i-1}\boldsymbol{Ap}_{i-1}\text{,}\\ \boldsymbol{p}_i&=\boldsymbol{r}_i+\beta_{i-1}\boldsymbol{p}_{i-1}\text{.} \end{align}$

作为简化的前提，我们现在推导 $\boldsymbol{r}_i$ 的正交性和 $\boldsymbol{p}_i$ 的共轭性，即对于 $i\neq j$,

$\begin{align} \boldsymbol{r}_i^\mathrm{T}\boldsymbol{r}_j&=0\text{,}\\ \boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_j&=0\text{.} \end{align}$

各个残量之间满足正交性的原因是 $\boldsymbol{r}_i$ 实际上可由 $\boldsymbol{v}_{i+1}$ 乘上一个系数而得。这是因为对于 $i=0$，$\boldsymbol{r}_0=\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{v}_1$，对于 $i>0$，

$\begin{align} \boldsymbol{r}_i&=\boldsymbol{b}-\boldsymbol{Ax}_i\\ &=\boldsymbol{b}-\boldsymbol{A}(\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{y}_i)\\ &=\boldsymbol{r}_0-\boldsymbol{AV}_i\boldsymbol{y}_i\\ &=\boldsymbol{r}_0-\boldsymbol{V}_{i+1}\boldsymbol{\tilde{H}}_i\boldsymbol{y}_i\\ &=\boldsymbol{r}_0-\boldsymbol{V}_i\boldsymbol{H}_i\boldsymbol{y}_i-h_{i+1,i}(\boldsymbol{e}_i^\mathrm{T}\boldsymbol{y}_i)\boldsymbol{v}_{i+1}\\ &=\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{v}_1-\boldsymbol{V}_i(\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{e}_1)-h_{i+1,i}(\boldsymbol{e}_i^\mathrm{T}\boldsymbol{y}_i)\boldsymbol{v}_{i+1}\\ &=-h_{i+1,i}(\boldsymbol{e}_i^\mathrm{T}\boldsymbol{y}_i)\boldsymbol{v}_{i+1}\text{.}\end{align}$

要得到 $\boldsymbol{p}_i$ 的共轭性，只需证明 $\boldsymbol{P}_i^\mathrm{T}\boldsymbol{AP}_i$ 是对角矩阵：

$\begin{align} \boldsymbol{P}_i^\mathrm{T}\boldsymbol{AP}_i&=\boldsymbol{U}_i^{-\mathrm{T}}\boldsymbol{V}_i^\mathrm{T}\boldsymbol{AV}_i\boldsymbol{U}_i^{-1}\\ &=\boldsymbol{U}_i^{-\mathrm{T}}\boldsymbol{H}_i\boldsymbol{U}_i^{-1}\\ &=\boldsymbol{U}_i^{-\mathrm{T}}\boldsymbol{L}_i\boldsymbol{U}_i\boldsymbol{U}_i^{-1}\\ &=\boldsymbol{U}_i^{-\mathrm{T}}\boldsymbol{L}_i \end{align}$

是对称的下三角矩阵，因而必然是对角矩阵。

现在我们可以单纯由$\boldsymbol{r}_i$ 的正交性和 $\boldsymbol{p}_i$ 的共轭性推导相对于缩放后的 $\boldsymbol{p}_i$ 的常数因子 $\alpha_i$ 和 $\beta_i$。

由于 $\boldsymbol{r}_i$ 的正交性，必然有 $\boldsymbol{r}_{i+1}^\mathrm{T}\boldsymbol{r}_i=(\boldsymbol{r}_i-\alpha_i\boldsymbol{Ap}_i)^\mathrm{T}\boldsymbol{r}_i=0$。于是

$\begin{align} \alpha_i&=\frac{\boldsymbol{r}_i^\mathrm{T}\boldsymbol{r}_i}{\boldsymbol{r}_i^\mathrm{T}\boldsymbol{Ap}_i}\\ &=\frac{\boldsymbol{r}_i^\mathrm{T}\boldsymbol{r}_i}{(\boldsymbol{p}_i-\beta_{i-1}\boldsymbol{p}_{i-1})^\mathrm{T}\boldsymbol{Ap}_i}\\ &=\frac{\boldsymbol{r}_i^\mathrm{T}\boldsymbol{r}_i}{\boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_i}\text{.} \end{align}$

类似地，由于 $\boldsymbol{p}_i$ 的共轭性，必然有 $\boldsymbol{p}_{i+1}^\mathrm{T}\boldsymbol{Ap}_i=(\boldsymbol{r}_{i+1}+\beta_i\boldsymbol{p}_i)^\mathrm{T}\boldsymbol{Ap}_i=0$。于是

$\begin{align} \beta_i&=-\frac{\boldsymbol{r}_{i+1}^\mathrm{T}\boldsymbol{Ap}_i}{\boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_i}\\ &=-\frac{\boldsymbol{r}_{i+1}^\mathrm{T}(\boldsymbol{r}_i-\boldsymbol{r}_{i+1})}{\alpha_i\boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_i}\\ &=\frac{\boldsymbol{r}_{i+1}^\mathrm{T}\boldsymbol{r}_{i+1}}{\boldsymbol{r}_i^\mathrm{T}\boldsymbol{r}_i}\text{.} \end{align}$

推导至此完成。

参考文献

1. Hestenes, M. R.; Stiefel, E. Methods of conjugate gradients for solving linear systems (PDF). Journal of Research of the National Bureau of Standards. 1952.12, 49 (6). 
2. Saad, Y. Chapter 6: Krylov Subspace Methods, Part I//Iterative methods for sparse linear systems 2nd. SIAM. 2003. ISBN 978-0898715347.