稀疏矩阵的对称逐次超松弛（SSOR）预处理共轭梯度法
我将为您讲解“稀疏矩阵的对称逐次超松弛（SSOR）预处理共轭梯度法”。这是一个结合了迭代法和预处理技术的高效算法，用于求解稀疏对称正定线性方程组。
题目描述
我们考虑求解大型稀疏对称正定（SPD）线性方程组：
\[Ax = b \]
其中 \(A \in \mathbb{R}^{n \times n}\) 是一个稀疏对称正定矩阵，\(b \in \mathbb{R}^n\) 是已知右端向量，\(x \in \mathbb{R}^n\) 是未知解向量。
共轭梯度法（CG）是求解此类问题的首选迭代法，但当矩阵 \(A\) 的条件数较大时，其收敛速度可能很慢。预处理技术通过求解一个等价但条件更好的方程组来加速收敛。SSOR预处理是一种基于矩阵分裂的简单而有效的预处理子，特别适用于对称正定矩阵。
解题过程循序渐进讲解
步骤1：理解问题背景与核心思想
-
共轭梯度法的局限性：
- CG法的收敛速度取决于矩阵 \(A\) 的条件数 \(\kappa(A)\)：条件数越大，收敛越慢。
- 对于病态问题（\(\kappa(A) \gg 1\)），需要很多迭代步才能达到满意精度。
-
预处理的基本思想：
- 寻找一个预处理矩阵 \(M\)，使得 \(M^{-1}A\) 的条件数远小于 \(A\) 的条件数，且 \(M\) 易于求逆。
- 求解等价方程组：\(M^{-1}Ax = M^{-1}b\)。
- 在CG法中，这体现为修改内积为基于 \(M\) 的内积。
-
SSOR预处理的来源：
- 基于矩阵分裂 \(A = D - L - U\)，其中 \(D\) 为 \(A\) 的对角部分，\(L\) 和 \(U\) 分别为严格下三角和严格上三角部分（且 \(U = L^T\)）。
- SOR（逐次超松弛）迭代法使用参数 \(\omega\) 加速收敛。SSOR是对称版本，即依次执行一次前向SOR和一次后向SOR，保证预处理子对称。
步骤2：推导SSOR预处理矩阵
-
矩阵分裂与SOR迭代矩阵：
- 前向SOR迭代格式：\(x^{(k+1)} = (D - \omega L)^{-1}[\omega U + (1-\omega)D] x^{(k)} + \omega (D - \omega L)^{-1} b\)。
- 后向SOR迭代格式（将 \(A\) 视为 \(D - U - L\)，角色互换）：\(x^{(k+1)} = (D - \omega U)^{-1}[\omega L + (1-\omega)D] x^{(k)} + \omega (D - \omega U)^{-1} b\)。
-
SSOR迭代与预处理矩阵：
- SSOR迭代是前后两次SOR的组合。其迭代矩阵的推导较复杂，但可以直接给出SSOR预处理矩阵 \(M_{\text{SSOR}}\)：
\[ M_{\text{SSOR}} = \frac{1}{\omega(2-\omega)} (D - \omega L) D^{-1} (D - \omega U) \]
其中 $ \omega \in (0, 2) $ 是松弛参数。
- 验证对称性：由于 \(U = L^T\)，有 \((D - \omega U) = (D - \omega L^T)\)，因此 \(M_{\text{SSOR}}\) 对称。
- 实际使用的形式：
- 定义下三角部分 \(C = D - \omega L\)。
- 则 \(M_{\text{SSOR}} = \frac{1}{\omega(2-\omega)} C D^{-1} C^T\)。
- 预处理步骤需要求解 \(M_{\text{SSOR}} z = r\)（其中 \(r\) 是残差），这可以通过两次三角求解完成：
\[ \text{解 } C y = r \cdot \sqrt{\omega(2-\omega)}, \quad \text{再解 } C^T z = D y \]
但通常更直接地，我们利用 $ M_{\text{SSOR}} $ 的分解形式嵌入CG算法。
步骤3：预处理共轭梯度法（PCG）框架
标准PCG算法步骤如下（给定初始猜测 \(x_0\)，预处理子 \(M\)）：
- 计算初始残差 \(r_0 = b - A x_0\)。
- 解 \(M z_0 = r_0\) 得到 \(z_0\)。
- 设初始搜索方向 \(p_0 = z_0\)。
- 对于 \(k = 0, 1, 2, \dots\) 直到收敛：
a. \(\alpha_k = \frac{r_k^T z_k}{p_k^T A p_k}\)。
b. \(x_{k+1} = x_k + \alpha_k p_k\)。
c. \(r_{k+1} = r_k - \alpha_k A p_k\)。
d. 解 \(M z_{k+1} = r_{k+1}\)。
e. \(\beta_k = \frac{r_{k+1}^T z_{k+1}}{r_k^T z_k}\)。
f. \(p_{k+1} = z_{k+1} + \beta_k p_k\)。
步骤4：SSOR-PCG的具体实现
对于SSOR预处理，步骤“解 \(M z = r\)”具体为：
- 前代求解：\((D - \omega L) y = r\)。
- 由于 \(D - \omega L\) 是下三角矩阵，且对角线为 \(D\) 的对角元（因为 \(L\) 对角为零），可以通过前代快速求解（复杂度与矩阵非零元数量成正比）。
- 计算中间量：\(t = D y\)（即对角缩放）。
- 后代求解：\((D - \omega U) z = t\)。
- 由于 \(D - \omega U = (D - \omega L)^T\)，这是上三角系统，可通过后代求解。
注意：实际中常将常数因子 \(\frac{1}{\omega(2-\omega)}\) 吸收到迭代中或忽略（因为CG步骤中的内积会自适应调整），因此常用简化形式 \(M = (D - \omega L) D^{-1} (D - \omega U)\)。
步骤5：参数选择与算法特性
-
松弛参数 \(\omega\)：
- 理论上 \(\omega \in (0, 2)\) 保证 \(M_{\text{SSOR}}\) 正定。
- 最优 \(\omega\) 难以理论确定，通常通过试验在 \([1.0, 1.5]\) 之间选取，或设为1（此时退化为对称Gauss-Seidel预处理）。
- 自适应策略：可根据残差下降情况动态调整。
-
算法优势：
- 保持对称正定性：\(M_{\text{SSOR}}\) 对称正定，可直接用于PCG框架。
- 易于实现：只需矩阵对角线、下三角部分，以及两次三角求解。
- 适合并行：稀疏三角求解可通过level-scheduling等技术并行化。
-
收敛性：
- SSOR预处理通常能显著降低条件数，尤其当 \(A\) 对角占优或具有某种结构时。
- 收敛速度优于无预处理CG，有时甚至接近不完全Cholesky预处理。
步骤6：数值示例（直观演示）
考虑一个简单二维Poisson方程离散得到的5点差分格式，矩阵 \(A\) 为对称正定、稀疏（每行最多5个非零元）。比较：
- 标准CG。
- SSOR-PCG（取 \(\omega = 1.2\)）。
可观察到SSOR-PCG以更少的迭代步达到相同残差精度，尤其当网格细化（条件数增大）时优势更明显。
总结
SSOR预处理共轭梯度法通过对称逐次超松弛分裂构造预处理子，将条件数改善与简单三角求解结合，是求解稀疏对称正定方程组的一种稳健高效方法。其核心在于利用矩阵的三角部分进行快速预处理，同时保持对称性以适配CG框架。该算法在科学计算中广泛应用，尤其适用于中等规模问题或作为更复杂预处理子的基础。