从理论到实践：线性与非线性最小二乘的算法演进与工程实现-尧图网络科技

1. 最小二乘法的数学本质

最小二乘法本质上是一种数学优化技术，它的核心思想是通过最小化误差平方和来寻找数据的最佳匹配函数。想象一下这样的场景：你有一堆散乱的数据点，想找一条直线或曲线尽可能靠近所有这些点。这就是最小二乘法要解决的问题。

从概率角度看，最小二乘解恰好对应着极大似然估计。假设误差服从正态分布，那么最大化观测数据的概率等价于最小化误差平方和。这种概率解释赋予了最小二乘法坚实的统计基础。

在实际应用中，我们通常会遇到两种形式的最小二乘问题：

线性最小二乘：模型关于参数是线性的，可以表示为矩阵方程
非线性最小二乘：模型关于参数是非线性的，需要迭代求解

这两种形式在工程实践中都非常常见。比如在传感器标定中常用线性最小二乘，而在SLAM（同步定位与地图构建）等复杂系统中则更多使用非线性最小二乘。

2. 线性最小二乘的求解方法

2.1 正规方程法

对于线性最小二乘问题 min||Ax-b||²，最直接的解法是求解正规方程：

AᵀAx = Aᵀb

这个方法的优点是理论简单，直接给出解析解。但在实际应用中存在几个问题：

计算AᵀA可能导致数值不稳定
当A条件数很大时，解可能不准确
对于大规模问题，矩阵求逆计算量太大

2.2 QR分解法

QR分解将矩阵A分解为正交矩阵Q和上三角矩阵R的乘积。利用正交矩阵的性质，最小二乘问题可以转化为：

Rx = Qᵀb

由于R是上三角矩阵，这个方程可以通过回代法高效求解。QR分解数值稳定性好，是求解中小规模线性最小二乘问题的首选方法。

在Eigen库中的实现示例：

#include <Eigen/Dense> using namespace Eigen; VectorXd qrSolver(const MatrixXd& A, const VectorXd& b) { HouseholderQR<MatrixXd> qr(A); return qr.solve(b); }

2.3 SVD分解法

奇异值分解(SVD)将矩阵A分解为UΣVᵀ，其中U和V是正交矩阵，Σ是对角矩阵。最小二乘解可以通过下式求得：

x = VΣ⁺Uᵀb

其中Σ⁺是Σ的伪逆。SVD方法最稳定可靠，特别适合处理秩亏或接近秩亏的矩阵，但计算成本也最高。

Eigen中的SVD求解示例：

VectorXd svdSolver(const MatrixXd& A, const VectorXd& b) { JacobiSVD<MatrixXd> svd(A, ComputeThinU | ComputeThinV); return svd.solve(b); }

3. 非线性最小二乘的迭代算法

3.1 高斯-牛顿法

高斯-牛顿法是非线性最小二乘最常用的求解方法。它通过在当前点对非线性函数进行一阶泰勒展开，将非线性问题转化为一系列线性最小二乘问题。

算法步骤：

初始化参数x₀
在每一步迭代中：
- 计算雅可比矩阵J和残差f(x)
- 求解线性方程 (JᵀJ)Δx = -Jᵀf(x)
- 更新参数 x = x + αΔx
直到收敛

高斯-牛顿法收敛速度快，但需要良好的初始值，且当JᵀJ接近奇异时可能不稳定。

3.2 莱文贝格-马夸特(LM)法

LM算法是高斯-牛顿法的改进版，通过引入阻尼因子来调节：

(JᵀJ + μI)Δx = -Jᵀf(x)

阻尼因子μ的调节策略：

当误差减小时，减小μ，使算法更接近高斯-牛顿法
当误差增大时，增大μ，使算法更接近梯度下降

这种自适应机制使LM算法兼具快速收敛和良好稳定性，成为非线性最小二乘的事实标准算法。

在Ceres Solver中的使用示例：

ceres::Problem problem; ceres::CostFunction* cost_function = new ceres::AutoDiffCostFunction<CostFunctor, 1, 2>(new CostFunctor); problem.AddResidualBlock(cost_function, nullptr, x); ceres::Solver::Options options; options.minimizer_progress_to_stdout = true; ceres::Solver::Summary summary; ceres::Solve(options, &problem, &summary);

4. 工程实践中的关键问题

4.1 鲁棒核函数的使用

在实际工程中，数据常包含异常值(outliers)，会严重影响最小二乘的结果。鲁棒核函数通过对大残差施加惩罚来减小异常值的影响。

常用核函数：

Huber核：对大残差施加线性惩罚
Cauchy核：对异常值更加鲁棒

Ceres中设置核函数的示例：

ceres::LossFunction* loss_function = new ceres::HuberLoss(1.0); problem.AddResidualBlock(cost_function, loss_function, parameters);

4.2 稀疏性问题处理

许多工程问题（如SLAM）会产生稀疏的雅可比矩阵。利用稀疏性可以大幅提升计算效率。

稀疏矩阵的存储方式：

压缩列存储(CSC)
压缩行存储(CSR)

Eigen中稀疏矩阵的使用：

SparseMatrix<double> J(m,n); J.insert(i,j) = value; // 填充非零元素

4.3 参数尺度归一化

当不同参数的尺度差异很大时，会导致数值问题。参数归一化可以改善条件数：

x_normalized = (x - mean)/std

在优化完成后，需要将结果转换回原始尺度。

5. 典型应用场景分析

5.1 曲线拟合问题

给定一组数据点{(xᵢ,yᵢ)}，寻找最佳拟合曲线y=f(x)。例如在传感器标定中，常用多项式拟合：

f(x) = a₀ + a₁x + a₂x² + ... + aₙxⁿ

实现要点：

多项式阶数选择：通过交叉验证确定
正则化：防止高阶多项式过拟合

5.2 SLAM中的位姿图优化

在SLAM系统中，最小二乘用于优化位姿图的约束：

min Σ||eᵢᵀΩᵢeᵢ||

其中eᵢ是约束误差，Ωᵢ是信息矩阵。这个问题通常具有稀疏性，可以使用g2o或GTSAM等专用库求解。

5.3 三维重建中的Bundle Adjustment

Bundle Adjustment通过最小化重投影误差来优化相机参数和三维点：

min Σ||π(Pᵢ,Xⱼ) - xⱼᵢ||²

其中π是投影函数，Pᵢ是相机参数，Xⱼ是三维点，xⱼᵢ是观测到的二维点。

6. 性能优化技巧

6.1 雅可比矩阵的计算优化

雅可比矩阵的计算是迭代算法的关键步骤。常用方法：

自动微分：Ceres等库提供支持
解析导数：手动推导更高效
数值差分：简单但精度低

6.2 线性求解器的选择

根据问题规模选择合适求解器：

稠密小矩阵：QR或SVD
稀疏中矩阵：Sparse QR
超大规模：共轭梯度(CG)

6.3 并行化计算

利用现代CPU/GPU的并行能力：

残差计算并行化
雅可比矩阵计算并行化
使用多线程线性代数库

7. 常见问题与调试技巧

7.1 算法不收敛的可能原因

初始值太差：尝试更好的初始化
步长过大：调整线搜索参数
数据有异常值：使用鲁棒核函数
参数尺度不一致：进行归一化

7.2 数值不稳定问题处理

添加正则化项
使用更稳定的分解方法(SVD)
增加浮点精度(double->long double)

7.3 结果验证方法

残差分析：检查分布是否随机
参数扰动测试：微小扰动应导致残差平滑变化
交叉验证：检查泛化性能

资讯详情