好的，根据你的要求，我们从“线性代数算法”领域中随机选取一个尚未讲解的算法题目。
题目：基于Krylov子空间方法的FOM（全正交化方法）算法解线性方程组
算法描述
我们要求解一个大型稀疏的非对称线性方程组：
\[A\mathbf{x} = \mathbf{b} \]
其中 \(A \in \mathbb{R}^{n \times n}\) 是一个大型稀疏矩阵（意味着 \(n\) 很大，且 \(A\) 中大部分元素为零），\(\mathbf{b} \in \mathbb{R}^n\) 是已知的右端向量，\(\mathbf{x}\) 是待求的未知向量。
对于此类问题，直接使用高斯消元法等直接法通常计算量和存储量巨大，不切实际。Krylov子空间迭代法是主流选择，它们通过在一个精心构造的低维子空间中寻找近似解来逼近真实解。
FOM (Full Orthogonalization Method)，即全正交化方法，是Arnoldi过程在线性方程组求解中的一个经典应用。其核心思想是：在由矩阵 \(A\) 和初始残差向量 \(\mathbf{r}_0\) 生成的Krylov子空间中，通过强制Galerkin条件（即让新的残差向量与当前的Krylov子空间正交）来构造近似解。
解题过程
FOM的求解过程可以清晰地分为几个步骤。我们假设从初始猜测解 \(\mathbf{x}_0\) 开始（通常可以取 \(\mathbf{x}_0 = \mathbf{0}\)），目标是构造一个在 \(m\) 维Krylov子空间中的近似解 \(\mathbf{x}_m\)，使得精度满足要求。
步骤 1：初始化与Arnoldi过程
- 计算初始残差：
\[ \mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0 \]
初始近似解 $ \mathbf{x}_0 $ 通常取零向量，此时 $ \mathbf{r}_0 = \mathbf{b} $。
-
初始化Arnoldi过程：
- 令 \(\beta = \|\mathbf{r}_0\|_2\) （2-范数，即向量长度）。
- 取第一个标准正交基向量：\(\mathbf{v}_1 = \mathbf{r}_0 / \beta\)。
Arnoldi过程的目标是为第 \(m\) 个Krylov子空间 \(\mathcal{K}_m(A, \mathbf{r}_0) = \text{span}\{\mathbf{r}_0, A\mathbf{r}_0, ..., A^{m-1}\mathbf{r}_0\}\) 构造一组标准正交基 \(\{\mathbf{v}_1, \mathbf{v}_2, ..., \mathbf{v}_m\}\)，并得到一个上Hessenberg矩阵 \(H_m\)。
-
执行标准Arnoldi迭代（对于 j=1, 2, ..., m）：
a. 矩阵向量乘：\(\mathbf{w} = A \mathbf{v}_j\)。
b. 正交化：对 \(i = 1\) 到 \(j\) 执行，计算内积 \(h_{i,j} = \langle \mathbf{v}_i, \mathbf{w} \rangle\)，并从 \(\mathbf{w}\) 中减去 \(\mathbf{v}_i\) 方向的投影：
\[ \mathbf{w} := \mathbf{w} - h_{i,j} \mathbf{v}_i \]
这一步保证了新向量 $ \mathbf{w} $ 与所有已生成的基向量 $ \mathbf{v}_i $ 正交。
c. **标准化**：计算 $ h_{j+1,j} = \|\mathbf{w}\|_2 $。如果 $ h_{j+1,j} = 0 $，则说明Krylov子空间已“崩溃”，过程可以提前终止（此时已找到精确解所在的子空间）。否则，令：
\[ \mathbf{v}_{j+1} = \mathbf{w} / h_{j+1,j} \]
这就得到了一个新的正交基向量。
经过 $ m $ 步后，我们得到了正交矩阵 $ V_m = [\mathbf{v}_1, \mathbf{v}_2, ..., \mathbf{v}_m] \in \mathbb{R}^{n \times m} $，以及一个上Hessenberg矩阵 $ H_m \in \mathbb{R}^{m \times m} $（其元素就是上面计算的 $ h_{i,j} $）。它们满足一个至关重要的**Arnoldi关系式**：
\[ A V_m = V_m H_m + h_{m+1, m} \mathbf{v}_{m+1} \mathbf{e}_m^T \]
其中 $ \mathbf{e}_m $ 是 $ m $ 维单位矩阵的第 $ m $ 列。
步骤 2：在Krylov子空间中构造近似解
FOM的核心是Galerkin条件：寻找近似解 \(\mathbf{x}_m \in \mathbf{x}_0 + \mathcal{K}_m(A, \mathbf{r}_0)\)，使得新的残差 \(\mathbf{r}_m = \mathbf{b} - A\mathbf{x}_m\) 与当前的Krylov子空间 \(\mathcal{K}_m(A, \mathbf{r}_0)\) 正交。即：
\[\mathbf{r}_m \perp \mathcal{K}_m(A, \mathbf{r}_0) \quad \text{或等价地} \quad V_m^T \mathbf{r}_m = \mathbf{0} \]
- 解的表示：
由于 \(\mathbf{x}_m \in \mathbf{x}_0 + \mathcal{K}_m\)，它可以表示为：
\[ \mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m \]
其中 $ \mathbf{y}_m \in \mathbb{R}^m $ 是一个待求的系数向量。
- 利用Galerkin条件求解 \(\mathbf{y}_m\)：
计算残差：
\[ \mathbf{r}_m = \mathbf{b} - A\mathbf{x}_m = \mathbf{r}_0 - A V_m \mathbf{y}_m \]
代入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 = V_m (\beta \mathbf{e}_1) $（其中 $ \mathbf{e}_1 $ 是 $ m $ 维单位矩阵的第一列），我们得到：
\[ \mathbf{r}_m = V_m (\beta \mathbf{e}_1 - H_m \mathbf{y}_m) - h_{m+1,m} (\mathbf{e}_m^T \mathbf{y}_m) \mathbf{v}_{m+1} \]
现在，将 $ \mathbf{r}_m $ 与 $ V_m $ 作内积（即要求 $ V_m^T \mathbf{r}_m = \mathbf{0} $）。由于 $ V_m^T V_m = I $ 且 $ V_m^T \mathbf{v}_{m+1} = \mathbf{0} $（因为基向量相互正交），Galerkin条件简化为：
\[ \beta \mathbf{e}_1 - H_m \mathbf{y}_m = \mathbf{0} \]
即：
\[ H_m \mathbf{y}_m = \beta \mathbf{e}_1 \]
-
求解小型稠密线性系统：
现在，我们只需要求解一个 \(m \times m\) 的上Hessenberg线性方程组 \(H_m \mathbf{y}_m = \beta \mathbf{e}_1\)。由于 \(m\) 通常远小于 \(n\)（比如几十到几百），这是一个小型稠密问题，可以用高斯消元法（或QR分解）稳定高效地求解。 -
得到近似解：
解出 \(\mathbf{y}_m\) 后，代回公式：
\[ \mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m \]
就得到了第 $ m $ 步的近似解。
步骤 3：检查收敛性与重启策略
- 计算残差范数以检查收敛：
我们并不需要显式计算 \(\mathbf{r}_m\) 来检查收敛。从 \(\mathbf{r}_m\) 的表达式可以看出，因为 \(V_m^T \mathbf{r}_m = \mathbf{0}\)，残差 \(\mathbf{r}_m\) 平行于 \(\mathbf{v}_{m+1}\)。因此，其2-范数为：
\[ \|\mathbf{r}_m\|_2 = |h_{m+1, m} (\mathbf{e}_m^T \mathbf{y}_m)| \]
在求解小型方程组 $ H_m \mathbf{y}_m = \beta \mathbf{e}_1 $ 后，我们可以轻易计算出 $ \mathbf{y}_m $ 的最后一个分量 $ \eta_m = \mathbf{e}_m^T \mathbf{y}_m $，从而得到残差范数 $ \| \mathbf{r}_m \|_2 = |h_{m+1, m} \eta_m| $。
- 判断与重启：
- 如果 \(\|\mathbf{r}_m\|_2\) 小于预设的容差（Tolerance），则算法收敛，输出 \(\mathbf{x}_m\) 作为最终解。
- 如果未收敛，且子空间维数 \(m\) 已达到预设的最大值（例如50或100），则需要进行重启（Restart）：
- 以当前的 \(\mathbf{x}_m\) 作为新的初始猜测 \(\mathbf{x}_0^{\text{new}}\)。
- 计算新的初始残差 \(\mathbf{r}_0^{\text{new}} = \mathbf{b} - A\mathbf{x}_0^{\text{new}}\)。
- 清空基向量集合，从步骤1开始，基于新的 \(\mathbf{r}_0^{\text{new}}\) 重新构建Krylov子空间并迭代。
- 重启是为了控制内存开销（存储所有基向量 \(V_m\)）和计算量（正交化成本随 \(m\) 增加而增长）。通过多次重启，FOM可以逐渐逼近真实解。
总结与要点
- 核心思想：在由 \(A\) 和初始残差生成的Krylov子空间中，通过强制Galerkin正交条件来寻找最佳近似解。
- 关键技术：使用Arnoldi过程构造Krylov子空间的标准正交基，并得到小型稠密的Hessenberg矩阵 \(H_m\)。
- 核心计算：将一个庞大的 \(n\) 维问题，约化为一个很小的 \(m\) 维稠密线性方程组 \(H_m \mathbf{y}_m = \beta \mathbf{e}_1\) 来求解。
- 优点：对于非对称矩阵，能产生一组正交基，数值稳定性较好。
- 缺点：
- 需要存储所有基向量 \(V_m\)，内存开销为 \(O(n \times m)\)。
- 正交化过程（步骤3.b）的计算成本随 \(m\) 呈二次方增长（\(O(m^2 n)\)）。
- 与GMRES的关系：FOM是GMRES（广义最小残差法）的“兄弟”算法。GMRES在同样的Krylov子空间中，采用Petrov-Galerkin条件（即最小化残差2-范数），这导致它需要求解一个最小二乘问题 \(\min \|\beta \mathbf{e}_1 - H_m \mathbf{y}_m\|_2\)，而非线性方程组。GMRES通常比FOM更鲁棒，但FOM在某些情况下（当 \(H_m\) 非奇异时）计算量略小。
通过以上步骤，FOM算法将求解大规模稀疏非对称线性方程组的问题，系统地转化为一系列小规模稠密问题的求解，并通过迭代和重启策略逼近最终答案。