别再死记硬背公式了！用Python（NumPy/SciPy）手撕QR分解的三种算法（附代码对比）

weixin_33728268

416人浏览 · 2026-06-03 14:12:11

weixin_33728268 · 2026-06-03 14:12:11 发布

用Python实战QR分解：从数学公式到高效代码的三种实现路径

QR分解作为线性代数中的核心算法，在机器学习、数据科学和工程计算中扮演着关键角色。但大多数教程仅停留在理论推导，缺乏实际编码指导。本文将打破这一常规，带你用NumPy和SciPy实现三种QR分解算法，并通过性能对比揭示它们在不同场景下的优劣。

1. QR分解为何值得你掌握？

当你处理最小二乘问题、特征值计算或线性方程组求解时，QR分解就像瑞士军刀般多功能。它将任意矩阵A分解为正交矩阵Q和上三角矩阵R的乘积，这种结构为数值计算带来了稳定性优势。

传统教学中，我们常被要求死记硬背施密特正交化的公式，却很少思考：为什么需要多种QR算法？每种方法在计算机中实际表现如何？让我们用Python解开这些疑问。

先看一个典型应用场景：求解超定方程组Ax=b的最小二乘解。QR分解解法比直接计算(A^T A)^-1 A^T b更稳定：

import numpy as np
from scipy.linalg import qr

# 超定方程组示例
A = np.random.rand(100, 80)  # 100个方程，80个未知数
b = np.random.rand(100)

# 传统法 (数值不稳定)
x_normal = np.linalg.inv(A.T @ A) @ A.T @ b

# QR分解法
Q, R = qr(A, mode='economic')
x_qr = np.linalg.solve(R, Q.T @ b)

print(f"传统方法误差: {np.linalg.norm(A @ x_normal - b):.3e}")
print(f"QR方法误差: {np.linalg.norm(A @ x_qr - b):.3e}")

2. 施密特正交化：最直观的QR实现

施密特正交化是教科书最爱介绍的方法，其数学过程直观：

将矩阵A的列向量α₁,α₂,...αₙ正交化为β₁,β₂,...βₙ
单位化得到Q的列向量
通过变换矩阵B求得R

Python实现揭示了其计算效率问题：

def gram_schmidt_qr(A):
    m, n = A.shape
    Q = np.zeros((m, n))
    R = np.zeros((n, n))
    
    for j in range(n):
        v = A[:, j]
        for i in range(j):
            R[i, j] = Q[:, i].T @ 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

# 测试病态矩阵
A = np.array([[1, 1, 1],
              [1e-10, 0, 0],
              [0, 1e-10, 0],
              [0, 0, 1e-10]])
Q, R = gram_schmidt_qr(A)
print("施密特法的正交性检验:")
print(Q.T @ Q)  # 应接近单位矩阵

注意：当矩阵列向量接近线性相关时，经典施密特法会丧失正交性。改进版可通过一次计算所有投影来提升稳定性。

3. 吉文斯旋转：适合稀疏矩阵的优雅方案

吉文斯旋转通过一系列平面旋转逐步引入零元素，特别适合稀疏矩阵。其核心思想是：每次消去一个非对角线元素，而不影响已处理的零元素。

def givens_rotation(a, b):
    if b == 0:
        c, s = 1, 0
    else:
        if abs(b) > abs(a):
            r = a / b
            s = 1 / np.sqrt(1 + r**2)
            c = s * r
        else:
            r = b / a
            c = 1 / np.sqrt(1 + r**2)
            s = c * r
    return c, s

def givens_qr(A):
    m, n = A.shape
    Q = np.eye(m)
    R = A.copy()
    
    for j in range(n):
        for i in range(m-1, j, -1):
            if R[i, j] != 0:
                c, s = givens_rotation(R[i-1, j], R[i, j])
                G = np.eye(m)
                G[[i-1, i], [i-1, i]] = c
                G[i-1, i], G[i, i-1] = s, -s
                R = G.T @ R
                Q = Q @ G
    return Q[:n].T, R[:n]

# 创建三对角矩阵示例
n = 5
A = np.diag(np.full(n, 2)) + np.diag(np.full(n-1, -1), 1) + np.diag(np.full(n-1, -1), -1)
Q, R = givens_qr(A)
print("\n吉文斯法的Q正交性:")
print(np.round(Q.T @ Q, 10))

性能对比表：

方法	时间复杂度	适合矩阵类型	数值稳定性
施密特正交化	O(mn²)	稠密	中等
吉文斯旋转	O(n³)	稀疏	高
豪斯霍尔德	O(2mn²-2n³/3)	通用	最高

4. 豪斯霍尔德反射：LAPACK采用的工业标准

豪斯霍尔德反射是专业数值库的首选，它通过镜面反射实现正交化，具有最优的数值稳定性。SciPy的qr函数默认采用此方法：

def householder_qr(A):
    m, n = A.shape
    R = A.copy()
    Q = np.eye(m)
    
    for k in range(n):
        x = R[k:, k]
        e = np.zeros_like(x)
        e[0] = np.sign(x[0]) * np.linalg.norm(x)
        v = (x - e) / np.linalg.norm(x - e)
        
        H = np.eye(m)
        H[k:, k:] -= 2 * np.outer(v, v)
        
        R = H @ R
        Q = Q @ H.T
    return Q, R

# 性能对比
A_large = np.random.rand(500, 400)
%timeit gram_schmidt_qr(A_large)  # 约1.2秒
%timeit givens_qr(A_large)       # 约8.7秒 
%timeit householder_qr(A_large)  # 约0.3秒

豪斯霍尔德的秘诀在于：

每次处理一列，通过反射将其下方元素归零
反射矩阵的特殊结构使得计算可以优化
避免了施密特法的正交性累积误差

5. 算法实战：如何根据场景选择最佳方法？

三种方法各有胜负：

教学演示 ：施密特法最直观
稀疏矩阵 ：吉文斯旋转最节省计算
生产环境 ：豪斯霍尔德反射是稳妥选择

病态矩阵测试揭示关键差异：

# 构建希尔伯特矩阵(经典病态矩阵)
n = 10
H = np.array([[1/(i+j+1) for j in range(n)] for i in range(n)])

methods = {
    "施密特": gram_schmidt_qr,
    "吉文斯": givens_qr,
    "豪斯霍尔德": householder_qr
}

for name, method in methods.items():
    Q, R = method(H)
    orth_error = np.linalg.norm(Q.T @ Q - np.eye(n))
    res_error = np.linalg.norm(Q @ R - H)
    print(f"{name}法 - 正交误差: {orth_error:.3e}, 分解误差: {res_error:.3e}")

输出结果通常显示豪斯霍尔德在病态情况下表现最优。对于特定场景，混合策略可能更佳——例如对稀疏部分用吉文斯，稠密部分用豪斯霍尔德。

在实现细节上，专业库会做更多优化：