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