4.5 带位移的隐式QR迭代法

QR迭代法中需要考虑的另一个重要问题就是运算量：每一步迭代都需要做一次QR分解和矩阵乘积，运算量为 $\mathcal{O}(n^3)$ 。即使每计算一个特征值只需迭代一步，则计算所有特征值也需要 $\mathcal{O}(n^4)$ 的运算量。因此QR迭代算法？？虽然简单易用，但离实际应用还很远。下面我们就想办法将总运算量从 $\mathcal{O}(n^4)$ 减小到 $\mathcal{O}(n^3)$ 。

为了实现这个目标, 我们需要利用 Hessenberg 矩阵. 具体步骤如下: 首先通过相似变化将 $A$ 转化成一个上 Hessenberg 矩阵, 然后再对这个 Hessenberg 矩阵实施隐式 QR 迭代. 所谓隐式 QR 迭代, 就是在 QR 迭代中, 我们不需要进行显式的 QR 分解, 以及 $R$ 和 $Q$ 的乘积. 这样就可以将 QR 迭代的每一步运算量从 $\mathcal{O}(n^{3})$ 降低到 $\mathcal{O}(n^{2})$ . 从而将总的运算量降低到 $\mathcal{O}(n^{3})$ .

带位移的隐式QR迭代法是当前计算中小规模矩阵的所有特征值的标准方法.

4.5.1 上Hessenberg矩阵

设 $H = [h_{ij}] \in \mathbb{R}^{n \times n}$ , 若当 $i > j + 1$ 时, 有 $h_{ij} = 0$ , 则称 $H$ 为上Hessenberg矩阵.

定理4.4设 $A\in \mathbb{R}^{n\times n}$ ，则存在正交矩阵 $Q\in \mathbb{R}^{n\times n}$ ，使得 $QAQ^{\mathrm{T}}$ 是上Hessenberg矩阵

下面我们以一个 $5 \times 5$ 的矩阵 $A$ 为例, 给出具体的上 Hessenberg 化过程, 所采用的工具仍然是 Householder 变换.

第一步：令 $Q_{1} = \mathrm{diag}(I_{1\times 1},H_{1})$ ，其中 $H_{1}$ 是对应于向量 $A(2:5,1)$ 的Householder矩阵.于是可得

由于用 $Q_1^{\mathsf{T}}$ 右乘 $Q_{1}A$ 时，不会改变 $Q_{1}A$ 第一列元素的值，故

第二步：令 $Q_{2} = \mathrm{diag}(I_{2\times 2},H_{2})$ ，其中 $H_{2}$ 是对应于向量 $A_{1}(3:5,2)$ 的Householder矩阵，则用 $Q_{2}$ 左乘 $A_{1}$ 时，不会改变 $A_{1}$ 的第一列元素的值.用 $Q_{2}^{\mathsf{T}}$ 右乘 $Q_{2}A_{1}$ 时，不会改变 $Q_{2}A_{1}$ 前两列元素的值.因此，

第三步：令 $Q_{3} = \mathrm{diag}(I_{3\times 3},H_{3})$ ，其中 $H_{3}$ 是对应于向量 $A_{2}(4:5,3)$ 的Householder矩阵，则有

这时, 我们就将 $A$ 转化成一个上Hessenberg矩阵, 即 $QAQ^{\mathsf{T}} = A_{3}$ 其中 $Q = Q_{3}Q_{2}Q_{1}$ 是正交矩阵, $A_{3}$ 是上Hessenberg矩阵.

下面是将任意一个矩阵转化成上Hessenberg矩阵的算法

算法4.7.上Hessenberg化(UpperHessenbergReduction)

1: Set $Q = I$
2: for $k = 1$ to $n - 2$ do
3: compute Householder matrix $H_{k} = I - \beta_{k}v_{k}v_{k}^{\mathsf{T}}$ with respect to $A(k + 1:n,k)$
4: $w_{k}^{\mathsf{T}} = v_{k}^{\mathsf{T}}A(k + 1:n,k:n)$
5: $A(k + 1:n,k:n) = H_kA(k + 1:n,k:n) = A(k + 1:n,k:n) - \beta_kv_kw^{\mathsf{T}}_k$
6: $w_{k} = A(1:n,k + 1:n)v_{k}$
7: $A(k:n,k + 1:n) = A(k:n,k + 1:n)H_k^\top = A(k:n,k + 1:n) - \beta_kw_kv_k^\top$
8: end for
9: for $k = n - 2$ to 1 do

10: % 用向后累积法计算 $Q$

11: end for

在实际计算时, 我们不需要显式地形成 Householder 矩阵 $H_{k}$ .
可以先存储所有 Householder 向量, 最后用向前累积法计算 $Q$ .

上述算法的运算量大约为 $\frac{14}{3} n^3 + \mathcal{O}(n^2)$ . 如果不需要计算特征向量, 则正交矩阵 $Q$ 也不用计算, 此时运算量大约为 $\frac{10}{3} n^3 + \mathcal{O}(n^2)$ .

上Hessenberg矩阵的一个很重要的性质就是在QR迭代中能保持形状不变.

定理4.5 设 $A \in \mathbb{R}^{n \times n}$ 是非奇异上Hessenberg矩阵, 其QR分解为 $A = QR$ , 则 $\tilde{A} \triangleq RQ$ 也是上Hessenberg矩阵. (板书)

证明. 设 $A = QR$ 是 $A$ 的 QR 分解, 则 $Q = AR^{-1}$ . 由于 $R$ 是一个上三角矩阵, 所以 $R^{-1}$ 也是一个上三角矩阵. 因此 $Q$ 的第 $j$ 列是 $A$ 的前 $j$ 列的线性组合. 又 $A$ 是上 Hessenberg 矩阵, 所以 $Q$ 也是一个上 Hessenberg 矩阵.

相类似地, 我们很容易验证 $RQ$ 也是一个上Hessenberg矩阵. 所以结论成立.

若 $A$ 是奇异的, 也可以通过选取适当的 $Q$ , 使得定理 4.5 中的结论成立. 事实上, 由基于 Gram-Schmidt 过程的 QR 分解可知, 如果 $A$ 是奇异的上 Hessenberg 矩阵, 则 $Q$ 仍可取为上 Hessenberg 矩阵.

由这个性质可知, 如果 $A \in \mathbb{R}^{n \times n}$ 是上Hessenberg矩阵, 则QR迭代中的每一个 $A_{k}$ 都是上Hessenberg矩阵. 这样, 在进行QR分解时, 运算量可大大降低.

Hessenberg 矩阵还有一个重要性质, 就是在 QR 迭代过程中能保持下次对角线元素非零.

定理4.6 设 $A \in \mathbb{R}^{n \times n}$ 是非奇异的上Hessenberg矩阵, 且下次对角线元素均非零, 即 $a_{i+1, i} \neq 0$ , $i = 1, 2, \ldots, n-1$ . 设其QR分解为 $A = QR$ , 则 $\tilde{A} \triangleq RQ$ 的下次对角线元素也都非零.

(留作练习)

易知, 如果上Hessenberg矩阵 $A$ 存在某个下次对角线元素为零, 则 $A$ 一定可约, 此时计算 $A$ 的特征值就转化为计算两个更小规模的矩阵的特征值. 因此, 我们只需考虑下次对角线均非零的情形.
需要指出的是, 如果上 Hessenberg 矩阵 $A$ 存在某个下次对角线元素为零, 则 $A$ 一定可约, 但反之却不成立, 即如果 $A$ 是可约的上 Hessenberg 矩阵, $A$ 的下次对角线元素仍可能均非零.

推论4.7 设 $A \in \mathbb{R}^{n \times n}$ 是非奇异的上Hessenberg矩阵, 且下次对角线元素均非零, 则在带位移的QR迭代中, 所有的 $A_{k}$ 的下次对角线元素均非零.

4.5.2 隐式QR迭代

在QR迭代中,我们需要先做QR分解 $A_{k} = Q_{k}R_{k}$ ，然后再计算 $A_{k + 1} = R_kQ_k$ .但事实上，我们可以将这个过程进行简化，即在不对 $A_{k}$ 进行QR分解的前提下，直接计算出 $A_{k + 1}$ .这就是隐式QR迭代.

我们这里考虑不可约的上Hessenberg矩阵，即 $A$ 的下次对角线元素都不为0.事实上，若 $A$ 是可约的，则 $A$ 就是一个块上三角矩阵，这时 $A$ 的特征值计算问题就转化成计算两个对角块的特征值问题.隐式QR迭代的理论基础就是下面的隐式 $Q$ 定理

定理4.8 (Implicit $Q$ Theorem) 设 $H = Q^{\mathsf{T}}AQ \in \mathbb{R}^{n \times n}$ 是一个不可约上Hessenberg矩阵, 其中 $Q \in \mathbb{R}^{n \times n}$ 是正交矩阵, 则 $Q$ 的第2至第 $n$ 列均由 $Q$ 的第一列所唯一确定 (可相差一个符号).

(板书)

证明. 设 $H = Q^{\mathsf{T}}AQ$ 和 $G = V^{\mathsf{T}}AV$ 都是不可约上Hessenberg矩阵, 其中 $Q = [q_{1}, q_{2}, \ldots, q_{n}], V = [v_{1}, v_{2}, \ldots, v_{n}]$ 都是正交矩阵, 且 $q_{1} = v_{1}$ . 下面我们只需证明 $q_{i} = v_{i}$ 或 $q_{i} = -v_{i}, i = 2, 3, \ldots, n,$ 即证明

记 $W = [w_{1},w_{2},\ldots ,w_{n}],H = [h_{ij}]$ ，则有

即

所以

因为 $q_{1} = v_{1}$ , 所以 $w_{1} = [1,0,\dots,0]^{\mathsf{T}}$ . 又 $G$ 是上Hessenberg矩阵, 利用归纳法, 我们可以证明 $w_{i}$ 的第 $i + 1$ 到第 $n$ 个元素均为 0. 所以 $W$ 是一个上三角矩阵. 又 $W$ 是正交矩阵, 所以 $W = \mathrm{diag}(1,\pm 1,\ldots ,\pm 1)$ . 由此, 定理结论成立. □

由于 $Q_{k}$ 的其它列都是由 $Q_{k}$ 的第一列唯一确定 (至多相差一个符号), 所以我们只要找到一个正交矩阵 $\tilde{Q}_{k}$ 使得其第一列与 $Q_{k}$ 的第一列相等, 且 $\tilde{Q}_{k}^{\mathsf{T}}A_{k}\tilde{Q}_{k}$ 为上Hessenberg矩阵, 则由隐式Q定理可知 $\tilde{Q}_k = WQ_k$ , 其中 $W = \mathrm{diag}(1,\pm 1,\ldots ,\pm 1)$ . 于是

又 $W^{\mathsf{T}}A_{k + 1}W$ 与 $A_{k + 1}$ 相似,且对角线元素相等,而其它元素也至多相差一个符号,所以不会影响 $A_{k + 1}$ 的收敛性,即下三角元素收敛到0,对角线元素收敛到 $A$ 的特征值.换句话说,在QR迭代法中,如果我们用 $\tilde{Q}_k^\top A_k\tilde{Q}_k$ 代替 $Q_{k}A_{k}Q_{k}^{\top}$ ,即直接令 $A_{k + 1} = \tilde{Q}_k^\top A_k\tilde{Q}_k$ ，则其收敛性与原QR迭代法没有任何区别!这就是隐式QR迭代的基本思想.

由于 $A$ 是上Hessenberg矩阵，因此在计算中可以使用Givens变换，即 $\tilde{Q}_k$ 是一系列Givens变换的乘积.下面我们举一个例子，具体说明如何利用隐式Q定理，由 $A_{1}$ 得到 $A_{2}$

例4.5 设 $A \in \mathbb{R}^{5 \times 5}$ 是一个不可约上Hessenberg矩阵, 即

第一步：构造一个Givens变换 $G_{1}$ ，其转置如下

这里 $G_{1}$ 的第一列 $[c_{1},s_{1},0,\ldots ,0]^{\mathsf{T}}$ 就是 $A_{1} - \sigma_{1}I$ 的第一列 $[a_{11} - \sigma_1,a_{21},0,\dots ,0]^{\mathsf{T}}$ 的单位化后的列向量,其中 $\sigma_{1}$ 是位移.于是有

与 $A_{1}$ 相比较, $A^{(1)}$ 在 (3,1) 位置上多出一个非零元, 我们把它记为 “+”, 并称之为 bulge. 在下面的计算过程中, 我们的目标就是将其“赶”出矩阵, 从而得到一个新的上 Hessenberg 矩阵, 即 $A_{2}$ .

第二步：为了消去这个 bulge, 我们可以构造 Givens 变换

为了保持与原矩阵的相似性, 需要再右乘 $G_{2}$ , 所以

此时, bugle 从 (3,1) 位置被 “赶” 到 (4,2) 位置.

第三步：与第二步类似，构造Givens变换

这时

于是, bugle 又从 (4,2) 位置又被“赶”到 (5,3) 位置.

第四步：再次构造Givens变换

于是

现在, bulge 已经被 “赶” 出矩阵, 且

其中 $\tilde{Q}_1 = G_1G_2G_3G_4$ . 通过直接计算可知, $\tilde{Q}_1$ 的第一列为 $[c_1, s_1, 0, 0, 0]^{\mathsf{T}}$ , 即 $A_1 - \sigma_1I$ 的第一列的单位化. 根据隐式 $\mathbf{Q}$ 定理, $A_2 \triangleq A^{(4)} = \tilde{Q}_1^{\mathsf{T}}A_1\tilde{Q}_1$ 就是我们所需要的矩阵.

如果 $A \in \mathbb{R}^{n \times n}$ 是上Hessenberg矩阵, 则使用上面的算法, 带位移的隐式QR迭代中每一步的运算量为 $6n^{2} + \mathcal{O}(n)$ , 如果需要计算 $Q$ 的话, 则总运算量为 $12n^{2} + \mathcal{O}(n)$ .

4.5.3 位移的选取

在带位移的QR迭代法中, 位移的选取非常重要. 通常, 位移越离某个特征值越近, 收敛速度就越快. 由习题4.5可知, 如果位移 $\sigma$ 与某个特征值非常接近, 则 $A_{k}(n,n) - \sigma$ 就非常接近于0. 这说明 $A_{k}(n,n)$ 通常会首先收敛到 $A$ 的一个特征值, 所以 $\sigma = A_{k}(n,n)$ 是一个不错的选择. 但是, 如果这个特征值是复数, 这种位移选取方法就可能失效.

下面我们介绍一种针对共轭复特征值的位移选取方法, 即双位移策略.

设 $\sigma \in \mathbb{C}$ 是 $A$ 的某个复特征值 $\lambda$ 的一个很好的近似, 则其共轭 $\bar{\sigma}$ 也应该是 $\bar{\lambda}$ 的一个很好的

近似. 因此我们可以考虑双位移策略, 即先以 $\sigma$ 为位移迭代一次, 然后再以 $\bar{\sigma}$ 为位移迭代一次, 如此不断交替进行迭代. 这样就有

容易验证

其中 $Q = Q_{1}Q_{2}$ .我们注意到 $\sigma$ 是复的,所以 $Q_{1}$ 和 $Q_{2}$ 都可能是复矩阵.但我们却可以选取适当的 $Q_{1}$ 和 $Q_{2}$ ，使的 $Q = Q_{1}Q_{2}$ 是实正交矩阵

引理4.9 在双位移QR迭代(4.3)中，我们可以选取酉矩阵 $Q_{1}$ 和 $Q_{2}$ 使得 $Q = Q_{1}Q_{2}$ 是实矩阵. (板书)

证明. 由于

所以

又 $Q_{1}Q_{2}$ 是酉矩阵, $R_{2}R_{1}$ 是上三角矩阵, 故 $(Q_{1}Q_{2})(R_{2}R_{1})$ 是实矩阵 $A_{1}^{2} - 2\mathrm{Re}(\sigma)A_{1} + |\sigma|^{2}I = (A_{1} - \sigma I)(A_{1} - \bar{\sigma}I)$ 的 QR 分解. 所以 $Q_{1}Q_{2}$ 和 $R_{2}R_{1}$ 都可以是实矩阵.

由这个引理可知, 存在 $Q_{1}$ 和 $Q_{2}$ , 使得 $Q = Q_{1}Q_{2}$ 是实正交矩阵, 从而 $A_{3} = Q^{\mathsf{T}}A_{1}Q$ 也是实矩阵. 这时, 在迭代过程 (4.3) 中, 我们无需计算 $A_{2}$ , 可直接由 $A_{1}$ 计算出 $A_{3}$ .

具体计算过程仍然是根据隐式 $Q$ 定理：我们只要找到一个实正交矩阵 $Q$ ，使得其第一列与 $A_1^2 - 2\mathrm{Re}(\sigma)A_1 + |\sigma|^2I$ 的第一列平行，并且 $A_3 = Q^\top A_1Q$ 是上Hessenberg矩阵即可.

由于 $A_1^2 - 2\mathrm{Re}(\sigma)A_1 + |\sigma|^2 I$ 的第一列为

所以 $Q$ 的第一列是上述向量的单位化. 其它过程可以通过隐式 QR 迭代来实现. 但此时的“bulge”是一个 $2 \times 2$ 的小矩阵. 因此, 在双位移隐式 QR 迭代过程中, 我们需要使用 Householder 变换.

需要指出的是, 双位移 QR 迭代法中的运算都是实数运算.

下面我们举一个例子, 具体说明如何在实数运算下实现双位移隐式 QR 迭代法.

例4.6 设 $A \in \mathbb{R}^{6 \times 6}$ 是一个不可约上Hessenberg矩阵, 即

第一步：构造一个正交矩阵 $H_{1} = \left[ \begin{array}{cc}\tilde{H}_{1}^{\mathsf{T}} & 0\\ 0 & I_{3\times 3} \end{array} \right],$ 其中 $\tilde{H}_1\in \mathbb{R}^{3\times 3}$ ，使得其第一列与向量(4.4)平行.于是有

与 $A_{1}$ 相比较, $A^{(1)}$ 在 $(3,1)$ , $(4,1)$ 和 $(4,2)$ 位置上出现 bulge. 在下面的计算过程中, 我们的目标就是把它们“赶”出矩阵, 从而得到一个新的上 Hessenberg 矩阵.

第二步：令 $H_{2} = \left[ \begin{array}{ccc}I_{1\times 1} & 0 & 0\\ 0 & \tilde{H}_{2}^{\mathsf{T}} & 0\\ 0 & 0 & I_{2\times 2} \end{array} \right],$ 其中 $\tilde{H}_2\in \mathbb{R}^{3\times 3}$ 是对应于 $A(2:4,1)$ 的Householder变换，使得

这时, 我们将bugle向右下角方向“赶”了一个位置.

第三步：与第二步类似，令 $H_{3} = \left[ \begin{array}{ccc}I_{2\times 2} & 0 & 0\\ 0 & \tilde{H}_{3}^{\mathsf{T}} & 0\\ 0 & 0 & I_{1\times 1} \end{array} \right],$ 其中 $\tilde{H}_3\in \mathbb{R}^{3\times 3}$ 是对应于 $A(3:5,2)$

的 Householder 变换, 使得

此时, bugle 又被向右下角方向 “赶” 了一个位置.

第四步：令 $H_{4} = \left[ \begin{array}{cc}I_{3\times 3} & 0\\ 0 & \tilde{H}_{4}^{\mathsf{T}} \end{array} \right],$ 其中 $\tilde{H}_4\in \mathbb{R}^{3\times 3}$ 是对应于 $A(4:6,3)$ 的Householder变换，使得

第五步：此时，只需构造一个Givens变换 $G_{5} = \left[ \begin{array}{cc}I_{4\times 4} & 0\\ 0 & G(4,5,\theta)^{\mathsf{T}} \end{array} \right]$ ，使得

现在, bulge 已经被全部消除, 且

其中 $Q = H_{1}H_{2}H_{3}H_{4}G_{5}$ . 通过直接计算可知, $Q$ 的第一列即为 $H_{1}$ 的第一列, 也就是向量 (4.4) 的单位化. 根据隐式 $\mathbf{Q}$ 定理, 可以直接令 $A_{3} \triangleq A^{(5)} = Q^{\mathrm{T}}AQ$ .

最后要考虑位移 $\sigma$ 的具体选取问题.在单位移QR迭代法中，如果 $A$ 的特征值都是实数的话，我们可以取 $\sigma = A_{k}(n,n)$ 推广到复共轭特征值上，我们可以取 $A_{k}$ 的右下角矩阵

的复共轭特征值作为双位移. 这样选取的位移就是 Francis 位移.

如果上述矩阵的两个特征值都是实的，则选取其中模较小的特征值做单位移QR迭代

一般来说, 采用 Francis 位移的 QR 迭代法会使得迭代矩阵的右下角收敛到一个上三角矩阵 (两个实特征值) 或一个 $2 \times 2$ 的矩阵 (一对复共轭特征值), 即 $A_{k}(n, n - 1)$ 或 $A_{k}(n - 1, n - 2)$ 趋向于 0. 带 Francis 位移的隐式 QR 迭代法通常具有二次渐进收敛性, 在实际计算中, 计算一个特征值大约平均迭代两步.

A 作为一种实用的算法, 需要给出一种有效的收敛判别准则, 一种简单实用的准则是

其中 $\varepsilon_{u}$ 是机器精度或给定的阈值

有时也可能是中间某个下次对角线元素 (不是最后两个) 首先收敛到 0, 此时可以转化为计算两个子矩阵 (对角块) 的特征值. 因此在迭代中需要判断所有下次对角线元素是否满足精度要求.
如果同时有多个下次对角线元素满足精度要求, 则可转化为多个子矩阵.

但需要指出的是, QR 迭代并不是对所有矩阵都收敛. 例如:

对于上面的矩阵, 采用Francis位移的QR迭代法无效

另外，也可以考虑多重位移策略，参见[132].

4.5.4 收缩

在隐式QR迭代过程中, 当矩阵 $A_{k + 1}$ 的某个下次对角线元素 $a_{i + 1,i}$ 很小时, 我们可以将其设为0.由于 $A_{k + 1}$ 是上Hessenberg矩阵, 这时 $A_{k + 1}$ 就可以写成分块上三角形式, 其中两个对角块都是上Hessenberg矩阵.因此我们可以将隐式QR迭代作用在这两个规模相对较小的矩阵上, 从而可以大大节约运算量.