求导函数
优先选用自动微分算法,某些情况可能需要用到解析微分算法,尽量避免数值微分算法
- 自动微分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); //向问题中添加误差项,本问题比较简单,添加一个就行。
例子: m i n a , b = ∑ i = 1 N ∣ ∣ y i − ( a x i + b ) ∣ ∣ 2 min_{a,b}=\sum_{i=1}^{N}||y_i-(ax_i+b)||^2 mina,b=i=1∑N∣∣yi−(axi+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 ) = x c o s θ + y s i n θ + c 偏 导 数 : ∂ f ∂ θ = − x s i n θ + y c o s θ ∂ f ∂ c = 1 f(x,y)=xcos\theta +ysin\theta+c\\ 偏导数:\\ \frac{\partial f}{\partial \theta}=-x sin\theta+ycos\theta\\ \frac{\partial f}{\partial c}=1 f(x,y)=xcosθ+ysinθ+c偏导数:∂θ∂f=−xsinθ+ycosθ∂c∂f=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;
}