Let EIGEN_SPARSE + SPARSE_NORMAL_CHOLESKY use block AMD.
Modify SparseNormalCholeskySolver to use a pre-ordered Jacobian
matrix.
Change-Id: Ib4d725d7a2d7bb94ea76dbb3a9b172784dbc8ea0
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a262c5b..d0f7c77 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -251,6 +251,15 @@
MESSAGE(" this option will result in an LGPL licensed version of ")
MESSAGE(" Ceres Solver as the Simplicial Cholesky factorization in Eigen")
MESSAGE(" is licensed under the LGPL. ")
+
+ IF (EIGEN_VERSION VERSION_LESS 3.2.2)
+ MESSAGE(" WARNING:")
+ MESSAGE("")
+ MESSAGE(" Your version of Eigen is older than version 3.2.2.")
+ MESSAGE(" The performance of SPARSE_NORMAL_CHOLESKY and SPARSE_SCHUR")
+ MESSAGE(" linear solvers will suffer. ")
+ ENDIF (EIGEN_VERSION VERSION_LESS 3.2.2)
+
ELSE (EIGENSPARSE)
MESSAGE(" Disabling the use of Eigen as a sparse linear algebra library.")
MESSAGE(" This does not affect the covariance estimation algorithm ")
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index 0940815..e7da6a8 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -28,9 +28,6 @@
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
-// This include must come before any #ifndef check on Ceres compile options.
-#include "ceres/internal/port.h"
-
#include "ceres/sparse_normal_cholesky_solver.h"
#include <algorithm>
@@ -48,15 +45,71 @@
#include "ceres/wall_time.h"
#include "Eigen/SparseCore"
+#ifdef CERES_USE_EIGEN_SPARSE
+#include "Eigen/SparseCholesky"
+#endif
namespace ceres {
namespace internal {
+namespace {
+
+#ifdef CERES_USE_EIGEN_SPARSE
+// A templated factorized and solve function, which allows us to use
+// the same code independent of whether a AMD or a Natural ordering is
+// used.
+template <typename SimplicialCholeskySolver>
+LinearSolver::Summary SimplicialLDLTSolve(
+ Eigen::MappedSparseMatrix<double, Eigen::ColMajor>& 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;
+ summary.message = "Success.";
+
+ if (do_symbolic_analysis) {
+ solver->analyzePattern(lhs);
+ 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.";
+ return summary;
+ }
+ }
+
+ solver->factorize(lhs);
+ event_logger->AddEvent("Factorize");
+ if (solver->info() != Eigen::Success) {
+ summary.termination_type = LINEAR_SOLVER_FAILURE;
+ summary.message = "Eigen failure. Unable to find numeric factorization.";
+ return summary;
+ }
+
+ const Vector rhs = VectorRef(rhs_and_solution, lhs.cols());
+
+ VectorRef(rhs_and_solution, lhs.cols()) = solver->solve(rhs);
+ event_logger->AddEvent("Solve");
+ if (solver->info() != Eigen::Success) {
+ summary.termination_type = LINEAR_SOLVER_FAILURE;
+ summary.message = "Eigen failure. Unable to do triangular solve.";
+ return summary;
+ }
+
+ return summary;
+}
+
+#endif // CERES_USE_EIGEN_SPARSE
+
+} // namespace
SparseNormalCholeskySolver::SparseNormalCholeskySolver(
const LinearSolver::Options& options)
: factor_(NULL),
cxsparse_factor_(NULL),
- options_(options){
+ options_(options) {
}
void SparseNormalCholeskySolver::FreeFactorization() {
@@ -141,12 +194,6 @@
#else
EventLogger event_logger("SparseNormalCholeskySolver::Eigen::Solve");
-
- LinearSolver::Summary summary;
- summary.num_iterations = 1;
- summary.termination_type = LINEAR_SOLVER_SUCCESS;
- summary.message = "Success.";
-
// Compute the normal equations. J'J delta = J'f and solve them
// using a sparse Cholesky factorization. Notice that when compared
// to SuiteSparse we have to explicitly compute the normal equations
@@ -178,42 +225,45 @@
outer_product_->mutable_cols(),
outer_product_->mutable_values());
- const Vector b = VectorRef(rhs_and_solution, outer_product_->num_rows());
- if (simplicial_ldlt_.get() == NULL || options_.dynamic_sparsity) {
- simplicial_ldlt_.reset(new SimplicialLDLT);
- // This is a crappy way to be doing this. But right now Eigen does
- // not expose a way to do symbolic analysis with a given
- // permutation pattern, so we cannot use a block analysis of the
- // Jacobian.
- simplicial_ldlt_->analyzePattern(AtA);
- if (simplicial_ldlt_->info() != Eigen::Success) {
- summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
- summary.message =
- "Eigen failure. Unable to find symbolic factorization.";
- return summary;
+ // Dynamic sparsity means that we cannot depend on a static analysis
+ // of sparsity structure of the jacobian, so we compute a new
+ // symbolic factorization every time.
+ if (options_.dynamic_sparsity) {
+ amd_ldlt_.reset(NULL);
+ }
+
+ bool do_symbolic_analysis = false;
+
+ // If using post ordering or dynamic sparsity, 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 ||
+ options_.dynamic_sparsity ||
+ !EIGEN_VERSION_AT_LEAST(3, 2, 2)) {
+ if (amd_ldlt_.get() == NULL) {
+ amd_ldlt_.reset(new SimplicialLDLTWithAMDOrdering);
+ do_symbolic_analysis = true;
}
- }
- event_logger.AddEvent("Analysis");
- simplicial_ldlt_->factorize(AtA);
- if(simplicial_ldlt_->info() != Eigen::Success) {
- summary.termination_type = LINEAR_SOLVER_FAILURE;
- summary.message =
- "Eigen failure. Unable to find numeric factorization.";
- return summary;
+ return SimplicialLDLTSolve(AtA,
+ do_symbolic_analysis,
+ amd_ldlt_.get(),
+ rhs_and_solution,
+ &event_logger);
}
- VectorRef(rhs_and_solution, outer_product_->num_rows()) =
- simplicial_ldlt_->solve(b);
- if(simplicial_ldlt_->info() != Eigen::Success) {
- summary.termination_type = LINEAR_SOLVER_FAILURE;
- summary.message =
- "Eigen failure. Unable to do triangular solve.";
- return summary;
+ // The common case
+ if (natural_ldlt_.get() == NULL) {
+ natural_ldlt_.reset(new SimplicialLDLTWithNaturalOrdering);
+ do_symbolic_analysis = true;
}
- event_logger.AddEvent("Solve");
- return summary;
+ return SimplicialLDLTSolve(AtA,
+ do_symbolic_analysis,
+ natural_ldlt_.get(),
+ rhs_and_solution,
+ &event_logger);
+
#endif // EIGEN_USE_EIGEN_SPARSE
}
@@ -291,7 +341,9 @@
summary.termination_type = LINEAR_SOLVER_FATAL_ERROR;
summary.message =
"CXSparse failure. Unable to find symbolic factorization.";
- } else if (!cxsparse_.SolveCholesky(AtA, cxsparse_factor_, rhs_and_solution)) {
+ } else if (!cxsparse_.SolveCholesky(AtA,
+ cxsparse_factor_,
+ rhs_and_solution)) {
summary.termination_type = LINEAR_SOLVER_FAILURE;
summary.message = "CXSparse::SolveCholesky failed.";
}
@@ -341,7 +393,8 @@
if (options_.dynamic_sparsity) {
factor_ = ss_.AnalyzeCholesky(&lhs, &summary.message);
} else {
- factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs, &summary.message);
+ factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs,
+ &summary.message);
}
}
}
@@ -359,7 +412,9 @@
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");
diff --git a/internal/ceres/sparse_normal_cholesky_solver.h b/internal/ceres/sparse_normal_cholesky_solver.h
index 6572835..d1efe00 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.h
+++ b/internal/ceres/sparse_normal_cholesky_solver.h
@@ -34,6 +34,8 @@
#ifndef CERES_INTERNAL_SPARSE_NORMAL_CHOLESKY_SOLVER_H_
#define CERES_INTERNAL_SPARSE_NORMAL_CHOLESKY_SOLVER_H_
+#include <vector>
+
// This include must come before any #ifndef check on Ceres compile options.
#include "ceres/internal/port.h"
@@ -76,7 +78,7 @@
const LinearSolver::PerSolveOptions& options,
double* rhs_and_solution);
- // Crashes if CERES_USE_LGPGL_CODE is not defined.
+ // Crashes if CERES_USE_EIGEN_SPARSE is not defined.
LinearSolver::Summary SolveImplUsingEigen(
CompressedRowSparseMatrix* A,
const LinearSolver::PerSolveOptions& options,
@@ -94,8 +96,16 @@
#ifdef CERES_USE_EIGEN_SPARSE
typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>,
- Eigen::Upper> SimplicialLDLT;
- scoped_ptr<SimplicialLDLT> simplicial_ldlt_;
+ Eigen::Upper,
+ Eigen::NaturalOrdering<int> >
+ SimplicialLDLTWithNaturalOrdering;
+ typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>,
+ Eigen::Upper,
+ Eigen::AMDOrdering<int> >
+ SimplicialLDLTWithAMDOrdering;
+
+ scoped_ptr<SimplicialLDLTWithNaturalOrdering> natural_ldlt_;
+ scoped_ptr<SimplicialLDLTWithAMDOrdering> amd_ldlt_;
#endif
scoped_ptr<CompressedRowSparseMatrix> outer_product_;
