Ceres-Solver 官网教程翻译与学习

最新推荐文章于 2023-08-25 12:31:41 发布

guchaono1

最新推荐文章于 2023-08-25 12:31:41 发布

阅读量1.4k

点赞数 2

分类专栏： c++ 文章标签： c++

本文链接：https://blog.csdn.net/guchaono1/article/details/116780172

版权

c++ 专栏收录该内容

13 篇文章 2 订阅

订阅专栏

ceres-solver 参考1这是一个系列教程
ceres-solver 参考2这是一个入门教程
官网教程推荐

仿函数（预备知识）

仿函数参考
ceres当中代价函数（cost function）的构建当中用到了仿函数（functor），使用方法是在costfunction函数中重载了（）运算符（operator）
示例程序：

class Functor 
{
private:
	char* ss;
	int count;
public:
	explicit Functor(char* str) :ss(str),count(0)
	{
		std::cout << "Functor: " << ss << "\n"; 
	}
	void operator() (std::string str)  {
		count++;
		std::cout << "The Functor is called " << count << " times " <<str<< "\n";
	}
};
int main()
{
	Functor functor("construct");
	functor("first");
	functor("second");
	Functor functor1("construct1");
	functor1("first1");
	functor1("second1");
}

结果：

Functor: construct
The Functor is called 1 times first
The Functor is called 2 times second
Functor: construct1
The Functor is called 1 times first1
The Functor is called 2 times second1

Ceres

ceres用来解决边界约束鲁棒非线性最小二乘问题（bounds constrained robustified non-linear least squares problems）

最小二乘问题

边界约束非线性最小二成问题

$\begin{aligned} \underset{x}{\min}&\enspace\frac{1}{2}{\sum\limits_i}\rho_i\Big(\|f_i(x_{i1},\dots,x_{ik})\|^2\Big) \\\mathrm{s.t.}&\enspace{}l_j\le{}x_j\le{}u_j \end{aligned}\tag{1}$

目标函数 $\min\limits_x\thinspace{}f(x)$ 表示使 $f (x)$ 最小的 $x$ 取值
约束条件 $\mathrm{s.t.}$ 是subject to 的简写

$\rho_i\Big(\|f_i(x_{i1},\dots,x_{ik})\|^2\Big)$ 称为残差块(ResidualBlock)
$f_i(\cdot)$ 称为代价函数(CostFuction)
$[x_{i1},\dots,x_{ik}]$ 称为参数块(ParameterBlock) 其中 $x_{ik}$ 是标量(Scalar),参数块可以是一个标量群，也可以是单个标量
$l_j$ , $u_j$ 是参数块(Parameter Block) $x_j$ 的边界条件
$\rho_i$ 称为损失函数(LossFunction) 是一个标量函数(Scalar Function)，用来减小异常值（outliers）对结果的非线性最小二乘求解结果的影响。当 $\rho_i\big(f_i(x_i)\big)=f_i(x_i)$ 且 $l_i=-\infty$ , $u_i=\infty$ 时，我们就得到通常意义上的非线性最小二乘问题
$\frac{1}{2}\thinspace\sum\limits_i\|f(x_{i1},\dots,x_{ik})\|^2\tag{2}$

Hello World

首先从一个查找函数最小值的问题起步
$\dfrac{1}{2}(10-x)^2$
这是一个很简单的问题，很容易得到最小值在 $x = 10$ 处，但是这是一个很好的起点来说明使用Ceres来解决问题的基本用法。
首先写一个仿函数（functor）来评价函数
$f (x) = 10 - x$

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

operator()是模板方法(Template Method)，这种方法假设输入与输出都是某个类型T，模板的使用允许ceres按照下面两种方法调用CostFunctor::operator()：

当只需要计算残差值时T=double
当需要Jacobian矩阵时T=Jet，Jet是Ceres定义的一个特殊类型
一旦我们有办法计算残差方程(Residual Function)，是时候构建一个非线性最小二乘问题(non-linear least squares problem)并利用Ceres来求解

void TestCeres()
{
	ceres::Problem problem;
	double initial_x = 5.0;
	double x = initial_x;
	//建立唯一的代价函数（cost function）也可以称为 残差（residual）
	//Auto-Differentiation用来获得微分(derivative	)或者 Jacobian矩阵
	ceres::CostFunction* cost_function =
		new ceres::AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
	problem.AddResidualBlock(cost_function, nullptr, &x);
	//run the solver!
	ceres::Solver::Options options;
	options.linear_solver_type = ceres::DENSE_QR;
	options.minimizer_progress_to_stdout = true;
	ceres::Solver::Summary summary;
	ceres::Solve(options, &problem, &summary);

	std::cout << summary.BriefReport() << "\n";
	std::cout << "x : " << initial_x << " -> " << x << "\n";
}

AutoDiffCostFucntion将 CostFunctor的一个实例地址作为输入，自动的对其进行微分，并赋予其一个 CostFunction接口
结果: 我的结果跟官网的结果稍有不同，我的目标值（cost）比官网小，不知道为什么

iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time
   0  1.250000e+01    0.00e+00    5.00e+00   0.00e+00   0.00e+00  1.00e+04        0    8.90e-04    1.96e-03
   1  1.249750e-07    1.25e+01    5.00e-04   5.00e+00   1.00e+00  3.00e+04        1    1.18e-03    7.61e-03
   2  1.388518e-16    1.25e-07    1.67e-08   5.00e-04   1.00e+00  9.00e+04        1    3.60e-04    8.09e-03
Ceres Solver Report: Iterations: 3, Initial cost: 1.250000e+01, Final cost: 1.388518e-16, Termination: CONVERGENCE
x : 5 -> 10

从 $x = 5$ 开始，求解器（solver）经过两次迭代得到最终结果10，仔细的读者应该发现这是一个线性问题，采用线性求解器应该就会求得优化结果。求解器的默认配置是面向非线性问题的，为了简化过程例子中并没有修改配置。使用Ceres-Solver确实有可能通过一次迭代得到这个问题的结果，在讲到Ceres收敛与参数设置时我们会详细的讨论这个问题。

导数（derivatives）

参考
 导数与微分
Ceres-Solver与大多数优化包一样，依赖于能够在任意参数下评估（evaluate）目标函数中每一项的函数值与导数，能够正确并有效的做到这一点对获得好的结果至关重要。Ceres Solver提供了几种不同方法，在上面的例子中你已经见到了其中的一种-Automatic Differentiation（自动微分），我们现在考虑其他两种可能性（possibilities）：解析求导（Analytic Derivatives）与数值求导（Numerical Derivatives）

数值微分（Numerical Derivatives）

在某些情况下，我们无法定义一个模板代价函数，比如计算残差（residual）时调用了一个无法控制的库函数，这种情况下就可以使用数值求导。用户定义一个仿函数（functor）来计算残差（residual），并构建一个NumericDiffCostFunction，以f(x)=10-x为例，对应的仿函数（functor）为：

struct NumericalDiffCostFunctor
{
	bool operator()(const double* const x, double* residual) const
	{
		residual[0] = 10.0 - x[0];
		return true;
	}
};

按照如下方法加入Problem中

void TestCeresNumerical()
{
	ceres::Problem problem;
	double initial_x = 5;
	double x = initial_x;
	ceres::CostFunction* cost_function =
		new ceres::NumericDiffCostFunction<NumericalDiffCostFunctor, ceres::CENTRAL, 1, 1>
		(new NumericalDiffCostFunctor);
	problem.AddResidualBlock(cost_function, nullptr, &x);
	ceres::Solver::Options option;
	option.linear_solver_type = ceres::DENSE_QR;
	option.minimizer_progress_to_stdout = true;
	ceres::Solver::Summary summary;
	ceres::Solve(option, &problem, &summary);

	std::cout << summary.BriefReport() << "\n";
	std::cout << "x: " << initial_x << " -> " << x << "\n";
}

与Automatic Differentiation中我们使用的代码做对比

	ceres::CostFunction* cost_function =
		new ceres::AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
	problem.AddResidualBlock(cost_function, nullptr, &x);

与 Automatic Differentiation中的代码结构几乎一样（identical），但是增加了一个额外的参数ceres::CENTRAL用来表明采用哪种有限差分方案来计算数值导数。ceres::CENTRAL定义如下：NumericDiffCostFunction

enum NumericDiffMethodType {
  // Compute central finite difference: f'(x) ~ (f(x+h) - f(x-h)) / 2h.
  CENTRAL,
  // Compute forward finite difference: f'(x) ~ (f(x+h) - f(x)) / h.
  FORWARD,
  // Adaptive numerical differentiation using Ridders' method. Provides more
  // accurate and robust derivatives at the expense of additional cost
  // function evaluations.
  RIDDERS
};

总的来说我们推荐自动微分（Automatic Differentiation）来取代数值微分（Numerical Differentiation）。C++的模板类让自动微分十分高效，而数值微分计算量大，容易产生数值误差（numeric error），并且收敛（convergence）更慢。

解析法求导（Analytic Derivative）

有些情况下可能无法使用自动微分。比如，在某些情况下使用封闭解比自动微分中使用的链式法则（chain rule）的计算导数更有效。
在这种情况下，需要提供你自己的残差（residual）与雅可比（Jacobian）计算代码。定义一个CostFunction的子类（subclass）,如果你知道参数（parameter）与残差（residual）在编译时的大小也可以定义SizedCostFunction的子类。下面以实现 $f (x) = 10 - x$ 的SimpleCostFunction为例

class QuadraticCostFunction : public ceres::SizedCostFunction<1, 1>
{
public:
	virtual ~QuadraticCostFunction() {}
	virtual bool Evaluate(double const* const* parameters,
		double* residuals,
		double** jacobians) const
	{
		const double x = parameters[0][0];
		residuals[0] = 10 - x;
		//compute the Jacobian if asked for
		if (jacobians != nullptr && jacobians[0] != nullptr)
		{
			jacobians[0][0] = -1;
		}
		return true;
	}
};

下面是我自己写的，教程中没有

void testQuadraticCostFunction()
{
	ceres::Problem problem;
	double initial_x = 5;
	double x = initial_x;
	ceres::CostFunction* cost_function = new QuadraticCostFunction;
	problem.AddResidualBlock(cost_function,nullptr,&x);
	ceres::Solver::Options option;
	option.linear_solver_type = ceres::DENSE_QR;
	option.minimizer_progress_to_stdout = true;
	ceres::Solver::Summary summary;
	ceres::Solve(option, &problem, &summary);
	std::cout << summary.BriefReport() << "\n";
	std::cout << "x: " << initial_x << " -> " << x << '\n';
}

这里代码应该是写错了,将SimpleCOstFunction写成了QuadraticCostFunction
QuadraticCostFunction::Evaluate有一个参数（parameter）数组，一个输出的残差（residual）数组，和一个输出的雅可比（Jacobian）数组，雅可比数组是可选的，evaluate函数会检查jacobians是否为非空（non-null），如果是非空的就填入残差方程的导数值（value of the derivative of residual function）。在本例当中残差方程是线性的，所以雅可比矩阵是常数（constant）
从上面的例子中我们可以看出，实现CostFunction对象的过程有些枯燥，除非你有非常好的理由自己实现Jacobian矩阵，否则我们还是推荐你使用自动微分或者数值微分来构建残差块（residual blocks）
计算导数目前为止是使用Ceres最复杂的部分，根据情况使用者可能需要更复杂的导数计算方法。这一部分只是浅显的介绍了如何向Ceres提供导数，一旦你适应了使用NumericDiffCostFunctionandAutoDiffCostFunction，我们建议看一看DynamicAutoDiffCostFunction, CostFunctionToFunctor, NumericDiffFunctorandConditionedCostFunction来获取更多构建与计算代价函数的方法。

Powell’s Function

鲍威尔方程
现在考虑一个稍微复杂点的例子—鲍威尔方程的最小化，令 $x=[x_1,x_2,x_3,x_4]$
$\begin{aligned} f_1(x)&=x_1+10x_2\\ f_2(x)&=\sqrt{5}(x_3-x_4)\\ f_3(x)&=(x_2-2x_3)^2\\ f_4(x)&=\sqrt{10}(x_1-x_4)^2\\ F(x)&=[f_1(x),f_2(x),f_3(x),f_4(x)] \end{aligned}$
$F (x)$ 是一个拥有四个参数，四组残差的方程，我们希望找到一个令 $\dfrac{1}{2}\|F(x)\|^2$ 最小化的 $x$
第一步，目标函数中的每一项都定义仿函数，下面是 $f_1(x_1,x_2)$ , $f_2(x_3,x_4)$ , $f_3(x_2,x_3)$ , $f_4(x_1,x_4)$ 的代码

struct F1
{
   template <typename T>
   bool operator()(const T* const x1, const T* const x2, T* residual) const
   {
   	residual[0] = x1[0] + 10 * x2[0];
   	return true;
   }
};
struct F2
{
   template <typename T>
   bool operator()(const T* const x3, const T* const x4, T* residual) const
   {
   	residual[0] = sqrt(5) * (x3[0] - x4[0]);
   }
};
struct F3
{
   template <typename T>
   bool operator()(const T* const x2, const T* const x3, T* residual) const
   {
   	residual[0] = (x2[0] - 2 * x3[0]) * (x2[0] - 2 * x3[0]);
   }
};
struct F4 
{
   template <typename T>
   bool operator()(const T* const x1, const T* const x4, T* residual)const
   {
   	residual = sqrt(10.0) * (x1[0] - x4[0]) * (x1[0] - x4[0]);
   	return true;
   }
};

利用上面定义的结构体按如下方法构建problem：

	double x1 = 3.0; double x2 = -1.0; double x3 = 0.0; double x4 = 1.0;
	ceres::Problem problem;
	problem.AddResidualBlock(
		new ceres::AutoDiffCostFunction<F1, 1, 1, 1>(new F1), nullptr, &x1, &x2);
	problem.AddResidualBlock(
		new ceres::AutoDiffCostFunction<F2, 1, 1, 1>(new F2), nullptr, &x3, &x4);
	problem.AddResidualBlock(
		new ceres::AutoDiffCostFunction<F3, 1, 1, 1>(new F3), nullptr, &x2, &x3);
	problem.AddResidualBlock(
		new ceres::AutoDiffCostFunction<F4, 1, 1, 1>(new F4), nullptr, &x1, &x4);
	ceres::Solver::Summary summary;
	ceres::Solver::Options option;
	option.max_num_iterations = 100;
	option.linear_solver_type = ceres::DENSE_QR;
	option.minimizer_progress_to_stdout = true;
	ceres::Solve(option, &problem, &summary);
	std::cout << "Initial: " << "x1 =" << x1 << ", x2= " << x2 << ", x3= " << x3 << ", x4= " << x4 << '\n';
	std::cout << summary.FullReport() << '\n';
	std::cout << "Final: " << "x1 =" << x1 << ", x2= " << x2 << ", x3= " << x3 << ", x4= " << x4 << '\n';

注意到每一个ResidualBlock只依赖两个参数，相应的残差对象也不依赖所有的四个参数。结果略。

Curve Fitting

到目前为止我们所见到的例子都是无数据的简单优化问题。最小方差与非线性最小方差分析的最初目的是对数据进行曲线拟合。现在正适合我们考虑这个问题。问题包含从曲线 $e^{0.3x+0.1}$ 采样的数据，并向数据添加带有标准差（standard deviation） $\sigma=2$ 的高斯（Gaussian）噪声。现在我们用数据拟合曲线 $y=e^{mx+c}$
我们先来定义一个模板类来计算残差。对每一个观测（observation）都有一个残差。

struct ExponentialResidual
{
	ExponentialResidual(double x, double y) :x_(x), y_(y) {}
	template <typename T>
	bool operator()(const T* const m, const T* const c, T* residual) const
	{
		residual[0] = y_ - exp(m[0] * x_ + c[0]);
		return ture;
	}
private:
	const double x_;
	const double y_;
};

假设观测数据存在一个 $2 n$ 大小的数组data当中, problem 的构建就是简单的对每一个观测数据创建一个CostFunction。

    const int kNumObservations = 67;
    // clang-format off
    const double data[] = {  0.000000e+00, 1.133898e+00,
... }；//这里有67组数 一组有一个x，一个y
	double m = 0.0;
	double c = 0.0;
	ceres::Problem problem;
	for (int i = 0; i <kNumObservations; i++)
	{
		ceres::CostFunction* cost_function =
			new ceres::AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
				new ExponentialResidual(data[2 * i], data[2 * i + 1]));
		problem.AddResidualBlock(cost_function, nullptr, &m, &c);
	}
	ceres::Solver::Options option;
	option.max_num_iterations = 100;
	option.linear_solver_type = ceres::DENSE_QR;
	option.minimizer_progress_to_stdout = true;
	ceres::Solver::Summary summary;
	ceres::Solve(option, &problem, &summary);
	std::cout << summary.BriefReport() << '\n';
	std::cout << "Initial m: " << 0.0 << " c: " << 0.0 << '\n';
	std::cout << "Final   m: " << m << " c: " << c << '\n';

iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time
   0  1.211734e+02    0.00e+00    3.61e+02   0.00e+00   0.00e+00  1.00e+04        0    1.80e-03    2.13e-03
   1  2.334822e+03   -2.21e+03    3.61e+02   7.52e-01  -1.87e+01  5.00e+03        1    8.37e-04    3.48e-03
   2  2.331438e+03   -2.21e+03    3.61e+02   7.51e-01  -1.86e+01  1.25e+03        1    4.32e-04    4.01e-03
   3  2.311313e+03   -2.19e+03    3.61e+02   7.48e-01  -1.85e+01  1.56e+02        1    4.24e-04    4.52e-03
   4  2.137268e+03   -2.02e+03    3.61e+02   7.22e-01  -1.70e+01  9.77e+00        1    4.26e-04    5.02e-03
   5  8.553131e+02   -7.34e+02    3.61e+02   5.78e-01  -6.32e+00  3.05e-01        1    4.28e-04    5.54e-03
   6  3.306595e+01    8.81e+01    4.10e+02   3.18e-01   1.37e+00  9.16e-01        1    2.11e-03    7.77e-03
   7  6.426770e+00    2.66e+01    1.81e+02   1.29e-01   1.10e+00  2.75e+00        1    2.01e-03    9.90e-03
   8  3.344546e+00    3.08e+00    5.51e+01   3.05e-02   1.03e+00  8.24e+00        1    2.24e-03    1.22e-02
   9  1.987485e+00    1.36e+00    2.33e+01   8.87e-02   9.94e-01  2.47e+01        1    2.01e-03    1.44e-02
  10  1.211585e+00    7.76e-01    8.22e+00   1.05e-01   9.89e-01  7.42e+01        1    2.38e-03    1.69e-02
  11  1.063265e+00    1.48e-01    1.44e+00   6.06e-02   9.97e-01  2.22e+02        1    2.05e-03    1.91e-02
  12  1.056795e+00    6.47e-03    1.18e-01   1.47e-02   1.00e+00  6.67e+02        1    2.02e-03    2.12e-02
  13  1.056751e+00    4.39e-05    3.79e-03   1.28e-03   1.00e+00  2.00e+03        1    2.02e-03    2.34e-02
Ceres Solver Report: Iterations: 14, Initial cost: 1.211734e+02, Final cost: 1.056751e+00, Termination: CONVERGENCE
Initial m: 0 c: 0
Final   m:0.291861 c: 0.131439

从参数值 $m = 0$ , $c = 0$ 开始，初始的目标函数值是121.173。Ceres找到了一组结果 $m = 0.291861, c = 0.131439$ 对应的目标函数值是1.05675。这些值与原始的模型参数 $y=e^{0.3x+0.1}$ 稍有不同，但是这个是可以预料到的。当从有噪声的数据中拟合曲线时，我们希望看到这个偏差（deviation）

Rubust Curve Fitting 鲁棒曲线拟合

如果假设数据中有一些离群值（outlier），也就是说有一些点不符合噪声模式，如果我们使用上面的代码来拟合这样的数据就会得到下面的曲线（摘自官网），注意拟合的曲线如何偏离真值。
在这里插入图片描述
为了处理离群值outlier，标准的方法是利用损失函数LossFunction, 损失函数减小了交大残差对残差块的影响，通常较大残差是由离群值影响的。为了使一个损失函数与残差块关联起来，我们将

problem.AddResidualBlock(Cost_function, nullptr, &m, &c)

替换为

problem.AddResidualBlock(const_function, new CauchyLoss(0.5), &m, &c);

CauchyLoss是Ceres Solver的的一个损失函数。参数0.5制定了损失函数的幅度。这样就可以得到下面的较好结果：略

Bundle Adjustment

光束平差法1
光束平差法2 推荐
 李群和李代数
 自然常数 $e$ ,工程中的自然数 $1$
编写Ceres的一个重要原因是我们需要解决大规模的光束平差问题。
给出一组测量的图像特征坐标和其他相应数据，光束平差法的目标是找到使重投影误差最小的3D点坐标和相机参数。这个优化问题通常被规划为非线性最小方差问题，其中误差是观测道德特征坐标与相应3D点在相机相平面投影之间的 $L_2$ 范数normal的平方square。Ceres 广泛的支持光束平差问题。

Reference
$L_0$ 范数：向量中非0元素的个数
$L_1$ 范数：向量中各个元素绝对值之和 $\|x\|_1=|x_1|+|x_2|+\dots+|x_n|$
$L_2$ 范数：向量各元素平方和后求平方根 $\|x_2\|_2=\sqrt{|x_1|^2+|x_2|^2+\dots+|x_n|^2}$

让我们来求解BAL数据集
此处未完待续

On Derivatives

Spivak Notation (斯皮瓦克标记法)

Michael Spivak – 微分几何的数学家
为了保持集体理智，我们使用Spivak关于导数的解释。这是一个实用的解释，使包含导数的表达式理解起来变得简单易懂
对于一个单变量的函数 $f$ , $f (a)$ 代表它再 $a$ 点的值， $D f$ 代表它的一次导数， $D f (a)$ 是在 $a$ 点的导数，也就是： $Df(a)=\frac{d}{dx}f(x)\Big|_{x=a}$
$D^k$ 代表 $f$ 的第 $k^{th}$ 次导数。
对于一个二元函数 $g (x, y)$ . $D_{1g}$ 与 $D_{2g}$ 分别代表 $g$ 的第一变量与第二变量的偏微分. 与经典标记法等效：
$D_{1g}=\frac{\partial}{\partial{x}}g(x,y)\\ D_{2g}=\frac{\partial}{\partial{y}}g(x,y)$
$D_g$ 代表 $g$ 的Jacobian矩阵：
$D_g=[D_{1g},D_{2g}]$
对于一个多变量函数 $g:\mathbb{R}^{n}\to\mathbb{R}^{m}$ , $D_g$ 代表 $m\times{n}$ 的Jacobian矩阵. $D_{ig}$ 是 $g$ 的第 $i$ 个坐标的偏微分，是 $D_g$ 的第 $i$ 列。
最后， $D_{1}^{2}g$ , $D_{1}D_{2}g$ 代表高阶偏微分。
For more see Michael Spivak’s book Calculus on Manifolds

Analytic Derivatives

考虑下面这个将曲线拟合到数据上的问题:
$\frac{b_{1}}{(1+e^{b_{2}-b_{3}x_{1}})^{1/b_{4}}}$
也就是说给定一些数据 $\{x_i,y_i\},\forall{i}=1,2,\dots,n$ , 确定最适合数据的参数 $b_1,b_2,b_3,b_4$ 。
这个问题也就是找到使下面目标函数最小的参数 $b_1,b_2,b_3,b_4$
$\begin{aligned}E(b_1,b_2,b_3,b_4)&=\sum_{i=1}^{n}{f^2(b_1,b_2,b_3,b_4;x_i,y_j)}\\ &=\sum_{i=1}^{n}\bigg(\frac{b_1}{(1+E^{b_2-b_3x_1})}-y_i\bigg)^2 \end{aligned}$
利用Ceres解决这个问题，我们需要定义一个代价函数来计算给定 $x$ 与 $y$ 关于 $b_1,b_2,b_3,b_4$ 的残差Residual $f$ 和导数.
使用初级微积分，我们可以得到：
$\begin{aligned} D_1f(b_1,b_2,b_3,b_4;x,y)&=\frac{1}{(1+e^{b_2-b_3x})^{1/b_4}}\\ D_2f(b_1,b_2,b_3,b_4;x,y)&=\frac{-b_1e^{b2-b_3x}}{b_4(1+e^{b_2-b_3x})^{1/b_4+1}}\\ D_3f(b_1,b_2,b_3,b_4;x,y)&=\frac{b_1xe^{b_2-b_3x}}{b_4(1+e^{b_2-b_3x})^{1/b_4+1}}\\ D_4f(b_1,b_2,b_3,b_4;x,y)&=\frac{b_1\log(1+e^{b_2-b_3x})}{b_4^2(1+e^{b_2-b_3x})^{1/b_4}} \end{aligned}$

Tat43AnalyticOptimized比Rat43Analytic快了2.8倍。这种运行时的差异十分普遍。为了获得解析导数的最佳性能，需要对通用子表达式进行优化。

什么时候需要使用解析微分？

1、表达式简
单，例如线性的
2、一个计算机代数软件像Maple, Mathematica, SymPy可以用来对目标函数计算符号微分，并生成C++代码。
3、对效率要求极高，可以利用一些代数结构来获得比自动微分更好的性能。
也就是说为了获得更好的性能，需要做大量的复杂工作。在走这条路之前，评估计算Jacobian矩阵在整个求解过程中所占的时间非常有必要。记住Amdahl’s Law对你很有帮助。
4、没有其他方法计算导数，例如要计算下面这个多项式Polyomial根关于 $x$ , $y$ 的导数：
$a_3(x,y)z^3+a_2(x,y)z^2+a_1(x,y)z+a_0(x,y)=0$
这个方程的导数求解需要使用Inverse Function Theorem
5、你就是喜欢手动计算导数。

附注Footnotes

最佳拟合的概念依赖用于评价拟合质量目标函数的选择，目标函数反过来也依赖于产生观测的基本噪声处理。当噪声是高斯噪声时，最小化方差和是正确的方法。在这种情况下，参数的优化值是极大似然估计

Numeric Derivatives

Forward Differences

Central Differences

Ridders’ Methord

Automatic Derivatives

Interfacing With Automatic Differentiation

Reference
当代价函数是一个明确的表达式时，可以直接使用自动微分。但是并不总是这样，有时代价函数需要与外部程序或数据进行交互。在这一章中我们考虑几种不同的方法。
我们考虑下面这个问题，通过查找参数 $\theta$ 与 $t$ 来优化下面这个优化问题：
$\begin{aligned} \min\quad&\sum{\big\|y_i-f(\|q_i\|^2)q_i\big\|^2}\\ \mathrm{such\ {}that}\quad&q_i=R(\theta)x_i+t \end{aligned}$
这里 $R$ 是一个二维的旋转矩阵，使用 $\theta$ 进行参数化。 $t$ 是一个二维向量。 $f$ 是一个外部的畸变函数distortion function
先考虑这个情况，首先有一个模板函数TemplatedComputeDistortion来计算函数 $f$ 。那么对应的残差仿函数就很简单，如下所示：

template <typename T> T TemplatedComputeDistortion(const T r2) {
  const double k1 = 0.0082;
  const double k2 = 0.000023;
  return 1.0 + k1 * y2 + k2 * r2 * r2;//此处y2 应为 r2
}

struct Affine2DWithDistortion {
  Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
    x[0] = x_in[0];
    x[1] = x_in[1];
    y[0] = y_in[0];
    y[1] = y_in[1];
  }

  template <typename T>
  bool operator()(const T* theta,
                  const T* t,
                  T* residuals) const {
    const T q_0 =  cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
    const T q_1 =  sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
    const T f = TemplatedComputeDistortion(q_0 * q_0 + q_1 * q_1); // !!!
    residuals[0] = y[0] - f * q_0;
    residuals[1] = y[1] - f * q_1;
    return true;
  }

  double x[2];
  double y[2];
};

现在我们考虑三种不能直接应用自动微分的 $f$ :

1.非模板函数
2.计算自身值与导数的函数
3.函数被定义成一个插值的函数表
下面依次考虑这三种情况.

A function that return its value(返回值的非模板函数)

假设给定一个方程ComputeDistortionValue声明signature如下

double ComputeDistortionValue(double r2)

这个函数计算了 $f$ 的值，函数内部的实现不重要，将这个函数与Affine2DWithDistortion对接总共分三步：

1、将ComputeDistortionValue封装warp成仿函数ComputeDistortionValueFunctor
2、对ComputeDistortionValueFunctor使用’NumericDiffCostFunction’进行数值微分
3、使用CostFunctionToFunctor对CostFunction对象进行封装，得到一个含有模板函数Operator()的仿函数对象，这个仿函数对象将NumericDiffCostFunction计算的Jacobian矩阵传送到Jet对象中

struct ComputeDistortionValueFunctor {
  bool operator()(const double* r2, double* value) const {
    *value = ComputeDistortionValue(r2[0]);
    return true;
  }
};

struct Affine2DWithDistortion {
  Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
    x[0] = x_in[0];
    x[1] = x_in[1];
    y[0] = y_in[0];
    y[1] = y_in[1];

    compute_distortion.reset(new ceres::CostFunctionToFunctor<1, 1>(
         new ceres::NumericDiffCostFunction<ComputeDistortionValueFunctor,
                                            ceres::CENTRAL,
                                            1,
                                            1>(
            new ComputeDistortionValueFunctor)));
  }

  template <typename T>
  bool operator()(const T* theta, const T* t, T* residuals) const {
    const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
    const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
    const T r2 = q_0 * q_0 + q_1 * q_1;
    T f;
    (*compute_distortion)(&r2, &f);
    residuals[0] = y[0] - f * q_0;
    residuals[1] = y[1] - f * q_1;
    return true;
  }

  double x[2];
  double y[2];
  std::unique_ptr<ceres::CostFunctionToFunctor<1, 1> > compute_distortion;
};

未完待续

guchaono1

关注

2
点赞
踩
19

收藏

觉得还不错? 一键收藏
1
评论
Ceres-Solver 官网教程翻译与学习

Ceres-Sovler仿函数欢迎使用Markdown编辑器新的改变功能快捷键合理的创建标题，有助于目录的生成如何改变文本的样式插入链接与图片如何插入一段漂亮的代码片生成一个适合你的列表创建一个表格设定内容居中、居左、居右SmartyPants创建一个自定义列表如何创建一个注脚注释也是必不可少的KaTeX数学公式新的甘特图功能，丰富你的文章UML 图表FLowchart流程图导出与导入导出导入仿函数仿函数参考ceres当中代价函数（cost function）的构建当中用到了仿函数（functor）
复制链接

扫一扫