mirror of
https://github.com/ceres-solver/ceres-solver.git
synced 2026-08-29 16:40:38 +08:00
Add an article on interfacing with automatic differentiation.
Doing this also necessitated some re-organization of the derivatives article into chapters and some minor edits. Change-Id: Ic08e83af138817173caa80a52a9e72707cd57512
This commit is contained in:
Sameer Agarwal
@@ -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 differentition 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 implict. 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.
|
||||
Reference in New Issue
Block a user