带位移反转的Lanczos算法

在许多实际应用中，位移反转变体是一种自然选择。 它对应于以下选择：

C=B(A-\sigma B)^{-1}. \tag{5.16}

这给出了接近所选位移点\sigma的特征值，并且我们可以在更少的步骤j后获得收敛。即使带有位移矩阵(A-\sigma B)的系统比直接变体中所需的B更难以解决，所需的较少步骤j通常会弥补这一点。

基本递归公式（5.9）被替换为

B(A-\sigma B)^{-1}V_j=V_{j}T_{j}+re_j^{\ast}, \tag{5.17}

现在使用B^{-1}-正交（5.10）基V_j和B正交辅助基W_j，其中V_j=BW_j。这两个基是双正交的， V_j^{\ast}W_j=I_{j}，与非对称双侧Lanczos算法（如第7.8节所述）相比。

矩阵T是对称的，

T_{j}=V_j^{\ast}(A-\sigma B)^{-1}V_j. \tag{5.18}

T_{j}的特征解，

T_{j}s_i^{(j)}=s_i^{(j)}\theta_i^{(j)},

用于获得近似特征值

\lambda_i^{(j)}=\sigma+\frac{1}{\theta_i^{(j)}} \tag{5.19}

和近似特征向量

x_i^{(j)}=W_js_i^{(j)}. \tag{5.20}

其残差为

\begin{aligned} r_i^{(j)} &= Ax_i^{(j)}-Bx_i^{(j)}\lambda_i^{(j)}=(A-\sigma B)x_i^{(j)}-\frac{1}{\theta^{(j)}_i}Bx_i^{(j)} \\ &= \frac{1}{\theta^{(j)}_i}\left((A-\sigma B)W_js_i^{(j)}\theta^{(j)}_i-V_j s^{(j)}\right) \\ &= -\frac{1}{\theta^{(j)}_i}(A-\sigma B)w_{j+1}\beta_j s^{(j)}_{j,i}, \end{aligned}

其范数在以下量小时也小

\beta_js_{j,i}^{(j)}/\theta_i^{(j)} \tag{5.21}

这个近似不是标准的，而是原特征问题的调和Ritz值（参见第3.2节的讨论）。

与标准情况一样，我们可以监控 \beta_js_{j,i}^{(j)}/\theta_i^{(j)} 并在其变小时标记特征值\lambda_i为收敛，而无需实际执行矩阵向量乘法（5.20）。我们保存这一操作，直到步骤j时估计指示收敛。

我们得到以下算法。

算法 5.5 GHEP的位移-逆 Lanczos 方法
(1)  从初始向量 r = x 开始。计算 q = B r 和 \beta_0 = \left\| q^* r \right\|^{1/2}
(2)  对于 j = 1, 2, \ldots，直到收敛，
(3)      w_j = r / \beta_{j-1}, v_j = q / \beta_{j-1}
(4)      计算 r = (A - \sigma B)^{-1} v_j
(5)      r = r - w_{j-1} B_{j-1}
(6)      a_j = v_j^* r
(7)      r = r - w_j a_j
(8)      必要时进行再正交化
(9)      计算 q = B r
(10)     B_j = \left\| q^* r \right\|^{1/2}
(11)     计算近似特征值 T_j = S \Theta^{(j)} S^*
(12)     测试收敛性
(13) 结束循环
(14) 计算近似特征向量 X = W_j S

我们注意到，只有少数步骤与之前的直接迭代算法不同。在实际实现中，可以使用相同的程序进行这两种计算。让我们评论那些不同的步骤：

(8)

现在我们期望对于相当少的向量j收敛，建议使用完全再正交化，即使选择性再正交化给出的结果足够准确。应用Gram-Schmidt正交化过程，

r=r-W_j(W_j^{\ast}(Br)),

直到r对基W_j是B-正交的。我们需要每次再正交化额外进行一次B的矩阵向量乘法，并且必须使用经典的Gram-Schmidt正交化过程。

我们可以保存并使用辅助基V_j以节省额外的B操作，

r=r-W_j(V_j^{\ast}r).

(14)

现在使用基W_j来获取特征向量（5.20），

x_i^{(j)}=W_js_i^{(j)},

对于每个标记为收敛的i。

当我们运行位移反转算子（5.16）时，我们因子化

LDL^{\ast}=P^T(A-\sigma B)P \tag{5.22}

对于适当的稀疏保持排列P，在开始时使用稀疏高斯消去法。

在实际迭代过程中，我们使用因子L和U，计算

r=P(L^{-\ast}(D^{-1}(L^{-1}(P^Tv_j)))) \tag{5.23}

在算法5.5的步骤(4)中。

如果特征值\lambda在位移\sigma的两侧，矩阵A-\sigma B是不定的，我们不能使用简单的Cholesky分解，而必须进行对称不定分解，例如Duff和Reid的MA47（参见第10.3节）。这种类型的算法使D块对角化，包含1x1和2x2块。很容易确定这种矩阵的惯性，这给出了关于位移\sigma两侧有多少特征值的信息。我们可以通过在端点进行两次因子化来确定区间内的特征值数量。这在基于[162]和[206]的软件中用于确保计算了区间内的所有特征值。

我们注意到，算法5.5不需要对矩阵B进行任何因子化，并且可以在B奇异时应用。这在实际情况下并不罕见，但必须采取特殊预防措施，以防止B的零空间分量进入计算；参见[161]或[340]。