算法

下面介绍的非厄米Lanczos方法（参见算法7.13）是一种双侧迭代算法，起始向量为p_1和q_1。它可以看作是通过双侧Gram-Schmidt过程对两个Krylov序列进行双正交化：

\{ q_1, Aq_1, A^2q_1, \ldots \} \quad \text{和} \quad \{ p_1, A^{\ast}p_1, (A^{\ast})^2 p_1, \ldots \}.

生成的两个向量序列\{q_i\}和\{p_i\}使用三项递归关系：

\beta_{j+1}q_{j+1} = A q_j - \alpha_j q_j - \gamma_{j} q_{j-1}, \tag{7.33}

\bar{\gamma}_{j+1} p_{j+1} = A^{\ast} p_j - \bar{\alpha}_j p_j - \bar{\beta}_{j} p_{j-1}. \tag{7.34}

向量\{q_i\}和\{p_i\}被称为Lanczos向量，分别张成{\mathcal K}^j(A, q_1)和{\mathcal K}^j(A^{\ast}, p_1)，并且是双正交的，即p^{\ast}_{k} q_{\ell} = 0如果k \neq \ell，以及w^{\ast}_{k} v_{k} = 1。用矩阵表示，在第j步，Lanczos方法生成两个n \times j矩阵Q_j和P_j：

Q_j = [ q_1, q_2, \ldots, q_j ], \quad P_j = [ p_1, p_2, \ldots, p_j ]

以及

T_j = \begin{bmatrix} \alpha_1 & \gamma_2 & & & \\ \beta_2 & \alpha_2 & \gamma_3 & & \\ & \beta_3 & \alpha_3 & \ddots & \\ & & \ddots & \ddots & \gamma_j & \\ & & & \beta_j & \alpha_j \end{bmatrix}.

计算得到的量满足称为Lanczos分解的控制关系：

AQ_j = Q_j T_j + \beta_{j+1} q_{j+1} e_{j}^{\ast}, \tag{7.35}

A^{\ast} P_j = P_j T_j^{\ast} + \bar{\gamma}_{j+1} p_{j+1} e_j^{\ast}, \tag{7.36}

P_j^{\ast} Q_j = I_j. \tag{7.37}

此外，P_j^{\ast} q_{j+1} = 0和p_{j+1}^{\ast} Q_j = 0。关系式(7.37)表明Lanczos向量（基）是双正交的。但注意Q_j和P_j都不是酉矩阵。在Lanczos基下，矩阵A由一个非厄米三对角矩阵表示，

P_j^{\ast} AQ_j = T_j. \tag{7.38}

在任何步骤j，我们可以计算T_j的特征解，

T_j z_i^{(j)} = \theta_i^{(j)} z_i^{(j)} \quad \text{和} \quad (w_i^{(j)})^{\ast} T_j = \theta_i^{(j)} (w_i^{(j)})^{\ast}.

矩阵A的特征值由T_j的特征值\theta_i^{(j)}近似，这些特征值称为Ritz值。

对于每个Ritz值\theta_i^{(j)}，存在相应的右和左Ritz向量，

x_i^{(j)} = Q_j z_i^{(j)} \quad \text{和} \quad y_i^{(j)} = P_j w_i^{(j)}. \tag{7.39}

Ritz值和向量收敛到原始矩阵A的特征值和特征向量的程度可以通过比较残差的范数与x_i^{(j)}和y_i^{(j)}的范数来评估，

r_i^{(j)} = A x_i^{(j)} - \theta_i^{(j)} x_i^{(j)}, \tag{7.40}

(s_i^{(j)})^{\ast} = (y_i^{(j)})^{\ast} A - \theta_i^{(j)} (y_i^{(j)})^{\ast}. \tag{7.41}

此外，根据方程(7.35)，右残差向量变为

r_i^{(j)} = \beta_{j+1} q_{j+1} (e_j^{\ast} z_i^{(j)}), \tag{148}

根据方程(7.36)，左残差向量变为

(s_i^{(j)})^{\ast} = \gamma_{j+1} p_{j+1}^{\ast} ((w_i^{(j)})^{\ast} e_j). \tag{149}

因此，与厄米情况（参见§4.4）一样，残差范数可以在不显式计算Ritz向量x_i^{(j)}和y_i^{(j)}的情况下获得，尽管\|x_i^{(j)}\|_2和\|y_i^{(j)}\|_2是不可用的。有关此主题的更详细讨论，请参见§7.8.2。

算法 7.13：NHEP的Lanczos方法

(1)  选择起始向量 q_1 和 p_1，使得 p_1^* q_1 = 1。
(2)  r = A q_1
(3)  s = A^* p_1
(4)  对于 j = 1, 2, ..., 直到收敛，
(5)      \alpha_j = p_j^* r
(6)      r := r - q_j \alpha_j
(7)      s := s - p_j \alpha_j
(8)      如果 (\|r\| = 0 和/或 \|s\| = 0)，则停止
(9)      w_j = r^* s
(10)     如果 (w_j = 0)，则停止
(11)     B_{j+1} = |w_j|^{1/2}
(12)     Y_{j+1} = w_j / B_{j+1}
(13)     q_{j+1} = r / \beta_{j+1}
(14)     P_{j+1} = s / Y_{j+1}
(15)     计算 T_j 的特征三元组 \left(\theta_i^{(j)}, z_i^{(j)}, w_i^{(j)}\right)。
(16)     测试收敛性。
(17)     如有必要，重新正交化。
(18)     r = A q_{j+1}
(19)     s = A^* p_{j+1}
(20)     r := r - q_j Y_{j+1}
(21)     s := s - p_j \beta_{j+1}
(22) 结束循环
(23) 计算近似特征向量 x_i^{(j)} 和 y_i^{(j)}。

我们现在对算法7.13的某些步骤进行评论：

(1)

初始起始向量p_1和q_1最好由用户根据关于A的所需特征向量的任何可用知识来选择。在没有此类知识的情况下，可以选择具有随机分布条目的q_1，并令p_1 = q_1。

(2), (3), (18), (19)

在这些步骤中需要矩阵-向量乘法例程来乘以A和A^{\ast}与任意向量。这通常是计算瓶颈。有关移位-逆情况下的实现注释，请参见下面的收敛性讨论。

(8)

这是方法可能崩溃的两种情况之一。实际上，这是一种理想的崩溃。如果r为零，则Lanczos向量\{q_1, q_2, \ldots, q_j\}张成A的一个（右）不变子空间；每个Ritz值和相应的右Ritz向量是A的一个精确特征值和特征向量。如果需要，可以通过选择任何满足P_j^{\ast} q_{j+1} = 0的q_{j+1}并设置\beta_{j+1} = 0来继续Lanczos算法。当s或r和s都消失时，可以采取类似的行动。

在实践中，精确的零向量很少见。r和/或s的范数可能会变得很小。应该给出一个容差值来检测相对于\|Q\|_2 \|T\|_2或\|P\|_2 \|T\|_2的微小\|r\|_2，以及相对于\|P\|_2 \|T\|_2的微小\|s\|_2。默认容差值是机器精度\epsilon的小倍数。

(10)

如果\omega_j = r^{\ast} s = 0在r或s消失之前，方法会出现本质崩溃。在大多数情况下，我们可以继续在Krylov子空间{\mathcal K}^{j+k}(A, r)和{\mathcal K}^{j+k}(A^{\ast}, s)中寻找新向量，对于某个整数k > 0，并在T_j的三对角线之外添加一个块；这种所谓的向前看过程在[178]中有描述，并在QMRPACK中实现（参见§7.8.3）。然而，如果我们的起始向量q_1和p_1具有不同的最小多项式（或者，例如，p_1和q_1是不同特征向量集的复合），即使这也不起作用，我们会有一个不匹配，也称为不可治愈的崩溃。有关进一步讨论，请参见[354]。在§7.9中讨论了崩溃的不同处理方法。

在实践中，精确的崩溃很少见。接近崩溃的情况更常见；即，\omega_j是非零的，但绝对值非常小。接近崩溃会导致停滞和不稳定。任何检测接近崩溃的准则要么在某些情况下过早停止，要么在其他情况下过晚停止。检测特征值问题中接近崩溃的合理折中准则是，如果\vert\omega_j\vert \leq \sqrt{\epsilon} \|r\|_2 \|s\|_2，则停止。

(15)

通过QR算法（参见§7.3）计算三对角矩阵T_j的特征分解的成本大约是每次迭代30j^3次浮点运算，累积成本从步骤1到j是15j^4次浮点运算。定期解决特征问题，例如每10步一次，可以降低这一成本。

尚未找到在j^2次浮点运算中近似一般三对角矩阵特征值的稳定算法，尽管最近提出了条件稳定算法，如[94, 189]。据作者所知，没有Lanczos算法的软件实现使用快速条件稳定特征求解器。

(16)

一旦确定了基Q_j和P_j，使得T_j的特征值\theta_i^{(j)}（参见(7.38)）近似于A的所有所需特征值，且残差足够小，计算就会停止，这些残差根据方程(7.42)和(7.43)计算。有关收敛性的更详细讨论，请参见§7.8.2。

如果没有重新双正交化（参见[105]），那么在有限精度算术中，一旦Ritz值收敛到A的特征值，这个Ritz值将在后续的Lanczos步骤中出现副本。例如，缩减三对角矩阵T_j的Ritz值集群可能近似于原始矩阵A的单个特征值。伪值[93]是一个简单的Ritz值，也是从T_j中删除第一行和第一列得到的j-1阶矩阵的特征值。这样的伪值应从考虑中剔除。T_j的非伪特征值被识别为原始矩阵A的特征值近似，并进行收敛性测试。

(17)

与厄米情况一样，在有限精度算术中，计算的Lanczos向量\{q_i\}和\{p_i\}的双正交性会恶化。可以使用双侧修正Gram-Schmidt（TSMGS）过程[354]重新双正交化第(j+1)个Lanczos向量；

for i=1,2,\ldots,j
q_{j+1} = q_{j+1} - q_i (p_i^{\ast} q_{j+1})
p_{j+1} = p_{j+1} - p_i (q_i^{\ast} p_{j+1})
end for

在每一步应用TSMGS过程非常耗时，并成为计算瓶颈。在[105]中，提出了一种有效的替代方案。这一主题在§7.9中重新讨论。