QR算法

QR 算法（QR algorithm）通过反复作 QR 分解和交换两个因子的次序，构造与原矩阵相似的一列矩阵，以计算特征值。实际特征值软件通常先化为 Hessenberg 或三对角形式，再结合位移、隐式变换和消去小耦合。基础的无位移迭代适合说明原理，却不代表现代软件的全部实现。

一次分解与交换

设 $A_0=\left(\begin{array}{cc}2& 1\\ 1& 2\end{array}\right)$。对第一列归一化得 $q_1=(2,1{)}^T/\sqrt{5}$，取与它正交的 $q_2=(-1,2{)}^T/\sqrt{5}$，于是 $$Q_0=\frac{1}{\sqrt{5}}\left(\begin{array}{cc}2& -1\\ 1& 2\end{array}\right),R_0=\left(\begin{array}{cc}\sqrt{5}& 4/\sqrt{5}\\ 0& 3/\sqrt{5}\end{array}\right),A_0=Q_0{R}_0.$$ 交换顺序后得到 $$A_1=R_0{Q}_0=\left(\begin{array}{cc}14/5& 3/5\\ 3/5& 6/5\end{array}\right).$$ 非对角项从一降为 $0.6$，对角项从两个二变成 $2.8,1.2$，接近真实特征值三和一。但它们在第一步还不是特征值；例如迹仍为四、行列式仍为三，正是矩阵整体保留的量。

相似变换保留特征值

对一般实方阵，每轮求 $A_k=Q_k{R}_k$，其中 $Q_k^T{Q}_k=I$，再令 $A_{k+1}=R_k{Q}_k$。由 $R_k=Q_k^T{A}_k$ 得 $$A_{k+1}=Q_k^T{A}_k{Q}_k.$$ 这是正交相似变换，所以 $$det(\lambda I-A_{k+1})=det(Q_k^T(\lambda I-A_k)Q_k)=det(\lambda I-A_k).$$ 两侧特征多项式完全相同，因而所有特征值及代数重数均保留。若初始矩阵对称，各轮仍对称；累积 $Z_k=Q_0\cdots Q_{k-1}$，有 $A_k=Z_k^T{A}_0{Z}_k$。若极限对角化成立，$Z_k$ 的列给出特征向量。

这个证明保证每轮没有改变特征值，却不保证任意矩阵的无位移迭代都收敛为对角矩阵。LAPACK 的文档将一般实矩阵的目标表述为实 Schur 形：一阶和二阶对角块分别表示实特征值与复共轭特征值。

本例可以写出精确迭代公式，因而图中的收敛速度可逐项复算；它不作为所有矩阵的统一收敛率。

本例的精确收敛公式

矩阵 $A_0$ 的正交特征方向为 $(1,1{)}^T$ 与 $(1,-1{)}^T$，特征值分别为三与一。正交迭代的第一方向与 $A_0^k{e}_1$ 平行，而 $$A_0^k{e}_1=\frac{1}{2}\left(\begin{array}{c}{3}^k+1\\ 3^k-1\end{array}\right).$$ 这一关系也可由各轮 QR 相乘得到：$A_0^k=Z_k{T}_k$，其中 $T_k$ 上三角，其第一列只含首个非零分量，因此第一列方向由 $Z_k$ 的第一列给出。

令 $r=3^k$，把归一化后的 $(r+1,r-1{)}^T$ 及其正交方向代入 $Z_k^T{A}_0{Z}_k$，可得 $$A_k=\left(\begin{array}{cc}3-2/(9^k+1)& 2\cdot 3^k/(9^k+1)\\ 2\cdot 3^k/(9^k+1)& 1+2/(9^k+1)\end{array}\right)$$ （采用本例连续的正交列符号约定）。因此非对角项按 $O(3^{-k})$ 衰减，对角误差按 $O(9^{-k})$ 衰减。QR 分解列符号的另一种合法选择可能改变非对角项符号，而不改变其模和特征值。

有解变式取旋转矩阵 $\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)$。它的特征值为 $\pm i$，不可能用实正交相似变换得到实对角矩阵；无位移 QR 也不会把这个二阶耦合消成两个实数。正确的实 Schur 目标应保留二阶块。

位移、结构化计算与消去耦合

位移 QR 先分解 $A_k-{\mu }_{k}I=Q_k{R}_k$，再令 $A_{k+1}=R_k{Q}_k+{\mu }_{k}I$，相似变换证明仍成立。合适的位移可加快尾部非对角项趋零；对称三对角问题常用 Wilkinson 位移。接近零的耦合经尺度检验后可分裂成两个更小的特征值问题，称为 deflation。

实际计算不会反复对完整稠密矩阵作普通 QR。一般矩阵先作 Hessenberg 化，对称矩阵先作三对角化，再用保持结构的隐式迭代降低每轮代价。舍入下检查分裂阈值和最终残差，较只查看某几个对角数字更能反映计算是否可靠。NumPy 的 QR 接口只提供分解本身，并不等同于完整特征值驱动程序。

对称矩阵的可运行演示

代码只用于实对称矩阵的无位移演示，返回对角近似、累积正交矩阵和轮数，不修改输入。以非对角 Frobenius 范数作为停止条件，轮数耗尽则报告失败。稠密矩阵每轮 $O(n^3)$、存储 $O(n^2)$；这是教学实现的成本，不是结构化生产 QR 的成本。

Python 3

import numpy as np


def symmetric_qr_demo(A, tol=1e-12, maxiter=500):
    A = np.array(A, dtype=float, copy=True)
    if A.ndim != 2 or A.shape[0] != A.shape[1] or not A.size or not np.isfinite(A).all():
        raise ValueError("nonempty finite square matrix required")
    if not np.array_equal(A, A.T):
        raise ValueError("real symmetric matrix required")
    if not np.isfinite(tol) or tol <= 0 or not isinstance(maxiter, int) or maxiter < 0:
        raise ValueError("positive tolerance and nonnegative integer maxiter required")
    original = A.copy(); Z = np.eye(A.shape[0])
    scale = max(1.0, np.linalg.norm(A))
    if not np.isfinite(scale):
        raise ArithmeticError("matrix norm outside floating-point range")
    threshold = tol*scale
    residual_threshold = 10*threshold
    if not np.isfinite(threshold) or not np.isfinite(residual_threshold):
        raise ArithmeticError("stopping threshold outside floating-point range")
    for k in range(maxiter+1):
        off = A-np.diag(np.diag(A))
        if np.linalg.norm(off) <= threshold:
            values = np.diag(A).copy()
            if np.linalg.norm(original@Z-Z*values) <= residual_threshold:
                return values, Z, k
            raise RuntimeError("eigenpair residual failed")
        if k == maxiter:
            break
        Q, R = np.linalg.qr(A)
        A = R@Q; Z = Z@Q
        if not np.isfinite(A).all() or not np.isfinite(Z).all():
            raise ArithmeticError("nonfinite QR iterate")
    raise RuntimeError("unshifted QR did not converge within iteration limit")


values, vectors, steps = symmetric_qr_demo([[2, 1], [1, 2]])
print(values)  # [3. 1.] up to rounding
print(np.linalg.norm(np.array([[2, 1], [1, 2]])@vectors-vectors*values))