迭代的交易
直接法变换 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 的关键。