QR 算法、Shift、Hessenberg 与 SVD

前面讨论 Krylov subspace methods 时，我们关心的是如何用

$$\mathbf{b},\ \mathbf{A}\mathbf{b},\ \mathbf{A}^2\mathbf{b},\dots$$

构造一个低维子空间，然后在这个子空间里近似求解线性系统或者近似捕捉谱信息

这一章换一个角度：如果目标不是解

$$\mathbf{A}\mathbf{x}=\mathbf{b}$$

而是求矩阵 $\mathbf{A}$ 的 eigenvalues，那么 QR 分解、Hessenberg 形式、shift 和 Krylov 方法之间会自然连在一起

核心主线可以先写成

$$\boxed{ \mathbf{A} \quad \longrightarrow \quad \mathbf{H} \quad \longrightarrow \quad \text{shifted QR iteration} \quad \longrightarrow \quad \text{eigenvalues} }$$

其中 $\mathbf{H}$ 是一个和 $\mathbf{A}$ 相似的 Hessenberg 矩阵。对于对称矩阵，$\mathbf{H}$ 会进一步简化成三对角矩阵

这条路线背后的思想是：我们不断改变坐标系，但不改变特征值；坐标系越来越接近特征向量方向时，矩阵就越来越接近上三角形式，最后特征值会出现在对角线上

从 QR 分解到 QR iteration

普通 QR 分解只是把一个矩阵分解成

$$\mathbf{A}_0=\mathbf{Q}_0\mathbf{R}_0$$

其中 $\mathbf{Q}_0$ 是正交矩阵，$\mathbf{R}_0$ 是上三角矩阵

这里要特别注意一点：QR 分解的数学定义虽然可以用 Gram-Schmidt 正交化来理解，但实际数值计算里一般不会直接使用最朴素的 Gram-Schmidt

原因是 Gram-Schmidt 在数值上可能损失正交性，尤其当列向量接近线性相关时，舍入误差会让得到的 $\mathbf{Q}$ 不够正交

实际算法更常用的是

• Householder reflections
• Givens rotations

Householder reflections 更适合一次性把一整列下面的元素消掉，常用于稠密矩阵的 QR 分解和 Hessenberg reduction；Givens rotations 则适合一次消掉一个指定元素，常用于稀疏结构、局部更新，以及隐式 QR step 中追逐非零元素

所以在理解 QR iteration 时，可以把 Gram-Schmidt 当作概念入口，但不要把它当成实际算法的默认实现

如果只停在这一步，它只是一个分解。但 QR eigenvalue algorithm 的关键是把顺序反过来，定义

$$\mathbf{A}_1=\mathbf{R}_0\mathbf{Q}_0$$

然后继续重复

$$\mathbf{A}_k=\mathbf{Q}_k\mathbf{R}_k$$ $$\mathbf{A}_{k+1}=\mathbf{R}_k\mathbf{Q}_k$$

这个过程叫 QR iteration

表面上看，它只是不断拆开再乘回去；但乘回去的顺序变了，所以矩阵元素会变化

关键问题是：这样做会不会改变特征值？

由于

$$\mathbf{A}_k=\mathbf{Q}_k\mathbf{R}_k$$

所以

$$\mathbf{R}_k=\mathbf{Q}_k^\top\mathbf{A}_k$$

代入下一步定义：

$$\begin{align*} \mathbf{A}_{k+1} &= \mathbf{R}_k\mathbf{Q}_k \\ &= \mathbf{Q}_k^\top\mathbf{A}_k\mathbf{Q}_k \end{align*}$$

因此 $\mathbf{A}_{k+1}$ 和 $\mathbf{A}_k$ 是相似矩阵

$$\mathbf{A}_{k+1} = \mathbf{Q}_k^\top\mathbf{A}_k\mathbf{Q}_k$$

相似矩阵有相同的特征值，所以整个序列

$$\mathbf{A}_0,\mathbf{A}_1,\mathbf{A}_2,\dots$$

虽然元素不断变化，但特征值始终不变

这就是 QR iteration 可以用来求特征值的第一个原因：

$$\boxed{ \text{每一步都是正交相似变换，因此不改变 eigenvalues} }$$

为什么最后看对角线

如果矩阵 $\mathbf{A}$ 是实对称矩阵，那么它可以被正交对角化：

$$\mathbf{A} = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^\top$$

其中

$$\mathbf{\Lambda} = \operatorname{diag}(\lambda_1,\lambda_2,\dots,\lambda_n)$$

如果 QR iteration 逐渐把当前坐标系调整到特征向量方向，那么矩阵在这个坐标系下就会越来越接近对角矩阵

$$\mathbf{A}_k \longrightarrow \mathbf{\Lambda}$$

这时对角线上的元素就是特征值

对于一般非对称矩阵，情况稍微不同。一般矩阵不一定能被正交对角化，但可以被化到 Schur form

$$\mathbf{A} = \mathbf{U}\mathbf{T}\mathbf{U}^\top$$

其中 $\mathbf{T}$ 是上三角矩阵

$$\mathbf{T} = \begin{bmatrix} \lambda_1 & * & * \\ 0 & \lambda_2 & * \\ 0 & 0 & \lambda_3 \end{bmatrix}$$

上三角矩阵的特征值就是对角线元素，因为

$$\det(\mathbf{T}-\lambda\mathbf{I}) = \prod_i (t_{ii}-\lambda)$$

所以更准确地说，QR iteration 的目标不是一定把矩阵变成对角矩阵，而是把它推向上三角形式

$$\boxed{ \text{一般矩阵趋向 Schur form；对称矩阵进一步趋向对角形式} }$$

非对角元素为什么会变小

QR iteration 可以理解成一种同时追踪多个方向的 power method

普通 power method 从一个向量出发，不断做

$$\mathbf{x},\ \mathbf{A}\mathbf{x},\ \mathbf{A}^2\mathbf{x},\dots$$

如果

$$|\lambda_1|>|\lambda_2|\geq \cdots$$

并且初始向量在 $\mathbf{v}_1$ 方向上有非零分量，那么反复乘 $\mathbf{A}$ 以后，最大的特征方向会逐渐占主导

QR iteration 做的是更完整的事情：它不是只追踪一个向量，而是追踪一组正交向量，并且每一步都重新正交化

可以把这个过程粗略理解为

$$\mathbf{A}^k \quad \text{的信息被不断重新整理到一组正交基里}$$

当这组正交基越来越接近特征向量基时，矩阵在这组基底下就越来越少发生“方向混合”

在特征向量基底中有

$$\mathbf{A}\mathbf{v}_i=\lambda_i\mathbf{v}_i$$

这表示每个方向只缩放自己，不混到别的方向。因此矩阵表示会接近对角或者上三角形式，非对角元素就会变小

所以非对角元素可以理解为不同特征方向之间的混合量。当坐标轴越来越贴近特征方向，这种混合量就逐渐消失

Shifted QR iteration

朴素 QR iteration 理论上很漂亮，但实际收敛可能很慢。尤其当两个特征值大小很接近时，类似 power method 里的收敛比率会很不理想

为了加速收敛，引入 shift

给定一个数 $\mu_k$，不直接对 $\mathbf{A}_k$ 做 QR 分解，而是对

$$\mathbf{A}_k-\mu_k\mathbf{I}$$

做 QR 分解：

$$\mathbf{A}_k-\mu_k\mathbf{I} = \mathbf{Q}_k\mathbf{R}_k$$

然后定义

$$\mathbf{A}_{k+1} = \mathbf{R}_k\mathbf{Q}_k+\mu_k\mathbf{I}$$

看起来好像只是先减去 $\mu_k\mathbf{I}$，最后又加回来，但重点在于 QR 分解得到的 $\mathbf{Q}_k$ 已经变了

它仍然保持相似关系。由

$$\mathbf{A}_k-\mu_k\mathbf{I} = \mathbf{Q}_k\mathbf{R}_k$$

得到

$$\mathbf{R}_k = \mathbf{Q}_k^\top(\mathbf{A}_k-\mu_k\mathbf{I})$$

因此

$$\begin{align*} \mathbf{A}_{k+1} &= \mathbf{R}_k\mathbf{Q}_k+\mu_k\mathbf{I} \\ &= \mathbf{Q}_k^\top(\mathbf{A}_k-\mu_k\mathbf{I})\mathbf{Q}_k +\mu_k\mathbf{I} \\ &= \mathbf{Q}_k^\top\mathbf{A}_k\mathbf{Q}_k - \mu_k\mathbf{Q}_k^\top\mathbf{Q}_k + \mu_k\mathbf{I} \\ &= \mathbf{Q}_k^\top\mathbf{A}_k\mathbf{Q}_k \end{align*}$$

所以 shifted QR 仍然不改变特征值

shift 的作用不是改变答案，而是改变每一步选择的正交变换，使它更快地把某个特征值分离出来

如果 $\mu_k$ 很接近某个特征值 $\lambda_j$，那么

$$\lambda_j-\mu_k$$

会变得很小。对 $\mathbf{A}_k-\mu_k\mathbf{I}$ 做 QR 分解时，这个方向会被更有针对性地处理，相关的次对角元素会更快下降

这就是为什么实际算法通常使用 shifted QR，而不是最朴素的 QR iteration

对于对称三对角矩阵，常用的 Wilkinson shift 往往会带来非常快的局部收敛。直观上，如果某个角度误差本来像

$$\sin\theta$$

一次迭代后可能变成类似

$$\sin^3\theta$$

那么当误差已经小的时候，它会极快地被压下去

Deflation

QR iteration 过程中，我们通常观察次对角线元素

$$a_{i+1,i}$$

如果某个次对角元素已经非常小，比如

$$|a_{i+1,i}|\approx 0$$

那么矩阵就近似分裂成块上三角形式

$$\mathbf{A}_k \approx \begin{bmatrix} \mathbf{B} & * \\ \mathbf{0} & \mathbf{C} \end{bmatrix}$$

块上三角矩阵的特征值等于对角块的特征值合在一起，所以可以把问题拆开

$$\lambda(\mathbf{A}_k) = \lambda(\mathbf{B})\cup\lambda(\mathbf{C})$$

这一步叫 deflation

如果右下角已经近似变成一个 $1\times 1$ 块

$$\begin{bmatrix} \mathbf{B} & * \\ \mathbf{0} & \lambda_n \end{bmatrix}$$

那么 $\lambda_n$ 就可以认为已经被求出来了。接下来只需要继续对左上角的 $\mathbf{B}$ 做 QR iteration

因此 QR 算法不是一口气把所有特征值同时“算完”，而是不断让某些次对角元素变小，然后把已经收敛的部分切出来

为什么先化成 Hessenberg

如果每一步都对完整的稠密矩阵做 QR 分解，代价很高

实际算法通常先用 Householder 变换把一般矩阵化成 Hessenberg 形式

$$\mathbf{H} = \begin{bmatrix} * & * & * & * \\ * & * & * & * \\ 0 & * & * & * \\ 0 & 0 & * & * \end{bmatrix}$$

也就是主对角线下一条以外的元素全为 $0$

这个化简通过正交相似变换完成

$$\mathbf{H} = \mathbf{Q}^\top\mathbf{A}\mathbf{Q}$$

所以 $\mathbf{H}$ 和 $\mathbf{A}$ 有相同特征值

Hessenberg 形式的好处是：它足够简单，但又能在 QR step 中被保持。也就是说，如果 $\mathbf{H}_k$ 是 Hessenberg 矩阵，那么一次 QR iteration 后得到的 $\mathbf{H}_{k+1}$ 仍然是 Hessenberg 矩阵

于是后续每一步 shifted QR 都可以利用这个结构，避免处理大量本来就应该为 $0$ 的元素

对于实对称矩阵，Hessenberg 形式会进一步退化成三对角矩阵

$$\mathbf{T} = \begin{bmatrix} * & * & 0 & 0 \\ * & * & * & 0 \\ 0 & * & * & * \\ 0 & 0 & * & * \end{bmatrix}$$

这是因为对称性要求

$$\mathbf{T}^\top=\mathbf{T}$$

而 Hessenberg 结构只允许主对角线下一条非零。两者合起来，就只剩主对角线和上下两条副对角线

所以实际求对称矩阵特征值的流程通常是

$$\boxed{ \mathbf{A} \xrightarrow{\mathbf{Q}^\top\mathbf{A}\mathbf{Q}} \mathbf{T} \xrightarrow{\text{shifted QR}} \mathbf{\Lambda} }$$

和 Krylov 方法的连接

Hessenberg 形式也正是 Krylov 方法里自然出现的结构

从一个向量 $\mathbf{b}$ 出发，构造 Krylov 子空间

$$K_m(\mathbf{A},\mathbf{b}) = \operatorname{span} \{\mathbf{b},\mathbf{A}\mathbf{b},\mathbf{A}^2\mathbf{b},\dots,\mathbf{A}^{m-1}\mathbf{b}\}$$

Arnoldi 过程会构造一组正交基

$$\mathbf{Q}_m=[\mathbf{q}_1,\mathbf{q}_2,\dots,\mathbf{q}_m]$$

使得

$$\operatorname{span}(\mathbf{Q}_m) = K_m(\mathbf{A},\mathbf{b})$$

并得到关系

$$\mathbf{A}\mathbf{Q}_m = \mathbf{Q}_{m+1}\bar{\mathbf{H}}_m$$

等价地，可以写成

$$\mathbf{A}\mathbf{Q}_m = \mathbf{Q}_m\mathbf{H}_m +h_{m+1,m}\mathbf{q}_{m+1}\mathbf{e}_m^\top$$

其中 $\mathbf{H}_m$ 是一个 $m\times m$ 的小 Hessenberg 矩阵

如果忽略最后那个残差项，那么可以把

$$\mathbf{H}_m \approx \mathbf{Q}_m^\top\mathbf{A}\mathbf{Q}_m$$

理解为 $\mathbf{A}$ 在 Krylov 子空间上的投影

于是求大矩阵 $\mathbf{A}$ 的部分特征值，可以变成求小矩阵 $\mathbf{H}_m$ 的特征值

$$\mathbf{H}_m\mathbf{y} = \theta\mathbf{y}$$

再令

$$\mathbf{x} = \mathbf{Q}_m\mathbf{y}$$

那么 $\theta$ 就是 $\mathbf{A}$ 的近似特征值，$\mathbf{x}$ 是对应的近似特征向量。这些 $\theta$ 叫 Ritz values

所以 QR 和 Krylov 的关系可以概括为

完整 QR：先把整个 $\mathbf{A}$ 化成 Hessenberg，再求全部特征值。

Krylov / Arnoldi：只在低维子空间里得到小 Hessenberg，再求部分特征值。

二者都离不开 Hessenberg 结构，只是规模不同

奇异值和 bidiagonal QR

特征值问题讨论的是

$$\mathbf{A}\mathbf{x} = \lambda\mathbf{x}$$

而奇异值来自 SVD

$$\mathbf{A} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\top$$

其中

$$\mathbf{\Sigma} = \operatorname{diag}(\sigma_1,\sigma_2,\dots)$$

奇异值和特征值的联系是

$$\mathbf{A}^\top\mathbf{A} = \mathbf{V}\mathbf{\Sigma}^2\mathbf{V}^\top$$

所以

$$\lambda_i(\mathbf{A}^\top\mathbf{A}) = \sigma_i^2$$

也就是

$$\boxed{ \sigma_i = \sqrt{\lambda_i(\mathbf{A}^\top\mathbf{A})} }$$

从理论上看，可以先求 $\mathbf{A}^\top\mathbf{A}$ 的特征值，再开平方得到奇异值

但实际数值计算通常不显式形成 $\mathbf{A}^\top\mathbf{A}$，因为这会把条件数平方

$$\kappa(\mathbf{A}^\top\mathbf{A}) = \kappa(\mathbf{A})^2$$

小奇异值方向的误差会因此被放大

更稳定的做法是先通过左右正交变换把 $\mathbf{A}$ 化成双对角矩阵

$$\mathbf{B} = \mathbf{U}^\top\mathbf{A}\mathbf{V}$$

其中

$$\mathbf{B} = \begin{bmatrix} d_1 & e_1 & 0 & 0 \\ 0 & d_2 & e_2 & 0 \\ 0 & 0 & d_3 & e_3 \\ 0 & 0 & 0 & d_4 \end{bmatrix}$$

注意这里不是普通的相似变换，因为左右两边的正交矩阵一般不同

普通特征值相似变换是

$$\mathbf{Q}^\top\mathbf{A}\mathbf{Q}$$

它保持特征值

而 SVD 里的左右正交变换

$$\mathbf{U}^\top\mathbf{A}\mathbf{V}$$

保持的是奇异值

原因是

$$\begin{align*} \mathbf{B}^\top\mathbf{B} &= (\mathbf{U}^\top\mathbf{A}\mathbf{V})^\top (\mathbf{U}^\top\mathbf{A}\mathbf{V}) \\ &= \mathbf{V}^\top\mathbf{A}^\top\mathbf{U} \mathbf{U}^\top\mathbf{A}\mathbf{V} \\ &= \mathbf{V}^\top\mathbf{A}^\top\mathbf{A}\mathbf{V} \end{align*}$$

所以 $\mathbf{B}^\top\mathbf{B}$ 和 $\mathbf{A}^\top\mathbf{A}$ 相似，具有相同特征值

因此

$$\sigma_i(\mathbf{B}) = \sigma_i(\mathbf{A})$$

接下来实际算法并不显式对 $\mathbf{B}^\top\mathbf{B}$ 做 QR，而是在 $\mathbf{B}$ 上做隐式的 shifted QR，也就是 bidiagonal QR / Golub-Kahan SVD step

可以这样理解：

理论目标：求 $\mathbf{B}^\top\mathbf{B}$ 的特征值。

实际操作：直接迭代 $\mathbf{B}$，避免形成 $\mathbf{B}^\top\mathbf{B}$。

每一步都是

$$\mathbf{B}_{k+1} = \mathbf{U}_k^\top\mathbf{B}_k\mathbf{V}_k$$

因此奇异值不变；但双对角线上的非对角元素会逐渐被消掉

$$\mathbf{B}_k \longrightarrow \mathbf{\Sigma}$$

最后对角线上的非负数就是奇异值

总结

QR eigenvalue algorithm 的关键不是 QR 分解本身，而是 QR 分解之后把顺序反过来

$$\mathbf{A}_k=\mathbf{Q}_k\mathbf{R}_k, \qquad \mathbf{A}_{k+1}=\mathbf{R}_k\mathbf{Q}_k$$

每一步都等价于

$$\mathbf{A}_{k+1} = \mathbf{Q}_k^\top\mathbf{A}_k\mathbf{Q}_k$$

所以特征值保持不变

随着迭代进行，矩阵逐渐接近 Schur form；如果矩阵对称，则逐渐接近对角形式。于是对角线元素给出特征值

实际算法会先把矩阵化成 Hessenberg 形式，或者在对称情形下化成三对角形式，再做 shifted QR。Hessenberg 让每一步便宜，shift 让收敛更快，deflation 让已经收敛的特征值可以被分离出去

Krylov / Arnoldi 方法也会产生 Hessenberg 矩阵。区别在于：完整 QR 处理整个空间，目标是全部特征值；Krylov 方法只处理低维子空间，目标通常是大型稀疏矩阵的一部分重要特征值

而对于奇异值，理论上可以通过 $\mathbf{A}^\top\mathbf{A}$ 的特征值求平方根；实际算法则更倾向于先做双对角化，再在双对角矩阵上做隐式 shifted QR，从而避免显式形成条件数更差的 $\mathbf{A}^\top\mathbf{A}$