不精确矩阵变换

这里仅对Cayley变换进行了分析，尽管移位-反演方法也可以使用。原因有三：首先，移位-反演Arnoldi方法和Cayley变换Arnoldi方法产生的Ritz向量相同。其次，这使得与(Jacobi)Davidson方法的联系更容易建立。第三，Cayley变换为问题提供了一种更自然的处理方式。

在应用于Cayley变换的Arnoldi方法中，T_{\mathrm{C}}，我们必须解一系列线性系统

(A - \mu B) w_j = (A - \nu B) v_j - s_j

对于j=1,\ldots,k，其中s_j是残差，v_j是最后一个Krylov基向量。结果w_j相对于v_1,\ldots,v_j进行正交化，得到v_{j+1}。也可以将Cayley变换应用于Krylov空间中的另一个向量，而不是v_j。在这种情况下，我们解线性系统

(A - \mu B) w_j = (A - \nu B) y_j - s_j, \tag{11.2}

其中y_j选择如下：它可以是一个Ritz向量(Jacobi-Davidson)或最后一个基向量(Arnoldi)或连续向量(有理Krylov)。残差范数\Vert s_j\Vert是后向误差的度量，因为我们可以将(11.2)重写为

(A - \mu B + s_j w_j^{\ast} / \Vert w_j\Vert^2 ) w_j= (A - \nu B) y_j .

当使用直接方法时，通常会发生 s_j/ \Vert A-\mu B\Vert \Vert w_j\Vert是机器精度的适度倍数。我们称直接方法计算了线性系统的后向稳定解。由于\Vert s_j\Vert非常小，我们谈论的是精确的特征值求解器。通过迭代方法获得小的残差范数也是可能的，但通常非常昂贵。在本节中，我们研究\Vert s_j\Vert较大的情况。为了说明我们所说的“大”是什么意思，需要对迭代线性系统求解器说几句话。线性系统Cx=b被称为以相对残差容差\tau求解，当解x满足\Vert b-Cx\Vert \leq \tau\Vert b\Vert时。Krylov方法[388]通常被使用。GMRES[389]，\mathrm{BiCGSTAB}(\ell)[410]，和QMR[179]是其中最广泛使用的。参见[41]以获取所有这些求解器的模板。当使用合适的预处理器时，它们的性能显著提高。因此，我们所说的“大误差”是指 10^{-8} \leq \tau \leq 10^{-2}.

通过将所有s_j对于j=1, \ldots,k-1放在一起， S_{k-1}=[s_{1},\ldots,s_{k-1}]，w_j在W_{k-1}中，y_j在Y_{k-1}中，我们得到

(A - \mu B)W_{k-1} = (A -\nu B) Y_{k-1} - S_{k-1} . \tag{11.3}

令W_{k-1}^{\dagger}为W_{k-1}的广义Moore-Penrose逆，即W_{k-1}^{\dagger} W_{k-1} = I，并定义 E_{k-1} = S_{k-1} W_{k-1}^{\dagger}。我们可以将关系重写为

(A - \mu B + E_{k-1}) W_{k-1} = (A - \nu B) Y_{k-1} . \tag{11.4}

这是我们在{k-1}步后应得到的“非精确”Cayley变换的关系

T_{\mathrm{IC}} = (A - \mu B + E_{k-1})^{-1} (A - \nu B).

注意，对于固定的\nu和\mu，T_{\mathrm{IC}}依赖于k以及每个线性系统(11.2)的求解方式。当(11.2)通过预处理器M或静态线性系统求解器求解时，即 w_j = M^{-1} (A-\nu B) y_j对于j=1, \ldots,k-1，则 E_{k-1} = E \equiv M - (A-\mu B)独立于k和y_j，因此T_{\mathrm{IC}}也是。

特征对从{\rm span}\{v_1,v_2,\ldots,v_k\}计算。在Arnoldi方法中，这是通过W_{k-1}的正交化产生的Hessenberg矩阵发生的；在Davidson方法中，使用A和B的投影，例如 V_k^{\ast} A V_k z = \theta V_k^{\ast} B V_k z。

{\rm Span}(v_1,\ldots,v_k)包含与T_{\mathrm{IC}}的分离极值特征值相关的特征向量。这是一个扰动的T_{\mathrm{C}}，因此与Ax = \lambda Bx的关系部分丧失，准确计算Ax = \lambda Bx的特征对可能很困难。为了了解哪些参数在扰动中起作用，从[323]中的定理4.3和4.4可以证明，如果T_{\mathrm{IC}}有不同的特征值，那么对于T_{\mathrm{C}}的每个特征对(\zeta,x)，存在T_{\mathrm{IC}}的特征对(\eta,u)，使得

\begin{aligned} \vert\zeta - \eta\vert \leq \kappa_1 \vert\zeta\vert \Vert E_{k-1}\Vert, \\ \Vert(I-uu^*/(u^*u)) x\Vert \leq \kappa_2 \vert\zeta\vert \Vert E_{k-1}\Vert \, , \end{aligned}

其中\kappa_1,\kappa_2 > 0。这些不等式表明，如果\Vert E_{k-1}\Vert和/或\zeta很小，(\zeta,x)和(\eta,u)对应得非常好。在这种情况下，(\zeta,x)可以通过T_{\mathrm{IC}}的特征对计算。由于\zeta = (\lambda-\mu)^{-1}(\lambda-\nu)，\zeta\approx 0，简化为\lambda \approx \nu。

本节可以总结如下。

• 非精确Cayley变换表现得像精确Cayley变换。这意味着仍然针对极点\mu附近的特征值。

• 然而，特征值受到扰动。扰动的程度取决于E_{k-1}。对于任何E_{k-1}，变换能够准确计算\nu附近的特征值。在实践中，这意味着只能非常准确地计算一个（简单）特征值。这也意味着\nu需要不时更新，因为\lambda是未知的。其他特征值仅根据误差水平\Vert E_{k-1}\Vert近似计算。

• 序列(11.2)对于j=1, \ldots,k-1定义了一个Krylov子空间\mathcal{K}^k(T_{\mathrm{IC}}, y_1)。

剩下的问题是现在如何计算T_{\mathrm{IC}}的特征对，或者如何利用它们计算Ax = \lambda Bx的特征对。在第11.2.3节和第11.2.4节中，我们讨论Rayleigh-Ritz技术；即，特征对从Ax = \lambda Bx在{\rm span}(v_1,\ldots,v_k)上的正交投影计算。在第11.2.7节中，特征对直接通过使用有理Krylov递推关系的T_{\mathrm{IC}}的特征对计算。第11.2.6节介绍了一种使用特征向量的递推关系和特征值的Rayleigh-Ritz投影的Lanczos算法。

注意，\mu和\nu不能选择得太远。假设需要特征值\lambda。从上面的结论来看，当\mu选择接近\lambda时，收敛更快，并且计算的特征值只有在\nu接近\lambda时才能准确。理论上，可以选择\mu=\nu，这通常在Davidson方法中使用。由于在这种情况下T_{\mathrm{C}}=I，后一种方法存在停滞的高风险。在第11.2.5节中讨论的Jacobi-Davidson方法中使用了鲁棒选择。