收缩

我们现在考虑引入收缩过程的以下实现方式。迄今为止，我们已经描述了仅计算一个特征对（eigenpair）的算法。若需计算多个特征对，则有两种可能的选择。

第一种选择是在重启时将v_1取为近似特征向量的线性组合。例如，如果我们需要计算p个最右侧的特征向量，我们可以取

\hat v_1 = \sum_{i=1}^p \rho_i \tilde u_i,

其中特征值按其实部的递减顺序编号。然后通过归一化\hat v_1得到向量v_1。系数\rho_i的最简单选择是取\rho_i = 1, i=1,\ldots,p。这种方法有几个缺点，其中最重要的是无法以系统的方式选择系数\rho_i。结果是对于难题，收敛性难以实现。

一个更可靠的替代方法是逐个计算特征对并使用收缩法。矩阵A可以通过逐步构建前k个舒尔向量来显式收缩。如果已经计算出不变子空间的前k-1个正交基U_{k-1} = [u_1,\ldots,u_{k-1}]，那么在计算特征值\lambda_{k}时，我们可以使用矩阵

\tilde A = A - U_{k-1} \Sigma U_{k-1}^{\ast},

其中\Sigma = \mathrm{diag}(\sigma_i)是一个位移的对角矩阵。\tilde A的特征值由两组组成。与舒尔向量u_1, \ldots, u_{k-1}相关的特征值将位移至\tilde \lambda_i = \lambda_i - \sigma_i，而其余特征值保持不变。如果寻求具有最大实部的特征值，则选择位移使得\lambda_k成为\tilde A中具有最大实部的下一个特征值。也可以通过简单地投影出与由U_{k-1}张成的不变子空间相关的分量来进行收缩；这将导致操作矩阵

\tilde A = A (I - U_{k-1} U_{k-1}^{\ast}).

注意，如果A U_{k-1} = U_{k-1} R_{k-1}是与前k-1个里兹值相关的部分舒尔分解，那么\tilde A = A - U_{k-1} R_{k-1} U_{k-1}^{\ast}。与舒尔向量u_1, \ldots, u_{k-1}相关的特征值现在都将移动到零。

一个更好的收缩实现方式，与Arnoldi过程很好地结合，是使用单个基v_1, v_2, \ldots , v_m，其前几个向量是已经收敛的舒尔向量。假设已有k-1个这样的向量收敛，记为v_1, v_2,\ldots,v_{k-1}。然后我们选择一个与v_1,\ldots,v_{k-1}正交且范数为1的向量v_k。接下来执行m-k步Arnoldi过程，其中向量v_j相对于所有之前的v_i，包括v_1,\ldots,v_{k-1}保持正交。这将生成子空间的正交基

{\rm span}\{ v_1,\ldots, v_{k-1} , v_k, A v_k, \ldots, A^{m-k} v_k \}. \tag{7.9}

因此，这个修改的Krylov子空间的维数通常是恒定的，等于m。结合Arnoldi方法的隐式收缩过程的概要如下。

算法 7.5：显式重启带消去的Arnoldi方法用于NHEP
(1)  设置 k := 1
(2)  循环
(3)      对于 j = k, k+1, \ldots, m 执行
(4)      w := Av_j
(5)      计算一组 j 个系数 h_{ij}，使得 w := w - \sum_{i=1}^{j} h_{ij} v_i 正交于所有先前的 v_i，i = 1, 2, \ldots, j。
(6)      h_{j+1, j} = \|w\|_2
(7)      v_{j+1} = \frac{w}{h_{j+1, j}}
(8)      计算与特征值 \tilde{\lambda}_k 相关的 A 的近似特征向量及其相关的残差范数估计 \rho_k。
(9)      将此特征向量正交化以得到近似的Schur向量 \widetilde{u}_k 并定义 v_k := \widetilde{u}_k。
(10)     如果 \rho_k 足够小，则（接受特征值）
(11)         h_{i,k} = v_i^* Av_k，i = 1, \ldots, k
(12)         设置 k := k + 1，
(13)         如果 k \geq \text{nev}，则停止，否则跳转到 (2)
(14)     否则
(15)         跳转到 (3)
(16)     结束如果

注意，在循环中，与特征值\lambda_1,\ldots,\lambda_{k-1}相关的舒尔向量在后续步骤中不会被触及。它们有时被称为“锁定向量”。同样，与这些向量对应的上三角矩阵也被锁定。

\underbrace{\left[v_1, v_2, \ldots, v_{k-1}\right.}_{Locked},\underbrace{\left. v_k, v_{k+1}, \ldots v_m \right] }_{Active}

当一个新的舒尔向量收敛时，步骤(10)计算与这个新基向量相关的R的第k列。在后续步骤中，近似特征值是算法中定义的m \times m海森堡矩阵H_m的特征值，其k \times k主子矩阵是上三角的。例如，当m=6且第二个舒尔向量k=2收敛后，矩阵H_m将具有以下形式

H_m = \begin{bmatrix} * & * & * & * & * & * \\ & * & * & * & * & * \\ & & * & * & * & * \\ & & * & * & * & * \\ & & & * & * & * \\ & & & & * & * \end{bmatrix} \tag{7.10}

在后续步骤中，只需考虑与2 \times 2上三角矩阵无关的特征值。

可以证明，在精确算术中，下部(2 \times 2)块中的(n-k) \times (n-k)海森堡矩阵与应用于矩阵

\tilde A = (I-U_{k-1} U_{k-1}^{\ast} ) A.

的Arnoldi运行所得到的矩阵相同。因此，我们隐式地从A的范围中投影出已经计算出的不变子空间。