跳到主要内容

QR 算法

一、直观理解：用正交相似变换把矩阵“三角化”​

QR 算法是计算矩阵全部（或大量）特征值的经典现代方法之一。基本思想：反复把矩阵做 QR 分解 Ak=QkRkA_k=Q_k R_k，再反过来乘 Ak+1=RkQkA_{k+1}=R_k Q_k。由于 Ak+1=QkTAkQkA_{k+1}=Q_k^\mathsf{T}A_k Q_k，每步都是正交相似变换，特征值不变；在适当条件下 AkA_k 趋向上三角（或实舒尔型），对角元逼近特征值。

可以把过程理解为：不断用正交变换显露谱信息，而不是直接对特征多项式暴力求根。

二、适用场景与优缺点​

类型说明
适用场景中小规模稠密矩阵的全部特征值；教学理解 eig 背后思想；对称矩阵三对角化后的 QR 迭代。
核心优点数值稳定（正交变换）；可求全部特征值；有位移加速等成熟变体。
局限与缺点朴素实现 O(n3)O(n^3) 每步且步数多；大稀疏矩阵应使用 eigs/迭代法；手写高效 QR 算法工程量大。

实际工作中直接调用 eig/schur，但理解 QR 有助于解释对称矩阵对角化、收敛与稳定性。

三、核心原理​

1. 基本 QR 迭代​

Ak=QkRk,Ak+1=RkQkA_k=Q_k R_k,\quad A_{k+1}=R_k Q_k

2. 对称情形​

对称矩阵可先 Householder 三对角化，再对三对角做隐式 QR，效率大增，且保持对称。

3. 位移（shift）​

用 Ak−μIA_k-\mu I 做 QR 可加速某个特征值收敛，Wilkinson 位移等策略使算法实用化。

4. 实非对称​

收敛到实 Schur 型（分块上三角，对角块 1×1 或 2×2 对应复共轭对）。

四、实践要点​

  • 小矩阵用 eig；
  • 只需少数极端特征值用 eigs；
  • 对称正定可与 Cholesky/svd 联系；
  • 检查正交性损失与残差 ∥AV−VD∥\|AV-VD\|。

五、MATLAB 示例​

%% 朴素 QR 迭代演示（教学用，非生产实现）
A0 = [ 2  1  0;
       1  3  1;
       0  1  4];
A = A0;
nIt = 30;
diagHist = zeros(nIt, 3);

for k = 1:nIt
    [Q, R] = qr(A);
    A = R * Q;
    diagHist(k, :) = diag(A)';
end

lamQR = sort(diag(A), "descend");
lamEig = sort(eig(A0), "descend");
fprintf("QR 迭代对角元 ≈ "); disp(lamQR');
fprintf("eig 特征值     = "); disp(lamEig');
fprintf("差异 = %.3e\n", norm(lamQR - lamEig));

figure;
plot(diagHist, "LineWidth", 1.2);
xlabel("迭代"); ylabel("对角元");
title("朴素 QR 迭代下对角元演化");
legend("a_{11}","a_{22}","a_{33}"); grid on;

% 残差检查：用 eig 的特征向量
[V,D] = eig(A0);
fprintf("||A V - V D|| = %.3e\n", norm(A0*V - V*D));

对称矩阵与正交对角化：

S = A0' * A0;   % SPD
[V,D] = eig(S);
fprintf("正交性 ||V'V-I|| = %.3e\n", norm(V'*V - eye(3)));
fprintf("重构 ||V D V'-S|| = %.3e\n", norm(V*D*V' - S));

六、如何解释结果​

报告特征值列表（按模或大小排序）、是否对称、残差 ∥AV−VD∥\|AV-VD\|、计算所用函数。若解释模态，把特征向量与实际变量对应，并注意归一化与符号。

七、常见误区​

1. 把朴素 QR 迭代当成工业实现​

真实 eig 含平衡、Hessenberg 化、隐式位移等。

2. 大稀疏矩阵调用稠密 eig​

内存与时间都可能不可接受。

3. 忽略复特征值成对出现​

实矩阵的复谱共轭对称。

4. 特征向量未归一/未对齐就比较​

应比较子空间或夹角。

附：建模流程中的落地清单（QR 算法）​

在竞赛或课程论文中，仅写出定义式往往不够。建议按下列顺序组织材料，既能保证可复现，也便于评委快速抓住贡献。

1. 问题对齐​

先用两三句话说明：题目中的哪个子问题对应本方法，输入是什么，输出要支持什么决策。避免“因为常见所以使用”。若存在更简单的替代方法，应说明为何不够用。

2. 符号与单位​

列出核心符号表，强调单位一致。许多错误并非算法失败，而是日/年、米/千米、概率/百分比混用。对无量纲参数，说明取值来源（估计、文献、情景设定）。

3. 关键公式回顾​

基本迭代 A=QR, A←RQA=QR,\ A\leftarrow RQ 保持相似，对角元趋向特征值（教学版）。

4. 计算与实现要点​

  • 固定随机种子（若涉及随机）；
  • 记录中间量以便手算核对小例子；
  • 用极限情形或解析特例做冒烟测试；
  • 保存主要输入输出到表格，便于敏感性分析复用。

5. 补充 MATLAB 片段​

A=hilb(4); for k=1:40, [Q,R]=qr(A); A=R*Q; end; disp(sort(diag(A))); disp(sort(eig(hilb(4))));

6. 结果报告模板​

建议至少报告：特征值列表、残差 ∥AV−VD∥\|AV-VD\|、是否对称、所用库函数。

7. 何时不该使用​

把朴素 QR 迭代当作工业实现细节全部；或对大稀疏矩阵调用稠密 eig。

完成以上清单后，再进入“小结”中的一句话金句，通常论文的方法段与验证段就会扎实许多。

八、小结​

QR 算法通过反复 QR 分解与正交相似变换显露矩阵的特征值结构，是稠密特征值计算的理论与实践支柱。建模中多数时候调用可靠库函数，但理解 QR 有助于正确使用 eig、诊断残差，并在谱方法、PCA、振动与稳定性分析中准确解释结果。