JOVANA
Library Glossary Getting Started Three Levels Fields How it works Mission
Join the mission
All guides

向解迭代逼近:从分裂法到 Krylov

当矩阵大到无法分解时,你只能猜测再改进。本篇从简单的定常迭代——Jacobi 与 Gauss-Seidel——建立直觉,然后揭示每一种现代求解器背后的结构性思想:Krylov 子空间。

迭代的交易

直接法变换 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).
Jacobi 每遍只用旧值;Gauss-Seidel 立即复用新值,通常使迭代次数减半。

飞跃: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 的关键。