用Python实现QR分解：从数学公式到可视化编程实践

线性代数中的QR分解是机器学习、数据科学和工程计算中的基础工具，但很多学习者却被抽象的数学符号和理论推导劝退。本文将带你用Python的NumPy和Matplotlib库，通过代码实现三种经典QR分解算法，并可视化每一步的计算过程，让抽象的矩阵运算变得直观可见。

1. QR分解为何重要：从理论到实践

QR分解是将一个矩阵A分解为正交矩阵Q和上三角矩阵R的乘积（A=QR）。这种分解在数值计算中有着广泛应用：

• 最小二乘问题 ：求解Ax=b时，当A不是方阵或条件数较大时，QR分解比直接求逆更稳定
• 特征值计算 ：QR算法是计算矩阵特征值的标准方法
• 机器学习 ：在主成分分析(PCA)和线性回归中都有重要应用

传统教学中，QR分解通常以数学公式形式呈现，例如施密特正交化的投影公式：

proj_u(v) = (u·v)/(u·u) * u

但对于编程实践者来说，更关心的是如何用代码实现这些数学概念。下面我们将用Python实现三种QR分解算法，并比较它们的特性。

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

施密特正交化是最容易理解的QR分解方法，它通过逐步正交化矩阵的列向量来构造Q矩阵。

import numpy as np

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

可视化关键步骤 ：我们可以绘制向量在正交化过程中的变化：

import matplotlib.pyplot as plt

A = np.array([[1, 1], [1, 0]], dtype=float)
Q, R = gram_schmidt_qr(A)

plt.quiver(0, 0, A[0,0], A[1,0], angles='xy', scale_units='xy', scale=1, color='r')
plt.quiver(0, 0, A[0,1], A[1,1], angles='xy', scale_units='xy', scale=1, color='b')
plt.quiver(0, 0, Q[0,1], Q[1,1], angles='xy', scale_units='xy', scale=1, color='g')
plt.xlim(-1, 2)
plt.ylim(-1, 2)
plt.grid()
plt.show()

施密特正交化的数值稳定性较差，当矩阵条件数较大时，正交化的向量可能会逐渐失去正交性。在实际应用中，通常使用下面两种更稳定的方法。

3. 吉文斯旋转：逐个元素消去的艺术

吉文斯旋转通过一系列平面旋转将矩阵转化为上三角形式。每个旋转只影响矩阵的两行，将一个特定元素变为零。

def givens_rotation(a, b):
    if b == 0:
        c = 1
        s = 0
    else:
        if abs(b) > abs(a):
            t = -a / b
            s = 1 / np.sqrt(1 + t**2)
            c = s * t
        else:
            t = -b / a
            c = 1 / np.sqrt(1 + t**2)
            s = c * t
    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+1, i-1:i+1] = np.array([[c, -s], [s, c]])
                R = G.T @ R
                Q = Q @ G
                
    return Q[:n].T, R[:n]

可视化旋转过程 ：我们可以绘制每个吉文斯旋转对矩阵的影响：

# 创建一个动画展示吉文斯旋转过程
from matplotlib.animation import FuncAnimation

A = np.random.rand(3, 3)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

def update(frame):
    ax.clear()
    # 这里简化展示，实际应跟踪每一步旋转
    Q, R = givens_qr(A)
    ax.quiver(0, 0, 0, A[0,0], A[1,0], A[2,0], color='r')
    ax.quiver(0, 0, 0, Q[0,0], Q[1,0], Q[2,0], color='g')
    ax.set_xlim(-1, 1)
    ax.set_ylim(-1, 1)
    ax.set_zlim(-1, 1)

ani = FuncAnimation(fig, update, frames=10, interval=500)
plt.close()

吉文斯旋转特别适合稀疏矩阵，因为它可以只对非零元素进行操作。但在处理稠密矩阵时，豪斯霍尔德变换通常更高效。

4. 豪斯霍尔德变换：镜像反射的威力

豪斯霍尔德变换使用镜像反射将矩阵的列向量逐步转化为上三角形式。每个变换可以将一个向量的多个元素同时置零。

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

稳定性分析 ：豪斯霍尔德变换是三种方法中数值稳定性最好的，特别适合处理病态矩阵。我们可以通过条件数来比较三种方法的稳定性：

np.random.seed(42)
A = np.random.rand(5, 5) * 1e-8
A[:, 0] = A[:, 1]  # 制造病态条件

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

for name, method in methods.items():
    Q, R = method(A)
    orthogonality_error = np.linalg.norm(Q.T @ Q - np.eye(5))
    reconstruction_error = np.linalg.norm(Q @ R - A)
    print(f"{name}: Orthogonality error={orthogonality_error:.2e}, Reconstruction error={reconstruction_error:.2e}")

输出结果可能显示豪斯霍尔德的误差最小，验证了其数值稳定性。

5. 三种方法的比较与选择指南

在实际应用中，我们需要根据矩阵特性和计算需求选择合适的方法。以下是三种方法的对比：

特性	施密特正交化	吉文斯旋转	豪斯霍尔德变换
数值稳定性	较差	中等	最好
计算复杂度	O(mn²)	O(mn²)～O(mn³)	O(mn²)
并行化难度	中等	较难	较易
适合矩阵类型	小规模、良态矩阵	稀疏矩阵	稠密矩阵
实现复杂度	最简单	中等	较复杂

何时选择哪种方法 ：

• 教学和理解QR分解原理：施密特正交化
• 处理稀疏矩阵或并行计算：吉文斯旋转
• 一般数值计算和稳定性要求高：豪斯霍尔德变换

对于Python用户，NumPy和SciPy已经提供了优化的QR分解实现：

import numpy as np
from scipy.linalg import qr

A = np.random.rand(5, 3)
Q, R = qr(A)  # 默认使用豪斯霍尔德方法
Q_gs, R_gs = qr(A, mode='economic')  # 经济型QR分解

理解这些底层实现不仅能帮助我们更好地使用库函数，还能在需要自定义分解时提供思路。比如，在处理特殊结构矩阵时，我们可能需要结合这些基本算法的变种。

通过代码实现和可视化，抽象的QR分解变得具体而直观。这种"手撸"算法的练习是理解线性代数计算本质的最佳途径，远比死记硬背公式有效。当你下次遇到QR分解的应用场景时，不妨回想这些可视化图像，它们会帮助你更深入地理解矩阵运算的几何意义。