基于Krylov子空间的FOM（全正交化方法）算法解线性方程组
字数 4261 2025-12-22 21:06:00

基于Krylov子空间的FOM（全正交化方法）算法解线性方程组

题目描述

给定大型稀疏非对称线性方程组 \(Ax = b\)，其中 \(A \in \mathbb{R}^{n \times n}\) 非奇异，\(b \in \mathbb{R}^n\)，且 \(n\) 很大，矩阵 \(A\) 稀疏。FOM（Full Orthogonalization Method，全正交化方法）是一种基于Krylov子空间的投影方法，通过在Krylov子空间上构造正交基，并将原方程投影到该子空间求解近似解。FOM的特点是每次迭代在子空间上求解一个小型稠密线性系统（通常为上Hessenberg系统），适用于非对称矩阵，是GMRES的前身之一。题目要求详细讲解FOM算法的推导、步骤、实现细节及关键性质。

解题过程循序渐进讲解

步骤1：理解投影方法的基本思想

对于线性方程组 \(Ax = b\)，投影方法的核心是：

  1. 选择一个维度为 \(m\) 的子空间 \(\mathcal{K}_m\)（搜索空间）。
  2. 选择一个约束子空间 \(\mathcal{L}_m\)（通常与搜索空间相同或不同）。
  3. 寻找近似解 \(x_m \in x_0 + \mathcal{K}_m\) 满足 Petrov–Galerkin 条件：

\[ b - A x_m \perp \mathcal{L}_m \]

如果 \(\mathcal{L}_m = \mathcal{K}_m\)，称为正交投影（Galerkin投影）；如果 \(\mathcal{L}_m = A\mathcal{K}_m\)，称为斜投影（如GMRES使用残差极小化）。

FOM采用 Galerkin 投影，即 \(\mathcal{L}_m = \mathcal{K}_m\)。

步骤2：构造Krylov子空间

给定初始猜测解 \(x_0\)，计算初始残差 \(r_0 = b - A x_0\)。
定义 \(m\) 维 Krylov 子空间：

\[\mathcal{K}_m(A, r_0) = \text{span}\{ r_0, A r_0, A^2 r_0, \dots, A^{m-1} r_0 \}. \]

FOM的搜索空间为 \(x_0 + \mathcal{K}_m(A, r_0)\)。

步骤3：Arnoldi过程生成正交基

为了数值稳定，用Arnoldi过程构造 \(\mathcal{K}_m\) 的一组标准正交基 \(\{v_1, v_2, \dots, v_m\}\)，其中 \(v_1 = r_0 / \|r_0\|_2\)。
Arnoldi过程递推生成：

\[A V_m = V_{m+1} \tilde{H}_m, \]

其中 \(V_m = [v_1, \dots, v_m] \in \mathbb{R}^{n \times m}\) 的列正交，\(\tilde{H}_m \in \mathbb{R}^{(m+1) \times m}\) 是上Hessenberg矩阵（几乎上三角，但多一次对角线）。
具体递推（对于 \(j = 1, \dots, m\)）：

  1. \(w = A v_j\)。
  2. 对 \(i = 1, \dots, j\)：\(h_{i,j} = v_i^T w\)，\(w = w - h_{i,j} v_i\)。
  3. \(h_{j+1,j} = \|w\|_2\)。
  4. 如果 \(h_{j+1,j} = 0\)，则子空间不变，提前终止；否则 \(v_{j+1} = w / h_{j+1,j}\)。

记 \(H_m\) 为 \(\tilde{H}_m\) 的前 \(m\) 行（即 \(m \times m\) 上Hessenberg矩阵），则有：

\[A V_m = V_m H_m + h_{m+1,m} v_{m+1} e_m^T, \]

其中 \(e_m\) 是第 \(m\) 个单位坐标向量。

步骤4：在Krylov子空间中表示近似解

设近似解为 \(x_m = x_0 + V_m y_m\)，其中 \(y_m \in \mathbb{R}^m\) 是待求系数向量。
残差向量：

\[r_m = b - A x_m = r_0 - A V_m y_m. \]

代入 Arnoldi 关系：\(A V_m = V_m H_m + h_{m+1,m} v_{m+1} e_m^T\)，得：

\[r_m = r_0 - V_m H_m y_m - h_{m+1,m} (e_m^T y_m) v_{m+1}. \]

注意 \(r_0 = \|r_0\|_2 v_1 = \|r_0\|_2 V_m e_1\)（因为 \(v_1\) 是 \(V_m\) 的第一列），所以：

\[r_m = V_m \big( \|r_0\|_2 e_1 - H_m y_m \big) - h_{m+1,m} (e_m^T y_m) v_{m+1}. \]

步骤5：施加Galerkin条件

Galerkin条件要求残差与子空间 \(\mathcal{K}_m\) 正交，即 \(V_m^T r_m = 0\)。
由于 \(V_m^T V_m = I_m\)，且 \(V_m^T v_{m+1} = 0\)，代入上式得：

\[V_m^T r_m = \|r_0\|_2 e_1 - H_m y_m = 0. \]

因此，系数向量 \(y_m\) 满足小型上Hessenberg线性系统：

\[H_m y_m = \|r_0\|_2 e_1. \]

步骤6：求解小型系统得到近似解

因为 \(H_m\) 是 \(m \times m\) 上Hessenberg矩阵，且通常 \(m \ll n\)，可用稠密矩阵解法（如QR分解、Givens旋转）稳定求解 \(y_m\)。
一旦得到 \(y_m\)，近似解为：

\[x_m = x_0 + V_m y_m. \]

实际计算中不需要存储所有 \(v_j\) 向量（除非需要后续迭代），但若只求最终解，只需在 Arnoldi 过程中同时更新 \(x_m\) 的递推形式。

步骤7：残差范数的简便计算

Galerkin条件不保证残差范数单调下降，但残差范数可以从解 \(y_m\) 方便计算。由前面表达式：

\[r_m = - h_{m+1,m} (e_m^T y_m) v_{m+1}. \]

因此残差范数为：

\[\|r_m\|_2 = |h_{m+1,m}| \cdot |e_m^T y_m|. \]

这提供了停机准则：若 \(\|r_m\|_2 < \text{tol}\)，则终止迭代。

步骤8：算法流程总结

  1. 输入 \(A, b, x_0\)，容差 \(\epsilon\)，最大迭代步数 \(m_{\max}\)。
  2. 计算 \(r_0 = b - A x_0\)，\(\beta = \|r_0\|_2\)，若 \(\beta < \epsilon\) 则输出 \(x_0\)。
  3. 设 \(v_1 = r_0 / \beta\)。
  4. 对 \(j = 1, 2, \dots, m_{\max}\) 做 Arnoldi 过程：
    • 计算 \(w = A v_j\)。
    • 对 \(i = 1, \dots, j\)：\(h_{i,j} = v_i^T w\)，\(w = w - h_{i,j} v_i\)。
    • \(h_{j+1,j} = \|w\|_2\)，如果 \(h_{j+1,j} = 0\) 则转到第6步。
    • \(v_{j+1} = w / h_{j+1,j}\)。
  5. 求解 \(H_j y_j = \beta e_1\) 得 \(y_j\)。
  6. 计算近似解 \(x_j = x_0 + V_j y_j\)。
  7. 计算残差范数 \(\|r_j\|_2 = |h_{j+1,j}| \cdot |e_j^T y_j|\)，若 \(\|r_j\|_2 < \epsilon\) 则终止并输出 \(x_j\)。
  8. 若未收敛且 \(j = m_{\max}\)，可重新以当前解为初值重启（类似GMRES）。

步骤9：关键性质与注意点

  • 收敛性：若 \(A\) 对称正定，FOM 等价于共轭梯度法（CG）；对非对称矩阵，FOM 不一定保证残差单调下降，但若 \(H_m\) 非奇异，解存在唯一。
  • 存储与计算成本：每步需存储所有基向量 \(v_1, \dots, v_{m+1}\)（\(O(n m)\) 内存），Arnoldi 内积计算量 \(O(m^2 n)\)。因此需限制 \(m\) 或采用重启策略（FOM(m)）。
  • 与GMRES关系：GMRES 在相同子空间上最小残差范数，需解最小二乘问题（\(\|\beta e_1 - \tilde{H}_m y\|_2\) 极小），而FOM解小型方阵系统，计算略简但收敛性可能不如GMRES稳定。
  • 中断处理：若 Arnoldi 过程中 \(h_{j+1,j} = 0\)，则 Krylov 子空间不变，此时 \(\mathcal{K}_j\) 是不变子空间，解精确。

总结

FOM 是经典的 Krylov 子空间投影方法，通过 Arnoldi 过程构造正交基，利用 Galerkin 条件导出小型上Hessenberg系统，迭代求解原方程。它适合非对称稀疏大系统，但需注意内存增长和可能的收敛振荡。实际应用中常与预处理技术结合加速收敛。

基于Krylov子空间的FOM（全正交化方法）算法解线性方程组 题目描述 给定大型稀疏非对称线性方程组 \( Ax = b \)，其中 \( A \in \mathbb{R}^{n \times n} \) 非奇异，\( b \in \mathbb{R}^n \)，且 \( n \) 很大，矩阵 \( A \) 稀疏。FOM（Full Orthogonalization Method，全正交化方法）是一种基于Krylov子空间的投影方法，通过在Krylov子空间上构造正交基，并将原方程投影到该子空间求解近似解。FOM的特点是每次迭代在子空间上求解一个小型稠密线性系统（通常为上Hessenberg系统），适用于非对称矩阵，是GMRES的前身之一。题目要求详细讲解FOM算法的推导、步骤、实现细节及关键性质。 解题过程循序渐进讲解 步骤1：理解投影方法的基本思想 对于线性方程组 \( Ax = b \)，投影方法的核心是： 选择一个维度为 \( m \) 的子空间 \( \mathcal{K}_ m \)（搜索空间）。 选择一个约束子空间 \( \mathcal{L}_ m \)（通常与搜索空间相同或不同）。 寻找近似解 \( x_ m \in x_ 0 + \mathcal{K}_ m \) 满足 Petrov–Galerkin 条件： \[ b - A x_ m \perp \mathcal{L}_ m \] 如果 \( \mathcal{L}_ m = \mathcal{K}_ m \)，称为正交投影（Galerkin投影）；如果 \( \mathcal{L}_ m = A\mathcal{K}_ m \)，称为斜投影（如GMRES使用残差极小化）。 FOM采用 Galerkin 投影，即 \( \mathcal{L}_ m = \mathcal{K}_ m \)。 步骤2：构造Krylov子空间 给定初始猜测解 \( x_ 0 \)，计算初始残差 \( r_ 0 = b - A x_ 0 \)。 定义 \( m \) 维 Krylov 子空间： \[ \mathcal{K}_ m(A, r_ 0) = \text{span}\{ r_ 0, A r_ 0, A^2 r_ 0, \dots, A^{m-1} r_ 0 \}. \] FOM的搜索空间为 \( x_ 0 + \mathcal{K}_ m(A, r_ 0) \)。 步骤3：Arnoldi过程生成正交基 为了数值稳定，用Arnoldi过程构造 \( \mathcal{K}_ m \) 的一组标准正交基 \( \{v_ 1, v_ 2, \dots, v_ m\} \)，其中 \( v_ 1 = r_ 0 / \|r_ 0\| 2 \)。 Arnoldi过程递推生成： \[ A V_ m = V {m+1} \tilde{H}_ m, \] 其中 \( V_ m = [ v_ 1, \dots, v_ m] \in \mathbb{R}^{n \times m} \) 的列正交，\( \tilde{H}_ m \in \mathbb{R}^{(m+1) \times m} \) 是上Hessenberg矩阵（几乎上三角，但多一次对角线）。 具体递推（对于 \( j = 1, \dots, m \)）： \( w = A v_ j \)。 对 \( i = 1, \dots, j \)：\( h_ {i,j} = v_ i^T w \)，\( w = w - h_ {i,j} v_ i \)。 \( h_ {j+1,j} = \|w\|_ 2 \)。 如果 \( h_ {j+1,j} = 0 \)，则子空间不变，提前终止；否则 \( v_ {j+1} = w / h_ {j+1,j} \)。 记 \( H_ m \) 为 \( \tilde{H} m \) 的前 \( m \) 行（即 \( m \times m \) 上Hessenberg矩阵），则有： \[ A V_ m = V_ m H_ m + h {m+1,m} v_ {m+1} e_ m^T, \] 其中 \( e_ m \) 是第 \( m \) 个单位坐标向量。 步骤4：在Krylov子空间中表示近似解 设近似解为 \( x_ m = x_ 0 + V_ m y_ m \)，其中 \( y_ m \in \mathbb{R}^m \) 是待求系数向量。 残差向量： \[ r_ m = b - A x_ m = r_ 0 - A V_ m y_ m. \] 代入 Arnoldi 关系：\( A V_ m = V_ m H_ m + h_ {m+1,m} v_ {m+1} e_ m^T \)，得： \[ r_ m = r_ 0 - V_ m H_ m y_ m - h_ {m+1,m} (e_ m^T y_ m) v_ {m+1}. \] 注意 \( r_ 0 = \|r_ 0\|_ 2 v_ 1 = \|r_ 0\| 2 V_ m e_ 1 \)（因为 \( v_ 1 \) 是 \( V_ m \) 的第一列），所以： \[ r_ m = V_ m \big( \|r_ 0\| 2 e_ 1 - H_ m y_ m \big) - h {m+1,m} (e_ m^T y_ m) v {m+1}. \] 步骤5：施加Galerkin条件 Galerkin条件要求残差与子空间 \( \mathcal{K} m \) 正交，即 \( V_ m^T r_ m = 0 \)。 由于 \( V_ m^T V_ m = I_ m \)，且 \( V_ m^T v {m+1} = 0 \)，代入上式得： \[ V_ m^T r_ m = \|r_ 0\|_ 2 e_ 1 - H_ m y_ m = 0. \] 因此，系数向量 \( y_ m \) 满足小型上Hessenberg线性系统： \[ H_ m y_ m = \|r_ 0\|_ 2 e_ 1. \] 步骤6：求解小型系统得到近似解 因为 \( H_ m \) 是 \( m \times m \) 上Hessenberg矩阵，且通常 \( m \ll n \)，可用稠密矩阵解法（如QR分解、Givens旋转）稳定求解 \( y_ m \)。 一旦得到 \( y_ m \)，近似解为： \[ x_ m = x_ 0 + V_ m y_ m. \] 实际计算中不需要存储所有 \( v_ j \) 向量（除非需要后续迭代），但若只求最终解，只需在 Arnoldi 过程中同时更新 \( x_ m \) 的递推形式。 步骤7：残差范数的简便计算 Galerkin条件不保证残差范数单调下降，但残差范数可以从解 \( y_ m \) 方便计算。由前面表达式： \[ r_ m = - h_ {m+1,m} (e_ m^T y_ m) v_ {m+1}. \] 因此残差范数为： \[ \|r_ m\| 2 = |h {m+1,m}| \cdot |e_ m^T y_ m|. \] 这提供了停机准则：若 \( \|r_ m\|_ 2 < \text{tol} \)，则终止迭代。 步骤8：算法流程总结 输入 \( A, b, x_ 0 \)，容差 \( \epsilon \)，最大迭代步数 \( m_ {\max} \)。 计算 \( r_ 0 = b - A x_ 0 \)，\( \beta = \|r_ 0\|_ 2 \)，若 \( \beta < \epsilon \) 则输出 \( x_ 0 \)。 设 \( v_ 1 = r_ 0 / \beta \)。 对 \( j = 1, 2, \dots, m_ {\max} \) 做 Arnoldi 过程： 计算 \( w = A v_ j \)。 对 \( i = 1, \dots, j \)：\( h_ {i,j} = v_ i^T w \)，\( w = w - h_ {i,j} v_ i \)。 \( h_ {j+1,j} = \|w\| 2 \)，如果 \( h {j+1,j} = 0 \) 则转到第6步。 \( v_ {j+1} = w / h_ {j+1,j} \)。 求解 \( H_ j y_ j = \beta e_ 1 \) 得 \( y_ j \)。 计算近似解 \( x_ j = x_ 0 + V_ j y_ j \)。 计算残差范数 \( \|r_ j\| 2 = |h {j+1,j}| \cdot |e_ j^T y_ j| \)，若 \( \|r_ j\|_ 2 < \epsilon \) 则终止并输出 \( x_ j \)。 若未收敛且 \( j = m_ {\max} \)，可重新以当前解为初值重启（类似GMRES）。 步骤9：关键性质与注意点 收敛性 ：若 \( A \) 对称正定，FOM 等价于共轭梯度法（CG）；对非对称矩阵，FOM 不一定保证残差单调下降，但若 \( H_ m \) 非奇异，解存在唯一。 存储与计算成本 ：每步需存储所有基向量 \( v_ 1, \dots, v_ {m+1} \)（\( O(n m) \) 内存），Arnoldi 内积计算量 \( O(m^2 n) \)。因此需限制 \( m \) 或采用重启策略（FOM(m)）。 与GMRES关系 ：GMRES 在相同子空间上最小残差范数，需解最小二乘问题（\( \|\beta e_ 1 - \tilde{H}_ m y\|_ 2 \) 极小），而FOM解小型方阵系统，计算略简但收敛性可能不如GMRES稳定。 中断处理 ：若 Arnoldi 过程中 \( h_ {j+1,j} = 0 \)，则 Krylov 子空间不变，此时 \( \mathcal{K}_ j \) 是不变子空间，解精确。 总结 FOM 是经典的 Krylov 子空间投影方法，通过 Arnoldi 过程构造正交基，利用 Galerkin 条件导出小型上Hessenberg系统，迭代求解原方程。它适合非对称稀疏大系统，但需注意内存增长和可能的收敛振荡。实际应用中常与预处理技术结合加速收敛。