Ceres 数值导数 （NumericDerivatives）进阶

  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​)​=i∑​f2(b1​,b2​,b3​,b4​;xi​,yi​)=i∑​((1+eb2​−b3​xi​)1/b4​b1​​−yi​)2​
为了用ceres来解决这个问题，我们需要定义一个CostFunction(代价函数)来计算给定的$x$和$y$的残差$f$和导数。利用数值导数的形式对这个代价函数进行求解，求解的方式有前向查分、中值差分和Ridder差分。

二、代价函数定义

1. Forward Differences

  我们不能在计算机上进行数值极限运算，所以我们做的事情，就是选择一个很小的值，并将导数近似为:
$$Df(x)\approx \frac{f(x+h)-f(x)}{h}$$
上面的公式是最简单最基本的数值微分形式。它被称为前向差分公式。我们定义代价函数，并用前向差分的方式进行求解，如下：

struct Rat43CostFunctor {
  Rat43CostFunctor(const double x, const double y) : x_(x), y_(y) {}

  bool operator()(const double* parameters, double* residuals) const {
    const double b1 = parameters[0];
    const double b2 = parameters[1];
    const double b3 = parameters[2];
    const double b4 = parameters[3];
    residuals[0] = b1 * pow(1.0 + exp(b2 -  b3 * x_), -1.0 / b4) - y_;
    return true;
  }

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

ceres::CostFunction* cost_function =
    new ceres::NumericDiffCostFunction<Rat43CostFunctor, ceres::FORWARD, 1, 4>(
        new Rat43CostFunctor(x_data[i],y_data[i]));
1
2
3

2. Central Differences

  中值差分的表达式如下：
$$Df(x)\approx \frac{f(x+h)-f(x-h)}{2h}$$
定义如下：

ceres::CostFunction* cost_function =
    new ceres::NumericDiffCostFunction<Rat43CostFunctor, ceres::CENTRAL, 1, 4>(
        new Rat43CostFunctor(x_data[i],y_data[i]));
1
2
3

3. Ridders Differences

  Ridders差分的表达式如下：
D f ( x ) = f ( x + h ) − f ( x − h ) 2 h + h 2 3 ! D 3 f ( x ) + h 4 5 ! D 5 f ( x ) + ⋯ = f ( x + h ) − f ( x − h ) 2 h + K 2 h 2 + K 4 h 4 + ⋯

Df(x)​=2hf(x+h)−f(x−h)​+3!h2​D3f(x)+5!h4​D5f(x)+⋯=2hf(x+h)−f(x−h)​+K2​h2+K4​h4+⋯​
定义如下：

ceres::CostFunction* cost_function =
    new ceres::NumericDiffCostFunction<Rat43CostFunctor, ceres::RIDDERS, 1, 4>(
        new Rat43CostFunctor(x_data[i],y_data[i]));
1
2
3

三、完整代码

#include 

#include "ceres/ceres.h"
#include "glog/logging.h"

#include 


struct Rat43CostFunctor {
  Rat43CostFunctor(const double x, const double y) : x_(x), y_(y) {}

  bool operator()(const double* parameters, double* residuals) const {
    const double b1 = parameters[0];
    const double b2 = parameters[1];
    const double b3 = parameters[2];
    const double b4 = parameters[3];
    residuals[0] = b1 * pow(1.0 + exp(b2 -  b3 * x_), -1.0 / b4) - y_;
    return true;
  }

  const double x_;
  const double y_;
};



int main(int argc, char** 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=20;//num of data

    double w_sigma=0.02;//sigma of noise

    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};//give initial value
  const double parameters_init[4] ={b1,b2,b3,b4};//give initial value
  // 构建 problem.
  ceres::Problem problem;

  // 建立优化方程.
  ceres::LossFunction* loss_function = new ceres::HuberLoss(1.0);
  for (int i = 0; i < N; i++)
  {
    ceres::CostFunction* cost_function =
        new ceres::NumericDiffCostFunction<Rat43CostFunctor, ceres::FORWARD, 1, 4>(
            new Rat43CostFunctor(x_data[i],y_data[i]));
    // ceres::CostFunction* cost_function = new Rat43Analytic(x_data[i],y_data[i]);
    problem.AddResidualBlock(cost_function, loss_function, 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() << "\n";
                            
  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
106

四、输出结果

parameter_true : 1 2 2 2 
parameter_init : 2 -1 0 1.5 
parameter_final : 1.00002 2.0147 2.00489 2.01526 
parameter_error : 1.51419e-05 0.0147027 0.0048945 0.0152647 
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(NumericDerivatives NumericDerivatives.cpp)
target_link_libraries(NumericDerivatives ${CERES_LIBRARIES} ${OpenCV_LIBS})
1
2
3
4
5
6
7
8
9
10

六、适用情况

  当无法分析或使用自动微分计算导数时，应使用数值微分。通常情况下，当您用一个外部库或函数时，我们不知道它的分析形式，或者即使知道，也无法以使用自动导数所需的方式重新编写它。当使用数值微分时，至少使用中值差分，如果执行时间不是问题，或者目标函数难以确定良好的静态相对步长，建议使用Ridders方法。