计算内部特征值

如果只寻找最小或最大的特征值对应的特征对，那么显然的重启方法效果相当好，但通常如果感兴趣的是内部特征值，这种方法表现不佳。问题在于，里兹值会单调地向外部特征值收敛，而一个接近目标值的里兹值可能正朝着其他外部特征值发展。甚至可能出现对应的里兹向量在所需特征向量方向上的分量很小的情况。显然，这样的里兹向量作为重启的候选者质量很差，因此问题来了，什么样的向量更适合重启？一个答案是由所谓的调和里兹向量提供的，这在第§3.2节中讨论过；参见[331,349,411]。

如我们所见，雅可比-戴维森方法生成子空间{\mathcal V}_{m}的基向量{v}_1, \ldots, {v}_{m}。对于A在这个子空间上的投影，我们计算向量{w}_j\equiv A{v}_j。调和里兹值是A^{-1}的里兹值的倒数，相对于由{w}_j张成的子空间。它们可以在不实际求逆A的情况下计算，因为一个调和里兹对(\widetilde{\theta}_j^{(m)},\widetilde{u}_j^{(m)})满足

A\widetilde{u}_j^{(m)} -\widetilde{\theta}_j^{(m)} \widetilde{u}_j^{(m)} \perp{\mathcal W}_{m} \equiv \mathrm{span}({w}_1,\ldots , {w}_{m}) \tag{4.51}

对于\widetilde{u}_j^{(m)} \in {V}_{m}且\widetilde{u}_j^{(m)}\neq0。这意味着调和里兹值是矩阵束({W}_{m}^\ast A{V}_{m}, {W}_{m}^\ast {V}_{m})的特征值，或者由于{W}_{m}=A{V}_{m}：

{W}_{m}^\ast {W}_{m} {s}_j^{(m)} - \widetilde{\theta}_j^{(m)}{W}_{m}^\ast {V}_{m} {s}_j^{(m)} =0 \, .\,

出于稳定性考虑，我们对{W}_{m}的列进行正交化处理，并相应地变换{V}_{m}的列。这进一步简化了方程：我们看到调和里兹值是对称矩阵{W}_{m}^\ast {V}_{m}的特征值的倒数。

在[349]中表明，对于厄米矩阵A，调和里兹值单调地向绝对值最小的非零特征值收敛。注意，调和里兹值无法识别A的零特征值，因为那对应于A^{-1}的无穷大特征值。同样，移位矩阵A-\tau I的调和里兹值单调地向最接近目标值\tau的特征值\lambda\neq \tau收敛。幸运的是，移位矩阵和未移位矩阵的搜索子空间{\mathcal V}_m相同，这便于计算任意移位的调和里兹对。移位矩阵对应的调和里兹向量，其调和里兹值最接近\tau，可以解释为最大化(A-\tau I)^{-1}的瑞利商。它代表了渐近意义上可获得的关于所需特征值的最佳信息，因此渐近意义上是重启时作为起始向量的最佳候选，前提是\tau\neq \lambda。

对于调和里兹值，校正方程需要考虑与A\tilde u^{(m)}_j的正交性，这导致斜投影。我们可以通过以下方式使用正交投影。如果\tilde u=\tilde u^{(m)}_j是被选中的特征向量近似，瑞利商\theta = \tilde u^\ast A\tilde u/\tilde u^\ast \tilde u导致具有最小范数的残差；即，对于r\equiv A\tilde u-\theta \tilde u，我们有\Vert r\Vert _2\leq \Vert A\tilde u-\mu\tilde u\Vert _2，对于任何标量\mu，包括调和里兹值\mu=\tilde\theta^{(m)}_j。此外，瑞利商的残差r与\tilde u正交。这使得r与校正方程中的算子(I- u u^\ast)(A-\theta I)(I-u u^\ast)“兼容”。这里u\equiv \tilde u/\Vert\tilde u\Vert _2。

基于调和里兹值和向量的雅可比-戴维森方法，结合重启和消减，在算法4.19中给出。该算法可用于计算目标值\tau右侧的一系列连续特征值。

要应用此算法，我们需要指定一个起始向量{v}_0，一个容差\epsilon，一个目标值\tau，以及一个指定应计算多少个接近\tau的特征对的数值{k}_{\max}。{m}_{\max}表示搜索子空间的最大维度。如果超过此值，将进行重启，使用指定维度{m}_{\min}的子空间。

完成后，将提供最接近\tau右侧的{k}_{\max}个特征值。计算的特征对(\widetilde\lambda_j,\widetilde{x}_{j})，\Vert\widetilde{x}_{j}\Vert _2=1，满足\Vert A \widetilde{x}_{j}-\widetilde\lambda_j\widetilde{x}_{j}\Vert _2\leq j\epsilon，其中\widetilde{x}_j表示\widetilde{X}的第j列。

对于外部特征值，已在第§4.7.3节中描述了一个更简单的算法。现在我们将根据前几小节的讨论，对算法的一些部分进行评论。

(1)

初始化阶段。

(3)-(7)

向量w = (A-\tau I)t通过修正的Gram-Schmidt方法与当前测试子空间{\mathcal W}_{m}正交。为了提高数值稳定性，可以用算法4.14中给出的模板（对于向量t）进行替代。

(8)-(10)

数值{M}_{i,j}表示方阵{M}\equiv {W}^\ast {V}的元素，其中{V}表示列向量为{v}_j的n乘{m}矩阵，{W}同理。因为{M}是厄米的，所以只计算该矩阵的上三角部分。

(11)-(13)

此时应计算问题{M}{s}=\widetilde\theta {s}的特征对。这可以用LAPACK中适用于厄米稠密矩阵的合适例程来完成。注意，调和里兹值是{M}的特征值的倒数。我们需要计算{V}{s}_1的瑞利商\theta，然后归一化{V}{s}_1，以便计算一个适当的残差r \perp {V}{s}_1。我们使用了{s}_1^\ast {V}^\ast (A-\tau I) {V}{s}_1 = {s}_1^\ast{M}^\ast {s}_1=\widetilde\theta_1。

向量s_j是m乘m矩阵S的列，\tilde\Theta =\mathrm{diag}(\widetilde\theta_1,\ldots,\widetilde\theta_m)。

(14)

停止准则是在归一化特征向量近似的残差范数低于\epsilon时接受一个特征向量近似。这意味着我们接受计算特征值的误差在\epsilon^2数量级，以及特征向量的角度误差在{\mathcal O}({\epsilon})数量级，前提是相关特征值是简单的且与其他特征值充分分离；参见(4.4)。

无法保证检测到所有想要的特征值；参见算法4.13的注释(14)和算法4.17的注释(13)。

(17)

这是在接受近似特征对后的重启。重启稍微复杂一些，因为涉及两个子空间。从这些子空间重新计算生成向量在(18)-(21)中完成。

(24)

此时进行重启，当子空间维度超过{m}_{\max}时。重启后，雅可比-戴维森迭代以维度{m}_{\min}的子空间继续进行。该子空间的选择基于最接近目标\tau的调和里兹值。

(31)-(32)

通过\widetilde{X}的因子表示计算特征向量的消减。矩阵\widetilde{X}的列是计算的特征向量。如果存在左预条件子{K}用于算子A-\theta I$，则可以像处理外部特征值的情况一样，通过类似Krylov求解器实现类似的简化。处理左预条件算子的有效模板在算法4.18中给出。