摘要：深入剖析QR算法的历史渊源、数学原理及其在现代工程中的核心应用，揭示这一被誉为"20世纪数值线性代数最伟大成就"的迭代方法如何让人类精确计算任意矩阵的全部特征值。
影响人类的50个最重要算法系列：5/50
在上一篇中，我们讲了SVD如何把矩阵"拆骨"成旋转和拉伸的组合。但有一个更基本的问题：一个矩阵的"固有频率"是什么？
这不是比喻。当工程师设计一座桥梁时，他们最担心的事情之一就是"共振"——如果风的频率恰好等于桥梁的固有振动频率，桥就会像被无形之手越推越猛，直到崩塌。1940年塔科马海峡大桥的坍塌，就是这个噩梦的真实上演。
要计算一个结构的固有频率，在数学上就是求解一个矩阵的特征值（Eigenvalue）。而QR算法，就是人类找到的计算任意矩阵全部特征值的最强武器。
一、什么是特征值？为什么它如此重要？
在解释QR算法之前，先说清楚它要解决的问题。
当一个矩阵A作用于某个特殊向量v时，如果结果只是把v放大或缩小了λ倍，而方向完全不变：
$$Av = \lambda v$$
那么v就叫做A的特征向量，λ就叫做对应的特征值。
“特征"这个词翻译得很精准——特征值就是矩阵的"本征属性”，是矩阵最内在的DNA。不管你怎么旋转坐标系、怎么换基底，特征值都不会变。
特征值的物理意义无处不在：
- 结构工程：建筑物、桥梁、飞机机翼的固有振动频率就是刚度矩阵的特征值。如果外力频率接近某个特征值，就会发生共振。
- 量子力学：原子的能级就是哈密顿算子的特征值。电子只能存在于特定能级上，这些能级由特征值方程决定。
- 稳定性分析：一个动力系统是否稳定，取决于其雅可比矩阵的特征值是否都在复平面的左半部分。
- Google PageRank：网页的重要性排名，本质上是互联网链接矩阵的主特征向量。
- 主成分分析：数据的主要变化方向，就是协方差矩阵的特征向量，对应的特征值表示该方向上的方差大小。
所以，计算特征值不是一个抽象的数学练习，而是关乎桥会不会塌、飞机会不会散架、核反应堆会不会失控的生死问题。
二、为什么计算特征值这么难？
你可能会想：特征值不就是解方程 det(A - λI) = 0 吗？把行列式展开成多项式，然后求根不就行了？
理论上没错。对于一个n×n的矩阵，det(A - λI) = 0 是一个n次多项式方程（叫做"特征多项式"），它的n个根就是n个特征值。
但这条路在实践中是死路。
第一个问题：高次多项式求根本身就是噩梦。 阿贝尔-鲁菲尼定理证明了，五次及以上的多项式方程没有通用的代数求根公式。你只能用数值方法近似求解。
第二个问题：特征多项式对系数极其敏感。 矩阵元素的微小扰动，可能导致特征多项式系数的巨大变化，进而导致求出的根完全错误。这在数值计算中叫做"病态问题"。
1961年，数值分析大师威尔金森（James Wilkinson）用一个著名的反例证明了这一点：一个20×20的矩阵，特征值本来是1, 2, 3, …, 20，但如果你先算出特征多项式再求根，由于浮点数舍入误差的放大，得到的"特征值"会偏离真实值几十甚至上百。
所以，所有严肃的数值计算教科书都会告诉你：永远不要通过特征多项式来计算特征值。
那该怎么办？
三、QR算法的诞生：一个28岁研究生的天才之作
1959年，英国国家物理实验室（NPL）的一位28岁研究生约翰·弗朗西斯（John G.F. Francis）提出了QR算法。几乎同时，苏联数学家库布兰诺夫斯卡娅（Vera Kublanovskaya）也独立发现了本质相同的方法。
弗朗西斯的想法优美得令人窒息：不去求解特征多项式，而是通过反复"旋转"矩阵，让它自己慢慢"露出"特征值。
具体来说，QR算法的每一步做两件事：
- QR分解：把当前矩阵A分解为一个正交矩阵Q和一个上三角矩阵R的乘积，即 A = QR
- 反向相乘：把Q和R反过来乘，得到新矩阵 A’ = RQ
然后对A’重复这个过程。
神奇的是：经过足够多次迭代后，矩阵会逐渐收敛为一个（准）上三角矩阵，而对角线上的元素就是原矩阵的特征值！
四、为什么这个"旋转"能逼出特征值？
QR算法的收敛性背后有深刻的数学原因，但核心直觉可以这样理解：
每一次QR迭代，本质上是对矩阵做了一次正交相似变换：
$$A’ = RQ = Q^T A Q$$
相似变换不改变特征值（这是线性代数的基本定理），但它改变了矩阵的"外观"。QR算法的精妙之处在于，它选择的这种特定的相似变换，会让矩阵的下三角部分越来越小，最终趋近于零。
当下三角部分变成零时，矩阵就变成了上三角矩阵——而上三角矩阵的特征值就是对角线元素，直接读出来就行。
打个比方：想象你有一块浑浊的水晶球，里面藏着宝石（特征值）。QR算法就像是一遍遍地打磨水晶球，每磨一次，水晶就更透明一点，宝石的轮廓就更清晰一点。磨到最后，宝石完全暴露出来。
五、工程加速：从理论到实用
原始的QR算法虽然正确，但速度太慢。对于一个n×n的矩阵，每次QR分解就需要O(n³)的计算量，而收敛可能需要很多次迭代。
工程实践中，有三个关键的加速技巧让QR算法变得真正实用：
1. 先化为Hessenberg形式
在开始QR迭代之前，先用Householder变换（这是我们系列中第6个算法的主题）把矩阵化为"上Hessenberg形式"——一种几乎是上三角、只在次对角线上有非零元素的特殊形式。
这一步只需要O(n³)的一次性计算，但之后每次QR迭代的成本从O(n³)降到了O(n²)。这是一个巨大的加速。
2. 带位移的QR迭代
原始QR算法的收敛速度取决于相邻特征值的比值。如果两个特征值很接近，收敛就会极慢。
解决方案是"位移"（Shift）：不对A直接做QR分解，而是对 A - σI 做QR分解，其中σ是对某个特征值的估计（通常取矩阵右下角的元素）。这相当于把目标特征值"拉"到零附近，让收敛速度从线性变成二次甚至三次。
Wilkinson位移策略（取右下角2×2子矩阵的特征值中更接近末尾元素的那个作为位移）几乎总能保证三次收敛，这意味着每次迭代就能让误差缩小到原来的立方。
3. 隐式QR（Francis步）
弗朗西斯还发现，你不需要显式地计算Q和R，而是可以通过一系列Givens旋转或Householder反射来"隐式地"完成QR步。这不仅节省了存储空间，还提高了数值稳定性。
结合这三个技巧，现代QR算法对于一个n×n的矩阵，总计算量大约是O(n³)——和高斯消元法解线性方程组的复杂度相当。这意味着，计算一个1000×1000矩阵的全部1000个特征值，在现代计算机上只需要几秒钟。
六、QR算法的应用：从桥梁到量子
结构工程：防止共振灾难
每一座桥梁、每一栋高楼、每一架飞机在设计阶段，都需要计算其结构的固有频率。工程师建立有限元模型，得到刚度矩阵和质量矩阵，然后用QR算法（或其变体）计算广义特征值问题的全部解。
如果某个固有频率落在风荷载、地震波或发动机振动的频率范围内，设计就必须修改。QR算法让这种"虚拟风洞测试"成为可能。
量子化学：计算分子能级
在量子力学中，分子的电子能级由哈密顿矩阵的特征值决定。对于小分子，这个矩阵可能只有几百维；对于蛋白质或纳米材料，可能达到几万维。QR算法及其变体是计算这些能级的标准工具。
控制系统：判断稳定性
一个线性控制系统是否稳定，完全取决于其系统矩阵的特征值。如果所有特征值的实部都是负的，系统就是稳定的；如果有任何一个特征值的实部为正，系统就会发散失控。
从自动驾驶到工业机器人，从电网调度到航天器姿态控制，QR算法在幕后默默守护着这些系统的稳定性。
Google PageRank：互联网的秩序
虽然PageRank通常用幂迭代法计算（因为只需要主特征向量），但QR算法在分析互联网链接结构的完整谱特性时仍然不可替代。通过分析转移矩阵的全部特征值分布，研究者可以理解网络的社区结构、信息传播模式和脆弱性。
七、一个被遗忘的天才
QR算法的发明者约翰·弗朗西斯的故事，是科学史上最令人唏嘘的篇章之一。
1959年发表那两篇改变数值计算历史的论文时，弗朗西斯只有28岁，还是一个没有博士学位的研究生。发表论文后不久，他就离开了学术界，进入工业界从事与数值计算无关的工作，此后几十年几乎完全消失在公众视野中。
直到2007年，数值分析界的学者们才重新找到了已经退休的弗朗西斯，并在一次学术会议上向他致敬。当时的弗朗西斯已经76岁，对自己年轻时的工作竟然产生了如此深远的影响感到惊讶。
今天，全世界每一台计算机上的LAPACK库（线性代数标准库）中，计算特征值的核心程序仍然是QR算法。每一次有限元分析、每一次振动模态计算、每一次控制系统设计，都在使用弗朗西斯60多年前的发明。
这可能是数学史上"影响力与知名度最不匹配"的算法之一。
八、QR算法与SVD的关系
细心的读者可能已经注意到：上一篇讲的SVD也是矩阵分解，这一篇讲的QR算法也是矩阵分解。它们是什么关系？
答案是：计算SVD的标准方法，底层就是QR算法。
Golub和Kahan在1965年提出的SVD计算方法，核心步骤是：先把矩阵化为双对角形式，然后对双对角矩阵反复做QR迭代，直到收敛。所以QR算法不仅自己能计算特征值，还是计算SVD的基础引擎。
可以说，QR算法是数值线性代数的"发动机"，而SVD是建立在这个发动机之上的"整车"。
📥 系列回顾：