1// Ceres Solver - A fast non-linear least squares minimizer 2// Copyright 2010, 2011, 2012 Google Inc. All rights reserved. 3// http://code.google.com/p/ceres-solver/ 4// 5// Redistribution and use in source and binary forms, with or without 6// modification, are permitted provided that the following conditions are met: 7// 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. 16// 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. 28// 29// Author: sameeragarwal@google.com (Sameer Agarwal) 30// 31// Abstract interface for objects solving linear systems of various 32// kinds. 33 34#ifndef CERES_INTERNAL_LINEAR_SOLVER_H_ 35#define CERES_INTERNAL_LINEAR_SOLVER_H_ 36 37#include <cstddef> 38#include <map> 39#include <string> 40#include <vector> 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" 49 50namespace ceres { 51namespace internal { 52 53class LinearOperator; 54 55// Abstract base class for objects that implement algorithms for 56// solving linear systems 57// 58// Ax = b 59// 60// It is expected that a single instance of a LinearSolver object 61// maybe used multiple times for solving multiple linear systems with 62// the same sparsity structure. This allows them to cache and reuse 63// information across solves. This means that calling Solve on the 64// same LinearSolver instance with two different linear systems will 65// result in undefined behaviour. 66// 67// Subclasses of LinearSolver use two structs to configure themselves. 68// The Options struct configures the LinearSolver object for its 69// lifetime. The PerSolveOptions struct is used to specify options for 70// a particular Solve call. 71class LinearSolver { 72 public: 73 struct Options { 74 Options() 75 : type(SPARSE_NORMAL_CHOLESKY), 76 preconditioner_type(JACOBI), 77 dense_linear_algebra_library_type(EIGEN), 78 sparse_linear_algebra_library_type(SUITE_SPARSE), 79 use_postordering(false), 80 min_num_iterations(1), 81 max_num_iterations(1), 82 num_threads(1), 83 residual_reset_period(10), 84 row_block_size(Eigen::Dynamic), 85 e_block_size(Eigen::Dynamic), 86 f_block_size(Eigen::Dynamic) { 87 } 88 89 LinearSolverType type; 90 91 PreconditionerType preconditioner_type; 92 93 DenseLinearAlgebraLibraryType dense_linear_algebra_library_type; 94 SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type; 95 96 // See solver.h for information about this flag. 97 bool use_postordering; 98 99 // Number of internal iterations that the solver uses. This 100 // parameter only makes sense for iterative solvers like CG. 101 int min_num_iterations; 102 int max_num_iterations; 103 104 // If possible, how many threads can the solver use. 105 int num_threads; 106 107 // Hints about the order in which the parameter blocks should be 108 // eliminated by the linear solver. 109 // 110 // For example if elimination_groups is a vector of size k, then 111 // the linear solver is informed that it should eliminate the 112 // parameter blocks 0 ... elimination_groups[0] - 1 first, and 113 // then elimination_groups[0] ... elimination_groups[1] - 1 and so 114 // on. Within each elimination group, the linear solver is free to 115 // choose how the parameter blocks are ordered. Different linear 116 // solvers have differing requirements on elimination_groups. 117 // 118 // The most common use is for Schur type solvers, where there 119 // should be at least two elimination groups and the first 120 // elimination group must form an independent set in the normal 121 // equations. The first elimination group corresponds to the 122 // num_eliminate_blocks in the Schur type solvers. 123 vector<int> elimination_groups; 124 125 // Iterative solvers, e.g. Preconditioned Conjugate Gradients 126 // maintain a cheap estimate of the residual which may become 127 // inaccurate over time. Thus for non-zero values of this 128 // parameter, the solver can be told to recalculate the value of 129 // the residual using a |b - Ax| evaluation. 130 int residual_reset_period; 131 132 // If the block sizes in a BlockSparseMatrix are fixed, then in 133 // some cases the Schur complement based solvers can detect and 134 // specialize on them. 135 // 136 // It is expected that these parameters are set programmatically 137 // rather than manually. 138 // 139 // Please see schur_complement_solver.h and schur_eliminator.h for 140 // more details. 141 int row_block_size; 142 int e_block_size; 143 int f_block_size; 144 }; 145 146 // Options for the Solve method. 147 struct PerSolveOptions { 148 PerSolveOptions() 149 : D(NULL), 150 preconditioner(NULL), 151 r_tolerance(0.0), 152 q_tolerance(0.0) { 153 } 154 155 // This option only makes sense for unsymmetric linear solvers 156 // that can solve rectangular linear systems. 157 // 158 // Given a matrix A, an optional diagonal matrix D as a vector, 159 // and a vector b, the linear solver will solve for 160 // 161 // | A | x = | b | 162 // | D | | 0 | 163 // 164 // If D is null, then it is treated as zero, and the solver returns 165 // the solution to 166 // 167 // A x = b 168 // 169 // In either case, x is the vector that solves the following 170 // optimization problem. 171 // 172 // arg min_x ||Ax - b||^2 + ||Dx||^2 173 // 174 // Here A is a matrix of size m x n, with full column rank. If A 175 // does not have full column rank, the results returned by the 176 // solver cannot be relied on. D, if it is not null is an array of 177 // size n. b is an array of size m and x is an array of size n. 178 double * D; 179 180 // This option only makes sense for iterative solvers. 181 // 182 // In general the performance of an iterative linear solver 183 // depends on the condition number of the matrix A. For example 184 // the convergence rate of the conjugate gradients algorithm 185 // is proportional to the square root of the condition number. 186 // 187 // One particularly useful technique for improving the 188 // conditioning of a linear system is to precondition it. In its 189 // simplest form a preconditioner is a matrix M such that instead 190 // of solving Ax = b, we solve the linear system AM^{-1} y = b 191 // instead, where M is such that the condition number k(AM^{-1}) 192 // is smaller than the conditioner k(A). Given the solution to 193 // this system, x = M^{-1} y. The iterative solver takes care of 194 // the mechanics of solving the preconditioned system and 195 // returning the corrected solution x. The user only needs to 196 // supply a linear operator. 197 // 198 // A null preconditioner is equivalent to an identity matrix being 199 // used a preconditioner. 200 LinearOperator* preconditioner; 201 202 203 // The following tolerance related options only makes sense for 204 // iterative solvers. Direct solvers ignore them. 205 206 // Solver terminates when 207 // 208 // |Ax - b| <= r_tolerance * |b|. 209 // 210 // This is the most commonly used termination criterion for 211 // iterative solvers. 212 double r_tolerance; 213 214 // For PSD matrices A, let 215 // 216 // Q(x) = x'Ax - 2b'x 217 // 218 // be the cost of the quadratic function defined by A and b. Then, 219 // the solver terminates at iteration i if 220 // 221 // i * (Q(x_i) - Q(x_i-1)) / Q(x_i) < q_tolerance. 222 // 223 // This termination criterion is more useful when using CG to 224 // solve the Newton step. This particular convergence test comes 225 // from Stephen Nash's work on truncated Newton 226 // methods. References: 227 // 228 // 1. Stephen G. Nash & Ariela Sofer, Assessing A Search 229 // Direction Within A Truncated Newton Method, Operation 230 // Research Letters 9(1990) 219-221. 231 // 232 // 2. Stephen G. Nash, A Survey of Truncated Newton Methods, 233 // Journal of Computational and Applied Mathematics, 234 // 124(1-2), 45-59, 2000. 235 // 236 double q_tolerance; 237 }; 238 239 // Summary of a call to the Solve method. We should move away from 240 // the true/false method for determining solver success. We should 241 // let the summary object do the talking. 242 struct Summary { 243 Summary() 244 : residual_norm(0.0), 245 num_iterations(-1), 246 termination_type(FAILURE) { 247 } 248 249 double residual_norm; 250 int num_iterations; 251 LinearSolverTerminationType termination_type; 252 }; 253 254 virtual ~LinearSolver(); 255 256 // Solve Ax = b. 257 virtual Summary Solve(LinearOperator* A, 258 const double* b, 259 const PerSolveOptions& per_solve_options, 260 double* x) = 0; 261 262 // The following two methods return copies instead of references so 263 // that the base class implementation does not have to worry about 264 // life time issues. Further, these calls are not expected to be 265 // frequent or performance sensitive. 266 virtual map<string, int> CallStatistics() const { 267 return map<string, int>(); 268 } 269 270 virtual map<string, double> TimeStatistics() const { 271 return map<string, double>(); 272 } 273 274 // Factory 275 static LinearSolver* Create(const Options& options); 276}; 277 278// This templated subclass of LinearSolver serves as a base class for 279// other linear solvers that depend on the particular matrix layout of 280// the underlying linear operator. For example some linear solvers 281// need low level access to the TripletSparseMatrix implementing the 282// LinearOperator interface. This class hides those implementation 283// details behind a private virtual method, and has the Solve method 284// perform the necessary upcasting. 285template <typename MatrixType> 286class TypedLinearSolver : public LinearSolver { 287 public: 288 virtual ~TypedLinearSolver() {} 289 virtual LinearSolver::Summary Solve( 290 LinearOperator* A, 291 const double* b, 292 const LinearSolver::PerSolveOptions& per_solve_options, 293 double* x) { 294 ScopedExecutionTimer total_time("LinearSolver::Solve", &execution_summary_); 295 CHECK_NOTNULL(A); 296 CHECK_NOTNULL(b); 297 CHECK_NOTNULL(x); 298 return SolveImpl(down_cast<MatrixType*>(A), b, per_solve_options, x); 299 } 300 301 virtual map<string, int> CallStatistics() const { 302 return execution_summary_.calls(); 303 } 304 305 virtual map<string, double> TimeStatistics() const { 306 return execution_summary_.times(); 307 } 308 309 private: 310 virtual LinearSolver::Summary SolveImpl( 311 MatrixType* A, 312 const double* b, 313 const LinearSolver::PerSolveOptions& per_solve_options, 314 double* x) = 0; 315 316 ExecutionSummary execution_summary_; 317}; 318 319// Linear solvers that depend on acccess to the low level structure of 320// a SparseMatrix. 321typedef TypedLinearSolver<BlockSparseMatrix> BlockSparseMatrixSolver; // NOLINT 322typedef TypedLinearSolver<CompressedRowSparseMatrix> CompressedRowSparseMatrixSolver; // NOLINT 323typedef TypedLinearSolver<DenseSparseMatrix> DenseSparseMatrixSolver; // NOLINT 324typedef TypedLinearSolver<TripletSparseMatrix> TripletSparseMatrixSolver; // NOLINT 325 326} // namespace internal 327} // namespace ceres 328 329#endif // CERES_INTERNAL_LINEAR_SOLVER_H_ 330