junit单元测试 mybatis 单独测试 springboot-CFANZ编程社区

现在我们将讨论自动微分算法。它是一种可以快速计算精确导数的算法，同时用户只要做与数值微分法类似的工作。下面的代码片段实现了对Rat43的CostFunction。

struct Rat43CostFunctor {
  Rat43CostFunctor(const double x, const double y) : x_(x), y_(y) {}

  template <typename T>
  bool operator()(const T* parameters, T* residuals) const {
    const T b1 = parameters[0];
    const T b2 = parameters[1];
    const T b3 = parameters[2];
    const T b4 = parameters[3];
    residuals[0] = b1 * pow(1.0 + exp(b2 -  b3 * x_), -1.0 / b4) - y_;
    return true;
  }

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


CostFunction* cost_function =
      new AutoDiffCostFunction<Rat43CostFunctor, 1, 4>(
        new Rat43CostFunctor(x, y));

注意，与数值微分法相比，在定义自动微分的Functor时，唯一的区别是对操作符operator()的设置。
PS：数值微分的代码

			struct Rat43CostFunctor {
			  Rat43CostFunctor(const double x, const double y) : x_(x), y_(y) {}
			
			  bool operator()(const double* parameters, double* residuals) const {
			    const double b1 = parameters[0];
			    const double b2 = parameters[1];
			    const double b3 = parameters[2];
			    const double b4 = parameters[3];
			    residuals[0] = b1 * pow(1.0 + exp(b2 -  b3 * x_), -1.0 / b4) - y_;
			    return true;
			  }
			
			  const double x_;
			  const double y_;
			}
			
			CostFunction* cost_function =
			  new NumericDiffCostFunction<Rat43CostFunctor, FORWARD, 1, 4>(
			    new Rat43CostFunctor(x, y));

在数值微分的情况下，它是

bool operator()(const double* parameters, double* residuals) const;

对于自动微分，它是一个模板化的函数格式

template <typename T> bool operator()(const T* parameters, T* residuals) const;

这个变化有什么影响呢？下表比较了采用不同方法计算Rat43的残差和雅可比矩阵所需的时间。

CostFunction	Time (ns)
Rat43Analytic	255
Rat43AnalyticOptimized	92
Rat43NumericDiffForward	262
Rat43NumericDiffCentral	517
Rat43NumericDiffRidders	3760
Rat43AutomaticDiff	129

我们可以使用自动微分(Rat43AutomaticDiff)得到精确的微分，这与编写数值微分代码的工作量差不多，但只比手动优化的解析微分慢40%。

那么它是如何工作的呢?为此，我们将学习 Dual Numbers 和 Jets 。

Dual Numbers对偶数 & Jets

对偶数是实数的延伸，类似于复数:而复数通过引入一个虚单位 $i^2=-1$ 来扩充实数，对偶数引入一个无穷小单位的 $ϵ$ ，从而 $ϵ^2=0$ 。
对偶数 $a + v ϵ$ 有两个分量，实分量a和无穷小分量v。
PS: 对偶数参考 https://zhuanlan.zhihu.com/p/380140763

令人惊讶的是，这个简单的改变导致了一个方便的方法来计算精确的导数，而不需要操作复杂的符号表达式。例如，考虑下面的函数
$f(x) = x^2 ,$
然后，
$\begin{split} f(10 + \epsilon) &= (10 + \epsilon)^2\\ &= 100 + 20 \epsilon + \epsilon^2\\ &= 100 + 20 \epsilon \end{split}$

观察ϵ的系数，我们发现Df(10)=20。事实上，这可以推广到非多项式的函数。考虑一个任意可微函数f(x)
然后我们可以通过考虑f在x附近的泰勒展开式来计算 $\epsilon)$ ，它给出了无穷级数
$\begin{split} f(x + \epsilon) &= f(x) + Df(x) \epsilon + D^2f(x) \frac{\epsilon^2}{2} + D^3f(x) \frac{\epsilon^3}{6} + \cdots\\ f(x + \epsilon) &= f(x) + Df(x) \epsilon \end{split}$

注意， $\epsilon^2 = 0$ 。

Jet是一个n维dual number。其中我们用n个无穷小单位 $\epsilon_i,\ i=1,...,n$ 来增加实数，并且存在性质 $\forall i, j\ :\epsilon_i\epsilon_j = 0$ 。那么Jet由实部a和n维无穷小部分v组成，即，
$\sum_j v_{j} \epsilon_j$

求和符号很乏味，所以我们还是写下来
$\mathbf{v}.$

上面这个式子， $\epsilon_i$ 是隐含的。然后，使用与上面相同的泰勒级数展开，我们可以看到:
$\mathbf{v}) = f(a) + Df(a) \mathbf{v}.$
类似地，对于一个多元函数 $f:\mathbb{R}^{n}\rightarrow \mathbb{R}^m$ ，对于 $x_i = a_i + \mathbf{v}_i,\ \forall i = 1,...,n$

对于每一个 $\mathbf{v}_i = e_i$ 是第i个标准基向量。然后，将上面的表达式简化为
$f(x_1,..., x_n) = f(a_1, ..., a_n) + \sum_i D_i f(a_1, ..., a_n) \epsilon_i$

我们可以通过检查的系数来提取雅可比矩阵的坐标 $\epsilon_i$ 。

Implementing Jets

为了在实践中发挥作用，我们需要有能力计算任意函数f，不仅在实数上，而且在对偶数上，但我们通常不通过泰勒展开来计算函数，

这就是c++模板和操作符重载发挥作用的地方。下面的代码片段有一个Jet的简单实现和一些操作它们的操作符/函数。

template<int N> struct Jet {
  double a;
  Eigen::Matrix<double, 1, N> v;
};

template<int N> Jet<N> operator+(const Jet<N>& f, const Jet<N>& g) {
  return Jet<N>(f.a + g.a, f.v + g.v);
}

template<int N> Jet<N> operator-(const Jet<N>& f, const Jet<N>& g) {
  return Jet<N>(f.a - g.a, f.v - g.v);
}

template<int N> Jet<N> operator*(const Jet<N>& f, const Jet<N>& g) {
  return Jet<N>(f.a * g.a, f.a * g.v + f.v * g.a);
}

template<int N> Jet<N> operator/(const Jet<N>& f, const Jet<N>& g) {
  return Jet<N>(f.a / g.a, f.v / g.a - f.a * g.v / (g.a * g.a));
}

template <int N> Jet<N> exp(const Jet<N>& f) {
  return Jet<T, N>(exp(f.a), exp(f.a) * f.v);
}

// This is a simple implementation for illustration purposes, the
// actual implementation of pow requires careful handling of a number
// of corner cases.
template <int N>  Jet<N> pow(const Jet<N>& f, const Jet<N>& g) {
  return Jet<N>(pow(f.a, g.a),
                g.a * pow(f.a, g.a - 1.0) * f.v +
                pow(f.a, g.a) * log(f.a); * g.v);
}

有了这些重载函数，我们现在可以用一个Jets数组而不是double双精度调用Rat43CostFunctor。将其与适当初始化的Jets组合在一起，我们可以计算雅可比矩阵，如下所示:

class Rat43Automatic : public ceres::SizedCostFunction<1,4> {
 public:
  Rat43Automatic(const Rat43CostFunctor* functor) : functor_(functor) {}
  virtual ~Rat43Automatic() {}
  virtual bool Evaluate(double const* const* parameters,
                        double* residuals,
                        double** jacobians) const {
    // Just evaluate the residuals if Jacobians are not required.
    if (!jacobians) return (*functor_)(parameters[0], residuals);

    // Initialize the Jets
    ceres::Jet<4> jets[4];
    for (int i = 0; i < 4; ++i) {
      jets[i].a = parameters[0][i];
      jets[i].v.setZero();
      jets[i].v[i] = 1.0;
    }

    ceres::Jet<4> result;
    (*functor_)(jets, &result);

    // Copy the values out of the Jet.
    residuals[0] = result.a;
    for (int i = 0; i < 4; ++i) {
      jacobians[0][i] = result.v[i];
    }
    return true;
  }

 private:
  std::unique_ptr<const Rat43CostFunctor> functor_;
};

实际上，这就是AutoDiffCostFunction的工作原理。

陷阱

自动微分将用户从计算和推理雅可比矩阵符号表达式的负担中解放出来，但是这种自由是有代价的。例如，考虑下面的简单仿函数(functor):

struct Functor {
  template <typename T> bool operator()(const T* x, T* residual) const {
    residual[0] = 1.0 - sqrt(x[0] * x[0] + x[1] * x[1]);
    return true;
  }
};

查看代码的残差计算，没有预见任何问题。但是，如果我们看一下雅可比矩阵的解析表达式
$\begin{split} y &= 1 - \sqrt{x_0^2 + x_1^2}\\ D_1y &= -\frac{x_0}{\sqrt{x_0^2 + x_1^2}},\ D_2y = -\frac{x_1}{\sqrt{x_0^2 + x_1^2}}\end{split}$