好的,已经仔细核对过您提供的列表,我接下来为您讲解一个在已列题目中尚未出现的重要线性代数迭代算法。
题目:基于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\) 是未知解向量。由于矩阵规模巨大,直接解法(如高斯消元)计算量和存储需求都过高,因此我们寻求一种基于投影的迭代方法,在Krylov子空间中寻找近似解。
FOM (Full Orthogonalization Method) 正是这样一类方法。它的核心思想是:在第 \(m\) 步迭代 (\(m \ll n\)),在由 \(A\) 和初始残差 \(\mathbf{r}_0\) 生成的 \(m\) 维Krylov子空间 \(\mathcal{K}_m(A, \mathbf{r}_0)\) 中,寻找一个解 \(\mathbf{x}_m\),使得新的残差 \(\mathbf{b} - A\mathbf{x}_m\) 与这个Krylov子空间正交。
循序渐进解题过程:
我们从初始猜测解 \(\mathbf{x}_0\) 开始。通常,为简单起见,可以设 \(\mathbf{x}_0 = \mathbf{0}\)。
步骤 1:初始化与Krylov子空间的生成
- 计算初始残差:
\[ \mathbf{r}_0 = \mathbf{b} - A \mathbf{x}_0 \]
令 $ \beta = \|\mathbf{r}_0\|_2 $ (2-范数)。如果 $ \beta = 0 $,则 $ \mathbf{x}_0 $ 已是精确解,算法终止。
-
构建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{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_m\}\) 来张成这个子空间。这通过 Arnoldi 过程 实现。- 令第一个基向量: \(\mathbf{v}_1 = \mathbf{r}_0 / \beta\)。
- 对于 \(j = 1, 2, \dots, m\):
a. 计算向量 \(\mathbf{w}_j = A \mathbf{v}_j\)。
b. 对于 \(i = 1\) 到 \(j\),执行 完全正交化:
\[ h_{i,j} = (\mathbf{w}_j, \mathbf{v}_i) \quad \text{(向量内积)} \]
\[ \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 $,则子空间是 $ A $-不变的,过程可以提前终止。
e. 否则,令 $ \mathbf{v}_{j+1} = \mathbf{w}_j / h_{j+1,j} $ 作为下一个标准正交基向量。
这个过程结束后,我们得到了一个 $ n \times m $ 的矩阵 $ V_m = [\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_m] $,其列是Krylov子空间的标准正交基,以及一个 $ (m+1) \times m $ 的上Hessenberg矩阵 $ \bar{H}_m $,其非零元素就是上面计算的 $ h_{i,j} $。
步骤 2:构建投影方程(Galerkin条件)
FOM方法施加的 Galerkin正交条件 是:
\[\mathbf{r}_m = \mathbf{b} - A\mathbf{x}_m \perp \mathcal{K}_m(A, \mathbf{r}_0) \]
由于 \(\mathcal{K}_m\) 的基是 \(V_m\),这个条件等价于:
\[V_m^T (\mathbf{b} - A\mathbf{x}_m) = \mathbf{0} \]
步骤 3:在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}_0 + V_m \mathbf{y}_m) = \mathbf{r}_0 - A V_m \mathbf{y}_m \]
步骤 4:利用Arnoldi关系简化
从Arnoldi过程,我们有一个关键关系:
\[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 \]
其中:
- \(\bar{H}_m\) 是 \((m+1) \times m\) 的Hessenberg矩阵。
- \(H_m\) 是 \(\bar{H}_m\) 去掉最后一行得到的 \(m \times m\) 上Hessenberg矩阵。
- \(\mathbf{e}_m\) 是 \(\mathbb{R}^m\) 中的第 \(m\) 个单位基向量。
同时,由于 \(\mathbf{v}_1 = \mathbf{r}_0 / \beta\),我们有 \(\mathbf{r}_0 = \beta \mathbf{v}_1 = V_{m+1} (\beta \mathbf{e}_1)\),其中 \(\mathbf{e}_1\) 是 \(\mathbb{R}^{m+1}\) 中的第一个单位基向量。
将这两个关系代入 \(\mathbf{r}_m\):
\[\mathbf{r}_m = V_{m+1} (\beta \mathbf{e}_1) - V_{m+1} \bar{H}_m \mathbf{y}_m = V_{m+1} (\beta \mathbf{e}_1 - \bar{H}_m \mathbf{y}_m) \]
步骤 5:应用Galerkin条件求解系数向量 \(\mathbf{y}_m\)
Galerkin条件 \(V_m^T \mathbf{r}_m = \mathbf{0}\) 意味着残差的前 \(m\) 个分量(在 \(V_m\) 坐标系下)为零。由于 \(V_m^T V_{m+1} = [I_m, \mathbf{0}]\),我们将条件应用到上面 \(\mathbf{r}_m\) 的表达式:
\[V_m^T \mathbf{r}_m = [I_m, \mathbf{0}] (\beta \mathbf{e}_1 - \bar{H}_m \mathbf{y}_m) = \beta \mathbf{e}_1' - H_m \mathbf{y}_m = \mathbf{0} \]
这里 \(\mathbf{e}_1'\) 是 \(\mathbb{R}^m\) 中的第一个单位基向量。于是我们得到了一个 小型稠密线性方程组:
\[H_m \mathbf{y}_m = \beta \mathbf{e}_1' \]
由于 \(m\) 很小(例如几十到几百),我们可以用直接法(如QR分解或Givens旋转)稳定地求解这个 \(m \times m\) 的上Hessenberg方程组。
步骤 6:形成近似解
求解出 \(\mathbf{y}_m\) 后,我们得到第 \(m\) 步的近似解:
\[\mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m \]
我们无需显式计算残差范数来判断收敛。根据 Arnoldi 关系,残差向量为 \(\mathbf{r}_m = -h_{m+1, m} (\mathbf{e}_m^T \mathbf{y}_m) \mathbf{v}_{m+1}\),因此其范数可以方便地计算为:
\[\| \mathbf{r}_m \|_2 = |h_{m+1, m} (\mathbf{e}_m^T \mathbf{y}_m)| \]
如果这个值小于预设的收敛容差,我们就接受 \(\mathbf{x}_m\) 作为近似解。否则,可以增加 \(m\)(重启或继续扩展子空间)并重复上述过程。
总结:
FOM算法通过在Krylov子空间中施加Galerkin条件,将大规模稀疏问题简化为一个小型稠密问题(求解 \(H_m \mathbf{y}_m = \beta \mathbf{e}_1\))。它需要存储所有 \(m\) 个基向量 \(\mathbf{v}_j\)(内存开销为 \(O(nm)\)),并在每一步进行完全正交化(计算开销为 \(O(nm^2)\))。虽然计算成本比不完全正交化的GMRES高,但它在每个子空间维数 \(m\) 内提供了最优的(在2-范数意义下)残差,并且解 \(\mathbf{x}_m\) 位于一个真正的Krylov子空间中。为了处理大 \(m\) 时的内存和计算成本,通常采用重启策略(FOM(m)),即每进行 \(m\) 步迭代,就用当前近似解作为新的初始猜测重新开始整个过程。