基于Krylov子空间的块FOM（Block Full Orthogonalization Method）算法解多重右端项线性方程组
题目描述
考虑求解多重右端项的线性方程组：
\[A X = B \]
其中：
- \(A \in \mathbb{R}^{n \times n}\) 是一个大规模稀疏矩阵（非奇异），
- \(B \in \mathbb{R}^{n \times s}\) 是右端项矩阵（\(s \ll n\)），
- \(X \in \mathbb{R}^{n \times s}\) 是待求的解矩阵。
块FOM（Block FOM）是经典FOM（全正交化方法）的推广，旨在通过一次Krylov子空间迭代同时逼近多个线性方程组的解。该方法尤其适用于右端项具有相关性的情况（例如，在参数化PDE求解、模型降阶中），可大幅减少计算量。
解题过程
步骤1：问题转化与基本思想
- 多重右端项方程组可视为 \(s\) 个独立方程：
\[ A \mathbf{x}_j = \mathbf{b}_j, \quad j=1,\dots,s \]
经典FOM每次只能解一个方程，需重复 \(s\) 次，计算成本高。
- 块FOM的核心思想：为所有右端项构建一个共同的块Krylov子空间：
\[ \mathcal{K}_m(A, B) = \text{span}\{B, AB, A^2B, \dots, A^{m-1}B\} \]
其中 \(m\) 是子空间维数（通常 \(m \ll n\)）。
- 在该子空间中寻找近似解 \(X_m\)，满足：
\[ X_m \in \mathcal{K}_m(A, B) \]
且残差矩阵 \(R_m = B - A X_m\) 与子空间正交。
步骤2：构建块Krylov子空间的正交基
块Krylov子空间的基由一系列块向量组成，需进行块正交化（类似块Arnoldi过程）：
-
初始化：
- 对 \(B\) 进行QR分解：\(B = V_1 R_0\)，其中 \(V_1 \in \mathbb{R}^{n \times s}\) 列正交（\(V_1^T V_1 = I_s\)），\(R_0 \in \mathbb{R}^{s \times s}\) 是上三角矩阵。
- 令 \(V_1\) 作为正交基的第一块。
-
迭代生成后续基（\(j=1,2,\dots,m\)）：
a. 计算块矩阵乘法：\(W = A V_j\)。
b. 对 \(W\) 与已生成的正交基 \(\{V_1,\dots,V_j\}\) 进行正交化：
\[ H_{i,j} = V_i^T W, \quad i=1,\dots,j \]
\[ \tilde{W} = W - \sum_{i=1}^{j} V_i H_{i,j} \]
c. 对 \(\tilde{W}\) 进行QR分解：\(\tilde{W} = V_{j+1} H_{j+1,j}\)，其中 \(V_{j+1} \in \mathbb{R}^{n \times s}\) 列正交，\(H_{j+1,j} \in \mathbb{R}^{s \times s}\) 是上三角矩阵。
d. 保存 \(H_{i,j}\) 和 \(H_{j+1,j}\) 构成块上Hessenberg矩阵 \(\mathcal{H}_m\)。
步骤3：构建近似解
- 经过 \(m\) 步迭代，得到正交基矩阵：
\[ \mathcal{V}_m = [V_1, V_2, \dots, V_m] \in \mathbb{R}^{n \times ms} \]
和块上Hessenberg矩阵：
\[ \mathcal{H}_m = \begin{pmatrix} H_{1,1} & H_{1,2} & \cdots & H_{1,m} \\ H_{2,1} & H_{2,2} & \cdots & H_{2,m} \\ 0 & H_{3,2} & \cdots & H_{3,m} \\ \vdots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & H_{m+1,m} \end{pmatrix} \in \mathbb{R}^{(m+1)s \times ms} \]
满足块Arnoldi关系：
\[ A \mathcal{V}_m = \mathcal{V}_{m+1} \mathcal{H}_m \]
-
在子空间中设近似解为 \(X_m = \mathcal{V}_m Y_m\)，其中 \(Y_m \in \mathbb{R}^{ms \times s}\) 待定。
-
代入原方程：
\[ A X_m = A \mathcal{V}_m Y_m = \mathcal{V}_{m+1} \mathcal{H}_m Y_m \]
残差矩阵：
\[ R_m = B - A X_m = V_1 R_0 - \mathcal{V}_{m+1} \mathcal{H}_m Y_m \]
由于 \(V_1 R_0 = \mathcal{V}_{m+1} (E_1 R_0)\)，其中 \(E_1 = [I_s, 0, \dots, 0]^T \in \mathbb{R}^{(m+1)s \times s}\)，有：
\[ R_m = \mathcal{V}_{m+1} (E_1 R_0 - \mathcal{H}_m Y_m) \]
- 强制残差与子空间 \(\mathcal{K}_m(A, B)\) 正交（即与 \(\mathcal{V}_m\) 的所有列正交）：
\[ \mathcal{V}_m^T R_m = 0 \]
因为 \(\mathcal{V}_{m+1}\) 列正交，等价于：
\[ E_1 R_0 - \bar{\mathcal{H}}_m Y_m = 0 \]
其中 \(\bar{\mathcal{H}}_m \in \mathbb{R}^{ms \times ms}\) 是 \(\mathcal{H}_m\) 的前 \(m\) 个块行（去掉最后一行 \(H_{m+1,m}\)）。这得到一个块上Hessenberg线性方程组。
步骤4：求解约化系统并得到解
- 解块方程组：
\[ \bar{\mathcal{H}}_m Y_m = E_1 R_0 \]
由于 \(\bar{\mathcal{H}}_m\) 是块上三角矩阵，可通过块前向回代快速求解。
- 得到近似解：
\[ X_m = \mathcal{V}_m Y_m \]
- 若残差范数 \(\|R_m\|_F\) 未满足容忍度，可增加 \(m\) 并继续扩展Krylov子空间（类似重启策略）。
关键点与注意事项
- 块正交化稳定性：需使用稳定的QR分解（如Householder或Modified Gram-Schmidt with reorthogonalization）避免基向量丢失正交性。
- 存储与计算：基矩阵 \(\mathcal{V}_m\) 大小为 \(n \times ms\)，当 \(ms\) 很大时需重启或使用稀疏存储。
- 收敛性：若右端项 \(B\) 的列张成的子空间与 \(A\) 的某些特征向量方向高度相关，块FOM可比单独求解每个方程更快收敛。
- 与块GMRES区别：块FOM要求残差与子空间正交（Galerkin条件），块GMRES则最小化残差范数（Petrov-Galerkin条件）。
总结
块FOM通过块Arnoldi过程构建共享的Krylov子空间，将多重右端项问题转化为低维块上Hessenberg方程组，显著减少了迭代次数和矩阵-向量乘法的开销。它特别适合右端项相关或需要同时求解多个类似方程组的场景，是大规模科学计算中高效的多重右端项求解器之一。