Update ceres to the latest version in google3.
Change-Id: I0165fffa55f60714f23e0096eac89fa68df75a05
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
new file mode 100644
index 0000000..7dad673
--- /dev/null
+++ b/docs/source/least_squares_fit.png
Binary files differ
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
new file mode 100644
index 0000000..9f98d00
--- /dev/null
+++ b/docs/source/loss.png
Binary files differ
diff --git a/docs/source/modeling.rst b/docs/source/modeling.rst
new file mode 100644
index 0000000..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
new file mode 100644
index 0000000..643d162
--- /dev/null
+++ b/docs/source/non_robust_least_squares_fit.png
Binary files differ
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
new file mode 100644
index 0000000..89003c9
--- /dev/null
+++ b/docs/source/robust_least_squares_fit.png
Binary files differ
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.
