ceres-solver / ceres-solver / e269b64f55f30f5e9802632b79479aa24d9d07de / . / include / ceres / solver.h

// 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_ |