Refactor SparseNormalCholeskySolver
Now that there is a single piece of code doing the outer product
computation for all three sparse linear algebra backends, move
this code one level up the call stack and there by make the actual
per-library solver code shorter and simpler.
Also fix a minor omission in the outer product computation code
where row/column blocks were not being copied over to the
outer product matrix.
Change-Id: I22a7967bdc659385b741901afefa7af312e676e5
diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc
index c5f64dd..9c438f6 100644
--- a/internal/ceres/compressed_row_sparse_matrix.cc
+++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -561,6 +561,9 @@
new CompressedRowSparseMatrix(num_cols, num_cols, num_nonzeros);
matrix->set_storage_type(storage_type);
+ *(matrix->mutable_row_blocks()) = blocks;
+ *(matrix->mutable_col_blocks()) = blocks;
+
// Compute block offsets for outer product matrix, which is used in
// ComputeOuterProduct.
vector<int>* block_offsets = matrix->mutable_block_offsets();
diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc
index ddaa8cc..66e2390 100644
--- a/internal/ceres/compressed_row_sparse_matrix_test.cc
+++ b/internal/ceres/compressed_row_sparse_matrix_test.cc
@@ -450,6 +450,9 @@
CompressedRowSparseMatrix::ComputeOuterProduct(
*random_matrix, program, outer_product.get());
+ EXPECT_EQ(outer_product->row_blocks(), random_matrix->col_blocks());
+ EXPECT_EQ(outer_product->col_blocks(), random_matrix->col_blocks());
+
Matrix actual_outer_product =
Eigen::MappedSparseMatrix<double, Eigen::ColMajor>(
outer_product->num_rows(),
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index 0a35658..9c725ff 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -35,6 +35,7 @@
#include <ctime>
#include <sstream>
+#include "Eigen/SparseCore"
#include "ceres/compressed_row_sparse_matrix.h"
#include "ceres/cxsparse.h"
#include "ceres/internal/eigen.h"
@@ -44,7 +45,6 @@
#include "ceres/triplet_sparse_matrix.h"
#include "ceres/types.h"
#include "ceres/wall_time.h"
-#include "Eigen/SparseCore"
#ifdef CERES_USE_EIGEN_SPARSE
#include "Eigen/SparseCholesky"
@@ -52,6 +52,34 @@
namespace ceres {
namespace internal {
+
+// Different sparse linear algebra libraries prefer different storage
+// orders for the input matrix. This trait class helps choose the
+// ordering based on the sparse linear algebra backend being used.
+//
+// The storage order is lower-triangular by default. It is only
+// SuiteSparse which prefers an upper triangular matrix. Saves a whole
+// matrix copy in the process.
+//
+// Note that this is the storage order for a compressed row sparse
+// matrix. All the sparse linear algebra libraries take compressed
+// column sparse matrices as input. We map these matrices to into
+// compressed column sparse matrices before calling them and in the
+// process, transpose them.
+//
+// TODO(sameeragarwal): This does not account for post ordering, where
+// the optimal storage order maybe different. Either get rid of post
+// ordering support entirely, or investigate making this trait class
+// richer.
+
+CompressedRowSparseMatrix::StorageType StorageTypeForSparseLinearAlgebraLibrary(
+ SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type) {
+ if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
+ return CompressedRowSparseMatrix::UPPER_TRIANGULAR;
+ }
+ return CompressedRowSparseMatrix::LOWER_TRIANGULAR;
+}
+
namespace {
#ifdef CERES_USE_EIGEN_SPARSE
@@ -59,12 +87,11 @@
// the same code independent of whether a AMD or a Natural ordering is
// used.
template <typename SimplicialCholeskySolver, typename SparseMatrixType>
-LinearSolver::Summary SimplicialLDLTSolve(
- const SparseMatrixType& lhs,
- const bool do_symbolic_analysis,
- SimplicialCholeskySolver* solver,
- double* rhs_and_solution,
- EventLogger* event_logger) {
+LinearSolver::Summary SimplicialLDLTSolve(const SparseMatrixType& lhs,
+ const bool do_symbolic_analysis,
+ SimplicialCholeskySolver* solver,
+ double* rhs_and_solution,
+ EventLogger* event_logger) {
LinearSolver::Summary summary;
summary.num_iterations = 1;
summary.termination_type = LINEAR_SOLVER_SUCCESS;
@@ -75,14 +102,12 @@
if (VLOG_IS_ON(2)) {
std::stringstream ss;
solver->dumpMemory(ss);
- VLOG(2) << "Symbolic Analysis\n"
- << ss.str();
+ VLOG(2) << "Symbolic Analysis\n" << ss.str();
}
event_logger->AddEvent("Analyze");
if (solver->info() != Eigen::Success) {
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
- summary.message =
- "Eigen failure. Unable to find symbolic factorization.";
+ summary.message = "Eigen failure. Unable to find symbolic factorization.";
return summary;
}
}
@@ -114,10 +139,7 @@
SparseNormalCholeskySolver::SparseNormalCholeskySolver(
const LinearSolver::Options& options)
- : factor_(NULL),
- cxsparse_factor_(NULL),
- options_(options) {
-}
+ : factor_(NULL), cxsparse_factor_(NULL), options_(options) {}
void SparseNormalCholeskySolver::FreeFactorization() {
if (factor_ != NULL) {
@@ -139,8 +161,7 @@
CompressedRowSparseMatrix* A,
const double* b,
const LinearSolver::PerSolveOptions& per_solve_options,
- double * x) {
-
+ double* x) {
const int num_cols = A->num_cols();
VectorRef(x, num_cols).setZero();
A->LeftMultiply(b, x);
@@ -151,24 +172,36 @@
scoped_ptr<CompressedRowSparseMatrix> regularizer;
if (A->col_blocks().size() > 0) {
regularizer.reset(CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(
- per_solve_options.D, A->col_blocks()));
+ per_solve_options.D, A->col_blocks()));
} else {
- regularizer.reset(new CompressedRowSparseMatrix(
- per_solve_options.D, num_cols));
+ regularizer.reset(
+ new CompressedRowSparseMatrix(per_solve_options.D, num_cols));
}
A->AppendRows(*regularizer);
}
+ if (outer_product_.get() == NULL) {
+ outer_product_.reset(
+ CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
+ *A,
+ StorageTypeForSparseLinearAlgebraLibrary(
+ options_.sparse_linear_algebra_library_type),
+ &pattern_));
+ }
+
+ CompressedRowSparseMatrix::ComputeOuterProduct(
+ *A, pattern_, outer_product_.get());
+
LinearSolver::Summary summary;
switch (options_.sparse_linear_algebra_library_type) {
case SUITE_SPARSE:
- summary = SolveImplUsingSuiteSparse(A, x);
+ summary = SolveImplUsingSuiteSparse(x);
break;
case CX_SPARSE:
- summary = SolveImplUsingCXSparse(A, x);
+ summary = SolveImplUsingCXSparse(x);
break;
case EIGEN_SPARSE:
- summary = SolveImplUsingEigen(A, x);
+ summary = SolveImplUsingEigen(x);
break;
default:
LOG(FATAL) << "Unknown sparse linear algebra library : "
@@ -183,8 +216,7 @@
}
LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingEigen(
- CompressedRowSparseMatrix* A,
- double * rhs_and_solution) {
+ double* rhs_and_solution) {
#ifndef CERES_USE_EIGEN_SPARSE
LinearSolver::Summary summary;
@@ -200,23 +232,8 @@
#else
EventLogger event_logger("SparseNormalCholeskySolver::Eigen::Solve");
- // Compute the normal equations. J'J delta = J'f and solve them
- // using a sparse Cholesky factorization.
- // Compute outer product as a compressed row lower triangular
- // matrix, because after mapping to a column major matrix, this will
- // become a compressed column upper triangular matrix. Which is the
- // default that Eigen uses.
- if (outer_product_.get() == NULL) {
- outer_product_.reset(
- CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
- *A, CompressedRowSparseMatrix::LOWER_TRIANGULAR, &pattern_));
- }
-
- CompressedRowSparseMatrix::ComputeOuterProduct(
- *A, pattern_, outer_product_.get());
-
- // Map to an upper triangular column major matrix.
+ // Map outer_product_ to an upper triangular column major matrix.
//
// outer_product_ is a compressed row sparse matrix and in lower
// triangular form, when mapped to a compressed column sparse
@@ -234,8 +251,7 @@
// If using post ordering or an old version of Eigen, we cannot
// depend on a preordered jacobian, so we work with a SimplicialLDLT
// decomposition with AMD ordering.
- if (options_.use_postordering ||
- !EIGEN_VERSION_AT_LEAST(3, 2, 2)) {
+ if (options_.use_postordering || !EIGEN_VERSION_AT_LEAST(3, 2, 2)) {
if (amd_ldlt_.get() == NULL) {
amd_ldlt_.reset(new SimplicialLDLTWithAMDOrdering);
do_symbolic_analysis = true;
@@ -248,7 +264,7 @@
&event_logger);
}
-#if EIGEN_VERSION_AT_LEAST(3,2,2)
+#if EIGEN_VERSION_AT_LEAST(3, 2, 2)
// The common case
if (natural_ldlt_.get() == NULL) {
natural_ldlt_.reset(new SimplicialLDLTWithNaturalOrdering);
@@ -266,8 +282,7 @@
}
LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse(
- CompressedRowSparseMatrix* A,
- double * rhs_and_solution) {
+ double* rhs_and_solution) {
#ifdef CERES_NO_CXSPARSE
LinearSolver::Summary summary;
@@ -288,20 +303,11 @@
summary.termination_type = LINEAR_SOLVER_SUCCESS;
summary.message = "Success.";
-
- // Compute outer product as a compressed row lower triangular
- // matrix, which would be mapped to a compressed column upper
- // triangular matrix, which is the representation used by CXSparse's
- // sparse Cholesky factorization.
- if (outer_product_.get() == NULL) {
- outer_product_.reset(
- CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
- *A, CompressedRowSparseMatrix::LOWER_TRIANGULAR, &pattern_));
- }
-
- CompressedRowSparseMatrix::ComputeOuterProduct(
- *A, pattern_, outer_product_.get());
-
+ // Map outer_product_ to an upper triangular column major matrix.
+ //
+ // outer_product_ is a compressed row sparse matrix and in lower
+ // triangular form, when mapped to a compressed column sparse
+ // matrix, it becomes an upper triangular matrix.
cs_di lhs = cxsparse_.CreateSparseMatrixTransposeView(outer_product_.get());
event_logger.AddEvent("Setup");
@@ -309,9 +315,8 @@
// Compute symbolic factorization if not available.
if (cxsparse_factor_ == NULL) {
if (options_.use_postordering) {
- cxsparse_factor_ = cxsparse_.BlockAnalyzeCholesky(&lhs,
- A->col_blocks(),
- A->col_blocks());
+ cxsparse_factor_ = cxsparse_.BlockAnalyzeCholesky(
+ &lhs, outer_product_->col_blocks(), outer_product_->col_blocks());
} else {
cxsparse_factor_ = cxsparse_.AnalyzeCholeskyWithNaturalOrdering(&lhs);
}
@@ -322,9 +327,8 @@
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
summary.message =
"CXSparse failure. Unable to find symbolic factorization.";
- } else if (!cxsparse_.SolveCholesky(&lhs,
- cxsparse_factor_,
- rhs_and_solution)) {
+ } else if (!cxsparse_.SolveCholesky(
+ &lhs, cxsparse_factor_, rhs_and_solution)) {
summary.termination_type = LINEAR_SOLVER_FAILURE;
summary.message = "CXSparse::SolveCholesky failed.";
}
@@ -335,8 +339,7 @@
}
LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse(
- CompressedRowSparseMatrix* A,
- double * rhs_and_solution) {
+ double* rhs_and_solution) {
#ifdef CERES_NO_SUITESPARSE
LinearSolver::Summary summary;
@@ -356,33 +359,25 @@
summary.num_iterations = 1;
summary.message = "Success.";
- // Compute outer product to compressed row upper triangular matrix,
- // this will be mapped to a compressed-column lower triangular
- // matrix, which is the fastest option for our default natural
- // ordering (see comment in cholmod_factorize.c:205 in SuiteSparse).
- if (outer_product_.get() == NULL) {
- outer_product_.reset(
- CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
- *A, CompressedRowSparseMatrix::UPPER_TRIANGULAR, &pattern_));
- }
-
- CompressedRowSparseMatrix::ComputeOuterProduct(
- *A, pattern_, outer_product_.get());
-
- const int num_cols = A->num_cols();
+ // Map outer_product_ to an lower triangular column major matrix.
+ //
+ // outer_product_ is a compressed row sparse matrix and in upper
+ // triangular form, when mapped to a compressed column sparse
+ // matrix, it becomes an lower triangular matrix.
+ const int num_cols = outer_product_->num_cols();
cholmod_sparse lhs =
ss_.CreateSparseMatrixTransposeView(outer_product_.get());
event_logger.AddEvent("Setup");
if (factor_ == NULL) {
if (options_.use_postordering) {
- factor_ = ss_.BlockAnalyzeCholesky(&lhs,
- A->col_blocks(),
- A->col_blocks(),
- &summary.message);
+ factor_ = ss_.BlockAnalyzeCholesky(
+ &lhs,
+ outer_product_->col_blocks(),
+ outer_product_->col_blocks(),
+ &summary.message);
} else {
- factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs,
- &summary.message);
+ factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs, &summary.message);
}
}
@@ -400,9 +395,8 @@
return summary;
}
- cholmod_dense* rhs = ss_.CreateDenseVector(rhs_and_solution,
- num_cols,
- num_cols);
+ cholmod_dense* rhs =
+ ss_.CreateDenseVector(rhs_and_solution, num_cols, num_cols);
cholmod_dense* solution = ss_.Solve(factor_, rhs, &summary.message);
event_logger.AddEvent("Solve");
@@ -421,5 +415,5 @@
#endif
}
-} // namespace internal
-} // namespace ceres
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/sparse_normal_cholesky_solver.h b/internal/ceres/sparse_normal_cholesky_solver.h
index 2a93bc5..ef9d7be 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.h
+++ b/internal/ceres/sparse_normal_cholesky_solver.h
@@ -67,18 +67,9 @@
const LinearSolver::PerSolveOptions& options,
double* x);
- LinearSolver::Summary SolveImplUsingSuiteSparse(
- CompressedRowSparseMatrix* A,
- double* rhs_and_solution);
-
-
- LinearSolver::Summary SolveImplUsingCXSparse(
- CompressedRowSparseMatrix* A,
- double* rhs_and_solution);
-
- LinearSolver::Summary SolveImplUsingEigen(
- CompressedRowSparseMatrix* A,
- double* rhs_and_solution);
+ LinearSolver::Summary SolveImplUsingSuiteSparse(double* rhs_and_solution);
+ LinearSolver::Summary SolveImplUsingCXSparse(double* rhs_and_solution);
+ LinearSolver::Summary SolveImplUsingEigen(double* rhs_and_solution);
void FreeFactorization();
