Add Nested Dissection based fill reducing ordering

With this change, the user can now choose between Approximate Minimum Degree and Nested Dissection as a fill reducing algorithm when using a sparse direct factorization based linear solver like SPARSE_NORMAL_CHOLESKY or SPARSE_SCHUR. Currenly only SUITE_SPARSE is supported. It requires that SuiteSparse be compiled with Metis support enabled. On most problems AMD is still the better choice, but in some cases like the grid3D dataset from https://lucacarlone.mit.edu/datasets/ the solution time with AMD is 57s and with NESDIS 38 on my M1 Mac. On some other problems at Google we have observed speedups of 10x, there is also a corresponding decrease in the total amount of memory used. This patch is based on the original work done by NeroBurner in https://ceres-solver-review.googlesource.com/c/ceres-solver/+/20580 1. Add a new enum to the public api LinearSolverOrderingType and a setting Solver::Options::linear_solver_ordering_type. 2. TrustRegionPreprocessor had some complicated logic which determined when linear solvers should reorder their matrices on their own and not this has been refactored into a more readable function that lives inside reorder_program.h/cc. 3. Plumbing in reorder_program.cc and trust_region_processor.cc to use nested dissection. 4. Update bundle_adjuster.cc to use nested dissection. Change-Id: I388b027934f86c58b4da2b65a4fa5204ea73bf40
2026-08-29 16:40:38 +08:00 · 2022-05-14 08:05:49 -07:00
parent aa62dd86a8
commit 39ec5e8f99
10 changed files with 387 additions and 114 deletions
@@ -85,7 +85,6 @@ static int MinParameterBlock(const ResidualBlock* residual_block,
  return min_parameter_block_position;
 }

-#if defined(CERES_USE_EIGEN_SPARSE)
 Eigen::SparseMatrix<int> CreateBlockJacobian(
    const TripletSparseMatrix& block_jacobian_transpose) {
  using SparseMatrix = Eigen::SparseMatrix<int>;
@@ -105,9 +104,9 @@ Eigen::SparseMatrix<int> CreateBlockJacobian(
  block_jacobian.setFromTriplets(triplets.begin(), triplets.end());
  return block_jacobian;
 }
-#endif

 void OrderingForSparseNormalCholeskyUsingSuiteSparse(
+    const LinearSolverOrderingType linear_solver_ordering_type,
    const TripletSparseMatrix& tsm_block_jacobian_transpose,
    const vector<ParameterBlock*>& parameter_blocks,
    const ParameterBlockOrdering& parameter_block_ordering,
@@ -120,40 +119,55 @@ void OrderingForSparseNormalCholeskyUsingSuiteSparse(
  cholmod_sparse* block_jacobian_transpose = ss.CreateSparseMatrix(
      const_cast<TripletSparseMatrix*>(&tsm_block_jacobian_transpose));

-  // If the user did not supply a useful ordering, then just use
-  // regular AMD.
-  if (parameter_block_ordering.NumGroups() <= 1) {
-    ss.Ordering(block_jacobian_transpose, OrderingType::AMD, &ordering[0]);
-  } else {
-    vector<int> constraints;
-    for (auto* parameter_block : parameter_blocks) {
-      constraints.push_back(parameter_block_ordering.GroupId(
-          parameter_block->mutable_user_state()));
+  if (linear_solver_ordering_type == ceres::AMD) {
+    if (parameter_block_ordering.NumGroups() <= 1) {
+      // The user did not supply a useful ordering so just go ahead
+      // and use AMD.
+      ss.Ordering(block_jacobian_transpose, OrderingType::AMD, &ordering[0]);
+    } else {
+      // The user supplied an ordering, so use CAMD.
+      vector<int> constraints;
+      constraints.reserve(parameter_blocks.size());
+      for (auto* parameter_block : parameter_blocks) {
+        constraints.push_back(parameter_block_ordering.GroupId(
+            parameter_block->mutable_user_state()));
+      }
+
+      // Renumber the entries of constraints to be contiguous integers
+      // as CAMD requires that the group ids be in the range [0,
+      // parameter_blocks.size() - 1].
+      MapValuesToContiguousRange(constraints.size(), &constraints[0]);
+      ss.ConstrainedApproximateMinimumDegreeOrdering(
+          block_jacobian_transpose, &constraints[0], ordering);
    }
-
-    // Renumber the entries of constraints to be contiguous integers
-    // as CAMD requires that the group ids be in the range [0,
-    // parameter_blocks.size() - 1].
-    MapValuesToContiguousRange(constraints.size(), &constraints[0]);
-    ss.ConstrainedApproximateMinimumDegreeOrdering(
-        block_jacobian_transpose, &constraints[0], ordering);
+  } else if (linear_solver_ordering_type == ceres::NESDIS) {
+    // If nested dissection is chosen as an ordering algorithm, then
+    // ignore any user provided linear_solver_ordering.
+    CHECK(SuiteSparse::IsNestedDissectionAvailable())
+        << "Congratulations, you found a Ceres bug! "
+        << "Please report this error to the developers.";
+    ss.Ordering(block_jacobian_transpose, OrderingType::NESDIS, &ordering[0]);
+  } else {
+    LOG(FATAL) << "Congratulations, you found a Ceres bug! "
+               << "Please report this error to the developers.";
  }

-  VLOG(2) << "Block ordering stats: "
-          << " flops: " << ss.mutable_cc()->fl
-          << " lnz  : " << ss.mutable_cc()->lnz
-          << " anz  : " << ss.mutable_cc()->anz;
-
  ss.Free(block_jacobian_transpose);
 #endif  // CERES_NO_SUITESPARSE
 }

 void OrderingForSparseNormalCholeskyUsingCXSparse(
-    const TripletSparseMatrix& tsm_block_jacobian_transpose, int* ordering) {
+    const LinearSolverOrderingType linear_solver_ordering_type,
+    const TripletSparseMatrix& tsm_block_jacobian_transpose,
+    int* ordering) {
 #ifdef CERES_NO_CXSPARSE
  LOG(FATAL) << "Congratulations, you found a Ceres bug! "
             << "Please report this error to the developers.";
 #else
+  CHECK_NE(linear_solver_ordering_type, NESDIS)
+      << "Congratulations, you found a Ceres bug! "
+      << "Please report this error to the developers.";
+
  // CXSparse works with J'J instead of J'. So compute the block
  // sparsity for J'J and compute an approximate minimum degree
  // ordering.
@@ -173,13 +187,18 @@ void OrderingForSparseNormalCholeskyUsingCXSparse(
 }

 void OrderingForSparseNormalCholeskyUsingEigenSparse(
-    const TripletSparseMatrix& tsm_block_jacobian_transpose, int* ordering) {
+    const LinearSolverOrderingType linear_solver_ordering_type,
+    const TripletSparseMatrix& tsm_block_jacobian_transpose,
+    int* ordering) {
 #ifndef CERES_USE_EIGEN_SPARSE
  LOG(FATAL) << "SPARSE_NORMAL_CHOLESKY cannot be used with EIGEN_SPARSE "
                "because Ceres was not built with support for "
                "Eigen's SimplicialLDLT decomposition. "
                "This requires enabling building with -DEIGENSPARSE=ON.";
 #else
+  CHECK_NE(linear_solver_ordering_type, NESDIS)
+      << "Congratulations, you found a Ceres bug! "
+      << "Please report this error to the developers.";

  // This conversion from a TripletSparseMatrix to a Eigen::Triplet
  // matrix is unfortunate, but unavoidable for now. It is not a
@@ -326,7 +345,7 @@ bool LexicographicallyOrderResidualBlocks(

 // Pre-order the columns corresponding to the Schur complement if
 // possible.
-static void MaybeReorderSchurComplementColumnsUsingSuiteSparse(
+static void ReorderSchurComplementColumnsUsingSuiteSparse(
    const ParameterBlockOrdering& parameter_block_ordering, Program* program) {
 #ifndef CERES_NO_SUITESPARSE
  SuiteSparse ss;
@@ -365,7 +384,7 @@ static void MaybeReorderSchurComplementColumnsUsingSuiteSparse(
 #endif
 }

-static void MaybeReorderSchurComplementColumnsUsingEigen(
+static void ReorderSchurComplementColumnsUsingEigen(
    const int size_of_first_elimination_group,
    const ProblemImpl::ParameterMap& parameter_map,
    Program* program) {
@@ -420,6 +439,7 @@ static void MaybeReorderSchurComplementColumnsUsingEigen(
 bool ReorderProgramForSchurTypeLinearSolver(
    const LinearSolverType linear_solver_type,
    const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+    const LinearSolverOrderingType linear_solver_ordering_type,
    const ProblemImpl::ParameterMap& parameter_map,
    ParameterBlockOrdering* parameter_block_ordering,
    Program* program,
@@ -485,12 +505,16 @@ bool ReorderProgramForSchurTypeLinearSolver(
  const int size_of_first_elimination_group =
      parameter_block_ordering->group_to_elements().begin()->second.size();

-  if (linear_solver_type == SPARSE_SCHUR) {
+  // Pre-ordering of the columns of the Schur complement only works if
+  // we are using approximate mininmum degree based ordering and
+  // SUITE_SPARSE or EIGEN_SPARSE.
+  if (linear_solver_type == SPARSE_SCHUR &&
+      linear_solver_ordering_type == ceres::AMD) {
    if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
-      MaybeReorderSchurComplementColumnsUsingSuiteSparse(
-          *parameter_block_ordering, program);
+      ReorderSchurComplementColumnsUsingSuiteSparse(*parameter_block_ordering,
+                                                    program);
    } else if (sparse_linear_algebra_library_type == EIGEN_SPARSE) {
-      MaybeReorderSchurComplementColumnsUsingEigen(
+      ReorderSchurComplementColumnsUsingEigen(
          size_of_first_elimination_group, parameter_map, program);
    }
  }
@@ -503,6 +527,7 @@ bool ReorderProgramForSchurTypeLinearSolver(

 bool ReorderProgramForSparseCholesky(
    const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+    const LinearSolverOrderingType linear_solver_ordering_type,
    const ParameterBlockOrdering& parameter_block_ordering,
    int start_row_block,
    Program* program,
@@ -526,12 +551,14 @@ bool ReorderProgramForSparseCholesky(

  if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
    OrderingForSparseNormalCholeskyUsingSuiteSparse(
+        linear_solver_ordering_type,
        *tsm_block_jacobian_transpose,
        parameter_blocks,
        parameter_block_ordering,
        &ordering[0]);
  } else if (sparse_linear_algebra_library_type == CX_SPARSE) {
-    OrderingForSparseNormalCholeskyUsingCXSparse(*tsm_block_jacobian_transpose,
+    OrderingForSparseNormalCholeskyUsingCXSparse(linear_solver_ordering_type,
+                                                 *tsm_block_jacobian_transpose,
                                                 &ordering[0]);
  } else if (sparse_linear_algebra_library_type == ACCELERATE_SPARSE) {
    // Accelerate does not provide a function to perform reordering without
@@ -544,7 +571,9 @@ bool ReorderProgramForSparseCholesky(

  } else if (sparse_linear_algebra_library_type == EIGEN_SPARSE) {
    OrderingForSparseNormalCholeskyUsingEigenSparse(
-        *tsm_block_jacobian_transpose, &ordering[0]);
+        linear_solver_ordering_type,
+        *tsm_block_jacobian_transpose,
+        &ordering[0]);
  }

  // Apply ordering.
@@ -569,4 +598,54 @@ int ReorderResidualBlocksByPartition(
  return it - residual_blocks->begin();
 }

+bool AreJacobianColumnsOrdered(
+    const LinearSolverType linear_solver_type,
+    const PreconditionerType preconditioner_type,
+    const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
+    const LinearSolverOrderingType linear_solver_ordering_type,
+    const bool use_postordering,
+    const bool dynamic_sparsity) {
+  if (use_postordering || dynamic_sparsity) {
+    return false;
+  }
+
+  if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
+    if (linear_solver_type == SPARSE_NORMAL_CHOLESKY ||
+        (linear_solver_type == CGNR && preconditioner_type == SUBSET)) {
+      return true;
+    }
+    if (linear_solver_type == SPARSE_SCHUR &&
+        linear_solver_ordering_type == ceres::AMD) {
+      return true;
+    }
+  }
+
+  // For all sparse linear algebra libraries other than SuiteSparse,
+  // nested dissection is not used for pre-ordering.
+  if (linear_solver_ordering_type == ceres::NESDIS) {
+    return false;
+  }
+
+  if (sparse_linear_algebra_library_type == ceres::EIGEN_SPARSE) {
+    if (linear_solver_type == SPARSE_NORMAL_CHOLESKY ||
+        linear_solver_type == SPARSE_SCHUR ||
+        (linear_solver_type == CGNR && preconditioner_type == SUBSET)) {
+      return true;
+    }
+  }
+
+  if (sparse_linear_algebra_library_type == ceres::CX_SPARSE) {
+    if (linear_solver_type == SPARSE_NORMAL_CHOLESKY ||
+        (linear_solver_type == CGNR && preconditioner_type == SUBSET)) {
+      return true;
+    }
+  }
+
+  if (sparse_linear_algebra_library_type == ceres::ACCELERATE_SPARSE) {
+    return false;
+  }
+
+  return false;
+}
+
 }  // namespace ceres::internal