简介这是Yousef Saad所著《Iterative Methods for Sparse Linear Systems》第二版PDF面向数值计算、科学计算与工程仿真领域的研究生、教师和科研人员系统讲解大规模稀疏线性方程组的迭代求解理论与算法。书中先由Jacobi、Gauss-Seidel、SOR等经典方法切入说明从初值出发通过逐步迭代逼近精确解的基本思想再重点分析CG、GMRES、BiCGSTAB等Krylov子空间方法的构造原理、收敛性质与并行计算优势随后围绕块Jacobi、ILU、AMG等预条件器讨论如何降低系数矩阵条件数、加速收敛并提升计算稳定性。正文包含大量算法推导、收敛性分析、数值实验及实际应用示例同时涉及内存受限条件下的迭代实现与预处理选择等实用问题。资源为单个PDF文件压缩包约447MB结构完整便于课程教学、自学阅读和科研检索。目前已有268人学习下载适合需要系统掌握稀疏线性系统迭代解法的高年级本科生、研究生及从事大规模数值计算的工程师。1. 从直接法到迭代法稀疏线性系统为什么要换求解思路大规模稀疏线性系统无处不在。有限元结构分析、油藏数值模拟、流体力学离散最后都落在一个 Axb 上矩阵动辄百万阶非零元占比却不到千分之一。直接法在这里遇到两个硬问题高斯消元过程会产生大量填充内存和计算量随维度快速膨胀稀疏 LU 分解的并行扩展性也远不如想象中理想。迭代法的思路完全不同——不追求一次性精确分解而是从初始解出发用稀疏矩阵-向量乘积反复逼近真解这正是 Iterative Methods for Sparse Linear Systems 这本书解决的核心问题把收敛性分析、算法实现和预条件策略串成一套可以落地的选型方法论。适合数值仿真、科学计算和所有需要处理大规模线性系统的工程师。2. 稀疏矩阵存储与经典迭代格式从Jacobi到SOR2.1 稀疏矩阵存储CSR格式与矩阵-向量乘积迭代方法的每次迭代都围绕稀疏矩阵-向量乘积SpMV展开所以存储格式直接决定性能上限。CSRCompressed Sparse Row是工业界最常用的格式之一它用三个数组描述整个矩阵values 存非零元数值col_indices 存每个非零元对应的列号row_ptr 存每行第一个非零元在 values 中的起始偏移。这样存储总量是 O(nnz n 1)与零元数量无关百万阶矩阵也能轻松装进内存。下面是一个 CSR 格式的 SpMV 实现展示了为什么迭代法能利用稀疏结构def csr_matvec(values, col_indices, row_ptr, x): n len(row_ptr) - 1 y [0.0] * n for i in range(n): start row_ptr[i] end row_ptr[i 1] s 0.0 for j in range(start, end): s values[j] * x[col_indices[j]] y[i] s return y这段代码按行遍历矩阵通过 row_ptr 拿到当前行在 values 中的区间再对区间内的非零元做乘加。col_indices[j] 提供列号确保 x 的对应分量被正确取到。相比稠密矩阵的 n×n 循环这里的内层迭代次数等于该行非零元个数完全避开零元上的无效运算。对五点差分矩阵每行通常只有 5 个非零元一次 SpMV 的耗时几乎就是读取 values 数组的时间这是迭代法能扛住百万规模问题的底层原因。2.2 Jacobi、Gauss-Seidel 与 SOR 的迭代格式Saad 在书中用矩阵分裂统一描述这三类经典方法。把 A 写成 A D - L - UD 是对角部分L 是严格下三角U 是严格上三角。Jacobi 取 M DGauss-Seidel 取 M D - LSOR 则在此基础上引入松弛因子 omega迭代矩阵变为 M(omega) (D - omega L)^(-1)[(1-omega)D omega U]。一个二维泊松方程五点差分离散后的 SOR 求解器实现如下import numpy as np def sor_solver(A, b, omega1.5, tol1e-8, max_iter1000): n len(b) x np.zeros(n) D np.diag(A) for it in range(max_iter): x_new x.copy() for i in range(n): # 用最新分量计算残差等价于Gauss-Seidel的逐次更新 r_i b[i] - A[i, :] x_new x_new[i] x_new[i] omega * r_i / D[i] if np.linalg.norm(b - A x_new) tol: return x_new, it 1 x x_new return x, max_iter这里的残差更新形式与标准 SOR 等价x_i 的新值等于旧值加上 omega 乘以当前残差除以对角元。omega 是松弛因子小于 1 是欠松弛大于 1 是超松弛。对泊松方程五点差分这类模型问题理论上可以解析求出最优 omega但工程中更常见的做法是扫几个候选值对比迭代步数后取最优。需要注意 omega 必须落在 (0, 2) 区间否则方法必然发散这是 SOR 收敛的必要条件。2.3 收敛性判断与迭代矩阵谱半径任意定常迭代格式 x^(k1) M^(-1)N x^k M^(-1)b 的收敛充要条件是迭代矩阵 G M^(-1)N 的谱半径 rho(G) 1。谱半径越小收敛越快。三类方法的迭代矩阵和收敛条件如下方法迭代矩阵 G收敛条件JacobiD^(-1)(LU)rho(G) 1严格对角占优时必然成立Gauss-Seidel(D-L)^(-1)Urho(G) 1A 对称正定时成立SOR(D-omega L)^(-1)[(1-omega)D omega U]0 omega 2A 对称正定时成立实际中直接算谱半径对大矩阵不现实工程上更常用的是观察残差范数的下降曲线。如果残差每步只下降一个固定比例说明谱半径接近 1需要考虑预条件或换方法。对五点差分离散的泊松方程Jacobi 迭代的谱半径约为 cos(pi/(n1))网格越细越接近 1。这就是经典迭代法在大规模问题上的效率瓶颈也是 Krylov 子空间方法登场的直接原因。3. Krylov子空间方法CG与GMRES的算法机制3.1 Krylov子空间的构造与Arnoldi过程Krylov 子空间 K_m(A, v) span{v, Av, A²v, …, A^(m-1)v} 的核心想法是用矩阵反复作用一个初始向量张成一个低维子空间再在其中寻找近似解。Arnoldi 过程是这个框架的基石它通过 Gram-Schmidt 正交化生成该子空间的一组正交基同时把原矩阵压缩成一个小的上 Hessenberg 矩阵。import numpy as np def arnoldi(A, v0, m): n len(v0) v v0 / np.linalg.norm(v0) V np.zeros((n, m 1)) H np.zeros((m 1, m)) V[:, 0] v for j in range(m): w A V[:, j] # 每步一次SpMV主导开销 for i in range(j 1): H[i, j] np.dot(w, V[:, i]) w w - H[i, j] * V[:, i] # 修正正交性 H[j 1, j] np.linalg.norm(w) if H[j 1, j] 1e-14: break # 子空间已不变可提前终止 V[:, j 1] w / H[j 1, j] return V, H返回的 V 和 H 满足关系 AV_m V_(m1)H_m。H 只在主对角线和次对角线下方有非零元原本 n 维的问题被压缩成一个 (m1)×m 的小规模最小二乘问题。H[j1, j] 接近 0 时说明 Krylov 子空间已经不变此时方法会在有限步内精确收敛这是理解 GMRES 和 CG 有限步终止性质的关键。3.2 CG算法对称正定矩阵的Lanczos路线当 A 对称正定时Arnoldi 过程退化为 Lanczos 过程H 变成对称三对角矩阵只需要保存三条对角线和两组基向量内存开销从 O(n·m) 降到 O(n)。CG 正是在这个框架下通过极小化 A-范数误差得到的算法。def conjugate_gradient(A, b, x0, tol1e-8, max_iter1000): x x0.copy() r b - A x p r.copy() rs_old np.dot(r, r) for it in range(max_iter): Ap A p alpha rs_old / np.dot(p, Ap) # 步长极小化A-范数误差 x x alpha * p r r - alpha * Ap rs_new np.dot(r, r) if np.sqrt(rs_new) tol: return x, it 1 p r (rs_new / rs_old) * p # 更新共轭搜索方向 rs_old rs_new return x, max_iteralpha 的表达式来自一维极小化p 序列在 A-内积下相互共轭。收敛速度由条件数 kappa lambda_max / lambda_min 决定理论误差上界为 (sqrt(kappa)-1)/(sqrt(kappa)1) 的平方。kappa 越大收敛越慢这就是为什么 CG 必须配合预条件使用。另外 CG 对舍入误差比较敏感残差可能短暂上升但只要矩阵严格对称正定整体迭代不会发散。若矩阵非对称直接用 CG 会产生完全不可信的结果。3.3 GMRES非对称问题的极小残差框架对非对称矩阵短递推关系不再成立GMRES 的做法是直接在 m 维 Krylov 子空间中极小化残差范数。借助 Arnoldi 分解原问题转化为一个 (m1)×m 的最小二乘问题只需在 Arnoldi 循环结束后做一次小规模求解。def gmres(A, b, x0, m30, tol1e-8, max_outer100): x x0.copy() for outer in range(max_outer): r b - A x beta np.linalg.norm(r) V np.zeros((A.shape[0], m 1)) H np.zeros((m 1, m)) V[:, 0] r / beta for j in range(m): w A V[:, j] for i in range(j 1): H[i, j] np.dot(w, V[:, i]) w w - H[i, j] * V[:, i] H[j 1, j] np.linalg.norm(w) if H[j 1, j] 1e-14: break V[:, j 1] w / H[j 1, j] k j 1 rhs np.zeros(k 1) rhs[0] beta y, _, _, _ np.linalg.lstsq(H[:k 1, :k], rhs, rcondNone) x x0 V[:, :k] y if np.linalg.norm(b - A x) tol: return x, outer * m k x0 x.copy() return x, max_outer * mGMRES 每一步都要保存一组正交基 V到第 m 步需要 n×m 的存储这是它的主要空间开销。m 是内层迭代维数即重启前的步数工程上通常取 20 到 50。重启会丢弃已构建的子空间信息导致收敛曲线出现平台甚至停滞。对特征值分布复杂的矩阵单纯调大 m 不一定有效更可靠的做法是配合预条件。3.4 FOM与GMRES的关系及实际选型书中 6.5.7 节详细讨论了 FOM 与 GMRES 的关系两者在相同的 Krylov 子空间上工作FOM 使用 Ritz-Galerkin 条件要求残差与子空间正交GMRES 则直接极小化残差范数。GMRES 的残差单调不增数值稳定性更强而 FOM 的残差可能波动。实际工程中对收敛行为要求严格时优先选 GMRES资源紧张且矩阵接近对称时可以考虑 CG 或 BiCGSTAB。BiCGSTAB 通过双正交化构造每步只需两次 SpMV 和有限个短递推向量存储占用 O(n)是非对称问题上存储与收敛之间的折中方案。方法适用矩阵存储每步SpMV数残差行为CG对称正定O(n)1单调下降数值上可微扰GMRES一般非奇异O(n·m)1单调不增BiCGSTAB一般非奇异O(n)2可波动需配合预条件4. 预条件器的设计与实现ILU与相关问题4.1 预条件为什么有效从条件数到谱聚集预条件的本质是构造一个与 A 相近的矩阵 M把原方程转化为 M^(-1)Ax M^(-1)b。理想情况下 M^(-1)A 的特征值聚集在 1 附近而不是散布在很宽的区间。因为 Krylov 方法的收敛速度由特征值分布决定条件数小只是一方面特征值簇的紧凑程度往往更关键。预条件器设计的目标就是让矩阵的谱分布更利于 Krylov 迭代。书中 10.2 节从最简单的 Jacobi 预条件讲起取 M D对对角占优矩阵有效但对泊松方程这类问题帮助有限。SSOR 预条件利用 SOR 分裂构造对称形式M (1/(2-omega))(D/omega L)(D/omega)^(-1)(D/omega U)^T对对称正定矩阵可以保持预条件后的对称性从而与 CG 无缝配合。4.2 ILU分解从ILU(0)到ILUTILU 的核心是做一个受限的 LU 分解。ILU(0) 只允许 L 和 U 在与原矩阵相同的非零图案上存在非零元其余位置直接填零。它对结构规则的问题效果不错但对强非对称或病态问题往往不够。更通用的是 ILUT引入两个阈值参数drop_tol 控制小元素的丢弃阈值fill_factor 控制每行最多保留的填充元素数量。from scipy.sparse.linalg import spilu, LinearOperator def setup_ilut_preconditioner(A, fill_factor10, drop_tol1e-4): A_csc A.tocsc() ilu spilu(A_csc, fill_factorfill_factor, drop_toldrop_tol) M LinearOperator( A_csc.shape, matveclambda x: ilu.solve(x), rmatveclambda x: ilu.solve(x, T), ) return Mspilu 返回一个不完全 LU 分解对象matvec 调用前代和回代完成 M^(-1)x。fill_factor10 表示每行最多保留原非零元数量 10 倍的填充drop_tol1e-4 表示绝对值小于该阈值的填充元素直接丢弃。这两个参数直接影响预条件效果和构建成本。一般经验是先固定 drop_tol1e-4从 fill_factor5 开始逐步加大观察 GMRES 迭代次数与总耗时的变化曲线找到拐点即可。fill_factor 调得过大预条件器构建时间和内存会爆炸但迭代次数并不会无限下降。4.3 左预条件、右预条件与灵活GMRES预条件可以作用在系统左侧或右侧。左预条件求解 M^(-1)Ax M^(-1)b但 GMRES 此时极小化的是预条件后残差的范数不是原始残差这可能导致收敛判据失真。右预条件求解 AM^(-1)u b再令 x M^(-1)u此时原始残差范数可以直接监控工程上更推荐这种做法。书中 9.4 节的 FGMRES 进一步放宽限制允许每一步使用不同的预条件器 M_j。这对多网格、非线性问题嵌套迭代等场景很有价值。使用 FGMRES 时需要额外保存每一步的 z_j M_j^(-1)v_j存储开销比标准 GMRES 多一组向量但换来的是预条件选择的灵活性。4.4 预条件器选型建议面对一个实际的稀疏系统预条件器的选择往往决定求解成败。下面是不同预条件器的适用场景对比预条件器适用场景内存开销构建代价并行友好度Jacobi对角占优矩阵快速验证O(n)极低高SSOR对称正定矩阵配合CGO(n)低一般ILU(0)结构规则问题O(nnz)中低ILUT一般稀疏矩阵效果好O(nnz·fill)中高低近似逆需要并行求解的场景O(nnz)中高我的做法是先跑一次无预条件 GMRES 建立基线如果上百步还不收敛直接上 ILUT 并调 fill_factor。二维 PDE 离散矩阵通常 ILU(0) 就够三维问题或强各向异性介质建议用 AMG 或 ILUT。千万自由度以上还要关注预条件器构建时的内存峰值ILU 的填充往往比预期多得多必要时用近似逆预条件器换取更好的并行扩展性。5. 残差曲线诊断复现收敛实验时的一个关键技巧在复现书中数值实验时有一个很实用的操作用残差曲线判断预条件器是否真正起了作用而不是只看最终迭代次数。先构造一个二维泊松方程五点差分矩阵分别用无预条件 GMRES 和 ILU 预条件 GMRES 求解记录每一步的残差范数。import numpy as np from scipy.sparse import diags from scipy.sparse.linalg import spilu, LinearOperator, gmres def poisson2d(n): N n * n diagonals [ 4.0 * np.ones(N), -np.ones(N - 1), -np.ones(N - 1), -np.ones(N - n), -np.ones(N - n), ] offsets [0, 1, -1, n, -n] A diags(diagonals, offsets, formatcsr) for k in range(1, n): A[k * n - 1, k * n] 0.0 A[k * n, k * n - 1] 0.0 return A.tocsr() n 50 A poisson2d(n) b np.ones(n * n) res_plain [] def cb_plain(xk): res_plain.append(np.linalg.norm(b - A xk)) gmres(A, b, tol1e-10, restart30, maxiter200, callbackcb_plain) ilu spilu(A.tocsc(), drop_tol1e-4, fill_factor10) M LinearOperator(A.shape, matveclambda x: ilu.solve(x), dtypeA.dtype) res_pre [] def cb_pre(xk): res_pre.append(np.linalg.norm(b - A xk)) gmres(A, b, tol1e-10, restart30, maxiter200, MM, callbackcb_pre)这里的关键技巧是用 callback 在每个迭代步获取当前解 xk再重新计算残差范数。scipy 的 gmres 不同版本对 callback 的传参有差异有的传残差有的传解向量统一按解向量处理再自行计算残差能保证实验代码在多个版本下一致运行。poisson2d 构造矩阵后把每个 n 行块末尾指向下一行块开头的两个非零元清零防止 1 和 -1 对角线产生跨行连接。对比 res_plain 和 res_pre 两条曲线时重点看两件事一是下降到 1e-10 所需的总步数无预条件 GMRES 通常要上百步ILU 预条件后大约十几步二是曲线斜率是否稳定。如果预条件后的残差曲线在某个阶段突然变平说明 ILU 对部分特征值没有起到聚集作用此时应提高 fill_factor 或改用其他预条件器。评测预条件器时不能只看迭代次数必须把 ILU 构建时间和每次迭代的 SpMV 成本一起算进墙钟时间某些场景下 SSOR 预条件 CG 的总耗时反而低于 ILU 预条件 GMRES。把 restart 和 fill_factor 做成二维扫描观察残差曲线下包络的形态是定位这类问题最直接的手段。本文还有配套的精品资源点击获取