-
Ceres学习笔记003--使用Ceres进行曲线拟合
使用Ceres进行曲线拟合
1 生成数据
离散数据点通过下述方法生成：
y = e 0.3 x + 0.1 + N ( 0 , σ 2 ) (1) y = e^{0.3x+0.1} + N(0,\sigma ^2) \tag{1} y=e0.3x+0.1+N(0,σ2)(1)
式中， N ( 0 , σ 2 ) N(0,\sigma^2) N(0,σ2) 表示标准差为 σ \sigma σ 的零均值高斯分布噪声。注：Ceres 官方文档中说的是高斯噪声，但样例代码中给出的数据实际上是随机噪声，因此，本文给出的数据与官方样例代码中的数据会不太一样，不过这都是细节，几乎没有影响。
2 构建最小二乘问题
对于 Ceres 使用基本步骤，本文不再赘述，可以参考之前的文章。
未知参数的曲线函数表达式如下：
y = e m x + c (2) y = e^{mx + c}\tag{2} y=emx+c(2)
用通过式（1）生成的数据点来拟合上述曲线，因此构建如下非线性最小二乘问题：
m i n ∑ i 1 2 ∣ ∣ y i − e m x i + c ∣ ∣ 2 (3) min\sum_i \frac{1}{2}||y_i - e^{mx_i + c}||^2 \tag{3} mini∑21∣∣yi−emxi+c∣∣2(3)3 代码实践
以下分别根据自动微分法和解析解法来实现
3.1 自动微分法
完整代码如下：
#include "ceres/ceres.h" #include "glog/logging.h" #include "opencv2/core.hpp" // 用户自定义残差计算模型 struct ExponentialResidual { // 样本数据观测值通过构造函数传入 ExponentialResidual(double x, double y) : x_(x), y_(y) {} templatebool operator()(const T* const m, const T* const c, T* residual) const { // 输出残差维度为1 // 输入2个参数块，每个参数块的维度为1 residual[0] = T(y_) - exp(m[0] * T(x_) + c[0]); return true; } private: const double x_; const double y_; }; int main(int argc, char** argv) { google::InitGoogleLogging(argv[0]); // 生成数据 const int kNumObservations = 67; // 生成67个点 double sigma = 0.2; // 标准差 cv::RNG rng; // Opencv随机数产生器 double data[2 * kNumObservations]; // 数据容器[x1,y1,x2,y2,...] for(int i = 0; i < kNumObservations; ++i) { double x = 0.075 * i; double y = exp(0.3 * x + 0.1); double noise = rng.gaussian(sigma); data[2 * i] = x; data[2 * i + 1] = y + noise; } // 设置参数初始值 // 输入2个参数块，每个参数块的维度为1 double m = 0.3; double c = 0.1; // 构建非线性最小二乘问题 ceres::Problem problem; for(int i = 0; i < kNumObservations; ++i) { // 添加残差块，需要一次指定代价函数、损失函数、参数块 // 本例中损失函数为单位函数 problem.AddResidualBlock( // 输出残差维度为1，输出参数块有2个，每个参数块维度都为1 new ceres::AutoDiffCostFunction - 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
输出结果如下：
Initial m: 0 c: 0 iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time 0 1.204244e+02 0.00e+00 3.60e+02 0.00e+00 0.00e+00 1.00e+04 0 2.31e-03 2.90e-03 1 2.425961e+03 -2.31e+03 3.60e+02 8.02e-01 -1.95e+01 5.00e+03 1 1.26e-03 4.42e-03 2 2.422258e+03 -2.30e+03 3.60e+02 8.02e-01 -1.95e+01 1.25e+03 1 4.97e-04 5.03e-03 3 2.400245e+03 -2.28e+03 3.60e+02 7.98e-01 -1.93e+01 1.56e+02 1 4.80e-04 5.60e-03 4 2.210383e+03 -2.09e+03 3.60e+02 7.66e-01 -1.77e+01 9.77e+00 1 4.73e-04 6.15e-03 5 8.483095e+02 -7.28e+02 3.60e+02 5.71e-01 -6.32e+00 3.05e-01 1 4.75e-04 6.70e-03 6 3.404435e+01 8.64e+01 4.10e+02 3.12e-01 1.37e+00 9.16e-01 1 2.15e-03 8.93e-03 7 7.242644e+00 2.68e+01 1.84e+02 1.27e-01 1.11e+00 2.75e+00 1 2.17e-03 1.12e-02 8 3.933925e+00 3.31e+00 5.81e+01 3.45e-02 1.03e+00 8.24e+00 1 2.15e-03 1.34e-02 9 2.333679e+00 1.60e+00 2.52e+01 9.77e-02 9.90e-01 2.47e+01 1 3.09e-03 1.77e-02 10 1.419436e+00 9.14e-01 8.73e+00 1.16e-01 9.83e-01 7.42e+01 1 2.33e-03 2.02e-02 11 1.245322e+00 1.74e-01 1.43e+00 6.69e-02 9.89e-01 2.22e+02 1 2.20e-03 2.25e-02 12 1.237976e+00 7.35e-03 9.21e-02 1.61e-02 9.91e-01 6.67e+02 1 2.74e-03 2.54e-02 13 1.237935e+00 4.11e-05 9.96e-04 1.28e-03 9.90e-01 2.00e+03 1 2.48e-03 2.81e-02 Solver Summary (v 2.0.0-eigen-(3.3.8)-no_lapack-eigensparse-no_openmp) Original Reduced Parameter blocks 2 2 Parameters 2 2 Residual blocks 67 67 Residuals 67 67 Minimizer TRUST_REGION Dense linear algebra library EIGEN Trust region strategy LEVENBERG_MARQUARDT Given Used Linear solver DENSE_QR DENSE_QR Threads 1 1 Linear solver ordering AUTOMATIC 2 Cost: Initial 1.204244e+02 Final 1.237935e+00 Change 1.191865e+02 Minimizer iterations 14 Successful steps 9 Unsuccessful steps 5 Time (in seconds): Preprocessor 0.000587 Residual only evaluation 0.002542 (14) Jacobian & residual evaluation 0.016143 (9) Linear solver 0.004534 (14) Minimizer 0.028694 Postprocessor 0.000090 Total 0.029372 Termination: CONVERGENCE (Function tolerance reached. |cost_change|/cost: 1.219979e-08 <= 1.000000e-06) Final m: 0.301681 c: 0.0907709- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
根据输出结果可以看出，参数迭代初始值为 m = 0, c = 0 ，初始目标函数残差为 120.4244，最终参数收敛于 m = 0.301681, c = 0.0907709，残差为 1.237935，最终估计出的参数并不完全为 m = 0.3, c = 0.1，这样的偏差是符合预期的，因为生成的数据中包含了噪声，事实上，如果使用 m = 0.3, c = 0.1 计算残差的话，残差值为 1.241615，比 m = 0.301681, c = 0.0907709 时的残差值 1.237935 更大。
3.2 解析解法
对于 f = y − e m x + c f = y - e^{mx + c} f=y−emx+c，偏导函数如下：
{ ∂ f ∂ m = − x ∗ e m x + c ∂ f ∂ c = − e m x + c (4) \left\{\right. \tag{4} {∂m∂f=−x∗emx+c∂c∂f=−emx+c(4)∂ f ∂ m = − x ∗ e m x + c ∂ f ∂ c = − e m x + c
完整代码如下：#include "ceres/ceres.h" #include "glog/logging.h" #include "opencv2/core.hpp" class ExponentialResidual : public ceres::SizedCostFunction<1, /*输出（Residual）维度大小*/\ 1, /*第1个输入参数块维度大小*/\ 1 /*第2个输入参数块维度大小*/> { public: ExponentialResidual(double x, double y) : x_(x), y_(y) {} virtual ~ExponentialResidual() {} // 用户自定义残差计算方法 virtual bool Evaluate(double const* const* x /*输入参数块*/, double* residuals /*输出残差*/, double** jacobians /*输出雅可比矩阵*/) const { // 本例中有两个输入参数块，每个参数块中有1个参数 double m = x[0][0]; double c = x[1][0]; // 本例中输出残差维度为1 double y0 = exp(m * x_ + c); residuals[0] = y_ - y0; if(jacobians == NULL) { return true; } // 残差对第1个参数块中的参数依次求偏导，即对m求偏导 if(jacobians[0] != NULL) { jacobians[0][0] = -y0 * x_; } // 残差对第2个参数块中的参数依次求偏导，即对c求偏导 if(jacobians[1] != NULL) { jacobians[1][0] = -y0; } return true; } private: const double x_; const double y_; }; int main(int argc, char** argv) { google::InitGoogleLogging(argv[0]); // 生成数据 const int kNumObservations = 67; // 67个点 double sigma = 0.2; // 标准差 cv::RNG rng; // Opencv随机数产生器 double data[2 * kNumObservations]; // 数据容器[x1,y1,x2,y2,......] for(int i = 0; i < kNumObservations; ++i) { double x = 0.075 * i; double y = exp(0.3 * x + 0.1); double noise = rng.gaussian(sigma); data[2 * i] = x; data[2 * i + 1] = y + noise; } // 设置参数初始值 // 输入2个参数块，每个参数块的维度为1 double m = 0.0; double c = 0.0; // 构建非线性最小二乘问题 ceres::Problem problem; for (int i = 0; i < kNumObservations; ++i) { // 添加残差块，需要依次指定代价函数、损失函数、参数块 // 本例中损失函数为单位函数 problem.AddResidualBlock(new ExponentialResidual(data[2 * i], data[2 * i + 1]), NULL, &m, &c); } // 配置求解器参数 ceres::Solver::Options options; // 设置最大迭代次数 options.max_num_iterations = 25; // 指定线性求解器来求解问题 options.linear_solver_type = ceres::DENSE_QR; // 输出每次迭代的信息 options.minimizer_progress_to_stdout = true; // 输出日志内容 ceres::Solver::Summary summary; // 开始优化求解 ceres::Solve(options, &problem, &summary); // 输出优化过程及结果 std::cout << summary.FullReport() << "\n"; std::cout << "Final m: " << m << " c: " << c << "\n"; std::system("pause"); return 0; }- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
结果如下：
iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time 0 1.204244e+02 0.00e+00 3.60e+02 0.00e+00 0.00e+00 1.00e+04 0 7.01e-04 1.22e-03 1 2.425961e+03 -2.31e+03 3.60e+02 8.02e-01 -1.95e+01 5.00e+03 1 1.37e-03 3.17e-03 2 2.422258e+03 -2.30e+03 3.60e+02 8.02e-01 -1.95e+01 1.25e+03 1 4.39e-04 3.82e-03 3 2.400245e+03 -2.28e+03 3.60e+02 7.98e-01 -1.93e+01 1.56e+02 1 3.70e-04 4.27e-03 4 2.210383e+03 -2.09e+03 3.60e+02 7.66e-01 -1.77e+01 9.77e+00 1 3.63e-04 4.69e-03 5 8.483095e+02 -7.28e+02 3.60e+02 5.71e-01 -6.32e+00 3.05e-01 1 3.62e-04 5.11e-03 6 3.404435e+01 8.64e+01 4.10e+02 3.12e-01 1.37e+00 9.16e-01 1 7.48e-04 5.92e-03 7 7.242644e+00 2.68e+01 1.84e+02 1.27e-01 1.11e+00 2.75e+00 1 7.51e-04 6.73e-03 8 3.933925e+00 3.31e+00 5.81e+01 3.45e-02 1.03e+00 8.24e+00 1 7.44e-04 7.54e-03 9 2.333679e+00 1.60e+00 2.52e+01 9.77e-02 9.90e-01 2.47e+01 1 7.65e-04 8.37e-03 10 1.419436e+00 9.14e-01 8.73e+00 1.16e-01 9.83e-01 7.42e+01 1 7.47e-04 9.18e-03 11 1.245322e+00 1.74e-01 1.43e+00 6.69e-02 9.89e-01 2.22e+02 1 7.44e-04 9.99e-03 12 1.237976e+00 7.35e-03 9.21e-02 1.61e-02 9.91e-01 6.67e+02 1 7.44e-04 1.08e-02 13 1.237935e+00 4.11e-05 9.96e-04 1.28e-03 9.90e-01 2.00e+03 1 7.43e-04 1.16e-02 Solver Summary (v 2.0.0-eigen-(3.3.8)-no_lapack-eigensparse-no_openmp) Original Reduced Parameter blocks 2 2 Parameters 2 2 Residual blocks 67 67 Residuals 67 67 Minimizer TRUST_REGION Dense linear algebra library EIGEN Trust region strategy LEVENBERG_MARQUARDT Given Used Linear solver DENSE_QR DENSE_QR Threads 1 1 Linear solver ordering AUTOMATIC 2 Cost: Initial 1.204244e+02 Final 1.237935e+00 Change 1.191865e+02 Minimizer iterations 14 Successful steps 9 Unsuccessful steps 5 Time (in seconds): Preprocessor 0.000520 Residual only evaluation 0.001796 (14) Jacobian & residual evaluation 0.003180 (9) Linear solver 0.003432 (14) Minimizer 0.011550 Postprocessor 0.000050 Total 0.012121 Termination: CONVERGENCE (Function tolerance reached. |cost_change|/cost: 1.219979e-08 <= 1.000000e-06) Final m: 0.301681 c: 0.0907709- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
4 鲁棒最小二乘拟合曲线
假设数据存在一些外点（完全不符合曲线模型的噪点），如果还是使用上述代码进行最小二乘拟合，拟合得到的曲线将如图所示：
为了减少异常值或外点（outliers）对非线性最小二乘优化问题的影响，常用的技巧是使用损失函数，损失函数能够抑制很大的残差，而通常只有外点才会有很大的残差值，相当于减少外点的权重，也就是鲁棒拟合。Ceres自带的损失函数有 TrivialLoss, HuberLoss, SoftLOneLoss, CauchyLoss, ArctanLoss, TolerantLoss, TukeyLoss, ComposedLoss, ScaledLoss。
以 CauchyLoss 为例，使用方法如下：
problem.AddResidualBlock( new ExponentialResidual(data[2 * i], data[2 * i + 1]), //Y用户自定义成本函数 new ceres::CauchyLoss(0.5), // 损失函数，CauchyLoss函数,0.5表示尺度 &m, // 第一个参数块 &c); // 第二个参数块- 1
- 2
- 3
- 4
- 5
加入损失函数后的曲线拟合效果如下图：
-
相关阅读:
0143__c++filt(1) command
第一类曲面积分：曲面微元dσ与其投影面积微元dxdy之间的关系推导
C++ 背包问题——01背包
告别单调，Django后台主页改造 - 使用AdminLTE组件
LeetCode199:二叉树的右视图
基于JavaSwing开发打小鸟游戏+实验报告 课程设计 大作业 毕业设计
document.load和document.ready之间的区别
【Python】OS 模块简介
【Jenkins+K8s】持续集成与交付 （二十）：K8s集群通过Deployment方式部署安装Jenkins
如何关闭Google Chrome的跨域限制（cross-origin）
- 原文地址：https://blog.csdn.net/qq_45006390/article/details/127647998
