Refactor SolverImpl::CreateReducedProgram. Break up CreateReducedProgram into smaller functions in preparation for more sophisticated ordering strategies. Change-Id: Ic3897522574fde770646d747fe383f5dbd7a6619
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc index ffc347e..add28ae 100644 --- a/internal/ceres/solver_impl.cc +++ b/internal/ceres/solver_impl.cc
@@ -38,8 +38,8 @@ #include "ceres/gradient_checking_cost_function.h" #include "ceres/iteration_callback.h" #include "ceres/levenberg_marquardt_strategy.h" -#include "ceres/linear_solver.h" #include "ceres/line_search_minimizer.h" +#include "ceres/linear_solver.h" #include "ceres/map_util.h" #include "ceres/minimizer.h" #include "ceres/ordered_groups.h" @@ -505,19 +505,6 @@ summary->trust_region_strategy_type = options.trust_region_strategy_type; summary->dogleg_type = options.dogleg_type; - // Only Schur types require the lexicographic reordering. - if (IsSchurType(options.linear_solver_type)) { - const int num_eliminate_blocks = - options.linear_solver_ordering - ->group_to_elements().begin() - ->second.size(); - if (!LexicographicallyOrderResidualBlocks(num_eliminate_blocks, - reduced_program.get(), - &summary->error)) { - return; - } - } - scoped_ptr<Evaluator> evaluator(CreateEvaluator(options, problem_impl->parameter_map(), reduced_program.get(), @@ -953,11 +940,14 @@ parameter_blocks->resize(j); } - CHECK(((program->NumResidualBlocks() == 0) && + if (!(((program->NumResidualBlocks() == 0) && (program->NumParameterBlocks() == 0)) || ((program->NumResidualBlocks() != 0) && - (program->NumParameterBlocks() != 0))) - << "Congratulations, you found a bug in Ceres. Please report it."; + (program->NumParameterBlocks() != 0)))) { + *error = "Congratulations, you found a bug in Ceres. Please report it."; + return false; + } + return true; } @@ -965,19 +955,14 @@ ProblemImpl* problem_impl, double* fixed_cost, string* error) { - EventLogger event_logger("CreateReducedProgram"); - CHECK_NOTNULL(options->linear_solver_ordering); Program* original_program = problem_impl->mutable_program(); scoped_ptr<Program> transformed_program(new Program(*original_program)); - event_logger.AddEvent("TransformedProgram"); ParameterBlockOrdering* linear_solver_ordering = options->linear_solver_ordering; - const int min_group_id = linear_solver_ordering->group_to_elements().begin()->first; - const int original_num_groups = linear_solver_ordering->NumGroups(); if (!RemoveFixedBlocksFromProgram(transformed_program.get(), linear_solver_ordering, @@ -986,97 +971,39 @@ return NULL; } - event_logger.AddEvent("RemoveFixedBlocks"); - if (transformed_program->NumParameterBlocks() == 0) { - if (transformed_program->NumResidualBlocks() > 0) { - *error = "Zero parameter blocks but non-zero residual blocks" - " in the reduced program. Congratulations, you found a " - "Ceres bug! Please report this error to the developers."; - return NULL; - } - LOG(WARNING) << "No varying parameter blocks to optimize; " << "bailing early."; return transformed_program.release(); } - // If the user supplied an linear_solver_ordering with just one - // group, it is equivalent to the user supplying NULL as - // ordering. Ceres is completely free to choose the parameter block - // ordering as it sees fit. For Schur type solvers, this means that - // the user wishes for Ceres to identify the e_blocks, which we do - // by computing a maximal independent set. - if (original_num_groups == 1 && IsSchurType(options->linear_solver_type)) { - vector<ParameterBlock*> schur_ordering; - const int num_eliminate_blocks = ComputeSchurOrdering(*transformed_program, - &schur_ordering); - CHECK_EQ(schur_ordering.size(), transformed_program->NumParameterBlocks()) - << "Congratulations, you found a Ceres bug! Please report this error " - << "to the developers."; - - for (int i = 0; i < schur_ordering.size(); ++i) { - linear_solver_ordering->AddElementToGroup( - schur_ordering[i]->mutable_user_state(), - (i < num_eliminate_blocks) ? 0 : 1); - } - } - event_logger.AddEvent("SchurOrdering"); - - if (!ApplyUserOrdering(problem_impl->parameter_map(), - linear_solver_ordering, - transformed_program.get(), - error)) { - return NULL; - } - event_logger.AddEvent("ApplyOrdering"); - - // If the user requested the use of a Schur type solver, and - // supplied a non-NULL linear_solver_ordering object with more than - // one elimination group, then it can happen that after all the - // parameter blocks which are fixed or unused have been removed from - // the program and the ordering, there are no more parameter blocks - // in the first elimination group. - // - // In such a case, the use of a Schur type solver is not possible, - // as they assume there is at least one e_block. Thus, we - // automatically switch to one of the other solvers, depending on - // the user's indicated preferences. if (IsSchurType(options->linear_solver_type) && - original_num_groups > 1 && linear_solver_ordering->GroupSize(min_group_id) == 0) { - string msg = "No e_blocks remaining. Switching from "; - if (options->linear_solver_type == SPARSE_SCHUR) { - options->linear_solver_type = SPARSE_NORMAL_CHOLESKY; - msg += "SPARSE_SCHUR to SPARSE_NORMAL_CHOLESKY."; - } else if (options->linear_solver_type == DENSE_SCHUR) { - // TODO(sameeragarwal): This is probably not a great choice. - // Ideally, we should have a DENSE_NORMAL_CHOLESKY, that can - // take a BlockSparseMatrix as input. - options->linear_solver_type = DENSE_QR; - msg += "DENSE_SCHUR to DENSE_QR."; - } else if (options->linear_solver_type == ITERATIVE_SCHUR) { - msg += StringPrintf("ITERATIVE_SCHUR with %s preconditioner " - "to CGNR with JACOBI preconditioner.", - PreconditionerTypeToString( - options->preconditioner_type)); - options->linear_solver_type = CGNR; - if (options->preconditioner_type != IDENTITY) { - // CGNR currently only supports the JACOBI preconditioner. - options->preconditioner_type = JACOBI; - } - } - - LOG(WARNING) << msg; + // If the user requested the use of a Schur type solver, and + // supplied a non-NULL linear_solver_ordering object with more than + // one elimination group, then it can happen that after all the + // parameter blocks which are fixed or unused have been removed from + // the program and the ordering, there are no more parameter blocks + // in the first elimination group. + // + // In such a case, the use of a Schur type solver is not possible, + // as they assume there is at least one e_block. Thus, we + // automatically switch to the closest solver to the one indicated + // by the user. + AlternateLinearSolverForSchurTypeLinearSolver(options); } - event_logger.AddEvent("AlternateSolver"); + if (IsSchurType(options->linear_solver_type)) { + if (!ReorderProgramForSchurTypeLinearSolver(problem_impl->parameter_map(), + linear_solver_ordering, + transformed_program.get(), + error)) { + return NULL; + } + return transformed_program.release(); + } - // Since the transformed program is the "active" program, and it is - // mutated, update the parameter offsets and indices. transformed_program->SetParameterOffsetsAndIndex(); - - event_logger.AddEvent("SetOffsets"); return transformed_program.release(); } @@ -1398,5 +1325,96 @@ return inner_iteration_minimizer.release(); } +void SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver( + Solver::Options* options) { + if (!IsSchurType(options->linear_solver_type)) { + return; + } + + string msg = "No e_blocks remaining. Switching from "; + if (options->linear_solver_type == SPARSE_SCHUR) { + options->linear_solver_type = SPARSE_NORMAL_CHOLESKY; + msg += "SPARSE_SCHUR to SPARSE_NORMAL_CHOLESKY."; + } else if (options->linear_solver_type == DENSE_SCHUR) { + // TODO(sameeragarwal): This is probably not a great choice. + // Ideally, we should have a DENSE_NORMAL_CHOLESKY, that can + // take a BlockSparseMatrix as input. + options->linear_solver_type = DENSE_QR; + msg += "DENSE_SCHUR to DENSE_QR."; + } else if (options->linear_solver_type == ITERATIVE_SCHUR) { + options->linear_solver_type = CGNR; + if (options->preconditioner_type != IDENTITY) { + msg += StringPrintf("ITERATIVE_SCHUR with %s preconditioner " + "to CGNR with JACOBI preconditioner.", + PreconditionerTypeToString( + options->preconditioner_type)); + // CGNR currently only supports the JACOBI preconditioner. + options->preconditioner_type = JACOBI; + } else { + msg += StringPrintf("ITERATIVE_SCHUR with IDENTITY preconditioner " + "to CGNR with IDENTITY preconditioner."); + } + } + LOG(WARNING) << msg; +} + +bool SolverImpl::ReorderProgramForSchurTypeLinearSolver( + const ProblemImpl::ParameterMap& parameter_map, + ParameterBlockOrdering* ordering, + Program* program, + string* error) { + // At this point one of two things is true. + // + // 1. The user did not specify an ordering - ordering has one + // group containined all the parameter blocks. + + // 2. The user specified an ordering, and the first group has + // non-zero elements. + // + // We handle these two cases in turn. + if (ordering->NumGroups() == 1) { + // If the user supplied an ordering with just one + // group, it is equivalent to the user supplying NULL as an + // ordering. Ceres is completely free to choose the parameter + // block ordering as it sees fit. For Schur type solvers, this + // means that the user wishes for Ceres to identify the e_blocks, + // which we do by computing a maximal independent set. + vector<ParameterBlock*> schur_ordering; + const int num_eliminate_blocks = ComputeSchurOrdering(*program, + &schur_ordering); + + CHECK_EQ(schur_ordering.size(), program->NumParameterBlocks()) + << "Congratulations, you found a Ceres bug! Please report this error " + << "to the developers."; + + // Update the ordering object. + for (int i = 0; i < schur_ordering.size(); ++i) { + double* parameter_block = schur_ordering[i]->mutable_user_state(); + const int group_id = (i < num_eliminate_blocks) ? 0 : 1; + ordering->AddElementToGroup(parameter_block, group_id); + } + + // Apply the parameter block re-ordering. Technically we could + // call ApplyUserOrdering, but this is cheaper and simpler. + swap(*program->mutable_parameter_blocks(), schur_ordering); + } else { + // The user supplied an ordering. + if (!ApplyUserOrdering(parameter_map, ordering, program, error)) { + return false; + } + } + + program->SetParameterOffsetsAndIndex(); + + const int num_eliminate_blocks = + ordering->group_to_elements().begin()->second.size(); + + // Schur type solvers also require that their residual blocks be + // lexicographically ordered. + return LexicographicallyOrderResidualBlocks(num_eliminate_blocks, + program, + error); +} + } // namespace internal } // namespace ceres
diff --git a/internal/ceres/solver_impl.h b/internal/ceres/solver_impl.h index c054830..6377e84 100644 --- a/internal/ceres/solver_impl.h +++ b/internal/ceres/solver_impl.h
@@ -151,6 +151,32 @@ const Program& program, const ProblemImpl::ParameterMap& parameter_map, Solver::Summary* summary); + + // If the linear solver is of Schur type, then replace it with the + // closest equivalent linear solver. This is done when the user + // requested a Schur type solver but the problem structure makes it + // impossible to use one. + // + // If the linear solver is not of Schur type, the function is a + // no-op. + static void AlternateLinearSolverForSchurTypeLinearSolver( + Solver::Options* options); + + // Schur type solvers require that all parameter blocks eliminated + // by the Schur eliminator occur before others and the residuals be + // sorted in lexicographic order of their parameter blocks. + // + // If ordering has atleast two groups, then apply the ordering, + // otherwise compute a new ordering using a Maximal Independent Set + // algorithm and apply it. + // + // Upon return, ordering contains the parameter block ordering that + // was used to order the program. + static bool ReorderProgramForSchurTypeLinearSolver( + const ProblemImpl::ParameterMap& parameter_map, + ParameterBlockOrdering* ordering, + Program* program, + string* error); }; } // namespace internal
diff --git a/internal/ceres/solver_impl_test.cc b/internal/ceres/solver_impl_test.cc index 24860d2..a752eff 100644 --- a/internal/ceres/solver_impl_test.cc +++ b/internal/ceres/solver_impl_test.cc
@@ -378,11 +378,6 @@ expected_residual_blocks.push_back(residual_blocks[3]); expected_residual_blocks.push_back(residual_blocks[2]); - EXPECT_TRUE(SolverImpl::LexicographicallyOrderResidualBlocks( - 2, - reduced_program.get(), - &error)); - EXPECT_EQ(reduced_program->residual_blocks().size(), expected_residual_blocks.size()); for (int i = 0; i < expected_residual_blocks.size(); ++i) { @@ -794,5 +789,62 @@ EXPECT_EQ(summary.final_cost, 1.0 / 2.0); } +TEST(SolverImpl, AlternateLinearSolverForSchurTypeLinearSolver) { + Solver::Options options; + + options.linear_solver_type = DENSE_QR; + SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver(&options); + EXPECT_EQ(options.linear_solver_type, DENSE_QR); + + options.linear_solver_type = DENSE_NORMAL_CHOLESKY; + SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver(&options); + EXPECT_EQ(options.linear_solver_type, DENSE_NORMAL_CHOLESKY); + + options.linear_solver_type = SPARSE_NORMAL_CHOLESKY; + SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver(&options); + EXPECT_EQ(options.linear_solver_type, SPARSE_NORMAL_CHOLESKY); + + options.linear_solver_type = CGNR; + SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver(&options); + EXPECT_EQ(options.linear_solver_type, CGNR); + + options.linear_solver_type = DENSE_SCHUR; + SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver(&options); + EXPECT_EQ(options.linear_solver_type, DENSE_QR); + + options.linear_solver_type = SPARSE_SCHUR; + SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver(&options); + EXPECT_EQ(options.linear_solver_type, SPARSE_NORMAL_CHOLESKY); + + options.linear_solver_type = ITERATIVE_SCHUR; + options.preconditioner_type = IDENTITY; + SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver(&options); + EXPECT_EQ(options.linear_solver_type, CGNR); + EXPECT_EQ(options.preconditioner_type, IDENTITY); + + options.linear_solver_type = ITERATIVE_SCHUR; + options.preconditioner_type = JACOBI; + SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver(&options); + EXPECT_EQ(options.linear_solver_type, CGNR); + EXPECT_EQ(options.preconditioner_type, JACOBI); + + options.linear_solver_type = ITERATIVE_SCHUR; + options.preconditioner_type = SCHUR_JACOBI; + SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver(&options); + EXPECT_EQ(options.linear_solver_type, CGNR); + EXPECT_EQ(options.preconditioner_type, JACOBI); + + options.linear_solver_type = ITERATIVE_SCHUR; + options.preconditioner_type = CLUSTER_JACOBI; + SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver(&options); + EXPECT_EQ(options.linear_solver_type, CGNR); + EXPECT_EQ(options.preconditioner_type, JACOBI); + + options.linear_solver_type = ITERATIVE_SCHUR; + options.preconditioner_type = CLUSTER_TRIDIAGONAL; + SolverImpl::AlternateLinearSolverForSchurTypeLinearSolver(&options); + EXPECT_EQ(options.linear_solver_type, CGNR); + EXPECT_EQ(options.preconditioner_type, JACOBI); +} } // namespace internal } // namespace ceres
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc index 7ece965..0cd064b 100644 --- a/internal/ceres/sparse_normal_cholesky_solver.cc +++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -205,11 +205,11 @@ if (factor_ == NULL) { if (options_.use_block_amd) { - factor_ = ss_.BlockAnalyzeCholesky(lhs.get(), + factor_ = ss_.BlockAnalyzeCholesky(&lhs, A->col_blocks(), A->row_blocks()); } else { - factor_ = ss_.AnalyzeCholesky(lhs.get()); + factor_ = ss_.AnalyzeCholesky(&lhs); } } event_logger.AddEvent("Analysis");