Squashed 'third_party/ceres/' content from commit e51e9b4
Change-Id: I763587619d57e594d3fa158dc3a7fe0b89a1743b
git-subtree-dir: third_party/ceres
git-subtree-split: e51e9b46f6ca88ab8b2266d0e362771db6d98067
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h
new file mode 100644
index 0000000..24c245d
--- /dev/null
+++ b/internal/ceres/linear_solver.h
@@ -0,0 +1,339 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2015 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)
+//
+// Abstract interface for objects solving linear systems of various
+// kinds.
+
+#ifndef CERES_INTERNAL_LINEAR_SOLVER_H_
+#define CERES_INTERNAL_LINEAR_SOLVER_H_
+
+#include <cstddef>
+#include <map>
+#include <string>
+#include <vector>
+#include "ceres/block_sparse_matrix.h"
+#include "ceres/casts.h"
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/context_impl.h"
+#include "ceres/dense_sparse_matrix.h"
+#include "ceres/execution_summary.h"
+#include "ceres/triplet_sparse_matrix.h"
+#include "ceres/types.h"
+#include "glog/logging.h"
+
+namespace ceres {
+namespace internal {
+
+enum LinearSolverTerminationType {
+ // Termination criterion was met.
+ LINEAR_SOLVER_SUCCESS,
+
+ // Solver ran for max_num_iterations and terminated before the
+ // termination tolerance could be satisfied.
+ LINEAR_SOLVER_NO_CONVERGENCE,
+
+ // Solver was terminated due to numerical problems, generally due to
+ // the linear system being poorly conditioned.
+ LINEAR_SOLVER_FAILURE,
+
+ // Solver failed with a fatal error that cannot be recovered from,
+ // e.g. CHOLMOD ran out of memory when computing the symbolic or
+ // numeric factorization or an underlying library was called with
+ // the wrong arguments.
+ LINEAR_SOLVER_FATAL_ERROR
+};
+
+// This enum controls the fill-reducing ordering a sparse linear
+// algebra library should use before computing a sparse factorization
+// (usually Cholesky).
+enum OrderingType {
+ NATURAL, // Do not re-order the matrix. This is useful when the
+ // matrix has been ordered using a fill-reducing ordering
+ // already.
+ AMD // Use the Approximate Minimum Degree algorithm to re-order
+ // the matrix.
+};
+
+class LinearOperator;
+
+// Abstract base class for objects that implement algorithms for
+// solving linear systems
+//
+// Ax = b
+//
+// It is expected that a single instance of a LinearSolver object
+// maybe used multiple times for solving multiple linear systems with
+// the same sparsity structure. This allows them to cache and reuse
+// information across solves. This means that calling Solve on the
+// same LinearSolver instance with two different linear systems will
+// result in undefined behaviour.
+//
+// Subclasses of LinearSolver use two structs to configure themselves.
+// The Options struct configures the LinearSolver object for its
+// lifetime. The PerSolveOptions struct is used to specify options for
+// a particular Solve call.
+class LinearSolver {
+ public:
+ struct Options {
+ LinearSolverType type = SPARSE_NORMAL_CHOLESKY;
+ PreconditionerType preconditioner_type = JACOBI;
+ VisibilityClusteringType visibility_clustering_type = CANONICAL_VIEWS;
+ DenseLinearAlgebraLibraryType dense_linear_algebra_library_type = EIGEN;
+ SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type =
+ SUITE_SPARSE;
+
+ // See solver.h for information about these flags.
+ bool use_postordering = false;
+ bool dynamic_sparsity = false;
+ bool use_explicit_schur_complement = false;
+
+ // Number of internal iterations that the solver uses. This
+ // parameter only makes sense for iterative solvers like CG.
+ int min_num_iterations = 1;
+ int max_num_iterations = 1;
+
+ // If possible, how many threads can the solver use.
+ int num_threads = 1;
+
+ // Hints about the order in which the parameter blocks should be
+ // eliminated by the linear solver.
+ //
+ // For example if elimination_groups is a vector of size k, then
+ // the linear solver is informed that it should eliminate the
+ // parameter blocks 0 ... elimination_groups[0] - 1 first, and
+ // then elimination_groups[0] ... elimination_groups[1] - 1 and so
+ // on. Within each elimination group, the linear solver is free to
+ // choose how the parameter blocks are ordered. Different linear
+ // solvers have differing requirements on elimination_groups.
+ //
+ // The most common use is for Schur type solvers, where there
+ // should be at least two elimination groups and the first
+ // elimination group must form an independent set in the normal
+ // equations. The first elimination group corresponds to the
+ // num_eliminate_blocks in the Schur type solvers.
+ std::vector<int> elimination_groups;
+
+ // Iterative solvers, e.g. Preconditioned Conjugate Gradients
+ // maintain a cheap estimate of the residual which may become
+ // inaccurate over time. Thus for non-zero values of this
+ // parameter, the solver can be told to recalculate the value of
+ // the residual using a |b - Ax| evaluation.
+ int residual_reset_period = 10;
+
+ // If the block sizes in a BlockSparseMatrix are fixed, then in
+ // some cases the Schur complement based solvers can detect and
+ // specialize on them.
+ //
+ // It is expected that these parameters are set programmatically
+ // rather than manually.
+ //
+ // Please see schur_complement_solver.h and schur_eliminator.h for
+ // more details.
+ int row_block_size = Eigen::Dynamic;
+ int e_block_size = Eigen::Dynamic;
+ int f_block_size = Eigen::Dynamic;
+
+ bool use_mixed_precision_solves = false;
+ int max_num_refinement_iterations = 0;
+ ContextImpl* context = nullptr;
+ };
+
+ // Options for the Solve method.
+ struct PerSolveOptions {
+ // This option only makes sense for unsymmetric linear solvers
+ // that can solve rectangular linear systems.
+ //
+ // Given a matrix A, an optional diagonal matrix D as a vector,
+ // and a vector b, the linear solver will solve for
+ //
+ // | A | x = | b |
+ // | D | | 0 |
+ //
+ // If D is null, then it is treated as zero, and the solver returns
+ // the solution to
+ //
+ // A x = b
+ //
+ // In either case, x is the vector that solves the following
+ // optimization problem.
+ //
+ // arg min_x ||Ax - b||^2 + ||Dx||^2
+ //
+ // Here A is a matrix of size m x n, with full column rank. If A
+ // does not have full column rank, the results returned by the
+ // solver cannot be relied on. D, if it is not null is an array of
+ // size n. b is an array of size m and x is an array of size n.
+ double* D = nullptr;
+
+ // This option only makes sense for iterative solvers.
+ //
+ // In general the performance of an iterative linear solver
+ // depends on the condition number of the matrix A. For example
+ // the convergence rate of the conjugate gradients algorithm
+ // is proportional to the square root of the condition number.
+ //
+ // One particularly useful technique for improving the
+ // conditioning of a linear system is to precondition it. In its
+ // simplest form a preconditioner is a matrix M such that instead
+ // of solving Ax = b, we solve the linear system AM^{-1} y = b
+ // instead, where M is such that the condition number k(AM^{-1})
+ // is smaller than the conditioner k(A). Given the solution to
+ // this system, x = M^{-1} y. The iterative solver takes care of
+ // the mechanics of solving the preconditioned system and
+ // returning the corrected solution x. The user only needs to
+ // supply a linear operator.
+ //
+ // A null preconditioner is equivalent to an identity matrix being
+ // used a preconditioner.
+ LinearOperator* preconditioner = nullptr;
+
+
+ // The following tolerance related options only makes sense for
+ // iterative solvers. Direct solvers ignore them.
+
+ // Solver terminates when
+ //
+ // |Ax - b| <= r_tolerance * |b|.
+ //
+ // This is the most commonly used termination criterion for
+ // iterative solvers.
+ double r_tolerance = 0.0;
+
+ // For PSD matrices A, let
+ //
+ // Q(x) = x'Ax - 2b'x
+ //
+ // be the cost of the quadratic function defined by A and b. Then,
+ // the solver terminates at iteration i if
+ //
+ // i * (Q(x_i) - Q(x_i-1)) / Q(x_i) < q_tolerance.
+ //
+ // This termination criterion is more useful when using CG to
+ // solve the Newton step. This particular convergence test comes
+ // from Stephen Nash's work on truncated Newton
+ // methods. References:
+ //
+ // 1. Stephen G. Nash & Ariela Sofer, Assessing A Search
+ // Direction Within A Truncated Newton Method, Operation
+ // Research Letters 9(1990) 219-221.
+ //
+ // 2. Stephen G. Nash, A Survey of Truncated Newton Methods,
+ // Journal of Computational and Applied Mathematics,
+ // 124(1-2), 45-59, 2000.
+ //
+ double q_tolerance = 0.0;
+ };
+
+ // Summary of a call to the Solve method. We should move away from
+ // the true/false method for determining solver success. We should
+ // let the summary object do the talking.
+ struct Summary {
+ double residual_norm = -1.0;
+ int num_iterations = -1;
+ LinearSolverTerminationType termination_type = LINEAR_SOLVER_FAILURE;
+ std::string message;
+ };
+
+ // If the optimization problem is such that there are no remaining
+ // e-blocks, a Schur type linear solver cannot be used. If the
+ // linear solver is of Schur type, this function implements a policy
+ // to select an alternate nearest linear solver to the one selected
+ // by the user. The input linear_solver_type is returned otherwise.
+ static LinearSolverType LinearSolverForZeroEBlocks(
+ LinearSolverType linear_solver_type);
+
+ virtual ~LinearSolver();
+
+ // Solve Ax = b.
+ virtual Summary Solve(LinearOperator* A,
+ const double* b,
+ const PerSolveOptions& per_solve_options,
+ double* x) = 0;
+
+ // This method returns copies instead of references so that the base
+ // class implementation does not have to worry about life time
+ // issues. Further, this calls are not expected to be frequent or
+ // performance sensitive.
+ virtual std::map<std::string, CallStatistics> Statistics() const {
+ return std::map<std::string, CallStatistics>();
+ }
+
+ // Factory
+ static LinearSolver* Create(const Options& options);
+};
+
+// This templated subclass of LinearSolver serves as a base class for
+// other linear solvers that depend on the particular matrix layout of
+// the underlying linear operator. For example some linear solvers
+// need low level access to the TripletSparseMatrix implementing the
+// LinearOperator interface. This class hides those implementation
+// details behind a private virtual method, and has the Solve method
+// perform the necessary upcasting.
+template <typename MatrixType>
+class TypedLinearSolver : public LinearSolver {
+ public:
+ virtual ~TypedLinearSolver() {}
+ virtual LinearSolver::Summary Solve(
+ LinearOperator* A,
+ const double* b,
+ const LinearSolver::PerSolveOptions& per_solve_options,
+ double* x) {
+ ScopedExecutionTimer total_time("LinearSolver::Solve", &execution_summary_);
+ CHECK(A != nullptr);
+ CHECK(b != nullptr);
+ CHECK(x != nullptr);
+ return SolveImpl(down_cast<MatrixType*>(A), b, per_solve_options, x);
+ }
+
+ virtual std::map<std::string, CallStatistics> Statistics() const {
+ return execution_summary_.statistics();
+ }
+
+ private:
+ virtual LinearSolver::Summary SolveImpl(
+ MatrixType* A,
+ const double* b,
+ const LinearSolver::PerSolveOptions& per_solve_options,
+ double* x) = 0;
+
+ ExecutionSummary execution_summary_;
+};
+
+// Linear solvers that depend on acccess to the low level structure of
+// a SparseMatrix.
+typedef TypedLinearSolver<BlockSparseMatrix> BlockSparseMatrixSolver; // NOLINT
+typedef TypedLinearSolver<CompressedRowSparseMatrix> CompressedRowSparseMatrixSolver; // NOLINT
+typedef TypedLinearSolver<DenseSparseMatrix> DenseSparseMatrixSolver; // NOLINT
+typedef TypedLinearSolver<TripletSparseMatrix> TripletSparseMatrixSolver; // NOLINT
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_LINEAR_SOLVER_H_