Jacobi-Davidson方法
 G. Sleijpen 和 H. van der Vorst

我们考虑将Jacobi-Davidson方法应用于GHEP（5.1）。类似于上一节中对Lanczos方法的处理（参见[444]），我们可以将Jacobi-Davidson方法应用于（5.1），其中搜索子空间的基是B-正交的，即

{V}_{m}^\ast B{V}_{m}=I_{m},

如果我们用{V}_{m}表示列向量为{v}_j的矩阵。

在这个子空间中，向量{u}\equiv {V}_{m}{s}的Ritz-Galerkin条件导致

({V}_{m}^\ast A {V}_{m}) {s} - \theta ({V}_{m}^\ast B {V}_{m}) {s} = 0, \tag{5.24}

或者，由于{V}_{m}的B-正交性：

({V}_{m}^\ast A {V}_{m}) {s} - \theta {s} = 0 .

这导致了Ritz向量{u}_j^{(m)}\equiv V_{m}{s}_j和Ritz值\theta_j^{(m)}。我们假设这些Ritz向量是相对于B-内积归一化的。

对于广义特征问题的特征向量分量t_j^{({m})}\perp B{u}_j^{(m)}的校正方程可以写成

\begin{aligned} \left(I - B{u}_j^{(m)} {{u}_j^{(m)}}^\ast\right)(A-\theta^{(m)}_jB)\left(I - {u}_j^{(m)} {{u}_j^{(m)}}^\ast B\right) t^{(m)}_j \\ = -r_j^{({m})}\equiv - (A-\theta_j^{(m)} B){u}_j^{(m)} . \end{aligned}.\tag{89}

如果矩阵A-\theta_j^{(m)} B的线性系统可以高效求解，那么我们可以计算精确解t_j^{({m})}。在其他情况下，我们可以像通常对Jacobi-Davidson方法所做的那样，使用应用于（5.25）的Krylov求解器计算t_j^{({m})}的近似解。注意，算子

\left( I - B{u}_j^{(m)} {{u}_j^{(m)}}^\ast \right)(A-\theta_j^{(m)} B)\left( I - {u}_j^{(m)} {{u}_j^{(m)}}^\ast B \right)

将空间(B{u}_j^{(m)})^\perp映射到空间({u}_j^{(m)})^\perp，因此如果我们在(B{u}_j^{(m)})^\perp和(B{u}_j^{(m)})^\perp之间使用Krylov求解器，则总是需要预处理（参见算法5.6中的项(31)）。

对于标准厄米特情况，所得到的方案可以与重启和消减相结合。如果我们希望在消减中使用正交算子，那么我们必须使用将给定广义系统化为Schur形式的B-正交矩阵：

A{Q}_{k}={Z}_{k}{D}_{k},

其中{Z}_{k}=B{Q}_{k}且{Q}_{k}是B-正交的。矩阵{D}_{k}是一个对角矩阵，其对角线上是计算出的{k}个特征值；{Q}_{k}的列是A的特征向量。这导致了对消减的前{k}个特征向量的斜投影：

\left( I - {Z}_{k} {Q}_{k}^\ast \right)(A-\lambda B)\left( I - {Q}_{k} {Z}_{k}^\ast \right) \, .

很容易验证，相对于空间(B{u}_j^{(m)})^\perp，消减后的算子B仍然是关于该空间的对称正定。我们可以简单地在该空间中使用B-内积，因为B和消减后的B在该空间上一致。

如果B不是良态的，即如果x^\ast By导致高度扭曲的内积，那么我们建议使用Jacobi-Davidson方法的QZ方法（参见第8.4节）。QZ方法不利用所涉及矩阵的对称性。

算法5.6表示一个用于外部特征值的重启和消减的Jacobi-Davidson模板。算法5.7给出了一个左预处理的Krylov求解器的模板。

算法 5.6: 用于GHEP的k_{\text{max}}外部特征值的Jacobi-Davidson方法

(1)  t = v_0, k = 0, m = 0; Q = [], Z = []
(2)  当 k < k_{\text{max}} 时
(3)      对于 i = 1, \ldots, m
(4)          t = t - (v_i^B * t) v_i
(5)      结束 for 循环
(6)      m = m + 1; \tilde{t} = B t, \text{norm} = \sqrt{t^* \tilde{t}},
(7)      v_m = t / \text{norm}, v_m^A = A v_m, v_m^B = \tilde{t} / \text{norm}
(8)      对于 i = 1, \ldots, m
(9)          M[i, m] = v_i^* v_m^A
(10)     结束 for 循环
(11)     计算 M 的所有特征对 (\theta_i, s_i) (其中 |s_i|_2 = 1)，并按 |\theta_i - \tau| \geq |\theta_{i-1} - \tau| 排序
(12)     u = V s_1, p = V^B s_1, u^A = V^A s_1, r = u^A - \theta_1 p
(13)     当 \frac{|r|_2}{|u|_2} \leq \epsilon 时
(14)         Q = [Q, u], Z = [Z, p], \tilde{\lambda}_{k+1} = \theta_1, k = k + 1
(15)         如果 (k = k_{\text{max}}), \tilde{X} = Q, 则停止
(16)         m = m - 1, M = 0
(17)         对于 i = 1, \ldots, m
(18)             v_i = V s_{i+1}, v_i^A = V^A s_{i+1}, v_i^B = V^B s_{i+1}
(19)             M[i, i] = \theta_{i+1}, s_i = e_i, \theta_i = \theta_{i+1}
(20)         结束 for 循环
(21)         u = v_1, p = v_1^B, r = v_1^A - \theta_1 p
(22)     结束 while 循环
(23)     如果 m > m_{\text{max}} 则
(24)         M = 0
(25)         对于 i = 2, \ldots, m_{\text{min}}
(26)             v_i = V s_i, v_i^A = V^A s_i, v_i^B = V^B s_i, M[i, i] = \theta_i
(27)         结束 for 循环
(28)         v_1 = u, v_1^B = p, M[1, 1] = \theta_1, m = m_{\text{min}}
(29)     结束 if 语句
(30)     \theta = \theta_1, \tilde{Q} = [Q, u], \tilde{Z} = [Z, p]
(31)     从 (I - \tilde{Z} \tilde{Q}^*)(A - \theta B)(I - \tilde{Q} \tilde{Z}^*) t = -r 中近似求解 t \perp \tilde{Z}
(32) 结束 while 循环

算法 5.7: 降阶 Jacobi-Davidson GHEP 修正方程的近似解法

使用左预条件器 \tilde{K}\equiv\left(I-\tilde{Z}\tilde{Q}^{*}\right) K\left(I-\tilde{Q}\tilde{Z}^{*}\right)，求解 \tilde{A}\equiv\left(I-\tilde{Z}\tilde{Q}^{*}\right)(A-\theta B)\left(I-\tilde{Q}\tilde{Z}^{*}\right)：

(1)  从 K\widehat{Z}=\widetilde{Z} 中求解 \widehat{Z}
(2)  计算 \mathcal{M}=\widetilde{Z}^{*}\widehat{Z}
(3)  分解 \mathcal{M}=\mathcal{L}\mathcal{U}
(4)  计算 \tilde{r}\equiv\tilde{K}_{j}^{-1} r：
(5)      (b') 从 K\widehat{r}=r 中求解 \widehat{r}
(6)      (c') \vec{\gamma}=\tilde{Z}^{*}\hat{r}
(7)             从 \mathcal{L}\vec{\beta}=\vec{\gamma} 中求解 \vec{\beta}
(8)             从 \mathcal{U}\vec{\alpha}=\vec{\beta} 中求解 \vec{\alpha}
(9)      (d') \tilde{r}=\widehat{r}-\widehat{Z}\vec{\alpha}
(10) 使用初始值 t_{0}=0，算子 \widetilde{K}_{j}^{-1}\widetilde{A}，和右端项 -\widetilde{r} 应用 Krylov 子空间方法，计算 z=\widetilde{K}_{j}^{-1}\widetilde{A} v：
(11)     (a) y=(A-\theta B) v
(12)     (b) 从 K\widehat{y}=y 中求解 \widehat{y}
(13)     (c) \vec{\gamma}=\widetilde{Z}^{*} y
(14)           从 \mathcal{L}\vec{\beta}=\vec{\gamma} 中求解 \vec{\beta}
(15)           从 \mathcal{U}\vec{\alpha}=\vec{\beta} 中求解 \vec{\alpha}
(16)     (d) z=\widehat{y}-\widehat{Z}\vec{\alpha}

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

完成后，当\tau选择大于\lambda_{\max}(A)时，将提供{k}_{\max}个最大特征值；如果\tau选择小于\lambda_{\min}，则提供{k}_{\max}个最小特征值。计算出的特征对(\widetilde\lambda_j,\widetilde{x}_{j})，\Vert\widetilde{x}_{j}\Vert _B=1，满足\Vert A \widetilde{x}_{j}- \widetilde\lambda_jB\widetilde{x}_{j}\Vert _2\leq j\epsilon，其中\widetilde{x}_j表示\widetilde{X}的第j列。特征向量是B-正交的：\widetilde{x}_i^\ast B\widetilde{x}_j=0，对于i \neq j。

现在让我们讨论算法5.6的不同步骤。

(1)

初始化阶段。

(3)-(5)

通过修正的Gram-Schmidt方法，使新向量{t}相对于当前搜索子空间B-正交。为了提高数值稳定性，可以用算法4.14中的模板替换，其中所有内积和范数应解释为B-内积，B-范数。

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

(7)

我们扩展n乘m矩阵{V}，{V}^A\equiv A{V}，和{V}^B\equiv B{V}，{V}表示包含当前搜索子空间的基向量{v}_i作为其列的矩阵；同样{V}^A和{V}^B。

(8)-(10)

计算对称矩阵{M}\equiv {V}^\ast A {V}={V}^\ast V^A（阶数{m}）的第{m}列。

(11)

可以用LAPACK中适用于厄米特稠密矩阵的适当例程来解决{m}乘{m}矩阵{M}的特征问题。我们选择计算标准Ritz值，这使得该算法适用于计算接近指定\tau的A - \lambda B的{k}_{\max}个外部特征值。如果需要计算频谱内部部分的特征值，则建议计算调和Petrov值；参见第8.4节。

(13)

停止准则是在归一化向量近似的残差范数低于\epsilon时接受一个特征向量近似。这意味着我们接受在计算的特征值中误差为\epsilon^2，在特征向量中的误差（角度）为O({\epsilon})（前提是所关心的特征值是简单的且与其他特征值很好地分离，并且B不是病态的；参见第5.7节中的（5.33）和（5.34））。

无法保证检测到所有想要的特征值；参见算法4.17（第页）中的项(13)。

(17)-(20)

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

(23)-(29)

一旦当前特征向量的搜索空间维度超过{m}_{\max}，我们就进行重启。该过程以由最接近目标值\tau的{m}_{\min}个Ritz向量张成的子空间重启。该子空间的构建在(25)-(27)中完成。

(31)

我们在{Q}中收集了锁定的（计算出的）特征向量，矩阵\tilde{Q}是Q扩展了当前特征向量近似u。这样做是为了获得更紧凑的公式；算法5.6中的步骤(31)的校正方程等价于方程（5.25）中的校正方程。新的校正{t}必须与{Z}=B{Q}的列以及{p}=B{u}正交。

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

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

有关该方法的完整理论背景，请参见[172]。有关使用特征向量的消减技术的详细信息，请参见第4.7.3节。