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