Home | History | Annotate | Download | only in ceres
      1 // Ceres Solver - A fast non-linear least squares minimizer
      2 // Copyright 2014 Google Inc. All rights reserved.
      3 // http://code.google.com/p/ceres-solver/
      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 (at) 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.
     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): 773782.
    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.
    966 CERES_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