solver.h revision 79397c21138f54fcff6ec067b44b847f1f7e0e98
15f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// Ceres Solver - A fast non-linear least squares minimizer
25f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// Copyright 2014 Google Inc. All rights reserved.
35f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// http://code.google.com/p/ceres-solver/
45f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)//
55f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// Redistribution and use in source and binary forms, with or without
65f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// modification, are permitted provided that the following conditions are met:
75f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)//
85f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// * Redistributions of source code must retain the above copyright notice,
95f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)//   this list of conditions and the following disclaimer.
105f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// * Redistributions in binary form must reproduce the above copyright notice,
115f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)//   this list of conditions and the following disclaimer in the documentation
125f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)//   and/or other materials provided with the distribution.
135f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// * Neither the name of Google Inc. nor the names of its contributors may be
145f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)//   used to endorse or promote products derived from this software without
155f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)//   specific prior written permission.
165f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)//
175f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
185f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
195f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
205f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27// POSSIBILITY OF SUCH DAMAGE.
28//
29// Author: sameeragarwal@google.com (Sameer Agarwal)
30
31#ifndef CERES_PUBLIC_SOLVER_H_
32#define CERES_PUBLIC_SOLVER_H_
33
34#include <cmath>
35#include <string>
36#include <vector>
37#include "ceres/crs_matrix.h"
38#include "ceres/internal/macros.h"
39#include "ceres/internal/port.h"
40#include "ceres/iteration_callback.h"
41#include "ceres/ordered_groups.h"
42#include "ceres/types.h"
43#include "ceres/internal/disable_warnings.h"
44
45namespace ceres {
46
47class Problem;
48
49// Interface for non-linear least squares solvers.
50class CERES_EXPORT Solver {
51 public:
52  virtual ~Solver();
53
54  // The options structure contains, not surprisingly, options that control how
55  // the solver operates. The defaults should be suitable for a wide range of
56  // problems; however, better performance is often obtainable with tweaking.
57  //
58  // The constants are defined inside types.h
59  struct CERES_EXPORT Options {
60    // Default constructor that sets up a generic sparse problem.
61    Options() {
62      minimizer_type = TRUST_REGION;
63      line_search_direction_type = LBFGS;
64      line_search_type = WOLFE;
65      nonlinear_conjugate_gradient_type = FLETCHER_REEVES;
66      max_lbfgs_rank = 20;
67      use_approximate_eigenvalue_bfgs_scaling = false;
68      line_search_interpolation_type = CUBIC;
69      min_line_search_step_size = 1e-9;
70      line_search_sufficient_function_decrease = 1e-4;
71      max_line_search_step_contraction = 1e-3;
72      min_line_search_step_contraction = 0.6;
73      max_num_line_search_step_size_iterations = 20;
74      max_num_line_search_direction_restarts = 5;
75      line_search_sufficient_curvature_decrease = 0.9;
76      max_line_search_step_expansion = 10.0;
77      trust_region_strategy_type = LEVENBERG_MARQUARDT;
78      dogleg_type = TRADITIONAL_DOGLEG;
79      use_nonmonotonic_steps = false;
80      max_consecutive_nonmonotonic_steps = 5;
81      max_num_iterations = 50;
82      max_solver_time_in_seconds = 1e9;
83      num_threads = 1;
84      initial_trust_region_radius = 1e4;
85      max_trust_region_radius = 1e16;
86      min_trust_region_radius = 1e-32;
87      min_relative_decrease = 1e-3;
88      min_lm_diagonal = 1e-6;
89      max_lm_diagonal = 1e32;
90      max_num_consecutive_invalid_steps = 5;
91      function_tolerance = 1e-6;
92      gradient_tolerance = 1e-10;
93      parameter_tolerance = 1e-8;
94
95#if defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CXSPARSE) && !defined(CERES_ENABLE_LGPL_CODE)
96      linear_solver_type = DENSE_QR;
97#else
98      linear_solver_type = SPARSE_NORMAL_CHOLESKY;
99#endif
100
101      preconditioner_type = JACOBI;
102      visibility_clustering_type = CANONICAL_VIEWS;
103      dense_linear_algebra_library_type = EIGEN;
104
105      // Choose a default sparse linear algebra library in the order:
106      //
107      //   SUITE_SPARSE > CX_SPARSE > EIGEN_SPARSE
108#if !defined(CERES_NO_SUITESPARSE)
109      sparse_linear_algebra_library_type = SUITE_SPARSE;
110#else
111  #if !defined(CERES_NO_CXSPARSE)
112      sparse_linear_algebra_library_type = CX_SPARSE;
113  #else
114    #if defined(CERES_USE_EIGEN_SPARSE)
115      sparse_linear_algebra_library_type = EIGEN_SPARSE;
116    #endif
117  #endif
118#endif
119
120      num_linear_solver_threads = 1;
121      use_postordering = false;
122      dynamic_sparsity = false;
123      min_linear_solver_iterations = 1;
124      max_linear_solver_iterations = 500;
125      eta = 1e-1;
126      jacobi_scaling = true;
127      use_inner_iterations = false;
128      inner_iteration_tolerance = 1e-3;
129      logging_type = PER_MINIMIZER_ITERATION;
130      minimizer_progress_to_stdout = false;
131      trust_region_problem_dump_directory = "/tmp";
132      trust_region_problem_dump_format_type = TEXTFILE;
133      check_gradients = false;
134      gradient_check_relative_precision = 1e-8;
135      numeric_derivative_relative_step_size = 1e-6;
136      update_state_every_iteration = false;
137    }
138
139    // Returns true if the options struct has a valid
140    // configuration. Returns false otherwise, and fills in *error
141    // with a message describing the problem.
142    bool IsValid(string* error) const;
143
144    // Minimizer options ----------------------------------------
145
146    // Ceres supports the two major families of optimization strategies -
147    // Trust Region and Line Search.
148    //
149    // 1. The line search approach first finds a descent direction
150    // along which the objective function will be reduced and then
151    // computes a step size that decides how far should move along
152    // that direction. The descent direction can be computed by
153    // various methods, such as gradient descent, Newton's method and
154    // Quasi-Newton method. The step size can be determined either
155    // exactly or inexactly.
156    //
157    // 2. The trust region approach approximates the objective
158    // function using using a model function (often a quadratic) over
159    // a subset of the search space known as the trust region. If the
160    // model function succeeds in minimizing the true objective
161    // function the trust region is expanded; conversely, otherwise it
162    // is contracted and the model optimization problem is solved
163    // again.
164    //
165    // Trust region methods are in some sense dual to line search methods:
166    // trust region methods first choose a step size (the size of the
167    // trust region) and then a step direction while line search methods
168    // first choose a step direction and then a step size.
169    MinimizerType minimizer_type;
170
171    LineSearchDirectionType line_search_direction_type;
172    LineSearchType line_search_type;
173    NonlinearConjugateGradientType nonlinear_conjugate_gradient_type;
174
175    // The LBFGS hessian approximation is a low rank approximation to
176    // the inverse of the Hessian matrix. The rank of the
177    // approximation determines (linearly) the space and time
178    // complexity of using the approximation. Higher the rank, the
179    // better is the quality of the approximation. The increase in
180    // quality is however is bounded for a number of reasons.
181    //
182    // 1. The method only uses secant information and not actual
183    // derivatives.
184    //
185    // 2. The Hessian approximation is constrained to be positive
186    // definite.
187    //
188    // So increasing this rank to a large number will cost time and
189    // space complexity without the corresponding increase in solution
190    // quality. There are no hard and fast rules for choosing the
191    // maximum rank. The best choice usually requires some problem
192    // specific experimentation.
193    //
194    // For more theoretical and implementation details of the LBFGS
195    // method, please see:
196    //
197    // Nocedal, J. (1980). "Updating Quasi-Newton Matrices with
198    // Limited Storage". Mathematics of Computation 35 (151): 773–782.
199    int max_lbfgs_rank;
200
201    // As part of the (L)BFGS update step (BFGS) / right-multiply step (L-BFGS),
202    // the initial inverse Hessian approximation is taken to be the Identity.
203    // However, Oren showed that using instead I * \gamma, where \gamma is
204    // chosen to approximate an eigenvalue of the true inverse Hessian can
205    // result in improved convergence in a wide variety of cases. Setting
206    // use_approximate_eigenvalue_bfgs_scaling to true enables this scaling.
207    //
208    // It is important to note that approximate eigenvalue scaling does not
209    // always improve convergence, and that it can in fact significantly degrade
210    // performance for certain classes of problem, which is why it is disabled
211    // by default.  In particular it can degrade performance when the
212    // sensitivity of the problem to different parameters varies significantly,
213    // as in this case a single scalar factor fails to capture this variation
214    // and detrimentally downscales parts of the jacobian approximation which
215    // correspond to low-sensitivity parameters. It can also reduce the
216    // robustness of the solution to errors in the jacobians.
217    //
218    // Oren S.S., Self-scaling variable metric (SSVM) algorithms
219    // Part II: Implementation and experiments, Management Science,
220    // 20(5), 863-874, 1974.
221    bool use_approximate_eigenvalue_bfgs_scaling;
222
223    // Degree of the polynomial used to approximate the objective
224    // function. Valid values are BISECTION, QUADRATIC and CUBIC.
225    //
226    // BISECTION corresponds to pure backtracking search with no
227    // interpolation.
228    LineSearchInterpolationType line_search_interpolation_type;
229
230    // If during the line search, the step_size falls below this
231    // value, it is truncated to zero.
232    double min_line_search_step_size;
233
234    // Line search parameters.
235
236    // Solving the line search problem exactly is computationally
237    // prohibitive. Fortunately, line search based optimization
238    // algorithms can still guarantee convergence if instead of an
239    // exact solution, the line search algorithm returns a solution
240    // which decreases the value of the objective function
241    // sufficiently. More precisely, we are looking for a step_size
242    // s.t.
243    //
244    //   f(step_size) <= f(0) + sufficient_decrease * f'(0) * step_size
245    //
246    double line_search_sufficient_function_decrease;
247
248    // In each iteration of the line search,
249    //
250    //  new_step_size >= max_line_search_step_contraction * step_size
251    //
252    // Note that by definition, for contraction:
253    //
254    //  0 < max_step_contraction < min_step_contraction < 1
255    //
256    double max_line_search_step_contraction;
257
258    // In each iteration of the line search,
259    //
260    //  new_step_size <= min_line_search_step_contraction * step_size
261    //
262    // Note that by definition, for contraction:
263    //
264    //  0 < max_step_contraction < min_step_contraction < 1
265    //
266    double min_line_search_step_contraction;
267
268    // Maximum number of trial step size iterations during each line search,
269    // if a step size satisfying the search conditions cannot be found within
270    // this number of trials, the line search will terminate.
271    int max_num_line_search_step_size_iterations;
272
273    // Maximum number of restarts of the line search direction algorithm before
274    // terminating the optimization. Restarts of the line search direction
275    // algorithm occur when the current algorithm fails to produce a new descent
276    // direction. This typically indicates a numerical failure, or a breakdown
277    // in the validity of the approximations used.
278    int max_num_line_search_direction_restarts;
279
280    // The strong Wolfe conditions consist of the Armijo sufficient
281    // decrease condition, and an additional requirement that the
282    // step-size be chosen s.t. the _magnitude_ ('strong' Wolfe
283    // conditions) of the gradient along the search direction
284    // decreases sufficiently. Precisely, this second condition
285    // is that we seek a step_size s.t.
286    //
287    //   |f'(step_size)| <= sufficient_curvature_decrease * |f'(0)|
288    //
289    // Where f() is the line search objective and f'() is the derivative
290    // of f w.r.t step_size (d f / d step_size).
291    double line_search_sufficient_curvature_decrease;
292
293    // During the bracketing phase of the Wolfe search, the step size is
294    // increased until either a point satisfying the Wolfe conditions is
295    // found, or an upper bound for a bracket containing a point satisfying
296    // the conditions is found.  Precisely, at each iteration of the
297    // expansion:
298    //
299    //   new_step_size <= max_step_expansion * step_size.
300    //
301    // By definition for expansion, max_step_expansion > 1.0.
302    double max_line_search_step_expansion;
303
304    TrustRegionStrategyType trust_region_strategy_type;
305
306    // Type of dogleg strategy to use.
307    DoglegType dogleg_type;
308
309    // The classical trust region methods are descent methods, in that
310    // they only accept a point if it strictly reduces the value of
311    // the objective function.
312    //
313    // Relaxing this requirement allows the algorithm to be more
314    // efficient in the long term at the cost of some local increase
315    // in the value of the objective function.
316    //
317    // This is because allowing for non-decreasing objective function
318    // values in a princpled manner allows the algorithm to "jump over
319    // boulders" as the method is not restricted to move into narrow
320    // valleys while preserving its convergence properties.
321    //
322    // Setting use_nonmonotonic_steps to true enables the
323    // non-monotonic trust region algorithm as described by Conn,
324    // Gould & Toint in "Trust Region Methods", Section 10.1.
325    //
326    // The parameter max_consecutive_nonmonotonic_steps controls the
327    // window size used by the step selection algorithm to accept
328    // non-monotonic steps.
329    //
330    // Even though the value of the objective function may be larger
331    // than the minimum value encountered over the course of the
332    // optimization, the final parameters returned to the user are the
333    // ones corresponding to the minimum cost over all iterations.
334    bool use_nonmonotonic_steps;
335    int max_consecutive_nonmonotonic_steps;
336
337    // Maximum number of iterations for the minimizer to run for.
338    int max_num_iterations;
339
340    // Maximum time for which the minimizer should run for.
341    double max_solver_time_in_seconds;
342
343    // Number of threads used by Ceres for evaluating the cost and
344    // jacobians.
345    int num_threads;
346
347    // Trust region minimizer settings.
348    double initial_trust_region_radius;
349    double max_trust_region_radius;
350
351    // Minimizer terminates when the trust region radius becomes
352    // smaller than this value.
353    double min_trust_region_radius;
354
355    // Lower bound for the relative decrease before a step is
356    // accepted.
357    double min_relative_decrease;
358
359    // For the Levenberg-Marquadt algorithm, the scaled diagonal of
360    // the normal equations J'J is used to control the size of the
361    // trust region. Extremely small and large values along the
362    // diagonal can make this regularization scheme
363    // fail. max_lm_diagonal and min_lm_diagonal, clamp the values of
364    // diag(J'J) from above and below. In the normal course of
365    // operation, the user should not have to modify these parameters.
366    double min_lm_diagonal;
367    double max_lm_diagonal;
368
369    // Sometimes due to numerical conditioning problems or linear
370    // solver flakiness, the trust region strategy may return a
371    // numerically invalid step that can be fixed by reducing the
372    // trust region size. So the TrustRegionMinimizer allows for a few
373    // successive invalid steps before it declares NUMERICAL_FAILURE.
374    int max_num_consecutive_invalid_steps;
375
376    // Minimizer terminates when
377    //
378    //   (new_cost - old_cost) < function_tolerance * old_cost;
379    //
380    double function_tolerance;
381
382    // Minimizer terminates when
383    //
384    //   max_i |x - Project(Plus(x, -g(x))| < gradient_tolerance
385    //
386    // This value should typically be 1e-4 * function_tolerance.
387    double gradient_tolerance;
388
389    // Minimizer terminates when
390    //
391    //   |step|_2 <= parameter_tolerance * ( |x|_2 +  parameter_tolerance)
392    //
393    double parameter_tolerance;
394
395    // Linear least squares solver options -------------------------------------
396
397    LinearSolverType linear_solver_type;
398
399    // Type of preconditioner to use with the iterative linear solvers.
400    PreconditionerType preconditioner_type;
401
402    // Type of clustering algorithm to use for visibility based
403    // preconditioning. This option is used only when the
404    // preconditioner_type is CLUSTER_JACOBI or CLUSTER_TRIDIAGONAL.
405    VisibilityClusteringType visibility_clustering_type;
406
407    // Ceres supports using multiple dense linear algebra libraries
408    // for dense matrix factorizations. Currently EIGEN and LAPACK are
409    // the valid choices. EIGEN is always available, LAPACK refers to
410    // the system BLAS + LAPACK library which may or may not be
411    // available.
412    //
413    // This setting affects the DENSE_QR, DENSE_NORMAL_CHOLESKY and
414    // DENSE_SCHUR solvers. For small to moderate sized probem EIGEN
415    // is a fine choice but for large problems, an optimized LAPACK +
416    // BLAS implementation can make a substantial difference in
417    // performance.
418    DenseLinearAlgebraLibraryType dense_linear_algebra_library_type;
419
420    // Ceres supports using multiple sparse linear algebra libraries
421    // for sparse matrix ordering and factorizations. Currently,
422    // SUITE_SPARSE and CX_SPARSE are the valid choices, depending on
423    // whether they are linked into Ceres at build time.
424    SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type;
425
426    // Number of threads used by Ceres to solve the Newton
427    // step. Currently only the SPARSE_SCHUR solver is capable of
428    // using this setting.
429    int num_linear_solver_threads;
430
431    // The order in which variables are eliminated in a linear solver
432    // can have a significant of impact on the efficiency and accuracy
433    // of the method. e.g., when doing sparse Cholesky factorization,
434    // there are matrices for which a good ordering will give a
435    // Cholesky factor with O(n) storage, where as a bad ordering will
436    // result in an completely dense factor.
437    //
438    // Ceres allows the user to provide varying amounts of hints to
439    // the solver about the variable elimination ordering to use. This
440    // can range from no hints, where the solver is free to decide the
441    // best possible ordering based on the user's choices like the
442    // linear solver being used, to an exact order in which the
443    // variables should be eliminated, and a variety of possibilities
444    // in between.
445    //
446    // Instances of the ParameterBlockOrdering class are used to
447    // communicate this information to Ceres.
448    //
449    // Formally an ordering is an ordered partitioning of the
450    // parameter blocks, i.e, each parameter block belongs to exactly
451    // one group, and each group has a unique non-negative integer
452    // associated with it, that determines its order in the set of
453    // groups.
454    //
455    // Given such an ordering, Ceres ensures that the parameter blocks in
456    // the lowest numbered group are eliminated first, and then the
457    // parmeter blocks in the next lowest numbered group and so on. Within
458    // each group, Ceres is free to order the parameter blocks as it
459    // chooses.
460    //
461    // If NULL, then all parameter blocks are assumed to be in the
462    // same group and the solver is free to decide the best
463    // ordering.
464    //
465    // e.g. Consider the linear system
466    //
467    //   x + y = 3
468    //   2x + 3y = 7
469    //
470    // There are two ways in which it can be solved. First eliminating x
471    // from the two equations, solving for y and then back substituting
472    // for x, or first eliminating y, solving for x and back substituting
473    // for y. The user can construct three orderings here.
474    //
475    //   {0: x}, {1: y} - eliminate x first.
476    //   {0: y}, {1: x} - eliminate y first.
477    //   {0: x, y}      - Solver gets to decide the elimination order.
478    //
479    // Thus, to have Ceres determine the ordering automatically using
480    // heuristics, put all the variables in group 0 and to control the
481    // ordering for every variable, create groups 0..N-1, one per
482    // variable, in the desired order.
483    //
484    // Bundle Adjustment
485    // -----------------
486    //
487    // A particular case of interest is bundle adjustment, where the user
488    // has two options. The default is to not specify an ordering at all,
489    // the solver will see that the user wants to use a Schur type solver
490    // and figure out the right elimination ordering.
491    //
492    // But if the user already knows what parameter blocks are points and
493    // what are cameras, they can save preprocessing time by partitioning
494    // the parameter blocks into two groups, one for the points and one
495    // for the cameras, where the group containing the points has an id
496    // smaller than the group containing cameras.
497    shared_ptr<ParameterBlockOrdering> linear_solver_ordering;
498
499    // Sparse Cholesky factorization algorithms use a fill-reducing
500    // ordering to permute the columns of the Jacobian matrix. There
501    // are two ways of doing this.
502
503    // 1. Compute the Jacobian matrix in some order and then have the
504    //    factorization algorithm permute the columns of the Jacobian.
505
506    // 2. Compute the Jacobian with its columns already permuted.
507
508    // The first option incurs a significant memory penalty. The
509    // factorization algorithm has to make a copy of the permuted
510    // Jacobian matrix, thus Ceres pre-permutes the columns of the
511    // Jacobian matrix and generally speaking, there is no performance
512    // penalty for doing so.
513
514    // In some rare cases, it is worth using a more complicated
515    // reordering algorithm which has slightly better runtime
516    // performance at the expense of an extra copy of the Jacobian
517    // matrix. Setting use_postordering to true enables this tradeoff.
518    bool use_postordering;
519
520    // Some non-linear least squares problems are symbolically dense but
521    // numerically sparse. i.e. at any given state only a small number
522    // of jacobian entries are non-zero, but the position and number of
523    // non-zeros is different depending on the state. For these problems
524    // it can be useful to factorize the sparse jacobian at each solver
525    // iteration instead of including all of the zero entries in a single
526    // general factorization.
527    //
528    // If your problem does not have this property (or you do not know),
529    // then it is probably best to keep this false, otherwise it will
530    // likely lead to worse performance.
531
532    // This settings affects the SPARSE_NORMAL_CHOLESKY solver.
533    bool dynamic_sparsity;
534
535    // Some non-linear least squares problems have additional
536    // structure in the way the parameter blocks interact that it is
537    // beneficial to modify the way the trust region step is computed.
538    //
539    // e.g., consider the following regression problem
540    //
541    //   y = a_1 exp(b_1 x) + a_2 exp(b_3 x^2 + c_1)
542    //
543    // Given a set of pairs{(x_i, y_i)}, the user wishes to estimate
544    // a_1, a_2, b_1, b_2, and c_1.
545    //
546    // Notice here that the expression on the left is linear in a_1
547    // and a_2, and given any value for b_1, b_2 and c_1, it is
548    // possible to use linear regression to estimate the optimal
549    // values of a_1 and a_2. Indeed, its possible to analytically
550    // eliminate the variables a_1 and a_2 from the problem all
551    // together. Problems like these are known as separable least
552    // squares problem and the most famous algorithm for solving them
553    // is the Variable Projection algorithm invented by Golub &
554    // Pereyra.
555    //
556    // Similar structure can be found in the matrix factorization with
557    // missing data problem. There the corresponding algorithm is
558    // known as Wiberg's algorithm.
559    //
560    // Ruhe & Wedin (Algorithms for Separable Nonlinear Least Squares
561    // Problems, SIAM Reviews, 22(3), 1980) present an analyis of
562    // various algorithms for solving separable non-linear least
563    // squares problems and refer to "Variable Projection" as
564    // Algorithm I in their paper.
565    //
566    // Implementing Variable Projection is tedious and expensive, and
567    // they present a simpler algorithm, which they refer to as
568    // Algorithm II, where once the Newton/Trust Region step has been
569    // computed for the whole problem (a_1, a_2, b_1, b_2, c_1) and
570    // additional optimization step is performed to estimate a_1 and
571    // a_2 exactly.
572    //
573    // This idea can be generalized to cases where the residual is not
574    // linear in a_1 and a_2, i.e., Solve for the trust region step
575    // for the full problem, and then use it as the starting point to
576    // further optimize just a_1 and a_2. For the linear case, this
577    // amounts to doing a single linear least squares solve. For
578    // non-linear problems, any method for solving the a_1 and a_2
579    // optimization problems will do. The only constraint on a_1 and
580    // a_2 is that they do not co-occur in any residual block.
581    //
582    // This idea can be further generalized, by not just optimizing
583    // (a_1, a_2), but decomposing the graph corresponding to the
584    // Hessian matrix's sparsity structure in a collection of
585    // non-overlapping independent sets and optimizing each of them.
586    //
587    // Setting "use_inner_iterations" to true enables the use of this
588    // non-linear generalization of Ruhe & Wedin's Algorithm II.  This
589    // version of Ceres has a higher iteration complexity, but also
590    // displays better convergence behaviour per iteration. Setting
591    // Solver::Options::num_threads to the maximum number possible is
592    // highly recommended.
593    bool use_inner_iterations;
594
595    // If inner_iterations is true, then the user has two choices.
596    //
597    // 1. Let the solver heuristically decide which parameter blocks
598    //    to optimize in each inner iteration. To do this leave
599    //    Solver::Options::inner_iteration_ordering untouched.
600    //
601    // 2. Specify a collection of of ordered independent sets. Where
602    //    the lower numbered groups are optimized before the higher
603    //    number groups. Each group must be an independent set. Not
604    //    all parameter blocks need to be present in the ordering.
605    shared_ptr<ParameterBlockOrdering> inner_iteration_ordering;
606
607    // Generally speaking, inner iterations make significant progress
608    // in the early stages of the solve and then their contribution
609    // drops down sharply, at which point the time spent doing inner
610    // iterations is not worth it.
611    //
612    // Once the relative decrease in the objective function due to
613    // inner iterations drops below inner_iteration_tolerance, the use
614    // of inner iterations in subsequent trust region minimizer
615    // iterations is disabled.
616    double inner_iteration_tolerance;
617
618    // Minimum number of iterations for which the linear solver should
619    // run, even if the convergence criterion is satisfied.
620    int min_linear_solver_iterations;
621
622    // Maximum number of iterations for which the linear solver should
623    // run. If the solver does not converge in less than
624    // max_linear_solver_iterations, then it returns MAX_ITERATIONS,
625    // as its termination type.
626    int max_linear_solver_iterations;
627
628    // Forcing sequence parameter. The truncated Newton solver uses
629    // this number to control the relative accuracy with which the
630    // Newton step is computed.
631    //
632    // This constant is passed to ConjugateGradientsSolver which uses
633    // it to terminate the iterations when
634    //
635    //  (Q_i - Q_{i-1})/Q_i < eta/i
636    double eta;
637
638    // Normalize the jacobian using Jacobi scaling before calling
639    // the linear least squares solver.
640    bool jacobi_scaling;
641
642    // Logging options ---------------------------------------------------------
643
644    LoggingType logging_type;
645
646    // By default the Minimizer progress is logged to VLOG(1), which
647    // is sent to STDERR depending on the vlog level. If this flag is
648    // set to true, and logging_type is not SILENT, the logging output
649    // is sent to STDOUT.
650    bool minimizer_progress_to_stdout;
651
652    // List of iterations at which the minimizer should dump the trust
653    // region problem. Useful for testing and benchmarking. If empty
654    // (default), no problems are dumped.
655    vector<int> trust_region_minimizer_iterations_to_dump;
656
657    // Directory to which the problems should be written to. Should be
658    // non-empty if trust_region_minimizer_iterations_to_dump is
659    // non-empty and trust_region_problem_dump_format_type is not
660    // CONSOLE.
661    string trust_region_problem_dump_directory;
662    DumpFormatType trust_region_problem_dump_format_type;
663
664    // Finite differences options ----------------------------------------------
665
666    // Check all jacobians computed by each residual block with finite
667    // differences. This is expensive since it involves computing the
668    // derivative by normal means (e.g. user specified, autodiff,
669    // etc), then also computing it using finite differences. The
670    // results are compared, and if they differ substantially, details
671    // are printed to the log.
672    bool check_gradients;
673
674    // Relative precision to check for in the gradient checker. If the
675    // relative difference between an element in a jacobian exceeds
676    // this number, then the jacobian for that cost term is dumped.
677    double gradient_check_relative_precision;
678
679    // Relative shift used for taking numeric derivatives. For finite
680    // differencing, each dimension is evaluated at slightly shifted
681    // values; for the case of central difference, this is what gets
682    // evaluated:
683    //
684    //   delta = numeric_derivative_relative_step_size;
685    //   f_initial  = f(x)
686    //   f_forward  = f((1 + delta) * x)
687    //   f_backward = f((1 - delta) * x)
688    //
689    // The finite differencing is done along each dimension. The
690    // reason to use a relative (rather than absolute) step size is
691    // that this way, numeric differentation works for functions where
692    // the arguments are typically large (e.g. 1e9) and when the
693    // values are small (e.g. 1e-5). It is possible to construct
694    // "torture cases" which break this finite difference heuristic,
695    // but they do not come up often in practice.
696    //
697    // TODO(keir): Pick a smarter number than the default above! In
698    // theory a good choice is sqrt(eps) * x, which for doubles means
699    // about 1e-8 * x. However, I have found this number too
700    // optimistic. This number should be exposed for users to change.
701    double numeric_derivative_relative_step_size;
702
703    // If true, the user's parameter blocks are updated at the end of
704    // every Minimizer iteration, otherwise they are updated when the
705    // Minimizer terminates. This is useful if, for example, the user
706    // wishes to visualize the state of the optimization every
707    // iteration.
708    bool update_state_every_iteration;
709
710    // Callbacks that are executed at the end of each iteration of the
711    // Minimizer. An iteration may terminate midway, either due to
712    // numerical failures or because one of the convergence tests has
713    // been satisfied. In this case none of the callbacks are
714    // executed.
715
716    // Callbacks are executed in the order that they are specified in
717    // this vector. By default, parameter blocks are updated only at
718    // the end of the optimization, i.e when the Minimizer
719    // terminates. This behaviour is controlled by
720    // update_state_every_variable. If the user wishes to have access
721    // to the update parameter blocks when his/her callbacks are
722    // executed, then set update_state_every_iteration to true.
723    //
724    // The solver does NOT take ownership of these pointers.
725    vector<IterationCallback*> callbacks;
726  };
727
728  struct CERES_EXPORT Summary {
729    Summary();
730
731    // A brief one line description of the state of the solver after
732    // termination.
733    string BriefReport() const;
734
735    // A full multiline description of the state of the solver after
736    // termination.
737    string FullReport() const;
738
739    bool IsSolutionUsable() const;
740
741    // Minimizer summary -------------------------------------------------
742    MinimizerType minimizer_type;
743
744    TerminationType termination_type;
745
746    // Reason why the solver terminated.
747    string message;
748
749    // Cost of the problem (value of the objective function) before
750    // the optimization.
751    double initial_cost;
752
753    // Cost of the problem (value of the objective function) after the
754    // optimization.
755    double final_cost;
756
757    // The part of the total cost that comes from residual blocks that
758    // were held fixed by the preprocessor because all the parameter
759    // blocks that they depend on were fixed.
760    double fixed_cost;
761
762    // IterationSummary for each minimizer iteration in order.
763    vector<IterationSummary> iterations;
764
765    // Number of minimizer iterations in which the step was
766    // accepted. Unless use_non_monotonic_steps is true this is also
767    // the number of steps in which the objective function value/cost
768    // went down.
769    int num_successful_steps;
770
771    // Number of minimizer iterations in which the step was rejected
772    // either because it did not reduce the cost enough or the step
773    // was not numerically valid.
774    int num_unsuccessful_steps;
775
776    // Number of times inner iterations were performed.
777    int num_inner_iteration_steps;
778
779    // All times reported below are wall times.
780
781    // When the user calls Solve, before the actual optimization
782    // occurs, Ceres performs a number of preprocessing steps. These
783    // include error checks, memory allocations, and reorderings. This
784    // time is accounted for as preprocessing time.
785    double preprocessor_time_in_seconds;
786
787    // Time spent in the TrustRegionMinimizer.
788    double minimizer_time_in_seconds;
789
790    // After the Minimizer is finished, some time is spent in
791    // re-evaluating residuals etc. This time is accounted for in the
792    // postprocessor time.
793    double postprocessor_time_in_seconds;
794
795    // Some total of all time spent inside Ceres when Solve is called.
796    double total_time_in_seconds;
797
798    // Time (in seconds) spent in the linear solver computing the
799    // trust region step.
800    double linear_solver_time_in_seconds;
801
802    // Time (in seconds) spent evaluating the residual vector.
803    double residual_evaluation_time_in_seconds;
804
805    // Time (in seconds) spent evaluating the jacobian matrix.
806    double jacobian_evaluation_time_in_seconds;
807
808    // Time (in seconds) spent doing inner iterations.
809    double inner_iteration_time_in_seconds;
810
811    // Number of parameter blocks in the problem.
812    int num_parameter_blocks;
813
814    // Number of parameters in the probem.
815    int num_parameters;
816
817    // Dimension of the tangent space of the problem (or the number of
818    // columns in the Jacobian for the problem). This is different
819    // from num_parameters if a parameter block is associated with a
820    // LocalParameterization
821    int num_effective_parameters;
822
823    // Number of residual blocks in the problem.
824    int num_residual_blocks;
825
826    // Number of residuals in the problem.
827    int num_residuals;
828
829    // Number of parameter blocks in the problem after the inactive
830    // and constant parameter blocks have been removed. A parameter
831    // block is inactive if no residual block refers to it.
832    int num_parameter_blocks_reduced;
833
834    // Number of parameters in the reduced problem.
835    int num_parameters_reduced;
836
837    // Dimension of the tangent space of the reduced problem (or the
838    // number of columns in the Jacobian for the reduced
839    // problem). This is different from num_parameters_reduced if a
840    // parameter block in the reduced problem is associated with a
841    // LocalParameterization.
842    int num_effective_parameters_reduced;
843
844    // Number of residual blocks in the reduced problem.
845    int num_residual_blocks_reduced;
846
847    //  Number of residuals in the reduced problem.
848    int num_residuals_reduced;
849
850    //  Number of threads specified by the user for Jacobian and
851    //  residual evaluation.
852    int num_threads_given;
853
854    // Number of threads actually used by the solver for Jacobian and
855    // residual evaluation. This number is not equal to
856    // num_threads_given if OpenMP is not available.
857    int num_threads_used;
858
859    //  Number of threads specified by the user for solving the trust
860    // region problem.
861    int num_linear_solver_threads_given;
862
863    // Number of threads actually used by the solver for solving the
864    // trust region problem. This number is not equal to
865    // num_threads_given if OpenMP is not available.
866    int num_linear_solver_threads_used;
867
868    // Type of the linear solver requested by the user.
869    LinearSolverType linear_solver_type_given;
870
871    // Type of the linear solver actually used. This may be different
872    // from linear_solver_type_given if Ceres determines that the
873    // problem structure is not compatible with the linear solver
874    // requested or if the linear solver requested by the user is not
875    // available, e.g. The user requested SPARSE_NORMAL_CHOLESKY but
876    // no sparse linear algebra library was available.
877    LinearSolverType linear_solver_type_used;
878
879    // Size of the elimination groups given by the user as hints to
880    // the linear solver.
881    vector<int> linear_solver_ordering_given;
882
883    // Size of the parameter groups used by the solver when ordering
884    // the columns of the Jacobian.  This maybe different from
885    // linear_solver_ordering_given if the user left
886    // linear_solver_ordering_given blank and asked for an automatic
887    // ordering, or if the problem contains some constant or inactive
888    // parameter blocks.
889    vector<int> linear_solver_ordering_used;
890
891    // True if the user asked for inner iterations to be used as part
892    // of the optimization.
893    bool inner_iterations_given;
894
895    // True if the user asked for inner iterations to be used as part
896    // of the optimization and the problem structure was such that
897    // they were actually performed. e.g., in a problem with just one
898    // parameter block, inner iterations are not performed.
899    bool inner_iterations_used;
900
901    // Size of the parameter groups given by the user for performing
902    // inner iterations.
903    vector<int> inner_iteration_ordering_given;
904
905    // Size of the parameter groups given used by the solver for
906    // performing inner iterations. This maybe different from
907    // inner_iteration_ordering_given if the user left
908    // inner_iteration_ordering_given blank and asked for an automatic
909    // ordering, or if the problem contains some constant or inactive
910    // parameter blocks.
911    vector<int> inner_iteration_ordering_used;
912
913    //  Type of preconditioner used for solving the trust region
914    //  step. Only meaningful when an iterative linear solver is used.
915    PreconditionerType preconditioner_type;
916
917    // Type of clustering algorithm used for visibility based
918    // preconditioning. Only meaningful when the preconditioner_type
919    // is CLUSTER_JACOBI or CLUSTER_TRIDIAGONAL.
920    VisibilityClusteringType visibility_clustering_type;
921
922    //  Type of trust region strategy.
923    TrustRegionStrategyType trust_region_strategy_type;
924
925    //  Type of dogleg strategy used for solving the trust region
926    //  problem.
927    DoglegType dogleg_type;
928
929    //  Type of the dense linear algebra library used.
930    DenseLinearAlgebraLibraryType dense_linear_algebra_library_type;
931
932    // Type of the sparse linear algebra library used.
933    SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type;
934
935    // Type of line search direction used.
936    LineSearchDirectionType line_search_direction_type;
937
938    // Type of the line search algorithm used.
939    LineSearchType line_search_type;
940
941    //  When performing line search, the degree of the polynomial used
942    //  to approximate the objective function.
943    LineSearchInterpolationType line_search_interpolation_type;
944
945    // If the line search direction is NONLINEAR_CONJUGATE_GRADIENT,
946    // then this indicates the particular variant of non-linear
947    // conjugate gradient used.
948    NonlinearConjugateGradientType nonlinear_conjugate_gradient_type;
949
950    // If the type of the line search direction is LBFGS, then this
951    // indicates the rank of the Hessian approximation.
952    int max_lbfgs_rank;
953  };
954
955  // Once a least squares problem has been built, this function takes
956  // the problem and optimizes it based on the values of the options
957  // parameters. Upon return, a detailed summary of the work performed
958  // by the preprocessor, the non-linear minmizer and the linear
959  // solver are reported in the summary object.
960  virtual void Solve(const Options& options,
961                     Problem* problem,
962                     Solver::Summary* summary);
963};
964
965// Helper function which avoids going through the interface.
966CERES_EXPORT void Solve(const Solver::Options& options,
967           Problem* problem,
968           Solver::Summary* summary);
969
970}  // namespace ceres
971
972#include "ceres/internal/reenable_warnings.h"
973
974#endif  // CERES_PUBLIC_SOLVER_H_
975