1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
|
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2015 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
#ifndef CERES_INTERNAL_SCHUR_COMPLEMENT_SOLVER_H_
#define CERES_INTERNAL_SCHUR_COMPLEMENT_SOLVER_H_
#include <memory>
#include <set>
#include <utility>
#include <vector>
#include "ceres/block_random_access_diagonal_matrix.h"
#include "ceres/block_random_access_matrix.h"
#include "ceres/block_sparse_matrix.h"
#include "ceres/block_structure.h"
#include "ceres/internal/port.h"
#include "ceres/linear_solver.h"
#include "ceres/schur_eliminator.h"
#include "ceres/types.h"
#ifdef CERES_USE_EIGEN_SPARSE
#include "Eigen/OrderingMethods"
#include "Eigen/SparseCholesky"
#endif
namespace ceres {
namespace internal {
class BlockSparseMatrix;
class SparseCholesky;
// Base class for Schur complement based linear least squares
// solvers. It assumes that the input linear system Ax = b can be
// partitioned into
//
// E y + F z = b
//
// Where x = [y;z] is a partition of the variables. The paritioning
// of the variables is such that, E'E is a block diagonal
// matrix. Further, the rows of A are ordered so that for every
// variable block in y, all the rows containing that variable block
// occur as a vertically contiguous block. i.e the matrix A looks like
//
// E F
// A = [ y1 0 0 0 | z1 0 0 0 z5]
// [ y1 0 0 0 | z1 z2 0 0 0]
// [ 0 y2 0 0 | 0 0 z3 0 0]
// [ 0 0 y3 0 | z1 z2 z3 z4 z5]
// [ 0 0 y3 0 | z1 0 0 0 z5]
// [ 0 0 0 y4 | 0 0 0 0 z5]
// [ 0 0 0 y4 | 0 z2 0 0 0]
// [ 0 0 0 y4 | 0 0 0 0 0]
// [ 0 0 0 0 | z1 0 0 0 0]
// [ 0 0 0 0 | 0 0 z3 z4 z5]
//
// This structure should be reflected in the corresponding
// CompressedRowBlockStructure object associated with A. The linear
// system Ax = b should either be well posed or the array D below
// should be non-null and the diagonal matrix corresponding to it
// should be non-singular.
//
// SchurComplementSolver has two sub-classes.
//
// DenseSchurComplementSolver: For problems where the Schur complement
// matrix is small and dense, or if CHOLMOD/SuiteSparse is not
// installed. For structure from motion problems, this is solver can
// be used for problems with upto a few hundred cameras.
//
// SparseSchurComplementSolver: For problems where the Schur
// complement matrix is large and sparse. It requires that Ceres be
// build with at least one sparse linear algebra library, as it
// computes a sparse Cholesky factorization of the Schur complement.
//
// This solver can be used for solving structure from motion problems
// with tens of thousands of cameras, though depending on the exact
// sparsity structure, it maybe better to use an iterative solver.
//
// The two solvers can be instantiated by calling
// LinearSolver::CreateLinearSolver with LinearSolver::Options::type
// set to DENSE_SCHUR and SPARSE_SCHUR
// respectively. LinearSolver::Options::elimination_groups[0] should
// be at least 1.
class CERES_EXPORT_INTERNAL SchurComplementSolver
: public BlockSparseMatrixSolver {
public:
explicit SchurComplementSolver(const LinearSolver::Options& options)
: options_(options) {
CHECK_GT(options.elimination_groups.size(), 1);
CHECK_GT(options.elimination_groups[0], 0);
CHECK(options.context != NULL);
}
SchurComplementSolver(const SchurComplementSolver&) = delete;
void operator=(const SchurComplementSolver&) = delete;
// LinearSolver methods
virtual ~SchurComplementSolver() {}
LinearSolver::Summary SolveImpl(
BlockSparseMatrix* A,
const double* b,
const LinearSolver::PerSolveOptions& per_solve_options,
double* x) override;
protected:
const LinearSolver::Options& options() const { return options_; }
const BlockRandomAccessMatrix* lhs() const { return lhs_.get(); }
void set_lhs(BlockRandomAccessMatrix* lhs) { lhs_.reset(lhs); }
const double* rhs() const { return rhs_.get(); }
void set_rhs(double* rhs) { rhs_.reset(rhs); }
private:
virtual void InitStorage(const CompressedRowBlockStructure* bs) = 0;
virtual LinearSolver::Summary SolveReducedLinearSystem(
const LinearSolver::PerSolveOptions& per_solve_options,
double* solution) = 0;
LinearSolver::Options options_;
std::unique_ptr<SchurEliminatorBase> eliminator_;
std::unique_ptr<BlockRandomAccessMatrix> lhs_;
std::unique_ptr<double[]> rhs_;
};
// Dense Cholesky factorization based solver.
class DenseSchurComplementSolver : public SchurComplementSolver {
public:
explicit DenseSchurComplementSolver(const LinearSolver::Options& options)
: SchurComplementSolver(options) {}
DenseSchurComplementSolver(const DenseSchurComplementSolver&) = delete;
void operator=(const DenseSchurComplementSolver&) = delete;
virtual ~DenseSchurComplementSolver() {}
private:
void InitStorage(const CompressedRowBlockStructure* bs) final;
LinearSolver::Summary SolveReducedLinearSystem(
const LinearSolver::PerSolveOptions& per_solve_options,
double* solution) final;
};
// Sparse Cholesky factorization based solver.
class SparseSchurComplementSolver : public SchurComplementSolver {
public:
explicit SparseSchurComplementSolver(const LinearSolver::Options& options);
SparseSchurComplementSolver(const SparseSchurComplementSolver&) = delete;
void operator=(const SparseSchurComplementSolver&) = delete;
virtual ~SparseSchurComplementSolver();
private:
void InitStorage(const CompressedRowBlockStructure* bs) final;
LinearSolver::Summary SolveReducedLinearSystem(
const LinearSolver::PerSolveOptions& per_solve_options,
double* solution) final;
LinearSolver::Summary SolveReducedLinearSystemUsingConjugateGradients(
const LinearSolver::PerSolveOptions& per_solve_options, double* solution);
// Size of the blocks in the Schur complement.
std::vector<int> blocks_;
std::unique_ptr<SparseCholesky> sparse_cholesky_;
std::unique_ptr<BlockRandomAccessDiagonalMatrix> preconditioner_;
};
} // namespace internal
} // namespace ceres
#endif // CERES_INTERNAL_SCHUR_COMPLEMENT_SOLVER_H_
|
