基于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\) 是未知解向量。
直接解法（如高斯消元）因内存和时间消耗过大而不适用，而迭代法（如GMRES）是更合适的选择。FOM（Full Orthogonalization Method）是一种基于Krylov子空间的投影方法，它通过构造一组标准正交基，将原始问题投影到较低维的子空间上求解，再回代得到近似解。
解题过程
FOM的核心思想是：在Krylov子空间 \(\mathcal{K}_m(A, \mathbf{r}_0) = \text{span}\{\mathbf{r}_0, A\mathbf{r}_0, \dots, A^{m-1}\mathbf{r}_0\}\) 中寻找近似解 \(\mathbf{x}_m\)，其中 \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\) 是初始残差（通常取初始猜测 \(\mathbf{x}_0 = \mathbf{0}\)），\(m \ll n\) 是子空间维度。
步骤1：构造Krylov子空间的标准正交基
使用Arnoldi过程生成一组标准正交基向量 \(\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_m\)，它们张成Krylov子空间 \(\mathcal{K}_m(A, \mathbf{r}_0)\)。具体步骤如下：
- 初始化：
- 计算初始残差 \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\)。
- 令 \(\mathbf{v}_1 = \mathbf{r}_0 / \|\mathbf{r}_0\|_2\)。
- 对 \(j = 1, 2, \dots, m\) 执行Arnoldi迭代：
a. 计算 \(\mathbf{w}_j = A \mathbf{v}_j\)。
b. 对 \(i = 1, \dots, j\) 执行正交化（使用修正的Gram-Schmidt过程）：- 计算内积 \(h_{i,j} = \mathbf{v}_i^T \mathbf{w}_j\)。
- 更新 \(\mathbf{w}_j := \mathbf{w}_j - h_{i,j} \mathbf{v}_i\)。
c. 计算 \(h_{j+1,j} = \|\mathbf{w}_j\|_2\)。
d. 如果 \(h_{j+1,j} = 0\)，则提前终止（此时子空间已不变）。
e. 否则，令 \(\mathbf{v}_{j+1} = \mathbf{w}_j / h_{j+1,j}\)。
这个过程可以写成矩阵形式：
\[ A V_m = V_{m+1} \bar{H}_m, \]
其中 \(V_m = [\mathbf{v}_1, \dots, \mathbf{v}_m] \in \mathbb{R}^{n \times m}\) 的列是标准正交基，\(\bar{H}_m \in \mathbb{R}^{(m+1) \times m}\) 是一个上Hessenberg矩阵（其元素为 \(h_{i,j}\)）。
步骤2：将原问题投影到Krylov子空间
FOM寻找近似解 \(\mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m\)，其中 \(\mathbf{y}_m \in \mathbb{R}^m\) 是子空间中的系数向量。
将 \(A \mathbf{x}_m = \mathbf{b}\) 代入，并利用 Arnoldi 关系式：
\[A (\mathbf{x}_0 + V_m \mathbf{y}_m) = \mathbf{b} \quad \Rightarrow \quad A V_m \mathbf{y}_m = \mathbf{r}_0. \]
将 \(A V_m = V_{m+1} \bar{H}_m\) 和 \(\mathbf{r}_0 = \|\mathbf{r}_0\|_2 \mathbf{v}_1 = \|\mathbf{r}_0\|_2 V_{m+1} \mathbf{e}_1\)（其中 \(\mathbf{e}_1 = [1, 0, \dots, 0]^T \in \mathbb{R}^{m+1}\)）代入，得到：
\[V_{m+1} \bar{H}_m \mathbf{y}_m = \|\mathbf{r}_0\|_2 V_{m+1} \mathbf{e}_1. \]
由于 \(V_{m+1}\) 的列正交，可左乘 \(V_{m+1}^T\)，简化成：
\[\bar{H}_m \mathbf{y}_m = \|\mathbf{r}_0\|_2 \mathbf{e}_1. \]
注意 \(\bar{H}_m\) 是 \((m+1) \times m\) 矩阵，这是一个超定方程组。FOM的关键在于强制残差正交于Krylov子空间，即要求残差 \(\mathbf{r}_m = \mathbf{b} - A \mathbf{x}_m\) 满足 \(\mathbf{r}_m \perp \mathcal{K}_m(A, \mathbf{r}_0)\)。由于 \(\mathbf{r}_m = \mathbf{r}_0 - A V_m \mathbf{y}_m = V_{m+1} (\|\mathbf{r}_0\|_2 \mathbf{e}_1 - \bar{H}_m \mathbf{y}_m)\)，正交条件等价于残差在子空间基向量上的投影为零，即：
\[V_m^T \mathbf{r}_m = \mathbf{0} \quad \Rightarrow \quad V_m^T V_{m+1} (\|\mathbf{r}_0\|_2 \mathbf{e}_1 - \bar{H}_m \mathbf{y}_m) = \mathbf{0}. \]
由于 \(V_m^T V_{m+1} = [I_m \ | \ \mathbf{0}]\)（前 \(m\) 行单位矩阵，最后一列零向量），上述条件简化为：
\[H_m \mathbf{y}_m = \|\mathbf{r}_0\|_2 \mathbf{e}_1, \]
其中 \(H_m \in \mathbb{R}^{m \times m}\) 是 \(\bar{H}_m\) 的前 \(m\) 行（即去掉最后一行的上Hessenberg矩阵）。
步骤3：求解投影方程组
现在我们需要求解一个 \(m \times m\) 的上Hessenberg线性方程组：
\[H_m \mathbf{y}_m = \|\mathbf{r}_0\|_2 \mathbf{e}_1. \]
由于 \(m\) 很小，我们可以用直接法（如QR分解或高斯消元）高效求解。注意：
- 如果 \(H_m\) 是奇异的，说明找到了精确解（在Krylov子空间中）或方法失效，但在实际中可通过重启策略处理。
- 解出 \(\mathbf{y}_m\) 后，近似解为：
\[ \mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m. \]
步骤4：计算残差与收敛判断
FOM的残差范数可以通过递推公式计算，无需显式计算 \(A \mathbf{x}_m\)。实际上，残差满足：
\[\mathbf{r}_m = \mathbf{b} - A \mathbf{x}_m = -h_{m+1,m} (\mathbf{e}_m^T \mathbf{y}_m) \mathbf{v}_{m+1}, \]
其中 \(\mathbf{e}_m\) 是第 \(m\) 个单位向量，\(h_{m+1,m}\) 来自 \(\bar{H}_m\) 的最后一行的次对角线元素。
因此，残差范数为：
\[\|\mathbf{r}_m\|_2 = |h_{m+1,m}| \cdot |y_{m}^{(m)}|, \]
这里 \(y_{m}^{(m)}\) 是 \(\mathbf{y}_m\) 的最后一个分量。
检查 \(\|\mathbf{r}_m\|_2\) 是否小于预设的容差 \(\epsilon\)。如果满足，则输出 \(\mathbf{x}_m\) 作为近似解；否则，可以增加 \(m\) 或采用重启策略（即用当前 \(\mathbf{x}_m\) 作为新的初始猜测，重新运行FOM）。
总结
FOM算法通过Arnoldi过程构造正交基，将大规模问题投影到一个小型的上Hessenberg系统，从而高效求解。与GMRES（它最小化残差范数）不同，FOM强制残差正交于Krylov子空间，这使得它在每一步的残差范数不一定单调下降，但通常收敛性与GMRES相似。FOM适合与预处理技术结合，以加速收敛。