Fix MathJax and reStructuredText rendering errors across Sphinx docs

Escape underscores inside \text and \textrm macros for MathJax 3 compatibility, fix directive and list indentation that caused accidental blockquotes and definition lists, add missing math role markup, and fix malformed C++ domain cross-references.

Change-Id: I4e4bc8c8ff9dedd9c4bd5b841eafb100e791ecd3
diff --git a/docs/source/automatic_derivatives.rst b/docs/source/automatic_derivatives.rst
index 10a957c..d0eee5f 100644
--- a/docs/source/automatic_derivatives.rst
+++ b/docs/source/automatic_derivatives.rst
@@ -59,7 +59,7 @@
 
 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.
+``Rat43`` using various methods.
 
 ==========================   =========
 CostFunction                 Time (ns)
diff --git a/docs/source/conf.py b/docs/source/conf.py
index dd1352f..671c2f8 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -79,6 +79,7 @@
     ('cpp:identifier', 'Ownership'),
     ('cpp:identifier', 'PreconditionerType'),
     ('cpp:identifier', 'ResidualBlockId'),
+    ('cpp:identifier', 'Solver'),
     ('cpp:identifier', 'SparseLinearAlgebraLibraryType'),
     ('cpp:identifier', 'TerminationType'),
     ('cpp:identifier', 'TrustRegionStrategyType'),
diff --git a/docs/source/gradient_solver.rst b/docs/source/gradient_solver.rst
index 7deef3c..fc89963 100644
--- a/docs/source/gradient_solver.rst
+++ b/docs/source/gradient_solver.rst
@@ -186,11 +186,11 @@
    approximation. The increase in quality is however is bounded for a
    number of reasons.
 
-     1. The method only uses secant information and not actual
-        derivatives.
+   1. The method only uses secant information and not actual
+      derivatives.
 
-     2. The Hessian approximation is constrained to be positive
-        definite.
+   2. The Hessian approximation is constrained to be positive
+      definite.
 
    So increasing this rank to a large number will cost time and space
    complexity without the corresponding increase in solution
@@ -218,22 +218,22 @@
 
    With:
 
-  .. math:: y_k = \nabla f_{k+1} - \nabla f_k
-  .. math:: s_k = x_{k+1} - x_k
+   .. math:: y_k = \nabla f_{k+1} - \nabla f_k
+   .. math:: s_k = x_{k+1} - x_k
 
-  Where :math:`f()` is the line search objective and :math:`x` the
-  vector of parameter values [NocedalWright]_.
+   Where :math:`f()` is the line search objective and :math:`x` the
+   vector of parameter values [NocedalWright]_.
 
-  It is important to note that approximate eigenvalue scaling does
-  **not** *always* improve convergence, and that it can in fact
-  *significantly* degrade performance for certain classes of problem,
-  which is why it is disabled by default.  In particular it can
-  degrade performance when the sensitivity of the problem to different
-  parameters varies significantly, as in this case a single scalar
-  factor fails to capture this variation and detrimentally downscales
-  parts of the Jacobian approximation which correspond to
-  low-sensitivity parameters. It can also reduce the robustness of the
-  solution to errors in the Jacobians.
+   It is important to note that approximate eigenvalue scaling does
+   **not** *always* improve convergence, and that it can in fact
+   *significantly* degrade performance for certain classes of problem,
+   which is why it is disabled by default.  In particular it can
+   degrade performance when the sensitivity of the problem to different
+   parameters varies significantly, as in this case a single scalar
+   factor fails to capture this variation and detrimentally downscales
+   parts of the Jacobian approximation which correspond to
+   low-sensitivity parameters. It can also reduce the robustness of the
+   solution to errors in the Jacobians.
 
 .. member:: LineSearchInterpolationType GradientProblemSolver::Options::line_search_interpolation_type
 
@@ -247,7 +247,7 @@
 
    The line search terminates if:
 
-   .. math:: \|\Delta x_k\|_\infty < \text{min_line_search_step_size}
+   .. math:: \|\Delta x_k\|_\infty < \text{min\_line\_search\_step\_size}
 
    where :math:`\|\cdot\|_\infty` refers to the max norm, and
    :math:`\Delta x_k` is the step change in the parameter values at
@@ -264,7 +264,7 @@
    value of the objective function sufficiently. More precisely, we
    are looking for a step size s.t.
 
-   .. math:: f(\text{step_size}) \le f(0) + \text{sufficient_decrease} * [f'(0) * \text{step_size}]
+   .. math:: f(\text{step\_size}) \le f(0) + \text{sufficient\_decrease} * [f'(0) * \text{step\_size}]
 
    This condition is known as the Armijo condition.
 
@@ -274,11 +274,11 @@
 
    In each iteration of the line search,
 
-   .. math:: \text{new_step_size} \geq \text{max_line_search_step_contraction} * \text{step_size}
+   .. math:: \text{new\_step\_size} \geq \text{max\_line\_search\_step\_contraction} * \text{step\_size}
 
    Note that by definition, for contraction:
 
-   .. math:: 0 < \text{max_step_contraction} < \text{min_step_contraction} < 1
+   .. math:: 0 < \text{max\_step\_contraction} < \text{min\_step\_contraction} < 1
 
 .. member:: double GradientProblemSolver::Options::min_line_search_step_contraction
 
@@ -286,11 +286,11 @@
 
    In each iteration of the line search,
 
-   .. math:: \text{new_step_size} \leq \text{min_line_search_step_contraction} * \text{step_size}
+   .. math:: \text{new\_step\_size} \leq \text{min\_line\_search\_step\_contraction} * \text{step\_size}
 
    Note that by definition, for contraction:
 
-   .. math:: 0 < \text{max_step_contraction} < \text{min_step_contraction} < 1
+   .. math:: 0 < \text{max\_step\_contraction} < \text{min\_step\_contraction} < 1
 
 .. member:: int GradientProblemSolver::Options::max_num_line_search_step_size_iterations
 
@@ -332,7 +332,7 @@
    decreases sufficiently. Precisely, this second condition
    is that we seek a step size s.t.
 
-   .. math:: \|f'(\text{step_size})\| \leq \text{sufficient_curvature_decrease} * \|f'(0)\|
+   .. math:: \|f'(\text{step\_size})\| \leq \text{sufficient\_curvature\_decrease} * \|f'(0)\|
 
    Where :math:`f()` is the line search objective and :math:`f'()` is the derivative
    of :math:`f` with respect to the step size: :math:`\frac{d f}{d~\text{step size}}`.
@@ -347,11 +347,11 @@
    satisfying the conditions is found.  Precisely, at each iteration
    of the expansion:
 
-   .. math:: \text{new_step_size} \leq \text{max_step_expansion} * \text{step_size}
+   .. math:: \text{new\_step\_size} \leq \text{max\_step\_expansion} * \text{step\_size}
 
    By definition for expansion
 
-   .. math:: \text{max_step_expansion} > 1.0
+   .. math:: \text{max\_step\_expansion} > 1.0
 
 .. member:: int GradientProblemSolver::Options::max_num_iterations
 
@@ -370,7 +370,7 @@
 
    Solver terminates if
 
-   .. math:: \frac{|\Delta \text{cost}|}{\text{cost}} \leq \text{function_tolerance}
+   .. math:: \frac{|\Delta \text{cost}|}{\text{cost}} \leq \text{function\_tolerance}
 
    where, :math:`\Delta \text{cost}` is the change in objective
    function value (up or down) in the current iteration of the line search.
@@ -381,7 +381,7 @@
 
    Solver terminates if
 
-   .. math:: \|x - \Pi \boxplus(x, -g(x))\|_\infty \leq \text{gradient_tolerance}
+   .. math:: \|x - \Pi \boxplus(x, -g(x))\|_\infty \leq \text{gradient\_tolerance}
 
    where :math:`\|\cdot\|_\infty` refers to the max norm, :math:`\Pi`
    is projection onto the bounds constraints and :math:`\boxplus` is
@@ -394,7 +394,7 @@
 
    Solver terminates if
 
-   .. math:: \|\Delta x\| \leq (\|x\| + \text{parameter_tolerance}) * \text{parameter_tolerance}
+   .. math:: \|\Delta x\| \leq (\|x\| + \text{parameter\_tolerance}) * \text{parameter\_tolerance}
 
    where :math:`\Delta x` is the step computed by the linear solver in
    the current iteration of the line search.
@@ -477,8 +477,8 @@
 
    Whether the solution returned by the optimization algorithm can be
    relied on to be numerically sane. This will be the case if
-   `GradientProblemSolver::Summary:termination_type` is set to `CONVERGENCE`,
-   `USER_SUCCESS` or `NO_CONVERGENCE`, i.e., either the solver
+   :member:`GradientProblemSolver::Summary::termination_type` is set to ``CONVERGENCE``,
+   ``USER_SUCCESS`` or ``NO_CONVERGENCE``, i.e., either the solver
    converged by meeting one of the convergence tolerances or because
    the user indicated that it had converged or it ran to the maximum
    number of iterations or time.
@@ -550,11 +550,11 @@
 
 .. member:: NonlinearConjugateGradientType GradientProblemSolver::Summary::nonlinear_conjugate_gradient_type
 
-   If the line search direction is `NONLINEAR_CONJUGATE_GRADIENT`,
+   If the line search direction is ``NONLINEAR_CONJUGATE_GRADIENT``,
    then this indicates the particular variant of non-linear conjugate
    gradient used.
 
 .. member:: int GradientProblemSolver::Summary::max_lbfgs_rank
 
-   If the type of the line search direction is `LBFGS`, then this
+   If the type of the line search direction is ``LBFGS``, then this
    indicates the rank of the Hessian approximation.
diff --git a/docs/source/gradient_tutorial.rst b/docs/source/gradient_tutorial.rst
index 66174f0..0d9c823 100644
--- a/docs/source/gradient_tutorial.rst
+++ b/docs/source/gradient_tutorial.rst
@@ -169,7 +169,7 @@
 And finally, if you would rather compute the derivatives by hand (say
 because the size of the parameter vector is too large to be
 automatically differentiated). Then you should define an instance of
-`FirstOrderFunction`, which is the analog of :class:`CostFunction` for
+:class:`FirstOrderFunction`, which is the analog of :class:`CostFunction` for
 non-linear least squares problems [#f3]_.
 
 .. code::
diff --git a/docs/source/installation.rst b/docs/source/installation.rst
index 3e0b438..84c61f9 100644
--- a/docs/source/installation.rst
+++ b/docs/source/installation.rst
@@ -9,10 +9,10 @@
 Dependencies
 ============
 
- .. note ::
+.. note::
 
-    Starting with version 2.2, Ceres Solver requires a **fully
-    C++17-compliant** compiler.
+   Starting with version 2.2, Ceres Solver requires a **fully
+   C++17-compliant** compiler.
 
 Ceres relies on a number of open source libraries, some of which are
 optional. For details on customizing the build process, see
@@ -338,12 +338,12 @@
    <https://github.com/microsoft/vcpkg#quick-start-windows>`_, e.g., using
    Visual Studio 2022 community edition, or simply run
 
-    .. code:: bat
+   .. code:: bat
 
-        git clone https://github.com/Microsoft/vcpkg.git
-        cd vcpkg
-        .\bootstrap-vcpkg.bat
-        .\vcpkg integrate install
+      git clone https://github.com/Microsoft/vcpkg.git
+      cd vcpkg
+      .\bootstrap-vcpkg.bat
+      .\vcpkg integrate install
 
 #. Use vcpkg to install and build Ceres and all its dependencies, e.g., for 64
    bit Windows
@@ -856,13 +856,13 @@
 ``Ceres_DIR`` should be the path to the exported Ceres build
 directory.
 
-  .. NOTE ::
+.. note::
 
-     You do not need to call include_directories(${CERES_INCLUDE_DIRS})
-     as the exported Ceres CMake target already contains the definitions
-     of its public include directories which will be automatically
-     included by CMake when compiling a target that links against Ceres.
-     In fact, since v2.0 ``CERES_INCLUDE_DIRS`` is not even set.
+   You do not need to call ``include_directories(${CERES_INCLUDE_DIRS})``
+   as the exported Ceres CMake target already contains the definitions
+   of its public include directories which will be automatically
+   included by CMake when compiling a target that links against Ceres.
+   In fact, since v2.0 ``CERES_INCLUDE_DIRS`` is not even set.
 
 Specify Ceres components
 -------------------------------------
@@ -961,30 +961,30 @@
 
 - The `official CMake tutorial <https://cmake.org/cmake/help/latest/guide/tutorial/index.html>`_
 
-   Provides a tour of the core features of CMake.
+  Provides a tour of the core features of CMake.
 
 - `cmake-packages documentation
   <https://cmake.org/cmake/help/latest/manual/cmake-packages.7.html>`_
 
-   Covers how to write a ``ProjectConfig.cmake`` file, discussed below,
-   for your own project when installing or exporting it using CMake,
-   and how these processes in conjunction with ``find_package()`` are
-   handled by CMake.
+  Covers how to write a ``ProjectConfig.cmake`` file, discussed below,
+  for your own project when installing or exporting it using CMake,
+  and how these processes in conjunction with ``find_package()`` are
+  handled by CMake.
 
-  .. NOTE :: **Targets in CMake.**
+.. note:: **Targets in CMake.**
 
-    All libraries and executables built using CMake are represented as
-    *targets* created using `add_library()
-    <https://cmake.org/cmake/help/latest/command/add_library.html>`_
-    and `add_executable()
-    <https://cmake.org/cmake/help/latest/command/add_executable.html>`_.
-    Targets encapsulate the rules and dependencies (which can be other
-    targets) required to build or link against an object.  This allows
-    CMake to implicitly manage dependency chains.  Thus it is
-    sufficient to tell CMake that a library target: ``B`` depends on a
-    previously declared library target ``A``, and CMake will
-    understand that this means that ``B`` also depends on all of the
-    public dependencies of ``A``.
+   All libraries and executables built using CMake are represented as
+   *targets* created using `add_library()
+   <https://cmake.org/cmake/help/latest/command/add_library.html>`_
+   and `add_executable()
+   <https://cmake.org/cmake/help/latest/command/add_executable.html>`_.
+   Targets encapsulate the rules and dependencies (which can be other
+   targets) required to build or link against an object.  This allows
+   CMake to implicitly manage dependency chains.  Thus it is
+   sufficient to tell CMake that a library target: ``B`` depends on a
+   previously declared library target ``A``, and CMake will
+   understand that this means that ``B`` also depends on all of the
+   public dependencies of ``A``.
 
 When a project like Ceres is installed using CMake, or its build
 directory is exported (see :ref:`section-install-vs-export`), in addition to
@@ -1002,19 +1002,19 @@
 - ``<PROJECT_NAME>Config.cmake`` (or
   ``<lower_case_project_name>-config.cmake``)
 
-   Which is written by the developers of the project, and is
-   configured with the selected options and installed locations when
-   the project is built and imports the project targets and/or defines
-   the legacy CMake variables: ``<PROJECT_NAME>_INCLUDE_DIRS`` &
-   ``<PROJECT_NAME>_LIBRARIES`` which are used by the caller.
+  Which is written by the developers of the project, and is
+  configured with the selected options and installed locations when
+  the project is built and imports the project targets and/or defines
+  the legacy CMake variables: ``<PROJECT_NAME>_INCLUDE_DIRS`` &
+  ``<PROJECT_NAME>_LIBRARIES`` which are used by the caller.
 
 The ``<PROJECT_NAME>Config.cmake`` typically includes a second file
 installed to the same location:
 
 - ``<PROJECT_NAME>Targets.cmake``
 
-   Which is autogenerated by CMake as part of the install process and defines
-   **imported targets** for the project in the caller's CMake scope.
+  Which is autogenerated by CMake as part of the install process and defines
+  **imported targets** for the project in the caller's CMake scope.
 
 An **imported target** contains the same information about a library
 as a CMake target that was declared locally in the current CMake
@@ -1126,7 +1126,7 @@
   ``git clone --recurse-submodules`` so that CMake automatically builds the
   bundled submodules in ``third_party/``.
 - **CMake minimum version**: Building Ceres now requires CMake 3.22 or later.
-- **CMake feature options use the ``WITH_`` prefix**:
+- **CMake feature options use the** ``WITH_`` **prefix**:
   The Ceres-specific CMake options have been standardized to use the ``WITH_``
   prefix:
 
@@ -1144,7 +1144,7 @@
     ``WITH_STRIPPED_DEBUG_SYMBOLS``
   - ``PROVIDE_UNINSTALL_TARGET`` :math:`\rightarrow` ``WITH_UNINSTALL_TARGET``
 
-- **SuiteSparse is now opt-in (``WITH_SUITESPARSE=OFF`` by default)**:
+- **SuiteSparse is now opt-in** (``WITH_SUITESPARSE=OFF`` **by default**):
   Because ``SuiteSparse``'s supernodal ``CHOLMOD`` and ``SPQR`` components are
   licensed under the GPL, ``WITH_SUITESPARSE`` now defaults to ``OFF``. Pass
   ``-DWITH_SUITESPARSE=ON`` to CMake to build Ceres with ``SuiteSparse``
diff --git a/docs/source/inverse_and_implicit_function_theorems.rst b/docs/source/inverse_and_implicit_function_theorems.rst
index a9173a6..f281eb0 100644
--- a/docs/source/inverse_and_implicit_function_theorems.rst
+++ b/docs/source/inverse_and_implicit_function_theorems.rst
@@ -109,7 +109,7 @@
 \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?
+:math:`y` with respect to :math:`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
diff --git a/docs/source/nnls_covariance.rst b/docs/source/nnls_covariance.rst
index 8b82887..b8ce5af 100644
--- a/docs/source/nnls_covariance.rst
+++ b/docs/source/nnls_covariance.rst
@@ -112,23 +112,24 @@
 numerical values is more complicated. Here again there are two
 cases.
 
-  a. The rank deficiency arises from overparameterization. e.g., a
-     four dimensional quaternion used to parameterize :math:`SO(3)`,
-     which is a three dimensional manifold. In cases like this, the
-     user should use an appropriate
-     :class:`Manifold`. Not only will this lead to better
-     numerical behaviour of the Solver, it will also expose the rank
-     deficiency to the :class:`Covariance` object so that it can
-     handle it correctly.
+a. The rank deficiency arises from overparameterization. e.g., a
+   four dimensional quaternion used to parameterize :math:`SO(3)`,
+   which is a three dimensional manifold. In cases like this, the
+   user should use an appropriate
+   :class:`Manifold`. Not only will this lead to better
+   numerical behaviour of the Solver, it will also expose the rank
+   deficiency to the :class:`Covariance` object so that it can
+   handle it correctly.
 
-  b. More general numerical rank deficiency in the Jacobian requires
-     the computation of the so called Singular Value Decomposition
-     (SVD) of :math:`J'J`. We do not know how to do this for large
-     sparse matrices efficiently. For small and moderate sized
-     problems this is done using dense linear algebra.
+b. More general numerical rank deficiency in the Jacobian requires
+   the computation of the so called Singular Value Decomposition
+   (SVD) of :math:`J'J`. We do not know how to do this for large
+   sparse matrices efficiently. For small and moderate sized
+   problems this is done using dense linear algebra.
 
 
 :class:`Covariance::Options`
+============================
 
 .. class:: Covariance::Options
 
@@ -157,10 +158,10 @@
    1. ``SPARSE_QR`` uses the sparse QR factorization algorithm to
       compute the decomposition
 
-       .. math::
+      .. math::
 
-          QR &= J\\
-          \left(J^\top J\right)^{-1} &= \left(R^\top R\right)^{-1}
+         QR &= J\\
+         \left(J^\top J\right)^{-1} &= \left(R^\top R\right)^{-1}
 
       The speed of this algorithm depends on the sparse linear algebra
       library being used. ``Eigen``'s sparse QR factorization is a
@@ -192,17 +193,17 @@
 
    Default: :math:`-1`
 
-    During QR factorization, if a column with Euclidean norm less than
-    ``column_pivot_threshold`` is encountered it is treated as zero.
+   During QR factorization, if a column with Euclidean norm less than
+   ``column_pivot_threshold`` is encountered it is treated as zero.
 
-    If ``column_pivot_threshold < 0``, then an automatic default value
-    of `20*(m+n)*eps*sqrt(max(diag(J’*J)))` is used.  Here `m` and `n`
-    are the number of rows and columns of the Jacobian (`J`)
-    respectively.
+   If ``column_pivot_threshold < 0``, then an automatic default value
+   of ``20*(m+n)*eps*sqrt(max(diag(J'*J)))`` is used.  Here :math:`m` and :math:`n`
+   are the number of rows and columns of the Jacobian (:math:`J`)
+   respectively.
 
-    This is an advanced option meant for users who know enough about
-    their Jacobian matrices that they can determine a value better
-    than the default.
+   This is an advanced option meant for users who know enough about
+   their Jacobian matrices that they can determine a value better
+   than the default.
 
 
 .. member:: int Covariance::Options::min_reciprocal_condition_number
@@ -236,7 +237,7 @@
 
    1. ``DENSE_SVD``
 
-      .. math:: \frac{\sigma_{\text{min}}}{\sigma_{\text{max}}}  < \sqrt{\text{min_reciprocal_condition_number}}
+      .. math:: \frac{\sigma_{\text{min}}}{\sigma_{\text{max}}}  < \sqrt{\text{min\_reciprocal\_condition\_number}}
 
       where :math:`\sigma_{\text{min}}` and
       :math:`\sigma_{\text{max}}` are the minimum and maximum
@@ -244,51 +245,51 @@
 
    2. ``SPARSE_QR``
 
-       .. math:: \operatorname{rank}(J) < \operatorname{num\_col}(J)
+      .. math:: \operatorname{rank}(J) < \operatorname{num\_col}(J)
 
-       Here :math:`\operatorname{rank}(J)` is the estimate of the rank
-       of :math:`J` returned by the sparse QR factorization
-       algorithm. It is a fairly reliable indication of rank
-       deficiency.
+      Here :math:`\operatorname{rank}(J)` is the estimate of the rank
+      of :math:`J` returned by the sparse QR factorization
+      algorithm. It is a fairly reliable indication of rank
+      deficiency.
 
 .. member:: int Covariance::Options::null_space_rank
 
-    When using ``DENSE_SVD``, the user has more control in dealing
-    with singular and near singular covariance matrices.
+   When using ``DENSE_SVD``, the user has more control in dealing
+   with singular and near singular covariance matrices.
 
-    As mentioned above, when the covariance matrix is near singular,
-    instead of computing the inverse of :math:`J'J`, the Moore-Penrose
-    pseudoinverse of :math:`J'J` should be computed.
+   As mentioned above, when the covariance matrix is near singular,
+   instead of computing the inverse of :math:`J'J`, the Moore-Penrose
+   pseudoinverse of :math:`J'J` should be computed.
 
-    If :math:`J'J` has the eigen decomposition :math:`(\lambda_i,
-    e_i)`, where :math:`\lambda_i` is the :math:`i^\textrm{th}`
-    eigenvalue and :math:`e_i` is the corresponding eigenvector, then
-    the inverse of :math:`J'J` is
+   If :math:`J'J` has the eigen decomposition :math:`(\lambda_i,
+   e_i)`, where :math:`\lambda_i` is the :math:`i^\textrm{th}`
+   eigenvalue and :math:`e_i` is the corresponding eigenvector, then
+   the inverse of :math:`J'J` is
 
-    .. math:: (J'J)^{-1} = \sum_i \frac{1}{\lambda_i} e_i e_i'
+   .. math:: (J'J)^{-1} = \sum_i \frac{1}{\lambda_i} e_i e_i'
 
-    and computing the pseudo inverse involves dropping terms from this
-    sum that correspond to small eigenvalues.
+   and computing the pseudo inverse involves dropping terms from this
+   sum that correspond to small eigenvalues.
 
-    How terms are dropped is controlled by
-    `min_reciprocal_condition_number` and `null_space_rank`.
+   How terms are dropped is controlled by
+   ``min_reciprocal_condition_number`` and ``null_space_rank``.
 
-    If `null_space_rank` is non-negative, then the smallest
-    `null_space_rank` eigenvalue/eigenvectors are dropped irrespective
-    of the magnitude of :math:`\lambda_i`. If the ratio of the
-    smallest non-zero eigenvalue to the largest eigenvalue in the
-    truncated matrix is still below min_reciprocal_condition_number,
-    then the `Covariance::Compute()` will fail and return `false`.
+   If ``null_space_rank`` is non-negative, then the smallest
+   ``null_space_rank`` eigenvalue/eigenvectors are dropped irrespective
+   of the magnitude of :math:`\lambda_i`. If the ratio of the
+   smallest non-zero eigenvalue to the largest eigenvalue in the
+   truncated matrix is still below ``min_reciprocal_condition_number``,
+   then :func:`Covariance::Compute` will fail and return ``false``.
 
-    Setting `null_space_rank = -1` drops all terms for which
+   Setting ``null_space_rank = -1`` drops all terms for which
 
-    .. math::  \frac{\lambda_i}{\lambda_{\textrm{max}}} < \textrm{min_reciprocal_condition_number}
+   .. math::  \frac{\lambda_i}{\lambda_{\textrm{max}}} < \textrm{min\_reciprocal\_condition\_number}
 
-    This option has no effect on ``SPARSE_QR``.
+   This option has no effect on ``SPARSE_QR``.
 
 .. member:: bool Covariance::Options::apply_loss_function
 
-   Default: `true`
+   Default: ``true``
 
    Even though the residual blocks in the problem may contain loss
    functions, setting ``apply_loss_function`` to false will turn off
@@ -329,7 +330,7 @@
    :class:`Covariance::Options` for more on the conditions under which
    this function returns ``false``.
 
-.. function:: bool GetCovarianceBlock(const double* parameter_block1, const double* parameter_block2, double* covariance_block) const
+.. function:: bool Covariance::GetCovarianceBlock(const double* parameter_block1, const double* parameter_block2, double* covariance_block) const
 
    Return the block of the cross-covariance matrix corresponding to
    ``parameter_block1`` and ``parameter_block2``.
@@ -344,12 +345,12 @@
    a ``parameter_block1_size x parameter_block2_size`` matrix. The
    returned covariance will be a row-major matrix.
 
-.. function:: bool GetCovarianceBlockInTangentSpace(const double* parameter_block1, const double* parameter_block2, double* covariance_block) const
+.. function:: bool Covariance::GetCovarianceBlockInTangentSpace(const double* parameter_block1, const double* parameter_block2, double* covariance_block) const
 
    Return the block of the cross-covariance matrix corresponding to
    ``parameter_block1`` and ``parameter_block2``.
-   Returns cross-covariance in the tangent space if a local
-   parameterization is associated with either parameter block;
+   Returns cross-covariance in the tangent space if a :class:`Manifold`
+   is associated with either parameter block;
    else returns cross-covariance in the ambient space.
 
    Compute must be called before the first call to ``GetCovarianceBlock``
diff --git a/docs/source/nnls_modeling.rst b/docs/source/nnls_modeling.rst
index 4b5a8ee..150ec13 100644
--- a/docs/source/nnls_modeling.rst
+++ b/docs/source/nnls_modeling.rst
@@ -269,7 +269,7 @@
    2-dimensional.
 
    By default :class:`AutoDiffCostFunction` will take ownership of the cost
-   functor pointer passed to it, ie. will call `delete` on the cost functor
+   functor pointer passed to it, ie. will call ``delete`` on the cost functor
    when the :class:`AutoDiffCostFunction` itself is deleted. However, this may
    be undesirable in certain cases, therefore it is also possible to specify
    ``DO_NOT_TAKE_OWNERSHIP`` as a second argument in the constructor,
@@ -320,7 +320,7 @@
    applications, this is not enough e.g., Bezier curve fitting, Neural
    Network training etc.
 
-     .. code-block:: c++
+   .. code-block:: c++
 
       template <typename CostFunctor, int Stride = 4>
       class DynamicAutoDiffCostFunction : public CostFunction {
@@ -331,24 +331,24 @@
    templated functor, but the signature of the functor differs
    slightly. The expected interface for the cost functors is:
 
-     .. code-block:: c++
+   .. code-block:: c++
 
-       struct MyCostFunctor {
-         template<typename T>
-         bool operator()(T const* const* parameters, T* residuals) const {
-         }
-       }
+      struct MyCostFunctor {
+        template<typename T>
+        bool operator()(T const* const* parameters, T* residuals) const {
+        }
+      }
 
    Since the sizing of the parameters is done at runtime, you must
    also specify the sizes after creating the dynamic autodiff cost
    function. For example:
 
-     .. code-block:: c++
+   .. code-block:: c++
 
-       auto* cost_function = new DynamicAutoDiffCostFunction<MyCostFunctor, 4>();
-       cost_function->AddParameterBlock(5);
-       cost_function->AddParameterBlock(10);
-       cost_function->SetNumResiduals(21);
+      auto* cost_function = new DynamicAutoDiffCostFunction<MyCostFunctor, 4>();
+      cost_function->AddParameterBlock(5);
+      cost_function->AddParameterBlock(10);
+      cost_function->SetNumResiduals(21);
 
    Under the hood, the implementation evaluates the cost function
    multiple times, computing a small set of the derivatives (four by
@@ -460,7 +460,7 @@
                                   Dimension of y -------------------------+
 
   In this example, there is usually an instance for each measurement
-  of `k`.
+  of ``k``.
 
   In the instantiation above, the template parameters following
   ``MyScalarCostFunctor``, ``1, 2, 2``, describe the functor as
@@ -470,7 +470,7 @@
   NumericDiffCostFunction also supports cost functions with a
   runtime-determined number of residuals. For example:
 
-   .. code-block:: c++
+  .. code-block:: c++
 
      auto functor = std::make_unique<CostFunctorWithDynamicNumResiduals>(1.0);
      auto* cost_function
@@ -498,7 +498,7 @@
   estimating :math:`f'(x)` by computing
   :math:`\frac{f(x+h)-f(x-h)}{2h}`.
 
-  The ``RIDDERS`` difference method[Ridders]_ is an adaptive scheme
+  The ``RIDDERS`` difference method [Ridders]_ is an adaptive scheme
   that estimates derivatives by performing multiple central
   differences at varying scales. Specifically, the algorithm starts at
   a certain :math:`h` and as the derivative is estimated, this step
@@ -512,6 +512,15 @@
   results, either try forward difference to improve performance or
   Ridders' method to improve accuracy.
 
+  .. warning::
+      A common beginner's error when first using
+      :class:`NumericDiffCostFunction` is to get the sizing wrong. In
+      particular, there is a tendency to set the template parameters to
+      (dimension of residual, number of parameters) instead of passing a
+      dimension parameter for *every parameter*. In the example above, that
+      would be ``<MyScalarCostFunctor, 1, 2>``, which is missing the last ``2``
+      argument. Please be careful when setting the size parameters.
+
 .. class:: NumericDiffOptions
 
    Options pertaining to numeric differentiation (e.g., step sizes and
@@ -553,77 +562,68 @@
    The factor by which to shrink the step size with each extrapolation in
    Ridders' method.
 
-  .. warning::
-      A common beginner's error when first using
-      :class:`NumericDiffCostFunction` is to get the sizing wrong. In
-      particular, there is a tendency to set the template parameters to
-      (dimension of residual, number of parameters) instead of passing a
-      dimension parameter for *every parameter*. In the example above, that
-      would be ``<MyScalarCostFunctor, 1, 2>``, which is missing the last ``2``
-      argument. Please be careful when setting the size parameters.
-
 
 Numeric Differentiation & Manifolds
 -----------------------------------
 
-   If your cost function depends on a parameter block that must lie on
-   a manifold and the functor cannot be evaluated for values of that
-   parameter block not on the manifold then you may have problems
-   numerically differentiating such functors.
+If your cost function depends on a parameter block that must lie on
+a manifold and the functor cannot be evaluated for values of that
+parameter block not on the manifold then you may have problems
+numerically differentiating such functors.
 
-   This is because numeric differentiation in Ceres is performed by
-   perturbing the individual coordinates of the parameter blocks that
-   a cost functor depends on. This perturbation assumes that the
-   parameter block lives on a Euclidean Manifold rather than the
-   actual manifold associated with the parameter block. As a result
-   some of the perturbed points may not lie on the manifold anymore.
+This is because numeric differentiation in Ceres is performed by
+perturbing the individual coordinates of the parameter blocks that
+a cost functor depends on. This perturbation assumes that the
+parameter block lives on a Euclidean Manifold rather than the
+actual manifold associated with the parameter block. As a result
+some of the perturbed points may not lie on the manifold anymore.
 
-   For example consider a four dimensional parameter block that is
-   interpreted as a unit Quaternion. Perturbing the coordinates of
-   this parameter block will violate the unit norm property of the
-   parameter block.
+For example consider a four dimensional parameter block that is
+interpreted as a unit Quaternion. Perturbing the coordinates of
+this parameter block will violate the unit norm property of the
+parameter block.
 
-   Fixing this problem requires that :class:`NumericDiffCostFunction`
-   be aware of the :class:`Manifold` associated with each
-   parameter block and only generate perturbations in the local
-   tangent space of each parameter block.
+Fixing this problem requires that :class:`NumericDiffCostFunction`
+be aware of the :class:`Manifold` associated with each
+parameter block and only generate perturbations in the local
+tangent space of each parameter block.
 
-   For now this is not considered to be a serious enough problem to
-   warrant changing the :class:`NumericDiffCostFunction` API. Further,
-   in most cases it is relatively straightforward to project a point
-   off the manifold back onto the manifold before using it in the
-   functor. For example in case of the Quaternion, normalizing the
-   4-vector before using it does the trick.
+For now this is not considered to be a serious enough problem to
+warrant changing the :class:`NumericDiffCostFunction` API. Further,
+in most cases it is relatively straightforward to project a point
+off the manifold back onto the manifold before using it in the
+functor. For example in case of the Quaternion, normalizing the
+4-vector before using it does the trick.
 
-   **Alternate Interface**
+**Alternate Interface**
 
-   For a variety of reasons, including compatibility with legacy code,
-   :class:`NumericDiffCostFunction` can also take
-   :class:`CostFunction` objects as input. The following describes
-   how.
+For a variety of reasons, including compatibility with legacy code,
+:class:`NumericDiffCostFunction` can also take
+:class:`CostFunction` objects as input. The following describes
+how.
 
-   To get a numerically differentiated cost function, define a
-   subclass of :class:`CostFunction` such that the
-   :func:`CostFunction::Evaluate` function ignores the ``jacobians``
-   parameter. The numeric differentiation wrapper will fill in the
-   jacobian parameter if necessary by repeatedly calling the
-   :func:`CostFunction::Evaluate` with small changes to the
-   appropriate parameters, and computing the slope. For performance,
-   the numeric differentiation wrapper class is templated on the
-   concrete cost function, even though it could be implemented only in
-   terms of the :class:`CostFunction` interface.
+To get a numerically differentiated cost function, define a
+subclass of :class:`CostFunction` such that the
+:func:`CostFunction::Evaluate` function ignores the ``jacobians``
+parameter. The numeric differentiation wrapper will fill in the
+jacobian parameter if necessary by repeatedly calling the
+:func:`CostFunction::Evaluate` with small changes to the
+appropriate parameters, and computing the slope. For performance,
+the numeric differentiation wrapper class is templated on the
+concrete cost function, even though it could be implemented only in
+terms of the :class:`CostFunction` interface.
 
-   The numerically differentiated version of a cost function for a
-   cost function can be constructed as follows:
+The numerically differentiated version of a cost function for a
+cost function can be constructed as follows:
 
-   .. code-block:: c++
+.. code-block:: c++
 
-     auto* cost_function
-         = new NumericDiffCostFunction<MyCostFunction, CENTRAL, 1, 4, 8>(...);
+   auto* cost_function
+       = new NumericDiffCostFunction<MyCostFunction, CENTRAL, 1, 4, 8>(...);
 
-   where ``MyCostFunction`` has 1 residual and 2 parameter blocks with
-   sizes 4 and 8 respectively. Look at the tests for a more detailed
-   example.
+where ``MyCostFunction`` has 1 residual and 2 parameter blocks with
+sizes 4 and 8 respectively. Look at the tests for a more detailed
+example.
 
 :class:`DynamicNumericDiffCostFunction`
 =======================================
@@ -634,7 +634,7 @@
    requires that the number of parameter blocks and their sizes be
    known at compile time. In a number of applications, this is not enough.
 
-     .. code-block:: c++
+   .. code-block:: c++
 
       template <typename CostFunctor, NumericDiffMethodType method = CENTRAL>
       class DynamicNumericDiffCostFunction : public CostFunction {
@@ -647,23 +647,23 @@
    functor, but the signature of the functor differs slightly. The
    expected interface for the cost functors is:
 
-     .. code-block:: c++
+   .. code-block:: c++
 
-       struct MyCostFunctor {
-         bool operator()(double const* const* parameters, double* residuals) const {
-         }
-       }
+      struct MyCostFunctor {
+        bool operator()(double const* const* parameters, double* residuals) const {
+        }
+      }
 
    Since the sizing of the parameters is done at runtime, you must
    also specify the sizes after creating the dynamic numeric diff cost
    function. For example:
 
-     .. code-block:: c++
+   .. code-block:: c++
 
-       auto cost_function = std::make_unique<DynamicNumericDiffCostFunction<MyCostFunctor>>();
-       cost_function->AddParameterBlock(5);
-       cost_function->AddParameterBlock(10);
-       cost_function->SetNumResiduals(21);
+      auto cost_function = std::make_unique<DynamicNumericDiffCostFunction<MyCostFunctor>>();
+      cost_function->AddParameterBlock(5);
+      cost_function->AddParameterBlock(10);
+      cost_function->SetNumResiduals(21);
 
    As a rule of thumb, try using :class:`NumericDiffCostFunction` before
    you use :class:`DynamicNumericDiffCostFunction`.
@@ -770,34 +770,34 @@
     };
 
 
-  Here ``ThirdPartyProjectionFunction`` is some third party library
-  function that we have no control over. So this function can compute
-  its value and we would like to use numeric differentiation to
-  compute its derivatives. In this case we can use a combination of
-  ``NumericDiffCostFunction`` and ``CostFunctionToFunctor`` to get the
-  job done.
+   Here ``ThirdPartyProjectionFunction`` is some third party library
+   function that we have no control over. So this function can compute
+   its value and we would like to use numeric differentiation to
+   compute its derivatives. In this case we can use a combination of
+   ``NumericDiffCostFunction`` and ``CostFunctionToFunctor`` to get the
+   job done.
 
-  .. code-block:: c++
+   .. code-block:: c++
 
-   struct CameraProjection {
-     explicit CameraProjection(double* observation)
-        : intrinsic_projection_(
-              std::make_unique<NumericDiffCostFunction<IntrinsicProjection, CENTRAL, 2, 5, 3>>()) {}
+    struct CameraProjection {
+      explicit CameraProjection(double* observation)
+         : intrinsic_projection_(
+               std::make_unique<NumericDiffCostFunction<IntrinsicProjection, CENTRAL, 2, 5, 3>>()) {}
 
-     template <typename T>
-     bool operator()(const T* rotation,
-                     const T* translation,
-                     const T* intrinsics,
-                     const T* point,
-                     T* residuals) const {
-       T transformed_point[3];
-       RotateAndTranslatePoint(rotation, translation, point, transformed_point);
-       return intrinsic_projection_(intrinsics, transformed_point, residuals);
-     }
+      template <typename T>
+      bool operator()(const T* rotation,
+                      const T* translation,
+                      const T* intrinsics,
+                      const T* point,
+                      T* residuals) const {
+        T transformed_point[3];
+        RotateAndTranslatePoint(rotation, translation, point, transformed_point);
+        return intrinsic_projection_(intrinsics, transformed_point, residuals);
+      }
 
-    private:
-     CostFunctionToFunctor<2, 5, 3> intrinsic_projection_;
-   };
+     private:
+      CostFunctionToFunctor<2, 5, 3> intrinsic_projection_;
+    };
 
 
 :class:`DynamicCostFunctionToFunctor`
@@ -915,20 +915,20 @@
 
 .. class:: GradientChecker
 
-    This class compares the Jacobians returned by a cost function
-    against derivatives estimated using finite differencing. It is
-    meant as a tool for unit testing, giving you more fine-grained
-    control than the check_gradients option in the solver options.
+   This class compares the Jacobians returned by a cost function
+   against derivatives estimated using finite differencing. It is
+   meant as a tool for unit testing, giving you more fine-grained
+   control than the check_gradients option in the solver options.
 
-    The condition enforced is that
+   The condition enforced is that
 
-    .. math:: \forall{i,j}: \frac{J_{ij} - J'_{ij}}{max_{ij}(J_{ij} - J'_{ij})} < r
+   .. math:: \forall{i,j}: \frac{J_{ij} - J'_{ij}}{max_{ij}(J_{ij} - J'_{ij})} < r
 
-    where :math:`J_{ij}` is the jacobian as computed by the supplied
-    cost function multiplied by the `Manifold::PlusJacobian`,
-    :math:`J'_{ij}` is the jacobian as computed by finite differences,
-    multiplied by the `Manifold::PlusJacobian` as well, and :math:`r`
-    is the relative precision.
+   where :math:`J_{ij}` is the jacobian as computed by the supplied
+   cost function multiplied by the :func:`Manifold::PlusJacobian`,
+   :math:`J'_{ij}` is the jacobian as computed by finite differences,
+   multiplied by the :func:`Manifold::PlusJacobian` as well, and :math:`r`
+   is the relative precision.
 
    Usage:
 
@@ -988,14 +988,14 @@
    is the variable. In case the user is interested in implementing a cost
    function of the form
 
-  .. math::  cost(x) = (x - \mu)^T S^{-1} (x - \mu)
+   .. math::  cost(x) = (x - \mu)^T S^{-1} (x - \mu)
 
-  where, :math:`\mu` is a vector and :math:`S` is a covariance matrix,
-  then, :math:`A = S^{-1/2}`, i.e the matrix :math:`A` is the square
-  root of the inverse of the covariance, also known as the stiffness
-  matrix. There are however no restrictions on the shape of
-  :math:`A`. It is free to be rectangular, which would be the case if
-  the covariance matrix :math:`S` is rank deficient.
+   where, :math:`\mu` is a vector and :math:`S` is a covariance matrix,
+   then, :math:`A = S^{-1/2}`, i.e the matrix :math:`A` is the square
+   root of the inverse of the covariance, also known as the stiffness
+   matrix. There are however no restrictions on the shape of
+   :math:`A`. It is free to be rectangular, which would be the case if
+   the covariance matrix :math:`S` is rank deficient.
 
 
 
@@ -1532,7 +1532,7 @@
    EuclideanManifold<3> manifold;
 
 If the ambient space dimensions is not known at compile time the
-template parameter needs to be set to `ceres::DYNAMIC` and the actual
+template parameter needs to be set to ``ceres::DYNAMIC`` and the actual
 dimension needs to be provided as a constructor argument:
 
 .. code-block:: c++
@@ -1712,7 +1712,7 @@
    SphereManifold<3> manifold;
 
 If the ambient space dimensions is not known at compile time the
-template parameter needs to be set to `ceres::DYNAMIC` and the actual
+template parameter needs to be set to ``ceres::DYNAMIC`` and the actual
 dimension needs to be provided as a constructor argument:
 
 .. code-block:: c++
@@ -2010,7 +2010,7 @@
 
 .. member:: bool Problem::Options::disable_all_safety_checks
 
-    Default: `false`
+    Default: ``false``
 
     By default, Ceres performs a variety of safety checks when
     constructing the problem. There is a small but measurable
@@ -2030,7 +2030,7 @@
 
     A Ceres global context to use for solving this problem. This may
     help to reduce computation time as Ceres can reuse expensive
-    objects to create.  The context object can be `nullptr`, in which
+    objects to create.  The context object can be ``nullptr``, in which
     case Ceres may create one.
 
     Ceres does NOT take ownership of the pointer.
@@ -2125,9 +2125,9 @@
    (unless :member:`Problem::Options::disable_all_safety_checks` is set to true).
 
    Repeated calls with the same double pointer and size but different
-   :class:`Manifold` is equivalent to calling `SetManifold(manifold)`,
+   :class:`Manifold` is equivalent to calling ``SetManifold(manifold)``,
    i.e., any previously associated :class:`Manifold` object will be replaced
-   with the `manifold`.
+   with the ``manifold``.
 
 .. function:: void Problem::AddParameterBlock(double* values, int size)
 
@@ -2223,27 +2223,27 @@
 
 .. function:: void Problem::SetParameterLowerBound(double* values, int index, double lower_bound)
 
-   Set the lower bound for the parameter at position `index` in the
-   parameter block corresponding to `values`. By default the lower
+   Set the lower bound for the parameter at position ``index`` in the
+   parameter block corresponding to ``values``. By default the lower
    bound is ``-std::numeric_limits<double>::max()``, which is treated
    by the solver as the same as :math:`-\infty`.
 
 .. function:: void Problem::SetParameterUpperBound(double* values, int index, double upper_bound)
 
-   Set the upper bound for the parameter at position `index` in the
-   parameter block corresponding to `values`. By default the value is
+   Set the upper bound for the parameter at position ``index`` in the
+   parameter block corresponding to ``values``. By default the value is
    ``std::numeric_limits<double>::max()``, which is treated by the
    solver as the same as :math:`\infty`.
 
 .. function:: double Problem::GetParameterLowerBound(const double* values, int index)
 
-   Get the lower bound for the parameter with position `index`. If the
+   Get the lower bound for the parameter with position ``index``. If the
    parameter is not bounded by the user, then its lower bound is
    ``-std::numeric_limits<double>::max()``.
 
 .. function:: double Problem::GetParameterUpperBound(const double* values, int index)
 
-   Get the upper bound for the parameter with position `index`. If the
+   Get the upper bound for the parameter with position ``index``. If the
    parameter is not bounded by the user, then its upper bound is
    ``std::numeric_limits<double>::max()``.
 
@@ -2289,9 +2289,9 @@
 
 .. function:: void Problem::GetResidualBlocks(std::vector<ResidualBlockId>* residual_blocks) const
 
-   Fills the passed `residual_blocks` vector with pointers to the
+   Fills the passed ``residual_blocks`` vector with pointers to the
    residual blocks currently in the problem. After this call,
-   `residual_blocks.size() == NumResidualBlocks`.
+   ``residual_blocks.size() == NumResidualBlocks``.
 
 .. function:: void Problem::GetParameterBlocksForResidualBlock(const ResidualBlockId residual_block, std::vector<double*>* parameter_blocks) const
 
@@ -2346,8 +2346,8 @@
    .. NOTE:: If an :class:`EvaluationCallback` is associated with the
       problem, then its
       :func:`EvaluationCallback::PrepareForEvaluation` method will be
-      called every time this method is called with `new_point =
-      true`. This conservatively assumes that the user may have
+      called every time this method is called with ``new_point =
+      true``. This conservatively assumes that the user may have
       changed the parameter values since the previous call to evaluate
       / solve.  For improved efficiency, and only if you know that the
       parameter values have not changed between calls, see
@@ -2396,8 +2396,8 @@
         double cost = 0.0;
         problem.Evaluate(Problem::EvaluateOptions(), &cost, nullptr, nullptr, nullptr);
 
-      The cost is evaluated at `x = 1`. If you wish to evaluate the
-      problem at `x = 2`, then
+      The cost is evaluated at ``x = 1``. If you wish to evaluate the
+      problem at ``x = 2``, then
 
       .. code-block:: c++
 
@@ -2651,11 +2651,11 @@
 
 .. function:: template <typename T> void CrossProduct(const T x[3], const T y[3], T x_cross_y[3])
 
-   .. math:: \text{x_cross_y} = x \times y
+   .. math:: \text{x\_cross\_y} = x \times y
 
 .. function:: template <typename T> void AngleAxisRotatePoint(const T angle_axis[3], const T pt[3], T result[3])
 
-   .. math:: y = R(\text{angle_axis}) x
+   .. math:: y = R(\text{angle\_axis}) x
 
 
 Cubic Interpolation
diff --git a/docs/source/nnls_solving.rst b/docs/source/nnls_solving.rst
index 5403f80..feb6bfb 100644
--- a/docs/source/nnls_solving.rst
+++ b/docs/source/nnls_solving.rst
@@ -86,21 +86,21 @@
 
 The basic trust region algorithm looks something like this.
 
-   1. Given an initial point :math:`x` and a trust region radius :math:`\mu`.
-   2. Solve
+1. Given an initial point :math:`x` and a trust region radius :math:`\mu`.
+2. Solve
 
-      .. math::
-         \arg \min_{\Delta x}& \frac{1}{2}\|J(x)\Delta x + F(x)\|^2 \\
-         \text{such that} &\|D(x)\Delta x\|^2 \le \mu\\
-         &L \le x + \Delta x \le U.
+   .. math::
+      \arg \min_{\Delta x}& \frac{1}{2}\|J(x)\Delta x + F(x)\|^2 \\
+      \text{such that} &\|D(x)\Delta x\|^2 \le \mu\\
+      &L \le x + \Delta x \le U.
 
-   3. :math:`\rho = \frac{\displaystyle \|F(x + \Delta x)\|^2 -
-      \|F(x)\|^2}{\displaystyle \|J(x)\Delta x + F(x)\|^2 -
-      \|F(x)\|^2}`
-   4. if :math:`\rho > \epsilon` then  :math:`x = x + \Delta x`.
-   5. if :math:`\rho > \eta_1` then :math:`\mu = 2  \mu`
-   6. else if :math:`\rho < \eta_2` then :math:`\mu = 0.5 * \mu`
-   7. Go to 2.
+3. :math:`\rho = \frac{\displaystyle \|F(x + \Delta x)\|^2 -
+   \|F(x)\|^2}{\displaystyle \|J(x)\Delta x + F(x)\|^2 -
+   \|F(x)\|^2}`
+4. if :math:`\rho > \epsilon` then  :math:`x = x + \Delta x`.
+5. if :math:`\rho > \eta_1` then :math:`\mu = 2  \mu`
+6. else if :math:`\rho < \eta_2` then :math:`\mu = 0.5 * \mu`
+7. Go to 2.
 
 Here, :math:`\mu` is the trust region radius, :math:`D(x)` is some
 matrix used to define a metric on the domain of :math:`F(x)` and
@@ -322,8 +322,8 @@
 .. math:: y = f_1(a_1, e^{b_1 x}) + f_2(a_2, e^{b_3 x^2 + c_1})
 
 In this case, we solve for the trust region step for the full problem,
-and then use it as the starting point to further optimize just `a_1`
-and `a_2`. For the linear case, this amounts to doing a single linear
+and then use it as the starting point to further optimize just :math:`a_1`
+and :math:`a_2`. For the linear case, this amounts to doing a single linear
 least squares solve. For non-linear problems, any method for solving
 the :math:`a_1` and :math:`a_2` optimization problems will do. The
 only constraint on :math:`a_1` and :math:`a_2` (if they are two
@@ -389,18 +389,18 @@
 
 The basic line search algorithm looks something like this:
 
-   1. Given an initial point :math:`x`
-   2. :math:`\Delta x = -H^{-1}(x) g(x)`
-   3. :math:`\arg \min_\mu \frac{1}{2} \| F(x + \mu \Delta x) \|^2`
-   4. :math:`x = x + \mu \Delta x`
-   5. Goto 2.
+1. Given an initial point :math:`x`
+2. :math:`\Delta x = -H^{-1}(x) g(x)`
+3. :math:`\arg \min_\mu \frac{1}{2} \| F(x + \mu \Delta x) \|^2`
+4. :math:`x = x + \mu \Delta x`
+5. Goto 2.
 
 Here :math:`H(x)` is some approximation to the Hessian of the
 objective function, and :math:`g(x)` is the gradient at
 :math:`x`. Depending on the choice of :math:`H(x)` we get a variety of
 different search directions :math:`\Delta x`.
 
-Step 4, which is a one dimensional optimization or `Line Search` along
+Step 4, which is a one dimensional optimization or *Line Search* along
 :math:`\Delta x` is what gives this class of methods its name.
 
 Different line search algorithms differ in their choice of the search
@@ -429,7 +429,7 @@
    are currently the best known general quasi-Newton algorithm.
 
 4. ``LBFGS`` A limited memory approximation to the full ``BFGS``
-   method in which the last `M` iterations are used to approximate the
+   method in which the last :math:`M` iterations are used to approximate the
    inverse Hessian used to compute a quasi-Newton step [Nocedal]_,
    [ByrdNocedal]_.
 
@@ -493,7 +493,7 @@
 implementation, then the user can also choose to use LAPACK's
 `DGEQRF`_ routine by setting
 :member:`Solver::Options::dense_linear_algebra_library_type` to
-``LAPACK``. Depending on the `LAPACK` and the underlying `BLAS`
+``LAPACK``. Depending on the ``LAPACK`` and the underlying ``BLAS``
 implementation this may perform better than using Eigen's Householder
 QR factorization.
 
@@ -560,7 +560,7 @@
 implementation, then the user can also choose to use LAPACK's
 `DPOTRF`_ routine by setting
 :member:`Solver::Options::dense_linear_algebra_library_type` to
-``LAPACK``. Depending on the `LAPACK` and the underlying `BLAS`
+``LAPACK``. Depending on the ``LAPACK`` and the underlying ``BLAS``
 implementation this may perform better than using Eigen's Cholesky
 factorization.
 
@@ -599,16 +599,16 @@
 The use of this linear solver requires that Ceres is compiled with
 support for at least one of:
 
- 1. `SuiteSparse <https://people.engr.tamu.edu/davis/suitesparse.html>`_ (``SUITE_SPARSE``).
- 2. `Apple's Accelerate framework
-    <https://developer.apple.com/documentation/accelerate/sparse-solvers-library?language=objc>`_
-    (``ACCELERATE_SPARSE``).
- 3. `Eigen's sparse linear solvers
-    <https://libeigen.gitlab.io/eigen/docs-nightly/group__SparseCholesky__Module.html>`_
-    (``EIGEN_SPARSE``).
- 4. `NVIDIA's cuDSS library
-    <https://developer.nvidia.com/cudss>`_
-    (``CUDA_SPARSE``).
+1. `SuiteSparse <https://people.engr.tamu.edu/davis/suitesparse.html>`_ (``SUITE_SPARSE``).
+2. `Apple's Accelerate framework
+   <https://developer.apple.com/documentation/accelerate/sparse-solvers-library?language=objc>`_
+   (``ACCELERATE_SPARSE``).
+3. `Eigen's sparse linear solvers
+   <https://libeigen.gitlab.io/eigen/docs-nightly/group__SparseCholesky__Module.html>`_
+   (``EIGEN_SPARSE``).
+4. `NVIDIA's cuDSS library
+   <https://developer.nvidia.com/cudss>`_
+   (``CUDA_SPARSE``).
 
 SuiteSparse, Accelerate, and cuDSS offer high performance sparse Cholesky
 factorization routines (using level-3 BLAS routines on the CPU or GPU
@@ -636,7 +636,7 @@
 Cholesky factorization or a sparse linear algebra library is not
 linked into Ceres, another option is the ``CGNR`` solver. This solver
 uses the `Conjugate Gradients
-<https://en.wikipedia.org/wiki/Conjugate_gradient_method>_` method on
+<https://en.wikipedia.org/wiki/Conjugate_gradient_method>`_ method on
 the *normal equations*, but without forming the normal equations
 explicitly. It exploits the relation
 
@@ -681,7 +681,7 @@
 ... ,z_{q}]`. Where, :math:`y` and :math:`z` correspond to camera and
 point parameters respectively.  Further, let the camera blocks be of
 size :math:`c` and the point blocks be of size :math:`s` (for most
-problems :math:`c` = :math:`6`--`9` and :math:`s = 3`). Ceres does not
+problems :math:`c` = :math:`6`--:math:`9` and :math:`s = 3`). Ceres does not
 impose any constancy requirement on these block sizes, but choosing
 them to be constant simplifies the exposition.
 
@@ -1004,7 +1004,7 @@
 structure of the scene. The similarity between a pair of cameras
 :math:`i` and :math:`j` is given by:
 
-  .. math:: S_{ij} = \frac{|V_i \cap V_j|}{|V_i| |V_j|}
+.. math:: S_{ij} = \frac{|V_i \cap V_j|}{|V_i| |V_j|}
 
 Here :math:`V_i` is the set of scene points visible in camera
 :math:`i`. This idea was first exploited by [KushalAgarwal]_ to create
@@ -1045,7 +1045,7 @@
 
 Suppose the Jacobian :math:`J` has been horizontally partitioned as
 
-  .. math:: J = \begin{bmatrix} P \\ Q \end{bmatrix}
+.. math:: J = \begin{bmatrix} P \\ Q \end{bmatrix}
 
 Where, :math:`Q` is the set of rows corresponding to the residual
 blocks in
@@ -1287,22 +1287,22 @@
 
    With:
 
-  .. math:: y_k = \nabla f_{k+1} - \nabla f_k
-  .. math:: s_k = x_{k+1} - x_k
+   .. math:: y_k = \nabla f_{k+1} - \nabla f_k
+   .. math:: s_k = x_{k+1} - x_k
 
-  Where :math:`f()` is the line search objective and :math:`x` the
-  vector of parameter values [NocedalWright]_.
+   Where :math:`f()` is the line search objective and :math:`x` the
+   vector of parameter values [NocedalWright]_.
 
-  It is important to note that approximate eigenvalue scaling does
-  **not** *always* improve convergence, and that it can in fact
-  *significantly* degrade performance for certain classes of problem,
-  which is why it is disabled by default.  In particular it can
-  degrade performance when the sensitivity of the problem to different
-  parameters varies significantly, as in this case a single scalar
-  factor fails to capture this variation and detrimentally downscales
-  parts of the Jacobian approximation which correspond to
-  low-sensitivity parameters. It can also reduce the robustness of the
-  solution to errors in the Jacobians.
+   It is important to note that approximate eigenvalue scaling does
+   **not** *always* improve convergence, and that it can in fact
+   *significantly* degrade performance for certain classes of problem,
+   which is why it is disabled by default.  In particular it can
+   degrade performance when the sensitivity of the problem to different
+   parameters varies significantly, as in this case a single scalar
+   factor fails to capture this variation and detrimentally downscales
+   parts of the Jacobian approximation which correspond to
+   low-sensitivity parameters. It can also reduce the robustness of the
+   solution to errors in the Jacobians.
 
 .. member:: LineSearchInterpolationType Solver::Options::line_search_interpolation_type
 
@@ -1318,7 +1318,7 @@
 
    The line search terminates if:
 
-   .. math:: \|\Delta x_k\|_\infty < \text{min_line_search_step_size}
+   .. math:: \|\Delta x_k\|_\infty < \text{min\_line\_search\_step\_size}
 
    where :math:`\|\cdot\|_\infty` refers to the max norm, and
    :math:`\Delta x_k` is the step change in the parameter values at
@@ -1335,7 +1335,7 @@
    value of the objective function sufficiently. More precisely, we
    are looking for a step size s.t.
 
-   .. math:: f(\text{step_size}) \le f(0) + \text{sufficient_decrease} * [f'(0) * \text{step_size}]
+   .. math:: f(\text{step\_size}) \le f(0) + \text{sufficient\_decrease} * [f'(0) * \text{step\_size}]
 
    This condition is known as the Armijo condition.
 
@@ -1345,11 +1345,11 @@
 
    In each iteration of the line search,
 
-   .. math:: \text{new_step_size} >= \text{max_line_search_step_contraction} * \text{step_size}
+   .. math:: \text{new\_step\_size} \ge \text{max\_line\_search\_step\_contraction} * \text{step\_size}
 
    Note that by definition, for contraction:
 
-   .. math:: 0 < \text{max_step_contraction} < \text{min_step_contraction} < 1
+   .. math:: 0 < \text{max\_step\_contraction} < \text{min\_step\_contraction} < 1
 
 .. member:: double Solver::Options::min_line_search_step_contraction
 
@@ -1357,11 +1357,11 @@
 
    In each iteration of the line search,
 
-   .. math:: \text{new_step_size} <= \text{min_line_search_step_contraction} * \text{step_size}
+   .. math:: \text{new\_step\_size} \le \text{min\_line\_search\_step\_contraction} * \text{step\_size}
 
    Note that by definition, for contraction:
 
-   .. math:: 0 < \text{max_step_contraction} < \text{min_step_contraction} < 1
+   .. math:: 0 < \text{max\_step\_contraction} < \text{min\_step\_contraction} < 1
 
 .. member:: int Solver::Options::max_num_line_search_step_size_iterations
 
@@ -1379,7 +1379,7 @@
    As this is an 'artificial' constraint (one imposed by the user, not
    the underlying math), if ``WOLFE`` line search is being used, *and*
    points satisfying the Armijo sufficient (function) decrease
-   condition have been found during the current search (in :math:`<=`
+   condition have been found during the current search (in :math:`\le`
    ``max_num_line_search_step_size_iterations``).  Then, the step size
    with the lowest function value which satisfies the Armijo condition
    will be returned as the new valid step, even though it does *not*
@@ -1408,7 +1408,7 @@
    decreases sufficiently. Precisely, this second condition
    is that we seek a step size s.t.
 
-   .. math:: \|f'(\text{step_size})\| <= \text{sufficient_curvature_decrease} * \|f'(0)\|
+   .. math:: \|f'(\text{step\_size})\| \le \text{sufficient\_curvature\_decrease} * \|f'(0)\|
 
    Where :math:`f()` is the line search objective and :math:`f'()` is the derivative
    of :math:`f` with respect to the step size: :math:`\frac{d f}{d~\text{step size}}`.
@@ -1423,11 +1423,11 @@
    satisfying the conditions is found.  Precisely, at each iteration
    of the expansion:
 
-   .. math:: \text{new_step_size} <= \text{max_step_expansion} * \text{step_size}
+   .. math:: \text{new\_step\_size} \le \text{max\_step\_expansion} * \text{step\_size}
 
    By definition for expansion
 
-   .. math:: \text{max_step_expansion} > 1.0
+   .. math:: \text{max\_step\_expansion} > 1.0
 
 .. member:: TrustRegionStrategyType Solver::Options::trust_region_strategy_type
 
@@ -1541,7 +1541,7 @@
 
    Solver terminates if
 
-   .. math:: \frac{|\Delta \text{cost}|}{\text{cost}} <= \text{function_tolerance}
+   .. math:: \frac{|\Delta \text{cost}|}{\text{cost}} \le \text{function\_tolerance}
 
    where, :math:`\Delta \text{cost}` is the change in objective
    function value (up or down) in the current iteration of
@@ -1553,7 +1553,7 @@
 
    Solver terminates if
 
-   .. math:: \|x - \Pi \boxplus(x, -g(x))\|_\infty <= \text{gradient_tolerance}
+   .. math:: \|x - \Pi \boxplus(x, -g(x))\|_\infty \le \text{gradient\_tolerance}
 
    where :math:`\|\cdot\|_\infty` refers to the max norm, :math:`\Pi`
    is projection onto the bounds constraints and :math:`\boxplus` is
@@ -1566,7 +1566,7 @@
 
    Solver terminates if
 
-   .. math:: \|\Delta x\| <= (\|x\| + \text{parameter_tolerance}) * \text{parameter_tolerance}
+   .. math:: \|\Delta x\| \le (\|x\| + \text{parameter\_tolerance}) * \text{parameter\_tolerance}
 
    where :math:`\Delta x` is the step computed by the linear solver in
    the current iteration.
@@ -1624,7 +1624,7 @@
 
    Suppose the Jacobian :math:`J` has been horizontally partitioned as
 
-       .. math:: J = \begin{bmatrix} P \\ Q \end{bmatrix}
+   .. math:: J = \begin{bmatrix} P \\ Q \end{bmatrix}
 
    Where, :math:`Q` is the set of rows corresponding to the residual
    blocks in
@@ -1687,22 +1687,22 @@
 
    Default: ``AMD``
 
-    The order in which variables are eliminated in a linear solver can
-    have a significant impact on the efficiency and accuracy of the
-    method. e.g., when doing sparse Cholesky factorization, there are
-    matrices for which a good ordering will give a Cholesky factor
-    with :math:`O(n)` storage, where as a bad ordering will result in
-    an completely dense factor.
+   The order in which variables are eliminated in a linear solver can
+   have a significant impact on the efficiency and accuracy of the
+   method. e.g., when doing sparse Cholesky factorization, there are
+   matrices for which a good ordering will give a Cholesky factor
+   with :math:`O(n)` storage, where as a bad ordering will result in
+   an completely dense factor.
 
-    Sparse direct solvers like ``SPARSE_NORMAL_CHOLESKY`` and
-    ``SPARSE_SCHUR`` use a fill reducing ordering of the columns and
-    rows of the matrix being factorized before computing the numeric
-    factorization.
+   Sparse direct solvers like ``SPARSE_NORMAL_CHOLESKY`` and
+   ``SPARSE_SCHUR`` use a fill reducing ordering of the columns and
+   rows of the matrix being factorized before computing the numeric
+   factorization.
 
-    This enum controls the type of algorithm used to compute this fill
-    reducing ordering. There is no single algorithm that works on all
-    matrices, so determining which algorithm works better is a matter
-    of empirical experimentation.
+   This enum controls the type of algorithm used to compute this fill
+   reducing ordering. There is no single algorithm that works on all
+   matrices, so determining which algorithm works better is a matter
+   of empirical experimentation.
 
 .. member:: std::shared_ptr<ParameterBlockOrdering> Solver::Options::linear_solver_ordering
 
@@ -1758,7 +1758,7 @@
    then it is probably best to keep this false, otherwise it will
    likely lead to worse performance.
 
-   This setting only affects the `SPARSE_NORMAL_CHOLESKY` solver.
+   This setting only affects the ``SPARSE_NORMAL_CHOLESKY`` solver.
 
 .. member:: bool Solver::Options::use_mixed_precision_solves
 
@@ -1801,7 +1801,7 @@
 
 .. member:: int Solver::Options::max_num_spse_iterations
 
-   Default: `5`
+   Default: ``5``
 
    Maximum number of iterations performed by
    ``SCHUR_POWER_SERIES_EXPANSION``. Each iteration corresponds to one
@@ -1810,7 +1810,7 @@
    whether it is used as a preconditioner or just to initialize the
    solution for ``ITERATIVE_SCHUR``.
 
-.. member:: bool Solver:Options::use_spse_initialization
+.. member:: bool Solver::Options::use_spse_initialization
 
    Default: ``false``
 
@@ -1820,7 +1820,7 @@
 
 .. member:: double Solver::Options::spse_tolerance
 
-   Default: `0.1`
+   Default: ``0.1``
 
    When ``use_spse_initialization`` is ``true``, this parameter along
    with ``max_num_spse_iterations`` controls the number of
@@ -1968,12 +1968,12 @@
 
    Default: ``/tmp``
 
-    Directory to which the problems should be written to. Should be
-    non-empty if
-    :member:`Solver::Options::trust_region_minimizer_iterations_to_dump` is
-    non-empty and
-    :member:`Solver::Options::trust_region_problem_dump_format_type` is not
-    ``CONSOLE``.
+   Directory to which the problems should be written to. Should be
+   non-empty if
+   :member:`Solver::Options::trust_region_minimizer_iterations_to_dump` is
+   non-empty and
+   :member:`Solver::Options::trust_region_problem_dump_format_type` is not
+   ``CONSOLE``.
 
 .. member:: DumpFormatType Solver::Options::trust_region_problem_dump_format_type
 
@@ -1984,10 +1984,10 @@
    is non-empty.  There are three options:
 
    * ``CONSOLE`` prints the linear least squares problem in a human
-      readable format to ``stderr``. The Jacobian is printed as a
-      dense matrix. The vectors :math:`D`, :math:`x` and :math:`f` are
-      printed as dense vectors. This should only be used for small
-      problems.
+     readable format to ``stderr``. The Jacobian is printed as a
+     dense matrix. The vectors :math:`D`, :math:`x` and :math:`f` are
+     printed as dense vectors. This should only be used for small
+     problems.
 
    * ``TEXTFILE`` Write out the linear least squares problem to the
      directory pointed to by
@@ -2005,7 +2005,7 @@
      used to parse and load the problem into memory.
 
      .. NOTE::
-        **Understanding the dumped ``x.txt`` vector:**
+        **Understanding the dumped** ``x.txt`` **vector:**
 
         The content of ``x.txt`` depends on the
         :member:`Solver::Options::trust_region_strategy_type` being
@@ -2014,12 +2014,12 @@
         matrix:
 
         * If the strategy is ``LEVENBERG_MARQUARDT``, ``x.txt``
-            contains the **damped trust region step** (the solution
-            to the scaled problem).
+          contains the **damped trust region step** (the solution
+          to the scaled problem).
         * If the strategy is ``DOGLEG``, ``x.txt`` contains the
-            **pure, un-damped Gauss-Newton step**. This is *not*
-            necessarily the final step taken by the optimizer, as it
-            has not yet been limited by the trust region radius.
+          **pure, un-damped Gauss-Newton step**. This is *not*
+          necessarily the final step taken by the optimizer, as it
+          has not yet been limited by the trust region radius.
 
 .. member:: bool Solver::Options::check_gradients
 
@@ -2048,7 +2048,7 @@
 
       This option only applies to the numeric differentiation used for
       checking the user provided derivatives when
-      `Solver::Options::check_gradients` is true. If you are using
+      :member:`Solver::Options::check_gradients` is ``true``. If you are using
       :class:`NumericDiffCostFunction` and are interested in changing
       the step size for numeric differentiation in your cost function,
       please have a look at :class:`NumericDiffOptions`.
@@ -2056,7 +2056,7 @@
    Relative initial step size used for taking numeric derivatives (via
    Ridders' adaptive extrapolation method; see
    :member:`NumericDiffOptions::ridders_relative_initial_step_size`)
-   when `Solver::Options::check_gradients` is `true`.
+   when :member:`Solver::Options::check_gradients` is ``true``.
 
    If parameters are close to zero, Ridders' method sets the initial
    step size directly to this value; otherwise it is multiplied by the
@@ -2183,18 +2183,18 @@
    Step was numerically valid, i.e., all values are finite and the
    step reduces the value of the linearized model.
 
-    **Note**: :member:`IterationSummary::step_is_valid` is `false`
-    when :member:`IterationSummary::iteration` = 0.
+   **Note**: :member:`IterationSummary::step_is_valid` is ``false``
+   when :member:`IterationSummary::iteration` = 0.
 
 .. member::  bool IterationSummary::step_is_nonmonotonic
 
-    Step did not reduce the value of the objective function
-    sufficiently, but it was accepted because of the relaxed
-    acceptance criterion used by the non-monotonic trust region
-    algorithm.
+   Step did not reduce the value of the objective function
+   sufficiently, but it was accepted because of the relaxed
+   acceptance criterion used by the non-monotonic trust region
+   algorithm.
 
-    **Note**: :member:`IterationSummary::step_is_nonmonotonic` is
-    `false` when when :member:`IterationSummary::iteration` = 0.
+   **Note**: :member:`IterationSummary::step_is_nonmonotonic` is
+   ``false`` when :member:`IterationSummary::iteration` = 0.
 
 .. member:: bool IterationSummary::step_is_successful
 
@@ -2204,12 +2204,12 @@
    relative reduction in the objective function value was greater than
    :member:`Solver::Options::min_relative_decrease`. However, if the
    non-monotonic trust region algorithm is used
-   (:member:`Solver::Options::use_nonmonotonic_steps` = `true`), then
+   (:member:`Solver::Options::use_nonmonotonic_steps` = ``true``), then
    even if the relative decrease is not sufficient, the algorithm may
    accept the step and the step is declared successful.
 
-   **Note**: :member:`IterationSummary::step_is_successful` is `false`
-   when when :member:`IterationSummary::iteration` = 0.
+   **Note**: :member:`IterationSummary::step_is_successful` is ``false``
+   when :member:`IterationSummary::iteration` = 0.
 
 .. member:: double IterationSummary::cost
 
@@ -2303,72 +2303,72 @@
       };
 
 
-  The solver uses the return value of ``operator()`` to decide whether
-  to continue solving or to terminate. The user can return three
-  values.
+   The solver uses the return value of ``operator()`` to decide whether
+   to continue solving or to terminate. The user can return three
+   values.
 
-  #. ``SOLVER_ABORT`` indicates that the callback detected an abnormal
-     situation. The solver returns without updating the parameter
-     blocks (unless ``Solver::Options::update_state_every_iteration`` is
-     set true). Solver returns with ``Solver::Summary::termination_type``
-     set to ``USER_FAILURE``.
+   #. ``SOLVER_ABORT`` indicates that the callback detected an abnormal
+      situation. The solver returns without updating the parameter
+      blocks (unless ``Solver::Options::update_state_every_iteration`` is
+      set true). Solver returns with ``Solver::Summary::termination_type``
+      set to ``USER_FAILURE``.
 
-  #. ``SOLVER_TERMINATE_SUCCESSFULLY`` indicates that there is no need
-     to optimize anymore (some user specified termination criterion
-     has been met). Solver returns with
-     ``Solver::Summary::termination_type``` set to ``USER_SUCCESS``.
+   #. ``SOLVER_TERMINATE_SUCCESSFULLY`` indicates that there is no need
+      to optimize anymore (some user specified termination criterion
+      has been met). Solver returns with
+      ``Solver::Summary::termination_type`` set to ``USER_SUCCESS``.
 
-  #. ``SOLVER_CONTINUE`` indicates that the solver should continue
-     optimizing.
+   #. ``SOLVER_CONTINUE`` indicates that the solver should continue
+      optimizing.
 
-  The return values can be used to implement custom termination
-  criterion that supercede the iteration/time/tolerance based
-  termination implemented by Ceres.
+   The return values can be used to implement custom termination
+   criterion that supercede the iteration/time/tolerance based
+   termination implemented by Ceres.
 
-  For example, the following :class:`IterationCallback` is used
-  internally by Ceres to log the progress of the optimization.
+   For example, the following :class:`IterationCallback` is used
+   internally by Ceres to log the progress of the optimization.
 
-  .. code-block:: c++
+   .. code-block:: c++
 
-    class LoggingCallback : public IterationCallback {
-     public:
-      explicit LoggingCallback(bool log_to_stdout)
-          : log_to_stdout_(log_to_stdout) {}
+     class LoggingCallback : public IterationCallback {
+      public:
+       explicit LoggingCallback(bool log_to_stdout)
+           : log_to_stdout_(log_to_stdout) {}
 
-      ~LoggingCallback() {}
+       ~LoggingCallback() {}
 
-      CallbackReturnType operator()(const IterationSummary& summary) {
-        const char* kReportRowFormat =
-            "% 4d: f:% 8e d:% 3.2e g:% 3.2e h:% 3.2e "
-            "rho:% 3.2e mu:% 3.2e eta:% 3.2e li:% 3d";
-        string output = StringPrintf(kReportRowFormat,
-                                     summary.iteration,
-                                     summary.cost,
-                                     summary.cost_change,
-                                     summary.gradient_max_norm,
-                                     summary.step_norm,
-                                     summary.relative_decrease,
-                                     summary.trust_region_radius,
-                                     summary.eta,
-                                     summary.linear_solver_iterations);
-        if (log_to_stdout_) {
-          cout << output << endl;
-        } else {
-          VLOG(1) << output;
-        }
-        return SOLVER_CONTINUE;
-      }
+       CallbackReturnType operator()(const IterationSummary& summary) {
+         const char* kReportRowFormat =
+             "% 4d: f:% 8e d:% 3.2e g:% 3.2e h:% 3.2e "
+             "rho:% 3.2e mu:% 3.2e eta:% 3.2e li:% 3d";
+         string output = StringPrintf(kReportRowFormat,
+                                      summary.iteration,
+                                      summary.cost,
+                                      summary.cost_change,
+                                      summary.gradient_max_norm,
+                                      summary.step_norm,
+                                      summary.relative_decrease,
+                                      summary.trust_region_radius,
+                                      summary.eta,
+                                      summary.linear_solver_iterations);
+         if (log_to_stdout_) {
+           cout << output << endl;
+         } else {
+           VLOG(1) << output;
+         }
+         return SOLVER_CONTINUE;
+       }
 
-     private:
-      const bool log_to_stdout_;
-    };
+      private:
+       const bool log_to_stdout_;
+     };
 
 
-  See `examples/iteration_callback_example.cc
-  <https://github.com/ceres-solver/ceres-solver/blob/master/examples/iteration_callback_example.cc>`_
-  for another example that uses
-  :member:`Solver::Options::update_state_every_iteration` to log
-  changes to the parameter blocks over the course of the optimization.
+   See `examples/iteration_callback_example.cc
+   <https://github.com/ceres-solver/ceres-solver/blob/master/examples/iteration_callback_example.cc>`_
+   for another example that uses
+   :member:`Solver::Options::update_state_every_iteration` to log
+   changes to the parameter blocks over the course of the optimization.
 
 
 :class:`CRSMatrix`
@@ -2449,8 +2449,8 @@
 
    Whether the solution returned by the optimization algorithm can be
    relied on to be numerically sane. This will be the case if
-   `Solver::Summary:termination_type` is set to `CONVERGENCE`,
-   `USER_SUCCESS` or `NO_CONVERGENCE`, i.e., either the solver
+   :member:`Solver::Summary::termination_type` is set to ``CONVERGENCE``,
+   ``USER_SUCCESS`` or ``NO_CONVERGENCE``, i.e., either the solver
    converged by meeting one of the convergence tolerances or because
    the user indicated that it had converged or it ran to the maximum
    number of iterations or time.
@@ -2491,7 +2491,7 @@
 
    Number of minimizer iterations in which the step was
    accepted. Unless :member:`Solver::Options::use_nonmonotonic_steps`
-   is `true` this is also the number of steps in which the objective
+   is ``true`` this is also the number of steps in which the objective
    function value/cost went down.
 
 .. member:: int Solver::Summary::num_unsuccessful_steps
@@ -2629,7 +2629,7 @@
    determines that the problem structure is not compatible with the
    linear solver requested or if the linear solver requested by the
    user is not available, e.g. The user requested
-   `SPARSE_NORMAL_CHOLESKY` but no sparse linear algebra library was
+   ``SPARSE_NORMAL_CHOLESKY`` but no sparse linear algebra library was
    available.
 
 .. member:: std::vector<int> Solver::Summary::linear_solver_ordering_given
@@ -2665,12 +2665,12 @@
 
 .. member:: bool Solver::Summary::inner_iterations_given
 
-   `True` if the user asked for inner iterations to be used as part of
+   ``true`` if the user asked for inner iterations to be used as part of
    the optimization.
 
 .. member:: bool Solver::Summary::inner_iterations_used
 
-   `True` if the user asked for inner iterations to be used as part of
+   ``true`` if the user asked for inner iterations to be used as part of
    the optimization and the problem structure was such that they were
    actually performed. For example, in a problem with just one parameter
    block, inner iterations are not performed.
@@ -2739,11 +2739,11 @@
 
 .. member:: NonlinearConjugateGradientType Solver::Summary::nonlinear_conjugate_gradient_type
 
-   If the line search direction is `NONLINEAR_CONJUGATE_GRADIENT`,
+   If the line search direction is ``NONLINEAR_CONJUGATE_GRADIENT``,
    then this indicates the particular variant of non-linear conjugate
    gradient used.
 
 .. member:: int Solver::Summary::max_lbfgs_rank
 
-   If the type of the line search direction is `LBFGS`, then this
+   If the type of the line search direction is ``LBFGS``, then this
    indicates the rank of the Hessian approximation.
diff --git a/docs/source/numerical_derivatives.rst b/docs/source/numerical_derivatives.rst
index 95b7271..3aaac21 100644
--- a/docs/source/numerical_derivatives.rst
+++ b/docs/source/numerical_derivatives.rst
@@ -36,11 +36,11 @@
 <https://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``.
+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++
 
@@ -204,22 +204,22 @@
 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`.
+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
diff --git a/docs/source/spivak_notation.rst b/docs/source/spivak_notation.rst
index 3ac56ba..4f10f1a 100644
--- a/docs/source/spivak_notation.rst
+++ b/docs/source/spivak_notation.rst
@@ -31,7 +31,7 @@
    D_1 g = \frac{\partial}{\partial x}g(x,y) \text{ and }  D_2 g  = \frac{\partial}{\partial y}g(x,y).
 
 
-:math:`Dg` denotes the Jacobian of `g`, i.e.,
+:math:`Dg` denotes the Jacobian of :math:`g`, i.e.,
 
 .. math::
 
diff --git a/docs/source/version_history.rst b/docs/source/version_history.rst
index 740d785..71207bb 100644
--- a/docs/source/version_history.rst
+++ b/docs/source/version_history.rst
@@ -403,7 +403,7 @@
 #. Fix unused arguments of ``Make1stOrderPerturbation`` (Dmitriy
    Korchemkin)
 #. Fix ``SuiteSparse`` path and version reporting (Sergiu Deitsch)
-#. Enable `GitHub` workflows and deprecate ``TravisCI`` (Sergiu
+#. Enable ``GitHub`` workflows and deprecate ``TravisCI`` (Sergiu
    Deitsch)
 #. Add missing includes (Sergiu Deitsch, Sameer Agarwal)
 #. Fix path for ``cuda-memcheck`` tests (Joydeep Biswas)
@@ -1304,7 +1304,7 @@
 #. ``Solver::Options::gradient_tolerance`` used to be a relative
    gradient tolerance. i.e., The solver converged when
 
-   .. math:: \|g(x)\|_\infty < \text{gradient_tolerance} *
+   .. math:: \|g(x)\|_\infty < \text{gradient\_tolerance} *
       \|g(x_0)\|_\infty
 
    where :math:`g(x)` is the gradient of the objective function at
@@ -1314,7 +1314,7 @@
    This has changed to an absolute tolerance, i.e. the solver
    converges when
 
-   .. math:: \|g(x)\|_\infty < \text{gradient_tolerance}
+   .. math:: \|g(x)\|_\infty < \text{gradient\_tolerance}
 
 #. Ceres cannot be built without the line search minimizer
    anymore. Thus the preprocessor define
@@ -1324,7 +1324,7 @@
 ---------
 #. Disabled warning C4251. (Björn Piltz)
 #. Do not propagate 3d party libs through
-   `IMPORTED_LINK_INTERFACE_LIBRARIES_[DEBUG/RELEASE]` mechanism when
+   ``IMPORTED_LINK_INTERFACE_LIBRARIES_[DEBUG/RELEASE]`` mechanism when
    building shared libraries. (Björn Piltz)
 #. Fixed errant verbose levels (Björn Piltz)
 #. Variety of code cleanups, optimizations and bug fixes to the line
@@ -1336,8 +1336,8 @@
 #. Fixed a memory leak in ``SchurComplementSolver``.
 #. Added ``epsilon()`` method to ``NumTraits<ceres::Jet<T, N>
    >``. (Filippo Basso)
-#. Fixed a bug in `CompressedRowSparseMatrix::AppendRows`` and
-   ``DeleteRows``.q
+#. Fixed a bug in ``CompressedRowSparseMatrix::AppendRows`` and
+   ``DeleteRows``.
 #. Handle empty problems consistently.
 #. Restore the state of the ``Problem`` after a call to
    ``Problem::Evaluate``. (Stefan Leutenegger)
@@ -1865,7 +1865,7 @@
 #. Diagonal scaling allows for equal upper and lower bounds. (Markus Moll)
 #. TrustRegionStrategy does not misuse LinearSolver:Summary anymore.
 #. Fix Eigen3 Row/Column Major storage issue. (Lena Gieseke)
-#. QuaternionToAngleAxis now guarantees an angle in $[-\pi, \pi]$. (Guoxuan Zhang)
+#. QuaternionToAngleAxis now guarantees an angle in :math:`[-\pi, \pi]`. (Guoxuan Zhang)
 #. Added a workaround for a compiler bug in the Android NDK to the
    Schur eliminator.
 #. The sparse linear algebra library is only logged in