Convert documentation from LaTeX to Sphinx.
A live version of the doc can be found at
http://homes.cs.washington.edu/~sagarwal/ceres-solver/
As it stands, the documentation has better hyperlinking
and coverage than the latex documentation now.
Change-Id: I7ede3aa83b9b9ef25104caf331e5727b4f5beae5
diff --git a/docs/Makefile b/docs/Makefile
new file mode 100644
index 0000000..fe17b50
--- /dev/null
+++ b/docs/Makefile
@@ -0,0 +1,153 @@
+# Makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS =
+SPHINXBUILD = sphinx-build
+PAPER =
+BUILDDIR = build
+
+# Internal variables.
+PAPEROPT_a4 = -D latex_paper_size=a4
+PAPEROPT_letter = -D latex_paper_size=letter
+ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
+# the i18n builder cannot share the environment and doctrees with the others
+I18NSPHINXOPTS = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
+
+.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest gettext
+
+help:
+ @echo "Please use \`make <target>' where <target> is one of"
+ @echo " html to make standalone HTML files"
+ @echo " dirhtml to make HTML files named index.html in directories"
+ @echo " singlehtml to make a single large HTML file"
+ @echo " pickle to make pickle files"
+ @echo " json to make JSON files"
+ @echo " htmlhelp to make HTML files and a HTML help project"
+ @echo " qthelp to make HTML files and a qthelp project"
+ @echo " devhelp to make HTML files and a Devhelp project"
+ @echo " epub to make an epub"
+ @echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
+ @echo " latexpdf to make LaTeX files and run them through pdflatex"
+ @echo " text to make text files"
+ @echo " man to make manual pages"
+ @echo " texinfo to make Texinfo files"
+ @echo " info to make Texinfo files and run them through makeinfo"
+ @echo " gettext to make PO message catalogs"
+ @echo " changes to make an overview of all changed/added/deprecated items"
+ @echo " linkcheck to check all external links for integrity"
+ @echo " doctest to run all doctests embedded in the documentation (if enabled)"
+
+clean:
+ -rm -rf $(BUILDDIR)/*
+
+html:
+ $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
+ @echo
+ @echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
+
+dirhtml:
+ $(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
+ @echo
+ @echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."
+
+singlehtml:
+ $(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml
+ @echo
+ @echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml."
+
+pickle:
+ $(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
+ @echo
+ @echo "Build finished; now you can process the pickle files."
+
+json:
+ $(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
+ @echo
+ @echo "Build finished; now you can process the JSON files."
+
+htmlhelp:
+ $(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
+ @echo
+ @echo "Build finished; now you can run HTML Help Workshop with the" \
+ ".hhp project file in $(BUILDDIR)/htmlhelp."
+
+qthelp:
+ $(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
+ @echo
+ @echo "Build finished; now you can run "qcollectiongenerator" with the" \
+ ".qhcp project file in $(BUILDDIR)/qthelp, like this:"
+ @echo "# qcollectiongenerator $(BUILDDIR)/qthelp/CeresSolver.qhcp"
+ @echo "To view the help file:"
+ @echo "# assistant -collectionFile $(BUILDDIR)/qthelp/CeresSolver.qhc"
+
+devhelp:
+ $(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp
+ @echo
+ @echo "Build finished."
+ @echo "To view the help file:"
+ @echo "# mkdir -p $$HOME/.local/share/devhelp/CeresSolver"
+ @echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/CeresSolver"
+ @echo "# devhelp"
+
+epub:
+ $(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub
+ @echo
+ @echo "Build finished. The epub file is in $(BUILDDIR)/epub."
+
+latex:
+ $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
+ @echo
+ @echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
+ @echo "Run \`make' in that directory to run these through (pdf)latex" \
+ "(use \`make latexpdf' here to do that automatically)."
+
+latexpdf:
+ $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
+ @echo "Running LaTeX files through pdflatex..."
+ $(MAKE) -C $(BUILDDIR)/latex all-pdf
+ @echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."
+
+text:
+ $(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text
+ @echo
+ @echo "Build finished. The text files are in $(BUILDDIR)/text."
+
+man:
+ $(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man
+ @echo
+ @echo "Build finished. The manual pages are in $(BUILDDIR)/man."
+
+texinfo:
+ $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
+ @echo
+ @echo "Build finished. The Texinfo files are in $(BUILDDIR)/texinfo."
+ @echo "Run \`make' in that directory to run these through makeinfo" \
+ "(use \`make info' here to do that automatically)."
+
+info:
+ $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
+ @echo "Running Texinfo files through makeinfo..."
+ make -C $(BUILDDIR)/texinfo info
+ @echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo."
+
+gettext:
+ $(SPHINXBUILD) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale
+ @echo
+ @echo "Build finished. The message catalogs are in $(BUILDDIR)/locale."
+
+changes:
+ $(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes
+ @echo
+ @echo "The overview file is in $(BUILDDIR)/changes."
+
+linkcheck:
+ $(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck
+ @echo
+ @echo "Link check complete; look for any errors in the above output " \
+ "or in $(BUILDDIR)/linkcheck/output.txt."
+
+doctest:
+ $(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest
+ @echo "Testing of doctests in the sources finished, look at the " \
+ "results in $(BUILDDIR)/doctest/output.txt."
diff --git a/docs/make.bat b/docs/make.bat
new file mode 100644
index 0000000..005a432
--- /dev/null
+++ b/docs/make.bat
@@ -0,0 +1,190 @@
+@ECHO OFF
+
+REM Command file for Sphinx documentation
+
+if "%SPHINXBUILD%" == "" (
+ set SPHINXBUILD=sphinx-build
+)
+set BUILDDIR=build
+set ALLSPHINXOPTS=-d %BUILDDIR%/doctrees %SPHINXOPTS% source
+set I18NSPHINXOPTS=%SPHINXOPTS% source
+if NOT "%PAPER%" == "" (
+ set ALLSPHINXOPTS=-D latex_paper_size=%PAPER% %ALLSPHINXOPTS%
+ set I18NSPHINXOPTS=-D latex_paper_size=%PAPER% %I18NSPHINXOPTS%
+)
+
+if "%1" == "" goto help
+
+if "%1" == "help" (
+ :help
+ echo.Please use `make ^<target^>` where ^<target^> is one of
+ echo. html to make standalone HTML files
+ echo. dirhtml to make HTML files named index.html in directories
+ echo. singlehtml to make a single large HTML file
+ echo. pickle to make pickle files
+ echo. json to make JSON files
+ echo. htmlhelp to make HTML files and a HTML help project
+ echo. qthelp to make HTML files and a qthelp project
+ echo. devhelp to make HTML files and a Devhelp project
+ echo. epub to make an epub
+ echo. latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter
+ echo. text to make text files
+ echo. man to make manual pages
+ echo. texinfo to make Texinfo files
+ echo. gettext to make PO message catalogs
+ echo. changes to make an overview over all changed/added/deprecated items
+ echo. linkcheck to check all external links for integrity
+ echo. doctest to run all doctests embedded in the documentation if enabled
+ goto end
+)
+
+if "%1" == "clean" (
+ for /d %%i in (%BUILDDIR%\*) do rmdir /q /s %%i
+ del /q /s %BUILDDIR%\*
+ goto end
+)
+
+if "%1" == "html" (
+ %SPHINXBUILD% -b html %ALLSPHINXOPTS% %BUILDDIR%/html
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The HTML pages are in %BUILDDIR%/html.
+ goto end
+)
+
+if "%1" == "dirhtml" (
+ %SPHINXBUILD% -b dirhtml %ALLSPHINXOPTS% %BUILDDIR%/dirhtml
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The HTML pages are in %BUILDDIR%/dirhtml.
+ goto end
+)
+
+if "%1" == "singlehtml" (
+ %SPHINXBUILD% -b singlehtml %ALLSPHINXOPTS% %BUILDDIR%/singlehtml
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The HTML pages are in %BUILDDIR%/singlehtml.
+ goto end
+)
+
+if "%1" == "pickle" (
+ %SPHINXBUILD% -b pickle %ALLSPHINXOPTS% %BUILDDIR%/pickle
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished; now you can process the pickle files.
+ goto end
+)
+
+if "%1" == "json" (
+ %SPHINXBUILD% -b json %ALLSPHINXOPTS% %BUILDDIR%/json
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished; now you can process the JSON files.
+ goto end
+)
+
+if "%1" == "htmlhelp" (
+ %SPHINXBUILD% -b htmlhelp %ALLSPHINXOPTS% %BUILDDIR%/htmlhelp
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished; now you can run HTML Help Workshop with the ^
+.hhp project file in %BUILDDIR%/htmlhelp.
+ goto end
+)
+
+if "%1" == "qthelp" (
+ %SPHINXBUILD% -b qthelp %ALLSPHINXOPTS% %BUILDDIR%/qthelp
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished; now you can run "qcollectiongenerator" with the ^
+.qhcp project file in %BUILDDIR%/qthelp, like this:
+ echo.^> qcollectiongenerator %BUILDDIR%\qthelp\CeresSolver.qhcp
+ echo.To view the help file:
+ echo.^> assistant -collectionFile %BUILDDIR%\qthelp\CeresSolver.ghc
+ goto end
+)
+
+if "%1" == "devhelp" (
+ %SPHINXBUILD% -b devhelp %ALLSPHINXOPTS% %BUILDDIR%/devhelp
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished.
+ goto end
+)
+
+if "%1" == "epub" (
+ %SPHINXBUILD% -b epub %ALLSPHINXOPTS% %BUILDDIR%/epub
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The epub file is in %BUILDDIR%/epub.
+ goto end
+)
+
+if "%1" == "latex" (
+ %SPHINXBUILD% -b latex %ALLSPHINXOPTS% %BUILDDIR%/latex
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished; the LaTeX files are in %BUILDDIR%/latex.
+ goto end
+)
+
+if "%1" == "text" (
+ %SPHINXBUILD% -b text %ALLSPHINXOPTS% %BUILDDIR%/text
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The text files are in %BUILDDIR%/text.
+ goto end
+)
+
+if "%1" == "man" (
+ %SPHINXBUILD% -b man %ALLSPHINXOPTS% %BUILDDIR%/man
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The manual pages are in %BUILDDIR%/man.
+ goto end
+)
+
+if "%1" == "texinfo" (
+ %SPHINXBUILD% -b texinfo %ALLSPHINXOPTS% %BUILDDIR%/texinfo
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The Texinfo files are in %BUILDDIR%/texinfo.
+ goto end
+)
+
+if "%1" == "gettext" (
+ %SPHINXBUILD% -b gettext %I18NSPHINXOPTS% %BUILDDIR%/locale
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Build finished. The message catalogs are in %BUILDDIR%/locale.
+ goto end
+)
+
+if "%1" == "changes" (
+ %SPHINXBUILD% -b changes %ALLSPHINXOPTS% %BUILDDIR%/changes
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.The overview file is in %BUILDDIR%/changes.
+ goto end
+)
+
+if "%1" == "linkcheck" (
+ %SPHINXBUILD% -b linkcheck %ALLSPHINXOPTS% %BUILDDIR%/linkcheck
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Link check complete; look for any errors in the above output ^
+or in %BUILDDIR%/linkcheck/output.txt.
+ goto end
+)
+
+if "%1" == "doctest" (
+ %SPHINXBUILD% -b doctest %ALLSPHINXOPTS% %BUILDDIR%/doctest
+ if errorlevel 1 exit /b 1
+ echo.
+ echo.Testing of doctests in the sources finished, look at the ^
+results in %BUILDDIR%/doctest/output.txt.
+ goto end
+)
+
+:end
diff --git a/docs/source/acknowledgements.rst b/docs/source/acknowledgements.rst
new file mode 100644
index 0000000..049a380
--- /dev/null
+++ b/docs/source/acknowledgements.rst
@@ -0,0 +1,25 @@
+.. _chapter-acknowledgements:
+
+================
+Acknowledgements
+================
+
+A number of people have helped with the development and open sourcing
+of Ceres.
+
+Fredrik Schaffalitzky when he was at Google started the development of
+Ceres, and even though much has changed since then, many of the ideas
+from his original design are still present in the current code.
+
+Amongst Ceres' users at Google two deserve special mention: William
+Rucklidge and James Roseborough. William was the first user of
+Ceres. He bravely took on the task of porting production code to an
+as-yet unproven optimization library, reporting bugs and helping fix
+them along the way. James is perhaps the most sophisticated user of
+Ceres at Google. He has reported and fixed bugs and helped evolve the
+API for the better.
+
+Since the initial release of Ceres, a number of people have
+contributed to Ceres by porting it to new platforms, reporting bugs,
+fixing bugs and adding new functionality. We acknowledge all of these
+contributions in :ref:`chapter-version-history`.
diff --git a/docs/source/bibliography.rst b/docs/source/bibliography.rst
new file mode 100644
index 0000000..a1a9883
--- /dev/null
+++ b/docs/source/bibliography.rst
@@ -0,0 +1,103 @@
+.. _sec-bibliography:
+
+============
+Bibliography
+============
+
+.. [Agarwal] S. Agarwal, N. Snavely, S. M. Seitz and R. Szeliski,
+ **Bundle Adjustment in the Large**, *Proceedings of the European
+ Conference on Computer Vision*, pp. 29--42, 2010.
+
+.. [Bjorck] A. Bjorck, **Numerical Methods for Least Squares
+ Problems**, SIAM, 1996
+
+.. [Brown] D. C. Brown, **A solution to the general problem of
+ multiple station analytical stereo triangulation**, Technical
+ Report 43, Patrick Airforce Base, Florida, 1958.
+
+.. [Byrd] R.H. Byrd, R.B. Schnabel, and G.A. Shultz, **Approximate
+ solution of the trust region problem by minimization over
+ two dimensional subspaces**, *Mathematical programming*,
+ 40(1):247–263, 1988.
+
+.. [Chen] Y. Chen, T. A. Davis, W. W. Hager, and
+ S. Rajamanickam, **Algorithm 887: CHOLMOD, Supernodal Sparse
+ Cholesky Factorization and Update/Downdate**, *TOMS*, 35(3), 2008.
+
+.. [Conn] A.R. Conn, N.I.M. Gould, and P.L. Toint, **Trust region
+ methods**, *Society for Industrial Mathematics*, 2000.
+
+.. [GolubPereyra] G.H. Golub and V. Pereyra, **The differentiation of
+ pseudo-inverses and nonlinear least squares problems whose
+ variables separate**, *SIAM Journal on numerical analysis*,
+ 10(2):413–432, 1973.
+
+.. [HartleyZisserman] R.I. Hartley & A. Zisserman, **Multiview
+ Geometry in Computer Vision**, Cambridge University Press, 2004.
+
+.. [KushalAgarwal] A. Kushal and S. Agarwal, **Visibility based
+ preconditioning for bundle adjustment**, *In Proceedings of the
+ IEEE Conference on Computer Vision and Pattern Recognition*, 2012.
+
+.. [Levenberg] K. Levenberg, **A method for the solution of certain
+ nonlinear problems in least squares**, *Quart. Appl. Math*,
+ 2(2):164–168, 1944.
+
+.. [LiSaad] Na Li and Y. Saad, **MIQR: A multilevel incomplete qr
+ preconditioner for large sparse least squares problems**, *SIAM
+ Journal on Matrix Analysis and Applications*, 28(2):524–550, 2007.
+
+.. [Madsen] K. Madsen, H.B. Nielsen, and O. Tingleff, **Methods for
+ nonlinear least squares problems**, 2004.
+
+.. [Mandel] J. Mandel, **On block diagonal and Schur complement
+ preconditioning**, *Numer. Math.*, 58(1):79–93, 1990.
+
+.. [Marquardt] D.W. Marquardt, **An algorithm for least squares
+ estimation of nonlinear parameters**, *J. SIAM*, 11(2):431–441,
+ 1963.
+
+.. [Mathew] T.P.A. Mathew, **Domain decomposition methods for the
+ numerical solution of partial differential equations**, Springer
+ Verlag, 2008.
+
+.. [NashSofer] S.G. Nash and A. Sofer, **Assessing a search direction
+ within a truncated newton method**, *Operations Research Letters*,
+ 9(4):219–221, 1990.
+
+.. [NocedalWright] J. Nocedal & S. Wright, **Numerical Optimization**,
+ Springer, 2004.
+
+.. [RuheWedin] A. Ruhe and P.Å. Wedin, **Algorithms for separable
+ nonlinear least squares problems**, Siam Review, 22(3):318–337,
+ 1980.
+
+.. [Saad] Y. Saad, **Iterative methods for sparse linear
+ systems**, SIAM, 2003.
+
+.. [Stigler] S. M. Stigler, **Gauss and the invention of least
+ squares**, *The Annals of Statistics*, 9(3):465-474, 1981.
+
+.. [TenenbaumDirector] J. Tenenbaum & B. Director, **How Gauss
+ Determined the Orbit of Ceres**.
+
+.. [TrefethenBau] L.N. Trefethen and D. Bau, **Numerical Linear
+ Algebra**, SIAM, 1997.
+
+.. [Triggs] B. Triggs, P. F. Mclauchlan, R. I. Hartley &
+ A. W. Fitzgibbon, **Bundle Adjustment: A Modern Synthesis**,
+ Proceedings of the International Workshop on Vision Algorithms:
+ Theory and Practice, pp. 298-372, 1999.
+
+.. [Wiberg] T. Wiberg, **Computation of principal components when data
+ are missing**, In Proc. *Second Symp. Computational Statistics*,
+ pages 229–236, 1976.
+
+.. [WrightHolt] S. J. Wright and J. N. Holt, **An Inexact
+ Levenberg Marquardt Method for Large Sparse Nonlinear Least
+ Squares**, *Journal of the Australian Mathematical Society Series
+ B*, 26(4):387–403, 1985.
+
+
+
+
diff --git a/docs/source/building.rst b/docs/source/building.rst
new file mode 100644
index 0000000..3edcac8
--- /dev/null
+++ b/docs/source/building.rst
@@ -0,0 +1,374 @@
+.. _chapter-building:
+
+========
+Building
+========
+
+Ceres source code and documentation are hosted at
+http://code.google.com/p/ceres-solver/ .
+
+.. _section-dependencies:
+
+Dependencies
+============
+
+Ceres relies on a number of open source libraries, some of which are
+optional. For details on customizing the build process, see
+:ref:`section-customizing` .
+
+1. `CMake <http://www.cmake.org>`_ is a cross platform build
+system. Ceres needs a relatively recent version of CMake (version
+2.8.0 or better).
+
+2. `eigen3 <http://eigen.tuxfamily.org/index.php?title=Main_Page>`_ is
+used for doing all the low level matrix and linear algebra operations.
+
+3. `google-glog <http://http://code.google.com/p/google-glog>`_ is
+used for error checking and logging. Ceres needs glog version 0.3.1 or
+later. Version 0.3 (which ships with Fedora 16) has a namespace bug
+which prevents Ceres from building.
+
+4. `gflags <http://code.google.com/p/gflags>`_ is a library for
+processing command line flags. It is used by some of the examples and
+tests. While it is not strictly necessary to build the library, we
+strongly recommend building the library with gflags.
+
+
+5. `SuiteSparse
+<http://www.cise.ufl.edu/research/sparse/SuiteSparse/>`_ is used for
+sparse matrix analysis, ordering and factorization. In particular
+Ceres uses the AMD, COLAMD and CHOLMOD libraries. This is an optional
+dependency.
+
+6. `CXSparse <http://www.cise.ufl.edu/research/sparse/CXSparse/>`_ is
+used for sparse matrix analysis, ordering and factorization. While it
+is similar to SuiteSparse in scope, its performance is a bit worse but
+is a much simpler library to build and does not have any other
+dependencies. This is an optional dependency.
+
+7. `BLAS <http://www.netlib.org/blas/>`_ and `LAPACK
+<http://www.netlib.org/lapack/>`_ routines are needed by
+SuiteSparse. We recommend either `GotoBlas2
+<http://www.tacc.utexas.edu/tacc- projects/gotoblas2>`_ or `ATLAS
+<http://math- atlas.sourceforge.net/>`_ , both of which ship with BLAS
+and LAPACK routines.
+
+8. `protobuf <http://code.google.com/p/protobuf/>`_ is used for
+serializing and deserializing linear least squares problems to
+disk. This is useful for debugging and testing. It is an optional
+depdendency and without it some of the tests will be disabled.
+
+.. _section-linux:
+
+Building on Linux
+=================
+We will use `Ubuntu <http://www.ubuntu.com>`_ as our example platform.
+
+#. ``CMake``
+
+ .. code-block:: bash
+
+ sudo apt-get install cmake
+
+#. ``gflags`` can either be installed from source via the ``autoconf``
+ invocation
+
+ .. code-block:: bash
+
+ tar -xvzf gflags-2.0.tar.gz
+ cd gflags-2.0
+ ./configure --prefix=/usr/local
+ make
+ sudo make install.
+
+
+ or via the ``deb`` or ``rpm`` packages available on the ``gflags`` website.
+
+#. ``google-glog`` must be configured to use the previously installed
+ ``gflags``, rather than the stripped down version that is bundled
+ with ``google-glog``. Assuming you have it installed in ``/usr/local`` the
+ following ``autoconf`` invocation installs it.
+
+ .. code-block:: bash
+
+ tar -xvzf glog-0.3.2.tar.gz
+ cd glog-0.3.2
+ ./configure --with-gflags=/usr/local/
+ make
+ sudo make install
+
+#. ``Eigen``
+
+ .. code-block:: bash
+
+ sudo apt-get install libeigen3-dev
+
+#. ``SuiteSparse`` and ``CXSparse``
+
+ .. code-block:: bash
+
+ sudo apt-get install libsuitesparse-dev
+
+ This should automatically bring in the necessary ``BLAS`` and
+ ``LAPACK`` dependencies. By co-incidence on Ubuntu, this also
+ installs ``CXSparse``.
+
+#. ``protobuf``
+
+ .. code-block:: bash
+
+ sudo apt-get install libprotobuf-dev
+
+
+We are now ready to build and test Ceres. Note that ``CMake`` requires
+the exact path to the ``libglog.a`` and ``libgflag.a``.
+
+.. code-block:: bash
+
+ tar zxf ceres-solver-1.2.1.tar.gz
+ mkdir ceres-bin
+ cd ceres-bin
+ cmake ../ceres-solver-1.2.1
+ make -j3
+ make test
+
+You can also try running the command line bundling application with one of the
+included problems, which comes from the University of Washington's BAL
+dataset [Agarwal]_.
+
+.. code-block:: bash
+
+ bin/simple_bundle_adjuster \
+ ../ceres-solver-1.2.1/data/problem-16-22106-pre.txt \
+
+This runs Ceres for a maximum of 10 iterations using the
+``DENSE_SCHUR`` linear solver. The output should look something like
+this.
+
+.. code-block:: bash
+
+ 0: f: 4.185660e+06 d: 0.00e+00 g: 1.09e+08 h: 0.00e+00 rho: 0.00e+00 mu: 1.00e+04 li: 0 it: 1.16e-01 tt: 3.39e-01
+ 1: f: 1.062590e+05 d: 4.08e+06 g: 8.99e+06 h: 5.36e+02 rho: 9.82e-01 mu: 3.00e+04 li: 1 it: 3.90e-01 tt: 7.29e-01
+ 2: f: 4.992817e+04 d: 5.63e+04 g: 8.32e+06 h: 3.19e+02 rho: 6.52e-01 mu: 3.09e+04 li: 1 it: 3.52e-01 tt: 1.08e+00
+ 3: f: 1.899774e+04 d: 3.09e+04 g: 1.60e+06 h: 1.24e+02 rho: 9.77e-01 mu: 9.26e+04 li: 1 it: 3.60e-01 tt: 1.44e+00
+ 4: f: 1.808729e+04 d: 9.10e+02 g: 3.97e+05 h: 6.39e+01 rho: 9.51e-01 mu: 2.78e+05 li: 1 it: 3.62e-01 tt: 1.80e+00
+ 5: f: 1.803399e+04 d: 5.33e+01 g: 1.48e+04 h: 1.23e+01 rho: 9.99e-01 mu: 8.33e+05 li: 1 it: 3.54e-01 tt: 2.16e+00
+ 6: f: 1.803390e+04 d: 9.02e-02 g: 6.35e+01 h: 8.00e-01 rho: 1.00e+00 mu: 2.50e+06 li: 1 it: 3.59e-01 tt: 2.52e+00
+
+ Ceres Solver Report
+ -------------------
+ Original Reduced
+ Parameter blocks 22122 22122
+ Parameters 66462 66462
+ Residual blocks 83718 83718
+ Residual 167436 167436
+ Trust Region Strategy LEVENBERG_MARQUARDT
+
+ Given Used
+ Linear solver DENSE_SCHUR DENSE_SCHUR
+ Preconditioner N/A N/A
+ Threads: 1 1
+ Linear solver threads 1 1
+ Linear solver ordering AUTOMATIC 22106,16
+
+ Cost:
+ Initial 4.185660e+06
+ Final 1.803390e+04
+ Change 4.167626e+06
+
+ Number of iterations:
+ Successful 6
+ Unsuccessful 0
+ Total 6
+
+ Time (in seconds):
+ Preprocessor 2.229e-01
+
+ Evaluator::Residuals 7.438e-02
+ Evaluator::Jacobians 6.790e-01
+ Linear Solver 1.681e+00
+ Minimizer 2.547e+00
+
+ Postprocessor 1.920e-02
+ Total 2.823e+00
+
+ Termination: FUNCTION_TOLERANCE
+
+.. section-osx:
+
+Building on Mac OS X
+====================
+
+On OS X, we recommend using the `homebrew
+<http://mxcl.github.com/homebrew/>`_ package manager.
+
+
+#. ``CMake``
+
+ .. code-block:: bash
+
+ brew install cmake
+
+#. ``google-glog`` and ``gflags``
+
+Installing ``google-glog`` takes also brings in ``gflags`` as a dependency.
+
+ .. code-block:: bash
+
+ brew install glog
+
+#. ``Eigen3``
+
+ .. code-block:: bash
+
+ brew install eigen
+
+#. ``SuiteSparse`` and ``CXSparse``
+
+ .. code-block:: bash
+
+ brew install suite-sparse
+
+#. ``protobuf``
+
+ .. code-block:: bash
+
+ brew install protobuf
+
+
+We are now ready to build and test Ceres.
+
+.. code-block:: bash
+
+ tar zxf ceres-solver-1.2.1.tar.gz
+ mkdir ceres-bin
+ cd ceres-bin
+ cmake ../ceres-solver-1.2.1
+ make -j3
+ make test
+
+
+Like the Linux build, you should now be able to run
+``bin/simple_bundle_adjuster``.
+
+.. _section-windows:
+
+Building on Windows with Visual Studio
+======================================
+
+On Windows, we support building with Visual Studio 2010 or newer. Note
+that the Windows port is less featureful and less tested than the
+Linux or Mac OS X versions due to the unavaliability of SuiteSparse
+and ``CXSparse``. Building is also more involved since there is no
+automated way to install the dependencies.
+
+#. Make a toplevel directory for deps & build & src somewhere: ``ceres/``
+#. Get dependencies; unpack them as subdirectories in ``ceres/``
+ (``ceres/eigen``, ``ceres/glog``, etc)
+
+ #. ``Eigen`` 3.1 (needed on Windows; 3.0.x will not work). There is
+ no need to build anything; just unpack the source tarball.
+
+ #. ``google-glog`` Open up the Visual Studio solution and build it.
+ #. ``gflags`` Open up the Visual Studio solution and build it.
+
+#. Unpack the Ceres tarball into ``ceres``. For the tarball, you
+ should get a directory inside ``ceres`` similar to
+ ``ceres-solver-1.3.0``. Alternately, checkout Ceres via ``git`` to
+ get ``ceres-solver.git`` inside ``ceres``.
+
+#. Install ``CMake``,
+
+#. Make a dir ``ceres/ceres-bin`` (for an out-of-tree build)
+
+#. Run ``CMake``; select the ``ceres-solver-X.Y.Z`` or
+ ``ceres-solver.git`` directory for the CMake file. Then select the
+ ``ceres-bin`` for the build dir.
+
+#. Try running "Configure". It won't work. It'll show a bunch of options.
+ You'll need to set:
+
+ #. ``GLOG_INCLUDE``
+ #. ``GLOG_LIB``
+ #. ``GFLAGS_LIB``
+ #. ``GFLAGS_INCLUDE``
+
+ to the appropriate place where you unpacked/built them.
+
+#. You may have to tweak some more settings to generate a MSVC
+ project. After each adjustment, try pressing Configure & Generate
+ until it generates successfully.
+
+#. Open the solution and build it in MSVC
+
+
+To run the tests, select the ``RUN_TESTS`` target and hit **Build
+RUN_TESTS** from the build menu.
+
+Like the Linux build, you should now be able to run ``bin/simple_bundle_adjuster``.
+
+Notes:
+
+#. The default build is Debug; consider switching it to release mode.
+#. Currently ``system_test`` is not working properly.
+#. Building Ceres as a DLL is not supported; patches welcome.
+#. CMake puts the resulting test binaries in ``ceres-bin/examples/Debug``
+ by default.
+#. The solvers supported on Windows are ``DENSE_QR``, ``DENSE_SCHUR``,
+ ``CGNR``, and ``ITERATIVE_SCHUR``.
+#. We're looking for someone to work with upstream ``SuiteSparse`` to
+ port their build system to something sane like ``CMake``, and get a
+ supported Windows port.
+
+
+.. _section-android:
+
+Building on Android
+===================
+
+
+Download the ``Android NDK``. Run ``ndk-build`` from inside the
+``jni`` directory. Use the ``libceres.a`` that gets created.
+
+.. _section-customizing:
+
+Customizing the build
+=====================
+
+It is possible to reduce the libraries needed to build Ceres and
+customize the build process by passing appropriate flags to
+``CMake``. Use these flags only if you really know what you are doing.
+
+#. ``-DPROTOBUF=OFF`` : ``protobuf`` is a big dependency and if you do not
+ care for the tests that depend on it and the logging support it
+ enables, you can use this flag to turn it off.
+
+#. ``-DSUITESPARSE=OFF`` : By default, Ceres will link to
+ ``SuiteSparse`` if all its dependencies are present. Use this flag
+ to buils Ceres without ``SuiteSparse``. This will also disable
+ dependency checking for ``LAPACK`` and ``BLAS`` This saves on
+ binary size, but the resulting version of Ceres is not suited to
+ large scale problems due to the lack of a sparse Cholesky solver.
+ This will reduce Ceres' dependencies down to ``Eigen``, ``gflags``
+ and ``google-glog``.
+
+#. ``-DCXSPARSE=OFF`` : By default, Ceres will link to ``CXSparse`` if all
+ its dependencies are present. Use this flag to buils Ceres without
+ ``CXSparse``. This saves on binary size, but the resulting version
+ of Ceres is not suited to large scale problems due to the lack of a
+ sparse Cholesky solver. This will reduce Ceres' dependencies down
+ to ``Eigen``, ``gflags`` and ``google-glog``.
+
+#. ``-DGFLAGS=OFF`` : Use this flag to build Ceres without
+ ``gflags``. This will also prevent some of the example code from
+ building.
+
+#. ``-DSCHUR_SPECIALIZATIONS=OFF`` : If you are concerned about binary
+ size/compilation time over some small (10-20%) performance gains in
+ the ``SPARSE_SCHUR`` solver, you can disable some of the template
+ specializations by using this flag.
+
+#. ``-DOPENMP=OFF`` : On certain platforms like Android,
+ multithreading with ``OpenMP`` is not supported. Use this flag to
+ disable multithreading.
+
diff --git a/docs/source/conf.py b/docs/source/conf.py
new file mode 100644
index 0000000..b63f914
--- /dev/null
+++ b/docs/source/conf.py
@@ -0,0 +1,242 @@
+# -*- coding: utf-8 -*-
+#
+# Ceres Solver documentation build configuration file, created by
+# sphinx-quickstart on Sun Jan 20 20:34:07 2013.
+#
+# This file is execfile()d with the current directory set to its containing dir.
+#
+# Note that not all possible configuration values are present in this
+# autogenerated file.
+#
+# All configuration values have a default; values that are commented out
+# serve to show the default.
+
+import sys, os
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#sys.path.insert(0, os.path.abspath('.'))
+
+# -- General configuration -----------------------------------------------------
+
+# If your documentation needs a minimal Sphinx version, state it here.
+#needs_sphinx = '1.0'
+
+# Add any Sphinx extension module names here, as strings. They can be extensions
+# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
+extensions = ['sphinx.ext.todo', 'sphinx.ext.mathjax', 'sphinx.ext.ifconfig']
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# The suffix of source filenames.
+source_suffix = '.rst'
+
+# The encoding of source files.
+#source_encoding = 'utf-8-sig'
+
+# The master toctree document.
+master_doc = 'index'
+
+# General information about the project.
+project = u'Ceres Solver'
+copyright = u'2013, Google Inc.'
+
+# The version info for the project you're documenting, acts as replacement for
+# |version| and |release|, also used in various other places throughout the
+# built documents.
+#
+# The short X.Y version.
+version = '1.5.0'
+# The full version, including alpha/beta/rc tags.
+release = '1.5.0'
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#language = None
+
+# There are two options for replacing |today|: either, you set today to some
+# non-false value, then it is used:
+#today = ''
+# Else, today_fmt is used as the format for a strftime call.
+#today_fmt = '%B %d, %Y'
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+exclude_patterns = []
+
+# The reST default role (used for this markup: `text`) to use for all documents.
+#default_role = None
+
+# If true, '()' will be appended to :func: etc. cross-reference text.
+#add_function_parentheses = True
+
+# If true, the current module name will be prepended to all description
+# unit titles (such as .. function::).
+#add_module_names = True
+
+# If true, sectionauthor and moduleauthor directives will be shown in the
+# output. They are ignored by default.
+#show_authors = False
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx'
+
+# A list of ignored prefixes for module index sorting.
+#modindex_common_prefix = []
+
+
+# -- Options for HTML output ---------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages. See the documentation for
+# a list of builtin themes.
+html_theme = 'nature'
+
+# Theme options are theme-specific and customize the look and feel of a theme
+# further. For a list of options available for each theme, see the
+# documentation.
+#html_theme_options = {}
+
+# Add any paths that contain custom themes here, relative to this directory.
+html_theme_path = ["_themes",]
+
+# The name for this set of Sphinx documents. If None, it defaults to
+# "<project> v<release> documentation".
+html_title = "The Ceres Solver Manual"
+
+# A shorter title for the navigation bar. Default is the same as html_title.
+#html_short_title = None
+
+# The name of an image file (relative to this directory) to place at the top
+# of the sidebar.
+#html_logo = None
+
+# The name of an image file (within the static path) to use as favicon of the
+# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32
+# pixels large.
+#html_favicon = None
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
+
+# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
+# using the given strftime format.
+#html_last_updated_fmt = '%b %d, %Y'
+
+# If true, SmartyPants will be used to convert quotes and dashes to
+# typographically correct entities.
+html_use_smartypants = True
+
+# Custom sidebar templates, maps document names to template names.
+#html_sidebars = {}
+
+# Additional templates that should be rendered to pages, maps page names to
+# template names.
+#html_additional_pages = {}
+
+# If false, no module index is generated.
+html_domain_indices = True
+
+# If false, no index is generated.
+html_use_index = True
+
+# If true, the index is split into individual pages for each letter.
+html_split_index = False
+
+# If true, links to the reST sources are added to the pages.
+html_show_sourcelink = True
+
+# If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
+html_show_sphinx = False
+
+# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
+html_show_copyright = True
+
+# If true, an OpenSearch description file will be output, and all pages will
+# contain a <link> tag referring to it. The value of this option must be the
+# base URL from which the finished HTML is served.
+#html_use_opensearch = ''
+
+# This is the file name suffix for HTML files (e.g. ".xhtml").
+#html_file_suffix = None
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'CeresSolverdoc'
+
+
+# -- Options for LaTeX output --------------------------------------------------
+
+latex_elements = {
+# The paper size ('letterpaper' or 'a4paper').
+#'papersize': 'letterpaper',
+
+# The font size ('10pt', '11pt' or '12pt').
+#'pointsize': '10pt',
+
+# Additional stuff for the LaTeX preamble.
+#'preamble': '',
+}
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title, author, documentclass [howto/manual]).
+latex_documents = [
+ ('index', 'CeresSolver.tex', u'The Ceres Solver Manual',
+ u'Sameer Agarwal \\& Keir Mierle', 'manual'),
+]
+
+# The name of an image file (relative to this directory) to place at the top of
+# the title page.
+#latex_logo = None
+
+# For "manual" documents, if this is true, then toplevel headings are parts,
+# not chapters.
+#latex_use_parts = False
+
+# If true, show page references after internal links.
+#latex_show_pagerefs = False
+
+# If true, show URL addresses after external links.
+#latex_show_urls = False
+
+# Documents to append as an appendix to all manuals.
+#latex_appendices = []
+
+# If false, no module index is generated.
+#latex_domain_indices = True
+
+
+# -- Options for manual page output --------------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+man_pages = [
+ ('index', 'ceressolver', u'The Ceres Solver Manual',
+ [u'Sameer Agarwal & Keir Mierle'], 1)
+]
+
+# If true, show URL addresses after external links.
+#man_show_urls = False
+
+
+# -- Options for Texinfo output ------------------------------------------------
+
+# Grouping the document tree into Texinfo files. List of tuples
+# (source start file, target name, title, author,
+# dir menu entry, description, category)
+texinfo_documents = [
+ ('index', 'CeresSolver', u'The Ceres Solver Manual',
+ u'Sameer Agarwal & Keir Mierle', 'CeresSolver', 'One line description of project.',
+ 'Miscellaneous'),
+]
+
+# Documents to append as an appendix to all manuals.
+#texinfo_appendices = []
+
+# If false, no module index is generated.
+#texinfo_domain_indices = True
+
+# How to display URL addresses: 'footnote', 'no', or 'inline'.
+#texinfo_show_urls = 'footnote'
diff --git a/docs/source/contributing.rst b/docs/source/contributing.rst
new file mode 100644
index 0000000..5e0347e
--- /dev/null
+++ b/docs/source/contributing.rst
@@ -0,0 +1,29 @@
+.. _chapter-contributing:
+
+=============
+Contributions
+=============
+
+
+We welcome contributions to Ceres, whether they are new features, bug
+fixes or tests. The Ceres mailing list [#f1]_ is the best place for
+all development related discussions. Please consider joining it. If
+you have ideas on how you would like to contribute to Ceres, it is a
+good idea to let us know on the mailinglist before you start
+development. We may have suggestions that will save effort when trying
+to merge your work into the main branch. If you are looking for ideas,
+please let us know about your interest and skills and we will be happy
+to make a suggestion or three.
+
+We follow Google's C++ Style Guide [#f2]_ and use ``git`` for version
+control.
+
+
+Gerrit Instructions
+===================
+
+TBD
+.. rubric:: Footnotes
+
+.. [#f1] http://groups.google.com/group/ceres-solver
+.. [#f2] http://google-styleguide.googlecode.com/svn/trunk/cppguide.xml
diff --git a/docs/source/fit.png b/docs/source/fit.png
new file mode 100644
index 0000000..b80e46b
--- /dev/null
+++ b/docs/source/fit.png
Binary files differ
diff --git a/docs/source/index.rst b/docs/source/index.rst
new file mode 100644
index 0000000..eaf9c90
--- /dev/null
+++ b/docs/source/index.rst
@@ -0,0 +1,94 @@
+.. Ceres Solver documentation master file, created by
+ sphinx-quickstart on Sat Jan 19 00:07:33 2013.
+ You can adapt this file completely to your liking, but it should at least
+ contain the root `toctree` directive.
+
+Introduction
+============
+
+Solving nonlinear least squares problems [#f1]_ comes up in a broad
+range of areas across science and engineering - from fitting curves in
+statistics, to constructing 3D models from photographs in computer
+vision. Ceres Solver [#f2]_ [#f3]_ is a portable C++ library for
+solving non-linear least squares problems. It is designed to solve
+small and large sparse problems accurately and efficiently.
+
+At Google, Ceres Solver has been used for solving a variety of
+problems in computer vision and machine learning. e.g., it is used to
+to estimate the pose of Street View cars, aircrafts, and satellites;
+to build 3D models for PhotoTours; to estimate satellite image sensor
+characteristics, and more.
+
+
+Features:
+
+#. A friendly :ref:`chapter-modeling`.
+
+#. Automatic and numeric differentiation.
+
+#. Robust loss functions and Local parameterizations.
+
+#. Multithreading.
+
+#. Trust-Region (Levenberg-Marquardt and Dogleg) and Line Search
+ (Nonlinear CG and L-BFGS) solvers.
+
+#. Variety of linear solvers.
+
+ a. Dense QR and Cholesky factorization (using `Eigen
+ <http://eigen.tuxfamily.org/index.php?title=Main_Page>`_) for
+ small problems.
+
+ b. Sparse Cholesky factorization (using `SuiteSparse
+ <http://www.cise.ufl.edu/research/sparse/SuiteSparse/>`_ and
+ `CXSparse <http://www.cise.ufl.edu/research/sparse/CSparse/>`_) for
+ large sparse problems.
+
+ c. Specialized solvers for bundle adjustment problems in computer
+ vision.
+
+ d. Iterative linear solvers with perconditioners for general sparse
+ and bundle adjustment problems.
+
+#. Portable: Runs on Linux, Windows, Mac OS X and Android. An iOS port is
+ underway.
+
+
+
+Contents
+========
+
+.. toctree::
+ :maxdepth: 1
+
+ building
+ tutorial
+ modeling
+ solving
+ reading
+ contributing
+ acknowledgements
+ version_history
+ bibliography
+ license
+
+.. rubric:: Footnotes
+
+.. [#f1] For a gentle but brief introduction to non-linear least
+ squares problems, please start by reading the
+ :ref:`chapter-tutorial`.
+
+.. [#f2] While there is some debate as to who invented of the method
+ of Least Squares [Stigler]_. There is no debate that it was
+ Carl Friedrich Gauss's prediction of the orbit of the newly
+ discovered asteroid Ceres based on just 41 days of
+ observations that brought it to the attention of the world
+ [TenenbaumDirector]_. We named our solver after Ceres to
+ celebrate this seminal event in the history of astronomy,
+ statistics and optimization.
+
+.. [#f3] For brevity, in the rest of this document we will just use
+ the term Ceres.
+
+
+
diff --git a/docs/source/license.rst b/docs/source/license.rst
new file mode 100644
index 0000000..3c962e1
--- /dev/null
+++ b/docs/source/license.rst
@@ -0,0 +1,42 @@
+=======
+License
+=======
+
+Ceres Solver is licensed under the New BSD license, whose terms are as follows.
+
+Copyright (c) 2010, 2011, 2012, 2013 Google Inc. All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice,
+ this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright notice,
+ this list of conditions and the following disclaimer in the documentation
+ and/or other materials provided with the distribution.
+3. Neither the name of Google Inc., nor the names of its contributors may
+ be used to endorse or promote products derived from this software without
+ specific prior written permission.
+
+This software is provided by the copyright holders and contributors "AS IS" and
+any express or implied warranties, including, but not limited to, the implied
+warranties of merchantability and fitness for a particular purpose are
+disclaimed. In no event shall Google Inc. be liable for any direct, indirect,
+incidental, special, exemplary, or consequential damages (including, but not
+limited to, procurement of substitute goods or services; loss of use, data, or
+profits; or business interruption) however caused and on any theory of
+liability, whether in contract, strict liability, or tort (including negligence
+or otherwise) arising in any way out of the use of this software, even if
+advised of the possibility of such damage.
+
+Citing Ceres Solver
+===================
+
+If you use Ceres Solver for an academic publication, please cite this
+manual. e.g., ::
+
+ @manual{ceres-manual,
+ Author = {Sameer Agarwal and Keir Mierle},
+ Title = {Ceres Solver: Tutorial \& Reference},
+ Organization = {Google Inc.}
+ }
diff --git a/docs/source/loss.png b/docs/source/loss.png
new file mode 100644
index 0000000..9f98d00
--- /dev/null
+++ b/docs/source/loss.png
Binary files differ
diff --git a/docs/source/modeling.rst b/docs/source/modeling.rst
new file mode 100644
index 0000000..ae5cfdf
--- /dev/null
+++ b/docs/source/modeling.rst
@@ -0,0 +1,1116 @@
+.. default-domain:: cpp
+
+
+.. cpp:namespace:: ceres
+
+
+.. _`chapter-modeling`:
+
+============
+Modeling API
+============
+
+
+Introduction
+------------
+
+Ceres solves robustified non-linear least squares problems of the form
+
+.. math:: \frac{1}{2}\sum_{i=1} \rho_i\left(\left\|f_i\left(x_{i_1}, ... ,x_{i_k}\right)\right\|^2\right).
+ :label: ceresproblem
+
+The term
+:math:`\rho_i\left(\left\|f_i\left(x_{i_1},...,x_{i_k}\right)\right\|^2\right)`
+is known as a ``ResidualBlock``, where :math:`f_i(\cdot)` is a
+:class:`CostFunction` that depends on the parameter blocks
+:math:`\left[x_{i_1},... , x_{i_k}\right]` and :math:`\rho_i` is a
+:class:`LossFunction`. In most optimization problems small groups of
+scalars occur together. For example the three components of a
+translation vector and the four components of the quaternion that
+define the pose of a camera. We refer to such a group of small scalars
+as a ``ParameterBlock``. Of course a ``ParameterBlock`` can just have
+a single parameter.
+
+:class:`CostFunction`
+---------------------
+
+.. class:: CostFunction
+
+ .. code-block:: c++
+
+ class CostFunction {
+ public:
+ virtual bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) = 0;
+ const vector<int16>& parameter_block_sizes();
+ int num_residuals() const;
+
+ protected:
+ vector<int16>* mutable_parameter_block_sizes();
+ void set_num_residuals(int num_residuals);
+ };
+
+ Given parameter blocks :math:`\left[x_{i_1}, ... , x_{i_k}\right]`,
+ a :class:`CostFunction` is responsible for computing a vector of
+ residuals and if asked a vector of Jacobian matrices, i.e., given
+ :math:`\left[x_{i_1}, ... , x_{i_k}\right]`, compute the vector
+ :math:`f_i\left(x_{i_1},...,x_{i_k}\right)` and the matrices
+
+ .. math:: J_{ij} = \frac{\partial}{\partial x_{i_j}}f_i\left(x_{i_1},...,x_{i_k}\right),\quad \forall j \in \{i_1,..., i_k\}
+
+ The signature of the class:`CostFunction` (number and sizes of
+ input parameter blocks and number of outputs) is stored in
+ :member:`CostFunction::parameter_block_sizes_` and
+ :member:`CostFunction::num_residuals_` respectively. User code
+ inheriting from this class is expected to set these two members
+ with the corresponding accessors. This information will be verified
+ by the :class:`Problem` when added with
+ :func:`Problem::AddResidualBlock`.
+
+.. function:: bool CostFunction::Evaluate(double const* const* parameters, double* residuals, double** jacobians)
+
+ This is the key methods. It implements the residual and Jacobian
+ computation.
+
+ ``parameters`` is an array of pointers to arrays containing the
+ various parameter blocks. parameters has the same number of
+ elements as :member:`CostFunction::parameter_block_sizes_`.
+ Parameter blocks are in the same order as
+ :member:`CostFunction::parameter_block_sizes_`.
+
+ ``residuals`` is an array of size ``num_residuals_``.
+
+
+ ``jacobians`` is an array of size
+ :member:`CostFunction::parameter_block_sizes_` containing pointers
+ to storage for Jacobian matrices corresponding to each parameter
+ block. The Jacobian matrices are in the same order as
+ :member:`CostFunction::parameter_block_sizes_`. ``jacobians[i]`` is
+ an array that contains :member:`CostFunction::num_residuals_` x
+ :member:`CostFunction::parameter_block_sizes_` ``[i]``
+ elements. Each Jacobian matrix is stored in row-major order, i.e.,
+ ``jacobians[i][r * parameter_block_size_[i] + c]`` =
+ :math:`\frac{\partial residual[r]}{\partial parameters[i][c]}`
+
+
+ If ``jacobians`` is ``NULL``, then no derivatives are returned;
+ this is the case when computing cost only. If ``jacobians[i]`` is
+ ``NULL``, then the Jacobian matrix corresponding to the
+ :math:`i^{\textrm{th}}` parameter block must not be returned, this
+ is the case when the a parameter block is marked constant.
+
+
+
+:class:`SizedCostFunction`
+--------------------------
+
+.. class:: SizedCostFunction
+
+ If the size of the parameter blocks and the size of the residual
+ vector is known at compile time (this is the common case), Ceres
+ provides :class:`SizedCostFunction`, where these values can be
+ specified as template parameters. In this case the user only needs
+ to implement the :func:`CostFunction::Evaluate`.
+
+ .. code-block:: c++
+
+ template<int kNumResiduals,
+ int N0 = 0, int N1 = 0, int N2 = 0, int N3 = 0, int N4 = 0,
+ int N5 = 0, int N6 = 0, int N7 = 0, int N8 = 0, int N9 = 0>
+ class SizedCostFunction : public CostFunction {
+ public:
+ virtual bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const = 0;
+ };
+
+
+:class:`AutoDiffCostFunction`
+-----------------------------
+
+.. class:: AutoDiffCostFunction
+
+ But even defining the :class:`SizedCostFunction` can be a tedious
+ affair if complicated derivative computations are involved. To this
+ end Ceres provides automatic differentiation.
+
+ To get an auto differentiated cost function, you must define a
+ class with a templated ``operator()`` (a functor) that computes the
+ cost function in terms of the template parameter ``T``. The
+ autodiff framework substitutes appropriate ``Jet`` objects for
+ ``T`` in order to compute the derivative when necessary, but this
+ is hidden, and you should write the function as if ``T`` were a
+ scalar type (e.g. a double-precision floating point number).
+
+ The function must write the computed value in the last argument
+ (the only non-``const`` one) and return true to indicate success.
+
+ For example, consider a scalar error :math:`e = k - x^\top y`,
+ where both :math:`x` and :math:`y` are two-dimensional vector
+ parameters and :math:`k` is a constant. The form of this error,
+ which is the difference between a constant and an expression, is a
+ common pattern in least squares problems. For example, the value
+ :math:`x^\top y` might be the model expectation for a series of
+ measurements, where there is an instance of the cost function for
+ each measurement :math:`k`.
+
+ The actual cost added to the total problem is :math:`e^2`, or
+ :math:`(k - x^\top y)^2`; however, the squaring is implicitly done
+ by the optimization framework.
+
+ To write an auto-differentiable cost function for the above model,
+ first define the object
+
+ .. code-block:: c++
+
+ class MyScalarCostFunctor {
+ MyScalarCostFunctor(double k): k_(k) {}
+
+ template <typename T>
+ bool operator()(const T* const x , const T* const y, T* e) const {
+ e[0] = T(k_) - x[0] * y[0] - x[1] * y[1];
+ return true;
+ }
+
+ private:
+ double k_;
+ };
+
+
+ Note that in the declaration of ``operator()`` the input parameters
+ ``x`` and ``y`` come first, and are passed as const pointers to arrays
+ of ``T``. If there were three input parameters, then the third input
+ parameter would come after ``y``. The output is always the last
+ parameter, and is also a pointer to an array. In the example above,
+ ``e`` is a scalar, so only ``e[0]`` is set.
+
+ Then given this class definition, the auto differentiated cost
+ function for it can be constructed as follows.
+
+ .. code-block:: c++
+
+ CostFunction* cost_function
+ = new AutoDiffCostFunction<MyScalarCostFunctor, 1, 2, 2>(
+ new MyScalarCostFunctor(1.0)); ^ ^ ^
+ | | |
+ Dimension of residual ------+ | |
+ Dimension of x ----------------+ |
+ Dimension of y -------------------+
+
+
+ In this example, there is usually an instance for each measurement
+ of ``k``.
+
+ In the instantiation above, the template parameters following
+ ``MyScalarCostFunction``, ``<1, 2, 2>`` describe the functor as
+ computing a 1-dimensional output from two arguments, both
+ 2-dimensional.
+
+ The framework can currently accommodate cost functions of up to 6
+ independent variables, and there is no limit on the dimensionality of
+ each of them.
+
+ **WARNING 1** Since the functor will get instantiated with
+ different types for ``T``, you must convert from other numeric
+ types to ``T`` before mixing computations with other variables
+ oftype ``T``. In the example above, this is seen where instead of
+ using ``k_`` directly, ``k_`` is wrapped with ``T(k_)``.
+
+ **WARNING 2** A common beginner's error when first using
+ :class:`AutoDiffCostFunction` 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 block*. In the example above, that
+ would be ``<MyScalarCostFunction, 1, 2>``, which is missing the 2
+ as the last template argument.
+
+
+:class:`NumericDiffCostFunction`
+--------------------------------
+
+.. class:: NumericDiffCostFunction
+
+ .. code-block:: c++
+
+ template <typename CostFunctionNoJacobian,
+ NumericDiffMethod method = CENTRAL, int M = 0,
+ int N0 = 0, int N1 = 0, int N2 = 0, int N3 = 0, int N4 = 0,
+ int N5 = 0, int N6 = 0, int N7 = 0, int N8 = 0, int N9 = 0>
+ class NumericDiffCostFunction
+ : public SizedCostFunction<M, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9> {
+ };
+
+
+ Create a :class:`CostFunction` as needed by the least squares
+ framework with jacobians computed via numeric (a.k.a. finite)
+ differentiation. For more details see
+ http://en.wikipedia.org/wiki/Numerical_differentiation.
+
+ To get an numerically differentiated :class:`CostFunction`, you
+ must define a class with a ``operator()`` (a functor) that computes
+ the residuals. The functor must write the computed value in the
+ last argument (the only non-``const`` one) and return ``true`` to
+ indicate success. e.g., an object of the form
+
+ .. code-block:: c++
+
+ struct ScalarFunctor {
+ public:
+ bool operator()(const double* const x1,
+ const double* const x2,
+ double* residuals) const;
+ }
+
+ For example, consider a scalar error :math:`e = k - x'y`, where
+ both :math:`x` and :math:`y` are two-dimensional column vector
+ parameters, the prime sign indicates transposition, and :math:`k`
+ is a constant. The form of this error, which is the difference
+ between a constant and an expression, is a common pattern in least
+ squares problems. For example, the value :math:`x'y` might be the
+ model expectation for a series of measurements, where there is an
+ instance of the cost function for each measurement :math:`k`.
+
+ To write an numerically-differentiable class:`CostFunction` for the
+ above model, first define the object
+
+ .. code-block:: c++
+
+ class MyScalarCostFunctor {
+ MyScalarCostFunctor(double k): k_(k) {}
+
+ bool operator()(const double* const x,
+ const double* const y,
+ double* residuals) const {
+ residuals[0] = k_ - x[0] * y[0] + x[1] * y[1];
+ return true;
+ }
+
+ private:
+ double k_;
+ };
+
+ Note that in the declaration of ``operator()`` the input parameters
+ ``x`` and ``y`` come first, and are passed as const pointers to
+ arrays of ``double`` s. If there were three input parameters, then
+ the third input parameter would come after ``y``. The output is
+ always the last parameter, and is also a pointer to an array. In
+ the example above, the residual is a scalar, so only
+ ``residuals[0]`` is set.
+
+ Then given this class definition, the numerically differentiated
+ :class:`CostFunction` with central differences used for computing
+ the derivative can be constructed as follows.
+
+ .. code-block:: c++
+
+ CostFunction* cost_function
+ = new NumericDiffCostFunction<MyScalarCostFunctor, CENTRAL, 1, 2, 2>(
+ new MyScalarCostFunctor(1.0)); ^ ^ ^
+ | | | |
+ Finite Differencing Scheme -+ | | |
+ Dimension of residual ----------+ | |
+ Dimension of x --------------------+ |
+ Dimension of y -----------------------+
+
+ In this example, there is usually an instance for each measumerent of `k`.
+
+ In the instantiation above, the template parameters following
+ ``MyScalarCostFunctor``, ``1, 2, 2``, describe the functor as
+ computing a 1-dimensional output from two arguments, both
+ 2-dimensional.
+
+ The framework can currently accommodate cost functions of up to 10
+ independent variables, and there is no limit on the dimensionality
+ of each of them.
+
+ The ``CENTRAL`` difference method is considerably more accurate at
+ the cost of twice as many function evaluations than forward
+ difference. Consider using central differences begin with, and only
+ after that works, trying forward difference to improve performance.
+
+ **WARNING** A common beginner's error when first using
+ 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.
+
+
+ **Alternate Interface**
+
+ For a variety of reason, 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 nececssary 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:
+
+ .. code-block:: c++
+
+ CostFunction* cost_function
+ = new NumericDiffCostFunction<MyCostFunction, CENTRAL, 1, 4, 8>(
+ new MyCostFunction(...), TAKE_OWNERSHIP);
+
+ 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:`NormalPrior`
+--------------------
+
+.. class:: NormalPrior
+
+ .. code-block:: c++
+
+ class NormalPrior: public CostFunction {
+ public:
+ // Check that the number of rows in the vector b are the same as the
+ // number of columns in the matrix A, crash otherwise.
+ NormalPrior(const Matrix& A, const Vector& b);
+
+ virtual bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const;
+ };
+
+ Implements a cost function of the form
+
+ .. math:: cost(x) = ||A(x - b)||^2
+
+ where, the matrix A and the vector b are fixed and x 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)
+
+ 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.
+
+
+:class:`ConditionedCostFunction`
+--------------------------------
+
+.. class:: ConditionedCostFunction
+
+ This class allows you to apply different conditioning to the residual
+ values of a wrapped cost function. An example where this is useful is
+ where you have an existing cost function that produces N values, but you
+ want the total cost to be something other than just the sum of these
+ squared values - maybe you want to apply a different scaling to some
+ values, to change their contribution to the cost.
+
+ Usage:
+
+ .. code-block:: c++
+
+ // my_cost_function produces N residuals
+ CostFunction* my_cost_function = ...
+ CHECK_EQ(N, my_cost_function->num_residuals());
+ vector<CostFunction*> conditioners;
+
+ // Make N 1x1 cost functions (1 parameter, 1 residual)
+ CostFunction* f_1 = ...
+ conditioners.push_back(f_1);
+
+ CostFunction* f_N = ...
+ conditioners.push_back(f_N);
+ ConditionedCostFunction* ccf =
+ new ConditionedCostFunction(my_cost_function, conditioners);
+
+
+ Now ``ccf`` 's ``residual[i]`` (i=0..N-1) will be passed though the
+ :math:`i^{\text{th}}` conditioner.
+
+ .. code-block:: c++
+
+ ccf_residual[i] = f_i(my_cost_function_residual[i])
+
+ and the Jacobian will be affected appropriately.
+
+:class:`CostFunctionToFunctor`
+------------------------------
+
+.. class:: CostFunctionToFunctor
+
+ :class:`CostFunctionToFunctor` is an adapter class that allows users to use
+ :class:`CostFunction` objects in templated functors which are to be used for
+ automatic differentiation. This allows the user to seamlessly mix
+ analytic, numeric and automatic differentiation.
+
+ For example, let us assume that
+
+ .. code-block:: c++
+
+ class IntrinsicProjection : public SizedCostFunction<2, 5, 3> {
+ public:
+ IntrinsicProjection(const double* observations);
+ virtual bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const;
+ };
+
+ is a :class:`CostFunction` that implements the projection of a
+ point in its local coordinate system onto its image plane and
+ subtracts it from the observed point projection. It can compute its
+ residual and either via analytic or numerical differentiation can
+ compute its jacobians.
+
+ Now we would like to compose the action of this
+ :class:`CostFunction` with the action of camera extrinsics, i.e.,
+ rotation and translation. Say we have a templated function
+
+ .. code-block:: c++
+
+ template<typename T>
+ void RotateAndTranslatePoint(const T* rotation,
+ const T* translation,
+ const T* point,
+ T* result);
+
+
+ Then we can now do the following,
+
+ .. code-block:: c++
+
+ struct CameraProjection {
+ CameraProjection(double* observation) {
+ intrinsic_projection_.reset(
+ new CostFunctionToFunctor<2, 5, 3>(new IntrinsicProjection(observation_)));
+ }
+ template <typename T>
+ bool operator(const T* rotation,
+ const T* translation,
+ const T* intrinsics,
+ const T* point,
+ T* residual) const {
+ T transformed_point[3];
+ RotateAndTranslatePoint(rotation, translation, point, transformed_point);
+
+ // Note that we call intrinsic_projection_, just like it was
+ // any other templated functor.
+ return (*intrinsic_projection_)(intrinsics, transformed_point, residual);
+ }
+
+ private:
+ scoped_ptr<CostFunctionToFunctor<2,5,3> > intrinsic_projection_;
+ };
+
+
+:class:`NumericDiffFunctor`
+---------------------------
+
+.. class:: NumericDiffFunctor
+
+ A wrapper class that takes a variadic functor evaluating a
+ function, numerically differentiates it and makes it available as a
+ templated functor so that it can be easily used as part of Ceres'
+ automatic differentiation framework.
+
+ For example, let us assume that
+
+ .. code-block:: c++
+
+ struct IntrinsicProjection
+ IntrinsicProjection(const double* observations);
+ bool operator()(const double* calibration,
+ const double* point,
+ double* residuals);
+ };
+
+ is a functor that implements the projection of a point in its local
+ coordinate system onto its image plane and subtracts it from the
+ observed point projection.
+
+ Now we would like to compose the action of this functor with the
+ action of camera extrinsics, i.e., rotation and translation, which
+ is given by the following templated function
+
+ .. code-block:: c++
+
+ template<typename T>
+ void RotateAndTranslatePoint(const T* rotation,
+ const T* translation,
+ const T* point,
+ T* result);
+
+ To compose the extrinsics and intrinsics, we can construct a
+ ``CameraProjection`` functor as follows.
+
+ .. code-block:: c++
+
+ struct CameraProjection {
+ typedef NumericDiffFunctor<IntrinsicProjection, CENTRAL, 2, 5, 3>
+ IntrinsicProjectionFunctor;
+
+ CameraProjection(double* observation) {
+ intrinsic_projection_.reset(
+ new IntrinsicProjectionFunctor(observation)) {
+ }
+
+ 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, residual);
+ }
+
+ private:
+ scoped_ptr<IntrinsicProjectionFunctor> intrinsic_projection_;
+ };
+
+ Here, we made the choice of using ``CENTRAL`` differences to compute
+ the jacobian of ``IntrinsicProjection``.
+
+ Now, we are ready to construct an automatically differentiated cost
+ function as
+
+ .. code-block:: c++
+
+ CostFunction* cost_function =
+ new AutoDiffCostFunction<CameraProjection, 2, 3, 3, 5>(
+ new CameraProjection(observations));
+
+ ``cost_function`` now seamlessly integrates automatic
+ differentiation of ``RotateAndTranslatePoint`` with a numerically
+ differentiated version of ``IntrinsicProjection``.
+
+
+:class:`LossFunction`
+---------------------
+
+.. class:: LossFunction
+
+ For least squares problems where the minimization may encounter
+ input terms that contain outliers, that is, completely bogus
+ measurements, it is important to use a loss function that reduces
+ their influence.
+
+ Consider a structure from motion problem. The unknowns are 3D
+ points and camera parameters, and the measurements are image
+ coordinates describing the expected reprojected position for a
+ point in a camera. For example, we want to model the geometry of a
+ street scene with fire hydrants and cars, observed by a moving
+ camera with unknown parameters, and the only 3D points we care
+ about are the pointy tippy-tops of the fire hydrants. Our magic
+ image processing algorithm, which is responsible for producing the
+ measurements that are input to Ceres, has found and matched all
+ such tippy-tops in all image frames, except that in one of the
+ frame it mistook a car's headlight for a hydrant. If we didn't do
+ anything special the residual for the erroneous measurement will
+ result in the entire solution getting pulled away from the optimum
+ to reduce the large error that would otherwise be attributed to the
+ wrong measurement.
+
+ Using a robust loss function, the cost for large residuals is
+ reduced. In the example above, this leads to outlier terms getting
+ down-weighted so they do not overly influence the final solution.
+
+ .. code-block:: c++
+
+ class LossFunction {
+ public:
+ virtual void Evaluate(double s, double out[3]) const = 0;
+ };
+
+
+ The key method is :func:`LossFunction::Evaluate`, which given a
+ non-negative scalar ``s``, computes
+
+ .. math:: out = \begin{bmatrix}\rho(s), & \rho'(s), & \rho''(s)\end{bmatrix}
+
+ Here the convention is that the contribution of a term to the cost
+ function is given by :math:`\frac{1}{2}\rho(s)`, where :math:`s
+ =\|f_i\|^2`. Calling the method with a negative value of :math:`s`
+ is an error and the implementations are not required to handle that
+ case.
+
+ Most sane choices of :math:`\rho` satisfy:
+
+ .. math::
+
+ \rho(0) &= 0\\
+ \rho'(0) &= 1\\
+ \rho'(s) &< 1 \text{ in the outlier region}\\
+ \rho''(s) &< 0 \text{ in the outlier region}
+
+ so that they mimic the squared cost for small residuals.
+
+ **Scaling**
+
+
+ Given one robustifier :math:`\rho(s)` one can change the length
+ scale at which robustification takes place, by adding a scale
+ factor :math:`a > 0` which gives us :math:`\rho(s,a) = a^2 \rho(s /
+ a^2)` and the first and second derivatives as :math:`\rho'(s /
+ a^2)` and :math:`(1 / a^2) \rho''(s / a^2)` respectively.
+
+
+ The reason for the appearance of squaring is that :math:`a` is in
+ the units of the residual vector norm whereas :math:`s` is a squared
+ norm. For applications it is more convenient to specify :math:`a` than
+ its square.
+
+Instances
+^^^^^^^^^
+
+Ceres includes a number of other loss functions. For simplicity we
+described their unscaled versions. The figure below illustrates their
+shape graphically. More details can be found in
+``include/ceres/loss_function.h``.
+
+.. figure:: loss.png
+ :figwidth: 500px
+ :height: 400px
+ :align: center
+
+ Shape of the various common loss functions.
+
+.. class:: TrivialLoss
+
+ .. math:: \rho(s) = s
+
+.. class:: HuberLoss
+
+ .. math:: \rho(s) = \begin{cases} s & s \le 1\\ 2 \sqrt{s} - 1 & s > 1 \end{cases}
+
+.. class:: SoftLOneLoss
+
+ .. math:: \rho(s) = 2 (\sqrt{1+s} - 1)
+
+.. class:: CauchyLoss
+
+ .. math:: \rho(s) = \log(1 + s)
+
+.. class:: ArctanLoss
+
+ .. math:: \rho(s) = \arctan(s)
+
+.. class:: TolerantLoss
+
+ .. math:: \rho(s,a,b) = b \log(1 + e^{(s - a) / b}) - b \log(1 + e^{-a / b})
+
+.. class:: ComposedLoss
+
+.. class:: ScaledLoss
+
+.. class:: LossFunctionWrapper
+
+
+Theory
+^^^^^^
+
+Let us consider a problem with a single problem and a single parameter
+block.
+
+.. math::
+
+ \min_x \frac{1}{2}\rho(f^2(x))
+
+
+Then, the robustified gradient and the Gauss-Newton Hessian are
+
+.. math::
+
+ g(x) &= \rho'J^\top(x)f(x)\\
+ H(x) &= J^\top(x)\left(\rho' + 2 \rho''f(x)f^\top(x)\right)J(x)
+
+where the terms involving the second derivatives of :math:`f(x)` have
+been ignored. Note that :math:`H(x)` is indefinite if
+:math:`\rho''f(x)^\top f(x) + \frac{1}{2}\rho' < 0`. If this is not
+the case, then its possible to re-weight the residual and the Jacobian
+matrix such that the corresponding linear least squares problem for
+the robustified Gauss-Newton step.
+
+
+Let :math:`\alpha` be a root of
+
+.. math:: \frac{1}{2}\alpha^2 - \alpha - \frac{\rho''}{\rho'}\|f(x)\|^2 = 0.
+
+
+Then, define the rescaled residual and Jacobian as
+
+.. math::
+
+ \tilde{f}(x) &= \frac{\sqrt{\rho'}}{1 - \alpha} f(x)\\
+ \tilde{J}(x) &= \sqrt{\rho'}\left(1 - \alpha
+ \frac{f(x)f^\top(x)}{\left\|f(x)\right\|^2} \right)J(x)
+
+
+In the case :math:`2 \rho''\left\|f(x)\right\|^2 + \rho' \lesssim 0`,
+we limit :math:`\alpha \le 1- \epsilon` for some small
+:math:`\epsilon`. For more details see [Triggs]_.
+
+With this simple rescaling, one can use any Jacobian based non-linear
+least squares algorithm to robustifed non-linear least squares
+problems.
+
+
+:class:`LocalParameterization`
+------------------------------
+
+.. class:: LocalParameterization
+
+ .. code-block:: c++
+
+ class LocalParameterization {
+ public:
+ virtual ~LocalParameterization() {}
+ virtual bool Plus(const double* x,
+ const double* delta,
+ double* x_plus_delta) const = 0;
+ virtual bool ComputeJacobian(const double* x, double* jacobian) const = 0;
+ virtual int GlobalSize() const = 0;
+ virtual int LocalSize() const = 0;
+ };
+
+ Sometimes the parameters :math:`x` can overparameterize a
+ problem. In that case it is desirable to choose a parameterization
+ to remove the null directions of the cost. More generally, if
+ :math:`x` lies on a manifold of a smaller dimension than the
+ ambient space that it is embedded in, then it is numerically and
+ computationally more effective to optimize it using a
+ parameterization that lives in the tangent space of that manifold
+ at each point.
+
+ For example, a sphere in three dimensions is a two dimensional
+ manifold, embedded in a three dimensional space. At each point on
+ the sphere, the plane tangent to it defines a two dimensional
+ tangent space. For a cost function defined on this sphere, given a
+ point :math:`x`, moving in the direction normal to the sphere at
+ that point is not useful. Thus a better way to parameterize a point
+ on a sphere is to optimize over two dimensional vector
+ :math:`\Delta x` in the tangent space at the point on the sphere
+ point and then "move" to the point :math:`x + \Delta x`, where the
+ move operation involves projecting back onto the sphere. Doing so
+ removes a redundant dimension from the optimization, making it
+ numerically more robust and efficient.
+
+ More generally we can define a function
+
+ .. math:: x' = \boxplus(x, \Delta x),
+
+ where :math:`x` has the same size as :math:`x`, and :math:`\Delta
+ x` is of size less than or equal to :math:`x`. The function
+ :math:`\boxplus`, generalizes the definition of vector
+ addition. Thus it satisfies the identity
+
+ .. math:: \boxplus(x, 0) = x,\quad \forall x.
+
+ Instances of :class:`LocalParameterization` implement the
+ :math:`\boxplus` operation and its derivative with respect to
+ :math:`\Delta x` at :math:`\Delta x = 0`.
+
+
+.. function:: int LocalParameterization::GlobalSize()
+
+ The dimension of the ambient space in which the parameter block
+ :math:`x` lives.
+
+.. function:: int LocalParamterization::LocaLocalSize()
+
+ The size of the tangent space
+ that :math:`\Delta x` lives in.
+
+.. function:: bool LocalParameterization::Plus(const double* x, const double* delta, double* x_plus_delta) const
+
+ :func:`LocalParameterization::Plus` implements :math:`\boxplus(x,\Delta x)`.
+
+.. function:: bool LocalParameterization::ComputeJacobian(const double* x, double* jacobian) const
+
+ Computes the Jacobian matrix
+
+ .. math:: J = \left . \frac{\partial }{\partial \Delta x} \boxplus(x,\Delta x)\right|_{\Delta x = 0}
+
+ in row major form.
+
+Instances
+^^^^^^^^^
+
+.. class:: IdentityParameterization
+
+ A trivial version of :math:`\boxplus` is when :math:`\Delta x` is
+ of the same size as :math:`x` and
+
+ .. math:: \boxplus(x, \Delta x) = x + \Delta x
+
+.. class:: SubsetParameterization
+
+ A more interesting case if :math:`x` is a two dimensional vector,
+ and the user wishes to hold the first coordinate constant. Then,
+ :math:`\Delta x` is a scalar and :math:`\boxplus` is defined as
+
+ .. math::
+
+ \boxplus(x, \Delta x) = x + \left[ \begin{array}{c} 0 \\ 1
+ \end{array} \right] \Delta x
+
+ :class:`SubsetParameterization` generalizes this construction to
+ hold any part of a parameter block constant.
+
+.. class:: QuaternionParameterization
+
+ Another example that occurs commonly in Structure from Motion
+ problems is when camera rotations are parameterized using a
+ quaternion. There, it is useful only to make updates orthogonal to
+ that 4-vector defining the quaternion. One way to do this is to let
+ :math:`\Delta x` be a 3 dimensional vector and define
+ :math:`\boxplus` to be
+
+ .. math:: \boxplus(x, \Delta x) = \left[ \cos(|\Delta x|), \frac{\sin\left(|\Delta x|\right)}{|\Delta x|} \Delta x \right] * x
+ :label: quaternion
+
+ The multiplication between the two 4-vectors on the right hand side
+ is the standard quaternion
+ product. :class:`QuaternionParameterization` is an implementation
+ of :eq:`quaternion`.
+
+
+:class:`Problem`
+----------------
+
+.. class:: Problem
+
+ :class:`Problem` holds the robustified non-linear least squares
+ problem :eq:`ceresproblem`. To create a least squares problem, use
+ the :func:`Problem::AddResidualBlock` and
+ :func:`Problem::AddParameterBlock` methods.
+
+ For example a problem containing 3 parameter blocks of sizes 3, 4
+ and 5 respectively and two residual blocks of size 2 and 6:
+
+ .. code-block:: c++
+
+ double x1[] = { 1.0, 2.0, 3.0 };
+ double x2[] = { 1.0, 2.0, 3.0, 5.0 };
+ double x3[] = { 1.0, 2.0, 3.0, 6.0, 7.0 };
+
+ Problem problem;
+ problem.AddResidualBlock(new MyUnaryCostFunction(...), x1);
+ problem.AddResidualBlock(new MyBinaryCostFunction(...), x2, x3);
+
+ :func:`Problem::AddResidualBlock` as the name implies, adds a
+ residual block to the problem. It adds a :class:`CostFunction` , an
+ optional :class:`LossFunction` and connects the
+ :class:`CostFunction` to a set of parameter block.
+
+ The cost function carries with it information about the sizes of
+ the parameter blocks it expects. The function checks that these
+ match the sizes of the parameter blocks listed in
+ ``parameter_blocks``. The program aborts if a mismatch is
+ detected. ``loss_function`` can be ``NULL``, in which case the cost
+ of the term is just the squared norm of the residuals.
+
+ The user has the option of explicitly adding the parameter blocks
+ using :func:`Problem::AddParameterBlock`. This causes additional correctness
+ checking; however, :func:`Problem::AddResidualBlock` implicitly adds the
+ parameter blocks if they are not present, so calling
+ :func:`Problem::AddParameterBlock` explicitly is not required.
+
+
+ :class:`Problem` by default takes ownership of the ``cost_function`` and
+ ``loss_function`` pointers. These objects remain live for the life of
+ the :class:`Problem` object. If the user wishes to keep control over the
+ destruction of these objects, then they can do this by setting the
+ corresponding enums in the ``Problem::Options`` struct.
+
+
+ Note that even though the Problem takes ownership of ``cost_function``
+ and ``loss_function``, it does not preclude the user from re-using
+ them in another residual block. The destructor takes care to call
+ delete on each ``cost_function`` or ``loss_function`` pointer only
+ once, regardless of how many residual blocks refer to them.
+
+ :func:`Problem::AddParameterBlock` explicitly adds a parameter
+ block to the :class:`Problem`. Optionally it allows the user to
+ associate a :class:`LocalParameterization` object with the parameter
+ block too. Repeated calls with the same arguments are
+ ignored. Repeated calls with the same double pointer but a
+ different size results in undefined behaviour.
+
+ You can set any parameter block to be constant using
+ :func:`Problem::SetParameterBlockConstant` and undo this using
+ :func:`SetParameterBlockVariable`.
+
+ In fact you can set any number of parameter blocks to be constant,
+ and Ceres is smart enough to figure out what part of the problem
+ you have constructed depends on the parameter blocks that are free
+ to change and only spends time solving it. So for example if you
+ constructed a problem with a million parameter blocks and 2 million
+ residual blocks, but then set all but one parameter blocks to be
+ constant and say only 10 residual blocks depend on this one
+ non-constant parameter block. Then the computational effort Ceres
+ spends in solving this problem will be the same if you had defined
+ a problem with one parameter block and 10 residual blocks.
+
+ **Ownership**
+
+ :class:`Problem` by default takes ownership of the
+ ``cost_function``, ``loss_function`` and ``local_parameterization``
+ pointers. These objects remain live for the life of the
+ :class:`Problem`. If the user wishes to keep control over the
+ destruction of these objects, then they can do this by setting the
+ corresponding enums in the :class:`Problem::Options` struct.
+
+ Even though :class:`Problem` takes ownership of these pointers, it
+ does not preclude the user from re-using them in another residual
+ or parameter block. The destructor takes care to call delete on
+ each pointer only once.
+
+
+.. function:: ResidualBlockId Problem::AddResidualBlock(CostFunction* cost_function, LossFunction* loss_function, const vector<double*> parameter_blocks)
+
+
+.. function:: void Problem::AddParameterBlock(double* values, int size, LocalParameterization* local_parameterization)
+ void Problem::AddParameterBlock(double* values, int size)
+
+
+.. function:: void Problem::SetParameterBlockConstant(double* values)
+
+.. function:: void Problem::SetParameterBlockVariable(double* values)
+
+
+.. function:: void Problem::SetParameterization(double* values, LocalParameterization* local_parameterization)
+
+
+.. function:: int Problem::NumParameterBlocks() const
+
+
+.. function:: int Problem::NumParameters() const
+
+
+.. function:: int Problem::NumResidualBlocks() const
+
+
+.. function:: int Problem::NumResiduals() const
+
+
+
+``rotation.h``
+--------------
+
+Many applications of Ceres Solver involve optimization problems where
+some of the variables correspond to rotations. To ease the pain of
+work with the various representations of rotations (angle-axis,
+quaternion and matrix) we provide a handy set of templated
+functions. These functions are templated so that the user can use them
+within Ceres Solver's automatic differentiation framework.
+
+.. function:: template<typename T> void AngleAxisToQuaternion(T const* angle_axis, T* quaternion);
+
+ Convert a value in combined axis-angle representation to a
+ quaternion.
+
+ The value ``angle_axis`` is a triple whose norm is an angle in radians,
+ and whose direction is aligned with the axis of rotation, and
+ ``quaternion`` is a 4-tuple that will contain the resulting quaternion.
+
+.. function:: template<typename T> void QuaternionToAngleAxis(T const* quaternion, T* angle_axis);
+
+ Convert a quaternion to the equivalent combined axis-angle
+ representation.
+
+ The value ``quaternion`` must be a unit quaternion - it is not
+ normalized first, and ``angle_axis`` will be filled with a value
+ whose norm is the angle of rotation in radians, and whose direction
+ is the axis of rotation.
+
+.. function:: template <typename T> void RotationMatrixToAngleAxis(T const * R, T * angle_axis);
+.. function:: template <typename T> void AngleAxisToRotationMatrix(T const * angle_axis, T * R);
+
+ Conversions between 3x3 rotation matrix (in column major order) and
+ axis-angle rotation representations.
+
+.. function:: template <typename T> void EulerAnglesToRotationMatrix(const T* euler, int row_stride, T* R);
+
+ Conversions between 3x3 rotation matrix (in row major order) and
+ Euler angle (in degrees) rotation representations.
+
+ The {pitch,roll,yaw} Euler angles are rotations around the {x,y,z}
+ axes, respectively. They are applied in that same order, so the
+ total rotation R is Rz * Ry * Rx.
+
+.. function:: template <typename T> inline void QuaternionToScaledRotation(const T q[4], T R[3 * 3]);
+
+ Convert a 4-vector to a 3x3 scaled rotation matrix.
+
+ The choice of rotation is such that the quaternion
+ :math:`\begin{bmatrix} 1 &0 &0 &0\end{bmatrix}` goes to an identity
+ matrix and for small :math:`a, b, c` the quaternion
+ :math:`\begin{bmatrix}1 &a &b &c\end{bmatrix}` goes to the matrix
+
+ .. math::
+
+ I + 2 \begin{bmatrix} 0 & -c & b \\ c & 0 & -a\\ -b & a & 0
+ \end{bmatrix} + O(q^2)
+
+ which corresponds to a Rodrigues approximation, the last matrix
+ being the cross-product matrix of :math:`\begin{bmatrix} a& b&
+ c\end{bmatrix}`. Together with the property that :math:`R(q1 * q2)
+ = R(q1) * R(q2)` this uniquely defines the mapping from :math:`q` to
+ :math:`R`.
+
+ The rotation matrix ``R`` is row-major.
+
+ No normalization of the quaternion is performed, i.e.
+ :math:`R = \|q\|^2 Q`, where :math:`Q` is an orthonormal matrix
+ such that :math:`\det(Q) = 1` and :math:`Q*Q' = I`.
+
+
+.. function:: template <typename T> inline void QuaternionToRotation(const T q[4], T R[3 * 3]);
+
+ Same as above except that the rotation matrix is normalized by the
+ Frobenius norm, so that :math:`R R' = I` (and :math:`\det(R) = 1`).
+
+.. function:: template <typename T> inline void UnitQuaternionRotatePoint(const T q[4], const T pt[3], T result[3]);
+
+ Rotates a point pt by a quaternion q:
+
+ .. math:: \text{result} = R(q) \text{pt}
+
+ Assumes the quaternion is unit norm. If you pass in a quaternion
+ with :math:`|q|^2 = 2` then you WILL NOT get back 2 times the
+ result you get for a unit quaternion.
+
+
+.. function:: template <typename T> inline void QuaternionRotatePoint(const T q[4], const T pt[3], T result[3]);
+
+ With this function you do not need to assume that q has unit norm.
+ It does assume that the norm is non-zero.
+
+.. function:: template<typename T> inline void QuaternionProduct(const T z[4], const T w[4], T zw[4]);
+
+ .. math:: zw = z * w
+
+ where :math`*` is the Quaternion product between 4-vectors.
+
+
+.. function:: template<typename T> inline void CrossProduct(const T x[3], const T y[3], T x_cross_y[3]);
+
+ .. math:: \text{x_cross_y} = x \times y
+
+.. function:: template<typename T> inline void AngleAxisRotatePoint(const T angle_axis[3], const T pt[3], T result[3]);
+
+ .. math:: y = R(\text{angle_axis}) x
diff --git a/docs/source/reading.rst b/docs/source/reading.rst
new file mode 100644
index 0000000..0612aa0
--- /dev/null
+++ b/docs/source/reading.rst
@@ -0,0 +1,10 @@
+===============
+Further Reading
+===============
+
+For a short but informative introduction to the subject we recommend
+the booklet by [Madsen]_ . For a general introduction to non-linear
+optimization we recommend [NocedalWright]_. [Bjorck]_ remains the
+seminal reference on least squares problems. [TrefethenBau]_ book is
+our favourite text on introductory numerical linear algebra. [Triggs]_
+provides a thorough coverage of the bundle adjustment problem.
diff --git a/docs/source/solving.rst b/docs/source/solving.rst
new file mode 100644
index 0000000..43eedee
--- /dev/null
+++ b/docs/source/solving.rst
@@ -0,0 +1,1203 @@
+
+.. default-domain:: cpp
+
+.. cpp:namespace:: ceres
+
+.. _chapter-solving:
+
+==========
+Solver API
+==========
+
+Effective use of Ceres requires some familiarity with the basic
+components of a nonlinear least squares solver, so before we describe
+how to configure the solver, we will begin by taking a brief look at
+how some of the core optimization algorithms in Ceres work and the
+various linear solvers and preconditioners that power it.
+
+.. _section-trust-region-methods:
+
+Trust Region Methods
+--------------------
+
+Let :math:`x \in \mathbb{R}^n` be an :math:`n`-dimensional vector of
+variables, and
+:math:`F(x) = \left[f_1(x), ... , f_{m}(x) \right]^{\top}` be a
+:math:`m`-dimensional function of :math:`x`. We are interested in
+solving the following optimization problem [#f1]_ .
+
+.. math:: \arg \min_x \frac{1}{2}\|F(x)\|^2\ .
+ :label: nonlinsq
+
+Here, the Jacobian :math:`J(x)` of :math:`F(x)` is an :math:`m\times
+n` matrix, where :math:`J_{ij}(x) = \partial_j f_i(x)` and the
+gradient vector :math:`g(x) = \nabla \frac{1}{2}\|F(x)\|^2 = J(x)^\top
+F(x)`. Since the efficient global minimization of :eq:`nonlinsq` for general
+:math:`F(x)` is an intractable problem, we will have to settle for
+finding a local minimum.
+
+
+The general strategy when solving non-linear optimization problems is
+to solve a sequence of approximations to the original problem
+[NocedalWright]_. At each iteration, the approximation is solved to
+determine a correction :math:`\Delta x` to the vector :math:`x`. For
+non-linear least squares, an approximation can be constructed by using
+the linearization :math:`F(x+\Delta x) \approx F(x) + J(x)\Delta x`,
+which leads to the following linear least squares problem:
+
+.. math:: \min_{\Delta x} \frac{1}{2}\|J(x)\Delta x + F(x)\|^2
+ :label: linearapprox
+
+Unfortunately, naively solving a sequence of these problems and
+updating :math:`x \leftarrow x+ \Delta x` leads to an algorithm that may not
+converge. To get a convergent algorithm, we need to control the size
+of the step :math:`\Delta x`. And this is where the idea of a trust-region
+comes in.
+
+.. Algorithm~\ref{alg:trust-region} describes the basic trust-region
+.. loop for non-linear least squares problems.
+
+.. \begin{algorithm} \caption{The basic trust-region
+ algorithm.\label{alg:trust-region}} \begin{algorithmic} \REQUIRE
+ Initial point `x` and a trust region radius `\mu`. \LOOP
+ \STATE{Solve `\arg \min_{\Delta x} \frac{1}{2}\|J(x)\Delta x +
+ F(x)\|^2` s.t. `\|D(x)\Delta x\|^2 \le \mu`} \STATE{`\rho =
+ \frac{\displaystyle \|F(x + \Delta x)\|^2 -
+ \|F(x)\|^2}{\displaystyle \|J(x)\Delta x + F(x)\|^2 - \|F(x)\|^2}`}
+ \IF {`\rho > \epsilon`} \STATE{`x = x + \Delta x`} \ENDIF \IF {`\rho
+ > \eta_1`} \STATE{`\rho = 2 * \rho`} \ELSE \IF {`\rho < \eta_2`}
+ \STATE {`\rho = 0.5 * \rho`} \ENDIF \ENDIF \ENDLOOP
+ \end{algorithmic} \end{algorithm}
+
+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
+:math:`\rho` measures the quality of the step :math:`\Delta x`, i.e.,
+how well did the linear model predict the decrease in the value of the
+non-linear objective. The idea is to increase or decrease the radius
+of the trust region depending on how well the linearization predicts
+the behavior of the non-linear objective, which in turn is reflected
+in the value of :math:`\rho`.
+
+The key computational step in a trust-region algorithm is the solution
+of the constrained optimization problem
+
+.. math:: \arg\min_{\Delta x} \frac{1}{2}\|J(x)\Delta x + F(x)\|^2\quad \text{such that}\quad \|D(x)\Delta x\|^2 \le \mu
+ :label: trp
+
+There are a number of different ways of solving this problem, each
+giving rise to a different concrete trust-region algorithm. Currently
+Ceres, implements two trust-region algorithms - Levenberg-Marquardt
+and Dogleg. The user can choose between them by setting
+:member:`Solver::Options::trust_region_strategy_type`.
+
+.. rubric:: Footnotes
+
+.. [#f1] At the level of the non-linear solver, the block and
+ structure is not relevant, therefore our discussion here is
+ in terms of an optimization problem defined over a state
+ vector of size :math:`n`.
+
+.. _section-levenberg-marquardt:
+
+Levenberg-Marquardt
+^^^^^^^^^^^^^^^^^^^
+
+The Levenberg-Marquardt algorithm [Levenberg]_ [Marquardt]_ is the
+most popular algorithm for solving non-linear least squares problems.
+It was also the first trust region algorithm to be developed
+[Levenberg]_ [Marquardt]_. Ceres implements an exact step [Madsen]_
+and an inexact step variant of the Levenberg-Marquardt algorithm
+[WrightHolt]_ [NashSofer]_.
+
+It can be shown, that the solution to :eq:`trp` can be obtained by
+solving an unconstrained optimization of the form
+
+.. math:: \arg\min_{\Delta x}& \frac{1}{2}\|J(x)\Delta x + F(x)\|^2 +\lambda \|D(x)\Delta x\|^2
+
+Where, :math:`\lambda` is a Lagrange multiplier that is inverse
+related to :math:`\mu`. In Ceres, we solve for
+
+.. math:: \arg\min_{\Delta x}& \frac{1}{2}\|J(x)\Delta x + F(x)\|^2 + \frac{1}{\mu} \|D(x)\Delta x\|^2
+ :label: lsqr
+
+The matrix :math:`D(x)` is a non-negative diagonal matrix, typically
+the square root of the diagonal of the matrix :math:`J(x)^\top J(x)`.
+
+Before going further, let us make some notational simplifications. We
+will assume that the matrix :math:`\sqrt{\mu} D` has been concatenated
+at the bottom of the matrix :math:`J` and similarly a vector of zeros
+has been added to the bottom of the vector :math:`f` and the rest of
+our discussion will be in terms of :math:`J` and :math:`f`, i.e, the
+linear least squares problem.
+
+.. math:: \min_{\Delta x} \frac{1}{2} \|J(x)\Delta x + f(x)\|^2 .
+ :label: simple
+
+For all but the smallest problems the solution of :eq:`simple` in
+each iteration of the Levenberg-Marquardt algorithm is the dominant
+computational cost in Ceres. Ceres provides a number of different
+options for solving :eq:`simple`. There are two major classes of
+methods - factorization and iterative.
+
+The factorization methods are based on computing an exact solution of
+:eq:`lsqr` using a Cholesky or a QR factorization and lead to an exact
+step Levenberg-Marquardt algorithm. But it is not clear if an exact
+solution of :eq:`lsqr` is necessary at each step of the LM algorithm
+to solve :eq:`nonlinsq`. In fact, we have already seen evidence
+that this may not be the case, as :eq:`lsqr` is itself a regularized
+version of :eq:`linearapprox`. Indeed, it is possible to
+construct non-linear optimization algorithms in which the linearized
+problem is solved approximately. These algorithms are known as inexact
+Newton or truncated Newton methods [NocedalWright]_.
+
+An inexact Newton method requires two ingredients. First, a cheap
+method for approximately solving systems of linear
+equations. Typically an iterative linear solver like the Conjugate
+Gradients method is used for this
+purpose [NocedalWright]_. Second, a termination rule for
+the iterative solver. A typical termination rule is of the form
+
+.. math:: \|H(x) \Delta x + g(x)\| \leq \eta_k \|g(x)\|.
+ :label: inexact
+
+Here, :math:`k` indicates the Levenberg-Marquardt iteration number and
+:math:`0 < \eta_k <1` is known as the forcing sequence. [WrightHolt]_
+prove that a truncated Levenberg-Marquardt algorithm that uses an
+inexact Newton step based on :eq:`inexact` converges for any
+sequence :math:`\eta_k \leq \eta_0 < 1` and the rate of convergence
+depends on the choice of the forcing sequence :math:`\eta_k`.
+
+Ceres supports both exact and inexact step solution strategies. When
+the user chooses a factorization based linear solver, the exact step
+Levenberg-Marquardt algorithm is used. When the user chooses an
+iterative linear solver, the inexact step Levenberg-Marquardt
+algorithm is used.
+
+.. _section-dogleg:
+
+Dogleg
+^^^^^^
+
+Another strategy for solving the trust region problem :eq:`trp` was
+introduced by M. J. D. Powell. The key idea there is to compute two
+vectors
+
+.. math::
+
+ \Delta x^{\text{Gauss-Newton}} &= \arg \min_{\Delta x}\frac{1}{2} \|J(x)\Delta x + f(x)\|^2.\\
+ \Delta x^{\text{Cauchy}} &= -\frac{\|g(x)\|^2}{\|J(x)g(x)\|^2}g(x).
+
+Note that the vector :math:`\Delta x^{\text{Gauss-Newton}}` is the
+solution to :eq:`linearapprox` and :math:`\Delta
+x^{\text{Cauchy}}` is the vector that minimizes the linear
+approximation if we restrict ourselves to moving along the direction
+of the gradient. Dogleg methods finds a vector :math:`\Delta x`
+defined by :math:`\Delta x^{\text{Gauss-Newton}}` and :math:`\Delta
+x^{\text{Cauchy}}` that solves the trust region problem. Ceres
+supports two variants that can be chose by setting
+:member:`Solver::Options::dogleg_type`.
+
+``TRADITIONAL_DOGLEG`` as described by Powell, constructs two line
+segments using the Gauss-Newton and Cauchy vectors and finds the point
+farthest along this line shaped like a dogleg (hence the name) that is
+contained in the trust-region. For more details on the exact reasoning
+and computations, please see Madsen et al [Madsen]_.
+
+``SUBSPACE_DOGLEG`` is a more sophisticated method that considers the
+entire two dimensional subspace spanned by these two vectors and finds
+the point that minimizes the trust region problem in this
+subspace [Byrd]_.
+
+The key advantage of the Dogleg over Levenberg Marquardt is that if
+the step computation for a particular choice of :math:`\mu` does not
+result in sufficient decrease in the value of the objective function,
+Levenberg-Marquardt solves the linear approximation from scratch with
+a smaller value of :math:`\mu`. Dogleg on the other hand, only needs
+to compute the interpolation between the Gauss-Newton and the Cauchy
+vectors, as neither of them depend on the value of :math:`\mu`.
+
+The Dogleg method can only be used with the exact factorization based
+linear solvers.
+
+.. _section-inner-iterations:
+
+Inner Iterations
+^^^^^^^^^^^^^^^^
+
+Some non-linear least squares problems have additional structure in
+the way the parameter blocks interact that it is beneficial to modify
+the way the trust region step is computed. e.g., consider the
+following regression problem
+
+.. math:: y = a_1 e^{b_1 x} + a_2 e^{b_3 x^2 + c_1}
+
+
+Given a set of pairs :math:`\{(x_i, y_i)\}`, the user wishes to estimate
+:math:`a_1, a_2, b_1, b_2`, and :math:`c_1`.
+
+Notice that the expression on the left is linear in :math:`a_1` and
+:math:`a_2`, and given any value for :math:`b_1, b_2` and :math:`c_1`,
+it is possible to use linear regression to estimate the optimal values
+of :math:`a_1` and :math:`a_2`. It's possible to analytically
+eliminate the variables :math:`a_1` and :math:`a_2` from the problem
+entirely. Problems like these are known as separable least squares
+problem and the most famous algorithm for solving them is the Variable
+Projection algorithm invented by Golub & Pereyra [GolubPereyra]_.
+
+Similar structure can be found in the matrix factorization with
+missing data problem. There the corresponding algorithm is known as
+Wiberg's algorithm [Wiberg]_.
+
+Ruhe & Wedin present an analysis of various algorithms for solving
+separable non-linear least squares problems and refer to *Variable
+Projection* as Algorithm I in their paper [RuheWedin]_.
+
+Implementing Variable Projection is tedious and expensive. Ruhe &
+Wedin present a simpler algorithm with comparable convergence
+properties, which they call Algorithm II. Algorithm II performs an
+additional optimization step to estimate :math:`a_1` and :math:`a_2`
+exactly after computing a successful Newton step.
+
+
+This idea can be generalized to cases where the residual is not
+linear in :math:`a_1` and :math:`a_2`, i.e.,
+
+.. 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
+least squares solve. For non-linear problems, any method for solving
+the `a_1` and `a_2` optimization problems will do. The only constraint
+on `a_1` and `a_2` (if they are two different parameter block) is that
+they do not co-occur in a residual block.
+
+This idea can be further generalized, by not just optimizing
+:math:`(a_1, a_2)`, but decomposing the graph corresponding to the
+Hessian matrix's sparsity structure into a collection of
+non-overlapping independent sets and optimizing each of them.
+
+Setting :member:`Solver::Options::use_inner_iterations` to ``true``
+enables the use of this non-linear generalization of Ruhe & Wedin's
+Algorithm II. This version of Ceres has a higher iteration
+complexity, but also displays better convergence behavior per
+iteration.
+
+Setting :member:`Solver::Options::num_threads` to the maximum number
+possible is highly recommended.
+
+.. _section-non-monotonic-steps:
+
+Non-monotonic Steps
+^^^^^^^^^^^^^^^^^^^
+
+Note that the basic trust-region algorithm described in
+Algorithm~\ref{alg:trust-region} is a descent algorithm in that they
+only accepts a point if it strictly reduces the value of the objective
+function.
+
+Relaxing this requirement allows the algorithm to be more efficient in
+the long term at the cost of some local increase in the value of the
+objective function.
+
+This is because allowing for non-decreasing objective function values
+in a princpled manner allows the algorithm to *jump over boulders* as
+the method is not restricted to move into narrow valleys while
+preserving its convergence properties.
+
+Setting :member:`Solver::Options::use_nonmonotonic_steps` to ``true``
+enables the non-monotonic trust region algorithm as described by Conn,
+Gould & Toint in [Conn]_.
+
+Even though the value of the objective function may be larger
+than the minimum value encountered over the course of the
+optimization, the final parameters returned to the user are the
+ones corresponding to the minimum cost over all iterations.
+
+The option to take non-monotonic is available for all trust region
+strategies.
+
+
+.. _section-linear-solver:
+
+LinearSolver
+------------
+
+Recall that in both of the trust-region methods described above, the
+key computational cost is the solution of a linear least squares
+problem of the form
+
+.. math:: \min_{\Delta x} \frac{1}{2} \|J(x)\Delta x + f(x)\|^2 .
+ :label: simple2
+
+Let :math:`H(x)= J(x)^\top J(x)` and :math:`g(x) = -J(x)^\top
+f(x)`. For notational convenience let us also drop the dependence on
+:math:`x`. Then it is easy to see that solving :eq:`simple2` is
+equivalent to solving the *normal equations*.
+
+.. math:: H \Delta x = g
+ :label: normal
+
+Ceres provides a number of different options for solving :eq:`normal`.
+
+.. _section-qr:
+
+``DENSE_QR``
+^^^^^^^^^^^^
+
+For small problems (a couple of hundred parameters and a few thousand
+residuals) with relatively dense Jacobians, ``DENSE_QR`` is the method
+of choice [Bjorck]_. Let :math:`J = QR` be the QR-decomposition of
+:math:`J`, where :math:`Q` is an orthonormal matrix and :math:`R` is
+an upper triangular matrix [TrefethenBau]_. Then it can be shown that
+the solution to :eq:`normal` is given by
+
+.. math:: \Delta x^* = -R^{-1}Q^\top f
+
+
+Ceres uses ``Eigen`` 's dense QR factorization routines.
+
+.. _section-cholesky:
+
+``DENSE_NORMAL_CHOLESKY`` & ``SPARSE_NORMAL_CHOLESKY``
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Large non-linear least square problems are usually sparse. In such
+cases, using a dense QR factorization is inefficient. Let :math:`H =
+R^\top R` be the Cholesky factorization of the normal equations, where
+:math:`R` is an upper triangular matrix, then the solution to
+:eq:`normal` is given by
+
+.. math::
+
+ \Delta x^* = R^{-1} R^{-\top} g.
+
+
+The observant reader will note that the :math:`R` in the Cholesky
+factorization of :math:`H` is the same upper triangular matrix
+:math:`R` in the QR factorization of :math:`J`. Since :math:`Q` is an
+orthonormal matrix, :math:`J=QR` implies that :math:`J^\top J = R^\top
+Q^\top Q R = R^\top R`. There are two variants of Cholesky
+factorization -- sparse and dense.
+
+``DENSE_NORMAL_CHOLESKY`` as the name implies performs a dense
+Cholesky factorization of the normal equations. Ceres uses
+``Eigen`` 's dense LDLT factorization routines.
+
+``SPARSE_NORMAL_CHOLESKY``, as the name implies performs a sparse
+Cholesky factorization of the normal equations. This leads to
+substantial savings in time and memory for large sparse
+problems. Ceres uses the sparse Cholesky factorization routines in
+Professor Tim Davis' ``SuiteSparse`` or ``CXSparse`` packages [Chen]_.
+
+.. _section-schur:
+
+``DENSE_SCHUR`` & ``SPARSE_SCHUR``
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+While it is possible to use ``SPARSE_NORMAL_CHOLESKY`` to solve bundle
+adjustment problems, bundle adjustment problem have a special
+structure, and a more efficient scheme for solving :eq:`normal`
+can be constructed.
+
+Suppose that the SfM problem consists of :math:`p` cameras and
+:math:`q` points and the variable vector :math:`x` has the block
+structure :math:`x = [y_{1}, ... ,y_{p},z_{1}, ... ,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 impose any constancy
+requirement on these block sizes, but choosing them to be constant
+simplifies the exposition.
+
+A key characteristic of the bundle adjustment problem is that there is
+no term :math:`f_{i}` that includes two or more point blocks. This in
+turn implies that the matrix :math:`H` is of the form
+
+.. math:: H = \left[ \begin{matrix} B & E\\ E^\top & C \end{matrix} \right]\ ,
+ :label: hblock
+
+where, :math:`B \in \mathbb{R}^{pc\times pc}` is a block sparse matrix
+with :math:`p` blocks of size :math:`c\times c` and :math:`C \in
+\mathbb{R}^{qs\times qs}` is a block diagonal matrix with :math:`q` blocks
+of size :math:`s\times s`. :math:`E \in \mathbb{R}^{pc\times qs}` is a
+general block sparse matrix, with a block of size :math:`c\times s`
+for each observation. Let us now block partition :math:`\Delta x =
+[\Delta y,\Delta z]` and :math:`g=[v,w]` to restate :eq:`normal`
+as the block structured linear system
+
+.. math:: \left[ \begin{matrix} B & E\\ E^\top & C \end{matrix}
+ \right]\left[ \begin{matrix} \Delta y \\ \Delta z
+ \end{matrix} \right] = \left[ \begin{matrix} v\\ w
+ \end{matrix} \right]\ ,
+ :label: linear2
+
+and apply Gaussian elimination to it. As we noted above, :math:`C` is
+a block diagonal matrix, with small diagonal blocks of size
+:math:`s\times s`. Thus, calculating the inverse of :math:`C` by
+inverting each of these blocks is cheap. This allows us to eliminate
+:math:`\Delta z` by observing that :math:`\Delta z = C^{-1}(w - E^\top
+\Delta y)`, giving us
+
+.. math:: \left[B - EC^{-1}E^\top\right] \Delta y = v - EC^{-1}w\ .
+ :label: schur
+
+The matrix
+
+.. math:: S = B - EC^{-1}E^\top
+
+is the Schur complement of :math:`C` in :math:`H`. It is also known as
+the *reduced camera matrix*, because the only variables
+participating in :eq:`schur` are the ones corresponding to the
+cameras. :math:`S \in \mathbb{R}^{pc\times pc}` is a block structured
+symmetric positive definite matrix, with blocks of size :math:`c\times
+c`. The block :math:`S_{ij}` corresponding to the pair of images
+:math:`i` and :math:`j` is non-zero if and only if the two images
+observe at least one common point.
+
+
+Now, eq-linear2 can be solved by first forming :math:`S`, solving for
+:math:`\Delta y`, and then back-substituting :math:`\Delta y` to
+obtain the value of :math:`\Delta z`. Thus, the solution of what was
+an :math:`n\times n`, :math:`n=pc+qs` linear system is reduced to the
+inversion of the block diagonal matrix :math:`C`, a few matrix-matrix
+and matrix-vector multiplies, and the solution of block sparse
+:math:`pc\times pc` linear system :eq:`schur`. For almost all
+problems, the number of cameras is much smaller than the number of
+points, :math:`p \ll q`, thus solving :eq:`schur` is
+significantly cheaper than solving :eq:`linear2`. This is the
+*Schur complement trick* [Brown]_.
+
+This still leaves open the question of solving :eq:`schur`. The
+method of choice for solving symmetric positive definite systems
+exactly is via the Cholesky factorization [TrefethenBau]_ and
+depending upon the structure of the matrix, there are, in general, two
+options. The first is direct factorization, where we store and factor
+:math:`S` as a dense matrix [TrefethenBau]_. This method has
+:math:`O(p^2)` space complexity and :math:`O(p^3)` time complexity and
+is only practical for problems with up to a few hundred cameras. Ceres
+implements this strategy as the ``DENSE_SCHUR`` solver.
+
+
+But, :math:`S` is typically a fairly sparse matrix, as most images
+only see a small fraction of the scene. This leads us to the second
+option: Sparse Direct Methods. These methods store :math:`S` as a
+sparse matrix, use row and column re-ordering algorithms to maximize
+the sparsity of the Cholesky decomposition, and focus their compute
+effort on the non-zero part of the factorization [Chen]_. Sparse
+direct methods, depending on the exact sparsity structure of the Schur
+complement, allow bundle adjustment algorithms to significantly scale
+up over those based on dense factorization. Ceres implements this
+strategy as the ``SPARSE_SCHUR`` solver.
+
+.. _section-cgnr:
+
+``CGNR``
+^^^^^^^^
+
+For general sparse problems, if the problem is too large for
+``CHOLMOD`` or a sparse linear algebra library is not linked into
+Ceres, another option is the ``CGNR`` solver. This solver uses the
+Conjugate Gradients solver on the *normal equations*, but without
+forming the normal equations explicitly. It exploits the relation
+
+.. math::
+ H x = J^\top J x = J^\top(J x)
+
+
+When the user chooses ``ITERATIVE_SCHUR`` as the linear solver, Ceres
+automatically switches from the exact step algorithm to an inexact
+step algorithm.
+
+.. _section-iterative_schur:
+
+``ITERATIVE_SCHUR``
+^^^^^^^^^^^^^^^^^^^
+
+Another option for bundle adjustment problems is to apply PCG to the
+reduced camera matrix :math:`S` instead of :math:`H`. One reason to do
+this is that :math:`S` is a much smaller matrix than :math:`H`, but
+more importantly, it can be shown that :math:`\kappa(S)\leq
+\kappa(H)`. Cseres implements PCG on :math:`S` as the
+``ITERATIVE_SCHUR`` solver. When the user chooses ``ITERATIVE_SCHUR``
+as the linear solver, Ceres automatically switches from the exact step
+algorithm to an inexact step algorithm.
+
+The cost of forming and storing the Schur complement :math:`S` can be
+prohibitive for large problems. Indeed, for an inexact Newton solver
+that computes :math:`S` and runs PCG on it, almost all of its time is
+spent in constructing :math:`S`; the time spent inside the PCG
+algorithm is negligible in comparison. Because PCG only needs access
+to :math:`S` via its product with a vector, one way to evaluate
+:math:`Sx` is to observe that
+
+.. math:: x_1 &= E^\top x
+.. math:: x_2 &= C^{-1} x_1
+.. math:: x_3 &= Ex_2\\
+.. math:: x_4 &= Bx\\
+.. math:: Sx &= x_4 - x_3
+ :label: schurtrick1
+
+Thus, we can run PCG on :math:`S` with the same computational effort
+per iteration as PCG on :math:`H`, while reaping the benefits of a
+more powerful preconditioner. In fact, we do not even need to compute
+:math:`H`, :eq:`schurtrick1` can be implemented using just the columns
+of :math:`J`.
+
+Equation :eq:`schurtrick1` is closely related to *Domain
+Decomposition methods* for solving large linear systems that arise in
+structural engineering and partial differential equations. In the
+language of Domain Decomposition, each point in a bundle adjustment
+problem is a domain, and the cameras form the interface between these
+domains. The iterative solution of the Schur complement then falls
+within the sub-category of techniques known as Iterative
+Sub-structuring [Saad]_ [Mathew]_.
+
+.. _section-preconditioner:
+
+Preconditioner
+--------------
+
+The convergence rate of Conjugate Gradients for
+solving :eq:`normal` depends on the distribution of eigenvalues
+of :math:`H` [Saad]_. A useful upper bound is
+:math:`\sqrt{\kappa(H)}`, where, :math:`\kappa(H)` is the condition
+number of the matrix :math:`H`. For most bundle adjustment problems,
+:math:`\kappa(H)` is high and a direct application of Conjugate
+Gradients to :eq:`normal` results in extremely poor performance.
+
+The solution to this problem is to replace :eq:`normal` with a
+*preconditioned* system. Given a linear system, :math:`Ax =b` and a
+preconditioner :math:`M` the preconditioned system is given by
+:math:`M^{-1}Ax = M^{-1}b`. The resulting algorithm is known as
+Preconditioned Conjugate Gradients algorithm (PCG) and its worst case
+complexity now depends on the condition number of the *preconditioned*
+matrix :math:`\kappa(M^{-1}A)`.
+
+The computational cost of using a preconditioner :math:`M` is the cost
+of computing :math:`M` and evaluating the product :math:`M^{-1}y` for
+arbitrary vectors :math:`y`. Thus, there are two competing factors to
+consider: How much of :math:`H`'s structure is captured by :math:`M`
+so that the condition number :math:`\kappa(HM^{-1})` is low, and the
+computational cost of constructing and using :math:`M`. The ideal
+preconditioner would be one for which :math:`\kappa(M^{-1}A)
+=1`. :math:`M=A` achieves this, but it is not a practical choice, as
+applying this preconditioner would require solving a linear system
+equivalent to the unpreconditioned problem. It is usually the case
+that the more information :math:`M` has about :math:`H`, the more
+expensive it is use. For example, Incomplete Cholesky factorization
+based preconditioners have much better convergence behavior than the
+Jacobi preconditioner, but are also much more expensive.
+
+
+The simplest of all preconditioners is the diagonal or Jacobi
+preconditioner, i.e., :math:`M=\operatorname{diag}(A)`, which for
+block structured matrices like :math:`H` can be generalized to the
+block Jacobi preconditioner.
+
+For ``ITERATIVE_SCHUR`` there are two obvious choices for block
+diagonal preconditioners for :math:`S`. The block diagonal of the
+matrix :math:`B` [Mandel]_ and the block diagonal :math:`S`, i.e, the
+block Jacobi preconditioner for :math:`S`. Ceres's implements both of
+these preconditioners and refers to them as ``JACOBI`` and
+``SCHUR_JACOBI`` respectively.
+
+For bundle adjustment problems arising in reconstruction from
+community photo collections, more effective preconditioners can be
+constructed by analyzing and exploiting the camera-point visibility
+structure of the scene [KushalAgarwal]. Ceres implements the two
+visibility based preconditioners described by Kushal & Agarwal as
+``CLUSTER_JACOBI`` and ``CLUSTER_TRIDIAGONAL``. These are fairly new
+preconditioners and Ceres' implementation of them is in its early
+stages and is not as mature as the other preconditioners described
+above.
+
+.. _section-ordering:
+
+Ordering
+--------
+
+The order in which variables are eliminated in a linear solver can
+have a significant of impact on the efficiency and accuracy of the
+method. For example 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.
+
+Ceres allows the user to provide varying amounts of hints to the
+solver about the variable elimination ordering to use. This can range
+from no hints, where the solver is free to decide the best ordering
+based on the user's choices like the linear solver being used, to an
+exact order in which the variables should be eliminated, and a variety
+of possibilities in between.
+
+Instances of the :class:`ParameterBlockOrdering` class are used to
+communicate this information to Ceres.
+
+Formally an ordering is an ordered partitioning of the parameter
+blocks. Each parameter block belongs to exactly one group, and each
+group has a unique integer associated with it, that determines its
+order in the set of groups. We call these groups *Elimination Groups*
+
+Given such an ordering, Ceres ensures that the parameter blocks in the
+lowest numbered elimination group are eliminated first, and then the
+parameter blocks in the next lowest numbered elimination group and so
+on. Within each elimination group, Ceres is free to order the
+parameter blocks as it chooses. e.g. Consider the linear system
+
+.. math::
+ x + y &= 3\\
+ 2x + 3y &= 7
+
+There are two ways in which it can be solved. First eliminating
+:math:`x` from the two equations, solving for y and then back
+substituting for :math:`x`, or first eliminating :math:`y`, solving
+for :math:`x` and back substituting for :math:`y`. The user can
+construct three orderings here.
+
+1. :math:`\{0: x\}, \{1: y\}` : Eliminate :math:`x` first.
+2. :math:`\{0: y\}, \{1: x\}` : Eliminate :math:`y` first.
+3. :math:`\{0: x, y\}` : Solver gets to decide the elimination order.
+
+Thus, to have Ceres determine the ordering automatically using
+heuristics, put all the variables in the same elimination group. The
+identity of the group does not matter. This is the same as not
+specifying an ordering at all. To control the ordering for every
+variable, create an elimination group per variable, ordering them in
+the desired order.
+
+If the user is using one of the Schur solvers (``DENSE_SCHUR``,
+``SPARSE_SCHUR``, ``ITERATIVE_SCHUR``) and chooses to specify an
+ordering, it must have one important property. The lowest numbered
+elimination group must form an independent set in the graph
+corresponding to the Hessian, or in other words, no two parameter
+blocks in in the first elimination group should co-occur in the same
+residual block. For the best performance, this elimination group
+should be as large as possible. For standard bundle adjustment
+problems, this corresponds to the first elimination group containing
+all the 3d points, and the second containing the all the cameras
+parameter blocks.
+
+If the user leaves the choice to Ceres, then the solver uses an
+approximate maximum independent set algorithm to identify the first
+elimination group [LiSaad]_.
+
+.. _section-solver-options:
+
+:class:`Solver::Options`
+------------------------
+
+.. class:: Solver::Options
+
+ :class:`Solver::Options` controls the overall behavior of the
+ solver. We list the various settings and their default values below.
+
+.. member:: TrustRegionStrategyType Solver::Options::trust_region_strategy_type
+
+ Default: ``LEVENBERG_MARQUARDT``
+
+ The trust region step computation algorithm used by
+ Ceres. Currently ``LEVENBERG_MARQUARDT`` and ``DOGLEG`` are the two
+ valid choices. See :ref:`section-levenberg-marquardt` and
+ :ref:`section-dogleg` for more details.
+
+.. member:: DoglegType Solver::Options::dogleg_type
+
+ Default: ``TRADITIONAL_DOGLEG``
+
+ Ceres supports two different dogleg strategies.
+ ``TRADITIONAL_DOGLEG`` method by Powell and the ``SUBSPACE_DOGLEG``
+ method described by [Byrd]_. See :ref:`section-dogleg` for more
+ details.
+
+.. member:: bool Solver::Options::use_nonmonotonic_steps
+
+ Default: ``false``
+
+ Relax the requirement that the trust-region algorithm take strictly
+ decreasing steps. See :ref:`section-non-monotonic-steps` for more
+ details.
+
+.. member:: int Solver::Options::max_consecutive_nonmonotonic_steps
+
+ Default: ``5``
+
+ The window size used by the step selection algorithm to accept
+ non-monotonic steps.
+
+.. member:: int Solver::Options::max_num_iterations
+
+ Default: ``50``
+
+ Maximum number of iterations for which the solver should run.
+
+.. member:: double Solver::Options::max_solver_time_in_seconds
+
+ Default: ``1e6``
+ Maximum amount of time for which the solver should run.
+
+.. member:: int Solver::Options::num_threads
+
+ Default: ``1``
+
+ Number of threads used by Ceres to evaluate the Jacobian.
+
+.. member:: double Solver::Options::initial_trust_region_radius
+
+ Default: ``1e4``
+
+ The size of the initial trust region. When the
+ ``LEVENBERG_MARQUARDT`` strategy is used, the reciprocal of this
+ number is the initial regularization parameter.
+
+.. member:: double Solver::Options::max_trust_region_radius
+
+ Default: ``1e16``
+
+ The trust region radius is not allowed to grow beyond this value.
+
+.. member:: double Solver::Options::min_trust_region_radius
+
+ Default: ``1e-32``
+
+ The solver terminates, when the trust region becomes smaller than
+ this value.
+
+.. member:: double Solver::Options::min_relative_decrease
+
+ Default: ``1e-3``
+
+ Lower threshold for relative decrease before a trust-region step is
+ acceped.
+
+.. member:: double Solver::Options::lm_min_diagonal
+
+ Default: ``1e6``
+
+ The ``LEVENBERG_MARQUARDT`` strategy, uses a diagonal matrix to
+ regularize the the trust region step. This is the lower bound on
+ the values of this diagonal matrix.
+
+.. member:: double Solver::Options::lm_max_diagonal
+
+ Default: ``1e32``
+
+ The ``LEVENBERG_MARQUARDT`` strategy, uses a diagonal matrix to
+ regularize the the trust region step. This is the upper bound on
+ the values of this diagonal matrix.
+
+.. member:: int Solver::Options::max_num_consecutive_invalid_steps
+
+ Default: ``5``
+
+ The step returned by a trust region strategy can sometimes be
+ numerically invalid, usually because of conditioning
+ issues. Instead of crashing or stopping the optimization, the
+ optimizer can go ahead and try solving with a smaller trust
+ region/better conditioned problem. This parameter sets the number
+ of consecutive retries before the minimizer gives up.
+
+.. member:: double Solver::Options::function_tolerance
+
+ Default: ``1e-6``
+
+ Solver terminates if
+
+ .. math:: \frac{|\Delta \text{cost}|}{\text{cost} < \text{function_tolerance}}
+
+ where, :math:`\Delta \text{cost}` is the change in objective function
+ value (up or down) in the current iteration of Levenberg-Marquardt.
+
+.. member:: double Solver::Options::gradient_tolerance
+
+ Default: ``1e-10``
+
+ Solver terminates if
+
+ .. math:: \frac{\|g(x)\|_\infty}{\|g(x_0)\|_\infty} < \text{gradient_tolerance}
+
+ where :math:`\|\cdot\|_\infty` refers to the max norm, and :math:`x_0` is
+ the vector of initial parameter values.
+
+.. member:: double Solver::Options::parameter_tolerance
+
+ Default: ``1e-8``
+
+ Solver terminates if
+
+ .. math:: \|\Delta x\| < (\|x\| + \text{parameter_tolerance}) * \text{parameter_tolerance}
+
+ where :math:`\Delta x` is the step computed by the linear solver in the
+ current iteration of Levenberg-Marquardt.
+
+.. member:: LinearSolverType Solver::Options::linear_solver_type
+
+ Default: ``SPARSE_NORMAL_CHOLESKY`` / ``DENSE_QR``
+
+ Type of linear solver used to compute the solution to the linear
+ least squares problem in each iteration of the Levenberg-Marquardt
+ algorithm. If Ceres is build with ``SuiteSparse`` linked in then
+ the default is ``SPARSE_NORMAL_CHOLESKY``, it is ``DENSE_QR``
+ otherwise.
+
+.. member:: PreconditionerType Solver::Options::preconditioner_type
+
+ Default: ``JACOBI``
+
+ The preconditioner used by the iterative linear solver. The default
+ is the block Jacobi preconditioner. Valid values are (in increasing
+ order of complexity) ``IDENTITY``, ``JACOBI``, ``SCHUR_JACOBI``,
+ ``CLUSTER_JACOBI`` and ``CLUSTER_TRIDIAGONAL``. See
+ :ref:`section-preconditioner` for more details.
+
+.. member:: SparseLinearAlgebraLibrary Solver::Options::sparse_linear_algebra_library
+
+ Default:``SUITE_SPARSE``
+
+ Ceres supports the use of two sparse linear algebra libraries,
+ ``SuiteSparse``, which is enabled by setting this parameter to
+ ``SUITE_SPARSE`` and ``CXSparse``, which can be selected by setting
+ this parameter to ```CX_SPARSE``. ``SuiteSparse`` is a
+ sophisticated and complex sparse linear algebra library and should
+ be used in general. If your needs/platforms prevent you from using
+ ``SuiteSparse``, consider using ``CXSparse``, which is a much
+ smaller, easier to build library. As can be expected, its
+ performance on large problems is not comparable to that of
+ ``SuiteSparse``.
+
+.. member:: int Solver::Options::num_linear_solver_threads
+
+ Default: ``1``
+
+ Number of threads used by the linear solver.
+
+.. member:: bool Solver::Options::use_inner_iterations
+
+ Default: ``false``
+
+ Use a non-linear version of a simplified variable projection
+ algorithm. Essentially this amounts to doing a further optimization
+ on each Newton/Trust region step using a coordinate descent
+ algorithm. For more details, see :ref:`section-inner-iterations`.
+
+.. member:: ParameterBlockOrdering* Solver::Options::inner_iteration_ordering
+
+ Default: ``NULL``
+
+ If :member:`Solver::Options::use_inner_iterations` true, then the user has
+ two choices.
+
+ 1. Let the solver heuristically decide which parameter blocks to
+ optimize in each inner iteration. To do this, set
+ :member:`Solver::Options::inner_iteration_ordering` to ``NULL``.
+
+ 2. Specify a collection of of ordered independent sets. The lower
+ numbered groups are optimized before the higher number groups
+ during the inner optimization phase. Each group must be an
+ independent set.
+
+ See :ref:`section-ordering` for more details.
+
+.. member:: ParameterBlockOrdering* Solver::Options::linear_solver_ordering
+
+ Default: ``NULL``
+
+ An instance of the ordering object informs the solver about the
+ desired order in which parameter blocks should be eliminated by the
+ linear solvers. See section~\ref{sec:ordering`` for more details.
+
+ If ``NULL``, the solver is free to choose an ordering that it
+ thinks is best. Note: currently, this option only has an effect on
+ the Schur type solvers, support for the ``SPARSE_NORMAL_CHOLESKY``
+ solver is forth coming.
+
+ See :ref:`section-ordering` for more details.
+
+.. member:: bool Solver::Options::use_block_amd
+
+ Default: ``true``
+
+ By virtue of the modeling layer in Ceres being block oriented, all
+ the matrices used by Ceres are also block oriented. When doing
+ sparse direct factorization of these matrices, the fill-reducing
+ ordering algorithms can either be run on the block or the scalar
+ form of these matrices. Running it on the block form exposes more
+ of the super-nodal structure of the matrix to the Cholesky
+ factorization routines. This leads to substantial gains in
+ factorization performance. Setting this parameter to true, enables
+ the use of a block oriented Approximate Minimum Degree ordering
+ algorithm. Settings it to ``false``, uses a scalar AMD
+ algorithm. This option only makes sense when using
+ :member:`Solver::Options::sparse_linear_algebra_library` = ``SUITE_SPARSE``
+ as it uses the ``AMD`` package that is part of ``SuiteSparse``.
+
+.. member:: int Solver::Options::linear_solver_min_num_iterations
+
+ Default: ``1``
+
+ Minimum number of iterations used by the linear solver. This only
+ makes sense when the linear solver is an iterative solver, e.g.,
+ ``ITERATIVE_SCHUR`` or ``CGNR``.
+
+.. member:: int Solver::Options::linear_solver_max_num_iterations
+
+ Default: ``500``
+
+ Minimum number of iterations used by the linear solver. This only
+ makes sense when the linear solver is an iterative solver, e.g.,
+ ``ITERATIVE_SCHUR`` or ``CGNR``.
+
+.. member:: double Solver::Options::eta
+
+ Default: ``1e-1``
+
+ Forcing sequence parameter. The truncated Newton solver uses this
+ number to control the relative accuracy with which the Newton step
+ is computed. This constant is passed to
+ ``ConjugateGradientsSolver`` which uses it to terminate the
+ iterations when
+
+ .. math:: \frac{Q_i - Q_{i-1}}{Q_i} < \frac{\eta}{i}
+
+.. member:: bool Solver::Options::jacobi_scaling
+
+ Default: ``true``
+
+ ``true`` means that the Jacobian is scaled by the norm of its
+ columns before being passed to the linear solver. This improves the
+ numerical conditioning of the normal equations.
+
+.. member:: LoggingType Solver::Options::logging_type
+
+ Default: ``PER_MINIMIZER_ITERATION``
+
+.. member:: bool Solver::Options::minimizer_progress_to_stdout
+
+ Default: ``false``
+
+ By default the :class:`Minimizer` progress is logged to ``STDERR``
+ depending on the ``vlog`` level. If this flag is set to true, and
+ :member:`Solver::Options::logging_type` is not ``SILENT``, the logging
+ output is sent to ``STDOUT``.
+
+.. member:: bool Solver::Options::return_initial_residuals
+
+ Default: ``false``
+
+.. member:: bool Solver::Options::return_final_residuals
+
+ Default: ``false``
+
+ If true, the vectors :member:`Solver::Summary::initial_residuals` and
+ :member:`Solver::Summary::final_residuals` are filled with the residuals
+ before and after the optimization. The entries of these vectors are
+ in the order in which ResidualBlocks were added to the Problem
+ object.
+
+.. member:: bool Solver::Options::return_initial_gradient
+
+ Default: ``false``
+
+.. member:: bool Solver::Options::return_final_gradient
+
+ Default: ``false``
+
+ If true, the vectors :member:`Solver::Summary::initial_gradient` and
+ :member:`Solver::Summary::final_gradient` are filled with the gradient
+ before and after the optimization. The entries of these vectors are
+ in the order in which ParameterBlocks were added to the Problem
+ object.
+
+ Since :member:`Problem::AddResidualBlock` adds ParameterBlocks to
+ the :class:`Problem` automatically if they do not already exist,
+ if you wish to have explicit control over the ordering of the
+ vectors, then use :member:`Problem::AddParameterBlock` to
+ explicitly add the ParameterBlocks in the order desired.
+
+.. member:: bool Solver::Options::return_initial_jacobian
+
+ Default: ``false``
+
+.. member:: bool Solver::Options::return_initial_jacobian
+
+ Default: ``false``
+
+ If ``true``, the Jacobian matrices before and after the
+ optimization are returned in
+ :member:`Solver::Summary::initial_jacobian` and
+ :member:`Solver::Summary::final_jacobian` respectively.
+
+ The rows of these matrices are in the same order in which the
+ ResidualBlocks were added to the Problem object. The columns are in
+ the same order in which the ParameterBlocks were added to the
+ Problem object.
+
+ Since :member:`Problem::AddResidualBlock` adds ParameterBlocks to
+ the :class:`Problem` automatically if they do not already exist,
+ if you wish to have explicit control over the ordering of the
+ vectors, then use :member:`Problem::AddParameterBlock` to
+ explicitly add the ParameterBlocks in the order desired.
+
+ The Jacobian matrices are stored as compressed row sparse
+ matrices. Please see ``include/ceres/crs_matrix.h`` for more
+ details of the format.
+
+.. member:: vector<int> Solver::Options::lsqp_iterations_to_dump
+
+ Default: ``empty``
+
+ List of iterations at which the optimizer should dump the linear
+ least squares problem to disk. Useful for testing and
+ benchmarking. If ``empty``, no problems are dumped.
+
+.. member:: string Solver::Options::lsqp_dump_directory
+
+ Default: ``/tmp``
+
+ If :member:`Solver::Options::lsqp_iterations_to_dump` is non-empty, then
+ this setting determines the directory to which the files containing
+ the linear least squares problems are written to.
+
+.. member:: DumpFormatType Solver::Options::lsqp_dump_format
+
+ Default: ``TEXTFILE``
+
+ The format in which linear least squares problems should be logged
+ when :member:`Solver::Options::lsqp_iterations_to_dump` 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.
+
+ * ``PROTOBUF`` Write out the linear least squares problem to the
+ directory pointed to by :member:`Solver::Options::lsqp_dump_directory` as
+ a protocol buffer. ``linear_least_squares_problems.h/cc``
+ contains routines for loading these problems. For details on the
+ on disk format used, see ``matrix.proto``. The files are named
+ ``lm_iteration_???.lsqp``. This requires that ``protobuf`` be
+ linked into Ceres Solver.
+
+ * ``TEXTFILE`` Write out the linear least squares problem to the
+ directory pointed to by member::`Solver::Options::lsqp_dump_directory` as
+ text files which can be read into ``MATLAB/Octave``. The Jacobian
+ is dumped as a text file containing :math:`(i,j,s)` triplets, the
+ vectors :math:`D`, `x` and `f` are dumped as text files
+ containing a list of their values.
+
+ A ``MATLAB/Octave`` script called ``lm_iteration_???.m`` is also
+ output, which can be used to parse and load the problem into memory.
+
+.. member:: bool Solver::Options::check_gradients
+
+ Default: ``false``
+
+ Check all Jacobians computed by each residual block with finite
+ differences. This is expensive since it involves computing the
+ derivative by normal means (e.g. user specified, autodiff, etc),
+ then also computing it using finite differences. The results are
+ compared, and if they differ substantially, details are printed to
+ the log.
+
+.. member:: double Solver::Options::gradient_check_relative_precision
+
+ Default: ``1e08``
+
+ Precision to check for in the gradient checker. If the relative
+ difference between an element in a Jacobian exceeds this number,
+ then the Jacobian for that cost term is dumped.
+
+.. member:: double Solver::Options::numeric_derivative_relative_step_size
+
+ Default: ``1e-6``
+
+ Relative shift used for taking numeric derivatives. For finite
+ differencing, each dimension is evaluated at slightly shifted
+ values, e.g., for forward differences, the numerical derivative is
+
+ .. math::
+
+ \delta &= numeric\_derivative\_relative\_step\_size\\
+ \Delta f &= \frac{f((1 + \delta) x) - f(x)}{\delta x}
+
+ The finite differencing is done along each dimension. The reason to
+ use a relative (rather than absolute) step size is that this way,
+ numeric differentiation works for functions where the arguments are
+ typically large (e.g. :math:`10^9`) and when the values are small
+ (e.g. :math:`10^{-5}`). It is possible to construct *torture cases*
+ which break this finite difference heuristic, but they do not come
+ up often in practice.
+
+.. member:: vector<IterationCallback> Solver::Options::callbacks
+
+ Callbacks that are executed at the end of each iteration of the
+ :class:`Minimizer`. They are executed in the order that they are
+ specified in this vector. By default, parameter blocks are updated
+ only at the end of the optimization, i.e when the
+ :class:`Minimizer` terminates. This behavior is controlled by
+ :member:`Solver::Options::update_state_every_variable`. If the user wishes
+ to have access to the update parameter blocks when his/her
+ callbacks are executed, then set
+ :member:`Solver::Options::update_state_every_iteration` to true.
+
+ The solver does NOT take ownership of these pointers.
+
+.. member:: bool Solver::Options::update_state_every_iteration
+
+ Default: ``false``
+
+ Normally the parameter blocks are only updated when the solver
+ terminates. Setting this to true update them in every
+ iteration. This setting is useful when building an interactive
+ application using Ceres and using an :class:`IterationCallback`.
+
+.. member:: string Solver::Options::solver_log
+
+ Default: ``empty``
+
+ If non-empty, a summary of the execution of the solver is recorded
+ to this file. This file is used for recording and Ceres'
+ performance. Currently, only the iteration number, total time and
+ the objective function value are logged. The format of this file is
+ expected to change over time as the performance evaluation
+ framework is fleshed out.
+
+:class:`ParameterBlockOrdering`
+-------------------------------
+
+.. class:: ParameterBlockOrdering
+
+ TBD
+
+:class:`IterationCallback`
+--------------------------
+
+.. class:: IterationCallback
+
+ TBD
+
+:class:`CRSMatrix`
+------------------
+
+.. class:: CRSMatrix
+
+ TBD
+
+:class:`Solver::Summary`
+------------------------
+
+.. class:: Solver::Summary
+
+ TBD
+
+:class:`GradientChecker`
+------------------------
+
+.. class:: GradientChecker
+
+
+
+
+
diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst
new file mode 100644
index 0000000..df910ef
--- /dev/null
+++ b/docs/source/tutorial.rst
@@ -0,0 +1,554 @@
+.. _chapter-tutorial:
+
+========
+Tutorial
+========
+
+.. highlight:: c++
+
+.. _section-hello-world:
+
+Hello World!
+============
+
+To get started, let us consider the problem of finding the minimum of
+the function
+
+.. math:: \frac{1}{2}(10 -x)^2.
+
+This is a trivial problem, whose minimum is located at :math:`x = 10`,
+but it is a good place to start to illustrate the basics of solving a
+problem with Ceres [#f1]_.
+
+Let us write this problem as a non-linear least squares problem by
+defining the scalar residual function :math:`f_1(x) = 10 - x`. Then
+:math:`F(x) = [f_1(x)]` is a residual vector with exactly one
+component.
+
+When solving a problem with Ceres, the first thing to do is to define
+a subclass of CostFunction. It is responsible for computing
+the value of the residual function and its derivative (also known as
+the Jacobian) with respect to :math:`x`.
+
+.. code-block:: c++
+
+ class SimpleCostFunction : public ceres::SizedCostFunction<1, 1> {
+ public:
+ virtual ~SimpleCostFunction() {}
+ virtual bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const {
+ const double x = parameters[0][0];
+ residuals[0] = 10 - x;
+
+ // Compute the Jacobian if asked for.
+ if (jacobians != NULL && jacobians[0] != NULL) {
+ jacobians[0][0] = -1;
+ }
+ return true;
+ }
+ };
+
+
+SimpleCostFunction is provided with an input array of
+parameters, an output array for residuals and an optional output array
+for Jacobians. In our example, there is just one parameter and one
+residual and this is known at compile time, therefore we can save some
+code and instead of inheriting from CostFunction, we can
+instaed inherit from the templated SizedCostFunction class.
+
+
+The jacobians array is optional, Evaluate is expected to check when it
+is non-null, and if it is the case then fill it with the values of the
+derivative of the residual function. In this case since the residual
+function is linear, the Jacobian is constant.
+
+Once we have a way of computing the residual vector, it is now time to
+construct a Non-linear least squares problem using it and have Ceres
+solve it.
+
+.. code-block:: c++
+
+ int main(int argc, char** argv) {
+ double x = 5.0;
+ ceres::Problem problem;
+
+ // The problem object takes ownership of the newly allocated
+ // SimpleCostFunction and uses it to optimize the value of x.
+ problem.AddResidualBlock(new SimpleCostFunction, NULL, &x);
+
+ // Run the solver!
+ Solver::Options options;
+ options.max_num_iterations = 10;
+ options.linear_solver_type = ceres::DENSE_QR;
+ options.minimizer_progress_to_stdout = true;
+ Solver::Summary summary;
+ Solve(options, &problem, &summary);
+ std::cout << summary.BriefReport() << "\n";
+ std::cout << "x : 5.0 -> " << x << "\n";
+ return 0;
+ }
+
+
+Compiling and running the program gives us
+
+.. code-block:: bash
+
+ 0: f: 1.250000e+01 d: 0.00e+00 g: 5.00e+00 h: 0.00e+00 rho: 0.00e+00 mu: 1.00e+04 li: 0 it: 0.00e+00 tt: 0.00e+00
+ 1: f: 1.249750e-07 d: 1.25e+01 g: 5.00e-04 h: 5.00e+00 rho: 1.00e+00 mu: 3.00e+04 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 2: f: 1.388518e-16 d: 1.25e-07 g: 1.67e-08 h: 5.00e-04 rho: 1.00e+00 mu: 9.00e+04 li: 1 it: 0.00e+00 tt: 0.00e+00
+ Ceres Solver Report: Iterations: 2, Initial cost: 1.250000e+01, Final cost: 1.388518e-16, Termination: PARAMETER_TOLERANCE.
+ x : 5.0 -> 10
+
+
+Starting from a :math:`x=5`, the solver in two iterations goes to 10
+[#f2]_. The careful reader will note that this is a linear problem and
+one linear solve should be enough to get the optimal value. The
+default configuration of the solver is aimed at non-linear problems,
+and for reasons of simplicity we did not change it in this example. It
+is indeed possible to obtain the solution to this problem using Ceres
+in one iteration. Also note that the solver did get very close to the
+optimal function value of 0 in the very first iteration. We will
+discuss these issues in greater detail when we talk about convergence
+and parameter settings for Ceres.
+
+.. rubric:: Footnotes
+
+.. [#f1] Full working code for this and other
+ examples in this manual can be found in the examples directory. Code
+ for this example can be found in ``examples/quadratic.cc``.
+
+.. [#f2] Actually the solver ran for three iterations, and it was
+ by looking at the value returned by the linear solver in the third
+ iteration, it observed that the update to the parameter block was too
+ small and declared convergence. Ceres only prints out the display at
+ the end of an iteration, and terminates as soon as it detects
+ convergence, which is why you only see two iterations here and not
+ three.
+
+
+.. _section-powell:
+
+Powell's Function
+=================
+
+Consider now a slightly more complicated example -- the minimization
+of Powell's function. Let :math:`x = \left[x_1, x_2, x_3, x_4 \right]`
+and
+
+.. math::
+
+ \begin{align}
+ f_1(x) &= x_1 + 10x_2 \\
+ f_2(x) &= \sqrt{5} (x_3 - x_4)\\
+ f_3(x) &= (x_2 - 2x_3)^2\\
+ f_4(x) &= \sqrt{10} (x_1 - x_4)^2\\
+ F(x) & = \left[f_1(x),\ f_2(x),\ f_3(x),\ f_4(x) \right]
+ \end{align}
+
+
+:math:`F(x)` is a function of four parameters, and has four
+residuals. Now, one way to solve this problem would be to define four
+CostFunction objects that compute the residual and Jacobians. e.g. the
+following code shows the implementation for :math:`f_4(x)`.
+
+.. code-block:: c++
+
+ class F4 : public ceres::SizedCostFunction<1, 4> {
+ public:
+ virtual ~F4() {}
+ virtual bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** jacobians) const {
+ double x1 = parameters[0][0];
+ double x4 = parameters[0][3];
+
+ residuals[0] = sqrt(10.0) * (x1 - x4) * (x1 - x4)
+
+ if (jacobians != NULL) {
+ jacobians[0][0] = 2.0 * sqrt(10.0) * (x1 - x4);
+ jacobians[0][1] = 0.0;
+ jacobians[0][2] = 0.0;
+ jacobians[0][3] = -2.0 * sqrt(10.0) * (x1 - x4);
+ }
+ return true;
+ }
+ };
+
+
+But this can get painful very quickly, especially for residuals
+involving complicated multi-variate terms. Ceres provides two ways
+around this problem. Numeric and automatic symbolic differentiation.
+
+Automatic Differentiation
+-------------------------
+
+With its automatic differentiation support, Ceres allows you to define
+templated objects/functors that will compute the ``residual`` and it
+takes care of computing the Jacobians as needed and filling the
+``jacobians`` arrays with them. For example, for :math:`f_4(x)` we
+define
+
+.. code-block:: c++
+
+ class F4 {
+ public:
+ template <typename T> bool operator()(const T* const x1,
+ const T* const x4,
+ T* residual) const {
+ residual[0] = T(sqrt(10.0)) * (x1[0] - x4[0]) * (x1[0] - x4[0]);
+ return true;
+ }
+ };
+
+
+The important thing to note here is that ``operator()`` is a templated
+method, which assumes that all its inputs and outputs are of some type
+``T``. The reason for using templates here is because Ceres will call
+``F4::operator<T>()``, with ``T=double`` when just the residual is
+needed, and with a special type ``T=Jet`` when the Jacobians are
+needed.
+
+Note also that the parameters are not packed
+into a single array, they are instead passed as separate arguments to
+``operator()``. Similarly we can define classes ``F1``,``F2``
+and ``F4``. Then let us consider the construction and solution
+of the problem. For brevity we only describe the relevant bits of
+code [#f3]_ .
+
+
+.. code-block:: c++
+
+ double x1 = 3.0; double x2 = -1.0; double x3 = 0.0; double x4 = 1.0;
+ // Add residual terms to the problem using the using the autodiff
+ // wrapper to get the derivatives automatically.
+ problem.AddResidualBlock(
+ new ceres::AutoDiffCostFunction<F1, 1, 1, 1>(new F1), NULL, &x1, &x2);
+ problem.AddResidualBlock(
+ new ceres::AutoDiffCostFunction<F2, 1, 1, 1>(new F2), NULL, &x3, &x4);
+ problem.AddResidualBlock(
+ new ceres::AutoDiffCostFunction<F3, 1, 1, 1>(new F3), NULL, &x2, &x3)
+ problem.AddResidualBlock(
+ new ceres::AutoDiffCostFunction<F4, 1, 1, 1>(new F4), NULL, &x1, &x4);
+
+
+A few things are worth noting in the code above. First, the object
+being added to the ``Problem`` is an ``AutoDiffCostFunction`` with
+``F1``, ``F2``, ``F3`` and ``F4`` as template parameters. Second, each
+``ResidualBlock`` only depends on the two parameters that the
+corresponding residual object depends on and not on all four
+parameters.
+
+Compiling and running ``powell.cc`` gives us:
+
+.. code-block:: bash
+
+ Initial x1 = 3, x2 = -1, x3 = 0, x4 = 1
+ 0: f: 1.075000e+02 d: 0.00e+00 g: 1.55e+02 h: 0.00e+00 rho: 0.00e+00 mu: 1.00e+04 li: 0 it: 0.00e+00 tt: 0.00e+00
+ 1: f: 5.036190e+00 d: 1.02e+02 g: 2.00e+01 h: 2.16e+00 rho: 9.53e-01 mu: 3.00e+04 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 2: f: 3.148168e-01 d: 4.72e+00 g: 2.50e+00 h: 6.23e-01 rho: 9.37e-01 mu: 9.00e+04 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 3: f: 1.967760e-02 d: 2.95e-01 g: 3.13e-01 h: 3.08e-01 rho: 9.37e-01 mu: 2.70e+05 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 4: f: 1.229900e-03 d: 1.84e-02 g: 3.91e-02 h: 1.54e-01 rho: 9.37e-01 mu: 8.10e+05 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 5: f: 7.687123e-05 d: 1.15e-03 g: 4.89e-03 h: 7.69e-02 rho: 9.37e-01 mu: 2.43e+06 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 6: f: 4.804625e-06 d: 7.21e-05 g: 6.11e-04 h: 3.85e-02 rho: 9.37e-01 mu: 7.29e+06 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 7: f: 3.003028e-07 d: 4.50e-06 g: 7.64e-05 h: 1.92e-02 rho: 9.37e-01 mu: 2.19e+07 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 8: f: 1.877006e-08 d: 2.82e-07 g: 9.54e-06 h: 9.62e-03 rho: 9.37e-01 mu: 6.56e+07 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 9: f: 1.173223e-09 d: 1.76e-08 g: 1.19e-06 h: 4.81e-03 rho: 9.37e-01 mu: 1.97e+08 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 10: f: 7.333425e-11 d: 1.10e-09 g: 1.49e-07 h: 2.40e-03 rho: 9.37e-01 mu: 5.90e+08 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 11: f: 4.584044e-12 d: 6.88e-11 g: 1.86e-08 h: 1.20e-03 rho: 9.37e-01 mu: 1.77e+09 li: 1 it: 0.00e+00 tt: 0.00e+00
+ Ceres Solver Report: Iterations: 12, Initial cost: 1.075000e+02, Final cost: 4.584044e-12, Termination: GRADIENT_TOLERANCE.
+ Final x1 = 0.00116741, x2 = -0.000116741, x3 = 0.000190535, x4 = 0.000190535
+
+It is easy to see that the optimal solution to this problem is at
+:math:`x_1=0, x_2=0, x_3=0, x_4=0` with an objective function value of
+:math:`0`. In 10 iterations, Ceres finds a solution with an objective
+function value of :math:`4\times 10^{-12}`.
+
+Numeric Differentiation
+-----------------------
+
+In some cases, its not possible to define a templated cost functor. In
+such a situation, numerical differentiation can be used. The user
+defines a functor which computes the residual value and construct a
+``NumericDiffCostFunction`` using it. e.g., for ``F4``, the
+corresponding functor would be
+
+.. code-block:: c++
+
+ class F4 {
+ public:
+ bool operator()(const double* const x1,
+ const double* const x4,
+ double* residual) const {
+ residual[0] = sqrt(10.0) * (x1[0] - x4[0]) * (x1[0] - x4[0]);
+ return true;
+ }
+ };
+
+
+Which can then be wrapped ``NumericDiffCostFunction`` and added to the
+``Problem`` as follows
+
+.. code-block:: c++
+
+ problem.AddResidualBlock(
+ new ceres::NumericDiffCostFunction<F4, ceres::CENTRAL, 1, 1, 1>(new F4), NULL, &x1, &x4);
+
+
+The construction looks almost identical to the used for automatic
+differentiation, except for an extra template parameter that indicates
+the kind of finite differencing scheme to be used for computing the
+numerical derivatives. ``examples/quadratic_numeric_diff.cc`` shows a
+numerically differentiated implementation of
+``examples/quadratic.cc``.
+
+**We recommend that if possible, automatic differentiation should be
+used. The use of C++ templates makes automatic differentiation
+extremely efficient, whereas numeric differentiation can be quite
+expensive, prone to numeric errors and leads to slower convergence.**
+
+
+.. rubric:: Footnotes
+
+.. [#f3] The full source code for this example can be found in ``examples/powell.cc``.
+
+.. _section-fitting:
+
+Fitting a Curve to Data
+=======================
+
+
+The examples we have seen until now are simple optimization problems
+with no data. The original purpose of least squares and non-linear
+least squares analysis was fitting curves to data. It is only
+appropriate that we now consider an example of such a problem
+[#f4]_. It contains data generated by sampling the curve :math:`y =
+e^{0.3x + 0.1}` and adding Gaussian noise with standard deviation
+:math:`\sigma = 0.2`.}. Let us fit some data to the curve
+
+.. math:: y = e^{mx + c}.
+
+We begin by defining a templated object to evaluate the
+residual. There will be a residual for each observation.
+
+.. code-block:: c++
+
+ class ExponentialResidual {
+ public:
+ ExponentialResidual(double x, double y)
+ : x_(x), y_(y) {}
+
+ template <typename T> bool operator()(const T* const m,
+ const T* const c,
+ T* residual) const {
+ residual[0] = T(y_) - exp(m[0] * T(x_) + c[0]);
+ return true;
+ }
+
+ private:
+ // Observations for a sample.
+ const double x_;
+ const double y_;
+ };
+
+Assuming the observations are in a :math:`2n` sized array called ``data``
+the problem construction is a simple matter of creating a
+``CostFunction`` for every observation.
+
+
+.. code-block: c++
+
+ double m = 0.0;
+ double c = 0.0;
+
+ Problem problem;
+ for (int i = 0; i < kNumObservations; ++i) {
+ problem.AddResidualBlock(
+ new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
+ new ExponentialResidual(data[2 * i], data[2 * i + 1])),
+ NULL,
+ &m, &c);
+ }
+
+Compiling and running ``data_fitting.cc`` gives us:
+
+.. code-block:: bash
+
+ 0: f: 1.211734e+02 d: 0.00e+00 g: 3.61e+02 h: 0.00e+00 rho: 0.00e+00 mu: 1.00e+04 li: 0 it: 0.00e+00 tt: 0.00e+00
+ 1: f: 1.211734e+02 d:-2.21e+03 g: 3.61e+02 h: 7.52e-01 rho:-1.87e+01 mu: 5.00e+03 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 2: f: 1.211734e+02 d:-2.21e+03 g: 3.61e+02 h: 7.51e-01 rho:-1.86e+01 mu: 1.25e+03 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 3: f: 1.211734e+02 d:-2.19e+03 g: 3.61e+02 h: 7.48e-01 rho:-1.85e+01 mu: 1.56e+02 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 4: f: 1.211734e+02 d:-2.02e+03 g: 3.61e+02 h: 7.22e-01 rho:-1.70e+01 mu: 9.77e+00 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 5: f: 1.211734e+02 d:-7.34e+02 g: 3.61e+02 h: 5.78e-01 rho:-6.32e+00 mu: 3.05e-01 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 6: f: 3.306595e+01 d: 8.81e+01 g: 4.10e+02 h: 3.18e-01 rho: 1.37e+00 mu: 9.16e-01 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 7: f: 6.426770e+00 d: 2.66e+01 g: 1.81e+02 h: 1.29e-01 rho: 1.10e+00 mu: 2.75e+00 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 8: f: 3.344546e+00 d: 3.08e+00 g: 5.51e+01 h: 3.05e-02 rho: 1.03e+00 mu: 8.24e+00 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 9: f: 1.987485e+00 d: 1.36e+00 g: 2.33e+01 h: 8.87e-02 rho: 9.94e-01 mu: 2.47e+01 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 10: f: 1.211585e+00 d: 7.76e-01 g: 8.22e+00 h: 1.05e-01 rho: 9.89e-01 mu: 7.42e+01 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 11: f: 1.063265e+00 d: 1.48e-01 g: 1.44e+00 h: 6.06e-02 rho: 9.97e-01 mu: 2.22e+02 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 12: f: 1.056795e+00 d: 6.47e-03 g: 1.18e-01 h: 1.47e-02 rho: 1.00e+00 mu: 6.67e+02 li: 1 it: 0.00e+00 tt: 0.00e+00
+ 13: f: 1.056751e+00 d: 4.39e-05 g: 3.79e-03 h: 1.28e-03 rho: 1.00e+00 mu: 2.00e+03 li: 1 it: 0.00e+00 tt: 0.00e+00
+ Ceres Solver Report: Iterations: 13, Initial cost: 1.211734e+02, Final cost: 1.056751e+00, Termination: FUNCTION_TOLERANCE.
+ Initial m: 0 c: 0
+ Final m: 0.291861 c: 0.131439
+
+
+Starting from parameter values :math:`m = 0, c=0` with an initial
+objective function value of :math:`121.173` Ceres finds a solution
+:math:`m= 0.291861, c = 0.131439` with an objective function value of
+:math:`1.05675`. These values are a a bit different than the
+parameters of the original model :math:`m=0.3, c= 0.1`, but this is
+expected. When reconstructing a curve from noisy data, we expect to
+see such deviations. Indeed, if you were to evaluate the objective
+function for :math:`m=0.3, c=0.1`, the fit is worse with an objective
+function value of :math:`1.082425`. The figure below illustrates the fit.
+
+.. figure:: fit.png
+ :figwidth: 500px
+ :height: 400px
+ :align: center
+
+ Least squares data fitting to the curve :math:`y = e^{0.3x +
+ 0.1}`. Observations were generated by sampling this curve uniformly
+ in the interval :math:`x=(0,5)` and adding Gaussian noise with
+ :math:`\sigma = 0.2`.
+
+.. rubric:: Footnotes
+
+.. [#f4] The full source code for this example can be found in ``examples/data_fitting.cc``.
+
+
+Bundle Adjustment
+=================
+
+One of the main reasons for writing Ceres was our need to solve large
+scale bundle adjustment
+problems [HartleyZisserman]_, [Triggs]_.
+
+Given a set of measured image feature locations and correspondences,
+the goal of bundle adjustment is to find 3D point positions and camera
+parameters that minimize the reprojection error. This optimization
+problem is usually formulated as a non-linear least squares problem,
+where the error is the squared :math:`L_2` norm of the difference between
+the observed feature location and the projection of the corresponding
+3D point on the image plane of the camera. Ceres has extensive support
+for solving bundle adjustment problems.
+
+Let us consider the solution of a problem from the `BAL <http://grail.cs.washington.edu/projects/bal/>`_ dataset [#f5]_.
+
+The first step as usual is to define a templated functor that computes
+the reprojection error/residual. The structure of the functor is
+similar to the ``ExponentialResidual``, in that there is an
+instance of this object responsible for each image observation.
+
+
+Each residual in a BAL problem depends on a three dimensional point
+and a nine parameter camera. The nine parameters defining the camera
+can are: Three for rotation as a Rodriquez axis-angle vector, three
+for translation, one for focal length and two for radial distortion.
+The details of this camera model can be found on Noah Snavely's
+`Bundler homepage <http://phototour.cs.washington.edu/bundler/>`_
+and the `BAL homepage <http://grail.cs.washington.edu/projects/bal/>`_.
+
+.. code-block:: c++
+
+ struct SnavelyReprojectionError {
+ SnavelyReprojectionError(double observed_x, double observed_y)
+ : observed_x(observed_x), observed_y(observed_y) {}
+ template <typename T>
+ bool operator()(const T* const camera,
+ const T* const point,
+ T* residuals) const {
+ // camera[0,1,2] are the angle-axis rotation.
+ T p[3];
+ ceres::AngleAxisRotatePoint(camera, point, p);
+ // camera[3,4,5] are the translation.
+ p[0] += camera[3]; p[1] += camera[4]; p[2] += camera[5];
+
+ // Compute the center of distortion. The sign change comes from
+ // the camera model that Noah Snavely's Bundler assumes, whereby
+ // the camera coordinate system has a negative z axis.
+ T xp = - p[0] / p[2];
+ T yp = - p[1] / p[2];
+
+ // Apply second and fourth order radial distortion.
+ const T& l1 = camera[7];
+ const T& l2 = camera[8];
+ T r2 = xp*xp + yp*yp;
+ T distortion = T(1.0) + r2 * (l1 + l2 * r2);
+
+ // Compute final projected point position.
+ const T& focal = camera[6];
+ T predicted_x = focal * distortion * xp;
+ T predicted_y = focal * distortion * yp;
+
+ // The error is the difference between the predicted and observed position.
+ residuals[0] = predicted_x - T(observed_x);
+ residuals[1] = predicted_y - T(observed_y);
+ return true;
+ }
+ double observed_x;
+ double observed_y;
+ } ;
+
+
+Note that unlike the examples before this is a non-trivial function
+and computing its analytic Jacobian is a bit of a pain. Automatic
+differentiation makes our life very simple here. The function
+``AngleAxisRotatePoint`` and other functions for manipulating
+rotations can be found in ``include/ceres/rotation.h``.
+
+Given this functor, the bundle adjustment problem can be constructed
+as follows:
+
+.. code-block:: c++
+
+ // Create residuals for each observation in the bundle adjustment problem. The
+ // parameters for cameras and points are added automatically.
+ ceres::Problem problem;
+ for (int i = 0; i < bal_problem.num_observations(); ++i) {
+ // Each Residual block takes a point and a camera as input and outputs a 2
+ // dimensional residual. Internally, the cost function stores the observed
+ // image location and compares the reprojection against the observation.
+ ceres::CostFunction* cost_function =
+ new ceres::AutoDiffCostFunction<SnavelyReprojectionError, 2, 9, 3>(
+ new SnavelyReprojectionError(
+ bal_problem.observations()[2 * i + 0],
+ bal_problem.observations()[2 * i + 1]));
+ problem.AddResidualBlock(cost_function,
+ NULL /* squared loss */,
+ bal_problem.mutable_camera_for_observation(i),
+ bal_problem.mutable_point_for_observation(i));
+ }
+
+
+Again note that that the problem construction for bundle adjustment is
+very similar to the curve fitting example.
+
+One way to solve this problem is to set
+``Solver::Options::linear_solver_type`` to
+``SPARSE_NORMAL_CHOLESKY`` and call ``Solve``. And while
+this is a reasonable thing to do, bundle adjustment problems have a
+special sparsity structure that can be exploited to solve them much
+more efficiently. Ceres provides three specialized solvers
+(collectively known as Schur based solvers) for this task. The example
+code uses the simplest of them ``DENSE_SCHUR``.
+
+.. code-block:: c++
+
+ ceres::Solver::Options options;
+ options.linear_solver_type = ceres::DENSE_SCHUR;
+ options.minimizer_progress_to_stdout = true;
+ ceres::Solver::Summary summary;
+ ceres::Solve(options, &problem, &summary);
+ std::cout << summary.FullReport() << "\n";
+
+
+For a more sophisticated bundle adjustment example which demonstrates
+the use of Ceres' more advanced features including its various linear
+solvers, robust loss functions and local parameterizations see
+``examples/bundle_adjuster.cc``.
+
+.. rubric:: Footnotes
+
+.. [#f5] The full source code for this example can be found in ``examples/simple_bundle_adjuster.cc``.
diff --git a/docs/source/version_history.rst b/docs/source/version_history.rst
new file mode 100644
index 0000000..6cfb80c
--- /dev/null
+++ b/docs/source/version_history.rst
@@ -0,0 +1,446 @@
+.. _chapter-version-history:
+
+===============
+Version History
+===============
+
+1.5.0
+=====
+
+New Features
+------------
+
+#. Ceres now supports Line search based optimization algorithms in
+ addition to trust region algorithms. Currently there is support for
+ gradient descent, non-linear conjugate gradient and LBFGS search
+ directions.
+
+#. Speedup the robust loss function correction logic when residual is
+ one dimensional.
+
+#. Changed ``NumericDiffCostFunction`` to take functors like
+ ``AutoDiffCostFunction``.
+
+#. Added support for mixing automatic, analytic and numeric
+ differentiation. This is done by adding ``CostFunctionToFunctor``
+ and ``NumericDiffFunctor`` objects.
+
+
+Bug Fixes
+---------
+
+#. Fixed varidic evaluation bug in ``AutoDiff``.
+
+#. Fixed ``SolverImpl`` tests.
+
+#. Fixed a bug in ``DenseSparseMatrix::ToDenseMatrix()``.
+
+#. Fixed an initialization bug in ``ProgramEvaluator``.
+
+
+1.4.0
+=====
+
+
+API Changes
+-----------
+
+The new ordering API breaks existing code. Here the common case fixes.
+
+**Before**
+
+.. code-block:: c++
+
+ options.linear_solver_type = ceres::DENSE_SCHUR
+ options.ordering_type = ceres::SCHUR
+
+**After**
+
+
+.. code-block:: c++
+
+ options.linear_solver_type = ceres::DENSE_SCHUR
+
+
+**Before**
+
+.. code-block:: c++
+
+ options.linear_solver_type = ceres::DENSE_SCHUR;
+ options.ordering_type = ceres::USER;
+ for (int i = 0; i < num_points; ++i) {
+ options.ordering.push_back(my_points[i])
+ }
+ for (int i = 0; i < num_cameras; ++i) {
+ options.ordering.push_back(my_cameras[i])
+ }
+ options.num_eliminate_blocks = num_points;
+
+
+**After**
+
+.. code-block:: c++
+
+ options.linear_solver_type = ceres::DENSE_SCHUR;
+ options.ordering = new ceres::ParameterBlockOrdering;
+ for (int i = 0; i < num_points; ++i) {
+ options.linear_solver_ordering->AddElementToGroup(my_points[i], 0);
+ }
+ for (int i = 0; i < num_cameras; ++i) {
+ options.linear_solver_ordering->AddElementToGroup(my_cameras[i], 1);
+ }
+
+
+New Features
+------------
+
+#. A new richer, more expressive and consistent API for ordering
+ parameter blocks.
+
+#. A non-linear generalization of Ruhe & Wedin's Algorithm II. This
+ allows the user to use variable projection on separable and
+ non-separable non-linear least squares problems. With
+ multithreading, this results in significant improvements to the
+ convergence behavior of the solver at a small increase in run time.
+
+#. An image denoising example using fields of experts. (Petter
+ Strandmark)
+
+#. Defines for Ceres version and ABI version.
+
+#. Higher precision timer code where available. (Petter Strandmark)
+
+#. Example Makefile for users of Ceres.
+
+#. IterationSummary now informs the user when the step is a
+ non-monotonic step.
+
+#. Fewer memory allocations when using ``DenseQRSolver``.
+
+#. GradientChecker for testing CostFunctions (William Rucklidge)
+
+#. Add support for cost functions with 10 parameter blocks in
+ ``Problem``. (Fisher)
+
+#. Add support for 10 parameter blocks in ``AutoDiffCostFunction``.
+
+
+Bug Fixes
+---------
+
+#. static cast to force Eigen::Index to long conversion
+
+#. Change LOG(ERROR) to LOG(WARNING) in ``schur_complement_solver.cc``.
+
+#. Remove verbose logging from ``DenseQRSolve``.
+
+#. Fix the Android NDK build.
+
+#. Better handling of empty and constant Problems.
+
+#. Remove an internal header that was leaking into the public API.
+
+#. Memory leak in ``trust_region_minimizer.cc``
+
+#. Schur ordering was operating on the wrong object (Ricardo Martin)
+
+#. MSVC fixes (Petter Strandmark)
+
+#. Various fixes to ``nist.cc`` (Markus Moll)
+
+#. Fixed a jacobian scaling bug.
+
+#. Numerically robust computation of ``model_cost_change``.
+
+#. Signed comparison compiler warning fixes (Ricardo Martin)
+
+#. Various compiler warning fixes all over.
+
+#. Inclusion guard fixes (Petter Strandmark)
+
+#. Segfault in test code (Sergey Popov)
+
+#. Replaced ``EXPECT/ASSERT_DEATH`` with the more portable
+ ``EXPECT_DEATH_IF_SUPPORTED`` macros.
+
+#. Fixed the camera projection model in Ceres' implementation of
+ Snavely's camera model. (Ricardo Martin)
+
+
+1.3.0
+=====
+
+New Features
+------------
+
+#. Android Port (Scott Ettinger also contributed to the port)
+
+#. Windows port. (Changchang Wu and Pierre Moulon also contributed to the port)
+
+#. New subspace Dogleg Solver. (Markus Moll)
+
+#. Trust region algorithm now supports the option of non-monotonic steps.
+
+#. New loss functions ``ArcTanLossFunction``, ``TolerantLossFunction``
+ and ``ComposedLossFunction``. (James Roseborough).
+
+#. New ``DENSE_NORMAL_CHOLESKY`` linear solver, which uses Eigen's
+ LDLT factorization on the normal equations.
+
+#. Cached symbolic factorization when using ``CXSparse``.
+ (Petter Strandark)
+
+#. New example ``nist.cc`` and data from the NIST non-linear
+ regression test suite. (Thanks to Douglas Bates for suggesting this.)
+
+#. The traditional Dogleg solver now uses an elliptical trust
+ region (Markus Moll)
+
+#. Support for returning initial and final gradients & Jacobians.
+
+#. Gradient computation support in the evaluators, with an eye
+ towards developing first order/gradient based solvers.
+
+#. A better way to compute ``Solver::Summary::fixed_cost``. (Markus Moll)
+
+#. ``CMake`` support for building documentation, separate examples,
+ installing and uninstalling the library and Gerrit hooks (Arnaud
+ Gelas)
+
+#. ``SuiteSparse4`` support (Markus Moll)
+
+#. Support for building Ceres without ``TR1`` (This leads to
+ slightly slower ``DENSE_SCHUR`` and ``SPARSE_SCHUR`` solvers).
+
+#. ``BALProblem`` can now write a problem back to disk.
+
+#. ``bundle_adjuster`` now allows the user to normalize and perturb the
+ problem before solving.
+
+#. Solver progress logging to file.
+
+#. Added ``Program::ToString`` and ``ParameterBlock::ToString`` to
+ help with debugging.
+
+#. Ability to build Ceres as a shared library (MacOS and Linux only),
+ associated versioning and build release script changes.
+
+#. Portable floating point classification API.
+
+
+Bug Fixes
+---------
+
+#. Fix how invalid step evaluations are handled.
+
+#. Change the slop handling around zero for model cost changes to use
+ relative tolerances rather than absolute tolerances.
+
+#. Fix an inadvertant integer to bool conversion. (Petter Strandmark)
+
+#. Do not link to ``libgomp`` when building on
+ windows. (Petter Strandmark)
+
+#. Include ``gflags.h`` in ``test_utils.cc``. (Petter
+ Strandmark)
+
+#. Use standard random number generation routines. (Petter Strandmark)
+
+#. ``TrustRegionMinimizer`` does not implicitly negate the
+ steps that it takes. (Markus Moll)
+
+#. 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)
+
+#. Added a workaround for a compiler bug in the Android NDK to the
+ Schur eliminator.
+
+#. The sparse linear algebra library is only logged in
+ Summary::FullReport if it is used.
+
+#. Rename the macro ``CERES_DONT_HAVE_PROTOCOL_BUFFERS``
+ to ``CERES_NO_PROTOCOL_BUFFERS`` for consistency.
+
+#. Fix how static structure detection for the Schur eliminator logs
+ its results.
+
+#. Correct example code in the documentation. (Petter Strandmark)
+
+#. Fix ``fpclassify.h`` to work with the Android NDK and STLport.
+
+#. Fix a memory leak in the ``levenber_marquardt_strategy_test.cc``
+
+#. Fix an early return bug in the Dogleg solver. (Markus Moll)
+
+#. Zero initialize Jets.
+#. Moved ``internal/ceres/mock_log.h`` to ``internal/ceres/gmock/mock-log.h``
+
+#. Unified file path handling in tests.
+
+#. ``data_fitting.cc`` includes ``gflags``
+
+#. Renamed Ceres' Mutex class and associated macros to avoid
+ namespace conflicts.
+
+#. Close the BAL problem file after reading it (Markus Moll)
+
+#. Fix IsInfinite on Jets.
+
+#. Drop alignment requirements for Jets.
+
+#. Fixed Jet to integer comparison. (Keith Leung)
+
+#. Fix use of uninitialized arrays. (Sebastian Koch & Markus Moll)
+
+#. Conditionally compile gflag dependencies.(Casey Goodlett)
+
+#. Add ``data_fitting.cc`` to the examples ``CMake`` file.
+
+
+1.2.3
+=====
+
+Bug Fixes
+---------
+
+#. ``suitesparse_test`` is enabled even when ``-DSUITESPARSE=OFF``.
+
+#. ``FixedArray`` internal struct did not respect ``Eigen``
+ alignment requirements (Koichi Akabe & Stephan Kassemeyer).
+
+#. Fixed ``quadratic.cc`` documentation and code mismatch
+ (Nick Lewycky).
+
+1.2.2
+=====
+
+Bug Fixes
+---------
+
+#. Fix constant parameter blocks, and other minor fixes (Markus Moll)
+
+#. Fix alignment issues when combining ``Jet`` and
+ ``FixedArray`` in automatic differeniation.
+
+#. Remove obsolete ``build_defs`` file.
+
+1.2.1
+=====
+
+New Features
+------------
+
+#. Powell's Dogleg solver
+
+#. Documentation now has a brief overview of Trust Region methods and
+ how the Levenberg-Marquardt and Dogleg methods work.
+
+Bug Fixes
+---------
+
+#. Destructor for ``TrustRegionStrategy`` was not virtual (Markus Moll)
+
+#. Invalid ``DCHECK`` in ``suitesparse.cc`` (Markus Moll)
+
+#. Iteration callbacks were not properly invoked (Luis Alberto Zarrabeiti)
+
+#. Logging level changes in ConjugateGradientsSolver
+
+#. VisibilityBasedPreconditioner setup does not account for skipped camera pairs. This was debugging code.
+
+#. Enable SSE support on MacOS
+
+#. ``system_test`` was taking too long and too much memory (Koichi Akabe)
+
+1.2.0
+=====
+
+New Features
+------------
+
+#. ``CXSparse`` support.
+
+#. Block oriented fill reducing orderings. This reduces the
+ factorization time for sparse ``CHOLMOD`` significantly.
+
+#. New Trust region loop with support for multiple trust region step
+ strategies. Currently only Levenberg-Marquardt is supported, but
+ this refactoring opens the door for Dog-leg, Stiehaug and others.
+
+#. ``CMake`` file restructuring. Builds in ``Release`` mode by
+ default, and now has platform specific tuning flags.
+
+#. Re-organized documentation. No new content, but better
+ organization.
+
+
+Bug Fixes
+---------
+
+#. Fixed integer overflow bug in ``block_random_access_sparse_matrix.cc``.
+
+#. Renamed some macros to prevent name conflicts.
+
+#. Fixed incorrent input to ``StateUpdatingCallback``.
+
+#. Fixes to AutoDiff tests.
+
+#. Various internal cleanups.
+
+
+1.1.1
+=====
+
+Bug Fixes
+---------
+
+#. Fix a bug in the handling of constant blocks. (Louis Simard)
+
+#. Add an optional lower bound to the Levenberg-Marquardt regularizer
+ to prevent oscillating between well and ill posed linear problems.
+
+#. Some internal refactoring and test fixes.
+
+1.1.0
+=====
+
+New Features
+------------
+
+#. New iterative linear solver for general sparse problems - ``CGNR``
+ and a block Jacobi preconditioner for it.
+
+#. Changed the semantics of how ``SuiteSparse`` dependencies are
+ checked and used. Now ``SuiteSparse`` is built by default, only if
+ all of its dependencies are present.
+
+#. Automatic differentiation now supports dynamic number of residuals.
+
+#. Support for writing the linear least squares problems to disk in
+ text format so that they can loaded into ``MATLAB``.
+
+#. Linear solver results are now checked for nan and infinities.
+
+#. Added ``.gitignore`` file.
+
+#. A better more robust build system.
+
+
+Bug Fixes
+---------
+
+#. Fixed a strict weak ordering bug in the schur ordering.
+
+#. Grammar and typos in the documents and code comments.
+
+#. Fixed tests which depended on exact equality between floating point values.
+
+1.0.0
+=====
+
+Initial Release. Nathan Wiegand contributed to the Mac OSX port.
