⼀⽂助你Ceres⼊门——CeresSolver新⼿向全攻略
Ceres solver 是⾕歌开发的⼀款⽤于⾮线性优化的库,在⾕歌的开源激光雷达slam项⽬cartographer中被⼤量使⽤。Ceres官⽹上的⽂档⾮常详细地介绍了其具体使⽤⽅法,相⽐于另外⼀个在slam中被⼴泛使⽤的图优化库G2O,ceres的⽂档可谓相当丰富详细(没有对⽐就没有伤害,主要是G2O资料太少了,对⽐起来就显得ceres的很多),下⾯我就介绍下如何使⽤ceres库进⾏简单的⾮线性优化,给各位ceres的初学者⼀个低门槛的⼊门教程,能⽐较直观地对使⽤ceres库求解优化问题的过程有⼀个清晰的了解。话不多说,开整。
新⼿村-Ceres简易例程
使⽤Ceres求解⾮线性优化问题,⼀共分为三个部分:
1、 第⼀部分:构建cost fuction,即代价函数,也就是寻优的⽬标式。这个部分需要使⽤仿函数(functor)这⼀技巧来实现,做法是定义⼀个cost function的结构体,在结构体内重载()运算符,具体实现⽅法后续介绍。
2、 第⼆部分:通过代价函数构建待求解的优化问题。
3、 第三部分:配置求解器参数并求解问题,这个步骤就是设置⽅程怎么求解、求解过程是否输出等,
然后调⽤⼀下Solve⽅法。
好了,此时你应该对ceres的⼤概使⽤流程有了⼀个基本的认识。下⾯我就基于ceres官⽹上的教程中的⼀个例程来详细介绍⼀下ceres的⽤法。
Ceres官⽹教程给出的例程中,求解的问题是求x使得取到最⼩值。(很容易⼼算出x的解应该是10)
好,来看代码:
#include<iostream>
#include<ceres/ceres.h>
using namespace std;
using namespace ceres;
//第⼀部分:构建代价函数,重载()符号,仿函数的⼩技巧
struct CostFunctor {
template <typename T>
bool operator()(const T* const x, T* residual) const {
residual[0] = T(10.0) - x[0];
return true;
}
};
//主函数
int main(int argc, char** argv) {
google::InitGoogleLogging(argv[0]);
// 寻优参数x的初始值,为5
double initial_x = 5.0;
double x = initial_x;
/
/ 第⼆部分:构建寻优问题
Problem problem;
CostFunction* cost_function =
new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor); //使⽤⾃动求导,将之前的代价函数结构体传⼊,第⼀个1是输出维度,即残差的维度,第 problem.AddResidualBlock(cost_function, NULL, &x); //向问题中添加误差项,本问题⽐较简单,添加⼀个就⾏。
//第三部分:配置并运⾏求解器
Solver::Options options;
options.linear_solver_type = ceres::DENSE_QR; //配置增量⽅程的解法
options.minimizer_progress_to_stdout = true;//输出到cout
Solver::Summary summary;//优化信息
Solve(options, &problem, &summary);//求解
std::cout << summary.BriefReport() << "\n";//输出优化的简要信息
//最终结果
std::cout << "x : " << initial_x
<< " -> " << x << "\n";
return0;
}
第⼀部分:构造代价函数结构体
这⾥的使⽤了仿函数的技巧,即在CostFunction结构体内,对()进⾏重载,这样的话,该结构体的⼀个实例就能具有类似⼀个函数的性
质,在代码编写过程中就能当做⼀个函数⼀样来使⽤。
关于仿函数,这⾥再多说⼏句,对结构体、类的⼀个实例,⽐如Myclass类的⼀个实例Obj1,如果Myclass⾥对()进⾏了重载,那Obj1
被创建之后,就可以将Obj1这个实例当做函数来⽤,⽐如Obj(x)这样,为了⽅便读者理解,下⾯随便编⼀段简单的⽰例代码,凑活看看
吧。
//仿函数的⽰例代码
#include<iostream>
using namespace std;
class Myclass
{
public:
Myclass(int x):_x(x){};
int operator()(const int n)const{
return n*_x;
}
private:
int _x;
};
int main()
{
Myclass Obj1(5);
cout<<Obj1(3)<<endl;
system("pau");
return0;
}
在我随便写的⽰教代码中,可以看到我将Myclass的()符号的功能定义成了将括号内的数n乘以隐藏参数x倍,其中x是Obj1对象的⼀个私有成员变量,是是在构造Obj1时候赋予的。因为重载了()符号,所以在主函数中Obj1这个对象就可以当做⼀个函数来使⽤,使⽤⽅法为Obj1(n),如果Obj1的内部成员变量_x是5,则此函数功能就是将输⼊参数扩⼤5倍,如果这个成员变量是50,Obj1()函数的功能就是将输⼊n扩⼤50倍,这也是仿函数技巧的⼀个优点,它能利⽤对象的成员变量来储存更多的函数内部参数。
了解了仿函数技巧的使⽤⽅法后,再回过头来看看ceres使⽤中构造CostFuction 的具体⽅法:
CostFunction结构体中,对括号符号重载的函数中,传⼊参数有两个,⼀个是待优化的变量x,另⼀个是残差residual,也就是代价函数的输出。
重载了()符号之后,CostFunction就可以传⼊AutoDiffCostFunction⽅法来构建寻优问题了。
第⼆部分:通过代价函数构建待求解的优化问题
Problem problem;
CostFunction* cost_function =
new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
problem.AddResidualBlock(cost_function, NULL, &x);
这⼀部分就是待求解的优化问题的构建过程,使⽤之前结构体创建⼀个实例,由于使⽤了仿函数技巧,该实例在使⽤上可以当做⼀个函数。基于该实例new了⼀个CostFunction结构体,这⾥使⽤的⾃动求导,将之前的代价函数结构体传⼊,第⼀个1是输出维度,即残差的维度,第⼆个1是输⼊维度,即待寻优参数x的维度。分别对应之前结构体中的residual和x。
向问题中添加误差项,本问题⽐较简单,添加⼀次就⾏(有的问题要不断多次添加ResidualBlock以构建最⼩⼆乘求解问题)。这⾥的参数NULL是指不使⽤核函数,&x表⽰x是待寻优参数。
第三部分:配置问题并求解问题
Solver::Options options;
options.linear_solver_type = ceres::DENSE_QR;
options.minimizer_progress_to_stdout = true;
Solver::Summary summary;
Solve(options, &problem, &summary);
std::cout << summary.BriefReport() << "\n";
std::cout << "x : " << initial_x
<< " -> " << x << "\n";
这⼀部分很好理解,创建⼀个Option,配置⼀下求解器的配置,创建⼀个Summary。最后调⽤Solve⽅法,求解。
最后输出结果:
iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time
04.512500e+010.00e+009.50e+000.00e+000.00e+001.00e+040 5.33e-04 3.46e-03
14.511598e-07 4.51e+019.50e-049.50e+00 1.00e+003.00e+041 5.00e-04 4.05e-03
25.012552e-16 4.51e-07 3.17e-089.50e-04 1.00e+009.00e+041 1.60e-05 4.09e-03
Ceres Solver Report: Iterations: 2, Initial cost: 4.512500e+01, Final cost: 5.012552e-16, Termination: CONVERGENCE
x : 0.5 -> 10
读者们看到这⾥相信已经对Ceres库的使⽤已经有了⼀个⼤概的认识,现在可以试着将代码实际运⾏⼀下来感受⼀下,加深⼀下理解。
博主的使⽤环境为Ubuntu 16.04,所以在此附上,怎么样,是不是很贴⼼:)。
附:代码:
//:
cmake_minimum_required(VERSION 2.8)
project(ceres)
find_package(Ceres REQUIRED)
include_directories(${CERES_INCLUDE_DIRS})
add_executable(u_ceres u_ceres.cpp)
target_link_libraries(u_ceres ${CERES_LIBRARIES})
进阶-更多的求导法
在上⾯的例⼦中,使⽤的是⾃动求导法(AutoDiffCostFunction),Ceres库中其实还有更多的求导⽅法可供选择(虽然⾃动求导的确是最省⼼的,⽽且⼀般情况下也是最快的。。。)。这⾥就简要介绍⼀下其他的求导⽅法:
数值求导法(⼀般⽐⾃动求导法收敛更慢,且更容易出现数值错误):
数值求导法的代价函数结构体构建和⾃动求导中的没有区别,只是在第⼆部分的构建求解问题中稍有区别,下⾯是官⽹给出的数值求导法的问题构建部分代码:
CostFunction* cost_function =
new NumericDiffCostFunction<NumericDiffCostFunctor, ceres::CENTRAL, 1, 1>(
new NumericDiffCostFunctor);
problem.AddResidualBlock(cost_function, NULL, &x);
乍⼀看和⾃动求导法中的代码没区别,除了代价函数结构体的名字定义得稍有不同,使⽤的是Numeri
cDiffCostFunction⽽⾮AutoDiffCostFunction,改动的地⽅只有在模板参数设置输⼊输出维度前⾯加了⼀个模板参数ceres::CENTRAL,表明使⽤的是数值求导法。
还有其他⼀些更多更复杂的求导法,不详述。
再进阶-曲线拟合
趁热打铁,阅读到这⾥想必读者们应该已经对Ceres库的使⽤已经⽐较了解了(如果前⾯认真看了的话),现在就来尝试解决⼀个更加复杂的问题来检验⼀下成果,顺便进阶⼀下。
问题:
拟合⾮线性函数的曲线(和官⽹上的例⼦不⼀样,稍微复杂⼀丢丢):
依然,先上代码:
代码之前先啰嗦⼏句,整个代码的思路还是先构建代价函数结构体,然后在[0,1]之间均匀⽣成待拟合曲线的1000个数据点,加上⽅差为1的⽩噪声,数据点⽤两个vector储存(x_data和y_data),然后构建待求解优化问题,最后求解,拟合曲线参数。
(PS. 本段代码中使⽤OpenCV的随机数产⽣器,要跑代码的同学可能要先装⼀下OpenCV)
#include<iostream>
#include<opencv2/core/core.hpp>
#include<ceres/ceres.h>
using namespace std;
using namespace cv;
//构建代价函数结构体,abc为待优化参数,residual为残差。
struct CURVE_FITTING_COST
{
CURVE_FITTING_COST(double x,double y):_x(x),_y(y){}
template <typename T>
bool operator()(const T* const abc,T* residual)const
{
residual[0]=_y-ceres::exp(abc[0]*_x*_x+abc[1]*_x+abc[2]);
return true;
}
const double _x,_y;
};
//主函数
int main()
{
//参数初始化设置,abc初始化为0,⽩噪声⽅差为1(使⽤OpenCV的随机数产⽣器)。double a=3,b=2,c=1;
double w=1;
RNG rng;
double abc[3]={0,0,0};
//⽣成待拟合曲线的数据散点,储存在Vector⾥,x_data,y_data。
vector<double> x_data,y_data;
for(int i=0;i<1000;i++)
{
double x=i/1000.0;
x_data.push_back(x);
y_data.push_back(exp(a*x*x+b*x+c)+rng.gaussian(w));
}
//反复使⽤AddResidualBlock⽅法(逐个散点,反复1000次)
/
/将每个点的残差累计求和构建最⼩⼆乘优化式
//不使⽤核函数,待优化参数是abc
ceres::Problem problem;
for(int i=0;i<1000;i++)
{
problem.AddResidualBlock(
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_QR;
options.minimizer_progress_to_stdout=true;
ceres::Solver::Summary summary;
ceres::Solve(options,&problem,&summary);
cout<<"a= "<<abc[0]<<endl;
cout<<"b= "<<abc[1]<<endl;
cout<<"c= "<<abc[2]<<endl;
return0;
}
}