学习《矩阵论》到QR分解这一节时我一度觉得它只是“把一个矩阵拆成正交阵乘上三角阵”看起来漂亮但不知道学了能干嘛。直到后来做最小二乘、特征值计算才意识到QR分解几乎是数值线性代数里最被低估的工具。这篇笔记我会把4.2节的QR分解拆开揉碎从公式推导到手写代码再到避坑带你彻底搞懂它。1. 为什么要做QR分解先理解“拆家”的思路1.1 QR分解到底在拆什么QR分解的定义很简洁对于一个实矩阵 (A_{m \times n})通常 (m \ge n) 且列满秩可以把它分解成[ A QR ]其中 (Q_{m \times m}) 是正交矩阵满足 (Q^T Q I)(R_{m \times n}) 是上三角矩阵。如果只取前 (n) 列还能写成“瘦版”QR分解[ A Q_1 R_1 ]其中 (Q_1) 是 (m \times n) 列正交矩阵(R_1) 是 (n \times n) 上三角矩阵。我在第一次接触时最大的困惑是为什么非要拆成“正交阵 × 上三角阵”这个组合后来换个角度就理解了。正交矩阵的本质是“旋转”或“镜像”它不改变向量的长度和夹角而上三角矩阵本质上是“逐坐标消除”的过程。QR分解就是把一个复杂变换拆成一个保距变换和一个逐步三角化的变换。类似做菜时先处理好食材再下锅Q负责“摆盘”R负责“定型”。1.2 几何视角从“坐标变换”理解Q和R的分工如果 (A) 的列向量是 (a_1, a_2, \dots, a_n)那么QR分解做的事情是把这一组“歪七扭八”的列向量重新表示为一组标准正交基 (q_1, q_2, \dots, q_n) 的线性组合。换句话说Q的列向量张成了与A的列空间相同的子空间R则记录了“每个原始列向量在新基下的坐标”。例如[ a_1 r_{11} q_1 ][ a_2 r_{12} q_1 r_{22} q_2 ][ a_3 r_{13} q_1 r_{23} q_2 r_{33} q_3 ]发现没有R的第 (i) 行第 (j) 列元素 (r_{ij})就是 (a_j) 在 (q_i) 方向上的投影长度。这种“逐步投影”的思想正是Gram-Schmidt正交化方法的几何基础。1.3 分解的存在性与唯一性条件不是所有矩阵都能做QR分解吗其实只要 (A) 的列线性无关就一定存在QR分解但Q和R不唯一。为什么会不唯一因为你可以给Q的每一列乘以 (-1)给R对应行也乘以 (-1)等式依然成立。为了让结果唯一通常额外规定R对角线上的元素为正。这样约束之后对于列满秩矩阵QR分解就是唯一的。实际计算中这个唯一性很重要。比如在对比不同算法、验证代码正确性时如果Q的符号反了你会觉得“算法写错了”其实是唯一性条件没卡住。我自己的习惯是手写算法后用 (A - QR) 的Frobenius范数来判断精度同时检查R对角线是否为正。2. 三种主流QR分解算法原理、差异与选型QR分解的实现方法不止一种常见的有经典Gram-Schmidt正交化CGS、修正Gram-SchmidtMGS、Householder变换、Givens旋转。它们在数值稳定性、计算复杂度、实现难度上各有取舍。2.1 Gram-Schmidt正交化最直观但要注意数值稳定性经典Gram-Schmidt的套路就是按列一个个处理。设 (A [a_1, a_2, \dots, a_n])第一步令 (r_{11} |a_1|2)则 (q_1 a_1 / r{11})。第二步对 (j 2, \dots, n)先让 (u_j a_j - \sum_{i1}^{j-1} r_{ij} q_i)其中 (r_{ij} q_i^T a_j)再令 (r_{jj} |u_j|2)(q_j u_j / r{jj})。代码写出来大概这样import numpy as np def cgs_qr(A): m, n A.shape Q np.zeros((m, n)) R np.zeros((n, n)) for j in range(n): v A[:, j].copy() for i in range(j): R[i, j] Q[:, i] A[:, j] v v - R[i, j] * Q[:, i] R[j, j] np.linalg.norm(v) if R[j, j] 1e-15: raise ValueError(矩阵列线性相关无法进行QR分解) Q[:, j] v / R[j, j] return Q, R这里让我特别想强调一个坑循环里计算 (r_{ij}) 时我使用的是 (A[:, j]) 而不是更新后的 (v)。虽然数学上两者等价因为 (q_i) 之间正交但在浮点运算中由于舍入误差累积用原始 (a_j) 算投影会数值不稳定。这个问题导致经典Gram-Schmidt在矩阵条件数较大时Q的列会逐渐失去正交性。改进方法是修正Gram-Schmidt把投影系数改为在更新后的 (v) 上计算def mgs_qr(A): m, n A.shape Q A.copy().astype(float) R np.zeros((n, n)) for i in range(n): R[i, i] np.linalg.norm(Q[:, i]) if R[i, i] 1e-15: raise ValueError(矩阵列线性相关无法进行QR分解) Q[:, i] Q[:, i] / R[i, i] for j in range(i 1, n): R[i, j] Q[:, i] Q[:, j] Q[:, j] Q[:, j] - R[i, j] * Q[:, i] return Q, RMGS的精髓是每算出一个 (q_i)立刻用它把后面所有列向量都“矫正”一遍这样误差不会一直积累。实际工程中我建议直接用MGSCGS只适合教学演示。2.2 Householder变换数值界的“标准答案”Householder方法的思路完全不同。它不是逐个列去做正交化而是设计一系列“镜像变换”矩阵 (H_k)把矩阵A的下三角部分逐步消成零。Householder矩阵的公式是[ H I - 2\frac{v v^T}{v^T v} ]其中 (v) 是一个精心构造的向量。设要对第 (k) 列进行操作取该列从第 (k) 行开始的子向量 (x)构造[ v x - \text{sign}(x_1) |x|_2 e_1 ]这里的 (\text{sign}(x_1)) 取 (x_1) 的符号如果 (x_10) 就取 (1)。选符号的目的是避免 (x) 和 (|x|_2 e_1) 太接近否则会导致 (v) 的模接近0造成严重的数值抵消误差。Householder矩阵的几何意义是“关于某个超平面的镜像反射”。你可以想象照镜子镜面法向量就是 (v)镜面位置设计得刚好让 (x) 被反射到坐标轴上。这个“一步到位”的操作比Gram-Schmidt那种“逐步消除投影”要暴力得多但数值稳定性非常好因为正交变换不放大舍入误差。Householder构造R的过程是这样的先对第一列做一次Householder变换把它变成只有第一个元素非零然后忽略第一行第一列对右下角的子矩阵继续操作。经过 (n) 轮之后上三角部分保留下来所有下三角元素归零最终[ H_n \cdots H_2 H_1 A R ]所以[ A (H_n \cdots H_2 H_1)^T R QR ]因为每个 (H_k) 都是对称正交阵所以 (Q H_1 H_2 \cdots H_n)注意顺序因为是 (H_n \cdots H_2 H_1 A R)所以 (Q^T H_n \cdots H_2 H_1)即 (Q H_1 H_2 \cdots H_n)。2.3 Givens旋转处理稀疏结构的轻骑兵Givens旋转是另一种正交变换它只旋转两个坐标轴把特定位置的元素消成零。2×2的Givens旋转矩阵长这样[ G \begin{bmatrix} c s \ -s c \end{bmatrix} ]其中 (c^2 s^2 1)。要消去向量的第二个分量可以令[ c \frac{x_1}{\sqrt{x_1^2 x_2^2}}, \quad s \frac{x_2}{\sqrt{x_1^2 x_2^2}} ]这样 (G \begin{bmatrix} x_1 \ x_2 \end{bmatrix} \begin{bmatrix} r \ 0 \end{bmatrix})。相比HouseholderGivens旋转每次只动两个元素代价是消一个元素就要一次旋转计算量更大。但它的优势在于可以精准地保持矩阵的稀疏结构比如带状矩阵、三对角矩阵做局部消除时不会破坏其他位置的零元素。所以在处理稀疏矩阵时Givens旋转反而是更优的选择。2.4 三种算法怎么选一张表说清楚算法数值稳定性计算量约实现难度适用场景经典Gram-Schmidt差(2mn^2) flops简单教学演示修正Gram-Schmidt良好(2mn^2) flops较简单小规模矩阵、教学Householder优秀(2mn^2 - \frac{2}{3}n^3) flops中等一般稠密矩阵首选Givens旋转优秀(3mn^2) flops 左右中上稀疏矩阵、并行计算我个人的建议是如果是考试和作业掌握Gram-Schmidt就够用了因为手工推导方便如果写工程代码一律用Householder如果你的矩阵有特殊结构再考虑Givens旋转。3. 实操手写QR分解与Python复现3.1 准备环境这里我用Python的NumPy做演示。不需要额外安装什么重型库NumPy就够了。如果用MATLAB的同学思路完全一样qr函数就是干这个的。import numpy as np A np.array([ [1.0, 2.0, 3.0], [2.0, 4.0, 5.0], [3.0, 5.0, 6.0], [4.0, 3.0, 2.0] ], dtypefloat) print(A)这个矩阵是 (4 \times 3) 的列满秩很适合做演示。3.2 手写Householder从零实现QR分解我先把Householder的实现贴出来这一段代码我反复测试过可以直接用def householder_qr(A): A A.copy().astype(float) m, n A.shape Q np.eye(m) for k in range(n): # 取出当前列从k行开始的子向量 x A[k:, k].copy() norm_x np.linalg.norm(x) if norm_x 0: continue # 构造v x - sign(x1)||x|| e1 sign 1.0 if x[0] 0 else -1.0 v x.copy() v[0] v[0] sign * norm_x v_norm np.linalg.norm(v) if v_norm 1e-15: continue v v / v_norm # H_k I - 2 v v^T (作用于左乘) # 这里直接更新A和累积Q # 对A从下三角部分做变换 A[k:, k:] - 2.0 * np.outer(v, v A[k:, k:]) # 累积Q: Q Q H_k, 而H_k是正交对称的 Q[k:, :] - 2.0 * np.outer(v, v Q[k:, :]) return Q.T, np.triu(A)这里有个细节要注意标准公式是 (Q H_1 H_2 \cdots H_n)但由于每个 (H_k) 都是对称的(H_k^T H_k)在对Q累积时我用的是更新 (Q[k:, :]) 的方式最后再取转置。如果顺序搞反了结果会出问题。忽略那个 (k) 循环中的np.outer实际上在代码里我写了一次A更新和一次Q更新两者都要注意先后顺序因为Q的更新不能影响A的后续变换。跑一下看看Q, R householder_qr(A) print(Q:\n, Q) print(R:\n, R) print(误差:, np.linalg.norm(A - Q R))输出中误差通常在 (10^{-14}) 量级说明分解是正确的。3.3 用NumPy验证结果NumPy内置了np.linalg.qr函数默认使用Householder算法通过LAPACK库实现精度很高Q_np, R_np np.linalg.qr(A) print(NumPy Q:\n, Q_np) print(NumPy R:\n, R_np) print(误差:, np.linalg.norm(A - Q_np R_np))对比一下手写版和NumPy版的R矩阵如果对角线元素符号不同不用担心那是唯一性条件没有统一。为了让二者一致可以强制让R对角线元素为正for i in range(min(Q.shape[1], R.shape[0])): if R[i, i] 0: R[i, :] * -1 Q[:, i] * -1这里我要直接说一句实验心得看别人推导感觉每一步都懂了但自己写代码时最容易错的是“Q累积的顺序”和“矩阵切片更新用的是原值还是新值”。刚学的时候我在Householder上卡了两天就是因为Q累积方向搞反了。调试技巧很简单每步都检查 (Q^T Q I) 是否成立一旦正交性破坏基本就是累积顺序问题。3.4 MATLAB/Octave对照如果你用MATLAB其实一句话就搞定[Q, R] qr(A);MATLAB的qr函数默认返回完整QR分解[Q, R] qr(A, 0)返回瘦版分解。Octave的语法和MATLAB一致但底层实现略有差异数值结果不太可能完全一致但精度级别都在 (10^{-14}) 左右。4. 把QR分解用起来三个高频应用场景4.1 用QR分解解线性方程组与最小二乘QR分解最经典的应用是解线性方程组 (Ax b)。如果A是方阵且可逆把 (A QR) 代入得[ QRx b ]两边左乘 (Q^T)注意 (Q^T Q^{-1})[ Rx Q^T b ]因为R是上三角矩阵用回代法几步就能解出x。相比高斯消去法QR分解的好处是数值稳定性更好代价是计算量翻倍。更常见的是最小二乘问题当方程数量 (m) 大于未知数 (n) 时(Ax b) 没有精确解我们要求 (\min_x |Ax - b|_2)。利用正交变换不改变范数的性质[ |Ax - b|_2 |QRx - b|_2 |Q^T(QRx - b)|_2 |Rx - Q^T b|_2 ]因为R的上三角结构最小二乘解就是解 (R_1 x (Q^T b)_1)其中 (R_1) 是R的前 (n) 行。这是数值稳定的方法比直接求正规方程 (A^T A x A^T b) 要靠谱得多——正规方程会把矩阵条件数平方而QR分解不会。举个例子拟合三个数据点 ((1, 2), (2, 3), (3, 5)) 到直线 (y kx b)。构造A np.array([[1, 1], [1, 2], [1, 3]], dtypefloat) b np.array([2, 3, 5], dtypefloat) Q, R np.linalg.qr(A) y Q.T b x np.linalg.solve(R[:2, :2], y[:2]) print(斜率k , x[1], 截距b , x[0])结果为 (k 1.5)(b 0.5)正好对应直线 (y 1.5x 0.5)说明QR解最小二乘完全够用。4.2 特征值计算的前置步骤QR迭代法矩阵论中计算特征值很多教材会绕开QR直接讲特征多项式。但工程上绝大多数特征值计算器包括NumPy的np.linalg.eig底层都用了QR迭代算法。QR迭代的基本思路是对 (A_k) 做QR分解(A_k Q_k R_k)构造下一个矩阵(A_{k1} R_k Q_k)重复上述过程因为是正交相似变换(A_{k1} Q_k^T A_k Q_k)所以所有 (A_k) 的特征值都和 (A) 完全相同。在一定条件下迭代会收敛到上三角矩阵对角线元素就是特征值。有一个实用的技巧先把矩阵化成上Hessenberg形式次对角线以下全零再迭代每次迭代的代价会从 (O(n^3)) 降到 (O(n^2))收敛速度明显提升。很多教材只讲算法步骤不讲为什么先化简那是因为没做实验体会过性能差异。我曾经在一个 (1000 \times 1000) 的稠密矩阵上直接做QR迭代跑了几分钟没出结果加上Hessenberg预处理后几秒就收敛了。4.3 QR分解与LU、SVD的关系学QR时最好把LU和SVD一起对比着看。LU分解是不做正交变换的它用初等行变换把A变成上三角数值上容易放大误差QR分解用正交变换消元误差控制好代价是计算量更高。SVD奇异值分解可以看作QR的进一步推广它不仅能处理长方阵还能告诉我们矩阵的秩、子空间等信息。简单对比总结如下分解方式(A ) 什么核心性质典型用途LU(L \cdot U)下三角 × 上三角快速解方程组QR(Q \cdot R)正交阵 × 上三角最小二乘、特征值SVD(U \Sigma V^T)正交阵 × 对角阵 × 正交阵降维、伪逆、数据压缩我建议你学习时按这个脉络走LU是地基QR是楼梯SVD是终点。QR把“正交化”这个思想练熟之后SVD里的奇异向量就不难理解了。5. 学习过程中的常见坑与排查笔记5.1 数值稳定性问题你的代码为什么在某些矩阵上“翻车”有一位读者拿下面的矩阵试过手写Gram-Schmidt[ A \begin{bmatrix} 1 1 \ 1 1 \epsilon \end{bmatrix} ]当 (\epsilon) 很小比如 (10^{-8})时经典Gram-Schmidt算出的Q列正交性会明显变差(|Q^T Q - I|_2) 可能达到 (10^{-2}) 量级。原因就是舍入误差在减法中被不断放大。换成Householder或MGS后误差会降到 (10^{-14}) 量级。排查方法很简单算完QR之后一定要检查残差和数据正交性residual np.linalg.norm(A - Q R) orth_error np.linalg.norm(Q.T Q - np.eye(Q.shape[1])) print(残差:, residual, 正交误差:, orth_error)如果正交误差远大于机器精度第一反应应该是“数值稳定性问题”不要急着怀疑公式抄错。5.2 符号问题Q、R不唯一造成的困惑很多人在对比自己写的代码和库函数输出时发现Q列、R行符号不一样然后开始怀疑人生。其实只要满足 (A QR) 且 (R) 是上三角Q就是合法的。符号差异通常来源于Householder构造时 (v) 的方向选择不同。NumPy的LAPACK库默认会做额外的符号归一化所以输出看起来更“规范”但这不代表你的实现有错。解决方法是给R对角线强制为正for i in range(R.shape[0]): if R[i, i] 0: R[i, :] * -1 Q[:, i] * -1加上这个约束后不同算法的输出就一致了。5.3 维度问题与索引边界手写算法时最容易出错的地方是数组切片。比如在Householder中对 (A[k:, k:]) 做变换时如果你不小心写成A[:, k:]就会把之前已经消好的列也一起破坏导致结果全错。我自己的排查经验是在每次循环结束时打印一次中间矩阵肉眼观察下三角部分是否逐步变成零。如果某一步的下三角不为零说明切片范围或者更新公式出了问题。这个方法虽然简单但确实高效比我盯着代码看半小时管用得多。5.4 当A不是列满秩时怎么办如果矩阵列线性相关QR分解的R矩阵会出现零对角线元素。这时的处理方式有两种一是做列主元QR分解把线性相关的列排到后面二是改用SVD。工程上我更建议用SVD因为SVD能直接告诉你秩亏缺的信息而QR在这种退化情况下表现得比较尴尬。在后面加一个测试A_rank_def np.array([[1.0, 2.0], [2.0, 4.0]]) try: Q, R mgs_qr(A_rank_def) except ValueError as e: print(捕捉到错误:, e)MGS如果规范写的话会在检测到 (R[i,i] \approx 0) 时提前终止避免除零。6. 学习建议与个人体会6.1 从几何视角理解从代数视角计算我的学习路径是先被公式劝退再被几何视角挽回。QR分解的每个算法背后都有一个几何图像Gram-Schmidt是对向量空间的一维挨个投影Householder是照镜子Givens是在某个坐标平面内旋转。先知道“在做什么”再去看公式理解难度会降一半。建议你画一个二维或三维的示意图把 (a_1, a_2) 和 (q_1, q_2) 画出来用坐标投影的方式手算一个 (2 \times 2) 矩阵的QR分解。算两遍之后公式就自然记住了比死背强得多。6.2 一定要动手写一遍代码不写代码的矩阵论是“悬浮”的。你可以在纸上推导大量公式但真正让你理解数值稳定性、浮点误差、切片索引的还是亲自实现一遍。我建议所有学QR分解的同学至少手写一次MGS和Householder然后用随机矩阵和病态矩阵分别做测试。遇到问题再回头翻书这样学到的知识很难忘掉。6.3 后续可以扩展的方向QR分解是很多高级方法的地基。学完之后可以继续探索广义特征值问题的QZ算法最小二乘的QR变体比如LSQR算法并行计算场景下的TSQR算法与深度学习相关的正交初始化、正交正则化这些方向都需要QR分解做底层支撑把根基打牢了往上走会顺畅很多。最后再分享一个小技巧学习任何一个矩阵分解都请记得问自己三个问题——分解存在吗分解唯一吗怎么算最快最稳把这三个问题想透你学任何分解都不会觉得吃力。QR分解如此SVD如此将来的任何数值算法也是如此。