數值分析系列（三）：大型線性系統的迭代法

封存文章

本文由舊站技術筆記整理而成，保留原有主題，但已在 2026 年重寫若干公式、程式碼與複雜度說明。它適合作為概念導論，不應取代針對特定偏微分方程或矩陣結構的求解器分析。

大型科學計算很少真的把矩陣當成一張密集表格。以三維 Poisson 方程在 $100\times100\times100$ 網格上的有限差分離散為例，未知數有一百萬個；若把矩陣密集儲存，雙精度資料約需 8 TB，密集 Gaussian elimination 的運算量亦達 $O(n^3)$ 。然而 Poisson 矩陣每列通常只有約七個非零元素，所以真正的問題不是「所有直接法都不可行」，而是要利用稀疏性、對稱性、正定性與網格結構。稀疏直接法的成本取決於填充，迭代法則主要以矩陣向量乘法和預條件器來換取逐步收斂。

這篇文章集中回答三個問題：

固定點迭代何時收斂？
Krylov 子空間法為何比單純矩陣分裂有效？
預條件如何改變「難解」的線性系統？

1. 矩陣分裂與固定點迭代

考慮

Ax=b,\qquad A=M-N,

其中 $M$ 容易求解。於是

Mx=Nx+b

可寫成迭代

x_{k+1}=M^{-1}Nx_k+M^{-1}b =Tx_k+c.

實作時通常不會顯式計算 $M^{-1}$ ，而是每一步解 $Mz=Nx_k+b$ 。若 $x^\star$ 為精確解，誤差滿足

e_{k+1}=Te_k,\qquad e_k=T^ke_0.

因此對所有初始值都收斂的充要條件是

\rho(T)<1,

其中 $\rho(T)$ 是譜半徑。這個條件比「某一個矩陣範數小於 1」更基本：若能找到一致矩陣範數使 $\lVert T\rVert<1$ ，它是方便的充分條件，但不是必要條件。

2. Jacobi 與 Gauss–Seidel

採用符號

A=D+L+U,

其中 $D$ 是對角部分， $L$ 、 $U$ 分別是嚴格下、上三角部分。Jacobi 法為

x_{k+1}=D^{-1}\bigl(b-(L+U)x_k\bigr),

迭代矩陣是

T_J=-D^{-1}(L+U).

一個簡潔的 NumPy 版本如下：

import numpy as np

def jacobi(A, b, x0=None, rtol=1e-8, maxiter=1_000):
    A = np.asarray(A, dtype=float)
    b = np.asarray(b, dtype=float)
    x = np.zeros_like(b) if x0 is None else np.array(x0, dtype=float)
    d = np.diag(A)

    if np.any(d == 0):
        raise ValueError("Jacobi 法要求所有對角元素非零")

    R = A - np.diag(d)
    bnorm = max(np.linalg.norm(b), 1.0)

    history = []
    for _ in range(maxiter):
        x_new = (b - R @ x) / d
        residual = np.linalg.norm(b - A @ x_new) / bnorm
        history.append(residual)
        if residual <= rtol:
            return x_new, history
        x = x_new

    return x, history

嚴格對角優勢是 Jacobi 收斂的常用充分條件，但

\rho(T_J)

必須由 $T_J$ 的特徵值定義；它不能一般地寫成逐列係數和的某個最大值。而且停止準則宜檢查相對殘差

\frac{\lVert b-Ax_k\rVert}{\max(\lVert b\rVert,1)},

因為相鄰兩次迭代很接近，不一定代表 $Ax_k\approx b$ 。

Gauss–Seidel 法會即時使用本輪已更新的分量：

(D+L)x_{k+1}=b-Ux_k.

它在不少橢圓型問題中較 Jacobi 快，但資料相依性令平行化較困難。兩者今天仍很重要，尤其作為 multigrid smoother 或較複雜預條件器的組件，而非大型問題的唯一求解器。

3. 共軛梯度法

當 $A$ 是實對稱正定矩陣（SPD）時， $Ax=b$ 等價於最小化

\phi(x)=\frac12x^\mathsf{T}Ax-b^\mathsf{T}x.

因為

\nabla\phi(x)=Ax-b=-r(x),

而 Hessian 正是 $A$ 。共軛梯度法（Conjugate Gradient, CG）不只沿最陡下降方向前進，而是建立互相 $A$ -共軛的搜尋方向：

p_i^\mathsf{T}Ap_j=0,\qquad i\ne j.

標準遞推是

\begin{aligned} \alpha_k&=\frac{r_k^\mathsf{T}r_k}{p_k^\mathsf{T}Ap_k},\\ x_{k+1}&=x_k+\alpha_kp_k,\\ r_{k+1}&=r_k-\alpha_kAp_k,\\ \beta_k&=\frac{r_{k+1}^\mathsf{T}r_{k+1}}{r_k^\mathsf{T}r_k},\\ p_{k+1}&=r_{k+1}+\beta_kp_k. \end{aligned}

要分清兩種正交性：在精確算術中，殘差 $r_i$ 與 $r_j$ 是 Euclidean 正交；搜尋方向 $p_i$ 與 $p_j$ 才是 $A$ -共軛。原筆記曾把殘差稱為「 $A$ -orthogonal」，這裡已更正。

def conjugate_gradient(A, b, x0=None, rtol=1e-10, maxiter=None):
    b = np.asarray(b, dtype=float)
    x = np.zeros_like(b) if x0 is None else np.array(x0, dtype=float)
    maxiter = len(b) if maxiter is None else maxiter

    r = b - A @ x
    p = r.copy()
    rr = r @ r
    bnorm = max(np.linalg.norm(b), 1.0)
    history = [np.sqrt(rr) / bnorm]

    for _ in range(maxiter):
        Ap = A @ p
        curvature = p @ Ap
        if curvature <= 0:
            raise ValueError("A 可能不是對稱正定矩陣")

        alpha = rr / curvature
        x = x + alpha * p
        r_new = r - alpha * Ap
        rr_new = r_new @ r_new
        history.append(np.sqrt(rr_new) / bnorm)

        if history[-1] <= rtol:
            return x, history

        p = r_new + (rr_new / rr) * p
        r, rr = r_new, rr_new

    return x, history

這裡只記錄可直接計算的殘差。若不知道 $x^\star$ ，就不能由 $\sqrt{x^\mathsf{T}Ax-2b^\mathsf{T}x}$ 得到誤差的 energy norm；該式還欠一個與 $x$ 無關但必需的常數 $b^\mathsf{T}A^{-1}b$ ，而根號內甚至可能為負。

CG 在 $x_0+\mathcal K_k(A,r_0)$ 中最小化 $A$ -範數誤差：

\lVert x_k-x^\star\rVert_A \le 2\left( \frac{\sqrt{\kappa_2(A)}-1} {\sqrt{\kappa_2(A)}+1} \right)^k \lVert x_0-x^\star\rVert_A.

在精確算術中，CG 最多 $n$ 步便會終止；浮點運算會破壞正交性，因此實際表現主要取決於特徵值分佈、條件數與預條件器。對稀疏矩陣，每一步通常是 $O(\operatorname{nnz}(A))$ ，而不是籠統的 $O(n^2)$ 。

4. GMRES 與 Arnoldi 過程

對一般非對稱矩陣，GMRES（Generalized Minimal Residual）在 Krylov 子空間

\mathcal K_k(A,r_0) =\operatorname{span}\{r_0,Ar_0,\ldots,A^{k-1}r_0\}

中尋找使二範數殘差最小的近似：

x_k=x_0+V_ky_k,\qquad y_k=\arg\min_y\lVert\beta e_1-\bar H_ky\rVert_2.

Arnoldi 過程建立正交基 $V_{k+1}$ 與上 Hessenberg 矩陣 $\bar H_k$ ，使

AV_k=V_{k+1}\bar H_k.

以下是未重啟的教學版本；大型實務問題通常應直接使用 SciPy，並加入 restart 與 preconditioner：

def gmres(A, b, x0=None, maxiter=None, rtol=1e-8):
    b = np.asarray(b, dtype=float)
    x0 = np.zeros_like(b) if x0 is None else np.array(x0, dtype=float)
    n = len(b)
    m = n if maxiter is None else min(maxiter, n)

    r0 = b - A @ x0
    beta = np.linalg.norm(r0)
    if beta == 0:
        return x0, [0.0]

    V = np.zeros((n, m + 1))
    H = np.zeros((m + 1, m))
    V[:, 0] = r0 / beta
    rhs = np.zeros(m + 1)
    rhs[0] = beta
    history = []

    for j in range(m):
        w = A @ V[:, j]
        for i in range(j + 1):
            H[i, j] = V[:, i] @ w
            w -= H[i, j] * V[:, i]

        H[j + 1, j] = np.linalg.norm(w)
        if H[j + 1, j] > 0:
            V[:, j + 1] = w / H[j + 1, j]

        y, *_ = np.linalg.lstsq(
            H[:j + 2, :j + 1],
            rhs[:j + 2],
            rcond=None,
        )
        x = x0 + V[:, :j + 1] @ y
        relres = np.linalg.norm(b - A @ x) / max(np.linalg.norm(b), 1.0)
        history.append(relres)

        if relres <= rtol or H[j + 1, j] == 0:
            return x, history

    return x, history

GMRES 的多項式描述要包括約束 $p(0)=1$ ：

\lVert r_k\rVert_2 = \min_{\substack{p\in\Pi_k\\p(0)=1}} \lVert p(A)r_0\rVert_2.

若 $A$ 為 normal matrix，便有譜上的上界

\lVert r_k\rVert_2 \le \min_{\substack{p\in\Pi_k\\p(0)=1}} \max_{\lambda\in\sigma(A)}|p(\lambda)| \lVert r_0\rVert_2.

對非奇異 $n\times n$ 矩陣，完整 GMRES 在精確算術下最多 $n$ 步收斂，更精確地說是由相對於 $r_0$ 的 minimal polynomial 次數控制。但儲存所有 Arnoldi 向量的成本隨步數增加；重啟 GMRES( $m$ ) 雖節省記憶體，卻可能停滯，所以「一般矩陣最多 $n$ 步」不能直接套用到重啟版本。

5. 預條件不是附加選項

我們不一定直接解 $Ax=b$ ，而會選一個容易作用、又近似 $A$ 的矩陣 $P$ ：

P^{-1}Ax=P^{-1}b.

理想的 $P^{-1}A$ 具有較集中的特徵值、較小的有效條件數，而且套用 $P^{-1}$ 的成本遠低於解原系統。常見選擇包括：

Jacobi 或 block-Jacobi；
incomplete LU / incomplete Cholesky；
sparse approximate inverse；
geometric 或 algebraic multigrid；
由 PDE、網格或 saddle-point 結構推導的 block preconditioner。

Multigrid 的 V-cycle 需要明確的平滑器、restriction、prolongation、各層 operator 與 coarse-grid solve；幾行沒有狀態更新的偽代碼並不是可執行預條件器。在 Poisson 類橢圓問題中，multigrid 的價值在於同時處理高頻與低頻誤差，理想情況下可做到接近線性的工作量。

6. 方法如何選擇

方法	典型矩陣	每步主要成本	記憶體	重要限制
Jacobi	對角優勢或適合作 smoother	$O(\operatorname{nnz}(A))$	$O(n)$	通常慢；需 $\rho(T_J)<1$
Gauss–Seidel	結構良好的稀疏系統	$O(\operatorname{nnz}(A))$	$O(n)$	平行化較難
CG	對稱正定	$O(\operatorname{nnz}(A))$	$O(n)$	非 SPD 時不可直接使用
GMRES	一般非對稱	matvec 加正交化	隨 Krylov 維度增加	重啟可能停滯

這張表不是「最佳方法排行榜」。Jacobi 並非對角優勢矩陣的普遍最優法； GMRES 的「最優」只指在當前 affine Krylov 子空間內最小化二範數殘差； CG 的最小化性亦建立在 SPD 假設上。

典型應用

計算流體力學： pressure-correction、projection method 與 Navier–Stokes 離散會產生大型稀疏或 saddle-point 系統。
結構力學： Newton 線性化、接觸問題與 modal analysis 反覆需要解線性系統。
擴散與熱傳： 橢圓／拋物型 PDE 的離散特別適合 multigrid 與 CG 類方法。
資料同化與機器學習： Hessian-vector product 與 matrix-free Krylov 方法可避免顯式形成大型矩陣。

7. 小結

大型線性系統的核心並不是把密集直接法與迭代法作簡單對比，而是先問：

矩陣是否稀疏、對稱、正定或具有 block 結構？
能否只提供 $v\mapsto Av$ 的 matrix-free 操作？
應以誤差、殘差，還是物理量作停止準則？
哪一個預條件器能把結構轉化為更快收斂？

矩陣分裂提供固定點觀點；CG 與 GMRES 把問題轉化為 Krylov 子空間中的多項式近似；預條件則把抽象收斂界帶回實際計算。三者合起來，才是現代迭代線性代數的基本框架。

參考資料

Kincaid, D., & Cheney, W. (2009). Numerical Analysis: Mathematics of Scientific Computing (3rd ed.). AMS.
Saad, Y. (2003). Iterative Methods for Sparse Linear Systems (2nd ed.). SIAM.
Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press.
Trefethen, L. N., & Bau, D. III. (1997). Numerical Linear Algebra. SIAM.
Hackbusch, W. (2016). Iterative Solution of Large Sparse Systems of Equations (2nd ed.). Springer.
Greenbaum, A. (1997). Iterative Methods for Solving Linear Systems. SIAM.
Davis, T. A. (2006). Direct Methods for Sparse Linear Systems. SIAM.
Kelley, C. T. (1995). Iterative Methods for Linear and Nonlinear Equations. SIAM.
Saad, Y., & van der Vorst, H. A. (2000). Iterative solution of linear systems in the 20th century. Journal of Computational and Applied Mathematics, 123(1–2), 1–33.