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

最新推荐文章于 2024-07-24 11:32:27 发布
最新推荐文章于 2024-07-24 11:32:27 发布
阅读量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 官方教程学习笔记（四）——曲线拟合Curve Fitting
小政哥的私房博客
03-26 5057
本节主要依据Ceres官方指南翻译而成。 最小二乘法（least square）和非线性最小二乘分析的本来目的就是对一组数据进行曲线拟合。本节将介绍曲线拟合的问题。本节所用的采样点根据y=e0.3x+0.1y=e0.3x+0.1y =e^{0.3x + 0.1}生成，并且加入标准差为σ=0.2σ=0.2\sigma = 0.2高斯噪声。这2n2n2n个数据，存入data[]当中。下面我们...
ceres-python:ceres-solver的Python绑定
05-17
**正文** ...通过理解和熟练使用ceres-python，用户可以解决复杂的数据拟合问题，提升算法性能，同时保持代码的简洁性和易维护性。无论是在学术研究还是工业应用中，ceres-python都是一个宝贵的资源。
Ceres学习笔记应用篇003_使用Ceres进行曲线拟合
wangchenxue的博客
07-19 746
Ceres学习笔记003_使用Ceres进行曲线拟合
Ceres学习笔记004--使用Ceres进行2D圆拟合
不学会C++不改名
11-02 753
本文将按照官方历程介绍，如何使用Ceres进行2D圆拟合
Ceres学习笔记001--初识Ceres
不学会C++不改名
11-02 1156
​ Ceres Solver 是由 Google 开发的开源 C++ 库，用于解决具有边界约束和一般无约束的非线性最小二乘鲁棒优化问题。
ceres学习笔记（三）：曲线拟合
weixin_43835047的博客
03-08 169
ceres学习笔记
Ceres学习笔记应用篇004_2D圆拟合
wangchenxue的博客
07-20 586
Ceres学习笔记004_2D圆拟合
【TOOL】ceres学习笔记(一) —— 教程练习
gene's blog
06-22 624
是由Google开发的开源C++库，用于解决具有边界约束的非线性最小二乘优化和一般无约束优化问题，成熟、功能丰富、高性能。与一般优化问题不同的是，非线性最小二乘优化问题的目标函数具有明确的物理意义——残差。
ceres学习笔记3——曲线拟合powell.cc
qq_42731705的博客
10-09 221
1.官方程序 对于曲线拟合问题，对于每个数据都会有一个残差，定义的残差对象中还需要包含观测到的数据，所以相比前面的有一些不同。由于最后是要计算残差的平方和最小，所以说对于每一组观测数据（对应一个残差），都需要添加到Problem的残差块中。 这里带估计的变量有两个，分别是m和c，官方程序还是把他们分开来了，每个优化变量都是一维的，实际同上一节的笔记记录，这里也是可以写成一个数组，把多个一维变量变成一个多维变量。 struct ExponentialResidual { ExponentialResid.
Ceres学习笔记应用篇001_Ceres Solver介绍、基本使用步骤及一个简单例子
wangchenxue的博客
07-19 1695
Ceres学习笔记001_Ceres Solver介绍、基本使用步骤及一个简单例子
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版的，加载后即可编译。
C语言 之 理解指针（1）
fallzzzzz的博客
07-20 1184
就是先获取该数组的首元素的地址，从前面我们知道，指针的类型的作用是指针+1 所跳过的字节数，int为整形，大小为4个字节，所以p+1就跳过4个字节，p+2就跳过8个字节，计算机的CPU（中央处理器）在处理数据的时候，需要的数据是在内存中读取的，处理后数据也会放回内存中，如果想对这些空间进行高效的管理，就可以。并且我们知道，数组的元素存放是随着下标增长，内存地址变大的。上面的代码创建了一个整型变量n，因为一个整形的大小是4个字节，所以内存中申请了4个字节，⽤于存放整数10，其中每个字节都有地址。
【C++修炼之路 第五章】string 类（补充）：string 扩展空间规律 + reserve 的缩容 + resize 的使用
唯有忍受别人不能忍受的，才能享受别人不能享受的
07-23 1127
【详细剖析】string 扩展空间规律 + reserve 的缩容性质与只改变capacity的性质 + resize 的使用 + shink_to_fit 的缩容函数
VS2019编译和使用gtest测试（C++）
最新发布
m0_65635427的博客
07-24 267
下载下来解压到文件夹，再在文件夹里面新建一个build文件夹，如下：再安装cmake，可以先检查一下是否安装了cmake，打开命令窗口，输入cmake -version：然后双击就对了，我这里已经安装好了，无法提供截图，注意一点就是在让选择是否建立系统环境变量（path）的时候，记得勾上第二个建立系统环境变量！下载好之后，记得检验一下是否安装成功，如果在普通用户下显示没有，建议用管理员身份打开命令窗口，再次尝试。之后命令窗口进入创建的build文件夹路径里面执行cmake ../
C++ | Leetcode C++题解之第275题H指数II
Ddddddd_158的博客
07-23 367
C++ | Leetcode C++题解之第275题H指数II
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、付费专栏及课程。

余额充值