Squashed 'third_party/ceres/' content from commit e51e9b4
Change-Id: I763587619d57e594d3fa158dc3a7fe0b89a1743b
git-subtree-dir: third_party/ceres
git-subtree-split: e51e9b46f6ca88ab8b2266d0e362771db6d98067
diff --git a/docs/source/automatic_derivatives.rst b/docs/source/automatic_derivatives.rst
new file mode 100644
index 0000000..0c48c80
--- /dev/null
+++ b/docs/source/automatic_derivatives.rst
@@ -0,0 +1,307 @@
+.. default-domain:: cpp
+
+.. cpp:namespace:: ceres
+
+.. _chapter-automatic_derivatives:
+
+=====================
+Automatic Derivatives
+=====================
+
+We will now consider automatic differentiation. It is a technique that
+can compute exact derivatives, fast, while requiring about the same
+effort from the user as is needed to use numerical differentiation.
+
+Don't believe me? Well here goes. The following code fragment
+implements an automatically differentiated ``CostFunction`` for `Rat43
+<http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml>`_.
+
+.. code-block:: c++
+
+ 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));
+
+Notice that compared to numeric differentiation, the only difference
+when defining the functor for use with automatic differentiation is
+the signature of the ``operator()``.
+
+In the case of numeric differentiation it was
+
+.. code-block:: c++
+
+ bool operator()(const double* parameters, double* residuals) const;
+
+and for automatic differentiation it is a templated function of the
+form
+
+.. code-block:: c++
+
+ template <typename T> bool operator()(const T* parameters, T* residuals) const;
+
+
+So what does this small change buy us? The following table compares
+the time it takes to evaluate the residual and the Jacobian for
+`Rat43` using various methods.
+
+========================== =========
+CostFunction Time (ns)
+========================== =========
+Rat43Analytic 255
+Rat43AnalyticOptimized 92
+Rat43NumericDiffForward 262
+Rat43NumericDiffCentral 517
+Rat43NumericDiffRidders 3760
+Rat43AutomaticDiff 129
+========================== =========
+
+We can get exact derivatives using automatic differentiation
+(``Rat43AutomaticDiff``) with about the same effort that is required
+to write the code for numeric differentiation but only :math:`40\%`
+slower than hand optimized analytical derivatives.
+
+So how does it work? For this we will have to learn about **Dual
+Numbers** and **Jets** .
+
+
+Dual Numbers & Jets
+===================
+
+.. NOTE::
+
+ Reading this and the next section on implementing Jets is not
+ necessary to use automatic differentiation in Ceres Solver. But
+ knowing the basics of how Jets work is useful when debugging and
+ reasoning about the performance of automatic differentiation.
+
+Dual numbers are an extension of the real numbers analogous to complex
+numbers: whereas complex numbers augment the reals by introducing an
+imaginary unit :math:`\iota` such that :math:`\iota^2 = -1`, dual
+numbers introduce an *infinitesimal* unit :math:`\epsilon` such that
+:math:`\epsilon^2 = 0` . A dual number :math:`a + v\epsilon` has two
+components, the *real* component :math:`a` and the *infinitesimal*
+component :math:`v`.
+
+Surprisingly, this simple change leads to a convenient method for
+computing exact derivatives without needing to manipulate complicated
+symbolic expressions.
+
+For example, consider the function
+
+.. math::
+
+ f(x) = x^2 ,
+
+Then,
+
+.. math::
+
+ \begin{align}
+ f(10 + \epsilon) &= (10 + \epsilon)^2\\
+ &= 100 + 20 \epsilon + \epsilon^2\\
+ &= 100 + 20 \epsilon
+ \end{align}
+
+Observe that the coefficient of :math:`\epsilon` is :math:`Df(10) =
+20`. Indeed this generalizes to functions which are not
+polynomial. Consider an arbitrary differentiable function
+:math:`f(x)`. Then we can evaluate :math:`f(x + \epsilon)` by
+considering the Taylor expansion of :math:`f` near :math:`x`, which
+gives us the infinite series
+
+.. math::
+ \begin{align}
+ 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{align}
+
+Here we are using the fact that :math:`\epsilon^2 = 0`.
+
+A `Jet <https://en.wikipedia.org/wiki/Jet_(mathematics)>`_ is a
+:math:`n`-dimensional dual number, where we augment the real numbers
+with :math:`n` infinitesimal units :math:`\epsilon_i,\ i=1,...,n` with
+the property that :math:`\forall i, j\ :\epsilon_i\epsilon_j = 0`. Then
+a Jet consists of a *real* part :math:`a` and a :math:`n`-dimensional
+*infinitesimal* part :math:`\mathbf{v}`, i.e.,
+
+.. math::
+ x = a + \sum_j v_{j} \epsilon_j
+
+The summation notation gets tedious, so we will also just write
+
+.. math::
+ x = a + \mathbf{v}.
+
+where the :math:`\epsilon_i`'s are implicit. Then, using the same
+Taylor series expansion used above, we can see that:
+
+.. math::
+
+ f(a + \mathbf{v}) = f(a) + Df(a) \mathbf{v}.
+
+Similarly for a multivariate function
+:math:`f:\mathbb{R}^{n}\rightarrow \mathbb{R}^m`, evaluated on
+:math:`x_i = a_i + \mathbf{v}_i,\ \forall i = 1,...,n`:
+
+.. math::
+ f(x_1,..., x_n) = f(a_1, ..., a_n) + \sum_i D_i f(a_1, ..., a_n) \mathbf{v}_i
+
+So if each :math:`\mathbf{v}_i = e_i` were the :math:`i^{\text{th}}`
+standard basis vector, then, the above expression would simplify to
+
+.. math::
+ f(x_1,..., x_n) = f(a_1, ..., a_n) + \sum_i D_i f(a_1, ..., a_n) \epsilon_i
+
+and we can extract the coordinates of the Jacobian by inspecting the
+coefficients of :math:`\epsilon_i`.
+
+Implementing Jets
+-----------------
+
+In order for the above to work in practice, we will need the ability
+to evaluate an arbitrary function :math:`f` not just on real numbers
+but also on dual numbers, but one does not usually evaluate functions
+by evaluating their Taylor expansions,
+
+This is where C++ templates and operator overloading comes into
+play. The following code fragment has a simple implementation of a
+``Jet`` and some operators/functions that operate on them.
+
+.. code-block:: c++
+
+ 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);
+ }
+
+
+With these overloaded functions in hand, we can now call
+``Rat43CostFunctor`` with an array of Jets instead of doubles. Putting
+that together with appropriately initialized Jets allows us to compute
+the Jacobian as follows:
+
+.. code-block:: c++
+
+ 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_;
+ };
+
+Indeed, this is essentially how :class:`AutoDiffCostFunction` works.
+
+
+Pitfalls
+========
+
+Automatic differentiation frees the user from the burden of computing
+and reasoning about the symbolic expressions for the Jacobians, but
+this freedom comes at a cost. For example consider the following
+simple functor:
+
+.. code-block:: c++
+
+ 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;
+ }
+ };
+
+Looking at the code for the residual computation, one does not foresee
+any problems. However, if we look at the analytical expressions for
+the Jacobian:
+
+.. math::
+
+ 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}}
+
+we find that it is an indeterminate form at :math:`x_0 = 0, x_1 =
+0`.
+
+There is no single solution to this problem. In some cases one needs
+to reason explicitly about the points where indeterminacy may occur
+and use alternate expressions using `L'Hopital's rule
+<https://en.wikipedia.org/wiki/L'H%C3%B4pital's_rule>`_ (see for
+example some of the conversion routines in `rotation.h
+<https://github.com/ceres-solver/ceres-solver/blob/master/include/ceres/rotation.h>`_. In
+other cases, one may need to regularize the expressions to eliminate
+these points.