Ceres入门笔记

介绍

Ceres可以解决下列形式的边界约束鲁棒非线性最小二乘问题

（1）

$minlimits_{x}quad frac{1}{2} sumlimits_{i} ho_{i}left( left| f_{i}left( x_{i1},ldots,x_{ik} ight) ight|^{2} ight)$

$s.t. quad l_{j} leqslant x_{j} leqslant u_{j}$

这种形式的问题广泛地出现在科学工程领域——从拟合统计曲线到计算机视觉中通过图像重建3D模型。

我们将通过Ceres求解器解决问题（1），本章给所有示例都提供了完整可行的代码。

表达式$ ho_{i}left( left| f_{i}left( x_{i1},ldots,x_{ik} ight) ight|^{2} ight)$作为ResidualBlock，其中$f_{i}left(cdot ight)$是取决于参数块$left[ x_{i1},ldots,x_{ik} ight]$的CostFunction。在大多数优化问题中，小群体的标量一起出现。比如，平移向量的三个组成部分和定义相机姿态四元数的四个分量。我们将这样一组小标量称为ParameterBlock。当然ParameterBlock也可以只有一个参数。$l_{j}$和$u_{j}$是参数块$x_{j}$的边界。

$ ho_{i}$是LossFunction,LossFunction是一个标量函数，用于减少异常值对非线性最小二乘解的影响。

一个特殊情况，当$ ho_{i}left( x ight)=x$，i.e.恒等函数,并且$l_{j}=-infty$和$u_{j}=infty$，我们得到了一更熟悉的非线性最小二乘问题。

（2）

$frac{1}{2} sumlimits_{i} left| f_{i}left( x_{i1},ldots,x_{ik} ight) ight|^{2}$

Hello World！

首先，考虑函数最小值的问题

$frac{1}{2}left(10-x ight)^{2}$

这是一个很简单的问题，其最小值为x=10，但是它非常适合用来解释如何通过Ceres解决问题。

第一步是编写一个函数来评估函数$fleft(x ight)=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() << "
";
  std::cout << "x : " << initial_x
            << " -> " << x << "
";
  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

 

一开始X =5，两次迭代中的求解器将得到了10 。仔细的读者会注意到这是一个线性问题，一个线性的解决方案应该足以获得最优值。而求解器的默认配置是针对非线性问题的，为了简单起见，我们在本例中没有改变它。 确实有可能在一次迭代中就可以使用Ceres来解决这个问题。 还要注意的时，求解器在第一次迭代中确实已经得到非常接近0的最优函数值。 当我们谈论Ceres的收敛和参数设置时，我们将更详细地讨论这些问题。

 

求导

像大多数用于优化的软件包一样，Ceres求解器依赖于能够在任意参数下评估目标函数中每个项的值和导数。 正确而有效地做到这一点对于取得好的结果至关重要。 Ceres Solver提供了许多方法。 

我们现在考虑其他两种可能性。 也就是解析求导和数值求导。

数值求导

在某些情况下，无法定义模板损失函数，例如，当残差评估涉及对您无法控制的库函数的调用时。 在这种情况下，可以使用数值区分。 用户定义了一个计算残差的函数，并用它构造一个NumericDiffCostFunction。 例如，$fleft(x ight)=10-x$对应的函数是

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的子类。 这里例如是实现$fleft(x ight)=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;
  }
};

SimpleCostFunction :: Evaluate提供输入parameters数组，residuals输出数组残差和jacobians输出数组雅可比。 jacobians数组是可选的，Evaluate需要检查它是否为非空值，如果是这种情况，则用残差函数的导数值填充它。 在这种情况下，由于残差函数是线性的，因此雅可比行列式是恒定的。

从上面的代码片段可以看出，实现CostFunction对象有点乏味。 我们建议，除非您有充分理由来自己管理雅可比计算，否则使用AutoDiffCostFunction或NumericDiffCostFunction来构建残差块。解析求导程序在examples/helloworld_analytic_diff.cc。

#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() << "
";
  std::cout << "x : " << initial_x
            << " -> " << x << "
";

  return 0;
}