基于Chebyshev加速的对称逐次超松弛迭代法（Chebyshev-SSOR）求解对称正定线性方程组
算法描述
本题目介绍Chebyshev加速的对称逐次超松弛迭代法（Chebyshev-SSOR），用于高效求解大规模稀疏对称正定线性方程组 \(A\mathbf{x} = \mathbf{b}\)，其中 \(A \in \mathbb{R}^{n \times n}\) 对称正定， \(\mathbf{x}, \mathbf{b} \in \mathbb{R}^{n}\)。该方法结合了对称逐次超松弛（SSOR）作为内迭代（预处理/平滑器），以及基于切比雪夫多项式的最优参数选择的外迭代（加速），以显著提高经典SSOR迭代法的收敛速度。与共轭梯度法（CG）不同，该方法不依赖于正交性条件，在某些并行计算场景或当矩阵条件数可以较好估计时具有优势。
我们的目标是求解 \(\mathbf{x}_* = A^{-1} \mathbf{b}\)。记初始猜测为 \(\mathbf{x}^{(0)}\)，迭代格式将生成序列 \(\{\mathbf{x}^{(k)}\}\) 以逼近 \(\mathbf{x}_*\)。
解题过程
步骤1：回顾对称逐次超松弛（SSOR）迭代
SSOR是求解 \(A\mathbf{x} = \mathbf{b}\) 的一种经典定常迭代法，适用于对称矩阵。它由两个半步组成：
首先，将矩阵 \(A\) 分解为 \(A = D - L - L^T\)，其中 \(D\) 是对角部分（假设 \(D\) 的对角元非零）， \(-L\) 是严格下三角部分。
给定第 \(k\) 次迭代的近似解 \(\mathbf{x}^{(k)}\)，SSOR迭代计算 \(\mathbf{x}^{(k+1)}\) 如下：
- 前向SOR步（Forward SOR）：对 \(i = 1, 2, ..., n\) 顺序计算：
\[ x_i^{(k+\frac{1}{2})} = (1-\omega)x_i^{(k)} + \frac{\omega}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+\frac{1}{2})} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k)} \right) \]
这里 $ \omega \in (0, 2) $ 是松弛参数。此步用到了刚更新出的分量 $ x_j^{(k+\frac{1}{2})} $（当 $ j < i $）。
- 后向SOR步（Backward SOR）：对 \(i = n, n-1, ..., 1\) 逆序计算：
\[ x_i^{(k+1)} = (1-\omega)x_i^{(k+\frac{1}{2})} + \frac{\omega}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+\frac{1}{2})} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k+1)} \right) \]
此步用到了刚更新出的分量 $ x_j^{(k+1)} $（当 $ j > i $）。
这两个步骤合起来称为一次SSOR迭代。将其写成矩阵形式更为简洁。定义 \(M_{SSOR}\) 为SSOR迭代的迭代矩阵，其计算格式可以等价地写为：
\[\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + M_{SSOR}^{-1} (\mathbf{b} - A\mathbf{x}^{(k)}) \]
其中，
\[M_{SSOR} = \frac{1}{\omega(2-\omega)} (D - \omega L) D^{-1} (D - \omega L^T) \]
这里 \(M_{SSOR}\) 对称正定（当 \(A\) 对称正定且 \(0 < \omega < 2\)）。因此，SSOR迭代也可视为用 \(M_{SSOR}\) 作为预处理子的不动点迭代。
步骤2：引入Chebyshev加速（多项式加速）
定常迭代法（如SSOR）的收敛速度由其迭代矩阵的谱半径决定。Chebyshev加速是一种非定常迭代法，它通过组合前几次迭代的信息，以在每次迭代中最小化某种误差范数，从而达到比原始定常迭代更快的收敛速度。
核心思想：将第 \(k\) 次迭代的解表示为初始残差的多项式形式，并选择多项式系数使得误差在某种意义下最小化。
设我们有一个基本的一阶定常迭代格式（如SSOR）：
\[\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + M^{-1} (\mathbf{b} - A\mathbf{x}^{(k)}) = \mathbf{x}^{(k)} + M^{-1} \mathbf{r}^{(k)} \]
其中 \(M\) 是迭代矩阵（对于SSOR， \(M = M_{SSOR}\)）， \(\mathbf{r}^{(k)} = \mathbf{b} - A\mathbf{x}^{(k)}\) 是残差。
这个格式可以改写为误差方程。令 \(\mathbf{e}^{(k)} = \mathbf{x}_* - \mathbf{x}^{(k)}\) 为误差向量。则
\[\mathbf{e}^{(k+1)} = (I - M^{-1}A) \mathbf{e}^{(k)} = G \mathbf{e}^{(k)} \]
其中 \(G = I - M^{-1}A\) 是误差传播矩阵（iteration matrix）。经过 \(k\) 步迭代后，误差满足：
\[\mathbf{e}^{(k)} = G^k \mathbf{e}^{(0)} \]
Chebyshev加速的目标：不直接使用 \(G^k\)，而是寻找一个 \(k\) 次多项式 \(P_k(\lambda)\) 满足 \(P_k(0) = 1\)，并用 \(P_k(G) \mathbf{e}^{(0)}\) 来近似误差，使得在某种范数下 \(\| P_k(G) \mathbf{e}^{(0)} \|\) 尽可能小。由于我们不知道 \(\mathbf{e}^{(0)}\)，这通常通过构造一个关于残差的递推格式来实现。
在实际应用中，假设矩阵 \(M^{-1}A\) 的特征值位于正实数区间 \([ \alpha, \beta ]\) 内，其中 \(0 < \alpha \le \beta\)。对于SSOR预处理，当 \(A\) 对称正定时， \(M_{SSOR}^{-1}A\) 的特征值是正的，且区间 \([\alpha, \beta]\) 可以估计（有时取 \(\alpha = \lambda_{\min}(M^{-1}A)\), \(\beta = \lambda_{\max}(M^{-1}A)\)）。
步骤3：构建Chebyshev加速迭代格式
Chebyshev加速迭代的标准三项递推格式如下：
给定初始猜测 \(\mathbf{x}^{(0)}\)，计算 \(\mathbf{x}^{(1)}\) 通常先进行一次基本迭代（如SSOR）：
\[\mathbf{x}^{(1)} = \mathbf{x}^{(0)} + M^{-1} \mathbf{r}^{(0)}, \quad \mathbf{r}^{(0)} = \mathbf{b} - A\mathbf{x}^{(0)} \]
然后，对于 \(k \ge 1\)，迭代格式为：
\[\mathbf{x}^{(k+1)} = \rho_k \left[ \mathbf{x}^{(k)} + M^{-1} \mathbf{r}^{(k)} - \mathbf{x}^{(k-1)} \right] + \mathbf{x}^{(k-1)} \]
其中系数 \(\rho_k\) 由以下公式给出：
\[\rho_k = \left( 1 - \frac{\mu}{2} C_{k-1}(\mu) / C_{k}(\mu) \right)^{-1} \]
这里 \(C_k(t)\) 是 \(k\) 次第一类切比雪夫多项式， \(\mu = \frac{\beta - \alpha}{\beta + \alpha}\)。
更常用的等价形式（避免直接计算切比雪夫多项式比值）是利用三项递推关系计算系数。定义：
\[\delta = \frac{\beta - \alpha}{\beta + \alpha}, \quad \gamma = \frac{2}{\beta + \alpha} \]
则迭代格式可写为：
- 初始化：
\[ \mathbf{x}^{(1)} = \mathbf{x}^{(0)} + \gamma M^{-1} \mathbf{r}^{(0)} \]
注意这里使用了初始缩放因子 \(\gamma\)，它对应于将特征值区间映射到 \([-1, 1]\) 相关的参数。
- 对 \(k = 1, 2, ...\) 进行迭代，直到收敛：
a. 计算残差： \(\mathbf{r}^{(k)} = \mathbf{b} - A\mathbf{x}^{(k)}\)
b. 计算预处理残差： \(\mathbf{z}^{(k)} = M^{-1} \mathbf{r}^{(k)}\)
c. 更新解：
\[ \mathbf{x}^{(k+1)} = \frac{4}{(2 - \mu^2)\rho_k - 2} \left( \gamma \mathbf{z}^{(k)} + \mathbf{x}^{(k)} - \mathbf{x}^{(k-1)} \right) + \mathbf{x}^{(k-1)} \]
其中 \(\rho_k = 1 / (1 - \frac{\mu^2}{4} \rho_{k-1})\)， \(\rho_0 = 2\)， \(\mu = \delta\)。这个系数递推来自切比雪夫多项式在 \([-1,1]\) 上的最小最大性质。
在实际编程中，更常见的实现是以下简化形式（通过变量替换得到）：
设 \(\theta = \frac{\beta + \alpha}{2}\), \(\delta = \frac{\beta - \alpha}{2}\)。则令 \(\rho_0 = 1/\theta\)。对 \(k=0,1,2,...\)：
- 如果 \(k=0\)：
\[ \mathbf{x}^{(1)} = \mathbf{x}^{(0)} + \rho_0 M^{-1} \mathbf{r}^{(0)} \]
- 否则（\(k \ge 1\)）：
\[ \rho_k = \frac{1}{\theta - \frac{\delta^2}{4} \rho_{k-1}} \]
\[ \mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \rho_k \left( M^{-1} \mathbf{r}^{(k)} + \frac{\rho_{k-1} \delta^2}{4} (\mathbf{x}^{(k)} - \mathbf{x}^{(k-1)}) \right) \]
这个形式更容易实现，它直接由切比雪夫多项式的三项递推导出。
步骤4：将SSOR作为内迭代（预处理子）
在Chebyshev-SSOR中，每一步迭代都需要计算 \(\mathbf{z} = M^{-1} \mathbf{r}\)，即求解 \(M \mathbf{z} = \mathbf{r}\)。由于我们选择 \(M = M_{SSOR}\)，而 \(M_{SSOR} = \frac{1}{\omega(2-\omega)} (D - \omega L) D^{-1} (D - \omega L^T)\)，求解 \(M \mathbf{z} = \mathbf{r}\) 等价于执行一次SSOR迭代（从零初始解开始，右端项为 \(\mathbf{r}\)）。
具体来说，给定向量 \(\mathbf{r}\)，计算 \(\mathbf{z} = M_{SSOR}^{-1} \mathbf{r}\) 的步骤如下（这本质上是解 \(M_{SSOR} \mathbf{z} = \mathbf{r}\)）：
- 前向代入：解 \((D - \omega L) \mathbf{y} = \omega(2-\omega) \mathbf{r}\) 得到 \(\mathbf{y}\)。
- 因为 \(D - \omega L\) 是下三角矩阵，可以通过前向代入快速求解：
\[ y_i = \omega(2-\omega) r_i + \omega \sum_{j=1}^{i-1} a_{ij} y_j, \quad i=1,...,n \]
-
缩放：计算 \(\mathbf{v} = D^{-1} \mathbf{y}\)，即 \(v_i = y_i / a_{ii}\)。
-
后向代入：解 \((D - \omega L^T) \mathbf{z} = D \mathbf{v}\) 得到 \(\mathbf{z}\)。
- 因为 \(D - \omega L^T\) 是上三角矩阵，可以通过后向代入求解：
\[ z_i = v_i + \frac{\omega}{a_{ii}} \sum_{j=i+1}^{n} a_{ij} z_j, \quad i=n,...,1 \]
注意：在实际迭代中，残差向量 \(\mathbf{r}^{(k)}\) 是已知的，我们通过以上三步计算出 \(\mathbf{z}^{(k)} = M_{SSOR}^{-1} \mathbf{r}^{(k)}\)，然后代入Chebyshev加速的更新公式。
步骤5：完整算法流程与参数选择
- 输入：对称正定矩阵 \(A \in \mathbb{R}^{n \times n}\)，右端项 \(\mathbf{b} \in \mathbb{R}^n\)，初始猜测 \(\mathbf{x}^{(0)}\)（可设为零向量），松弛参数 \(\omega\)（通常取 \(1 < \omega < 2\)，如1.2, 1.5等，需实验调整），迭代次数上限 \(K_{max}\)，容差 \(\epsilon\)。
- 估计特征值边界：需要估计 \(M_{SSOR}^{-1}A\) 的特征值区间 \([\alpha, \beta]\)。这是一项挑战。常用方法包括：
- 通过若干步Lanczos迭代估计最大、最小特征值。
- 利用矩阵的性质给出理论界。例如，对于某些模型问题（如泊松方程离散），区间是已知的。
- 经验取值：有时取 \(\alpha = 0.1 \lambda_{\min}(A)\), \(\beta = \lambda_{\max}(A)\) 的粗略估计，然后通过试验微调。
- 如果难以估计，可退化为使用固定参数（此时退化为某种多项式加速，但非最优）。
- 计算常数：
\[ \theta = \frac{\beta + \alpha}{2}, \quad \delta = \frac{\beta - \alpha}{2}, \quad \rho_0 = 1/\theta \]
- 初始化：
- 计算初始残差： \(\mathbf{r}^{(0)} = \mathbf{b} - A\mathbf{x}^{(0)}\)
- 用SSOR作为预处理子求解： \(\mathbf{z}^{(0)} = M_{SSOR}^{-1} \mathbf{r}^{(0)}\)（即执行步骤4中的三步计算）。
- 更新解： \(\mathbf{x}^{(1)} = \mathbf{x}^{(0)} + \rho_0 \mathbf{z}^{(0)}\)
- 计算新残差： \(\mathbf{r}^{(1)} = \mathbf{b} - A\mathbf{x}^{(1)}\)
- 设置 \(k = 1\)。
- 迭代主循环：当 \(\|\mathbf{r}^{(k)}\| / \|\mathbf{b}\| > \epsilon\) 且 \(k < K_{max}\) 时：
a. 计算预处理残差： \(\mathbf{z}^{(k)} = M_{SSOR}^{-1} \mathbf{r}^{(k)}\)。
b. 计算Chebyshev系数： \(\rho_k = \frac{1}{\theta - \frac{\delta^2}{4} \rho_{k-1}}\)。
c. 更新解：
\[ \mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \rho_k \left( \mathbf{z}^{(k)} + \frac{\rho_{k-1} \delta^2}{4} (\mathbf{x}^{(k)} - \mathbf{x}^{(k-1)}) \right) \]
d. 计算新残差： $ \mathbf{r}^{(k+1)} = \mathbf{b} - A\mathbf{x}^{(k+1)} $。
e. $ k = k + 1 $。
- 输出：近似解 \(\mathbf{x}^{(k)}\)。
步骤6：算法特点与小结
- 加速效果：相比经典SSOR迭代，Chebyshev加速能显著减少迭代次数，特别是当 \(M^{-1}A\) 的特征值分布较为集中时（即 \(\alpha\) 和 \(\beta\) 接近）。其收敛速度与 \(\frac{\sqrt{\beta/\alpha} - 1}{\sqrt{\beta/\alpha} + 1}\) 相关，类似于预处理共轭梯度法（PCG），但无需计算内积，在某些并行架构上可能更有优势。
- 参数依赖：算法性能严重依赖于特征值区间 \([\alpha, \beta]\) 估计的准确性。过高估计 \(\beta/\alpha\) 比值会减慢收敛；而低估可能导致算法不稳定（参数 \(\rho_k\) 可能发散）。通常需要保守估计，确保区间包含所有特征值。
- 与CG对比：对于对称正定问题，预处理共轭梯度法（PCG）通常收敛更快且对特征值区间估计不敏感。但Chebyshev加速无需计算向量内积，避免了全局通信，在极端大规模并行计算中有时被采用。
- 适用性：特别适用于特征值边界已知或容易估计的问题（如来自规则网格离散的偏微分方程），且当SSOR是一个有效的平滑器或预处理子时。
通过以上步骤，我们将SSOR的稳定性和Chebyshev多项式的最优逼近性质结合，得到了一个高效、适合于求解大规模稀疏对称正定线性方程组的迭代算法。