From 9189f4ea4bb2d71ea7f6b9d9bd3290415aee323d Mon Sep 17 00:00:00 2001 From: Sameer Agarwal Date: Fri, 19 Apr 2013 17:09:49 -0700 Subject: [PATCH] Enable pre-ordering for SPARSE_NORMAL_CHOLESKY. Sparse Cholesky factorization algorithms use a fill-reducing ordering to permute the columns of the Jacobian matrix. There are two ways of doing this. 1. Compute the Jacobian matrix in some order and then have the factorization algorithm permute the columns of the Jacobian. 2. Compute the Jacobian with its columns already permuted. The first option incurs a significant memory penalty. The factorization algorithm has to make a copy of the permuted Jacobian matrix. Starting with this change Ceres pre-permutes the columns of the Jacobian matrix and generally speaking, there is no performance penalty for doing so. 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 matrix. Setting Solver::Options::use_postordering to true enables this tradeoff. This change also removes Solver::Options::use_block_amd as an option. All matrices are ordered using their block structure. The ability to order them by their scalar sparsity structure has been removed. Here is what performance on looks like on some BAL problems. Memory ====== HEAD pre-ordering 16-22106 137957376.0 113516544.0 49-7776 56688640.0 46628864.0 245-198739 1718005760.0 1383550976.0 257-65132 387715072.0 319512576.0 356-226730 2014826496.0 1626087424.0 744-543562 4903358464.0 3957878784.0 1024-110968 968626176.0 822071296.0 Time ==== HEAD pre-ordering 16-22106 3.8 3.7 49-7776 1.9 1.8 245-198739 82.6 81.9 257-65132 14.0 13.4 356-226730 98.8 95.8 744-543562 325.2 301.6 1024-110968 42.1 37.1 Change-Id: I6b2e25f3fed7310f88905386a7898ac94d37467e --- examples/bundle_adjuster.cc | 5 --- include/ceres/solver.h | 43 +++++++++++-------- .../iterative_schur_complement_solver.cc | 1 - internal/ceres/linear_solver.h | 6 +-- internal/ceres/preconditioner.h | 4 -- internal/ceres/schur_complement_solver.cc | 6 +-- internal/ceres/solver_impl.cc | 6 +-- .../ceres/sparse_normal_cholesky_solver.cc | 5 ++- .../ceres/unsymmetric_linear_solver_test.cc | 1 - .../ceres/visibility_based_preconditioner.cc | 6 +-- 10 files changed, 35 insertions(+), 48 deletions(-) diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index 8fff9f0bc..59420c91a 100644 --- a/examples/bundle_adjuster.cc +++ b/examples/bundle_adjuster.cc @@ -96,9 +96,6 @@ DEFINE_double(eta, 1e-2, "Default value for eta. Eta determines the " "accuracy of each linear solve of the truncated newton step. " "Changing this parameter can affect solve performance."); -DEFINE_bool(use_block_amd, true, "Use a block oriented fill reducing " - "ordering."); - DEFINE_int32(num_threads, 1, "Number of threads."); DEFINE_int32(num_iterations, 5, "Number of iterations."); DEFINE_double(max_solver_time, 1e32, "Maximum solve time in seconds."); @@ -139,8 +136,6 @@ void SetOrdering(BALProblem* bal_problem, Solver::Options* options) { const int camera_block_size = bal_problem->camera_block_size(); double* cameras = bal_problem->mutable_cameras(); - options->use_block_amd = FLAGS_use_block_amd; - if (options->use_inner_iterations) { if (FLAGS_blocks_for_inner_iterations == "cameras") { LOG(INFO) << "Camera blocks for inner iterations"; diff --git a/include/ceres/solver.h b/include/ceres/solver.h index 8c2ff32d8..97d082d80 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h @@ -95,13 +95,8 @@ class Solver { #endif num_linear_solver_threads = 1; - -#if defined(CERES_NO_SUITESPARSE) - use_block_amd = false; -#else - use_block_amd = true; -#endif linear_solver_ordering = NULL; + use_postordering = false; use_inner_iterations = false; inner_iteration_ordering = NULL; linear_solver_min_num_iterations = 1; @@ -356,19 +351,29 @@ class Solver { // deallocate the memory when destroyed. ParameterBlockOrdering* linear_solver_ordering; - // By virtue of the modeling layer in Ceres being block oriented, - // all the matrices used by Ceres are also block oriented. When - // doing sparse direct factorization of these matrices (for - // SPARSE_NORMAL_CHOLESKY, SPARSE_SCHUR and ITERATIVE in - // conjunction with CLUSTER_TRIDIAGONAL AND CLUSTER_JACOBI - // preconditioners), the fill-reducing ordering algorithms can - // either be run on the block or the scalar form of these matrices. - // Running it on the block form exposes more of the super-nodal - // structure of the matrix to the factorization routines. Setting - // this parameter to true runs the ordering algorithms in block - // form. Currently this option only makes sense with - // sparse_linear_algebra_library = SUITE_SPARSE. - bool use_block_amd; + // Note: This option only applies to the SPARSE_NORMAL_CHOLESKY + // solver when used with SUITE_SPARSE. + + // Sparse Cholesky factorization algorithms use a fill-reducing + // ordering to permute the columns of the Jacobian matrix. There + // are two ways of doing this. + + // 1. Compute the Jacobian matrix in some order and then have the + // factorization algorithm permute the columns of the Jacobian. + + // 2. Compute the Jacobian with its columns already permuted. + + // The first option incurs a significant memory penalty. The + // factorization algorithm has to make a copy of the permuted + // Jacobian matrix, thus Ceres pre-permutes the columns of the + // Jacobian matrix and generally speaking, there is no performance + // penalty for doing so. + + // 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 + // // matrix. Setting use_postordering to true enables this tradeoff. + bool use_postordering; // Some non-linear least squares problems have additional // structure in the way the parameter blocks interact that it is diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc index cf5562ac7..15e0bdcd8 100644 --- a/internal/ceres/iterative_schur_complement_solver.cc +++ b/internal/ceres/iterative_schur_complement_solver.cc @@ -99,7 +99,6 @@ LinearSolver::Summary IterativeSchurComplementSolver::SolveImpl( preconditioner_options.type = options_.preconditioner_type; preconditioner_options.sparse_linear_algebra_library = options_.sparse_linear_algebra_library; - preconditioner_options.use_block_amd = options_.use_block_amd; preconditioner_options.num_threads = options_.num_threads; preconditioner_options.row_block_size = options_.row_block_size; preconditioner_options.e_block_size = options_.e_block_size; diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h index f4bd0fb6f..ca10faa24 100644 --- a/internal/ceres/linear_solver.h +++ b/internal/ceres/linear_solver.h @@ -74,7 +74,7 @@ class LinearSolver { : type(SPARSE_NORMAL_CHOLESKY), preconditioner_type(JACOBI), sparse_linear_algebra_library(SUITE_SPARSE), - use_block_amd(true), + use_postordering(false), min_num_iterations(1), max_num_iterations(1), num_threads(1), @@ -90,8 +90,8 @@ class LinearSolver { SparseLinearAlgebraLibraryType sparse_linear_algebra_library; - // See solver.h for explanation of this option. - bool use_block_amd; + // See solver.h for information about this flag. + bool use_postordering; // Number of internal iterations that the solver uses. This // parameter only makes sense for iterative solvers like CG. diff --git a/internal/ceres/preconditioner.h b/internal/ceres/preconditioner.h index bfc8464db..d7c882936 100644 --- a/internal/ceres/preconditioner.h +++ b/internal/ceres/preconditioner.h @@ -47,7 +47,6 @@ class Preconditioner : public LinearOperator { Options() : type(JACOBI), sparse_linear_algebra_library(SUITE_SPARSE), - use_block_amd(true), num_threads(1), row_block_size(Eigen::Dynamic), e_block_size(Eigen::Dynamic), @@ -58,9 +57,6 @@ class Preconditioner : public LinearOperator { SparseLinearAlgebraLibraryType sparse_linear_algebra_library; - // See solver.h for explanation of this option. - bool use_block_amd; - // If possible, how many threads the preconditioner can use. int num_threads; diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc index 875f45510..8afb12150 100644 --- a/internal/ceres/schur_complement_solver.cc +++ b/internal/ceres/schur_complement_solver.cc @@ -286,11 +286,7 @@ bool SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse( // Symbolic factorization is computed if we don't already have one handy. if (factor_ == NULL) { - if (options().use_block_amd) { - factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs, blocks_, blocks_); - } else { - factor_ = ss_.AnalyzeCholesky(cholmod_lhs); - } + factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs, blocks_, blocks_); } cholmod_dense* cholmod_solution = diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc index 0d7b83c40..8217484e6 100644 --- a/internal/ceres/solver_impl.cc +++ b/internal/ceres/solver_impl.cc @@ -1004,7 +1004,8 @@ Program* SolverImpl::CreateReducedProgram(Solver::Options* options, return transformed_program.release(); } - if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY) { + if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY && + options->sparse_linear_algebra_library == SUITE_SPARSE) { ReorderProgramForSparseNormalCholesky(transformed_program.get()); return transformed_program.release(); } @@ -1091,11 +1092,10 @@ LinearSolver* SolverImpl::CreateLinearSolver(Solver::Options* options, linear_solver_options.preconditioner_type = options->preconditioner_type; linear_solver_options.sparse_linear_algebra_library = options->sparse_linear_algebra_library; - + linear_solver_options.use_postordering = options->use_postordering; linear_solver_options.num_threads = options->num_linear_solver_threads; options->num_linear_solver_threads = linear_solver_options.num_threads; - linear_solver_options.use_block_amd = options->use_block_amd; const map >& groups = options->linear_solver_ordering->group_to_elements(); for (map >::const_iterator it = groups.begin(); diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc index 0cd064b55..bc1f98334 100644 --- a/internal/ceres/sparse_normal_cholesky_solver.cc +++ b/internal/ceres/sparse_normal_cholesky_solver.cc @@ -204,14 +204,15 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( event_logger.AddEvent("Setup"); if (factor_ == NULL) { - if (options_.use_block_amd) { + if (options_.use_postordering) { factor_ = ss_.BlockAnalyzeCholesky(&lhs, A->col_blocks(), A->row_blocks()); } else { - factor_ = ss_.AnalyzeCholesky(&lhs); + factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs); } } + event_logger.AddEvent("Analysis"); cholmod_dense* sol = ss_.SolveCholesky(&lhs, factor_, rhs); diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/unsymmetric_linear_solver_test.cc index 0b0d593cf..c8a15c026 100644 --- a/internal/ceres/unsymmetric_linear_solver_test.cc +++ b/internal/ceres/unsymmetric_linear_solver_test.cc @@ -62,7 +62,6 @@ class UnsymmetricLinearSolverTest : public ::testing::Test { LinearSolver::Options options; options.type = linear_solver_type; options.sparse_linear_algebra_library = sparse_linear_algebra_library; - options.use_block_amd = false; scoped_ptr solver(LinearSolver::Create(options)); LinearSolver::PerSolveOptions per_solve_options; diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc index f5c46bfcc..4b1e26af6 100644 --- a/internal/ceres/visibility_based_preconditioner.cc +++ b/internal/ceres/visibility_based_preconditioner.cc @@ -421,11 +421,7 @@ bool VisibilityBasedPreconditioner::Factorize() { // Symbolic factorization is computed if we don't already have one handy. if (factor_ == NULL) { - if (options_.use_block_amd) { - factor_ = ss_.BlockAnalyzeCholesky(lhs, block_size_, block_size_); - } else { - factor_ = ss_.AnalyzeCholesky(lhs); - } + factor_ = ss_.BlockAnalyzeCholesky(lhs, block_size_, block_size_); } bool status = ss_.Cholesky(lhs, factor_);