From 39ec5e8f99c595226ab3cd760d4945e61c31fc57 Mon Sep 17 00:00:00 2001 From: Sameer Agarwal Date: Sat, 14 May 2022 08:05:49 -0700 Subject: [PATCH] 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 --- examples/bundle_adjuster.cc | 43 +++--- include/ceres/solver.h | 130 +++++++++++++--- include/ceres/types.h | 29 ++++ internal/ceres/reorder_program.cc | 145 ++++++++++++++---- internal/ceres/reorder_program.h | 13 ++ internal/ceres/reorder_program_test.cc | 1 + internal/ceres/solver.cc | 83 ++++++---- internal/ceres/trust_region_preprocessor.cc | 38 +++-- internal/ceres/types.cc | 17 ++ .../ceres/visibility_based_preconditioner.cc | 2 +- 10 files changed, 387 insertions(+), 114 deletions(-) diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index d03c2ff59..87e18e110 100644 --- a/examples/bundle_adjuster.cc +++ b/examples/bundle_adjuster.cc @@ -95,7 +95,7 @@ DEFINE_string(sparse_linear_algebra_library, "suite_sparse", "Options are: suite_sparse and cx_sparse."); DEFINE_string(dense_linear_algebra_library, "eigen", "Options are: eigen, lapack, and cuda"); -DEFINE_string(ordering, "automatic", "Options are: automatic, user."); +DEFINE_string(ordering, "amd", "Options are: amd, nesdis and user"); DEFINE_bool(use_quaternions, false, "If true, uses quaternions to represent " "rotations. If false, angle axis is used."); @@ -128,7 +128,6 @@ DEFINE_int32(max_num_refinement_iterations, 0, "Iterative refinement iterations" DEFINE_string(initial_ply, "", "Export the BAL file data as a PLY file."); DEFINE_string(final_ply, "", "Export the refined BAL file data as a PLY " "file."); - // clang-format on namespace ceres::examples { @@ -228,24 +227,30 @@ void SetOrdering(BALProblem* bal_problem, Solver::Options* options) { // the right ParameterBlock ordering, or by manually specifying a // suitable ordering vector and defining // Options::num_eliminate_blocks. - if (CERES_GET_FLAG(FLAGS_ordering) == "automatic") { - return; + if (CERES_GET_FLAG(FLAGS_ordering) == "user") { + auto* ordering = new ceres::ParameterBlockOrdering; + + // The points come before the cameras. + for (int i = 0; i < num_points; ++i) { + ordering->AddElementToGroup(points + point_block_size * i, 0); + } + + for (int i = 0; i < num_cameras; ++i) { + // When using axis-angle, there is a single parameter block for + // the entire camera. + ordering->AddElementToGroup(cameras + camera_block_size * i, 1); + } + + options->linear_solver_ordering.reset(ordering); + } else { + if (CERES_GET_FLAG(FLAGS_ordering) == "amd") { + options->linear_solver_ordering_type = AMD; + } else if (CERES_GET_FLAG(FLAGS_ordering) == "nesdis") { + options->linear_solver_ordering_type = NESDIS; + } else { + LOG(FATAL) << "Unknown ordering type: " << CERES_GET_FLAG(FLAGS_ordering); + } } - - auto* ordering = new ceres::ParameterBlockOrdering; - - // The points come before the cameras. - for (int i = 0; i < num_points; ++i) { - ordering->AddElementToGroup(points + point_block_size * i, 0); - } - - for (int i = 0; i < num_cameras; ++i) { - // When using axis-angle, there is a single parameter block for - // the entire camera. - ordering->AddElementToGroup(cameras + camera_block_size * i, 1); - } - - options->linear_solver_ordering.reset(ordering); } void SetMinimizerOptions(Solver::Options* options) { diff --git a/include/ceres/solver.h b/include/ceres/solver.h index 026fc1c08..43146b110 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h @@ -395,12 +395,33 @@ class CERES_EXPORT Solver { #endif // The order in which variables are eliminated in a linear solver - // can have a significant of impact on the efficiency and accuracy - // of the method. e.g., when doing sparse Cholesky factorization, + // can have a significant impact on the efficiency and accuracy of + // the method. e.g., when doing sparse Cholesky factorization, // there are matrices for which a good ordering will give a // Cholesky factor with O(n) storage, where as a bad ordering will // result in an completely dense factor. // + // Sparse direct solvers like SPARSE_NORMAL_CHOLESKY and + // SPARSE_SCHUR use a fill reducing ordering of the columns and + // rows of the matrix being factorized before computing the + // numeric factorization. + // + // This enum controls the type of algorithm used to compute + // this fill reducing ordering. There is no single algorithm + // that works on all matrices, so determining which algorithm + // works better is a matter of empirical experimentation. + // + // Implementation status: + // + // AMD works for SUITE_SPARSE, CX_SPARSE, EIGEN_SPARSE & + // ACCELERATE_SPARSE. + // + // NESDIS currently works for SUITE_SPARSE when using + // SPARSE_NORMAL_CHOLESKY, SPARSE_SCHUR, CGNR + SUBSET, + // ITERATIVE_SCHUR + CLUSTER_JACOBI, ITERATIVE_SCHUR + + // CLUSTER_TRIDIAGONAL linear solvers. + LinearSolverOrderingType linear_solver_ordering_type = AMD; + // Ceres allows the user to provide varying amounts of hints to // the solver about the variable elimination ordering to use. This // can range from no hints, where the solver is free to decide the @@ -418,15 +439,47 @@ class CERES_EXPORT Solver { // associated with it, that determines its order in the set of // groups. // - // Given such an ordering, Ceres ensures that the parameter blocks in - // the lowest numbered group are eliminated first, and then the - // parameter blocks in the next lowest numbered group and so on. Within - // each group, Ceres is free to order the parameter blocks as it - // chooses. + // The exact interpretation of this information depends on the + // values of linear_solver_ordering_type and + // linear_solver_type/preconditioner_type and + // sparse_linear_algebra_type. // - // If nullptr, then all parameter blocks are assumed to be in the - // same group and the solver is free to decide the best - // ordering. + // sparse_linear_algebra_library_type != SUITE_SPARSE + // ================================================== + // + // The value of linear_solver_ordering_type does not matter and: + // + // a. linear_solver_type = SPARSE_NORMAL_CHOLESKY or + // linear_solver_type = CGNR and preconditioner_type = SUBSET + // + // then the value of linear_solver_ordering is ignored. + // + // b. linear_solver_type = SPARSE_SCHUR/DENSE_SCHUR/ITERATIVE_SCHUR + // + // then if linear_solver_ordering is non-null, then the lowest + // number group is used as the first elimination group to compute + // the Schur complement. All other groups are ignored. + // + // sparse_linear_algebra_library_type == SUITE_SPARSE + // ================================================== + // + // linear_solver_ordering_type = AMD + // --------------------------------- + // + // linear_solver_type = SPARSE_NORMAL_CHOLESKY or + // linear_solver_type = CGNR and preconditioner_type = SUBSET + // + // if linear_solver_ordering = nullptr, then Ceres assumes that + // all parameter blocks are in the same elimination group and uses + // the Approximate Minimum Degree algorithm to compute a fill + // reducing ordering. + // + // If linear_solver_order is not null, then a Constrained + // Approximate Minimum Degree (CAMD) ordering used where the + // parameter blocks in the lowest numbered group are eliminated + // first, and then the parameter blocks in the next lowest + // numbered group and so on. Within each group, CAMD free to order + // the parameter blocks as it chooses. // // e.g. Consider the linear system // @@ -442,18 +495,53 @@ class CERES_EXPORT Solver { // {0: y}, {1: x} - eliminate y first. // {0: x, y} - Solver gets to decide the elimination order. // - // Thus, to have Ceres determine the ordering automatically using - // heuristics, put all the variables in group 0 and to control the - // ordering for every variable, create groups 0..N-1, one per - // variable, in the desired order. + // Thus, to have Ceres determine the ordering automatically, put + // all the variables in group 0 and to control the ordering for + // every variable, create groups 0..N-1, one per variable, in the + // desired order. + // + // b. linear_solver_type = SPARSE_SCHUR/DENSE_SCHUR/ITERATIVE_SCHUR + // + // If linear_solver_ordering is null then an greedy maximum + // independent set algorithm is used to compute the first + // elimination group, and AMD is applied to the corresponding + // Schur complement matrix or a subset thereof used that is used + // as a preconditioner. + // + // If linear_solver_ordering is not null, then CAMD is used to + // compute a fill reducing ordering. + // + // linear_solver_ordering = NESDIS + // ------------------------------- + // + // a. linear_solver_type = SPARSE_NORMAL_CHOLESKY or + // linear_solver_type = CGNR and preconditioner_type = SUBSET + // + // then the value of linear_solver_ordering is ignored and a + // Nested Dissection algorithm is used to compute a fill reducing + // ordering. + // + // b. linear_solver_type = SPARSE_SCHUR/DENSE_SCHUR/ITERATIVE_SCHUR + // + // If linear_solver_ordering is null then an greedy maximum + // independent set algorithm is used to compute the first + // elimination group, and Nested Dissection is applied to the + // corresponding Schur complement matrix or a subset thereof used + // that is used as a preconditioner. + // + // If linear_solver_ordering is not null, the parameter blocks in + // the lowest numbered group are eliminated to compute the Schur + // complement. All other groups are ignored and Nested Dissection + // is applied to the corresponding Schur complement matrix or a + // subset thereof used that is used as a preconditioner. // // Bundle Adjustment // ----------------- // - // A particular case of interest is bundle adjustment, where the user - // has two options. The default is to not specify an ordering at all, - // the solver will see that the user wants to use a Schur type solver - // and figure out the right elimination ordering. + // A particular case of interest is bundle adjustment, where the + // user has two options. The default is not specifying ordering at + // all; in this case the solver will see that the user wants to + // use a Schur type solver and figure out an elimination ordering. // // But if the user already knows what parameter blocks are points and // what are cameras, they can save preprocessing time by partitioning @@ -500,6 +588,10 @@ class CERES_EXPORT Solver { // Jacobian matrix and generally speaking, there is no performance // penalty for doing so. + // TODO(sameeragarwal): Remove this option. It is too obscure and + // there is no clear way of figuring out when this is a useful + // thing to do. + // // In some rare cases, it is worth using a more complicated // reordering algorithm which has slightly better runtime // performance at the expense of an extra copy of the Jacobian @@ -946,6 +1038,8 @@ class CERES_EXPORT Solver { SPARSE_NORMAL_CHOLESKY; #endif + LinearSolverOrderingType linear_solver_ordering_type; + // Size of the elimination groups given by the user as hints to // the linear solver. std::vector linear_solver_ordering_given; diff --git a/include/ceres/types.h b/include/ceres/types.h index e52242381..2fc5abccd 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h @@ -183,6 +183,30 @@ enum SparseLinearAlgebraLibraryType { NO_SPARSE }; +// The order in which variables are eliminated in a linear solver +// can have a significant of impact on the efficiency and accuracy +// of the method. e.g., when doing sparse Cholesky factorization, +// there are matrices for which a good ordering will give a +// Cholesky factor with O(n) storage, where as a bad ordering will +// result in an completely dense factor. +// +// So sparse direct solvers like SPARSE_NORMAL_CHOLESKY and +// SPARSE_SCHUR and preconditioners like SUBSET, CLUSTER_JACOBI & +// CLUSTER_TRIDIAGONAL use a fill reducing ordering of the columns and +// rows of the matrix being factorized before actually the numeric +// factorization. +// +// This enum controls the class of algorithm used to compute this +// fill reducing ordering. There is no single algorithm that works +// on all matrices, so determining which algorithm works better is a +// matter of empirical experimentation. +enum LinearSolverOrderingType { + // Approximate Minimum Degree. + AMD, + // Nested Dissection. + NESDIS +}; + enum DenseLinearAlgebraLibraryType { EIGEN, LAPACK, @@ -467,6 +491,11 @@ CERES_EXPORT const char* SparseLinearAlgebraLibraryTypeToString( CERES_EXPORT bool StringToSparseLinearAlgebraLibraryType( std::string value, SparseLinearAlgebraLibraryType* type); +CERES_EXPORT const char* LinearSolverOrderingTypeToString( + LinearSolverOrderingType type); +CERES_EXPORT bool StringToLinearSolverOrderingType( + std::string value, LinearSolverOrderingType* type); + CERES_EXPORT const char* DenseLinearAlgebraLibraryTypeToString( DenseLinearAlgebraLibraryType type); CERES_EXPORT bool StringToDenseLinearAlgebraLibraryType( diff --git a/internal/ceres/reorder_program.cc b/internal/ceres/reorder_program.cc index c6765fef7..1f50a39ef 100644 --- a/internal/ceres/reorder_program.cc +++ b/internal/ceres/reorder_program.cc @@ -85,7 +85,6 @@ static int MinParameterBlock(const ResidualBlock* residual_block, return min_parameter_block_position; } -#if defined(CERES_USE_EIGEN_SPARSE) Eigen::SparseMatrix CreateBlockJacobian( const TripletSparseMatrix& block_jacobian_transpose) { using SparseMatrix = Eigen::SparseMatrix; @@ -105,9 +104,9 @@ Eigen::SparseMatrix 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& parameter_blocks, const ParameterBlockOrdering& parameter_block_ordering, @@ -120,40 +119,55 @@ void OrderingForSparseNormalCholeskyUsingSuiteSparse( cholmod_sparse* block_jacobian_transpose = ss.CreateSparseMatrix( const_cast(&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 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 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 diff --git a/internal/ceres/reorder_program.h b/internal/ceres/reorder_program.h index 64f409504..8ae456ec4 100644 --- a/internal/ceres/reorder_program.h +++ b/internal/ceres/reorder_program.h @@ -35,6 +35,7 @@ #include "ceres/internal/disable_warnings.h" #include "ceres/internal/export.h" +#include "ceres/linear_solver.h" #include "ceres/parameter_block_ordering.h" #include "ceres/problem_impl.h" #include "ceres/types.h" @@ -75,6 +76,7 @@ CERES_NO_EXPORT bool LexicographicallyOrderResidualBlocks( CERES_NO_EXPORT bool ReorderProgramForSchurTypeLinearSolver( LinearSolverType linear_solver_type, SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type, + LinearSolverOrderingType linear_solver_ordering_type, const ProblemImpl::ParameterMap& parameter_map, ParameterBlockOrdering* parameter_block_ordering, Program* program, @@ -92,6 +94,7 @@ CERES_NO_EXPORT bool ReorderProgramForSchurTypeLinearSolver( // ordering will take it into account, otherwise it will be ignored. CERES_NO_EXPORT bool ReorderProgramForSparseCholesky( SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type, + LinearSolverOrderingType linear_solver_ordering_type, const ParameterBlockOrdering& parameter_block_ordering, int start_row_block, Program* program, @@ -111,6 +114,16 @@ CERES_NO_EXPORT int ReorderResidualBlocksByPartition( const std::unordered_set& bottom_residual_blocks, Program* program); +// The return value of this function indicates whether the columns of +// the Jacobian can be reordered using a fill reducing ordering. +CERES_NO_EXPORT bool AreJacobianColumnsOrdered( + LinearSolverType linear_solver_type, + PreconditionerType preconditioner_type, + SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type, + LinearSolverOrderingType linear_solver_ordering_type, + bool use_postordering, + bool dynamic_sparsity); + } // namespace ceres::internal #include "ceres/internal/reenable_warnings.h" diff --git a/internal/ceres/reorder_program_test.cc b/internal/ceres/reorder_program_test.cc index 0fbab4e10..872f3812c 100644 --- a/internal/ceres/reorder_program_test.cc +++ b/internal/ceres/reorder_program_test.cc @@ -185,6 +185,7 @@ class ReorderProgramForSparseCholeskyUsingSuiteSparseTest std::string error; EXPECT_TRUE(ReorderProgramForSparseCholesky(ceres::SUITE_SPARSE, + ceres::AMD, linear_solver_ordering, 0, /* use all rows */ program, diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc index 150c5550f..dbaba3e89 100644 --- a/internal/ceres/solver.cc +++ b/internal/ceres/solver.cc @@ -50,6 +50,7 @@ #include "ceres/schur_templates.h" #include "ceres/solver_utils.h" #include "ceres/stringprintf.h" +#include "ceres/suitesparse.h" #include "ceres/types.h" #include "ceres/wall_time.h" @@ -107,6 +108,11 @@ bool CommonOptionsAreValid(const Solver::Options& options, string* error) { return true; } +bool IsNestedDissectionAvailable(SparseLinearAlgebraLibraryType type) { + return (type == SUITE_SPARSE) && + internal::SuiteSparse::IsNestedDissectionAvailable(); +} + bool TrustRegionOptionsAreValid(const Solver::Options& options, string* error) { OPTION_GT(initial_trust_region_radius, 0.0); OPTION_GT(min_trust_region_radius, 0.0); @@ -182,8 +188,10 @@ bool TrustRegionOptionsAreValid(const Solver::Options& options, string* error) { name, sparse_linear_algebra_library_name); return false; - } else if (!IsSparseLinearAlgebraLibraryTypeAvailable( - options.sparse_linear_algebra_library_type)) { + } + + if (!IsSparseLinearAlgebraLibraryTypeAvailable( + options.sparse_linear_algebra_library_type)) { *error = StringPrintf( "Can't use %s with " "Solver::Options::sparse_linear_algebra_library_type = %s, " @@ -192,6 +200,22 @@ bool TrustRegionOptionsAreValid(const Solver::Options& options, string* error) { sparse_linear_algebra_library_name); return false; } + + if (options.linear_solver_ordering_type == ceres::NESDIS && + !IsNestedDissectionAvailable( + options.sparse_linear_algebra_library_type)) { + if (options.sparse_linear_algebra_library_type == SUITE_SPARSE) { + *error = StringPrintf( + "Can't use NESDIS with SUITE_SPARSE because SuiteSparse was " + "compiled without support for Metis."); + } else { + *error = StringPrintf( + "Can't use NESDIS with " + "Solver::Options::sparse_linear_algebra_library_type = %s.", + sparse_linear_algebra_library_name); + } + return false; + } } } @@ -339,7 +363,7 @@ void PreSolveSummarize(const Solver::Options& options, &(summary->inner_iteration_ordering_given)); // clang-format off - summary->dense_linear_algebra_library_type = options.dense_linear_algebra_library_type; // NOLINT + summary->dense_linear_algebra_library_type = options.dense_linear_algebra_library_type; summary->dogleg_type = options.dogleg_type; summary->inner_iteration_time_in_seconds = 0.0; summary->num_line_search_steps = 0; @@ -348,18 +372,19 @@ void PreSolveSummarize(const Solver::Options& options, summary->line_search_polynomial_minimization_time_in_seconds = 0.0; summary->line_search_total_time_in_seconds = 0.0; summary->inner_iterations_given = options.use_inner_iterations; - summary->line_search_direction_type = options.line_search_direction_type; // NOLINT - summary->line_search_interpolation_type = options.line_search_interpolation_type; // NOLINT + summary->line_search_direction_type = options.line_search_direction_type; + summary->line_search_interpolation_type = options.line_search_interpolation_type; summary->line_search_type = options.line_search_type; summary->linear_solver_type_given = options.linear_solver_type; summary->max_lbfgs_rank = options.max_lbfgs_rank; summary->minimizer_type = options.minimizer_type; - summary->nonlinear_conjugate_gradient_type = options.nonlinear_conjugate_gradient_type; // NOLINT + summary->nonlinear_conjugate_gradient_type = options.nonlinear_conjugate_gradient_type; summary->num_threads_given = options.num_threads; summary->preconditioner_type_given = options.preconditioner_type; - summary->sparse_linear_algebra_library_type = options.sparse_linear_algebra_library_type; // NOLINT - summary->trust_region_strategy_type = options.trust_region_strategy_type; // NOLINT - summary->visibility_clustering_type = options.visibility_clustering_type; // NOLINT + summary->sparse_linear_algebra_library_type = options.sparse_linear_algebra_library_type; + summary->linear_solver_ordering_type = options.linear_solver_ordering_type; + summary->trust_region_strategy_type = options.trust_region_strategy_type; + summary->visibility_clustering_type = options.visibility_clustering_type; // clang-format on } @@ -367,11 +392,14 @@ void PostSolveSummarize(const internal::PreprocessedProblem& pp, Solver::Summary* summary) { internal::OrderingToGroupSizes(pp.options.linear_solver_ordering.get(), &(summary->linear_solver_ordering_used)); + // TODO(sameeragarwal): Update the preprocessor to collapse the + // second and higher groups into one group when nested dissection is + // used. internal::OrderingToGroupSizes(pp.options.inner_iteration_ordering.get(), &(summary->inner_iteration_ordering_used)); // clang-format off - summary->inner_iterations_used = pp.inner_iteration_minimizer.get() != nullptr; // NOLINT + summary->inner_iterations_used = pp.inner_iteration_minimizer != nullptr; summary->linear_solver_type_used = pp.linear_solver_options.type; summary->num_threads_used = pp.options.num_threads; summary->preconditioner_type_used = pp.options.preconditioner_type; @@ -379,7 +407,7 @@ void PostSolveSummarize(const internal::PreprocessedProblem& pp, internal::SetSummaryFinalCost(summary); - if (pp.reduced_program.get() != nullptr) { + if (pp.reduced_program != nullptr) { SummarizeReducedProgram(*pp.reduced_program, summary); } @@ -389,7 +417,7 @@ void PostSolveSummarize(const internal::PreprocessedProblem& pp, // case if the preprocessor failed, or if the reduced problem did // not contain any parameter blocks. Thus, only extract the // evaluator statistics if one exists. - if (pp.evaluator.get() != nullptr) { + if (pp.evaluator != nullptr) { const map& evaluator_statistics = pp.evaluator->Statistics(); { @@ -411,7 +439,7 @@ void PostSolveSummarize(const internal::PreprocessedProblem& pp, // Again, like the evaluator, there may or may not be a linear // solver from which we can extract run time statistics. In // particular the line search solver does not use a linear solver. - if (pp.linear_solver.get() != nullptr) { + if (pp.linear_solver != nullptr) { const map& linear_solver_statistics = pp.linear_solver->Statistics(); const CallStatistics& call_stats = FindWithDefault( @@ -666,17 +694,6 @@ string Solver::Summary::FullReport() const { dense_linear_algebra_library_type)); } - if (linear_solver_type_used == SPARSE_NORMAL_CHOLESKY || - linear_solver_type_used == SPARSE_SCHUR || - (linear_solver_type_used == ITERATIVE_SCHUR && - (preconditioner_type_used == CLUSTER_JACOBI || - preconditioner_type_used == CLUSTER_TRIDIAGONAL))) { - StringAppendF(&report, - "\nSparse linear algebra library %15s\n", - SparseLinearAlgebraLibraryTypeToString( - sparse_linear_algebra_library_type)); - } - StringAppendF(&report, "Trust region strategy %19s", TrustRegionStrategyTypeToString(trust_region_strategy_type)); @@ -687,9 +704,23 @@ string Solver::Summary::FullReport() const { StringAppendF(&report, " (SUBSPACE)"); } } - StringAppendF(&report, "\n"); - StringAppendF(&report, "\n"); + if (linear_solver_type_used == SPARSE_NORMAL_CHOLESKY || + linear_solver_type_used == SPARSE_SCHUR || + (linear_solver_type_used == CGNR && + preconditioner_type_used == SUBSET) || + (linear_solver_type_used == ITERATIVE_SCHUR && + (preconditioner_type_used == CLUSTER_JACOBI || + preconditioner_type_used == CLUSTER_TRIDIAGONAL))) { + StringAppendF( + &report, + "\nSparse linear algebra library %15s + %s\n", + SparseLinearAlgebraLibraryTypeToString( + sparse_linear_algebra_library_type), + LinearSolverOrderingTypeToString(linear_solver_ordering_type)); + } + + StringAppendF(&report, "\n"); StringAppendF(&report, "%45s %21s\n", "Given", "Used"); StringAppendF(&report, "Linear solver %25s%25s\n", diff --git a/internal/ceres/trust_region_preprocessor.cc b/internal/ceres/trust_region_preprocessor.cc index 9467b5ac6..ed62aa9d8 100644 --- a/internal/ceres/trust_region_preprocessor.cc +++ b/internal/ceres/trust_region_preprocessor.cc @@ -113,6 +113,7 @@ bool ReorderProgram(PreprocessedProblem* pp) { return ReorderProgramForSchurTypeLinearSolver( options.linear_solver_type, options.sparse_linear_algebra_library_type, + options.linear_solver_ordering_type, pp->problem->parameter_map(), options.linear_solver_ordering.get(), pp->reduced_program.get(), @@ -123,6 +124,7 @@ bool ReorderProgram(PreprocessedProblem* pp) { !options.dynamic_sparsity) { return ReorderProgramForSparseCholesky( options.sparse_linear_algebra_library_type, + options.linear_solver_ordering_type, *options.linear_solver_ordering, 0, /* use all the rows of the jacobian */ pp->reduced_program.get(), @@ -138,6 +140,7 @@ bool ReorderProgram(PreprocessedProblem* pp) { return ReorderProgramForSparseCholesky( options.sparse_linear_algebra_library_type, + options.linear_solver_ordering_type, *options.linear_solver_ordering, pp->linear_solver_options.subset_preconditioner_start_row_block, pp->reduced_program.get(), @@ -200,6 +203,7 @@ bool SetupLinearSolver(PreprocessedProblem* pp) { options.visibility_clustering_type; pp->linear_solver_options.sparse_linear_algebra_library_type = options.sparse_linear_algebra_library_type; + pp->linear_solver_options.dense_linear_algebra_library_type = options.dense_linear_algebra_library_type; pp->linear_solver_options.use_explicit_schur_complement = @@ -210,10 +214,6 @@ bool SetupLinearSolver(PreprocessedProblem* pp) { pp->linear_solver_options.max_num_refinement_iterations = options.max_num_refinement_iterations; pp->linear_solver_options.num_threads = options.num_threads; - pp->linear_solver_options.ordering_type = OrderingType::NATURAL; - if (options.use_postordering) { - pp->linear_solver_options.ordering_type = OrderingType::AMD; - } pp->linear_solver_options.context = pp->problem->context(); if (IsSchurType(pp->linear_solver_options.type)) { @@ -227,20 +227,24 @@ bool SetupLinearSolver(PreprocessedProblem* pp) { if (pp->linear_solver_options.elimination_groups.size() == 1) { pp->linear_solver_options.elimination_groups.push_back(0); } + } - if (options.linear_solver_type == SPARSE_SCHUR) { - // When using SPARSE_SCHUR, we ignore the user's postordering - // preferences if CX_SPARSE is the sparse linear algebra - // backend. - // - // This ensures that the linear solver does not assume that a - // fill-reducing pre-ordering has been done. - // - // TODO(sameeragarwal): Implement the reordering of parameter - // blocks for CX_SPARSE. - if (options.sparse_linear_algebra_library_type == CX_SPARSE) { - pp->linear_solver_options.ordering_type = OrderingType::AMD; - } + if (AreJacobianColumnsOrdered(options.linear_solver_type, + options.preconditioner_type, + options.sparse_linear_algebra_library_type, + options.linear_solver_ordering_type, + options.use_postordering, + options.dynamic_sparsity)) { + pp->linear_solver_options.ordering_type = OrderingType::NATURAL; + } else { + if (options.linear_solver_ordering_type == ceres::AMD) { + pp->linear_solver_options.ordering_type = OrderingType::AMD; + } else if (options.linear_solver_ordering_type == ceres::NESDIS) { + pp->linear_solver_options.ordering_type = OrderingType::NESDIS; + } else { + LOG(FATAL) << "Congratulations you have found a bug in Ceres Solver." + << " Please report this to the maintainers. : " + << options.linear_solver_ordering_type; } } diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc index 48242678b..eea031237 100644 --- a/internal/ceres/types.cc +++ b/internal/ceres/types.cc @@ -124,6 +124,23 @@ bool StringToSparseLinearAlgebraLibraryType( return false; } +const char* LinearSolverOrderingTypeToString(LinearSolverOrderingType type) { + switch (type) { + CASESTR(AMD); + CASESTR(NESDIS); + default: + return "UNKNOWN"; + } +} + +bool StringToLinearSolverOrderingType(string value, + LinearSolverOrderingType* type) { + UpperCase(&value); + STRENUM(AMD); + STRENUM(NESDIS); + return false; +} + const char* DenseLinearAlgebraLibraryTypeToString( DenseLinearAlgebraLibraryType type) { switch (type) { diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc index bde281406..0a7b19f97 100644 --- a/internal/ceres/visibility_based_preconditioner.cc +++ b/internal/ceres/visibility_based_preconditioner.cc @@ -106,7 +106,7 @@ VisibilityBasedPreconditioner::VisibilityBasedPreconditioner( LinearSolver::Options sparse_cholesky_options; sparse_cholesky_options.sparse_linear_algebra_library_type = options_.sparse_linear_algebra_library_type; - sparse_cholesky_options.ordering_type = options.ordering_type; + sparse_cholesky_options.ordering_type = options_.ordering_type; sparse_cholesky_ = SparseCholesky::Create(sparse_cholesky_options); const time_t init_time = time(nullptr);