Ceres入门笔记
介绍
Ceres可以解决下列形式的边界约束鲁棒非线性最小二乘问题
（1）
$\min\limits_{x}\quad \frac{1}{2} \sum\limits_{i}\rho_{i}\left( \left\| f_{i}\left( x_{i1},\ldots,x_{ik}\right)\right\|^{2} \right)$
$s.t. \quad l_{j} \leqslant x_{j} \leqslant u_{j}$
这种形式的问题广泛地出现在科学工程领域——从拟合统计曲线到计算机视觉中通过图像重建3D模型。
我们将通过Ceres求解器解决问题（1），本章给所有示例都提供了完整可行的代码。
表达式$\rho_{i}\left( \left\| f_{i}\left( x_{i1},\ldots,x_{ik}\right)\right\|^{2} \right)$作为ResidualBlock，其中$f_{i}\left(\cdot\right)$是取决于参数块$\left[ x_{i1},\ldots,x_{ik} \right]$的CostFunction。在大多数优化问题中，小群体的标量一起出现。比如，平移向量的三个组成部分和定义相机姿态四元数的四个分量。我们将这样一组小标量称为ParameterBlock。当然ParameterBlock也可以只有一个参数。$l_{j}$和$u_{j}$是参数块$x_{j}$的边界。
$\rho_{i}$是LossFunction,LossFunction是一个标量函数，用于减少异常值对非线性最小二乘解的影响。
一个特殊情况，当$\rho_{i}\left( x \right)=x$，i.e.恒等函数,并且$l_{j}=-\infty$和$u_{j}=\infty$，我们得到了一更熟悉的非线性最小二乘问题。
（2）
$\frac{1}{2} \sum\limits_{i} \left\| f_{i}\left( x_{i1},\ldots,x_{ik}\right)\right\|^{2}$
Hello World！
首先，考虑函数最小值的问题
$\frac{1}{2}\left(10-x\right)^{2}$
这是一个很简单的问题，其最小值为x=10，但是它非常适合用来解释如何通过Ceres解决问题。
第一步是编写一个函数来评估函数$f\left(x\right)=10-x$；
struct CostFunctor { template <typename T> bool operator()(const T* const x, T* residual) const { residual[0] = T(10.0) - x[0]; return true; } };
这里需要注意的一点是operator()是一个模板方法，它假定所有的输入和输出都是T类型的。此处使用模板允许调用CostFunction::operator<T>()。当只有残差被用到时T=double，当用到雅可比时T=Jet。
一旦我们有了计算残差函数的方法，就可以用它构造一个非线性最小二乘问题，并让Ceres解决它。
int main(int argc, char** argv) { google::InitGoogleLogging(argv[0]); // The variable to solve for with its initial value. double initial_x = 5.0; double x = initial_x; // Build the problem. Problem problem; // Set up the only cost function (also known as residual). This uses // auto-differentiation to obtain the derivative (jacobian). CostFunction* cost_function = new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);//CostFunctor结构 problem.AddResidualBlock(cost_function, NULL, &x); // Run the solver! Solver::Options options; options.linear_solver_type = ceres::DENSE_QR; options.minimizer_progress_to_stdout = true; Solver::Summary summary; Solve(options, &problem, &summary); std::cout << summary.BriefReport() << "\n"; std::cout << "x : " << initial_x << " -> " << x << "\n"; return 0; }
AutoDiffCostFunction将CostFunction作为输入，并且自动求导，并为其提供CostFunction接口。
编译运行examples/helloworld.cc（在下载的Ceres的examples文件夹下）有以下结果
iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time 0 4.512500e+01 0.00e+00 9.50e+00 0.00e+00 0.00e+00 1.00e+04 0 5.33e-04 3.46e-03 1 4.511598e-07 4.51e+01 9.50e-04 9.50e+00 1.00e+00 3.00e+04 1 5.00e-04 4.05e-03 2 5.012552e-16 4.51e-07 3.17e-08 9.50e-04 1.00e+00 9.00e+04 1 1.60e-05 4.09e-03 Ceres Solver Report: Iterations: 2, Initial cost: 4.512500e+01, Final cost: 5.012552e-16, Termination: CONVERGENCE x : 5.0 -> 10
struct NumericDiffCostFunctor { bool operator()(const double* const x, double* residual) const { residual[0] = 10.0 - x[0]; return true; } };
添加到problem中：
CostFunction* cost_function = new NumericDiffCostFunction<NumericDiffCostFunctor, ceres::CENTRAL, 1, 1>(new NumericDiffCostFunctor); problem.AddResidualBlock(cost_function, NULL, &x);
注意，自动求导时我们用的是
CostFunction* cost_function =
new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
problem.AddResidualBlock(cost_function, NULL, &x);
除了一个额外的模板参数，该参数表明用于计算数值导数（该示例在examples/helloworld_numeric_diff.cc）的有限差分格式的种类之外，其结构看起来与用于自动求导的格式几乎完全相同。
一般来说，我们建议自动求导而不是数值求导。 C ++模板的使用使得自动求导变得高效，而数值求导很expensive，容易出现数值错误，并导致收敛速度变慢。
解析求导
在某些情况下，使用自动区分是不可能的。 例如，以封闭形式计算导数比依赖自动求导代码所使用的链式规则更有效。
在这种情况下，可以提供自己写的残差和雅可比计算代码。 为此，如果您知道编译时参数和残差的大小，请定义CostFunction或SizedCostFunction的子类。 这里例如是实现$f\left(x\right)=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 != NULL && jacobians[0] != NULL) { jacobians[0][0] = -1;//残差的导数为-1，残差函数是线性的。 } return true; } };
#include <vector> #include "ceres/ceres.h" #include "glog/logging.h" using ceres::CostFunction; using ceres::SizedCostFunction; using ceres::Problem; using ceres::Solver; using ceres::Solve; // A CostFunction implementing analytically derivatives for the // function f(x) = 10 - x. class QuadraticCostFunction : public SizedCostFunction<1 /* number of residuals */, 1 /* size of first parameter */> { public: virtual ~QuadraticCostFunction() {} virtual bool Evaluate(double const* const* parameters, double* residuals, double** jacobians) const { double x = parameters[0][0]; // f(x) = 10 - x. residuals[0] = 10 - x; // f'(x) = -1. Since there's only 1 parameter and that parameter // has 1 dimension, there is only 1 element to fill in the // jacobians. // // Since the Evaluate function can be called with the jacobians // pointer equal to NULL, the Evaluate function must check to see // if jacobians need to be computed. // // For this simple problem it is overkill to check if jacobians[0] // is NULL, but in general when writing more complex // CostFunctions, it is possible that Ceres may only demand the // derivatives w.r.t. a subset of the parameter blocks. if (jacobians != NULL && jacobians[0] != NULL) { jacobians[0][0] = -1; } return true; } }; int main(int argc, char** argv) { google::InitGoogleLogging(argv[0]); // The variable to solve for with its initial value. It will be // mutated in place by the solver. double x = 0.5; const double initial_x = x; // Build the problem. Problem problem; // Set up the only cost function (also known as residual). CostFunction* cost_function = new QuadraticCostFunction; problem.AddResidualBlock(cost_function, NULL, &x); //x待估计参数 // Run the solver! Solver::Options options; options.minimizer_progress_to_stdout = true; Solver::Summary summary; Solve(options, &problem, &summary); std::cout << summary.BriefReport() << "\n"; std::cout << "x : " << initial_x << " -> " << x << "\n"; return 0; }
