预处理同时迭代法

利用同时或块迭代是一个众所周知的方法，它比单向量方法有重要改进，并允许我们计算一个(m > 1)维的不变子空间，而不是一次计算一个特征向量。它还可以作为并行计算机上单向量方法的加速技术，因为对于极值特征值的收敛通常随着块的大小增加而增加，并且每一步都可以自然地在多种多处理器计算机上实现。

块算法是单向量方法的直接推广，通常与Rayleigh-Ritz过程结合使用。对于前面讨论的每个预处理特征求解器，可以很容易地提出几种同时迭代法的变体；例如，参见带有移位的预处理幂法的块版本[266,149,63,130]。为了简洁起见，我们这里只描述一种方法——局部最优PCG方法的同时版本，建议在[268]中——并且仅针对铅笔B - \mu A。

算法 11.11：用于广义特征值问题的同步PCG方法

(1) 选择一组初始向量 x_{j}^{(0)}, j=1,\ldots, m
(2) 对于 i=0,\ldots，直到收敛执行：
(3) 对于 j=1,\ldots, m，计算 \mu_{j}^{(i)}:=\left(x_{j}^{(i)}, B x_{j}^{(i)}\right)/\left(x_{j}^{(i)}, A x_{j}^{(i)}\right)
(4) 对于 j=1,\ldots, m，计算 r_{j}:=B x_{j}^{(i)}-\mu_{j}^{(i)} A x_{j}^{(i)}
(5) 对于 j=1,\ldots, m，计算 w_{j}:=\operatorname{Tr}_{j}
(6) 对 B-\mu A 使用Rayleigh-Ritz方法在
\operatorname{span}\left\{w_{1},\ldots, w_{m}, x_{1}^{(i)},\ldots, x_{m}^{(i)}, x_{1}^{(i-1)},\ldots, x_{m}^{(i-1)}\right\},
(7) 对于 j=1,\ldots, m，计算 x_{j}^{(i+1)}:= 对应于从顶部数第j个Ritz值的第j个Ritz向量

与单向量算法一样，需要采取特殊措施来克服算法11.11中试验子空间中几乎线性相关基的问题。

目前还没有理论可以准确预测算法11.11的收敛速度。然而，通过类比已知的PCG系统求解器的收敛理论，我们预计\mu^{(i)}_j向\mu_j的收敛是线性的，其比率为

\left( \frac{1 - \sqrt{\xi_j} }{1+\sqrt{\xi_j}} \right)^2,\quad \xi_j = \frac {\delta_0}{\delta_1}\frac{\mu_j - \mu_{m+1}}{\mu_j - \mu_{\min}}. \tag{11.18}

我们通过数值比较（参见[268]）了块变体的最速上升算法11.6和算法11.11，块大小m=3。然而，我们只绘制了两个顶部特征值的误差，省略了第三个特征值。两种红色阴影代表块SA方法，两种蓝色阴影对应于算法11.11。在黑白打印上也很容易区分这些方法，因为算法11.11总是收敛得更快。两条直线对应于(11.18)预测的线性收敛。

在所有测试中，B=I，A是一个对角矩阵，最小条目为1, 2, 3，最大条目为10^{10}，我们测量特征值误差为

\lambda_i^{(i)} - \lambda_j, \ j=1,2.

问题的大小在100到2000之间变化。初始猜测对每次新的运行都是随机的。预处理器T是一个随机的对称正定矩阵，固定值为\kappa (TA) = \delta_1 / \delta_0（参见11.9）。

我们的随机初始猜测导致初始误差非常大，因为矩阵A条件很差，\kappa_2(A)=10^{10}。我们在以下所有图中观察到许多初始误差在10^{10}水平，但两种测试方法在几次迭代后成功地将误差降低到个位数。

我们可以看到，A的巨大条件数、问题的大小、特征值在不需要部分的分布以及预处理器T的特定选择都不会显著影响算法11.11的收敛，因为在所有图中，收敛历史线的元素彼此非常接近。此外，我们的理论预测（11.18）算法11.11的收敛速率相对准确，尽管有些悲观。\delta_1 / \delta_0的10倍增加导致迭代次数增加——块SA方法为10倍，而算法11.11仅为约3倍——正如我们所预期的那样。我们观察到，深色中的第一个特征值的收敛通常比浅色和虚线的第二个特征值快。最后，我们还观察到当迭代次数接近解决问题大小的一半时，超线性加速。然而，迭代次数和问题大小之间的这种大比率在实践中是不典型的，因此我们通过增加问题的大小来避免它；例如，在图11.3中，前两次运行是在大小为200的问题上进行的；注意到超线性收敛后，我们立即将大小增加到1000。

图11.2：共轭梯度与最速上升，\delta _1 / \delta _0=10

图11.3：共轭梯度与最速上升，\delta _1 / \delta _0=100

图11.4：共轭梯度与最速上升，\delta _1 / \delta _0=1000

图11.2到11.4清楚地说明了在§11.1中提出的两个陈述。首先，预处理求解器的性能在很大程度上取决于所使用的预处理器的质量。一个好的预处理器导致图11.2上的快速收敛。一个糟糕的预处理器显著减慢了图11.4上的收敛速度。我们注意到问题的大小并不真正影响收敛速度。其次，即使使用相同的预处理器，每种方法的具体实现也会产生很大的差异。这在质量较差的预处理器上尤为明显，例如在图11.4上。算法11.11，局部最优PCG方法，比算法11.6收敛大约快一百倍，尽管这两种算法使用相同的预处理器并且每次迭代的成本相似。

总之，我们的数值结果表明，算法11.11是一种真正的共轭梯度方法。我们最近的数值结果[269]表明，算法11.11实际上是对称正定特征值问题预处理求解器整个类别的最优方法。

已知[91]有块预处理方法用于同时计算谱的两端，但它们并没有比单独计算最小和最大特征值的标准方法提供太多改进。

另一个已知的构建块特征求解器的想法是使用非线性优化方法来最小化或最大化Rayleigh-Ritz方法中的投影矩阵的迹（参见，例如，[91,16,366]和§9.4]，这里T=I。如果使用共轭梯度方法进行优化，这会导致与算法11.11类似的方法，但算法11.11有一个使用Rayleigh-Ritz方法的优势。

使用锁定，一种降阶形式，利用不同特征向量的不均匀收敛率，提高了上述预处理同时迭代方法的性能，与没有预处理的经典方法类似；参见§4.3。由于同时迭代计算的每个近似特征向量的不同收敛率，通常的做法是提取那些已经收敛的特征向量并进行降阶。我们冻结提取的近似值并将其从后续迭代中移除，这样就不需要继续将它们与任何矩阵相乘。然而，我们仍然需要将它们包含在Rayleigh-Ritz方法的试验子空间的基中，或者在需要时对冻结向量进行后续正交化。

前一种可能性似乎更自然，并且比后者更容易编程。利用已知的关于Rayleigh-Ritz方法精度的结果，分析其对迭代方法的影响也相当简单；参见，例如，[387,267]。然而，正交化的成本可能更低，因为它减少了试验子空间的维数。

我们注意到，当尝试一次计算一组特征向量时，也应该在单向量迭代方法中使用锁定。

不幸的是，如果我们想要能够在进一步迭代过程中研究由此产生的误差的传播，那么在预处理特征求解器中使用正交化进行锁定可能并不简单。

作为一个例子，让我们考虑我们最简单的预处理特征求解器，算法11.5，带有由正交投影器P^\perp定义的附加正交化，该投影器位于已经计算的特征向量所张成的子空间的正交补上，可以写成

x^{(i+1)} = P^\perp w^{(i)} + \tau^{(i)} x^{(i)}, \qquad w^{(i)} = T (B x^{(i)} - \mu^{(i)} A x^{(i)}). \tag{11.19}

我们需要弄清楚在我们的正交投影器P^\perp中要使用哪些标量积。

首先，我们需要为已经计算的特征向量所张成的子空间的正交补选择一个标量积。基于A的标量积在这里似乎是一个自然的选择。当B是正定矩阵时，通常也使用基于B的标量积。

其次，我们需要定义一个标量积，其中投影器P^\perp是正交的。传统的方法是使用与第一步相同的标量积。在代码中实现这一点也很简单。

不幸的是，通过这样的选择，方法（11.19）中的迭代算子失去了相对于基于T^{-1}的标量积的对称性。这使得理论研究正交化对近似计算的特征向量的影响变得相当复杂（参见[147,149]），其中进行了直接的扰动分析。

为了保持对称性，我们必须使用基于T^{-1}的正交投影器P^\perp，尽管我们在第一步中使用了不同的标量积，例如基于A的标量积来定义正交补。通过这种选择，我们可以使用标准的和简单的后向误差分析[264,265]，而不是直接分析[147,149]，但给定w的P^\perp w的实际计算需要特别注意。

遵循[264,265]，我们取由冻结的近似特征向量张成的原始子空间（我们称之为\tilde X），并找到子空间TA \tilde X的基于T^{-1}的正交基。从数学上讲，后者的基于T^{-1}的正交补与原始子空间\tilde X的基于A的正交补一致。现在，我们可以使用一个标准的基于T^{-1}的投影器到TA \tilde X的基于T^{-1}的正交补上。显然，必须能够计算基于T^{-1}的标量积才能使用这个技巧。

我们还注意到，使用基于T^{-1}的标量积，以及任何基于病态矩阵的标量积，可能会导致不稳定的计算。