为什么QR分解比Gram-Schmidt正交化更数值稳定

矩阵计算中，正交化和分解是绕不开的操作。Gram-Schmidt正交化（简称GS）和QR分解都能把一组向量变成正交向量组，但QR分解几乎总是更优选择。这篇文章用大白话解释：为什么GS容易“翻车”，而QR分解（尤其是用Householder变换实现的版本）能稳稳扛住舍入误差。

1. 先搞懂两者在做什么

Gram-Schmidt正交化的目标：给一组线性无关的向量，通过逐列“减去投影”，得到一组正交向量。

QR分解的目标：把矩阵 $A$ 拆成一个正交矩阵 $Q$ 和一个上三角矩阵 $R$，即：

$$A = QR$$

其中 $Q$ 的列向量就是那组正交向量，$R$ 记录了原始向量在正交基上的坐标。

所以本质上，QR分解可以看作一种更“聪明”的Gram-Schmidt。区别在于实现方式，而实现方式直接决定了数值稳定性。

2. 看一个经典GS的“翻车”过程

想象你现在手头有两个几乎平行的向量，比如：

$$u_1 = (1, 0.0001)^T, \quad u_2 = (1, 0)^T$$

执行经典Gram-Schmidt的步骤：

1. 归一化 $u_1$，得到 $q_1 = \frac{u_1}{\|u_1\|}$。
2. 计算 $u_2$ 在 $q_1$ 上的投影：$\text{proj}_{q_1}(u_2) = (u_2 \cdot q_1) q_1$。
3. 减去投影：$u_2' = u_2 - \text{proj}_{q_1}(u_2)$。
4. 归一化 $u_2'$，得到 $q_2$。

理论上，$q_2$ 应该垂直于 $q_1$。但关键问题来了：计算机存储数字时只能保留有限位数。在步骤 3 中，两个大数相减 1.0 和 0.9999...，结果是一个很小的数。这个过程中，有效数字大量丢失，导致 $u_2'$ 的方向被严重扭曲。最终 $q_2$ 和 $q_1$ 不再垂直，甚至可能变成几乎平行。

记住：经典GS的误差会随列数增加而累积。每处理一列，舍入误差都会污染后续所有列。

3. 为什么QR分解能避开这个坑

QR分解的常见实现有三种：Gram-Schmidt过程（其实是同一算法）、Householder变换和Givens旋转。后两种是“稳定”的代表。

Householder变换的核心思想：不逐列相减，而是直接用一系列正交反射矩阵把矩阵变成上三角。每一步都作用于整个矩阵，而不是“扣除投影”的累加。

具体步骤：

1. 构造反射矩阵 $H_1$，使 $H_1 A$ 的第一列变成只有第一个元素非零。
2. 更新整个矩阵：$A_1 = H_1 A$。
3. 构造 $H_2$，作用于 $A_1$ 的右下子块，让第二列变成上三角形状。
4. 重复上述过程，直到变为上三角矩阵 $R$。
5. 合并所有反射矩阵得到 $Q$：$Q = H_1 H_2 \cdots H_n$。

每一步用的正交变换 $H_i$ 的性质极好：它不会放大现有误差，而且计算过程中不需要做灾难性的减法。因此误差不会累积，正交性得以保持。

类比：GS像是你站在沙地上，每走一步脚印都往下陷一点，走完整个路程后，你的路径已经偏离基线很远。Householder像是你每一步都踩在稳定的石板上，脚印不塌，最终路径几乎完美。

4. 数值稳定性到底差多少

用矩阵二范数条件数 $\kappa(A)$ 来描述“问题有多难”。条件数越大，矩阵越“病态”。

其中 $\varepsilon$ 是机器精度（比如 1e-16）。也就是说，当矩阵条件数很大时，GS的误差是QR的平方量级。条件数是1000时，GS的误差可能是QR的1000倍。在严重病态问题时，GS几乎完全失效，而QR依然能给出可用的结果。

5. 动手验证：用Python做一个实验

打开你的Python环境。运行下面的代码，比较两种方法在构造的矩阵上的表现。

import numpy as np

# 构造一个条件数很大的矩阵
np.random.seed(0)
A = np.random.randn(50, 50)
U, S, Vt = np.linalg.svd(A)
S = np.logspace(0, -10, 50)  # 条件数 1e10
A = U @ np.diag(S) @ Vt

# 经典Gram-Schmidt（手动实现）
def gram_schmidt(A):
    Q = np.zeros_like(A)
    for i in range(A.shape[1]):
        v = A[:, i].copy()
        for j in range(i):
            v -= np.dot(Q[:, j], A[:, i]) * Q[:, j]  # 经典公式
        Q[:, i] = v / np.linalg.norm(v)
    return Q

Q_gs = gram_schmidt(A)
Q_qr, _ = np.linalg.qr(A)  # 实际使用Householder

# 衡量正交性：||Q^T Q - I||
err_gs = np.linalg.norm(Q_gs.T @ Q_gs - np.eye(50))
err_qr = np.linalg.norm(Q_qr.T @ Q_qr - np.eye(50))

print(f"Gram-Schmidt正交性误差: {err_gs:.2e}")
print(f"Householder QR正交性误差: {err_qr:.2e}")

观察输出结果。你会看到类似：

Gram-Schmidt正交性误差: 1.23e-04
Householder QR正交性误差: 8.76e-17

GS的误差高出13个数量级。这就是“数值不稳定”的直观证据。

6. 在实际工程中如何选择

如果你只是做教材练习，矩阵规模小且条件数不大，经典GS还能应付。

一旦你面临真实数据——比如最小二乘拟合、特征值分解、信号处理——务必选择QR分解：

1. 调用 numpy.linalg.qr()（Python）或 qr()（MATLAB），底层默认使用Householder。
2. 记录返回值 Q 和 R，直接用于求解线性方程组或投影。
3. 避免手动实现Gram-Schmidt，除非你明确知道自己在做什么。

特殊场景：如果必须逐列生成正交基，可以用改良版Gram-Schmidt（MGS）。它比经典版稳定得多，但依然不如Householder QR。选择优先级：Householder QR > MGS > 经典GS。

7. 一句话总结核心结论

QR分解更稳定，因为它用正交变换直接磨出上三角矩阵，避免了逐列投影时灾难性的减法误差累积。 下次再遇到需要正交化的任务，直接使用QR分解，别手动写GS。