ceres学习三（导数）

最新推荐文章于 2023-04-01 16:25:44 发布

磊磊哈哈

最新推荐文章于 2023-04-01 16:25:44 发布

阅读量233

点赞数

原文链接：http://ceres-solver.org/numerical_derivatives.html

版权

SLAM十四讲专栏收录该内容

37 篇文章 8 订阅

订阅专栏

求导函数

优先选用自动微分算法，某些情况可能需要用到解析微分算法，尽量避免数值微分算法

自动微分Automatic Derivatives

struct CostFunctor {
   template <typename T>
   bool operator()(const T* const x, T* residual) const {
     residual[0] = T(10.0) - x[0];
     return true;
   }
};
Problem problem;
  CostFunction* cost_function =
      new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor); //使用自动求导，将之前的代价函数结构体传入，第一个1是输出维度，即残差的维度，第二个1是输入维度，即待寻优参数x的维度。
  problem.AddResidualBlock(cost_function, NULL, &x); //向问题中添加误差项，本问题比较简单，添加一个就行。

例子： $min_{a,b}=\sum_{i=1}^{N}||y_i-(ax_i+b)||^2$

#include <ceres/ceres.h>
#include <opencv/cv.hpp>
using namespace std;
using namespace ceres;
using namespace cv;
//y=ax+b+w
//构建代价函数
struct costfunction
{
  costfunction(double x,double y): _x (x),_y (y){}
  template <typename T> bool operator()(const T* const ab,T* residual)const
  {
    residual[0]=-T(_y)+(ab[0]*T(_x)+ab[1]);
   return true; 
  }
  const double _x,_y;
};
int main(int argc, const char** argv) {
    
    vector<double> x_data,y_data;
    int N=100;//数据个数
    double w_sigma=0.1;//高斯模糊
    cv::RNG rng;//opencv随机数
    double ab[2]={0,0};
    double a=1.0,b=0;//a,b真实值

    //构建观测数据
    for(int i=0;i<N;i++){
      x_data.push_back(i);
      y_data.push_back(a*i+b+rng.gaussian(w_sigma));
    }
	//构建问题
    Problem problem;
    for(int i=0;i<N;i++){
      CostFunction* cost_function=new AutoDiffCostFunction<costfunction,1,2>(new costfunction(x_data[i],y_data[i]));
      problem.AddResidualBlock(cost_function,NULL,ab);
    }
	//求解
    Solver::Options options;
    
    options.minimizer_progress_to_stdout=true;

    Solver::Summary summary;
    Solve(options,&problem,&summary);
    cout<<summary.BriefReport()<<endl;
    cout<<"ab:"<<ab[0]<<ab[1]<<endl;
    return 0;
}

数值法求导 Numeric Derivatives
在某些情况下，无法定义模板化的成本函子，例如，当残差的评估涉及对您无法控制的库函数的调用时。在这种情况下，可以使用数值微分，
无法通过解析或自动微分计算导数时，应使用数字微分。

struct NumericDiffCostFunctor {
  bool operator()(const double*  x, double* residual) const {//这里不使用模板
    residual[0] = 10.0 - x[0];
    return true;
  }
};
CostFunction* cost_function =
  new NumericDiffCostFunction<NumericDiffCostFunctor, ceres::CENTRAL, 1, 1>(
      new NumericDiffCostFunctor);
problem.AddResidualBlock(cost_function, nullptr, &x);

例子：

#include <ceres/ceres.h>
#include <opencv/cv.hpp>
using namespace std;
using namespace ceres;
using namespace cv;
//y-ax-b
struct costfunction
{
  costfunction(double x,double y): _x (x),_y (y){}
  bool operator()(const double* const ab,double* residual)const{//这里不使用模板
    residual[0]=(ab[0]*_x)+ab[1]-_y;
    return true;
  };
  const double _x,_y;
};



int main(int argc, const char** argv) {
    
  double a=1.0, b=0;         // 真实参数值
  int N=100;                          // 数据点
  double w_sigma=0.1;                 // 噪声Sigma值
  cv::RNG rng;                        // OpenCV随机数产生器
  double ab[2] = {0,0};            // abc参数的估计值

  vector<double> x_data, y_data;      // 数据
  cout<<"generating data: "<<endl;
  for ( int i=0; i<N; i++ )
  {
      double x = i/100.0;
      x_data.push_back ( x );
      y_data.push_back (a*x + b + rng.gaussian ( w_sigma ));
      cout<<x_data[i]<<" "<<y_data[i]<<endl;
  }
  Problem problem;
  for ( int i=0; i<N; i++ ){
    CostFunction* cost_function=new NumericDiffCostFunction<costfunction,CENTRAL,1,2>(new costfunction(x_data[i],y_data[i]));
    problem.AddResidualBlock(cost_function,NULL,ab);
  }
    Solver::Options options;
  options.minimizer_progress_to_stdout = true;
  Solver::Summary summary;
  Solve(options, &problem, &summary);
  std::cout << summary.BriefReport() << "\n";
  std::cout << "output a: " << ab[0]
            << "output b: " << ab[1] << "\n";
    return 0;
}

解析微分法Analytic Derivatives
多元函数
$f(x,y)=xcos\theta +ysin\theta+c\\ 偏导数:\\ \frac{\partial f}{\partial \theta}=-x sin\theta+ycos\theta\\ \frac{\partial f}{\partial c}=1$
什么时候使用解析法：
（1）表达式很简单，例如大部分是线性的。
（2）可以使用诸如Maple，Mathematica或SymPy之类的计算机代数系统来象征性地区分目标函数并生成C ++对其进行评估。
（3）性能是最重要的问题，并且存在代数结构，可以用来获得比自动微分更好的性能。
也就是说，要想从分析衍生工具中获得最佳性能，就需要大量的工作。在走这条路之前，有必要测量花费在评估Jacobian上的时间占总求解时间的一部分，使用Amdahl定律。
（4）复杂的求导函数，其他方法可以计算导数。
没有其他方法可以计算导数


#include "ceres/ceres.h"
#include <opencv2/core/core.hpp>
using namespace std;
using ceres::AutoDiffCostFunction;
using ceres::CostFunction;
using ceres::Problem;
using ceres::Solver;

class AnalyticCostFunctor : public ceres::SizedCostFunction<1//输出的残差维度
                                                        ,2> //第一个参数的纬度大小
                                                        {
   public:
     AnalyticCostFunctor(const double x, const double y) : x_(x), y_(y) {}
     virtual ~AnalyticCostFunctor() {}
     virtual bool Evaluate(double const* const* parameters,
                           double* residuals,
                           double** jacobians) const {
       const double theta = parameters[0][0];//读取第一个参数，的第一个纬度值
       const double c = parameters[0][1];//读取第一个参数，的第二个纬度值
       residuals[0] =cos(theta)*x_+sin(theta)*y_+c;//计算残差

       if (!jacobians) return true;
       double* jacobian = jacobians[0];
       if (!jacobian) return true;

       jacobian[0] = -sin(theta)*x_+cos(theta)*y_;//求雅克比矩阵，即求偏导
       jacobian[1] = 1.0;//求偏导
       return true;
     }
   private:
     const double x_;
     const double y_;
 };
int main(int argc, char** argv)
{
    double theta=-M_PI/4;
    double c1=1;
    double k=-tan(theta);
    double b1=-c1/sin(theta);         // a is k b is b真实参数值
    cout<<"ture theta="<<theta<<"c1="<<c1<<endl;

    int N1=100;                          // 数据点
    double w_sigma=0.01;                 // 噪声Sigma值
    cv::RNG rng;                        // OpenCV随机数产生器
    double thetac[2] = {0,0};            // abc参数的估计值

    vector<double> l1_x_data, l1_y_data;      // 数据
    cout<<"generating data: "<<endl;
    for ( int i=0; i<N1; i++ )
    {
        double x = i/100.0;
        l1_x_data.push_back ( x );
        l1_y_data.push_back (k*x + b1 + rng.gaussian ( w_sigma ));
     //   cout<<x_data[i]<<" "<<y_data[i]<<endl;
    }
    // Build the problem.
    Problem problem;
    for ( int i=0; i<N1; i++ )
    {
        CostFunction* cost_function =new AnalyticCostFunctor( l1_x_data[i], l1_y_data[i] );//这里直接导入参数
        problem.AddResidualBlock(cost_function, NULL, thetac);
    }

    // Run the solver!
    Solver::Options options;
    options.minimizer_progress_to_stdout = true;
    Solver::Summary summary;
    Solve(options, &problem, &summary);
  //  std::cout << summary.BriefReport() << "\n";
    double opt_theta =thetac[0];
    double opt_c1 =thetac[1];
    std::cout << "opt_theta=: " << opt_theta
              << "opt_c1 =: " << opt_c1;
    return 0;
}

磊磊哈哈

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
ceres学习三（导数）

求导函数优先选用自动微分算法，某些情况可能需要用到解析微分算法，尽量避免数值微分算法自动微分Automatic Derivativesstruct CostFunctor { template <typename T> bool operator()(const T* const x, T* residual) const { residual[0] = T(10.0) - x[0]; return true; }};Problem probl
复制链接

扫一扫

专栏目录