1. Ceres 简介

Ceres 是一个广泛使用的最小二乘问题求解库。在 Ceres 中，我们作为用户，只需按照一定步骤定义待解的优化问题，然后交给求解器计算。Ceres 求解的最小二乘问题最一般的形式如下（带边界的核函数最小二乘）：
min⁡x12∑iρi(∥fi(xi1,…xin)∥2)s.t. lj≤xj≤uj\begin{array}{ll} \min _x & \frac{1}{2} \sum_i \rho_i\left(\left\|f_i\left(x_{i_1}, \ldots x_{i_n}\right)\right\|^2\right) \\ \text { s.t. } & l_j \leq x_j \leq u_j \end{array} minx s.t. 21∑iρi(∥fi(xi1,…xin)∥2)lj≤xj≤uj

在这个问题中，x1,...xnx_1,...x_nx1,...xn 为优化变量，又称参数块(Parameter blocks)，fif_ifi 称为代价函数(Cost function)，也称为残差块(Residual blocks)，在 SLAM 中也可理解为误差项。ljl_jlj 和 uju_juj 为第 jjj 个优化变量的上限和下限。在最简单的情况下，取 lj=−∞l_j=-\inftylj=−∞，uj=∞u_j=\inftyuj=∞ (不限制优化变量的边界)。此时，目标函数由许多平方项经过一个核函数 ρ(⋅)\rho(\cdot)ρ(⋅) 之后求和组成。同样，可以取 ρ\rhoρ 为恒等函数，那么目标函数即为许多项的平方和，我们就得到了无约束的最小二乘问题。

为了让 Ceres 帮我们求解这个问题，我们需要做以下几件事：

定义每个参数块。参数块通常为平凡的向量，但是在 SLAM 里也可以定义成四元数、李代数这种特殊的结构。如果是向量，那么我们需要为每个参数块分配一个 double 数组来存储变量的值。
定义残差块的计算方式。残差块通常关联若干个参数块，对它们进行一些自定义的计算，然后返回残差值。Ceres 对它们求平方和之后，作为目标函数的值。
残差块往往也需要定义雅克比的计算方式。在 Ceres 中，可以使用它提供的“自动求导”功能，也可以手动指定雅克比的计算过程。如果要使用自动求导，那么残差块需要按照特定的写法书写：残差的计算过程应该是一个带模板的括号运算符。
把所有的参数块和残差块加入 Ceres 定义的 Problem 对象中，调用 Solve 函数求解即可。求解之前，我们可以传入一些配置信息，例如迭代次数、终止条件等，也可以使用默认的配置。

2. 使用 Ceres 拟合曲线

下面的代码演示了如何使用 Ceres 求解问题：

#include <iostream>
#include <opencv2/core/core.hpp>
#include <ceres/ceres.h>
#include <chrono>using namespace std;// 代价函数的计算模型
struct CURVE_FITTING_COST {CURVE_FITTING_COST(double x, double y) : _x(x), _y(y) {}// 残差的计算template<typename T>bool operator()(const T *const abc, // 模型参数，有3维T *residual) const {// 这里残差的计算公式：y-exp(ax^2+bx+c) residual[0] = T(_y) - ceres::exp(abc[0] * T(_x) * T(_x) + abc[1] * T(_x) + abc[2]); return true;}const double _x, _y;    // x,y数据
};int main(int argc, char **argv) {double ar = 1.0, br = 2.0, cr = 1.0;         // 真实参数值double ae = 2.0, be = -1.0, ce = 5.0;        // 估计参数值int N = 100;                                 // 数据点double w_sigma = 1.0;                        // 噪声Sigma值double inv_sigma = 1.0 / w_sigma;cv::RNG rng;                                 // OpenCV随机数产生器vector<double> x_data, y_data;      // 数据/** @code 生成一组带噪声的数据*/  for (int i = 0; i < N; i++) {double x = i / 100.0;x_data.push_back(x);y_data.push_back(exp(ar * x * x + br * x + cr) + rng.gaussian(w_sigma * w_sigma));}/** @endcode */double abc[3] = {ae, be, ce};// 构建最小二乘问题ceres::Problem problem;for (int i = 0; i < N; i++) {problem.AddResidualBlock(     // 向问题中添加误差项// 使用自动求导，模板参数：误差类型，输出维度，输入维度，维数要与前面struct中一致new ceres::AutoDiffCostFunction<CURVE_FITTING_COST, 1, 3>(new CURVE_FITTING_COST(x_data[i], y_data[i])),nullptr,            // 核函数，这里不使用，为空abc                 // 待估计参数);}// 配置求解器ceres::Solver::Options options;     // 这里有很多配置项可以填options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;  // 增量方程如何求解options.minimizer_progress_to_stdout = true;   // 输出到coutceres::Solver::Summary summary;                // 优化信息chrono::steady_clock::time_point t1 = chrono::steady_clock::now();ceres::Solve(options, &problem, &summary);  // 开始优化chrono::steady_clock::time_point t2 = chrono::steady_clock::now();chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);cout << "solve time cost = " << time_used.count() << " seconds. " << endl;// 输出结果cout << summary.BriefReport() << endl;cout << "estimated a,b,c = ";for (auto a:abc) cout << a << " ";cout << endl;return 0;
}

程序需要说明的地方均已加注释。可以看到，我们利用 OpenCV 的噪声生成器生成了 100100100 个带高斯噪声的数据，随后利用 Ceres 进行拟合。这里演示的 Ceres 用法有如下几项：

定义残差块的类。方法是书写一个类（或结构体），并在类中定义带模板参数的 () 运算符，这样该类就成为了一个拟函数（Functor）。这种定义方式使得 Ceres 可以像调用函数一样，对该类的某个对象（比如说 a）调用 a() 方法。事实上，Ceres 会把雅克比矩阵作为类型参数传入此函数，从而实现自动求导的功能。
程序中的 double abc[3] 即参数块，而对于残差块，我们对每一个数据构造 CURVE_FITTING_COST 对象，然后调用 AddResidualBlock 将误差项添加到目标函数中。由于优化需要梯度，我们有若干种选择：
（1）使用 Ceres 的自动求导（Auto Diff）；
（2）使用数值求导（Numeric Diff）；
（3）自行推导解析的导数形式，提供给 Ceres。
因为自动求导在编码上是最方便的，于是我们使用自动求导(1)。
自动求导需要指定误差项和优化变量的维度(即AutoDiffCostFunction需要传入输入维度和输出维度)。这里的误差项是标量，维度为 1；优化的是 a, b, c 三个量，维度为 3。于是，在自动求导类 AutoDiffCostFunction 的模板参数中设定变量维度为 1、3。
设定好问题后，调用 Solve 函数进行求解。你可以在 options 里配置（非常详细的）优化选项。例如，我们可以选择使用 Line Search 还是 Trust Region，迭代次数，步长等等。读者可以查看 Options 的定义，看看有哪些优化方法可选，当然默认的配置已经可以用在很广泛的问题上了。

2.1 数值求导与自动求导的区别

We will now consider automatic differentiation. It is a technique that can compute exact derivatives, fast, while requiring about the same effort from the user as is needed to use numerical differentiation.

3. 数值解的使用方法

// 1.这里使用的是ceres的自动求导
struct FITTING_COST{FITTING_COST(double molecule,double denominator,double delta_combine):_molecule(molecule),_denominator(denominator),_delta_distance(delta_combine){}// 残差的计算template<typename T>bool operator()(const T *const res, T *residual) const{residual[0] = (res[0] + _delta_distance) * (_denominator - res[1]) - _molecule;return true;}const double _molecule, _denominator, _delta_distance;
};// 2.这里使用的是解析解求导方法
// A CostFunction implementing analytically derivatives for the function
class FITTING_COST_EXTRA: public ceres::SizedCostFunction<1 /* number of residuals */,2 /* size of first parameter */> {public:FITTING_COST_EXTRA(double molecule,double denominator,double delta_combine):_molecule(molecule),_denominator(denominator),_delta_distance(delta_combine){}virtual ~FITTING_COST_EXTRA() {}virtual bool Evaluate(double const* const* res,double* residual,double** jacobians) const {residual[0] = (res[0][0] + _delta_distance) * (_denominator - res[0][1]) - _molecule;// Since the Evaluate function can be called with the jacobians// pointer equal to NULL, the Evaluate function must check to see// if jacobians need to be computed.//// For this simple problem it is overkill to check if jacobians[0]// is NULL, but in general when writing more complex// CostFunctions, it is possible that Ceres may only demand the// derivatives w.r.t. a subset of the parameter blocks.if (jacobians != NULL && jacobians[0] != NULL) {jacobians[0][0] = _denominator - res[0][1];jacobians[0][1] = -res[0][0] - _delta_distance;}return true;}const double _molecule, _denominator, _delta_distance;
};int main(){......................................................................double res[2] = {res_distance, res_delta_x};//build the least square problemceres::Problem problem;// 1.自动求导方法// for (size_t j = 0; j < bucket.size(); j++)// {//     problem.AddResidualBlock(   //add error terms//         new ceres::AutoDiffCostFunction<FITTING_COST,1,2>(//             new FITTING_COST(//                 bucket[j]._pre_molecule,//                 bucket[j]._pre_denominator,//                 bucket[j]._pre_delta_distance)//         ),//         nullptr,//         res//     );// }// 2.自行推导导数方式for (size_t j = 0; j < bucket.size(); j++){problem.AddResidualBlock(   //add error termsnew FITTING_COST_EXTRA(bucket[j]._pre_molecule,bucket[j]._pre_denominator,bucket[j]._pre_delta_distance),nullptr,res);}problem.SetParameterUpperBound(res,1,0.6);problem.SetParameterLowerBound(res,1,-0.5);//  config solversceres::Solver::Options options;options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;options.minimizer_progress_to_stdout = false;//options.min_line_search_step_size = 0.5;//options.max_num_iterations = 10;ceres::Solver::Summary summary; //optimize infoceres::Solve(options, &problem, &summary); //start optimization......................................................................
}

Ceres 曲线拟合相关推荐

一夜爆火的SLAM技术即将颠覆哪些领域
SLAM的英文全程是 Simultaneous Localization and Mapping,中文称作「同时定位与地图创建」.SLAM试图解决这样的问题:一个机器人在未知的环境中运动,如何通过对环 ...
深蓝视觉SLAM课程第四讲--相机模型，非线性优化（G2O）
课程Github地址:https://github.com/wrk666/VSLAM-Course/tree/master 0. 内容对应于十四讲中的第5讲和第6讲回顾十四讲P24-26 1. 针 ...
SLAM专题（6）-- 非线性优化
0. 摘要: SLAM中经常遇到非线性优化问题:多个误差项平方和组成的最小二乘问题,两个最常见的梯度下降方法-非线性优化方案:高斯牛顿法.裂纹伯格-马夸尔特方法. 两个基于c++编写的优化库:来自谷歌 ...
视觉SLAM理论与实践学习笔记
需要学习注释起来方便过后忘记留存代码示例在尾部 You Know BOOL State::XueXi(StateOutFile& file, const char* name, Obj ...
曲线拟合（高斯牛顿法，Ceres,g2o)
实践:曲线拟合问题我们先通过高斯牛顿法来求解最小二乘问题,然后介绍使用优化库(Ceres/g2o)来求解此问题. 例题: 考虑一条满足如下方程的曲线: y = exp(ax2x^2x2+ bx + ...
非线性优化库Ceres学习笔记7（鲁棒的曲线拟合）
现在假设给定的数据有一些异常值,即我们有一些不服从噪声模型的点.如果我们使用上面的代码来拟合这些数据,我们将得到如下所示的拟合. 在这个时候需要应用损失函数(Loss function)来对异常数据进 ...
Ceres入门——Ceres的基本使用方法
Ceres入门--Ceres的基本使用方法 1.使用流程 2.实例分析--HelloWorld 2.1 构建代价函数(cost function) 2.2 构建寻优问题 2.3 配置并运行求解器 2. ...
Ceres和g2o的配置和使用
上文非线性优化介绍了非线性优化的基本求解方法,并使用C++手动实现了曲线拟合实例.本文介绍ceres和g2o库的配置方法,并通过曲线拟合实例介绍其使用方法. Ceres安装 Google Ceres ...
视觉SLAM十四讲学习笔记-第六讲-非线性优化的实践-高斯牛顿法和曲线拟合
专栏系列文章如下: 视觉SLAM十四讲学习笔记-第一讲_goldqiu的博客-CSDN博客视觉SLAM十四讲学习笔记-第二讲-初识SLAM_goldqiu的博客-CSDN博客视觉SLAM十四讲学习 ...
视觉SLAM——ceres非线性最小二乘求解器
ceres求解器曲线拟合代码分解一二 2.1 2.2 2.3

Ceres 曲线拟合