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.Å. 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.