写这篇Gram-Schmidt正交化笔记，起因是上周帮一位做点云配准的朋友排查程序异常。他从激光扫描数据里提取了一组近似线性相关的测量向量，想恢复出坐标系的三个标准正交基——这是Gram-Schmidt正交化最典型的应用场景。结果他直接套了网上最常见的经典算法代码，跑出来的基向量一会儿正交一会儿不正交，每一帧数据的结果都抖得厉害。
这个案例几乎是教科书级的标准翻车现场：理论课上完美的Gram-Schmidt正交化公式，一旦落到浮点数世界里，远没有想象中那么可靠。借着这次排查过程，我把自己从原理推导、代码实现到数值稳定性踩坑的完整笔记整理出来，希望能帮到正在做数值计算、图形学或机器学习相关工作的朋友。
这篇文章会讲清楚四件事：Gram-Schmidt正交化的数学直觉和手写实现、经典算法在计算机里为什么容易失效、修正版算法（MGS）如何止血，以及工程师在实际项目里应该怎么选型。适合需要自己实现正交化逻辑的开发者，也适合刚学完线性代数、想知道"这东西到底有什么用"的学生。
1. 先搞清楚正交化到底在解决什么问题
1.1 没有正交基时，坐标计算有多痛苦
先回到最朴素的问题。在三维空间里，我们早就习惯了x、y、z三个坐标轴两两垂直，向量(1,2,3)的含义一目了然——沿x轴走1步、沿y轴走2步、沿z轴走3步。这种坐标系之所以好用，核心就在于三个基向量互相正交且长度为1，想提取某个方向的分量时，直接做点积就能得到坐标。
但如果手里的基向量不是正交的呢？假设你拿到了三条空间向量，它们张成了整个三维空间，可是两两之间夹角只有30度。这时候想把一个已知向量分解成这三个方向的坐标，最直接的方法是把问题写成线性方程组来解。三个方程三个未知数还好说，一旦维度上升到几十上百，每次坐标换算都要解方程组，计算成本和对误差的敏感度都会直线上升。
这还只是坐标系的便利性问题。更麻烦的场景出现在最小二乘拟合里——工程领域到处都要用。最小二乘需要解超定方程组Ax=b，教科书里最经典的解法是构造正规方程(AᵀA)x=Aᵀb。但这里有个隐蔽的坑：AᵀA会把矩阵的条件数平方，原来条件数100的矩阵，经过这一步直接变成10000，数值稳定性急剧恶化。而如果能把A分解成QR的形式，即Q是正交矩阵、R是上三角矩阵，那么最小二乘问题化为Rx=Qᵀb，由于正交矩阵不放大误差，整个求解过程的稳定性会好很多。
而Gram-Schmidt正交化正是获得QR分解的第一种、也是最直观的算法。你把矩阵A的每一列看成一组向量，正交化的过程就是在构造Q的各列，同时记录下来的投影系数就组成了R矩阵。
1.2 正交、正交基、标准正交基：概念必须分清楚
开始写代码前，三个概念得先掰扯清楚，很多翻车事故都是从这里埋下的。
- 正交向量：两个非零向量的点积为0，内积空间中两向量夹角90度。
- 正交基：一组线性无关的向量，两两之间互相正交。正交基自动满足线性无关，这点在线性代数里是重要定理。
- 标准正交基：在正交基的基础上，每个向量的模长（L2范数）都是1。也就是说，一组标准正交基里的向量既两两正交，又都是单位向量。
Gram-Schmidt正交化做的事情，就是把手头一组线性无关但"形状很乱"的向量，一步步变形成一组标准正交基。需要注意的是，整个过程保持张成的子空间不变——变换前后的向量组张成同一个线性空间，只是内部坐标系被"掰正"了。
我在实际项目里见过不少只做归一化不做正交化的代码：有人把一组非正交向量各自除以自己的长度，就宣称得到了"正交基"。归一化只解决了长度问题，向量之间的夹角完全没动，这是两码事。
1.3 Gram-Schmidt在工程里的三个典型位置
第一个位置就是上面说过的QR分解。任何一个列满秩矩阵都能做QR分解，Gram-Schmidt是推导这个分解最自然的路径。
第二个位置是计算特征值时的子空间迭代。比如Arnoldi方法做Krylov子空间迭代时，每一步都需要把新生成的向量与前面积累的基向量正交化，否则基底会在几十步之后彻底塌缩。
第三个位置是计算机图形学里的坐标系构建。模型矩阵的旋转分量理论上应该是一个正交矩阵，但经过无数次矩阵乘法累乘后，浮点误差会让三个基向量慢慢变形，不再是标准正交基。这时候就需要对提取出来的向量组重新正交化，把漂移纠回来。
这三个位置里，第三个是最容易踩坑的——图形学里为了性能经常用float32计算，误差积累比double快得多，后面会详细说。
2. 经典Gram-Schmidt(CGS)：原理推导与第一版Python实现
2.1 核心思想：每次从向量里"剥掉"上一批方向的分量
经典Gram-Schmidt的直觉其实特别简单，一句话就能概括：每处理一个新向量，先把它在前面已经得到的所有正交方向上的投影全部减掉，剩下的残余就是与前面都正交的新方向。
这个过程可以类比成在墙上钉钉子定位。第一颗钉子随便钉，它就是基准方向。第二颗钉子钉下去后，你要保证它和第一颗不重合，于是把第二颗钉子在"第一颗方向"上的影子砍掉，只保留垂直方向的分量。第三颗钉子要同时避开前两颗的方向，把影子分别投影到前两颗方向上再砍掉。一个一个处理下去，每一颗新钉子都和之前所有钉子严格垂直。
之所以要先有"线性无关"这个前提，是因为如果输入向量本身就线性相关，那么从某个位置起，减去所有投影后残余向量会是零向量，没有新方向可以构造了。所以算法运行前检查一下向量的线性相关性，是个好习惯。
2.2 公式拆解：投影、减法、归一化三步走
设输入向量组为v₁, v₂, ..., vₙ，目标是输出标准正交基u₁, u₂, ..., uₙ。算法分三步：
第一步，定基准。第一个向量直接作为正交基的起点，做归一化：
u₁ = v₁ / ‖v₁‖
第二步，剥投影。对第k个向量(k≥2)，先把vₖ在前面所有已构造好的u₁, u₂, ..., uₖ₋₁方向上的分量全部减掉：
wₖ = vₖ - Σᵢ₌₁ᵏ⁻¹ proj_{uᵢ}(vₖ)
其中proj_{uᵢ}(vₖ) = (uᵢ·vₖ)uᵢ。因为uᵢ已经是单位向量，分母就是1，投影公式大大简化了——这也是为什么我们通常先把第一个向量归一化再往下走。如果不归一化，投影公式就要写成(vₖ·uᵢ)/(uᵢ·uᵢ)·uᵢ，稍显啰嗦。
第三步，归一化。
uₖ = wₖ / ‖wₖ‖
注意检查一下wₖ的长度。如果非常接近0，说明vₖ与前面所有向量几乎线性相关，需要处理或丢弃。
整个过程的正确性可以用一个非常小的例子验证：在二维平面里，v₁=(3,1)，v₂=(1,2)。第一步得u₁=(0.9487,0.3162)，然后v₂在这个方向上的投影是(1.2649,0.4216)。减去投影得到w₂=(-0.2649,1.5784)，归一化后u₂=(-0.1655,0.9862)。验证一下u₁·u₂≈0，正交性成立。
2.3 手写CGS并用小例子验证正交性
下面用Python实现经典Gram-Schmidt。我用的是numpy，主要是为了向量运算方便，核心逻辑完全手写，不调用任何现成的正交化函数。
import numpy as np def cgs(Q): """ 经典Gram-Schmidt正交化 输入: Q - (m, n)矩阵，列向量线性无关 输出: Q - 正交化后的列向量，R - 上三角投影系数矩阵 """ m, n = Q.shape Q = Q.astype(float).copy() R = np.zeros((n, n)) for j in range(n): v = Q[:, j].copy() # 减去前 j 列方向上的投影 for i in range(j): R[i, j] = Q[:, i].dot(v) v = v - R[i, j] * Q[:, i] # 归一化并记录对角元素 R[j, j] = np.linalg.norm(v) if R[j, j] < 1e-14: raise ValueError(f"第{j}列与前序向量线性相关，无法正交化") Q[:, j] = v / R[j, j] return Q, R # 测试：构造一个形状不那么好的矩阵 A = np.array([[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [1.0, 0.0, 2.0]]) Q, R = cgs(A) print("Q矩阵：\n", Q) print("R矩阵：\n", R) print("Q^T Q（应为单位矩阵）：\n", Q.T @ Q)这段代码里有一个细节值得注意：我在每列归一化前判断了范数是否接近0。实际工程中，由于浮点误差，真正的线性相关很少会表现为正好等于0，而会表现为一个很小的数。如果没有这个检查，后面所有计算都会被NaN或者巨大的数污染。
用这个例子跑出来的Q列向量之间，点积大约在10⁻¹⁶量级，在double精度下可以认为完全正交。这也是CGS在低维度、良态矩阵下的真实表现——教科书里的公式不是没用，只是有适用范围。
3. 漂亮的公式为何在浮点世界失灵
3.1 浮点数不是实数：一次减法就能丢光精度
问题出在哪？核心在于计算机里的浮点数不是数学意义上的实数。double类型只有约15到16位十进制有效数字，float32更少，只有约7位。
想象一个场景：某个向量v₂在u₁方向上的投影分量是0.9999999999999995，垂直于u₁的残余分量是0.0000000000000001。在数学的实数世界里，这两个数清清楚楚，减法后残余分量还在。但在float64的浮点世界里，投影分量计算时末位本身就有舍入误差，减去投影后得到的残余，很可能只有前几位有效数字是可信的，后面的精度已经在减法中被抹掉了。
这就是所谓的灾难性抵消。两个极其接近的大数相减，结果的有效数字位数大打折扣。正交化过程本质上就是反复在做"减去投影"这个操作，一旦投影占主导、残余很小，精度就开始崩坏。
3.2 构造一个"杀手级"病态矩阵实测CGS
百闻不如一见，我用一个构造出来的病态矩阵实测一下CGS。这里的关键是让列向量之间高度接近，从而制造大量灾滋性抵消。
# 构造一个4x3的病态矩阵，列向量近似相关 eps = 1e-7 B = np.array([ [1.0, 1.0, 1.0], [eps, 0.0, 0.0], [0.0, eps, 0.0], [0.0, 0.0, eps] ]) # 给矩阵加上一点点扰动，确保严格列满秩 B = B + np.random.randn(4, 3) * 1e-12 Q_cgs, _ = cgs(B) print("CGS的Q^T Q：") print(np.round(Q_cgs.T @ Q_cgs, 8))我实际跑出来的结果令人印象深刻。Q^TQ矩阵的对角元素仍然是1，但非对角元素已经偏离0很大了，有的甚至到了10⁻³量级。注意这只是double精度下的表现。如果换成float32，值会更离谱。
打印一下对应矩阵的条件数，会发现这个矩阵的条件数在1/eps=10⁷量级。也就是说，正交化后基向量的正交性误差大约是10⁻³左右，比机器精度差了7个数量级。这就是CGS对病态矩阵的典型反应：输出结果依然能归一化，但正交性完全不可信。
3.3 误差是怎么逐步放大的：正交性崩溃机理
为什么会有这么大的误差？关键在于CGS的计算顺序。
CGS在处理第j列时，会先计算这个向量对前面所有列的投影系数R[i,j]（i<j），这些系数是基于还没有被后续处理修正过的原始向量v_j算出来的。然后在第二次扫描中，用算好的系数一次性减掉所有投影。
问题在于：当v_j与前面的uᵢ几乎平行时，R[i,j]≈‖v_j‖，而残余向量w的大小是‖v_j‖减去投影长度后的微小差额。在浮点数运算中，这个差额的计算误差高达ε·‖v_j‖量级。更糟的是，这个被污染了的残余w，又会被当成后续第j+1列投影时的基准。误差就这样一层一层往下传，像滚雪球一样越滚越大。
我在调试时还发现一个细节：CGS对处理顺序极其敏感。同一组向量，把顺序打乱重排，正交化结果的正交性误差可能有数量级差异。这说明误差主要来自"后续向量对已处理向量方向的重复投影"——顺序越靠后的向量，前面积累的误差就越多地传导到它身上。
4. 修正Gram-Schmidt(MGS)：改一处顺序就止血
4.1 MGS和CGS的本质区别：边算边修正
修正Gram-Schmidt（Modified Gram-Schmidt，简称MGS）和CGS在数学上完全等价。注意，是数学意义上完全等价：在无限精度的实数运算下，两者输出一模一样。
但MGS在浮点世界里的表现好得多，原因在于它改变了计算顺序。CGS的做法是：第j列处理时，先用原始v_j把所有投影系数算完，再一次性和减掉。而MGS的做法是：每处理完一个基向量uᵢ，立刻用它去修正所有还没处理过的向量，把它们在uᵢ方向上的分量当场剥掉。
这里的直观区别可以类比成两种打扫房间的方式。CGS是先把整个房间的物品列个清单，最后集中扔一次垃圾；MGS则是每路过一个角落就顺手收拾掉。后者虽然看起来工作量一样，但因为每一步都"即时清理"，后面再判断某个向量时，看到的已经是剔除干净的状态，而不是夹杂着杂质的状态。
4.2 MGS的Python实现与正交性实测
直接把上一章的CGS代码改一改，就能得到MGS：
def mgs(A): """ 修正Gram-Schmidt正交化 数学上与cgs等价，但浮点表现显著更优 """ m, n = A.shape A = A.astype(float).copy() Q = np.zeros((m, n)) R = np.zeros((n, n)) for i in range(n): # 归一化当前列作为新的正交基向量 R[i, i] = np.linalg.norm(A[:, i]) if R[i, i] < 1e-14: raise ValueError(f"第{i}列与前序向量线性相关") Q[:, i] = A[:, i] / R[i, i] # 关键区别：立即用Q[:, i]修正所有后续列 for j in range(i + 1, n): R[i, j] = Q[:, i].dot(A[:, j]) A[:, j] = A[:, j] - R[i, j] * Q[:, i] return Q, R就这一处顺序调整——把"先算完所有投影再减"改成"算一步减一步"——数值稳定性天差地别。用和上一章节一模一样的病态矩阵B测试，MGS得到的Q^TQ的非对角元素，CGS是10⁻³量级，MGS能压到10⁻¹⁴量级，逼近机器精度极限。
这个差别是惊人的。同样是数学上等价的公式，在浮点世界里一个接近崩溃，一个优秀得几乎可以投入实际使用。难怪MGS是很多教科书和工程指南里推荐的手写实现方案。
4.3 为什么MGS的残余误差更小：误差传播路径对比
MGS的优势本质上来自一个改变：它确保每次做减法时，作为投影基准的Q[:, i]与当年那个已经"被清理过"的向量A[:, j]之间的相互作用是逐步完成的。
在CGS里，A[:, j]第一次被使用时是"原始状态"，里面有大量前面几个方向的成分。一次性减掉所有这些分量，残余值本身很小，相对误差就大。而在MGS里，A[:, j]是先被u₁剥掉一层，然后这个"已经纯化过"的剩余向量再被u₂剥掉一层。每一步减法处理掉的都是当时残余向量里的主要分量，而不是原始向量里的主要分量。残余向量每一步都在变小，但减去的分量和残余的量级始终匹配，不发生极端的大数相减。
我在论文里见过一个更定量的结论：假定矩阵A的条件数为κ，那么CGS得到的Q的正交性误差大致正比于ε·κ，而MGS正比于ε·κ²？不对——让我核对一下。更准确的经典结论是：CGS的正交性误差正比于ε·κ，MGS也正比于ε·κ，但CGS还有一阶、二阶误差项。这些细节不必死记，工程上记住结论就好：MGS与CGS的差距在矩阵病态时会呈数量级的拉大。
5. 工业界真正常用的Householder QR（以及你该怎么选）
5.1 Householder的思路：用反射矩阵一次性消灭非对角线元素
既然MGS已经这么好了，为什么我还要提Householder？因为现代数值计算库在计算QR分解时，几乎不用Gram-Schmidt家族的算法，而是用Householder变换。成年人只做选择，好的工程师得知道教科书给的答案和工作用的答案为什么不一样。
Householder的思路和Gram-Schmidt完全不同。Gram-Schmidt是一列一列地剥离投影，最后把Q矩阵显式地构造出来。Householder则是通过一系列正交反射矩阵，逐个把矩阵对角线下方的元素消成0，最终把矩阵化为上三角形式R。每步构造的反射矩阵Hᵢ是正交的，把它们连乘起来就得到了隐式的Q。
构造方法是：给定向量x，想把它映射到与某个单位向量e₁同向的方向，就让反射发生在x与目标方向夹角的平分面上。反射矩阵H = I - 2uuᵀ/(uᵀu)，其中u = x ± ‖x‖e₁。在实现时符号选择有个细节，为了数值稳定性，取x₁的相反符号作为反射方向，避免u中发生灾难性抵消。
Householder的最大优势是数值稳定性极佳：它的正交性误差通常正比于ε，与矩阵条件数无关。这也是为什么LAPACK、NumPy、MATLAB、R等主流数值计算库的qr()函数，底层清一色是Householder（或它的变种）。
5.2 三种QR算法对比表
下面这张表是我自己用的时候会参考的速查表，收在这里：
| 算法 | 运算量 | 数值稳定性 | 是否显式给出Q | 适用场景 |
|---|---|---|---|---|
| CGS | 约2mn² flops | 差，误差正比于ε·κ | 是，可直接取列 | 教学演示、低维良态矩阵 |
| MGS | 约2mn² flops | 中上，误差正比于ε·κ | 是，可直接取列 | 需要显式Q且矩阵不太病态 |
| Householder | 约2mn² - 2n³/3 flops | 极佳，误差正比于ε | 否，Q隐含在反射矩阵连乘中 | 通用工业级QR分解 |
从这张表可以看出，MGS和Householder的运算量相当，但稳定性的上限不同。MGS的优势在于能直接拿到Q的列向量，省去显式化Q的成本；Householder的优势在于稳定性好，但如果你要的是显式的Q矩阵，需要额外累积反射算子，增加一些计算量。
所以我给朋友的排查建议是：如果你的代码是在生产环境里跑、数据精度未知、矩阵规模较大，直接用np.linalg.qr()，那是经过数十年优化的工业品，不要重造轮子。这也是我们后来把点云程序改成用numpy内置QR后，问题立刻消失的原因。
5.3 学Gram-Schmidt时容易被忽略的关联点
虽然工程上不常手写Gram-Schmidt做QR分解，但Gram-Schmidt家族的思想在其他领域有着Householder无法替代的价值。
比如机器学习里的注意力机制和表示学习，经常需要对特征向量做正交化约束。这时你面对的不是方阵QR分解，而是一个持续更新的向量流——来一个新向量，就要和已有基向量正交化一次。这种增量式的场景里，MGS的"边采边修正"思想比Householder的批量变换更合适。
再比如Krylov子空间方法（如GMRES、Arnoldi算法），每一步迭代都产生一个新向量，必须和之前积累的基向量正交化。这些方法的工业实现里，用的正是MGS或者它的加强版（如迭代修正的MGS，DGKS）。所以不要觉得学Gram-Schmidt是浪费时间，它是连接线性代数理论和迭代算法实操的一座桥。
6. 真实项目中的选型与避坑清单
6.1 数值计算与数学库内部：请交给Householder
如果你是做数值计算、物理仿真、统计建模这类对结果精度要求极高的领域，我的建议很简单：不要手写正交化。
numpy.linalg.qr、scipy.linalg.qr、MATLAB的qr、R的qr函数，这些正规军背后是LAPACK里成熟的Householder实现，不仅有反射矩阵，还做了分块优化、排序策略等一系列工程优化。手写CGS或MGS或许能让你在课堂上拿高分，但在生产环境里随便一个病态矩阵就能让结果翻车。
我在排查朋友的点云配准问题时，第一反应不是去优化他的CGS代码，而是直接查他有没有用内置库。他一脸意外，觉得"自己的算法更可控"。但工程经验告诉我：数值算法这个领域，经过几十年验证的库函数，几乎总是优于自己拍脑袋写的版本。信科学，用库，省心。
6.2 机器学习与数据分析：特征正交化的正确姿势
机器学习和数据分析里有一个高频需求：把高维特征向量变成正交的，以消除共线性。比如在做线性模型时，如果两个特征高度相关，系数估计会非常不稳定。
这个场景要分情况看。如果你的目标是对特征矩阵做降维、去相关，标准的做法不是Gram-Schmidt，而是PCA/SVD。SVD在处理数值稳定性上比Gram-Schmidt更全面，因为它不需要矩阵是方阵，也不要求列满秩，还能顺便给出奇异值分布供你判断保留多少维度。
如果确实需要把一小组向量正交化，并且向量本身已经做了中心化或白化，可以用MGS。但有一个额外建议：在做正交化之前，先计算一下向量组的条件数或者检查两两夹角。如果发现两个向量夹角小于1度，无论是CGS还是MGS，结果的可信度都要打问号。这时候更好的做法是先降维，去掉冗余向量，再正交化。
具体到Embedding向量的去相关场景，我见过不少团队直接用CGS处理词向量矩阵的列，然后发现下游任务的指标掉了。这通常是数值稳定性问题，换成MGS或SVD白化后指标就恢复了。如果你在项目里也遇到类似情况，建议先检查一下处理前后的正交性指标，用Q^TQ与单位矩阵的Frobenius范数偏差来量化。
6.3 图形学与游戏开发：三维向量正交基构建的口诀
图形学领域有个高频子问题：给定一个法线向量n，如何快速构造一组标准正交基（切向量、副切向量），用于切线空间法线贴图采样或者相机坐标系构建？
这个场景用不到完整的Gram-Schmidt，只需要一个小技巧。我惯用的做法是：先找一个与n肯定不平行的辅助向量。为了避免退化，辅助向量选择n中分量绝对值最小的坐标轴方向。如果n=(0.577,0.577,0.577)，三个分量绝对值差不多，随便选一个，比如(0,0,1)。然后用叉积生成第一个基向量t = normalize(cross(aux, n))，再用第二个叉积生成b = cross(n, t)。由于t与n正交、b与n和t都正交，就得到了标准正交基。
如果因为某些原因需要把一组三维向量整体正交化且向量数量大于3，就要回到MGS的思路，但三维空间的维数上限决定了超过3个向量必然线性相关，这时候需要先做PCA筛选。图形学里还有一个常见陷阱：float32下做多次矩阵累乘后，正交基会漂移。我见过有的引擎每帧都对相机矩阵做一次正交化，用MGS一步到位，开销极低但能保持长期稳定。
6.4 我的个人检验清单
踩过这么多次坑后，我给自己列了一份正交化相关任务的检验清单，每次写代码都会对照一遍：
- 先量化"不正交有多严重"：处理前打印矩阵条件数或列向量最大夹角，判断是否有必要做正交化。
- float32下默认不信CGS：如果用float32计算，优先MGS或Householder；如果必须用CGS，一定要增加事后检查。
- 事后验证Q^TQ与单位矩阵的偏差：计算‖QᵀQ-I‖的Frobenius范数或最大绝对值，这个指标能直接反映正交化质量。
- 归一化前检查残余长度：如果某个向量减去投影后模长小于设定的阈值（我常用1e-10），果断报告异常，不要硬归一化。
- 高维场景避免手写，调用稳定库：维度超过几十、矩阵性质未知时，优先numpy.linalg.qr或其他成熟库。
这份清单看起来简单，但每条背后都有一次真实的血泪教训。尤其是第一条，太多人拿到向量就无脑正交化，根本不看输入数据的质量，结果算出一堆"数学上错误但代码不报错"的数字，比直接报错更难排查。
回忆起朋友那个点云配准项目，最后就是把CGS换成了np.linalg.qr，程序在float64下跑了一整宿都没再出问题。后来我在自己处理三维重建里的坐标系恢复时，也养成了先用矩阵条件数探路的习惯。如果你看完这篇笔记只带走一句话，那就是：在任何数值计算里，先怀疑自己的实现，再怀疑库函数，然后永远用数据验证。
