算法

Jacobi-Davidson算法在算法8.1中给出。该算法试图计算广义Schur对 ((\alpha,\beta),q)，其中比值\beta/\alpha在复平面上最接近指定的目标值\tau。 算法包含重启机制以限制搜索空间的维度，并通过已经收敛的左右Schur向量进行收缩。

算法 8.1：用于 GNHEP 中接近 \tau 的 k_{\max} 内部特征值的 Jacobi-Davidson QZ 方法
(1) t=v_{0}, k=0, \nu_{0}=1/\sqrt{1+\|\tau\|^{2}}, \mu_{0}=-\tau \nu_{0}, m=0;
(2) 2=[], Z=[], S=[], T=[]
(3) 当 k
(4) 对于 i=1,\ldots, m
(5) t=t-\left(v_{i}^{*} t\right) v_{i}
(6) 2=[], Z=[], S=[], T=[] m=m+1, v_m=t/\|t\|_2, v_m^A=A v_m, v_m^B=B v_m, w=\nu_0 v_m^A+\mu_0 v_m^B
(7) 对于 i=1,\ldots, k
(8) w=w-\left(z_{i}^{*} w\right) z_{i}
(9) 对于 i=1,\ldots, m-1
(10) w_{m}=w/\|w\|_{2}
(11) 对于 i=1,\ldots, m-1
(12) M_{i, m}^{A}=w_{i}^{*} v_{m}^{A}, M_{m, i}^{A}=w_{m}^{*} v_{i}^{A}, M_{i, m}^{B}=w_{i}^{*} v_{m}^{B}, M_{m, i}^{B}=w_{m}^{*} v_{i}^{B}
(13) M_{m, m}^{A}=w_{m}^{*} v_{m}^{A}, M_{m, m}^{B}=w_{m}^{*} v_{m}^{B}
(14) 计算 QZ 分解：M^{A} S^{R}=S^{L} T^{A}, M^{B} S^{R}=S^{L} T^{B} 使得：\left\|T_{i, i}^{A}/ T_{i, i}^{B}-\tau\right\|\leq\left\|T_{i+1, i+1}^{A}/ T_{i+1, i+1}^{B}-\tau\right\|
(15) u=V s_1^R, p=W s_1^L, u^A=V^A s_1^R, u^B=V^B s_1^R,\zeta=T_{1,1}^A,\eta=T_{1,1}^B
(16) r=\eta u^{A}-\zeta u^{B},\tilde{a}=Z^{*} u^{A},\tilde{b}=Z^{*} u^{B},\tilde{r}=r-Z(\eta\tilde{a}-\zeta\tilde{b})
(17) 当 \|\tilde{r}\|_{2}\leq\epsilon
(18) R^{A}=\left[\begin{array}{cc}R^{A}&\widetilde{a}\\ 0&\zeta\end{array}\right],\quad R^{B}=\left[\begin{array}{cc}R^{B}&\widetilde{b}\\ 0&\eta\end{array}\right]
(19) Q=[Q, u], Z=[Z, p], k=k+1
(20) 如果 k=k_{\max} 则停止
(21) m=m-1
(22) 对于 i=1,\ldots, m
(23)  v_i=V s_{i+1}^R, v_i^A=V^A s_{i+1}^R, v_i^B=V^B s_{i+1}^R,
(24) w_{i}=W s_{i+1}^{L}, s_{i}^{R}=s_{i}^{L}=e_{i}
(25) M^{A}, M^{B} 是 T^{A}, T^{B} 的下 m 阶 m 阶块，分别。
(26) u=v_{1}, p=w_{1}, u^{A}=v_{1}^{A}, u^{B}=v_{1}^{B},\zeta=T_{1,1}^{A},\eta=T_{1,1}^{B}
(27) r=\eta u^{A}-\zeta u^{B},\widetilde{a}=Z^{*} u^{A},\widetilde{b}=Z^{*} u^{B},\widetilde{r}=r-Z(\eta\widetilde{s}-\zeta\widetilde{t})
(28) 结束当
(29) 如果 m\geq m_{\max} 则
(30) 对于 i=2,\ldots, m_{\min}
(31) v_{i}=V s_{i}^{R}, v_{i}^{A}=V^{A} s_{i}^{R}, v_{i}^{B}=V^{B} s_{i}^{R}, w_{i}=W s_{i}^{L}
(32) M^{A}, M^{B} 是 T^{A}, T^{B} 的前 m_{\min} 阶 m_{\min} 阶块，分别。
(33) v_1=u, v_1^A=u^A, v_1^B=u^B, w_1=p, m=m_{\min}
(34) 结束如果
(35) \widetilde{Q}=[Q, u],\widetilde{Z}=[Z, p]
(36) 解 t(\perp\widetilde{Q})  (大致) 来自
(37) \left(I-\widetilde{Z}\widetilde{Z}^{*}\right)(\eta A-\zeta B)\left(I-\widetilde{Q}\widetilde{Q}^{*}\right) t=-\widetilde{r}
(38) 结束当

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

完成后，将交付k_{\max}个接近\tau的广义特征值，以及相应的缩减Schur形式AQ=ZR^A，BQ={Z}R^B，其中Q和Z是n乘k_{\max}的正交矩阵，R^A和R^B是k_{\max}乘k_{\max}的上三角矩阵。广义特征值是R^A和R^B的对角线元素。计算形式满足 \Vert A q_{j}- {Z} R^Ae_{j}\Vert _2=O(\epsilon)， \Vert B q_{j}- {Z} R^Be_{j}\Vert _2=O(\epsilon)， 其中q_j是Q的第j列。

计算的缩减Schur形式的精度取决于目标值\tau与特征值 (\alpha_j,\beta_j)\equiv(R^A_{j,j},R^B_{j,j})之间的距离。 如果我们忽略机器精度和\epsilon^2阶的项，那么我们有 \Vert A q_j- {Z} R^Ae_j\Vert _2\leq j\gamma_A \epsilon， \Vert B q_j- {Z} R^Be_j\Vert _2\leq j\gamma_B \epsilon， 其中常数 \gamma_A和\gamma_B由

\gamma_A\equiv\frac{\vert\mu_0\vert}{\vert\nu_0\alpha_j+\mu_0\beta_j} \quad\text{和}\quad \gamma_B = \frac{\vert\nu_0\vert}{\vert\nu_0\alpha_j+\mu_0\beta_j\vert}.

如果 \mu_0/\nu_0=-\tau，如算法步骤(1)中所示，那么 \gamma_A=\vert\tau\vert/\vert\alpha_j-\tau\beta_j\vert， \gamma_B=1/\vert\alpha_j-\tau\beta_j\vert。这些值可能很大 如果 \tau\approx \alpha_j/\beta_j。实际上，如果\tau接近检测到的特征值，也能达到\epsilon阶的精度。当使用值 (\mu_0,\nu_0)=(1,\bar\tau)进行额外的细化步骤时，可以保证\epsilon阶的精度。

现在我们将解释算法的主要连续阶段。

(1)

初始化阶段。 标量\nu_0和\mu_0的选择特别有效 如果\tau在谱的内部。 如果\tau是特征值（可以轻松测试），此选择会导致分解失败。

(4)-(5)

新向量t通过修改的Gram-Schmidt方法相对于当前搜索子空间V进行正交化。 同样，向量 w=(\nu_0A+\mu_0B)t相对于当前测试子空间W进行正交化。这两个正交化过程可以通过算法4.14（第页）中的模板替换，以提高数值稳定性。

我们扩展子空间V，V^A\equiv AV，V^B\equiv BV，以及W。 V表示包含当前搜索子空间基向量v_i作为列的矩阵。其他矩阵以类似明显的方式定义。

(12)-(13)

计算矩阵 M^A\equiv W^\ast A V 和 M^B\equiv W^\ast B V的第{m}行和列。

注意，标量M_{i,{m}}^B也可以从标量M_{i,{m}}^A以及步骤(10)中w_i^\ast w的正交化常数计算得出。

(15)

可以通过LAPACK中适用于密集矩阵束的合适例程计算{m}乘{m}矩阵对(M^A, M^B)的QZ分解。

我们选择计算广义Petrov对，这使得算法适合于计算k_{\max}个内部广义特征值 \beta A-\alpha B，其中\alpha/\beta接近指定的\tau。

有关重新排序广义Schur形式的算法，请参见[448,449,171]。

(18)

停止准则是在归一化右Schur向量近似残差的范数低于\epsilon时接受广义本征对近似。这意味着我们在计算的广义特征值中接受\epsilon阶的不精确性，以及Schur向量中O({\epsilon})（角度）的不精确性（前提是所关心的特征值是简单的且与其他特征值充分分离）。

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

(22)

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

(30)

当当前特征向量的搜索空间的维度超过{m}_{\max}时，我们进行重启。重启过程以{m}_{\min}个左右Ritz向量为基础，这些Ritz向量对应于最接近目标值\tau的广义Ritz对。

(36)

我们将已锁定的（计算的）右Schur向量收集在{Q}中，矩阵\tilde{Q}是{Q}扩展了当前右Schur特征向量近似{u}。同样，收敛的左Schur向量已被收集在{Z}中，该矩阵扩展了{p}。这样做是为了获得更紧凑的表述；算法8.1步骤(37)中的修正方程等价于方程(8.13)中关于收缩对(8.15)的修正方程。新修正{t}必须与{Q}的列以及{u}正交。

当然，修正方程可以通过任何合适的进程求解，例如，设计用于求解非对称系统的预处理Krylov子空间方法。然而，由于不同的投影，我们总是需要一个预处理器（如果没有任何其他可用，则可能是恒等算子），该预处理器通过相同的斜投影进行收缩，以便我们在 \widetilde{Q}^\perp和自身之间获得映射。 由于\widetilde{Q}和\widetilde{Z}的出现，必须谨慎使用矩阵 \eta A -\zeta B的预处理器。 预处理器的包含可以按照算法8.2进行。确保迭代求解器的初始向量t_0满足正交约束 \widetilde Q^\ast t_0=0。请注意，如果在几次Jacobi-Davidson迭代中保持{K}不变，算法8.2每步可以节省显著的计算量。在这种情况下，可以从先前步骤中保存\widehat{Z}的列。同样，矩阵{\mathcal M}及其 {\mathcal L}{\mathcal U}分解也可以从先前步骤更新。

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

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

有关此方法的完整理论背景以及Schur向量的收缩技术细节，请参见[172]。

算法 8.2：Deflated Jacobi-Davidson GNHEP修正方程的近似解

使用左预处理器 \tilde{K}\equiv\left(I-\tilde{Z}\tilde{Z}^{*}\right) K\left(I-\tilde{Q}\tilde{Q}^{*}\right) 来求解 \tilde{A}\equiv\left(I-\widetilde{Z}\tilde{Z}^{*}\right)(\eta A-\zeta B)\left(I-\widetilde{Q}\tilde{Q}^{*}\right) ：

(1) 求解 \widehat{Z} 使得 K\widehat{Z}=\widetilde{Z}，
(2) 计算 \mathcal{M}=\widetilde{Q}^{*}\widehat{Z}，
(3) 分解 \mathcal{M}=\mathcal{L}\mathcal{U}，
(4) 计算 \tilde{r}\equiv\tilde{K}^{-1} r 如下：
(5) \quad\left(b^{\prime}\right) 求解 \widehat{r} 使得 K\widehat{r}=r，
(6) \left(c^{\prime}\right)\vec{\gamma}=\widetilde{Q}^{*}\vec{r}
(7) 求解 \vec{\beta} 使得 \mathcal{L}\vec{\beta}=\vec{\gamma}
(8) 求解 \vec{\alpha} 使得 \mathcal{U}\vec{\alpha}=\vec{\beta}，
(9) (d') \tilde{r}=\widehat{r}-\widehat{Z}\vec{\alpha}
应用 Krylov 子空间方法，初始值 t_{0}=0，操作符 \widetilde{K}^{-1}\widetilde{A}，右侧为 -\tilde{r}，给定的 v 计算 z=\widetilde{K}^{-1}\widetilde{A} v 如下：
(a) y=(\eta A-\zeta B) v，
(b) 求解 \widehat{y} 使得 K\widehat{y}=y，
(c) \vec{\gamma}=\widetilde{Q}^{*}\widehat{y}
求解 \vec{\beta} 使得 \mathcal{L}\vec{\beta}=\vec{\gamma}，
求解 \vec{\alpha} 使得 \mathcal{U}\vec{\alpha}=\vec{\beta}，
(d) z=\widehat{y}-\widehat{Z}\vec{\alpha}。