基于Krylov子空间方法的FOM（全正交化方法）算法解线性方程组
题目描述
考虑求解大型稀疏非对称线性方程组：
\[A\mathbf{x} = \mathbf{b} \]
其中 \(A \in \mathbb{R}^{n \times n}\) 是一个大型稀疏方阵（可能非对称），\(\mathbf{b} \in \mathbb{R}^n\) 是已知右端向量，\(\mathbf{x} \in \mathbb{R}^n\) 是待求解向量。由于 \(A\) 规模大且稀疏，直接法（如高斯消元）在存储和计算上代价高，因此常采用Krylov子空间迭代法。FOM（Full Orthogonalization Method，全正交化方法）是其中一种基于Arnoldi过程构建正交基，并将原问题投影到Krylov子空间上求解的算法。它适用于非对称矩阵，并在Krylov子空间上构建最小残差解。
循序渐进讲解
1. 核心思想
FOM的核心思想是：在每一步迭代 \(m\) 中，构造一个Krylov子空间：
\[\mathcal{K}_m(A, \mathbf{r}_0) = \text{span}\{\mathbf{r}_0, A\mathbf{r}_0, A^2\mathbf{r}_0, \dots, A^{m-1}\mathbf{r}_0\} \]
其中 \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\) 是初始残差，\(\mathbf{x}_0\) 是初始猜测解（常取零向量）。FOM在该子空间中寻找近似解 \(\mathbf{x}_m\)，使得残差 \(\mathbf{r}_m = \mathbf{b} - A\mathbf{x}_m\) 与子空间正交，从而得到投影方程。
2. 算法步骤详解
步骤1：初始化
- 选择初始猜测解 \(\mathbf{x}_0\)（例如 \(\mathbf{x}_0 = \mathbf{0}\)）。
- 计算初始残差：\(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\)。
- 令 \(\beta = \|\mathbf{r}_0\|_2\)，并设置第一个基向量：\(\mathbf{v}_1 = \mathbf{r}_0 / \beta\)。
目的：构建Krylov子空间的第一组标准正交基。
步骤2：Arnoldi过程构建正交基
对 \(j = 1, 2, \dots, m\) 执行：
- 计算 \(\mathbf{w} = A\mathbf{v}_j\)。
- 对 \(i = 1, 2, \dots, j\) 执行正交化（采用修正Gram-Schmidt）：
\[ h_{i,j} = \mathbf{v}_i^T \mathbf{w}, \quad \mathbf{w} = \mathbf{w} - h_{i,j} \mathbf{v}_i \]
- 计算 \(h_{j+1,j} = \|\mathbf{w}\|_2\)。
- 如果 \(h_{j+1,j} = 0\)，则子空间已不变，提前终止。
- 否则，令 \(\mathbf{v}_{j+1} = \mathbf{w} / h_{j+1,j}\)。
结果：得到正交矩阵 \(V_m = [\mathbf{v}_1, \dots, \mathbf{v}_m] \in \mathbb{R}^{n \times m}\) 和上Hessenberg矩阵 \(H_m \in \mathbb{R}^{m \times m}\)（元素为 \(h_{i,j}\)），满足：
\[A V_m = V_m H_m + h_{m+1,m} \mathbf{v}_{m+1} \mathbf{e}_m^T \]
其中 \(\mathbf{e}_m\) 是第 \(m\) 个单位向量。
步骤3：构造投影方程
FOM在子空间 \(\mathcal{K}_m\) 中寻找近似解 \(\mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m\)，其中 \(\mathbf{y}_m \in \mathbb{R}^m\) 是系数向量。代入原方程：
\[A(\mathbf{x}_0 + V_m \mathbf{y}_m) \approx \mathbf{b} \]
利用 Arnoldi 关系 \(A V_m = V_m H_m + h_{m+1,m} \mathbf{v}_{m+1} \mathbf{e}_m^T\) 和 \(\mathbf{r}_0 = \beta \mathbf{v}_1\)，可得：
\[\mathbf{b} - A\mathbf{x}_0 = \beta \mathbf{v}_1 = V_m H_m \mathbf{y}_m + h_{m+1,m} \mathbf{e}_m^T \mathbf{y}_m \mathbf{v}_{m+1} \]
要求残差 \(\mathbf{r}_m = \mathbf{b} - A\mathbf{x}_m\) 与子空间 \(\mathcal{K}_m\) 正交，即 \(V_m^T \mathbf{r}_m = 0\)。将上式左乘 \(V_m^T\)，注意到 \(V_m^T V_m = I_m\) 且 \(V_m^T \mathbf{v}_{m+1} = 0\)，得到投影方程：
\[H_m \mathbf{y}_m = \beta \mathbf{e}_1 \]
其中 \(\mathbf{e}_1 = [1, 0, \dots, 0]^T \in \mathbb{R}^m\)。
物理意义：将原 \(n\) 维问题约化为 \(m\) 维上Hessenberg线性方程组。
步骤4：求解约化系统并更新解
- 求解 \(m\) 维线性方程组：\(H_m \mathbf{y}_m = \beta \mathbf{e}_1\)。由于 \(H_m\) 是上Hessenberg矩阵，可用Givens旋转或Householder变换稳定求解（类似稠密QR分解）。
- 得到近似解：\(\mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m\)。
注意：由于子空间维数 \(m\) 通常远小于 \(n\)，该方程组可通过直接法高效求解。
步骤5：检查收敛性与重启策略
- 计算残差范数：理论上，FOM的残差满足 \(\|\mathbf{r}_m\|_2 = h_{m+1,m} |\mathbf{e}_m^T \mathbf{y}_m|\)，可利用已计算的量直接得到，无需显式计算 \(A\mathbf{x}_m\)。
- 若 \(\|\mathbf{r}_m\|_2 < \epsilon\)（预设容差），则停止并输出 \(\mathbf{x}_m\)。
- 若未收敛且 \(m\) 达到最大允许维数（如 \(m=30\)），则采用重启策略：以当前解 \(\mathbf{x}_m\) 作为新的初始猜测 \(\mathbf{x}_0\)，重新计算残差并重复Arnoldi过程，避免子空间过大导致存储和计算成本增加。
3. 算法特性与注意事项
- 优点：适用于非对称矩阵；在Krylov子空间上残差正交，通常收敛稳定；无需存储整个矩阵 \(A\)，只需矩阵-向量乘法。
- 缺点：需存储全部基向量 \(V_m\)，存储成本为 \(O(mn)\)；每步需完全正交化，计算量为 \(O(m^2 n)\)。重启可缓解内存压力，但可能影响收敛速度。
- 与GMRES关系：FOM是GMRES的前身，区别在于GMRES最小化残差范数 \(\|\mathbf{b} - A\mathbf{x}_m\|_2\)，而FOM强制残差与子空间正交。对于对称正定矩阵，FOM退化为共轭梯度法（CG）。
- 实际应用：FOM常用于中等规模非对称问题，或作为理论分析工具。对于大规模问题，常使用重启策略或转向GMRES等更高效变体。
总结
FOM通过Arnoldi过程构建Krylov子空间的正交基，将原方程投影到低维上Hessenberg系统求解。其核心是正交性条件和投影技术，结合重启策略平衡存储与收敛性。理解FOM有助于掌握Krylov子空间方法的基础，并为学习GMRES、Arnoldi特征值算法等提供桥梁。