4.3 共軛梯度法#

概念說明#

Jacobi、Gauss-Seidel、SOR 都是把 \(A\mathbf{x} = \mathbf{b}\) 寫成不動點迭代，收斂速度被譜半徑鎖死。共軛梯度法（Conjugate Gradient, CG）換了一個世界觀：若 \(A\) 對稱正定（SPD），解線性系統等價於最小化二次函數

\[ \phi(\mathbf{x}) = \frac{1}{2}\mathbf{x}^T A \mathbf{x} - \mathbf{b}^T \mathbf{x}, \qquad \nabla \phi = A\mathbf{x} - \mathbf{b} \]

於是線性代數問題變成最佳化問題（與 5.1 節梯度下降相連）。但 CG 不做樸素梯度下降，而是修正方向：把當前負梯度對「已用過的方向」做正交化，使每一步的方向彼此 A-共軛（ \(\mathbf{p}_i^T A \mathbf{p}_j = 0\) ）。

這個性質帶來近乎魔術的保證：在精確算術下，** \(n\) 步內必達精確解**——因為 \(n\) 個兩兩 A-共軛方向構成基底，逐次沿它們最小化，等價於在 \(n\) 維空間做完全部一維搜尋。更妙的是：若 \(A\) 只有 \(m \ll n\) 個相異特徵值，CG 約 \(m\) 步就收斂；一般情形誤差滿足

\[ \|\mathbf{x}_k - \mathbf{x}^*\|_A \le 2 \left(\frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}\right)^k \|\mathbf{x}_0 - \mathbf{x}^*\|_A, \qquad \kappa = \frac{\lambda_{\max}}{\lambda_{\min}} \]

條件數 \(\kappa\) 越小越快。這催生了預處理（preconditioning）：找一個便宜求逆的 \(M \approx A\) ，解 \(M^{-1}A\mathbf{x} = M^{-1}\mathbf{b}\) 使 \(\kappa\) 大幅下降——預處理共軛梯度（PCG）是科學計算中最重要的演算法之一，入選「二十世紀十大演算法」。

CG 只需矩陣-向量乘與向量運算，完美契合稀疏矩陣；對非對稱矩陣則需換用 GMRES、BiCGSTAB 等親戚。

生活類比#

在橢圓形山谷找最低點。樸素梯度下降像沿著最陡方向走，卻因山谷狹長而不斷「之」字折返；CG 則像聰明的登山者：記住走過的方向並承諾「不再重複修正同一方向」，直線穿過山谷，幾步就到底。

演算法步驟#

初始化 \(\mathbf{x}_0\) ， \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\) ， \(\mathbf{p}_0 = \mathbf{r}_0\) 。
迭代 \(k = 0, 1, \dots\) ：
\(\alpha_k = \dfrac{\mathbf{r}_k^T \mathbf{r}_k}{\mathbf{p}_k^T A \mathbf{p}_k}\)
\(\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{p}_k\)
\(\mathbf{r}_{k+1} = \mathbf{r}_k - \alpha_k A\mathbf{p}_k\)
若 \(\|\mathbf{r}_{k+1}\|\) 足夠小，停止。
\(\beta_k = \dfrac{\mathbf{r}_{k+1}^T \mathbf{r}_{k+1}}{\mathbf{r}_k^T \mathbf{r}_k}\) ， \(\mathbf{p}_{k+1} = \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k\)
每輪只多一次矩陣-向量乘 \(A\mathbf{p}_k\) （結果可與下一步共用）。

Python 範例#

import numpy as np

def cg(A, b, x0=None, tol=1e-10, max_iter=None):
    n = len(b)
    x = np.zeros(n) if x0 is None else x0.copy()
    r = b - A @ x
    p = r.copy()
    rs = r @ r
    hist = [np.sqrt(rs)]
    for k in range(1, (max_iter or n) + 1):
        Ap = A @ p
        alpha = rs / (p @ Ap)
        x += alpha * p
        r -= alpha * Ap
        rs_new = r @ r
        hist.append(np.sqrt(rs_new))
        if np.sqrt(rs_new) < tol:
            break
        beta = rs_new / rs
        p = r + beta * p
        rs = rs_new
    return x, hist

np.random.seed(0)
Q, _ = np.linalg.qr(np.random.randn(6, 6))
A = Q @ np.diag([1, 2, 5, 10, 50, 100]) @ Q.T   # SPD
b = np.random.randn(6)
x, hist = cg(A, b)
print(f"收斂於 {len(hist)-1} 步，殘差 {hist[-1]:.2e}")
print("解驗證誤差:", f"{np.linalg.norm(A @ x - b):.2e}")
print("殘差下降:", [f"{v:.1e}" for v in hist])

複雜度與收斂性分析#

每輪成本 \(O(\text{nnz})\) （一次矩陣-向量乘）加 \(O(n)\) 向量運算，記憶體 \(O(n)\) ——只需存幾個向量，不像直接法需分解。理論上 \(n\) 步精確收斂，實務收斂與 \(\sqrt{\kappa}\) 指數相關： \(\kappa = 10^4\) 時每步誤差約縮 0.98，仍偏慢；配合好的預處理把 \(\kappa\) 壓到 \(10^2\) 以下，幾十步內收斂。CG 的收斂可視為在 Krylov 子空間 \(\text{span}\{\mathbf{r}_0, A\mathbf{r}_0, \dots\}\) 上做最佳逼近，譜分布越聚集收斂越快。使用守則：對稱正定→CG＋預處理；非對稱→GMRES/BiCGSTAB；小而稠密→直接法（LU）往往更划算。

練習#

用 CG 解 4.2 節的模型問題（ \(n = 100\) ），比較與 GS、最佳 SOR 的迭代次數，體會為何 CG 取代了 SOR 成為主流。
構造 \(\kappa = 1\) 的矩陣（對角陣 \(I\) ），驗證 CG 一步收斂；構造 \(\kappa = 10^8\) 的矩陣，觀察收斂停滯。
實作 Jacobi 預處理（ \(M = \text{diag}(A)\) ）的 PCG，對 \(\kappa\) 很大的對角系統比較迭代次數下降幅度。