14. 项目案例：实 Francis QR 算法

本章把正交相似变换、Hessenberg 结构、隐式位移和数值验收串成一条完整工程链。目标是计算实 Schur 分解

$$A\approx UT{U}^T,$$

其中 $U$ 正交，$T$ 为实准上三角矩阵：对角上允许出现代表共轭复特征值的 $2\times 2$ 块。

14.1 总体流程

实现可以分成五层：

1. 用 Householder 反射将 $A$ 化为上 Hessenberg 矩阵 $H$；
2. 找到尚未收敛的活动子块；
3. 对该子块执行 Francis 双位移步并追赶凸起；
4. 检测可忽略的次对角元并做消去；
5. 对已隔离的 $1\times 1$ 或 $2\times 2$ 块整理标准形。

每次左、右正交变换都必须同步累积到 $U$，否则最终只能得到特征值，无法验证 Schur 分解。

14.2 为什么先化 Hessenberg 形

若直接对稠密矩阵反复做 QR 分解，每步成本是 $O(n^3)$。Hessenberg 矩阵只有第一条次对角线可能非零，一次隐式 QR 扫描只需 $O(n^2)$，且相似变换保持结构。

一次性的 Hessenberg 化为 $O(n^3)$；随后若需要 $O(n)$ 次扫描，总体仍是典型的 $O(n^3)$ 算法。

14.3 隐式双位移

取活动子块尾部 $2\times 2$ 矩阵的两个特征值 ${\mu }_1,{\mu }_2$。即便它们是共轭复数，和与积仍为实数。定义

$$p(H)=(H-{\mu }_{1}I)(H-{\mu }_{2}I)=H^2-sH+tI,$$

其中 $s={\mu }_1+{\mu }_2$、$t={\mu }_1{\mu }_2$。

算法无需显式形成 $p(H)$。只取 $p(H)e_1$ 的前三个分量，构造一个 Householder 变换把它对准 $e_1$，便在 Hessenberg 带外制造一个很小的“凸起”。隐式 Q 定理保证这一过程等价于带双位移的 QR 步。

14.4 凸起追赶

凸起出现后，依次对局部的三维或二维向量构造 Householder 变换，从左边消去带外元素，再从右边施加同一正交变换。凸起沿次对角线向右下移动，最终离开活动子块。

实现难点不在反射公式，而在索引边界：

• 左更新只触及局部行与其右侧；
• 右更新要覆盖相应列和必要的上方行；
• 累积矩阵 $U$ 必须做完全一致的右更新；
• 每步后主动清理理论上应为零的带外小量。

14.5 消去与 $2\times 2$ 实块

当

$$|h_{i+1,i}|\le \tau (|h_{ii}|+|h_{i+1,i+1}|)$$

时，可将该次对角元置零并分裂问题。阈值 $\tau$ 应与机器精度及矩阵尺度相适应。

实矩阵若有非实特征值，不能被实正交相似变换化成真正的上三角矩阵。因此算法必须接受收敛的 $2\times 2$ 块，而不是强迫所有次对角元都变为零。

14.6 三类验收指标

项目实现不能只比较特征值。至少检查：

正交性

$$\parallel U^{T}U-I{\parallel }_F.$$

后向残差

$$\frac{\parallel U^{T}AU-T{\parallel }_F}{\parallel A{\parallel }_F}.$$

结构残差

检查 $T$ 的第一条次对角线以下是否接近零，并确认未消去的次对角元只属于合法的 $2\times 2$ 块。

在普通随机矩阵上，可靠实现通常可把分解残差控制到机器精度乘以温和维数因子的量级；随着规模和非正规性上升，正交性与结构误差会更敏感。这正是需要同时报告多种指标的原因。

14.7 从数学算法到可靠程序

• 把“寻找活动块”“单次双位移”“消去判断”拆成独立函数，便于单测。
• 对 $1\times 1$、$2\times 2$、已分块和近消去矩阵设置专门测试。
• 设置每个活动块的最大扫描次数，并在失败时报告位置和残差。
• 用尺度归一化后的误差比较不同矩阵，避免绝对误差误导。
• 以 NumPy 或 LAPACK 的 Schur 结果作参照，但验证重点应是自身分解恒等式。

14.8 双位移起始向量的显式公式

设活动块左上角从索引 $l$ 开始，右下 $2\times 2$ 块的迹与行列式为 $s,t$。无需形成 $H^2$，$p(H)e_l$ 的非零部分只涉及前三个分量：

$$x=h_{ll}^2+h_{l,l+1}{h}_{l+1,l}-s{h}_{ll}+t,$$$$y=h_{l+1,l}(h_{ll}+h_{l+1,l+1}-s),z=h_{l+1,l}{h}_{l+2,l+1}.$$

对 $(x,y,z{)}^T$ 构造 Householder 反射即可启动凸起。实现前应按 $|x|+|y|+|z|$ 缩放，避免极大或极小元素导致溢出/下溢。

14.9 一次凸起追赶的更新范围

在位置 $k$ 构造长度至多为 3 的反射 $P_k$。若活动块为 $l:r$，更新遵循

$$H_{k:k+2,,k-1:r}\leftarrow P_k^T{H}_{k:k+2,,k-1:r},$$$$H_{l:min(r,k+3),,k:k+2}\leftarrow H_{l:min(r,k+3),,k:k+2}{P}_k,$$

并累积

$$U_{:,,k:k+2}\leftarrow U_{:,,k:k+2}{P}_k.$$

边界处反射长度降为 2。右更新若截得过短，会破坏相似关系；截得过长虽仍正确，却浪费计算。带外理论零应在局部变换结束后清理，而不是在更新前随意截断。

14.10 $2\times 2$ 块与实 Schur 标准化

对实块

$$B=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right),$$

判别式 $\mathrm{\Delta }=(a-d{)}^2+4bc$ 决定特征值是否为实数。若 $\mathrm{\Delta }\ge 0$，可用稳定 Givens 旋转把块进一步三角化；若 $\mathrm{\Delta }<0$，应保留代表共轭对的 $2\times 2$ 块，并通过正交相似变换把它整理成一致的准三角形式。

直接使用求根公式可能发生相消。计算实特征值时，应先稳定求出较大根，再用行列式关系得到另一根。

14.11 收敛保护与异常位移

实际 QR 迭代可能暂时停滞。可靠实现通常需要：

• 局部尺度化的亏损检查；
• 对极小活动块的专门处理；
• 多次未亏损时使用 exceptional shift；
• 每个活动块的迭代上限和清晰失败信息；
• 对 NaN、无穷和过大缩放的前置检测。

异常位移不是改变目标特征值，而是暂时改变迭代多项式，帮助凸起摆脱不利的局部结构。

14.12 伪代码骨架

H, U = hessenberg_reduction(A)
r = n - 1
while r >= 0:
    寻找活动块 l:r
    若 r == l: 接受 1x1 块，r -= 1
    否则若块可作为 2x2 接受: 标准化并令 r -= 2
    否则:
        从尾部 2x2 块取得双位移的迹 s 与行列式 t
        构造 p(H)e_l 的前三个分量
        启动并追赶凸起，同时更新 U
        依据局部阈值执行亏损
return U, H

14.13 测试矩阵族

仅测试随机正态矩阵不够。至少覆盖：

• 已经是上三角或 Hessenberg 的矩阵；
• 含多个独立块和接近亏损次对角元的矩阵；
• 具有共轭复特征值对的实矩阵；
• 重根、近重根与高度非正规矩阵；
• $1\times 1$、$2\times 2$、零矩阵和尺度差异极大的边界输入。

测试结果要联合判断正交性、相似残差、结构残差和是否在有限扫描内完成。

14.14 这一算法串起了什么

Francis QR 同时使用了数值代数最核心的思想：正交变换保证稳定性，结构化约简降低复杂度，位移加速收敛，亏损把全局问题拆成局部问题，而后向误差把“算出一个答案”提升为“知道答案可信”。

自检

1. 为什么双位移可以始终在实数运算中完成？
2. Hessenberg 化为何能把一次 QR 扫描降到 $O(n^2)$？
3. 验证 Schur 分解为什么不能只比较特征值？
4. 起始向量为什么只需要三个分量？
5. 哪些测试矩阵能暴露随机矩阵不容易触发的错误？