# Analytic Derivatives

## Example of Analytic Derivatives

考虑通过样本数据对如下问题 ([Rat43](http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml)) 进行曲线拟合

$$
y=\frac{b\_{1}}{\left(1+e^{b\_{2}-b\_{3} x}\right)^{1 / b\_{4}}}
$$

给定一些样本数据 $$\left{x\_{i}, y\_{i}\right}, \forall i=1, \ldots, n$$ ，估计参数 $$b\_{1}, b\_{2}, b\_{3},b\_{4}$$ 。可以建模成寻找使得如下目标函数取得最小值时， $$b\_{1}, b\_{2}, b\_{3},b\_{4}$$ 的值：

$$
\begin{aligned} E\left(b\_{1}, b\_{2}, b\_{3}, b\_{4}\right) & =\sum\_{i} f^{2}\left(b\_{1}, b\_{2}, b\_{3}, b\_{4} ; x\_{i}, y\_{i}\right) \ & =\sum\_{i}\left(\frac{b\_{1}}{\left(1+e^{b\_{2}-b\_{3} x\_{i}}\right)^{1 / b\_{4}}}-y\_{i}\right)^{2} \end{aligned}
$$

使用 Ceres 解决该问题，我们需要定义一个 costFunction，用于计算给定 $$x$$ 和 $$y$$ 的残差 $$f$$ 及其相对于 $$b\_{1}、b\_{2}、b\_{3}$$ 和 $$b\_{4}$$ 的导数。

$$
\begin{array}{l} D\_{1} f\left(b\_{1}, b\_{2}, b\_{3}, b\_{4} ; x, y\right)=\frac{1}{\left(1+e^{b\_{2}-b\_{3} x}\right)^{1 / b\_{4}}} \ D\_{2} f\left(b\_{1}, b\_{2}, b\_{3}, b\_{4} ; x, y\right)=\frac{-b\_{1} e^{b\_{2}-b\_{3} x}}{b\_{4}\left(1+e^{b\_{2}-b\_{3} x}\right)^{1 / b\_{4}+1}} \ D\_{3} f\left(b\_{1}, b\_{2}, b\_{3}, b\_{4} ; x, y\right)=\frac{b\_{1} x e^{b\_{2}-b\_{3} x}}{b\_{4}\left(1+e^{b\_{2}-b\_{3} x}\right)^{1 / b\_{4}+1}} \ D\_{4} f\left(b\_{1}, b\_{2}, b\_{3}, b\_{4} ; x, y\right)=\frac{b\_{1} \log \left(1+e^{b\_{2}-b\_{3} x}\right)}{b\_{4}^{2}\left(1+e^{b\_{2}-b\_{3} x}\right)^{1 / b\_{4}}} \end{array}
$$

有了这些导数，我们现在可以将 CostFunction 实现为：

```cpp
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_;
 };
```

这些难以阅读并且有很多冗余。因此在实践中我们会用一些子表达式来提高其效率，这会给我们带来类似的结果：

```cpp
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_;
 };
```

两种实现的区别如下：

|  CostFunction | Times(ns) |
| :-----------: | :-------: |
| Rat43Analytic |    255    |
| Rat43Analytic |     92    |

## When should you use analytical derivatives?

* 表达式很简单，例如大多是线性的。
* 可以使用 [Maple](https://www.maplesoft.com/products/maple/) 、[Mathematica](https://www.wolfram.com/mathematica/) 或 [SymPy](http://www.sympy.org/en/index.html) 等计算机代数系统对目标函数进行微分，并生成 C++ 对其进行评估。
* 性能是最先考虑的，可以利用计算中的代数结构来获得比自动微分更好的性能。尽管如此，要想从 analytical derivatives 中获得最佳性能，还需要做大量的工作。在采用这种方法之前，最好先估算一下求解 jacobian 所花费的时间占总求解时间的比例，使用 [Amdahl’s Law](https://en.wikipedia.org/wiki/Amdahl's_law) 即可。
* 您喜欢链式法则，并且喜欢手工计算所有导数。
* There is no other way to compute the derivatives, e.g. you wish to compute the derivative of the root of a polynomial, $$a\_{3}(x, y) z^{3}+a\_{2}(x, y) z^{2}+a\_{1}(x, y) z+a\_{0}(x, y)=0$$, with respect to $$x,y$$, This requires the use of the [Inverse Function Theorem](https://en.wikipedia.org/wiki/Inverse_function_theorem).