基于Krylov子空间方法的FOM（全正交化方法）算法解线性方程组
题目描述
本题目讲解如何运用全正交化方法（Full Orthogonalization Method, FOM） 这一Krylov子空间投影算法，来求解大规模稀疏非对称线性方程组 \(A\mathbf{x} = \mathbf{b}\)。其中，\(A \in \mathbb{R}^{n \times n}\) 是一个大型稀疏矩阵（通常非对称），\(\mathbf{b} \in \mathbb{R}^n\) 是已知右端项，目标是求解未知向量 \(\mathbf{x} \in \mathbb{R}^n\)。
FOM 属于一类重要的迭代法，它不直接对 \(A\) 进行分解，而是通过在一个逐渐增大的Krylov子空间中寻找近似解来高效求解。我们将从算法思想、数学推导、实现步骤和收敛性等方面，循序渐进地进行讲解。
解题过程
第一步：理解Krylov子空间投影法的基本思想
对于大型稀疏线性方程组，直接法（如LU分解）因内存和计算量过大往往不可行。投影迭代法的核心思想是：在第 \(m\) 步迭代（\(m \ll n\)），在一个 \(m\) 维的子空间 \(\mathcal{K}_m\) 中寻找近似解 \(\mathbf{x}_m\)，并强制要求残差 \(\mathbf{r}_m = \mathbf{b} - A\mathbf{x}_m\) 与另一个 \(m\) 维子空间 \(\mathcal{L}_m\) 正交。这种正交条件称为 Petrov-Galerkin 条件。
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\) 是初始猜测（常取为零向量）。
- 约束空间：与搜索空间相同，即 \(\mathcal{L}_m = \mathcal{K}_m\)。
因此，FOM 的 Petrov-Galerkin 条件为：
\[\mathbf{r}_m = \mathbf{b} - A\mathbf{x}_m \perp \mathcal{K}_m(A, \mathbf{r}_0) \]
这意味着残差与整个当前的Krylov子空间正交。
第二步：建立Krylov子空间的标准正交基（Arnoldi过程）
为了在子空间 \(\mathcal{K}_m\) 中方便地表示向量并施加正交条件，我们需要为其构造一组标准正交基 \(\{\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_m\}\)。这通过 Arnoldi迭代 实现，该过程可看作是对矩阵 \(A\) 进行不完整的正交相似变换，得到一个上Hessenberg矩阵 \(H_m\)。
Arnoldi过程的详细步骤：
- 初始化：给定初始残差 \(\mathbf{r}_0\)。计算 \(\beta = \|\mathbf{r}_0\|_2\)，并设第一个基向量 \(\mathbf{v}_1 = \mathbf{r}_0 / \beta\)。
- 迭代构造：对于 \(j = 1, 2, \dots, m\) 执行：
a. 计算矩阵向量积：\(\mathbf{w} = A \mathbf{v}_j\)。
b. 全正交化：对 \(\mathbf{w}\) 关于所有已生成的基向量 \(\mathbf{v}_1, \dots, \mathbf{v}_j\) 进行Gram-Schmidt正交化：
\[ h_{ij} = (\mathbf{w}, \mathbf{v}_i), \quad \text{for } i = 1, \dots, j \]
\[ \mathbf{w} \leftarrow \mathbf{w} - \sum_{i=1}^{j} h_{ij} \mathbf{v}_i \]
c. 计算 $ h_{j+1, j} = \|\mathbf{w}\|_2 $。
d. 如果 $ h_{j+1, j} = 0 $，则Krylov子空间已达到不变子空间，迭代可提前终止。否则，设 $ \mathbf{v}_{j+1} = \mathbf{w} / h_{j+1, j} $。
这个过程产生一个 \(n \times (m+1)\) 的列正交矩阵 \(V_{m+1} = [\mathbf{v}_1, \dots, \mathbf{v}_{m+1}]\)，以及一个 \((m+1) \times m\) 的上Hessenberg矩阵 \(\bar{H}_m\)，其前 \(m\) 行构成方阵 \(H_m\)，满足核心关系式：
\[A V_m = V_{m+1} \bar{H}_m = V_m H_m + h_{m+1, m} \mathbf{v}_{m+1} \mathbf{e}_m^T \]
其中 \(\mathbf{e}_m\) 是第 \(m\) 个单位坐标向量。
第三步：在Krylov子空间中表示近似解并施加正交条件
我们设近似解 \(\mathbf{x}_m\) 可表示为：
\[\mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m \]
其中 \(\mathbf{y}_m \in \mathbb{R}^m\) 是系数向量。于是残差为：
\[\mathbf{r}_m = \mathbf{b} - A\mathbf{x}_m = \mathbf{r}_0 - A V_m \mathbf{y}_m \]
由于 \(\mathbf{r}_0 = \beta \mathbf{v}_1 = V_{m+1} (\beta \mathbf{e}_1)\)（这里 \(\beta = \|\mathbf{r}_0\|_2\)，且 \(\mathbf{e}_1\) 是第一个单位坐标向量），代入 Arnoldi 关系式得：
\[\mathbf{r}_m = V_{m+1} (\beta \mathbf{e}_1 - \bar{H}_m \mathbf{y}_m) \]
FOM 要求 \(\mathbf{r}_m \perp \mathcal{K}_m = \text{span}(V_m)\)，即 \(V_m^T \mathbf{r}_m = 0\)。由于 \(V_m^T V_{m+1} = [I_m, \mathbf{0}]\)，这个正交条件等价于：
\[\beta \mathbf{e}_1 - H_m \mathbf{y}_m = 0 \]
因此，我们得到了一个 \(m \times m\) 的上Hessenberg线性方程组：
\[H_m \mathbf{y}_m = \beta \mathbf{e}_1 \]
第四步：高效求解缩减后的线性方程组
上Hessenberg系统 \(H_m \mathbf{y}_m = \beta \mathbf{e}_1\) 的规模 \(m\) 远小于原问题规模 \(n\)。我们可以用 Givens旋转 对其进行QR分解，从而高效稳定地求解。
- 初始化：设 \(Q^{(0)} = I_m\)，\(R^{(0)} = H_m\)，\(\mathbf{g}^{(0)} = \beta \mathbf{e}_1\)。
- 应用Givens旋转：对于 \(i = 1, \dots, m-1\)，构造Givens旋转矩阵 \(G_i\)，作用于 \(R^{(i-1)}\) 的第 \(i\) 行和第 \(i+1\) 行，以消去 \(R^{(i-1)}\) 在第 \((i+1, i)\) 位置的元素。同时，对右端项做同样的旋转：\(\mathbf{g}^{(i)} = G_i \mathbf{g}^{(i-1)}\)。
- 回代求解：经过 \(m-1\) 次旋转后，我们得到上三角矩阵 \(R = R^{(m-1)}\) 和变换后的右端项 \(\mathbf{g} = \mathbf{g}^{(m-1)}\)。然后通过回代求解上三角系统 \(R \mathbf{y}_m = \mathbf{g}\)。
- 重构近似解：计算 \(\mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m\)。
优势：每次迭代增加一维时，只需对新的一列 \(H_m\) 应用一个额外的Givens旋转，更新是增量的，计算代价低。
第五步：收敛性判断与重新启动策略
- 收敛判断：近似解的残差范数可以通过公式高效计算，无需显式计算 \(\mathbf{r}_m\)：
\[ \|\mathbf{r}_m\|_2 = |h_{m+1, m}| \cdot |\mathbf{e}_m^T \mathbf{y}_m| \]
当 \(\|\mathbf{r}_m\|_2\) 小于预设的收敛容差 \(\epsilon\) 时，迭代终止，\(\mathbf{x}_m\) 即为满足精度要求的解。
- 重新启动的必要性：FOM 是“完全”正交化，每一步都需要与所有之前的基向量做内积，其计算量和存储量随迭代步数 \(m\) 线性增长（每次矩阵向量乘一次，内积 \(O(m)\) 个，存储 \(m\) 个基向量）。当 \(m\) 很大时，计算和存储变得不可行。此外，数值误差的积累也可能导致基向量失去正交性。
- 重新启动FOM：一种常用策略是设定一个最大子空间维数 \(m_{\max}\)。当迭代步数达到 \(m_{\max}\) 仍未收敛时，我们取当前近似解 \(\mathbf{x}_{m_{\max}}\) 作为新的初始猜测 \(\mathbf{x}_0\)，重新计算初始残差，然后重启一个新的 Arnoldi 过程。这个过程可以重复，直到满足收敛条件。
第六步：算法总结与特性
FOM 算法流程总结如下：
- 输入：矩阵 \(A\)，右端项 \(\mathbf{b}\)，初始解 \(\mathbf{x}_0\)，容差 \(\epsilon\)，最大内迭代步数 \(m_{\max}\)。
- 外循环（重新启动循环）：
a. 计算残差：\(\mathbf{r} = \mathbf{b} - A\mathbf{x}_0\)，\(\beta = \|\mathbf{r}\|_2\)。若 \(\beta < \epsilon\)，输出 \(\mathbf{x}_0\) 并终止。
b. 设 \(\mathbf{v}_1 = \mathbf{r} / \beta\)。
c. 内循环（Arnoldi过程与FOM求解）：对于 \(j = 1, 2, \dots, m_{\max}\)：
i. 执行 Arnoldi 过程一步，得到新的基向量 \(\mathbf{v}_{j+1}\) 和 \(H_j\) 的新一列。
ii. 更新对上Hessenberg系统 \(H_j \mathbf{y}_j = \beta \mathbf{e}_1\) 的QR分解（通过一个Givens旋转）。
iii. 计算当前残差范数的估计 \(\|\mathbf{r}_j\|_2\)。
iv. 如果 \(\|\mathbf{r}_j\|_2 < \epsilon\)，则计算 \(\mathbf{y}_j\)（回代），更新解 \(\mathbf{x} = \mathbf{x}_0 + V_j \mathbf{y}_j\)，输出 \(\mathbf{x}\) 并终止。
d. 若内循环完成（\(j = m_{\max}\)）仍未收敛，则计算 \(\mathbf{y}_{m_{\max}}\)，更新 \(\mathbf{x}_0 = \mathbf{x}_0 + V_{m_{\max}} \mathbf{y}_{m_{\max}}\)，然后返回步骤 a 进行重启。
FOM的特性：
- 它是精确的投影方法，在无限精度的算术中，第 \(m\) 步近似解 \(\mathbf{x}_m\) 满足残差与 \(\mathcal{K}_m\) 正交，这等价于在 \(\mathcal{K}_m\) 中最小化某个与 \(A\) 相关的范数（当 \(A\) 对称正定时即为能量范数）。
- 对于对称矩阵，FOM 退化为著名的 共轭梯度法（CG）。
- 相比于GMRES（另一种流行的Krylov方法，它在 \(\mathcal{K}_m\) 中直接最小化残差2-范数），FOM 不需要存储上Hessenberg矩阵的QR分解以外的额外信息来重构解，但其残差范数不是单调递减的，收敛曲线可能不如GMRES平滑。
通过以上步骤，FOM 为求解大规模稀疏非对称线性方程组提供了一个基于Krylov子空间投影的强大工具。其核心在于利用Arnoldi过程构造正交基，将大规模问题约化为一个小型上Hessenberg系统，并通过高效稳定的数值线性代数技术求解。