分块矩阵的预处理共轭梯度法解多重右端项线性方程组
问题描述
考虑对称正定分块矩阵 \(A \in \mathbb{R}^{n \times n}\) 和多个右端项 \(B = [b_1, b_2, \dots, b_k] \in \mathbb{R}^{n \times k}\)，需解线性方程组 \(AX = B\)，其中 \(X = [x_1, x_2, \dots, x_k]\) 为待求解矩阵。直接对每个右端项独立使用共轭梯度法（CG）计算成本高。分块预处理共轭梯度法通过同时处理所有右端项，并利用分块矩阵结构及预处理技术加速收敛。
解题过程
- 问题形式化
将原问题拆分为 \(k\) 个独立子问题：
\[ A x_i = b_i, \quad i = 1, 2, \dots, k. \]
分块CG法的核心思想是同步迭代求解所有子问题，共享Krylov子空间信息以减少总迭代次数。
-
分块CG算法框架
- 初始化：
设 \(X_0\) 为初始解估计（常取零矩阵），计算初始残差 \(R_0 = B - A X_0\)。 - 预处理步骤：
选取分块预处理子 \(M\)（如分块对角预处理或不完全Cholesky分解），使得 \(M \approx A^{-1}\) 且易求逆。对残差应用预处理：\(Z_0 = M R_0\)。 - 搜索方向初始化：
设初始搜索方向 \(P_0 = Z_0\)。
- 初始化：
-
迭代过程
对迭代步 \(j = 0, 1, 2, \dots\) 执行：- 步长计算：
计算标量步长 \(\alpha_j\) 需最小化所有残差的范数。通过投影到当前搜索方向：
- 步长计算：
\[ \alpha_j = (R_j^\top Z_j) / (P_j^\top A P_j). \]
此处分子为残差与预处理残差的内积矩阵，分母为搜索方向与矩阵 $ A $ 的作用量。
- 更新解与残差：
\[ X_{j+1} = X_j + P_j \alpha_j, \quad R_{j+1} = R_j - A P_j \alpha_j. \]
- 预处理残差更新：
\[ Z_{j+1} = M R_{j+1}. \]
- 搜索方向更新：
计算系数 \(\beta_j = (R_{j+1}^\top Z_{j+1}) / (R_j^\top Z_j)\)，并更新搜索方向：
\[ P_{j+1} = Z_{j+1} + P_j \beta_j. \]
此步骤确保搜索方向的 $ A $-共轭性。
-
收敛判定
检查残差范数 \(\|R_{j+1}\|_F\)（Frobenius范数）是否小于预设阈值。若满足则终止，否则继续迭代。 -
预处理技术细节
- 分块对角预处理：取 \(M\) 为 \(A\) 的分块对角部分，逆易计算。
- 不完全Cholesky预处理：对 \(A\) 进行不完全Cholesky分解 \(A \approx L L^\top\)，则 \(M = (L L^\top)^{-1}\)。
预处理目的为降低条件数，加速收敛。
-
优势与复杂度
- 通过分块处理，矩阵-矩阵乘法（如 \(A P_j\)）可利用BLAS-3级运算优化。
- 共享Krylov子空间减少总迭代次数，尤其当右端项相关时效果显著。
- 复杂度从 \(O(k \cdot \kappa(A) n^2)\) 降至 \(O(\kappa(M^{-1}A) \cdot n^2 k)\)，其中 \(\kappa\) 为条件数。
此方法结合分块结构、预处理技术与共轭梯度法，高效处理多重右端项问题，适用于大规模科学计算。