Ax=b稀疏线性方程组的解法介绍
稀疏线性方程组 Ax = b 的解法介绍
在科学计算、工程仿真、有限元分析、计算机图形学等领域，常需求解大规模线性方程组：
A x = b
其中：
- A ∈ ℝⁿˣⁿ 是系数矩阵
- x ∈ ℝⁿ 是未知向量
- b ∈ ℝⁿ 是右端项向量
当矩阵 A 大部分元素为零 时，称为稀疏矩阵（Sparse Matrix）。对这类系统，直接使用稠密矩阵算法（如高斯消元）效率极低，因此需采用专门的稀疏求解方法。
一、稀疏矩阵的特点
- 非零元素数量远小于 n²（例如每行仅有 O(1) 或 O(log n) 个非零元）
- 存储方式：仅存储非零元素及其位置（如 CSR、CSC 格式）
- 计算中避免对零元素进行操作，节省内存与时间
二、稀疏线性方程组的求解方法分类
| 类别 | 方法 | 特点 |
|---|---|---|
| 直接法 | LU 分解（稀疏）、Cholesky 分解 | 精确解，适合中小规模或病态系统 |
| 迭代法 | Jacobi, Gauss-Seidel, SOR | 简单，但收敛慢 |
| Krylov 子空间法 | CG, GMRES, BiCGSTAB | 主流方法，适合大规模稀疏系统 |
三、主要解法介绍
1. 直接法（Direct Methods）
（1）稀疏 LU 分解
若 A 非奇异，可分解为：
A = P L U Q
其中：
- L：单位下三角矩阵（Lower）
- U：上三角矩阵（Upper）
- P、Q：行/列置换矩阵（用于减少“填充”fill-in）
求解过程：
- 分解：A → P L U Q
- 前代：L y = Pᵀ b
- 回代：U z = y
- 重排：x = Q z
✅ 优点：稳定、精确
❌ 缺点：内存消耗大，填充效应严重
🔧 工具：UMFPACK、SuperLU、MKL PARDISO
（2）稀疏 Cholesky 分解（A 对称正定）
若 A 对称正定（SPD），则：
A = P L Lᵀ Pᵀ
求解步骤：
- 分解 A
- 前代：L y = Pᵀ b
- 回代：Lᵀ z = y
- 重排：x = P z
✅ 优点：比 LU 更快更稳定
🔧 工具：CHOLMOD、Eigen::SimplicialLLT
2. 迭代法（Iterative Methods）
（1）Jacobi 方法
迭代格式：
x^(k+1)i = (1/a(ii)) (b_i − Σ_(j≠i) a_(ij) x^(k)_j)
- 逐点更新，依赖前一次迭代值
- 收敛条件：A 严格对角占优
❌ 收敛慢，一般不单独使用
（2）Gauss-Seidel 方法
迭代格式：
x^(k+1)i = (1/a(ii)) (b_i − Σ_(j<i) a_(ij) x^(k+1)j − Σ(j>i) a_(ij) x^(k)_j)
- 使用最新更新值，通常比 Jacobi 快
- 串行性强，难以并行
（3）SOR（Successive Over-Relaxation）
Gauss-Seidel 的加速版本：
x^(k+1)i = (1−ω) x^(k)i + ω (1/a(ii)) (b_i − Σ(j<i) a_(ij) x^(k+1)j − Σ(j>i) a_(ij) x^(k)_j)
其中 ω ∈ (0,2) 为松弛因子：
- ω = 1：退化为 Gauss-Seidel
- ω > 1：超松弛（常用 1.2 ~ 1.8）
⚠️ 需调参，收敛速度依赖 ω
3. Krylov 子空间方法（主流）
Krylov 子空间定义为：
𝒦_k(A,b) = span{ b, A b, A² b, ..., A^(k−1) b }
在该子空间中寻找近似解。
（1）共轭梯度法（Conjugate Gradient, CG）
适用条件：A 对称正定（SPD）
特点：
- 最小化能量范数误差
- 理论上 n 步内收敛（实际因舍入误差需迭代至收敛）
- 每步只需矩阵-向量乘 A·v 和内积计算
迭代流程（简略）：
- r₀ = b − A x₀, p₀ = r₀
- For k = 0,1,... until convergence:
- α_k = (r_kᵀ r_k) / (p_kᵀ A p_k)
- x_(k+1) = x_k + α_k p_k
- r_(k+1) = r_k − α_k A p_k
- β_k = (r_(k+1)ᵀ r_(k+1)) / (r_kᵀ r_k)
- p_(k+1) = r_(k+1) + β_k p_k
✅ 广泛用于有限元、几何处理
🔧 工具：PETSc、SciPy (scipy.sparse.linalg.cg)
（2）GMRES（Generalized Minimal Residual）
适用条件：A 非对称、非奇异
特点：
- 最小化残差 ‖b − A x_k‖
- 使用 Arnoldi 过程构造正交基
- 每 m 步可重启（Restarted GMRES）
✅ 通用性强
❌ 存储开销随迭代步数增长
🔧 工具：PETSc、MATLABgmres
（3）BiCGSTAB（Biconjugate Gradient Stabilized）
适用条件：A 非对称
特点：
- 改进 BiCG 的振荡问题
- 收敛较稳定，每步计算量适中
- 不保证收敛，但实践中表现良好
✅ 常用于 CFD、电磁场仿真
🔧 工具：Eigen、SciPy
四、预条件技术（Preconditioning）
为加速迭代法收敛，引入预条件子 M ≈ A⁻¹，求解：
M⁻¹ A x = M⁻¹ b
理想预条件子：
- M⁻¹ A 的谱分布集中（条件数小）
- M 容易求解（如 M y = r 易解）
常见预条件子：
- 对角预条件（Jacobi）：M = diag(A)
- SSOR：基于 SOR 的预条件
- 不完全 LU（ILU）：稀疏近似 LU 分解
- 代数多重网格（AMG）：对某些 PDE 离散系统极高效
✅ 实践中，ILU + GMRES/BiCGSTAB 是通用组合
五、方法选择建议
| 场景 | 推荐方法 |
|---|---|
| 小到中等规模（n < 10⁵），A SPD | 稀疏 Cholesky |
| 小到中等规模，A 一般 | 稀疏 LU（带填充分解） |
| 大规模，A SPD | CG + ILU/AMG 预条件 |
| 大规模，A 非对称 | GMRES(m) 或 BiCGSTAB + ILU |
| 结构化网格（如PDE） | AMG 预条件 |
| 需多次求解同类系统 | 直接法（分解后可重用） |