学《矩阵论》学到4.2节QR分解的时候我突然意识到前面学的东西都可以换个角度看。LU分解是高斯消去法的矩阵写法特征值分解是找不变方向而QR分解本质上是给矩阵做了一次系统性正交化——把一组可能彼此“纠缠”的列向量拆成一个正交矩阵Q和一个上三角矩阵R。这一节内容是后续最小二乘、特征值迭代、数值线性代数许多算法的地基。这篇文章是我学习这一节时整理的完整笔记包含定义、三种构造方法、稳定性对比和几个实际应用适合正在啃矩阵论或者数值分析教材的读者参考。1. 从“解方程组”到“找正交基”QR分解在做什么先回顾一下前面学过的内容。LU分解解决的是“方阵、主元不为零”情形下的线性方程组本质上就是把高斯消去法写成了矩阵乘积的形式。但现实里我遇到更多的问题是两类一类是方程个数多于未知数的超定方程组另一类是矩阵本身接近病态、主元小得离谱的情况。这时候LU不一定不能用但至少不太好用。QR分解的处理思路完全不同。它不把精力花在“消元”上而是先把矩阵A重新“整理”一遍。具体来说对于列满秩的m×n实矩阵Am≥n存在m×m正交矩阵Q和m×n上三角矩阵R使得A Q R如果只需要n列也可以取Q的前n列得到所谓“瘦QR分解”reduced QR此时Q是m×n列正交矩阵R是n×n上三角矩阵。正交矩阵的最大好处是逆矩阵等于转置这给后续计算省了无数事。用一个小例子立刻能看出它在做什么。取A [[1, 1], [2, 1], [2, 0]]这个矩阵的第一列和第二列在三维空间里指向两个不同的方向。对它做QR分解可以得到Q [[1/3, 2/3], [2/3, 1/3], [2/3, -2/3]]R [[3, 1], [0, 1]]验证一下Q的第一列乘3等于(1,2,2)正好是A的第一列Q的第二列乘1再加上Q的第一列等于(1,1,0)正好是A的第二列。看到这里几何意义就很清楚了Q的两列是把A的两列张成的那个平面重新换了一套“标准正交坐标轴”R记录的则是原来两个向量在这套新坐标轴下的坐标。换句话说QR分解不是在解某个具体方程而是在改造矩阵本身——把一组可能方向很接近、数值上很“挤”的列向量变成了彼此垂直、长度规范的标准正交基以及一个描述“原列向量如何由这些基组合出来”的上三角系数矩阵。为什么这个改造有用因为方程Axb在分解之后变成QRxb两边左乘Q^T得到Rx Q^T bR是上三角矩阵解这个方程只需要从最后一行开始逐个回代复杂度只有O(n²)。而QR分解本身的复杂度对m×n矩阵大约是O(mn²)量级。也就是说我们把一个“难解”的一般线性系统变成了一次性付出稍高代价的分解加一个极便宜的回代。很多工程计算里同一个A会被反复用来解不同的b这时候把A分解一次、存好Q和R之后每次求解都只需要O(n²)的回代整体收益非常可观。还有一点容易被忽略QR分解并不要求A是方阵。mn的瘦长矩阵、列满秩的矩阵统统可以做。这正是它在线性回归、参数估计这类超定问题里大显身手的原因。2. 三种构造QR的路径Gram-Schmidt、Householder与Givens矩阵论教材通常用Gram-Schmidt正交化来证明QR分解的存在性但实际工程里大量使用的是Householder变换和Givens旋转。这三条路径各有各的适用场景我一个个说。2.1 Gram-Schmidt正交化最直观但数值上最脆弱Gram-Schmidt的想法非常朴素从第一列开始逐列做投影把当前列中与前面所有已正交化方向平行的分量减掉剩下的部分就是新的正交方向。还是用上面那个例子。A的第一列是a₁(1,2,2)长度为3于是q₁ (1/3, 2/3, 2/3)R的第一个对角元r₁₁3。接着处理第二列a₂(1,1,0)先算它在q₁方向上的投影长度r₁₂ q₁ᵀa₂ 1/3×1 2/3×1 2/3×0 1从a₂里减去这个投影得到v₂ a₂ - r₁₂q₁ (1,1,0) - (1/3,2/3,2/3) (2/3,1/3,-2/3)v₂的长度恰好是1所以q₂ (2/3,1/3,-2/3)r₂₂ 1整理一下Q的两列就是q₁和q₂R就是[[3,1],[0,1]]。整个过程很顺每一步都有明确的几何画面。把Gram-Schmidt写成矩阵语言就是经典Gram-Schmidt算法for j 1..n: v_j a_j for i 1..j-1: r_ij q_i^T a_j v_j v_j - r_ij q_i r_jj ||v_j|| q_j v_j / r_jj但这里有个隐藏很深的坑经典Gram-Schmidt在浮点运算下会逐渐丢失正交性。原因在于当a_j与前面某个q_i方向很接近时投影量r_ij很大做减法“a_j - r_ij q_i”时会出现两个大数相减有效数字被抵消掉误差就混进去了。这个误差会被后面每一步继续放大。改进办法是Modified Gram-Schmidt修正Gram-Schmidt它在算出q_i之后立刻把剩余所有列都减去对q_i的投影而不是等到处理第j列时才回头使用原始a_jfor i 1..n: r_ii ||v_i|| q_i v_i / r_ii for j i1..n: r_ij q_i^T v_j v_j v_j - r_ij q_i两者理论上是等价的但后者的数值稳定性明显更好因为每一步处理时v_j已经被前面的q_1,...,q_{i-1}修正过减去的投影不会累积那么多舍入误差。不过即使修正过Gram-Schmidt系列在强病态情况下依然不如Householder靠谱。作为学习路径我很推荐它但真要写生产代码我会优先选下一种。2.2 Householder变换反射一下就消掉一整列Householder变换的出发点是另一个几何直觉给定一个向量a能不能找一面“镜子”让a反射之后只保留第一个分量其余分量全部变成零这面镜子就是法向量为w的超平面对应的变换矩阵是H I - 2wwᵀ/(wᵀw)可以验证H是对称阵而且HᵀHI所以H是正交矩阵。它作用在任意向量x上效果是x关于法向量w的平面做镜像反射。关键问题是怎么选w。如果我想让Ha变成αe₁e₁是第一个标准基向量取w a - αe₁其中α ±||a||。为了让数值计算更稳取α与a的第一个分量同号避免w的两个分量接近抵消。还是用刚才的例子。a(1,2,2)||a||3a₁10取α3于是w a - 3e₁ (-2, 2, 2)算一下H₁ I - 2wwᵀ/12作用在a上确实得到3e₁。这个例子还出现了一个小巧合A的第二列v(1,1,0)而wᵀv -220 0所以H₁作用在第二列上完全不改变它。于是H₁A [[3, 1], [0, 1], [0, 0]]接下来对右下角2×2块里的向量(1,1)再做一次Householder变换H₂把它变成(√2,0)。重复同样的构造取w₂(1,1)-√2e₁(1-√2,1)H₂作用后得到R H₂H₁A [[3, 1], [0, √2], [0, 0]]因为每次乘的H都是正交矩阵它们的乘积也是正交矩阵于是A Q R其中Q H₁H₂H₂按3×3嵌入理解作用于前两行。Householder变换的整个流程是第一步用一个H₁消掉第一列主对角线以下的全部元素第二步对右下角子块再用一个H₂消掉第二列主对角线以下的元素如此反复。每步只处理一个右下角子块思路非常机械适合编程实现。它比Gram-Schmidt稳是因为每次变换本质上是正交反射不涉及“大数相减再重新归一化”这种容易损失精度的操作。实际LAPACK里的dgeqrf等标准QR分解实现底层基本就是Householder变换。2.3 Givens旋转定点消零稀疏矩阵的利器Givens旋转的思路更“局部”。二维平面里旋转矩阵G [[c, -s], [s, c]]可以把向量(x,y)旋转到x轴上其中c x/√(x²y²)s y/√(x²y²)。这个性质看着简单但好处是它只影响两个分量。对于n阶矩阵Givens旋转就是在一个坐标平面上做旋转其他坐标保持不变。具体来说G(i,j,θ)是把单位阵的第i行第i列、第i行第j列、第j行第i列、第j行第j列四个位置替换成相应的c和s。它同样是正交矩阵。做QR时Givens旋转的用法是从上到下、从左到右逐个消去主对角线以下的元素。比如对第一列先用G(1,2)消掉第2行元素再用G(1,3)消掉第3行元素逐一把第一列主对角线以下全部清零然后转向第二列用G(2,3)、G(2,4)等继续消。每步只改变两行所以对于稀疏矩阵来说这种局部性非常宝贵——你不会像Householder那样一个反射就把本来稀疏的矩阵变稠密。代价是计算量更大。Householder一次反射能清掉一整列下面的所有元素Givens一次只消掉一个元素而且每步都要算一次平方根。对那些本身就是上Hessenberg形式主对角线下一层之外全是零的矩阵Givens旋转的数量和每个旋转的flops都很少所以在QR迭代求特征值时Hessenberg矩阵配Givens几乎是标准搭配。三种方法放在一起对比方法主要思想计算量m×n列满秩数值稳定性典型场景经典Gram-Schmidt逐列投影约2mn²较差理论证明、教学手算修正Gram-Schmidt逐列投影并即时修正约2mn²中等小规模实现、Arnoldi迭代Householder镜像反射约2mn²−(2/3)n³优稠密矩阵通用首选Givens平面旋转约3mn²量级优稀疏矩阵、Hessenberg矩阵3. 数值稳定性这堂课为什么教材例题总是“恰好”这么干净矩阵论课上用Gram-Schmidt讲QR分解作业里也都是那种数字精巧的小矩阵算出来的Q列向量几乎是完美的正交基。但当我第一次拿经典Gram-Schmidt去解一个稍微病态一点的实际数据矩阵时发现QᵀQ离单位阵差得离谱——这堂课是数值稳定性给我上的。问题的根源在于浮点数只能保留有限位有效数字。假设某一列a_j几乎落在前面某些q_i张成的子空间里那么投影量r_ij会非常大减法a_j - r_ij q_i的结果是一个很小的向量。这个“小向量”的有效数字在大数相减中已经损失了相当一部分再拿去归一化得到的方向就与理想情况有明显偏差。更糟的是经典GS是逐列处理的前一步的误差会污染后续所有列导致Q的正交性随列数增加而加速恶化。有人可能会问条件数差的矩阵用精确算术算出来不是一样的吗对但现实中我们永远在用浮点算术。数值稳定性关注的正是舍入误差会不会被算法本身放大。Householder和Givens之所以稳定是因为它们每一步都是正交变换几何上相当于旋转/反射整个空间不会像Gram-Schmidt那样把“当前向量减去投影”这一步放大误差。我实际测试过一个3×3例子A的第一列和第二列故意取成几乎线性相关即A [[1, 1], [1ε, 1], [1, 1]]经典GS算出的q₁和q₂内积可能达到10⁻²级别QᵀQ的非对角元完全不可接受。换成修正GS能改善一个数量级但仍然不够理想。用Householder变换QᵀQ的非对角元能压到机器精度附近。这里不是说要完全抛弃Gram-Schmidt。修正Gram-Schmidt在Krylov子空间方法比如Arnoldi迭代里依然是核心工具因为那些场景下我们本来就要逐列生成正交基无法一次性用Householder整体分解。关键在于用对地方手算推导用小矩阵没问题但写代码做数值计算默认优先Householder。另一个容易踩的坑是“重新正交化”。有的老代码会在Gram-Schmidt后检查正交性发现不够再做一遍。这是应急手段不是好设计。重新正交化能改善结果但会掩盖算法本身的不稳定而且没法预测效果。与其这样不如一开始就用数值稳定的分解方法。4. 应用场景的展开最小二乘、特征值迭代与矩阵求逆QR分解在《矩阵论》里只占一小节但它在数值计算的各个角落反复出现。这里挑三个高频场景展开。4.1 最小二乘问题比“正规方程”稳定得多超定方程组Ax≈b在工程里太常见了。随便举一个曲线拟合的例子有三个数据点(0,1),(1,2),(2,2)想用直线y a bx去拟合。设计矩阵A [[1, 0], [1, 1], [1, 2]]b (1, 2, 2)这是3个方程、2个未知数通常不存在严格解。教科书上最经典的做法是解正规方程AᵀAx Aᵀb但正规方程有个致命缺陷AᵀA的条件数等于cond(A)²矩阵越病态条件数被放大得越厉害舍入误差对解的扰动呈平方级上升。用QR分解做最小二乘就不一样了。因为Q是正交矩阵它不改变向量的长度||Ax - b||² ||QRx - b||² ||Rx - Qᵀb||²这里Q和R取瘦QR形式R是n×n上三角。R的前n行是上三角矩阵后面的行全是零所以上式里能找到最小二乘解的部分只有前n行后m-n行对应的残差是固定的。于是问题退化成解一个n×n上三角方程Rx (Qᵀb)的前n行直接用上面例子算一下。对A做瘦QR可以得到R [[√3, √3], [0, √2]]Qᵀb (5/√3, 1/√2)只取两行解上三角方程√3a √3b 5/√3 √2b 1/√2得到b 1/2a 7/6。和正规方程解出来一模一样但整个过程不需要计算AᵀA也不会把条件数平方放大。用Python写也就几行import numpy as np A np.array([[1, 0], [1, 1], [1, 2]], dtypefloat) b np.array([1, 2, 2], dtypefloat) Q, R np.linalg.qr(A, modereduced) x np.linalg.solve(R, Q.T b) print(x) # [1.16666667 0.5]很多数值线性代数库里的最小二乘求解器底层用的就是带列主元的QR分解原因就在这里。4.2 特征值迭代QR分解最出人意料的应用第一次看到QR迭代算法时我是有点震惊的——原来求特征值还能用QR分解自己迭代出来。算法本身极其简单A_1 A for k 1, 2, ...: A_k Q_k R_k A_{k1} R_k Q_k每次迭代先把当前矩阵分解成Q_k R_k然后把两个因子交换顺序乘回去。为什么这个看似奇怪的迭代会收敛到特征值因为每次迭代都是相似变换A_{k1} R_k Q_k Q_kᵀ(Q_k R_k)Q_k Q_kᵀ A_k Q_k所以所有A_k的特征值都和原始矩阵A完全一样。而这个过程会让矩阵逐渐趋近于上三角形式实矩阵特征值不完全相同时对角元就是特征值遇到共轭复特征值则会留下2×2对角块。举一个2×2小例子A [[2,1],[1,2]]。对它做一次QR分解得到Q [[2/√5, -1/√5], [1/√5, 2/√5]]R [[√5, 4/√5], [0, 3/√5]]然后算RQA₁ RQ [[14/5, 3/5], [3/5, 6/5]]直观地看原来A的非对角元是1迭代一次后变成0.6对角元向3和1靠拢这个矩阵的特征值正好是3和1。继续迭代下去非对角元会越来越小最终逼近对角阵diag(3,1)。真实工业界当然不会直接对稠密矩阵做这个迭代因为每一步QR都很贵。标准做法是先用Householder变换把A化为上Hessenberg形式主对角线下方只有一层非零元再做带位移的QR迭代比如Wilkinson位移每一步QR再配合Givens旋转高效完成。但所有这些手段的地基还是QR分解本身。学4.2节时知道这层背景后面看特征值数值算法会非常受用。4.3 矩阵求逆利用R上三角结构省算力对一个可逆方阵A做QR分解AQR之后求逆变得非常直接A⁻¹ R⁻¹ Qᵀ因为Q是正交矩阵Qᵀ就是Q⁻¹R是上三角矩阵求逆只需要从最后一行开始逐个回代整个过程的计算量和数值稳定性都比老老实实的高斯消去法更可控。这里的核心思想还是那句先把A化简成容易处理的结构再利用这种结构做后续计算。实际工程里矩阵条件数偏高又必须反反复复求解Axb时我往往会把QR分解结果缓存下来之后每次只需要算Qᵀb和一次回代。相比每次重新做消去这样能节约大量重复计算而且精度有保障。5. 学习QR分解时容易踩的坑这一节聊聊我学习过程中碰到的几个具体问题都是教材不会特意提醒、但自己动手时很快会撞上的地方。第一个坑是把Gram-Schmidt和QR分解当成两个不相干的话题。它们其实是同一件事的两种说法Gram-Schmidt过程天然就产生了一组正交基和对应的系数把这些系数按列排成矩阵写成AQR就是QR分解。教材把它们放在不同章节讲考试时却默认你理解这一点刚学的时候容易晕。第二个坑是忽略分解不唯一。QR分解不是绝对的。如果Q的某一列取反同时R的对应行取反乘积QR仍然等于A。所以手算出来的Q和调库算出来的Q逐元素对不上是非常正常的。判断对错的标准是看QR是否等于A而不是比较Q本身。同样道理教材里有时说“列满秩时QR分解唯一”那是有附加条件的——通常要求R的对角线元素为正。去掉这个限制符号自由度就回来了。第三个坑是把“瘦QR”和“完整QR”搞混。完整QR里Q是m×m方阵R是m×n矩阵下半部分是零行瘦QR里Q只有n列R是n×n方阵。很多教材混用这两种形式做题时要注意题目到底要哪种。在NumPy里np.linalg.qr(A)默认给完整形式np.linalg.qr(A, modereduced)给瘦形式用错维度会直接报错或者得到一堆奇怪的零。第四个坑是手算和程序结果对不上时先怀疑R的符号再怀疑自己的手算过程。我最初用Householder手算一个3×3矩阵的QR和MATLAB结果差在第二列全是相反数折腾半天才意识到这是分解的符号自由度不是算错。第五个坑是在病态问题上继续用正规方程。很多人因为只会解线性最小二乘的正规方程遇到条件数大的设计矩阵解出来数值误差大到没法用。换个QR分解条件数不会像AᵀA那样被平方放大结果会好非常多。另外一个容易被忽略的细节是复矩阵。矩阵论里经常遇到Hermite矩阵这时QR分解里的转置要换成共轭转置Q满足QQI称为酉矩阵。Householder变换的公式变成H I - 2vv/(v*v)Givens旋转的构造也要带上共轭。公式形式几乎一样但漏掉共轭符号会导致完全错误的结果。我的学习建议是拿一个小矩阵先用Gram-Schmidt手算一遍完整的QR再用Householder手算一遍最后用np.linalg.qr验证一次。三次得到的Q可能不完全相同但QR的乘积都等于原矩阵。亲手走完这个过程你对分解不唯一性、瘦QR和完整QR的区别、两种算法的计算路径都会有一个比教材文字深刻得多的理解。之后再去读QR迭代求特征值的代码你会发现自己已经能看懂每一行在做什么了。