Krylov子空间方法在特征值问题中的应用
好的，我们将聚焦于Krylov子空间方法在线性代数特征值问题中的运用。我将从最基础的概念开始，逐步构建起对该方法的完整理解。
步骤 1：问题定义与目标
我们要解决的标准特征值问题是：
对于给定的一个 \(n \times n\) 矩阵 \(A\)，寻找标量 \(\lambda\)（特征值）和非零向量 \(\mathbf{v}\)（特征向量），使得：
\[A \mathbf{v} = \lambda \mathbf{v} \]
我们的目标通常是计算矩阵 \(A\) 的一部分特征值和对应的特征向量，例如模最大的几个、模最小的几个、或者实部在某个区间内的特征值。对于大型稀疏矩阵（即矩阵中绝大多数元素为零），直接使用QR算法等稠密矩阵方法是不可行的，因为其计算量和存储需求过高（\(O(n^3)\) 计算量和 \(O(n^2)\) 存储）。Krylov子空间方法正是为此类大规模稀疏特征值问题设计的迭代算法。
步骤 2：Krylov子空间的定义
Krylov子空间方法是基于一个核心的数学结构——Krylov子空间。
给定一个矩阵 \(A\) 和一个非零初始向量 \(\mathbf{q}_1\)（通常经过标准化，即 \(\|\mathbf{q}_1\| = 1\)），对应的 \(m\) 维 Krylov 子空间定义为：
\[\mathcal{K}_m(A, \mathbf{q}_1) = \text{span} \{ \mathbf{q}_1, A\mathbf{q}_1, A^2\mathbf{q}_1, \dots, A^{m-1}\mathbf{q}_1 \} \]
其中 \(m \ll n\)。这个空间由初始向量 \(\mathbf{q}_1\) 经过反复左乘矩阵 \(A\) 所张成。直观上，\(A\) 对这个空间的作用会“暴露”出与其“大”特征值相关的方向信息。
步骤 3：构建正交基——Arnoldi过程
为了在 Krylov 子空间中进行数值稳定的计算，我们需要为其构造一组标准正交基（正交归一化基）。这通过 Arnoldi迭代 过程实现，它本质上是一种改进的 Gram-Schmidt 正交化过程。
算法过程如下：
- 选择初始向量 \(\mathbf{q}_1\)，满足 \(\|\mathbf{q}_1\| = 1\)。
- 对于 \(j = 1, 2, \dots, m\) 执行：
a. 计算 \(\mathbf{v} = A \mathbf{q}_j\)。
b. 对于 \(i = 1\) 到 \(j\)，计算投影系数 \(h_{ij} = \mathbf{q}_i^T \mathbf{v}\)，并从 \(\mathbf{v}\) 中减去这个投影：\(\mathbf{v} = \mathbf{v} - h_{ij} \mathbf{q}_i\)。
c. 令 \(h_{j+1, j} = \|\mathbf{v}\|\)。
d. 如果 \(h_{j+1, j} = 0\) 则停止，否则令 \(\mathbf{q}_{j+1} = \mathbf{v} / h_{j+1, j}\)。
经过 \(m\) 步后，我们得到了一组标准正交向量 \(\mathbf{q}_1, \dots, \mathbf{q}_{m+1}\)，以及一个 \((m+1) \times m\) 的上 Hessenberg 矩阵 \(\tilde{H}_m\)（即除了主对角线和第一条次对角线外，其他元素为零）。它们满足一个关键关系式：
\[A Q_m = Q_{m+1} \tilde{H}_m \]
其中，\(Q_m = [\mathbf{q}_1, \dots, \mathbf{q}_m]\) 是 \(n \times m\) 矩阵，其列向量构成了 \(\mathcal{K}_m\) 的一组标准正交基；\(Q_{m+1}\) 则在 \(Q_m\) 基础上增加了一列 \(\mathbf{q}_{m+1}\)。如果我们只取前 \(m\) 行，则有：
\[A Q_m = Q_m H_m + h_{m+1, m} \mathbf{q}_{m+1} \mathbf{e}_m^T \]
这里 \(H_m\) 是 \(\tilde{H}_m\) 的前 \(m\) 行，一个 \(m \times m\) 的上 Hessenberg 矩阵。
步骤 4：将大规模问题投影为小规模问题
这是 Krylov 子空间方法的核心思想。我们寻找原问题 \(A\mathbf{v} = \lambda \mathbf{v}\) 在 Krylov 子空间 \(\mathcal{K}_m\) 中的“近似解”。具体地，我们寻找近似特征对 \((\theta, \mathbf{y})\)，其中 \(\theta\) 是标量，\(\mathbf{y}\) 是 \(m\) 维向量，使得近似特征向量 \(\mathbf{x} = Q_m \mathbf{y} \in \mathcal{K}_m\) 满足某种最佳条件。
常用的准则被称为 Galerkin条件：要求残差 \(\mathbf{r} = A \mathbf{x} - \theta \mathbf{x}\) 与 Krylov 子空间 \(\mathcal{K}_m\) 正交。
即：
\[\mathbf{r} \perp \mathcal{K}_m \quad \Leftrightarrow \quad Q_m^T \mathbf{r} = \mathbf{0} \]
将 \(\mathbf{x} = Q_m \mathbf{y}\) 和 \(A Q_m\) 的关系代入，我们得到：
\[Q_m^T (A Q_m \mathbf{y} - \theta Q_m \mathbf{y}) = (H_m \mathbf{y} - \theta \mathbf{y}) = \mathbf{0} \]
因此，投影特征值问题 是求解一个小规模的 \(m \times m\) 矩阵 \(H_m\) 的特征值问题：
\[H_m \mathbf{y} = \theta \mathbf{y} \]
这个小问题可以通过 QR 算法等稠密方法高效求解。解得的 \((\theta_i, \mathbf{y}_i)\) 被称为 Ritz对，其中 \(\theta_i\) 是 Ritz 值，\(\mathbf{x}_i = Q_m \mathbf{y}_i\) 是 Ritz 向量，它们是原问题在大 Krylov 子空间中的近似特征对。
步骤 5：重启技术与收敛性
随着子空间维数 \(m\) 的增大，存储 \(Q_m\) 和计算 \(H_m\) 特征值的开销会增加。但通常，远在 \(m\) 接近 \(n\) 之前，特征值近似值 \(\theta_i\) 中就会有几个已经收敛到 \(A\) 的某些极端特征值（如模最大或最小的特征值）。
为了计算更多特征值或加速收敛，我们需要 重启技术。基本思想是：
- 运行 Arnoldi 过程到指定的维数 \(m\)。
- 计算 \(H_m\) 的特征值，并选取 \(k (k < m)\) 个“想要的”Ritz 对（例如，实部最大的几个）。
- 以这 \(k\) 个 Ritz 向量的线性组合（或其中“最好”的一个）作为新的初始向量 \(\mathbf{q}_1\)，重新开始 Arnoldi 过程。
这样，我们不断地“刷新”Krylov 子空间，使其更聚焦于我们感兴趣的特征方向，同时控制了子空间的规模。
步骤 6：重要变体——Lanczos算法（针对对称矩阵）
当矩阵 \(A\) 是实对称或复 Hermitian 矩阵时，Arnoldi 过程会大幅简化。由于对称性，正交化过程中 \(\mathbf{v}\) 只需要与最近的两个基向量正交即可。这个简化的算法称为 Lanczos迭代。
它生成的三项递推关系为：
\[A Q_m = Q_m T_m + \beta_m \mathbf{q}_{m+1} \mathbf{e}_m^T \]
其中 \(T_m\) 是一个实对称的三对角矩阵（只有主对角线和相邻的两条次对角线非零）。这使得投影特征值问题 \(T_m \mathbf{y} = \theta \mathbf{y}\) 的求解和收敛性分析变得更为简单高效。Lanczos 方法是计算大型对称矩阵特征值问题最核心的工具之一。
总结
Krylov子空间方法解决特征值问题的精髓在于：
- 构建空间：通过矩阵 \(A\) 与一个初始向量的幂作用，构建出能反映 \(A\) 主要谱信息的 Krylov 子空间。
- 正交化：通过 Arnoldi（或对称情形的 Lanczos）过程，获得子空间的标准正交基，并得到一个投影矩阵 \(H_m\)（或 \(T_m\)）。
- 投影与求解：将大规模的原特征值问题，投影为一个小规模的、基于 \(H_m\) 的特征值问题，并用稠密方法求解，得到 Ritz 对作为近似解。
- 重启迭代：通过重启技术，控制内存消耗，并引导迭代过程收敛到我们感兴趣的特征对。
这种方法在科学计算中应用极广，是求解大规模稀疏矩阵部分特征值的基石算法。