docs/source/modeling_faqs.rst

   1 .. _chapter-modeling_faqs:
   2
   3 .. default-domain:: cpp
   4
   5 .. cpp:namespace:: ceres
   6
   7 ========
   8 Modeling
   9 ========
  10
  11 #. Use analytical/automatic derivatives.
  12
  13    This is the single most important piece of advice we can give to
  14    you. It is tempting to take the easy way out and use numeric
  15    differentiation. This is a bad idea. Numeric differentiation is
  16    slow, ill-behaved, hard to get right, and results in poor
  17    convergence behaviour.
  18
  19    Ceres allows the user to define templated functors which will
  20    be automatically differentiated. For most situations this is enough
  21    and we recommend using this facility. In some cases the derivatives
  22    are simple enough or the performance considerations are such that
  23    the overhead of automatic differentiation is too much. In such
  24    cases, analytic derivatives are recommended.
  25
  26    The use of numerical derivatives should be a measure of last
  27    resort, where it is simply not possible to write a templated
  28    implementation of the cost function.
  29
  30    In many cases it is not possible to do analytic or automatic
  31    differentiation of the entire cost function, but it is generally
  32    the case that it is possible to decompose the cost function into
  33    parts that need to be numerically differentiated and parts that can
  34    be automatically or analytically differentiated.
  35
  36    To this end, Ceres has extensive support for mixing analytic,
  37    automatic and numeric differentiation. See
  38    :class:`CostFunctionToFunctor`.
  39
  40 #. When using Quaternions,  consider using :class:`QuaternionParameterization`.
  41
  42    `Quaternions <https://en.wikipedia.org/wiki/Quaternion>`_ are a
  43    four dimensional parameterization of the space of three dimensional
  44    rotations :math:`SO(3)`.  However, the :math:`SO(3)` is a three
  45    dimensional set, and so is the tangent space of a
  46    Quaternion. Therefore, it is sometimes (not always) benefecial to
  47    associate a local parameterization with parameter blocks
  48    representing a Quaternion. Assuming that the order of entries in
  49    your parameter block is :math:`w,x,y,z`, you can use
  50    :class:`QuaternionParameterization`.
  51
  52    .. NOTE::
  53
  54      If you are using `Eigen's Quaternion
  55      <http://eigen.tuxfamily.org/dox/classEigen_1_1Quaternion.html>`_
  56      object, whose layout is :math:`x,y,z,w`, then you should use
  57      :class:`EigenQuaternionParameterization`.
  58
  59
  60 #. How do I solve problems with general linear & non-linear
  61    **inequality** constraints with Ceres Solver?
  62
  63    Currently, Ceres Solver only supports upper and lower bounds
  64    constraints on the parameter blocks.
  65
  66    A crude way of dealing with inequality constraints is have one or
  67    more of your cost functions check if the inequalities you are
  68    interested in are satisfied, and if not return false instead of
  69    true. This will prevent the solver from ever stepping into an
  70    infeasible region.
  71
  72    This requires that the starting point for the optimization be a
  73    feasible point.  You also risk pre-mature convergence using this
  74    method.
  75
  76 #. How do I solve problems with general linear & non-linear **equality**
  77    constraints with Ceres Solver?
  78
  79    There is no built in support in ceres for solving problems with
  80    equality constraints.  Currently, Ceres Solver only supports upper
  81    and lower bounds constraints on the parameter blocks.
  82
  83    The trick described above for dealing with inequality
  84    constraints will **not** work for equality constraints.
  85
  86 #. How do I set one or more components of a parameter block constant?
  87
  88    Using :class:`SubsetParameterization`.
  89
  90 #. Putting `Inverse Function Theorem
  91    <http://en.wikipedia.org/wiki/Inverse_function_theorem>`_ to use.
  92
  93    Every now and then we have to deal with functions which cannot be
  94    evaluated analytically. Computing the Jacobian in such cases is
  95    tricky. A particularly interesting case is where the inverse of the
  96    function is easy to compute analytically. An example of such a
  97    function is the Coordinate transformation between the `ECEF
  98    <http://en.wikipedia.org/wiki/ECEF>`_ and the `WGS84
  99    <http://en.wikipedia.org/wiki/World_Geodetic_System>`_ where the
 100    conversion from WGS84 to ECEF is analytic, but the conversion
 101    back to WGS84 uses an iterative algorithm. So how do you compute the
 102    derivative of the ECEF to WGS84 transformation?
 103
 104    One obvious approach would be to numerically
 105    differentiate the conversion function. This is not a good idea. For
 106    one, it will be slow, but it will also be numerically quite
 107    bad.
 108
 109    Turns out you can use the `Inverse Function Theorem
 110    <http://en.wikipedia.org/wiki/Inverse_function_theorem>`_ in this
 111    case to compute the derivatives more or less analytically.
 112
 113    The key result here is. If :math:`x = f^{-1}(y)`, and :math:`Df(x)`
 114    is the invertible Jacobian of :math:`f` at :math:`x`. Then the
 115    Jacobian :math:`Df^{-1}(y) = [Df(x)]^{-1}`, i.e., the Jacobian of
 116    the :math:`f^{-1}` is the inverse of the Jacobian of :math:`f`.
 117
 118    Algorithmically this means that given :math:`y`, compute :math:`x =
 119    f^{-1}(y)` by whatever means you can. Evaluate the Jacobian of
 120    :math:`f` at :math:`x`. If the Jacobian matrix is invertible, then
 121    its inverse is the Jacobian of :math:`f^{-1}(y)` at  :math:`y`.
 122
 123    One can put this into practice with the following code fragment.
 124
 125    .. code-block:: c++
 126
 127       Eigen::Vector3d ecef; // Fill some values
 128       // Iterative computation.
 129       Eigen::Vector3d lla = ECEFToLLA(ecef);
 130       // Analytic derivatives
 131       Eigen::Matrix3d lla_to_ecef_jacobian = LLAToECEFJacobian(lla);
 132       bool invertible;
 133       Eigen::Matrix3d ecef_to_lla_jacobian;
 134       lla_to_ecef_jacobian.computeInverseWithCheck(ecef_to_lla_jacobian, invertible);