Tim Shen1. 线性最小二乘问题概述
在工程计算和数据分析中，我们经常会遇到需要求解超定线性方程组的情况。这类问题可以表述为：给定一个m×n的矩阵A（m>n）和m维向量b，寻找n维向量θ使得Aθ尽可能接近b。由于方程数量多于未知数，通常没有精确解，因此我们转而寻求最小二乘解。
最小二乘问题的数学表述是：
min_θ ||Aθ - b||²
这个优化问题的解可以通过求解正规方程获得：
AᵀAθ = Aᵀb
理论上，当A列满秩时，解可以表示为：
θ* = (AᵀA)⁻¹Aᵀb
然而在实际数值计算中，直接计算逆矩阵(AᵀA)⁻¹既低效又不稳定。下面我将详细介绍三种更优的求解方法，并分析它们各自的适用场景。
需要模型API调用？ 免费领10W Token，多模型网关一键接入 Claude、DeepSeek 等主流模型。
2. 三种数值求解方法比较
2.1 直接求逆法的问题
在大学线性代数课程中，我们学过两种求逆方法：伴随矩阵法和高斯消元法。让我们分析它们的计算复杂度：
- 伴随矩阵法：
- 计算每个元素的余子式：需要计算(n-1)×(n-1)矩阵的行列式
- 行列式计算复杂度为O(n!)
- 总复杂度约为O(n·n!)
- 高斯消元法：
- 前向消元：O(n³)
- 回代求解：O(n³)
- 总复杂度为O(n³)
即使使用相对高效的高斯消元法，对于大型矩阵（如n=1000），n³=10⁹次运算仍然非常耗时。此外，计算机浮点运算带来的舍入误差会在求逆过程中不断累积，导致结果不精确。
提示：在实际工程中，当矩阵维度n>100时，就应当避免直接求逆，转而使用更稳定的数值方法。
2.2 QR分解法
QR分解将矩阵A分解为正交矩阵Q和上三角矩阵R的乘积：
A = QR
求解过程如下：
- 对A进行QR分解（常用Householder变换或Givens旋转）
- 正规方程简化为：Rθ = Qᵀb
- 回代求解上三角方程组
QR分解的优势：
- 计算复杂度：O(mn²)，比直接求逆更高效
- 数值稳定性好，适合条件数不是特别大的矩阵
- 当A列满秩时，能给出唯一解
我在实际项目中发现，对于中等规模（m<10⁴，n<10³）且良态的问题，QR分解通常是首选方法。一个典型的应用场景是传感器标定，其中我们需要处理数千个测量方程和几十个待标定参数。
2.3 SVD分解法
奇异值分解（SVD）将任意矩阵A分解为：
A = UΣVᵀ
求解最小二乘问题的步骤：
- 计算A的SVD分解
- 构建Σ的伪逆Σ⁺（非零奇异值取倒数）
- 解为：θ = VΣ⁺Uᵀb
SVD的特点：
- 能处理秩亏矩阵（A不满秩）
- 自动给出最小范数解
- 计算复杂度较高：O(mn²)到O(mn²+n³)
- 数值稳定性最好
在图像处理和信号处理领域，我经常遇到病态问题（如图像去模糊），此时SVD是必不可少的工具。通过观察奇异值的衰减情况，还可以判断问题的病态程度并决定截断策略。
3. 方法选择指南
根据我的工程实践经验，这三种方法的选择应考虑以下因素：
| 方法 | 适用场景 | 优点 | 缺点 |
|---|---|---|---|
| 直接求逆 | 小规模(n<100)理论分析 | 概念简单 | 效率低、不稳定 |
| QR分解 | 中大规模列满秩问题 | 效率较高稳定性好 | 不能处理秩亏 |
| SVD分解 | 病态或秩亏问题 | 最稳定处理秩亏 | 计算量最大 |
一个实用的建议是：先尝试QR分解，如果发现矩阵条件数很大（cond(A)>1e10）或明显秩亏，再转向SVD。在MATLAB或Python中，可以使用cond()函数估计矩阵条件数。
4. 实现细节与注意事项
4.1 QR分解的实现
在实际编程中，我推荐使用以下库函数：
Python (NumPy/SciPy):
python复制import scipy.linalg
# Householder QR
Q, R = scipy.linalg.qr(A, mode='economic')
# 求解最小二乘问题
theta = scipy.linalg.solve_triangular(R, Q.T @ b)
MATLAB:
matlab复制[Q,R] = qr(A,0); % 经济型QR分解
theta = R\(Q'*b); % 回代求解
注意事项：
- 使用"经济型"QR分解以避免计算不必要的零元素
- 确保使用稳定的QR算法（如Householder变换）
- 解三角方程组时优先使用专用函数（如solve_triangular）
4.2 SVD分解的实现
Python实现示例：
python复制import numpy as np
U, s, Vh = np.linalg.svd(A, full_matrices=False)
theta = Vh.T @ np.diag(1/s) @ U.T @ b
MATLAB实现：
matlab复制[U,S,V] = svd(A,'econ');
theta = V*diag(1./diag(S))*(U'*b);
使用技巧：
- 设置full_matrices=False或'econ'以提高计算效率
- 对于病态问题，可以截断小奇异值（Tikhonov正则化）
- 大矩阵可考虑使用随机SVD等近似算法
5. 常见问题与解决方案
5.1 数值不稳定问题
症状：解对数据微小扰动非常敏感，结果波动大
解决方法：
- 改用SVD分解
- 添加正则化项：min ||Aθ-b||² + λ||θ||²
- 检查数据是否需要标准化
5.2 秩亏问题
症状：矩阵条件数极大或计算中出现奇异警告
解决方法：
- 使用SVD获取最小范数解
- 分析问题物理意义，可能需添加约束条件
- 考虑变量之间存在精确线性关系
5.3 大规模问题
症状：内存不足或计算时间过长
解决方法：
- 使用迭代法（如LSQR）
- 考虑矩阵稀疏性，使用稀疏矩阵库
- 分布式计算（如Spark MLlib）
6. 实际案例分享
最近在一个机器人定位项目中，我们需要处理约10,000个测距方程来估计100个位置参数。最初尝试直接求逆导致数值不稳定，后改用QR分解成功解决。关键代码片段：
python复制# 数据预处理：标准化
A = (A - np.mean(A, axis=0)) / np.std(A, axis=0)
b = b - np.mean(b)
# 带列主元的QR分解
Q, R, piv = scipy.linalg.qr(A, mode='economic', pivoting=True)
theta_piv = np.zeros(A.shape[1])
theta_piv[piv] = scipy.linalg.solve_triangular(R, Q.T @ b)
# 结果后处理
theta = theta_piv / np.std(A, axis=0)
这个案例教会我：数据预处理（标准化）和适当的数值方法选择同样重要。通过列主元QR分解，我们不仅获得了稳定解，还识别出了对解影响最大的关键测量。