docs/source/interfacing_with_autodiff.rst - ceres-solver - Git at Google

 .. default-domain:: cpp

 .. cpp:namespace:: ceres

 .. _chapter-interfacing_with_automatic_differentiation:

 Interfacing with Automatic Differentiation
 ==========================================

 Automatic differentiation is straightforward to use in cases where an
 explicit expression for the cost function is available. But this is
 not always possible. Often one has to interface with external routines
 or data. In this chapter we will consider a number of different ways
 of doing so.

 To do this, we will consider the problem of finding parameters
 :math:`\theta` and :math:`t` that solve an optimization problem of the
 form:

 .. math::
    \min & \quad \sum_i \left \|y_i - f\left (\|q_{i}\|^2\right) q_i
    \right \|^2\\
    \text{such that} & \quad q_i = R(\theta) x_i + t

 Here, :math:`R` is a two dimensional rotation matrix parameterized
 using the angle :math:`\theta` and :math:`t` is a two dimensional
 vector. :math:`f` is an external distortion function.

 We begin by considering the case, where we have a templated function
 :code:`TemplatedComputeDistortion` that can compute the function
 :math:`f`. Then the implementation of the corresponding residual
 functor is straightforward and will look as follows:

 .. code-block:: c++
    :emphasize-lines: 21

    template <typename T> T TemplatedComputeDistortion(const T r2) {
      const double k1 = 0.0082;
      const double k2 = 0.000023;
      return 1.0 + k1 * r2 + k2 * r2 * r2;
    }

    struct Affine2DWithDistortion {
      Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
        x[0] = x_in[0];
        x[1] = x_in[1];
        y[0] = y_in[0];
        y[1] = y_in[1];
      }

      template <typename T>
      bool operator()(const T* theta,
                      const T* t,
                      T* residuals) const {
        const T q_0 =  cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
        const T q_1 =  sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
        const T f = TemplatedComputeDistortion(q_0 * q_0 + q_1 * q_1);
        residuals[0] = y[0] - f * q_0;
        residuals[1] = y[1] - f * q_1;
        return true;
      }

      double x[2];
      double y[2];
    };

 So far so good, but let us now consider three ways of defining
 :math:`f` which are not directly amenable to being used with automatic
 differentiation:

 #. A non-templated function that evaluates its value.
 #. A function that evaluates its value and derivative.
 #. A function that is defined as a table of values to be interpolated.

 We will consider them in turn below.

 A function that returns its value
 ----------------------------------

 Suppose we were given a function :code:`ComputeDistortionValue` with
 the following signature

 .. code-block:: c++

    double ComputeDistortionValue(double r2);

 that computes the value of :math:`f`. The actual implementation of the
 function does not matter. Interfacing this function with
 :code:`Affine2DWithDistortion` is a three step process:

 1. Wrap :code:`ComputeDistortionValue` into a functor
    :code:`ComputeDistortionValueFunctor`.
 2. Numerically differentiate :code:`ComputeDistortionValueFunctor`
    using :class:`NumericDiffCostFunction` to create a
    :class:`CostFunction`.
 3. Wrap the resulting :class:`CostFunction` object using
    :class:`CostFunctionToFunctor`. The resulting object is a functor
    with a templated :code:`operator()` method, which pipes the
    Jacobian computed by :class:`NumericDiffCostFunction` into the
    appropriate :code:`Jet` objects.

 An implementation of the above three steps looks as follows:

 .. code-block:: c++
    :emphasize-lines: 15,16,17,18,19,20, 29

    struct ComputeDistortionValueFunctor {
      bool operator()(const double* r2, double* value) const {
        *value = ComputeDistortionValue(r2[0]);
        return true;
      }
    };

    struct Affine2DWithDistortion {
      Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
        x[0] = x_in[0];
        x[1] = x_in[1];
        y[0] = y_in[0];
        y[1] = y_in[1];

        compute_distortion = std::make_unique<ceres::CostFunctionToFunctor<1, 1>>(
          std::make_unique<ceres::NumericDiffCostFunction<
                ComputeDistortionValueFunctor
              , ceres::CENTRAL, 1, 1
            >
          >()
        );
      }

      template <typename T>
      bool operator()(const T* theta, const T* t, T* residuals) const {
        const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
        const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
        const T r2 = q_0 * q_0 + q_1 * q_1;
        T f;
        (*compute_distortion)(&r2, &f);
        residuals[0] = y[0] - f * q_0;
        residuals[1] = y[1] - f * q_1;
        return true;
      }

      double x[2];
      double y[2];
      std::unique_ptr<ceres::CostFunctionToFunctor<1, 1>> compute_distortion;
    };


 A function that returns its value and derivative
 ------------------------------------------------

 Now suppose we are given a function :code:`ComputeDistortionValue`
 that is able to compute its value and optionally its Jacobian on demand
 and has the following signature:

 .. code-block:: c++

    void ComputeDistortionValueAndJacobian(double r2,
                                           double* value,
                                           double* jacobian);

 Again, the actual implementation of the function does not
 matter. Interfacing this function with :code:`Affine2DWithDistortion`
 is a two step process:

 1. Wrap :code:`ComputeDistortionValueAndJacobian` into a
    :class:`CostFunction` object which we call
    :code:`ComputeDistortionFunction`.
 2. Wrap the resulting :class:`ComputeDistortionFunction` object using
    :class:`CostFunctionToFunctor`. The resulting object is a functor
    with a templated :code:`operator()` method, which pipes the
    Jacobian computed by :class:`NumericDiffCostFunction` into the
    appropriate :code:`Jet` objects.

 The resulting code will look as follows:

 .. code-block:: c++
    :emphasize-lines: 21,22, 33

    class ComputeDistortionFunction : public ceres::SizedCostFunction<1, 1> {
     public:
      virtual bool Evaluate(double const* const* parameters,
                            double* residuals,
                            double** jacobians) const {
        if (!jacobians) {
          ComputeDistortionValueAndJacobian(parameters[0][0], residuals, nullptr);
        } else {
          ComputeDistortionValueAndJacobian(parameters[0][0], residuals, jacobians[0]);
        }
        return true;
      }
    };

    struct Affine2DWithDistortion {
      Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
        x[0] = x_in[0];
        x[1] = x_in[1];
        y[0] = y_in[0];
        y[1] = y_in[1];
        compute_distortion.reset(
            new ceres::CostFunctionToFunctor<1, 1>(new ComputeDistortionFunction));
      }

      template <typename T>
      bool operator()(const T* theta,
                      const T* t,
                      T* residuals) const {
        const T q_0 =  cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
        const T q_1 =  sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
        const T r2 = q_0 * q_0 + q_1 * q_1;
        T f;
        (*compute_distortion)(&r2, &f);
        residuals[0] = y[0] - f * q_0;
        residuals[1] = y[1] - f * q_1;
        return true;
      }

      double x[2];
      double y[2];
      std::unique_ptr<ceres::CostFunctionToFunctor<1, 1> > compute_distortion;
    };


 A function that is defined as a table of values
 -----------------------------------------------

 The third and final case we will consider is where the function
 :math:`f` is defined as a table of values on the interval :math:`[0,
 100)`, with a value for each integer.

 .. code-block:: c++

    vector<double> distortion_values;

 There are many ways of interpolating a table of values. Perhaps the
 simplest and most common method is linear interpolation. But it is not
 a great idea to use linear interpolation because the interpolating
 function is not differentiable at the sample points.

 A simple (well behaved) differentiable interpolation is the `Cubic
 Hermite Spline
 <http://en.wikipedia.org/wiki/Cubic_Hermite_spline>`_. Ceres Solver
 ships with routines to perform Cubic & Bi-Cubic interpolation that is
 automatic differentiation friendly.

 Using Cubic interpolation requires first constructing a
 :class:`Grid1D` object to wrap the table of values and then
 constructing a :class:`CubicInterpolator` object using it.

 The resulting code will look as follows:

 .. code-block:: c++
    :emphasize-lines: 10,11,12,13, 24, 32,33

    struct Affine2DWithDistortion {
      Affine2DWithDistortion(const double x_in[2],
                             const double y_in[2],
                             const std::vector<double>& distortion_values) {
        x[0] = x_in[0];
        x[1] = x_in[1];
        y[0] = y_in[0];
        y[1] = y_in[1];

        grid.reset(new ceres::Grid1D<double, 1>(
            &distortion_values[0], 0, distortion_values.size()));
        compute_distortion.reset(
            new ceres::CubicInterpolator<ceres::Grid1D<double, 1> >(*grid));
      }

      template <typename T>
      bool operator()(const T* theta,
                      const T* t,
                      T* residuals) const {
        const T q_0 =  cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
        const T q_1 =  sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
        const T r2 = q_0 * q_0 + q_1 * q_1;
        T f;
        compute_distortion->Evaluate(r2, &f);
        residuals[0] = y[0] - f * q_0;
        residuals[1] = y[1] - f * q_1;
        return true;
      }

      double x[2];
      double y[2];
      std::unique_ptr<ceres::Grid1D<double, 1> > grid;
      std::unique_ptr<ceres::CubicInterpolator<ceres::Grid1D<double, 1> > > compute_distortion;
    };

 In the above example we used :class:`Grid1D` and
 :class:`CubicInterpolator` to interpolate a one dimensional table of
 values. :class:`Grid2D` combined with :class:`CubicInterpolator` lets
 the user to interpolate two dimensional tables of values. Note that
 neither :class:`Grid1D` or :class:`Grid2D` are limited to scalar
 valued functions, they also work with vector valued functions.
	.. default-domain:: cpp

	.. cpp:namespace:: ceres

	.. _chapter-interfacing_with_automatic_differentiation:

	Interfacing with Automatic Differentiation
	==========================================

	Automatic differentiation is straightforward to use in cases where an
	explicit expression for the cost function is available. But this is
	not always possible. Often one has to interface with external routines
	or data. In this chapter we will consider a number of different ways
	of doing so.

	To do this, we will consider the problem of finding parameters
	:math:`\theta` and :math:`t` that solve an optimization problem of the
	form:

	.. math::
	\min & \quad \sum_i \left \\|y_i - f\left (\\|q_{i}\\|^2\right) q_i
	\right \\|^2\\
	\text{such that} & \quad q_i = R(\theta) x_i + t

	Here, :math:`R` is a two dimensional rotation matrix parameterized
	using the angle :math:`\theta` and :math:`t` is a two dimensional
	vector. :math:`f` is an external distortion function.

	We begin by considering the case, where we have a templated function
	:code:`TemplatedComputeDistortion` that can compute the function
	:math:`f`. Then the implementation of the corresponding residual
	functor is straightforward and will look as follows:

	.. code-block:: c++
	:emphasize-lines: 21

	template <typename T> T TemplatedComputeDistortion(const T r2) {
	const double k1 = 0.0082;
	const double k2 = 0.000023;
	return 1.0 + k1 * r2 + k2 * r2 * r2;
	}

	struct Affine2DWithDistortion {
	Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
	x[0] = x_in[0];
	x[1] = x_in[1];
	y[0] = y_in[0];
	y[1] = y_in[1];
	}

	template <typename T>
	bool operator()(const T* theta,
	const T* t,
	T* residuals) const {
	const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
	const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
	const T f = TemplatedComputeDistortion(q_0 * q_0 + q_1 * q_1);
	residuals[0] = y[0] - f * q_0;
	residuals[1] = y[1] - f * q_1;
	return true;
	}

	double x[2];
	double y[2];
	};

	So far so good, but let us now consider three ways of defining
	:math:`f` which are not directly amenable to being used with automatic
	differentiation:

	#. A non-templated function that evaluates its value.
	#. A function that evaluates its value and derivative.
	#. A function that is defined as a table of values to be interpolated.

	We will consider them in turn below.

	A function that returns its value
	----------------------------------

	Suppose we were given a function :code:`ComputeDistortionValue` with
	the following signature

	.. code-block:: c++

	double ComputeDistortionValue(double r2);

	that computes the value of :math:`f`. The actual implementation of the
	function does not matter. Interfacing this function with
	:code:`Affine2DWithDistortion` is a three step process:

	1. Wrap :code:`ComputeDistortionValue` into a functor
	:code:`ComputeDistortionValueFunctor`.
	2. Numerically differentiate :code:`ComputeDistortionValueFunctor`
	using :class:`NumericDiffCostFunction` to create a
	:class:`CostFunction`.
	3. Wrap the resulting :class:`CostFunction` object using
	:class:`CostFunctionToFunctor`. The resulting object is a functor
	with a templated :code:`operator()` method, which pipes the
	Jacobian computed by :class:`NumericDiffCostFunction` into the
	appropriate :code:`Jet` objects.

	An implementation of the above three steps looks as follows:

	.. code-block:: c++
	:emphasize-lines: 15,16,17,18,19,20, 29

	struct ComputeDistortionValueFunctor {
	bool operator()(const double* r2, double* value) const {
	*value = ComputeDistortionValue(r2[0]);
	return true;
	}
	};

	struct Affine2DWithDistortion {
	Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
	x[0] = x_in[0];
	x[1] = x_in[1];
	y[0] = y_in[0];
	y[1] = y_in[1];

	compute_distortion = std::make_unique<ceres::CostFunctionToFunctor<1, 1>>(
	std::make_unique<ceres::NumericDiffCostFunction<
	ComputeDistortionValueFunctor
	, ceres::CENTRAL, 1, 1
	>
	>()
	);
	}

	template <typename T>
	bool operator()(const T* theta, const T* t, T* residuals) const {
	const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
	const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
	const T r2 = q_0 * q_0 + q_1 * q_1;
	T f;
	(*compute_distortion)(&r2, &f);
	residuals[0] = y[0] - f * q_0;
	residuals[1] = y[1] - f * q_1;
	return true;
	}

	double x[2];
	double y[2];
	std::unique_ptr<ceres::CostFunctionToFunctor<1, 1>> compute_distortion;
	};


	A function that returns its value and derivative
	------------------------------------------------

	Now suppose we are given a function :code:`ComputeDistortionValue`
	that is able to compute its value and optionally its Jacobian on demand
	and has the following signature:

	.. code-block:: c++

	void ComputeDistortionValueAndJacobian(double r2,
	double* value,
	double* jacobian);

	Again, the actual implementation of the function does not
	matter. Interfacing this function with :code:`Affine2DWithDistortion`
	is a two step process:

	1. Wrap :code:`ComputeDistortionValueAndJacobian` into a
	:class:`CostFunction` object which we call
	:code:`ComputeDistortionFunction`.
	2. Wrap the resulting :class:`ComputeDistortionFunction` object using
	:class:`CostFunctionToFunctor`. The resulting object is a functor
	with a templated :code:`operator()` method, which pipes the
	Jacobian computed by :class:`NumericDiffCostFunction` into the
	appropriate :code:`Jet` objects.

	The resulting code will look as follows:

	.. code-block:: c++
	:emphasize-lines: 21,22, 33

	class ComputeDistortionFunction : public ceres::SizedCostFunction<1, 1> {
	public:
	virtual bool Evaluate(double const* const* parameters,
	double* residuals,
	double** jacobians) const {
	if (!jacobians) {
	ComputeDistortionValueAndJacobian(parameters[0][0], residuals, nullptr);
	} else {
	ComputeDistortionValueAndJacobian(parameters[0][0], residuals, jacobians[0]);
	}
	return true;
	}
	};

	struct Affine2DWithDistortion {
	Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
	x[0] = x_in[0];
	x[1] = x_in[1];
	y[0] = y_in[0];
	y[1] = y_in[1];
	compute_distortion.reset(
	new ceres::CostFunctionToFunctor<1, 1>(new ComputeDistortionFunction));
	}

	template <typename T>
	bool operator()(const T* theta,
	const T* t,
	T* residuals) const {
	const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
	const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
	const T r2 = q_0 * q_0 + q_1 * q_1;
	T f;
	(*compute_distortion)(&r2, &f);
	residuals[0] = y[0] - f * q_0;
	residuals[1] = y[1] - f * q_1;
	return true;
	}

	double x[2];
	double y[2];
	std::unique_ptr<ceres::CostFunctionToFunctor<1, 1> > compute_distortion;
	};


	A function that is defined as a table of values
	-----------------------------------------------

	The third and final case we will consider is where the function
	:math:`f` is defined as a table of values on the interval :math:`[0,
	100)`, with a value for each integer.

	.. code-block:: c++

	vector<double> distortion_values;

	There are many ways of interpolating a table of values. Perhaps the
	simplest and most common method is linear interpolation. But it is not
	a great idea to use linear interpolation because the interpolating
	function is not differentiable at the sample points.

	A simple (well behaved) differentiable interpolation is the `Cubic
	Hermite Spline
	<http://en.wikipedia.org/wiki/Cubic_Hermite_spline>`_. Ceres Solver
	ships with routines to perform Cubic & Bi-Cubic interpolation that is
	automatic differentiation friendly.

	Using Cubic interpolation requires first constructing a
	:class:`Grid1D` object to wrap the table of values and then
	constructing a :class:`CubicInterpolator` object using it.

	The resulting code will look as follows:

	.. code-block:: c++
	:emphasize-lines: 10,11,12,13, 24, 32,33

	struct Affine2DWithDistortion {
	Affine2DWithDistortion(const double x_in[2],
	const double y_in[2],
	const std::vector<double>& distortion_values) {
	x[0] = x_in[0];
	x[1] = x_in[1];
	y[0] = y_in[0];
	y[1] = y_in[1];

	grid.reset(new ceres::Grid1D<double, 1>(
	&distortion_values[0], 0, distortion_values.size()));
	compute_distortion.reset(
	new ceres::CubicInterpolator<ceres::Grid1D<double, 1> >(*grid));
	}

	template <typename T>
	bool operator()(const T* theta,
	const T* t,
	T* residuals) const {
	const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
	const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
	const T r2 = q_0 * q_0 + q_1 * q_1;
	T f;
	compute_distortion->Evaluate(r2, &f);
	residuals[0] = y[0] - f * q_0;
	residuals[1] = y[1] - f * q_1;
	return true;
	}

	double x[2];
	double y[2];
	std::unique_ptr<ceres::Grid1D<double, 1> > grid;
	std::unique_ptr<ceres::CubicInterpolator<ceres::Grid1D<double, 1> > > compute_distortion;
	};

	In the above example we used :class:`Grid1D` and
	:class:`CubicInterpolator` to interpolate a one dimensional table of
	values. :class:`Grid2D` combined with :class:`CubicInterpolator` lets
	the user to interpolate two dimensional tables of values. Note that
	neither :class:`Grid1D` or :class:`Grid2D` are limited to scalar
	valued functions, they also work with vector valued functions.