nnls_solving.rst 105 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742
  1. .. highlight:: c++
  2. .. default-domain:: cpp
  3. .. cpp:namespace:: ceres
  4. .. _chapter-nnls_solving:
  5. ================================
  6. Solving Non-linear Least Squares
  7. ================================
  8. Introduction
  9. ============
  10. Effective use of Ceres Solver requires some familiarity with the basic
  11. components of a non-linear least squares solver, so before we describe
  12. how to configure and use the solver, we will take a brief look at how
  13. some of the core optimization algorithms in Ceres Solver work.
  14. Let :math:`x \in \mathbb{R}^n` be an :math:`n`-dimensional vector of
  15. variables, and
  16. :math:`F(x) = \left[f_1(x), ... , f_{m}(x) \right]^{\top}` be a
  17. :math:`m`-dimensional function of :math:`x`. We are interested in
  18. solving the optimization problem [#f1]_
  19. .. math:: \arg \min_x \frac{1}{2}\|F(x)\|^2\ . \\
  20. L \le x \le U
  21. :label: nonlinsq
  22. Where, :math:`L` and :math:`U` are vector lower and upper bounds on
  23. the parameter vector :math:`x`. The inequality holds component-wise.
  24. Since the efficient global minimization of :eq:`nonlinsq` for
  25. general :math:`F(x)` is an intractable problem, we will have to settle
  26. for finding a local minimum.
  27. In the following, the Jacobian :math:`J(x)` of :math:`F(x)` is an
  28. :math:`m\times n` matrix, where :math:`J_{ij}(x) = D_j f_i(x)`
  29. and the gradient vector is :math:`g(x) = \nabla \frac{1}{2}\|F(x)\|^2
  30. = J(x)^\top F(x)`.
  31. The general strategy when solving non-linear optimization problems is
  32. to solve a sequence of approximations to the original problem
  33. [NocedalWright]_. At each iteration, the approximation is solved to
  34. determine a correction :math:`\Delta x` to the vector :math:`x`. For
  35. non-linear least squares, an approximation can be constructed by using
  36. the linearization :math:`F(x+\Delta x) \approx F(x) + J(x)\Delta x`,
  37. which leads to the following linear least squares problem:
  38. .. math:: \min_{\Delta x} \frac{1}{2}\|J(x)\Delta x + F(x)\|^2
  39. :label: linearapprox
  40. Unfortunately, naively solving a sequence of these problems and
  41. updating :math:`x \leftarrow x+ \Delta x` leads to an algorithm that
  42. may not converge. To get a convergent algorithm, we need to control
  43. the size of the step :math:`\Delta x`. Depending on how the size of
  44. the step :math:`\Delta x` is controlled, non-linear optimization
  45. algorithms can be divided into two major categories [NocedalWright]_.
  46. 1. **Trust Region** The trust region approach approximates the
  47. objective function using a model function (often a quadratic) over
  48. a subset of the search space known as the trust region. If the
  49. model function succeeds in minimizing the true objective function
  50. the trust region is expanded; conversely, otherwise it is
  51. contracted and the model optimization problem is solved again.
  52. 2. **Line Search** The line search approach first finds a descent
  53. direction along which the objective function will be reduced and
  54. then computes a step size that decides how far should move along
  55. that direction. The descent direction can be computed by various
  56. methods, such as gradient descent, Newton's method and Quasi-Newton
  57. method. The step size can be determined either exactly or
  58. inexactly.
  59. Trust region methods are in some sense dual to line search methods:
  60. trust region methods first choose a step size (the size of the trust
  61. region) and then a step direction while line search methods first
  62. choose a step direction and then a step size. Ceres Solver implements
  63. multiple algorithms in both categories.
  64. .. _section-trust-region-methods:
  65. Trust Region Methods
  66. ====================
  67. The basic trust region algorithm looks something like this.
  68. 1. Given an initial point :math:`x` and a trust region radius :math:`\mu`.
  69. 2. Solve
  70. .. math::
  71. \arg \min_{\Delta x}& \frac{1}{2}\|J(x)\Delta x + F(x)\|^2 \\
  72. \text{such that} &\|D(x)\Delta x\|^2 \le \mu\\
  73. &L \le x + \Delta x \le U.
  74. 3. :math:`\rho = \frac{\displaystyle \|F(x + \Delta x)\|^2 -
  75. \|F(x)\|^2}{\displaystyle \|J(x)\Delta x + F(x)\|^2 -
  76. \|F(x)\|^2}`
  77. 4. if :math:`\rho > \epsilon` then :math:`x = x + \Delta x`.
  78. 5. if :math:`\rho > \eta_1` then :math:`\mu = 2 \mu`
  79. 6. else if :math:`\rho < \eta_2` then :math:`\mu = 0.5 * \mu`
  80. 7. Go to 2.
  81. Here, :math:`\mu` is the trust region radius, :math:`D(x)` is some
  82. matrix used to define a metric on the domain of :math:`F(x)` and
  83. :math:`\rho` measures the quality of the step :math:`\Delta x`, i.e.,
  84. how well did the linear model predict the decrease in the value of the
  85. non-linear objective. The idea is to increase or decrease the radius
  86. of the trust region depending on how well the linearization predicts
  87. the behavior of the non-linear objective, which in turn is reflected
  88. in the value of :math:`\rho`.
  89. The key computational step in a trust-region algorithm is the solution
  90. of the constrained optimization problem
  91. .. math::
  92. \arg \min_{\Delta x}&\quad \frac{1}{2}\|J(x)\Delta x + F(x)\|^2 \\
  93. \text{such that} &\quad \|D(x)\Delta x\|^2 \le \mu\\
  94. &\quad L \le x + \Delta x \le U.
  95. :label: trp
  96. There are a number of different ways of solving this problem, each
  97. giving rise to a different concrete trust-region algorithm. Currently,
  98. Ceres implements two trust-region algorithms - Levenberg-Marquardt
  99. and Dogleg, each of which is augmented with a line search if bounds
  100. constraints are present [Kanzow]_. The user can choose between them by
  101. setting :member:`Solver::Options::trust_region_strategy_type`.
  102. .. rubric:: Footnotes
  103. .. [#f1] At the level of the non-linear solver, the block structure is
  104. not relevant, therefore our discussion here is in terms of an
  105. optimization problem defined over a state vector of size
  106. :math:`n`. Similarly the presence of loss functions is also
  107. ignored as the problem is internally converted into a pure
  108. non-linear least squares problem.
  109. .. _section-levenberg-marquardt:
  110. Levenberg-Marquardt
  111. -------------------
  112. The Levenberg-Marquardt algorithm [Levenberg]_ [Marquardt]_ is the
  113. most popular algorithm for solving non-linear least squares problems.
  114. It was also the first trust region algorithm to be developed
  115. [Levenberg]_ [Marquardt]_. Ceres implements an exact step [Madsen]_
  116. and an inexact step variant of the Levenberg-Marquardt algorithm
  117. [WrightHolt]_ [NashSofer]_.
  118. It can be shown, that the solution to :eq:`trp` can be obtained by
  119. solving an unconstrained optimization of the form
  120. .. math:: \arg\min_{\Delta x} \frac{1}{2}\|J(x)\Delta x + F(x)\|^2 +\lambda \|D(x)\Delta x\|^2
  121. :label: lsqr-naive
  122. Where, :math:`\lambda` is a Lagrange multiplier that is inversely
  123. related to :math:`\mu`. In Ceres, we solve for
  124. .. math:: \arg\min_{\Delta x} \frac{1}{2}\|J(x)\Delta x + F(x)\|^2 + \frac{1}{\mu} \|D(x)\Delta x\|^2
  125. :label: lsqr
  126. The matrix :math:`D(x)` is a non-negative diagonal matrix, typically
  127. the square root of the diagonal of the matrix :math:`J(x)^\top J(x)`.
  128. Before going further, let us make some notational simplifications.
  129. We will assume that the matrix :math:`\frac{1}{\sqrt{\mu}} D` has been
  130. concatenated at the bottom of the matrix :math:`J(x)` and a
  131. corresponding vector of zeroes has been added to the bottom of
  132. :math:`F(x)`, i.e.:
  133. .. math:: J(x) = \begin{bmatrix} J(x) \\ \frac{1}{\sqrt{\mu}} D
  134. \end{bmatrix},\quad F(x) = \begin{bmatrix} F(x) \\ 0
  135. \end{bmatrix}.
  136. This allows us to re-write :eq:`lsqr` as
  137. .. math:: \min_{\Delta x} \frac{1}{2} \|J(x)\Delta x + F(x)\|^2 .
  138. :label: simple
  139. and only talk about :math:`J(x)` and :math:`F(x)` going forward.
  140. For all but the smallest problems the solution of :eq:`simple` in each
  141. iteration of the Levenberg-Marquardt algorithm is the dominant
  142. computational cost. Ceres provides a number of different options for
  143. solving :eq:`simple`. There are two major classes of methods -
  144. factorization and iterative.
  145. The factorization methods are based on computing an exact solution of
  146. :eq:`lsqr` using a Cholesky or a QR factorization and lead to the so
  147. called exact step Levenberg-Marquardt algorithm. But it is not clear
  148. if an exact solution of :eq:`lsqr` is necessary at each step of the
  149. Levenberg-Mardquardt algorithm. We have already seen evidence that
  150. this may not be the case, as :eq:`lsqr` is itself a regularized
  151. version of :eq:`linearapprox`. Indeed, it is possible to construct
  152. non-linear optimization algorithms in which the linearized problem is
  153. solved approximately. These algorithms are known as inexact Newton or
  154. truncated Newton methods [NocedalWright]_.
  155. An inexact Newton method requires two ingredients. First, a cheap
  156. method for approximately solving systems of linear
  157. equations. Typically an iterative linear solver like the Conjugate
  158. Gradients method is used for this purpose [NocedalWright]_. Second, a
  159. termination rule for the iterative solver. A typical termination rule
  160. is of the form
  161. .. math:: \|H(x) \Delta x + g(x)\| \leq \eta_k \|g(x)\|.
  162. :label: inexact
  163. Here, :math:`k` indicates the Levenberg-Marquardt iteration number and
  164. :math:`0 < \eta_k <1` is known as the forcing sequence. [WrightHolt]_
  165. prove that a truncated Levenberg-Marquardt algorithm that uses an
  166. inexact Newton step based on :eq:`inexact` converges for any
  167. sequence :math:`\eta_k \leq \eta_0 < 1` and the rate of convergence
  168. depends on the choice of the forcing sequence :math:`\eta_k`.
  169. Ceres supports both exact and inexact step solution strategies. When
  170. the user chooses a factorization based linear solver, the exact step
  171. Levenberg-Marquardt algorithm is used. When the user chooses an
  172. iterative linear solver, the inexact step Levenberg-Marquardt
  173. algorithm is used.
  174. We will talk more about the various linear solvers that you can use in
  175. :ref:`section-linear-solver`.
  176. .. _section-dogleg:
  177. Dogleg
  178. ------
  179. Another strategy for solving the trust region problem :eq:`trp` was
  180. introduced by
  181. `M. J. D. Powell <https://en.wikipedia.org/wiki/Michael_J._D._Powell>`_. The
  182. key idea there is to compute two vectors
  183. .. math::
  184. \Delta x^{\text{Gauss-Newton}} &= \arg \min_{\Delta x}\frac{1}{2} \|J(x)\Delta x + f(x)\|^2.\\
  185. \Delta x^{\text{Cauchy}} &= -\frac{\|g(x)\|^2}{\|J(x)g(x)\|^2}g(x).
  186. Note that the vector :math:`\Delta x^{\text{Gauss-Newton}}` is the
  187. solution to :eq:`linearapprox` and :math:`\Delta
  188. x^{\text{Cauchy}}` is the vector that minimizes the linear
  189. approximation if we restrict ourselves to moving along the direction
  190. of the gradient. Dogleg methods finds a vector :math:`\Delta x`
  191. defined by :math:`\Delta x^{\text{Gauss-Newton}}` and :math:`\Delta
  192. x^{\text{Cauchy}}` that solves the trust region problem. Ceres
  193. supports two variants that can be chose by setting
  194. :member:`Solver::Options::dogleg_type`.
  195. ``TRADITIONAL_DOGLEG`` as described by Powell, constructs two line
  196. segments using the Gauss-Newton and Cauchy vectors and finds the point
  197. farthest along this line shaped like a dogleg (hence the name) that is
  198. contained in the trust-region. For more details on the exact reasoning
  199. and computations, please see Madsen et al [Madsen]_.
  200. ``SUBSPACE_DOGLEG`` is a more sophisticated method that considers the
  201. entire two dimensional subspace spanned by these two vectors and finds
  202. the point that minimizes the trust region problem in this subspace
  203. [ByrdSchnabel]_.
  204. The key advantage of the Dogleg over Levenberg-Marquardt is that if
  205. the step computation for a particular choice of :math:`\mu` does not
  206. result in sufficient decrease in the value of the objective function,
  207. Levenberg-Marquardt solves the linear approximation from scratch with
  208. a smaller value of :math:`\mu`. Dogleg on the other hand, only needs
  209. to compute the interpolation between the Gauss-Newton and the Cauchy
  210. vectors, as neither of them depend on the value of :math:`\mu`. As a
  211. result the Dogleg method only solves one linear system per successful
  212. step, while Levenberg-Marquardt may need to solve an arbitrary number
  213. of linear systems before it can make progress [LourakisArgyros]_.
  214. A disadvantage of the Dogleg implementation in Ceres Solver is that is
  215. can only be used with method can only be used with exact factorization
  216. based linear solvers.
  217. .. _section-inner-iterations:
  218. Inner Iterations
  219. ----------------
  220. Some non-linear least squares problems have additional structure in
  221. the way the parameter blocks interact that it is beneficial to modify
  222. the way the trust region step is computed. For example, consider the
  223. following regression problem
  224. .. math:: y = a_1 e^{b_1 x} + a_2 e^{b_3 x^2 + c_1}
  225. Given a set of pairs :math:`\{(x_i, y_i)\}`, the user wishes to estimate
  226. :math:`a_1, a_2, b_1, b_2`, and :math:`c_1`.
  227. Notice that the expression on the left is linear in :math:`a_1` and
  228. :math:`a_2`, and given any value for :math:`b_1, b_2` and :math:`c_1`,
  229. it is possible to use linear regression to estimate the optimal values
  230. of :math:`a_1` and :math:`a_2`. It's possible to analytically
  231. eliminate the variables :math:`a_1` and :math:`a_2` from the problem
  232. entirely. Problems like these are known as separable least squares
  233. problem and the most famous algorithm for solving them is the Variable
  234. Projection algorithm invented by Golub & Pereyra [GolubPereyra]_.
  235. Similar structure can be found in the matrix factorization with
  236. missing data problem. There the corresponding algorithm is known as
  237. Wiberg's algorithm [Wiberg]_.
  238. Ruhe & Wedin present an analysis of various algorithms for solving
  239. separable non-linear least squares problems and refer to *Variable
  240. Projection* as Algorithm I in their paper [RuheWedin]_.
  241. Implementing Variable Projection is tedious and expensive. Ruhe &
  242. Wedin present a simpler algorithm with comparable convergence
  243. properties, which they call Algorithm II. Algorithm II performs an
  244. additional optimization step to estimate :math:`a_1` and :math:`a_2`
  245. exactly after computing a successful Newton step.
  246. This idea can be generalized to cases where the residual is not
  247. linear in :math:`a_1` and :math:`a_2`, i.e.,
  248. .. math:: y = f_1(a_1, e^{b_1 x}) + f_2(a_2, e^{b_3 x^2 + c_1})
  249. In this case, we solve for the trust region step for the full problem,
  250. and then use it as the starting point to further optimize just `a_1`
  251. and `a_2`. For the linear case, this amounts to doing a single linear
  252. least squares solve. For non-linear problems, any method for solving
  253. the :math:`a_1` and :math:`a_2` optimization problems will do. The
  254. only constraint on :math:`a_1` and :math:`a_2` (if they are two
  255. different parameter block) is that they do not co-occur in a residual
  256. block.
  257. This idea can be further generalized, by not just optimizing
  258. :math:`(a_1, a_2)`, but decomposing the graph corresponding to the
  259. Hessian matrix's sparsity structure into a collection of
  260. non-overlapping independent sets and optimizing each of them.
  261. Setting :member:`Solver::Options::use_inner_iterations` to ``true``
  262. enables the use of this non-linear generalization of Ruhe & Wedin's
  263. Algorithm II. This version of Ceres has a higher iteration
  264. complexity, but also displays better convergence behavior per
  265. iteration.
  266. Setting :member:`Solver::Options::num_threads` to the maximum number
  267. possible is highly recommended.
  268. .. _section-non-monotonic-steps:
  269. Non-monotonic Steps
  270. -------------------
  271. Note that the basic trust-region algorithm described in
  272. :ref:`section-trust-region-methods` is a descent algorithm in that it
  273. only accepts a point if it strictly reduces the value of the objective
  274. function.
  275. Relaxing this requirement allows the algorithm to be more efficient in
  276. the long term at the cost of some local increase in the value of the
  277. objective function.
  278. This is because allowing for non-decreasing objective function values
  279. in a principled manner allows the algorithm to *jump over boulders* as
  280. the method is not restricted to move into narrow valleys while
  281. preserving its convergence properties.
  282. Setting :member:`Solver::Options::use_nonmonotonic_steps` to ``true``
  283. enables the non-monotonic trust region algorithm as described by Conn,
  284. Gould & Toint in [Conn]_.
  285. Even though the value of the objective function may be larger than the
  286. minimum value encountered over the course of the optimization, the
  287. final parameters returned to the user are the ones corresponding to
  288. the minimum cost over all iterations.
  289. The option to take non-monotonic steps is available for all trust
  290. region strategies.
  291. .. _section-line-search-methods:
  292. Line Search Methods
  293. ===================
  294. .. NOTE::
  295. The line search method in Ceres Solver cannot handle bounds
  296. constraints right now, so it can only be used for solving
  297. unconstrained problems.
  298. The basic line search algorithm looks something like this:
  299. 1. Given an initial point :math:`x`
  300. 2. :math:`\Delta x = -H^{-1}(x) g(x)`
  301. 3. :math:`\arg \min_\mu \frac{1}{2} \| F(x + \mu \Delta x) \|^2`
  302. 4. :math:`x = x + \mu \Delta x`
  303. 5. Goto 2.
  304. Here :math:`H(x)` is some approximation to the Hessian of the
  305. objective function, and :math:`g(x)` is the gradient at
  306. :math:`x`. Depending on the choice of :math:`H(x)` we get a variety of
  307. different search directions :math:`\Delta x`.
  308. Step 4, which is a one dimensional optimization or `Line Search` along
  309. :math:`\Delta x` is what gives this class of methods its name.
  310. Different line search algorithms differ in their choice of the search
  311. direction :math:`\Delta x` and the method used for one dimensional
  312. optimization along :math:`\Delta x`. The choice of :math:`H(x)` is the
  313. primary source of computational complexity in these
  314. methods. Currently, Ceres Solver supports four choices of search
  315. directions, all aimed at large scale problems.
  316. 1. ``STEEPEST_DESCENT`` This corresponds to choosing :math:`H(x)` to
  317. be the identity matrix. This is not a good search direction for
  318. anything but the simplest of the problems. It is only included here
  319. for completeness.
  320. 2. ``NONLINEAR_CONJUGATE_GRADIENT`` A generalization of the Conjugate
  321. Gradient method to non-linear functions. The generalization can be
  322. performed in a number of different ways, resulting in a variety of
  323. search directions. Ceres Solver currently supports
  324. ``FLETCHER_REEVES``, ``POLAK_RIBIERE`` and ``HESTENES_STIEFEL``
  325. directions.
  326. 3. ``BFGS`` A generalization of the Secant method to multiple
  327. dimensions in which a full, dense approximation to the inverse
  328. Hessian is maintained and used to compute a quasi-Newton step
  329. [NocedalWright]_. ``BFGS`` and its limited memory variant ``LBFGS``
  330. are currently the best known general quasi-Newton algorithm.
  331. 4. ``LBFGS`` A limited memory approximation to the full ``BFGS``
  332. method in which the last `M` iterations are used to approximate the
  333. inverse Hessian used to compute a quasi-Newton step [Nocedal]_,
  334. [ByrdNocedal]_.
  335. Currently Ceres Solver supports both a backtracking and interpolation
  336. based `Armijo line search algorithm
  337. <https://en.wikipedia.org/wiki/Backtracking_line_search>`_ (``ARMIJO``)
  338. , and a sectioning / zoom interpolation (strong) `Wolfe condition line
  339. search algorithm <https://en.wikipedia.org/wiki/Wolfe_conditions>`_
  340. (``WOLFE``).
  341. .. NOTE::
  342. In order for the assumptions underlying the ``BFGS`` and ``LBFGS``
  343. methods to be satisfied the ``WOLFE`` algorithm must be used.
  344. .. _section-linear-solver:
  345. Linear Solvers
  346. ==============
  347. Observe that for both of the trust-region methods described above, the
  348. key computational cost is the solution of a linear least squares
  349. problem of the form
  350. .. math:: \min_{\Delta x} \frac{1}{2} \|J(x)\Delta x + F(x)\|^2 .
  351. :label: simple2
  352. Let :math:`H(x)= J(x)^\top J(x)` and :math:`g(x) = -J(x)^\top
  353. F(x)`. For notational convenience let us also drop the dependence on
  354. :math:`x`. Then it is easy to see that solving :eq:`simple2` is
  355. equivalent to solving the *normal equations*.
  356. .. math:: H \Delta x = g
  357. :label: normal
  358. Ceres provides a number of different options for solving :eq:`normal`.
  359. .. _section-qr:
  360. DENSE_QR
  361. --------
  362. For small problems (a couple of hundred parameters and a few thousand
  363. residuals) with relatively dense Jacobians, QR-decomposition is the
  364. method of choice [Bjorck]_. Let :math:`J = QR` be the QR-decomposition
  365. of :math:`J`, where :math:`Q` is an orthonormal matrix and :math:`R`
  366. is an upper triangular matrix [TrefethenBau]_. Then it can be shown
  367. that the solution to :eq:`normal` is given by
  368. .. math:: \Delta x^* = -R^{-1}Q^\top f
  369. You can use QR-decomposition by setting
  370. :member:`Solver::Options::linear_solver_type` to ``DENSE_QR``.
  371. By default (``Solver::Options::dense_linear_algebra_library_type =
  372. EIGEN``) Ceres Solver will use `Eigen Householder QR factorization
  373. <https://eigen.tuxfamily.org/dox-devel/classEigen_1_1HouseholderQR.html>`_
  374. .
  375. If Ceres Solver has been built with an optimized LAPACK
  376. implementation, then the user can also choose to use LAPACK's
  377. `DGEQRF`_ routine by setting
  378. :member:`Solver::Options::dense_linear_algebra_library_type` to
  379. ``LAPACK``. Depending on the `LAPACK` and the underlying `BLAS`
  380. implementation this may perform better than using Eigen's Householder
  381. QR factorization.
  382. .. _DGEQRF: https://netlib.org/lapack/explore-html/df/dc5/group__variants_g_ecomputational_ga3766ea903391b5cf9008132f7440ec7b.html
  383. If an NVIDIA GPU is available and Ceres Solver has been built with
  384. CUDA support enabled, then the user can also choose to perform the
  385. QR-decomposition on the GPU by setting
  386. :member:`Solver::Options::dense_linear_algebra_library_type` to
  387. ``CUDA``. Depending on the GPU this can lead to a substantial
  388. speedup. Using CUDA only makes sense for moderate to large sized
  389. problems. This is because to perform the decomposition on the GPU the
  390. matrix :math:`J` needs to be transferred from the CPU to the GPU and
  391. this incurs a cost. So unless the speedup from doing the decomposition
  392. on the GPU is large enough to also account for the time taken to
  393. transfer the Jacobian to the GPU, using CUDA will not be better than
  394. just doing the decomposition on the CPU.
  395. .. _section-dense-normal-cholesky:
  396. DENSE_NORMAL_CHOLESKY
  397. ---------------------
  398. It is often the case that the number of rows in the Jacobian :math:`J`
  399. are much larger than the the number of columns. The complexity of QR
  400. factorization scales linearly with the number of rows, so beyond a
  401. certain size it is more efficient to solve :eq:`normal` using a dense
  402. `Cholesky factorization
  403. <https://en.wikipedia.org/wiki/Cholesky_decomposition>`_.
  404. Let :math:`H = R^\top R` be the Cholesky factorization of the normal
  405. equations, where :math:`R` is an upper triangular matrix, then the
  406. solution to :eq:`normal` is given by
  407. .. math::
  408. \Delta x^* = R^{-1} R^{-\top} g.
  409. The observant reader will note that the :math:`R` in the Cholesky
  410. factorization of :math:`H` is the same upper triangular matrix
  411. :math:`R` in the QR factorization of :math:`J`. Since :math:`Q` is an
  412. orthonormal matrix, :math:`J=QR` implies that :math:`J^\top J = R^\top
  413. Q^\top Q R = R^\top R`.
  414. Unfortunately, forming the matrix :math:`H = J'J` squares the
  415. condition number. As a result while the cost of forming :math:`H` and
  416. computing its Cholesky factorization is lower than computing the
  417. QR-factorization of :math:`J`, we pay the price in terms of increased
  418. numerical instability and potential failure of the Cholesky
  419. factorization for ill-conditioned Jacobians.
  420. You can use dense Cholesky factorization by setting
  421. :member:`Solver::Options::linear_solver_type` to
  422. ``DENSE_NORMAL_CHOLESKY``.
  423. By default (``Solver::Options::dense_linear_algebra_library_type =
  424. EIGEN``) Ceres Solver will use `Eigen's LLT factorization`_ routine.
  425. .. _Eigen's LLT Factorization: https://eigen.tuxfamily.org/dox/classEigen_1_1LLT.html
  426. If Ceres Solver has been built with an optimized LAPACK
  427. implementation, then the user can also choose to use LAPACK's
  428. `DPOTRF`_ routine by setting
  429. :member:`Solver::Options::dense_linear_algebra_library_type` to
  430. ``LAPACK``. Depending on the `LAPACK` and the underlying `BLAS`
  431. implementation this may perform better than using Eigen's Cholesky
  432. factorization.
  433. .. _DPOTRF: https://www.netlib.org/lapack/explore-html/d1/d7a/group__double_p_ocomputational_ga2f55f604a6003d03b5cd4a0adcfb74d6.html
  434. If an NVIDIA GPU is available and Ceres Solver has been built with
  435. CUDA support enabled, then the user can also choose to perform the
  436. Cholesky factorization on the GPU by setting
  437. :member:`Solver::Options::dense_linear_algebra_library_type` to
  438. ``CUDA``. Depending on the GPU this can lead to a substantial speedup.
  439. Using CUDA only makes sense for moderate to large sized problems. This
  440. is because to perform the decomposition on the GPU the matrix
  441. :math:`H` needs to be transferred from the CPU to the GPU and this
  442. incurs a cost. So unless the speedup from doing the decomposition on
  443. the GPU is large enough to also account for the time taken to transfer
  444. the Jacobian to the GPU, using CUDA will not be better than just doing
  445. the decomposition on the CPU.
  446. .. _section-sparse-normal-cholesky:
  447. SPARSE_NORMAL_CHOLESKY
  448. ----------------------
  449. Large non-linear least square problems are usually sparse. In such
  450. cases, using a dense QR or Cholesky factorization is inefficient. For
  451. such problems, Cholesky factorization routines which treat :math:`H`
  452. as a sparse matrix and computes a sparse factor :math:`R` are better
  453. suited [Davis]_. This can lead to substantial savings in memory and
  454. CPU time for large sparse problems.
  455. You can use dense Cholesky factorization by setting
  456. :member:`Solver::Options::linear_solver_type` to
  457. ``SPARSE_NORMAL_CHOLESKY``.
  458. The use of this linear solver requires that Ceres is compiled with
  459. support for at least one of:
  460. 1. `SuiteSparse <https://people.engr.tamu.edu/davis/suitesparse.html>`_ (``SUITE_SPARSE``).
  461. 2. `Apple's Accelerate framework
  462. <https://developer.apple.com/documentation/accelerate/sparse_solvers?language=objc>`_
  463. (``ACCELERATE_SPARSE``).
  464. 3. `Eigen's sparse linear solvers
  465. <https://eigen.tuxfamily.org/dox/group__SparseCholesky__Module.html>`_
  466. (``EIGEN_SPARSE``).
  467. SuiteSparse and Accelerate offer high performance sparse Cholesky
  468. factorization routines as they level-3 BLAS routines
  469. internally. Eigen's sparse Cholesky routines are *simplicial* and do
  470. not use dense linear algebra routines and as a result cannot compete
  471. with SuiteSparse and Accelerate, especially on large problems. As a
  472. result to get the best performance out of SuiteSparse it should be
  473. linked to high quality BLAS and LAPACK implementations e.g. `ATLAS
  474. <https://math-atlas.sourceforge.net/>`_, `OpenBLAS
  475. <https://www.openblas.net/>`_ or `Intel MKL
  476. <https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html>`_.
  477. A critical part of a sparse Cholesky factorization routine is the use
  478. a fill-reducing ordering. By default Ceres Solver uses the Approximate
  479. Minimum Degree (``AMD``) ordering, which usually performs well, but
  480. there are other options that may perform better depending on the
  481. actual sparsity structure of the Jacobian. See :ref:`section-ordering`
  482. for more details.
  483. .. _section-cgnr:
  484. CGNR
  485. ----
  486. For general sparse problems, if the problem is too large for sparse
  487. Cholesky factorization or a sparse linear algebra library is not
  488. linked into Ceres, another option is the ``CGNR`` solver. This solver
  489. uses the `Conjugate Gradients
  490. <https://en.wikipedia.org/wiki/Conjugate_gradient_method>_` method on
  491. the *normal equations*, but without forming the normal equations
  492. explicitly. It exploits the relation
  493. .. math::
  494. H x = J^\top J x = J^\top(J x)
  495. Because ``CGNR`` never solves the linear system exactly, when the user
  496. chooses ``CGNR`` as the linear solver, Ceres automatically switches
  497. from the exact step algorithm to an inexact step algorithm. This also
  498. means that ``CGNR`` can only be used with ``LEVENBERG_MARQUARDT`` and
  499. not with ``DOGLEG`` trust region strategy.
  500. ``CGNR`` by default runs on the CPU. However, if an NVIDIA GPU is
  501. available and Ceres Solver has been built with CUDA support enabled,
  502. then the user can also choose to run ``CGNR`` on the GPU by setting
  503. :member:`Solver::Options::sparse_linear_algebra_library_type` to
  504. ``CUDA_SPARSE``. The key complexity of ``CGNR`` comes from evaluating
  505. the two sparse-matrix vector products (SpMV) :math:`Jx` and
  506. :math:`J'y`. GPUs are particularly well suited for doing sparse
  507. matrix-vector products. As a result, for large problems using a GPU
  508. can lead to a substantial speedup.
  509. The convergence of Conjugate Gradients depends on the conditioner
  510. number :math:`\kappa(H)`. Usually :math:`H` is quite poorly
  511. conditioned and a `Preconditioner
  512. <https://en.wikipedia.org/wiki/Preconditioner>`_ must be used to get
  513. reasonable performance. See section on :ref:`section-preconditioner`
  514. for more details.
  515. .. _section-schur:
  516. DENSE_SCHUR & SPARSE_SCHUR
  517. --------------------------
  518. While it is possible to use ``SPARSE_NORMAL_CHOLESKY`` to solve bundle
  519. adjustment problems, they have a special sparsity structure that can
  520. be exploited to solve the normal equations more efficiently.
  521. Suppose that the bundle adjustment problem consists of :math:`p`
  522. cameras and :math:`q` points and the variable vector :math:`x` has the
  523. block structure :math:`x = [y_{1}, ... ,y_{p},z_{1},
  524. ... ,z_{q}]`. Where, :math:`y` and :math:`z` correspond to camera and
  525. point parameters respectively. Further, let the camera blocks be of
  526. size :math:`c` and the point blocks be of size :math:`s` (for most
  527. problems :math:`c` = :math:`6`--`9` and :math:`s = 3`). Ceres does not
  528. impose any constancy requirement on these block sizes, but choosing
  529. them to be constant simplifies the exposition.
  530. The key property of bundle adjustment problems which we will exploit
  531. is the fact that no term :math:`f_{i}` in :eq:`nonlinsq` includes two
  532. or more point blocks at the same time. This in turn implies that the
  533. matrix :math:`H` is of the form
  534. .. math:: H = \left[ \begin{matrix} B & E\\ E^\top & C \end{matrix} \right]\ ,
  535. :label: hblock
  536. where :math:`B \in \mathbb{R}^{pc\times pc}` is a block sparse matrix
  537. with :math:`p` blocks of size :math:`c\times c` and :math:`C \in
  538. \mathbb{R}^{qs\times qs}` is a block diagonal matrix with :math:`q` blocks
  539. of size :math:`s\times s`. :math:`E \in \mathbb{R}^{pc\times qs}` is a
  540. general block sparse matrix, with a block of size :math:`c\times s`
  541. for each observation. Let us now block partition :math:`\Delta x =
  542. [\Delta y,\Delta z]` and :math:`g=[v,w]` to restate :eq:`normal`
  543. as the block structured linear system
  544. .. math:: \left[ \begin{matrix} B & E\\ E^\top & C \end{matrix}
  545. \right]\left[ \begin{matrix} \Delta y \\ \Delta z
  546. \end{matrix} \right] = \left[ \begin{matrix} v\\ w
  547. \end{matrix} \right]\ ,
  548. :label: linear2
  549. and apply Gaussian elimination to it. As we noted above, :math:`C` is
  550. a block diagonal matrix, with small diagonal blocks of size
  551. :math:`s\times s`. Thus, calculating the inverse of :math:`C` by
  552. inverting each of these blocks is cheap. This allows us to eliminate
  553. :math:`\Delta z` by observing that :math:`\Delta z = C^{-1}(w - E^\top
  554. \Delta y)`, giving us
  555. .. math:: \left[B - EC^{-1}E^\top\right] \Delta y = v - EC^{-1}w\ .
  556. :label: schur
  557. The matrix
  558. .. math:: S = B - EC^{-1}E^\top
  559. is the Schur complement of :math:`C` in :math:`H`. It is also known as
  560. the *reduced camera matrix*, because the only variables
  561. participating in :eq:`schur` are the ones corresponding to the
  562. cameras. :math:`S \in \mathbb{R}^{pc\times pc}` is a block structured
  563. symmetric positive definite matrix, with blocks of size :math:`c\times
  564. c`. The block :math:`S_{ij}` corresponding to the pair of images
  565. :math:`i` and :math:`j` is non-zero if and only if the two images
  566. observe at least one common point.
  567. Now :eq:`linear2` can be solved by first forming :math:`S`, solving
  568. for :math:`\Delta y`, and then back-substituting :math:`\Delta y` to
  569. obtain the value of :math:`\Delta z`. Thus, the solution of what was
  570. an :math:`n\times n`, :math:`n=pc+qs` linear system is reduced to the
  571. inversion of the block diagonal matrix :math:`C`, a few matrix-matrix
  572. and matrix-vector multiplies, and the solution of block sparse
  573. :math:`pc\times pc` linear system :eq:`schur`. For almost all
  574. problems, the number of cameras is much smaller than the number of
  575. points, :math:`p \ll q`, thus solving :eq:`schur` is significantly
  576. cheaper than solving :eq:`linear2`. This is the *Schur complement
  577. trick* [Brown]_.
  578. This still leaves open the question of solving :eq:`schur`. As we
  579. discussed when considering the exact solution of the normal equations
  580. using Cholesky factorization, we have two options.
  581. 1. ``DENSE_SCHUR`` - The first is **dense Cholesky factorization**,
  582. where we store and factor :math:`S` as a dense matrix. This method has
  583. :math:`O(p^2)` space complexity and :math:`O(p^3)` time complexity and
  584. is only practical for problems with up to a few hundred cameras.
  585. 2. ``SPARSE_SCHUR`` - For large bundle adjustment problems :math:`S`
  586. is typically a fairly sparse matrix, as most images only see a small
  587. fraction of the scene. This leads us to the second option: **sparse
  588. Cholesky factorization** [Davis]_. Here we store :math:`S` as a
  589. sparse matrix, use row and column re-ordering algorithms to maximize
  590. the sparsity of the Cholesky decomposition, and focus their compute
  591. effort on the non-zero part of the factorization [Davis]_ [Chen]_
  592. . Sparse direct methods, depending on the exact sparsity structure of
  593. the Schur complement, allow bundle adjustment algorithms to scenes
  594. with thousands of cameras.
  595. .. _section-iterative_schur:
  596. ITERATIVE_SCHUR
  597. ---------------
  598. Another option for bundle adjustment problems is to apply Conjugate
  599. Gradients to the reduced camera matrix :math:`S` instead of
  600. :math:`H`. One reason to do this is that :math:`S` is a much smaller
  601. matrix than :math:`H`, but more importantly, it can be shown that
  602. :math:`\kappa(S)\leq \kappa(H)` [Agarwal]_.
  603. Ceres implements Conjugate Gradients on :math:`S` as the
  604. ``ITERATIVE_SCHUR`` solver. When the user chooses ``ITERATIVE_SCHUR``
  605. as the linear solver, Ceres automatically switches from the exact step
  606. algorithm to an inexact step algorithm.
  607. The key computational operation when using Conjuagate Gradients is the
  608. evaluation of the matrix vector product :math:`Sx` for an arbitrary
  609. vector :math:`x`. Because PCG only needs access to :math:`S` via its
  610. product with a vector, one way to evaluate :math:`Sx` is to observe
  611. that
  612. .. math:: x_1 &= E^\top x\\
  613. x_2 &= C^{-1} x_1\\
  614. x_3 &= Ex_2\\
  615. x_4 &= Bx\\
  616. Sx &= x_4 - x_3
  617. :label: schurtrick1
  618. Thus, we can run Conjugate Gradients on :math:`S` with the same
  619. computational effort per iteration as Conjugate Gradients on
  620. :math:`H`, while reaping the benefits of a more powerful
  621. preconditioner. In fact, we do not even need to compute :math:`H`,
  622. :eq:`schurtrick1` can be implemented using just the columns of
  623. :math:`J`.
  624. Equation :eq:`schurtrick1` is closely related to *Domain Decomposition
  625. methods* for solving large linear systems that arise in structural
  626. engineering and partial differential equations. In the language of
  627. Domain Decomposition, each point in a bundle adjustment problem is a
  628. domain, and the cameras form the interface between these domains. The
  629. iterative solution of the Schur complement then falls within the
  630. sub-category of techniques known as Iterative Sub-structuring [Saad]_
  631. [Mathew]_.
  632. While in most cases the above method for evaluating :math:`Sx` is the
  633. way to go, for some problems it is better to compute the Schur
  634. complemenent :math:`S` explicitly and then run Conjugate Gradients on
  635. it. This can be done by settin
  636. ``Solver::Options::use_explicit_schur_complement`` to ``true``. This
  637. option can only be used with the ``SCHUR_JACOBI`` preconditioner.
  638. .. _section-schur_power_series_expansion:
  639. SCHUR_POWER_SERIES_EXPANSION
  640. ----------------------------
  641. It can be shown that the inverse of the Schur complement can be
  642. written as an infinite power-series [Weber]_ [Zheng]_:
  643. .. math:: S &= B - EC^{-1}E^\top\\
  644. &= B(I - B^{-1}EC^{-1}E^\top)\\
  645. S^{-1} &= (I - B^{-1}EC^{-1}E^\top)^{-1} B^{-1}\\
  646. & = \sum_{i=0}^\infty \left(B^{-1}EC^{-1}E^\top\right)^{i} B^{-1}
  647. As a result a truncated version of this power series expansion can be
  648. used to approximate the inverse and therefore the solution to
  649. :eq:`schur`. Ceres allows the user to use Schur power series expansion
  650. in three ways.
  651. 1. As a linear solver. This is what [Weber]_ calls **Power Bundle
  652. Adjustment** and corresponds to using the truncated power series to
  653. approximate the inverse of the Schur complement. This is done by
  654. setting the following options.
  655. .. code-block:: c++
  656. Solver::Options::linear_solver_type = ITERATIVE_SCHUR
  657. Solver::Options::preconditioner_type = IDENTITY
  658. Solver::Options::use_spse_initialization = true
  659. Solver::Options::max_linear_solver_iterations = 0;
  660. // The following two settings are worth tuning for your application.
  661. Solver::Options::max_num_spse_iterations = 5;
  662. Solver::Options::spse_tolerance = 0.1;
  663. 2. As a preconditioner for ``ITERATIVE_SCHUR``. Any method for
  664. approximating the inverse of a matrix can also be used as a
  665. preconditioner. This is enabled by setting the following options.
  666. .. code-block:: c++
  667. Solver::Options::linear_solver_type = ITERATIVE_SCHUR
  668. Solver::Options::preconditioner_type = SCHUR_POWER_SERIES_EXPANSION;
  669. Solver::Options::use_spse_initialization = false;
  670. // This is worth tuning for your application.
  671. Solver::Options::max_num_spse_iterations = 5;
  672. 3. As initialization for ``ITERATIIVE_SCHUR`` with any
  673. preconditioner. This is a combination of the above two, where the
  674. Schur Power Series Expansion
  675. .. code-block:: c++
  676. Solver::Options::linear_solver_type = ITERATIVE_SCHUR
  677. Solver::Options::preconditioner_type = ... // Preconditioner of your choice.
  678. Solver::Options::use_spse_initialization = true
  679. Solver::Options::max_linear_solver_iterations = 0;
  680. // The following two settings are worth tuning for your application.
  681. Solver::Options::max_num_spse_iterations = 5;
  682. // This only affects the initialization but not the preconditioner.
  683. Solver::Options::spse_tolerance = 0.1;
  684. .. _section-mixed-precision:
  685. Mixed Precision Solves
  686. ======================
  687. Generally speaking Ceres Solver does all its arithmetic in double
  688. precision. Sometimes though, one can use single precision arithmetic
  689. to get substantial speedups. Currently, for linear solvers that
  690. perform Cholesky factorization (sparse or dense) the user has the
  691. option cast the linear system to single precision and then use
  692. single precision Cholesky factorization routines to solve the
  693. resulting linear system. This can be enabled by setting
  694. :member:`Solver::Options::use_mixed_precision_solves` to ``true``.
  695. Depending on the conditioning of the problem, the use of single
  696. precision factorization may lead to some loss of accuracy. Some of
  697. this accuracy can be recovered by performing `Iterative Refinement
  698. <https://en.wikipedia.org/wiki/Iterative_refinement>`_. The number of
  699. iterations of iterative refinement are controlled by
  700. :member:`Solver::Options::max_num_refinement_iterations`. The default
  701. value of this parameter is zero, which means if
  702. :member:`Solver::Options::use_mixed_precision_solves` is ``true``,
  703. then no iterative refinement is performed. Usually 2-3 refinement
  704. iterations are enough.
  705. Mixed precision solves are available in the following linear solver
  706. configurations:
  707. 1. ``DENSE_NORMAL_CHOLESKY`` + ``EIGEN``/ ``LAPACK`` / ``CUDA``.
  708. 2. ``DENSE_SCHUR`` + ``EIGEN``/ ``LAPACK`` / ``CUDA``.
  709. 3. ``SPARSE_NORMAL_CHOLESKY`` + ``EIGEN_SPARSE`` / ``ACCELERATE_SPARSE``
  710. 4. ``SPARSE_SCHUR`` + ``EIGEN_SPARSE`` / ``ACCELERATE_SPARSE``
  711. Mixed precision solves area not available when using ``SUITE_SPARSE``
  712. as the sparse linear algebra backend because SuiteSparse/CHOLMOD does
  713. not support single precision solves.
  714. .. _section-preconditioner:
  715. Preconditioners
  716. ===============
  717. The convergence rate of Conjugate Gradients for solving :eq:`normal`
  718. depends on the distribution of eigenvalues of :math:`H` [Saad]_. A
  719. useful upper bound is :math:`\sqrt{\kappa(H)}`, where,
  720. :math:`\kappa(H)` is the condition number of the matrix :math:`H`. For
  721. most non-linear least squares problems, :math:`\kappa(H)` is high and
  722. a direct application of Conjugate Gradients to :eq:`normal` results in
  723. extremely poor performance.
  724. The solution to this problem is to replace :eq:`normal` with a
  725. *preconditioned* system. Given a linear system, :math:`Ax =b` and a
  726. preconditioner :math:`M` the preconditioned system is given by
  727. :math:`M^{-1}Ax = M^{-1}b`. The resulting algorithm is known as
  728. Preconditioned Conjugate Gradients algorithm (PCG) and its worst case
  729. complexity now depends on the condition number of the *preconditioned*
  730. matrix :math:`\kappa(M^{-1}A)`.
  731. The computational cost of using a preconditioner :math:`M` is the cost
  732. of computing :math:`M` and evaluating the product :math:`M^{-1}y` for
  733. arbitrary vectors :math:`y`. Thus, there are two competing factors to
  734. consider: How much of :math:`H`'s structure is captured by :math:`M`
  735. so that the condition number :math:`\kappa(HM^{-1})` is low, and the
  736. computational cost of constructing and using :math:`M`. The ideal
  737. preconditioner would be one for which :math:`\kappa(M^{-1}A)
  738. =1`. :math:`M=A` achieves this, but it is not a practical choice, as
  739. applying this preconditioner would require solving a linear system
  740. equivalent to the unpreconditioned problem. It is usually the case
  741. that the more information :math:`M` has about :math:`H`, the more
  742. expensive it is use. For example, Incomplete Cholesky factorization
  743. based preconditioners have much better convergence behavior than the
  744. Jacobi preconditioner, but are also much more expensive.
  745. For a survey of the state of the art in preconditioning linear least
  746. squares problems with general sparsity structure see [GouldScott]_.
  747. Ceres Solver comes with an number of preconditioners suited for
  748. problems with general sparsity as well as the special sparsity
  749. structure encountered in bundle adjustment problems.
  750. IDENTITY
  751. --------
  752. This is equivalent to using an identity matrix as a preconditioner,
  753. i.e. no preconditioner at all.
  754. JACOBI
  755. ------
  756. The simplest of all preconditioners is the diagonal or Jacobi
  757. preconditioner, i.e., :math:`M=\operatorname{diag}(A)`, which for
  758. block structured matrices like :math:`H` can be generalized to the
  759. block Jacobi preconditioner. The ``JACOBI`` preconditioner in Ceres
  760. when used with :ref:`section-cgnr` refers to the block diagonal of
  761. :math:`H` and when used with :ref:`section-iterative_schur` refers to
  762. the block diagonal of :math:`B` [Mandel]_.
  763. For detailed performance data about the performance of ``JACOBI`` on
  764. bundle adjustment problems see [Agarwal]_.
  765. SCHUR_JACOBI
  766. ------------
  767. Another obvious choice for :ref:`section-iterative_schur` is the block
  768. diagonal of the Schur complement matrix :math:`S`, i.e, the block
  769. Jacobi preconditioner for :math:`S`. In Ceres we refer to it as the
  770. ``SCHUR_JACOBI`` preconditioner.
  771. For detailed performance data about the performance of
  772. ``SCHUR_JACOBI`` on bundle adjustment problems see [Agarwal]_.
  773. CLUSTER_JACOBI and CLUSTER_TRIDIAGONAL
  774. --------------------------------------
  775. For bundle adjustment problems arising in reconstruction from
  776. community photo collections, more effective preconditioners can be
  777. constructed by analyzing and exploiting the camera-point visibility
  778. structure of the scene.
  779. The key idea is to cluster the cameras based on the visibility
  780. structure of the scene. The similarity between a pair of cameras
  781. :math:`i` and :math:`j` is given by:
  782. .. math:: S_{ij} = \frac{|V_i \cap V_j|}{|V_i| |V_j|}
  783. Here :math:`V_i` is the set of scene points visible in camera
  784. :math:`i`. This idea was first exploited by [KushalAgarwal]_ to create
  785. the ``CLUSTER_JACOBI`` and the ``CLUSTER_TRIDIAGONAL`` preconditioners
  786. which Ceres implements.
  787. The performance of these two preconditioners depends on the speed and
  788. clustering quality of the clustering algorithm used when building the
  789. preconditioner. In the original paper, [KushalAgarwal]_ used the
  790. Canonical Views algorithm [Simon]_, which while producing high quality
  791. clusterings can be quite expensive for large graphs. So, Ceres
  792. supports two visibility clustering algorithms - ``CANONICAL_VIEWS``
  793. and ``SINGLE_LINKAGE``. The former is as the name implies Canonical
  794. Views algorithm of [Simon]_. The latter is the the classic `Single
  795. Linkage Clustering
  796. <https://en.wikipedia.org/wiki/Single-linkage_clustering>`_
  797. algorithm. The choice of clustering algorithm is controlled by
  798. :member:`Solver::Options::visibility_clustering_type`.
  799. SCHUR_POWER_SERIES_EXPANSION
  800. ----------------------------
  801. As explained in :ref:`section-schur_power_series_expansion`, the Schur
  802. complement matrix admits a power series expansion and a truncated
  803. version of this power series can be used as a preconditioner for
  804. ``ITERATIVE_SCHUR``. When used as a preconditioner
  805. :member:`Solver::Options::max_num_spse_iterations` controls the number
  806. of terms in the power series that are used.
  807. SUBSET
  808. ------
  809. This is a preconditioner for problems with general sparsity. Given a
  810. subset of residual blocks of a problem, it uses the corresponding
  811. subset of the rows of the Jacobian to construct a preconditioner
  812. [Dellaert]_.
  813. Suppose the Jacobian :math:`J` has been horizontally partitioned as
  814. .. math:: J = \begin{bmatrix} P \\ Q \end{bmatrix}
  815. Where, :math:`Q` is the set of rows corresponding to the residual
  816. blocks in
  817. :member:`Solver::Options::residual_blocks_for_subset_preconditioner`. The
  818. preconditioner is the matrix :math:`(Q^\top Q)^{-1}`.
  819. The efficacy of the preconditioner depends on how well the matrix
  820. :math:`Q` approximates :math:`J^\top J`, or how well the chosen
  821. residual blocks approximate the full problem.
  822. This preconditioner is NOT available when running ``CGNR`` using
  823. ``CUDA``.
  824. .. _section-ordering:
  825. Ordering
  826. ========
  827. The order in which variables are eliminated in a linear solver can
  828. have a significant of impact on the efficiency and accuracy of the
  829. method. For example when doing sparse Cholesky factorization, there
  830. are matrices for which a good ordering will give a Cholesky factor
  831. with :math:`O(n)` storage, whereas a bad ordering will result in an
  832. completely dense factor.
  833. Ceres allows the user to provide varying amounts of hints to the
  834. solver about the variable elimination ordering to use. This can range
  835. from no hints, where the solver is free to decide the best possible
  836. ordering based on the user's choices like the linear solver being
  837. used, to an exact order in which the variables should be eliminated,
  838. and a variety of possibilities in between.
  839. The simplest thing to do is to just set
  840. :member:`Solver::Options::linear_solver_ordering_type` to ``AMD``
  841. (default) or ``NESDIS`` based on your understanding of the problem or
  842. empirical testing.
  843. More information can be commmuniucated by using an instance
  844. :class:`ParameterBlockOrdering` class.
  845. Formally an ordering is an ordered partitioning of the
  846. parameter blocks, i.e, each parameter block belongs to exactly
  847. one group, and each group has a unique non-negative integer
  848. associated with it, that determines its order in the set of
  849. groups.
  850. e.g. Consider the linear system
  851. .. math::
  852. x + y &= 3 \\
  853. 2x + 3y &= 7
  854. There are two ways in which it can be solved. First eliminating
  855. :math:`x` from the two equations, solving for :math:`y` and then back
  856. substituting for :math:`x`, or first eliminating :math:`y`, solving
  857. for :math:`x` and back substituting for :math:`y`. The user can
  858. construct three orderings here.
  859. 1. :math:`\{0: x\}, \{1: y\}` - eliminate :math:`x` first.
  860. 2. :math:`\{0: y\}, \{1: x\}` - eliminate :math:`y` first.
  861. 3. :math:`\{0: x, y\}` - Solver gets to decide the elimination order.
  862. Thus, to have Ceres determine the ordering automatically, put all the
  863. variables in group 0 and to control the ordering for every variable,
  864. create groups :math:`0 \dots N-1`, one per variable, in the desired
  865. order.
  866. ``linear_solver_ordering == nullptr`` and an ordering where all the
  867. parameter blocks are in one elimination group mean the same thing -
  868. the solver is free to choose what it thinks is the best elimination
  869. ordering using the ordering algorithm (specified using
  870. :member:`Solver::Options::linear_solver_ordering_type`). Therefore in
  871. the following we will only consider the case where
  872. ``linear_solver_ordering != nullptr``.
  873. The exact interpretation of the ``linear_solver_ordering`` depends on
  874. the values of :member:`Solver::Options::linear_solver_ordering_type`,
  875. :member:`Solver::Options::linear_solver_type`,
  876. :member:`Solver::Options::preconditioner_type` and
  877. :member:`Solver::Options::sparse_linear_algebra_library_type` as we will
  878. explain below.
  879. Bundle Adjustment
  880. -----------------
  881. If the user is using one of the Schur solvers (``DENSE_SCHUR``,
  882. ``SPARSE_SCHUR``, ``ITERATIVE_SCHUR``) and chooses to specify an
  883. ordering, it must have one important property. The lowest numbered
  884. elimination group must form an independent set in the graph
  885. corresponding to the Hessian, or in other words, no two parameter
  886. blocks in the first elimination group should co-occur in the same
  887. residual block. For the best performance, this elimination group
  888. should be as large as possible. For standard bundle adjustment
  889. problems, this corresponds to the first elimination group containing
  890. all the 3d points, and the second containing the parameter blocks for
  891. all the cameras.
  892. If the user leaves the choice to Ceres, then the solver uses an
  893. approximate maximum independent set algorithm to identify the first
  894. elimination group [LiSaad]_.
  895. ``sparse_linear_algebra_library_type = SUITE_SPARSE``
  896. -----------------------------------------------------
  897. **linear_solver_ordering_type = AMD**
  898. A constrained Approximate Minimum Degree (CAMD) ordering is used where
  899. the parameter blocks in the lowest numbered group are eliminated
  900. first, and then the parameter blocks in the next lowest numbered group
  901. and so on. Within each group, CAMD is free to order the parameter blocks
  902. as it chooses.
  903. **linear_solver_ordering_type = NESDIS**
  904. a. ``linear_solver_type = SPARSE_NORMAL_CHOLESKY`` or
  905. ``linear_solver_type = CGNR`` and ``preconditioner_type = SUBSET``
  906. The value of ``linear_solver_ordering`` is ignored and a Nested
  907. Dissection algorithm is used to compute a fill reducing ordering.
  908. b. ``linear_solver_type = SPARSE_SCHUR/DENSE_SCHUR/ITERATIVE_SCHUR``
  909. ONLY the lowest group are used to compute the Schur complement, and
  910. Nested Dissection is used to compute a fill reducing ordering for
  911. the Schur Complement (or its preconditioner).
  912. ``sparse_linear_algebra_library_type = EIGEN_SPARSE/ACCELERATE_SPARSE``
  913. -----------------------------------------------------------------------
  914. a. ``linear_solver_type = SPARSE_NORMAL_CHOLESKY`` or
  915. ``linear_solver_type = CGNR`` and ``preconditioner_type = SUBSET``
  916. The value of ``linear_solver_ordering`` is ignored and ``AMD`` or
  917. ``NESDIS`` is used to compute a fill reducing ordering as requested
  918. by the user.
  919. b. ``linear_solver_type = SPARSE_SCHUR/DENSE_SCHUR/ITERATIVE_SCHUR``
  920. ONLY the lowest group are used to compute the Schur complement, and
  921. ``AMD`` or ``NESID`` is used to compute a fill reducing ordering
  922. for the Schur Complement (or its preconditioner) as requested by
  923. the user.
  924. .. _section-solver-options:
  925. :class:`Solver::Options`
  926. ========================
  927. .. class:: Solver::Options
  928. :class:`Solver::Options` controls the overall behavior of the
  929. solver. We list the various settings and their default values below.
  930. .. function:: bool Solver::Options::IsValid(std::string* error) const
  931. Validate the values in the options struct and returns true on
  932. success. If there is a problem, the method returns false with
  933. ``error`` containing a textual description of the cause.
  934. .. member:: MinimizerType Solver::Options::minimizer_type
  935. Default: ``TRUST_REGION``
  936. Choose between ``LINE_SEARCH`` and ``TRUST_REGION`` algorithms. See
  937. :ref:`section-trust-region-methods` and
  938. :ref:`section-line-search-methods` for more details.
  939. .. member:: LineSearchDirectionType Solver::Options::line_search_direction_type
  940. Default: ``LBFGS``
  941. Choices are ``STEEPEST_DESCENT``, ``NONLINEAR_CONJUGATE_GRADIENT``,
  942. ``BFGS`` and ``LBFGS``.
  943. See :ref:`section-line-search-methods` for more details.
  944. .. member:: LineSearchType Solver::Options::line_search_type
  945. Default: ``WOLFE``
  946. Choices are ``ARMIJO`` and ``WOLFE`` (strong Wolfe conditions).
  947. Note that in order for the assumptions underlying the ``BFGS`` and
  948. ``LBFGS`` line search direction algorithms to be satisfied, the
  949. ``WOLFE`` line search must be used.
  950. See :ref:`section-line-search-methods` for more details.
  951. .. member:: NonlinearConjugateGradientType Solver::Options::nonlinear_conjugate_gradient_type
  952. Default: ``FLETCHER_REEVES``
  953. Choices are ``FLETCHER_REEVES``, ``POLAK_RIBIERE`` and
  954. ``HESTENES_STIEFEL``.
  955. .. member:: int Solver::Options::max_lbfgs_rank
  956. Default: ``20``
  957. The LBFGS hessian approximation is a low rank approximation to
  958. the inverse of the Hessian matrix. The rank of the
  959. approximation determines (linearly) the space and time
  960. complexity of using the approximation. Higher the rank, the
  961. better is the quality of the approximation. The increase in
  962. quality is however is bounded for a number of reasons.
  963. 1. The method only uses secant information and not actual
  964. derivatives.
  965. 2. The Hessian approximation is constrained to be positive
  966. definite.
  967. So increasing this rank to a large number will cost time and
  968. space complexity without the corresponding increase in solution
  969. quality. There are no hard and fast rules for choosing the
  970. maximum rank. The best choice usually requires some problem
  971. specific experimentation.
  972. For more theoretical and implementation details of the LBFGS
  973. method, please see [Nocedal]_.
  974. .. member:: bool Solver::Options::use_approximate_eigenvalue_bfgs_scaling
  975. Default: ``false``
  976. As part of the ``BFGS`` update step / ``LBFGS`` right-multiply
  977. step, the initial inverse Hessian approximation is taken to be the
  978. Identity. However, [Oren]_ showed that using instead :math:`I *
  979. \gamma`, where :math:`\gamma` is a scalar chosen to approximate an
  980. eigenvalue of the true inverse Hessian can result in improved
  981. convergence in a wide variety of cases. Setting
  982. ``use_approximate_eigenvalue_bfgs_scaling`` to true enables this
  983. scaling in ``BFGS`` (before first iteration) and ``LBFGS`` (at each
  984. iteration).
  985. Precisely, approximate eigenvalue scaling equates to
  986. .. math:: \gamma = \frac{y_k' s_k}{y_k' y_k}
  987. With:
  988. .. math:: y_k = \nabla f_{k+1} - \nabla f_k
  989. .. math:: s_k = x_{k+1} - x_k
  990. Where :math:`f()` is the line search objective and :math:`x` the
  991. vector of parameter values [NocedalWright]_.
  992. It is important to note that approximate eigenvalue scaling does
  993. **not** *always* improve convergence, and that it can in fact
  994. *significantly* degrade performance for certain classes of problem,
  995. which is why it is disabled by default. In particular it can
  996. degrade performance when the sensitivity of the problem to different
  997. parameters varies significantly, as in this case a single scalar
  998. factor fails to capture this variation and detrimentally downscales
  999. parts of the Jacobian approximation which correspond to
  1000. low-sensitivity parameters. It can also reduce the robustness of the
  1001. solution to errors in the Jacobians.
  1002. .. member:: LineSearchIterpolationType Solver::Options::line_search_interpolation_type
  1003. Default: ``CUBIC``
  1004. Degree of the polynomial used to approximate the objective
  1005. function. Valid values are ``BISECTION``, ``QUADRATIC`` and
  1006. ``CUBIC``.
  1007. .. member:: double Solver::Options::min_line_search_step_size
  1008. Default: ``1e-9``
  1009. The line search terminates if:
  1010. .. math:: \|\Delta x_k\|_\infty < \text{min_line_search_step_size}
  1011. where :math:`\|\cdot\|_\infty` refers to the max norm, and
  1012. :math:`\Delta x_k` is the step change in the parameter values at
  1013. the :math:`k`-th iteration.
  1014. .. member:: double Solver::Options::line_search_sufficient_function_decrease
  1015. Default: ``1e-4``
  1016. Solving the line search problem exactly is computationally
  1017. prohibitive. Fortunately, line search based optimization algorithms
  1018. can still guarantee convergence if instead of an exact solution,
  1019. the line search algorithm returns a solution which decreases the
  1020. value of the objective function sufficiently. More precisely, we
  1021. are looking for a step size s.t.
  1022. .. math:: f(\text{step_size}) \le f(0) + \text{sufficient_decrease} * [f'(0) * \text{step_size}]
  1023. This condition is known as the Armijo condition.
  1024. .. member:: double Solver::Options::max_line_search_step_contraction
  1025. Default: ``1e-3``
  1026. In each iteration of the line search,
  1027. .. math:: \text{new_step_size} >= \text{max_line_search_step_contraction} * \text{step_size}
  1028. Note that by definition, for contraction:
  1029. .. math:: 0 < \text{max_step_contraction} < \text{min_step_contraction} < 1
  1030. .. member:: double Solver::Options::min_line_search_step_contraction
  1031. Default: ``0.6``
  1032. In each iteration of the line search,
  1033. .. math:: \text{new_step_size} <= \text{min_line_search_step_contraction} * \text{step_size}
  1034. Note that by definition, for contraction:
  1035. .. math:: 0 < \text{max_step_contraction} < \text{min_step_contraction} < 1
  1036. .. member:: int Solver::Options::max_num_line_search_step_size_iterations
  1037. Default: ``20``
  1038. Maximum number of trial step size iterations during each line
  1039. search, if a step size satisfying the search conditions cannot be
  1040. found within this number of trials, the line search will stop.
  1041. The minimum allowed value is 0 for trust region minimizer and 1
  1042. otherwise. If 0 is specified for the trust region minimizer, then
  1043. line search will not be used when solving constrained optimization
  1044. problems.
  1045. As this is an 'artificial' constraint (one imposed by the user, not
  1046. the underlying math), if ``WOLFE`` line search is being used, *and*
  1047. points satisfying the Armijo sufficient (function) decrease
  1048. condition have been found during the current search (in :math:`<=`
  1049. ``max_num_line_search_step_size_iterations``). Then, the step size
  1050. with the lowest function value which satisfies the Armijo condition
  1051. will be returned as the new valid step, even though it does *not*
  1052. satisfy the strong Wolfe conditions. This behaviour protects
  1053. against early termination of the optimizer at a sub-optimal point.
  1054. .. member:: int Solver::Options::max_num_line_search_direction_restarts
  1055. Default: ``5``
  1056. Maximum number of restarts of the line search direction algorithm
  1057. before terminating the optimization. Restarts of the line search
  1058. direction algorithm occur when the current algorithm fails to
  1059. produce a new descent direction. This typically indicates a
  1060. numerical failure, or a breakdown in the validity of the
  1061. approximations used.
  1062. .. member:: double Solver::Options::line_search_sufficient_curvature_decrease
  1063. Default: ``0.9``
  1064. The strong Wolfe conditions consist of the Armijo sufficient
  1065. decrease condition, and an additional requirement that the
  1066. step size be chosen s.t. the *magnitude* ('strong' Wolfe
  1067. conditions) of the gradient along the search direction
  1068. decreases sufficiently. Precisely, this second condition
  1069. is that we seek a step size s.t.
  1070. .. math:: \|f'(\text{step_size})\| <= \text{sufficient_curvature_decrease} * \|f'(0)\|
  1071. Where :math:`f()` is the line search objective and :math:`f'()` is the derivative
  1072. of :math:`f` with respect to the step size: :math:`\frac{d f}{d~\text{step size}}`.
  1073. .. member:: double Solver::Options::max_line_search_step_expansion
  1074. Default: ``10.0``
  1075. During the bracketing phase of a Wolfe line search, the step size
  1076. is increased until either a point satisfying the Wolfe conditions
  1077. is found, or an upper bound for a bracket containing a point
  1078. satisfying the conditions is found. Precisely, at each iteration
  1079. of the expansion:
  1080. .. math:: \text{new_step_size} <= \text{max_step_expansion} * \text{step_size}
  1081. By definition for expansion
  1082. .. math:: \text{max_step_expansion} > 1.0
  1083. .. member:: TrustRegionStrategyType Solver::Options::trust_region_strategy_type
  1084. Default: ``LEVENBERG_MARQUARDT``
  1085. The trust region step computation algorithm used by
  1086. Ceres. Currently ``LEVENBERG_MARQUARDT`` and ``DOGLEG`` are the two
  1087. valid choices. See :ref:`section-levenberg-marquardt` and
  1088. :ref:`section-dogleg` for more details.
  1089. .. member:: DoglegType Solver::Options::dogleg_type
  1090. Default: ``TRADITIONAL_DOGLEG``
  1091. Ceres supports two different dogleg strategies.
  1092. ``TRADITIONAL_DOGLEG`` method by Powell and the ``SUBSPACE_DOGLEG``
  1093. method described by [ByrdSchnabel]_ . See :ref:`section-dogleg`
  1094. for more details.
  1095. .. member:: bool Solver::Options::use_nonmonotonic_steps
  1096. Default: ``false``
  1097. Relax the requirement that the trust-region algorithm take strictly
  1098. decreasing steps. See :ref:`section-non-monotonic-steps` for more
  1099. details.
  1100. .. member:: int Solver::Options::max_consecutive_nonmonotonic_steps
  1101. Default: ``5``
  1102. The window size used by the step selection algorithm to accept
  1103. non-monotonic steps.
  1104. .. member:: int Solver::Options::max_num_iterations
  1105. Default: ``50``
  1106. Maximum number of iterations for which the solver should run.
  1107. .. member:: double Solver::Options::max_solver_time_in_seconds
  1108. Default: ``1e9``
  1109. Maximum amount of time for which the solver should run.
  1110. .. member:: int Solver::Options::num_threads
  1111. Default: ``1``
  1112. Number of threads used by Ceres to evaluate the Jacobian.
  1113. .. member:: double Solver::Options::initial_trust_region_radius
  1114. Default: ``1e4``
  1115. The size of the initial trust region. When the
  1116. ``LEVENBERG_MARQUARDT`` strategy is used, the reciprocal of this
  1117. number is the initial regularization parameter.
  1118. .. member:: double Solver::Options::max_trust_region_radius
  1119. Default: ``1e16``
  1120. The trust region radius is not allowed to grow beyond this value.
  1121. .. member:: double Solver::Options::min_trust_region_radius
  1122. Default: ``1e-32``
  1123. The solver terminates, when the trust region becomes smaller than
  1124. this value.
  1125. .. member:: double Solver::Options::min_relative_decrease
  1126. Default: ``1e-3``
  1127. Lower threshold for relative decrease before a trust-region step is
  1128. accepted.
  1129. .. member:: double Solver::Options::min_lm_diagonal
  1130. Default: ``1e-6``
  1131. The ``LEVENBERG_MARQUARDT`` strategy, uses a diagonal matrix to
  1132. regularize the trust region step. This is the lower bound on
  1133. the values of this diagonal matrix.
  1134. .. member:: double Solver::Options::max_lm_diagonal
  1135. Default: ``1e32``
  1136. The ``LEVENBERG_MARQUARDT`` strategy, uses a diagonal matrix to
  1137. regularize the trust region step. This is the upper bound on
  1138. the values of this diagonal matrix.
  1139. .. member:: int Solver::Options::max_num_consecutive_invalid_steps
  1140. Default: ``5``
  1141. The step returned by a trust region strategy can sometimes be
  1142. numerically invalid, usually because of conditioning
  1143. issues. Instead of crashing or stopping the optimization, the
  1144. optimizer can go ahead and try solving with a smaller trust
  1145. region/better conditioned problem. This parameter sets the number
  1146. of consecutive retries before the minimizer gives up.
  1147. .. member:: double Solver::Options::function_tolerance
  1148. Default: ``1e-6``
  1149. Solver terminates if
  1150. .. math:: \frac{|\Delta \text{cost}|}{\text{cost}} <= \text{function_tolerance}
  1151. where, :math:`\Delta \text{cost}` is the change in objective
  1152. function value (up or down) in the current iteration of
  1153. Levenberg-Marquardt.
  1154. .. member:: double Solver::Options::gradient_tolerance
  1155. Default: ``1e-10``
  1156. Solver terminates if
  1157. .. math:: \|x - \Pi \boxplus(x, -g(x))\|_\infty <= \text{gradient_tolerance}
  1158. where :math:`\|\cdot\|_\infty` refers to the max norm, :math:`\Pi`
  1159. is projection onto the bounds constraints and :math:`\boxplus` is
  1160. Plus operation for the overall manifold associated with the
  1161. parameter vector.
  1162. .. member:: double Solver::Options::parameter_tolerance
  1163. Default: ``1e-8``
  1164. Solver terminates if
  1165. .. math:: \|\Delta x\| <= (\|x\| + \text{parameter_tolerance}) * \text{parameter_tolerance}
  1166. where :math:`\Delta x` is the step computed by the linear solver in
  1167. the current iteration.
  1168. .. member:: LinearSolverType Solver::Options::linear_solver_type
  1169. Default: ``SPARSE_NORMAL_CHOLESKY`` / ``DENSE_QR``
  1170. Type of linear solver used to compute the solution to the linear
  1171. least squares problem in each iteration of the Levenberg-Marquardt
  1172. algorithm. If Ceres is built with support for ``SuiteSparse`` or
  1173. ``Accelerate`` or ``Eigen``'s sparse Cholesky factorization, the
  1174. default is ``SPARSE_NORMAL_CHOLESKY``, it is ``DENSE_QR``
  1175. otherwise.
  1176. .. member:: PreconditionerType Solver::Options::preconditioner_type
  1177. Default: ``JACOBI``
  1178. The preconditioner used by the iterative linear solver. The default
  1179. is the block Jacobi preconditioner. Valid values are (in increasing
  1180. order of complexity) ``IDENTITY``, ``JACOBI``, ``SCHUR_JACOBI``,
  1181. ``CLUSTER_JACOBI``, ``CLUSTER_TRIDIAGONAL``, ``SUBSET`` and
  1182. ``SCHUR_POWER_SERIES_EXPANSION``. See :ref:`section-preconditioner`
  1183. for more details.
  1184. .. member:: VisibilityClusteringType Solver::Options::visibility_clustering_type
  1185. Default: ``CANONICAL_VIEWS``
  1186. Type of clustering algorithm to use when constructing a visibility
  1187. based preconditioner. The original visibility based preconditioning
  1188. paper and implementation only used the canonical views algorithm.
  1189. This algorithm gives high quality results but for large dense
  1190. graphs can be particularly expensive. As its worst case complexity
  1191. is cubic in size of the graph.
  1192. Another option is to use ``SINGLE_LINKAGE`` which is a simple
  1193. thresholded single linkage clustering algorithm that only pays
  1194. attention to tightly coupled blocks in the Schur complement. This
  1195. is a fast algorithm that works well.
  1196. The optimal choice of the clustering algorithm depends on the
  1197. sparsity structure of the problem, but generally speaking we
  1198. recommend that you try ``CANONICAL_VIEWS`` first and if it is too
  1199. expensive try ``SINGLE_LINKAGE``.
  1200. .. member:: std::unordered_set<ResidualBlockId> Solver::Options::residual_blocks_for_subset_preconditioner
  1201. ``SUBSET`` preconditioner is a preconditioner for problems with
  1202. general sparsity. Given a subset of residual blocks of a problem,
  1203. it uses the corresponding subset of the rows of the Jacobian to
  1204. construct a preconditioner.
  1205. Suppose the Jacobian :math:`J` has been horizontally partitioned as
  1206. .. math:: J = \begin{bmatrix} P \\ Q \end{bmatrix}
  1207. Where, :math:`Q` is the set of rows corresponding to the residual
  1208. blocks in
  1209. :member:`Solver::Options::residual_blocks_for_subset_preconditioner`. The
  1210. preconditioner is the matrix :math:`(Q^\top Q)^{-1}`.
  1211. The efficacy of the preconditioner depends on how well the matrix
  1212. :math:`Q` approximates :math:`J^\top J`, or how well the chosen
  1213. residual blocks approximate the full problem.
  1214. If ``Solver::Options::preconditioner_type == SUBSET``, then
  1215. ``residual_blocks_for_subset_preconditioner`` must be non-empty.
  1216. .. member:: DenseLinearAlgebraLibrary Solver::Options::dense_linear_algebra_library_type
  1217. Default: ``EIGEN``
  1218. Ceres supports using multiple dense linear algebra libraries for
  1219. dense matrix factorizations. Currently ``EIGEN``, ``LAPACK`` and
  1220. ``CUDA`` are the valid choices. ``EIGEN`` is always available,
  1221. ``LAPACK`` refers to the system ``BLAS + LAPACK`` library which may
  1222. or may not be available. ``CUDA`` refers to Nvidia's GPU based
  1223. dense linear algebra library which may or may not be available.
  1224. This setting affects the ``DENSE_QR``, ``DENSE_NORMAL_CHOLESKY``
  1225. and ``DENSE_SCHUR`` solvers. For small to moderate sized problem
  1226. ``EIGEN`` is a fine choice but for large problems, an optimized
  1227. ``LAPACK + BLAS`` or ``CUDA`` implementation can make a substantial
  1228. difference in performance.
  1229. .. member:: SparseLinearAlgebraLibrary Solver::Options::sparse_linear_algebra_library_type
  1230. Default: The highest available according to: ``SUITE_SPARSE`` >
  1231. ``ACCELERATE_SPARSE`` > ``EIGEN_SPARSE`` > ``NO_SPARSE``
  1232. Ceres supports the use of three sparse linear algebra libraries,
  1233. ``SuiteSparse``, which is enabled by setting this parameter to
  1234. ``SUITE_SPARSE``, ``Acclerate``, which can be selected by setting
  1235. this parameter to ``ACCELERATE_SPARSE`` and ``Eigen`` which is
  1236. enabled by setting this parameter to ``EIGEN_SPARSE``. Lastly,
  1237. ``NO_SPARSE`` means that no sparse linear solver should be used;
  1238. note that this is irrespective of whether Ceres was compiled with
  1239. support for one.
  1240. ``SuiteSparse`` is a sophisticated sparse linear algebra library
  1241. and should be used in general. On MacOS you may want to use the
  1242. ``Accelerate`` framework.
  1243. If your needs/platforms prevent you from using ``SuiteSparse``,
  1244. consider using the sparse linear algebra routines in ``Eigen``. The
  1245. sparse Cholesky algorithms currently included with ``Eigen`` are
  1246. not as sophisticated as the ones in ``SuiteSparse`` and
  1247. ``Accelerate`` and as a result its performance is considerably
  1248. worse.
  1249. .. member:: LinearSolverOrderingType Solver::Options::linear_solver_ordering_type
  1250. Default: ``AMD``
  1251. The order in which variables are eliminated in a linear solver can
  1252. have a significant impact on the efficiency and accuracy of the
  1253. method. e.g., when doing sparse Cholesky factorization, there are
  1254. matrices for which a good ordering will give a Cholesky factor
  1255. with :math:`O(n)` storage, where as a bad ordering will result in
  1256. an completely dense factor.
  1257. Sparse direct solvers like ``SPARSE_NORMAL_CHOLESKY`` and
  1258. ``SPARSE_SCHUR`` use a fill reducing ordering of the columns and
  1259. rows of the matrix being factorized before computing the numeric
  1260. factorization.
  1261. This enum controls the type of algorithm used to compute this fill
  1262. reducing ordering. There is no single algorithm that works on all
  1263. matrices, so determining which algorithm works better is a matter
  1264. of empirical experimentation.
  1265. .. member:: std::shared_ptr<ParameterBlockOrdering> Solver::Options::linear_solver_ordering
  1266. Default: ``nullptr``
  1267. An instance of the ordering object informs the solver about the
  1268. desired order in which parameter blocks should be eliminated by the
  1269. linear solvers.
  1270. If ``nullptr``, the solver is free to choose an ordering that it
  1271. thinks is best.
  1272. See :ref:`section-ordering` for more details.
  1273. .. member:: bool Solver::Options::use_explicit_schur_complement
  1274. Default: ``false``
  1275. Use an explicitly computed Schur complement matrix with
  1276. ``ITERATIVE_SCHUR``.
  1277. By default this option is disabled and ``ITERATIVE_SCHUR``
  1278. evaluates evaluates matrix-vector products between the Schur
  1279. complement and a vector implicitly by exploiting the algebraic
  1280. expression for the Schur complement.
  1281. The cost of this evaluation scales with the number of non-zeros in
  1282. the Jacobian.
  1283. For small to medium sized problems there is a sweet spot where
  1284. computing the Schur complement is cheap enough that it is much more
  1285. efficient to explicitly compute it and use it for evaluating the
  1286. matrix-vector products.
  1287. .. NOTE::
  1288. This option can only be used with the ``SCHUR_JACOBI``
  1289. preconditioner.
  1290. .. member:: bool Solver::Options::dynamic_sparsity
  1291. Default: ``false``
  1292. Some non-linear least squares problems are symbolically dense but
  1293. numerically sparse. i.e. at any given state only a small number of
  1294. Jacobian entries are non-zero, but the position and number of
  1295. non-zeros is different depending on the state. For these problems
  1296. it can be useful to factorize the sparse jacobian at each solver
  1297. iteration instead of including all of the zero entries in a single
  1298. general factorization.
  1299. If your problem does not have this property (or you do not know),
  1300. then it is probably best to keep this false, otherwise it will
  1301. likely lead to worse performance.
  1302. This setting only affects the `SPARSE_NORMAL_CHOLESKY` solver.
  1303. .. member:: bool Solver::Options::use_mixed_precision_solves
  1304. Default: ``false``
  1305. If true, the Gauss-Newton matrix is computed in *double* precision, but
  1306. its factorization is computed in **single** precision. This can result in
  1307. significant time and memory savings at the cost of some accuracy in the
  1308. Gauss-Newton step. Iterative refinement is used to recover some
  1309. of this accuracy back.
  1310. If ``use_mixed_precision_solves`` is true, we recommend setting
  1311. ``max_num_refinement_iterations`` to 2-3.
  1312. See :ref:`section-mixed-precision` for more details.
  1313. .. member:: int Solver::Options::max_num_refinement_iterations
  1314. Default: ``0``
  1315. Number steps of the iterative refinement process to run when
  1316. computing the Gauss-Newton step, see
  1317. :member:`Solver::Options::use_mixed_precision_solves`.
  1318. .. member:: int Solver::Options::min_linear_solver_iterations
  1319. Default: ``0``
  1320. Minimum number of iterations used by the linear solver. This only
  1321. makes sense when the linear solver is an iterative solver, e.g.,
  1322. ``ITERATIVE_SCHUR`` or ``CGNR``.
  1323. .. member:: int Solver::Options::max_linear_solver_iterations
  1324. Default: ``500``
  1325. Minimum number of iterations used by the linear solver. This only
  1326. makes sense when the linear solver is an iterative solver, e.g.,
  1327. ``ITERATIVE_SCHUR`` or ``CGNR``.
  1328. .. member:: int Solver::Options::max_num_spse_iterations
  1329. Default: `5`
  1330. Maximum number of iterations performed by
  1331. ``SCHUR_POWER_SERIES_EXPANSION``. Each iteration corresponds to one
  1332. more term in the power series expansion od the inverse of the Schur
  1333. complement. This value controls the maximum number of iterations
  1334. whether it is used as a preconditioner or just to initialize the
  1335. solution for ``ITERATIVE_SCHUR``.
  1336. .. member:: bool Solver:Options::use_spse_initialization
  1337. Default: ``false``
  1338. Use Schur power series expansion to initialize the solution for
  1339. ``ITERATIVE_SCHUR``. This option can be set ``true`` regardless of
  1340. what preconditioner is being used.
  1341. .. member:: double Solver::Options::spse_tolerance
  1342. Default: `0.1`
  1343. When ``use_spse_initialization`` is ``true``, this parameter along
  1344. with ``max_num_spse_iterations`` controls the number of
  1345. ``SCHUR_POWER_SERIES_EXPANSION`` iterations performed for
  1346. initialization. It is not used to control the preconditioner.
  1347. .. member:: double Solver::Options::eta
  1348. Default: ``1e-1``
  1349. Forcing sequence parameter. The truncated Newton solver uses this
  1350. number to control the relative accuracy with which the Newton step
  1351. is computed. This constant is passed to
  1352. ``ConjugateGradientsSolver`` which uses it to terminate the
  1353. iterations when
  1354. .. math:: \frac{Q_i - Q_{i-1}}{Q_i} < \frac{\eta}{i}
  1355. .. member:: bool Solver::Options::jacobi_scaling
  1356. Default: ``true``
  1357. ``true`` means that the Jacobian is scaled by the norm of its
  1358. columns before being passed to the linear solver. This improves the
  1359. numerical conditioning of the normal equations.
  1360. .. member:: bool Solver::Options::use_inner_iterations
  1361. Default: ``false``
  1362. Use a non-linear version of a simplified variable projection
  1363. algorithm. Essentially this amounts to doing a further optimization
  1364. on each Newton/Trust region step using a coordinate descent
  1365. algorithm. For more details, see :ref:`section-inner-iterations`.
  1366. **Note** Inner iterations cannot be used with :class:`Problem`
  1367. objects that have an :class:`EvaluationCallback` associated with
  1368. them.
  1369. .. member:: std::shared_ptr<ParameterBlockOrdering> Solver::Options::inner_iteration_ordering
  1370. Default: ``nullptr``
  1371. If :member:`Solver::Options::use_inner_iterations` true, then the
  1372. user has two choices.
  1373. 1. Let the solver heuristically decide which parameter blocks to
  1374. optimize in each inner iteration. To do this, set
  1375. :member:`Solver::Options::inner_iteration_ordering` to ``nullptr``.
  1376. 2. Specify a collection of of ordered independent sets. The lower
  1377. numbered groups are optimized before the higher number groups
  1378. during the inner optimization phase. Each group must be an
  1379. independent set. Not all parameter blocks need to be included in
  1380. the ordering.
  1381. See :ref:`section-ordering` for more details.
  1382. .. member:: double Solver::Options::inner_iteration_tolerance
  1383. Default: ``1e-3``
  1384. Generally speaking, inner iterations make significant progress in
  1385. the early stages of the solve and then their contribution drops
  1386. down sharply, at which point the time spent doing inner iterations
  1387. is not worth it.
  1388. Once the relative decrease in the objective function due to inner
  1389. iterations drops below ``inner_iteration_tolerance``, the use of
  1390. inner iterations in subsequent trust region minimizer iterations is
  1391. disabled.
  1392. .. member:: LoggingType Solver::Options::logging_type
  1393. Default: ``PER_MINIMIZER_ITERATION``
  1394. Valid values are ``SILENT`` and ``PER_MINIMIZER_ITERATION``.
  1395. .. member:: bool Solver::Options::minimizer_progress_to_stdout
  1396. Default: ``false``
  1397. By default the Minimizer's progress is logged to ``STDERR``
  1398. depending on the ``vlog`` level. If this flag is set to true, and
  1399. :member:`Solver::Options::logging_type` is not ``SILENT``, the
  1400. logging output is sent to ``STDOUT``.
  1401. For ``TRUST_REGION_MINIMIZER`` the progress display looks like
  1402. .. code-block:: bash
  1403. iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time
  1404. 0 4.185660e+06 0.00e+00 1.09e+08 0.00e+00 0.00e+00 1.00e+04 0 7.59e-02 3.37e-01
  1405. 1 1.062590e+05 4.08e+06 8.99e+06 5.36e+02 9.82e-01 3.00e+04 1 1.65e-01 5.03e-01
  1406. 2 4.992817e+04 5.63e+04 8.32e+06 3.19e+02 6.52e-01 3.09e+04 1 1.45e-01 6.48e-01
  1407. Here
  1408. #. ``cost`` is the value of the objective function.
  1409. #. ``cost_change`` is the change in the value of the objective
  1410. function if the step computed in this iteration is accepted.
  1411. #. ``|gradient|`` is the max norm of the gradient.
  1412. #. ``|step|`` is the change in the parameter vector.
  1413. #. ``tr_ratio`` is the ratio of the actual change in the objective
  1414. function value to the change in the value of the trust
  1415. region model.
  1416. #. ``tr_radius`` is the size of the trust region radius.
  1417. #. ``ls_iter`` is the number of linear solver iterations used to
  1418. compute the trust region step. For direct/factorization based
  1419. solvers it is always 1, for iterative solvers like
  1420. ``ITERATIVE_SCHUR`` it is the number of iterations of the
  1421. Conjugate Gradients algorithm.
  1422. #. ``iter_time`` is the time take by the current iteration.
  1423. #. ``total_time`` is the total time taken by the minimizer.
  1424. For ``LINE_SEARCH_MINIMIZER`` the progress display looks like
  1425. .. code-block:: bash
  1426. 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
  1427. 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
  1428. 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
  1429. Here
  1430. #. ``f`` is the value of the objective function.
  1431. #. ``d`` is the change in the value of the objective function if
  1432. the step computed in this iteration is accepted.
  1433. #. ``g`` is the max norm of the gradient.
  1434. #. ``h`` is the change in the parameter vector.
  1435. #. ``s`` is the optimal step length computed by the line search.
  1436. #. ``it`` is the time take by the current iteration.
  1437. #. ``tt`` is the total time taken by the minimizer.
  1438. .. member:: std::vector<int> Solver::Options::trust_region_minimizer_iterations_to_dump
  1439. Default: ``empty``
  1440. List of iterations at which the trust region minimizer should dump
  1441. the trust region problem. Useful for testing and benchmarking. If
  1442. ``empty``, no problems are dumped.
  1443. .. member:: std::string Solver::Options::trust_region_problem_dump_directory
  1444. Default: ``/tmp``
  1445. Directory to which the problems should be written to. Should be
  1446. non-empty if
  1447. :member:`Solver::Options::trust_region_minimizer_iterations_to_dump` is
  1448. non-empty and
  1449. :member:`Solver::Options::trust_region_problem_dump_format_type` is not
  1450. ``CONSOLE``.
  1451. .. member:: DumpFormatType Solver::Options::trust_region_problem_dump_format_type
  1452. Default: ``TEXTFILE``
  1453. The format in which trust region problems should be logged when
  1454. :member:`Solver::Options::trust_region_minimizer_iterations_to_dump`
  1455. is non-empty. There are three options:
  1456. * ``CONSOLE`` prints the linear least squares problem in a human
  1457. readable format to ``stderr``. The Jacobian is printed as a
  1458. dense matrix. The vectors :math:`D`, :math:`x` and :math:`f` are
  1459. printed as dense vectors. This should only be used for small
  1460. problems.
  1461. * ``TEXTFILE`` Write out the linear least squares problem to the
  1462. directory pointed to by
  1463. :member:`Solver::Options::trust_region_problem_dump_directory` as
  1464. text files which can be read into ``MATLAB/Octave``. The Jacobian
  1465. is dumped as a text file containing :math:`(i,j,s)` triplets, the
  1466. vectors :math:`D`, `x` and `f` are dumped as text files
  1467. containing a list of their values.
  1468. A ``MATLAB/Octave`` script called
  1469. ``ceres_solver_iteration_???.m`` is also output, which can be
  1470. used to parse and load the problem into memory.
  1471. .. member:: bool Solver::Options::check_gradients
  1472. Default: ``false``
  1473. Check all Jacobians computed by each residual block with finite
  1474. differences. This is expensive since it involves computing the
  1475. derivative by normal means (e.g. user specified, autodiff, etc),
  1476. then also computing it using finite differences. The results are
  1477. compared, and if they differ substantially, the optimization fails
  1478. and the details are stored in the solver summary.
  1479. .. member:: double Solver::Options::gradient_check_relative_precision
  1480. Default: ``1e-8``
  1481. Precision to check for in the gradient checker. If the relative
  1482. difference between an element in a Jacobian exceeds this number,
  1483. then the Jacobian for that cost term is dumped.
  1484. .. member:: double Solver::Options::gradient_check_numeric_derivative_relative_step_size
  1485. Default: ``1e-6``
  1486. .. NOTE::
  1487. This option only applies to the numeric differentiation used for
  1488. checking the user provided derivatives when when
  1489. `Solver::Options::check_gradients` is true. If you are using
  1490. :class:`NumericDiffCostFunction` and are interested in changing
  1491. the step size for numeric differentiation in your cost function,
  1492. please have a look at :class:`NumericDiffOptions`.
  1493. Relative shift used for taking numeric derivatives when
  1494. `Solver::Options::check_gradients` is `true`.
  1495. For finite differencing, each dimension is evaluated at slightly
  1496. shifted values, e.g., for forward differences, the numerical
  1497. derivative is
  1498. .. math::
  1499. \delta &= gradient\_check\_numeric\_derivative\_relative\_step\_size\\
  1500. \Delta f &= \frac{f((1 + \delta) x) - f(x)}{\delta x}
  1501. The finite differencing is done along each dimension. The reason to
  1502. use a relative (rather than absolute) step size is that this way,
  1503. numeric differentiation works for functions where the arguments are
  1504. typically large (e.g. :math:`10^9`) and when the values are small
  1505. (e.g. :math:`10^{-5}`). It is possible to construct *torture cases*
  1506. which break this finite difference heuristic, but they do not come
  1507. up often in practice.
  1508. .. member:: bool Solver::Options::update_state_every_iteration
  1509. Default: ``false``
  1510. If ``update_state_every_iteration`` is ``true``, then Ceres Solver
  1511. will guarantee that at the end of every iteration and before any
  1512. user :class:`IterationCallback` is called, the parameter blocks are
  1513. updated to the current best solution found by the solver. Thus the
  1514. IterationCallback can inspect the values of the parameter blocks
  1515. for purposes of computation, visualization or termination.
  1516. If ``update_state_every_iteration`` is ``false`` then there is no
  1517. such guarantee, and user provided :class:`IterationCallback` s
  1518. should not expect to look at the parameter blocks and interpret
  1519. their values.
  1520. .. member:: std::vector<IterationCallback*> Solver::Options::callbacks
  1521. Default: ``empty``
  1522. Callbacks that are executed at the end of each iteration of the
  1523. minimizer. They are executed in the order that they are
  1524. specified in this vector.
  1525. By default, parameter blocks are updated only at the end of the
  1526. optimization, i.e., when the minimizer terminates. This means that
  1527. by default, if an :class:`IterationCallback` inspects the parameter
  1528. blocks, they will not see them changing in the course of the
  1529. optimization.
  1530. To tell Ceres to update the parameter blocks at the end of each
  1531. iteration and before calling the user's callback, set
  1532. :member:`Solver::Options::update_state_every_iteration` to
  1533. ``true``.
  1534. See `examples/iteration_callback_example.cc
  1535. <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/iteration_callback_example.cc>`_
  1536. for an example of an :class:`IterationCallback` that uses
  1537. :member:`Solver::Options::update_state_every_iteration` to log
  1538. changes to the parameter blocks over the course of the
  1539. optimization.
  1540. The solver does NOT take ownership of these pointers.
  1541. :class:`ParameterBlockOrdering`
  1542. ===============================
  1543. .. class:: ParameterBlockOrdering
  1544. ``ParameterBlockOrdering`` is a class for storing and manipulating
  1545. an ordered collection of groups/sets with the following semantics:
  1546. Group IDs are non-negative integer values. Elements are any type
  1547. that can serve as a key in a map or an element of a set.
  1548. An element can only belong to one group at a time. A group may
  1549. contain an arbitrary number of elements.
  1550. Groups are ordered by their group id.
  1551. .. function:: bool ParameterBlockOrdering::AddElementToGroup(const double* element, const int group)
  1552. Add an element to a group. If a group with this id does not exist,
  1553. one is created. This method can be called any number of times for
  1554. the same element. Group ids should be non-negative numbers. Return
  1555. value indicates if adding the element was a success.
  1556. .. function:: void ParameterBlockOrdering::Clear()
  1557. Clear the ordering.
  1558. .. function:: bool ParameterBlockOrdering::Remove(const double* element)
  1559. Remove the element, no matter what group it is in. If the element
  1560. is not a member of any group, calling this method will result in a
  1561. crash. Return value indicates if the element was actually removed.
  1562. .. function:: void ParameterBlockOrdering::Reverse()
  1563. Reverse the order of the groups in place.
  1564. .. function:: int ParameterBlockOrdering::GroupId(const double* element) const
  1565. Return the group id for the element. If the element is not a member
  1566. of any group, return -1.
  1567. .. function:: bool ParameterBlockOrdering::IsMember(const double* element) const
  1568. True if there is a group containing the parameter block.
  1569. .. function:: int ParameterBlockOrdering::GroupSize(const int group) const
  1570. This function always succeeds, i.e., implicitly there exists a
  1571. group for every integer.
  1572. .. function:: int ParameterBlockOrdering::NumElements() const
  1573. Number of elements in the ordering.
  1574. .. function:: int ParameterBlockOrdering::NumGroups() const
  1575. Number of groups with one or more elements.
  1576. :class:`IterationSummary`
  1577. =========================
  1578. .. class:: IterationSummary
  1579. :class:`IterationSummary` describes the state of the minimizer at
  1580. the end of each iteration.
  1581. .. member:: int IterationSummary::iteration
  1582. Current iteration number.
  1583. .. member:: bool IterationSummary::step_is_valid
  1584. Step was numerically valid, i.e., all values are finite and the
  1585. step reduces the value of the linearized model.
  1586. **Note**: :member:`IterationSummary::step_is_valid` is `false`
  1587. when :member:`IterationSummary::iteration` = 0.
  1588. .. member:: bool IterationSummary::step_is_nonmonotonic
  1589. Step did not reduce the value of the objective function
  1590. sufficiently, but it was accepted because of the relaxed
  1591. acceptance criterion used by the non-monotonic trust region
  1592. algorithm.
  1593. **Note**: :member:`IterationSummary::step_is_nonmonotonic` is
  1594. `false` when when :member:`IterationSummary::iteration` = 0.
  1595. .. member:: bool IterationSummary::step_is_successful
  1596. Whether or not the minimizer accepted this step or not.
  1597. If the ordinary trust region algorithm is used, this means that the
  1598. relative reduction in the objective function value was greater than
  1599. :member:`Solver::Options::min_relative_decrease`. However, if the
  1600. non-monotonic trust region algorithm is used
  1601. (:member:`Solver::Options::use_nonmonotonic_steps` = `true`), then
  1602. even if the relative decrease is not sufficient, the algorithm may
  1603. accept the step and the step is declared successful.
  1604. **Note**: :member:`IterationSummary::step_is_successful` is `false`
  1605. when when :member:`IterationSummary::iteration` = 0.
  1606. .. member:: double IterationSummary::cost
  1607. Value of the objective function.
  1608. .. member:: double IterationSummary::cost_change
  1609. Change in the value of the objective function in this
  1610. iteration. This can be positive or negative.
  1611. .. member:: double IterationSummary::gradient_max_norm
  1612. Infinity norm of the gradient vector.
  1613. .. member:: double IterationSummary::gradient_norm
  1614. 2-norm of the gradient vector.
  1615. .. member:: double IterationSummary::step_norm
  1616. 2-norm of the size of the step computed in this iteration.
  1617. .. member:: double IterationSummary::relative_decrease
  1618. For trust region algorithms, the ratio of the actual change in cost
  1619. and the change in the cost of the linearized approximation.
  1620. This field is not used when a linear search minimizer is used.
  1621. .. member:: double IterationSummary::trust_region_radius
  1622. Size of the trust region at the end of the current iteration. For
  1623. the Levenberg-Marquardt algorithm, the regularization parameter is
  1624. 1.0 / :member:`IterationSummary::trust_region_radius`.
  1625. .. member:: double IterationSummary::eta
  1626. For the inexact step Levenberg-Marquardt algorithm, this is the
  1627. relative accuracy with which the step is solved. This number is
  1628. only applicable to the iterative solvers capable of solving linear
  1629. systems inexactly. Factorization-based exact solvers always have an
  1630. eta of 0.0.
  1631. .. member:: double IterationSummary::step_size
  1632. Step sized computed by the line search algorithm.
  1633. This field is not used when a trust region minimizer is used.
  1634. .. member:: int IterationSummary::line_search_function_evaluations
  1635. Number of function evaluations used by the line search algorithm.
  1636. This field is not used when a trust region minimizer is used.
  1637. .. member:: int IterationSummary::linear_solver_iterations
  1638. Number of iterations taken by the linear solver to solve for the
  1639. trust region step.
  1640. Currently this field is not used when a line search minimizer is
  1641. used.
  1642. .. member:: double IterationSummary::iteration_time_in_seconds
  1643. Time (in seconds) spent inside the minimizer loop in the current
  1644. iteration.
  1645. .. member:: double IterationSummary::step_solver_time_in_seconds
  1646. Time (in seconds) spent inside the trust region step solver.
  1647. .. member:: double IterationSummary::cumulative_time_in_seconds
  1648. Time (in seconds) since the user called Solve().
  1649. :class:`IterationCallback`
  1650. ==========================
  1651. .. class:: IterationCallback
  1652. Interface for specifying callbacks that are executed at the end of
  1653. each iteration of the minimizer.
  1654. .. code-block:: c++
  1655. class IterationCallback {
  1656. public:
  1657. virtual ~IterationCallback() {}
  1658. virtual CallbackReturnType operator()(const IterationSummary& summary) = 0;
  1659. };
  1660. The solver uses the return value of ``operator()`` to decide whether
  1661. to continue solving or to terminate. The user can return three
  1662. values.
  1663. #. ``SOLVER_ABORT`` indicates that the callback detected an abnormal
  1664. situation. The solver returns without updating the parameter
  1665. blocks (unless ``Solver::Options::update_state_every_iteration`` is
  1666. set true). Solver returns with ``Solver::Summary::termination_type``
  1667. set to ``USER_FAILURE``.
  1668. #. ``SOLVER_TERMINATE_SUCCESSFULLY`` indicates that there is no need
  1669. to optimize anymore (some user specified termination criterion
  1670. has been met). Solver returns with
  1671. ``Solver::Summary::termination_type``` set to ``USER_SUCCESS``.
  1672. #. ``SOLVER_CONTINUE`` indicates that the solver should continue
  1673. optimizing.
  1674. The return values can be used to implement custom termination
  1675. criterion that supercede the iteration/time/tolerance based
  1676. termination implemented by Ceres.
  1677. For example, the following :class:`IterationCallback` is used
  1678. internally by Ceres to log the progress of the optimization.
  1679. .. code-block:: c++
  1680. class LoggingCallback : public IterationCallback {
  1681. public:
  1682. explicit LoggingCallback(bool log_to_stdout)
  1683. : log_to_stdout_(log_to_stdout) {}
  1684. ~LoggingCallback() {}
  1685. CallbackReturnType operator()(const IterationSummary& summary) {
  1686. const char* kReportRowFormat =
  1687. "% 4d: f:% 8e d:% 3.2e g:% 3.2e h:% 3.2e "
  1688. "rho:% 3.2e mu:% 3.2e eta:% 3.2e li:% 3d";
  1689. string output = StringPrintf(kReportRowFormat,
  1690. summary.iteration,
  1691. summary.cost,
  1692. summary.cost_change,
  1693. summary.gradient_max_norm,
  1694. summary.step_norm,
  1695. summary.relative_decrease,
  1696. summary.trust_region_radius,
  1697. summary.eta,
  1698. summary.linear_solver_iterations);
  1699. if (log_to_stdout_) {
  1700. cout << output << endl;
  1701. } else {
  1702. VLOG(1) << output;
  1703. }
  1704. return SOLVER_CONTINUE;
  1705. }
  1706. private:
  1707. const bool log_to_stdout_;
  1708. };
  1709. See `examples/evaluation_callback_example.cc
  1710. <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/iteration_callback_example.cc>`_
  1711. for another example that uses
  1712. :member:`Solver::Options::update_state_every_iteration` to log
  1713. changes to the parameter blocks over the course of the optimization.
  1714. :class:`CRSMatrix`
  1715. ==================
  1716. .. class:: CRSMatrix
  1717. A compressed row sparse matrix used primarily for communicating the
  1718. Jacobian matrix to the user.
  1719. .. member:: int CRSMatrix::num_rows
  1720. Number of rows.
  1721. .. member:: int CRSMatrix::num_cols
  1722. Number of columns.
  1723. .. member:: std::vector<int> CRSMatrix::rows
  1724. :member:`CRSMatrix::rows` is a :member:`CRSMatrix::num_rows` + 1
  1725. sized array that points into the :member:`CRSMatrix::cols` and
  1726. :member:`CRSMatrix::values` array.
  1727. .. member:: std::vector<int> CRSMatrix::cols
  1728. :member:`CRSMatrix::cols` contain as many entries as there are
  1729. non-zeros in the matrix.
  1730. For each row ``i``, ``cols[rows[i]]`` ... ``cols[rows[i + 1] - 1]``
  1731. are the indices of the non-zero columns of row ``i``.
  1732. .. member:: std::vector<double> CRSMatrix::values
  1733. :member:`CRSMatrix::values` contain as many entries as there are
  1734. non-zeros in the matrix.
  1735. For each row ``i``,
  1736. ``values[rows[i]]`` ... ``values[rows[i + 1] - 1]`` are the values
  1737. of the non-zero columns of row ``i``.
  1738. e.g., consider the 3x4 sparse matrix
  1739. .. code-block:: c++
  1740. 0 10 0 4
  1741. 0 2 -3 2
  1742. 1 2 0 0
  1743. The three arrays will be:
  1744. .. code-block:: c++
  1745. -row0- ---row1--- -row2-
  1746. rows = [ 0, 2, 5, 7]
  1747. cols = [ 1, 3, 1, 2, 3, 0, 1]
  1748. values = [10, 4, 2, -3, 2, 1, 2]
  1749. :class:`Solver::Summary`
  1750. ========================
  1751. .. class:: Solver::Summary
  1752. Summary of the various stages of the solver after termination.
  1753. .. function:: std::string Solver::Summary::BriefReport() const
  1754. A brief one line description of the state of the solver after
  1755. termination.
  1756. .. function:: std::string Solver::Summary::FullReport() const
  1757. A full multiline description of the state of the solver after
  1758. termination.
  1759. .. function:: bool Solver::Summary::IsSolutionUsable() const
  1760. Whether the solution returned by the optimization algorithm can be
  1761. relied on to be numerically sane. This will be the case if
  1762. `Solver::Summary:termination_type` is set to `CONVERGENCE`,
  1763. `USER_SUCCESS` or `NO_CONVERGENCE`, i.e., either the solver
  1764. converged by meeting one of the convergence tolerances or because
  1765. the user indicated that it had converged or it ran to the maximum
  1766. number of iterations or time.
  1767. .. member:: MinimizerType Solver::Summary::minimizer_type
  1768. Type of minimization algorithm used.
  1769. .. member:: TerminationType Solver::Summary::termination_type
  1770. The cause of the minimizer terminating.
  1771. .. member:: std::string Solver::Summary::message
  1772. Reason why the solver terminated.
  1773. .. member:: double Solver::Summary::initial_cost
  1774. Cost of the problem (value of the objective function) before the
  1775. optimization.
  1776. .. member:: double Solver::Summary::final_cost
  1777. Cost of the problem (value of the objective function) after the
  1778. optimization.
  1779. .. member:: double Solver::Summary::fixed_cost
  1780. The part of the total cost that comes from residual blocks that
  1781. were held fixed by the preprocessor because all the parameter
  1782. blocks that they depend on were fixed.
  1783. .. member:: std::vector<IterationSummary> Solver::Summary::iterations
  1784. :class:`IterationSummary` for each minimizer iteration in order.
  1785. .. member:: int Solver::Summary::num_successful_steps
  1786. Number of minimizer iterations in which the step was
  1787. accepted. Unless :member:`Solver::Options::use_nonmonotonic_steps`
  1788. is `true` this is also the number of steps in which the objective
  1789. function value/cost went down.
  1790. .. member:: int Solver::Summary::num_unsuccessful_steps
  1791. Number of minimizer iterations in which the step was rejected
  1792. either because it did not reduce the cost enough or the step was
  1793. not numerically valid.
  1794. .. member:: int Solver::Summary::num_inner_iteration_steps
  1795. Number of times inner iterations were performed.
  1796. .. member:: int Solver::Summary::num_line_search_steps
  1797. Total number of iterations inside the line search algorithm across
  1798. all invocations. We call these iterations "steps" to distinguish
  1799. them from the outer iterations of the line search and trust region
  1800. minimizer algorithms which call the line search algorithm as a
  1801. subroutine.
  1802. .. member:: double Solver::Summary::preprocessor_time_in_seconds
  1803. Time (in seconds) spent in the preprocessor.
  1804. .. member:: double Solver::Summary::minimizer_time_in_seconds
  1805. Time (in seconds) spent in the minimizer.
  1806. .. member:: double Solver::Summary::postprocessor_time_in_seconds
  1807. Time (in seconds) spent in the post processor.
  1808. .. member:: double Solver::Summary::total_time_in_seconds
  1809. Time (in seconds) spent in the solver.
  1810. .. member:: double Solver::Summary::linear_solver_time_in_seconds
  1811. Time (in seconds) spent in the linear solver computing the trust
  1812. region step.
  1813. .. member:: int Solver::Summary::num_linear_solves
  1814. Number of times the Newton step was computed by solving a linear
  1815. system. This does not include linear solves used by inner
  1816. iterations.
  1817. .. member:: double Solver::Summary::residual_evaluation_time_in_seconds
  1818. Time (in seconds) spent evaluating the residual vector.
  1819. .. member:: int Solver::Summary::num_residual_evaluations
  1820. Number of times only the residuals were evaluated.
  1821. .. member:: double Solver::Summary::jacobian_evaluation_time_in_seconds
  1822. Time (in seconds) spent evaluating the Jacobian matrix.
  1823. .. member:: int Solver::Summary::num_jacobian_evaluations
  1824. Number of times only the Jacobian and the residuals were evaluated.
  1825. .. member:: double Solver::Summary::inner_iteration_time_in_seconds
  1826. Time (in seconds) spent doing inner iterations.
  1827. .. member:: int Solver::Summary::num_parameter_blocks
  1828. Number of parameter blocks in the problem.
  1829. .. member:: int Solver::Summary::num_parameters
  1830. Number of parameters in the problem.
  1831. .. member:: int Solver::Summary::num_effective_parameters
  1832. Dimension of the tangent space of the problem (or the number of
  1833. columns in the Jacobian for the problem). This is different from
  1834. :member:`Solver::Summary::num_parameters` if a parameter block is
  1835. associated with a :class:`Manifold`.
  1836. .. member:: int Solver::Summary::num_residual_blocks
  1837. Number of residual blocks in the problem.
  1838. .. member:: int Solver::Summary::num_residuals
  1839. Number of residuals in the problem.
  1840. .. member:: int Solver::Summary::num_parameter_blocks_reduced
  1841. Number of parameter blocks in the problem after the inactive and
  1842. constant parameter blocks have been removed. A parameter block is
  1843. inactive if no residual block refers to it.
  1844. .. member:: int Solver::Summary::num_parameters_reduced
  1845. Number of parameters in the reduced problem.
  1846. .. member:: int Solver::Summary::num_effective_parameters_reduced
  1847. Dimension of the tangent space of the reduced problem (or the
  1848. number of columns in the Jacobian for the reduced problem). This is
  1849. different from :member:`Solver::Summary::num_parameters_reduced` if
  1850. a parameter block in the reduced problem is associated with a
  1851. :class:`Manifold`.
  1852. .. member:: int Solver::Summary::num_residual_blocks_reduced
  1853. Number of residual blocks in the reduced problem.
  1854. .. member:: int Solver::Summary::num_residuals_reduced
  1855. Number of residuals in the reduced problem.
  1856. .. member:: int Solver::Summary::num_threads_given
  1857. Number of threads specified by the user for Jacobian and residual
  1858. evaluation.
  1859. .. member:: int Solver::Summary::num_threads_used
  1860. Number of threads actually used by the solver for Jacobian and
  1861. residual evaluation.
  1862. .. member:: LinearSolverType Solver::Summary::linear_solver_type_given
  1863. Type of the linear solver requested by the user.
  1864. .. member:: LinearSolverType Solver::Summary::linear_solver_type_used
  1865. Type of the linear solver actually used. This may be different from
  1866. :member:`Solver::Summary::linear_solver_type_given` if Ceres
  1867. determines that the problem structure is not compatible with the
  1868. linear solver requested or if the linear solver requested by the
  1869. user is not available, e.g. The user requested
  1870. `SPARSE_NORMAL_CHOLESKY` but no sparse linear algebra library was
  1871. available.
  1872. .. member:: std::vector<int> Solver::Summary::linear_solver_ordering_given
  1873. Size of the elimination groups given by the user as hints to the
  1874. linear solver.
  1875. .. member:: std::vector<int> Solver::Summary::linear_solver_ordering_used
  1876. Size of the parameter groups used by the solver when ordering the
  1877. columns of the Jacobian. This maybe different from
  1878. :member:`Solver::Summary::linear_solver_ordering_given` if the user
  1879. left :member:`Solver::Summary::linear_solver_ordering_given` blank
  1880. and asked for an automatic ordering, or if the problem contains
  1881. some constant or inactive parameter blocks.
  1882. .. member:: std::string Solver::Summary::schur_structure_given
  1883. For Schur type linear solvers, this string describes the template
  1884. specialization which was detected in the problem and should be
  1885. used.
  1886. .. member:: std::string Solver::Summary::schur_structure_used
  1887. For Schur type linear solvers, this string describes the template
  1888. specialization that was actually instantiated and used. The reason
  1889. this will be different from
  1890. :member:`Solver::Summary::schur_structure_given` is because the
  1891. corresponding template specialization does not exist.
  1892. Template specializations can be added to ceres by editing
  1893. ``internal/ceres/generate_template_specializations.py``
  1894. .. member:: bool Solver::Summary::inner_iterations_given
  1895. `True` if the user asked for inner iterations to be used as part of
  1896. the optimization.
  1897. .. member:: bool Solver::Summary::inner_iterations_used
  1898. `True` if the user asked for inner iterations to be used as part of
  1899. the optimization and the problem structure was such that they were
  1900. actually performed. For example, in a problem with just one parameter
  1901. block, inner iterations are not performed.
  1902. .. member:: std::vector<int> Solver::Summary::inner_iteration_ordering_given
  1903. Size of the parameter groups given by the user for performing inner
  1904. iterations.
  1905. .. member:: std::vector<int> Solver::Summary::inner_iteration_ordering_used
  1906. Size of the parameter groups given used by the solver for
  1907. performing inner iterations. This maybe different from
  1908. :member:`Solver::Summary::inner_iteration_ordering_given` if the
  1909. user left :member:`Solver::Summary::inner_iteration_ordering_given`
  1910. blank and asked for an automatic ordering, or if the problem
  1911. contains some constant or inactive parameter blocks.
  1912. .. member:: PreconditionerType Solver::Summary::preconditioner_type_given
  1913. Type of the preconditioner requested by the user.
  1914. .. member:: PreconditionerType Solver::Summary::preconditioner_type_used
  1915. Type of the preconditioner actually used. This may be different
  1916. from :member:`Solver::Summary::linear_solver_type_given` if Ceres
  1917. determines that the problem structure is not compatible with the
  1918. linear solver requested or if the linear solver requested by the
  1919. user is not available.
  1920. .. member:: VisibilityClusteringType Solver::Summary::visibility_clustering_type
  1921. Type of clustering algorithm used for visibility based
  1922. preconditioning. Only meaningful when the
  1923. :member:`Solver::Summary::preconditioner_type_used` is
  1924. ``CLUSTER_JACOBI`` or ``CLUSTER_TRIDIAGONAL``.
  1925. .. member:: TrustRegionStrategyType Solver::Summary::trust_region_strategy_type
  1926. Type of trust region strategy.
  1927. .. member:: DoglegType Solver::Summary::dogleg_type
  1928. Type of dogleg strategy used for solving the trust region problem.
  1929. .. member:: DenseLinearAlgebraLibraryType Solver::Summary::dense_linear_algebra_library_type
  1930. Type of the dense linear algebra library used.
  1931. .. member:: SparseLinearAlgebraLibraryType Solver::Summary::sparse_linear_algebra_library_type
  1932. Type of the sparse linear algebra library used.
  1933. .. member:: LineSearchDirectionType Solver::Summary::line_search_direction_type
  1934. Type of line search direction used.
  1935. .. member:: LineSearchType Solver::Summary::line_search_type
  1936. Type of the line search algorithm used.
  1937. .. member:: LineSearchInterpolationType Solver::Summary::line_search_interpolation_type
  1938. When performing line search, the degree of the polynomial used to
  1939. approximate the objective function.
  1940. .. member:: NonlinearConjugateGradientType Solver::Summary::nonlinear_conjugate_gradient_type
  1941. If the line search direction is `NONLINEAR_CONJUGATE_GRADIENT`,
  1942. then this indicates the particular variant of non-linear conjugate
  1943. gradient used.
  1944. .. member:: int Solver::Summary::max_lbfgs_rank
  1945. If the type of the line search direction is `LBFGS`, then this
  1946. indicates the rank of the Hessian approximation.