计算内部特征值

为了重启目的，特别是当目标是内部特征值时，调和Ritz向量已被推荐用于对称矩阵[331]；参见第§4.7.4节。

调和Ritz值的概念[349]很容易推广到非对称矩阵[411]。如我们所见，Jacobi-Davidson方法生成子空间{\mathcal V}_{m}的基向量{v}_1, \ldots, {v}_{m}。为了将A投影到这个子空间上，我们需要计算向量{w}_j\equiv A{v}_j。相对于由{w}_j张成的子空间的A^{-1}的Ritz值的倒数被称为调和Ritz值。调和Ritz值可以在不反转A的情况下计算，因为调和Ritz对(\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{7.103}

对于\widetilde{u}_j^{({m})} \in {V}_{m}且\widetilde{u}_j^{({m})}\neq 0。这意味着调和Ritz值是矩阵束(W_{m}^\ast AV_{m}, W_{m}^\ast V_{m})的特征值：

W_{m}^\ast AV_{m} s_j^{({m})} -\widetilde\theta_j^{({m})} W_{m}^\ast V_{m} s_j^{({m})} =0 .

外部标准Ritz值通常收敛到A的外部特征值。同样，移位矩阵A-\tau I的内部调和Ritz值通常收敛到最接近移位\tau的特征值\lambda\neq \tau。幸运的是，移位矩阵和未移位矩阵的搜索子空间{\mathcal V}_{m}一致，这有助于计算调和Ritz对。出于稳定性考虑，我们构造了正交的{V}_{m}和{W}_{m}：{W}_{m}使得(A-\tau I){V}_{m}={W}_{m} M^A_{m}，其中M^A_{m}是{m}阶上三角矩阵。

在生成的方案中，我们计算的是（部分）Schur分解而不是部分特征分解。也就是说，我们希望计算向量{q}_1,\ldots, {q}_{k}，使得A{Q}_k={Q}_k{R}_k，其中{Q}_{k}^\ast {Q}_{k}=I_{k}，且{R}_{k}是{k}阶上三角矩阵。{R}_{k}的对角元素表示A的特征值，对应的特征向量可以从{Q}_{k}和{R}_{k}计算得到。

基于调和Ritz值和向量的Jacobi-Davidson算法，包括部分Schur形式的约简，在算法7.19中给出。该算法包括重启和缩减技术。它可以用于计算接近\tau的多个特征值。

算法 7.19：Jacobi-Davidson 方法用于 NHEP 的 k_{\max} 内部特征值

(1) t = v_0, k = 0, m = 0; Q = [], R = []
(2) 当 k < k_{\max}
(3)     对于 i = 1, \ldots, m
(4)         t = t - (v_i^* t) v_i
(5)     m = m + 1, v_m = t / \|t\|_2, v_m^A = A v_m - \tau v_m, w = v_m^A
(6)     对于 i = 1, \ldots, k
(7)         w = w - (q_i^* w) q_i
(8)     对于 i = 1, \ldots, m-1
(9)         M_{i, m}^A = w_i^* w, w = w - M_{i, m}^A w_i
(10)    M_{m, m}^A = \|w\|_2, w_m = w / M_{m, m}^A
(11)    对于 i = 1, \ldots, m-1
(12)        M_{i, m} = w_i^* v_m, M_{m, i} = w_m^* v_i
(13)    M_{m, m} = w_m^* v_m
(14)    进行QZ分解 M^A S^R = S^L T^A, M S^R = S^L T, 使得：\left\|T_{i, i}^A / T_{i, i}\right\| \leq \left\|T_{i+1, i+1}^A / T_{i+1, i+1}\right\|
(15)    u = V s_1^R, u^A = V^A s_1^R, \vartheta = \overline{T_{1,1}} \cdot T_{1,1}^A
(16)    r = u^A - \vartheta u, \tilde{a} = Q^* r, \tilde{r} = r - Q \tilde{a}
(17)    当 \|\tilde{r}\|_2 \leq \epsilon
(18)        R = \left[\begin{array}{cc} R & \tilde{a} \\ 0 & \vartheta + \tau \end{array}\right], Q = [Q, u], k = k + 1
(19)        如果 k = k_{\max} 则停止
(20)        m = m - 1
(21)        对于 i = 1, \ldots, m
(22)            v_i = V s_{i+1}^R, v_i^A = V^A s_{i+1}^R, w_i = W s_{i+1}^L, s_i^R = s_i^L = e_i
(23)        M^A, M 是 T^A, T 的下 m 块
(24)        u = v_1, u^A = v_1^A, \vartheta = \overline{M_{1,1}} \cdot M_{1,1}^A
(25)        r = u^A - \vartheta u, \tilde{a} = Q^* r, \tilde{r} = r - Q \tilde{a}
(26)    结束当
(27)    如果 m \geq m_{\max} 则
(28)        对于 i = 2, \ldots, m_{\min}
(29)            v_i = V s_i^R, v_i^A = V^A s_i^R, w_i = W s_i^L
(30)        M^A, M 是 T^A, T 的前 m_{\min} 块
(31)        v_1 = u, v_1^A = u^A, w_1 = W s_1^L, m = m_{\min}
(32)    结束如果
(33)    \theta = \vartheta + \tau, \tilde{Q} = [Q, u]
(34)    解 t (\perp \tilde{Q}) （近似地）从：\left(I - \tilde{Q} \tilde{Q}^*\right)(A - \theta I)\left(I - \tilde{Q} \tilde{Q}^*\right) t = -\tilde{r}
(35) 结束当

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

完成后，将提供接近\tau的k_{\max}个特征值，以及相应的缩减Schur形式AQ=QR，其中Q是n阶正交矩阵，R是k_{\max}阶上三角矩阵。特征值位于R的对角线上。计算的Schur形式满足\Vert A q_{j}- {Q} R e_{j}\Vert _2\leq j \epsilon，其中q_j是Q的第j列。

我们将根据前几小节的讨论对算法的某些部分进行评论。

(1)

初始化阶段。

(3)-(4)

新计算的修正通过修正的Gram-Schmidt方法与当前搜索子空间正交化。我们选择存储矩阵V^A\equiv A{V}-\tau {V}；{v}_{m}^A是该矩阵的扩展向量。{W}的扩展向量通过使{v}_{m}^A与检测到的Schur向量空间和当前测试子空间通过修正的Gram-Schmidt方法正交化得到。为了提高数值稳定性，可以将Gram-Schmidt步骤替换为算法4.14中给出的模板。

{V}表示列向量为{v}_j的n阶矩阵；同样，W和V^A。

(8)-(9)

{M}_{i,j}表示方阵{M}\equiv {W}^\ast {V}的元素。{M^A}_{i,j}表示上三角矩阵{M^A}\equiv W^\ast V^A的元素。

（注意，如果{F}\equiv{Q}^\ast V^A，则V^A={W} {M^A}+{Q}{F}。因此，可以从{W}、{M^A}、{Q}和{F}重建V^A。特别是，可以从这些量计算r。与其存储{n}维矩阵V^A，不如存储{k}阶矩阵{F}（在(3)-(4)中计算的元素{q}_i^\ast {w}）。这种方法节省内存空间。然而，为了避免不稳定性，缩减过程需要特别注意。）

(14)

此时需要计算矩阵束(M^A,M)的QZ分解（广义Schur形式）：{M^A} {S}^R={S}^L{T^A}和{M}{S}^R={S}^L {T}，其中{S}^R和{S}^L是酉矩阵，{T^A}和{T}是上三角矩阵。这可以通过LAPACK中的适合密集矩阵束的例程来完成。

在每一步中，需要对Schur形式进行排序，使得\vert T^A_{1,1}/{T}_{1,1}\vert最小。只有当{m}\geq {m}_{\max}时，Schur形式的排序才需要使得{T^A}和{T}的前{m}_{\max}个对角元素表示最接近0的{m}_{\max}个调和Ritz值。为了便于说明，这里我们对所有对角元素进行了排序。

有关广义Schur形式的排序算法，请参见[33,448,449]和[171, Chap. 6B]。

{S^R}是m阶矩阵，列向量为{s}_j^R；同样，{S^L}。

(15)

\theta的值需要一些注意。我们选择计算与调和Ritz向量{u}对应的Rayleigh商\vartheta（而不是调和Ritz值）（参见[412]）。Rayleigh商由(A-\tau I){u}-\vartheta {u}\perp {u}而不是\perp {W}的要求得出；然后r\perp {u}。 \overline{{T}_{1,1}}表示{{T}_{1,1}}的复共轭。

(17)

停止准则是在归一化Schur向量近似的残差范数低于\epsilon时接受该近似。这意味着我们接受计算的特征值中的\epsilon量级的误差，以及特征向量中的O(\epsilon)量级的误差（前提是所关注的特征值是简单的且与其他特征值很好地分离，且左右特征向量之间的角度很小）。

无法保证检测到所有所需的特征值；参见算法4.17（第页）的注释(13)。

(20)

这是一个在近似特征对被接受后的重启。

(27)

此时，当子空间的维度超过{m}_{\max}时进行重启。重启后，Jacobi-Davidson迭代以维度为{m}_{\min}的子空间继续进行。

(33)

通过{Q}的因子表示计算的特征向量的缩减。矩阵{Q}的列是收敛的特征向量。如果存在左预条件器{K}用于算子A-\theta I，则可以实现类似于外部特征值情况下的类似缩减。算法4.18给出了高效处理左预条件算子的模板。