算法

为了利用Lanczos过程解决广义特征值问题 (8.21)，我们首先将其转换为等价的标准特征值问题。有几种可能的变换方法；以下列出了最常见的几种。

(a)

如果需要求解矩阵对\{A,B\}中模最大的几个特征值，且矩阵B非奇异，我们可以通过乘以B^{-1}来得到

B^{-1} A x = \lambda x.

矩阵B^{-1}A关于A或B对称。

(b)

如果需要求解接近给定值\sigma的几个特征值，我们应用SI方法

\left( (A - \sigma B)^{-1} B \right) x = (\lambda - \sigma)^{-1} x.

特别是，如果需要求解模最小的特征值，可以选择移位值\sigma=0。注意，矩阵(A - \sigma B)^{-1} B关于B或(A-\sigma B)对称。

本节余下部分将重点讨论情况(b)，并将变换后的问题写为

H^{-1} B x = \mu x,

其中

H = A - \sigma B\quad \mathrm{且} \quad\mu = \frac{1}{\lambda-\sigma}.

注意，即使A和B是实数矩阵，\{A,B\}的特征值也可能是复数，因此可能需要使用复数\sigma。即使使用了复数移位，Lanczos递归仍可以通过实数运算实现；参见 [362]。

对称不定Lanczos算法的模板在算法 8.4中展示。生成的Lanczos向量\{q_j\}和标量\{ \alpha_j, \beta_j, \omega_j\}满足以下控制方程：

(H^{-1} B ) Q_j = Q_j (\Omega^{-1}_j T_j ) +\frac{\beta_{j+1}}{\omega_{j+1}} q_{j+1} e^T_j, \tag{8.22}

其中

T_j = \begin{bmatrix}\alpha_1 & \beta_2 & & \\ \beta_2 & \alpha_2 & \ddots & \\ & \ddots & \ddots & \\ & & \beta_j & \alpha_j \end{bmatrix} \quad\text{和}\quad \Omega_j = \mathrm{diag}(\omega_1,\omega_2, \ldots, \omega_j).

Lanczos向量 Q_j = \left[ \begin{array}{cccc}q_1 & q_2 & \cdots & q_j \\\end{array} \right]是B正交的，即，

Q^T_j B Q_j = \Omega_j, \quad\quad Q^T_j B q_{j+1} = 0,

并且被归一化为单位欧几里得范数，即\Vert q_j\Vert _2 = 1。

简化的特征值问题是

T_j u^{(j)}_i = \theta^{(j)}_i \Omega_j u^{(j)}_i.

它继承了与原问题相同的数值困难。例如，Ritz值\theta^{(j)}_i可能是复数，甚至是缺陷的。换句话说，它可能属于\Omega^{-1}_j T_j的非对角Jordan块。

一旦计算出\{\theta^{(j)}_i,u^{(j)}_i\}，右和左Ritz向量分别定义为

{x}^{(j)}_i := Q_j u^{(j)}_i \quad \mathrm{且} \quad{y}^{(j)}_i := B \bar{Q}_j \bar{u}^{(j)}_i = B \bar{x}^{(j)}_i,

相应的。量\{\theta^{(j)}_i, x^{(j)}_i, y^{(j)}_i\}被称为矩阵H^{-1} B的Ritz三元组。

与Ritz对\{\theta^{(j)}_i, x^{(j)}_i\}对应的右残差向量r^{(j)}_i是

r^{(j)}_i = (H^{-1}B)x^{(j)}_i - \theta^{(j)}_ix^{(j)}_i= \gamma^{(j)}_i q_{j+1}, \tag{8.23}

其中

\gamma^{(j)}_i = \frac{\beta_{j+1}}{\omega_{j+1}} (e_j^Tu^{(j)}_i).

注意，在 (8.23)中的最后一个等式是通过将方程 (8.22)右乘u^{(j)}_i并利用x^{(j)}_i的定义得到的。左残差向量(s^{(j)}_i)^{\ast}与右残差向量r^{(j)}_i的关系为(s^{(j)}_i)^{\ast} = (r^{(j)}_i)^T B。由此得出

\Vert r^{(j)}_i\Vert _2 = \vert\gamma^{(j)}_i\vert \tag{8.24}

和

\Vert s^{(j)}_i\Vert _2 = \vert\gamma^{(j)}_i\vert\, \Vert B q_{j+1}\Vert _2. \tag{8.25}

注意，向量Bq_{j+1}在对称不定Lanczos算法的每次迭代中都是直接可用的。不需要额外的B调用。因此，Lanczos算法的一个吸引人的特点是可以在不显式计算Ritz向量{x}^{(j)}_i和{y}^{(j)}_i的情况下计算残差向量。

算法 8.4：对称不定Lanczos方法
(1) 选择一个起始向量 q
(2) 计算 \beta_1 = \|q\|^2
(3) 缩放 g_1 = \frac{q}{\beta_1}
(4) w = Bq_1
(5) w_i = w^T g_1
(6) 对于 j = 1,2,... 直到收敛
(7) 解 Hr = w 得到 r
(8) 对于 j > 1，r := r - q_{j-1}(B_{j-1}^T w / w_{j-1})
(9) a_j = w^T r
(10) r := r - \frac{\alpha_j}{w_j} g_j
(11) T = \|r\|^2
(12) 如果 T = 0，则停止
(13) 如有需要则重新正交化
(14) \beta_{j+1} = \|r\|
(15)  w = Bq_{j+1} 
(16) w_{j+1} = w^T g_{j+1}
(17) 如果 w_{j+1} = 0，则停止
(18) B_{j+1} = \tau w_{j+1}
(19) 解简化问题并检查收敛性。
(20) 结束循环
(21) 计算近似特征向量 x_i^{(j)}

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

(1)

理想的初始向量是与所需特征值对应的特征向量的线性组合。如果无法获得此类向量的近似值，则可以使用随机向量。初始向量r通常通过乘以H^{-1} B进行预处理。这一预处理步骤在[161]中提出，并且在B奇异时至关重要；详见§8.6.4。

(4)和(15)

矩阵向量乘法Bq_j应通过用户提供的子程序执行，该子程序可以利用矩阵B的数据结构。

(7)

必须求解线性方程组H r = w以得到r。设

H = {\mathcal L}{\mathcal D}{\mathcal L}^T

为稀疏LDL^{\rm T}分解（参见§10.3），在for循环之前计算。然后向量r可以通过三个步骤高效计算： (7)^{\prime} 解{\mathcal L} u = w 得到u
解{\mathcal D} v = u得到v
解{\mathcal L}^T r = v 得到r

(12)

当\Vert r\Vert _2 = 0时，表明j\times j矩阵对T_j - \lambda \Omega_j的所有特征值也是H^{-1} B的特征值，且Q_j的列张成一个不变子空间。可以选择退出算法或继续算法，将\beta_j的值设为零并选择r为与Q_j的列B正交的随机向量。

(13)

在有限精度算术中，向量\{q_j\}不再是B正交的。标准对称Lanczos过程的再正交化方案可以扩展到对称不定算法；详见§4.4。

(14)

在对称Lanczos算法（参见§4.4）中，Lanczos向量被缩放为B正交，即q_j^TBq_j = 1。然而，在对称不定算法中，向量被缩放使得

\Vert q_j\Vert _2 = 1.

这种缩放变化有两个原因，都与B的不定性有关。第一个原因是，即使A和B是实数，将Lanczos向量缩放为q_j^TBq_j = 1可能涉及复数运算。（通过跟踪q_j^TBq_j的符号，可以避免复数运算，但仍使用\vert q_j^TBq_j \vert^{1/2}作为范数。然而，这种算法稳定性较差。）第二个原因是，当B不定性时，并非每个子空间都有B正交基，即使有，也可能高度病态[417,20]。这种缩放选择并不能解决Lanczos向量基可能病态的问题。然而，它提供了一种方便的方法来检测基何时变得病态，具体来说，当\vert\omega_{j+1}\vert值极小时[357]。

(17)

如果\omega_{j+1}=0，Lanczos过程会遭遇类似于非厄米Lanczos过程中发生的崩溃；参见§7.8。可以采用针对非厄米Lanczos过程开发的相同补救措施，如前瞻[364,178]和自适应阻塞[29,298]。

(19)

矩阵束A - \lambda B的近似特征值通过求解简化的j\times j广义特征值问题来确定

T_j u = \theta \Omega_j u.

有几种算法可用于求解此问题。如果j较小，可以应用标准QZ算法于\{T_j,\Omega_j\}或QR算法于\Omega_j^{-1}T_j。然而，这些算法都没有利用\{T_j,\Omega_j\}的结构。在[357]中包含了对利用对称性的算法的良好综述，一些更新的工作可以在[406,407]中找到。

对于中等大小的j（例如，j\approx 100），求解简化问题在计算上变得显著。对于较长的Lanczos运行，建议定期求解简化问题，可能每五到十次Lanczos迭代一次。

对称不定Lanczos方法内循环的适当停止准则是

\max\left\{ \Vert r^{(j)}_i\Vert _2, \Vert s^{(j)}_i\Vert_2 \right\} = \vert\gamma_i^{(j)}\vert \max \left\{ 1, \Vert Bq_{k+1}\Vert _2 \right\} \leq \epsilon,

其中\epsilon是用户给定的容差值。注意，Bq_{k+1}在Lanczos运行结束时可用，其欧几里得范数可以通过一次点积计算。无需额外调用B来计算该项。关于收敛性的详细测试程序在以下§8.6.3中讨论。

(21)

仅当对应的Ritz值根据以下描述的测试收敛时，才计算近似特征向量。