随机化Chebyshev迭代法求解大型稀疏线性方程组
1. 问题描述
我们考虑求解大型稀疏线性方程组
\[A\mathbf{x} = \mathbf{b}, \]
其中 \(A \in \mathbb{R}^{n \times n}\) 是对称正定（SPD）矩阵，\(\mathbf{b} \in \mathbb{R}^n\) 是已知向量，\(\mathbf{x} \in \mathbb{R}^n\) 是待求解向量。在许多科学与工程问题中（如偏微分方程离散化），\(A\) 通常非常大（\(n\) 可达数百万甚至更大）且稀疏。直接法（如 LU 分解）因内存和计算成本过高而不适用，因此迭代法是更可行的选择。
Chebyshev 迭代法是一种多项式迭代方法，它通过构造一系列 Chebyshev 多项式来加速收敛，特别适用于系数矩阵的特征值边界已知的情况。其基本思想是：将解向量表示为初始猜测与一系列修正项的线性组合，这些修正项由 Chebyshev 多项式生成，以最小化迭代误差在某种范数下的上界。
2. 核心思想
Chebyshev 迭代法属于多项式加速方法。它不依赖于 Krylov 子空间的正交性（如共轭梯度法），而是利用 Chebyshev 多项式在区间 \([-1,1]\) 上的最小最大性质，将矩阵 \(A\) 的特征值映射到该区间，从而设计迭代格式，使误差在谱半径上快速衰减。
3. 逐步推导与算法细节
步骤 1：将问题转化为多项式优化问题
设第 \(k\) 步迭代的解为 \(\mathbf{x}_k\)，误差为 \(\mathbf{e}_k = \mathbf{x}_* - \mathbf{x}_k\)（\(\mathbf{x}_*\) 为精确解）。由于 \(A\) 对称正定，其特征值满足 \(0 < \lambda_{\min} = \lambda_1 \leq \lambda_2 \leq \dots \leq \lambda_n = \lambda_{\max}\)。
我们希望通过一个 \(k\) 次多项式 \(P_k(\lambda)\) 满足 \(P_k(0) = 1\)，来构造迭代：
\[\mathbf{e}_k = P_k(A) \mathbf{e}_0, \]
其中 \(\mathbf{e}_0\) 是初始误差。我们希望选择 \(P_k\) 使得 \(\|P_k(A)\|_2\) 尽可能小，从而误差衰减快。因为 \(A\) 对称，有
\[\|P_k(A)\|_2 = \max_{\lambda \in [\lambda_{\min},\lambda_{\max}]} |P_k(\lambda)|. \]
于是问题转化为：在多项式集合 \(\{P_k: P_k(0)=1\}\) 中，寻找一个多项式使得
\[\max_{\lambda \in [\lambda_{\min},\lambda_{\max}]} |P_k(\lambda)| \]
最小。
步骤 2：Chebyshev多项式的最小最大性质
标准 Chebyshev 多项式 \(T_k(t)\) 定义在 \([-1,1]\) 上：
\[T_k(t) = \cos(k \arccos t), \quad |t| \leq 1. \]
它在 \([-1,1]\) 上满足 \(|T_k(t)| \leq 1\)，并且在首项系数为 \(2^{k-1}\) 的所有 \(k\) 次多项式中，\(T_k(t)/2^{k-1}\) 是零点的唯一多项式，在 \([-1,1]\) 上与零的偏差最小。
我们需要将区间 \([\lambda_{\min},\lambda_{\max}]\) 映射到 \([-1,1]\)。设
\[\lambda = \frac{\lambda_{\max} + \lambda_{\min}}{2} + \frac{\lambda_{\max} - \lambda_{\min}}{2} t, \quad t \in [-1,1]. \]
定义缩放平移后的 Chebyshev 多项式：
\[\tilde{P}_k(\lambda) = \frac{T_k\left( \frac{\lambda_{\max} + \lambda_{\min} - 2\lambda}{\lambda_{\max} - \lambda_{\min}} \right)}{T_k\left( \frac{\lambda_{\max} + \lambda_{\min}}{\lambda_{\max} - \lambda_{\min}} \right)}. \]
注意 \(\tilde{P}_k(0) = 1\)，且 \(|\tilde{P}_k(\lambda)| \leq 1 / |T_k(\gamma)|\) 对于 \(\lambda \in [\lambda_{\min},\lambda_{\max}]\)，其中
\[\gamma = \frac{\lambda_{\max} + \lambda_{\min}}{\lambda_{\max} - \lambda_{\min}} > 1. \]
由于 \(|T_k(\gamma)|\) 随 \(k\) 指数增长，误差会指数衰减。
步骤 3：递推公式的导出
利用 Chebyshev 多项式的三项递推关系：
\[T_{k+1}(t) = 2t T_k(t) - T_{k-1}(t), \quad T_0(t)=1, \ T_1(t)=t. \]
令
\[\alpha = \frac{\lambda_{\max} + \lambda_{\min}}{2}, \quad \beta = \frac{\lambda_{\max} - \lambda_{\min}}{2}, \quad t = \frac{\alpha - \lambda}{\beta}. \]
定义 \(C_k = T_k(\alpha / \beta)\)，则
\[\tilde{P}_k(\lambda) = \frac{T_k((\alpha - \lambda)/\beta)}{C_k}. \]
代入三项递推，可得到解向量迭代格式。经过推导（过程略），得到 Chebyshev迭代法 的递推形式：
初始化：选 \(\mathbf{x}_0\)，计算 \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\)。
设 \(\rho = \frac{\lambda_{\max} - \lambda_{\min}}{\lambda_{\max} + \lambda_{\min}}\)，并取
\[\omega_1 = \frac{2}{2\lambda_{\max} - \lambda_{\min}}, \quad \mathbf{x}_1 = \mathbf{x}_0 + \omega_1 \mathbf{r}_0. \]
对于 \(k \geq 1\)：
\[\omega_{k+1} = \frac{4}{4\lambda_{\max} - \lambda_{\min} - \left( \frac{\lambda_{\max} - \lambda_{\min}}{2} \right)^2 \omega_k}, \]
\[ \mathbf{x}_{k+1} = \mathbf{x}_k + \omega_{k+1} (\mathbf{b} - A\mathbf{x}_k) + \frac{\omega_{k+1}}{\omega_k} (1 - \omega_k \lambda_{\min}) (\mathbf{x}_k - \mathbf{x}_{k-1}). \]
注意：实际常用简化形式（见下一步）。
步骤 4：简化递推公式（实用形式）
令
\[\theta = \frac{\lambda_{\max} + \lambda_{\min}}{2}, \quad \delta = \frac{\lambda_{\max} - \lambda_{\min}}{2}, \quad \rho = \frac{\delta}{\theta}. \]
并设
\[\tau_k = T_k(1/\rho), \quad \sigma_k = \frac{2}{\delta} \cdot \frac{\tau_{k-1}}{\tau_k}. \]
可得更稳定的递推（避免直接使用 \(\lambda_{\min}, \lambda_{\max}\) 造成数值敏感）：
- 初始化：\(\mathbf{x}_0\) 任意，\(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\)，\(\mathbf{x}_1 = \mathbf{x}_0 + \frac{2}{2\lambda_{\max} - \lambda_{\min}} \mathbf{r}_0\)。
- 对 \(k=1,2,\dots\) 直到收敛：
\[\mathbf{r}_k = \mathbf{b} - A\mathbf{x}_k, \]
\[ \omega_{k+1} = \left[ \theta - \frac{\delta^2}{4} \omega_k \right]^{-1} \cdot 2, \]
\[ \mathbf{x}_{k+1} = \mathbf{x}_k + \omega_{k+1} \mathbf{r}_k + \frac{\omega_{k+1}}{\omega_k} (1 - \omega_k \lambda_{\min}) (\mathbf{x}_k - \mathbf{x}_{k-1}). \]
步骤 5：收敛性分析
误差满足：
\[\frac{\|\mathbf{e}_k\|_A}{\|\mathbf{e}_0\|_A} \leq 2 \left( \frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1} \right)^k, \]
其中 \(\kappa = \lambda_{\max}/\lambda_{\min}\) 是条件数，\(\|\mathbf{e}\|_A = \sqrt{\mathbf{e}^T A \mathbf{e}}\) 是能量范数。
与最速下降法相比，Chebyshev 迭代避免了最坏情况的方向，收敛速度由 \(\kappa\) 决定，且不需要计算内积（无正交化），适合并行计算。
步骤 6：与共轭梯度法（CG）的比较
- CG 在能量范数下是最优的，但需要全局内积（并行时通信成本高）。
- Chebyshev 迭代不需要内积，只需已知 \([\lambda_{\min},\lambda_{\max}]\)，适合分布式计算。
- 如果特征值边界估计不准，Chebyshev 可能收敛慢甚至发散，CG 则对谱分布不敏感（但仍受条件数影响）。
- Chebyshev 可作为预处理子的加速方法（多项式预处理）。
步骤 7：算法总结
输入：\(A, \mathbf{b}\)，初值 \(\mathbf{x}_0\)，\(\lambda_{\min}\)，\(\lambda_{\max}\)，容差 \(\epsilon\)，最大步数 \(K\)。
输出：近似解 \(\mathbf{x}_k\)。
- 计算 \(\theta = (\lambda_{\max} + \lambda_{\min})/2\)，\(\delta = (\lambda_{\max} - \lambda_{\min})/2\)，\(\rho = \delta/\theta\)。
- 计算 \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\)。
- 设 \(\omega_1 = 2 / (2\lambda_{\max} - \lambda_{\min})\)，\(\mathbf{x}_1 = \mathbf{x}_0 + \omega_1 \mathbf{r}_0\)。
- 对 \(k = 1,2,\dots,K\)：
- 计算残差 \(\mathbf{r}_k = \mathbf{b} - A\mathbf{x}_k\)。
- 如果 \(\|\mathbf{r}_k\|_2 < \epsilon\)，退出。
- 计算 \(\omega_{k+1} = 4 / (4\lambda_{\max} - \lambda_{\min} - \delta^2 \omega_k)\)。
- 更新：
\[ \mathbf{x}_{k+1} = \mathbf{x}_k + \omega_{k+1} \mathbf{r}_k + \frac{\omega_{k+1}}{\omega_k} (1 - \omega_k \lambda_{\min}) (\mathbf{x}_k - \mathbf{x}_{k-1}). \]
- 返回 \(\mathbf{x}_k\)。
8. 关键点
- Chebyshev 迭代法需要预先估计 \(\lambda_{\min}\) 和 \(\lambda_{\max}\)，可用 Gershgorin 圆定理或少量迭代（如幂法）得到。
- 该方法特别适合在已知谱边界且矩阵-向量乘易于并行的情况下，作为 CG 的替代。
- 若特征值分布不均匀，可结合多项式预处理（如 Chebyshev 多项式平滑）改善收敛。