Initial commit of Ceres Solver.
diff --git a/include/ceres/solver.h b/include/ceres/solver.h new file mode 100644 index 0000000..15fd733 --- /dev/null +++ b/include/ceres/solver.h
@@ -0,0 +1,379 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2010, 2011, 2012 Google Inc. All rights reserved. +// http://code.google.com/p/ceres-solver/ +// +// 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 <string> +#include <vector> + +#include "ceres/iteration_callback.h" +#include "ceres/internal/macros.h" +#include "ceres/internal/port.h" +#include "ceres/types.h" + +namespace ceres { + +class Problem; + +// Interface for non-linear least squares solvers. +class 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 Options { + // Default constructor that sets up a generic sparse problem. + Options() { + minimizer_type = LEVENBERG_MARQUARDT; + max_num_iterations = 50; + max_solver_time_sec = 1.0e9; + num_threads = 1; + tau = 1e-4; + min_relative_decrease = 1e-3; + function_tolerance = 1e-6; + gradient_tolerance = 1e-10; + parameter_tolerance = 1e-8; +#ifndef CERES_NO_SUITESPARSE + linear_solver_type = SPARSE_NORMAL_CHOLESKY; +#else + linear_solver_type = DENSE_QR; +#endif // CERES_NO_SUITESPARSE + preconditioner_type = JACOBI; + num_linear_solver_threads = 1; + num_eliminate_blocks = 0; + ordering_type = NATURAL; + linear_solver_min_num_iterations = 1; + linear_solver_max_num_iterations = 500; + eta = 1e-1; + jacobi_scaling = true; + logging_type = PER_MINIMIZER_ITERATION; + minimizer_progress_to_stdout = false; + return_initial_residuals = false; + return_final_residuals = false; + lsqp_dump_format = "lm_iteration_%03d.lsqp"; + crash_and_dump_lsqp_on_failure = false; + check_gradients = false; + gradient_check_relative_precision = 1e-8; + numeric_derivative_relative_step_size = 1e-6; + update_state_every_iteration = false; + } + + // Minimizer options ---------------------------------------- + + MinimizerType minimizer_type; + + // Maximum number of iterations for the minimizer to run for. + int max_num_iterations; + + // Maximum time for which the minimizer should run for. + double max_solver_time_sec; + + // Number of threads used by Ceres for evaluating the cost and + // jacobians. + int num_threads; + + // For Levenberg-Marquardt, the initial value for the + // regularizer. This is the inversely related to the size of the + // initial trust region. + double tau; + + // For trust region methods, this is lower threshold for the + // relative decrease before a step is accepted. + double min_relative_decrease; + + // Minimizer terminates when + // + // (new_cost - old_cost) < function_tolerance * old_cost; + // + double function_tolerance; + + // Minimizer terminates when + // + // max_i |gradient_i| < gradient_tolerance * max_i|initial_gradient_i| + // + // This value should typically be 1e-4 * function_tolerance. + double gradient_tolerance; + + // Minimizer terminates when + // + // |step|_2 <= parameter_tolerance * ( |x|_2 + parameter_tolerance) + // + double parameter_tolerance; + + // Linear least squares solver options ------------------------------------- + + LinearSolverType linear_solver_type; + + // Type of preconditioner to use with the iterative linear solvers. + PreconditionerType preconditioner_type; + + // Number of threads used by Ceres to solve the Newton + // step. Currently only the SPARSE_SCHUR solver is capable of + // using this setting. + int num_linear_solver_threads; + + // For Schur reduction based methods, the first 0 to num blocks are + // eliminated using the Schur reduction. For example, when solving + // traditional structure from motion problems where the parameters are in + // two classes (cameras and points) then num_eliminate_blocks would be the + // number of points. + // + // This parameter is used in conjunction with the ordering. + // Applies to: Preprocessor and linear least squares solver. + int num_eliminate_blocks; + + // Internally Ceres reorders the parameter blocks to help the + // various linear solvers. This parameter allows the user to + // influence the re-ordering strategy used. For structure from + // motion problems use SCHUR, for other problems NATURAL (default) + // is a good choice. In case you wish to specify your own ordering + // scheme, for example in conjunction with num_eliminate_blocks, + // use USER. + OrderingType ordering_type; + + // The ordering of the parameter blocks. The solver pays attention + // to it if the ordering_type is set to USER and the vector is + // non-empty. + vector<double*> ordering; + + + // Minimum number of iterations for which the linear solver should + // run, even if the convergence criterion is satisfied. + int linear_solver_min_num_iterations; + + // Maximum number of iterations for which the linear solver should + // run. If the solver does not converge in less than + // linear_solver_max_num_iterations, then it returns + // MAX_ITERATIONS, as its termination type. + int linear_solver_max_num_iterations; + + // 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; + + // Normalize the jacobian using Jacobi scaling before calling + // the linear least squares solver. + bool jacobi_scaling; + + // Logging options --------------------------------------------------------- + + LoggingType logging_type; + + // 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; + + bool return_initial_residuals; + bool return_final_residuals; + + // List of iterations at which the optimizer should dump the + // linear least squares problem to disk. Useful for testing and + // benchmarking. If empty (default), no problems are dumped. + // + // This is ignored if protocol buffers are disabled. + vector<int> lsqp_iterations_to_dump; + + // Format string for the file name used for dumping the least + // squares problem to disk. If the format is 'ascii', then the + // problem is logged to the screen; don't try this with large + // problems or expect a frozen terminal. + string lsqp_dump_format; + + // Dump the linear least squares problem to disk if the minimizer + // fails due to NUMERICAL_FAILURE and crash the process. This flag + // is useful for generating debugging information. The problem is + // dumped in a file whose name is determined by + // Solver::Options::lsqp_dump_format. + // + // Note: This requires a version of Ceres built with protocol buffers. + bool crash_and_dump_lsqp_on_failure; + + // 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; + + // 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; + + // Relative shift used for taking numeric derivatives. For finite + // differencing, each dimension is evaluated at slightly shifted + // values; for the case of central difference, this is what gets + // evaluated: + // + // delta = 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 differentation 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 numeric_derivative_relative_step_size; + + // If true, the user's parameter blocks are updated at the end of + // every Minimizer iteration, otherwise they are updated when the + // Minimizer terminates. This is useful if, for example, the user + // wishes to visualize the state of the optimization every + // iteration. + bool update_state_every_iteration; + + // Callbacks that are executed at the end of each iteration of the + // Minimizer. They 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_variable. If the user wishes to have access + // to the update 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. + vector<IterationCallback*> callbacks; + }; + + struct Summary { + Summary(); + + // A brief one line description of the state of the solver after + // termination. + string BriefReport() const; + + // A full multiline description of the state of the solver after + // termination. + string FullReport() const; + + // Minimizer summary ------------------------------------------------- + SolverTerminationType termination_type; + + // If the solver did not run, or there was a failure, a + // description of the error. + string error; + + // Cost of the problem before and after the optimization. See + // problem.h for definition of the cost of a problem. + double initial_cost; + double final_cost; + + // 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; + + // Residuals before and after the optimization. Each vector + // contains problem.NumResiduals() elements. Residuals are in the + // same order in which they were added to the problem object when + // constructing this problem. + vector<double> initial_residuals; + vector<double> final_residuals; + + vector<IterationSummary> iterations; + + int num_successful_steps; + int num_unsuccessful_steps; + + double preprocessor_time_in_seconds; + double minimizer_time_in_seconds; + double total_time_in_seconds; + + // Preprocessor summary. + int num_parameter_blocks; + int num_parameters; + int num_residual_blocks; + int num_residuals; + + int num_parameter_blocks_reduced; + int num_parameters_reduced; + int num_residual_blocks_reduced; + int num_residuals_reduced; + + int num_eliminate_blocks_given; + int num_eliminate_blocks_used; + + int num_threads_given; + int num_threads_used; + + int num_linear_solver_threads_given; + int num_linear_solver_threads_used; + + LinearSolverType linear_solver_type_given; + LinearSolverType linear_solver_type_used; + + PreconditionerType preconditioner_type; + OrderingType ordering_type; + }; + + // 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 minmizer 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. +void Solve(const Solver::Options& options, + Problem* problem, + Solver::Summary* summary); + +} // namespace ceres + +#endif // CERES_PUBLIC_SOLVER_H_