More pre-ordering support.
1. CX_SPARSE supports pre-ordering of the jacobian.
2. Add support for constrained approximate minimum degree ordering
for SuiteSparse versions >= 4.2.0
3. Using 2, support for pre-ordering for SPARSE_SCHUR when used
with SUITE_SPARSE.
4. Using 2, support for user orderings in SPARSE_NORMAL_CHOLESKY.
5. Minor cleanups in documentation and code all around.
6. Test update and refactoring.
Change-Id: Ibfe3ac95d59d54ab14d1d60a07f767688070f29f
diff --git a/internal/ceres/cxsparse.cc b/internal/ceres/cxsparse.cc
index b7f2520..3279697 100644
--- a/internal/ceres/cxsparse.cc
+++ b/internal/ceres/cxsparse.cc
@@ -93,6 +93,11 @@
return cs_schol(1, A);
}
+cs_dis* CXSparse::AnalyzeCholeskyWithNaturalOrdering(cs_di* A) {
+ // order = 0 for Natural ordering.
+ return cs_schol(0, A);
+}
+
cs_dis* CXSparse::BlockAnalyzeCholesky(cs_di* A,
const vector<int>& row_blocks,
const vector<int>& col_blocks) {
@@ -173,6 +178,20 @@
return cs_compress(&tsm_wrapper);
}
+void CXSparse::ApproximateMinimumDegreeOrdering(cs_di* A, int* ordering) {
+ int* cs_ordering = cs_amd(1, A);
+ copy(cs_ordering, cs_ordering + A->m, ordering);
+ cs_free(cs_ordering);
+}
+
+cs_di* CXSparse::TransposeMatrix(cs_di* A) {
+ return cs_di_transpose(A, 1);
+}
+
+cs_di* CXSparse::MatrixMatrixMultiply(cs_di* A, cs_di* B) {
+ return cs_di_multiply(A, B);
+}
+
void CXSparse::Free(cs_di* sparse_matrix) {
cs_di_spfree(sparse_matrix);
}
diff --git a/internal/ceres/cxsparse.h b/internal/ceres/cxsparse.h
index d34b635..6004301 100644
--- a/internal/ceres/cxsparse.h
+++ b/internal/ceres/cxsparse.h
@@ -70,14 +70,49 @@
// with Free. May return NULL if the compression or allocation fails.
cs_di* CreateSparseMatrix(TripletSparseMatrix* A);
+ // B = A'
+ //
+ // The returned matrix should be deallocated with Free when not used
+ // anymore.
+ cs_di* TransposeMatrix(cs_di* A);
+
+ // C = A * B
+ //
+ // The returned matrix should be deallocated with Free when not used
+ // anymore.
+ cs_di* MatrixMatrixMultiply(cs_di* A, cs_di* B);
+
// Computes a symbolic factorization of A that can be used in SolveCholesky.
+ //
// The returned matrix should be deallocated with Free when not used anymore.
cs_dis* AnalyzeCholesky(cs_di* A);
+ // Computes a symbolic factorization of A that can be used in
+ // SolveCholesky, but does not compute a fill-reducing ordering.
+ //
+ // The returned matrix should be deallocated with Free when not used anymore.
+ cs_dis* AnalyzeCholeskyWithNaturalOrdering(cs_di* A);
+
+ // Computes a symbolic factorization of A that can be used in
+ // SolveCholesky. The difference from AnalyzeCholesky is that this
+ // function first detects the block sparsity of the matrix using
+ // information about the row and column blocks and uses this block
+ // sparse matrix to find a fill-reducing ordering. This ordering is
+ // then used to find a symbolic factorization. This can result in a
+ // significant performance improvement AnalyzeCholesky on block
+ // sparse matrices.
+ //
+ // The returned matrix should be deallocated with Free when not used
+ // anymore.
cs_dis* BlockAnalyzeCholesky(cs_di* A,
const vector<int>& row_blocks,
const vector<int>& col_blocks);
+ // Compute an fill-reducing approximate minimum degree ordering of
+ // the matrix A. ordering should be non-NULL and should point to
+ // enough memory to hold the ordering for the rows of A.
+ void ApproximateMinimumDegreeOrdering(cs_di* A, int* ordering);
+
void Free(cs_di* sparse_matrix);
void Free(cs_dis* symbolic_factorization);
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc
index 0defcd6..09f61d7 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/internal/ceres/schur_complement_solver.cc
@@ -276,26 +276,42 @@
return true;
}
- cholmod_sparse* cholmod_lhs = ss_.CreateSparseMatrix(tsm);
- // The matrix is symmetric, and the upper triangular part of the
- // matrix contains the values.
- cholmod_lhs->stype = 1;
+ cholmod_sparse* cholmod_lhs = NULL;
+ if (options().use_postordering) {
+ // If we are going to do a full symbolic analysis of the schur
+ // complement matrix from scratch and not rely on the
+ // pre-ordering, then the fastest path in cholmod_factorize is the
+ // one corresponding to upper triangular matrices.
+
+ // Create a upper triangular symmetric matrix.
+ cholmod_lhs = ss_.CreateSparseMatrix(tsm);
+ cholmod_lhs->stype = 1;
+
+ if (factor_ == NULL) {
+ factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs, blocks_, blocks_);
+ }
+ } else {
+ // If we are going to use the natural ordering (i.e. rely on the
+ // pre-ordering computed by solver_impl.cc), then the fastest
+ // path in cholmod_factorize is the one corresponding to lower
+ // triangular matrices.
+
+ // Create a upper triangular symmetric matrix.
+ cholmod_lhs = ss_.CreateSparseMatrixTranspose(tsm);
+ cholmod_lhs->stype = -1;
+
+ if (factor_ == NULL) {
+ factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(cholmod_lhs);
+ }
+ }
cholmod_dense* cholmod_rhs =
ss_.CreateDenseVector(const_cast<double*>(rhs()), num_rows, num_rows);
-
- // Symbolic factorization is computed if we don't already have one handy.
- if (factor_ == NULL) {
- factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs, blocks_, blocks_);
- }
-
cholmod_dense* cholmod_solution =
ss_.SolveCholesky(cholmod_lhs, factor_, cholmod_rhs);
ss_.Free(cholmod_lhs);
- cholmod_lhs = NULL;
ss_.Free(cholmod_rhs);
- cholmod_rhs = NULL;
if (cholmod_solution == NULL) {
LOG(WARNING) << "CHOLMOD solve failed.";
@@ -339,7 +355,8 @@
// Compute symbolic factorization if not available.
if (cxsparse_factor_ == NULL) {
- cxsparse_factor_ = CHECK_NOTNULL(cxsparse_.BlockAnalyzeCholesky(lhs, blocks_, blocks_));
+ cxsparse_factor_ =
+ CHECK_NOTNULL(cxsparse_.BlockAnalyzeCholesky(lhs, blocks_, blocks_));
}
// Solve the linear system.
diff --git a/internal/ceres/schur_complement_solver_test.cc b/internal/ceres/schur_complement_solver_test.cc
index 1820bc9..57fd263 100644
--- a/internal/ceres/schur_complement_solver_test.cc
+++ b/internal/ceres/schur_complement_solver_test.cc
@@ -87,7 +87,8 @@
int problem_id,
bool regularization,
ceres::LinearSolverType linear_solver_type,
- ceres::SparseLinearAlgebraLibraryType sparse_linear_algebra_library) {
+ ceres::SparseLinearAlgebraLibraryType sparse_linear_algebra_library,
+ bool use_postordering) {
SetUpFromProblemId(problem_id);
LinearSolver::Options options;
options.elimination_groups.push_back(num_eliminate_blocks);
@@ -95,6 +96,7 @@
A->block_structure()->cols.size() - num_eliminate_blocks);
options.type = linear_solver_type;
options.sparse_linear_algebra_library = sparse_linear_algebra_library;
+ options.use_postordering = use_postordering;
scoped_ptr<LinearSolver> solver(LinearSolver::Create(options));
@@ -129,32 +131,49 @@
scoped_array<double> sol_d;
};
+TEST_F(SchurComplementSolverTest, DenseSchurWithSmallProblem) {
+ ComputeAndCompareSolutions(2, false, DENSE_SCHUR, SUITE_SPARSE, true);
+ ComputeAndCompareSolutions(2, true, DENSE_SCHUR, SUITE_SPARSE, true);
+}
+
+TEST_F(SchurComplementSolverTest, DenseSchurWithLargeProblem) {
+ ComputeAndCompareSolutions(3, false, DENSE_SCHUR, SUITE_SPARSE, true);
+ ComputeAndCompareSolutions(3, true, DENSE_SCHUR, SUITE_SPARSE, true);
+}
+
#ifndef CERES_NO_SUITESPARSE
-TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparse) {
- ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, SUITE_SPARSE);
- ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, SUITE_SPARSE);
- ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, SUITE_SPARSE);
- ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, SUITE_SPARSE);
+TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseSmallProblemNoPostOrdering) {
+ ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, SUITE_SPARSE, false);
+ ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, SUITE_SPARSE, false);
+}
+
+TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseSmallProblemPostOrdering) {
+ ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, SUITE_SPARSE, true);
+ ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, SUITE_SPARSE, true);
+}
+
+TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseLargeProblemNoPostOrdering) {
+ ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, SUITE_SPARSE, false);
+ ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, SUITE_SPARSE, false);
+}
+
+TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseLargeProblemPostOrdering) {
+ ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, SUITE_SPARSE, true);
+ ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, SUITE_SPARSE, true);
}
#endif // CERES_NO_SUITESPARSE
#ifndef CERES_NO_CXSPARSE
-TEST_F(SchurComplementSolverTest, SparseSchurWithCXSparse) {
- ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, CX_SPARSE);
- ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, CX_SPARSE);
- ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, CX_SPARSE);
- ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, CX_SPARSE);
+TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseSmallProblem) {
+ ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, SUITE_SPARSE, true);
+ ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, SUITE_SPARSE, true);
+}
+
+TEST_F(SchurComplementSolverTest, SparseSchurWithSuiteSparseLargeProblem) {
+ ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, SUITE_SPARSE, true);
+ ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, SUITE_SPARSE, true);
}
#endif // CERES_NO_CXSPARSE
-TEST_F(SchurComplementSolverTest, DenseSchur) {
- // The sparse linear algebra library type is ignored for
- // DENSE_SCHUR.
- ComputeAndCompareSolutions(2, false, DENSE_SCHUR, SUITE_SPARSE);
- ComputeAndCompareSolutions(3, false, DENSE_SCHUR, SUITE_SPARSE);
- ComputeAndCompareSolutions(2, true, DENSE_SCHUR, SUITE_SPARSE);
- ComputeAndCompareSolutions(3, true, DENSE_SCHUR, SUITE_SPARSE);
-}
-
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc
index 5779299..52056f7 100644
--- a/internal/ceres/solver_impl.cc
+++ b/internal/ceres/solver_impl.cc
@@ -33,7 +33,9 @@
#include <cstdio>
#include <iostream> // NOLINT
#include <numeric>
+#include <string>
#include "ceres/coordinate_descent_minimizer.h"
+#include "ceres/cxsparse.h"
#include "ceres/evaluator.h"
#include "ceres/gradient_checking_cost_function.h"
#include "ceres/iteration_callback.h"
@@ -995,7 +997,9 @@
}
if (IsSchurType(options->linear_solver_type)) {
- if (!ReorderProgramForSchurTypeLinearSolver(problem_impl->parameter_map(),
+ if (!ReorderProgramForSchurTypeLinearSolver(options->linear_solver_type,
+ options->sparse_linear_algebra_library,
+ problem_impl->parameter_map(),
linear_solver_ordering,
transformed_program.get(),
error)) {
@@ -1004,9 +1008,15 @@
return transformed_program.release();
}
- if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY &&
- options->sparse_linear_algebra_library == SUITE_SPARSE) {
- ReorderProgramForSparseNormalCholesky(transformed_program.get());
+ if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
+ if (!ReorderProgramForSparseNormalCholesky(
+ options->sparse_linear_algebra_library,
+ linear_solver_ordering,
+ transformed_program.get(),
+ error)) {
+ return NULL;
+ }
+
return transformed_program.release();
}
@@ -1093,6 +1103,18 @@
linear_solver_options.sparse_linear_algebra_library =
options->sparse_linear_algebra_library;
linear_solver_options.use_postordering = options->use_postordering;
+
+ // Ignore user's postordering preferences and force it to be true if
+ // cholmod_camd is not available. This ensures that the linear
+ // solver does not assume that a fill-reducing pre-ordering has been
+ // done.
+#if !defined(CERES_NO_SUITESPARSE) && defined(CERES_NO_CAMD)
+ if (IsSchurType(linear_solver_options.type) &&
+ linear_solver_options.sparse_linear_algebra_library == SUITE_SPARSE) {
+ linear_solver_options.use_postordering = true;
+ }
+#endif
+
linear_solver_options.num_threads = options->num_linear_solver_threads;
options->num_linear_solver_threads = linear_solver_options.num_threads;
@@ -1115,48 +1137,6 @@
return LinearSolver::Create(linear_solver_options);
}
-bool SolverImpl::ApplyUserOrdering(
- const ProblemImpl::ParameterMap& parameter_map,
- const ParameterBlockOrdering* ordering,
- Program* program,
- string* error) {
- if (ordering->NumElements() != program->NumParameterBlocks()) {
- *error = StringPrintf("User specified ordering does not have the same "
- "number of parameters as the problem. The problem"
- "has %d blocks while the ordering has %d blocks.",
- program->NumParameterBlocks(),
- ordering->NumElements());
- return false;
- }
-
- vector<ParameterBlock*>* parameter_blocks =
- program->mutable_parameter_blocks();
- parameter_blocks->clear();
-
- const map<int, set<double*> >& groups =
- ordering->group_to_elements();
-
- for (map<int, set<double*> >::const_iterator group_it = groups.begin();
- group_it != groups.end();
- ++group_it) {
- const set<double*>& group = group_it->second;
- for (set<double*>::const_iterator parameter_block_ptr_it = group.begin();
- parameter_block_ptr_it != group.end();
- ++parameter_block_ptr_it) {
- ProblemImpl::ParameterMap::const_iterator parameter_block_it =
- parameter_map.find(*parameter_block_ptr_it);
- if (parameter_block_it == parameter_map.end()) {
- *error = StringPrintf("User specified ordering contains a pointer "
- "to a double that is not a parameter block in "
- "the problem. The invalid double is in group: %d",
- group_it->first);
- return false;
- }
- parameter_blocks->push_back(parameter_block_it->second);
- }
- }
- return true;
-}
// Find the minimum index of any parameter block to the given residual.
// Parameter blocks that have indices greater than num_eliminate_blocks are
@@ -1364,64 +1344,51 @@
LOG(WARNING) << msg;
}
-bool SolverImpl::ReorderProgramForSchurTypeLinearSolver(
+bool SolverImpl::ApplyUserOrdering(
const ProblemImpl::ParameterMap& parameter_map,
- ParameterBlockOrdering* ordering,
+ const ParameterBlockOrdering* parameter_block_ordering,
Program* program,
string* error) {
- // At this point one of two things is true.
- //
- // 1. The user did not specify an ordering - ordering has one group
- // containing all the parameter blocks.
-
- // 2. The user specified an ordering, and the first group has
- // non-zero elements.
- //
- // We handle these two cases in turn.
- if (ordering->NumGroups() == 1) {
- // If the user supplied an ordering with just one
- // group, it is equivalent to the user supplying NULL as an
- // ordering. Ceres is completely free to choose the parameter
- // block ordering as it sees fit. For Schur type solvers, this
- // means that the user wishes for Ceres to identify the e_blocks,
- // which we do by computing a maximal independent set.
- vector<ParameterBlock*> schur_ordering;
- const int num_eliminate_blocks = ComputeSchurOrdering(*program,
- &schur_ordering);
-
- CHECK_EQ(schur_ordering.size(), program->NumParameterBlocks())
- << "Congratulations, you found a Ceres bug! Please report this error "
- << "to the developers.";
-
- // Update the ordering object.
- for (int i = 0; i < schur_ordering.size(); ++i) {
- double* parameter_block = schur_ordering[i]->mutable_user_state();
- const int group_id = (i < num_eliminate_blocks) ? 0 : 1;
- ordering->AddElementToGroup(parameter_block, group_id);
- }
-
- // Apply the parameter block re-ordering. Technically we could
- // call ApplyUserOrdering, but this is cheaper and simpler.
- swap(*program->mutable_parameter_blocks(), schur_ordering);
- } else {
- // The user supplied an ordering.
- if (!ApplyUserOrdering(parameter_map, ordering, program, error)) {
- return false;
- }
+ const int num_parameter_blocks = program->NumParameterBlocks();
+ if (parameter_block_ordering->NumElements() != num_parameter_blocks) {
+ *error = StringPrintf("User specified ordering does not have the same "
+ "number of parameters as the problem. The problem"
+ "has %d blocks while the ordering has %d blocks.",
+ num_parameter_blocks,
+ parameter_block_ordering->NumElements());
+ return false;
}
- program->SetParameterOffsetsAndIndex();
+ vector<ParameterBlock*>* parameter_blocks =
+ program->mutable_parameter_blocks();
+ parameter_blocks->clear();
- const int num_eliminate_blocks =
- ordering->group_to_elements().begin()->second.size();
+ const map<int, set<double*> >& groups =
+ parameter_block_ordering->group_to_elements();
- // Schur type solvers also require that their residual blocks be
- // lexicographically ordered.
- return LexicographicallyOrderResidualBlocks(num_eliminate_blocks,
- program,
- error);
+ for (map<int, set<double*> >::const_iterator group_it = groups.begin();
+ group_it != groups.end();
+ ++group_it) {
+ const set<double*>& group = group_it->second;
+ for (set<double*>::const_iterator parameter_block_ptr_it = group.begin();
+ parameter_block_ptr_it != group.end();
+ ++parameter_block_ptr_it) {
+ ProblemImpl::ParameterMap::const_iterator parameter_block_it =
+ parameter_map.find(*parameter_block_ptr_it);
+ if (parameter_block_it == parameter_map.end()) {
+ *error = StringPrintf("User specified ordering contains a pointer "
+ "to a double that is not a parameter block in "
+ "the problem. The invalid double is in group: %d",
+ group_it->first);
+ return false;
+ }
+ parameter_blocks->push_back(parameter_block_it->second);
+ }
+ }
+ return true;
}
+
TripletSparseMatrix* SolverImpl::CreateJacobianBlockSparsityTranspose(
const Program* program) {
@@ -1468,34 +1435,185 @@
return tsm;
}
-void SolverImpl::ReorderProgramForSparseNormalCholesky(Program* program) {
-#ifndef CERES_NO_SUITESPARSE
- // Set the offsets and index for CreateJacobianSparsityTranspose.
- program->SetParameterOffsetsAndIndex();
+bool SolverImpl::ReorderProgramForSchurTypeLinearSolver(
+ const LinearSolverType linear_solver_type,
+ const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+ const ProblemImpl::ParameterMap& parameter_map,
+ ParameterBlockOrdering* parameter_block_ordering,
+ Program* program,
+ string* error) {
+ if (parameter_block_ordering->NumGroups() == 1) {
+ // If the user supplied an parameter_block_ordering with just one
+ // group, it is equivalent to the user supplying NULL as an
+ // parameter_block_ordering. Ceres is completely free to choose the
+ // parameter block ordering as it sees fit. For Schur type solvers,
+ // this means that the user wishes for Ceres to identify the
+ // e_blocks, which we do by computing a maximal independent set.
+ vector<ParameterBlock*> schur_ordering;
+ const int num_eliminate_blocks = ComputeSchurOrdering(*program,
+ &schur_ordering);
- // Compute a block sparse presentation of J'.
- scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
- SolverImpl::CreateJacobianBlockSparsityTranspose(program));
+ CHECK_EQ(schur_ordering.size(), program->NumParameterBlocks())
+ << "Congratulations, you found a Ceres bug! Please report this error "
+ << "to the developers.";
- // Order rows using AMD.
- SuiteSparse ss;
- cholmod_sparse* block_jacobian_transpose =
- ss.CreateSparseMatrix(tsm_block_jacobian_transpose.get());
+ // Update the parameter_block_ordering object.
+ for (int i = 0; i < schur_ordering.size(); ++i) {
+ double* parameter_block = schur_ordering[i]->mutable_user_state();
+ const int group_id = (i < num_eliminate_blocks) ? 0 : 1;
+ parameter_block_ordering->AddElementToGroup(parameter_block, group_id);
+ }
- vector<int> ordering(program->NumParameterBlocks(), -1);
- ss.ApproximateMinimumDegreeOrdering(block_jacobian_transpose, &ordering[0]);
- ss.Free(block_jacobian_transpose);
+ // We could call ApplyUserOrdering but this is cheaper and
+ // simpler.
+ swap(*program->mutable_parameter_blocks(), schur_ordering);
+ } else {
+ // The user provided an ordering with more than one elimination
+ // group. Trust the user and apply the ordering.
+ if (!ApplyUserOrdering(parameter_map,
+ parameter_block_ordering,
+ program,
+ error)) {
+ return false;
+ }
+ }
- // Apply ordering.
- vector<ParameterBlock*>& parameter_blocks =
- *(program->mutable_parameter_blocks());
- const vector<ParameterBlock*> parameter_blocks_copy(parameter_blocks);
- for (int i = 0; i < program->NumParameterBlocks(); ++i) {
- parameter_blocks[i] = parameter_blocks_copy[ordering[i]];
+ // Pre-order the columns corresponding to the schur complement if
+ // possible.
+#if !defined(CERES_NO_SUITESPARSE) && !defined(CERES_NO_CAMD)
+ if (linear_solver_type == SPARSE_SCHUR &&
+ sparse_linear_algebra_library_type == SUITE_SPARSE) {
+ vector<int> constraints;
+ vector<ParameterBlock*>& parameter_blocks =
+ *(program->mutable_parameter_blocks());
+
+ for (int i = 0; i < parameter_blocks.size(); ++i) {
+ constraints.push_back(
+ parameter_block_ordering->GroupId(
+ parameter_blocks[i]->mutable_user_state()));
+ }
+
+ // Set the offsets and index for CreateJacobianSparsityTranspose.
+ program->SetParameterOffsetsAndIndex();
+ // Compute a block sparse presentation of J'.
+ scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
+ SolverImpl::CreateJacobianBlockSparsityTranspose(program));
+
+ SuiteSparse ss;
+ cholmod_sparse* block_jacobian_transpose =
+ ss.CreateSparseMatrix(tsm_block_jacobian_transpose.get());
+
+ vector<int> ordering(parameter_blocks.size(), 0);
+ ss.ConstrainedApproximateMinimumDegreeOrdering(block_jacobian_transpose,
+ &constraints[0],
+ &ordering[0]);
+ ss.Free(block_jacobian_transpose);
+
+ const vector<ParameterBlock*> parameter_blocks_copy(parameter_blocks);
+ for (int i = 0; i < program->NumParameterBlocks(); ++i) {
+ parameter_blocks[i] = parameter_blocks_copy[ordering[i]];
+ }
}
#endif
program->SetParameterOffsetsAndIndex();
+ // Schur type solvers also require that their residual blocks be
+ // lexicographically ordered.
+ const int num_eliminate_blocks =
+ parameter_block_ordering->group_to_elements().begin()->second.size();
+ return LexicographicallyOrderResidualBlocks(num_eliminate_blocks,
+ program,
+ error);
+}
+
+bool SolverImpl::ReorderProgramForSparseNormalCholesky(
+ const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+ const ParameterBlockOrdering* parameter_block_ordering,
+ Program* program,
+ string* error) {
+ // Set the offsets and index for CreateJacobianSparsityTranspose.
+ program->SetParameterOffsetsAndIndex();
+ // Compute a block sparse presentation of J'.
+ scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
+ SolverImpl::CreateJacobianBlockSparsityTranspose(program));
+
+ vector<int> ordering(program->NumParameterBlocks(), 0);
+ vector<ParameterBlock*>& parameter_blocks =
+ *(program->mutable_parameter_blocks());
+
+ if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
+#ifdef CERES_NO_SUITESPARSE
+ *error = "Can't use SPARSE_NORMAL_CHOLESKY with SUITE_SPARSE because "
+ "SuiteSparse was not enabled when Ceres was built.";
+ return false;
+#else
+ SuiteSparse ss;
+ cholmod_sparse* block_jacobian_transpose =
+ ss.CreateSparseMatrix(tsm_block_jacobian_transpose.get());
+
+# ifdef CERES_NO_CAMD
+ // No cholmod_camd, so ignore user's parameter_block_ordering and
+ // use plain old AMD.
+ ss.ApproximateMinimumDegreeOrdering(block_jacobian_transpose, &ordering[0]);
+# else
+ if (parameter_block_ordering->NumGroups() > 1) {
+ // If the user specified more than one elimination groups use them
+ // to constrain the ordering.
+ vector<int> constraints;
+ for (int i = 0; i < parameter_blocks.size(); ++i) {
+ constraints.push_back(
+ parameter_block_ordering->GroupId(
+ parameter_blocks[i]->mutable_user_state()));
+ }
+ ss.ConstrainedApproximateMinimumDegreeOrdering(
+ block_jacobian_transpose,
+ &constraints[0],
+ &ordering[0]);
+ } else {
+ ss.ApproximateMinimumDegreeOrdering(block_jacobian_transpose,
+ &ordering[0]);
+ }
+# endif // CERES_NO_CAMD
+
+ ss.Free(block_jacobian_transpose);
+#endif // CERES_NO_SUITESPARSE
+
+ } else if (sparse_linear_algebra_library_type == CX_SPARSE) {
+#ifndef CERES_NO_CXSPARSE
+
+ // CXSparse works with J'J instead of J'. So compute the block
+ // sparsity for J'J and compute an approximate minimum degree
+ // ordering.
+ CXSparse cxsparse;
+ cs_di* block_jacobian_transpose;
+ block_jacobian_transpose =
+ cxsparse.CreateSparseMatrix(tsm_block_jacobian_transpose.get());
+ cs_di* block_jacobian = cxsparse.TransposeMatrix(block_jacobian_transpose);
+ cs_di* block_hessian =
+ cxsparse.MatrixMatrixMultiply(block_jacobian_transpose, block_jacobian);
+ cxsparse.Free(block_jacobian);
+ cxsparse.Free(block_jacobian_transpose);
+
+ cxsparse.ApproximateMinimumDegreeOrdering(block_hessian, &ordering[0]);
+ cxsparse.Free(block_hessian);
+#else // CERES_NO_CXSPARSE
+ *error = "Can't use SPARSE_NORMAL_CHOLESKY with CX_SPARSE because "
+ "CXSparse was not enabled when Ceres was built.";
+ return false;
+#endif // CERES_NO_CXSPARSE
+ } else {
+ *error = "Unknown sparse linear algebra library.";
+ return false;
+ }
+
+ // Apply ordering.
+ const vector<ParameterBlock*> parameter_blocks_copy(parameter_blocks);
+ for (int i = 0; i < program->NumParameterBlocks(); ++i) {
+ parameter_blocks[i] = parameter_blocks_copy[ordering[i]];
+ }
+
+ program->SetParameterOffsetsAndIndex();
+ return true;
}
} // namespace internal
diff --git a/internal/ceres/solver_impl.h b/internal/ceres/solver_impl.h
index 60fd2a2..ebfb813 100644
--- a/internal/ceres/solver_impl.h
+++ b/internal/ceres/solver_impl.h
@@ -103,15 +103,6 @@
static LinearSolver* CreateLinearSolver(Solver::Options* options,
string* error);
- // Reorder the parameter blocks in program using the ordering. A
- // return value of true indicates success and false indicates an
- // error was encountered whose cause is logged to LOG(ERROR).
- static bool ApplyUserOrdering(const ProblemImpl::ParameterMap& parameter_map,
- const ParameterBlockOrdering* ordering,
- Program* program,
- string* error);
-
-
// Reorder the residuals for program, if necessary, so that the
// residuals involving e block (i.e., the first num_eliminate_block
// parameter blocks) occur together. This is a necessary condition
@@ -163,29 +154,6 @@
static void AlternateLinearSolverForSchurTypeLinearSolver(
Solver::Options* options);
- // Schur type solvers require that all parameter blocks eliminated
- // by the Schur eliminator occur before others and the residuals be
- // sorted in lexicographic order of their parameter blocks.
- //
- // If ordering has at least two groups, then apply the ordering,
- // otherwise compute a new ordering using a Maximal Independent Set
- // algorithm and apply it.
- //
- // Upon return, ordering contains the parameter block ordering that
- // was used to order the program.
- static bool ReorderProgramForSchurTypeLinearSolver(
- const ProblemImpl::ParameterMap& parameter_map,
- ParameterBlockOrdering* ordering,
- Program* program,
- string* error);
-
- // CHOLMOD when doing the sparse cholesky factorization of the
- // Jacobian matrix, reorders its columns to reduce the
- // fill-in. Compute this permutation and re-order the parameter
- // blocks.
- //
- static void ReorderProgramForSparseNormalCholesky(Program* program);
-
// Create a TripletSparseMatrix which contains the zero-one
// structure corresponding to the block sparsity of the transpose of
// the Jacobian matrix.
@@ -193,6 +161,53 @@
// Caller owns the result.
static TripletSparseMatrix* CreateJacobianBlockSparsityTranspose(
const Program* program);
+
+ // Reorder the parameter blocks in program using the ordering
+ static bool ApplyUserOrdering(
+ const ProblemImpl::ParameterMap& parameter_map,
+ const ParameterBlockOrdering* parameter_block_ordering,
+ Program* program,
+ string* error);
+
+ // Sparse cholesky factorization routines when doing the sparse
+ // cholesky factorization of the Jacobian matrix, reorders its
+ // columns to reduce the fill-in. Compute this permutation and
+ // re-order the parameter blocks.
+ //
+ // If the parameter_block_ordering contains more than one
+ // elimination group and support for constrained fill-reducing
+ // ordering is available in the sparse linear algebra library
+ // (SuiteSparse version >= 4.2.0) then the fill reducing
+ // ordering will take it into account, otherwise it will be ignored.
+ static bool ReorderProgramForSparseNormalCholesky(
+ const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+ const ParameterBlockOrdering* parameter_block_ordering,
+ Program* program,
+ string* error);
+
+ // Schur type solvers require that all parameter blocks eliminated
+ // by the Schur eliminator occur before others and the residuals be
+ // sorted in lexicographic order of their parameter blocks.
+ //
+ // If the parameter_block_ordering only contains one elimination
+ // group then a maximal independent set is computed and used as the
+ // first elimination group, otherwise the user's ordering is used.
+ //
+ // If the linear solver type is SPARSE_SCHUR and support for
+ // constrained fill-reducing ordering is available in the sparse
+ // linear algebra library (SuiteSparse version >= 4.2.0) then
+ // columns of the schur complement matrix are ordered to reduce the
+ // fill-in the Cholesky factorization.
+ //
+ // Upon return, ordering contains the parameter block ordering that
+ // was used to order the program.
+ static bool ReorderProgramForSchurTypeLinearSolver(
+ const LinearSolverType linear_solver_type,
+ const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+ const ProblemImpl::ParameterMap& parameter_map,
+ ParameterBlockOrdering* parameter_block_ordering,
+ Program* program,
+ string* error);
};
} // namespace internal
diff --git a/internal/ceres/solver_impl_test.cc b/internal/ceres/solver_impl_test.cc
index a752eff..e99a3de 100644
--- a/internal/ceres/solver_impl_test.cc
+++ b/internal/ceres/solver_impl_test.cc
@@ -846,5 +846,83 @@
EXPECT_EQ(options.linear_solver_type, CGNR);
EXPECT_EQ(options.preconditioner_type, JACOBI);
}
+
+TEST(SolverImpl, CreateJacobianBlockSparsityTranspose) {
+ ProblemImpl problem;
+ double x[2];
+ double y[3];
+ double z;
+
+ problem.AddParameterBlock(x, 2);
+ problem.AddParameterBlock(y, 3);
+ problem.AddParameterBlock(&z, 1);
+
+ problem.AddResidualBlock(new MockCostFunctionBase<2, 2, 0, 0>(), NULL, x);
+ problem.AddResidualBlock(new MockCostFunctionBase<3, 1, 2, 0>(), NULL, &z, x);
+ problem.AddResidualBlock(new MockCostFunctionBase<4, 1, 3, 0>(), NULL, &z, y);
+ problem.AddResidualBlock(new MockCostFunctionBase<5, 1, 3, 0>(), NULL, &z, y);
+ problem.AddResidualBlock(new MockCostFunctionBase<1, 2, 1, 0>(), NULL, x, &z);
+ problem.AddResidualBlock(new MockCostFunctionBase<2, 1, 3, 0>(), NULL, &z, y);
+ problem.AddResidualBlock(new MockCostFunctionBase<2, 2, 1, 0>(), NULL, x, &z);
+ problem.AddResidualBlock(new MockCostFunctionBase<1, 3, 0, 0>(), NULL, y);
+
+ TripletSparseMatrix expected_block_sparse_jacobian(3, 8, 14);
+ {
+ int* rows = expected_block_sparse_jacobian.mutable_rows();
+ int* cols = expected_block_sparse_jacobian.mutable_cols();
+ double* values = expected_block_sparse_jacobian.mutable_values();
+ rows[0] = 0;
+ cols[0] = 0;
+
+ rows[1] = 2;
+ cols[1] = 1;
+ rows[2] = 0;
+ cols[2] = 1;
+
+ rows[3] = 2;
+ cols[3] = 2;
+ rows[4] = 1;
+ cols[4] = 2;
+
+ rows[5] = 2;
+ cols[5] = 3;
+ rows[6] = 1;
+ cols[6] = 3;
+
+ rows[7] = 0;
+ cols[7] = 4;
+ rows[8] = 2;
+ cols[8] = 4;
+
+ rows[9] = 2;
+ cols[9] = 5;
+ rows[10] = 1;
+ cols[10] = 5;
+
+ rows[11] = 0;
+ cols[11] = 6;
+ rows[12] = 2;
+ cols[12] = 6;
+
+ rows[13] = 1;
+ cols[13] = 7;
+ fill(values, values + 14, 1.0);
+ expected_block_sparse_jacobian.set_num_nonzeros(14);
+ }
+
+ Program* program = problem.mutable_program();
+ program->SetParameterOffsetsAndIndex();
+
+ scoped_ptr<TripletSparseMatrix> actual_block_sparse_jacobian(
+ SolverImpl::CreateJacobianBlockSparsityTranspose(program));
+
+ Matrix expected_dense_jacobian;
+ expected_block_sparse_jacobian.ToDenseMatrix(&expected_dense_jacobian);
+
+ Matrix actual_dense_jacobian;
+ actual_block_sparse_jacobian->ToDenseMatrix(&actual_dense_jacobian);
+ EXPECT_EQ((expected_dense_jacobian - actual_dense_jacobian).norm(), 0.0);
+}
+
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index bc1f983..9601142 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -133,34 +133,37 @@
// factorized. CHOLMOD/SuiteSparse on the other hand can just work
// off of Jt to compute the Cholesky factorization of the normal
// equations.
- cs_di* A2 = cs_transpose(&At, 1);
- cs_di* AtA = cs_multiply(&At, A2);
+ cs_di* A2 = cxsparse_.TransposeMatrix(&At);
+ cs_di* AtA = cxsparse_.MatrixMatrixMultiply(&At, A2);
cxsparse_.Free(A2);
if (per_solve_options.D != NULL) {
A->DeleteRows(num_cols);
}
-
event_logger.AddEvent("Setup");
// Compute symbolic factorization if not available.
if (cxsparse_factor_ == NULL) {
- cxsparse_factor_ = CHECK_NOTNULL(cxsparse_.AnalyzeCholesky(AtA));
+ if (options_.use_postordering) {
+ cxsparse_factor_ =
+ CHECK_NOTNULL(cxsparse_.BlockAnalyzeCholesky(AtA,
+ A->col_blocks(),
+ A->col_blocks()));
+ } else {
+ cxsparse_factor_ =
+ CHECK_NOTNULL(cxsparse_.AnalyzeCholeskyWithNaturalOrdering(AtA));
+ }
}
-
event_logger.AddEvent("Analysis");
-
// Solve the linear system.
if (cxsparse_.SolveCholesky(AtA, cxsparse_factor_, Atb.data())) {
VectorRef(x, Atb.rows()) = Atb;
summary.termination_type = TOLERANCE;
}
-
event_logger.AddEvent("Solve");
cxsparse_.Free(AtA);
-
event_logger.AddEvent("Teardown");
return summary;
}
@@ -205,11 +208,13 @@
if (factor_ == NULL) {
if (options_.use_postordering) {
- factor_ = ss_.BlockAnalyzeCholesky(&lhs,
- A->col_blocks(),
- A->row_blocks());
+ factor_ =
+ CHECK_NOTNULL(ss_.BlockAnalyzeCholesky(&lhs,
+ A->col_blocks(),
+ A->row_blocks()));
} else {
- factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs);
+ factor_ =
+ CHECK_NOTNULL(ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs));
}
}
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc
index fe2edd3..57d12a1 100644
--- a/internal/ceres/suitesparse.cc
+++ b/internal/ceres/suitesparse.cc
@@ -323,6 +323,21 @@
cholmod_amd(matrix, NULL, 0, ordering, &cc_);
}
+void SuiteSparse::ConstrainedApproximateMinimumDegreeOrdering(
+ cholmod_sparse* matrix,
+ int* constraints,
+ int* ordering) {
+#ifndef CERES_NO_CAMD
+ cholmod_camd(matrix, NULL, 0, constraints, ordering, &cc_);
+#else
+ LOG(FATAL) << "Congratulations you have found a bug in Ceres."
+ << "Ceres Solver was compiled with SuiteSparse "
+ << "version 4.1.0 or less. Calling this function "
+ << "in that case is a bug. Please contact the"
+ << "the Ceres Solver developers".
+#endif
+}
+
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h
index e138623..27182b8 100644
--- a/internal/ceres/suitesparse.h
+++ b/internal/ceres/suitesparse.h
@@ -33,6 +33,7 @@
#ifndef CERES_INTERNAL_SUITESPARSE_H_
#define CERES_INTERNAL_SUITESPARSE_H_
+
#ifndef CERES_NO_SUITESPARSE
#include <cstring>
@@ -43,6 +44,20 @@
#include "cholmod.h"
#include "glog/logging.h"
+// Before SuiteSparse version 4.2.0, cholmod_camd was only enabled
+// if SuiteSparse was compiled with Metis support. This makes
+// calling and linking into cholmod_camd problematic even though it
+// has nothing to do with Metis. This has been fixed reliably in
+// 4.2.0.
+//
+// The fix was actually committed in 4.1.0, but there is
+// some confusion about a silent update to the tar ball, so we are
+// being conservative and choosing the next minor version where
+// things are stable.
+#if (SUITESPARSE_VERSION<4002)
+#define CERES_NO_CAMD
+#endif
+
namespace ceres {
namespace internal {
@@ -189,6 +204,38 @@
// ordering.
void ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, int* ordering);
+
+ // Before SuiteSparse version 4.2.0, cholmod_camd was only enabled
+ // if SuiteSparse was compiled with Metis support. This makes
+ // calling and linking into cholmod_camd problematic even though it
+ // has nothing to do with Metis. This has been fixed reliably in
+ // 4.2.0.
+ //
+ // The fix was actually committed in 4.1.0, but there is
+ // some confusion about a silent update to the tar ball, so we are
+ // being conservative and choosing the next minor version where
+ // things are stable.
+ static bool IsConstrainedApproximateMinimumDegreeOrderingAvailable() {
+ return (SUITESPARSE_VERSION>4001);
+ }
+
+ // Find a fill reducing approximate minimum degree
+ // ordering. constraints is an array which associates with each
+ // column of the matrix an elimination group. i.e., all columns in
+ // group 0 are eliminated first, all columns in group 1 are
+ // eliminated next etc. This function finds a fill reducing ordering
+ // that obeys these constraints.
+ //
+ // Calling ApproximateMinimumDegreeOrdering is equivalent to calling
+ // ConstrainedApproximateMinimumDegreeOrdering with a constraint
+ // array that puts all columns in the same elimination group.
+ //
+ // If CERES_NO_CAMD is defined then calling this function will
+ // result in a crash.
+ void ConstrainedApproximateMinimumDegreeOrdering(cholmod_sparse* matrix,
+ int* constraints,
+ int* ordering);
+
void Free(cholmod_sparse* m) { cholmod_free_sparse(&m, &cc_); }
void Free(cholmod_dense* m) { cholmod_free_dense(&m, &cc_); }
void Free(cholmod_factor* m) { cholmod_free_factor(&m, &cc_); }
diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/unsymmetric_linear_solver_test.cc
index c8a15c0..34b03be 100644
--- a/internal/ceres/unsymmetric_linear_solver_test.cc
+++ b/internal/ceres/unsymmetric_linear_solver_test.cc
@@ -56,12 +56,7 @@
sol_regularized_.reset(problem->x_D.release());
}
- void TestSolver(
- LinearSolverType linear_solver_type,
- SparseLinearAlgebraLibraryType sparse_linear_algebra_library) {
- LinearSolver::Options options;
- options.type = linear_solver_type;
- options.sparse_linear_algebra_library = sparse_linear_algebra_library;
+ void TestSolver(const LinearSolver::Options& options) {
scoped_ptr<LinearSolver> solver(LinearSolver::Create(options));
LinearSolver::PerSolveOptions per_solve_options;
@@ -72,13 +67,22 @@
scoped_ptr<SparseMatrix> transformed_A;
- if (linear_solver_type == DENSE_QR ||
- linear_solver_type == DENSE_NORMAL_CHOLESKY) {
+ if (options.type == DENSE_QR ||
+ options.type == DENSE_NORMAL_CHOLESKY) {
transformed_A.reset(new DenseSparseMatrix(*A_));
- } else if (linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
- transformed_A.reset(new CompressedRowSparseMatrix(*A_));
+ } else if (options.type == SPARSE_NORMAL_CHOLESKY) {
+ CompressedRowSparseMatrix* crsm = new CompressedRowSparseMatrix(*A_);
+ // Add row/column blocks structure.
+ for (int i = 0; i < A_->num_rows(); ++i) {
+ crsm->mutable_row_blocks()->push_back(1);
+ }
+
+ for (int i = 0; i < A_->num_cols(); ++i) {
+ crsm->mutable_col_blocks()->push_back(1);
+ }
+ transformed_A.reset(crsm);
} else {
- LOG(FATAL) << "Unknown linear solver : " << linear_solver_type;
+ LOG(FATAL) << "Unknown linear solver : " << options.type;
}
// Unregularized
unregularized_solve_summary =
@@ -115,22 +119,50 @@
};
TEST_F(UnsymmetricLinearSolverTest, DenseQR) {
- TestSolver(DENSE_QR, SUITE_SPARSE);
+ LinearSolver::Options options;
+ options.type = DENSE_QR;
+ TestSolver(options);
}
TEST_F(UnsymmetricLinearSolverTest, DenseNormalCholesky) {
- TestSolver(DENSE_NORMAL_CHOLESKY, SUITE_SPARSE);
+ LinearSolver::Options options;
+ options.type = DENSE_NORMAL_CHOLESKY;
+ TestSolver(options);
}
#ifndef CERES_NO_SUITESPARSE
-TEST_F(UnsymmetricLinearSolverTest, SparseNormalCholeskyUsingSuiteSparse) {
- TestSolver(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE);
+TEST_F(UnsymmetricLinearSolverTest, SparseNormalCholeskyUsingSuiteSparsePreOrdering) {
+ LinearSolver::Options options;
+ options.sparse_linear_algebra_library = SUITE_SPARSE;
+ options.type = SPARSE_NORMAL_CHOLESKY;
+ options.use_postordering = false;
+ TestSolver(options);
+}
+
+TEST_F(UnsymmetricLinearSolverTest, SparseNormalCholeskyUsingSuiteSparsePostOrdering) {
+ LinearSolver::Options options;
+ options.sparse_linear_algebra_library = SUITE_SPARSE;
+ options.type = SPARSE_NORMAL_CHOLESKY;
+ options.use_postordering = true;
+ TestSolver(options);
}
#endif
#ifndef CERES_NO_CXSPARSE
-TEST_F(UnsymmetricLinearSolverTest, SparseNormalCholeskyUsingCXSparse) {
- TestSolver(SPARSE_NORMAL_CHOLESKY, CX_SPARSE);
+TEST_F(UnsymmetricLinearSolverTest, SparseNormalCholeskyUsingCXSparsePreOrdering) {
+ LinearSolver::Options options;
+ options.sparse_linear_algebra_library = CX_SPARSE;
+ options.type = SPARSE_NORMAL_CHOLESKY;
+ options.use_postordering = false;
+ TestSolver(options);
+}
+
+TEST_F(UnsymmetricLinearSolverTest, SparseNormalCholeskyUsingCXSparsePostOrdering) {
+ LinearSolver::Options options;
+ options.sparse_linear_algebra_library = CX_SPARSE;
+ options.type = SPARSE_NORMAL_CHOLESKY;
+ options.use_postordering = true;
+ TestSolver(options);
}
#endif
