Krylov 子空間法──線性方程的數值解法 (三)：共軛梯度法

本文的閱讀等級：高級

共軛梯度法 (conjugate gradient method) 是一個適用於實對稱正定矩陣的線性方程數值解法。顧名思義，共軛梯度法的核心是共軛 (conjugacy) 和梯度 (一階導數)。共軛能夠加快收斂，梯度則提供正交基底。因為這兩個特性，共軛梯度法的結構簡單優美，儲存量及運算量少，並且無須設定參數。對於大尺寸矩陣，我們往往無法使用直接法求解，譬如 Cholesky 分解 (見“Cholesky 分解”)，這時候可以採用以迭代方式計算的共軛梯度法。此外，對於大型非線性最佳化問題，共軛梯度法也是最有效的數值算法之一。

令 $A$ 為一個 $n\times n$ 階實對稱矩陣， $\mathbf{b}\in\mathbb{R}^n$ 為非零向量。假設 $\{\mathbf{p}_1,\ldots,\mathbf{p}_n\}$ 是 $\mathbb{R}^n$ 的一組基底。令 $P=\begin{bmatrix} \mathbf{p}_1&\cdots&\mathbf{p}_n \end{bmatrix}$ 。明顯地， $P$ 是 $n\times n$ 階可逆矩陣。線性方程 $A\mathbf{x}=\mathbf{b}$ 的精確解 $\mathbf{x}_\ast=A^{-1}\mathbf{b}$ 可表示為

$\displaystyle \mathbf{x}_\ast=\alpha_1\mathbf{p}_1+\cdots+\alpha_n\mathbf{p}_n=P\boldsymbol{\alpha}$ ，

其中 $\boldsymbol{\alpha}=(\alpha_1,\ldots,\alpha_n)^T$ 。線性方程 $A\mathbf{x}=\mathbf{b}$ 等號兩邊左乘 $P^T$ ，並代入上式可得等價式

$P^TAP\boldsymbol{\alpha}=P^T\mathbf{b}$ ，

其中係數矩陣 $P^TAP$ 仍為對稱矩陣， $\boldsymbol{\alpha}$ 是待解的向量。若 $P^TAP$ 是對角矩陣， $(P^TAP)_{ij}=\mathbf{p}_i^TA\mathbf{p}_j=0$ ， $i\neq j$ ，則

$\displaystyle \begin{bmatrix} \mathbf{p}_1^TA\mathbf{p}_1&&\\ &\ddots&\\ &&\mathbf{p}_n^TA\mathbf{p}_n \end{bmatrix}\begin{bmatrix} \alpha_1\\ \vdots\\ \alpha_n \end{bmatrix}=\begin{bmatrix} \mathbf{p}_1^T\mathbf{b}\\ \vdots\\ \mathbf{p}_n^T\mathbf{b} \end{bmatrix}$ 。

因為 $\mathbf{p}_i\neq\mathbf{0}$ ，若 $A$ 是正定矩陣，則 $\mathbf{p}_i^TA\mathbf{p}_i>0$ (見“特殊矩陣 (6)：正定矩陣”)。換句話說，正定矩陣保證可逆，故

$\displaystyle \alpha_i=\frac{\mathbf{p}_i^T\mathbf{b}}{\mathbf{p}_i^TA\mathbf{p}_i},~~i=1,\ldots,n$ ，

精確解可表示為

$\displaystyle \mathbf{x}_\ast=\sum_{i=1}^n\frac{\mathbf{p}_i^T\mathbf{b}}{\mathbf{p}_i^TA\mathbf{p}_i}\mathbf{p}_i$ 。

前述推演引申出下列定義：對於一個實對稱矩陣 $A$ ，若 $\mathbf{p}_i^TA\mathbf{p}_j=0$ ，我們說 $\mathbf{p}_i$ 和 $\mathbf{p}_j$ 對於 $A$ 共軛 (conjugate) 或 $A$ ─正交 ( $A$ -orthogonal)。本文考慮 $A$ 是正定矩陣，但此性質並不納入共軛的基本定義。若 $A=0$ ，則任何兩個向量皆為共軛。若 $A=I$ ，則共軛等價於正交。我們稱 $\{\mathbf{p}_1,\ldots,\mathbf{p}_k\}$ 對於 $A$ 為共軛向量集，若 $\mathbf{p}_i^TA\mathbf{p}_j=0$ ， $i\neq j$ (在不造成混淆的情況下，以下簡稱共軛向量集)。對於實對稱正定矩陣 $A$ ，若 $\{\mathbf{p}_1,\ldots,\mathbf{p}_k\}$ 為非零共軛向量集，則它們是線性獨立向量。證明於下：假設數組 $c_1,\ldots,c_k$ 使得 $c_1\mathbf{p}_1+\cdots+c_k\mathbf{p}_k=\mathbf{0}$ 。等號兩邊左乘 $\mathbf{p}_j^TA$ ， $\sum_{i=1}^kc_i\mathbf{p}_j^TA\mathbf{p}_i=c_j\mathbf{p}_j^TA\mathbf{p}_j=0$ 。因為 $\mathbf{p}_j^TA\mathbf{p}_j>0$ ，證明 $c_j=0$ 。在向量空間 $\mathbb{R}^n$ ，一旦獲得完整的非零共軛向量集 $\{\mathbf{p}_1,\ldots,\mathbf{p}_n\}$ ，透過矩陣─向量乘法與內積運算即可求出精確解。

如欲採用迭代法計算近似解，我們必須改寫精確解的表達式。給定初始猜測解 $\mathbf{x}_0$ ，以及非零共軛向量集 $\{\mathbf{p}_0,\mathbf{p}_1,\ldots,\mathbf{p}_{n-1}\}$ ，令

$\displaystyle \mathbf{x}_\ast-\mathbf{x}_0=\alpha_0\mathbf{p}_0+\alpha_1\mathbf{p}_1+\cdots+\alpha_{n-1}\mathbf{p}_{n-1}$ ，

且

$\mathbf{x}_j-\mathbf{x}_0=\alpha_0\mathbf{p}_0+\alpha_1\mathbf{p}_1+\cdots+\alpha_{j-1}\mathbf{p}_{j-1}$ ， $j>0$ 。

因此， $\mathbf{x}_n=\mathbf{x}_\ast$ 。下列遞迴式成立：

$\mathbf{x}_{j+1}=\mathbf{x}_j+\alpha_j\mathbf{p}_j,~~0\le j\le n-1$ 。

利用共軛性很容易解出 $\alpha_j$ 。考慮

$\displaystyle \mathbf{p}_j^TA(\mathbf{x}_\ast-\mathbf{x}_j)=\mathbf{p}_j^TA(\alpha_j\mathbf{p}_j+\cdots+\alpha_{n-1}\mathbf{p}_{n-1})=\alpha_j\mathbf{p}_j^TA\mathbf{p}_j$ ，

即得

$\displaystyle \alpha_j=\frac{\mathbf{p}_j^TA(\mathbf{x}_\ast-\mathbf{x}_j)}{\mathbf{p}_j^TA\mathbf{p}_j}=\frac{\mathbf{p}_j^T(\mathbf{b}-A\mathbf{x}_j)}{\mathbf{p}_j^TA\mathbf{p}_j}=\frac{\mathbf{p}_j^T\mathbf{r}_j}{\mathbf{p}_j^TA\mathbf{p}_j}$ ，

其中 $\mathbf{r}_j=\mathbf{b}-A\mathbf{x}_j$ ，稱為殘差 (residual)。

下面這個性質說明共軛向量集的價值所在。令 $\mathcal{P}_j=\hbox{span}\{\mathbf{p}_0,\mathbf{p}_1,\ldots,\mathbf{p}_{j-1}\}$ 。考慮二次型

$\displaystyle E(\mathbf{x})=\frac{1}{2}(\mathbf{x}_\ast-\mathbf{x})^TA(\mathbf{x}_\ast-\mathbf{x})$ 。

當 $\mathbf{x}=\mathbf{x}_\ast$ ， $E(\mathbf{x})$ 有極小值 $0$ 。限定 $\mathbf{x}-\mathbf{x}_0\in\mathcal{P}_j$ ，我們可以證明遞迴式產生的 $\mathbf{x}_j$ ， $\mathbf{x}_j-\mathbf{x}_0\in\mathcal{P}_j$ ，最小化 $E(\mathbf{x})$ 。令 $n\times j$ 階矩陣 $P_j=\begin{bmatrix} \mathbf{p}_0&\mathbf{p}_1&\cdots&\mathbf{p}_{j-1} \end{bmatrix}$ 。設 $\mathbf{x}=\mathbf{x}_0+P_j\mathbf{y}$ ，上面的約束優化問題等價於尋找 $\mathbf{y}\in\mathbb{R}^j$ 使最小化

$\displaystyle E(\mathbf{y})=\frac{1}{2}(\mathbf{x}_\ast-(\mathbf{x}_0+P_j\mathbf{y}))^TA(\mathbf{x}_\ast-(\mathbf{x}_0+P_j\mathbf{y}))$ 。

因為 $A^T=A$ ，求導可得 (見“矩陣導數”，SV-10)

$\displaystyle \frac{\partial E}{\partial\mathbf{y}}=-P_j^TA(\mathbf{x}_\ast-(\mathbf{x}_0+P_j\mathbf{y}))=-P_j^T(\mathbf{b}-A(\mathbf{x}_0+P_j\mathbf{y}))$ 。

最小化 $E(\mathbf{y})$ 的充要條件為 $P_j^T(\mathbf{b}-A(\mathbf{x}_0+P_j\mathbf{y}))=P_j^T(\mathbf{b}-A\mathbf{x})=0$ ，即 $\mathbf{b}-A\mathbf{x}$ 屬於 $P_j^T$ 的零空間 (nullspace) $N(P_j^T)$ ，也就是說， $\mathbf{b}-A\mathbf{x}\perp\mathcal{P}_j$ ，因為 $N(P_j^T)=\mathcal{P}_j^\perp$ (見“線性代數基本定理 (二)”)。我們用歸納法來證明 $\mathbf{r}_j\perp\mathcal{P}_j$ (即 $\mathbf{b}-A\mathbf{x}_j\perp\mathcal{P}_j$ )， $0\le j\le n-1$ 。若 $j=0$ ， $\mathcal{P}_0$ 是空集合，命題自然成立。假設 $\mathbf{r}_j\perp\mathcal{P}_j$ 。寫出 $\mathbf{r}_j$ 的遞迴式，

$\mathbf{r}_{j+1}=\mathbf{b}-A\mathbf{x}_{j+1}=\mathbf{b}-A(\mathbf{x}_j+\alpha_j\mathbf{p}_j)=\mathbf{r}_j-\alpha_jA\mathbf{p}_j$ 。

等號兩邊左乘 $\mathbf{p}_j^T$ ，再套用 $\alpha_j$ 的公式，

$\mathbf{p}_j^T\mathbf{r}_{j+1}=\mathbf{p}_j^T\mathbf{r}_j-\alpha_j\mathbf{p}_j^TA\mathbf{p}_j=0$ 。

對於 $i<j$ ， $\mathbf{p}_i\in\mathcal{P}_j$ ，使用先前假設 $\mathbf{r}_j\perp\mathbf{p}_i$ ，以及共軛性，可得

$\mathbf{p}_i^T\mathbf{r}_{j+1}=\mathbf{p}_i^T\mathbf{r}_j-\alpha_j\mathbf{p}_i^TA\mathbf{p}_j=0$ 。

合併以上結果，證明 $\mathbf{r}_{j+1}\perp\mathcal{P}_{j+1}$ 。

因為 $\mathcal{P}_1\subset\mathcal{P}_2\subset\cdots\subset\mathcal{P}_n=\mathbb{R}^n$ ，由上述性質可推論 $E(\mathbf{x}_0)\ge E(\mathbf{x}_1)\ge\cdots\ge E(\mathbf{x}_n)=0$ 。理論上，使用至多 $n$ 次計算可得精確解 $\mathbf{x}_\ast$ 。接下來的工作是找出適當的非零共軛向量集 $\{\mathbf{p}_0,\mathbf{p}_1,\ldots,\mathbf{p}_{j-1}\}$ 使得子空間 $\mathcal{P}_j$ 包含較接近 $\mathbf{x}_\ast$ 的近似解。最佳化理論提供了一些線索。將目標函數改寫如下：

$\displaystyle \begin{aligned} E(\mathbf{x})&=\frac{1}{2}(\mathbf{x}_\ast-\mathbf{x})^TA(\mathbf{x}_\ast-\mathbf{x})\\ &=\frac{1}{2}\mathbf{x}^TA\mathbf{x}-\mathbf{x}^TA\mathbf{x}_\ast+\frac{1}{2}\mathbf{x}_\ast^TA\mathbf{x}_\ast\\ &=\frac{1}{2}\mathbf{x}^TA\mathbf{x}-\mathbf{x}^T\mathbf{b}+\frac{1}{2}\mathbf{x}_\ast^TA\mathbf{x}_\ast. \end{aligned}$

最小化 $E(\mathbf{x})$ 等價於最小化

$\displaystyle f(\mathbf{x})=\frac{1}{2}\mathbf{x}^TA\mathbf{x}-\mathbf{x}^T\mathbf{b}$ 。

求導可得

$\displaystyle \nabla f=\frac{\partial f}{\partial \mathbf{x}}=A\mathbf{x}-\mathbf{b}$ ，

故線性方程 $A\mathbf{x}=\mathbf{b}$ 的唯一解等於二次優化問題 $\min_{\mathbf{x}} f(\mathbf{x})$ 的唯一解。給定初始猜測解 $\mathbf{x}_0$ ，梯度下降法 (gradient descent) 使用下列遞迴公式更新近似解 (見“最佳化理論與正定矩陣”)：

$\displaystyle \mathbf{x}_{j+1}=\mathbf{x}_j-\eta_j\nabla f(\mathbf{x}_j)=\mathbf{x}_j+\eta_j\mathbf{r}_j,~~j=0,1,2,\ldots$

其中 $-\nabla f(\mathbf{x}_j)=\mathbf{b}-A\mathbf{x}_j=\mathbf{r}_j$ 表示朝精確解的前進方向， $\eta_j$ 是一個夠小的正數，稱為步長 (step size)，使得 $f(\mathbf{x}_0)\ge f(\mathbf{x}_1)\ge f(\mathbf{x}_2)\ge\cdots$ 。若 $\mathbf{r}_j=\mathbf{0}$ ，則 $\mathbf{x}_j$ 即為線性方程的解。我們可以設 $\eta_j=\arg\min_{\eta>0}f(\mathbf{x}_j+\eta\mathbf{r}_j)$ ^[1]，這個方法稱為最陡下降法 (steepest descent)。然而，當 $A$ 的條件數 (condition number) 大時， $\mathbf{r}_j$ 未必指向靠近精確解的方向，最陡下降法的收斂速度因此變得很慢 (見“條件數”)。為了克服這個問題，共軛梯度法提出一個非常精巧的想法：以共軛向量 $\mathbf{p}_1,\mathbf{p}_2,\ldots$ 取代殘差 (負梯度) $\mathbf{r}_1,\mathbf{r}_2,\ldots$ ，這麼做的原因是限定 $\mathbf{x}-\mathbf{x}_0\in\mathcal{P}_j$ ，遞迴式 $\mathbf{x}_j=\mathbf{x}_{j-1}+\alpha_{j-1}\mathbf{p}_{j-1}$ 得到的 $\mathbf{x}_j$ 可使 $f(\mathbf{x})$ 有最小值。

下面介紹共軛梯度法。給定初始猜測解 $\mathbf{x}_0$ ，令 $\mathbf{p}_0=\mathbf{r}_0=\mathbf{b}-A\mathbf{x}_0$ 。重抄近似解和殘差的遞迴公式：

$\begin{aligned} \mathbf{x}_{j+1}&=\mathbf{x}_j+\alpha_j\mathbf{p}_j\\ \mathbf{r}_{j+1}&=\mathbf{r}_j-\alpha_jA\mathbf{p}_j,\end{aligned}$

其中 $\alpha_j={\mathbf{p}_j^T\mathbf{r}_j}/{\mathbf{p}_j^TA\mathbf{p}_j}$ 。設想共軛向量 $\mathbf{p}_{j+1}$ 可表示為 $\mathbf{r}_0,\mathbf{r}_1,\ldots,\mathbf{r}_{j+1}$ 的線性組合，但為了減少計算，我們採用下列遞迴算式：

$\mathbf{p}_{j+1}=\mathbf{r}_{j+1}+\beta_j\mathbf{p}_j$ ，

其中 $\beta_j$ 稱為動量 (momentum)，滿足共軛性 $0=\mathbf{p}_{j+1}^TA\mathbf{p}_j=\mathbf{r}_{j+1}^TA\mathbf{p}_j+\beta_j\mathbf{p}_j^TA\mathbf{p}_j$ ，由此解出 $\beta_j=-{\mathbf{r}_{j+1}^TA\mathbf{p}_j}/{\mathbf{p}_j^TA\mathbf{p}_j}$ 。如果 $\beta_j=0$ ，則共軛梯度法退化為最陡下降法。

定義 Krylov 子空間

$\displaystyle \mathcal{K}_j=\hbox{span}\{\mathbf{r}_0,A\mathbf{r}_0,\ldots,A^{j-1}\mathbf{r}_0\}$

且 $A(\mathcal{K}_j)=\{A\mathbf{y}\vert\mathbf{y}\in\mathcal{K}_j\}$ 。若共軛梯度法於 $\mathbf{x}_{j-1}$ 尚未停止，下列性質成立：

$\hbox{span}\{\mathbf{r}_0,\mathbf{r}_1,\ldots,\mathbf{r}_{j-1}\}=\mathcal{K}_j$
$\mathcal{P}_j=\hbox{span}\{\mathbf{p}_0,\mathbf{p}_1,\ldots,\mathbf{p}_{j-1}\}=\mathcal{K}_j$
$\mathbf{p}_i^TA\mathbf{p}_{j-1}=0$ ， $i<j-1$ 。

根據性質(1)和(2)， $\mathbf{r}_j\perp\mathcal{P}_j$ 可推得 $\mathbf{r}_j\perp\mathbf{r}_i$ ， $i<j$ ，也就是說，殘差 $\{\mathbf{r}_0,\mathbf{r}_1,\ldots,\mathbf{r}_{j-1}\}$ 為正交集。性質(3)表明 $\{\mathbf{p}_0,\mathbf{p}_1,\ldots,\mathbf{p}_{j-1}\}$ 為共軛向量集。

我們用歸納法同時證明這三個性質。當 $j=1$ ，性質(1-3)顯然成立。假設對於 $j$ ，性質(1-3)成立。考慮 $\mathbf{r}_j=\mathbf{r}_{j-1}+\alpha_{j-1}A\mathbf{p}_{j-1}$ 。根據假設，性質(1)給出 $\mathbf{r}_{j-1}\in\mathcal{K}_j\subseteq\mathcal{K}_{j+1}$ ，由性質(2)可知 $A\mathbf{p}_{j-1}\in A(\mathcal{K}_j)\subseteq\mathcal{K}_{j+1}$ ，因此 $\mathbf{r}_j\in\mathcal{K}_{j+1}$ 。此外，因為 $\mathbf{r}_j\perp\mathcal{P}_j$ ，即 $\mathbf{r}_j\perp\mathcal{K}_j$ ，推論 $\mathbf{r}_j\notin\mathcal{K}_j$ ，否則 $\mathbf{r}_{j}=\mathbf{0}$ 。所以， $\hbox{span}\{\mathbf{r}_0,\mathbf{r}_1,\ldots,\mathbf{r}_j\}=\mathcal{K}_{j+1}$ 。

接著證明性質(2)，考慮 $\mathbf{p}_j=\mathbf{r}_j+\beta_{j-1}\mathbf{p}_{j-1}$ 。性質(1)表明 $\mathbf{r}_j\in\mathcal{K}_{j+1}$ ，性質(2)假設 $\mathbf{p}_{j-1}\in\mathcal{K}_j\subseteq\mathcal{K}_{j+1}$ ，故 $\mathbf{p}_j\in\mathcal{K}_{j+1}$ ，也就證明 $\mathcal{P}_{j+1}=\mathcal{K}_{j+1}$ 。

欲證明性質(3)，考慮 $\mathbf{p}_j^TA\mathbf{p}_i=\mathbf{r}_j^TA\mathbf{p}_i+\beta_{j-1}\mathbf{p}_{j-1}^TA\mathbf{p}_i$ 。若 $i=j-1$ ，由 $\beta_{j-1}$ 的公式可知 $\mathbf{p}_j^TA\mathbf{p}_{j-1}=0$ 。若 $i<j-1$ ， $\mathbf{r}_j^TA\mathbf{p}_{i}=0$ ，因為 $\mathbf{r}_j\perp\mathcal{K}_j$ 且 $A\mathbf{p}_i\in\mathcal{K}_{i+2}\subseteq\mathcal{K}_j$ ； $\mathbf{p}_{j-1}^TA\mathbf{p}_i=0$ ，源自性質(3)的假設。所以， $\mathbf{p}_i^TA\mathbf{p}_j=0$ ， $i<j$ 。

為減少計算量，利用以上性質化簡 $\alpha_j$ 和 $\beta_j$ 。因為 $\mathbf{r}_j\perp\mathcal{P}_j$ ，即 $\mathbf{r}_j\perp\mathbf{p}_{i}$ ， $i<j$ ，可得 $\mathbf{p}_j^T\mathbf{r}_j=(\mathbf{r}_j+\beta_{j-1}\mathbf{p}_{j-1})^T\mathbf{r}_j=\mathbf{r}_j^T\mathbf{r}_j$ ，故

$\displaystyle \alpha_j=\frac{\mathbf{r}_j^T\mathbf{r}_j}{\mathbf{p}_j^TA\mathbf{p}_j}$ 。

因為 $\mathbf{p}_0=\mathbf{r}_0$ ，推知 $\alpha_0=\eta_0$ (見附註^[1])，共軛梯度法的第一次調整量與最陡下降法相同。使用 $A\mathbf{p}_j=-\frac{1}{\alpha_j}(\mathbf{r}_{j+1}-\mathbf{r}_j)$ ， $\beta_j$ 的分子為 $\mathbf{r}_{j+1}^TA\mathbf{p}_j=-\frac{1}{\alpha_j}\mathbf{r}_{j+1}^T(\mathbf{r}_{j+1}-\mathbf{r}_j)=-\frac{1}{\alpha_j}\mathbf{r}_{j+1}^T\mathbf{r}_{j+1}$ ，分母為 $\mathbf{p}_j^TA\mathbf{p}_j=\frac{1}{\alpha_j}\mathbf{r}_j^T\mathbf{r}_j$ ，故

$\displaystyle \beta_j=\frac{\mathbf{r}_{j+1}^T\mathbf{r}_{j+1}}{\mathbf{r}_j^T\mathbf{r}_j}$ 。

我們將共軛梯度法整理於下：對於一 $n\times n$ 階實對稱正定矩陣 $A$ 和非零向量 $\mathbf{b}\in\mathbb{R}^n$ ，設定初始解 $\mathbf{x}_0$ (一般設為零)。

共軛梯度法

$\mathbf{r}_0\leftarrow\mathbf{b}-A\mathbf{x}_0$ ， $\mathbf{p}_0\leftarrow\mathbf{r}_0$
對於 $j=0,1,\ldots,n-1$
$\alpha_j\leftarrow {\mathbf{r}_j^T\mathbf{r}_j}/{\mathbf{p}_j^TA\mathbf{p}_j}$
$\mathbf{x}_{j+1}\leftarrow\mathbf{x}_j+\alpha_j\mathbf{p}_j$
$\mathbf{r}_{j+1}\leftarrow\mathbf{r}_j-\alpha_jA\mathbf{p}_j$
若 $\Vert\mathbf{r}_{j+1}\Vert/\Vert\mathbf{b}\Vert<\epsilon$ ，停止計算。
$\beta_j\leftarrow\mathbf{r}_{j+1}^T\mathbf{r}_{j+1}/\mathbf{r}_j^T\mathbf{r}_j$
$\mathbf{p}_{j+1}\leftarrow \mathbf{r}_{j+1}+\beta_j\mathbf{p}_j$

若不考慮數值計算引進的捨入誤差，共軛梯度法於第 $k$ 次迭代後得到精確解， $k\le n$ 。為控制收斂，設相對殘量 $\Vert\mathbf{r}_{j+1}\Vert/\Vert\mathbf{b}\Vert$ 的容許誤差為 $\epsilon$ 。一般而言，共軛梯度法可快速逼近至精確解，收斂速度取決於 $A$ 的條件數。如果 $A$ 的條件數過大 (即接近病態)，通過前處理 (preconditioning) 可改善收斂性能。具體地說，將 $A\mathbf{x}=\mathbf{b}$ 以等價的 $M^{-1}A\mathbf{x}=M^{-1}\mathbf{b}$ 取代，其中 $M$ 是實對稱正定矩陣使得 $M^{-1}A$ 有較小的條件數^[2]。下圖顯示共軛梯度法 (紅色線段) 與最陡下降法 (綠色線段) 的收斂軌跡，兩個算法產生相同的 $\mathbf{x}_1$ ，但共軛梯度法迭代二次 ( $n=2$ ) 便得到解。

共軛梯度法與最陡下降法的收斂軌跡 From Wikimedia

共軛梯度法具有兩個標誌性質：(1) 朝向精確解的推進方向有共軛性， $\mathbf{p}_i^TA\mathbf{p}_j=0$ ， $i\neq j$ ，(2) 殘差有正交性， $\mathbf{r}_i^T\mathbf{r}_j=0$ ， $i\neq j$ 。最後我們說明何以共軛梯度法被歸於 Krylov 子空間法。共軛梯度法設定近似解 $\mathbf{x}_j=\mathbf{x}_0+\alpha_0\mathbf{p}_0+\alpha_1\mathbf{p}_1+\cdots+\alpha_{j-1}\mathbf{p}_{j-1}$ ，可知 $\mathbf{x}_j-\mathbf{x}_0\in\mathcal{K}_j$ ，且殘差滿足 $\mathbf{r}_j\perp\mathcal{K}_j$ ，此即 Petrov-Galerkin 條件 $\mathbf{b}-A\mathbf{x}_j\perp\mathcal{K}_j$ (見“Krylov 子空間法──線性方程的數值解法 (二)：GMRES 與 FOM”)。共軛梯度法與 FOM 有相同的定義條件，差異在於前者的 Krylov 子空間基底為共軛向量集並以遞迴式計算近似解 $\mathbf{x}_j$ ，後者則以 Lanczos 算法取代 Arnoldi 算法 (因為 $A^T=A$ ) 求得 Krylov 子空間的一組單範正交基底 (orthonormal basis) 之後再計算線性方程解 $\mathbf{x}_j-\mathbf{x}_0$ 。

附註：
[1] 寫出

$\displaystyle f(\mathbf{x}_j+\eta\mathbf{r}_j)=\frac{\eta^2}{2}\mathbf{r}_j^TA\mathbf{r}_j+\eta\mathbf{r}_j^TA\mathbf{x}_j+\frac{1}{2}\mathbf{x}_j^T\mathbf{x}_j-\mathbf{x}_j^T\mathbf{b}-\eta\mathbf{r}_j^T\mathbf{b}$ ，

對 $\eta$ 求導可得

$\displaystyle \frac{df}{d\eta}=\eta\mathbf{r}_j^TA\mathbf{r}_j+\mathbf{r}_j^TA\mathbf{x}_j-\mathbf{r}_j^T\mathbf{b}=\eta\mathbf{r}_j^TA\mathbf{r}_j-\mathbf{r}_j^T\mathbf{r}_j$ 。

設上式為零，解出

$\displaystyle \eta_j=\frac{\mathbf{r}_j^T\mathbf{r}_j}{\mathbf{r}_j^TA\mathbf{r}_j}$ 。

[2] 維基百科：Preconditioner

繼續閱讀：

從滿正交法 (FOM) 推導共軛梯度法 (CG)

2 Responses to Krylov 子空間法──線性方程的數值解法 (三)：共軛梯度法

TZ says:

08/02/2017 at 11:42 pm

感謝, 十分有幫助！

Jing says:

12/26/2018 at 11:31 am

谢谢，也对我帮助很大！

	輕鬆談如何教學二項式定理？… on 牛頓的二項式定理 (上)
	madhouse on 高斯消去法
	WishMobile on 翻轉 LU 分解
	周子傑 on Karush-Kuhn-Tucker (KKT) 條件
	Cloud Huang on 線性泛函與伴隨
	jeff85898 on 利用分塊矩陣證明 det(AB)=(det A)(det…