分数Krylov子空间方法
 A. Ruhe

分数Krylov子空间方法用于计算正则束的特征值和特征向量，

(A-\lambda B)x=0.\tag{8.17}

它是位移-反演Arnoldi算法的推广，在一次运行中使用多个不同位移的因子分解。这样，我们可以很好地近似于所有选定位移周围区域的特征值；参见[378]。一个典型的应用是线性动态系统的模型降阶，其中希望在预定频率范围内获得响应；例如，参见[379]或[184]。

分数Krylov开始时是一个带有位移\sigma_1和初始向量v_1的位移和反演Arnoldi迭代。它使用基本递归计算正交基V_j，

(A-\sigma_1B)^{-1}BV_j=V_{j+1}H_{j+1,j}\;. \tag{8.18}

（下标表示矩阵的大小，H_{j+1,j}是(j+1) \times j，而V_j是一个每列长度为n的j列矩阵。）我们可以从H_{j,j}计算Ritz值，并以通常的方式从V_j计算Ritz近似特征向量。解这个小规模的特征值问题

H_{j,j}Z^{(j)}=Z^{(j)}\Theta^{(j)} \tag{8.19}

并得到Ritz值\lambda^{(j)}_i=\sigma_1+1/\theta_i^{(j)}和Ritz向量y^{(j)}_i=V_jz_{.,i}^{(j)}，其中i=1,\dots,j。经过适度的步数j后，我们可以很好地近似于最接近位移\sigma_1的那些特征值\lambda_i。

分数Krylov算法的思想是在我们希望获得接近新位移的特征值近似时，继续使用一个新的位移\sigma_2。棘手的是避免丢弃在旧位移\sigma_1下计算得到的基V_j中收集的所有信息。如果我们将基V_{j+1}替换为一个新的基W_{j+1}，它跨越与V_{j+1}相同的子空间，但可以解释为带有新位移\sigma_2的Arnoldi递归(8.18)的正交基，这是可能的。同时，梯形（海森堡）矩阵H_{j+1,j}被替换为一个相同形式的新矩阵\tilde{H}_{j+1,j}。

将递归(8.18)重写为

BV_j=(A-\sigma_1B)V_{j+1}H_{j+1,j}\;.

在左侧添加一个矩阵，使得右侧的\sigma_1被新位移\sigma_2替换，

(\sigma_1-\sigma_2)BV_{j+1}H_{j+1,j}+BV_j=(A-\sigma_2B)V_{j+1}H_{j+1,j}\;,

并注意到B是左侧唯一的大矩阵，

BV_{j+1}(I_{j+1,j}+(\sigma_1-\sigma_2)H_{j+1,j})=(A-\sigma_2B)V_{j+1}H_{j+1,j}\;.

左侧的矩阵K_{j+1,j}=I_{j+1,j}+(\sigma_1-\sigma_2)H_{j+1,j}是梯形的，与H_{j+1,j}具有相同的非零模式，我们可以通过预乘(A-\sigma_2B)^{-1}回到类似于递归(8.18)的表达式，

(A-\sigma_2B)^{-1}BV_{j+1}K_{j+1,j}=V_{j+1}H_{j+1,j}\;.

剩下的就是去掉左侧的因子K_{j+1,j}，这是通过以下方式完成的。对K_{j+1,j}进行QR分解并得到

(A-\sigma_2B)^{-1}BV_{j+1}Q_{j+1,j+1}\left[\begin{array}{c}{R_{j,j}}\\ 0\end{array}\right]=V_{j+1}H_{j+1,j}\;.

从右侧乘以三角矩阵R_{j,j}的逆，我们知道如果海森堡矩阵是未简化的，它是非奇异的，

(A-\sigma_2B)^{-1}BV_{j+1}Q_{j+1,j}=V_{j+1}H_{j+1,j}R_{j,j}^{-1}\;. \tag{8.20}

引入新的正交基V_{j+1}Q_{j+1,j+1}，并注意到束(8.17)在新的基下由全矩阵表示，

L_{j+1,j}=Q_{j+1,j+1}^{\ast}H_{j+1,j}R_{j,j}^{-1},

我们可以通过从底部向上应用Householder算法将这个L_{j+1,j}转换为梯形（上部海森堡）形式（参见[377]），

L_{j+1,j}= \left[\begin{array}{cc}P_j & 0 \\ 0&1\end{array}\right]\tilde{H}_{j+1,j}P_j^{\ast}\;.

从右侧乘以正交的P_j并得到

(A-\sigma_2B)^{-1}BV_{j+1}Q_{j+1,j}P_j=V_{j+1}Q_{j+1}\left[\begin{array}{cc}P_j & 0 \\ 0& 1\end{array}\right]\tilde{H}_{j+1,j}

或

(A-\sigma_2B)^{-1}BW_{j}=W_{j+1}\tilde{H}_{j+1,j},

一个带有新位移\sigma_2的Arnoldi递归(8.18)，新的正交基，

W_{j+1}=V_{j+1}Q_{j+1,j+1}\left[\begin{array}{cc}P_j & 0 \\ 0& 1\end{array}\right]\;,

和新的梯形（上部海森堡）矩阵\tilde{H}_{j+1,j}。

我们再次强调，所有这些变换都是在不实际对大n \times n矩阵A和B进行任何操作的情况下实现的。我们甚至可以积累所有的Q和P因子，避免显式形成W，从而避免对长向量的额外工作。还要注意，即使我们得到一个Arnoldi递归，它也不从原始初始向量v_1开始，而是从w_1开始，后者是整个分数Krylov算法中计算的所有向量的线性组合。

在实际实现中，这与锁定、清除和隐式重启相结合。首先用第一个位移\sigma_1运行位移和反演Arnoldi。当\sigma_1周围适当数量的特征值收敛时，比如经过m步，锁定这些收敛的特征值并清除那些完全在感兴趣区域之外的特征值，留下一个(j \leq m)步的Arnoldi递归(8.18)。然后引入新的位移\sigma_2，并按照本节的描述获取新的基W_{j+1}，替换V_{j+1}。从新的位移开始，对这个新基的最后一个向量进行操作

r:=(A-\sigma_2B)^{-1}Bv_{j+1}

并得到带有新位移\sigma_2的Arnoldi递归(8.18)中的下一个基向量v_{j+2}。继续直到我们在\sigma_2周围的特征值集合收敛，并重复新的位移过程，直到所有感兴趣的特征值都收敛或预定的频率范围内的所有位移都已使用。

算法 8.3：GNHEP的有理Krylov子空间方法
(1) 以 \sigma_{1} 作为起始位移，v_{1} 作为单位起始向量，基础大小 j=1
(2) 对于 i=1,2,... 直到收敛
(3) 将 j 步Arnoldi递归扩展到 m 步：(A-\sigma_i B)^{-1} B V_m=V_{m+1} H_{m+1, m}
(4) 计算Ritz值，锁定和清除
(5) 在感兴趣的区域内确定新的位移 \sigma_{i+1}
(6) 分解 I_{j+1, j}+\left(\sigma_{i+1}-\sigma_{i}\right) H_{j+1, j}=Q_{j+1, j+1} R_{j+1, j}
(7) H_{j+1, j}:=\left[\begin{array}{cc} P_j^*& 0\\ 0& 1\end{array}\right] Q_{j+1, j+1}^* H_{j+1, j} R_{j, j}^{-1} P_j
(8) V_{j+1}:=V_{j+1} Q_{j+1, j+1}\left[\begin{array}{cc}P_{j}& 0\\ 0& 1\end{array}\right]

注意在步骤(3)中，我们只对位移和反演矩阵进行m-j次操作，前j个基向量是从前一个位移\sigma_{i-1}保存的。步骤(4)中的锁定、清除和分解操作与隐式重启Arnoldi中描述的操作非常相似，如§7.6所述；我们实际上使用了非厄米版本的厚重启[463]。在步骤(5)中，我们选择新的位移\sigma_{i+1}作为复平面上我们感兴趣的值，以研究矩阵束(8.17)的行为。然而，太接近特征值的选择会导致矩阵A-\sigma_{i+1} B接近奇异。此外，新的位移接近已收敛的Ritz值将使步骤(6)中的上三角因子R_{j,j}接近奇异。

当原始束(8.17)中的第二个矩阵B是正定的，我们可以选择使基V_j对B正交。如果A也是厄米的，那么H_{j,j}也是厄米的，我们可以使用三对角矩阵；参见[322]。

这种分数Krylov的表述与早期[378]使用的不同，在那里整个过程中使用相同的基向量集 V_j，并且束(8.17)被转换为两个海森堡矩阵。我们发现，这种新的表述提供了一种更自然的方式来继续使用新的位移，并在我们面临因线性依赖而失去精度的风险时发出信号。

参见报告[379]中的模型降阶的数值示例。