用Python实战QR分解：从数学公式到代码实现的三种算法对比

QR分解作为线性代数中的核心算法，在机器学习、数据科学和工程计算中扮演着重要角色。本文将带你用Python(NumPy/SciPy)实现施密特正交化、吉文斯旋转和豪斯霍尔德反射这三种QR分解算法，并通过可视化展示其工作原理。

1. QR分解基础与准备工作

QR分解是将矩阵A分解为正交矩阵Q和上三角矩阵R的过程，即A=QR。这种分解在求解线性方程组、最小二乘问题和特征值计算中非常有用。

首先设置Python环境并导入必要库：

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

创建一个用于演示的5×5随机矩阵：

np.random.seed(42)
A = np.random.randn(5, 5)
print("原始矩阵A:\n", A)

QR分解的核心数学性质包括：

• Q是正交矩阵(Q^TQ=I)
• R是上三角矩阵
• 对于可逆矩阵A，QR分解存在且当R对角元为正时唯一

2. 施密特正交化实现

施密特正交化是最直观的QR分解方法，通过对矩阵列向量逐步正交化实现。

算法步骤 ：

1. 取矩阵A的第一列作为q₁的初始向量并归一化
2. 对后续每一列，减去在前面的所有q向量上的投影
3. 对结果向量归一化得到新的q向量

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] = np.dot(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

验证分解结果：

Q, R = gram_schmidt_qr(A)
print("Q的正交性检验:\n", np.round(Q.T @ Q, 10))
print("重构误差:", np.linalg.norm(A - Q @ R))

注意：经典施密特正交化在数值计算中可能不稳定，改进版本会使用额外的重新正交化步骤。

3. 吉文斯旋转算法实现

吉文斯旋转通过一系列平面旋转将矩阵转化为上三角形式，特别适合稀疏矩阵。

旋转矩阵原理 ： 对于要消去的元素a[j,i]，构造旋转矩阵G(i,j,θ)使得：

c = a[i,i]/sqrt(a[i,i]² + a[j,i]²)
s = a[j,i]/sqrt(a[i,i]² + a[j,i]²)

Python实现代码：

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

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

可视化旋转过程：

def visualize_givens():
    fig, ax = plt.subplots()
    x = np.array([3, 1])
    line, = ax.plot([0, x[0]], [0, x[1]], 'b-', lw=2)
    
    def update(theta):
        c, s = np.cos(theta), np.sin(theta)
        G = np.array([[c, s], [-s, c]])
        new_x = G @ x
        line.set_data([0, new_x[0]], [0, new_x[1]])
        return line,
    
    ani = FuncAnimation(fig, update, frames=np.linspace(0, np.pi/2, 100),
                        interval=50, blit=True)
    plt.show()

4. 豪斯霍尔德反射算法实现

豪斯霍尔德反射使用镜面反射将矩阵列向量映射到坐标轴方向，是数值计算中最稳定的QR分解方法。

反射向量计算 ： 对于向量x，反射向量v和对应的Householder矩阵H为：

v = x - sign(x₁)||x||e₁
H = I - 2vvᵀ/(vᵀv)

Python实现代码：

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).reshape(-1, 1)
        
        if np.linalg.norm(v) > 1e-10:
            v = v / np.linalg.norm(v)
            H = np.eye(m)
            H[k:, k:] -= 2 * v @ v.T
            R = H @ R
            Q = Q @ H.T
    
    return Q[:, :n], R[:n]

性能对比实验：

import time

methods = {
    "Gram-Schmidt": gram_schmidt_qr,
    "Givens": givens_qr,
    "Householder": householder_qr
}

sizes = [10, 50, 100, 200]
results = {name: [] for name in methods}

for size in sizes:
    A = np.random.randn(size, size)
    for name, method in methods.items():
        start = time.time()
        Q, R = method(A)
        elapsed = time.time() - start
        error = np.linalg.norm(A - Q @ R)
        results[name].append((size, elapsed, error))

5. 三种算法对比与应用建议

数值稳定性对比 ：

• 施密特正交化：稳定性最差，特别是当矩阵接近秩亏损时
• 吉文斯旋转：稳定性好，适合稀疏矩阵
• 豪斯霍尔德反射：数值稳定性最好，是LAPACK等库的默认选择

计算复杂度比较 ：

算法	浮点运算次数	适用场景
施密特	≈2mn²	教学理解
吉文斯	≈3mn² - n³	稀疏矩阵
豪斯霍尔德	≈2mn² - 2n³/3	通用场景

实际应用建议 ：

1. 对于小规模或教学目的，施密特正交化最直观
2. 处理稀疏矩阵或需要并行化时，吉文斯旋转更优
3. 在大多数数值计算场景中，豪斯霍尔德反射是首选

SciPy中的优化实现：

# LAPACK中的高效实现
Q, R = la.qr(A, mode='economic')  # 经济型QR分解

QR分解在最小二乘问题中的应用示例：

def least_squares_qr(A, b):
    Q, R = la.qr(A, mode='economic')
    y = Q.T @ b
    return la.solve_triangular(R[:A.shape[1]], y[:A.shape[1]])

在实现这些算法时，我发现在处理病态矩阵时，豪斯霍尔德方法的鲁棒性确实明显优于其他两种方法。特别是在处理接近线性相关的列向量时，经典施密特正交化会产生显著的正交性误差，而豪斯霍尔德方法仍能保持较好的数值稳定性。