Q^{\ast } H Q的稳定性

如算法 7.8 所示的收缩正交变换构建过程显然是稳定的（即，分量相对精确地表示了在精确算术中可获得的变化）。毫无疑问，相似变换 Q^* H Q 在数值上保留了 H 的特征值。然而，关于这些变换在清除过程中数值上保持海森堡形式的效果存在严重疑问。需要对基本算法进行修改，以确保如果 y^*H = \theta y^*，那么 H_+ \equiv Q^* H Q 在数值上是海森堡的（即，次对角线以下的所有元素相对于 \Vert H \Vert 都是微小的）。

H_+ = Q^* H Q 为海森堡形式取决于表达式中的项 g y^* H R 消失，即

Q^* H Q = L^*HR + g y^*H R + \theta e_1 e_1^*.

然而，进一步检查发现，

e_1^* Q = e_1^* L + (e_1^* y) g^* = \eta_1 g^*,

其中 \eta_1 是 y 的第一个分量。因此，

\Vert g\Vert = \frac{1}{\vert \eta_1 \vert}, \tag{7.22}

因此当 y 的第一个分量很小时，可能会出现数值困难。具体来说，y^*H = \theta y^* 并且因此在精确算术中 y^* H R = 0。然而，在有限精度下，计算得到的 fl(y^* H) = \theta y^* + z^*。误差 z 相对于 \Vert H \Vert 约为 \epsilon_M，但

\Vert g y^*H R \Vert = \Vert g\Vert \cdot \Vert z^* R \Vert =\frac{1}{\vert \eta_1 \vert} \Vert z^* R \Vert,

因此这一项可能相当大。如果 \eta_1 = O(\epsilon_M)，它可能达到 O(1) 量级。这是一个严重的问题，并且在实践中如果不进行我们提出的修改，这种情况就会发生。

解决办法是引入逐步可接受的扰动和向量 y 的重缩放，同时强制满足以下条件：

Q^* y = e_1 \ \ \mathrm{和}\ \ y^* H Q = \theta e_1^*

在有限精度下具有足够的准确性。为了实现这一点，我们将设计一个方案，以数值方式实现

y^* H q_j = 0 \ \ \mathrm{对于} \ \ j > 1

如文献 [420] 所示，这一修改足以在数值上建立 q_i^* H q_j = 0，相对于 \Vert H \Vert，对于 i > j + 1。

基本思想如下：如果在第 j 步，计算得到的 fl(y^* H q_j ) 不够小，则通过缩放向量 y_j 和分量 \eta_{j+1} 来调整使其足够小。在计算 q_{j+1} 之前进行这种重缩放，我们有 y_j \leftarrow y_j \phi 和 \hat{y}_j \leftarrow \hat{y}_j \psi，其中 y^* = [ y_j^* , \hat{y}_j^*]。当然，\Vert y \Vert 不应因这种缩放而改变，因此这是必要的。这给出了以下方程组来确定 \phi 和 \psi：如果 \vert y_j^* H_j \hat{q}_j + \rho_j \beta_j \eta_{j+1} \vert >\epsilon_M \tau_{j+1}，

\begin{aligned} (\tau_j \phi)^2 + (1 - \tau_j^2)\psi^2 &= 1, \\ y_j^* H_j \hat{q}_j\phi + \rho_j \beta_j \eta_{j+1}\psi &= \pm\epsilon_M \tau_{j+1}. \end{aligned}

如果 \rho_j 在 \sqrt{\epsilon_M} 量级，那么这种缩放可以被吸收到 \rho_j 中，而不改变 y 并且不影响 Q 的数值正交性。当 y 被修改时，之前计算的 q_i（2 \le i < j）不需要改变。在第 j 步之后，向量 y_j 仅在后续步骤中进行重缩放，并且定义 q_i（2 \le i \le j）的公式对于 y_j 的缩放是不变的。有关详细信息，请参阅文献 [420]。

算法 7.9 所示的代码实现了这一方案，以计算可接受的 Q。^[1] 在实践中，这种变换可以在不存储 Q 的情况下就地计算并应用于 H := Q^* H Q 和 V := VQ。然而，这种实现非常微妙，构造 Q 的细节要求避免了额外的存储。该代码与上述描述略有不同，因为向量 y 仅在确定所有列 2 到 n 之后才被重缩放以具有单位范数。

算法 7.9：IRAM的稳定正交消去变换
(1)  从满足 \|y\|=1 的n维向量y开始，\eta=e_n^* y  以及满足 T y=\theta y 的n阶对称三对角矩阵T
(2)  Q=0
(3)  当 y(i)=0, y(i)=\epsilon_M / n, i:=i+1，结束循环
(4)  \sigma=y(1)^2, \quad \tau_0=|y(1)|
(5)  对于 j=2,\ldots, n
(6)      \sigma=\sigma+y(j)^2, \quad \tau=\sqrt{\sigma}
(7)      \gamma=-(y(j)/\tau)/\tau_0
(8)      \rho=\tau_0/\tau
(9)      Q(1: j-1, j)=y(1: j-1)\gamma, Q(j, j)=\rho
(10)     如果 (j\lt n 且 \tau\lt.05)
(11)         \tau_p=\sqrt{\sigma+y(j+1)^2}
(12)         \alpha=y(1: j)^* H(1: j, 1: j) Q(1: j, j)
(13)         \delta=y(j+1) H(j, j+1)\rho
(14)         \sigma_0=-1
(15)         如果 \operatorname{sign}(\alpha)\cdot\operatorname{sign}(\delta)\lt 0, \sigma_0=1，结束如果
(16)         如果 |\delta+\alpha|>\epsilon_M\tau_p
(17)             如果 \rho<\sqrt{\epsilon_M} / 100
(18)                 Q(j, j)=Q(j, j)\sigma_0(\epsilon_M\tau_p+|\alpha|)/|\delta|
(19)             否则
(20)                 如果 |\alpha|>|\delta|
(21)                     \phi=\sigma_0(\epsilon_M\tau_p+|\delta|)/|\alpha|
(22)                     y(1: j)=y(1: j)\phi
(23)                     \sigma=\sigma\phi^2, \quad \tau=\tau|\phi|
(24)                 否则
(25)                     \psi=\sigma_0(\epsilon_M\tau_p+|\alpha|)/|\delta|
(26)                     y(j+1: n)=y(j+1: n)\psi
(27)                 结束如果
(28)             结束如果
(29)         结束如果
(30)     结束如果
(31)     \tau_0=\tau
(32) 结束循环
(33) Q(:, 1)=y/\|y\|

有几个实现问题需要注意：

(3)

这里显示的扰动避免了特征向量 y 中初始零项的问题。理论上，当 H 未简化时，这种情况不应该发生，但在数值上，当 H 的对角线元素很小但非零时，这种情况可能发生。有一种更简洁的实现方法可以不修改零项。这是最简单（也是最粗糙）的修正。

(10)

一旦 \tau_j = \Vert y_j\Vert 足够大，就不需要进一步的修正。删除这个 if-clause 将恢复到算法 7.8 所示的未修改的 Q 计算。

(16)

这显示了几种可能的修改 y 的方法之一，以实现 Q^* H Q 的次对角线以下元素数值上微小的目标。更复杂的策略会将尽可能多的缩放吸收到对角线元素 Q(j,j) = \rho 中。这里，不缩放 \rho 而是缩放 y 的分支被设计成使得 \mathrm{fl}(1 + (\rho \psi)^2) = \mathrm{fl}(1 + \rho^2) = 1。

1. 现在有一种更好的算法用于构建和应用这些变换；参见D. Sorensen的《隐式重启Arnoldi方法的去膨胀化》，CAAM技术报告，TR 98-12，莱斯大学，1998年（2000年8月修订）。