不精确有理Krylov方法

在本节中，我们将非精确Cayley变换、(Jacobi) Davidson方法和有理Krylov方法的思想结合起来。粗略地说，我们的方法与非精确Cayley变换Arnoldi方法或预处理Lanczos方法基本相同，唯一的区别在于零点\nu和极点\mu可能在每次迭代中更新。与预处理Lanczos方法一样，Ritz向量是从海森堡矩阵中计算出来的。此外，Ritz值也是从递推关系中计算出来的。

我们现在详细讨论这种方法，使用下面给出的针对非精确Cayley变换设计的有理Krylov方法算法。

算法 11.4：用于广义非线性特征值问题的近似有理Krylov方法

(1) 给定 v_{1}，满足 \left\|v_{1}\right\|=1，和一个零值 \theta_{0}。令残差为 r_{0}=A v_{1}-\theta_{0} B v_{1}。
(2) 对于 j=1,2,\ldots 执行
(3)     选择一个极点 \mu_{j}，令 \nu_{j}=\theta_{j-1}，并解 w_{j}：
          \left(A-\mu_j B\right) w_j=r_{j-1}
(4)     将 w_j 相对于 v_1,\ldots, v_j 正交化，得到 v_{j+1}。在Hessenberg矩阵 H_{j+1, j} 的第j列中保留Gram-Schmidt系数 h_{j}。
(5)     构建 L_{j+1, j} 和 K_{j+1, j} 的第j列：
          \left(h_j-\left[\begin{array}{c} t_j\\ 0\end{array}\right]\right),\left(h_j\mu_j-\left[\begin{array}{c} t_j\\ 0\end{array}\right]\nu_j\right)
        其中 t_1=1，对于 j>1，t_j=L_{j, j-1} z_{j-1}。
(6)    解第j个 j \times j 特征值问题
          L_{j+1, j}^{\dagger} K_{j+1, j} z_{j}=\theta_{j} z_{j}
          并计算Ritz对 \left(\theta_{j}, x_{j}=V_{j+1} L_{j+1, j} z_{j}\right)。
(7)     计算残差 r_{j}=A x_{j}-\theta_{j} B x_{j}。收敛性检查。

让我们逐步分析这个算法。线性系统的解导致关系式（11.2），其中y_j = x_{j-1}是一个Ritz向量，\nu_j是前一次迭代中的相关Ritz值\theta_{j-1}，使得右端是残差 r_{j-1} = A x_{j-1} - \theta_{j-1} B x_{j-1}。 由于 x_{j-1} = V_j L_{j,j-1} z_{j-1}，我们也可以写成 x_{j-1}=V_jt_j，其中 t_j被称为延续向量。经过Gram-Schmidt正交化后，我们有 w_j = V_{j+1} h_j，因此我们可以将（11.2）重写为

A V_{j+1}\left( h_j - \left[\begin{array}{c} t_j \\ 0 \end{array}\right] \nu_j \right)- s_j.

将所有方程收集起来，对于j=1,\ldots,k，我们有

A V_{k+1} L_{k+1,k} = B V_{k+1} K_{k+1,k} - S_k \ , \tag{11.6}

这是另一种写法（11.3）。注意，当从该方程中省略S_k时，我们得到通常的有理Krylov子空间（RKS）递推关系。同样在这种情况下，L_{k+1,k}和K_{k+1,k}是上部海森堡矩阵，并且L_{k+1,k}具有满秩。令 L_{k+1,k}^{\dagger}表示L_{k+1,k}的广义Moore-Penrose逆。利用 E_k = S_k L_{k+1,k}^{\dagger} V_{k+1}^{\ast}，我们也可以写成

(A + E_k) V_{k+1} L_{k+1,k} = B V_{k+1} K_{k+1,k}, \tag{11.7}

这对应于（11.4）。换句话说，RKS系数矩阵L_{k+1,k}和K_{k+1,k}以及向量V_{k+1}可以被视为通过对扰动问题应用精确有理Krylov方法计算得到的结果 (A + E_k) u = \eta B u。类似于“非精确”Cayley变换的术语，我们谈论“非精确”有理Krylov方法。

在我们继续之前，必须给出精确有理Krylov方法的一些性质。以下引理解释了如何在有理Krylov方法中计算Ritz值。

引理 11.1 设 V_{k+1}, L_{k+1, k} 和 K_{k+1, k} 使得 V_{k+1} 和 L_{k+1, k} 具有满秩。设 A V_{k+1} L_{k+1, k}=B V_{k+1} K_{k+1, k} 成立。假设 B 可逆。则 (\theta, x) 是 B^{-1} A x=\lambda x 相对于范围 \left(V_{k+1} L_{k+1, k}\right) 的Ritz对，当且仅当 x_{k}=V_{k+1} L_{k+1, k} z_{k} 且

L_{k+1, k}^{\dagger} K_{k+1, k} z=\theta z.

此外，残差是 A x-\theta B x=f_{k} 其中

f_{k}=B V_{k+1}\left(K_{k+1, k}-\theta L_{k+1, k}\right) z.

与(Jacobi) Davidson方法一样，RKS方法对向量应用非精确Cayley变换。区别在于计算Ritz对的方式。在(Jacobi) Davidson方法中，Ritz对来自于Ax = \lambda Bx在子空间上的Galerkin投影。在RKS方法中，Ritz对是利用引理11.1从递推关系中计算出来的，假设E_k=0。

通过非精确Cayley变换，变换可能对精确Cayley变换造成较大扰动，但只要\nu选择得当，仍可用于计算一个特征对。这里的情况也是如此。不精确有理Krylov方法提供的S_k（或E_k）不是随机的，而是在所需特征对的方向上很小，只要各种参数设置得当。

由算法11.4计算的Ritz对(\theta,x)产生的残差

(A + E_k) x - \theta B x = f_k

由引理11.1定义。与原始问题相关的残差可以写成

A x - \theta B x = f_k - E_k x

或者，根据（11.6），

A x - \theta B x = f_k -S_k z \ .

为了使残差变小，E_k x = S_k z必须趋向于零。关于不精确有理Krylov方法收敛性的数学论证已由Lehoucq和Meerbergen[291]以及De Samblanx[109]给出。我们仅概述理论和一些数值结果。在[291]和[109]中，给出了数学和启发式论证，即 S_j z_j \approx s_j，使得 \Vert s_j\Vert \leq \tau \Vert r_{j-1}\Vert（见第11.2.2节）， 残差满足递推关系

\Vert r_j\Vert \approx \Vert f_j\Vert + \tau \Vert r_{j-1}\Vert \ .

当项f_j对S_j z_j占优势时，j=1,\ldots,k，则 收敛由精确有理Krylov过程的收敛主导： A x - \theta B x \approx f_j。当第二项占优势时，我们有 \Vert r_j\Vert \approx \tau^j\Vert r_0\Vert；即， r_j以收敛比率\tau趋向于零。

• 如果底层精确有理Krylov方法以收敛比率\rho线性收敛，使得 \Vert f_{k+1}\Vert \leq \rho \Vert f_k\Vert，并且 如果使用相对残差容差\tau来求解线性系统，则有效收敛比率大约为 \max(\rho,\tau)。当\mu_j对所有j保持恒定时，获得线性和超线性收敛。

• 当 \mu_j=\theta_{j-1}在每次迭代中，如在Davidson或Jacobi-Davidson方法中，并且 \nu_j=\theta_{j-1}+\epsilon_{j-1} （见第11.2.5节），则收敛是二次的。项f_j以二次方式收敛到零，而第二项在\tau不依赖于j时线性收敛。通常，收敛具有二次行为，直到 S_j z_j成为r_j中的主导项。从那时起，收敛变为线性，收敛比率为\tau。在实践中，建议选择\tau与\Vert r_j\Vert成比例，以保持二次收敛。通过 \tau_j = \varrho \Vert r_{j-1}\Vert，我们得到
  \Vert r_j\Vert \approx \Vert f_j\Vert + \varrho \Vert r_{j-1}\Vert^2,
  这显示了第二项的二次收敛。