9.Ceres官方教程-On Derivatives~Analytic Derivatives

考虑以下曲线(Rat43 https://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml)的拟合问题:

也就是说，给定一些数据，确定最适合该数据的参数
我们面临的问题就是求解使下列表达式的取值最小：

最佳拟合的概念取决于用来衡量拟合质量的目标函数的选择，而该目标函数又取决于产生观测结果的潜在噪声过程。当噪声是高斯噪声时，将差的平方和最小化是正确的做法。在这种情况下，参数的最优值是最大似然估计

为了使用Ceres求解器来解决这个问题，我们需要定义一个CostFunction来计算给定x和y的残差f及其对b1,b2,b3和b4的导数。

根据高等数学中微积分的知识，我们可以算出f的一系列导数：

根据这些手动计算出的导数，我们现在可以实现CostFunction了:

class Rat43Analytic : public SizedCostFunction<1,4> {
   public:
     Rat43Analytic(const double x, const double y) : x_(x), y_(y) {}
     virtual ~Rat43Analytic() {}
     virtual bool Evaluate(double const* const* parameters,
                           double* residuals,
                           double** jacobians) const {
       const double b1 = parameters[0][0];
       const double b2 = parameters[0][1];
       const double b3 = parameters[0][2];
       const double b4 = parameters[0][3];

       residuals[0] = b1 *  pow(1 + exp(b2 -  b3 * x_), -1.0 / b4) - y_;

       if (!jacobians) return true;
       double* jacobian = jacobians[0];
       if (!jacobian) return true;

       jacobian[0] = pow(1 + exp(b2 - b3 * x_), -1.0 / b4);
       jacobian[1] = -b1 * exp(b2 - b3 * x_) *
                     pow(1 + exp(b2 - b3 * x_), -1.0 / b4 - 1) / b4;
       jacobian[2] = x_ * b1 * exp(b2 - b3 * x_) *
                     pow(1 + exp(b2 - b3 * x_), -1.0 / b4 - 1) / b4;
       jacobian[3] = b1 * log(1 + exp(b2 - b3 * x_)) *
                     pow(1 + exp(b2 - b3 * x_), -1.0 / b4) / (b4 * b4);
       return true;
     }

    private:
     const double x_;
     const double y_;
 };

这是一段冗长乏味的代码，难以阅读，并且有很多冗余。所以在实践中，我们会缓存一些子表达式来提高它的效率，这将会得到如下结果:

class Rat43AnalyticOptimized : public SizedCostFunction<1,4> {
   public:
     Rat43AnalyticOptimized(const double x, const double y) : x_(x), y_(y) {}
     virtual ~Rat43AnalyticOptimized() {}
     virtual bool Evaluate(double const* const* parameters,
                           double* residuals,
                           double** jacobians) const {
       const double b1 = parameters[0][0];
       const double b2 = parameters[0][1];
       const double b3 = parameters[0][2];
       const double b4 = parameters[0][3];

       const double t1 = exp(b2 -  b3 * x_);
       const double t2 = 1 + t1;
       const double t3 = pow(t2, -1.0 / b4);
       residuals[0] = b1 * t3 - y_;

       if (!jacobians) return true;
       double* jacobian = jacobians[0];
       if (!jacobian) return true;

       const double t4 = pow(t2, -1.0 / b4 - 1);
       jacobian[0] = t3;
       jacobian[1] = -b1 * t1 * t4 / b4;
       jacobian[2] = -x_ * jacobian[1];
       jacobian[3] = b1 * log(t2) * t3 / (b4 * b4);
       return true;
     }

   private:
     const double x_;
     const double y_;
 };

这两种实现在性能上有什么不同?

Rat43AnalyticOptimized比Rat43Analytic快2.8倍。这种运行时上的差异并不少见。为了从分析计算的导数中获得最佳性能，通常需要优化代码以考虑公共子表达式。

什么时候应该使用analytical derivatives?

1.表达式很简单，例如大部分是线性的
2.像Maple、Mathematica或symy这样的计算机代数系统可以用来符号化地区分目标函数，并生成c++来计算它们。
3.式子中有一些代数结构可以实现比自动微分有更好的性能。
也就是说, 获得在计算倒数之外的最大性能需要大量的工作.在沿着这条路径走下去之前，评估雅可比矩阵的计算花费是整个求解时间的一小部分是很有用的，，记住Amdahl法则是你的朋友。
4.没有其他方法来计算导数，例如你想要计算一个多项式的根的导数:

对于x，y，这就需要用到反函数定理

5.你喜欢链式法则也喜欢手动计算代数运算。