diff options
Diffstat (limited to 'docs')
36 files changed, 7672 insertions, 107 deletions
diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt index b21a472..1c88478 100644 --- a/docs/CMakeLists.txt +++ b/docs/CMakeLists.txt @@ -1,5 +1,5 @@ # Ceres Solver - A fast non-linear least squares minimizer -# Copyright 2012 Google Inc. All rights reserved. +# Copyright 2013 Google Inc. All rights reserved. # http://code.google.com/p/ceres-solver/ # # Redistribution and use in source and binary forms, with or without @@ -25,42 +25,5 @@ # 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. -# -# Author: arnaudgelas@gmail.com (Arnaud Gelas) - -FIND_PROGRAM(PYGMENTIZE NAMES pygmentize - PATHS - /usr/bin - /usr/local/bin - /opt/bin - /opt/local/bin) - -IF (${PYGMENTIZE} MATCHES PYGMENTIZE-NOTFOUND) - MESSAGE(SEND_ERROR "Pygmentize is not installed") -ENDIF (${PYGMENTIZE} MATCHES PYGMENTIZE-NOTFOUND) - -FIND_PACKAGE(LATEX REQUIRED) - -FILE(DOWNLOAD http://minted.googlecode.com/files/minted.sty - ${CMAKE_BINARY_DIR}/docs/minted.sty - SHOW_PROGRESS) - -ADD_CUSTOM_TARGET(UserGuide ALL - COMMAND ${CMAKE_COMMAND} -E copy_directory - ${CMAKE_SOURCE_DIR}/docs ${CMAKE_BINARY_DIR}/docs - COMMAND ${PDFLATEX_COMPILER} - -shell-escape ${CMAKE_BINARY_DIR}/docs/ceres-solver.tex - -output-directory ${CMAKE_BINARY_DIR}/docs - WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/docs - COMMAND ${BIBTEX_COMPILER} ceres-solver.aux - COMMAND ${PDFLATEX_COMPILER} - -shell-escape ${CMAKE_BINARY_DIR}/docs/ceres-solver.tex - -output-directory ${CMAKE_BINARY_DIR}/docs - WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/docs - COMMAND ${PDFLATEX_COMPILER} - -shell-escape ${CMAKE_BINARY_DIR}/docs/ceres-solver.tex - -output-directory ${CMAKE_BINARY_DIR}/docs - WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/docs) -INSTALL(FILES ${CMAKE_BINARY_DIR}/docs/ceres-solver.pdf - DESTINATION share/ceres/docs) +ADD_SUBDIRECTORY(source) diff --git a/docs/ceres-solver.pdf b/docs/ceres-solver.pdf Binary files differdeleted file mode 100644 index 8f9618c..0000000 --- a/docs/ceres-solver.pdf +++ /dev/null diff --git a/docs/ceres-solver.tex b/docs/ceres-solver.tex index cf5e52d..05cc021 100644 --- a/docs/ceres-solver.tex +++ b/docs/ceres-solver.tex @@ -1,7 +1,7 @@ %%% Build instructions %%% pdflatex -shell-escape ceres && bibtex ceres && pdflatex -shell-escape ceres && pdflatex -shell-escape ceres -\documentclass[11pt,letterpaper,oneside]{memoir} +\documentclass[10pt,letterpaper,oneside]{memoir} \usepackage{fouriernc} \usepackage[T1]{fontenc} \usepackage{minted,amsmath,amssymb,amsthm,url,booktabs} @@ -63,9 +63,9 @@ \MakeLowercase{Ceres Solver: Tutorial \& Reference} } \author{ -\scshape\MakeLowercase{Sameer Agarwal} \\ \texttt{sameeragarwal@google.com} -\and -\scshape\MakeLowercase{Keir Mierle} \\ \texttt{ keir@google.com} +\scshape\MakeLowercase{Sameer Agarwal} \\ \texttt{sameeragarwal@google.com} +\and +\scshape\MakeLowercase{Keir Mierle} \\ \texttt{ mierle@gmail.com} } \checkandfixthelayout @@ -104,9 +104,9 @@ Building this pdf from source requires a relatively recent installation of \texttt{LaTeX}~\footnote{\url{http://www.tug.org/texlive/}}, \texttt{minted.sty}\footnote{\url{http://code.google.com/p/minted/}} and \texttt{pygments}\footnote{\url{http://pygments.org/}}. Despite our best efforts, this manual remains a work in progress and the source code for Ceres Solver remains the ultimate reference. - \input{changes} \input{introduction} +\input{license} \input{build} %% Tutorial diff --git a/docs/ceres.bib b/docs/ceres.bib new file mode 100644 index 0000000..b26a984 --- /dev/null +++ b/docs/ceres.bib @@ -0,0 +1,219 @@ +@article{stigler1981gauss, + title={Gauss and the invention of least squares}, + author={Stigler, S.M.}, + journal={The Annals of Statistics}, + volume={9}, + number={3}, + pages={465--474}, + year={1981}, + publisher={Institute of Mathematical Statistics} +} + +@inproceedings{kushal2012, + title={Visibility Based Preconditioning for Bundle Adjustment}, + author={Kushal, A. and Agarwal, S.}, + booktitle={Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition}, + year={2012}, +} + +@article{nash1990assessing, + title={Assessing a search direction within a truncated-Newton method}, + author={Nash, S.G. and Sofer, A.}, + journal={Operations Research Letters}, + volume={9}, + number={4}, + pages={219--221}, + year={1990} +} + +@inproceedings{wu2011multicore, + title={Multicore bundle adjustment}, + author={Wu, C. and Agarwal, S. and Curless, B. and Seitz, S.M.}, + booktitle={Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition}, + pages={3057--3064}, + year={2011}, +} + +@inproceedings{Agarwal10bal, + author = {Agarwal, S. and Snavely, N. and Seitz, S. M. and Szeliski, R.}, + title = {Bundle adjustment in the large}, +booktitle = {Proceedings of the European Conference on Computer Vision}, + year = {2010}, + pages = {29--42}, +} + +@article{li2007miqr, + title={MIQR: A multilevel incomplete QR preconditioner for large sparse least-squares problems}, + author={Li, Na and Saad, Y.}, + journal={SIAM Journal on Matrix Analysis and Applications}, + volume={28}, + number={2}, + pages={524--550}, + year={2007}, + publisher={Society for Industrial and Applied Mathematics} +} + +@article{wright1985inexact, + Author = {Wright, S. J. and Holt, J. N.}, + Journal = {Journal of the Australian Mathematical Society Series B}, + Pages = "387--403", + Title = {{An Inexact Levenberg-Marquardt Method for Large Sparse Nonlinear Least Squares}}, + Volume = 26, + Number = 4, + Year = 1985 +} + +@article{elsner1984note, + Author = {Elsner, L.}, + Journal = {Numer. Math.}, + Number = {1}, + Pages = {127--128}, + Publisher = {Springer}, + Title = {{A note on optimal block-scaling of matrices}}, + Volume = {44}, + Year = {1984}} + +@article{hestenes1952methods, + Author = {Hestenes, M.R. and Stiefel, E.}, + Journal = {Journal of research of the National Bureau of Standards}, + Number = {6}, + Pages = {409--436}, + Title = {{Methods of conjugate gradients for solving linear systems}}, + Volume = {49}, + Year = {1952}} + +@book{mathew2008domain, + Author = {Mathew, T.P.A.}, + Publisher = {Springer Verlag}, + Title = {{Domain decomposition methods for the numerical solution of partial differential equations}}, + Year = {2008}} + +@book{smith2004domain, + Author = {Smith, B.F. and Bjorstad, P.E. and Gropp, W.}, + Publisher = {Cambridge University Press}, + Title = {{Domain decomposition}}, + Year = {2004}} + +@article{demmel1983condition, + Author = {Demmel, J.}, + Journal = {SINUM}, + Number = {3}, + Pages = {599--610}, + Title = {{The condition number of equivalence transformations that block diagonalize matrix pencils}}, + Volume = {20}, + Year = {1983}} + +@article{eisenstat1982optimal, + Author = {Eisenstat, S.C. and Lewis, J.W. and Schultz, M.H.}, + Journal = {Linear Algebra Appl.}, + Pages = {181--186}, + Title = {{Optimal Block Diagonal Scaling of Block 2-Cyclic Matrices.}}, + Volume = {44}, + Year = {1982}} + +@article{mandel1990block, + Author = {Mandel, J.}, + Journal = {Numer. Math.}, + Number = {1}, + Pages = {79--93}, + Publisher = {Springer}, + Title = {{On block diagonal and Schur complement preconditioning}}, + Volume = {58}, + Year = {1990}} + +@book{davis2006direct, + Author = {Davis, T.A.}, + Publisher = {SIAM}, + Title = {{Direct methods for sparse linear systems}}, + Year = {2006}} + +@TechReport{brown-58, + author = {D. C. Brown}, + title = {A solution to the general problem of multiple station analytical stereo triangulation}, + institution = {Patrick Airforce Base}, + year = {1958}, + Tkey = {AFMTC TR 58-8}, + number = {43}, + address = {Florida}, +} + +@book{hartley-zisserman-book-2004, + Author = {Hartley, R.~I. and Zisserman, A.}, + Publisher = {Cambridge University Press}, + Title = {Multiple View Geometry in Computer Vision}, + Year = {2003}} + + +@book{trefethen1997numerical, + Author = {Trefethen, L.N. and Bau, D.}, + Publisher = {SIAM}, + Title = {{Numerical Linear Algebra}}, + Year = {1997}} + +@book{saad2003iterative, + Author = {Saad, Y.}, + Publisher = {SIAM}, + Title = {{Iterative methods for sparse linear systems}}, + Year = {2003}} + + +@book{nocedal2000numerical, + Author = {Nocedal, J. and Wright, S. J.}, + Publisher = {Springer}, + Title = {{Numerical Optimization}}, + Year = {2000}} + +@book{bjorck1996numerical, + Author = {Bj{\"o}rck, A.}, + Publisher = {SIAM}, + Title = {{Numerical methods for least squares problems}}, + Year = {1996}} + +@book{madsen2004methods, + Author = {Madsen, K. and Nielsen, H.B. and Tingleff, O.}, + Title = {{Methods for non-linear least squares problems}}, + Year = {2004}} + + +@article{marquardt1963algorithm, + Author = {Marquardt, D.W.}, + Journal = {J. SIAM}, + Number = {2}, + Pages = {431--441}, + Publisher = {SIAM}, + Title = {{An algorithm for least-squares estimation of nonlinear parameters}}, + Volume = {11}, + Year = {1963}} + +@article{levenberg1944method, + Author = {Levenberg, K.}, + Journal = {Quart. Appl. Math}, + Number = {2}, + Pages = {164--168}, + Title = {{A method for the solution of certain nonlinear problems in least squares}}, + Volume = {2}, + Year = {1944}} + +@article{chen2006acs, + Article = {22}, + Author = {Chen, Y. and Davis, T. A. and Hager, W. W. and Rajamanickam, S.}, + Journal = {TOMS}, + Number = {3}, + Title = {{Algorithm 887: CHOLMOD, Supernodal Sparse {Cholesky} Factorization and Update/Downdate}}, + Volume = {35}, + Year = {2008}} + + +@inproceedings{triggs-etal-1999, + Author = {Triggs, B. and McLauchlan, P. F. and Hartley, R. I. and Fitzgibbon, A. W.}, + Booktitle = {Vision Algorithms}, + Pages = {298-372}, + Title = {{Bundle Adjustment - A Modern Synthesis}}, + Year = {1999}} + + +@article{tennenbaum-director, +Author = {Tennenbaum, J. and Director, B.}, +Title = {{How Gauss Determined the Orbit of Ceres}} +} + diff --git a/docs/changes.tex b/docs/changes.tex index 159340a..72f7950 100644 --- a/docs/changes.tex +++ b/docs/changes.tex @@ -1,6 +1,23 @@ %!TEX root = ceres-solver.tex \chapter{Version History} +\section*{1.5.0} +\subsection{New Features} +\begin{itemize} +\item 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. +\item Speedup the robust loss function correction logic when residual is one dimensional. +\item Changed \texttt{NumericDiffCostFunction} to take functors like \texttt{AutoDiffCostFunction}. +\item Added support for mixing automatic, analytic and numeric differentiation. This is done by adding \texttt{CostFunctionToFunctor} and \texttt{NumericDiffFunctor} objects. +\end{itemize} + +\subsection{Bug Fixes} +\begin{itemize} +\item Fixed varidic evaluation bug in \texttt{AutoDiff}. +\item Fixed \texttt{SolverImpl} tests. +\item Fixed a bug in \texttt{DenseSparseMatrix::ToDenseMatrix()}. +\item Fixed an initialization bug in \texttt{ProgramEvaluator}. +\end{itemize} + \section*{1.4.0} \subsection{API Changes} The new ordering API breaks existing code. Here the common case fixes. diff --git a/docs/curvefitting.tex b/docs/curvefitting.tex index 07ccd48..c4bacc2 100644 --- a/docs/curvefitting.tex +++ b/docs/curvefitting.tex @@ -17,7 +17,8 @@ class ExponentialResidual { 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]); // $y - e^{mx + c}$ + // $y - e^{mx + c}$ + residual[0] = T(y_) - exp(m[0] * T(x_) + c[0]); return true; } diff --git a/docs/helloworld.tex b/docs/helloworld.tex index 04b677a..938ecba 100644 --- a/docs/helloworld.tex +++ b/docs/helloworld.tex @@ -5,7 +5,7 @@ To get started, let us consider the problem of finding the minimum of the functi \begin{equation} \frac{1}{2}(10 -x)^2. \end{equation} -This is a trivial problem, whose minimum is easy to see is located at $x = 10$, but it is a good place to start to illustrate the basics of solving a problem with Ceres\footnote{Full working code for this and other examples in this manual can be found in the \texttt{examples} directory. Code for this example can be found in \texttt{examples/quadratic.cc}}. +This is a trivial problem, whose minimum is located at $x = 10$, but it is a good place to start to illustrate the basics of solving a problem with Ceres\footnote{Full working code for this and other examples in this manual can be found in the \texttt{examples} directory. Code for this example can be found in \texttt{examples/quadratic.cc}}. Let us write this problem as a non-linear least squares problem by defining the scalar residual function $f_1(x) = 10 - x$. Then $F(x) = [f_1(x)]$ is a residual vector with exactly one component. diff --git a/docs/introduction.tex b/docs/introduction.tex index 65afef3..e9b845b 100644 --- a/docs/introduction.tex +++ b/docs/introduction.tex @@ -38,37 +38,3 @@ Amongst Ceres' users at Google two deserve special mention: William Rucklidge an Nathan Wiegand contributed the MacOS port. -\chapter{License} -Ceres Solver is licensed under the New BSD license, whose terms are as follows. - -\begin{quotation} - -\noindent -Copyright (c) 2010, 2011, 2012, Google Inc. All rights reserved. - -\noindent -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are met: -\begin{enumerate} -\item Redistributions of source code must retain the above copyright notice, - this list of conditions and the following disclaimer. -\item 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. -\item 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. -\end{enumerate} - -\noindent -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. -\end{quotation} diff --git a/docs/license.tex b/docs/license.tex new file mode 100644 index 0000000..602f6f8 --- /dev/null +++ b/docs/license.tex @@ -0,0 +1,35 @@ +%!TEX root = ceres-solver.tex +\chapter{License} +Ceres Solver is licensed under the New BSD license, whose terms are as follows. + +\begin{quotation} + +\noindent +Copyright (c) 2010, 2011, 2012, 2013 Google Inc. All rights reserved. + +\noindent +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: +\begin{enumerate} +\item Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. +\item 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. +\item 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. +\end{enumerate} + +\noindent +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. +\end{quotation} diff --git a/docs/modeling.tex b/docs/modeling.tex index e6dce1f..40fc844 100644 --- a/docs/modeling.tex +++ b/docs/modeling.tex @@ -1,13 +1,33 @@ %!TEX root = ceres-solver.tex \chapter{Modeling} \label{chapter:api} + +\section{Introduction} +Ceres solves robustified non-linear least squares problems of the form +\begin{equation} +\frac{1}{2}\sum_{i=1} \rho_i\left(\left\|f_i\left(x_{i_1},\hdots,x_{i_k}\right)\right\|^2\right). +\label{eq:ceresproblem} +\end{equation} +The term +$\rho_i\left(\left\|f_i\left(x_{i_1},\hdots,x_{i_k}\right)\right\|^2\right)$ +is known as a Residual Block, where $f_i(\cdot)$ is a cost function +that depends on the parameter blocks $\left[x_{i_1}, \hdots , + x_{i_k}\right]$ and $\rho_i$ is a loss function. 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 +\texttt{ParameterBlock}. Of course a \texttt{ParameterBlock} can just +have a single parameter. + \section{\texttt{CostFunction}} +\label{sec:costfunction} Given parameter blocks $\left[x_{i_1}, \hdots , x_{i_k}\right]$, a \texttt{CostFunction} is responsible for computing a vector of residuals and if asked a vector of Jacobian matrices, i.e., given $\left[x_{i_1}, \hdots , x_{i_k}\right]$, compute the vector $f_i\left(x_{i_1},\hdots,x_{i_k}\right)$ and the matrices \begin{equation} -J_{ij} = \frac{\partial}{\partial x_{j}}f_i\left(x_{i_1},\hdots,x_{i_k}\right),\quad \forall j \in \{i_1,\hdots, i_k\} +J_{ij} = \frac{\partial}{\partial x_{i_j}}f_i\left(x_{i_1},\hdots,x_{i_k}\right),\quad \forall j \in \{i_1,\hdots, i_k\} \end{equation} \begin{minted}{c++} class CostFunction { @@ -460,19 +480,16 @@ The cost \texttt{AddParameterBlock} explicitly adds a parameter block to the \texttt{Problem}. Optionally it allows the user to associate a 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 - -\texttt{Problem::SetParameterBlockConstant} - -and undo this using - -\texttt{Problem::SetParameterBlockVariable}. +You can set any parameter block to be constant using \texttt{SetParameterBlockConstant} and undo this using \texttt{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. +\subsubsection{Ownership} \texttt{Problem} by default takes ownership of the \texttt{cost\_function}, \texttt{loss\_function} and \\ \texttt{local\_parameterization} pointers. These objects remain live for the life of the \texttt{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 \texttt{Options} struct. Even though \texttt{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 + do this by setting the corresponding enums in the \texttt{Options} struct. + +Even though \texttt{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. diff --git a/docs/powell.tex b/docs/powell.tex index 5ff6ddc..ad86a42 100644 --- a/docs/powell.tex +++ b/docs/powell.tex @@ -39,7 +39,7 @@ But this can get painful very quickly, especially for residuals involving compli \section{Automatic Differentiation} \label{sec:tutorial:autodiff} 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 \texttt{jacobians} arrays with them. For example, for $f_4(x)$ we define -\begin{minted}[frame=lines,mathescape]{c++} +\begin{minted}[mathescape]{c++} class F4 { public: template <typename T> bool operator()(const T* const x1, @@ -63,7 +63,7 @@ and \texttt{F4}. Then let us consider the construction and solution of the prob \begin{minted}[mathescape]{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. +// wrapper to get the derivatives automatically. problem.AddResidualBlock( new ceres::AutoDiffCostFunction<F1, 1, 1, 1>(new F1), NULL, &x1, &x2); problem.AddResidualBlock( @@ -100,9 +100,30 @@ Final x1 = 0.000583994, x2 = -5.83994e-05, x3 = 9.55401e-05, x4 = 9.55401e-05 It is easy to see that the optimal solution to this problem is at $x_1=0, x_2=0, x_3=0, x_4=0$ with an objective function value of $0$. In 10 iterations, Ceres finds a solution with an objective function value of $4\times 10^{-12}$. \section{Numeric Differentiation} -If a templated implementation is not possible then a \texttt{NumericDiffCostFunction} object can be used. The user defines a \texttt{CostFunction} object whose \texttt{Evaluate} method is only computes the residuals. A wrapper object \texttt{NumericDiffCostFunction} then uses it to compute the residuals and the Jacobian using finite differencing. \texttt{examples/quadratic\_numeric\_diff.cc} shows a numerically differentiated implementation of \texttt{examples/quadratic.cc}. +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 \texttt{NumericDiffCostFunction} using it. e.g., for \texttt{F4}, the corresponding functor would be +\begin{minted}[mathescape]{c++} +class F4 { + public: + bool operator()(const double* const x1, + const double* const x4, + double* residual) const { + // $f_4 = \sqrt{10} * (x_1 - x_4)^2$ + residual[0] = sqrt(10.0) * (x1[0] - x4[0]) * (x1[0] - x4[0]); + return true; + } +}; +\end{minted} + +Which can then be wrapped \texttt{NumericDiffCostFunction} and added to the \texttt{Problem} object as follows + +\begin{minted}[mathescape]{c++} +problem.AddResidualBlock( + new ceres::NumericDiffCostFunction<F4, ceres::CENTRAL, 1, 1, 1>(new F4), NULL, &x1, &x4); +\end{minted} + +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. \texttt{examples/quadratic\_numeric\_diff.cc} shows a numerically differentiated implementation of \texttt{examples/quadratic.cc}. -We recommend that if possible, automatic differentiation should be used. The use of -C++ templates makes automatic differentiation extremely efficient, +\textbf{We recommend that if possible, automatic differentiation should be used. The use of +\texttt{C++} templates makes automatic differentiation extremely efficient, whereas numeric differentiation can be quite expensive, prone to -numeric errors and leads to slower convergence.
\ No newline at end of file +numeric errors and leads to slower convergence.} diff --git a/docs/reference-overview.tex b/docs/reference-overview.tex index dac0ef2..23ab82b 100644 --- a/docs/reference-overview.tex +++ b/docs/reference-overview.tex @@ -1,18 +1,18 @@ %!TEX root = ceres-solver.tex \chapter{Overview} \label{chapter:overview} -Ceres solves robustified non-linear least squares problems of the form -\begin{equation} -\frac{1}{2}\sum_{i=1} \rho_i\left(\left\|f_i\left(x_{i_1},\hdots,x_{i_k}\right)\right\|^2\right). -\label{eq:ceresproblem} -\end{equation} -Where $f_i(\cdot)$ is a cost function that depends on the parameter blocks $\left[x_{i_1}, \hdots , x_{i_k}\right]$ and $\rho_i$ is a loss function. 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 Parameter Block. Of course a parameter block can just have a single parameter. -The term $ \rho_i\left(\left\|f_i\left(x_{i_1},\hdots,x_{i_k}\right)\right\|^2\right)$ is known as a Residual Block. A Ceres problem is a collection of residual blocks, each of which depends on a subset of the parameter blocks. Solving problems using Ceres consists of two steps. -\begin{enumerate} -\item{Modeling} Define parameter blocks and residual blocks and build a \texttt{Problem} object containing them. -\item{Solving} Configure and run the solver. -\end{enumerate} +\begin{description} +\item{\textbf{Modeling}} Constructing an optimization problem by + specifying its parameters and the terms in the objective function. +\item{\textbf{Solving}} Configuring and running the solver. +\end{description} -These two steps are mostly independent of each other. This is by design. Modeling the optimization problem should not depend on how the solver and the user should be able to switch between various solver settings and strategies without changing the way the problem is modeled. In the next two chapters we will consider each of these steps in detail.
\ No newline at end of file +The two steps are mostly independent of each other. This is by +design. Modeling the optimization problem should not depend on how the +solver works. The user should be able model the problem once, and then +switch between various solver settings and strategies without touching +the problem. + +In the next two chapters we will consider each of these steps in detail. diff --git a/docs/source/CMakeLists.txt b/docs/source/CMakeLists.txt new file mode 100644 index 0000000..91b76df --- /dev/null +++ b/docs/source/CMakeLists.txt @@ -0,0 +1,19 @@ +FIND_PACKAGE(Sphinx REQUIRED) + +# HTML output directory +SET(SPHINX_HTML_DIR "${CMAKE_BINARY_DIR}/docs/html") + +# Install documentation +INSTALL(DIRECTORY ${SPHINX_HTML_DIR} + DESTINATION share/doc/ceres + COMPONENT Doc + PATTERN "${SPHINX_HTML_DIR}/*") + +# Building using 'make_docs.py' python script +ADD_CUSTOM_TARGET(ceres_docs ALL + python + "${CMAKE_SOURCE_DIR}/scripts/make_docs.py" + "${CMAKE_SOURCE_DIR}" + "${CMAKE_BINARY_DIR}/docs" + "${SPHINX_EXECUTABLE}" + COMMENT "Building HTML documentation with Sphinx") diff --git a/docs/source/_themes/armstrong/LICENSE b/docs/source/_themes/armstrong/LICENSE new file mode 100644 index 0000000..894aa01 --- /dev/null +++ b/docs/source/_themes/armstrong/LICENSE @@ -0,0 +1,26 @@ +Copyright (c) 2011 Bay Citizen & Texas Tribune + +Original ReadTheDocs.org code +Copyright (c) 2010 Charles Leifer, Eric Holscher, Bobby Grace + +Permission is hereby granted, free of charge, to any person +obtaining a copy of this software and associated documentation +files (the "Software"), to deal in the Software without +restriction, including without limitation the rights to use, +copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the +Software is furnished to do so, subject to the following +conditions: + +The above copyright notice and this permission notice shall be +included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +OTHER DEALINGS IN THE SOFTWARE. + diff --git a/docs/source/_themes/armstrong/globaltoc.html b/docs/source/_themes/armstrong/globaltoc.html new file mode 100644 index 0000000..20d8641 --- /dev/null +++ b/docs/source/_themes/armstrong/globaltoc.html @@ -0,0 +1,11 @@ +{# + basic/globaltoc.html + ~~~~~~~~~~~~~~~~~~~~ + + Sphinx sidebar template: global table of contents. + + :copyright: Copyright 2007-2011 by the Sphinx team, see AUTHORS. + :license: BSD, see LICENSE for details. +#} +<h3><a href="{{ pathto(master_doc) }}">{{ _('Ceres Solver') }}</a></h3> +{{ toctree() }} diff --git a/docs/source/_themes/armstrong/layout.html b/docs/source/_themes/armstrong/layout.html new file mode 100644 index 0000000..3faa34c --- /dev/null +++ b/docs/source/_themes/armstrong/layout.html @@ -0,0 +1,80 @@ +{% extends "basic/layout.html" %} + +{% set script_files = script_files + [pathto("_static/searchtools.js", 1)] %} + +{% block htmltitle %} +{{ super() }} + +<meta name="viewport" content="width=device-width; initial-scale=1.0; maximum-scale=1.0; user-scalable=0;"/> + +{% endblock %} + + +{%- macro sidebar() %} + {%- if render_sidebar %} + <div class="sphinxsidebar"> + <div class="sphinxsidebarwrapper"> + {%- block sidebarlogo %} + {%- if logo %} + <p class="logo"><a href="{{ pathto(master_doc) }}"> + <img class="logo" src="{{ pathto('_static/' + logo, 1) }}" alt="Logo"/> + </a></p> + {%- endif %} + {%- endblock %} + {%- if sidebars != None %} + {#- new style sidebar: explicitly include/exclude templates #} + {%- for sidebartemplate in sidebars %} + {%- include sidebartemplate %} + {%- endfor %} + {%- else %} + {#- old style sidebars: using blocks -- should be deprecated #} + {%- block sidebartoc %} + {%- include "globaltoc.html" %} + {%- endblock %} + {%- block sidebarsourcelink %} + {%- include "sourcelink.html" %} + {%- endblock %} + {%- if customsidebar %} + {%- include customsidebar %} + {%- endif %} + {%- block sidebarsearch %} + {%- include "searchbox.html" %} + {%- endblock %} + {%- endif %} + </div> + </div> + {%- endif %} +{%- endmacro %} + + +{% block footer %} +<div class="footer"> +{%- if show_copyright %} + {%- if hasdoc('copyright') %} + {% trans path=pathto('copyright'), copyright=copyright|e %}© <a href="{{ path }}">Copyright</a> {{ copyright }}.{% endtrans %} + {%- else %} + {% trans copyright=copyright|e %}© Copyright {{ copyright }}.{% endtrans %} + {%- endif %} +{%- endif %} +{%- if last_updated %} + {% trans last_updated=last_updated|e %}Last updated on {{ last_updated }}.{% endtrans %} +{%- endif %} +</div> + + +{% if theme_analytics_code %} +<!-- Google Analytics Code --> +<script type="text/javascript"> + var _gaq = _gaq || []; + _gaq.push(['_setAccount', '{{ theme_analytics_code }}']); + _gaq.push(['_trackPageview']); + + (function() { + var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true; + ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js'; + var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s); + })(); +</script> +{% endif %} + +{% endblock %} diff --git a/docs/source/_themes/armstrong/rtd-themes.conf b/docs/source/_themes/armstrong/rtd-themes.conf new file mode 100644 index 0000000..5930488 --- /dev/null +++ b/docs/source/_themes/armstrong/rtd-themes.conf @@ -0,0 +1,65 @@ +[theme] +inherit = default +stylesheet = rtd.css +pygment_style = default +show_sphinx = False + +[options] +show_rtd = True + +white = #ffffff +almost_white = #f8f8f8 +barely_white = #f2f2f2 +dirty_white = #eeeeee +almost_dirty_white = #e6e6e6 +dirtier_white = #dddddd +lighter_gray = #cccccc +gray_a = #aaaaaa +gray_9 = #999999 +light_gray = #888888 +gray_7 = #777777 +gray = #666666 +dark_gray = #444444 +gray_2 = #222222 +black = #111111 +light_color = #e8ecef +light_medium_color = #DDEAF0 +medium_color = #8ca1af +medium_color_link = #86989b +medium_color_link_hover = #a6b8bb +dark_color = #465158 + +h1 = #000000 +h2 = #465158 +h3 = #6c818f + +link_color = #444444 +link_color_decoration = #CCCCCC + +medium_color_hover = #697983 +green_highlight = #8ecc4c + + +positive_dark = #609060 +positive_medium = #70a070 +positive_light = #e9ffe9 + +negative_dark = #900000 +negative_medium = #b04040 +negative_light = #ffe9e9 +negative_text = #c60f0f + +ruler = #abc + +viewcode_bg = #f4debf +viewcode_border = #ac9 + +highlight = #ffe080 + +code_background = #eeeeee + +background = #465158 +background_link = #ffffff +background_link_half = #ffffff +background_text = #eeeeee +background_text_link = #86989b diff --git a/docs/source/_themes/armstrong/static/rtd.css_t b/docs/source/_themes/armstrong/static/rtd.css_t new file mode 100644 index 0000000..90354c3 --- /dev/null +++ b/docs/source/_themes/armstrong/static/rtd.css_t @@ -0,0 +1,781 @@ +/* + * rtd.css + * ~~~~~~~~~~~~~~~ + * + * Sphinx stylesheet -- sphinxdoc theme. Originally created by + * Armin Ronacher for Werkzeug. + * + * Customized for ReadTheDocs by Eric Pierce & Eric Holscher + * + * :copyright: Copyright 2007-2010 by the Sphinx team, see AUTHORS. + * :license: BSD, see LICENSE for details. + * + */ + +/* RTD colors + * light blue: {{ theme_light_color }} + * medium blue: {{ theme_medium_color }} + * dark blue: {{ theme_dark_color }} + * dark grey: {{ theme_grey_color }} + * + * medium blue hover: {{ theme_medium_color_hover }}; + * green highlight: {{ theme_green_highlight }} + * light blue (project bar): {{ theme_light_color }} + */ + +@import url("basic.css"); + +/* PAGE LAYOUT -------------------------------------------------------------- */ + +body { + font: 100%/1.5 "ff-meta-web-pro-1","ff-meta-web-pro-2",Arial,"Helvetica Neue",sans-serif; + text-align: center; + color: black; + background-color: {{ theme_background }}; + padding: 0; + margin: 0; +} + +div.document { + text-align: left; + background-color: {{ theme_light_color }}; +} + +div.bodywrapper { + background-color: {{ theme_white }}; + border-left: 1px solid {{ theme_lighter_gray }}; + border-bottom: 1px solid {{ theme_lighter_gray }}; + margin: 0 0 0 16em; +} + +div.body { + margin: 0; + padding: 0.5em 1.3em; + max-width: 55em; + min-width: 20em; +} + +div.related { + font-size: 1em; + background-color: {{ theme_background }}; +} + +div.documentwrapper { + float: left; + width: 100%; + background-color: {{ theme_light_color }}; +} + + +/* HEADINGS --------------------------------------------------------------- */ + +h1 { + margin: 0; + padding: 0.7em 0 0.3em 0; + font-size: 1.5em; + line-height: 1.15; + color: {{ theme_h1 }}; + clear: both; +} + +h2 { + margin: 2em 0 0.2em 0; + font-size: 1.35em; + padding: 0; + color: {{ theme_h2 }}; +} + +h3 { + margin: 1em 0 -0.3em 0; + font-size: 1.2em; + color: {{ theme_h3 }}; +} + +div.body h1 a, div.body h2 a, div.body h3 a, div.body h4 a, div.body h5 a, div.body h6 a { + color: black; +} + +h1 a.anchor, h2 a.anchor, h3 a.anchor, h4 a.anchor, h5 a.anchor, h6 a.anchor { + display: none; + margin: 0 0 0 0.3em; + padding: 0 0.2em 0 0.2em; + color: {{ theme_gray_a }} !important; +} + +h1:hover a.anchor, h2:hover a.anchor, h3:hover a.anchor, h4:hover a.anchor, +h5:hover a.anchor, h6:hover a.anchor { + display: inline; +} + +h1 a.anchor:hover, h2 a.anchor:hover, h3 a.anchor:hover, h4 a.anchor:hover, +h5 a.anchor:hover, h6 a.anchor:hover { + color: {{ theme_gray_7 }}; + background-color: {{ theme_dirty_white }}; +} + + +/* LINKS ------------------------------------------------------------------ */ + +/* Normal links get a pseudo-underline */ +a { + color: {{ theme_link_color }}; + text-decoration: none; + border-bottom: 1px solid {{ theme_link_color_decoration }}; +} + +/* Links in sidebar, TOC, index trees and tables have no underline */ +.sphinxsidebar a, +.toctree-wrapper a, +.indextable a, +#indices-and-tables a { + color: {{ theme_dark_gray }}; + text-decoration: none; + border-bottom: none; +} + +/* Most links get an underline-effect when hovered */ +a:hover, +div.toctree-wrapper a:hover, +.indextable a:hover, +#indices-and-tables a:hover { + color: {{ theme_black }}; + text-decoration: none; + border-bottom: 1px solid {{ theme_black }}; +} + +/* Footer links */ +div.footer a { + color: {{ theme_background_text_link }}; + text-decoration: none; + border: none; +} +div.footer a:hover { + color: {{ theme_medium_color_link_hover }}; + text-decoration: underline; + border: none; +} + +/* Permalink anchor (subtle grey with a red hover) */ +div.body a.headerlink { + color: {{ theme_lighter_gray }}; + font-size: 1em; + margin-left: 6px; + padding: 0 4px 0 4px; + text-decoration: none; + border: none; +} +div.body a.headerlink:hover { + color: {{ theme_negative_text }}; + border: none; +} + + +/* NAVIGATION BAR --------------------------------------------------------- */ + +div.related ul { + height: 2.5em; +} + +div.related ul li { + margin: 0; + padding: 0.65em 0; + float: left; + display: block; + color: {{ theme_background_link_half }}; /* For the >> separators */ + font-size: 0.8em; +} + +div.related ul li.right { + float: right; + margin-right: 5px; + color: transparent; /* Hide the | separators */ +} + +/* "Breadcrumb" links in nav bar */ +div.related ul li a { + order: none; + background-color: inherit; + font-weight: bold; + margin: 6px 0 6px 4px; + line-height: 1.75em; + color: {{ theme_background_link }}; + text-shadow: 0 1px rgba(0, 0, 0, 0.5); + padding: 0.4em 0.8em; + border: none; + border-radius: 3px; +} +/* previous / next / modules / index links look more like buttons */ +div.related ul li.right a { + margin: 0.375em 0; + background-color: {{ theme_medium_color_hover }}; + text-shadow: 0 1px rgba(0, 0, 0, 0.5); + border-radius: 3px; + -webkit-border-radius: 3px; + -moz-border-radius: 3px; +} +/* All navbar links light up as buttons when hovered */ +div.related ul li a:hover { + background-color: {{ theme_medium_color }}; + color: {{ theme_white }}; + text-decoration: none; + border-radius: 3px; + -webkit-border-radius: 3px; + -moz-border-radius: 3px; +} +/* Take extra precautions for tt within links */ +a tt, +div.related ul li a tt { + background: inherit !important; + color: inherit !important; +} + + +/* SIDEBAR ---------------------------------------------------------------- */ + +div.sphinxsidebarwrapper { + padding: 0; +} + +div.sphinxsidebar { + margin: 0; + margin-left: -100%; + float: left; + top: 3em; + left: 0; + padding: 0 1em; + width: 14em; + font-size: 1em; + text-align: left; + background-color: {{ theme_light_color }}; +} + +div.sphinxsidebar img { + max-width: 12em; +} + +div.sphinxsidebar h3, div.sphinxsidebar h4 { + margin: 1.2em 0 0.3em 0; + font-size: 1em; + padding: 0; + color: {{ theme_gray_2 }}; + font-family: "ff-meta-web-pro-1", "ff-meta-web-pro-2", "Arial", "Helvetica Neue", sans-serif; +} + +div.sphinxsidebar h3 a { + color: {{ theme_grey_color }}; +} + +div.sphinxsidebar ul, +div.sphinxsidebar p { + margin-top: 0; + padding-left: 0; + line-height: 130%; + background-color: {{ theme_light_color }}; +} + +/* No bullets for nested lists, but a little extra indentation */ +div.sphinxsidebar ul ul { + list-style-type: none; + margin-left: 1.5em; + padding: 0; +} + +/* A little top/bottom padding to prevent adjacent links' borders + * from overlapping each other */ +div.sphinxsidebar ul li { + padding: 1px 0; +} + +/* A little left-padding to make these align with the ULs */ +div.sphinxsidebar p.topless { + padding-left: 0 0 0 1em; +} + +/* Make these into hidden one-liners */ +div.sphinxsidebar ul li, +div.sphinxsidebar p.topless { + white-space: nowrap; + overflow: hidden; +} +/* ...which become visible when hovered */ +div.sphinxsidebar ul li:hover, +div.sphinxsidebar p.topless:hover { + overflow: visible; +} + +/* Search text box and "Go" button */ +#searchbox { + margin-top: 2em; + margin-bottom: 1em; + background: {{ theme_dirtier_white }}; + padding: 0.5em; + border-radius: 6px; + -moz-border-radius: 6px; + -webkit-border-radius: 6px; +} +#searchbox h3 { + margin-top: 0; +} + +/* Make search box and button abut and have a border */ +input, +div.sphinxsidebar input { + border: 1px solid {{ theme_gray_9 }}; + float: left; +} + +/* Search textbox */ +input[type="text"] { + margin: 0; + padding: 0 3px; + height: 20px; + width: 144px; + border-top-left-radius: 3px; + border-bottom-left-radius: 3px; + -moz-border-radius-topleft: 3px; + -moz-border-radius-bottomleft: 3px; + -webkit-border-top-left-radius: 3px; + -webkit-border-bottom-left-radius: 3px; +} +/* Search button */ +input[type="submit"] { + margin: 0 0 0 -1px; /* -1px prevents a double-border with textbox */ + height: 22px; + color: {{ theme_dark_gray }}; + background-color: {{ theme_light_color }}; + padding: 1px 4px; + font-weight: bold; + border-top-right-radius: 3px; + border-bottom-right-radius: 3px; + -moz-border-radius-topright: 3px; + -moz-border-radius-bottomright: 3px; + -webkit-border-top-right-radius: 3px; + -webkit-border-bottom-right-radius: 3px; +} +input[type="submit"]:hover { + color: {{ theme_white }}; + background-color: {{ theme_green_highlight }}; +} + +div.sphinxsidebar p.searchtip { + clear: both; + padding: 0.5em 0 0 0; + background: {{ theme_dirtier_white }}; + color: {{ theme_gray }}; + font-size: 0.9em; +} + +/* Sidebar links are unusual */ +div.sphinxsidebar li a, +div.sphinxsidebar p a { + background: {{ theme_light_color }}; /* In case links overlap main content */ + border-radius: 3px; + -moz-border-radius: 3px; + -webkit-border-radius: 3px; + border: 1px solid transparent; /* To prevent things jumping around on hover */ + padding: 0 5px 0 5px; +} +div.sphinxsidebar li a:hover, +div.sphinxsidebar p a:hover { + color: {{ theme_black }}; + text-decoration: none; + border: 1px solid {{ theme_light_gray }}; +} + +/* Tweak any link appearing in a heading */ +div.sphinxsidebar h3 a { +} + + + + +/* OTHER STUFF ------------------------------------------------------------ */ + +cite, code, tt { + font-family: 'Consolas', 'Deja Vu Sans Mono', + 'Bitstream Vera Sans Mono', monospace; + font-size: 0.95em; + letter-spacing: 0.01em; +} + +tt { + background-color: {{ theme_code_background }}; + color: {{ theme_dark_gray }}; +} + +tt.descname, tt.descclassname, tt.xref { + border: 0; +} + +hr { + border: 1px solid {{ theme_ruler }}; + margin: 2em; +} + +pre, #_fontwidthtest { + font-family: 'Consolas', 'Deja Vu Sans Mono', + 'Bitstream Vera Sans Mono', monospace; + margin: 1em 2em; + font-size: 0.95em; + letter-spacing: 0.015em; + line-height: 120%; + padding: 0.5em; + border: 1px solid {{ theme_lighter_gray }}; + background-color: {{ theme_code_background }}; + border-radius: 6px; + -moz-border-radius: 6px; + -webkit-border-radius: 6px; +} + +pre a { + color: inherit; + text-decoration: underline; +} + +td.linenos pre { + padding: 0.5em 0; +} + +div.quotebar { + background-color: {{ theme_almost_white }}; + max-width: 250px; + float: right; + padding: 2px 7px; + border: 1px solid {{ theme_lighter_gray }}; +} + +div.topic { + background-color: {{ theme_almost_white }}; +} + +table { + border-collapse: collapse; + margin: 0 -0.5em 0 -0.5em; +} + +table td, table th { + padding: 0.2em 0.5em 0.2em 0.5em; +} + + +/* ADMONITIONS AND WARNINGS ------------------------------------------------- */ + +/* Shared by admonitions, warnings and sidebars */ +div.admonition, +div.warning, +div.sidebar { + font-size: 0.9em; + margin: 2em; + padding: 0; + /* + border-radius: 6px; + -moz-border-radius: 6px; + -webkit-border-radius: 6px; + */ +} +div.admonition p, +div.warning p, +div.sidebar p { + margin: 0.5em 1em 0.5em 1em; + padding: 0; +} +div.admonition pre, +div.warning pre, +div.sidebar pre { + margin: 0.4em 1em 0.4em 1em; +} +div.admonition p.admonition-title, +div.warning p.admonition-title, +div.sidebar p.sidebar-title { + margin: 0; + padding: 0.1em 0 0.1em 0.5em; + color: white; + font-weight: bold; + font-size: 1.1em; + text-shadow: 0 1px rgba(0, 0, 0, 0.5); +} +div.admonition ul, div.admonition ol, +div.warning ul, div.warning ol, +div.sidebar ul, div.sidebar ol { + margin: 0.1em 0.5em 0.5em 3em; + padding: 0; +} + + +/* Admonitions and sidebars only */ +div.admonition, div.sidebar { + border: 1px solid {{ theme_positive_dark }}; + background-color: {{ theme_positive_light }}; +} +div.admonition p.admonition-title, +div.sidebar p.sidebar-title { + background-color: {{ theme_positive_medium }}; + border-bottom: 1px solid {{ theme_positive_dark }}; +} + + +/* Warnings only */ +div.warning { + border: 1px solid {{ theme_negative_dark }}; + background-color: {{ theme_negative_light }}; +} +div.warning p.admonition-title { + background-color: {{ theme_negative_medium }}; + border-bottom: 1px solid {{ theme_negative_dark }}; +} + + +/* Sidebars only */ +div.sidebar { + max-width: 200px; +} + + + +div.versioninfo { + margin: 1em 0 0 0; + border: 1px solid {{ theme_lighter_gray }}; + background-color: {{ theme_light_medium_color }}; + padding: 8px; + line-height: 1.3em; + font-size: 0.9em; +} + +.viewcode-back { + font-family: 'Lucida Grande', 'Lucida Sans Unicode', 'Geneva', + 'Verdana', sans-serif; +} + +div.viewcode-block:target { + background-color: {{ theme_viewcode_bg }}; + border-top: 1px solid {{ theme_viewcode_border }}; + border-bottom: 1px solid {{ theme_viewcode_border }}; +} + +dl { + margin: 1em 0 2.5em 0; +} + +/* Highlight target when you click an internal link */ +dt:target { + background: {{ theme_highlight }}; +} +/* Don't highlight whole divs */ +div.highlight { + background: transparent; +} +/* But do highlight spans (so search results can be highlighted) */ +span.highlight { + background: {{ theme_highlight }}; +} + +div.footer { + background-color: {{ theme_background }}; + color: {{ theme_background_text }}; + padding: 0 2em 2em 2em; + clear: both; + font-size: 0.8em; + text-align: center; +} + +p { + margin: 0.8em 0 0.5em 0; +} + +.section p img { + margin: 1em 2em; +} + + +/* MOBILE LAYOUT -------------------------------------------------------------- */ + +@media screen and (max-width: 600px) { + + h1, h2, h3, h4, h5 { + position: relative; + } + + ul { + padding-left: 1.75em; + } + + div.bodywrapper a.headerlink, #indices-and-tables h1 a { + color: {{ theme_almost_dirty_white }}; + font-size: 80%; + float: right; + line-height: 1.8; + position: absolute; + right: -0.7em; + visibility: inherit; + } + + div.bodywrapper h1 a.headerlink, #indices-and-tables h1 a { + line-height: 1.5; + } + + pre { + font-size: 0.7em; + overflow: auto; + word-wrap: break-word; + white-space: pre-wrap; + } + + div.related ul { + height: 2.5em; + padding: 0; + text-align: left; + } + + div.related ul li { + clear: both; + color: {{ theme_dark_color }}; + padding: 0.2em 0; + } + + div.related ul li:last-child { + border-bottom: 1px dotted {{ theme_medium_color }}; + padding-bottom: 0.4em; + margin-bottom: 1em; + width: 100%; + } + + div.related ul li a { + color: {{ theme_dark_color }}; + padding-right: 0; + } + + div.related ul li a:hover { + background: inherit; + color: inherit; + } + + div.related ul li.right { + clear: none; + padding: 0.65em 0; + margin-bottom: 0.5em; + } + + div.related ul li.right a { + color: {{ theme_white }}; + padding-right: 0.8em; + } + + div.related ul li.right a:hover { + background-color: {{ theme_medium_color }}; + } + + div.body { + clear: both; + min-width: 0; + word-wrap: break-word; + } + + div.bodywrapper { + margin: 0 0 0 0; + } + + div.sphinxsidebar { + float: none; + margin: 0; + width: auto; + } + + div.sphinxsidebar input[type="text"] { + height: 2em; + line-height: 2em; + width: 70%; + } + + div.sphinxsidebar input[type="submit"] { + height: 2em; + margin-left: 0.5em; + width: 20%; + } + + div.sphinxsidebar p.searchtip { + background: inherit; + margin-bottom: 1em; + } + + div.sphinxsidebar ul li, div.sphinxsidebar p.topless { + white-space: normal; + } + + .bodywrapper img { + display: block; + margin-left: auto; + margin-right: auto; + max-width: 100%; + } + + div.documentwrapper { + float: none; + } + + div.admonition, div.warning, pre, blockquote { + margin-left: 0em; + margin-right: 0em; + } + + .body p img { + margin: 0; + } + + #searchbox { + background: transparent; + } + + .related:not(:first-child) li { + display: none; + } + + .related:not(:first-child) li.right { + display: block; + } + + div.footer { + padding: 1em; + } + + .rtd_doc_footer .badge { + float: none; + margin: 1em auto; + position: static; + } + + .rtd_doc_footer .badge.revsys-inline { + margin-right: auto; + margin-bottom: 2em; + } + + table.indextable { + display: block; + width: auto; + } + + .indextable tr { + display: block; + } + + .indextable td { + display: block; + padding: 0; + width: auto !important; + } + + .indextable td dt { + margin: 1em 0; + } + + ul.search { + margin-left: 0.25em; + } + + ul.search li div.context { + font-size: 90%; + line-height: 1.1; + margin-bottom: 1; + margin-left: 0; + } + +} diff --git a/docs/source/_themes/armstrong/theme.conf b/docs/source/_themes/armstrong/theme.conf new file mode 100644 index 0000000..5930488 --- /dev/null +++ b/docs/source/_themes/armstrong/theme.conf @@ -0,0 +1,65 @@ +[theme] +inherit = default +stylesheet = rtd.css +pygment_style = default +show_sphinx = False + +[options] +show_rtd = True + +white = #ffffff +almost_white = #f8f8f8 +barely_white = #f2f2f2 +dirty_white = #eeeeee +almost_dirty_white = #e6e6e6 +dirtier_white = #dddddd +lighter_gray = #cccccc +gray_a = #aaaaaa +gray_9 = #999999 +light_gray = #888888 +gray_7 = #777777 +gray = #666666 +dark_gray = #444444 +gray_2 = #222222 +black = #111111 +light_color = #e8ecef +light_medium_color = #DDEAF0 +medium_color = #8ca1af +medium_color_link = #86989b +medium_color_link_hover = #a6b8bb +dark_color = #465158 + +h1 = #000000 +h2 = #465158 +h3 = #6c818f + +link_color = #444444 +link_color_decoration = #CCCCCC + +medium_color_hover = #697983 +green_highlight = #8ecc4c + + +positive_dark = #609060 +positive_medium = #70a070 +positive_light = #e9ffe9 + +negative_dark = #900000 +negative_medium = #b04040 +negative_light = #ffe9e9 +negative_text = #c60f0f + +ruler = #abc + +viewcode_bg = #f4debf +viewcode_border = #ac9 + +highlight = #ffe080 + +code_background = #eeeeee + +background = #465158 +background_link = #ffffff +background_link_half = #ffffff +background_text = #eeeeee +background_text_link = #86989b diff --git a/docs/source/acknowledgements.rst b/docs/source/acknowledgements.rst new file mode 100644 index 0000000..36c1562 --- /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 the :ref:`chapter-version-history`. diff --git a/docs/source/bibliography.rst b/docs/source/bibliography.rst new file mode 100644 index 0000000..7ba435a --- /dev/null +++ b/docs/source/bibliography.rst @@ -0,0 +1,119 @@ +.. _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. + +.. [ByrdNocedal] R. H. Byrd, J. Nocedal, R. B. Schanbel, + **Representations of Quasi-Newton Matrices and their use in Limited + Memory Methods**, *Mathematical Programming* 63(4):129–-156, 1994. + +.. [ByrdSchnabel] 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. + +.. [KanataniMorris] K. Kanatani and D. D. Morris, **Gauges and gauge + transformations for uncertainty description of geometric structure + with indeterminacy**, *IEEE Transactions on Information Theory* + 47(5):2017-2028, 2001. + +.. [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. + +.. [Nocedal] J. Nocedal, **Updating Quasi-Newton Matrices with Limited + Storage**, *Mathematics of Computation*, 35(151): 773--782, 1980. + +.. [NocedalWright] J. Nocedal & S. Wright, **Numerical Optimization**, + Springer, 2004. + +.. [Oren] S. S. Oren, **Self-scaling Variable Metric (SSVM) Algorithms + Part II: Implementation and Experiments**, Management Science, + 20(5), 863-874, 1974. + +.. [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..c326fd1 --- /dev/null +++ b/docs/source/building.rst @@ -0,0 +1,392 @@ +.. _chapter-building: + +===================== +Building Ceres Solver +===================== + +Stable Ceres Solver releases are available for download at +`code.google.com <http://code.google.com/p/ceres-solver/>`_. For the +more adventurous, the git repository is hosted on `Gerrit +<https://ceres-solver-review.googlesource.com/>`_. + +.. _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://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, CAMD, COLAMD and CHOLMOD libraries. This is an optional +dependency. + +6. `CXSparse <http://www.cise.ufl.edu/research/sparse/CXSparse/>`_ is +a sparse matrix library similar in scope to ``SuiteSparse`` but with +no dependencies on ``LAPACK`` and ``BLAS``. This makes for a simpler +build process and a smaller binary. The simplicity comes at a cost -- +for all but the most trivial matrices, ``SuiteSparse`` is +significantly faster than ``CXSparse``. + +7. `BLAS <http://www.netlib.org/blas/>`_ and `LAPACK +<http://www.netlib.org/lapack/>`_ routines are needed by +SuiteSparse. We recommend `ATLAS +<http://math-atlas.sourceforge.net/>`_, which includes BLAS and LAPACK +routines. It is also possible to use `OpenBLAS +<https://github.com/xianyi/OpenBLAS>`_ . However, one needs to be +careful to `turn off the threading +<https://github.com/xianyi/OpenBLAS/wiki/faq#wiki-multi-threaded>`_ +inside ``OpenBLAS`` as it conflicts with use of threads in Ceres. + +.. _section-linux: + +Building on Linux +================= +We will use `Ubuntu <http://www.ubuntu.com>`_ as our example +platform. Start by installing all the dependencies. + +.. code-block:: bash + + # CMake + sudo apt-get install cmake + # gflags + tar -xvzf gflags-2.0.tar.gz + cd gflags-2.0 + ./configure --prefix=/usr/local + make + sudo make install. + # google-glog must be configured to use the previously installed gflags + tar -xvzf glog-0.3.2.tar.gz + cd glog-0.3.2 + ./configure --with-gflags=/usr/local/ + make + sudo make install + # BLAS & LAPACK + sudo apt-get install libatlas-base-dev + # Eigen3 + sudo apt-get install libeigen3-dev + # SuiteSparse and CXSparse + sudo apt-get install libsuitesparse-dev + +We are now ready to build and test Ceres. + +.. code-block:: bash + + tar zxf ceres-solver-1.7.0.tar.gz + mkdir ceres-bin + cd ceres-bin + cmake ../ceres-solver-1.7.0 + 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.7.0/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 to install the +dependencies. There is no need to install ``BLAS`` or ``LAPACK`` +separately as OS X ships with optimized ``BLAS`` and ``LAPACK`` +routines as part of the `vecLib +<https://developer.apple.com/library/mac/#documentation/Performance/Conceptual/vecLib/Reference/reference.html>`_ +framework. + +.. code-block:: bash + + # CMake + brew install cmake + # google-glog and gflags + brew install glog + # Eigen3 + brew install eigen + # SuiteSparse and CXSparse + brew install suite-sparse + + +We are now ready to build and test Ceres. + +.. code-block:: bash + + tar zxf ceres-solver-1.7.0.tar.gz + mkdir ceres-bin + cd ceres-bin + cmake ../ceres-solver-1.7.0 + 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 unavailability 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. + +#. ``-DSUITESPARSE=OFF``: By default, Ceres will link to + ``SuiteSparse`` if all its dependencies are present. Use this flag + to build Ceres without ``SuiteSparse``. This will also disable + dependency checking for ``LAPACK`` and ``BLAS``. 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 builds Ceres + without ``CXSparse``. 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. + +#. ``-DLINE_SEARCH_MINIMIZER=OFF``: The line search based minimizer is + mostly suitable for large scale optimization problems, or when sparse + linear algebra libraries are not available. You can further save on + some compile time and binary size by using this flag. + +#. ``-DOPENMP=OFF``: On certain platforms like Android, + multi-threading with ``OpenMP`` is not supported. Use this flag to + disable multithreading. + +#. ``-DBUILD_DOCUMENTATION=ON``: Use this flag to enable building the + documentation. In addition, ``make ceres_docs`` can be used to + build only the documentation. + +.. _section-using-ceres: + +Using Ceres with CMake +====================== + +Once the library is installed with ``make install``, it is possible to +use CMake with `FIND_PACKAGE() +<http://www.cmake.org/cmake/help/v2.8.10/cmake.html#command:find_package>`_ +in order to compile **user code** against Ceres. For example, for +`examples/helloworld.cc +<https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/helloworld.cc>`_ +the following CMakeList.txt can be used: + +.. code-block:: cmake + + CMAKE_MINIMUM_REQUIRED(VERSION 2.8) + + PROJECT(helloworld) + + FIND_PACKAGE(Ceres REQUIRED) + INCLUDE_DIRECTORIES(${CERES_INCLUDES}) + + # helloworld + ADD_EXECUTABLE(helloworld helloworld.cc) + TARGET_LINK_LIBRARIES(helloworld ${CERES_LIBRARIES}) + +Specify Ceres version +--------------------- + +Additionally, when CMake has found Ceres it can check the package +version, if it has been specified in the `FIND_PACKAGE() +<http://www.cmake.org/cmake/help/v2.8.10/cmake.html#command:find_package>`_ +call. For example: + +.. code-block:: cmake + + FIND_PACKAGE(Ceres 1.2.3 REQUIRED) + +The version is an optional argument. + +Local installations +------------------- + +If Ceres was installed in a non-standard path by specifying +-DCMAKE_INSTALL_PREFIX="/some/where/local", then the user should add +the **PATHS** option to the ``FIND_PACKAGE()`` command. e.g., + +.. code-block:: cmake + + FIND_PACKAGE(Ceres REQUIRED PATHS "/some/where/local/") + +Note that this can be used to have multiple versions of Ceres installed. + +Compiling against static or shared library +------------------------------------------ + +.. code-block:: cmake + + TARGET_LINK_LIBRARIES(helloworld ${CERES_LIBRARIES}) + +will result in a statically linked binary. Changing this line to + +.. code-block:: cmake + + TARGET_LINK_LIBRARIES(helloworld ${CERES_LIBRARIES_SHARED}) + +will result in a dynamically linked binary. diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100644 index 0000000..f5ffb6d --- /dev/null +++ b/docs/source/conf.py @@ -0,0 +1,241 @@ +# -*- 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.7' +# The full version, including alpha/beta/rc tags. +release = '1.7.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 = 'armstrong' + +# 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 = "Ceres Solver" + +# 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 = False + +# 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'Ceres Solver', + 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'Ceres Solver', + [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'Ceres Solver', + 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..20fe34d --- /dev/null +++ b/docs/source/contributing.rst @@ -0,0 +1,140 @@ +.. _chapter-contributing: + +============= +Contributions +============= + + +We welcome contributions to Ceres, whether they are new features, bug +fixes or tests. The Ceres `mailing +<http://groups.google.com/group/ceres-solver>`_ list 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 mailing list 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 +<http://google-styleguide.googlecode.com/svn/trunk/cppguide.xml>`_ and +use `git <http://git-scm.com/>`_ for version control. We use the +`Gerrit <https://ceres-solver-review.googlesource.com/>`_ to collaborate and +review changes to Ceres. Gerrit enables pre-commit reviews so that +Ceres can maintain a linear history with clean, reviewed commits, and +no merges. + +We now describe how to set up your development environment and submit +a change list for review via Gerrit. + +Setting up your Development Environment +======================================= + +1. Download and configure ``git``. + + * Mac ``brew install git``. + * Linux ``sudo apt-get install git``. + * Windows. Download `msysgit + <https://code.google.com/p/msysgit/>`_, which includes a minimal + `Cygwin <http://www.cygwin.com/>`_ install. + +2. Sign up for `Gerrit + <https://ceres-solver-review.googlesource.com/>`_. You will also + need to sign the Contributor License Agreement (CLA) with Google, + which gives Google a royalty-free unlimited license to use your + contributions. You retain copyright. + +3. Clone the Ceres Solver ``git`` repository from Gerrit. + + .. code-block:: bash + + git clone https://ceres-solver.googlesource.com/ceres-solver + + +4. Build Ceres, following the instructions in + :ref:`chapter-building`. + + On Mac and Linux, the ``CMake`` build will download and enable + the Gerrit pre-commit hook automatically. This pre-submit hook + creates `Change-Id: ...` lines in your commits. + + If this does not work OR you are on Windows, execute the + following in the root directory of the local ``git`` repository: + + .. code-block:: bash + + curl -o .git/hooks/commit-msg https://ceres-solver-review.googlesource.com/tools/hooks/commit-msg + chmod +x .git/hooks/commit-msg + +5. Configure your Gerrit password with a ``.netrc`` (Mac and Linux) + or ``_netrc`` (Windows) which allows pushing to Gerrit without + having to enter a very long random password every time: + + * Sign into `http://ceres-solver-review.googlesource.com + <http://ceres-solver-review.googlesource.com>`_. + + * Click ``Settings -> HTTP Password -> Obtain Password``. + + * (maybe) Select an account for multi-login. This should be the + same as your Gerrit login. + + * Click ``Allow access`` when the page requests access to your + ``git`` repositories. + + * Copy the contents of the ``netrc`` into the clipboard. + + - On Mac and Linux, paste the contents into ``~/.netrc``. + + - On Windows, by default users do not have a ``%HOME%`` + setting. + + + Executing ``setx HOME %USERPROFILE%`` in a terminal will set up + the ``%HOME%`` environment variable persistently, and is used + by ``git`` to find ``%HOME%\_netrc``. + + Then, create a new text file named ``_netrc`` and put it in + e.g. ``C:\Users\username`` where ``username`` is your user + name. + + +Submitting a change to Ceres Solver +=================================== + +1. Make your changes against master or whatever branch you + like. Commit your changes as one patch. When you commit, the Gerrit + hook will add a `Change-Id:` line as the last line of the commit. + +2. Push your changes to the Ceres Gerrit instance: + + .. code-block:: bash + + git push origin HEAD:refs/for/master + + When the push succeeds, the console will display a URL showing the + address of the review. Go to the URL and add reviewers; typically + this is Sameer or Keir at this point. + +3. Wait for a review. + +4. Once review comments come in, address them. Please reply to each + comment in Gerrit, which makes the re-review process easier. After + modifying the code in your ``git`` instance, *don't make a new + commit*. Instead, update the last commit using a command like the + following: + + .. code-block:: bash + + git commit --amend -a + + This will update the last commit, so that it has both the original + patch and your updates as a single commit. You will have a chance + to edit the commit message as well. Push the new commit to Gerrit + as before. + + Gerrit will use the ``Change-Id:`` to match the previous commit + with the new one. The review interface retains your original patch, + but also shows the new patch. + + Publish your responses to the comments, and wait for a new round + of reviews. diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100644 index 0000000..f20dad4 --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,51 @@ +.. 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. + +============ +Ceres Solver +============ + +Ceres Solver is a portable C++ library for solving non-linear least +squares problems. + +* Download the latest stable `release + <https://ceres-solver.googlecode.com/files/ceres-solver-1.6.0.tar.gz>`_ + or clone the `repository + <https://ceres-solver.googlesource.com/ceres-solver>`_ + +* Read the :ref:`chapter-tutorial` + +* Browse the :ref:`chapter-modeling` API and :ref:`chapter-solving` API. + +* Join the `mailing list + <https://groups.google.com/forum/?fromgroups#!forum/ceres-solver>`_ + and ask questions. + +* File bugs, feature requests in the `issue tracker + <https://code.google.com/p/ceres-solver/issues/list>`_. + +* If you use Ceres Solver for a publication, you must cite it as:: + + @misc{ceres-solver, + author = "Sameer Agarwal and Keir Mierle and Others", + title = "Ceres Solver", + howpublished = "\url{https://code.google.com/p/ceres-solver/}", + } + +.. toctree:: + :maxdepth: 1 + :hidden: + + introduction + building + tutorial + modeling + solving + reading + contributing + acknowledgements + version_history + bibliography + license diff --git a/docs/source/introduction.rst b/docs/source/introduction.rst new file mode 100644 index 0000000..19a6f2e --- /dev/null +++ b/docs/source/introduction.rst @@ -0,0 +1,81 @@ +.. _chapter-introduction: + +============ +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 accurately and efficiently. + +**Features** + +#. A friendly :ref:`chapter-modeling` API. + +#. 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 preconditioners for general sparse + and bundle adjustment problems. + +#. Portable: Runs on Linux, Windows, Mac OS X and Android. + + +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. + +`Blender <http://www.blender.org>`_ uses Ceres for `motion tracking +<http://mango.blender.org/development/planar-tracking-preview/>`_ and +`bundle adjustment +<http://wiki.blender.org/index.php/Dev:Ref/Release_Notes/2.67/Motion_Tracker>`_. + + +.. 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 the method of + Least Squares [Stigler]_, there is no debate that it was + `Carl Friedrich Gauss + <http://en.wikipedia.org/wiki/Carl_Friedrich_Gauss>`_ who + brought it to the attention of the world. Using just 22 + observations of the newly discovered asteroid `Ceres + <http://en.wikipedia.org/wiki/Ceres_(dwarf_planet)>`_, Gauss + used the method of least squares to correctly predict when + and where the asteroid will emerge from behind the Sun + [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/least_squares_fit.png b/docs/source/least_squares_fit.png Binary files differnew file mode 100644 index 0000000..7dad673 --- /dev/null +++ b/docs/source/least_squares_fit.png diff --git a/docs/source/license.rst b/docs/source/license.rst new file mode 100644 index 0000000..58d70df --- /dev/null +++ b/docs/source/license.rst @@ -0,0 +1,30 @@ +======= +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. diff --git a/docs/source/loss.png b/docs/source/loss.png Binary files differnew file mode 100644 index 0000000..9f98d00 --- /dev/null +++ b/docs/source/loss.png diff --git a/docs/source/modeling.rst b/docs/source/modeling.rst new file mode 100644 index 0000000..8e6de12 --- /dev/null +++ b/docs/source/modeling.rst @@ -0,0 +1,1530 @@ +.. default-domain:: cpp + +.. cpp:namespace:: ceres + +.. _`chapter-modeling`: + +======== +Modeling +======== + +Recall that 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 expression +: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]`. 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 be a single parameter. :math:`\rho_i` is a +:class:`LossFunction`. A :class:`LossFunction` is a scalar function +that is used to reduce the influence of outliers on the solution of +non-linear least squares problems. + +In this chapter we will describe the various classes that are part of +Ceres Solver's modeling API, and how they can be used to construct an +optimization problem. Once a problem has been constructed, various +methods for solving them will be discussed in +:ref:`chapter-solving`. It is by design that the modeling and the +solving APIs are orthogonal to each other. This enables +switching/tweaking of various solver parameters without having to +touch the problem once it has been successfully modeled. + +:class:`CostFunction` +--------------------- + +The single biggest task when modeling a problem is specifying the +residuals and their derivatives. This is done using +:class:`CostFunction` objects. + +.. 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) + + Compute the residual vector and the Jacobian matrices. + + ``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_` and the + 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 a parameter block is marked constant. + + **NOTE** The return value indicates whether the computation of the + residuals and/or jacobians was successful or not. + + This can be used to communicate numerical failures in Jacobian + computations for instance. + + A more interesting and common use is to impose constraints on the + parameters. If the initial values of the parameter blocks satisfy + the constraints, then returning false whenever the constraints are + not satisfied will prevent the solver from moving into the + infeasible region. This is not a very sophisticated mechanism for + enforcing constraints, but is often good enough for things like + non-negativity constraints. + + Note that it is important that the initial values of the parameter + block must be feasible, otherwise the solver will declare a + numerical problem at iteration 0. + + +: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), + :class:`SizeCostFunction` can be used where these values can be + specified as template parameters and the user only needs to + implement :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 + + Defining a :class:`CostFunction` or a :class:`SizedCostFunction` + can be a tedious and error prone especially when computing + derivatives. To this end Ceres provides `automatic differentiation + <http://en.wikipedia.org/wiki/Automatic_differentiation>`_. + + .. code-block:: c++ + + template <typename CostFunctor, + int M, // Number of residuals, or ceres::DYNAMIC. + int N0, // Number of parameters in block 0. + int N1 = 0, // Number of parameters in block 1. + int N2 = 0, // Number of parameters in block 2. + int N3 = 0, // Number of parameters in block 3. + int N4 = 0, // Number of parameters in block 4. + int N5 = 0, // Number of parameters in block 5. + int N6 = 0, // Number of parameters in block 6. + int N7 = 0, // Number of parameters in block 7. + int N8 = 0, // Number of parameters in block 8. + int N9 = 0> // Number of parameters in block 9. + class AutoDiffCostFunction : public + SizedCostFunction<M, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9> { + }; + + 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. + Please see :class:`CostFunction` for details on how the return + value may be used to impose simple constraints on the parameter + block. + + 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 10 + 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 + of type ``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:`DynamicAutoDiffCostFunction` +------------------------------------ + +.. class:: DynamicAutoDiffCostFunction + + :class:`AutoDiffCostFunction` requires that the number of parameter + blocks and their sizes be known at compile time, e.g., Bezier curve + fitting, Neural Network training etc. It also has an upper limit of + 10 parameter blocks. In a number of applications, this is not + enough. + + .. code-block:: c++ + + template <typename CostFunctor, int Stride = 4> + class DynamicAutoDiffCostFunction : public CostFunction { + }; + + In such cases :class:`DynamicAutoDiffCostFunction` can be + used. Like :class:`AutoDiffCostFunction` the user must define a + templated functor, but the signature of the functor differs + slightly. The expected interface for the cost functors is: + + .. code-block:: c++ + + struct MyCostFunctor { + template<typename T> + bool operator()(T const* const* parameters, T* residuals) const { + } + } + + Since the sizing of the parameters is done at runtime, you must + also specify the sizes after creating the dynamic autodiff cost + function. For example: + + .. code-block:: c++ + + DynamicAutoDiffCostFunction<MyCostFunctor, 4> cost_function( + new MyCostFunctor()); + cost_function.AddParameterBlock(5); + cost_function.AddParameterBlock(10); + cost_function.SetNumResiduals(21); + + Under the hood, the implementation evaluates the cost function + multiple times, computing a small set of the derivatives (four by + default, controlled by the ``Stride`` template parameter) with each + pass. There is a performance tradeoff with the size of the passes; + Smaller sizes are more cache efficient but result in larger number + of passes, and larger stride lengths can destroy cache-locality + while reducing the number of passes over the cost function. The + optimal value depends on the number and sizes of the various + parameter blocks. + + As a rule of thumb, try using :class:`AutoDiffCostFunction` before + you use :class:`DynamicAutoDiffCostFunction`. + +:class:`NumericDiffCostFunction` +-------------------------------- + +.. class:: NumericDiffCostFunction + + In some cases, its not possible to define a templated cost functor, + for example when the evaluation of the residual involves a call to a + library function that you do not have control over. In such a + situation, `numerical differentiation + <http://en.wikipedia.org/wiki/Numerical_differentiation>`_ can be + used. + + .. 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> { + }; + + To get a 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. Please see :class:`CostFunction` for details on + how the return value may be used to impose simple constraints on + the parameter block. 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 measurement + 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 necessary by repeatedly calling the + :func:`CostFunction::Evaluate` with small changes to the + appropriate parameters, and computing the slope. For performance, + the numeric differentiation wrapper class is templated on the + concrete cost function, even though it could be implemented only in + terms of the :class:`CostFunction` interface. + + The numerically differentiated version of a cost function for a + cost function can be constructed as follows: + + .. 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:`NumericDiffFunctor` +--------------------------- + +.. class:: NumericDiffFunctor + + Sometimes parts of a cost function can be differentiated + automatically or analytically but others require numeric + differentiation. :class:`NumericDiffFunctor` is 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:`CostFunctionToFunctor` +------------------------------ + +.. class:: CostFunctionToFunctor + + Just like :class:`NumericDiffFunctor` allows numeric + differentiation to be mixed with automatic differentiation, + :class:`CostFunctionToFunctor` provides an even more general + mechanism. :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:`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:`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:`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 predefined 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 + + Given two loss functions ``f`` and ``g``, implements the loss + function ``h(s) = f(g(s))``. + + .. code-block:: c++ + + class ComposedLoss : public LossFunction { + public: + explicit ComposedLoss(const LossFunction* f, + Ownership ownership_f, + const LossFunction* g, + Ownership ownership_g); + }; + +.. class:: ScaledLoss + + Sometimes you want to simply scale the output value of the + robustifier. For example, you might want to weight different error + terms differently (e.g., weight pixel reprojection errors + differently from terrain errors). + + Given a loss function :math:`\rho(s)` and a scalar :math:`a`, :class:`ScaledLoss` + implements the function :math:`a \rho(s)`. + + Since we treat the a ``NULL`` Loss function as the Identity loss + function, :math:`rho` = ``NULL``: is a valid input and will result + in the input being scaled by :math:`a`. This provides a simple way + of implementing a scaled ResidualBlock. + +.. class:: LossFunctionWrapper + + Sometimes after the optimization problem has been constructed, we + wish to mutate the scale of the loss function. For example, when + performing estimation from data which has substantial outliers, + convergence can be improved by starting out with a large scale, + optimizing the problem and then reducing the scale. This can have + better convergence behavior than just using a loss function with a + small scale. + + This templated class allows the user to implement a loss function + whose scale can be mutated after an optimization problem has been + constructed. e.g, + + .. code-block:: c++ + + Problem problem; + + // Add parameter blocks + + CostFunction* cost_function = + new AutoDiffCostFunction < UW_Camera_Mapper, 2, 9, 3>( + new UW_Camera_Mapper(feature_x, feature_y)); + + LossFunctionWrapper* loss_function(new HuberLoss(1.0), TAKE_OWNERSHIP); + problem.AddResidualBlock(cost_function, loss_function, parameters); + + Solver::Options options; + Solver::Summary summary; + Solve(options, &problem, &summary); + + loss_function->Reset(new HuberLoss(1.0), TAKE_OWNERSHIP); + Solve(options, &problem, &summary); + + +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 robustified 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:`AutoDiffLocalParameterization` +-------------------------------------- + +.. class:: AutoDiffLocalParameterization + + :class:`AutoDiffLocalParameterization` does for + :class:`LocalParameterization` what :class:`AutoDiffCostFunction` + does for :class:`CostFunction`. It allows the user to define a + templated functor that implements the + :func:`LocalParameterization::Plus` operation and it uses automatic + differentiation to implement the computation of the Jacobian. + + To get an auto differentiated local parameterization, you must + define a class with a templated operator() (a functor) that computes + + .. math:: x' = \boxplus(x, \Delta x), + + For example, Quaternions have a three dimensional local + parameterization. It's plus operation can be implemented as (taken + from `internal/ceres/auto_diff_local_parameterization_test.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/include/ceres/local_parameterization.h>`_ + ) + + .. code-block:: c++ + + struct QuaternionPlus { + template<typename T> + bool operator()(const T* x, const T* delta, T* x_plus_delta) const { + const T squared_norm_delta = + delta[0] * delta[0] + delta[1] * delta[1] + delta[2] * delta[2]; + + T q_delta[4]; + if (squared_norm_delta > T(0.0)) { + T norm_delta = sqrt(squared_norm_delta); + const T sin_delta_by_delta = sin(norm_delta) / norm_delta; + q_delta[0] = cos(norm_delta); + q_delta[1] = sin_delta_by_delta * delta[0]; + q_delta[2] = sin_delta_by_delta * delta[1]; + q_delta[3] = sin_delta_by_delta * delta[2]; + } else { + // We do not just use q_delta = [1,0,0,0] here because that is a + // constant and when used for automatic differentiation will + // lead to a zero derivative. Instead we take a first order + // approximation and evaluate it at zero. + q_delta[0] = T(1.0); + q_delta[1] = delta[0]; + q_delta[2] = delta[1]; + q_delta[3] = delta[2]; + } + + Quaternionproduct(q_delta, x, x_plus_delta); + return true; + } + }; + + Then given this struct, the auto differentiated local + parameterization can now be constructed as + + .. code-block:: c++ + + LocalParameterization* local_parameterization = + new AutoDiffLocalParameterization<QuaternionPlus, 4, 3>; + | | + Global Size ---------------+ | + Local Size -------------------+ + + **WARNING:** Since the functor will get instantiated with different + types for ``T``, you must to convert from other numeric types to + ``T`` before mixing computations with other variables of type + ``T``. In the example above, this is seen where instead of using + ``k_`` directly, ``k_`` is wrapped with ``T(k_)``. + + +: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. + + :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 behavior. + + 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. + + 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. + +.. function:: ResidualBlockId Problem::AddResidualBlock(CostFunction* cost_function, LossFunction* loss_function, const vector<double*> parameter_blocks) + + Add a residual block to the overall cost function. 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 AddParameterBlock. This causes additional correctness + checking; however, AddResidualBlock implicitly adds the parameter + blocks if they are not present, so calling AddParameterBlock + explicitly is not required. + + The Problem object by default takes ownership of the + cost_function and loss_function pointers. These objects remain + live for the life of the 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 Options struct. + + Note: 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. + + Example usage: + + .. code-block:: c++ + + double x1[] = {1.0, 2.0, 3.0}; + double x2[] = {1.0, 2.0, 5.0, 6.0}; + double x3[] = {3.0, 6.0, 2.0, 5.0, 1.0}; + + Problem problem; + + problem.AddResidualBlock(new MyUnaryCostFunction(...), NULL, x1); + problem.AddResidualBlock(new MyBinaryCostFunction(...), NULL, x2, x1); + + +.. function:: void Problem::AddParameterBlock(double* values, int size, LocalParameterization* local_parameterization) + + Add a parameter block with appropriate size to the problem. + Repeated calls with the same arguments are ignored. Repeated calls + with the same double pointer but a different size results in + undefined behavior. + +.. function:: void Problem::AddParameterBlock(double* values, int size) + + Add a parameter block with appropriate size and parameterization to + the problem. Repeated calls with the same arguments are + ignored. Repeated calls with the same double pointer but a + different size results in undefined behavior. + +.. function:: void Problem::RemoveResidualBlock(ResidualBlockId residual_block) + + Remove a residual block from the problem. Any parameters that the residual + block depends on are not removed. The cost and loss functions for the + residual block will not get deleted immediately; won't happen until the + problem itself is deleted. + + **WARNING:** Removing a residual or parameter block will destroy + the implicit ordering, rendering the jacobian or residuals returned + from the solver uninterpretable. If you depend on the evaluated + jacobian, do not use remove! This may change in a future release. + Hold the indicated parameter block constant during optimization. + +.. function:: void Problem::RemoveParameterBlock(double* values) + + Remove a parameter block from the problem. The parameterization of + the parameter block, if it exists, will persist until the deletion + of the problem (similar to cost/loss functions in residual block + removal). Any residual blocks that depend on the parameter are also + removed, as described above in RemoveResidualBlock(). If + Problem::Options::enable_fast_parameter_block_removal is true, then + the removal is fast (almost constant time). Otherwise, removing a + parameter block will incur a scan of the entire Problem object. + + **WARNING:** Removing a residual or parameter block will destroy + the implicit ordering, rendering the jacobian or residuals returned + from the solver uninterpretable. If you depend on the evaluated + jacobian, do not use remove! This may change in a future release. + +.. function:: void Problem::SetParameterBlockConstant(double* values) + + Hold the indicated parameter block constant during optimization. + +.. function:: void Problem::SetParameterBlockVariable(double* values) + + Allow the indicated parameter to vary during optimization. + +.. function:: void Problem::SetParameterization(double* values, LocalParameterization* local_parameterization) + + Set the local parameterization for one of the parameter blocks. + The local_parameterization is owned by the Problem by default. It + is acceptable to set the same parameterization for multiple + parameters; the destructor is careful to delete local + parameterizations only once. The local parameterization can only be + set once per parameter, and cannot be changed once set. + +.. function:: int Problem::NumParameterBlocks() const + + Number of parameter blocks in the problem. Always equals + parameter_blocks().size() and parameter_block_sizes().size(). + +.. function:: int Problem::NumParameters() const + + The size of the parameter vector obtained by summing over the sizes + of all the parameter blocks. + +.. function:: int Problem::NumResidualBlocks() const + + Number of residual blocks in the problem. Always equals + residual_blocks().size(). + +.. function:: int Problem::NumResiduals() const + + The size of the residual vector obtained by summing over the sizes + of all of the residual blocks. + +.. function int Problem::ParameterBlockSize(const double* values) const; + + The size of the parameter block. + +.. function int Problem::ParameterBlockLocalSize(const double* values) const; + + The size of local parameterization for the parameter block. If + there is no local parameterization associated with this parameter + block, then ``ParameterBlockLocalSize`` = ``ParameterBlockSize``. + + +.. function void Problem::GetParameterBlocks(vector<double*>* parameter_blocks) const; + + Fills the passed ``parameter_blocks`` vector with pointers to the + parameter blocks currently in the problem. After this call, + ``parameter_block.size() == NumParameterBlocks``. + +.. function:: bool Problem::Evaluate(const Problem::EvaluateOptions& options, double* cost, vector<double>* residuals, vector<double>* gradient, CRSMatrix* jacobian) + + Evaluate a :class:`Problem`. Any of the output pointers can be + `NULL`. Which residual blocks and parameter blocks are used is + controlled by the :class:`Problem::EvaluateOptions` struct below. + + .. code-block:: c++ + + Problem problem; + double x = 1; + problem.Add(new MyCostFunction, NULL, &x); + + double cost = 0.0; + problem.Evaluate(Problem::EvaluateOptions(), &cost, NULL, NULL, NULL); + + The cost is evaluated at `x = 1`. If you wish to evaluate the + problem at `x = 2`, then + + .. code-block:: c++ + + x = 2; + problem.Evaluate(Problem::EvaluateOptions(), &cost, NULL, NULL, NULL); + + is the way to do so. + + **NOTE** If no local parameterizations are used, then the size of + the gradient vector is the sum of the sizes of all the parameter + blocks. If a parameter block has a local parameterization, then + it contributes "LocalSize" entries to the gradient vector. + +.. class:: Problem::EvaluateOptions + + Options struct that is used to control :func:`Problem::Evaluate`. + +.. member:: vector<double*> Problem::EvaluateOptions::parameter_blocks + + The set of parameter blocks for which evaluation should be + performed. This vector determines the order in which parameter + blocks occur in the gradient vector and in the columns of the + jacobian matrix. If parameter_blocks is empty, then it is assumed + to be equal to a vector containing ALL the parameter + blocks. Generally speaking the ordering of the parameter blocks in + this case depends on the order in which they were added to the + problem and whether or not the user removed any parameter blocks. + + **NOTE** This vector should contain the same pointers as the ones + used to add parameter blocks to the Problem. These parameter block + should NOT point to new memory locations. Bad things will happen if + you do. + +.. member:: vector<ResidualBlockId> Problem::EvaluateOptions::residual_blocks + + The set of residual blocks for which evaluation should be + performed. This vector determines the order in which the residuals + occur, and how the rows of the jacobian are ordered. If + residual_blocks is empty, then it is assumed to be equal to the + vector containing all the parameter blocks. + +``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:: void AngleAxisToQuaternion<T>(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:: void QuaternionToAngleAxis<T>(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:: void RotationMatrixToAngleAxis<T, row_stride, col_stride>(const MatrixAdapter<const T, row_stride, col_stride>& R, T * angle_axis) +.. function:: void AngleAxisToRotationMatrix<T, row_stride, col_stride>(T const * angle_axis, const MatrixAdapter<T, row_stride, col_stride>& R) +.. function:: void RotationMatrixToAngleAxis<T>(T const * R, T * angle_axis) +.. function:: void AngleAxisToRotationMatrix<T>(T const * angle_axis, T * R) + + Conversions between 3x3 rotation matrix with given column and row strides and + axis-angle rotation representations. The functions that take a pointer to T instead + of a MatrixAdapter assume a column major representation with unit row stride and a column stride of 3. + +.. function:: void EulerAnglesToRotationMatrix<T, row_stride, col_stride>(const T* euler, const MatrixAdapter<T, row_stride, col_stride>& R) +.. function:: void EulerAnglesToRotationMatrix<T>(const T* euler, int row_stride, T* R) + + Conversions between 3x3 rotation matrix with given column and row strides 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. + + The function that takes a pointer to T as the rotation matrix assumes a row + major representation with unit column stride and a row stride of 3. + The additional parameter row_stride is required to be 3. + +.. function:: void QuaternionToScaledRotation<T, row_stride, col_stride>(const T q[4], const MatrixAdapter<T, row_stride, col_stride>& R) +.. function:: void QuaternionToScaledRotation<T>(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`. + + In the function that accepts a pointer to T instead of a MatrixAdapter, + the rotation matrix ``R`` is a row-major matrix with unit column stride + and a row stride of 3. + + 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:: void QuaternionToRotation<T>(const T q[4], const MatrixAdapter<T, row_stride, col_stride>& R) +.. function:: void QuaternionToRotation<T>(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:: void UnitQuaternionRotatePoint<T>(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:: void QuaternionRotatePoint<T>(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:: void QuaternionProduct<T>(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:: void CrossProduct<T>(const T x[3], const T y[3], T x_cross_y[3]) + + .. math:: \text{x_cross_y} = x \times y + +.. function:: void AngleAxisRotatePoint<T>(const T angle_axis[3], const T pt[3], T result[3]) + + .. math:: y = R(\text{angle_axis}) x diff --git a/docs/source/non_robust_least_squares_fit.png b/docs/source/non_robust_least_squares_fit.png Binary files differnew file mode 100644 index 0000000..643d162 --- /dev/null +++ b/docs/source/non_robust_least_squares_fit.png diff --git a/docs/source/reading.rst b/docs/source/reading.rst new file mode 100644 index 0000000..4e27567 --- /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 favorite text on introductory numerical linear algebra. [Triggs]_ +provides a thorough coverage of the bundle adjustment problem. diff --git a/docs/source/robust_least_squares_fit.png b/docs/source/robust_least_squares_fit.png Binary files differnew file mode 100644 index 0000000..89003c9 --- /dev/null +++ b/docs/source/robust_least_squares_fit.png diff --git a/docs/source/solving.rst b/docs/source/solving.rst new file mode 100644 index 0000000..d8b9f4a --- /dev/null +++ b/docs/source/solving.rst @@ -0,0 +1,2229 @@ + +.. default-domain:: cpp + +.. cpp:namespace:: ceres + +.. _chapter-solving: + +======= +Solving +======= + + +Introduction +============ + +Effective use of Ceres requires some familiarity with the basic +components of a nonlinear least squares solver, so before we describe +how to configure and use the solver, we will take a brief look at how +some of the core optimization algorithms in Ceres work. + +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`. Depending on how the size of +the step :math:`\Delta x` is controlled, non-linear optimization +algorithms can be divided into two major categories [NocedalWright]_. + +1. **Trust Region** The trust region approach approximates the + objective function using using a model function (often a quadratic) + over a subset of the search space known as the trust region. If the + model function succeeds in minimizing the true objective function + the trust region is expanded; conversely, otherwise it is + contracted and the model optimization problem is solved again. + +2. **Line Search** The line search approach first finds a descent + direction along which the objective function will be reduced and + then computes a step size that decides how far should move along + that direction. The descent direction can be computed by various + methods, such as gradient descent, Newton's method and Quasi-Newton + method. The step size can be determined either exactly or + inexactly. + +Trust region methods are in some sense dual to line search methods: +trust region methods first choose a step size (the size of the trust +region) and then a step direction while line search methods first +choose a step direction and then a step size. Ceres implements +multiple algorithms in both categories. + +.. _section-trust-region-methods: + +Trust Region Methods +==================== + +The basic trust region algorithm looks something like this. + + 1. Given an initial point :math:`x` and a trust region radius :math:`\mu`. + 2. :math:`\arg \min_{\Delta x} \frac{1}{2}\|J(x)\Delta + x + F(x)\|^2` s.t. :math:`\|D(x)\Delta x\|^2 \le \mu` + 3. :math:`\rho = \frac{\displaystyle \|F(x + \Delta x)\|^2 - + \|F(x)\|^2}{\displaystyle \|J(x)\Delta x + F(x)\|^2 - + \|F(x)\|^2}` + 4. if :math:`\rho > \epsilon` then :math:`x = x + \Delta x`. + 5. if :math:`\rho > \eta_1` then :math:`\rho = 2 \rho` + 6. else if :math:`\rho < \eta_2` then :math:`\rho = 0.5 * \rho` + 7. Goto 2. + +Here, :math:`\mu` is the trust region radius, :math:`D(x)` is some +matrix used to define a metric on the domain of :math:`F(x)` and +: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 + 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 +[ByrdSchnabel]_. + +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 principled 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 steps is available for all trust +region strategies. + + +.. _section-line-search-methods: + +Line Search Methods +=================== + +**The implementation of line search algorithms in Ceres Solver is +fairly new and not very well tested, so for now this part of the +solver should be considered beta quality. We welcome reports of your +experiences both good and bad on the mailinglist.** + +Line search algorithms + + 1. Given an initial point :math:`x` + 2. :math:`\Delta x = -H^{-1}(x) g(x)` + 3. :math:`\arg \min_\mu \frac{1}{2} \| F(x + \mu \Delta x) \|^2` + 4. :math:`x = x + \mu \Delta x` + 5. Goto 2. + +Here :math:`H(x)` is some approximation to the Hessian of the +objective function, and :math:`g(x)` is the gradient at +:math:`x`. Depending on the choice of :math:`H(x)` we get a variety of +different search directions -`\Delta x`. + +Step 4, which is a one dimensional optimization or `Line Search` along +:math:`\Delta x` is what gives this class of methods its name. + +Different line search algorithms differ in their choice of the search +direction :math:`\Delta x` and the method used for one dimensional +optimization along :math:`\Delta x`. The choice of :math:`H(x)` is the +primary source of computational complexity in these +methods. Currently, Ceres Solver supports three choices of search +directions, all aimed at large scale problems. + +1. ``STEEPEST_DESCENT`` This corresponds to choosing :math:`H(x)` to + be the identity matrix. This is not a good search direction for + anything but the simplest of the problems. It is only included here + for completeness. + +2. ``NONLINEAR_CONJUGATE_GRADIENT`` A generalization of the Conjugate + Gradient method to non-linear functions. The generalization can be + performed in a number of different ways, resulting in a variety of + search directions. Ceres Solver currently supports + ``FLETCHER_REEVES``, ``POLAK_RIBIRERE`` and ``HESTENES_STIEFEL`` + directions. + +3. ``BFGS`` A generalization of the Secant method to multiple + dimensions in which a full, dense approximation to the inverse + Hessian is maintained and used to compute a quasi-Newton step + [NocedalWright]_. BFGS is currently the best known general + quasi-Newton algorithm. + +4. ``LBFGS`` A limited memory approximation to the full ``BFGS`` + method in which the last `M` iterations are used to approximate the + inverse Hessian used to compute a quasi-Newton step [Nocedal]_, + [ByrdNocedal]_. + +Currently Ceres Solver supports both a backtracking and interpolation +based Armijo line search algorithm, and a sectioning / zoom +interpolation (strong) Wolfe condition line search algorithm. +However, note that in order for the assumptions underlying the +``BFGS`` and ``LBFGS`` methods to be guaranteed to be satisfied the +Wolfe line search algorithm should be used. + +.. _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:: MinimizerType Solver::Options::minimizer_type + + Default: ``TRUST_REGION`` + + Choose between ``LINE_SEARCH`` and ``TRUST_REGION`` algorithms. See + :ref:`section-trust-region-methods` and + :ref:`section-line-search-methods` for more details. + +.. member:: LineSearchDirectionType Solver::Options::line_search_direction_type + + Default: ``LBFGS`` + + Choices are ``STEEPEST_DESCENT``, ``NONLINEAR_CONJUGATE_GRADIENT``, + ``BFGS`` and ``LBFGS``. + +.. member:: LineSearchType Solver::Options::line_search_type + + Default: ``WOLFE`` + + Choices are ``ARMIJO`` and ``WOLFE`` (strong Wolfe conditions). + Note that in order for the assumptions underlying the ``BFGS`` and + ``LBFGS`` line search direction algorithms to be guaranteed to be + satisifed, the ``WOLFE`` line search should be used. + +.. member:: NonlinearConjugateGradientType Solver::Options::nonlinear_conjugate_gradient_type + + Default: ``FLETCHER_REEVES`` + + Choices are ``FLETCHER_REEVES``, ``POLAK_RIBIRERE`` and + ``HESTENES_STIEFEL``. + +.. member:: int Solver::Options::max_lbfs_rank + + Default: 20 + + The L-BFGS hessian approximation is a low rank approximation to the + inverse of the Hessian matrix. The rank of the approximation + determines (linearly) the space and time complexity of using the + approximation. Higher the rank, the better is the quality of the + approximation. The increase in quality is however is bounded for a + number of reasons. + + 1. The method only uses secant information and not actual + derivatives. + + 2. The Hessian approximation is constrained to be positive + definite. + + So increasing this rank to a large number will cost time and space + complexity without the corresponding increase in solution + quality. There are no hard and fast rules for choosing the maximum + rank. The best choice usually requires some problem specific + experimentation. + +.. member:: bool Solver::Options::use_approximate_eigenvalue_bfgs_scaling + + Default: ``false`` + + As part of the ``BFGS`` update step / ``LBFGS`` right-multiply + step, the initial inverse Hessian approximation is taken to be the + Identity. However, [Oren]_ showed that using instead :math:`I * + \gamma`, where :math:`\gamma` is a scalar chosen to approximate an + eigenvalue of the true inverse Hessian can result in improved + convergence in a wide variety of cases. Setting + ``use_approximate_eigenvalue_bfgs_scaling`` to true enables this + scaling in ``BFGS`` (before first iteration) and ``LBFGS`` (at each + iteration). + + Precisely, approximate eigenvalue scaling equates to + + .. math:: \gamma = \frac{y_k' s_k}{y_k' y_k} + + With: + + .. math:: y_k = \nabla f_{k+1} - \nabla f_k + .. math:: s_k = x_{k+1} - x_k + + Where :math:`f()` is the line search objective and :math:`x` the + vector of parameter values [NocedalWright]_. + + It is important to note that approximate eigenvalue scaling does + **not** *always* improve convergence, and that it can in fact + *significantly* degrade performance for certain classes of problem, + which is why it is disabled by default. In particular it can + degrade performance when the sensitivity of the problem to different + parameters varies significantly, as in this case a single scalar + factor fails to capture this variation and detrimentally downscales + parts of the Jacobian approximation which correspond to + low-sensitivity parameters. It can also reduce the robustness of the + solution to errors in the Jacobians. + +.. member:: LineSearchIterpolationType Solver::Options::line_search_interpolation_type + + Default: ``CUBIC`` + + Degree of the polynomial used to approximate the objective + function. Valid values are ``BISECTION``, ``QUADRATIC`` and + ``CUBIC``. + +.. member:: double Solver::Options::min_line_search_step_size + + The line search terminates if: + + .. math:: \|\Delta x_k\|_\infty < \text{min_line_search_step_size} + + where :math:`\|\cdot\|_\infty` refers to the max norm, and + :math:`\Delta x_k` is the step change in the parameter values at + the :math:`k`-th iteration. + +.. member:: double Solver::Options::line_search_sufficient_function_decrease + + Default: ``1e-4`` + + Solving the line search problem exactly is computationally + prohibitive. Fortunately, line search based optimization algorithms + can still guarantee convergence if instead of an exact solution, + the line search algorithm returns a solution which decreases the + value of the objective function sufficiently. More precisely, we + are looking for a step size s.t. + + .. math:: f(\text{step_size}) \le f(0) + \text{sufficient_decrease} * [f'(0) * \text{step_size}] + + This condition is known as the Armijo condition. + +.. member:: double Solver::Options::max_line_search_step_contraction + + Default: ``1e-3`` + + In each iteration of the line search, + + .. math:: \text{new_step_size} >= \text{max_line_search_step_contraction} * \text{step_size} + + Note that by definition, for contraction: + + .. math:: 0 < \text{max_step_contraction} < \text{min_step_contraction} < 1 + +.. member:: double Solver::Options::min_line_search_step_contraction + + Default: ``0.6`` + + In each iteration of the line search, + + .. math:: \text{new_step_size} <= \text{min_line_search_step_contraction} * \text{step_size} + + Note that by definition, for contraction: + + .. math:: 0 < \text{max_step_contraction} < \text{min_step_contraction} < 1 + +.. member:: int Solver::Options::max_num_line_search_step_size_iterations + + Default: ``20`` + + Maximum number of trial step size iterations during each line + search, if a step size satisfying the search conditions cannot be + found within this number of trials, the line search will stop. + + As this is an 'artificial' constraint (one imposed by the user, not + the underlying math), if ``WOLFE`` line search is being used, *and* + points satisfying the Armijo sufficient (function) decrease + condition have been found during the current search (in :math:`<=` + ``max_num_line_search_step_size_iterations``). Then, the step size + with the lowest function value which satisfies the Armijo condition + will be returned as the new valid step, even though it does *not* + satisfy the strong Wolfe conditions. This behaviour protects + against early termination of the optimizer at a sub-optimal point. + +.. member:: int Solver::Options::max_num_line_search_direction_restarts + + Default: ``5`` + + Maximum number of restarts of the line search direction algorithm + before terminating the optimization. Restarts of the line search + direction algorithm occur when the current algorithm fails to + produce a new descent direction. This typically indicates a + numerical failure, or a breakdown in the validity of the + approximations used. + +.. member:: double Solver::Options::line_search_sufficient_curvature_decrease + + Default: ``0.9`` + + The strong Wolfe conditions consist of the Armijo sufficient + decrease condition, and an additional requirement that the + step size be chosen s.t. the *magnitude* ('strong' Wolfe + conditions) of the gradient along the search direction + decreases sufficiently. Precisely, this second condition + is that we seek a step size s.t. + + .. math:: \|f'(\text{step_size})\| <= \text{sufficient_curvature_decrease} * \|f'(0)\| + + Where :math:`f()` is the line search objective and :math:`f'()` is the derivative + of :math:`f` with respect to the step size: :math:`\frac{d f}{d~\text{step size}}`. + +.. member:: double Solver::Options::max_line_search_step_expansion + + Default: ``10.0`` + + During the bracketing phase of a Wolfe line search, the step size + is increased until either a point satisfying the Wolfe conditions + is found, or an upper bound for a bracket containinqg a point + satisfying the conditions is found. Precisely, at each iteration + of the expansion: + + .. math:: \text{new_step_size} <= \text{max_step_expansion} * \text{step_size} + + By definition for expansion + + .. math:: \text{max_step_expansion} > 1.0 + +.. 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 [ByrdSchnabel]_ . 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 + accepted. + +.. member:: double Solver::Options::min_lm_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::max_lm_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:: 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. + + See :ref:`section-ordering` for more details. + +.. member:: bool Solver::Options::use_post_ordering + + Default: ``false`` + + Sparse Cholesky factorization algorithms use a fill-reducing + ordering to permute the columns of the Jacobian matrix. There are + two ways of doing this. + + 1. Compute the Jacobian matrix in some order and then have the + factorization algorithm permute the columns of the Jacobian. + + 2. Compute the Jacobian with its columns already permuted. + + The first option incurs a significant memory penalty. The + factorization algorithm has to make a copy of the permuted Jacobian + matrix, thus Ceres pre-permutes the columns of the Jacobian matrix + and generally speaking, there is no performance penalty for doing + so. + + In some rare cases, it is worth using a more complicated reordering + algorithm which has slightly better runtime performance at the + expense of an extra copy of the Jacobian matrix. Setting + ``use_postordering`` to ``true`` enables this tradeoff. + +.. member:: int Solver::Options::min_linear_solver_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::max_linear_solver_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:: 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:: double Solver::Options::inner_itearation_tolerance + + Default: ``1e-3`` + + Generally speaking, inner iterations make significant progress in + the early stages of the solve and then their contribution drops + down sharply, at which point the time spent doing inner iterations + is not worth it. + + Once the relative decrease in the objective function due to inner + iterations drops below ``inner_iteration_tolerance``, the use of + inner iterations in subsequent trust region minimizer iterations is + disabled. + +.. 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. Not all parameter blocks need to be included in + the ordering. + + See :ref:`section-ordering` for more details. + +.. 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``. + + For ``TRUST_REGION_MINIMIZER`` the progress display looks like + + .. 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: 6.91e-06 tt: 1.91e-03 + 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: 2.81e-05 tt: 1.99e-03 + 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: 1.00e-05 tt: 2.01e-03 + + Here + + #. ``f`` is the value of the objective function. + #. ``d`` is the change in the value of the objective function if + the step computed in this iteration is accepted. + #. ``g`` is the max norm of the gradient. + #. ``h`` is the change in the parameter vector. + #. ``rho`` is the ratio of the actual change in the objective + function value to the change in the the value of the trust + region model. + #. ``mu`` is the size of the trust region radius. + #. ``li`` is the number of linear solver iterations used to compute + the trust region step. For direct/factorization based solvers it + is always 1, for iterative solvers like ``ITERATIVE_SCHUR`` it + is the number of iterations of the Conjugate Gradients + algorithm. + #. ``it`` is the time take by the current iteration. + #. ``tt`` is the the total time taken by the minimizer. + + For ``LINE_SEARCH_MINIMIZER`` the progress display looks like + + .. code-block:: bash + + 0: f: 2.317806e+05 d: 0.00e+00 g: 3.19e-01 h: 0.00e+00 s: 0.00e+00 e: 0 it: 2.98e-02 tt: 8.50e-02 + 1: f: 2.312019e+05 d: 5.79e+02 g: 3.18e-01 h: 2.41e+01 s: 1.00e+00 e: 1 it: 4.54e-02 tt: 1.31e-01 + 2: f: 2.300462e+05 d: 1.16e+03 g: 3.17e-01 h: 4.90e+01 s: 2.54e-03 e: 1 it: 4.96e-02 tt: 1.81e-01 + + Here + + #. ``f`` is the value of the objective function. + #. ``d`` is the change in the value of the objective function if + the step computed in this iteration is accepted. + #. ``g`` is the max norm of the gradient. + #. ``h`` is the change in the parameter vector. + #. ``s`` is the optimal step length computed by the line search. + #. ``it`` is the time take by the current iteration. + #. ``tt`` is the the total time taken by the minimizer. + +.. member:: vector<int> Solver::Options::trust_region_minimizer_iterations_to_dump + + Default: ``empty`` + + List of iterations at which the trust region minimizer should dump + the trust region problem. Useful for testing and benchmarking. If + ``empty``, no problems are dumped. + +.. member:: string Solver::Options::trust_region_problem_dump_directory + + Default: ``/tmp`` + + Directory to which the problems should be written to. Should be + non-empty if + :member:`Solver::Options::trust_region_minimizer_iterations_to_dump` is + non-empty and + :member:`Solver::Options::trust_region_problem_dump_format_type` is not + ``CONSOLE``. + +.. member:: DumpFormatType Solver::Options::trust_region_problem_dump_format + + Default: ``TEXTFILE`` + + The format in which trust region problems should be logged when + :member:`Solver::Options::trust_region_minimizer_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. + + * ``TEXTFILE`` Write out the linear least squares problem to the + directory pointed to by + :member:`Solver::Options::trust_region_problem_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 + ``ceres_solver_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 + + ``ParameterBlockOrdering`` is a class for storing and manipulating + an ordered collection of groups/sets with the following semantics: + + Group IDs are non-negative integer values. Elements are any type + that can serve as a key in a map or an element of a set. + + An element can only belong to one group at a time. A group may + contain an arbitrary number of elements. + + Groups are ordered by their group id. + +.. function:: bool ParameterBlockOrdering::AddElementToGroup(const double* element, const int group) + + Add an element to a group. If a group with this id does not exist, + one is created. This method can be called any number of times for + the same element. Group ids should be non-negative numbers. Return + value indicates if adding the element was a success. + +.. function:: void ParameterBlockOrdering::Clear() + + Clear the ordering. + +.. function:: bool ParameterBlockOrdering::Remove(const double* element) + + Remove the element, no matter what group it is in. If the element + is not a member of any group, calling this method will result in a + crash. Return value indicates if the element was actually removed. + +.. function:: void ParameterBlockOrdering::Reverse() + + Reverse the order of the groups in place. + +.. function:: int ParameterBlockOrdering::GroupId(const double* element) const + + Return the group id for the element. If the element is not a member + of any group, return -1. + +.. function:: bool ParameterBlockOrdering::IsMember(const double* element) const + + True if there is a group containing the parameter block. + +.. function:: int ParameterBlockOrdering::GroupSize(const int group) const + + This function always succeeds, i.e., implicitly there exists a + group for every integer. + +.. function:: int ParameterBlockOrdering::NumElements() const + + Number of elements in the ordering. + +.. function:: int ParameterBlockOrdering::NumGroups() const + + Number of groups with one or more elements. + + +:class:`IterationCallback` +-------------------------- + +.. class:: IterationSummary + + :class:`IterationSummary` describes the state of the optimizer + after each iteration of the minimization. Note that all times are + wall times. + + .. code-block:: c++ + + struct IterationSummary { + // Current iteration number. + int32 iteration; + + // Step was numerically valid, i.e., all values are finite and the + // step reduces the value of the linearized model. + // + // Note: step_is_valid is false when iteration = 0. + bool step_is_valid; + + // Step did not reduce the value of the objective function + // sufficiently, but it was accepted because of the relaxed + // acceptance criterion used by the non-monotonic trust region + // algorithm. + // + // Note: step_is_nonmonotonic is false when iteration = 0; + bool step_is_nonmonotonic; + + // Whether or not the minimizer accepted this step or not. If the + // ordinary trust region algorithm is used, this means that the + // relative reduction in the objective function value was greater + // than Solver::Options::min_relative_decrease. However, if the + // non-monotonic trust region algorithm is used + // (Solver::Options:use_nonmonotonic_steps = true), then even if the + // relative decrease is not sufficient, the algorithm may accept the + // step and the step is declared successful. + // + // Note: step_is_successful is false when iteration = 0. + bool step_is_successful; + + // Value of the objective function. + double cost; + + // Change in the value of the objective function in this + // iteration. This can be positive or negative. + double cost_change; + + // Infinity norm of the gradient vector. + double gradient_max_norm; + + // 2-norm of the size of the step computed by the optimization + // algorithm. + double step_norm; + + // For trust region algorithms, the ratio of the actual change in + // cost and the change in the cost of the linearized approximation. + double relative_decrease; + + // Size of the trust region at the end of the current iteration. For + // the Levenberg-Marquardt algorithm, the regularization parameter + // mu = 1.0 / trust_region_radius. + double trust_region_radius; + + // For the inexact step Levenberg-Marquardt algorithm, this is the + // relative accuracy with which the Newton(LM) step is solved. This + // number affects only the iterative solvers capable of solving + // linear systems inexactly. Factorization-based exact solvers + // ignore it. + double eta; + + // Step sized computed by the line search algorithm. + double step_size; + + // Number of function evaluations used by the line search algorithm. + int line_search_function_evaluations; + + // Number of iterations taken by the linear solver to solve for the + // Newton step. + int linear_solver_iterations; + + // Time (in seconds) spent inside the minimizer loop in the current + // iteration. + double iteration_time_in_seconds; + + // Time (in seconds) spent inside the trust region step solver. + double step_solver_time_in_seconds; + + // Time (in seconds) since the user called Solve(). + double cumulative_time_in_seconds; + }; + +.. class:: IterationCallback + + .. code-block:: c++ + + class IterationCallback { + public: + virtual ~IterationCallback() {} + virtual CallbackReturnType operator()(const IterationSummary& summary) = 0; + }; + + Interface for specifying callbacks that are executed at the end of + each iteration of the Minimizer. The solver uses the return value of + ``operator()`` to decide whether to continue solving or to + terminate. The user can return three values. + + #. ``SOLVER_ABORT`` indicates that the callback detected an abnormal + situation. The solver returns without updating the parameter + blocks (unless ``Solver::Options::update_state_every_iteration`` is + set true). Solver returns with ``Solver::Summary::termination_type`` + set to ``USER_ABORT``. + + #. ``SOLVER_TERMINATE_SUCCESSFULLY`` indicates that there is no need + to optimize anymore (some user specified termination criterion + has been met). Solver returns with + ``Solver::Summary::termination_type``` set to ``USER_SUCCESS``. + + #. ``SOLVER_CONTINUE`` indicates that the solver should continue + optimizing. + + For example, the following ``IterationCallback`` is used internally + by Ceres to log the progress of the optimization. + + .. code-block:: c++ + + class LoggingCallback : public IterationCallback { + public: + explicit LoggingCallback(bool log_to_stdout) + : log_to_stdout_(log_to_stdout) {} + + ~LoggingCallback() {} + + CallbackReturnType operator()(const IterationSummary& summary) { + const char* kReportRowFormat = + "% 4d: f:% 8e d:% 3.2e g:% 3.2e h:% 3.2e " + "rho:% 3.2e mu:% 3.2e eta:% 3.2e li:% 3d"; + string output = StringPrintf(kReportRowFormat, + summary.iteration, + summary.cost, + summary.cost_change, + summary.gradient_max_norm, + summary.step_norm, + summary.relative_decrease, + summary.trust_region_radius, + summary.eta, + summary.linear_solver_iterations); + if (log_to_stdout_) { + cout << output << endl; + } else { + VLOG(1) << output; + } + return SOLVER_CONTINUE; + } + + private: + const bool log_to_stdout_; + }; + + + +:class:`CRSMatrix` +------------------ + +.. class:: CRSMatrix + + .. code-block:: c++ + + struct CRSMatrix { + int num_rows; + int num_cols; + vector<int> cols; + vector<int> rows; + vector<double> values; + }; + + A compressed row sparse matrix used primarily for communicating the + Jacobian matrix to the user. + + A compressed row matrix stores its contents in three arrays, + ``rows``, ``cols`` and ``values``. + + ``rows`` is a ``num_rows + 1`` sized array that points into the ``cols`` and + ``values`` array. For each row ``i``: + + ``cols[rows[i]]`` ... ``cols[rows[i + 1] - 1]`` are the indices of the + non-zero columns of row ``i``. + + ``values[rows[i]]`` ... ``values[rows[i + 1] - 1]`` are the values of the + corresponding entries. + + ``cols`` and ``values`` contain as many entries as there are + non-zeros in the matrix. + + e.g, consider the 3x4 sparse matrix + + .. code-block:: c++ + + 0 10 0 4 + 0 2 -3 2 + 1 2 0 0 + + The three arrays will be: + + .. code-block:: c++ + + -row0- ---row1--- -row2- + rows = [ 0, 2, 5, 7] + cols = [ 1, 3, 1, 2, 3, 0, 1] + values = [10, 4, 2, -3, 2, 1, 2] + + +:class:`Solver::Summary` +------------------------ + +.. class:: Solver::Summary + + Note that all times reported in this struct are wall times. + + .. code-block:: c++ + + struct Summary { + // A brief one line description of the state of the solver after + // termination. + string BriefReport() const; + + // A full multiline description of the state of the solver after + // termination. + string FullReport() const; + + // Minimizer summary ------------------------------------------------- + MinimizerType minimizer_type; + + SolverTerminationType termination_type; + + // If the solver did not run, or there was a failure, a + // description of the error. + string error; + + // Cost of the problem before and after the optimization. See + // problem.h for definition of the cost of a problem. + double initial_cost; + double final_cost; + + // The part of the total cost that comes from residual blocks that + // were held fixed by the preprocessor because all the parameter + // blocks that they depend on were fixed. + double fixed_cost; + + vector<IterationSummary> iterations; + + int num_successful_steps; + int num_unsuccessful_steps; + int num_inner_iteration_steps; + + // When the user calls Solve, before the actual optimization + // occurs, Ceres performs a number of preprocessing steps. These + // include error checks, memory allocations, and reorderings. This + // time is accounted for as preprocessing time. + double preprocessor_time_in_seconds; + + // Time spent in the TrustRegionMinimizer. + double minimizer_time_in_seconds; + + // After the Minimizer is finished, some time is spent in + // re-evaluating residuals etc. This time is accounted for in the + // postprocessor time. + double postprocessor_time_in_seconds; + + // Some total of all time spent inside Ceres when Solve is called. + double total_time_in_seconds; + + double linear_solver_time_in_seconds; + double residual_evaluation_time_in_seconds; + double jacobian_evaluation_time_in_seconds; + double inner_iteration_time_in_seconds; + + // Preprocessor summary. + int num_parameter_blocks; + int num_parameters; + int num_effective_parameters; + int num_residual_blocks; + int num_residuals; + + int num_parameter_blocks_reduced; + int num_parameters_reduced; + int num_effective_parameters_reduced; + int num_residual_blocks_reduced; + int num_residuals_reduced; + + int num_eliminate_blocks_given; + int num_eliminate_blocks_used; + + int num_threads_given; + int num_threads_used; + + int num_linear_solver_threads_given; + int num_linear_solver_threads_used; + + LinearSolverType linear_solver_type_given; + LinearSolverType linear_solver_type_used; + + vector<int> linear_solver_ordering_given; + vector<int> linear_solver_ordering_used; + + bool inner_iterations_given; + bool inner_iterations_used; + + vector<int> inner_iteration_ordering_given; + vector<int> inner_iteration_ordering_used; + + PreconditionerType preconditioner_type; + + TrustRegionStrategyType trust_region_strategy_type; + DoglegType dogleg_type; + + SparseLinearAlgebraLibraryType sparse_linear_algebra_library; + + LineSearchDirectionType line_search_direction_type; + LineSearchType line_search_type; + int max_lbfgs_rank; + }; + + +Covariance Estimation +===================== + +Background +---------- + +One way to assess the quality of the solution returned by a +non-linear least squares solve is to analyze the covariance of the +solution. + +Let us consider the non-linear regression problem + +.. math:: y = f(x) + N(0, I) + +i.e., the observation :math:`y` is a random non-linear function of the +independent variable :math:`x` with mean :math:`f(x)` and identity +covariance. Then the maximum likelihood estimate of :math:`x` given +observations :math:`y` is the solution to the non-linear least squares +problem: + +.. math:: x^* = \arg \min_x \|f(x)\|^2 + +And the covariance of :math:`x^*` is given by + +.. math:: C(x^*) = \left(J'(x^*)J(x^*)\right)^{-1} + +Here :math:`J(x^*)` is the Jacobian of :math:`f` at :math:`x^*`. The +above formula assumes that :math:`J(x^*)` has full column rank. + +If :math:`J(x^*)` is rank deficient, then the covariance matrix :math:`C(x^*)` +is also rank deficient and is given by the Moore-Penrose pseudo inverse. + +.. math:: C(x^*) = \left(J'(x^*)J(x^*)\right)^{\dagger} + +Note that in the above, we assumed that the covariance matrix for +:math:`y` was identity. This is an important assumption. If this is +not the case and we have + +.. math:: y = f(x) + N(0, S) + +Where :math:`S` is a positive semi-definite matrix denoting the +covariance of :math:`y`, then the maximum likelihood problem to be +solved is + +.. math:: x^* = \arg \min_x f'(x) S^{-1} f(x) + +and the corresponding covariance estimate of :math:`x^*` is given by + +.. math:: C(x^*) = \left(J'(x^*) S^{-1} J(x^*)\right)^{-1} + +So, if it is the case that the observations being fitted to have a +covariance matrix not equal to identity, then it is the user's +responsibility that the corresponding cost functions are correctly +scaled, e.g. in the above case the cost function for this problem +should evaluate :math:`S^{-1/2} f(x)` instead of just :math:`f(x)`, +where :math:`S^{-1/2}` is the inverse square root of the covariance +matrix :math:`S`. + +Gauge Invariance +---------------- + +In structure from motion (3D reconstruction) problems, the +reconstruction is ambiguous upto a similarity transform. This is +known as a *Gauge Ambiguity*. Handling Gauges correctly requires the +use of SVD or custom inversion algorithms. For small problems the +user can use the dense algorithm. For more details see the work of +Kanatani & Morris [KanataniMorris]_. + + +:class:`Covariance` +------------------- + +:class:`Covariance` allows the user to evaluate the covariance for a +non-linear least squares problem and provides random access to its +blocks. The computation assumes that the cost functions compute +residuals such that their covariance is identity. + +Since the computation of the covariance matrix requires computing the +inverse of a potentially large matrix, this can involve a rather large +amount of time and memory. However, it is usually the case that the +user is only interested in a small part of the covariance +matrix. Quite often just the block diagonal. :class:`Covariance` +allows the user to specify the parts of the covariance matrix that she +is interested in and then uses this information to only compute and +store those parts of the covariance matrix. + +Rank of the Jacobian +-------------------- + +As we noted above, if the Jacobian is rank deficient, then the inverse +of :math:`J'J` is not defined and instead a pseudo inverse needs to be +computed. + +The rank deficiency in :math:`J` can be *structural* -- columns +which are always known to be zero or *numerical* -- depending on the +exact values in the Jacobian. + +Structural rank deficiency occurs when the problem contains parameter +blocks that are constant. This class correctly handles structural rank +deficiency like that. + +Numerical rank deficiency, where the rank of the matrix cannot be +predicted by its sparsity structure and requires looking at its +numerical values is more complicated. Here again there are two +cases. + + a. The rank deficiency arises from overparameterization. e.g., a + four dimensional quaternion used to parameterize :math:`SO(3)`, + which is a three dimensional manifold. In cases like this, the + user should use an appropriate + :class:`LocalParameterization`. Not only will this lead to better + numerical behaviour of the Solver, it will also expose the rank + deficiency to the :class:`Covariance` object so that it can + handle it correctly. + + b. More general numerical rank deficiency in the Jacobian requires + the computation of the so called Singular Value Decomposition + (SVD) of :math:`J'J`. We do not know how to do this for large + sparse matrices efficiently. For small and moderate sized + problems this is done using dense linear algebra. + + +:class:`Covariance::Options` + +.. class:: Covariance::Options + +.. member:: int Covariance::Options::num_threads + + Default: ``1`` + + Number of threads to be used for evaluating the Jacobian and + estimation of covariance. + +.. member:: CovarianceAlgorithmType Covariance::Options::algorithm_type + + Default: ``SPARSE_QR`` or ``DENSE_SVD`` + + Ceres supports three different algorithms for covariance + estimation, which represent different tradeoffs in speed, accuracy + and reliability. + + 1. ``DENSE_SVD`` uses ``Eigen``'s ``JacobiSVD`` to perform the + computations. It computes the singular value decomposition + + .. math:: U S V^\top = J + + and then uses it to compute the pseudo inverse of J'J as + + .. math:: (J'J)^{\dagger} = V S^{\dagger} V^\top + + It is an accurate but slow method and should only be used for + small to moderate sized problems. It can handle full-rank as + well as rank deficient Jacobians. + + 2. ``SPARSE_CHOLESKY`` uses the ``CHOLMOD`` sparse Cholesky + factorization library to compute the decomposition : + + .. math:: R^\top R = J^\top J + + and then + + .. math:: \left(J^\top J\right)^{-1} = \left(R^\top R\right)^{-1} + + It a fast algorithm for sparse matrices that should be used when + the Jacobian matrix J is well conditioned. For ill-conditioned + matrices, this algorithm can fail unpredictabily. This is + because Cholesky factorization is not a rank-revealing + factorization, i.e., it cannot reliably detect when the matrix + being factorized is not of full + rank. ``SuiteSparse``/``CHOLMOD`` supplies a heuristic for + checking if the matrix is rank deficient (cholmod_rcond), but it + is only a heuristic and can have both false positive and false + negatives. + + Recent versions of ``SuiteSparse`` (>= 4.2.0) provide a much more + efficient method for solving for rows of the covariance + matrix. Therefore, if you are doing ``SPARSE_CHOLESKY``, we strongly + recommend using a recent version of ``SuiteSparse``. + + 3. ``SPARSE_QR`` uses the ``SuiteSparseQR`` sparse QR factorization + library to compute the decomposition + + .. math:: + + QR &= J\\ + \left(J^\top J\right)^{-1} &= \left(R^\top R\right)^{-1} + + It is a moderately fast algorithm for sparse matrices, which at + the price of more time and memory than the ``SPARSE_CHOLESKY`` + algorithm is numerically better behaved and is rank revealing, + i.e., it can reliably detect when the Jacobian matrix is rank + deficient. + + Neither ``SPARSE_CHOLESKY`` or ``SPARSE_QR`` are capable of computing + the covariance if the Jacobian is rank deficient. + +.. member:: int Covariance::Options::min_reciprocal_condition_number + + Default: :math:`10^{-14}` + + If the Jacobian matrix is near singular, then inverting :math:`J'J` + will result in unreliable results, e.g, if + + .. math:: + + J = \begin{bmatrix} + 1.0& 1.0 \\ + 1.0& 1.0000001 + \end{bmatrix} + + which is essentially a rank deficient matrix, we have + + .. math:: + + (J'J)^{-1} = \begin{bmatrix} + 2.0471e+14& -2.0471e+14 \\ + -2.0471e+14 2.0471e+14 + \end{bmatrix} + + + This is not a useful result. Therefore, by default + :func:`Covariance::Compute` will return ``false`` if a rank + deficient Jacobian is encountered. How rank deficiency is detected + depends on the algorithm being used. + + 1. ``DENSE_SVD`` + + .. math:: \frac{\sigma_{\text{min}}}{\sigma_{\text{max}}} < \sqrt{\text{min_reciprocal_condition_number}} + + where :math:`\sigma_{\text{min}}` and + :math:`\sigma_{\text{max}}` are the minimum and maxiumum + singular values of :math:`J` respectively. + + 2. ``SPARSE_CHOLESKY`` + + .. math:: \text{cholmod_rcond} < \text{min_reciprocal_conditioner_number} + + Here cholmod_rcond is a crude estimate of the reciprocal + condition number of :math:`J^\top J` by using the maximum and + minimum diagonal entries of the Cholesky factor :math:`R`. There + are no theoretical guarantees associated with this test. It can + give false positives and negatives. Use at your own risk. The + default value of ``min_reciprocal_condition_number`` has been + set to a conservative value, and sometimes the + :func:`Covariance::Compute` may return false even if it is + possible to estimate the covariance reliably. In such cases, the + user should exercise their judgement before lowering the value + of ``min_reciprocal_condition_number``. + + 3. ``SPARSE_QR`` + + .. math:: \operatorname{rank}(J) < \operatorname{num\_col}(J) + + Here :\math:`\operatorname{rank}(J)` is the estimate of the + rank of `J` returned by the ``SuiteSparseQR`` algorithm. It is + a fairly reliable indication of rank deficiency. + +.. member:: int Covariance::Options::null_space_rank + + When using ``DENSE_SVD``, the user has more control in dealing + with singular and near singular covariance matrices. + + As mentioned above, when the covariance matrix is near singular, + instead of computing the inverse of :math:`J'J`, the Moore-Penrose + pseudoinverse of :math:`J'J` should be computed. + + If :math:`J'J` has the eigen decomposition :math:`(\lambda_i, + e_i)`, where :math:`lambda_i` is the :math:`i^\textrm{th}` + eigenvalue and :math:`e_i` is the corresponding eigenvector, then + the inverse of :math:`J'J` is + + .. math:: (J'J)^{-1} = \sum_i \frac{1}{\lambda_i} e_i e_i' + + and computing the pseudo inverse involves dropping terms from this + sum that correspond to small eigenvalues. + + How terms are dropped is controlled by + `min_reciprocal_condition_number` and `null_space_rank`. + + If `null_space_rank` is non-negative, then the smallest + `null_space_rank` eigenvalue/eigenvectors are dropped irrespective + of the magnitude of :math:`\lambda_i`. If the ratio of the + smallest non-zero eigenvalue to the largest eigenvalue in the + truncated matrix is still below min_reciprocal_condition_number, + then the `Covariance::Compute()` will fail and return `false`. + + Setting `null_space_rank = -1` drops all terms for which + + .. math:: \frac{\lambda_i}{\lambda_{\textrm{max}}} < \textrm{min_reciprocal_condition_number} + + This option has no effect on ``SPARSE_QR`` and ``SPARSE_CHOLESKY`` + algorithms. + +.. member:: bool Covariance::Options::apply_loss_function + + Default: `true` + + Even though the residual blocks in the problem may contain loss + functions, setting ``apply_loss_function`` to false will turn off + the application of the loss function to the output of the cost + function and in turn its effect on the covariance. + +.. class:: Covariance + + :class:`Covariance::Options` as the name implies is used to control + the covariance estimation algorithm. Covariance estimation is a + complicated and numerically sensitive procedure. Please read the + entire documentation for :class:`Covariance::Options` before using + :class:`Covariance`. + +.. function:: bool Covariance::Compute(const vector<pair<const double*, const double*> >& covariance_blocks, Problem* problem) + + Compute a part of the covariance matrix. + + The vector ``covariance_blocks``, indexes into the covariance + matrix block-wise using pairs of parameter blocks. This allows the + covariance estimation algorithm to only compute and store these + blocks. + + Since the covariance matrix is symmetric, if the user passes + ``<block1, block2>``, then ``GetCovarianceBlock`` can be called with + ``block1``, ``block2`` as well as ``block2``, ``block1``. + + ``covariance_blocks`` cannot contain duplicates. Bad things will + happen if they do. + + Note that the list of ``covariance_blocks`` is only used to + determine what parts of the covariance matrix are computed. The + full Jacobian is used to do the computation, i.e. they do not have + an impact on what part of the Jacobian is used for computation. + + The return value indicates the success or failure of the covariance + computation. Please see the documentation for + :class:`Covariance::Options` for more on the conditions under which + this function returns ``false``. + +.. function:: bool GetCovarianceBlock(const double* parameter_block1, const double* parameter_block2, double* covariance_block) const + + Return the block of the covariance matrix corresponding to + ``parameter_block1`` and ``parameter_block2``. + + Compute must be called before the first call to ``GetCovarianceBlock`` + and the pair ``<parameter_block1, parameter_block2>`` OR the pair + ``<parameter_block2, parameter_block1>`` must have been present in the + vector covariance_blocks when ``Compute`` was called. Otherwise + ``GetCovarianceBlock`` will return false. + + ``covariance_block`` must point to a memory location that can store + a ``parameter_block1_size x parameter_block2_size`` matrix. The + returned covariance will be a row-major matrix. + +Example Usage +------------- + +.. code-block:: c++ + + double x[3]; + double y[2]; + + Problem problem; + problem.AddParameterBlock(x, 3); + problem.AddParameterBlock(y, 2); + <Build Problem> + <Solve Problem> + + Covariance::Options options; + Covariance covariance(options); + + vector<pair<const double*, const double*> > covariance_blocks; + covariance_blocks.push_back(make_pair(x, x)); + covariance_blocks.push_back(make_pair(y, y)); + covariance_blocks.push_back(make_pair(x, y)); + + CHECK(covariance.Compute(covariance_blocks, &problem)); + + double covariance_xx[3 * 3]; + double covariance_yy[2 * 2]; + double covariance_xy[3 * 2]; + covariance.GetCovarianceBlock(x, x, covariance_xx) + covariance.GetCovarianceBlock(y, y, covariance_yy) + covariance.GetCovarianceBlock(x, y, covariance_xy) + diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst new file mode 100644 index 0000000..1e5756a --- /dev/null +++ b/docs/source/tutorial.rst @@ -0,0 +1,717 @@ +.. highlight:: c++ + +.. default-domain:: cpp + +.. _chapter-tutorial: + +======== +Tutorial +======== +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 expression +: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]`. 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 be a single parameter. + +:math:`\rho_i` is a :class:`LossFunction`. A :class:`LossFunction` is +a scalar function that is used to reduce the influence of outliers on +the solution of non-linear least squares problems. As a special case, +when :math:`\rho_i(x) = x`, i.e., the identity function, we get the +more familiar `non-linear least squares problem +<http://en.wikipedia.org/wiki/Non-linear_least_squares>`_. + +.. math:: \frac{1}{2}\sum_{i=1} \left\|f_i\left(x_{i_1}, ... ,x_{i_k}\right)\right\|^2. + :label: ceresproblem2 + +In this chapter we will learn how to solve :eq:`ceresproblem` using +Ceres Solver. Full working code for all the examples described in this +chapter and more can be found in the `examples +<https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/>`_ +directory. + +.. _section-hello-world: + +Hello World! +============ + +To get started, 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]_. + +The first step is to write a functor that will evaluate this the +function :math:`f(x) = 10 - x`: + +.. code-block:: c++ + + struct CostFunctor { + template <typename T> + bool operator()(const T* const x, T* residual) const { + residual[0] = T(10.0) - x[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 +``CostFunctor::operator<T>()``, with ``T=double`` when just the +residual is needed, and with a special type ``T=Jet`` when the +Jacobians are needed. In :ref:`section-derivatives` we discuss the +various ways of supplying derivatives to Ceres in more detail. + +Once we have a way of computing the residual function, 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) { + google::InitGoogleLogging(argv[0]); + + // The variable to solve for with its initial value. + double initial_x = 5.0; + double x = initial_x; + + // Build the problem. + Problem problem; + + // Set up the only cost function (also known as residual). This uses + // auto-differentiation to obtain the derivative (jacobian). + CostFunction* cost_function = + new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor); + problem.AddResidualBlock(cost_function, NULL, &x); + + // Run the solver! + Solver::Options options; + 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 : " << initial_x + << " -> " << x << "\n"; + return 0; + } + +:class:`AutoDiffCostFunction` takes a ``CostFunctor`` as input, +automatically differentiates it and gives it a :class:`CostFunction` +interface. + +Compiling and running `examples/helloworld.cc +<https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/helloworld.cc>`_ +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: 6.91e-06 tt: 1.91e-03 + 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: 2.81e-05 tt: 1.99e-03 + 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: 1.00e-05 tt: 2.01e-03 + Ceres Solver Report: Iterations: 2, Initial cost: 1.250000e+01, Final cost: 1.388518e-16, Termination: PARAMETER_TOLERANCE. + x : 5 -> 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] `examples/helloworld.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/helloworld.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-derivatives: + + +Derivatives +=========== + +Ceres Solver like most optimization packages, depends on being able to +evaluate the value and the derivatives of each term in the objective +function at arbitrary parameter values. Doing so correctly and +efficiently is essential to getting good results. Ceres Solver +provides a number of ways of doing so. You have already seen one of +them in action -- +Automatic Differentiation in `examples/helloworld.cc +<https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/helloworld.cc>`_ + +We now consider the other two possibilities. Analytic and numeric +derivatives. + + +Numeric Derivatives +------------------- + +In some cases, its not possible to define a templated cost functor, +for example when the evaluation of the residual involves a call to a +library function that you do not have control over. In such a +situation, numerical differentiation can be used. The user defines a +functor which computes the residual value and construct a +:class:`NumericDiffCostFunction` using it. e.g., for :math:`f(x) = 10 - x` +the corresponding functor would be + +.. code-block:: c++ + + struct NumericDiffCostFunctor { + bool operator()(const double* const x, double* residual) const { + residual[0] = 10.0 - x[0]; + return true; + } + }; + +Which is added to the :class:`Problem` as: + +.. code-block:: c++ + + CostFunction* cost_function = + new NumericDiffCostFunction<NumericDiffCostFunctor, ceres::CENTRAL, 1, 1, 1>( + new NumericDiffCostFunctor) + problem.AddResidualBlock(cost_function, NULL, &x); + +Notice the parallel from when we were using automatic differentiation + +.. code-block:: c++ + + CostFunction* cost_function = + new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor); + problem.AddResidualBlock(cost_function, NULL, &x); + +The construction looks almost identical to the one 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 [#f3]_. For more details see the documentation +for :class:`NumericDiffCostFunction`. + +**Generally speaking we recommend automatic differentiation instead of +numeric differentiation. The use of C++ templates makes automatic +differentiation efficient, whereas numeric differentiation is +expensive, prone to numeric errors, and leads to slower convergence.** + + +Analytic Derivatives +-------------------- + +In some cases, using automatic differentiation is not possible. For +example, it may be the case that it is more efficient to compute the +derivatives in closed form instead of relying on the chain rule used +by the automatic differentiation code. + +In such cases, it is possible to supply your own residual and jacobian +computation code. To do this, define a subclass of +:class:`CostFunction` or :class:`SizedCostFunction` if you know the +sizes of the parameters and residuals at compile time. Here for +example is ``SimpleCostFunction`` that implements :math:`f(x) = 10 - +x`. + +.. code-block:: c++ + + class QuadraticCostFunction : public ceres::SizedCostFunction<1, 1> { + public: + virtual ~QuadraticCostFunction() {} + 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::Evaluate`` is provided with an input array of +``parameters``, an output array ``residuals`` for residuals and an +output array ``jacobians`` for Jacobians. 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 [#f4]_ . + +As can be seen from the above code fragments, implementing +:class:`CostFunction` objects is a bit tedious. We recommend that +unless you have a good reason to manage the jacobian computation +yourself, you use :class:`AutoDiffCostFunction` or +:class:`NumericDiffCostFunction` to construct your residual blocks. + +More About Derivatives +---------------------- + +Computing derivatives is by far the most complicated part of using +Ceres, and depending on the circumstance the user may need more +sophisticated ways of computing derivatives. This section just +scratches the surface of how derivatives can be supplied to +Ceres. Once you are comfortable with using +:class:`NumericDiffCostFunction` and :class:`AutoDiffCostFunction` we +recommend taking a look at :class:`DynamicAutoDiffCostFunction`, +:class:`CostFunctionToFunctor`, :class:`NumericDiffFunctor` and +:class:`ConditionedCostFunction` for more advanced ways of +constructing and computing cost functions. + +.. rubric:: Footnotes + +.. [#f3] `examples/helloworld_numeric_diff.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/helloworld_numeric_diff.cc>`_. + +.. [#f4] `examples/helloworld_analytic_diff.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/helloworld_analytic_diff.cc>`_. + + +.. _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, has four residuals +and we wish to find :math:`x` such that :math:`\frac{1}{2}\|F(x)\|^2` +is minimized. + +Again, the first step is to define functors that evaluate of the terms +in the objective functor. Here is the code for evaluating +:math:`f_4(x_1, x_4)`: + +.. code-block:: c++ + + struct F4 { + 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; + } + }; + + +Similarly, we can define classes ``F1``, ``F2`` and ``F4`` to evaluate +:math:`f_1(x_1, x_2)`, :math:`f_2(x_3, x_4)` and :math:`f_3(x_2, x_3)` +respectively. Using these, the problem can be constructed as follows: + + +.. code-block:: c++ + + double x1 = 3.0; double x2 = -1.0; double x3 = 0.0; double x4 = 1.0; + + Problem problem; + + // Add residual terms to the problem using the using the autodiff + // wrapper to get the derivatives automatically. + problem.AddResidualBlock( + new AutoDiffCostFunction<F1, 1, 1, 1>(new F1), NULL, &x1, &x2); + problem.AddResidualBlock( + new AutoDiffCostFunction<F2, 1, 1, 1>(new F2), NULL, &x3, &x4); + problem.AddResidualBlock( + new AutoDiffCostFunction<F3, 1, 1, 1>(new F3), NULL, &x2, &x3) + problem.AddResidualBlock( + new AutoDiffCostFunction<F4, 1, 1, 1>(new F4), NULL, &x1, &x4); + + +Note that each ``ResidualBlock`` only depends on the two parameters +that the corresponding residual object depends on and not on all four +parameters. Compiling and running `examples/powell.cc +<https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/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}`. + + +.. rubric:: Footnotes + +.. [#f5] `examples/powell.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/powell.cc>`_. + + +.. _section-fitting: + +Curve Fitting +============= + +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 +[#f6]_. 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++ + + struct ExponentialResidual { + 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 +:class:`CostFunction` for every observation. + + +.. code-block:: c++ + + double m = 0.0; + double c = 0.0; + + Problem problem; + for (int i = 0; i < kNumObservations; ++i) { + CostFunction* cost_function = + new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>( + new ExponentialResidual(data[2 * i], data[2 * i + 1])); + problem.AddResidualBlock(cost_function, NULL, &m, &c); + } + +Compiling and running `examples/curve_fitting.cc +<https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/curve_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:: least_squares_fit.png + :figwidth: 500px + :height: 400px + :align: center + + Least squares curve fitting. + + +.. rubric:: Footnotes + +.. [#f6] `examples/curve_fitting.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/curve_fitting.cc>`_ + + +Robust Curve Fitting +===================== + +Now suppose the data we are given has some outliers, i.e., we have +some points that do not obey the noise model. If we were to use the +code above to fit such data, we would get a fit that looks as +below. Notice how the fitted curve deviates from the ground truth. + +.. figure:: non_robust_least_squares_fit.png + :figwidth: 500px + :height: 400px + :align: center + +To deal with outliers, a standard technique is to use a +:class:`LossFunction`. Loss functions, reduce the influence of +residual blocks with high residuals, usually the ones corresponding to +outliers. To associate a loss function in a residual block, we change + +.. code-block:: c++ + + problem.AddResidualBlock(cost_function, NULL , &m, &c); + +to + +.. code-block:: c++ + + problem.AddResidualBlock(cost_function, new CauchyLoss(0.5) , &m, &c); + +:class:`CauchyLoss` is one of the loss functions that ships with Ceres +Solver. The argument :math:`0.5` specifies the scale of the loss +function. As a result, we get the fit below [#f7]_. Notice how the +fitted curve moves back closer to the ground truth curve. + +.. figure:: robust_least_squares_fit.png + :figwidth: 500px + :height: 400px + :align: center + + Using :class:`LossFunction` to reduce the effect of outliers on a + least squares fit. + + +.. rubric:: Footnotes + +.. [#f7] `examples/robust_curve_fitting.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/robust_curve_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 solve a problem from the `BAL +<http://grail.cs.washington.edu/projects/bal/>`_ dataset [#f8]_. + +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 +are: three for rotation as a Rodriques' axis-angle vector, three +for translation, one for focal length and two for radial distortion. +The details of this camera model can be found the `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; + } + + // Factory to hide the construction of the CostFunction object from + // the client code. + static ceres::CostFunction* Create(const double observed_x, + const double observed_y) { + return (new ceres::AutoDiffCostFunction<SnavelyReprojectionError, 2, 9, 3>( + new SnavelyReprojectionError(observed_x, observed_y))); + } + + 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 life much simpler. The function +:func:`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++ + + ceres::Problem problem; + for (int i = 0; i < bal_problem.num_observations(); ++i) { + 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)); + } + + +Notice that the problem construction for bundle adjustment is very +similar to the curve fitting example -- one term is added to the +objective function per observation. + +Since this large sparse problem (well large for ``DENSE_QR`` anyways), +one way to solve this problem is to set +:member:`Solver::Options::linear_solver_type` to +``SPARSE_NORMAL_CHOLESKY`` and call :member:`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 +<https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/bundle_adjuster.cc>`_ + + +.. rubric:: Footnotes + +.. [#f8] `examples/simple_bundle_adjuster.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/simple_bundle_adjuster.cc>`_ + + +Other Examples +============== + +Besides the examples in this chapter, the `example +<https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/>`_ +directory contains a number of other examples: + +#. `bundle_adjuster.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/bundle_adjuster.cc>`_ + shows how to use the various features of Ceres to solve bundle + adjustment problems. + +#. `circle_fit.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/circle_fit.cc>`_ + shows how to fit data to a circle. + +#. `denoising.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/denoising.cc>`_ + implements image denoising using the `Fields of Experts + <http://www.gris.informatik.tu-darmstadt.de/~sroth/research/foe/index.html>`_ + model. + +#. `nist.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/nist.cc>`_ + implements and attempts to solves the `NIST + <http://www.itl.nist.gov/div898/strd/nls/nls_main.shtm>`_ + non-linear regression problems. + +#. `libmv_bundle_adjuster.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/libmv_bundle_adjuster.cc>`_ + is the bundle adjustment algorithm used by `Blender <www.blender.org>`_/libmv. + + diff --git a/docs/source/version_history.rst b/docs/source/version_history.rst new file mode 100644 index 0000000..5e1a150 --- /dev/null +++ b/docs/source/version_history.rst @@ -0,0 +1,714 @@ +.. _chapter-version-history: + +=============== +Version History +=============== + +1.7.0 +===== + +New Features +------------ + +#. Sparse and dense covariance estimation. +#. A new Wolfe line search. (Alex Stewart) +#. ``BFGS`` line search direction. (Alex Stewart) +#. C API +#. Speeded up the use of loss functions > 17x. +#. Use of Inner iterations can now be adaptively stopped. Iteration + and runtime statistics for inner iterations are not reported in + ``Solver::Summary`` and ``Solver::Summary::FullReport``. +#. Add BlockRandomAccessCRSMatrix. +#. Speeded up automatic differentiation by 7\%. +#. Bundle adjustment example from libmv/Blender (Sergey Sharybin) +#. Add the ability to turn shared library compilation on and off +#. No more dependence on Protocol Buffers. +#. Incomplete LQ factorization. +#. Ability to write trust region problems to disk. +#. Add sinh, cosh, tanh and tan functions to automatic differentiation + (Johannes Schönberger) + +Bug Fixes +--------- + +#. Add documentation for minimizer progress output. +#. Lint and other cleanups (William Rucklidge and James Roseborough) +#. Collections port fix for MSC 2008 (Sergey Sharybin) +#. Various corrections and cleanups in the documentation. +#. Change the path where CeresConfig.cmake is installed (Pablo + Speciale) +#. Minor erros in documentation (Pablo Speciale) +#. Updated depend.cmake to follow CMake IF convention. (Joydeep + Biswas) +#. Stablize the schur ordering algorithm. +#. Update license header in split.h. +#. Enabling -O4 (link-time optimization) only if compiler/linker + support it. (Alex Stewart) +#. Consistent glog path across files. +#. ceres-solver.spec: Use cleaner, more conventional Release string + (Taylor Braun-Jones) +#. Fix compile bug on RHEL6 due to missing header (Taylor Braun-Jones) +#. CMake file is less verbose. +#. Use the latest upstream version of google-test and gmock. +#. Rationalize some of the variable names in ``Solver::Options``. +#. Improve Summary::FullReport when line search is used. +#. Expose line search parameters in ``Solver::Options``. +#. Fix update of L-BFGS history buffers after they become full. (Alex + Stewart) +#. Fix configuration error on systems without SuiteSparse installed + (Sergey Sharybin) +#. Enforce the read call returns correct value in ``curve_fitting_c.c`` + (Arnaud Gelas) +#. Fix DynamicAutoDiffCostFunction (Richard Stebbing) +#. Fix Problem::RemoveParameterBlock documentation (Johannes + Schönberger) +#. Fix a logging bug in parameter_block.h +#. Refactor the preconditioner class structure. +#. Fix an uninitialized variable warning when building with ``GCC``. +#. Fix a reallocation bug in + ``CreateJacobianBlockSparsityTranspose``. (Yuliy Schwartzburg) +#. Add a define for O_BINARY. + + +1.6.0 +===== + +New Features +------------ + +#. Major Performance improvements. + + a. Schur type solvers (``SPARSE_SCHUR``, ``DENSE_SCHUR``, + ``ITERATIVE_SCHUR``) are significantly faster due to custom BLAS + routines and fewer heap allocations. + + b. ``SPARSE_SCHUR`` when used with ``CX_SPARSE`` now uses a block + AMD for much improved factorization performance. + + c. The jacobian matrix is pre-ordered so that + ``SPARSE_NORMAL_CHOLESKY`` and ``SPARSE_SCHUR`` do not have to + make copies inside ``CHOLMOD``. + + d. Faster autodiff by replacing division by multplication by inverse. + + e. When compiled without threads, the schur eliminator does not pay + the penalty for locking and unlocking mutexes. + +#. Users can now use ``linear_solver_ordering`` to affect the + fill-reducing ordering used by ``SUITE_SPARSE`` for + ``SPARSE_NORMAL_CHOLESKY``. + +#. ``Problem`` can now report the set of parameter blocks it knows about. + +#. ``TrustRegionMinimizer`` uses the evaluator to compute the gradient + instead of a matrix vector multiply. + +#. On ``Mac OS``, whole program optimization is enabled. + +#. Users can now use automatic differentiation to define new + ``LocalParameterization`` objects. (Sergey Sharybin) + +#. Enable larger tuple sizes for Visual Studio 2012. (Petter Strandmark) + + +Bug Fixes +--------- + +#. Update the documentation for ``CostFunction``. +#. Fixed a typo in the documentation. (Pablo Speciale) +#. Fix a typo in suitesparse.cc. +#. Bugfix in ``NumericDiffCostFunction``. (Nicolas Brodu) +#. Death to BlockSparseMatrixBase. +#. Change Minimizer::Options::min_trust_region_radius to double. +#. Update to compile with stricter gcc checks. (Joydeep Biswas) +#. Do not modify cached CMAKE_CXX_FLAGS_RELEASE. (Sergey Sharybin) +#. ``<iterator>`` needed for back_insert_iterator. (Petter Strandmark) +#. Lint cleanup. (William Rucklidge) +#. Documentation corrections. (Pablo Speciale) + + +1.5.0 +===== + +Backward Incompatible API Changes +--------------------------------- + +#. Added ``Problem::Evaluate``. Now you can evaluate a problem or any + part of it without calling the solver. + + In light of this the following settings have been deprecated and + removed from the API. + + - ``Solver::Options::return_initial_residuals`` + - ``Solver::Options::return_initial_gradient`` + - ``Solver::Options::return_initial_jacobian`` + - ``Solver::Options::return_final_residuals`` + - ``Solver::Options::return_final_gradient`` + - ``Solver::Options::return_final_jacobian`` + + Instead we recommend using something like this. + + .. code-block:: c++ + + Problem problem; + // Build problem + + vector<double> initial_residuals; + problem.Evaluate(Problem::EvaluateOptions(), + NULL, /* No cost */ + &initial_residuals, + NULL, /* No gradient */ + NULL /* No jacobian */ ); + + Solver::Options options; + Solver::Summary summary; + Solver::Solve(options, &problem, &summary); + + vector<double> final_residuals; + problem.Evaluate(Problem::EvaluateOptions(), + NULL, /* No cost */ + &final_residuals, + NULL, /* No gradient */ + NULL /* No jacobian */ ); + + +New Features +------------ +#. Problem now supports removal of ParameterBlocks and + ResidualBlocks. There is a space/time tradeoff in doing this which + is controlled by + ``Problem::Options::enable_fast_parameter_block_removal``. + +#. 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. + +#. Added ``Problem::Evaluate``. Now you can evaluate a problem or any + part of it without calling the solver. In light of this the + following settings have been deprecated and removed from the API. + + - ``Solver::Options::return_initial_residuals`` + - ``Solver::Options::return_initial_gradient`` + - ``Solver::Options::return_initial_jacobian`` + - ``Solver::Options::return_final_residuals`` + - ``Solver::Options::return_final_gradient`` + - ``Solver::Options::return_final_jacobian`` + +#. New, much improved HTML documentation using Sphinx. + +#. 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 to the API. + +#. Sped up the robust loss function correction logic when residual is + one dimensional. + +#. Sped up ``DenseQRSolver`` by changing the way dense jacobians are + stored. This is a 200-500% improvement in linear solver performance + depending on the size of the problem. + +#. ``DENSE_SCHUR`` now supports multi-threading. + +#. Greatly expanded ``Summary::FullReport``: + + - Report the ordering used by the ``LinearSolver``. + - Report the ordering used by the inner iterations. + - Execution timing breakdown into evaluations and linear solves. + - Effective size of the problem solved by the solver, which now + accounts for the size of the tangent space when using a + ``LocalParameterization``. + +#. Ceres when run at the ``VLOG`` level 3 or higher will report + detailed timing information about its internals. + +#. Remove extraneous initial and final residual evaluations. This + speeds up the solver a bit. + +#. Automatic differenatiation with a dynamic number of parameter + blocks. (Based on an idea by Thad Hughes). + +#. Sped up problem construction and destruction. + +#. Added matrix adapters to ``rotation.h`` so that the rotation matrix + routines can work with row and column major matrices. (Markus Moll) + +#. ``SCHUR_JACOBI`` can now be used without ``SuiteSparse``. + +#. A ``.spec`` file for producing RPMs. (Taylor Braun-Jones) + +#. ``CMake`` can now build the sphinx documentation (Pablo Speciale) + +#. Add support for creating a CMake config file during build to make + embedding Ceres in other CMake-using projects easier. (Pablo + Speciale). + +#. Better error reporting in ``Problem`` for missing parameter blocks. + +#. A more flexible ``Android.mk`` and a more modular build. If binary + size and/or compile time is a concern, larger parts of the solver + can be disabled at compile time. + +Bug Fixes +--------- +#. Compilation fixes for MSVC2010 (Sergey Sharybin) + +#. Fixed "deprecated conversion from string constant to char*" + warnings. (Pablo Speciale) + +#. Correctly propagate ifdefs when building without Schur eliminator + template specializations. + +#. Correct handling of ``LIB_SUFFIX`` on Linux. (Yuliy Schwartzburg). + +#. Code and signature cleanup in ``rotation.h``. + +#. Make examples independent of internal code. + +#. Disable unused member in ``gtest`` which results in build error on + OS X with latest Xcode. (Taylor Braun-Jones) + +#. Pass the correct flags to the linker when using + ``pthreads``. (Taylor Braun-Jones) + +#. Only use ``cmake28`` macro when building on RHEL6. (Taylor + Braun-Jones) + +#. Remove ``-Wno-return-type-c-linkage`` when compiling with + GCC. (Taylor Braun-Jones) + +#. Fix ``No previous prototype`` warnings. (Sergey Sharybin) + +#. MinGW build fixes. (Sergey Sharybin) + +#. Lots of minor code and lint fixes. (William Rucklidge) + +#. Fixed a bug in ``solver_impl.cc`` residual evaluation. (Markus + Moll) + +#. Fixed varidic evaluation bug in ``AutoDiff``. + +#. Fixed ``SolverImpl`` tests. + +#. Fixed a bug in ``DenseSparseMatrix::ToDenseMatrix()``. + +#. Fixed an initialization bug in ``ProgramEvaluator``. + +#. Fixes to Android.mk paths (Carlos Hernandez) + +#. Modify ``nist.cc`` to compute accuracy based on ground truth + solution rather than the ground truth function value. + +#. Fixed a memory leak in ``cxsparse.cc``. (Alexander Mordvintsev). + +#. Fixed the install directory for libraries by correctly handling + ``LIB_SUFFIX``. (Taylor Braun-Jones) + +1.4.0 +===== + +Backward Incompatible 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. |
