用Python+NumPy实战QR分解：从施密特正交化到矩阵可视化

线性代数中那些抽象的概念，是否总让你感到困惑？今天我们不谈枯燥的数学推导，而是用Python代码和可视化手段，带你亲手实现QR分解的核心算法——施密特正交化及其改进版本。通过NumPy的矩阵运算和Matplotlib的动态演示，你将直观理解正交矩阵的构造过程、浮点运算误差的影响，以及为什么改进算法能提供更稳定的数值结果。这种"代码先行"的学习方式，特别适合那些喜欢通过实践来掌握理论的技术爱好者。

1. 环境准备与基础概念

在开始编码前，我们需要配置好Python环境并理解几个核心概念。推荐使用Anaconda创建独立环境，确保库版本一致：

# 创建并激活conda环境
conda create -n linear_algebra python=3.8
conda activate linear_algebra

# 安装必要库
pip install numpy matplotlib ipython

QR分解的核心是将矩阵A分解为Q(正交矩阵)和R(上三角矩阵)的乘积。理解这个分解的关键在于掌握正交矩阵的特性：

• 正交矩阵：列向量两两正交且长度为1，满足QᵀQ = I
• 上三角矩阵：主对角线下方元素全为零
• 施密特正交化：将线性无关向量组转化为正交向量组的过程

让我们先看一个简单的2D例子。假设有两个向量：

import numpy as np

v1 = np.array([1, 1])
v2 = np.array([1, 0])

这两个向量不正交，我们需要通过施密特过程找到一组正交基。

2. 经典施密特正交化实现

经典施密特算法的步骤清晰明了：

1. 第一个基向量直接归一化
2. 后续每个向量减去它在前面所有基向量上的投影
3. 将结果归一化得到新的基向量

用Python实现这个逻辑：

def classical_gram_schmidt(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

让我们测试这个实现：

A = np.array([[1, 1], [1, 0]], dtype=float)
Q, R = classical_gram_schmidt(A)
print("正交矩阵Q:\n", Q)
print("上三角矩阵R:\n", R)
print("验证QᵀQ:\n", Q.T @ Q)

输出应该显示Q的列向量正交，且QᵀQ接近单位矩阵。但当我们尝试接近线性相关的向量时，问题就出现了：

A = np.array([[1, 1.0001], [1, 1]], dtype=float)
Q, R = classical_gram_schmidt(A)
print("验证正交性:\n", Q.T @ Q)

你会发现非对角线元素不再接近零，这就是浮点误差累积的结果。

3. 改进施密特正交化算法

改进算法通过即时更新剩余向量来减少误差传播。具体实现：

def modified_gram_schmidt(A):
    """改进施密特正交化"""
    m, n = A.shape
    Q = np.zeros((m, n))
    R = np.zeros((n, n))
    V = A.copy()
    
    for i in range(n):
        R[i, i] = np.linalg.norm(V[:, i])
        Q[:, i] = V[:, i] / R[i, 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

对比两种算法的数值稳定性：

np.random.seed(42)
A = np.random.randn(10, 10)  # 随机生成10x10矩阵

# 经典算法
Q1, R1 = classical_gram_schmidt(A)
error1 = np.linalg.norm(Q1.T @ Q1 - np.eye(10))

# 改进算法
Q2, R2 = modified_gram_schmidt(A)
error2 = np.linalg.norm(Q2.T @ Q2 - np.eye(10))

print(f"经典算法正交性误差: {error1:.4e}")
print(f"改进算法正交性误差: {error2:.4e}")

改进算法通常能降低1-2个数量级的误差。为了更直观理解，我们可以可视化正交化过程：

import matplotlib.pyplot as plt

def plot_vectors(vectors, colors, title):
    plt.figure()
    ax = plt.gca()
    for v, c in zip(vectors, colors):
        ax.quiver(0, 0, v[0], v[1], angles='xy', scale_units='xy', scale=1, color=c)
    plt.xlim(-2, 2)
    plt.ylim(-2, 2)
    plt.grid()
    plt.title(title)
    plt.show()

# 原始向量
v1 = np.array([1, 0.5])
v2 = np.array([1, 1])
plot_vectors([v1, v2], ['r', 'b'], '原始向量')

# 经典施密特过程
u1 = v1 / np.linalg.norm(v1)
v2_ortho = v2 - np.dot(v2, u1) * u1
plot_vectors([u1, v2_ortho], ['g', 'm'], '经典施密特正交化结果')

# 改进算法过程
u1_m = v1 / np.linalg.norm(v1)
v2_m = v2 - np.dot(v2, u1_m) * u1_m
u2_m = v2_m / np.linalg.norm(v2_m)
plot_vectors([u1_m, u2_m], ['c', 'y'], '改进算法结果')

4. 完整QR分解与应用实例

现在我们将改进算法应用于实际问题的求解。考虑超定方程组Ax=b的最小二乘解：

# 生成测试数据
np.random.seed(42)
A = np.random.randn(50, 10)  # 50个方程，10个未知数
x_true = np.random.randn(10)
b = A @ x_true + 0.1 * np.random.randn(50)  # 添加噪声

# 使用QR分解求解
Q, R = modified_gram_schmidt(A)
x_qr = np.linalg.solve(R, Q.T @ b)

# 与正规方程解比较
x_normal = np.linalg.solve(A.T @ A, A.T @ b)

print("QR解与真实值误差:", np.linalg.norm(x_qr - x_true))
print("正规方程解与真实值误差:", np.linalg.norm(x_normal - x_true))

QR分解通常能提供更稳定的数值解，特别是当A条件数较大时。我们可以进一步分析不同算法的数值表现：

算法类型	正交性误差	计算复杂度	数值稳定性
经典施密特	较高	O(mn²)	较差
改进施密特	中等	O(mn²)	较好
Householder	低	O(mn²)	优秀
Givens旋转	低	O(mn²)	优秀

虽然我们主要讨论了施密特方法，但在实际应用中，Householder变换通常是更优的选择。NumPy的np.linalg.qr函数默认使用Householder算法：

Q_np, R_np = np.linalg.qr(A, mode='reduced')
error_np = np.linalg.norm(Q_np.T @ Q_np - np.eye(10))
print(f"NumPy QR正交性误差: {error_np:.4e}")

理解这些底层算法后，你就能根据具体问题选择合适的分解方法。比如在稀疏矩阵场景，Givens旋转可能更高效；而当需要逐步更新分解时，施密特方法可能更合适。