Ceres的使用

淩亂°似流年 2022-02-25 05:38 275阅读 0赞

[Ceres官网][Ceres]

### 文章目录 ###

*   *  Ceres介绍
     *  求解过程
     *   *  举例叙述
         *  求导方法
     *  曲线拟合
     *   *  举例
         *  outliers的处理

## Ceres介绍 ##

*  ceres 可以解**约束非线性最小二乘**问题：  
    ![在这里插入图片描述][20190331222705630.png_pic_center]

1.   f i f\_i fi：cost function 代价函数，SLAM中叫做误差项
2.   x j x\_j xj：优化变量
3.   ρ i \\rho\_\{i\} ρi： 核函数

![在这里插入图片描述][20190401222053334.png]

在最简单的情况下：核函数 ρ i = x    \\rho\_i=x\\space\\space ρi=x  约束条件： l j = − ∞    u j = ∞ l\_j=-\\infty \\space\\space u\_j=\\infty lj=−∞  uj=∞ ，我们得到最熟悉的非线性最小二乘问题：  
![在这里插入图片描述][20190331223555620.png_pic_center]

--------------------

## 求解过程 ##

一般利用Ceres求解非线性优化问题分为三个过程：

1.  构建代价函数(cost function)
2.  通过代价函数来构建优化问题
3.  配置求解器

### 举例叙述 ###

问题：求解  1 2 ( 10 − x ) 2 \\frac\{1\}\{2\}(10-x)^2 21(10−x)2的最小值时候的 x x x

1.  构建代价函数

struct CostFunctor
    { 
        template <typename T>
        bool operator()(const T* const x, T* residual) const
        { 
            residual[0] = T(10.0) - x[0];
            return true;
        }
    };

residual：残差（一维）  
x：模型参数（一维）

1.  构建优化问题

double initial_x = 5.0;
        double x = initial_x; //变量初始值
    
    
        ceres::Problem  problem;
        CostFunction* cost_function = new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
        problem.AddResidualBlock(cost_function, NULL, &x);

(1) 先定义一个求解问题：`ceres::Problem problem`  
(2) 使用自动求导，将之前定义的cost function传入，模板参数：误差类型（自定义的仿函数cost function）、输出维度（残差）、输入维度（变量）；顺序与cost function定义相反  
(3)　添加误差项：　（２）中的cost\_function的指针，核函数（NULL不使用），x优化变量

1.  配置求解器

Solver::Options options;
        options.linear_solver_type=DENSE_QR;　　//增量方程如何求解
        options.minimizer_progress_to_stdout=true;　//迭代过程信息输出到stdout
    
    
        Solver::Summary summary; //优化信息
        Solve(options, &problem, &summary);　//开始执行求解

### 求导方法 ###

*  Automatic Derivatives(自动求导)

效率比数值求导高

*  Numeric Derivatives(数值求导)
 *  Analytic Derivatives(解析求导)

--------------------

## 曲线拟合 ##

### 举例 ###

*  拟合 y = e x p ( a x 2 + b x + c ) + w y=exp(ax^2+bx+c)+w y=exp(ax2\+bx\+c)\+w  
     a   b   c a \\space b \\space c a b c是曲线参数， w w w是高斯噪声，假设我们有Ｎ（１００）个带噪声的数据点，拟合曲线，求曲线参数：  
     min ⁡ a b c 1 2 ∑ i = 1 N ∣ ∣ y i − e x p ( a x i 2 + b x i + c ) ∣ ∣ 2 \\min \\limits\_\{abc\} \\frac\{1\}\{2\}\\sum\_\{i=1\}^\{N\}||y\_i-exp(ax\_i^2+bx\_i+c)||^2 abcmin21i=1∑N∣∣yi−exp(axi2\+bxi\+c)∣∣2

1.  定义代价函数

#include <iostream>
    #include <opencv4/opencv2/core.hpp>
    #include <ceres/ceres.h>
    #include <chrono>
    #include <fstream>
    
    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, T* residual) const
        { 
            residual[0]=T(_y) - ceres::exp(abc[0]*T(_x) * T(_x) + abc[1]*T(_x) +abc[2]);
            return true;
        }
    
        const double _x, _y;
    };
    
    
    
    
    
    int main()
    { 
        string fp("../res/data.txt");
        ofstream outfile;
        outfile.open(fp,ios::trunc);
    
        double a=1.0, b=2.0, c=1.0; //真实的参数
        int N=100 ;  //数据点
        double w_sigma = 1.0;
        cv::RNG rng;
        double abc[3] = { 0,0,0};
    
    
    
    
        // 数据生成
        vector<double> x_data, y_data;
        cout << "gernerating data" << endl;
    
        for(int i=0; i<N; i++)
        { 
            double x = i/100.0, y=0;
            x_data.push_back(x);
            y = exp(a*x*x+b*x+c)+rng.gaussian(w_sigma);
            y_data.push_back(y);
    
            outfile << x << "\t"<< y << endl;
    
            cout << x_data[i] << " " << y_data[i] << endl;
        }
    
    
    
    
    
        // 构建最小二乘问题
        ceres::Problem problem;
        for(int i=0; i<N; 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;
        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()  << "miliseconds" << endl;
    
        cout << summary.BriefReport() << endl;
        cout << "estimate a,b,c = " ;
        for(auto a:abc)cout << a<< " ";
        cout << endl;
    
        outfile << abc[0] << "\t"<< abc[1] <<  "\t" << abc[2] << endl;
        outfile.close();
        return 0;
    }

![在这里插入图片描述][20190402214326865.png_pic_center]  
估计的参数 a:0.891943   b:2.17039   c:0.944142  
原始参数：a: 1   b: 2   b:1

### outliers的处理 ###

*  假如有少数极端值（outliers）,那么对参数的估计有很大影响，如：

![在这里插入图片描述][20190402215541871.png_pic_center]

这时候可以借助核函数（loss function）来解决，将

problem.AddResidualBlock(new ceres::AutoDiffCostFunction<CURVE_FITTING_COST, 1, 3>(new CURVE_FITTING_COST(x_data[i], y_data[i])), nullptr, abc);

nullptr　改为　new ceres::CauchyLoss(0.5)

problem.AddResidualBlock(new ceres::AutoDiffCostFunction<CURVE_FITTING_COST, 1, 3>(new CURVE_FITTING_COST(x_data[i], y_data[i])), new ceres::CauchyLoss(0.5), abc);

![在这里插入图片描述][2019040222101867.png_pic_center]

[Ceres]: http://www.ceres-solver.org/
[20190331222705630.png_pic_center]: /images/20220225/bbc8a45f53eb4a1fabfcec9b223a8ea3.png
[20190401222053334.png]: /images/20220225/21df509cc31546799850d0df34feb2d1.png
[20190331223555620.png_pic_center]: /images/20220225/092106c468db44e9ba7088617cdca618.png
[20190402214326865.png_pic_center]: /images/20220225/32225d06e0744d73bbd52e0bc4c56b1c.png
[20190402215541871.png_pic_center]: /images/20220225/e196c2cdcc564cc789a32017ff0d5be5.png
[2019040222101867.png_pic_center]: /images/20220225/1fba02e2f1a7486aa02d65d875e2ea46.png