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.
