非对称Lanczos双正交化算法在求解非对称矩阵部分特征对中的应用
我们来讲解一个在大型稀疏非对称矩阵特征值计算中非常重要的算法。这个算法是Arnoldi算法的非对称推广，但与Arnoldi算法构建正交基不同，它构建的是双正交基，从而将一个一般矩阵转化为三对角矩阵，大大降低了特征值问题的规模。
题目描述
给定一个大型稀疏的非对称矩阵 \(A \in \mathbb{C}^{n \times n}\)，我们希望计算其部分特征值（通常是模最大或最小的几个）及对应的特征向量。直接对A应用特征值分解（如QR算法）在n很大时计算代价过高。非对称Lanczos双正交化算法 通过构建两个Krylov子空间，将A投影到一个较小的三对角矩阵上，然后通过求解这个三对角矩阵的特征值问题来近似原矩阵的部分特征对。
核心问题：如何高效、稳定地为非对称矩阵构建双正交的Krylov子空间基，并将其转化为三对角形式，进而计算特征值近似解。
解题过程循序渐进讲解
步骤1：算法核心思想与目标
对于一个非对称矩阵A，其特征向量本身不一定正交。Lanczos双正交化算法的核心思想是同时构建两个Krylov子空间：
\[\mathcal{K}_m(A, v_1) = \text{span}\{v_1, Av_1, A^2v_1, \dots, A^{m-1}v_1\} \]
\[ \mathcal{K}_m(A^H, w_1) = \text{span}\{w_1, A^H w_1, (A^H)^2 w_1, \dots, (A^H)^{m-1} w_1\} \]
其中 \(v_1\) 和 \(w_1\) 是两个给定的、满足 \(w_1^H v_1 = 1\) 的初始向量（通常随机选取并标准化）。目标是生成两组向量 \(\{v_1, v_2, \dots, v_m\}\) 和 \(\{w_1, w_2, \dots, w_m\}\)，使得：
- 它们分别构成两个Krylov子空间的基。
- 它们满足双正交性：\(w_i^H v_j = 0\) 当 \(i \neq j\)，且 \(w_i^H v_i = 1\)。
- 矩阵A在这两组基下的投影是一个三对角矩阵 \(T_m\)。
这样，大型稀疏矩阵A的特征值问题就转化为小型稠密三对角矩阵 \(T_m\) 的特征值问题。
步骤2：三对角矩阵的形式与投影关系
我们希望得到如下形式的分解：
\[A V_m = V_m T_m + \beta_{m+1} v_{m+1} e_m^T \]
\[ A^H W_m = W_m T_m^H + \bar{\gamma}_{m+1} w_{m+1} e_m^T \]
并且满足：
\[W_m^H V_m = I_m, \quad w_{m+1}^H v_{m+1} = 1, \quad W_m^H v_{m+1} = 0, \quad V_m^H w_{m+1} = 0 \]
其中：
- \(V_m = [v_1, v_2, \dots, v_m]\)， \(W_m = [w_1, w_2, \dots, w_m]\)
- \(e_m\) 是第m个标准单位向量
- \(T_m\) 是一个非对称的三对角矩阵：
\[T_m = \begin{pmatrix} \alpha_1 & \gamma_2 & & \\ \beta_2 & \alpha_2 & \gamma_3 & \\ & \beta_3 & \alpha_3 & \ddots \\ & & \ddots & \ddots & \gamma_m \\ & & & \beta_m & \alpha_m \end{pmatrix} \]
- \(\beta_{m+1}\) 和 \(\gamma_{m+1}\) 是标量，用于控制向后迭代。
关键：如果算法精确执行，则 \(T_m\) 的特征值（称为Ritz值）是A的特征值的近似，\(V_m y\)（其中y是 \(T_m\) 的特征向量）是A对应特征向量的近似。
步骤3：非对称Lanczos双正交化过程（核心递推公式）
算法从随机初始向量 \(v_1, w_1\) 开始，满足 \(\|v_1\|_2 = \|w_1\|_2 = 1\) 且 \(w_1^H v_1 = 1\)。然后进行迭代，对于 \(j = 1, 2, \dots, m\)：
- 计算矩阵乘积：
\[ \hat{v} = A v_j \]
\[ \hat{w} = A^H w_j \]
这是算法中计算量最大的步骤，但由于A是稀疏的，可以高效计算。
- 计算对角元 \(\alpha_j\) ：
\[ \alpha_j = w_j^H \hat{v} \quad (\text{或等价地 } w_j^H A v_j) \]
注意，由于双正交性，\(\alpha_j\) 也可以从 \(\hat{v}\) 和 \(w_j\) 的内积得到。
- 正交化（双正交化）：
- 对于向量 \(v\)：
\[ v_{temp} = \hat{v} - \alpha_j v_j - \beta_j v_{j-1} \]
这里我们约定 $ \beta_1 v_0 = 0 $。
- 对于向量 \(w\)：
\[ w_{temp} = \hat{w} - \bar{\alpha}_j w_j - \bar{\gamma}_j w_{j-1} \]
同样约定 $ \gamma_1 w_0 = 0 $。
- 注意，这里的 \(\beta_j\) 和 \(\gamma_j\) 是上一步（j-1步）已经计算出的次对角元。
- 确保双正交性并计算新系数：
理论上，\(v_{temp}\) 应该与 \(w_1, \dots, w_j\) 正交，\(w_{temp}\) 应该与 \(v_1, \dots, v_j\) 正交。但实际上，由于数值误差，可能需要进行一次“完全正交化”或选择性重正交化来保持稳定性（这是一个实际实现中的重要细节，但基础算法描述中常假设精确算术）。
然后，计算：
\[ \delta = w_{temp}^H v_{temp} \]
注意，在精确算术和理想双正交条件下，\(\delta\) 应等于 \(\beta_{j+1} \gamma_{j+1}\)。我们通常先计算一个标量积：
\[ \beta_{j+1} = \|v_{temp}\|_2, \quad \gamma_{j+1} = \frac{\delta}{\beta_{j+1}} \]
但更常见的标准做法是：
\[ \beta_{j+1} = \|v_{temp}\|_2, \quad \gamma_{j+1} = \frac{w_{temp}^H v_{temp}}{\beta_{j+1}} \]
并检查 \(\beta_{j+1} \gamma_{j+1} \approx \delta\)。为了保证 \( w_{j+1}^H v_{j+1} = 1 \，我们设置：
\[ v_{j+1} = v_{temp} / \beta_{j+1} \]
\[ w_{j+1} = w_{temp} / \bar{\gamma}_{j+1} \]
这里除法的目的是标准化，使得 \(w_{j+1}^H v_{j+1} = 1\)。
递推公式总结（标准的三项递归形式）：
\[\beta_{j+1} v_{j+1} = A v_j - \alpha_j v_j - \gamma_j v_{j-1} \]
\[ \bar{\gamma}_{j+1} w_{j+1} = A^H w_j - \bar{\alpha}_j w_j - \bar{\beta}_j w_{j-1} \]
其中 \(\gamma_1 v_0 = 0\), \(\beta_1 w_0 = 0\)。注意，在实际算法中，我们通常让 \(\beta_{j+1}\) 和 \(\gamma_{j+1}\) 为实数（取模），从而 \(v_{j+1}\) 和 \(w_{j+1}\) 的范数为1。
步骤4：算法步骤伪代码
- 输入：矩阵A，初始向量 \(v_1, w_1\)（满足 \(\|v_1\|=\|w_1\|=1, w_1^H v_1=1\)），迭代步数m。
- 初始化：\(\beta_1 = 0, \gamma_1 = 0, v_0 = 0, w_0 = 0\)。
- For \(j = 1\) to \(m\) do：
a. 计算 \(\alpha_j = w_j^H (A v_j)\)。
b. 计算 \(\hat{v} = A v_j - \alpha_j v_j - \beta_j v_{j-1}\)。
c. 计算 \(\hat{w} = A^H w_j - \bar{\alpha}_j w_j - \bar{\gamma}_j w_{j-1}\)。
d. 计算 \(\beta_{j+1} = \|\hat{v}\|_2\)。如果 \(\beta_{j+1} = 0\)，则算法提前终止（子空间不变）。
e. 计算 \(\gamma_{j+1} = \hat{w}^H \hat{v} / \beta_{j+1}\)。
f. 计算 \(v_{j+1} = \hat{v} / \beta_{j+1}\)。
g. 计算 \(w_{j+1} = \hat{w} / \bar{\gamma}_{j+1}\)。 - 输出：三对角矩阵 \(T_m\)（由 \(\alpha_1, \dots, \alpha_m\), \(\beta_2, \dots, \beta_m\), \(\gamma_2, \dots, \gamma_m\) 构成），以及矩阵 \(V_m, W_m\)。
步骤5：计算特征值近似（Ritz值）
得到 \(T_m\) 后，我们通过标准的稠密特征值算法（如QR算法）计算其全部特征值 \(\theta_i\) 和右特征向量 \(s_i\)（满足 \(T_m s_i = \theta_i s_i\)）。
- 近似特征值（Ritz值）：\(\theta_i\)。
- 近似右特征向量（Ritz向量）：\(u_i = V_m s_i\)。
- 近似左特征向量：可以通过 \(W_m t_i\) 得到，其中 \(t_i\) 是 \(T_m^H\) 的特征向量。
残差估计：可以证明，近似特征对 \((\theta_i, u_i)\) 的残差范数满足：
\[\|A u_i - \theta_i u_i\|_2 = |\beta_{m+1}| \cdot |e_m^T s_i| \]
由于 \(\beta_{m+1}\) 已知，\(|e_m^T s_i|\) 可以从 \(s_i\) 的最后一个分量得到。这个量可以用来评估近似的精度。当残差足够小时，\((\theta_i, u_i)\) 就是原矩阵一个很好的近似特征对。
步骤6：算法的性质、优势与挑战
-
优势：
- 低存储：与Arnoldi算法需要存储所有基向量不同，Lanczos双正交化理论上只需要存储最近的几个向量（三项递归），但为了计算Ritz向量，通常仍需存储 \(V_m\) 和 \(W_m\)。
- 高效：每次迭代只需两次矩阵-向量乘（A和A^H）。
- 降维：将大规模问题转化为小规模三对角特征值问题。
-
挑战与注意事项：
- 数值不稳定性：在有限精度计算中，双正交性会逐渐丧失，导致“幽灵”特征值（重复出现的近似特征值）等问题。通常需要采用重正交化策略，即定期用全部或部分已生成的基向量对当前向量进行重新正交化，但这会增加计算和存储成本。
- 提前终止：如果 \(\beta_{j+1}\) 或 \(\gamma_{j+1}\) 为零，则Krylov子空间是不变的，算法可以提前终止，并找到精确的不变子空间。
- A^H 的需求：需要计算 \(A^H w_j\)，对于某些应用，A^H的快速计算可能不易实现。
步骤7：总结
非对称Lanczos双正交化算法通过构建双正交基，为非对称矩阵特征值问题提供了一个强大的Krylov子空间框架。它将原问题投影到一个小型三对角矩阵上，使得计算部分特征值成为可能。尽管存在数值稳定性的挑战，但通过重正交化等技术，它在实际的大规模科学计算问题中（如流体力学、量子化学等领域的线性稳定性分析）有着广泛的应用。算法的核心在于巧妙的三项递归关系，它同时推进两个Krylov子空间，并保持它们之间的双正交性。