1 // Ceres Solver - A fast non-linear least squares minimizer
2 // Copyright 2015 Google Inc. All rights reserved.
3 // http://ceres-solver.org/
5 // Redistribution and use in source and binary forms, with or without
6 // modification, are permitted provided that the following conditions are met:
8 // * Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright notice,
11 // this list of conditions and the following disclaimer in the documentation
12 // and/or other materials provided with the distribution.
13 // * Neither the name of Google Inc. nor the names of its contributors may be
14 // used to endorse or promote products derived from this software without
15 // specific prior written permission.
17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 // POSSIBILITY OF SUCH DAMAGE.
29 // Author: sameeragarwal@google.com (Sameer Agarwal)
31 // Abstract interface for objects solving linear systems of various
34 #ifndef CERES_INTERNAL_LINEAR_SOLVER_H_
35 #define CERES_INTERNAL_LINEAR_SOLVER_H_
41 #include "ceres/block_sparse_matrix.h"
42 #include "ceres/casts.h"
43 #include "ceres/compressed_row_sparse_matrix.h"
44 #include "ceres/dense_sparse_matrix.h"
45 #include "ceres/execution_summary.h"
46 #include "ceres/triplet_sparse_matrix.h"
47 #include "ceres/types.h"
48 #include "glog/logging.h"
53 enum LinearSolverTerminationType {
54 // Termination criterion was met.
55 LINEAR_SOLVER_SUCCESS,
57 // Solver ran for max_num_iterations and terminated before the
58 // termination tolerance could be satisfied.
59 LINEAR_SOLVER_NO_CONVERGENCE,
61 // Solver was terminated due to numerical problems, generally due to
62 // the linear system being poorly conditioned.
63 LINEAR_SOLVER_FAILURE,
65 // Solver failed with a fatal error that cannot be recovered from,
66 // e.g. CHOLMOD ran out of memory when computing the symbolic or
67 // numeric factorization or an underlying library was called with
68 // the wrong arguments.
69 LINEAR_SOLVER_FATAL_ERROR
72 // This enum controls the fill-reducing ordering a sparse linear
73 // algebra library should use before computing a sparse factorization
74 // (usually Cholesky).
76 NATURAL, // Do not re-order the matrix. This is useful when the
77 // matrix has been ordered using a fill-reducing ordering
79 AMD // Use the Approximate Minimum Degree algorithm to re-order
85 // Abstract base class for objects that implement algorithms for
86 // solving linear systems
90 // It is expected that a single instance of a LinearSolver object
91 // maybe used multiple times for solving multiple linear systems with
92 // the same sparsity structure. This allows them to cache and reuse
93 // information across solves. This means that calling Solve on the
94 // same LinearSolver instance with two different linear systems will
95 // result in undefined behaviour.
97 // Subclasses of LinearSolver use two structs to configure themselves.
98 // The Options struct configures the LinearSolver object for its
99 // lifetime. The PerSolveOptions struct is used to specify options for
100 // a particular Solve call.
105 : type(SPARSE_NORMAL_CHOLESKY),
106 preconditioner_type(JACOBI),
107 visibility_clustering_type(CANONICAL_VIEWS),
108 dense_linear_algebra_library_type(EIGEN),
109 sparse_linear_algebra_library_type(SUITE_SPARSE),
110 use_postordering(false),
111 dynamic_sparsity(false),
112 use_explicit_schur_complement(false),
113 min_num_iterations(1),
114 max_num_iterations(1),
116 residual_reset_period(10),
117 row_block_size(Eigen::Dynamic),
118 e_block_size(Eigen::Dynamic),
119 f_block_size(Eigen::Dynamic) {
122 LinearSolverType type;
123 PreconditionerType preconditioner_type;
124 VisibilityClusteringType visibility_clustering_type;
125 DenseLinearAlgebraLibraryType dense_linear_algebra_library_type;
126 SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type;
128 // See solver.h for information about these flags.
129 bool use_postordering;
130 bool dynamic_sparsity;
131 bool use_explicit_schur_complement;
133 // Number of internal iterations that the solver uses. This
134 // parameter only makes sense for iterative solvers like CG.
135 int min_num_iterations;
136 int max_num_iterations;
138 // If possible, how many threads can the solver use.
141 // Hints about the order in which the parameter blocks should be
142 // eliminated by the linear solver.
144 // For example if elimination_groups is a vector of size k, then
145 // the linear solver is informed that it should eliminate the
146 // parameter blocks 0 ... elimination_groups[0] - 1 first, and
147 // then elimination_groups[0] ... elimination_groups[1] - 1 and so
148 // on. Within each elimination group, the linear solver is free to
149 // choose how the parameter blocks are ordered. Different linear
150 // solvers have differing requirements on elimination_groups.
152 // The most common use is for Schur type solvers, where there
153 // should be at least two elimination groups and the first
154 // elimination group must form an independent set in the normal
155 // equations. The first elimination group corresponds to the
156 // num_eliminate_blocks in the Schur type solvers.
157 std::vector<int> elimination_groups;
159 // Iterative solvers, e.g. Preconditioned Conjugate Gradients
160 // maintain a cheap estimate of the residual which may become
161 // inaccurate over time. Thus for non-zero values of this
162 // parameter, the solver can be told to recalculate the value of
163 // the residual using a |b - Ax| evaluation.
164 int residual_reset_period;
166 // If the block sizes in a BlockSparseMatrix are fixed, then in
167 // some cases the Schur complement based solvers can detect and
168 // specialize on them.
170 // It is expected that these parameters are set programmatically
171 // rather than manually.
173 // Please see schur_complement_solver.h and schur_eliminator.h for
180 // Options for the Solve method.
181 struct PerSolveOptions {
184 preconditioner(NULL),
189 // This option only makes sense for unsymmetric linear solvers
190 // that can solve rectangular linear systems.
192 // Given a matrix A, an optional diagonal matrix D as a vector,
193 // and a vector b, the linear solver will solve for
198 // If D is null, then it is treated as zero, and the solver returns
203 // In either case, x is the vector that solves the following
204 // optimization problem.
206 // arg min_x ||Ax - b||^2 + ||Dx||^2
208 // Here A is a matrix of size m x n, with full column rank. If A
209 // does not have full column rank, the results returned by the
210 // solver cannot be relied on. D, if it is not null is an array of
211 // size n. b is an array of size m and x is an array of size n.
214 // This option only makes sense for iterative solvers.
216 // In general the performance of an iterative linear solver
217 // depends on the condition number of the matrix A. For example
218 // the convergence rate of the conjugate gradients algorithm
219 // is proportional to the square root of the condition number.
221 // One particularly useful technique for improving the
222 // conditioning of a linear system is to precondition it. In its
223 // simplest form a preconditioner is a matrix M such that instead
224 // of solving Ax = b, we solve the linear system AM^{-1} y = b
225 // instead, where M is such that the condition number k(AM^{-1})
226 // is smaller than the conditioner k(A). Given the solution to
227 // this system, x = M^{-1} y. The iterative solver takes care of
228 // the mechanics of solving the preconditioned system and
229 // returning the corrected solution x. The user only needs to
230 // supply a linear operator.
232 // A null preconditioner is equivalent to an identity matrix being
233 // used a preconditioner.
234 LinearOperator* preconditioner;
237 // The following tolerance related options only makes sense for
238 // iterative solvers. Direct solvers ignore them.
240 // Solver terminates when
242 // |Ax - b| <= r_tolerance * |b|.
244 // This is the most commonly used termination criterion for
245 // iterative solvers.
248 // For PSD matrices A, let
250 // Q(x) = x'Ax - 2b'x
252 // be the cost of the quadratic function defined by A and b. Then,
253 // the solver terminates at iteration i if
255 // i * (Q(x_i) - Q(x_i-1)) / Q(x_i) < q_tolerance.
257 // This termination criterion is more useful when using CG to
258 // solve the Newton step. This particular convergence test comes
259 // from Stephen Nash's work on truncated Newton
260 // methods. References:
262 // 1. Stephen G. Nash & Ariela Sofer, Assessing A Search
263 // Direction Within A Truncated Newton Method, Operation
264 // Research Letters 9(1990) 219-221.
266 // 2. Stephen G. Nash, A Survey of Truncated Newton Methods,
267 // Journal of Computational and Applied Mathematics,
268 // 124(1-2), 45-59, 2000.
273 // Summary of a call to the Solve method. We should move away from
274 // the true/false method for determining solver success. We should
275 // let the summary object do the talking.
278 : residual_norm(0.0),
280 termination_type(LINEAR_SOLVER_FAILURE) {
283 double residual_norm;
285 LinearSolverTerminationType termination_type;
289 // If the optimization problem is such that there are no remaining
290 // e-blocks, a Schur type linear solver cannot be used. If the
291 // linear solver is of Schur type, this function implements a policy
292 // to select an alternate nearest linear solver to the one selected
293 // by the user. The input linear_solver_type is returned otherwise.
294 static LinearSolverType LinearSolverForZeroEBlocks(
295 LinearSolverType linear_solver_type);
297 virtual ~LinearSolver();
300 virtual Summary Solve(LinearOperator* A,
302 const PerSolveOptions& per_solve_options,
305 // The following two methods return copies instead of references so
306 // that the base class implementation does not have to worry about
307 // life time issues. Further, these calls are not expected to be
308 // frequent or performance sensitive.
309 virtual std::map<std::string, int> CallStatistics() const {
310 return std::map<std::string, int>();
313 virtual std::map<std::string, double> TimeStatistics() const {
314 return std::map<std::string, double>();
318 static LinearSolver* Create(const Options& options);
321 // This templated subclass of LinearSolver serves as a base class for
322 // other linear solvers that depend on the particular matrix layout of
323 // the underlying linear operator. For example some linear solvers
324 // need low level access to the TripletSparseMatrix implementing the
325 // LinearOperator interface. This class hides those implementation
326 // details behind a private virtual method, and has the Solve method
327 // perform the necessary upcasting.
328 template <typename MatrixType>
329 class TypedLinearSolver : public LinearSolver {
331 virtual ~TypedLinearSolver() {}
332 virtual LinearSolver::Summary Solve(
335 const LinearSolver::PerSolveOptions& per_solve_options,
337 ScopedExecutionTimer total_time("LinearSolver::Solve", &execution_summary_);
341 return SolveImpl(down_cast<MatrixType*>(A), b, per_solve_options, x);
344 virtual std::map<std::string, int> CallStatistics() const {
345 return execution_summary_.calls();
348 virtual std::map<std::string, double> TimeStatistics() const {
349 return execution_summary_.times();
353 virtual LinearSolver::Summary SolveImpl(
356 const LinearSolver::PerSolveOptions& per_solve_options,
359 ExecutionSummary execution_summary_;
362 // Linear solvers that depend on acccess to the low level structure of
364 typedef TypedLinearSolver<BlockSparseMatrix> BlockSparseMatrixSolver; // NOLINT
365 typedef TypedLinearSolver<CompressedRowSparseMatrix> CompressedRowSparseMatrixSolver; // NOLINT
366 typedef TypedLinearSolver<DenseSparseMatrix> DenseSparseMatrixSolver; // NOLINT
367 typedef TypedLinearSolver<TripletSparseMatrix> TripletSparseMatrixSolver; // NOLINT
369 } // namespace internal
372 #endif // CERES_INTERNAL_LINEAR_SOLVER_H_