Enable usage of schur power series expansion preconditioner. Add an option to use schur power series expansion for initialization of pcg solution in ITERATIVE_SCHUR linear solver. Change-Id: Ifb8bce02bc5f5ceebc74f961eefd3f6dd2ffab4a
diff --git a/include/ceres/solver.h b/include/ceres/solver.h index 2a1dd18..27198ac 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h
@@ -691,6 +691,23 @@ // as its termination type. int max_linear_solver_iterations = 500; + // Maximum number of iterations performed by SCHUR_POWER_SERIES_EXPANSION. + // This value controls the maximum number of iterations whether it is used + // as a preconditioner or just to initialize the solution for + // ITERATIVE_SCHUR. + int max_num_spse_iterations = 5; + + // Use SCHUR_POWER_SERIES_EXPANSION to initialize the solution for + // ITERATIVE_SCHUR. This option can be set true regardless of what + // preconditioner is being used. + bool use_power_series_expansion_initialization = false; + + // When use_power_series_expansion_initialization is true, this parameter + // along with max_num_spse_iterations controls the number of + // SCHUR_POWER_SERIES_EXPANSION iterations performed for initialization. It + // is not used to control the preconditioner. + double spse_tolerance = 0.1; + // Forcing sequence parameter. The truncated Newton solver uses // this number to control the relative accuracy with which the // Newton step is computed.
diff --git a/include/ceres/types.h b/include/ceres/types.h index 52e1e7b..1b07ab3 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h
@@ -107,9 +107,6 @@ // Use power series expansion to approximate the inversion of Schur complement // as a preconditioner. - // WARNING! Application of this preconditioner currently is not integrated - // into linear solvers, so failure to use it via public API is expected - // behaviour. SCHUR_POWER_SERIES_EXPANSION, // Visibility clustering based preconditioners.
diff --git a/internal/ceres/implicit_schur_complement.cc b/internal/ceres/implicit_schur_complement.cc index 92631de..8f7a2a0 100644 --- a/internal/ceres/implicit_schur_complement.cc +++ b/internal/ceres/implicit_schur_complement.cc
@@ -57,6 +57,7 @@ b_ = b; compute_ftf_inverse_ = + options_.use_power_series_expansion_initialization || options_.preconditioner_type == JACOBI || options_.preconditioner_type == SCHUR_POWER_SERIES_EXPANSION;
diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc index 5f56be2..8be33bc 100644 --- a/internal/ceres/iterative_schur_complement_solver.cc +++ b/internal/ceres/iterative_schur_complement_solver.cc
@@ -43,6 +43,7 @@ #include "ceres/implicit_schur_complement.h" #include "ceres/internal/eigen.h" #include "ceres/linear_solver.h" +#include "ceres/power_series_expansion_preconditioner.h" #include "ceres/preconditioner.h" #include "ceres/schur_jacobi_preconditioner.h" #include "ceres/triplet_sparse_matrix.h" @@ -90,9 +91,18 @@ return summary; } - // Initialize the solution to the Schur complement system to zero. + // Initialize the solution to the Schur complement system. reduced_linear_system_solution_.resize(schur_complement_->num_rows()); reduced_linear_system_solution_.setZero(); + if (options_.use_power_series_expansion_initialization) { + PowerSeriesExpansionPreconditioner pse_solver( + schur_complement_.get(), + options_.max_num_spse_iterations, + options_.spse_tolerance); + pse_solver.RightMultiplyAndAccumulate( + schur_complement_->rhs().data(), + reduced_linear_system_solution_.data()); + } CreatePreconditioner(A); if (preconditioner_ != nullptr) { @@ -168,6 +178,12 @@ preconditioner_ = std::make_unique<SparseMatrixPreconditionerWrapper>( schur_complement_->block_diagonal_FtF_inverse()); break; + case SCHUR_POWER_SERIES_EXPANSION: + // Ignoring the value of spse_tolerance to ensure preconditioner stays + // fixed during the iterations of cg. + preconditioner_ = std::make_unique<PowerSeriesExpansionPreconditioner>( + schur_complement_.get(), options_.max_num_spse_iterations, 0); + break; case SCHUR_JACOBI: preconditioner_ = std::make_unique<SchurJacobiPreconditioner>( *A->block_structure(), preconditioner_options);
diff --git a/internal/ceres/iterative_schur_complement_solver_test.cc b/internal/ceres/iterative_schur_complement_solver_test.cc index 80e388a..a98b295 100644 --- a/internal/ceres/iterative_schur_complement_solver_test.cc +++ b/internal/ceres/iterative_schur_complement_solver_test.cc
@@ -74,7 +74,9 @@ num_eliminate_blocks_ = problem->num_eliminate_blocks; } - AssertionResult TestSolver(double* D) { + AssertionResult TestSolver(double* D, + PreconditionerType preconditioner_type, + bool use_power_series_expansion_initialization) { TripletSparseMatrix triplet_A( A_->num_rows(), A_->num_cols(), A_->num_nonzeros()); A_->ToTripletSparseMatrix(&triplet_A); @@ -95,7 +97,10 @@ options.elimination_groups.push_back(num_eliminate_blocks_); options.elimination_groups.push_back(0); options.max_num_iterations = num_cols_; - options.preconditioner_type = SCHUR_JACOBI; + options.max_num_spse_iterations = 1; + options.use_power_series_expansion_initialization = + use_power_series_expansion_initialization; + options.preconditioner_type = preconditioner_type; IterativeSchurComplementSolver isc(options); Vector isc_sol(num_cols_); @@ -119,16 +124,30 @@ std::unique_ptr<double[]> D_; }; -TEST_F(IterativeSchurComplementSolverTest, NormalProblem) { +TEST_F(IterativeSchurComplementSolverTest, NormalProblemSchurJacobi) { SetUpProblem(2); - EXPECT_TRUE(TestSolver(nullptr)); - EXPECT_TRUE(TestSolver(D_.get())); + EXPECT_TRUE(TestSolver(nullptr, SCHUR_JACOBI, false)); + EXPECT_TRUE(TestSolver(D_.get(), SCHUR_JACOBI, false)); +} + +TEST_F(IterativeSchurComplementSolverTest, + NormalProblemSchurJacobiWithPowerSeriesExpansionInitialization) { + SetUpProblem(2); + EXPECT_TRUE(TestSolver(nullptr, SCHUR_JACOBI, true)); + EXPECT_TRUE(TestSolver(D_.get(), SCHUR_JACOBI, true)); +} + +TEST_F(IterativeSchurComplementSolverTest, + NormalProblemPowerSeriesExpansionPreconditioner) { + SetUpProblem(5); + EXPECT_TRUE(TestSolver(nullptr, SCHUR_POWER_SERIES_EXPANSION, false)); + EXPECT_TRUE(TestSolver(D_.get(), SCHUR_POWER_SERIES_EXPANSION, false)); } TEST_F(IterativeSchurComplementSolverTest, ProblemWithNoFBlocks) { SetUpProblem(3); - EXPECT_TRUE(TestSolver(nullptr)); - EXPECT_TRUE(TestSolver(D_.get())); + EXPECT_TRUE(TestSolver(nullptr, SCHUR_JACOBI, false)); + EXPECT_TRUE(TestSolver(D_.get(), SCHUR_JACOBI, false)); } } // namespace internal
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h index 4916e41..69e984f 100644 --- a/internal/ceres/linear_solver.h +++ b/internal/ceres/linear_solver.h
@@ -165,6 +165,23 @@ int min_num_iterations = 1; int max_num_iterations = 1; + // Maximum number of iterations performed by SCHUR_POWER_SERIES_EXPANSION. + // This value controls the maximum number of iterations whether it is used + // as a preconditioner or just to initialize the solution for + // ITERATIVE_SCHUR. + int max_num_spse_iterations = 5; + + // Use SCHUR_POWER_SERIES_EXPANSION to initialize the solution for + // ITERATIVE_SCHUR. This option can be set true regardless of what + // preconditioner is being used. + bool use_power_series_expansion_initialization = false; + + // When use_power_series_expansion_initialization is true, this parameter + // along with max_num_spse_iterations controls the number of + // SCHUR_POWER_SERIES_EXPANSION iterations performed for initialization. It + // is not used to control the preconditioner. + double spse_tolerance = 0.1; + // If possible, how many threads can the solver use. int num_threads = 1;
diff --git a/internal/ceres/power_series_expansion_preconditioner.cc b/internal/ceres/power_series_expansion_preconditioner.cc index 1b2dbb9..7a36d92 100644 --- a/internal/ceres/power_series_expansion_preconditioner.cc +++ b/internal/ceres/power_series_expansion_preconditioner.cc
@@ -34,13 +34,11 @@ PowerSeriesExpansionPreconditioner::PowerSeriesExpansionPreconditioner( const ImplicitSchurComplement* isc, - const double spse_tolerance, - const int min_num_iterations, - const int max_num_iterations) + const int max_num_spse_iterations, + const double spse_tolerance) : isc_(isc), - spse_tolerance_(spse_tolerance), - min_num_iterations_(min_num_iterations), - max_num_iterations_(max_num_iterations) {} + max_num_spse_iterations_(max_num_spse_iterations), + spse_tolerance_(spse_tolerance) {} PowerSeriesExpansionPreconditioner::~PowerSeriesExpansionPreconditioner() = default; @@ -66,8 +64,7 @@ isc_->InversePowerSeriesOperatorRightMultiplyAccumulate( previous_series_term.data(), series_term.data()); yref += series_term; - if (i >= min_num_iterations_ && - (i >= max_num_iterations_ || series_term.norm() < norm_threshold)) { + if (i >= max_num_spse_iterations_ || series_term.norm() < norm_threshold) { break; } std::swap(previous_series_term, series_term);
diff --git a/internal/ceres/power_series_expansion_preconditioner.h b/internal/ceres/power_series_expansion_preconditioner.h index c5cf32d..a8bb9a6 100644 --- a/internal/ceres/power_series_expansion_preconditioner.h +++ b/internal/ceres/power_series_expansion_preconditioner.h
@@ -46,9 +46,8 @@ : public Preconditioner { public: PowerSeriesExpansionPreconditioner(const ImplicitSchurComplement* isc, - const double spse_tolerance_, - const int min_num_iterations_, - const int max_num_iterations_); + const int max_num_spse_iterations, + const double spse_tolerance); PowerSeriesExpansionPreconditioner( const PowerSeriesExpansionPreconditioner&) = delete; void operator=(const PowerSeriesExpansionPreconditioner&) = delete; @@ -60,9 +59,8 @@ private: const ImplicitSchurComplement* isc_; + const int max_num_spse_iterations_; const double spse_tolerance_; - const int min_num_iterations_; - const int max_num_iterations_; }; } // namespace ceres::internal
diff --git a/internal/ceres/power_series_expansion_preconditioner_test.cc b/internal/ceres/power_series_expansion_preconditioner_test.cc index a101831..a291405 100644 --- a/internal/ceres/power_series_expansion_preconditioner_test.cc +++ b/internal/ceres/power_series_expansion_preconditioner_test.cc
@@ -81,10 +81,9 @@ TEST_F(PowerSeriesExpansionPreconditionerTest, InverseValidPreconditionerToleranceReached) { const double spse_tolerance = kEpsilon; - const int min_iterations = 1; - const int max_iterations = 50; + const int max_num_iterations = 50; PowerSeriesExpansionPreconditioner preconditioner( - isc_.get(), spse_tolerance, min_iterations, max_iterations); + isc_.get(), max_num_iterations, spse_tolerance); Vector x(num_f_cols_), y(num_f_cols_); for (int i = 0; i < num_f_cols_; i++) { @@ -104,10 +103,10 @@ TEST_F(PowerSeriesExpansionPreconditionerTest, InverseValidPreconditionerMaxIterations) { - const double spse_tolerance = 1 / kEpsilon; - const int num_iterations = 50; + const double spse_tolerance = 0; + const int max_num_iterations = 50; PowerSeriesExpansionPreconditioner preconditioner_fixed_n_iterations( - isc_.get(), spse_tolerance, num_iterations, num_iterations); + isc_.get(), max_num_iterations, spse_tolerance); Vector x(num_f_cols_), y(num_f_cols_); for (int i = 0; i < num_f_cols_; i++) { @@ -129,10 +128,9 @@ TEST_F(PowerSeriesExpansionPreconditionerTest, InverseInvalidBadPreconditionerTolerance) { const double spse_tolerance = 1 / kEpsilon; - const int min_iterations = 1; - const int max_iterations = 50; + const int max_num_iterations = 50; PowerSeriesExpansionPreconditioner preconditioner_bad_tolerance( - isc_.get(), spse_tolerance, min_iterations, max_iterations); + isc_.get(), max_num_iterations, spse_tolerance); Vector x(num_f_cols_), y(num_f_cols_); for (int i = 0; i < num_f_cols_; i++) { @@ -151,9 +149,9 @@ TEST_F(PowerSeriesExpansionPreconditionerTest, InverseInvalidBadPreconditionerMaxIterations) { const double spse_tolerance = kEpsilon; - const int num_iterations = 1; + const int max_num_iterations = 1; PowerSeriesExpansionPreconditioner preconditioner_bad_iterations_limit( - isc_.get(), spse_tolerance, num_iterations, num_iterations); + isc_.get(), max_num_iterations, spse_tolerance); Vector x(num_f_cols_), y(num_f_cols_); for (int i = 0; i < num_f_cols_; i++) {
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc index b07d4f0..71b4537 100644 --- a/internal/ceres/solver.cc +++ b/internal/ceres/solver.cc
@@ -259,12 +259,25 @@ return false; } - if (options.use_explicit_schur_complement && - options.preconditioner_type != SCHUR_JACOBI) { - *error = - "use_explicit_schur_complement only supports " - "SCHUR_JACOBI as the preconditioner."; - return false; + if (options.use_explicit_schur_complement) { + if (options.preconditioner_type != SCHUR_JACOBI) { + *error = + "use_explicit_schur_complement only supports " + "SCHUR_JACOBI as the preconditioner."; + return false; + } + if (options.use_power_series_expansion_initialization) { + *error = + "use_explicit_schur_complement does not support " + "use_power_series_expansion_initialization."; + return false; + } + } + + if (options.use_power_series_expansion_initialization || + options.preconditioner_type == SCHUR_POWER_SERIES_EXPANSION) { + OPTION_GE(max_num_spse_iterations, 1) + OPTION_GE(spse_tolerance, 0.0) } if (options.use_mixed_precision_solves) {
diff --git a/internal/ceres/trust_region_preprocessor.cc b/internal/ceres/trust_region_preprocessor.cc index a5a986f..e9406d1 100644 --- a/internal/ceres/trust_region_preprocessor.cc +++ b/internal/ceres/trust_region_preprocessor.cc
@@ -199,6 +199,11 @@ options.max_linear_solver_iterations; pp->linear_solver_options.type = options.linear_solver_type; pp->linear_solver_options.preconditioner_type = options.preconditioner_type; + pp->linear_solver_options.use_power_series_expansion_initialization = + options.use_power_series_expansion_initialization; + pp->linear_solver_options.spse_tolerance = options.spse_tolerance; + pp->linear_solver_options.max_num_spse_iterations = + options.max_num_spse_iterations; pp->linear_solver_options.visibility_clustering_type = options.visibility_clustering_type; pp->linear_solver_options.sparse_linear_algebra_library_type =
diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc index c0e3355..1d514de 100644 --- a/internal/ceres/types.cc +++ b/internal/ceres/types.cc
@@ -81,6 +81,7 @@ CASESTR(IDENTITY); CASESTR(JACOBI); CASESTR(SCHUR_JACOBI); + CASESTR(SCHUR_POWER_SERIES_EXPANSION); CASESTR(CLUSTER_JACOBI); CASESTR(CLUSTER_TRIDIAGONAL); CASESTR(SUBSET); @@ -94,6 +95,7 @@ STRENUM(IDENTITY); STRENUM(JACOBI); STRENUM(SCHUR_JACOBI); + STRENUM(SCHUR_POWER_SERIES_EXPANSION); STRENUM(CLUSTER_JACOBI); STRENUM(CLUSTER_TRIDIAGONAL); STRENUM(SUBSET);