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/
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_
39
#include <glog/logging.h>
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"
44
#include "ceres/triplet_sparse_matrix.h"
45
#include "ceres/types.h"
52
// Abstract base class for objects that implement algorithms for
53
// solving linear systems
57
// It is expected that a single instance of a LinearSolver object
58
// maybe used multiple times for solving multiple linear systems with
59
// the same sparsity structure. This allows them to cache and reuse
60
// information across solves. This means that calling Solve on the
61
// same LinearSolver instance with two different linear systems will
62
// result in undefined behaviour.
64
// Subclasses of LinearSolver use two structs to configure themselves.
65
// The Options struct configures the LinearSolver object for its
66
// lifetime. The PerSolveOptions struct is used to specify options for
67
// a particular Solve call.
72
: type(SPARSE_NORMAL_CHOLESKY),
73
preconditioner_type(JACOBI),
74
min_num_iterations(1),
75
max_num_iterations(1),
77
num_eliminate_blocks(0),
78
residual_reset_period(10),
79
row_block_size(Dynamic),
80
e_block_size(Dynamic),
81
f_block_size(Dynamic) {
84
LinearSolverType type;
86
PreconditionerType preconditioner_type;
88
// Number of internal iterations that the solver uses. This
89
// parameter only makes sense for iterative solvers like CG.
90
int min_num_iterations;
91
int max_num_iterations;
93
// If possible, how many threads can the solver use.
96
// Eliminate 0 to num_eliminate_blocks - 1 from the Normal
97
// equations to form a schur complement. Only used by the Schur
98
// complement based solver. The most common use for this parameter
99
// is in the case of structure from motion problems where we have
100
// camera blocks and point blocks. Then setting the
101
// num_eliminate_blocks to the number of points allows the solver
102
// to use the Schur complement trick. For more details see the
103
// description of this parameter in solver.h.
104
int num_eliminate_blocks;
106
// Iterative solvers, e.g. Preconditioned Conjugate Gradients
107
// maintain a cheap estimate of the residual which may become
108
// inaccurate over time. Thus for non-zero values of this
109
// parameter, the solver can be told to recalculate the value of
110
// the residual using a |b - Ax| evaluation.
111
int residual_reset_period;
113
// If the block sizes in a BlockSparseMatrix are fixed, then in
114
// some cases the Schur complement based solvers can detect and
115
// specialize on them.
117
// It is expected that these parameters are set programmatically
118
// rather than manually.
120
// Please see schur_complement_solver.h and schur_eliminator.h for
127
// Options for the Solve method.
128
struct PerSolveOptions {
131
preconditioner(NULL),
136
// This option only makes sense for unsymmetric linear solvers
137
// that can solve rectangular linear systems.
139
// Given a matrix A, an optional diagonal matrix D as a vector,
140
// and a vector b, the linear solver will solve for
145
// If D is null, then it is treated as zero, and the solver returns
150
// In either case, x is the vector that solves the following
151
// optimization problem.
153
// arg min_x ||Ax - b||^2 + ||Dx||^2
155
// Here A is a matrix of size m x n, with full column rank. If A
156
// does not have full column rank, the results returned by the
157
// solver cannot be relied on. D, if it is not null is an array of
158
// size n. b is an array of size m and x is an array of size n.
161
// This option only makes sense for iterative solvers.
163
// In general the performance of an iterative linear solver
164
// depends on the condition number of the matrix A. For example
165
// the convergence rate of the conjugate gradients algorithm
166
// is proportional to the square root of the condition number.
168
// One particularly useful technique for improving the
169
// conditioning of a linear system is to precondition it. In its
170
// simplest form a preconditioner is a matrix M such that instead
171
// of solving Ax = b, we solve the linear system AM^{-1} y = b
172
// instead, where M is such that the condition number k(AM^{-1})
173
// is smaller than the conditioner k(A). Given the solution to
174
// this system, x = M^{-1} y. The iterative solver takes care of
175
// the mechanics of solving the preconditioned system and
176
// returning the corrected solution x. The user only needs to
177
// supply a linear operator.
179
// A null preconditioner is equivalent to an identity matrix being
180
// used a preconditioner.
181
LinearOperator* preconditioner;
184
// The following tolerance related options only makes sense for
185
// iterative solvers. Direct solvers ignore them.
187
// Solver terminates when
189
// |Ax - b| <= r_tolerance * |b|.
191
// This is the most commonly used termination criterion for
192
// iterative solvers.
195
// For PSD matrices A, let
197
// Q(x) = x'Ax - 2b'x
199
// be the cost of the quadratic function defined by A and b. Then,
200
// the solver terminates at iteration i if
202
// i * (Q(x_i) - Q(x_i-1)) / Q(x_i) < q_tolerance.
204
// This termination criterion is more useful when using CG to
205
// solve the Newton step. This particular convergence test comes
206
// from Stephen Nash's work on truncated Newton
207
// methods. References:
209
// 1. Stephen G. Nash & Ariela Sofer, Assessing A Search
210
// Direction Within A Truncated Newton Method, Operation
211
// Research Letters 9(1990) 219-221.
213
// 2. Stephen G. Nash, A Survey of Truncated Newton Methods,
214
// Journal of Computational and Applied Mathematics,
215
// 124(1-2), 45-59, 2000.
220
// Summary of a call to the Solve method. We should move away from
221
// the true/false method for determining solver success. We should
222
// let the summary object do the talking.
225
: residual_norm(0.0),
227
termination_type(FAILURE) {
230
double residual_norm;
232
LinearSolverTerminationType termination_type;
235
virtual ~LinearSolver();
238
virtual Summary Solve(LinearOperator* A,
240
const PerSolveOptions& per_solve_options,
244
static LinearSolver* Create(const Options& options);
247
// This templated subclass of LinearSolver serves as a base class for
248
// other linear solvers that depend on the particular matrix layout of
249
// the underlying linear operator. For example some linear solvers
250
// need low level access to the TripletSparseMatrix implementing the
251
// LinearOperator interface. This class hides those implementation
252
// details behind a private virtual method, and has the Solve method
253
// perform the necessary upcasting.
254
template <typename MatrixType>
255
class TypedLinearSolver : public LinearSolver {
257
virtual ~TypedLinearSolver() {}
258
virtual LinearSolver::Summary Solve(
261
const LinearSolver::PerSolveOptions& per_solve_options,
266
return SolveImpl(down_cast<MatrixType*>(A), b, per_solve_options, x);
270
virtual LinearSolver::Summary SolveImpl(
273
const LinearSolver::PerSolveOptions& per_solve_options,
277
// Linear solvers that depend on acccess to the low level structure of
279
typedef TypedLinearSolver<BlockSparseMatrix> BlockSparseMatrixSolver; // NOLINT
280
typedef TypedLinearSolver<BlockSparseMatrixBase> BlockSparseMatrixBaseSolver; // NOLINT
281
typedef TypedLinearSolver<CompressedRowSparseMatrix> CompressedRowSparseMatrixSolver; // NOLINT
282
typedef TypedLinearSolver<DenseSparseMatrix> DenseSparseMatrixSolver; // NOLINT
283
typedef TypedLinearSolver<TripletSparseMatrix> TripletSparseMatrixSolver; // NOLINT
285
} // namespace internal
288
#endif // CERES_INTERNAL_LINEAR_SOLVER_H_
added
removed
internal/ceres/linear_solver.h
