1. QR分解的本质与数学定义

QR分解是将任意m×n矩阵A分解为一个正交矩阵Q和一个上三角矩阵R的乘积的过程。这种分解在数值线性代数中占据核心地位，因为它不仅提供了矩阵的结构化表示，更为诸多数值计算问题提供了稳定高效的解法。

从数学角度看，对于实矩阵A，其QR分解满足： A = QR 其中：

• Q是m×m正交矩阵（QᵀQ = I）
• R是m×n上三角矩阵（当m>n时，底部为全零行）

当矩阵A列满秩时，可以进一步得到"经济型"分解： A = Q₁R₁ 其中Q₁是m×n矩阵，R₁是n×n上三角矩阵

2. QR分解的三大经典算法实现

2.1 Householder变换法

Householder反射是QR分解最稳定的数值算法。通过构造一系列反射矩阵Hₖ，逐步将A的下三角部分归零：

1. 对第k列向量aₖ，计算反射向量： v = aₖ - sign(a₁ₖ)||aₖ||e₁

2. 构造Householder矩阵： H = I - 2vvᵀ/(vᵀv)

3. 应用变换： A ← HₖA

经过n步（m≥n时）后，A被转化为上三角矩阵R，而Q = H₁H₂...Hₙ

关键技巧：为避免数值误差，通常取sign(a₁ₖ)与a₁ₖ相反，并采用向量化的块操作提升性能

2.2 Givens旋转法

适用于稀疏矩阵或需要选择性归零的场景。通过平面旋转逐步消元：

1. 对元素a_{ij}（i>j），计算旋转角度： c = a_{jj}/r, s = -a_{ij}/r r = √(a_{jj}² + a_{ij}²)

2. 构造Givens矩阵G(i,j,θ)

3. 应用旋转： A ← GA

该方法需要约2n³/3次浮点运算，比Householder变换多50%计算量

2.3 Gram-Schmidt正交化

最直观但数值稳定性最差的方法。分为经典GS和改进GS两种形式：

经典GS过程： for j = 1:n v = A(:,j) for i = 1:j-1 v = v - (Q(:,i)'*A(:,j))*Q(:,i) end Q(:,j) = v/||v|| end

改进GS通过调整计算顺序提升稳定性： for j = 1:n v = A(:,j) for i = 1:j-1 R(i,j) = Q(:,i)'*v v = v - R(i,j)*Q(:,i) end R(j,j) = ||v|| Q(:,j) = v/R(j,j) end

3. QR分解的数值特性分析

3.1 稳定性对比

• Householder变换：向后稳定，误差界约O(mn²)ε
• Givens旋转：同Householder相当，但实现更复杂
• Gram-Schmidt：改进GS可达O(mn²)ε，经典GS可能O(mn²)κ(A)ε

3.2 计算复杂度

4. QR分解的典型应用场景

4.1 线性最小二乘问题

求解min||Ax-b||₂时，正规方程AᵀAx=Aᵀb的条件数为κ(A)²，而QR分解直接得到： x = R₁⁻¹(Q₁ᵀb) 条件数仅为κ(A)

4.2 特征值计算（QR算法）

通过迭代过程： Aₖ = QₖRₖ Aₖ₊₁ = RₖQₖ 最终收敛到Schur标准形

4.3 矩阵秩揭示

R矩阵的非零对角线元素个数即为矩阵秩，比SVD计算更高效

5. 实际编程实现要点

5.1 Householder实现技巧

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

5.2 数值稳定性增强措施

1. 列主元策略：每次选择范数最大的列优先处理
2. 迭代精化：对残差进行二次求解
3. 条件数估计：通过R矩阵判断问题适定性

6. 常见问题排查指南

7. 性能优化实践

对于大规模矩阵：

1. 分块算法：将矩阵划分为适合缓存的大小（通常64×64到256×256）
2. 并行化：利用Householder向量的树形归约
3. GPU加速：特别适合Givens旋转的批处理

实测数据表明，在双精度浮点下：

• 1000×1000矩阵：CPU约0.5秒，GPU加速后可达0.05秒
• 稀疏矩阵（密度5%）：专用算法可提速10倍以上