Eigen 库实现最小二乘算法(Least Squares)
一、最小二乘法基本原理
给定一个超定方程组 A x = b Ax = b Ax=b,当 A ∈ R m × n , m > n A \in \mathbb{R}^{m \times n}, m > n A∈Rm×n,m>n 时,一般无法精确解出 x x x。因此我们寻找一个使残差 ∥ A x − b ∥ 2 2 \|Ax - b\|_2^2 ∥Ax−b∥22 最小的解。
其解析解为:
x = ( A T A ) − 1 A T b x = (A^T A)^{-1} A^T b x=(ATA)−1ATb
或者使用更稳定的方式:
- QR分解
- SVD分解
二、基于 Eigen 的最小二乘解法
方法1:正规方程(Normal Equation)
#include <Eigen/Dense>
#include <iostream>
using namespace Eigen;int main() {MatrixXd A(4, 2); // 超定方程,4个点拟合一个一次函数 y = ax + bVectorXd b(4);A << 1, 1,2, 1,3, 1,4, 1;b << 6, 5, 7, 10;// 求解 x = (A^T A)^{-1} A^T bVector2d x = (A.transpose() * A).inverse() * A.transpose() * b;std::cout << "Least Squares Solution (a, b): " << x.transpose() << std::endl;return 0;
}
说明:适用于小规模数据,数值不稳定时不推荐。
方法2:QR 分解法(推荐)
Vector2d x = A.colPivHouseholderQr().solve(b);
更稳定可靠,效率也好。支持大型数据处理。
方法3:SVD 分解(最稳健)
Vector2d x = A.bdcSvd(ComputeThinU | ComputeThinV).solve(b);
适合矩阵奇异或接近奇异的情况(如点分布退化等)。
三、完整最小二乘拟合一元线性函数示例(y = ax + b)
#include <Eigen/Dense>
#include <iostream>
#include <vector>int main() {using namespace Eigen;std::vector<double> x_data = {1, 2, 3, 4, 5};std::vector<double> y_data = {2.2, 2.8, 3.6, 4.5, 5.1};const int N = x_data.size();MatrixXd A(N, 2);VectorXd b(N);for (int i = 0; i < N; ++i) {A(i, 0) = x_data[i];A(i, 1) = 1.0; // 常数项b(i) = y_data[i];}// 求解最小二乘(推荐 QR)Vector2d result = A.colPivHouseholderQr().solve(b);std::cout << "拟合直线: y = " << result(0) << " * x + " << result(1) << std::endl;return 0;
}
四、拓展:非线性最小二乘可通过迭代(如高斯-牛顿)
如拟合非线性模型 y = exp ( a x + b ) y = \exp(ax + b) y=exp(ax+b),需使用迭代优化,比如:
- 高斯-牛顿法
- Levenberg-Marquardt 法(Ceres/GTSAM 使用)
Eigen 本身不含自动求导工具,但可以配合手动 Jacobian 实现。