1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126 |
- // Ceres Solver - A fast non-linear least squares minimizer
- // Copyright 2023 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;
- // 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;
- // 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.
- //
- // This options is available when linear solver uses sparse or dense
- // cholesky factorization, except when sparse_linear_algebra_library_type =
- // SUITE_SPARSE.
- 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;
- // 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.
- // Each iteration corresponds to one more term in the power series expansion
- // od the inverse of the Schur complement. 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;
- // 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;
- 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_nonmonotonic_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.
- 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 = AMD;
- // 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_used
- // 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_
|