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