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/numerical_derivatives.rst b/docs/source/numerical_derivatives.rst
new file mode 100644
index 0000000..57b46bf
--- /dev/null
+++ b/docs/source/numerical_derivatives.rst
@@ -0,0 +1,403 @@
+.. default-domain:: cpp
+
+.. cpp:namespace:: ceres
+
+.. _chapter-numerical_derivatives:
+
+===================
+Numeric derivatives
+===================
+
+The other extreme from using analytic derivatives is to use numeric
+derivatives. The key observation here is that the process of
+differentiating a function :math:`f(x)` w.r.t :math:`x` can be written
+as the limiting process:
+
+.. math::
+ Df(x) = \lim_{h \rightarrow 0} \frac{f(x + h) - f(x)}{h}
+
+
+Forward Differences
+===================
+
+Now of course one cannot perform the limiting operation numerically on
+a computer so we do the next best thing, which is to choose a small
+value of :math:`h` and approximate the derivative as
+
+.. math::
+ Df(x) \approx \frac{f(x + h) - f(x)}{h}
+
+
+The above formula is the simplest most basic form of numeric
+differentiation. It is known as the *Forward Difference* formula.
+
+So how would one go about constructing a numerically differentiated
+version of ``Rat43Analytic`` (`Rat43
+<http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml>`_) in
+Ceres Solver. This is done in two steps:
+
+ 1. Define *Functor* that given the parameter values will evaluate the
+ residual for a given :math:`(x,y)`.
+ 2. Construct a :class:`CostFunction` by using
+ :class:`NumericDiffCostFunction` to wrap an instance of
+ ``Rat43CostFunctor``.
+
+.. code-block:: c++
+
+ 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));
+
+This is about the minimum amount of work one can expect to do to
+define the cost function. The only thing that the user needs to do is
+to make sure that the evaluation of the residual is implemented
+correctly and efficiently.
+
+Before going further, it is instructive to get an estimate of the
+error in the forward difference formula. We do this by considering the
+`Taylor expansion <https://en.wikipedia.org/wiki/Taylor_series>`_ of
+:math:`f` near :math:`x`.
+
+.. math::
+ \begin{align}
+ f(x+h) &= f(x) + h Df(x) + \frac{h^2}{2!} D^2f(x) +
+ \frac{h^3}{3!}D^3f(x) + \cdots \\
+ Df(x) &= \frac{f(x + h) - f(x)}{h} - \left [\frac{h}{2!}D^2f(x) +
+ \frac{h^2}{3!}D^3f(x) + \cdots \right]\\
+ Df(x) &= \frac{f(x + h) - f(x)}{h} + O(h)
+ \end{align}
+
+i.e., the error in the forward difference formula is
+:math:`O(h)` [#f4]_.
+
+
+Implementation Details
+----------------------
+
+:class:`NumericDiffCostFunction` implements a generic algorithm to
+numerically differentiate a given functor. While the actual
+implementation of :class:`NumericDiffCostFunction` is complicated, the
+net result is a :class:`CostFunction` that roughly looks something
+like the following:
+
+.. code-block:: c++
+
+ class Rat43NumericDiffForward : public SizedCostFunction<1,4> {
+ public:
+ Rat43NumericDiffForward(const Rat43Functor* functor) : functor_(functor) {}
+ virtual ~Rat43NumericDiffForward() {}
+ virtual bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const {
+ functor_(parameters[0], residuals);
+ if (!jacobians) return true;
+ double* jacobian = jacobians[0];
+ if (!jacobian) return true;
+
+ const double f = residuals[0];
+ double parameters_plus_h[4];
+ for (int i = 0; i < 4; ++i) {
+ std::copy(parameters, parameters + 4, parameters_plus_h);
+ const double kRelativeStepSize = 1e-6;
+ const double h = std::abs(parameters[i]) * kRelativeStepSize;
+ parameters_plus_h[i] += h;
+ double f_plus;
+ functor_(parameters_plus_h, &f_plus);
+ jacobian[i] = (f_plus - f) / h;
+ }
+ return true;
+ }
+
+ private:
+ std::unique_ptr<Rat43Functor> functor_;
+ };
+
+
+Note the choice of step size :math:`h` in the above code, instead of
+an absolute step size which is the same for all parameters, we use a
+relative step size of :math:`\text{kRelativeStepSize} = 10^{-6}`. This
+gives better derivative estimates than an absolute step size [#f2]_
+[#f3]_. This choice of step size only works for parameter values that
+are not close to zero. So the actual implementation of
+:class:`NumericDiffCostFunction`, uses a more complex step size
+selection logic, where close to zero, it switches to a fixed step
+size.
+
+
+Central Differences
+===================
+
+:math:`O(h)` error in the Forward Difference formula is okay but not
+great. A better method is to use the *Central Difference* formula:
+
+.. math::
+ Df(x) \approx \frac{f(x + h) - f(x - h)}{2h}
+
+Notice that if the value of :math:`f(x)` is known, the Forward
+Difference formula only requires one extra evaluation, but the Central
+Difference formula requires two evaluations, making it twice as
+expensive. So is the extra evaluation worth it?
+
+To answer this question, we again compute the error of approximation
+in the central difference formula:
+
+.. math::
+ \begin{align}
+ f(x + h) &= f(x) + h Df(x) + \frac{h^2}{2!}
+ D^2f(x) + \frac{h^3}{3!} D^3f(x) + \frac{h^4}{4!} D^4f(x) + \cdots\\
+ f(x - h) &= f(x) - h Df(x) + \frac{h^2}{2!}
+ D^2f(x) - \frac{h^3}{3!} D^3f(c_2) + \frac{h^4}{4!} D^4f(x) +
+ \cdots\\
+ Df(x) & = \frac{f(x + h) - f(x - h)}{2h} + \frac{h^2}{3!}
+ D^3f(x) + \frac{h^4}{5!}
+ D^5f(x) + \cdots \\
+ Df(x) & = \frac{f(x + h) - f(x - h)}{2h} + O(h^2)
+ \end{align}
+
+The error of the Central Difference formula is :math:`O(h^2)`, i.e.,
+the error goes down quadratically whereas the error in the Forward
+Difference formula only goes down linearly.
+
+Using central differences instead of forward differences in Ceres
+Solver is a simple matter of changing a template argument to
+:class:`NumericDiffCostFunction` as follows:
+
+.. code-block:: c++
+
+ CostFunction* cost_function =
+ new NumericDiffCostFunction<Rat43CostFunctor, CENTRAL, 1, 4>(
+ new Rat43CostFunctor(x, y));
+
+But what do these differences in the error mean in practice? To see
+this, consider the problem of evaluating the derivative of the
+univariate function
+
+.. math::
+ f(x) = \frac{e^x}{\sin x - x^2},
+
+at :math:`x = 1.0`.
+
+It is easy to determine that :math:`Df(1.0) =
+140.73773557129658`. Using this value as reference, we can now compute
+the relative error in the forward and central difference formulae as a
+function of the absolute step size and plot them.
+
+.. figure:: forward_central_error.png
+ :figwidth: 100%
+ :align: center
+
+Reading the graph from right to left, a number of things stand out in
+the above graph:
+
+ 1. The graph for both formulae have two distinct regions. At first,
+ starting from a large value of :math:`h` the error goes down as
+ the effect of truncating the Taylor series dominates, but as the
+ value of :math:`h` continues to decrease, the error starts
+ increasing again as roundoff error starts to dominate the
+ computation. So we cannot just keep on reducing the value of
+ :math:`h` to get better estimates of :math:`Df`. The fact that we
+ are using finite precision arithmetic becomes a limiting factor.
+ 2. Forward Difference formula is not a great method for evaluating
+ derivatives. Central Difference formula converges much more
+ quickly to a more accurate estimate of the derivative with
+ decreasing step size. So unless the evaluation of :math:`f(x)` is
+ so expensive that you absolutely cannot afford the extra
+ evaluation required by central differences, **do not use the
+ Forward Difference formula**.
+ 3. Neither formula works well for a poorly chosen value of :math:`h`.
+
+
+Ridders' Method
+===============
+
+So, can we get better estimates of :math:`Df` without requiring such
+small values of :math:`h` that we start hitting floating point
+roundoff errors?
+
+One possible approach is to find a method whose error goes down faster
+than :math:`O(h^2)`. This can be done by applying `Richardson
+Extrapolation
+<https://en.wikipedia.org/wiki/Richardson_extrapolation>`_ to the
+problem of differentiation. This is also known as *Ridders' Method*
+[Ridders]_.
+
+Let us recall, the error in the central differences formula.
+
+.. math::
+ \begin{align}
+ Df(x) & = \frac{f(x + h) - f(x - h)}{2h} + \frac{h^2}{3!}
+ D^3f(x) + \frac{h^4}{5!}
+ D^5f(x) + \cdots\\
+ & = \frac{f(x + h) - f(x - h)}{2h} + K_2 h^2 + K_4 h^4 + \cdots
+ \end{align}
+
+The key thing to note here is that the terms :math:`K_2, K_4, ...`
+are independent of :math:`h` and only depend on :math:`x`.
+
+Let us now define:
+
+.. math::
+
+ A(1, m) = \frac{f(x + h/2^{m-1}) - f(x - h/2^{m-1})}{2h/2^{m-1}}.
+
+Then observe that
+
+.. math::
+
+ Df(x) = A(1,1) + K_2 h^2 + K_4 h^4 + \cdots
+
+and
+
+.. math::
+
+ Df(x) = A(1, 2) + K_2 (h/2)^2 + K_4 (h/2)^4 + \cdots
+
+Here we have halved the step size to obtain a second central
+differences estimate of :math:`Df(x)`. Combining these two estimates,
+we get:
+
+.. math::
+
+ Df(x) = \frac{4 A(1, 2) - A(1,1)}{4 - 1} + O(h^4)
+
+which is an approximation of :math:`Df(x)` with truncation error that
+goes down as :math:`O(h^4)`. But we do not have to stop here. We can
+iterate this process to obtain even more accurate estimates as
+follows:
+
+.. math::
+
+ A(n, m) = \begin{cases}
+ \frac{\displaystyle f(x + h/2^{m-1}) - f(x -
+ h/2^{m-1})}{\displaystyle 2h/2^{m-1}} & n = 1 \\
+ \frac{\displaystyle 4^{n-1} A(n - 1, m + 1) - A(n - 1, m)}{\displaystyle 4^{n-1} - 1} & n > 1
+ \end{cases}
+
+It is straightforward to show that the approximation error in
+:math:`A(n, 1)` is :math:`O(h^{2n})`. To see how the above formula can
+be implemented in practice to compute :math:`A(n,1)` it is helpful to
+structure the computation as the following tableau:
+
+.. math::
+ \begin{array}{ccccc}
+ A(1,1) & A(1, 2) & A(1, 3) & A(1, 4) & \cdots\\
+ & A(2, 1) & A(2, 2) & A(2, 3) & \cdots\\
+ & & A(3, 1) & A(3, 2) & \cdots\\
+ & & & A(4, 1) & \cdots \\
+ & & & & \ddots
+ \end{array}
+
+So, to compute :math:`A(n, 1)` for increasing values of :math:`n` we
+move from the left to the right, computing one column at a
+time. Assuming that the primary cost here is the evaluation of the
+function :math:`f(x)`, the cost of computing a new column of the above
+tableau is two function evaluations. Since the cost of evaluating
+:math:`A(1, n)`, requires evaluating the central difference formula
+for step size of :math:`2^{1-n}h`
+
+Applying this method to :math:`f(x) = \frac{e^x}{\sin x - x^2}`
+starting with a fairly large step size :math:`h = 0.01`, we get:
+
+.. math::
+ \begin{array}{rrrrr}
+ 141.678097131 &140.971663667 &140.796145400 &140.752333523 &140.741384778\\
+ &140.736185846 &140.737639311 &140.737729564 &140.737735196\\
+ & &140.737736209 &140.737735581 &140.737735571\\
+ & & &140.737735571 &140.737735571\\
+ & & & &140.737735571\\
+ \end{array}
+
+Compared to the *correct* value :math:`Df(1.0) = 140.73773557129658`,
+:math:`A(5, 1)` has a relative error of :math:`10^{-13}`. For
+comparison, the relative error for the central difference formula with
+the same stepsize (:math:`0.01/2^4 = 0.000625`) is :math:`10^{-5}`.
+
+The above tableau is the basis of Ridders' method for numeric
+differentiation. The full implementation is an adaptive scheme that
+tracks its own estimation error and stops automatically when the
+desired precision is reached. Of course it is more expensive than the
+forward and central difference formulae, but is also significantly
+more robust and accurate.
+
+Using Ridder's method instead of forward or central differences in
+Ceres is again a simple matter of changing a template argument to
+:class:`NumericDiffCostFunction` as follows:
+
+.. code-block:: c++
+
+ CostFunction* cost_function =
+ new NumericDiffCostFunction<Rat43CostFunctor, RIDDERS, 1, 4>(
+ new Rat43CostFunctor(x, y));
+
+The following graph shows the relative error of the three methods as a
+function of the absolute step size. For Ridders's method we assume
+that the step size for evaluating :math:`A(n,1)` is :math:`2^{1-n}h`.
+
+.. figure:: forward_central_ridders_error.png
+ :figwidth: 100%
+ :align: center
+
+Using the 10 function evaluations that are needed to compute
+:math:`A(5,1)` we are able to approximate :math:`Df(1.0)` about a 1000
+times better than the best central differences estimate. To put these
+numbers in perspective, machine epsilon for double precision
+arithmetic is :math:`\approx 2.22 \times 10^{-16}`.
+
+Going back to ``Rat43``, let us also look at the runtime cost of the
+various methods for computing numeric derivatives.
+
+========================== =========
+CostFunction Time (ns)
+========================== =========
+Rat43Analytic 255
+Rat43AnalyticOptimized 92
+Rat43NumericDiffForward 262
+Rat43NumericDiffCentral 517
+Rat43NumericDiffRidders 3760
+========================== =========
+
+As expected, Central Differences is about twice as expensive as
+Forward Differences and the remarkable accuracy improvements of
+Ridders' method cost an order of magnitude more runtime.
+
+Recommendations
+===============
+
+Numeric differentiation should be used when you cannot compute the
+derivatives either analytically or using automatic differentiation. This
+is usually the case when you are calling an external library or
+function whose analytic form you do not know or even if you do, you
+are not in a position to re-write it in a manner required to use
+:ref:`chapter-automatic_derivatives`.
+
+
+When using numeric differentiation, use at least Central Differences,
+and if execution time is not a concern or the objective function is
+such that determining a good static relative step size is hard,
+Ridders' method is recommended.
+
+.. rubric:: Footnotes
+
+.. [#f2] `Numerical Differentiation
+ <https://en.wikipedia.org/wiki/Numerical_differentiation#Practical_considerations_using_floating_point_arithmetic>`_
+.. [#f3] [Press]_ Numerical Recipes, Section 5.7
+.. [#f4] In asymptotic error analysis, an error of :math:`O(h^k)`
+ means that the absolute-value of the error is at most some
+ constant times :math:`h^k` when :math:`h` is close enough to
+ :math:`0`.