Gram-Schmidt正交化是线性代数里一个看起来公式很短、用起来却处处是坑的方法。我最早学它的时候老师写了几行投影公式就带过去了作业也能算但一到真正做数值实验发现结果总是不对劲明明理论上是正交的向量计算机算出来却歪得离谱。这篇笔记就是把“Gram-Schmidt正交化”这件事从原理到代码完整梳理一遍包括经典版本和修正版本的区别、手算细节、Python实现以及在最小二乘拟合里的实际用法。如果你是正在学线性代数的学生或者做数值计算时被“正交性丢失”坑过的工程师这篇笔记应该能帮你少走不少弯路。1. 正交化在解决什么问题先忘掉公式看几何直觉1.1 从直角坐标系说起为什么“正交”这么重要要理解Gram-Schmidt正交化得先理解“正交基”到底好在哪。想象你在三维空间里建坐标系如果三个坐标轴是两两垂直、长度都为1的那么任何一个向量都可以直接分解成三个方向上的独立分量互不干扰。反过来如果用三个歪歪扭扭的向量当坐标轴描述同一个点时要考虑它们之间的夹角关系计算量和出错概率都会上升。工程里也一样。求解线性方程组时如果系数矩阵的列向量接近“歪”在一起微小扰动就会被成倍放大这就是条件数大的本质。而把一组线性无关的向量转换成一组两两正交的单位向量等于把歪斜的坐标系“掰正”后续的最小二乘、特征值计算、信号处理都会变得干净很多。这里先建立一个基础概念表后面所有讨论都围绕这几项展开概念数学定义直观理解内积 ⟨a,b⟩对应分量相乘再求和反映两个向量的“重合程度”范数 ‖a‖内积开根号向量的长度正交⟨a,b⟩0两个向量相互垂直标准正交基两两正交且范数均为1一组单位长度的直角坐标轴Gram-Schmidt正交化做的事情就是把任意一组线性无关向量 a₁, a₂, …, aₙ变成一组标准正交基 q₁, q₂, …, qₙ并且保证前 k 个 q 张成的子空间和前 k 个 a 张成的子空间完全一致。这个“子空间保持不变”的性质很重要后续很多算法依赖它。1.2 三维空间里自己做一次“手工正交化”不用公式先看三步操作。手上有一组向量 a₁, a₂, a₃想把它们变成互相垂直的单位向量。第一步直接把 a₁ 当成第一个方向把它缩成长度为1的向量得到 q₁。第二步a₂ 肯定和 q₁ 有夹角把 a₂ 在 q₁ 方向上的投影减掉剩下的部分就是和 q₁ 垂直的分量记作 v₂再把 v₂ 归一化得到 q₂。第三步a₃ 分别减去它在 q₁ 方向上的投影、在 q₂ 方向上的投影剩下的分量同时垂直于 q₁ 和 q₂归一化得到 q₃。整个过程就像在一个房间里敲钉子先钉第一根再以它为基准调整第二根的方向让它和前一根垂直第三根再同时参考前两根的方向保证和它们都垂直。每一根新钉子都要把已经钉好的钉子的“影子”去掉这就是Gram-Schmidt正交化的核心直觉。这里的投影公式是向量 a 在单位向量 q 上的投影是 ⟨a,q⟩q。如果 q 不是单位向量投影是 (⟨a,q⟩/⟨q,q⟩)q但我们在算法里每一步都会先归一化所以用 ⟨a,q⟩q 就够。这一步“先归一化、再投影”是很多人在手算时容易忽略的地方。1.3 正交化与QR分解的关系把正交化写成分解形式就是 A QR。Q 的列就是前面得到的 q₁, q₂, …, qₙR 是一个上三角矩阵记录的是每一步投影时留下的系数。R 的具体含义是a_j Σ_{i1}^{j} R[i,j] * q_i。也就是说原始向量 a_j 在新坐标系 Q 下的坐标恰好就是 R 的第 j 列。由于 q₁, …, qₙ 两两正交这个分解把矩阵 A 的信息拆成“正交部分 Q”和“三角部分 R”在解方程、求最小二乘解、计算特征值时都能大显身手。QR分解和Gram-Schmidt本质上是一件事的两种说法理解了正交化就理解了QR分解的一半。2. 经典Gram-SchmidtCGS算法步骤拆解与手算实例2.1 CGS算法流程经典Gram-Schmidt的算法流程可以写成下面这份伪代码注意这里的 A 按列存放每一列是输入向量输入m×n 矩阵 A列向量线性无关 输出Qm×n 正交列、Rn×n 上三角 for j 1 to n: v A[:, j] for i 1 to j-1: R[i, j] q_i^T * A[:, j] v v - R[i, j] * q_i R[j, j] ||v|| q_j v / R[j, j]实现要点有三个。第一内层循环计算投影系数时用的是“原始列 A[:, j]”去和 q_i 做内积所以 R[i, j] q_i^T A[:, j]。第二每次算完一个 q_i就用它对 v 做一次减法把已经正交化过的方向从 v 中剔除。第三每一步结束后一定要做归一化否则后面再投影时系数会整体缩放算出来就乱了。这个算法的时间复杂度是 O(mn²)每一对 (i, j) 都要做一次长度为 m 的内积和一次向量减法。空间上只需要额外存一个 v 向量内存开销很小。2.2 一个三维手算完整流程光看伪代码容易飘我用手算一个例子完整走一遍。取三个向量a₁ (1, 1, 0)a₂ (1, 0, 1)a₃ (0, 1, 1)第一步算 q₁。v₁ a₁范数 ‖v₁‖ √(1² 1² 0²) √2所以q₁ (1/√2, 1/√2, 0)第二步算 q₂。先求投影系数R[1,2] q₁ᵀ a₂ (1/√2) × 1 (1/√2) × 0 0 × 1 1/√2于是v₂ a₂ - R[1,2] q₁ (1, 0, 1) - (1/√2)(1/√2, 1/√2, 0) (1, 0, 1) - (1/2, 1/2, 0) (1/2, -1/2, 1)计算范数‖v₂‖ √((1/2)² (-1/2)² 1²) √(1/4 1/4 1) √(3/2) √6/2归一化得到q₂ (1/√6, -1/√6, 2/√6)第三步算 q₃。需要减掉两个方向上的投影。先算系数R[1,3] q₁ᵀ a₃ (1/√2) × 0 (1/√2) × 1 0 × 1 1/√2R[2,3] q₂ᵀ a₃ (1/√6) × 0 (-1/√6) × 1 (2/√6) × 1 1/√6于是v₃ a₃ - R[1,3]q₁ - R[2,3]q₂ (0, 1, 1) - (1/2, 1/2, 0) - (1/6, -1/6, 2/6) (-2/3, 2/3, 2/3)范数‖v₃‖ √(4/9 4/9 4/9) √(12/9) 2/√3归一化q₃ (-1/√3, 1/√3, 1/√3)最终得到的 Q 和 R 是Q [ 1/√2 1/√6 -1/√3 ] [ 1/√2 -1/√6 1/√3 ] [ 0 2/√6 1/√3 ] R [ √2 1/√2 1/√2 ] [ 0 √6/2 1/√6 ] [ 0 0 2/√3 ]可以验证一下q₁ᵀq₂ 1/√12 - 1/√12 0 0q₁ᵀq₃ -1/√6 1/√6 0 0q₂ᵀq₃ -1/√18 - 1/√18 2/√18 0。三个向量两两正交范数也都是1。再把 Q 乘上 R能还原出原始的 A 列说明分解没问题。2.3 手算时容易忽略的三个细节第一个细节投影系数一定要用原始列 a_j 和 q_i 做内积不要用正在更新的 v 和 q_i 做内积。虽然在精确算术下两种写法数学上等价但语义上 R 记录的是 a_j 在 q_i 方向上的投影大小用原始列更符合定义后面讲修正算法时这个区别会放大成质变。第二个细节每次归一化之前先检查 v 的范数是不是接近0。范数接近0说明当前列和前几列几乎线性相关在调用方那边可能是数据重复、特征冗余先排查数据不要硬算。第三个细节算完分式后尽量保持根号形式到最终结果再化简中间步骤用分数可以避免小数误差积累。手算时用小数很容易在前几步看着正常最后验证正交性时差出万分之一。3. 修正Gram-SchmidtMGS一个改变命运的“顺序调整”3.1 经典版会出什么问题数值不稳定的根源理论上CGS在实数精确运算下完全正确但计算机里的浮点数有精度上限通常是15到16位有效数字。问题出在“refine”的链条上。回想CGS的做法计算 q₂ 时先用 v₂ a₂ - 投影₁再归一化计算 q₃ 时又用原始 a₃ 减掉两个投影。表面上看没问题但浮点环境下早期投影误差会残留下来。尤其当矩阵的列之间“接近相关”时比如第二列和第一列几乎平行v₂ 会变得很小此时 v₂ 里的相对误差会被放大到难以接受的程度后面所有步骤都建立在这个被污染的 v₂ 上误差像多米诺骨牌一样往后传。可以做个类比CGS是先量好所有距离再一次性修正MGS是每往前走一步立刻把后面所有的坐标都重新校准一次。后者的每一步都在“当前实际值”上操作误差不会积累得那么严重。3.2 MGS的算法流程每算出一个 q_i 就立刻更新剩余向量修正Gram-Schmidt的核心改动是把“用原始列 j 去减投影”改成“用当前残留向量 v_j 去减投影”并且每算出一个 q_i立刻更新它后面所有还没处理的列。伪代码如下输入m×n 矩阵 A 输出Q、R V A 的副本转成浮点数 for i 1 to n: R[i, i] ||V[:, i]|| Q[:, i] V[:, i] / R[i, i] for j i1 to n: R[i, j] Q[:, i]^T * V[:, j] V[:, j] V[:, j] - R[i, j] * Q[:, i]和CGS的区别主要在于内层循环的处理对象。CGS中每一列 j 在处理时独立地用自己的原始列 a_j 去减所有投影MGS中当 i1 算完 q₁ 后立刻把第2到第n列都减掉它们在 q₁ 方向上的分量此后第2列留下的 v 已经是“去掉 q₁ 方向后的残留向量”。等 i2 时再在这个残留向量上继续操作。这样做的数学原理是a_j 在 q₁ 方向上的投影可以先减掉之后计算 q₂ 时看到的 V[:,2] 已经和 q₁ 正交了继续减投影也不会把已经正交的分量重新加回来。在精确算术下MGS和CGS的结果完全一样但浮点误差的传播路径完全不同。对比项经典Gram-Schmidt (CGS)修正Gram-Schmidt (MGS)内层投影对象原始列 A[:, j]当前残留向量 V[:, j]更新时机处理每一列时才逐步减投影每得到一个 q_i 就立刻更新后面所有列浮点误差积累较明显明显改善计算复杂度O(mn²)O(mn²)教材出现频率高中等3.3 数值实验希尔伯特矩阵暴露的差距口说无凭直接上实验。希尔伯特矩阵是典型的病态矩阵它的元素是 H[i,j] 1/(ij1)矩阵规模一大条件数飙升。用它来对比CGS和MGS最有说服力。import numpy as np def cgs(A): m, n A.shape Q np.zeros((m, n)) R np.zeros((n, n)) for j in range(n): v A[:, j].astype(float).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) Q[:, j] v / R[j, j] return Q, R def mgs(A): m, n A.shape Q np.zeros((m, n)) R np.zeros((n, n)) V A.astype(float).copy() for i in range(n): R[i, i] np.linalg.norm(V[:, i]) Q[:, i] V[:, i] / R[i, i] for j in range(i 1, n): R[i, j] Q[:, i] V[:, j] V[:, j] V[:, j] - R[i, j] * Q[:, i] return Q, R n 10 A np.array([[1.0 / (i j 1) for j in range(n)] for i in range(n)]) Q_cgs, _ cgs(A) Q_mgs, _ mgs(A) err_cgs np.max(np.abs(Q_cgs.T Q_cgs - np.eye(n))) err_mgs np.max(np.abs(Q_mgs.T Q_mgs - np.eye(n))) print(CGS 正交性误差:, err_cgs) print(MGS 正交性误差:, err_mgs)在我的机器上n10 的希尔伯特矩阵CGS的正交性误差能达到10⁻⁷这个量级而MGS能到10⁻¹⁵附近。别小看这8个数量级的差距在特征值计算或深度学习里正交性误差会被后续迭代不断放大最后等于是拿着一把歪尺子在量数据。所以如果只能记一个结论那就是工程上用MGS别用CGS。4. 代码实战从零手写CGS与MGS并用于最小二乘拟合4.1 一份可以直接拿来用的numpy实现把上面两个函数整理成完整模块加上类型转换和边界检查可以在项目里直接调用。import numpy as np def gram_schmidt_cgs(A): 经典Gram-Schmidt正交化返回Q(列正交)和R(上三角)。 A np.asarray(A, dtypefloat) m, n A.shape if n m: raise ValueError(列数不能大于行数否则必然线性相关) 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 - R[i, j] * Q[:, i] norm_v np.linalg.norm(v) if norm_v 1e-12: raise ValueError(f第{j}列与前序列接近线性相关) R[j, j] norm_v Q[:, j] v / norm_v return Q, R def gram_schmidt_mgs(A): 修正Gram-Schmidt正交化数值稳定性更好。 A np.asarray(A, dtypefloat) m, n A.shape if n m: raise ValueError(列数不能大于行数否则必然线性相关) Q np.zeros((m, n)) R np.zeros((n, n)) V A.copy() for i in range(n): norm_v np.linalg.norm(V[:, i]) if norm_v 1e-12: raise ValueError(f第{i}列与前序列接近线性相关) Q[:, i] V[:, i] / norm_v R[i, i] norm_v for j in range(i 1, n): R[i, j] Q[:, i] V[:, j] V[:, j] - R[i, j] * Q[:, i] return Q, R这里做了一件CGS和MGS都需要的防护归一化之前检查范数是否小于一个阈值。这一步看起来多余但能省下大量调试时间尤其是自动处理数据管道时某列数据全为零或者完全重复程序会在这里直接报错而不是往下输出一堆NaN。4.2 验证正确性正交性、重构与误差度量写都写了必须验证。第一步验证基本的正交性用前面手算的例子A np.array([ [1.0, 1.0, 0.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0], ]) Q1, R1 gram_schmidt_cgs(A) Q2, R2 gram_schmidt_mgs(A) print(CGS 正交性误差:, np.max(np.abs(Q1.T Q1 - np.eye(3)))) print(MGS 正交性误差:, np.max(np.abs(Q2.T Q2 - np.eye(3)))) print(CGS 重构误差:, np.max(np.abs(Q1 R1 - A))) print(MGS 重构误差:, np.max(np.abs(Q2 R2 - A)))理论上的输出应该接近机器精度四个误差都在10⁻¹⁵量级。这证明了两个算法在良性矩阵上都能正确重建原矩阵。进一步用随机矩阵做压力测试重复100次每次生成100×30的随机矩阵记录两种算法的正交性误差最大值。我实测下来CGS的正交性误差通常在10⁻¹⁴到10⁻¹³之间MGS也在这个范围两者在随机良态矩阵上差别不大。这就印证了一个经验CGS和MGS的差距只有在病态矩阵上才会被放大平时小矩阵怎么算都对但工程数据谁也不知道什么时候会碰到病态情况所以默认用MGS更稳妥。4.3 实战应用用QR分解求解最小二乘拟合问题纸上谈兵到此为止看一个真实场景。假设有一组实验数据点 (x, y)我怀疑它们符合抛物线关系 y c₀ c₁x c₂x²用带噪声的数据做拟合。rng np.random.default_rng(42) x np.linspace(0, 1, 20) y_true 1.5 2.0 * x - 0.8 * x**2 y y_true 0.05 * rng.standard_normal(x.size) A np.vstack([np.ones_like(x), x, x**2]).T Q, R gram_schmidt_mgs(A) c_mgs np.linalg.solve(R, Q.T y) c_lstsq, _, _, _ np.linalg.lstsq(A, y, rcondNone) print(MGSQR 求解系数:, c_mgs) print(np.linalg.lstsq 系数:, c_lstsq)两种方法求出来的系数应该非常接近差距在10⁻¹⁵量级。这里的核心思想是最小二乘问题 min ‖Ax - y‖₂其解满足正规方程 AᵀAx Aᵀy但直接构造 AᵀA 会让条件数变成原来的平方数值上非常危险。用 QR 分解把 A 换成 QR目标函数变成 ‖Qy - Rx‖₂因为 Q 是正交矩阵范数保持不变问题归约为一个上三角方程组的求解避免了平方条件数的灾难。在实际项目中用这个方法做曲线拟合、传感器校准、信号去趋势都很常见。特别是数据量上千、特征数量几百的时候直接求 AᵀA 的逆很容易炸QR路径则稳定很多。4.4 性能与精度权衡什么时候用CGS也没问题既然MGS更稳是不是CGS就该扔进垃圾桶也不完全是。CGS有一个MGS没有的优势内层循环可以批量向量化。如果把所有投影系数一次性算出来用矩阵乘法完成减法在GPU或者优化良好的BLAS库上CGS的实际吞吐量可能比逐列更新的MGS更高。前提是矩阵的条件数足够好你能接受10⁻¹²量级的正交性误差。如果只是做教学演示、快速原型、或者矩阵本身非常良性CGS完全够用。更稳妥的选择其实是直接用householder变换实现QR分解这也是LAPACK标准库和numpy.linalg.qr内部采用的方法。Householder比MGS更稳数值性质极好性能也高。手写Gram-Schmidt的价值在于理解原理、理解数值误差从哪来以及在一些无法引入成熟库的嵌入式或教学环境里自己实现。5. 常见问题与排查技巧实录5.1 QᵀQ不等于单位矩阵先查条件数和数值精度这是出现频率最高的问题。很多人发现手写的QR分解里 QᵀQ 和单位矩阵差了10⁻⁸甚至10⁻⁵第一反应是自己写错了。我的排查顺序是先打印输入矩阵的条件数判断矩阵是不是病态再换用numpy.linalg.qr的标准库结果对比最后再检查自己的实现里是不是用了CGS版本。一个经验法则是如果 QᵀQ 的非对角元达到10⁻⁸量级并且用的是CGS先换成MGS试试。如果换成MGS之后仍然差再考虑更底层的问题比如数据本身存在重复列、或者矩阵本来就不是满秩的。不要在一开始就怀疑浮点环境99%的情况是算法或者数据的问题。5.2 归一化时出现NaN或无穷大这种问题多半出现在数据自动处理流程中。某列向量全部为零或者某列恰好是前面若干列的线性组合导致v的范数等于0或接近0归一化时就会出现除零或者数值溢出。解决方法是前置一个秩检查或者在循环里加一个范数阈值判断。我有一次处理传感器数据时某个通道因为硬件故障输出全零Gram-Schmidt直接返回NaN后面的模型全部崩掉。加了阈值检查后程序会在进入NaN之前报错排查效率高了很多。阈值取多少要根据数据量纲来我的经验是先用矩阵整体的量级乘以1e-12作为默认阈值。5.3 手写QR和numpy.linalg.qr结果为什么不一样经常有人发现自己手写MGS得到的Q、R和numpy.linalg.qr的结果对不上以为代码有问题。实际上QR分解的结果不唯一任意一个对角元素全为±1的对角矩阵DQ′ QDR′ DR仍然是合法的QR分解。也就是说某些列整体乘以-1仍然是正交基。验证方法很简单把两个Q的每一列对比看是否只差一个符号或者直接对比 QR 是否都等于原矩阵。只要重构误差小、正交性好就不用纠结具体列符号是否完全一致。如果非要和库函数对齐可以检查并让每个q的第一非零元素都变为正数但一般没必要。5.4 到底该不该自己写工具选型建议我把问题整理成一张速查表方便直接对号入座场景推荐做法原因学习原理、教学演示手写CGS和MGS理解误差来源和算法差异工程数据拟合numpy.linalg.qrHouseholder数值稳定、性能好嵌入式/无库环境手写MGS代码简单且稳定超大矩阵GPU加速库函数或专门的正交化模块手写循环难发挥硬件性能低精度实验快速原型CGS也行实现简单先跑通再说一个更进阶的技巧是给MGS加上“重正交化”在算出q_i后再检查一下新得到的q_i和之前已算出的q_j的内积如果超过阈值就再减一次投影。这个技巧在一些极端病态的问题里可以把正交性误差再压低几个数量级代价是计算量稍微增加但有时能救命。最后说点我自己的体会。刚开始学Gram-Schmidt时我也觉得MGS只是把CGS的顺序变了一下没必要较真。直到一次做矩阵特征值迭代用CGS预处理基底迭代十几轮后结果完全发散换成MGS立刻稳定下来才真正理解“顺序”两个字在数值计算里有多重要。这块内容值得每个做数值计算的人亲手写一遍、跑一遍、踩一遍坑比看十遍公式都管用。