1. Remove constant_sparsity from LinearSolver::Options. It introduces
unnecessarily complexity in the structure of linear solvers and preconditioners.
This is the first step towards cleaning up the Preconditioner interface.
2. Minor tweaks and cleanups to the various linear solvers.
diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc
index 95edf53..e006e10 100644
--- a/internal/ceres/compressed_row_sparse_matrix.cc
+++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -40,7 +40,13 @@
namespace {
// Helper functor used by the constructor for reordering the contents
-// of a TripletSparseMatrix.
+// of a TripletSparseMatrix. This comparator assumes thay there are no
+// duplicates in the pair of arrays rows and cols, i.e., there is no
+// indices i and j (not equal to each other) s.t.
+//
+// rows[i] == rows[j] && cols[i] == cols[j]
+//
+// If this is the case, this functor will not be a StrictWeakOrdering.
struct RowColLessThan {
RowColLessThan(const int* rows, const int* cols)
: rows(rows), cols(cols) {
diff --git a/internal/ceres/dense_qr_solver.cc b/internal/ceres/dense_qr_solver.cc
index 3285054..21c5fb5 100644
--- a/internal/ceres/dense_qr_solver.cc
+++ b/internal/ceres/dense_qr_solver.cc
@@ -62,7 +62,7 @@
}
// rhs = [b;0] to account for the additional rows in the lhs.
- Vector rhs(num_rows + ((per_solve_options.D !=NULL) ? num_cols : 0));
+ Vector rhs(num_rows + ((per_solve_options.D != NULL) ? num_cols : 0));
rhs.setZero();
rhs.head(num_rows) = ConstVectorRef(b, num_rows);
diff --git a/internal/ceres/detect_structure.h b/internal/ceres/detect_structure.h
index 8af4f23..5f8e1b4 100644
--- a/internal/ceres/detect_structure.h
+++ b/internal/ceres/detect_structure.h
@@ -49,7 +49,7 @@
// is known as compile time.
//
// For more details about e_blocks and f_blocks, see
-// schur_complement.h. This information is used to initialized an
+// schur_eliminator.h. This information is used to initialized an
// appropriate template specialization of SchurEliminator.
void DetectStructure(const CompressedRowBlockStructure& bs,
const int num_eliminate_blocks,
diff --git a/internal/ceres/implicit_schur_complement.cc b/internal/ceres/implicit_schur_complement.cc
index bd90884..b09d253 100644
--- a/internal/ceres/implicit_schur_complement.cc
+++ b/internal/ceres/implicit_schur_complement.cc
@@ -42,10 +42,8 @@
namespace internal {
ImplicitSchurComplement::ImplicitSchurComplement(int num_eliminate_blocks,
- bool constant_sparsity,
bool preconditioner)
: num_eliminate_blocks_(num_eliminate_blocks),
- constant_sparsity_(constant_sparsity),
preconditioner_(preconditioner),
A_(NULL),
D_(NULL),
@@ -62,7 +60,7 @@
const double* b) {
// Since initialization is reasonably heavy, perhaps we can save on
// constructing a new object everytime.
- if ((A_ == NULL) || !constant_sparsity_) {
+ if (A_ == NULL) {
A_.reset(new PartitionedMatrixView(A, num_eliminate_blocks_));
}
@@ -71,7 +69,7 @@
// Initialize temporary storage and compute the block diagonals of
// E'E and F'E.
- if ((!constant_sparsity_) || (block_diagonal_EtE_inverse_ == NULL)) {
+ if (block_diagonal_EtE_inverse_ == NULL) {
block_diagonal_EtE_inverse_.reset(A_->CreateBlockDiagonalEtE());
if (preconditioner_) {
block_diagonal_FtF_inverse_.reset(A_->CreateBlockDiagonalFtF());
@@ -92,17 +90,10 @@
// The block diagonals of the augmented linear system contain
// contributions from the diagonal D if it is non-null. Add that to
// the block diagonals and invert them.
- if (D_ != NULL) {
- AddDiagonalAndInvert(D_, block_diagonal_EtE_inverse_.get());
- if (preconditioner_) {
- AddDiagonalAndInvert(D_ + A_->num_cols_e(),
- block_diagonal_FtF_inverse_.get());
- }
- } else {
- AddDiagonalAndInvert(NULL, block_diagonal_EtE_inverse_.get());
- if (preconditioner_) {
- AddDiagonalAndInvert(NULL, block_diagonal_FtF_inverse_.get());
- }
+ AddDiagonalAndInvert(D_, block_diagonal_EtE_inverse_.get());
+ if (preconditioner_) {
+ AddDiagonalAndInvert((D_ == NULL) ? NULL : D_ + A_->num_cols_e(),
+ block_diagonal_FtF_inverse_.get());
}
// Compute the RHS of the Schur complement system.
diff --git a/internal/ceres/implicit_schur_complement.h b/internal/ceres/implicit_schur_complement.h
index 37a319f..b9ebaa4 100644
--- a/internal/ceres/implicit_schur_complement.h
+++ b/internal/ceres/implicit_schur_complement.h
@@ -91,20 +91,13 @@
// num_eliminate_blocks is the number of E blocks in the matrix
// A.
//
- // constant_sparsity indicates if across calls to Init, the sparsity
- // structure of the matrix A remains constant or not. This makes for
- // significant savings across multiple matrices A, e.g. when used in
- // conjunction with an optimization algorithm.
- //
// 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.
- ImplicitSchurComplement(int num_eliminate_blocks,
- bool constant_sparsity,
- bool preconditioner);
+ ImplicitSchurComplement(int num_eliminate_blocks, bool preconditioner);
virtual ~ImplicitSchurComplement();
// Initialize the Schur complement for a linear least squares
@@ -151,7 +144,6 @@
void UpdateRhs();
int num_eliminate_blocks_;
- bool constant_sparsity_;
bool preconditioner_;
scoped_ptr<PartitionedMatrixView> A_;
diff --git a/internal/ceres/implicit_schur_complement_test.cc b/internal/ceres/implicit_schur_complement_test.cc
index e5dd471..9cf17af 100644
--- a/internal/ceres/implicit_schur_complement_test.cc
+++ b/internal/ceres/implicit_schur_complement_test.cc
@@ -83,7 +83,6 @@
const int num_schur_rows = blhs.num_rows();
LinearSolver::Options options;
- options.constant_sparsity = false;
options.num_eliminate_blocks = num_eliminate_blocks_;
options.type = DENSE_SCHUR;
@@ -121,7 +120,7 @@
Vector reference_solution;
ReducedLinearSystemAndSolution(D, &lhs, &rhs, &reference_solution);
- ImplicitSchurComplement isc(num_eliminate_blocks_, false, true);
+ ImplicitSchurComplement isc(num_eliminate_blocks_, true);
isc.Init(*A_, D, b_.get());
int num_sc_cols = lhs.cols();
diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc
index 42d7e86..3d837b3 100644
--- a/internal/ceres/iterative_schur_complement_solver.cc
+++ b/internal/ceres/iterative_schur_complement_solver.cc
@@ -66,10 +66,9 @@
CHECK_NOTNULL(A->block_structure());
// Initialize a ImplicitSchurComplement object.
- if ((schur_complement_ == NULL) || (!options_.constant_sparsity)) {
+ if (schur_complement_ == NULL) {
schur_complement_.reset(
new ImplicitSchurComplement(options_.num_eliminate_blocks,
- options_.constant_sparsity,
options_.preconditioner_type == JACOBI));
}
schur_complement_->Init(*A, per_solve_options.D, b);
@@ -116,7 +115,7 @@
new VisibilityBasedPreconditioner(*A->block_structure(), options_));
}
is_preconditioner_good =
- visibility_based_preconditioner_->Compute(*A, per_solve_options.D);
+ visibility_based_preconditioner_->Update(*A, per_solve_options.D);
cg_per_solve_options.preconditioner =
visibility_based_preconditioner_.get();
break;
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h
index 5860ecc..299a541 100644
--- a/internal/ceres/linear_solver.h
+++ b/internal/ceres/linear_solver.h
@@ -55,10 +55,11 @@
// Ax = b
//
// It is expected that a single instance of a LinearSolver object
-// maybe used multiple times for solving different linear
-// systems. This allows them to cache and reuse information across
-// solves if for example the sparsity of the linear system remains
-// constant.
+// maybe used multiple times for solving multiple linear systems with
+// the same sparsity structure. This allows them to cache and reuse
+// information across solves. This means that calling Solve on the
+// same LinearSolver instance with two different linear systems will
+// result in undefined behaviour.
//
// Subclasses of LinearSolver use two structs to configure themselves.
// The Options struct configures the LinearSolver object for its
@@ -73,7 +74,6 @@
min_num_iterations(1),
max_num_iterations(1),
num_threads(1),
- constant_sparsity(false),
num_eliminate_blocks(0),
residual_reset_period(10),
row_block_size(Dynamic),
@@ -93,10 +93,6 @@
// If possible, how many threads can the solver use.
int num_threads;
- // If possible cache and reuse the symbolic factorization across
- // multiple calls.
- bool constant_sparsity;
-
// Eliminate 0 to num_eliminate_blocks - 1 from the Normal
// equations to form a schur complement. Only used by the Schur
// complement based solver. The most common use for this parameter
@@ -121,8 +117,8 @@
// It is expected that these parameters are set programmatically
// rather than manually.
//
- // Please see explicit_schur_complement_solver_impl.h for more
- // details.
+ // Please see schur_complement_solver.h and schur_eliminator.h for
+ // more details.
int row_block_size;
int e_block_size;
int f_block_size;
@@ -244,6 +240,7 @@
const PerSolveOptions& per_solve_options,
double* x) = 0;
+ // Factory
static LinearSolver* Create(const Options& options);
};
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc
index 2bc8cdd..4db7c31 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/internal/ceres/schur_complement_solver.cc
@@ -57,7 +57,7 @@
const LinearSolver::PerSolveOptions& per_solve_options,
double* x) {
const time_t start_time = time(NULL);
- if (!options_.constant_sparsity || (eliminator_.get() == NULL)) {
+ if (eliminator_.get() == NULL) {
InitStorage(A->block_structure());
DetectStructure(*A->block_structure(),
options_.num_eliminate_blocks,
@@ -262,13 +262,6 @@
ss_.Free(cholmod_rhs);
cholmod_rhs = NULL;
- // If sparsity is not constant across calls, then reset the symbolic
- // factorization.
- if (!options().constant_sparsity) {
- ss_.Free(symbolic_factor_);
- symbolic_factor_ = NULL;
- }
-
if (cholmod_solution == NULL) {
LOG(ERROR) << "CHOLMOD solve failed.";
return false;
diff --git a/internal/ceres/schur_complement_solver.h b/internal/ceres/schur_complement_solver.h
index 039bc09..7cd019a 100644
--- a/internal/ceres/schur_complement_solver.h
+++ b/internal/ceres/schur_complement_solver.h
@@ -159,7 +159,7 @@
SuiteSparse ss_;
// Symbolic factorization of the reduced linear system. Precomputed
- // once and reused if constant_sparsity_ is true.
+ // once and reused in subsequent calls.
cholmod_factor* symbolic_factor_;
DISALLOW_COPY_AND_ASSIGN(SparseSchurComplementSolver);
};
diff --git a/internal/ceres/schur_complement_solver_test.cc b/internal/ceres/schur_complement_solver_test.cc
index 16c4200..5e0cc40 100644
--- a/internal/ceres/schur_complement_solver_test.cc
+++ b/internal/ceres/schur_complement_solver_test.cc
@@ -65,8 +65,8 @@
sol_d.reset(new double[num_cols]);
LinearSolver::Options options;
- options.constant_sparsity = false;
options.type = DENSE_QR;
+
scoped_ptr<LinearSolver> qr(LinearSolver::Create(options));
TripletSparseMatrix triplet_A(A->num_rows(),
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc
index ed07d9d..bb59b5f 100644
--- a/internal/ceres/solver_impl.cc
+++ b/internal/ceres/solver_impl.cc
@@ -455,7 +455,6 @@
}
LinearSolver::Options linear_solver_options;
- linear_solver_options.constant_sparsity = true;
linear_solver_options.min_num_iterations =
options->linear_solver_min_num_iterations;
linear_solver_options.max_num_iterations =
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index 59222dc..d4c2f3d 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -100,11 +100,6 @@
A->DeleteRows(num_cols);
}
- if (!options_.constant_sparsity) {
- ss_.Free(symbolic_factor_);
- symbolic_factor_ = NULL;
- }
-
summary.num_iterations = 1;
if (sol != NULL) {
memcpy(x, sol->x, num_cols * sizeof(*x));
diff --git a/internal/ceres/symmetric_linear_solver_test.cc b/internal/ceres/symmetric_linear_solver_test.cc
index 8a7ce68..f33adb4 100644
--- a/internal/ceres/symmetric_linear_solver_test.cc
+++ b/internal/ceres/symmetric_linear_solver_test.cc
@@ -63,7 +63,6 @@
LinearSolver::Options options;
options.max_num_iterations = 10;
- options.constant_sparsity = false;
LinearSolver::PerSolveOptions per_solve_options;
per_solve_options.r_tolerance = 1e-9;
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index aca7752..a930595 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -341,7 +341,6 @@
LinearSolver::Options eliminator_options;
eliminator_options.num_eliminate_blocks = options_.num_eliminate_blocks;
eliminator_options.num_threads = options_.num_threads;
- eliminator_options.constant_sparsity = true;
DetectStructure(bs, options_.num_eliminate_blocks,
&eliminator_options.row_block_size,
@@ -352,9 +351,9 @@
eliminator_->Init(options_.num_eliminate_blocks, &bs);
}
-// Compute the values of the preconditioner matrix and factorize it.
-bool VisibilityBasedPreconditioner::Compute(const BlockSparseMatrixBase& A,
- const double* D) {
+// Update the values of the preconditioner matrix and factorize it.
+bool VisibilityBasedPreconditioner::Update(const BlockSparseMatrixBase& A,
+ const double* D) {
const time_t start_time = time(NULL);
const int num_rows = m_->num_rows();
CHECK_GT(num_rows, 0);
diff --git a/internal/ceres/visibility_based_preconditioner.h b/internal/ceres/visibility_based_preconditioner.h
index fa095ca..5f14842 100644
--- a/internal/ceres/visibility_based_preconditioner.h
+++ b/internal/ceres/visibility_based_preconditioner.h
@@ -133,7 +133,7 @@
// options.num_eliminate_blocks = num_points;
// VisibilityBasedPreconditioner preconditioner(
// *A.block_structure(), options);
-// preconditioner.Compute(A, NULL);
+// preconditioner.Update(A, NULL);
// preconditioner.RightMultiply(x, y);
//
@@ -160,7 +160,7 @@
const LinearSolver::Options& options);
virtual ~VisibilityBasedPreconditioner();
- // Compute the numerical value of the preconditioner for the linear
+ // Update the numerical value of the preconditioner for the linear
// system:
//
// | A | x = |b|
@@ -171,12 +171,12 @@
//
// D can be NULL, in which case its interpreted as a diagonal matrix
// of size zero.
- bool Compute(const BlockSparseMatrixBase& A,
- const double* D);
+ bool Update(const BlockSparseMatrixBase& A, const double* D);
+
// LinearOperator interface. Since the operator is symmetric,
// LeftMultiply and num_cols are just calls to RightMultiply and
- // num_rows respectively. Compute() must be called before
+ // num_rows respectively. Update() must be called before
// RightMultiply can be called.
virtual void RightMultiply(const double* x, double* y) const;
virtual void LeftMultiply(const double* x, double* y) const {
diff --git a/internal/ceres/visibility_based_preconditioner_test.cc b/internal/ceres/visibility_based_preconditioner_test.cc
index 390f081..244d0b0 100644
--- a/internal/ceres/visibility_based_preconditioner_test.cc
+++ b/internal/ceres/visibility_based_preconditioner_test.cc
@@ -132,7 +132,7 @@
}
AssertionResult PreconditionerValuesMatch() {
- preconditioner_->Compute(*A_, D_.get());
+ preconditioner_->Update(*A_, D_.get());
const HashSet<pair<int, int> >& cluster_pairs = get_cluster_pairs();
const BlockRandomAccessSparseMatrix* m = get_m();
Matrix preconditioner_matrix;
