线性最小二乘问题的数值解法：QR分解与SVD比较

1. 线性最小二乘问题回顾与求解困境

在上一篇文章中，我们详细讨论了线性最小二乘问题的基本形式。给定一个线性模型f(x;θ)=Aθ，其中A∈R^{m×n}（m≥n），b∈R^m，最小二乘问题可以表示为：

min_θ ||Aθ - b||²

这等价于求解超定线性方程组Aθ=b。从理论上讲，当A列满秩时，其解析解为：

θ* = (A^T A)^(-1) A^T b

然而在实际工程应用中，直接使用这个解析解公式会遇到几个严重问题：

注意：在数值计算中直接求逆矩阵不仅效率低下，还会引入显著的舍入误差，这是需要特别警惕的。

1.1 逆矩阵的计算复杂度问题

考虑使用伴随矩阵法求逆矩阵的计算过程：

1. 计算行列式det(A)：需要O(n!)次运算
2. 计算余子式矩阵：每个元素需要计算一个(n-1)×(n-1)行列式
3. 构造伴随矩阵并进行转置
4. 最终进行标量乘法

即使采用相对高效的高斯消元法，复杂度也达到O(n³)。对于大规模问题（如n>1000），这样的计算成本是难以承受的。

1.2 数值稳定性问题

计算机使用浮点数表示实数时存在固有局限：

• 舍入误差：每次浮点运算都可能引入微小误差
• 误差累积：矩阵求逆涉及大量加减乘除运算，误差会不断累积
• 病态问题：当矩阵条件数较大时，微小扰动会导致解的巨大变化

2. 避免直接求逆的数值解法

2.1 QR分解法

2.1.1 QR分解原理

对矩阵A∈R^{m×n}（m≥n）进行QR分解：

A = Q₁R

其中：

• Q₁∈R^{m×n}列正交（Q₁^T Q₁ = I_n）
• R∈R^{n×n}为上三角矩阵

当A列满秩时，R可逆。将分解代入正规方程：

(A^T A)θ = A^T b
⇒ R^T Q₁^T Q₁ R θ = R^T Q₁^T b
⇒ R^T R θ = R^T (Q₁^T b)

两边左乘(R^T)^(-1)得到：

Rθ = Q₁^T b

这样就转化为一个上三角系统的求解问题。

2.1.2 实际计算步骤

1. 对A进行QR分解（常用Householder变换或Givens旋转）
2. 计算c = Q₁^T b
3. 回代求解Rθ = c

提示：现代数值线性代数库（如LAPACK）提供了高效的QR分解实现，通常比直接求逆快一个数量级。

2.2 SVD分解法

2.2.1 SVD分解原理

对任意矩阵A∈R^{m×n}，存在奇异值分解：

A = UΣV^T

其中：

• U∈R^{m×m}为正交矩阵（左奇异向量）
• V∈R^{n×n}为正交矩阵（右奇异向量）
• Σ∈R^{m×n}为对角矩阵（奇异值σ₁≥σ₂≥...≥σ_r>0）

将SVD代入正规方程：

A^T A θ = A^T b
⇒ VΣ^T U^T UΣV^T θ = VΣ^T U^T b
⇒ V(Σ^T Σ)V^T θ = VΣ^T U^T b

令y=V^Tθ，c=U^Tb，得到：

(Σ^T Σ)y = Σ^T c

这是一个对角方程组，易于求解。

2.2.2 实际计算步骤

1. 计算A的SVD分解
2. 确定有效秩r（忽略小于阈值的奇异值）
3. 计算c = U^T b
4. 求解y_i = c_i/σ_i (i=1,...,r)
5. 取y_j=0 (j=r+1,...,n)
6. 恢复解θ = V y

2.3 两种方法的比较

3. 分解方法的数学基础

3.1 QR分解的几何解释

QR分解本质上是Gram-Schmidt正交化过程的矩阵表示。给定矩阵A=[a₁,...,a_n]，通过以下步骤构造Q和R：

1. q₁ = a₁ / ||a₁||
2. 对于j=2,...,n：
   • v_j = a_j - Σ_{i=1}^{j-1}(q_i^T a_j)q_i
   • q_j = v_j / ||v_j||

这些关系可以紧凑地表示为A=QR，其中R的元素r_ij=q_i^T a_j。

3.2 SVD的代数基础

SVD可以看作对称矩阵特征值分解的推广。关键步骤：

1. 计算A^T A的特征值和特征向量：
   A^T A v_i = σ_i² v_i
2. 奇异值σ_i = √λ_i
3. 左奇异向量u_i = A v_i / σ_i
4. 构造U,Σ,V矩阵

4. 实际应用中的注意事项

4.1 阈值选择

在SVD求解中，确定有效秩时需要设置奇异值阈值。常用策略：

• 相对阈值：忽略σ_i < ε·σ₁
• 绝对阈值：基于机器精度和问题规模确定

4.2 病态问题处理

当条件数cond(A)=σ₁/σ_r很大时：

1. 使用截断SVD（TSVD）：忽略小的奇异值
2. 添加正则化项（Tikhonov正则化）：
   min ||Aθ-b||² + λ||θ||²

4.3 实现建议

1. 优先使用成熟的数值库（如LAPACK、Eigen等）
2. 对于大型稀疏矩阵，考虑迭代方法（如LSQR）
3. 定期检查残差||Aθ-b||以验证解的质量

5. 性能优化技巧

5.1 内存访问优化

• 对于大型矩阵，使用分块算法提高缓存利用率
• 考虑矩阵的存储顺序（行优先/列优先）

5.2 并行计算

• QR分解的Householder变换可并行化
• SVD的双对角化过程也有并行算法

5.3 特殊情况处理

• 对于Toeplitz等特殊结构矩阵，存在快速算法
• 增量式求解：当新增数据时，不必重新分解整个矩阵

6. 常见问题排查

6.1 解不稳定的可能原因

1. 矩阵接近秩亏（小的奇异值被保留）
2. 数据尺度差异大（建议对输入数据进行标准化）
3. 算法实现中的数值误差累积

6.2 诊断方法

1. 计算残差范数||Aθ-b||
2. 检查解的范数||θ||是否异常大
3. 绘制奇异值衰减曲线

6.3 解决方案

1. 增加正则化项
2. 使用更稳定的算法（如SVD代替QR）
3. 提高计算精度（如使用双精度）

在实际工程应用中，我通常会先尝试QR分解，因为它在大规模问题中效率较高。当遇到收敛问题或结果不稳定时，再转向SVD方法。对于特别大的问题，迭代法可能是唯一可行的选择。记住，没有放之四海而皆准的最佳方法，关键是根据具体问题的特性选择最适合的算法。