# 从零实现QR分解：CGS、MGS和Householder算法对比与实战（附Python代码） 如果你曾经在机器学习、计算机视觉或者科学计算领域里处理过线性最小二乘问题，那么QR分解这个名字对你来说一定不陌生。它不仅是求解线性方程组、计算特征值的核心工具，更是许多数值算法的基石。但你是否真正理解过，那些封装在`numpy.linalg.qr`或`scipy.linalg.qr`背后的算法究竟是如何工作的？为什么有的算法在数值上更稳定，而有的则容易在浮点运算中“翻车”？ 今天，我们不依赖任何现成的库，从最基础的线性代数原理出发，亲手实现三种主流的QR分解算法：经典的Gram-Schmidt正交化（CGS）、改进的Gram-Schmidt正交化（MGS）以及Householder变换。我会带你深入算法的每一步，用Python代码将理论落地，并通过实际的数值实验，直观感受它们在精度、稳定性和效率上的差异。这篇文章适合那些不满足于“调包”、渴望理解底层机制，并希望在自己的项目中实现定制化矩阵运算的开发者。 ## 1. QR分解的核心思想与几何直观 在深入算法细节之前，我们有必要先厘清QR分解究竟要解决什么问题。给定一个**m × n**的实矩阵**A**（通常假设 m ≥ n，即行数不少于列数），QR分解的目标是将其分解为两个矩阵的乘积： **A = Q R** 其中： * **Q** 是一个 **m × m** 的正交矩阵（当进行“完全分解”时），满足 **QᵀQ = I**。这意味着它的列向量构成了一组标准正交基。 * **R** 是一个 **m × n** 的上三角矩阵。在简化形式中，我们通常只取 **R** 的前 n 行 n 列（一个 n × n 的上三角方阵），以及 **Q** 的前 n 列（一个 m × n 的列正交矩阵）。 从几何视角看，QR分解可以理解为对矩阵 **A** 的列空间进行了一次“重新标架”。**A** 的每一列原本是站在标准基下的向量。**Q** 的列提供了一组新的、彼此垂直且长度为1的坐标轴（标准正交基）。而 **R** 的上三角元素，则记录了 **A** 的每一列在这组新坐标轴下的“坐标”。因为新基是正交的，所以 **A** 的列向量在新基下的表示会变得非常简洁——每个向量只与“前面”的基向量有关，这就形成了上三角结构。 这种分解的巨大价值体现在多个方面： * **求解线性方程组**：对于系统 **A x = b**，将其转化为 **Q R x = b**。由于 **Q** 是正交的，方程变为 **R x = Qᵀ b**，而 **R** 是上三角矩阵，可以通过简单的回代法快速求解。 * **最小二乘问题**：在拟合数据时，我们经常需要最小化 **||A x - b||₂**。利用QR分解，目标函数变为 **||Qᵀ(A x - b)||₂ = ||R x - Qᵀ b||₂**。由于正交变换不改变向量长度，问题简化为求解一个上三角系统，数值上比直接计算 **AᵀA** 更稳定。 * **特征值计算**：著名的QR迭代算法就是通过反复进行QR分解来逼近矩阵的特征值。 为了量化不同算法的性能，我们首先定义一个简单的测试矩阵，并在后续章节中用它来验证我们的实现： ```python import numpy as np # 定义一个条件数较大的测试矩阵，以凸显数值稳定性问题 def generate_test_matrix(m, n, condition_number=1e6): """生成一个指定尺寸和条件数的测试矩阵""" np.random.seed(42) # 固定随机种子以便复现 U, _ = np.linalg.qr(np.random.randn(m, m)) # 随机正交矩阵 U V, _ = np.linalg.qr(np.random.randn(n, n)) # 随机正交矩阵 V # 构建奇异值矩阵 S s = np.logspace(0, np.log10(condition_number), n) S = np.zeros((m, n)) for i in range(n): S[i, i] = s[i] A = U @ S @ V.T return A.astype(np.float64) # 确保是双精度浮点数 # 生成一个 5x3 的测试矩阵 A_test = generate_test_matrix(5, 3, condition_number=1e4) print("测试矩阵 A (5x3):") print(A_test) print(f"\n矩阵 A 的条件数 (近似): {np.linalg.cond(A_test):.2e}") ``` > **提示**：我们特意生成了一个条件数较大的矩阵。条件数衡量了矩阵对输入误差的敏感程度，条件数越大，数值计算越容易不稳定，是检验算法鲁棒性的好工具。 ## 2. 经典Gram-Schmidt正交化 (CGS)：直观但脆弱 经典Gram-Schmidt算法（CGS）的思想直接源于线性代数教材中构造正交基的过程。给定矩阵 **A = [a₁, a₂, ..., aₙ]**，我们希望构造一组标准正交向量 **q₁, q₂, ..., qₙ**，使得每个 **aⱼ** 都可以由 **q₁** 到 **qⱼ** 线性表示。这个过程是逐列进行的： 1. 第一个向量直接归一化：**q₁ = a₁ / ||a₁||**，同时记录缩放因子 **r₁₁ = ||a₁||**。 2. 对于第 j 个向量 **aⱼ** (j > 1)： * 计算它在前 j-1 个正交方向上的投影分量：对于每个 i < j，计算投影系数 **rᵢⱼ = qᵢᵀ aⱼ**。 * 从 **aⱼ** 中减去所有这些投影，得到与前面所有 **qᵢ** 都正交的向量：**vⱼ = aⱼ - Σᵢ₌₁ʲ⁻¹ rᵢⱼ qᵢ**。 * 将这个正交向量归一化：**qⱼ = vⱼ / ||vⱼ||**，记录 **rⱼⱼ = ||vⱼ||**。 将所有 **rᵢⱼ** 按 i ≤ j 排列，就构成了上三角矩阵 **R**。所有 **qⱼ** 按列排列，就构成了矩阵 **Q**。 下面是用Python实现的CGS算法： ```python def qr_cgs(A): """ 使用经典Gram-Schmidt正交化进行QR分解。 参数: A: 输入矩阵，形状 (m, n)， m >= n 返回: Q: 正交矩阵，形状 (m, n) R: 上三角矩阵，形状 (n, n) """ m, n = A.shape Q = np.zeros((m, n), dtype=A.dtype) R = np.zeros((n, n), dtype=A.dtype) for j in range(n): # 复制当前列 v = A[:, j].copy().astype(np.float64) # 减去在前面的所有正交向量上的投影 for i in range(j): # 计算投影系数 R[i, j] = np.dot(Q[:, i], A[:, j]) # 减去投影分量 v = v - R[i, j] * Q[:, i] # 计算当前列的范数并归一化 R[j, j] = np.linalg.norm(v) if R[j, j] > 1e-15: # 避免除零 Q[:, j] = v / R[j, j] else: Q[:, j] = v # 零向量或线性相关，保持为零 # 在实际应用中，这里可能需要处理秩亏的情况 return Q, R ``` 让我们用测试矩阵来运行它，并检查分解的质量： ```python Q_cgs, R_cgs = qr_cgs(A_test) # 检查分解的正确性：A ≈ Q * R 吗？ reconstruction_error_cgs = np.linalg.norm(A_test - Q_cgs @ R_cgs, 'fro') print(f"CGS 重构误差 (Frobenius范数): {reconstruction_error_cgs:.2e}") # 检查 Q 的正交性：QᵀQ ≈ I 吗？ orthogonality_error_cgs = np.linalg.norm(Q_cgs.T @ Q_cgs - np.eye(3), 'fro') print(f"CGS 正交性误差 (QᵀQ - I): {orthogonality_error_cgs:.2e}") ``` 如果你运行这段代码，可能会发现即使对于一个中等条件数的矩阵，`orthogonality_error_cgs`（正交性误差）也可能远大于 `reconstruction_error_cgs`（重构误差）。**这正是CGS算法著名的数值缺陷**。 **CGS的数值不稳定性根源**： 问题出在 `v = v - R[i, j] * Q[:, i]` 这一步。在浮点运算中，当我们从一个向量中连续减去多个其他向量的分量时，由于舍入误差，计算出的新向量 **v** 可能无法完全正交于之前所有的 **Q[:, i]**。更糟糕的是，这种正交性误差会随着迭代累积和放大。在几何上，可以想象为由于计算精度的限制，我们每一步“剔除”投影都不够干净，残留的微小分量在后续步骤中不断污染新的正交方向。 为了更直观地对比，我们先看看后续更稳定的算法结果，但问题的严重性已经显现。当矩阵条件数很大或列向量接近线性相关时，CGS算法产生的 **Q** 矩阵可能严重偏离正交性，进而导致基于此分解的后续计算（如求解最小二乘问题）完全失败。 ## 3. 改进的Gram-Schmidt正交化 (MGS)：一个关键的顺序调整 改进的Gram-Schmidt算法（MGS）是针对CGS数值缺陷的一个巧妙而有效的修补。它的核心洞察在于：**计算顺序至关重要**。CGS是一次性计算向量 **aⱼ** 在所有先前方向上的投影系数，然后一次性减去。而MGS采用了一种“逐次投影”的策略。 MGS算法的过程如下，注意其内循环的不同： 1. 初始化：令 **V = A**（副本）。 2. 对于 i = 1 到 n： a. 对第 i 列进行归一化：**rᵢᵢ = ||V[:, i]||**, **Q[:, i] = V[:, i] / rᵢᵢ**。 b. **关键步骤**：对于每一个 j = i+1 到 n： * 计算当前正交向量 **Q[:, i]** 对后续向量 **V[:, j]** 的投影系数：**rᵢⱼ = Q[:, i]ᵀ V[:, j]**。 * **立即**从 **V[:, j]** 中减去这个投影分量：**V[:, j] = V[:, j] - rᵢⱼ Q[:, i]**。 注意，在步骤b中，一旦我们得到了 **Q[:, i]**，就立刻用它去“清理”所有尚未处理的列向量 **V[:, j]**。这样，当后续轮到处理第 j 列时，它里面已经不含 **Q[:, i]** 方向的分量了。这种“即时的减法”极大地减少了舍入误差的累积。 以下是MGS的Python实现： ```python def qr_mgs(A): """ 使用改进的Gram-Schmidt正交化进行QR分解。 参数: A: 输入矩阵，形状 (m, n)， m >= n 返回: Q: 正交矩阵，形状 (m, n) R: 上三角矩阵，形状 (n, n) """ m, n = A.shape # 创建V作为A的副本，我们将在V上直接操作 V = A.copy().astype(np.float64) Q = np.zeros((m, n), dtype=A.dtype) R = np.zeros((n, n), dtype=A.dtype) for i in range(n): # 计算当前列的范数并归一化，得到 q_i R[i, i] = np.linalg.norm(V[:, i]) if R[i, i] > 1e-15: Q[:, i] = V[:, i] / R[i, i] else: Q[:, i] = V[:, i] # 处理秩亏 # 立即用 q_i 去修正所有后续的列 for j in range(i+1, n): R[i, j] = np.dot(Q[:, i], V[:, j]) V[:, j] = V[:, j] - R[i, j] * Q[:, i] return Q, R ``` 现在，让我们对比MGS和CGS在同一个问题上的表现： ```python Q_mgs, R_mgs = qr_mgs(A_test) reconstruction_error_mgs = np.linalg.norm(A_test - Q_mgs @ R_mgs, 'fro') orthogonality_error_mgs = np.linalg.norm(Q_mgs.T @ Q_mgs - np.eye(3), 'fro') print("=== MGS 算法性能 ===") print(f"MGS 重构误差: {reconstruction_error_mgs:.2e}") print(f"MGS 正交性误差: {orthogonality_error_mgs:.2e}") print("\n=== 与 CGS 对比 (误差比) ===") print(f"重构误差比 (CGS/MGS): {reconstruction_error_cgs/reconstruction_error_mgs:.2f}") print(f"正交性误差比 (CGS/MGS): {orthogonality_error_cgs/orthogonality_error_mgs:.2f}") ``` 在我的测试中，MGS的正交性误差通常比CGS小好几个数量级。这个简单的顺序调整，带来了数值稳定性质的提升。MGS是许多需要中等精度QR分解场景下的可靠选择，它比CGS更稳定，同时又比接下来要介绍的Householder变换更易于理解（在某些并行化实现中也有其优势）。 为了更系统地比较，我们可以设计一个实验，测试在不同条件数下两种算法的正交性误差： | 矩阵条件数 | CGS 正交性误差 | MGS 正交性误差 | 改进倍数 (CGS/MGS) | | :--- | :--- | :--- | :--- | | 10² | ~1e-15 | ~1e-15 | ~1 | | 10⁴ | ~1e-11 | ~1e-15 | ~1e4 | | 10⁶ | ~1e-7 | ~1e-15 | ~1e8 | | 10⁸ | ~1e-3 | ~1e-14 | ~1e11 | > **注意**：上表为示意性数据，实际误差与矩阵的具体元素也有关，但趋势是明确的：随着条件数增大，CGS的误差急剧恶化，而MGS则能保持接近机器精度的优良正交性。 ## 4. Householder变换：基于反射的工业级算法 如果说Gram-Schmidt系列是“建设性”的（一步步构建正交基），那么Householder变换则是“破坏性”或“消元性”的。它不直接构造 **Q**，而是通过一系列精心设计的正交反射变换，直接将原矩阵 **A** “雕刻”成上三角矩阵 **R**。每一个Householder变换的目标是**将当前列向量下方的所有元素清零**。 **Householder变换的几何本质**： 给定一个向量 **x**，我们想找到一个正交变换 **H**，使得 **Hx** 与某个坐标轴（如第一个坐标轴）对齐，即除了第一个分量外，其他分量全为零。**H** 构造为一个反射矩阵： **H = I - 2 u uᵀ / (uᵀu)** 其中向量 **u** 是反射超平面的法向量。通过选择 **u = x ± ||x|| e₁**（其中 **e₁** 是第一个标准基向量），可以证明 **Hx** 就等于 **∓||x|| e₁**，从而实现了消元。符号的选择通常取 **-sign(x₁)||x||** 以避免数值上的相减消去。 **基于Householder的QR分解步骤**： 1. 对矩阵 **A** 的第一列，构造一个Householder矩阵 **H₁**，使得 **H₁A** 的第一列除了第一个元素外全为零。 2. 接着，忽略第一行第一列，对右下角的子矩阵重复此过程，构造 **H₂**。注意 **H₂** 需要嵌入到一个更大的单位矩阵中，以保证它不影响已被清零的部分。 3. 持续进行，直到将 **A** 化为上三角矩阵 **R**。即 **Hₙ ... H₂ H₁ A = R**。 4. 由于每个 **Hᵢ** 都是对称且正交的（**Hᵢ = Hᵢᵀ = Hᵢ⁻¹**），所以 **Q = H₁ H₂ ... Hₙ**。在实际存储时，我们通常不会显式构造出完整的 **Q**，而是保存每个 **Hᵢ** 对应的 **u** 向量，在需要时再进行运算。 下面是Householder QR分解的实现，它显式计算了 **Q** 矩阵： ```python def qr_householder(A): """ 使用Householder变换进行QR分解。 参数: A: 输入矩阵，形状 (m, n)， m >= n 返回: Q: 正交矩阵，形状 (m, m) R: 上三角矩阵，形状 (m, n) """ m, n = A.shape R = A.copy().astype(np.float64) Q = np.eye(m, dtype=np.float64) # 初始化为单位矩阵 for k in range(n): # 取当前列的下半部分 x = R[k:, k] # 计算范数，并选择符号以避免数值问题 norm_x = np.linalg.norm(x) if norm_x == 0: continue # 如果已经是零向量，跳过 # 构造Householder向量 u # 使用 -sign(x[0])*norm_x 来增强稳定性 alpha = -np.sign(x[0]) * norm_x u = x.copy() u[0] = u[0] - alpha beta = np.dot(u, u) # 如果 beta 太小，说明 u 几乎是零向量，跳过反射 if beta < 1e-20: continue # 应用Householder变换到 R 的剩余部分 # H = I - (2/beta) * u * uᵀ # 计算 w = (2/beta) * R[k:, k:].T @ u 会更高效 for j in range(k, n): # 计算 u 与 R 第 j 列下半部分的内积 gamma = np.dot(u, R[k:, j]) # 更新 R 的第 j 列 R[k:, j] = R[k:, j] - (2.0 * gamma / beta) * u # 同时将变换累积到 Q 矩阵上 # Q = Q * H， 由于 H 对称，等价于 Q = Q - (2/beta) * (Q[:, k:] @ u) * uᵀ for i in range(m): # 计算 Q 的第 i 行与 u 的内积（只涉及 k 行之后） gamma = np.dot(Q[i, k:], u) Q[i, k:] = Q[i, k:] - (2.0 * gamma / beta) * u # 通常我们返回“经济型”QR分解，即 Q 的前 n 列和 R 的前 n 行 return Q[:, :n], R[:n, :n] ``` 让我们评估Householder算法的表现： ```python Q_house, R_house = qr_householder(A_test) # 注意这里返回的Q是m x m，我们取前n列进行经济型分解对比 Q_house_econ = Q_house[:, :3] reconstruction_error_house = np.linalg.norm(A_test - Q_house_econ @ R_house, 'fro') orthogonality_error_house = np.linalg.norm(Q_house_econ.T @ Q_house_econ - np.eye(3), 'fro') print("=== Householder 算法性能 ===") print(f"Householder 重构误差: {reconstruction_error_house:.2e}") print(f"Householder 正交性误差: {orthogonality_error_house:.2e}") print("\n=== 三种算法正交性误差总结 ===") print(f"CGS: {orthogonality_error_cgs:.2e}") print(f"MGS: {orthogonality_error_mgs:.2e}") print(f"Householder: {orthogonality_error_house:.2e}") ``` 在绝大多数情况下，Householder变换能给出**最优的数值稳定性**。它的正交性误差通常接近机器精度（对于双精度浮点数约为1e-15），即使面对病态矩阵也是如此。这是因为Householder变换本质上是精确的反射操作，舍入误差的传播方式比Gram-Schmidt中的连续减法更可控。 **算法特性对比表**： | 特性 | 经典Gram-Schmidt (CGS) | 改进Gram-Schmidt (MGS) | Householder变换 | | :--- | :--- | :--- | :--- | | **核心思想** | 逐列正交化，一次性减去所有投影 | 逐列正交化，立即减去每个投影 | 逐列反射消元 | | **数值稳定性** | 差，误差易累积 | 好，比CGS有显著提升 | **优秀**，工业标准 | | **计算复杂度** | O(mn²) | O(mn²) | O(mn² - n³/3) (略优) | | **存储需求** | 需额外存储Q | 可原地操作（覆盖A） | 可原地存储反射向量 | | **显式Q矩阵** | 直接得到 | 直接得到 | 需额外计算或累积 | | **适用场景** | 教学、理解原理 | 中等精度需求、某些并行架构 | **高精度、高稳定性需求（默认选择）** | | **几何解释** | 逐步构建正交基 | 逐步构建正交基（顺序优化） | 一系列镜像反射 | ## 5. 实战应用：用自实现的QR分解求解最小二乘问题 理论最终要服务于实践。让我们用一个简单的线性回归例子，来验证我们自实现的QR分解能否正确工作，并比较不同算法在求解实际问题时的精度。 假设我们有一组数据点，想用一条直线 y = β₀ + β₁ x 来拟合。这对应于求解超定方程组 **A β ≈ b**，其中 **A** 的第一列是全1（对应截距），第二列是x值；**b** 是y值；**β** 是待求系数。 我们将使用Householder算法（最稳定）来求解，并与NumPy的权威实现`np.linalg.lstsq`进行对比。 ```python def solve_least_squares_qr(A, b, method='householder'): """ 使用指定的QR分解方法求解最小二乘问题 min ||Ax - b||_2 """ if method == 'householder': Q, R = qr_householder(A) Q = Q[:, :A.shape[1]] # 取经济型Q elif method == 'mgs': Q, R = qr_mgs(A) elif method == 'cgs': Q, R = qr_cgs(A) else: raise ValueError("方法必须是 'householder', 'mgs' 或 'cgs'") # 计算 c = Qᵀ b c = Q.T @ b # 回代求解 R β = c n = R.shape[1] beta = np.zeros(n) for i in range(n-1, -1, -1): # 从最后一行开始 beta[i] = (c[i] - np.dot(R[i, i+1:], beta[i+1:])) / R[i, i] return beta # 生成一些带噪声的线性数据 np.random.seed(123) n_samples = 50 x = np.linspace(0, 10, n_samples) true_beta = np.array([1.5, 0.8]) # 真实参数: 截距1.5，斜率0.8 y_true = true_beta[0] + true_beta[1] * x y_noisy = y_true + np.random.randn(n_samples) * 0.5 # 加入噪声 # 构造设计矩阵 A A_design = np.column_stack([np.ones_like(x), x]) b = y_noisy # 使用不同方法求解 beta_house = solve_least_squares_qr(A_design, b, 'householder') beta_mgs = solve_least_squares_qr(A_design, b, 'mgs') beta_cgs = solve_least_squares_qr(A_design, b, 'cgs') beta_numpy, residuals, rank, s = np.linalg.lstsq(A_design, b, rcond=None) print("=== 最小二乘拟合结果对比 ===") print(f"真实参数: {true_beta}") print(f"NumPy lstsq 求解: {beta_numpy}") print(f"Householder QR 求解: {beta_house}") print(f"MGS QR 求解: {beta_mgs}") print(f"CGS QR 求解: {beta_cgs}") print("\n=== 与NumPy结果的绝对误差 ===") print(f"Householder 误差: {np.abs(beta_numpy - beta_house)}") print(f"MGS 误差: {np.abs(beta_numpy - beta_mgs)}") print(f"CGS 误差: {np.abs(beta_numpy - beta_cgs)}") ``` 在这个例子中，由于问题本身是良态的，三种自实现算法可能都能给出与NumPy非常接近的结果。但你可以尝试修改数据，增加噪声或使设计矩阵 **A** 的列接近线性相关（例如，添加一个与现有列高度相关的列），这时CGS算法的解可能会明显偏离，而MGS和Householder则能保持稳健。 最后，我想分享一点在实现这些算法时容易踩的坑：**永远要注意浮点数的比较和除零问题**。在计算向量范数进行归一化前，检查其是否大于一个极小的阈值（如`1e-15`）。对于秩亏矩阵，需要有相应的处理逻辑，比如跳过零向量或引入列选主元（Column Pivoting）的QR分解，这能进一步稳定算法并处理秩不足的情况。这些细节，正是从“理解原理”到“实现可用代码”的关键跨越。