| // Ceres Solver - A fast non-linear least squares minimizer |
| // Copyright 2022 Google Inc. All rights reserved. |
| // http://ceres-solver.org/ |
| // |
| // Redistribution and use in source and binary forms, with or without |
| // modification, are permitted provided that the following conditions are met: |
| // |
| // * Redistributions of source code must retain the above copyright notice, |
| // this list of conditions and the following disclaimer. |
| // * Redistributions in binary form must reproduce the above copyright notice, |
| // this list of conditions and the following disclaimer in the documentation |
| // and/or other materials provided with the distribution. |
| // * Neither the name of Google Inc. nor the names of its contributors may be |
| // used to endorse or promote products derived from this software without |
| // specific prior written permission. |
| // |
| // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" |
| // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
| // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
| // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE |
| // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR |
| // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF |
| // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS |
| // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN |
| // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) |
| // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
| // POSSIBILITY OF SUCH DAMAGE. |
| // |
| // Author: sameeragarwal@google.com (Sameer Agarwal) |
| |
| #ifndef CERES_PUBLIC_SOLVER_H_ |
| #define CERES_PUBLIC_SOLVER_H_ |
| |
| #include <cmath> |
| #include <memory> |
| #include <string> |
| #include <unordered_set> |
| #include <vector> |
| |
| #include "ceres/crs_matrix.h" |
| #include "ceres/internal/config.h" |
| #include "ceres/internal/disable_warnings.h" |
| #include "ceres/internal/export.h" |
| #include "ceres/iteration_callback.h" |
| #include "ceres/ordered_groups.h" |
| #include "ceres/problem.h" |
| #include "ceres/types.h" |
| |
| namespace ceres { |
| |
| // Interface for non-linear least squares solvers. |
| class CERES_EXPORT Solver { |
| public: |
| virtual ~Solver(); |
| |
| // The options structure contains, not surprisingly, options that control how |
| // the solver operates. The defaults should be suitable for a wide range of |
| // problems; however, better performance is often obtainable with tweaking. |
| // |
| // The constants are defined inside types.h |
| struct CERES_EXPORT Options { |
| // Returns true if the options struct has a valid |
| // configuration. Returns false otherwise, and fills in *error |
| // with a message describing the problem. |
| bool IsValid(std::string* error) const; |
| |
| // Minimizer options ---------------------------------------- |
| |
| // Ceres supports the two major families of optimization strategies - |
| // Trust Region and Line Search. |
| // |
| // 1. The line search approach first finds a descent direction |
| // along which the objective function will be reduced and then |
| // computes a step size that decides how far should move along |
| // that direction. The descent direction can be computed by |
| // various methods, such as gradient descent, Newton's method and |
| // Quasi-Newton method. The step size can be determined either |
| // exactly or inexactly. |
| // |
| // 2. The trust region approach approximates the objective |
| // function using a model function (often a quadratic) over |
| // a subset of the search space known as the trust region. If the |
| // model function succeeds in minimizing the true objective |
| // function the trust region is expanded; conversely, otherwise it |
| // is contracted and the model optimization problem is solved |
| // again. |
| // |
| // Trust region methods are in some sense dual to line search methods: |
| // trust region methods first choose a step size (the size of the |
| // trust region) and then a step direction while line search methods |
| // first choose a step direction and then a step size. |
| MinimizerType minimizer_type = TRUST_REGION; |
| |
| LineSearchDirectionType line_search_direction_type = LBFGS; |
| LineSearchType line_search_type = WOLFE; |
| NonlinearConjugateGradientType nonlinear_conjugate_gradient_type = |
| FLETCHER_REEVES; |
| |
| // The LBFGS hessian approximation is a low rank approximation to |
| // the inverse of the Hessian matrix. The rank of the |
| // approximation determines (linearly) the space and time |
| // complexity of using the approximation. Higher the rank, the |
| // better is the quality of the approximation. The increase in |
| // quality is however is bounded for a number of reasons. |
| // |
| // 1. The method only uses secant information and not actual |
| // derivatives. |
| // |
| // 2. The Hessian approximation is constrained to be positive |
| // definite. |
| // |
| // So increasing this rank to a large number will cost time and |
| // space complexity without the corresponding increase in solution |
| // quality. There are no hard and fast rules for choosing the |
| // maximum rank. The best choice usually requires some problem |
| // specific experimentation. |
| // |
| // For more theoretical and implementation details of the LBFGS |
| // method, please see: |
| // |
| // Nocedal, J. (1980). "Updating Quasi-Newton Matrices with |
| // Limited Storage". Mathematics of Computation 35 (151): 773-782. |
| int max_lbfgs_rank = 20; |
| |
| // As part of the (L)BFGS update step (BFGS) / right-multiply step (L-BFGS), |
| // the initial inverse Hessian approximation is taken to be the Identity. |
| // However, Oren showed that using instead I * \gamma, where \gamma is |
| // chosen to approximate an eigenvalue of the true inverse Hessian can |
| // result in improved convergence in a wide variety of cases. Setting |
| // use_approximate_eigenvalue_bfgs_scaling to true enables this scaling. |
| // |
| // It is important to note that approximate eigenvalue scaling does not |
| // always improve convergence, and that it can in fact significantly degrade |
| // performance for certain classes of problem, which is why it is disabled |
| // by default. In particular it can degrade performance when the |
| // sensitivity of the problem to different parameters varies significantly, |
| // as in this case a single scalar factor fails to capture this variation |
| // and detrimentally downscales parts of the jacobian approximation which |
| // correspond to low-sensitivity parameters. It can also reduce the |
| // robustness of the solution to errors in the jacobians. |
| // |
| // Oren S.S., Self-scaling variable metric (SSVM) algorithms |
| // Part II: Implementation and experiments, Management Science, |
| // 20(5), 863-874, 1974. |
| bool use_approximate_eigenvalue_bfgs_scaling = false; |
| |
| // Degree of the polynomial used to approximate the objective |
| // function. Valid values are BISECTION, QUADRATIC and CUBIC. |
| // |
| // BISECTION corresponds to pure backtracking search with no |
| // interpolation. |
| LineSearchInterpolationType line_search_interpolation_type = CUBIC; |
| |
| // If during the line search, the step_size falls below this |
| // value, it is truncated to zero. |
| double min_line_search_step_size = 1e-9; |
| |
| // Line search parameters. |
| |
| // Solving the line search problem exactly is computationally |
| // prohibitive. Fortunately, line search based optimization |
| // algorithms can still guarantee convergence if instead of an |
| // exact solution, the line search algorithm returns a solution |
| // which decreases the value of the objective function |
| // sufficiently. More precisely, we are looking for a step_size |
| // s.t. |
| // |
| // f(step_size) <= f(0) + sufficient_decrease * f'(0) * step_size |
| // |
| double line_search_sufficient_function_decrease = 1e-4; |
| |
| // In each iteration of the line search, |
| // |
| // new_step_size >= max_line_search_step_contraction * step_size |
| // |
| // Note that by definition, for contraction: |
| // |
| // 0 < max_step_contraction < min_step_contraction < 1 |
| // |
| double max_line_search_step_contraction = 1e-3; |
| |
| // In each iteration of the line search, |
| // |
| // new_step_size <= min_line_search_step_contraction * step_size |
| // |
| // Note that by definition, for contraction: |
| // |
| // 0 < max_step_contraction < min_step_contraction < 1 |
| // |
| double min_line_search_step_contraction = 0.6; |
| |
| // Maximum number of trial step size iterations during each line |
| // search, if a step size satisfying the search conditions cannot |
| // be found within this number of trials, the line search will |
| // terminate. |
| |
| // The minimum allowed value is 0 for trust region minimizer and 1 |
| // otherwise. If 0 is specified for the trust region minimizer, |
| // then line search will not be used when solving constrained |
| // optimization problems. |
| int max_num_line_search_step_size_iterations = 20; |
| |
| // Maximum number of restarts of the line search direction algorithm before |
| // terminating the optimization. Restarts of the line search direction |
| // algorithm occur when the current algorithm fails to produce a new descent |
| // direction. This typically indicates a numerical failure, or a breakdown |
| // in the validity of the approximations used. |
| int max_num_line_search_direction_restarts = 5; |
| |
| // The strong Wolfe conditions consist of the Armijo sufficient |
| // decrease condition, and an additional requirement that the |
| // step-size be chosen s.t. the _magnitude_ ('strong' Wolfe |
| // conditions) of the gradient along the search direction |
| // decreases sufficiently. Precisely, this second condition |
| // is that we seek a step_size s.t. |
| // |
| // |f'(step_size)| <= sufficient_curvature_decrease * |f'(0)| |
| // |
| // Where f() is the line search objective and f'() is the derivative |
| // of f w.r.t step_size (d f / d step_size). |
| double line_search_sufficient_curvature_decrease = 0.9; |
| |
| // During the bracketing phase of the Wolfe search, the step size is |
| // increased until either a point satisfying the Wolfe conditions is |
| // found, or an upper bound for a bracket containing a point satisfying |
| // the conditions is found. Precisely, at each iteration of the |
| // expansion: |
| // |
| // new_step_size <= max_step_expansion * step_size. |
| // |
| // By definition for expansion, max_step_expansion > 1.0. |
| double max_line_search_step_expansion = 10.0; |
| |
| TrustRegionStrategyType trust_region_strategy_type = LEVENBERG_MARQUARDT; |
| |
| // Type of dogleg strategy to use. |
| DoglegType dogleg_type = TRADITIONAL_DOGLEG; |
| |
| // The classical trust region methods are descent methods, in that |
| // they only accept a point if it strictly reduces the value of |
| // the objective function. |
| // |
| // Relaxing this requirement allows the algorithm to be more |
| // efficient in the long term at the cost of some local increase |
| // in the value of the objective function. |
| // |
| // This is because allowing for non-decreasing objective function |
| // values in a principled manner allows the algorithm to "jump over |
| // boulders" as the method is not restricted to move into narrow |
| // valleys while preserving its convergence properties. |
| // |
| // Setting use_nonmonotonic_steps to true enables the |
| // non-monotonic trust region algorithm as described by Conn, |
| // Gould & Toint in "Trust Region Methods", Section 10.1. |
| // |
| // The parameter max_consecutive_nonmonotonic_steps controls the |
| // window size used by the step selection algorithm to accept |
| // non-monotonic steps. |
| // |
| // Even though the value of the objective function may be larger |
| // than the minimum value encountered over the course of the |
| // optimization, the final parameters returned to the user are the |
| // ones corresponding to the minimum cost over all iterations. |
| bool use_nonmonotonic_steps = false; |
| int max_consecutive_nonmonotonic_steps = 5; |
| |
| // Maximum number of iterations for the minimizer to run for. |
| int max_num_iterations = 50; |
| |
| // Maximum time for which the minimizer should run for. |
| double max_solver_time_in_seconds = 1e9; |
| |
| // Number of threads used by Ceres for evaluating the cost and |
| // jacobians. |
| int num_threads = 1; |
| |
| // Trust region minimizer settings. |
| double initial_trust_region_radius = 1e4; |
| double max_trust_region_radius = 1e16; |
| |
| // Minimizer terminates when the trust region radius becomes |
| // smaller than this value. |
| double min_trust_region_radius = 1e-32; |
| |
| // Lower bound for the relative decrease before a step is |
| // accepted. |
| double min_relative_decrease = 1e-3; |
| |
| // For the Levenberg-Marquadt algorithm, the scaled diagonal of |
| // the normal equations J'J is used to control the size of the |
| // trust region. Extremely small and large values along the |
| // diagonal can make this regularization scheme |
| // fail. max_lm_diagonal and min_lm_diagonal, clamp the values of |
| // diag(J'J) from above and below. In the normal course of |
| // operation, the user should not have to modify these parameters. |
| double min_lm_diagonal = 1e-6; |
| double max_lm_diagonal = 1e32; |
| |
| // Sometimes due to numerical conditioning problems or linear |
| // solver flakiness, the trust region strategy may return a |
| // numerically invalid step that can be fixed by reducing the |
| // trust region size. So the TrustRegionMinimizer allows for a few |
| // successive invalid steps before it declares NUMERICAL_FAILURE. |
| int max_num_consecutive_invalid_steps = 5; |
| |
| // Minimizer terminates when |
| // |
| // (new_cost - old_cost) < function_tolerance * old_cost; |
| // |
| double function_tolerance = 1e-6; |
| |
| // Minimizer terminates when |
| // |
| // max_i |x - Project(Plus(x, -g(x))| < gradient_tolerance |
| // |
| // This value should typically be 1e-4 * function_tolerance. |
| double gradient_tolerance = 1e-10; |
| |
| // Minimizer terminates when |
| // |
| // |step|_2 <= parameter_tolerance * ( |x|_2 + parameter_tolerance) |
| // |
| double parameter_tolerance = 1e-8; |
| |
| // Linear least squares solver options ------------------------------------- |
| |
| LinearSolverType linear_solver_type = |
| #if defined(CERES_NO_SPARSE) |
| DENSE_QR; |
| #else |
| SPARSE_NORMAL_CHOLESKY; |
| #endif |
| |
| // Type of preconditioner to use with the iterative linear solvers. |
| PreconditionerType preconditioner_type = JACOBI; |
| |
| // Type of clustering algorithm to use for visibility based |
| // preconditioning. This option is used only when the |
| // preconditioner_type is CLUSTER_JACOBI or CLUSTER_TRIDIAGONAL. |
| VisibilityClusteringType visibility_clustering_type = CANONICAL_VIEWS; |
| |
| // Subset preconditioner is a preconditioner for problems with |
| // general sparsity. Given a subset of residual blocks of a |
| // problem, it uses the corresponding subset of the rows of the |
| // Jacobian to construct a preconditioner. |
| // |
| // Suppose the Jacobian J has been horizontally partitioned as |
| // |
| // J = [P] |
| // [Q] |
| // |
| // Where, Q is the set of rows corresponding to the residual |
| // blocks in residual_blocks_for_subset_preconditioner. |
| // |
| // The preconditioner is the inverse of the matrix Q'Q. |
| // |
| // Obviously, the efficacy of the preconditioner depends on how |
| // well the matrix Q approximates J'J, or how well the chosen |
| // residual blocks approximate the non-linear least squares |
| // problem. |
| // |
| // If Solver::Options::preconditioner_type == SUBSET, then |
| // residual_blocks_for_subset_preconditioner must be non-empty. |
| std::unordered_set<ResidualBlockId> |
| residual_blocks_for_subset_preconditioner; |
| |
| // Ceres supports using multiple dense linear algebra libraries for dense |
| // matrix factorizations. Currently EIGEN, LAPACK and CUDA are the valid |
| // choices. EIGEN is always available, LAPACK refers to the system BLAS + |
| // LAPACK library which may or may not be available. CUDA refers to Nvidia's |
| // GPU based dense linear algebra library, which may or may not be |
| // available. |
| // |
| // This setting affects the DENSE_QR, DENSE_NORMAL_CHOLESKY and DENSE_SCHUR |
| // solvers. For small to moderate sized problem EIGEN is a fine choice but |
| // for large problems, an optimized LAPACK + BLAS or CUDA implementation can |
| // make a substantial difference in performance. |
| DenseLinearAlgebraLibraryType dense_linear_algebra_library_type = EIGEN; |
| |
| // Ceres supports using multiple sparse linear algebra libraries for sparse |
| // matrix ordering and factorizations. |
| SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type = |
| #if !defined(CERES_NO_SUITESPARSE) |
| SUITE_SPARSE; |
| #elif !defined(CERES_NO_ACCELERATE_SPARSE) |
| ACCELERATE_SPARSE; |
| #elif defined(CERES_USE_EIGEN_SPARSE) |
| EIGEN_SPARSE; |
| #else |
| NO_SPARSE; |
| #endif |
| |
| // The order in which variables are eliminated in a linear solver |
| // can have a significant impact on the efficiency and accuracy of |
| // the method. e.g., when doing sparse Cholesky factorization, |
| // there are matrices for which a good ordering will give a |
| // Cholesky factor with O(n) storage, where as a bad ordering will |
| // result in an completely dense factor. |
| // |
| // Sparse direct solvers like SPARSE_NORMAL_CHOLESKY and |
| // SPARSE_SCHUR use a fill reducing ordering of the columns and |
| // rows of the matrix being factorized before computing the |
| // numeric factorization. |
| // |
| // This enum controls the type of algorithm used to compute |
| // this fill reducing ordering. There is no single algorithm |
| // that works on all matrices, so determining which algorithm |
| // works better is a matter of empirical experimentation. |
| // |
| // The exact behaviour of this setting is affected by the value of |
| // linear_solver_ordering as described below. |
| LinearSolverOrderingType linear_solver_ordering_type = AMD; |
| |
| // Besides specifying the fill reducing ordering via |
| // linear_solver_ordering_type, Ceres allows the user to provide varying |
| // amounts of hints to the linear solver about the variable elimination |
| // ordering to use. This can range from no hints, where the solver is free |
| // to decide the best possible ordering based on the user's choices like the |
| // linear solver being used, to an exact order in which the variables should |
| // be eliminated, and a variety of possibilities in between. |
| // |
| // Instances of the ParameterBlockOrdering class are used to communicate |
| // this information to Ceres. |
| // |
| // Formally an ordering is an ordered partitioning of the parameter blocks, |
| // i.e, each parameter block belongs to exactly one group, and each group |
| // has a unique non-negative integer associated with it, that determines its |
| // order in the set of groups. |
| // |
| // e.g. Consider the linear system |
| // |
| // x + y = 3 |
| // 2x + 3y = 7 |
| // |
| // There are two ways in which it can be solved. First eliminating x from |
| // the two equations, solving for y and then back substituting for x, or |
| // first eliminating y, solving for x and back substituting for y. The user |
| // can construct three orderings here. |
| // |
| // {0: x}, {1: y} - eliminate x first. |
| // {0: y}, {1: x} - eliminate y first. |
| // {0: x, y} - Solver gets to decide the elimination order. |
| // |
| // Thus, to have Ceres determine the ordering automatically, put all the |
| // variables in group 0 and to control the ordering for every variable |
| // create groups 0 ... N-1, one per variable, in the desired |
| // order. |
| // |
| // linear_solver_ordering == nullptr and an ordering where all the parameter |
| // blocks are in one elimination group mean the same thing - the solver is |
| // free to choose what it thinks is the best elimination ordering. Therefore |
| // in the following we will only consider the case where |
| // linear_solver_ordering is nullptr. |
| // |
| // The exact interpretation of this information depends on the values of |
| // linear_solver_ordering_type and linear_solver_type/preconditioner_type |
| // and sparse_linear_algebra_type. |
| // |
| // Bundle Adjustment |
| // ================= |
| // |
| // If the user is using one of the Schur solvers (DENSE_SCHUR, |
| // SPARSE_SCHUR, ITERATIVE_SCHUR) and chooses to specify an |
| // ordering, it must have one important property. The lowest |
| // numbered elimination group must form an independent set in the |
| // graph corresponding to the Hessian, or in other words, no two |
| // parameter blocks in in the first elimination group should |
| // co-occur in the same residual block. For the best performance, |
| // this elimination group should be as large as possible. For |
| // standard bundle adjustment problems, this corresponds to the |
| // first elimination group containing all the 3d points, and the |
| // second containing the all the cameras parameter blocks. |
| // |
| // If the user leaves the choice to Ceres, then the solver uses an |
| // approximate maximum independent set algorithm to identify the first |
| // elimination group. |
| // |
| // sparse_linear_algebra_library_type = SUITE_SPARSE |
| // ================================================= |
| // |
| // linear_solver_ordering_type = AMD |
| // --------------------------------- |
| // |
| // A Constrained Approximate Minimum Degree (CAMD) ordering used where the |
| // parameter blocks in the lowest numbered group are eliminated first, and |
| // then the parameter blocks in the next lowest numbered group and so |
| // on. Within each group, CAMD free to order the parameter blocks as it |
| // chooses. |
| // |
| // linear_solver_ordering_type = NESDIS |
| // ------------------------------------- |
| // |
| // a. linear_solver_type = SPARSE_NORMAL_CHOLESKY or |
| // linear_solver_type = CGNR and preconditioner_type = SUBSET |
| // |
| // The value of linear_solver_ordering is ignored and a Nested Dissection |
| // algorithm is used to compute a fill reducing ordering. |
| // |
| // b. linear_solver_type = SPARSE_SCHUR/DENSE_SCHUR/ITERATIVE_SCHUR |
| // |
| // ONLY the lowest group are used to compute the Schur complement, and |
| // Nested Dissection is used to compute a fill reducing ordering for the |
| // Schur Complement (or its preconditioner). |
| // |
| // sparse_linear_algebra_library_type = EIGEN_SPARSE or ACCELERATE_SPARSE |
| // ====================================================================== |
| // |
| // a. linear_solver_type = SPARSE_NORMAL_CHOLESKY or |
| // linear_solver_type = CGNR and preconditioner_type = SUBSET |
| // |
| // then the value of linear_solver_ordering is ignored and AMD or NESDIS is |
| // used to compute a fill reducing ordering as requested by the user. |
| // |
| // b. linear_solver_type = SPARSE_SCHUR/DENSE_SCHUR/ITERATIVE_SCHUR |
| // |
| // ONLY the lowest group are used to compute the Schur complement, and AMD |
| // or NESDIS is used to compute a fill reducing ordering for the Schur |
| // Complement (or its preconditioner). |
| std::shared_ptr<ParameterBlockOrdering> linear_solver_ordering; |
| |
| // Use an explicitly computed Schur complement matrix with |
| // ITERATIVE_SCHUR. |
| // |
| // By default this option is disabled and ITERATIVE_SCHUR |
| // evaluates matrix-vector products between the Schur |
| // complement and a vector implicitly by exploiting the algebraic |
| // expression for the Schur complement. |
| // |
| // The cost of this evaluation scales with the number of non-zeros |
| // in the Jacobian. |
| // |
| // For small to medium sized problems there is a sweet spot where |
| // computing the Schur complement is cheap enough that it is much |
| // more efficient to explicitly compute it and use it for evaluating |
| // the matrix-vector products. |
| // |
| // Enabling this option tells ITERATIVE_SCHUR to use an explicitly |
| // computed Schur complement. |
| // |
| // NOTE: This option can only be used with the SCHUR_JACOBI |
| // preconditioner. |
| bool use_explicit_schur_complement = false; |
| |
| // Sparse Cholesky factorization algorithms use a fill-reducing |
| // ordering to permute the columns of the Jacobian matrix. There |
| // are two ways of doing this. |
| |
| // 1. Compute the Jacobian matrix in some order and then have the |
| // factorization algorithm permute the columns of the Jacobian. |
| |
| // 2. Compute the Jacobian with its columns already permuted. |
| |
| // The first option incurs a significant memory penalty. The |
| // factorization algorithm has to make a copy of the permuted |
| // Jacobian matrix, thus Ceres pre-permutes the columns of the |
| // Jacobian matrix and generally speaking, there is no performance |
| // penalty for doing so. |
| |
| // Some non-linear least squares problems are symbolically dense but |
| // numerically sparse. i.e. at any given state only a small number |
| // of jacobian entries are non-zero, but the position and number of |
| // non-zeros is different depending on the state. For these problems |
| // it can be useful to factorize the sparse jacobian at each solver |
| // iteration instead of including all of the zero entries in a single |
| // general factorization. |
| // |
| // If your problem does not have this property (or you do not know), |
| // then it is probably best to keep this false, otherwise it will |
| // likely lead to worse performance. |
| |
| // This settings only affects the SPARSE_NORMAL_CHOLESKY solver. |
| bool dynamic_sparsity = false; |
| |
| // TODO(sameeragarwal): Further expand the documentation for the |
| // following two options. |
| // TODO(joydeepbiswas): Update the documentation for the mixed precision |
| // option with CUDA. |
| |
| // NOTE1: EXPERIMENTAL FEATURE, UNDER DEVELOPMENT, USE AT YOUR OWN RISK. |
| // |
| // If use_mixed_precision_solves is true, the Gauss-Newton matrix |
| // is computed in double precision, but its factorization is |
| // computed in single precision. This can result in significant |
| // time and memory savings at the cost of some accuracy in the |
| // Gauss-Newton step. Iterative refinement is used to recover some |
| // of this accuracy back. |
| // |
| // If use_mixed_precision_solves is true, we recommend setting |
| // max_num_refinement_iterations to 2-3. |
| // |
| // NOTE2: The following two options are currently only applicable |
| // if sparse_linear_algebra_library_type is EIGEN_SPARSE or |
| // ACCELERATE_SPARSE, and linear_solver_type is SPARSE_NORMAL_CHOLESKY |
| // or SPARSE_SCHUR. |
| bool use_mixed_precision_solves = false; |
| |
| // Number steps of the iterative refinement process to run when |
| // computing the Gauss-Newton step. |
| int max_num_refinement_iterations = 0; |
| |
| // Some non-linear least squares problems have additional |
| // structure in the way the parameter blocks interact that it is |
| // beneficial to modify the way the trust region step is computed. |
| // |
| // e.g., consider the following regression problem |
| // |
| // y = a_1 exp(b_1 x) + a_2 exp(b_3 x^2 + c_1) |
| // |
| // Given a set of pairs{(x_i, y_i)}, the user wishes to estimate |
| // a_1, a_2, b_1, b_2, and c_1. |
| // |
| // Notice here that the expression on the left is linear in a_1 |
| // and a_2, and given any value for b_1, b_2 and c_1, it is |
| // possible to use linear regression to estimate the optimal |
| // values of a_1 and a_2. Indeed, its possible to analytically |
| // eliminate the variables a_1 and a_2 from the problem all |
| // together. Problems like these are known as separable least |
| // squares problem and the most famous algorithm for solving them |
| // is the Variable Projection algorithm invented by Golub & |
| // Pereyra. |
| // |
| // Similar structure can be found in the matrix factorization with |
| // missing data problem. There the corresponding algorithm is |
| // known as Wiberg's algorithm. |
| // |
| // Ruhe & Wedin (Algorithms for Separable Nonlinear Least Squares |
| // Problems, SIAM Reviews, 22(3), 1980) present an analysis of |
| // various algorithms for solving separable non-linear least |
| // squares problems and refer to "Variable Projection" as |
| // Algorithm I in their paper. |
| // |
| // Implementing Variable Projection is tedious and expensive, and |
| // they present a simpler algorithm, which they refer to as |
| // Algorithm II, where once the Newton/Trust Region step has been |
| // computed for the whole problem (a_1, a_2, b_1, b_2, c_1) and |
| // additional optimization step is performed to estimate a_1 and |
| // a_2 exactly. |
| // |
| // This idea can be generalized to cases where the residual is not |
| // linear in a_1 and a_2, i.e., Solve for the trust region step |
| // for the full problem, and then use it as the starting point to |
| // further optimize just a_1 and a_2. For the linear case, this |
| // amounts to doing a single linear least squares solve. For |
| // non-linear problems, any method for solving the a_1 and a_2 |
| // optimization problems will do. The only constraint on a_1 and |
| // a_2 is that they do not co-occur in any residual block. |
| // |
| // This idea can be further generalized, by not just optimizing |
| // (a_1, a_2), but decomposing the graph corresponding to the |
| // Hessian matrix's sparsity structure in a collection of |
| // non-overlapping independent sets and optimizing each of them. |
| // |
| // Setting "use_inner_iterations" to true enables the use of this |
| // non-linear generalization of Ruhe & Wedin's Algorithm II. This |
| // version of Ceres has a higher iteration complexity, but also |
| // displays better convergence behaviour per iteration. Setting |
| // Solver::Options::num_threads to the maximum number possible is |
| // highly recommended. |
| bool use_inner_iterations = false; |
| |
| // If inner_iterations is true, then the user has two choices. |
| // |
| // 1. Let the solver heuristically decide which parameter blocks |
| // to optimize in each inner iteration. To do this leave |
| // Solver::Options::inner_iteration_ordering untouched. |
| // |
| // 2. Specify a collection of of ordered independent sets. Where |
| // the lower numbered groups are optimized before the higher |
| // number groups. Each group must be an independent set. Not |
| // all parameter blocks need to be present in the ordering. |
| std::shared_ptr<ParameterBlockOrdering> inner_iteration_ordering; |
| |
| // Generally speaking, inner iterations make significant progress |
| // in the early stages of the solve and then their contribution |
| // drops down sharply, at which point the time spent doing inner |
| // iterations is not worth it. |
| // |
| // Once the relative decrease in the objective function due to |
| // inner iterations drops below inner_iteration_tolerance, the use |
| // of inner iterations in subsequent trust region minimizer |
| // iterations is disabled. |
| double inner_iteration_tolerance = 1e-3; |
| |
| // Minimum number of iterations for which the linear solver should |
| // run, even if the convergence criterion is satisfied. |
| int min_linear_solver_iterations = 0; |
| |
| // Maximum number of iterations for which the linear solver should |
| // run. If the solver does not converge in less than |
| // max_linear_solver_iterations, then it returns MAX_ITERATIONS, |
| // as its termination type. |
| int max_linear_solver_iterations = 500; |
| |
| // Maximum number of iterations performed by SCHUR_POWER_SERIES_EXPANSION. |
| // This value controls the maximum number of iterations whether it is used |
| // as a preconditioner or just to initialize the solution for |
| // ITERATIVE_SCHUR. |
| int max_num_spse_iterations = 5; |
| |
| // Use SCHUR_POWER_SERIES_EXPANSION to initialize the solution for |
| // ITERATIVE_SCHUR. This option can be set true regardless of what |
| // preconditioner is being used. |
| bool use_spse_initialization = false; |
| |
| // When use_spse_initialization is true, this parameter along with |
| // max_num_spse_iterations controls the number of |
| // SCHUR_POWER_SERIES_EXPANSION iterations performed for initialization. It |
| // is not used to control the preconditioner. |
| double spse_tolerance = 0.1; |
| |
| // Forcing sequence parameter. The truncated Newton solver uses |
| // this number to control the relative accuracy with which the |
| // Newton step is computed. |
| // |
| // This constant is passed to ConjugateGradientsSolver which uses |
| // it to terminate the iterations when |
| // |
| // (Q_i - Q_{i-1})/Q_i < eta/i |
| double eta = 1e-1; |
| |
| // Normalize the jacobian using Jacobi scaling before calling |
| // the linear least squares solver. |
| bool jacobi_scaling = true; |
| |
| // Logging options --------------------------------------------------------- |
| |
| LoggingType logging_type = PER_MINIMIZER_ITERATION; |
| |
| // By default the Minimizer progress is logged to VLOG(1), which |
| // is sent to STDERR depending on the vlog level. If this flag is |
| // set to true, and logging_type is not SILENT, the logging output |
| // is sent to STDOUT. |
| bool minimizer_progress_to_stdout = false; |
| |
| // List of iterations at which the minimizer should dump the trust |
| // region problem. Useful for testing and benchmarking. If empty |
| // (default), no problems are dumped. |
| std::vector<int> trust_region_minimizer_iterations_to_dump; |
| |
| // Directory to which the problems should be written to. Should be |
| // non-empty if trust_region_minimizer_iterations_to_dump is |
| // non-empty and trust_region_problem_dump_format_type is not |
| // CONSOLE. |
| std::string trust_region_problem_dump_directory = "/tmp"; |
| DumpFormatType trust_region_problem_dump_format_type = TEXTFILE; |
| |
| // Finite differences options ---------------------------------------------- |
| |
| // Check all jacobians computed by each residual block with finite |
| // differences. This is expensive since it involves computing the |
| // derivative by normal means (e.g. user specified, autodiff, |
| // etc), then also computing it using finite differences. The |
| // results are compared, and if they differ substantially, details |
| // are printed to the log. |
| bool check_gradients = false; |
| |
| // Relative precision to check for in the gradient checker. If the |
| // relative difference between an element in a jacobian exceeds |
| // this number, then the jacobian for that cost term is dumped. |
| double gradient_check_relative_precision = 1e-8; |
| |
| // WARNING: This option only applies to the to the numeric |
| // differentiation used for checking the user provided derivatives |
| // when when Solver::Options::check_gradients is true. If you are |
| // using NumericDiffCostFunction and are interested in changing |
| // the step size for numeric differentiation in your cost |
| // function, please have a look at |
| // include/ceres/numeric_diff_options.h. |
| // |
| // Relative shift used for taking numeric derivatives when |
| // Solver::Options::check_gradients is true. |
| // |
| // For finite differencing, each dimension is evaluated at |
| // slightly shifted values; for the case of central difference, |
| // this is what gets evaluated: |
| // |
| // delta = gradient_check_numeric_derivative_relative_step_size; |
| // f_initial = f(x) |
| // f_forward = f((1 + delta) * x) |
| // f_backward = f((1 - delta) * x) |
| // |
| // The finite differencing is done along each dimension. The |
| // reason to use a relative (rather than absolute) step size is |
| // that this way, numeric differentiation works for functions where |
| // the arguments are typically large (e.g. 1e9) and when the |
| // values are small (e.g. 1e-5). It is possible to construct |
| // "torture cases" which break this finite difference heuristic, |
| // but they do not come up often in practice. |
| // |
| // TODO(keir): Pick a smarter number than the default above! In |
| // theory a good choice is sqrt(eps) * x, which for doubles means |
| // about 1e-8 * x. However, I have found this number too |
| // optimistic. This number should be exposed for users to change. |
| double gradient_check_numeric_derivative_relative_step_size = 1e-6; |
| |
| // If update_state_every_iteration is true, then Ceres Solver will |
| // guarantee that at the end of every iteration and before any |
| // user provided IterationCallback is called, the parameter blocks |
| // are updated to the current best solution found by the |
| // solver. Thus the IterationCallback can inspect the values of |
| // the parameter blocks for purposes of computation, visualization |
| // or termination. |
| |
| // If update_state_every_iteration is false then there is no such |
| // guarantee, and user provided IterationCallbacks should not |
| // expect to look at the parameter blocks and interpret their |
| // values. |
| bool update_state_every_iteration = false; |
| |
| // Callbacks that are executed at the end of each iteration of the |
| // Minimizer. An iteration may terminate midway, either due to |
| // numerical failures or because one of the convergence tests has |
| // been satisfied. In this case none of the callbacks are |
| // executed. |
| |
| // Callbacks are executed in the order that they are specified in |
| // this vector. By default, parameter blocks are updated only at the |
| // end of the optimization, i.e when the Minimizer terminates. This |
| // behaviour is controlled by update_state_every_iteration. If the |
| // user wishes to have access to the updated parameter blocks when |
| // his/her callbacks are executed, then set |
| // update_state_every_iteration to true. |
| // |
| // The solver does NOT take ownership of these pointers. |
| std::vector<IterationCallback*> callbacks; |
| }; |
| |
| struct CERES_EXPORT Summary { |
| // A brief one line description of the state of the solver after |
| // termination. |
| std::string BriefReport() const; |
| |
| // A full multiline description of the state of the solver after |
| // termination. |
| std::string FullReport() const; |
| |
| bool IsSolutionUsable() const; |
| |
| // Minimizer summary ------------------------------------------------- |
| MinimizerType minimizer_type = TRUST_REGION; |
| |
| TerminationType termination_type = FAILURE; |
| |
| // Reason why the solver terminated. |
| std::string message = "ceres::Solve was not called."; |
| |
| // Cost of the problem (value of the objective function) before |
| // the optimization. |
| double initial_cost = -1.0; |
| |
| // Cost of the problem (value of the objective function) after the |
| // optimization. |
| double final_cost = -1.0; |
| |
| // The part of the total cost that comes from residual blocks that |
| // were held fixed by the preprocessor because all the parameter |
| // blocks that they depend on were fixed. |
| double fixed_cost = -1.0; |
| |
| // IterationSummary for each minimizer iteration in order. |
| std::vector<IterationSummary> iterations; |
| |
| // Number of minimizer iterations in which the step was |
| // accepted. Unless use_non_monotonic_steps is true this is also |
| // the number of steps in which the objective function value/cost |
| // went down. |
| int num_successful_steps = -1; |
| |
| // Number of minimizer iterations in which the step was rejected |
| // either because it did not reduce the cost enough or the step |
| // was not numerically valid. |
| int num_unsuccessful_steps = -1; |
| |
| // Number of times inner iterations were performed. |
| int num_inner_iteration_steps = -1; |
| |
| // Total number of iterations inside the line search algorithm |
| // across all invocations. We call these iterations "steps" to |
| // distinguish them from the outer iterations of the line search |
| // and trust region minimizer algorithms which call the line |
| // search algorithm as a subroutine. |
| int num_line_search_steps = -1; |
| |
| // All times reported below are wall times. |
| |
| // When the user calls Solve, before the actual optimization |
| // occurs, Ceres performs a number of preprocessing steps. These |
| // include error checks, memory allocations, and reorderings. This |
| // time is accounted for as preprocessing time. |
| double preprocessor_time_in_seconds = -1.0; |
| |
| // Time spent in the TrustRegionMinimizer. |
| double minimizer_time_in_seconds = -1.0; |
| |
| // After the Minimizer is finished, some time is spent in |
| // re-evaluating residuals etc. This time is accounted for in the |
| // postprocessor time. |
| double postprocessor_time_in_seconds = -1.0; |
| |
| // Some total of all time spent inside Ceres when Solve is called. |
| double total_time_in_seconds = -1.0; |
| |
| // Time (in seconds) spent in the linear solver computing the |
| // trust region step. |
| double linear_solver_time_in_seconds = -1.0; |
| |
| // Number of times the Newton step was computed by solving a |
| // linear system. This does not include linear solves used by |
| // inner iterations. |
| int num_linear_solves = -1; |
| |
| // Time (in seconds) spent evaluating the residual vector. |
| double residual_evaluation_time_in_seconds = -1.0; |
| |
| // Number of residual only evaluations. |
| int num_residual_evaluations = -1; |
| |
| // Time (in seconds) spent evaluating the jacobian matrix. |
| double jacobian_evaluation_time_in_seconds = -1.0; |
| |
| // Number of Jacobian (and residual) evaluations. |
| int num_jacobian_evaluations = -1; |
| |
| // Time (in seconds) spent doing inner iterations. |
| double inner_iteration_time_in_seconds = -1.0; |
| |
| // Cumulative timing information for line searches performed as part of the |
| // solve. Note that in addition to the case when the Line Search minimizer |
| // is used, the Trust Region minimizer also uses a line search when |
| // solving a constrained problem. |
| |
| // Time (in seconds) spent evaluating the univariate cost function as part |
| // of a line search. |
| double line_search_cost_evaluation_time_in_seconds = -1.0; |
| |
| // Time (in seconds) spent evaluating the gradient of the univariate cost |
| // function as part of a line search. |
| double line_search_gradient_evaluation_time_in_seconds = -1.0; |
| |
| // Time (in seconds) spent minimizing the interpolating polynomial |
| // to compute the next candidate step size as part of a line search. |
| double line_search_polynomial_minimization_time_in_seconds = -1.0; |
| |
| // Total time (in seconds) spent performing line searches. |
| double line_search_total_time_in_seconds = -1.0; |
| |
| // Number of parameter blocks in the problem. |
| int num_parameter_blocks = -1; |
| |
| // Number of parameters in the problem. |
| int num_parameters = -1; |
| |
| // Dimension of the tangent space of the problem (or the number of |
| // columns in the Jacobian for the problem). This is different |
| // from num_parameters if a parameter block is associated with a |
| // Manifold. |
| int num_effective_parameters = -1; |
| |
| // Number of residual blocks in the problem. |
| int num_residual_blocks = -1; |
| |
| // Number of residuals in the problem. |
| int num_residuals = -1; |
| |
| // Number of parameter blocks in the problem after the inactive |
| // and constant parameter blocks have been removed. A parameter |
| // block is inactive if no residual block refers to it. |
| int num_parameter_blocks_reduced = -1; |
| |
| // Number of parameters in the reduced problem. |
| int num_parameters_reduced = -1; |
| |
| // Dimension of the tangent space of the reduced problem (or the |
| // number of columns in the Jacobian for the reduced |
| // problem). This is different from num_parameters_reduced if a |
| // parameter block in the reduced problem is associated with a |
| // Manifold. |
| int num_effective_parameters_reduced = -1; |
| |
| // Number of residual blocks in the reduced problem. |
| int num_residual_blocks_reduced = -1; |
| |
| // Number of residuals in the reduced problem. |
| int num_residuals_reduced = -1; |
| |
| // Is the reduced problem bounds constrained. |
| bool is_constrained = false; |
| |
| // Number of threads specified by the user for Jacobian and |
| // residual evaluation. |
| int num_threads_given = -1; |
| |
| // Number of threads actually used by the solver for Jacobian and |
| // residual evaluation. This number is not equal to |
| // num_threads_given if OpenMP is not available. |
| int num_threads_used = -1; |
| |
| // Type of the linear solver requested by the user. |
| LinearSolverType linear_solver_type_given = |
| #if defined(CERES_NO_SPARSE) |
| DENSE_QR; |
| #else |
| SPARSE_NORMAL_CHOLESKY; |
| #endif |
| // Type of the linear solver actually used. This may be different |
| // from linear_solver_type_given if Ceres determines that the |
| // problem structure is not compatible with the linear solver |
| // requested or if the linear solver requested by the user is not |
| // available, e.g. The user requested SPARSE_NORMAL_CHOLESKY but |
| // no sparse linear algebra library was available. |
| LinearSolverType linear_solver_type_used = |
| #if defined(CERES_NO_SPARSE) |
| DENSE_QR; |
| #else |
| SPARSE_NORMAL_CHOLESKY; |
| #endif |
| |
| bool mixed_precision_solves_used = false; |
| |
| LinearSolverOrderingType linear_solver_ordering_type; |
| |
| // Size of the elimination groups given by the user as hints to |
| // the linear solver. |
| std::vector<int> linear_solver_ordering_given; |
| |
| // Size of the parameter groups used by the solver when ordering |
| // the columns of the Jacobian. This maybe different from |
| // linear_solver_ordering_given if the user left |
| // linear_solver_ordering_given blank and asked for an automatic |
| // ordering, or if the problem contains some constant or inactive |
| // parameter blocks. |
| std::vector<int> linear_solver_ordering_used; |
| |
| // For Schur type linear solvers, this string describes the |
| // template specialization which was detected in the problem and |
| // should be used. |
| std::string schur_structure_given; |
| |
| // This is the Schur template specialization that was actually |
| // instantiated and used. The reason this will be different from |
| // schur_structure_given is because the corresponding template |
| // specialization does not exist. |
| // |
| // Template specializations can be added to ceres by editing |
| // internal/ceres/generate_template_specializations.py |
| std::string schur_structure_used; |
| |
| // True if the user asked for inner iterations to be used as part |
| // of the optimization. |
| bool inner_iterations_given = false; |
| |
| // True if the user asked for inner iterations to be used as part |
| // of the optimization and the problem structure was such that |
| // they were actually performed. e.g., in a problem with just one |
| // parameter block, inner iterations are not performed. |
| bool inner_iterations_used = false; |
| |
| // Size of the parameter groups given by the user for performing |
| // inner iterations. |
| std::vector<int> inner_iteration_ordering_given; |
| |
| // Size of the parameter groups given used by the solver for |
| // performing inner iterations. This maybe different from |
| // inner_iteration_ordering_given if the user left |
| // inner_iteration_ordering_given blank and asked for an automatic |
| // ordering, or if the problem contains some constant or inactive |
| // parameter blocks. |
| std::vector<int> inner_iteration_ordering_used; |
| |
| // Type of the preconditioner requested by the user. |
| PreconditionerType preconditioner_type_given = IDENTITY; |
| |
| // Type of the preconditioner actually used. This may be different |
| // from linear_solver_type_given if Ceres determines that the |
| // problem structure is not compatible with the linear solver |
| // requested or if the linear solver requested by the user is not |
| // available. |
| PreconditionerType preconditioner_type_used = IDENTITY; |
| |
| // Type of clustering algorithm used for visibility based |
| // preconditioning. Only meaningful when the preconditioner_type |
| // is CLUSTER_JACOBI or CLUSTER_TRIDIAGONAL. |
| VisibilityClusteringType visibility_clustering_type = CANONICAL_VIEWS; |
| |
| // Type of trust region strategy. |
| TrustRegionStrategyType trust_region_strategy_type = LEVENBERG_MARQUARDT; |
| |
| // Type of dogleg strategy used for solving the trust region |
| // problem. |
| DoglegType dogleg_type = TRADITIONAL_DOGLEG; |
| |
| // Type of the dense linear algebra library used. |
| DenseLinearAlgebraLibraryType dense_linear_algebra_library_type = EIGEN; |
| |
| // Type of the sparse linear algebra library used. |
| SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type = |
| NO_SPARSE; |
| |
| // Type of line search direction used. |
| LineSearchDirectionType line_search_direction_type = LBFGS; |
| |
| // Type of the line search algorithm used. |
| LineSearchType line_search_type = WOLFE; |
| |
| // When performing line search, the degree of the polynomial used |
| // to approximate the objective function. |
| LineSearchInterpolationType line_search_interpolation_type = CUBIC; |
| |
| // If the line search direction is NONLINEAR_CONJUGATE_GRADIENT, |
| // then this indicates the particular variant of non-linear |
| // conjugate gradient used. |
| NonlinearConjugateGradientType nonlinear_conjugate_gradient_type = |
| FLETCHER_REEVES; |
| |
| // If the type of the line search direction is LBFGS, then this |
| // indicates the rank of the Hessian approximation. |
| int max_lbfgs_rank = -1; |
| }; |
| |
| // Once a least squares problem has been built, this function takes |
| // the problem and optimizes it based on the values of the options |
| // parameters. Upon return, a detailed summary of the work performed |
| // by the preprocessor, the non-linear minimizer and the linear |
| // solver are reported in the summary object. |
| virtual void Solve(const Options& options, |
| Problem* problem, |
| Solver::Summary* summary); |
| }; |
| |
| // Helper function which avoids going through the interface. |
| CERES_EXPORT void Solve(const Solver::Options& options, |
| Problem* problem, |
| Solver::Summary* summary); |
| |
| } // namespace ceres |
| |
| #include "ceres/internal/reenable_warnings.h" |
| |
| #endif // CERES_PUBLIC_SOLVER_H_ |