基本算法

给定一个n阶矩阵A和初始的n阶p块向量P_1和Q_1，使得P^{\ast}_1 Q_1 = I，块Lanczos方法生成一系列n\times p的右和左Lanczos块向量\{ Q_i \}和\{ P_i \}。这些向量是Krylov子空间的双正交基 {\mathcal K}^j(A,Q_1) 和 {\mathcal K}^j(A^{\ast},P_1)。 具体来说，经过j步后，Lanczos过程确定一个块三对角矩阵

T_{[j]} = \begin{bmatrix} A_1 & B_2 & & \\ C_2 & A_2 & \ddots \\ & \ddots & \ddots & B_j \\ & & C_j & A_j \end{bmatrix},

以及右和左Lanczos向量矩阵

Q_{[j]} = \left[\, Q_1, Q_2, \ldots , Q_j \,\right] \quad \text{和} \quad P_{[j]} = \left[\, P_1, P_2, \ldots , P_j \,\right]

满足双正交条件

P^{\ast}_{[j]} Q_{[j]} = I. \tag{7.49}

块Lanczos方法使用三项递归

Q_{j+1} C_{j+1} = AQ_j - Q_j A_j - Q_{j-1} B_j, \tag{7.50}

P_{j+1} B^{\ast}_{j+1} = A^{\ast} P_j - P_j A^{\ast}_j - P_{j-1}C^{\ast}_j \tag{7.51}

对于j = 1, 2, \ldots，其中Q_0 = P_0 = 0。在矩阵表示法中，这些量满足控制关系

AQ_{[j]} = Q_{[j]} T_{[j]} + Q_{j+1}C_{j+1}E^{\ast}_j, \tag{7.52}

A^{\ast} P_{[j]} = P_{[j]} T^{\ast}_{[j]} +P_{j+1} B^{\ast}_{j+1}E^{\ast}_j, \tag{7.53}

其中E_j是一个jp阶p矩阵，其底部方块是一个单位矩阵，其他地方为零。

为了使用Lanczos方法近似A的特征值和特征向量，我们解jp \times jp块三对角矩阵T_{[j]}的特征值问题。每个特征三元组( \theta^{(j)}_i, z^{(j)}_i, w^{(j)}_i ) of T_{[j]}，

T_{[j]} z^{(j)}_i = \theta^{(j)}_i z^{(j)}_i\quad \text{和} \quad (w^{(j)}_i)^{\ast} T_{[j]} = \theta^{(j)}_i (w^{(j)}_i)^{\ast}

确定一个Ritz三元组( \theta^{(j)}_i, x^{(j)}_i, y^{(j)}_i )，其中 x^{(j)}_i = Q_{[j]} z^{(j)}_i和 y^{(j)}_i = P_{[j]} w^{(j)}_i。 通常，Ritz三元组首先近似A的外特征三元组。

为了评估近似值， 设r^{(j)}_i和s^{(j)}_i表示相应的右和左残差向量：

r^{(j)}_i = Ax^{(j)}_i - \theta^{(j)}_i x^{(j)}_i, \tag{7.54}

(s^{(j)}_i)^{\ast} = (y^{(j)}_i)^{\ast} A-\theta^{(j)}_i (y^{(j)}_i)^{\ast}. \tag{7.55}

分别代入(7.52)和(7.53)，出现

r^{(j)}_i = Q_{j+1}C_{j+1} (E^{\ast}_j z^{(j)}_i), \tag{7.56}

(s^{(j)}_i)^{\ast} = ((w^{(j)}_i)^{\ast} E_j) B_{j+1} P^{\ast}_{j+1}. \tag{7.57}

一个Ritz值\theta^{(j)}_i被认为已经收敛，如果两个残差范数都足够小。同样，残差也可以用来确定三元组的后向误差界限，如同标准的非厄米Lanczos算法的情况（见§7.8）。

基本非厄米Lanczos方法的一个简单阻塞版本在算法7.14中呈现。

算法 7.14：NHEP的块Lanczos方法

(1)  选择 n \times p 起始块 P_1 和 Q_1，使得 P_1^* Q_1 = I
(2)  R = AQ_1
(3)  S = A^* P_1
(4)  对于 j = 1, 2, ..., 直到收敛
(5)      A_j = P_j R
(6)      R := R - Q_j A_j
(7)      S := S - P_j A
(8)      QR分解 R = Q_{j+1} C_{j+1}
(9)      QR分解 S = P_{j+1} B
(10)     如果 R 和/或 S 秩亏，停止
(11)     计算SVD：P_{j+1}^* Q_{j+1} = U_j \Sigma_j V_j^*
(12)     如果 \Sigma_j 是（近乎）奇异的，停止
(13)     B_{j+1} := B_{j+1} U_j E_j
(14)     C_{j+1} := E_j^{1/2} V^* C_{j+1}
(15)     Q_{j+1} := Q_{j+1} V_j E_j
(16)     P_{j+1} := P_{j+1} U_j E_j^{-1/2}
(17)     计算 T_{[j]} 的特征三元组 \left(\theta_i^{(j)}, z_i^{(j)}, w_i^{(j)}\right)
(18)     测试收敛性
(19)     如有必要，重新双正交化
(20)     R = AQ_{j+1}
(21)     S = A^* P_{j+1}
(22)     R := R - Q_j B_{j+1}
(23)     S := S - P_j C
(24) 结束循环
(25) 计算近似特征向量 x_i^{(j)} 和 y_i^{(j)}。

我们现在对算法7.14的某些步骤进行评论。

(1)

起始向量Q_1和P_1最好由用户选择，以体现任何关于A的所需特征向量的可用知识。在没有此类知识的情况下，选择Q_1为一个实正交化的随机n\times p矩阵，并令P_1 = Q_1。

块大小p应选择大于或等于任何所需特征值的多重性。如果多重性未知，典型的p值为2，3和6。

(2), (3), (20), (21)

用户需要提供子程序来计算任意n\times p矩阵Z的AZ和A^{\ast}Z乘积。这为用户提供了一次遍历定义A和A^{\ast}的数据结构的机会，可能提高了效率。

(8), (9), (10)

首先，计算R^{\ast}和S的QR分解。如果R^{\ast}和S都是满秩的，那么Q_{j+1}和P_{j+1}是n阶p矩阵，使得 Q^{\ast}_{j+1} Q_{j+1} = P^{\ast}_{j+1} P_{j+1} = I，并且两者 B^{\ast}_{j+1}和C_{j+1}都是p阶p右上三角矩阵。如果R或S不是满秩的，那么揭示了一个不变子空间。在秩亏情况下继续递归在§7.9.2中讨论。

(12)

如果\Sigma_j是奇异的或接近奇异的，即如果有一个零或非常小的奇异值，那么有一个中断。必须采取适当的行动继续。见§7.9.2可能的处理方法。

(17), (18)

可以使用§7.3中讨论的方法计算T_{[j]}的特征三元组。Ritz值的收敛可以通过残差的范数测试（7.56）和（7.57）。关于访问近似精度的更多细节，见§7.8.2和§7.13 和[29]。

当j变大时，解T_{[j]}的特征问题变得耗时。一个简单的降低成本的方法是定期这样做，比如每五步。

(19)

与未阻塞的Lanczos方法一样，在有限精度算术的情况下，块Lanczos向量\{ Q_i \}和\{ P_i \}的双正交性恶化。TSMGS过程可以用来对Q_{j+1}和P_{j+1}进行原地双正交化 Q_{[j]} =\left[\, Q_1, Q_2, \ldots, Q_j\, \right]和 P_{[j]} =\left[\, P_1, P_2, \ldots, P_j\, \right]。

PROCEDURE TSMGS
for i=1,2,\ldots,j
Q_{j+1} := Q_{j+1} - Q_i ( P^{\ast}_i Q_{j+1})
P_{j+1} := P_{j+1} - P_i ( Q^{\ast}_i P_{j+1})
end

(25)

这仅在Lanczos向量的双正交性保持良好时推荐。近似特征向量x^{(j)}_i 和y^{(j)}_i对应于收敛的Ritz 值\theta^{(j)}_i被计算，并且残差（7.54）和（7.55） 相对于\Vert x^{(j)}_i\Vert _2和\Vert y^{(j)}_i\Vert _2再次检查。注意，范数可能大于在步骤（17）处计算的估计值。这是Lanczos基向量病态的一个副作用。