Ceres 自动求导（AutomaticDerivatives）进阶

  ceres中有三种求导的方式，分别是自动求导、数值导数和解析导数，本节详细讲述自动求导的使用方式。

一、问题定义

  考虑一个直线拟合问题：
$$y=\frac{b_1}{(1+e^{b_2-b_{3}x}{)}^{1/b_4}}$$
我们给定一些数据 { x i , y i } ,   ∀ i = 1 , . . . , n \{x_i, y_i\},\ \forall i=1,... ,n {xi​,yi​}, ∀i=1,...,n来求解参数$b_1,b_2,b_3$和$b_4$，这个问题可以等价为，需要寻找一组参数$b_1,b_2,b_3$和$b_4$，使得如下方程的数值最小：
E ( b 1 , b 2 , b 3 , b 4 ) = ∑ i f 2 ( b 1 , b 2 , b 3 , b 4 ; x i , y i ) = ∑ i ( b 1 ( 1 + e b 2 − b 3 x i ) 1 / b 4 − y i ) 2 E(b1,b2,b3,b4)=∑if2(b1,b2,b3,b4;xi,yi)=∑i(b1(1+eb2−b3xi)1/b4−yi)2

E(b1​,b2​,b3​,b4​)​=i∑​f2(b1​,b2​,b3​,b4​;xi​,yi​)=i∑​((1+eb2​−b3​xi​)1/b4​b1​​−yi​)2​
为了用ceres来解决这个问题，我们需要定义一个CostFunction(损失函数)来计算给定的$x$和$y$的残差$f$和导数。由于是采用的自动求导的方式，因此我们的代价函数直接用原函数即可。

二、代价函数定义

struct CostFunctor
{
    CostFunctor (const double x,const double y): x_(x),y_(y){}
    template <typename T>
    bool operator()(const T* parameters, T* residuals) const{
        const T b1=parameters[0];
        const T b2=parameters[1];
        const T b3=parameters[2];
        const T b4=parameters[3];
        residuals[0]=b1*pow(1.0+exp(b2-b3*x_),-1.0/b4)-y_;
        return true;
    }

    private:
        const double x_;
        const double y_;
};
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

三、完整代码

#include 
#include 
#include 

struct CostFunctor
{
    CostFunctor (const double x,const double y): x_(x),y_(y){}
    template <typename T>
    bool operator()(const T* parameters, T* residuals) const{
        const T b1=parameters[0];
        const T b2=parameters[1];
        const T b3=parameters[2];
        const T b4=parameters[3];
        residuals[0]=b1*pow(1.0+exp(b2-b3*x_),-1.0/b4)-y_;
        return true;
    }

    private:
        const double x_;
        const double y_;
};


int main(int argc, char const *argv[])
{
    /*生成观测数据及初始化参数*/
    double B1=1.0,B2= 2.0,B3=2.0,B4=2.0;//true value
    double b1=2.0,b2=-1.0,b3=0.0,b4=1.5;//estimate value
    int N=200;//num of data
    double w_sigma=0.02;//sigma of noise
    /*此处仅用opencv库添加随机噪声，如果有其余添加噪声的方法，可以不用调用opencv*/
    cv::RNG rng((unsigned)time(NULL));//opencv random generator
    std::vector<double> x_data,y_data;
    for (int i = 0; i < N; i++)
    {
        double X=i;
        x_data.push_back(X);
        y_data.push_back(B1*pow(1+exp(B2-B3*X),-1/B4)+rng.gaussian(w_sigma*w_sigma));
    }

    /*定义需要优化的变量，并赋予初值*/
    double parameters[4]={b1,b2,b3,b4};
    const double parameters_init[4]={b1,b2,b3,b4};//赋予初值


    /*构建问题*/
    ceres::Problem problem;
    /*建立优化方程*/
    ceres::LossFunction* loss_function=new ceres::HuberLoss(1.0);

    for (int i = 0; i < N; i++)
    {
        //add the AddResidualBlock
        problem.AddResidualBlock(
            /*autoDiff:error type, output dim, input dim*/
            new ceres::AutoDiffCostFunction<CostFunctor,1,4>(new CostFunctor(x_data[i],y_data[i])),
            nullptr,
            parameters
        );
    }
    
    problem.SetParameterLowerBound(parameters,0,-3.5);
    problem.SetParameterUpperBound(parameters,0, 3.5);
    problem.SetParameterLowerBound(parameters,1,-3.5);
    problem.SetParameterUpperBound(parameters,1, 3.5);
    problem.SetParameterLowerBound(parameters,2,-3.5);
    problem.SetParameterUpperBound(parameters,2, 3.5);
    problem.SetParameterLowerBound(parameters,3,-3.5);
    problem.SetParameterUpperBound(parameters,3, 3.5);

    /*运行优化器*/
    ceres::Solver::Options options;
    options.linear_solver_type=ceres::DENSE_QR;
    options.max_num_iterations=5e2;
    options.function_tolerance=1e-10;
    options.minimizer_progress_to_stdout=true;
    ceres::Solver::Summary summary;
    ceres::Solve(options,&problem,&summary);
    std::cout<<summary.BriefReport()<<std::endl;
    std::cout << "parameter_true : " << B1
                                   <<" "<<B2
                                   <<" "<<B3
                                   <<" "<<B4
                                   <<" " <<std::endl;

    std::cout << "parameter_init : " << parameters_init[0]
                                    <<" "<<parameters_init[1]
                                    <<" "<<parameters_init[2]
                                    <<" "<<parameters_init[3]
                                    <<" " <<std::endl;

    std::cout << "parameter_final : " << parameters[0]
                                    <<" "<<parameters[1]
                                    <<" "<<parameters[2]
                                    <<" "<<parameters[3]
                                    <<" " <<std::endl;

    std::cout << "parameter_error : " << parameters[0]-B1
                                    <<" "<<parameters[1]-B2
                                    <<" "<<parameters[2]-B3
                                    <<" "<<parameters[3]-B4
                                    <<" " <<std::endl;
    
    return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

四、输出结果

parameter_true : 1 2 2 2 
parameter_init : 2 -1 0 1.5 
parameter_final : 0.999979 1.98831 1.99625 1.99256 
parameter_error : -2.13889e-05 -0.0116935 -0.0037504 -0.00744402 
1
2
3
4

五、配置文件

cmake_minimum_required(VERSION 2.8)
project(ceresCurveFitting)
set(CMAKE_CXX_STANDARD 14)
find_package(OpenCV 3 REQUIRED )
find_package(Ceres REQUIRED )
include_directories(${OpenCV_INCLUDE_DIRS})
include_directories( ${CERES_INCLUDE_DIRS})
include_directories("/usr/include/eigen3")

add_executable(AutomaticDerivatives AutomaticDerivatives.cpp)
target_link_libraries(AutomaticDerivatives ${CERES_LIBRARIES} ${OpenCV_LIBS})
1
2
3
4
5
6
7
8
9
10
11