Ceres学习笔记003--使用Ceres进行曲线拟合

最新推荐文章于 2024-04-25 16:59:43 发布
最新推荐文章于 2024-04-25 16:59:43 发布
阅读量1.1k 收藏 7
点赞数 3
分类专栏： Ceres学习笔记 文章标签： ceres c++ 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_45006390/article/details/127647998
版权

使用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) {}

	template<typename T>
	bool 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<ExponentialResidual, 1, 1, 1>(new ExponentialResidual(data[2 * i], data[2 * i + 1])),
			NULL,		//	损失函数、单位函数
			&m,		//	第1个参数块
			&c);		//	第2个参数块
	}

	//	配置求解器参数
	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 << " Initial m: " << 0.0 << " c: " << 0.0 << "\n";
	std::cout << " Final m: " << m << " c: " << c << "\n";

	std::system("pause");
	return 0;
}

输出结果如下：

 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

根据输出结果可以看出，参数迭代初始值为 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\{\begin{matrix} \frac{\partial f}{\partial m} = - x * e^{mx + c} \\ \frac{\partial f}{\partial c} = - e^{mx+c} \end{matrix}\right. \tag{4} {∂m∂f​=−x∗emx+c∂c∂f​=−emx+c​(4)
完整代码如下：

#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;
}

结果如下：

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

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);    // 第二个参数块

加入损失函数后的曲线拟合效果如下图：
增加损失函数实现鲁棒拟合

看到我请叫我学C++
关注
  • 3
    点赞
  • 踩
  • 7
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
Ceres solver 学习-求解器类型
zfjBIT的专栏
04-06 1137
Ceres solver 里面定义了8种线性求解器（ceres/include/ceres/types.h），分别是: **DENSE_QR:**用于小规模最小二乘问题的求解，基于Eigen； List item DENSE_NORMAL_CHOLESKY&SPARSE_NORMAL_CHOLESKY： cholesky分解，用于具有稀疏性的大规模非线性最小二乘问题求解；前者基于Eigen，后者基于SuiteSparse or CXSparse； DENSE_SCHUR&SPARSE_SCH
ceres之共轭梯度求解器
xiechaoyi123的博客
06-15 989
参考源码：https://github.com/ceres-solver/ceres-solver/blob/master/internal/ceres/conjugate_gradients_solver.cc 求解目标：Ax=B,已知系数矩阵A和残差项B，求解使得Ax=B的x 求解输入：（1）A,B （2）控制参数：r_tolerance(残差项下限，默认为0) q_tolerance(二次项下限.默认为0)： ...
Ceres Solver从零开始手把手教学使用
主要学习关于图像处理和SLAM领域
11-03 1804
笔者已经半年没有更新新的内容了，最近学习视觉SLAM的过程中发现自己之前学习的库基础不够扎实。Ceres作为一个优化库在视觉SLAM中使用非常广泛，笔者会跟随着Ceres Solver的官方教程一直学习下来，并且会在专栏中持续更新，如需转载请声明出处。
fatal error: ceres/ceres.h: 没有那个文件或目录
最新发布
weixin_58848887的博客
04-25 418
编译过程中遇到：“fatal error: ceres/ceres.h: 没有那个文件或目录”问题出现在你的Ceres安装出现问题。
Ceres 曲线拟合
泠山的博客
09-16 1143
Ceres 是一个广泛使用的最小二乘问题求解库。在 Ceres 中，我们作为用户，只需按照一定步骤定义待解的优化问题，然后交给求解器计算。min⁡x12∑iρi(∥fi(xi1,…xin)∥2)s.t.lj≤xj≤ujminx​s.t.​21​∑i​ρi​(∥fi​(xi1​​,…xin​​)∥2)lj​≤xj​≤uj​​在这个问题中，x。
Ceres入门(二)曲线拟合与BA求解
dushib的博客
01-26 1286
在一篇博客中，通过分析helloword的自动求导和节写求导简单例子，了解了Ceres的基本流程。本片博客在上一片基础之上，以高博十四讲内容为基础，分析Ceres两个使用案例 一、曲线拟合 1、问题描述 其中a,b,c为待估计的参数，w为噪声。在程序里利用模型生成x，y的数据，在给数据添加服从高斯分布的噪声。之后用ceres优化求解参数a,b,c。 2、求解代码 代码部分仍然与上一篇博客类似，分为三个部分 (1)、第一部分：构建cost fuction，即代价函数 struct CURV
ceres的常见用法总结
Enochzhu的博客
08-08 1811
ceres学习笔记用法
ceres-python:ceres-solver的Python绑定
05-17
**正文** ...通过理解和熟练使用ceres-python，用户可以解决复杂的数据拟合问题，提升算法性能，同时保持代码的简洁性和易维护性。无论是在学术研究还是工业应用中，ceres-python都是一个宝贵的资源。
ceres-solver-2.1.0.7z
03-11
在`ceres-solver-2.1.0`目录中，你通常会找到以下文件和子目录： - `include/`：包含Ceres库的头文件，用户需要引用这些头文件来使用Ceres的功能。 - `src/`：源代码目录，实现了所有的优化算法和功能。 - `...
ceres-solver-1.13.0.tar.gz
05-21
这个"ceres-solver-1.13.0.tar.gz"文件是Ceres Solver的1.13.0版本，包含了该库的所有源代码和其他相关文件，适用于Linux、Mac OS X和Windows等多种操作系统。以下将详细介绍Ceres Solver的关键知识点。 1. **非...
第三章作业代码及Ceres学习笔记-kaho1
08-03
**曲线拟合**\n在“第三章作业”中提到了曲线拟合，这是Ceres的一个典型应用场景。通过构建误差项并最小化这些误差，Ceres能够找到最佳的参数设置，使得模型尽可能地贴近实际数据。对于解析求导，Ceres默认使用加减...
ceres-solver-1.14.0编译库(windows)
06-29
win10下，通过cmake3.22.5和vs2015编译ceres1.14.0版本库，内含eigen3.4.0、gflags-2.2.2、glog-0.6.0、suitesparse-1.4.0和libcxsparse.lib，并附有vs属性表，包括debug和release版的，加载后即可编译。
ceres使用
qq_41485193的博客
08-09 1176
ceres
Ceres简介及示例（4）Curve Fitting(曲线拟合)
热爱生活，忠于自己
04-02 851
首先定义一个模板对象来计算残差。每一个观察值(采样点)都有一个残差，
Ceres入门——Ceres的基本使用方法
热门推荐
Mr.Andy的博客
05-02 1万+
Ceres入门——Ceres的基本使用方法1.使用流程2.实例分析2.1 构建代价函数（cost function） Ceres solver 是谷歌开发的一款用于非线性优化的库，在谷歌的开源激光雷达slam项目cartographer中被大量使用。Ceres官网上的文档非常详细地介绍了其具体使用方法，相比于另外一个在slam中被广泛使用的图优化库G2O，ceres的文档可谓相当丰富详细。本篇博客...
ceres优化库的使用
小博的博客
09-25 2106
谷歌ceres优化库在slam领域的使用
Ceres solver使用教程
zhaoliang38的博客
03-30 581
Ceres solver使用教程
Ceres 中的线性求解器类型（linear_solver_type）
一抹烟霞的博客
03-29 688
Ceres Solver中的线性求解器类型(linear_solver_type)有多个选项，包括： DENSE_QR：使用稠密QR分解方法求解线性方程组。适用于内存足够的小规模问题，求解速度较快。 DENSE_SCHUR：使用Schur补方法进行求解。适用于大规模问题，但内存消耗较大。 SPARSE_NORMAL_CHOLESKY：使用稀疏Cholesky分解方法求解正规方程。适用于具有稀...
Ceres 库：基础使用，以手写高斯-牛顿法为例
JasonLi的博客
02-17 2414
记录了使用Ceres库的基础方式，包括安装、编译、CMake工程建立，基础API的使用等内容。并以手写高斯-牛顿法为例，进行演示学习。
ceres-solver-1.14.0.tar.gz
10-03
ceres-solver-1.14.0.tar.gz 是一个开源的库文件，用于解决非线性最小二乘问题。它提供了一个高效、灵活和稳定的框架，可用于求解数学模型中的最优化问题。 这个版本的ceres-solver库是经过一系列改进和优化的结果...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值