blob: bdf7b569891d90f900b90714e2fa0891e2fbbcba [file] [log] [blame]
Keir Mierle8ebb0732012-04-30 23:09:08 -07001// 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 Agarwal509f68c2013-02-20 01:39:03 -080038#include <map>
Sameer Agarwal0c52f1e2012-09-17 11:30:14 -070039#include <vector>
Keir Mierle8ebb0732012-04-30 23:09:08 -070040#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 Agarwalbdd87c02013-01-29 16:24:31 -080044#include "ceres/execution_summary.h"
Keir Mierle8ebb0732012-04-30 23:09:08 -070045#include "ceres/triplet_sparse_matrix.h"
46#include "ceres/types.h"
Sameer Agarwal509f68c2013-02-20 01:39:03 -080047#include "glog/logging.h"
Keir Mierle8ebb0732012-04-30 23:09:08 -070048
49namespace ceres {
50namespace internal {
51
52class 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 Agarwala9d8ef82012-05-14 02:28:05 -070060// 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 Mierle8ebb0732012-04-30 23:09:08 -070065//
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.
70class LinearSolver {
71 public:
72 struct Options {
73 Options()
74 : type(SPARSE_NORMAL_CHOLESKY),
75 preconditioner_type(JACOBI),
Sameer Agarwalb0518732012-05-29 00:27:57 -070076 sparse_linear_algebra_library(SUITE_SPARSE),
Sameer Agarwal9189f4e2013-04-19 17:09:49 -070077 use_postordering(false),
Keir Mierle8ebb0732012-04-30 23:09:08 -070078 min_num_iterations(1),
79 max_num_iterations(1),
80 num_threads(1),
Keir Mierle8ebb0732012-04-30 23:09:08 -070081 residual_reset_period(10),
Sameer Agarwal31730ef2013-02-28 11:20:28 -080082 row_block_size(Eigen::Dynamic),
83 e_block_size(Eigen::Dynamic),
84 f_block_size(Eigen::Dynamic) {
Keir Mierle8ebb0732012-04-30 23:09:08 -070085 }
86
87 LinearSolverType type;
88
89 PreconditionerType preconditioner_type;
90
Sameer Agarwalb0518732012-05-29 00:27:57 -070091 SparseLinearAlgebraLibraryType sparse_linear_algebra_library;
92
Sameer Agarwal9189f4e2013-04-19 17:09:49 -070093 // See solver.h for information about this flag.
94 bool use_postordering;
Sameer Agarwal7a3c43b2012-06-05 23:10:59 -070095
Keir Mierle8ebb0732012-04-30 23:09:08 -070096 // 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 Agarwal0c52f1e2012-09-17 11:30:14 -0700104 // 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 Agarwal931c3092013-02-25 09:46:21 -0800109 // 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 Agarwal0c52f1e2012-09-17 11:30:14 -0700114 //
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 Mierle8ebb0732012-04-30 23:09:08 -0700121
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 Agarwala9d8ef82012-05-14 02:28:05 -0700136 // Please see schur_complement_solver.h and schur_eliminator.h for
137 // more details.
Keir Mierle8ebb0732012-04-30 23:09:08 -0700138 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 Mierlef7898fb2012-05-05 20:55:08 -0700169 // arg min_x ||Ax - b||^2 + ||Dx||^2
Keir Mierle8ebb0732012-04-30 23:09:08 -0700170 //
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 Agarwal42a84b82013-02-01 12:22:53 -0800259 // 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 Agarwala9d8ef82012-05-14 02:28:05 -0700271 // Factory
Keir Mierle8ebb0732012-04-30 23:09:08 -0700272 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.
282template <typename MatrixType>
283class 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 Agarwal42a84b82013-02-01 12:22:53 -0800291 ScopedExecutionTimer total_time("LinearSolver::Solve", &execution_summary_);
Keir Mierle8ebb0732012-04-30 23:09:08 -0700292 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 Agarwal42a84b82013-02-01 12:22:53 -0800298 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 Mierle8ebb0732012-04-30 23:09:08 -0700306 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 Agarwal42a84b82013-02-01 12:22:53 -0800312
313 ExecutionSummary execution_summary_;
Keir Mierle8ebb0732012-04-30 23:09:08 -0700314};
315
316// Linear solvers that depend on acccess to the low level structure of
317// a SparseMatrix.
318typedef TypedLinearSolver<BlockSparseMatrix> BlockSparseMatrixSolver; // NOLINT
Keir Mierle8ebb0732012-04-30 23:09:08 -0700319typedef TypedLinearSolver<CompressedRowSparseMatrix> CompressedRowSparseMatrixSolver; // NOLINT
320typedef TypedLinearSolver<DenseSparseMatrix> DenseSparseMatrixSolver; // NOLINT
321typedef TypedLinearSolver<TripletSparseMatrix> TripletSparseMatrixSolver; // NOLINT
322
323} // namespace internal
324} // namespace ceres
325
326#endif // CERES_INTERNAL_LINEAR_SOLVER_H_