2.Ceres官方教程-非线性最小二乘～Derivatives(导数)

像大多数优化软件包一样，Ceres求解器依赖其能够在任意参数值下评估目标函数中每一项的值和导数。 正确而高效地做到这一点是取得好结果的关键。Ceres提供了一系列解决方案，其中一个就是在Hello World中用到的Automatic Differentiation (自动微分算法)。我们将探讨另外两种可能性，即解析法（Analytic）和数值法（Numeric ）求导。

1.Numeric Derivatives(数值法求导)

在某些情况下，只定义一个仿函数是不行的，例如，当残差的求值涉及调用您无法控制的库函数时。在这种情况下，可以使用数值微分法。
用户定义一个仿函数(functor)来计算残差值，并且构建一个数值微分代价函数(NumericDiffCostFunction)。比如对于 f(x)=10−x 对应函数体如下：

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

/*
对比Hello world中的Functor的定义
可以发现，这次没用模板类。
*/

然后继续添加Problem

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

对比在Hello World中使用automatic算法

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

可以发现两种算法在构建Problem时候基本差不多。但是在用Nummeric算法时需要额外给定一个参数ceres::CENTRAL 。这个参数告诉计算机如何计算导数。
更多具体介绍可以参看NumericDiffCostFunction的文档(http://ceres-solver.org/numerical_derivatives.html?highlight=numericdiffcostfunction)

一般来说，Ceres官方更加推荐自动微分算法，因为C++模板类使自动算法有更高的效率。数值微分算法通常来说计算更复杂，收敛更缓慢。

代码在ceres-solver-1.14.0/examples/helloworld_numeric_diff.cc

#include "ceres/ceres.h"
#include "glog/logging.h"

using ceres::NumericDiffCostFunction;
using ceres::CENTRAL;
using ceres::CostFunction;
using ceres::Problem;
using ceres::Solver;
using ceres::Solve;

// A cost functor that implements the residual r = 10 - x.
struct CostFunctor {
  bool operator()(const double* const x, double* residual) const {
    residual[0] = 10.0 - x[0];
    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). This uses
  // numeric differentiation to obtain the derivative (jacobian).
  CostFunction* cost_function =
      new NumericDiffCostFunction<CostFunctor, CENTRAL, 1, 1> (new CostFunctor);
  problem.AddResidualBlock(cost_function, NULL, &x);

  // Run the solver!
  Solver::Options options;
  options.minimizer_progress_to_stdout = true;
  options.linear_solver_type = ceres::DENSE_QR;
  Solver::Summary summary;
  Solve(options, &problem, &summary);

  std::cout << summary.BriefReport() << "\n";
  std::cout << "initial x : " << initial_x << ", "
            << "result x : " << x << "\n";
  return 0;
}

结果

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    1.44e-04    3.64e-04
   1  4.511690e-07    4.51e+01    9.50e-04   9.50e+00   1.00e+00  3.00e+04        1    9.17e-05    5.92e-04
   2  5.012669e-16    4.51e-07    3.17e-08   9.50e-04   1.00e+00  9.00e+04        1    8.76e-06    6.08e-04
trust_region_minimizer.cc:687 Terminating: Parameter tolerance reached. Relative step_norm: 3.166246e-09 <= 1.000000e-08.
Ceres Solver Report: Iterations: 3, Initial cost: 4.512500e+01, Final cost: 5.012669e-16, Termination: CONVERGENCE
initial x : 0.5, result x : 10

2.Analytic Derivatives(解析法求导)

有些时候，应用自动求解算法时不可能的。比如在某些情况下，计算导数的时候，使用闭合解（closed form，也被称为解析解）会比使用自动微分算法中的链式法则（chain rule）更有效率。

解析解，数值解参考: https://blog.csdn.net/qianmianyuan/article/details/8998281
例如: 求解3÷7。解析解为3/7，数值解为0.429；

链式法则: 微积分中求导法则，用于求复合函数的导数
参考: https://blog.csdn.net/weixin_44010756/article/details/109662171

在这种情况下，你就可以自己编写残差和雅可比行列式的计算代码了。为此，定义一个CostFunction的子类，或者在编译时知道了参数的大小和残差，则可以定义SizedCostFunction的子类。

例如，为了实现f(x) = 10 - x，有一个简单的例子

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;
  }
};

parameters：输入数组
residuals/jacobians：输出数组
Evaluate函数用于检查jacobians是否为非空值。如果不为空，那么就把残差方程的导数值赋值给它。在这种情况下，由于残差函数是线性的，雅可比矩阵是常数。

从上面的代码片段可以看出，实现CostFunction对象有点繁琐。我们建议，除非有特殊的原因自己管理雅可比矩阵计算，否则您可以使用AutoDiffCostFunction或NumericDiffCostFunction来构造残差块。

代码在ceres-solver-1.14.0/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);

  // Run the solver!
  Solver::Options options;
  options.minimizer_progress_to_stdout = true;
  options.linear_solver_type = ceres::DENSE_QR;
  Solver::Summary summary;
  Solve(options, &problem, &summary);

  std::cout << summary.BriefReport() << "\n";
  std::cout << "initial x : " << initial_x << ", "
            << "result x : " << x << "\n";

  return 0;
}

结果

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    2.16e-05    7.73e-05
   1  4.511598e-07    4.51e+01    9.50e-04   9.50e+00   1.00e+00  3.00e+04        1    3.16e-05    1.39e-04
   2  5.012552e-16    4.51e-07    3.17e-08   9.50e-04   1.00e+00  9.00e+04        1    2.13e-05    1.88e-04
trust_region_minimizer.cc:687 Terminating: Parameter tolerance reached. Relative step_norm: 3.166209e-09 <= 1.000000e-08.
Ceres Solver Report: Iterations: 3, Initial cost: 4.512500e+01, Final cost: 5.012552e-16, Termination: CONVERGENCE
initial x : 0.5, result x : 10

3.More About Derivatives(其他导数计算方法)

到目前为止，计算导数其实是Ceres最复杂的部分了。根据需要，用户有时候还需要更复杂的导数计算算法。这一节仅仅是大体介绍了如何使用Ceres进行导数计算最浅显的部分。对NumericDiffCostFunction和AutoDiffCostFunction方法都很熟悉之后，还可以深入研究一下DynamicAutoDiffCostFunction, CostFunctionToFunctor, NumericDiffFunctor和ConditionedCostFunction，从而实现更高级的代价函数的计算方法。