saberbud/ceres-solver
1
0
Fork 0
mirror of https://github.com/ceres-solver/ceres-solver.git synced 2026-08-29 08:34:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add a section on implicit and inverse function theorems

Browse Source 
Change-Id: I0e6c7d2850a33d03aa629579f049ad44a7618621
This commit is contained in:
Sameer Agarwal
2021-09-24 21:44:16 -07:00
parent 686428f5cd
commit 7ef4a12218
3 changed files with 215 additions and 46 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/source/derivatives.rst
+1
View File
@@ -58,3 +58,4 @@ For the rest, read on.
numerical_derivatives
automatic_derivatives
interfacing_with_autodiff
inverse_and_implicit_function_theorems
docs/source/inverse_and_implicit_function_theorems.rst
+214
View File
@@ -0,0 +1,214 @@
.. default-domain:: cpp
.. cpp:namespace:: ceres
.. _chapter-inverse_function_theorem:
==========================================
Using Inverse & Implicit Function Theorems
==========================================
Until now we have considered methods for computing derivatives that
work directly on the function being differentiated. However, this is
not always possible. For example, if the function can only be computed
via an iterative algorithm, or there is no explicit definition of the
function available. In this section we will see how we can use two
basic results from calculus to get around these difficulties.
Inverse Function Theorem
========================
Suppose we wish to evaluate the derivative of a function :math:`f(x)`,
but evaluating :math:`f(x)` is not easy. Say it involves running an
iterative algorithm. You could try automatically differentiating the
iterative algorithm, but even if that is possible, it can become quite
expensive.
In some cases we get lucky, and computing the inverse of :math:`f(x)`
is an easy operation. In these cases, we can use the `Inverse Function
Theorem <http://en.wikipedia.org/wiki/Inverse_function_theorem>`_ to
compute the derivative exactly. Here is the key idea:
Assuming that :math:`y=f(x)` is continuously differentiable in a
neighborhood of a point :math:`x` and :math:`Df(x)` is the invertible
Jacobian of :math:`f` at :math:`x`, then by applying the chain rule to
the identity :math:`f^{-1}(f(x)) = x`, we have
:math:`Df^{-1}(f(x))Df(x) = I`, or :math:`Df^{-1}(y) = (Df(x))^{-1}`,
i.e., the Jacobian of :math:`f^{-1}` is the inverse of the Jacobian of
:math:`f`, or :math:`Df(x) = (Df^{-1}(y))^{-1}`.
For example, let :math:`f(x) = e^x`. Now of course we know that
:math:`Df(x) = e^x`, but let's try and compute it via the Inverse
Function Theorem. For :math:`x > 0`, we have :math:`f^{-1}(y) = \log
y`, so :math:`Df^{-1}(y) = \frac{1}{y}`, so :math:`Df(x) =
(Df^{-1}(y))^{-1} = y = e^x`.
You maybe wondering why the above is true. A smoothly differentiable
function in a small neighborhood is well approximated by a linear
function. Indeed this is a good way to think about the Jacobian, it is
the matrix that best approximates the function linearly. Once you do
that, it is straightforward to see that *locally* :math:`f^{-1}(y)` is
best approximated linearly by the inverse of the Jacobian of
:math:`f(x)`.
Let us now consider a more practical example.
Geodetic Coordinate System Conversion
-------------------------------------
When working with data related to the Earth, one can use two different
coordinate systems. The familiar (latitude, longitude, height)
Latitude-Longitude-Altitude coordinate system or the `ECEF
<http://en.wikipedia.org/wiki/ECEF>`_ coordinate systems. The former
is familiar but is not terribly convenient analytically. The latter is
a Cartesian system but not particularly intuitive. So systems that
process earth related data have to go back and forth between these
coordinate systems.
The conversion between the LLA and the ECEF coordinate system requires
a model of the Earth, the most commonly used one being `WGS84
<https://en.wikipedia.org/wiki/World_Geodetic_System#1984_version>`_.
Going from the spherical :math:`(\phi,\lambda,h)` to the ECEF
:math:`(x,y,z)` coordinates is easy.
.. math::
\chi &= \sqrt{1 - e^2 \sin^2 \phi}
X &= \left( \frac{a}{\chi} + h \right) \cos \phi \cos \lambda
Y &= \left( \frac{a}{\chi} + h \right) \cos \phi \sin \lambda
Z &= \left(\frac{a(1-e^2)}{\chi} +h \right) \sin \phi
Here :math:`a` and :math:`e^2` are constants defined by `WGS84
<https://en.wikipedia.org/wiki/World_Geodetic_System#1984_version>`_.
Going from ECEF to LLA coordinates requires an iterative algorithm. So
to compute the derivative of the this transformation we invoke the
Inverse Function Theorem as follows:
.. code-block:: c++
Eigen::Vector3d ecef; // Fill some values
// Iterative computation.
Eigen::Vector3d lla = ECEFToLLA(ecef);
// Analytic derivatives
Eigen::Matrix3d lla_to_ecef_jacobian = LLAToECEFJacobian(lla);
bool invertible;
Eigen::Matrix3d ecef_to_lla_jacobian;
lla_to_ecef_jacobian.computeInverseWithCheck(ecef_to_lla_jacobian, invertible);
Implicit Function Theorem
=========================
Consider now the problem where we have two variables :math:`x \in
\mathbb{R}^m` and :math:`y \in \mathbb{R}^n` and a function
:math:`F:\mathbb{R}^m \times \mathbb{R}^n \rightarrow \mathbb{R}^n`
such that :math:`F(x,y) = 0` and we wish to calculate the Jacobian of
:math:`y` with respect to `x`. How do we do this?
If for a given value of :math:`(x,y)`, the partial Jacobian
:math:`D_2F(x,y)` is full rank, then the `Implicit Function Theorem
<https://en.wikipedia.org/wiki/Implicit_function_theorem>`_ tells us
that there exists a neighborhood of :math:`x` and a function :math:`G`
such :math:`y = G(x)` in this neighborhood. Differentiating
:math:`F(x,G(x)) = 0` gives us
.. math::
D_1F(x,y) + D_2F(x,y)DG(x) &= 0
DG(x) &= -(D_2F(x,y))^{-1} D_1 F(x,y)
D y(x) &= -(D_2F(x,y))^{-1} D_1 F(x,y)
This means that we can compute the derivative of :math:`y` with
respect to :math:`x` by multiplying the Jacobian of :math:`F` w.r.t
:math:`x` by the inverse of the Jacobian of :math:`F` w.r.t :math:`y`.
Let's consider two examples.
Roots of a Polynomial
---------------------
The first example we consider is a classic. Let :math:`p(x) = a_0 +
a_1 x + \dots + a_n x^n` be a degree :math:`n` polynomial, and we wish
to compute the derivative of its roots with respect to its
coefficients. There is no closed form formula for computing the roots
of a general degree :math:`n` polynomial. `Galois
<https://en.wikipedia.org/wiki/%C3%89variste_Galois>`_ and `Abel
<https://en.wikipedia.org/wiki/Niels_Henrik_Abel>`_ proved that. There
are numerical algorithms like computing the eigenvalues of the
`Companion Matrix
<https://nhigham.com/2021/03/23/what-is-a-companion-matrix/>`_, but
differentiating an eigenvalue solver does not seem like fun. But the
Implicit Function Theorem offers us a simple path.
If :math:`x` is a root of :math:`p(x)`, then :math:`F(\mathbf{a}, x) =
a_0 + a_1 x + \dots + a_n x^n = 0`. So,
.. math::
D_1 F(\mathbf{a}, x) &= [1, x, x^2, \dots, x^n]
D_2 F(\mathbf{a}, x) &= \sum_{k=1}^n k a_k x^{k-1} = Dp(x)
Dx(a) &= \frac{-1}{Dp(x)} [1, x, x^2, \dots, x^n]
Differentiating the Solution to an Optimization Problem
-------------------------------------------------------
Sometimes we are required to solve optimization problems inside
optimization problems, and this requires computing the derivative of
the optimal solution (or a fixed point) of an optimization problem
w.r.t its parameters.
Let :math:`\theta \in \mathbb{R}^m` be a vector, :math:`A(\theta) \in
\mathbb{R}^{k\times n}` be a matrix whose entries are a function of
:math:`\theta` with :math:`k \ge n` and let :math:`b \in \mathbb{R}^k`
be a constant vector, then consider the linear least squares problem:
.. math::
x^* = \arg \min_x \|A(\theta) x - b\|_2^2
How do we compute :math:`D_\theta x^*(\theta)`?
One approach would be to observe that :math:`x^*(\theta) =
(A^\top(\theta)A(\theta))^{-1}A^\top(\theta)b` and then differentiate
this w.r.t :math:`\theta`. But this would require differentiating
through the inverse of the matrix
:math:`(A^\top(\theta)A(\theta))^{-1}`. Not exactly easy. Let's use
the Implicit Function Theorem instead.
The first step is to observe that :math:`x^*` satisfies the so called
*normal equations*.
.. math::
A^\top(\theta)A(\theta)x^* - A^\top(\theta)b = 0
We will compute :math:`D_\theta x^*` column-wise, treating
:math:`A(\theta)` as a function of one coordinate (:math:`\theta_i`)
of :math:`\theta` at a time. So using the normal equations, let's
define :math:`F(\theta_i, x^*) = A^\top(\theta_i)A(\theta_i)x^* -
A^\top(\theta_i)b = 0`. Using which can now compute:
.. math::
D_1F(\theta_i, x^*) &= D_{\theta_i}A^\top A + A^\top
D_{\theta_i}Ax^* - D_{\theta_i} A^\top b = g_i
D_2F(\theta_i, x^*) &= A^\top A
Dx^*(\theta_i) & = -(A^\top A)^{-1} g_i
Dx^*(\theta) & = -(A^\top A )^{-1} \left[g_1, \dots, g_m\right]
Observe that we only need to compute the inverse of :math:`A^\top A`,
to compute :math:`D x^*(\theta)`, which we needed anyways to compute
:math:`x^*`.
docs/source/modeling_faqs.rst
-46
View File
@@ -86,49 +86,3 @@ Modeling
#. How do I set one or more components of a parameter block constant?
Using :class:`SubsetParameterization`.
#. Putting `Inverse Function Theorem
<http://en.wikipedia.org/wiki/Inverse_function_theorem>`_ to use.
Every now and then we have to deal with functions which cannot be
evaluated analytically. Computing the Jacobian in such cases is
tricky. A particularly interesting case is where the inverse of the
function is easy to compute analytically. An example of such a
function is the Coordinate transformation between the `ECEF
<http://en.wikipedia.org/wiki/ECEF>`_ and the `WGS84
<http://en.wikipedia.org/wiki/World_Geodetic_System>`_ where the
conversion from WGS84 to ECEF is analytic, but the conversion
back to WGS84 uses an iterative algorithm. So how do you compute the
derivative of the ECEF to WGS84 transformation?
One obvious approach would be to numerically
differentiate the conversion function. This is not a good idea. For
one, it will be slow, but it will also be numerically quite
bad.
Turns out you can use the `Inverse Function Theorem
<http://en.wikipedia.org/wiki/Inverse_function_theorem>`_ in this
case to compute the derivatives more or less analytically.
The key result here is. If :math:`x = f^{-1}(y)`, and :math:`Df(x)`
is the invertible Jacobian of :math:`f` at :math:`x`. Then the
Jacobian :math:`Df^{-1}(y) = [Df(x)]^{-1}`, i.e., the Jacobian of
the :math:`f^{-1}` is the inverse of the Jacobian of :math:`f`.
Algorithmically this means that given :math:`y`, compute :math:`x =
f^{-1}(y)` by whatever means you can. Evaluate the Jacobian of
:math:`f` at :math:`x`. If the Jacobian matrix is invertible, then
its inverse is the Jacobian of :math:`f^{-1}(y)` at :math:`y`.
One can put this into practice with the following code fragment.
.. code-block:: c++
Eigen::Vector3d ecef; // Fill some values
// Iterative computation.
Eigen::Vector3d lla = ECEFToLLA(ecef);
// Analytic derivatives
Eigen::Matrix3d lla_to_ecef_jacobian = LLAToECEFJacobian(lla);
bool invertible;
Eigen::Matrix3d ecef_to_lla_jacobian;
lla_to_ecef_jacobian.computeInverseWithCheck(ecef_to_lla_jacobian, invertible);