迭代的交易
直接法變換 A 直到答案落出;迭代法從一個猜測 x_0 出發並精化它,產生收斂到真實 x 的 x_1, x_2, ...。對於巨型稀疏系統這筆交易很誘人:每次迭代的主要成本是一次矩陣-向量乘積 A v,對稀疏 A 而言僅需 O(nnz)——與非零元個數成正比,而非 n^2。你從不構造 L 或 U,從不承受填充。你只透過乘法接觸 A。
定常迭代:Jacobi 與 Gauss-Seidel
最古老的思想是矩陣分裂:寫 A = M - N,其中 M 易求逆,把 A x = b 化為 M x = N x + b,並迭代 x_{k+1} = M^-1 (N x_k + b)。在 Jacobi 迭代中,M 就是 A 的對角部分,於是每個新分量僅用舊值求解。Gauss-Seidel取 M 為下三角部分,在同一遍掃描內立即複用剛更新的分量——通常以約一半的步數收斂。
Jacobi sweep for A x = b (a_ii on the diagonal):
for i = 1..n:
x_new[i] = ( b[i] - sum_{j != i} a[i][j] * x_old[j] ) / a[i][i]
x_old <- x_new # all components use OLD values
Gauss-Seidel sweep (in place):
for i = 1..n:
x[i] = ( b[i] - sum_{j < i} a[i][j]*x[j] # already updated
- sum_{j > i} a[i][j]*x[j] # still old
) / a[i][i]
Convergence guaranteed if A is strictly diagonally dominant
(|a_ii| > sum_{j != i} |a_ij| for every row).飛躍:Krylov 子空間
定常方法浪費資訊:在第 k 步它只記住 x_k,丟棄了之前的每一個向量。現代洞見是把它們全部保留。從殘差 r_0 = b - A x_0 出發,反覆乘以 A 生成 r_0, A r_0, A^2 r_0, ...,它們的張成就是 Krylov 子空間 K_m = span{r_0, A r_0, ..., A^{m-1} r_0}。
為何這是正確的空間?因為你想要的答案 A^-1 b 可以寫成 A 的多項式作用於 b(Cayley-Hamilton 保證了這一點),而該多項式中 A 的各次冪恰好落在 Krylov 子空間內。於是 K_m 包含著對真解越來越好的逼近,且它僅用你已經負擔得起的稀疏矩陣-向量乘積建構。一個 Krylov 方法按某種準則在 K_m 內挑出最佳逼近——而準則的選擇正是下一篇中區分共軛梯度與 GMRES 的關鍵。