GEMV形式的Arnoldi过程

首先，我们将从稍有不同的视角来审视阿诺尔迪过程。阿诺尔迪过程在第§7.5节中已有推导，并在算法7.3中以修正格拉姆-施密特算法的形态呈现。本节中，我们将采用经典的格拉姆-施密特过程结合正交化细化[96]，以矩阵-向量（或GEMV）形式描述该过程。阿诺尔迪关系涉及矩阵A、基矩阵V_k及残差向量f_k，其形式为：

A V_k = V_k H_k + f_k e_k^{\ast},

其中，V_k \in {\mathcal C}^{n \times k}的列向量正交归一，V_k^{\ast} f_k = 0，且H_k \in {\mathcal C}^{k\times k}为上部海森堡矩阵，其次对角线元素非负。我们称此为A的k步阿诺尔迪分解。从构造过程可轻易看出，H_k = V_k^{\ast} A V_k为上部海森堡矩阵。当A为厄米矩阵时，这意味着H_k为实对称且三对角矩阵。V_k的列向量被称为阿诺尔迪向量。

算法7.6给出了计算k步阿诺尔迪分解的模板。

算法 7.6：k步Arnoldi分解
(1)  对起始向量 v，计算 v_1 = \frac{v}{\|v\|}
(2)  w = Av_1
(3)  h = v_1^* w; f_1 = w - v_1 h
(4)  V_1 = [v_1]; H_1 = [h]
(5)  对于 j = 1, 2, 3, \ldots, k-1
(6)      \beta_j = \|f_j\|
(7)      v_{j+1} = {f_j}/{\beta_j}
(8)      V_{j+1} := [V_j, v_{j+1}]; \hat{H}_j := \left[\begin{array}{c} H_j\\ \beta_j e_j^* \end{array}\right]
(9)      w = Av_{j+1}
(10)     h = V_{j+1}^* w
(11)     f_{j+1} = w - V_{j+1} h
(12)     如果 \|f_{j+1}\| \lt \eta \|h\|
(13)         s = V_{j+1}^* f_{j+1}
(14)         f_{j+1} := f_{j+1} - V_{j+1} s
(15)         h := h + s
(16)     结束如果
(17)     H_{j+1} := [\hat{H}_j, h]
(18) 结束循环

接下来，我们将详细描述一些实现细节，参照算法7.6中的各个阶段。

(1)

选择初始起始向量v并归一化。理想情况下，对于特征值计算，应尝试构建一个在感兴趣的特征向量方向上占主导地位的v。在没有其他考虑的情况下，随机向量是一个合理的选择。

(2)

初始化阿诺尔迪过程。在后续步骤中，f_1将成为新的阿诺尔迪基向量。此步骤应包含与其他步骤相同的再正交化过程。

(7)

归一化f_j以获得新的基向量v_{j+1}，部分更新H_j以得到(j+1) \times j矩阵\hat{H}_{j}，并计算新信息w \leftarrow Av_{j+1}。范数计算\beta_j = \Vert f_j \Vert在此处无需重新计算，已在步骤(12)中计算并可在此处重用。

(10)

这是经典的格拉姆-施密特步骤，用于将w相对于V_{j+1}的列向量进行正交化。此形式允许使用二级矩阵-向量运算。详细的讨论将在下文给出。

(12)

此if-语句中的序列确保新方向f_{j+1}在数值上与先前计算的方向（即V_{j+1}的列向量）正交。参数\eta必须指定（0 \le \eta < \infty）。此测试询问“w = Av_{j+1}是否接近已构建的空间？”比值\frac{\Vert f_{j+1}\Vert}{ \Vert h \Vert} =\frac{\sin(\theta)}{\cos(\theta)}，其中\theta是向量w与\mathrm{span}(V_{j+1})之间的角度（即w与其投影V_{j+1} h之间的角度）。\eta值越大，正交性要求越严格。当\eta = 0时，算法退化为经典格拉姆-施密特（CGS）。还需进一步细节以完全指定此过程。如果在更新f_{j+1}和h后，关系\Vert f_{j+1} \Vert < \eta \Vert h \Vert仍然成立，则在数值上我们有w \in \mathrm{span}(V_{j+1})。此时，阿诺尔迪过程应终止或生成一个随机向量与V_{j+1}的列向量正交以替换f_{j+1}，并将\beta_{j+1}设为零以继续分解。ARPACK实现中使用的\eta值为1/\sqrt{2}（参见第§7.6.9节）。关于此正交化装置的进一步讨论将在下文给出。

(17)

H_{j+1}现在是一个(j+1) \times (j+1)的上部海森堡矩阵。

密集矩阵-向量乘积V_{j+1}^{\ast} w和V_{j+1} h可以使用二级BLAS操作GEMV表示。这为从工作站到超级计算机的几乎所有平台提供了显著的性能优势。此外，ScaLAPACK项目[52]及多家高性能计算机供应商已为多种并行机器优化了这些内核。

用于将新信息Av_k相对于现有基V_k进行正交化的机制是CGS。它以不稳定著称，在此设置中若不加修改将彻底失败。一种补救措施是使用算法7.3中采用的修正格拉姆-施密特过程。遗憾的是，在我们即将讨论的重启情况下，这也无法产生正交向量，且在此设置中无法用二级BLAS表示。在步骤(5)中，我们引入了Daniel、Gragg、Kaufman和Stewart（DGKS）在[96]中提出的迭代细化技术。该方案提供了一种极佳的方法来构建在数值上与V_{j+1}正交的向量f_{j+1}。步骤(5)中指定的简单测试用于避免在不需要时进行此DGKS修正。

此机制以非常合理的成本维持了全工作精度的正交性。我们即将讨论的重启方案所施加的特殊情况使得这一修改对于获得准确的特征值和数值上正交的舒尔向量（厄米情况下的特征向量）至关重要。基于MGS的方案通常足以求解线性系统，因为它们确实构建了一个线性独立的基集，尽管可能不是数值上正交的。MGS正交性的质量依赖于原始向量集的条件数（线性独立性）。我们即将介绍的隐式重启机制在未保持数值正交性时效果较差，甚至可能失败。

已有充分文献证明，未保持正交性会导致兰索斯或阿诺尔迪过程中出现多种数值困难。在厄米情况下，Paige[347]表明，正交性丧失恰好在H_j的特征值接近 A 的特征值时发生。实际上，阿诺尔迪向量在相关近似特征向量的方向上失去正交性。此外，未能保持正交性会导致阿诺尔迪方法产生近似特征值的虚假副本。在厄米情况下，已有多种无需完全再正交化的补救措施被提出。这些在第§4.4节（第页）中进行了讨论。