舒尔形式与重启

如果我们希望结合重启和缩减技术，情况会变得更加复杂，因为非厄米矩阵通常不具备正交的特征系统。由于我们倾向于至少在测试子空间中使用正交基，因此我们计算缩减矩阵的舒尔形式。我们不是计算矩阵 A 的特征向量，而是计算其部分舒尔形式 A{Q}_{k}={Q}_{k}{R}_{k}，其中 {Q}_{k} 是一个 n 乘 {k} 的正交矩阵，而 {R}_{k} 是一个 {k} 乘 {k} 的上三对角矩阵。一个适应舒尔形式的方案在算法 7.18 中给出。该方案包括当当前子空间的维度超过 {m}_{\max} 时进行重启的可能性。

该方案计算复平面上接近目标 \tau 的 {k}_{\max} 个特征值。在这里我们必须采取不精确的方法，因为一般非厄米矩阵的特征值在复平面上是无序的。哪些 Ritz 值被选为接近的特征值，除了其他因素外，还取决于相应特征向量的角度。通常，如果 \tau 选在特征值分布的外部，从 Ritz 值中进行选择是合适的。如果我们希望计算接近某个内部点的 {M} 的特征值，那么该方案可能效果不佳，我们建议在这种情况下使用调和 Ritz 值。计算内部特征值的方案将在 §7.12.3 中介绍。

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

(1) t = v_{0}, k = 0, m = 0; Q = [], R = []
(2) 当 k
(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
(6)   对于 i = 1, \ldots, m-1
(7)     M_{i, m} = v_i^* v_m^A, M_{m, i} = v_m^* v_i^A
(8)   M_{m, m} = v_{m}^{*} v_{m}^{A}
(9) 进行 Schur 分解 M = S T S^{*}，使得 \left\|T_{i, i} - \tau\right\| \leq \left\|T_{i+1, i+1} - \tau\right\|
(10) u = V s_{1}, u^{A} = V^{A} s_{1}, \theta = T_{1,1}, r = u^{A} - \theta u, \tilde{a} = Q^{*} r, \tilde{r} = r - Q\tilde{a}
(11) 当 \|\tilde{r}\|_{2} \leq \epsilon_{M}
(12)   R = \begin{bmatrix}R & \tilde{a} \\ 0 & \theta\end{bmatrix}, Q = [Q, u], k = k + 1
(13)   如果 k = k_{\max} 则停止
(14)   m = m - 1
(15)   对于 i = 1, \ldots, m
(16)     v_{i} = V s_{i+1}, v_{i}^{A} = V^{A} s_{i+1}, s_{i} = e_{i}
(17)   M = \text{T的下m乘m块}
(18)   u = v_{1}, \theta = M_{1,1}, r = v_{1}^{A} - \theta u, \tilde{a} = Q^{*} r, \tilde{r} = r - Q\tilde{a}
(19) 结束当
(20) 如果 m \geq m_{\max} 则
(21)   对于 i = 2, \ldots, m_{\min}
(22)     v_{i} = V s_{i}, v_{i}^{A} = V^{A} s_{i}
(23)   M = T 的前 m_{\min}\times m_{\min}块
(24)   v_{1} = u, v_{1}^{A} = u^{A}, m = m_{\min}
(25) 结束如果
(26) \widetilde{Q} = [Q,u]
(27) 解 t(\perp \widetilde{Q}) （近似地）从：\left(I - \widetilde{Q} \widetilde{Q}^{*}\right)(A - \theta I)\left(I - \widetilde{Q} \widetilde{Q}^{*}\right) t = -\tilde{r}
(28) 结束当

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

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

现在我们将讨论算法的组成部分。

(1)

初始化阶段。

(3)-(4)

新向量 {t} 通过改进的 Gram-Schmidt 方法与当前搜索子空间正交。为了提高数值稳定性，可以使用算法 4.14 中给出的迭代改进的 Gram-Schmidt 模板进行替换。

如果 {m}=0，则这是一个空循环。

如果新向量 {t} 与搜索子空间之间的角度非常小，则生成的向量 {v}_j 几乎没有意义，方法实际上会停滞。一个随机向量 {t} 有助于克服这种停滞。更复杂的策略可以在 [190] 中找到。

(6)-(8)

我们计算稠密矩阵的最后一列和行 {M}\equiv V^\ast AV=V^\ast V^A（阶数为 m）； {V^A}\equiv A{V}。矩阵 {V} 表示列向量为 {v}_j 的 n 乘 {m} 矩阵，{V^A} 同样。

(9)

{m} 乘 {m} 矩阵 {M} 的舒尔约简可以通过标准稠密特征问题求解器解决。我们选择计算标准 Ritz 值，这使得算法适合计算接近指定 \tau 的 {k}_{\max} 个外部特征值。如果必须计算谱内部部分的特征值，则建议使用调和 Ritz 值；参见 §7.12.3。

在每一步中，舒尔形式必须排序，使得 {T}_{1,1} 最接近 \tau 。只有当 {m}\geq {m}_{\max} 时，舒尔形式的排序才必须使得 {T} 的所有 {m}_{\max} 个前导对角元素最接近 \tau 。为了便于表示，我们在这里对所有对角元素进行了排序。

有关舒尔形式排序算法的模板，请参见 [448,449,33,171, Chap. 6B]。

S 是一个 m 乘 m 的矩阵，列向量为 s_j。

(11)

停止准则是当残差范数（对于归一化舒尔向量近似）低于 \epsilon 时接受舒尔向量近似。这意味着我们接受计算特征值中的不精确度在 \epsilon 数量级，舒尔向量中的不精确度（在角度上）为 O(\epsilon)（前提是所关注的特征值是简单的且与其他特征值充分分离）。 对于更定量的分析和比例常数，参见 §7.13。

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

如果矩阵是实的，则所有特征对要么是实的，要么以复共轭对出现。如果检测到复特征对，则其共轭已知，也可以进行缩减。这使得算法更高效。

(15)-(16)

在接受一个 Ritz 对后，我们继续搜索下一个对，剩余的 Ritz 向量作为初始搜索空间的基。

(20)

一旦当前舒尔向量的搜索空间维度超过 {m}_{\max}，我们就进行重启。重启过程以最接近目标值 \tau 的 {m}_{\min} 个 Ritz 向量所张成的子空间进行。

(27)

我们将锁定的（已计算的）舒尔向量收集在 {Q} 中，矩阵 \tilde{Q} 是 {Q} 扩展了当前舒尔向量近似 {u}。这样做是为了获得更紧凑的公式；校正方程等价于 (4.50) 中的方程。新校正 {t} 必须与 {Q} 的列以及 {u} 正交。

当然，校正方程可以通过任何合适的流程来求解，例如设计用于求解非对称系统的预处理 Krylov 子空间方法。由于 \widetilde{Q} 的出现，在使用预处理矩阵 A-\theta I 时必须小心。预处理器的包含可以遵循与单向量 Jacobi-Davidson 算法相同的原则；参见算法 4.18 的模板。确保迭代求解器的初始向量 {t}_0 满足正交约束 \widetilde{Q}^\ast {t}_0=0。注意，如果 {K} 在几次 Jacobi-Davidson 迭代中保持不变，算法 4.18 可以在每一步中节省显著的计算量。在这种情况下，可以从先前步骤中保存 \widehat{Q} 的列。同样，矩阵 {\mathcal M} 及其 {\mathcal L}{\mathcal U} 分解也可以从先前步骤中更新。

没有必要非常精确地求解校正方程。一种常用于不精确牛顿方法的策略 [113]，在这里也很有效：随着 Jacobi-Davidson 迭代步骤增加精度，例如在校正方程中实现残差减少 2^{-\ell} 在第 \ell 次 Jacobi-Davidson 迭代中（当检测到舒尔向量时，\ell 重置为 0）。

特别是在最初的几次初始步骤中，近似特征值 \theta 可能非常不准确，此时精确求解校正方程没有意义。在这个阶段，暂时将 \theta 替换为 \tau 或在搜索子空间扩展中取 {t}=-r 可能更有效 [335,172]。

关于缩减，参见 §8.4.2，其中 Z_j 替换为 Q_j，B 替换为 I。 关于该方法的完整理论背景以及舒尔向量的缩减技术细节，参见 [172]。