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
|
// 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)
//
// An iterative solver for solving the Schur complement/reduced camera
// linear system that arise in SfM problems.
#ifndef CERES_INTERNAL_IMPLICIT_SCHUR_COMPLEMENT_H_
#define CERES_INTERNAL_IMPLICIT_SCHUR_COMPLEMENT_H_
#include <memory>
#include "ceres/internal/disable_warnings.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/export.h"
#include "ceres/linear_operator.h"
#include "ceres/linear_solver.h"
#include "ceres/partitioned_matrix_view.h"
#include "ceres/types.h"
namespace ceres {
namespace internal {
class BlockSparseMatrix;
// This class implements various linear algebraic operations related
// to the Schur complement without explicitly forming it.
//
//
// Given a reactangular linear system Ax = b, where
//
// A = [E F]
//
// The normal equations are given by
//
// A'Ax = A'b
//
// |E'E E'F||y| = |E'b|
// |F'E F'F||z| |F'b|
//
// and the Schur complement system is given by
//
// [F'F - F'E (E'E)^-1 E'F] z = F'b - F'E (E'E)^-1 E'b
//
// Now if we wish to solve Ax = b in the least squares sense, one way
// is to form this Schur complement system and solve it using
// Preconditioned Conjugate Gradients.
//
// The key operation in a conjugate gradient solver is the evaluation of the
// matrix vector product with the Schur complement
//
// S = F'F - F'E (E'E)^-1 E'F
//
// It is straightforward to see that matrix vector products with S can
// be evaluated without storing S in memory. Instead, given (E'E)^-1
// (which for our purposes is an easily inverted block diagonal
// matrix), it can be done in terms of matrix vector products with E,
// F and (E'E)^-1. This class implements this functionality and other
// auxilliary bits needed to implement a CG solver on the Schur
// complement using the PartitionedMatrixView object.
//
// THREAD SAFETY: This class is nqot thread safe. In particular, the
// RightMultiply (and the LeftMultiply) methods are not thread safe as
// they depend on mutable arrays used for the temporaries needed to
// compute the product y += Sx;
class CERES_NO_EXPORT ImplicitSchurComplement final : public LinearOperator {
public:
// num_eliminate_blocks is the number of E blocks in the matrix
// A.
//
// preconditioner indicates whether the inverse of the matrix F'F
// should be computed or not as a preconditioner for the Schur
// Complement.
//
// TODO(sameeragarwal): Get rid of the two bools below and replace
// them with enums.
explicit ImplicitSchurComplement(const LinearSolver::Options& options);
// Initialize the Schur complement for a linear least squares
// problem of the form
//
// |A | x = |b|
// |diag(D)| |0|
//
// If D is null, then it is treated as a zero dimensional matrix. It
// is important that the matrix A have a BlockStructure object
// associated with it and has a block structure that is compatible
// with the SchurComplement solver.
void Init(const BlockSparseMatrix& A, const double* D, const double* b);
// y += Sx, where S is the Schur complement.
void RightMultiply(const double* x, double* y) const final;
// The Schur complement is a symmetric positive definite matrix,
// thus the left and right multiply operators are the same.
void LeftMultiply(const double* x, double* y) const final {
RightMultiply(x, y);
}
// y = (E'E)^-1 (E'b - E'F x). Given an estimate of the solution to
// the Schur complement system, this method computes the value of
// the e_block variables that were eliminated to form the Schur
// complement.
void BackSubstitute(const double* x, double* y);
int num_rows() const final { return A_->num_cols_f(); }
int num_cols() const final { return A_->num_cols_f(); }
const Vector& rhs() const { return rhs_; }
const BlockSparseMatrix* block_diagonal_EtE_inverse() const {
return block_diagonal_EtE_inverse_.get();
}
const BlockSparseMatrix* block_diagonal_FtF_inverse() const {
return block_diagonal_FtF_inverse_.get();
}
private:
void AddDiagonalAndInvert(const double* D, BlockSparseMatrix* matrix);
void UpdateRhs();
const LinearSolver::Options& options_;
std::unique_ptr<PartitionedMatrixViewBase> A_;
const double* D_;
const double* b_;
std::unique_ptr<BlockSparseMatrix> block_diagonal_EtE_inverse_;
std::unique_ptr<BlockSparseMatrix> block_diagonal_FtF_inverse_;
Vector rhs_;
// Temporary storage vectors used to implement RightMultiply.
mutable Vector tmp_rows_;
mutable Vector tmp_e_cols_;
mutable Vector tmp_e_cols_2_;
mutable Vector tmp_f_cols_;
};
} // namespace internal
} // namespace ceres
#include "ceres/internal/reenable_warnings.h"
#endif // CERES_INTERNAL_IMPLICIT_SCHUR_COMPLEMENT_H_
|
