Simplify IterativeRefiner Change the loop structure of IterativeRefiner to unconditionally refine for max_num_iterations. This is done for two reasons. 1. We expect to use this refinement for a small number of iterations where the convergence test is useless. 2. Eliminating the convergence test means we can restructure the loop and save on a sparse matrix-vector multiply, saving precious compute. Change-Id: I6347f453a5d19d234af2a2eb1bce811048963e06
diff --git a/internal/ceres/iterative_refiner.cc b/internal/ceres/iterative_refiner.cc index b73c0b0..fb0e45b 100644 --- a/internal/ceres/iterative_refiner.cc +++ b/internal/ceres/iterative_refiner.cc
@@ -49,68 +49,25 @@ lhs_x_solution_.resize(num_cols); } -IterativeRefiner::Summary IterativeRefiner::Refine( - const SparseMatrix& lhs, - const double* rhs_ptr, - SparseCholesky* sparse_cholesky, - double* solution_ptr) { - Summary summary; +void IterativeRefiner::Refine(const SparseMatrix& lhs, + const double* rhs_ptr, + SparseCholesky* sparse_cholesky, + double* solution_ptr) { const int num_cols = lhs.num_cols(); Allocate(num_cols); - ConstVectorRef rhs(rhs_ptr, num_cols); VectorRef solution(solution_ptr, num_cols); - - summary.lhs_max_norm = ConstVectorRef(lhs.values(), lhs.num_nonzeros()) - .lpNorm<Eigen::Infinity>(); - summary.rhs_max_norm = rhs.lpNorm<Eigen::Infinity>(); - summary.solution_max_norm = solution.lpNorm<Eigen::Infinity>(); - - // residual = rhs - lhs * solution - lhs_x_solution_.setZero(); - lhs.RightMultiply(solution_ptr, lhs_x_solution_.data()); - residual_ = rhs - lhs_x_solution_; - summary.residual_max_norm = residual_.lpNorm<Eigen::Infinity>(); - - for (summary.num_iterations = 0; - summary.num_iterations < max_num_iterations_; - ++summary.num_iterations) { - // Check the current solution for convergence. - const double kTolerance = 5e-15; // From Hogg & Scott. - // residual_tolerance = (|A| |x| + |b|) * kTolerance; - const double residual_tolerance = - (summary.lhs_max_norm * summary.solution_max_norm + - summary.rhs_max_norm) * - kTolerance; - VLOG(3) << "Refinement:" - << " iter: " << summary.num_iterations - << " |A|: " << summary.lhs_max_norm - << " |b|: " << summary.rhs_max_norm - << " |x|: " << summary.solution_max_norm - << " |b - Ax|: " << summary.residual_max_norm - << " tol: " << residual_tolerance; - // |b - Ax| < (|A| |x| + |b|) * kTolerance; - if (summary.residual_max_norm < residual_tolerance) { - summary.converged = true; - break; - } - - // Solve for lhs * correction = residual - correction_.setZero(); - std::string ignored_message; - sparse_cholesky->Solve( - residual_.data(), correction_.data(), &ignored_message); - solution += correction_; - summary.solution_max_norm = solution.lpNorm<Eigen::Infinity>(); - + for (int i = 0; i < max_num_iterations_; ++i) { // residual = rhs - lhs * solution lhs_x_solution_.setZero(); lhs.RightMultiply(solution_ptr, lhs_x_solution_.data()); residual_ = rhs - lhs_x_solution_; - summary.residual_max_norm = residual_.lpNorm<Eigen::Infinity>(); + // solution += lhs^-1 residual + std::string ignored_message; + sparse_cholesky->Solve( + residual_.data(), correction_.data(), &ignored_message); + solution += correction_; } - - return summary; }; } // namespace internal
diff --git a/internal/ceres/iterative_refiner.h b/internal/ceres/iterative_refiner.h index 9021e81..f969935 100644 --- a/internal/ceres/iterative_refiner.h +++ b/internal/ceres/iterative_refiner.h
@@ -38,8 +38,8 @@ namespace ceres { namespace internal { -class SparseMatrix; class SparseCholesky; +class SparseMatrix; // Iterative refinement // (https://en.wikipedia.org/wiki/Iterative_refinement) is the process @@ -53,39 +53,18 @@ // IterativeRefiner implements this process for Symmetric Positive // Definite linear systems. // -// The above iterative loop is run until max_num_iterations is reached -// or the following convergence criterion is satisfied: -// -// |b - Ax| -// ------------- < 5e-15 -// |A| |x| + |b| -// -// All norms in the above expression are max-norms. The above -// expression is what is recommended and used by Hogg & Scott in "A -// fast and robust mixed-precision solver for the solution of sparse -// symmetric linear systems". -// -// For example usage, please see sparse_normal_cholesky_solver.cc +// The above iterative loop is run until max_num_iterations is reached. class IterativeRefiner { public: - struct Summary { - bool converged = false; - int num_iterations = -1; - double lhs_max_norm = -1; - double rhs_max_norm = -1; - double solution_max_norm = -1; - double residual_max_norm = -1; - }; - - // max_num_iterations is the maximum number of refinement iterations - // to perform. + // max_num_iterations is the number of refinement iterations to + // perform. IterativeRefiner(int max_num_iterations); // Needed for mocking. virtual ~IterativeRefiner(); // Given an initial estimate of the solution of lhs * x = rhs, use - // iterative refinement to improve it. + // max_num_iterations rounds of iterative refinement to improve it. // // sparse_cholesky is assumed to contain an already computed // factorization (or approximation thereof) of lhs. @@ -94,14 +73,10 @@ // to lhs * x = rhs. It can be zero. // // This method is virtual to facilitate mocking. - // - // TODO(sameeragarwal): Consider dropping the Summary object, and - // simplifying the internal implementation to improve efficiency, - // since we do not seem to be using the output at all. - virtual Summary Refine(const SparseMatrix& lhs, - const double* rhs, - SparseCholesky* sparse_cholesky, - double* solution); + virtual void Refine(const SparseMatrix& lhs, + const double* rhs, + SparseCholesky* sparse_cholesky, + double* solution); private: void Allocate(int num_cols);
diff --git a/internal/ceres/iterative_refiner_test.cc b/internal/ceres/iterative_refiner_test.cc index 0c9ef69..ecd9a36 100644 --- a/internal/ceres/iterative_refiner_test.cc +++ b/internal/ceres/iterative_refiner_test.cc
@@ -28,8 +28,9 @@ // // Author: sameeragarwal@google.com (Sameer Agarwal) -#include "Eigen/Dense" #include "ceres/iterative_refiner.h" + +#include "Eigen/Dense" #include "ceres/internal/eigen.h" #include "ceres/sparse_cholesky.h" #include "ceres/sparse_matrix.h" @@ -58,7 +59,6 @@ // y += Ax virtual void RightMultiply(const double* x, double* y) const { VectorRef(y, m_.cols()) += m_ * ConstVectorRef(x, m_.cols()); - } // y += A'x virtual void LeftMultiply(const double* x, double* y) const { @@ -70,7 +70,7 @@ virtual const double* values() const { return m_.data(); } virtual int num_rows() const { return m_.cols(); } virtual int num_cols() const { return m_.cols(); } - virtual int num_nonzeros() const {return m_.cols() * m_.cols(); } + virtual int num_nonzeros() const { return m_.cols() * m_.cols(); } // The following methods are not needed for tests in this file. virtual void SquaredColumnNorm(double* x) const DO_NOT_CALL; @@ -139,37 +139,19 @@ int num_cols_; int max_num_iterations_; Matrix lhs_; - Vector rhs_; - Vector solution_; + Vector rhs_, solution_; }; -TEST_F(IterativeRefinerTest, - ExactSolutionWithExactFactorizationReturnsInZeroIterations) { - FakeSparseMatrix lhs(lhs_); - FakeSparseCholesky<double> sparse_cholesky(lhs_); - IterativeRefiner refiner(max_num_iterations_); - Vector refined_solution = solution_; - auto summary = refiner.Refine( - lhs, rhs_.data(), &sparse_cholesky, refined_solution.data()); - EXPECT_EQ(summary.num_iterations, 0); - EXPECT_TRUE(summary.converged); - EXPECT_NEAR( - (refined_solution - solution_).norm() / solution_.norm(), 0.0, 5e-15); -} - -TEST_F(IterativeRefinerTest, - RandomSolutionWithExactFactorizationReturnsInOneIteration) { +TEST_F(IterativeRefinerTest, RandomSolutionWithExactFactorizationConverges) { FakeSparseMatrix lhs(lhs_); FakeSparseCholesky<double> sparse_cholesky(lhs_); IterativeRefiner refiner(max_num_iterations_); Vector refined_solution(num_cols_); refined_solution.setRandom(); - auto summary = refiner.Refine( - lhs, rhs_.data(), &sparse_cholesky, refined_solution.data()); - EXPECT_EQ(summary.num_iterations, 1); - EXPECT_TRUE(summary.converged); - EXPECT_NEAR( - (refined_solution - solution_).norm() / solution_.norm(), 0.0, 5e-15); + refiner.Refine(lhs, rhs_.data(), &sparse_cholesky, refined_solution.data()); + EXPECT_NEAR((lhs_ * refined_solution - rhs_).norm(), + 0.0, + std::numeric_limits<double>::epsilon()); } TEST_F(IterativeRefinerTest, @@ -181,11 +163,10 @@ IterativeRefiner refiner(max_num_iterations_); Vector refined_solution(num_cols_); refined_solution.setRandom(); - auto summary = refiner.Refine( - lhs, rhs_.data(), &sparse_cholesky, refined_solution.data()); - EXPECT_TRUE(summary.converged); - EXPECT_NEAR( - (refined_solution - solution_).norm() / solution_.norm(), 0.0, 5e-14); + refiner.Refine(lhs, rhs_.data(), &sparse_cholesky, refined_solution.data()); + EXPECT_NEAR((lhs_ * refined_solution - rhs_).norm(), + 0.0, + std::numeric_limits<double>::epsilon()); } } // namespace internal
diff --git a/internal/ceres/sparse_cholesky_test.cc b/internal/ceres/sparse_cholesky_test.cc index 1c79bc9..34098f3 100644 --- a/internal/ceres/sparse_cholesky_test.cc +++ b/internal/ceres/sparse_cholesky_test.cc
@@ -236,10 +236,10 @@ public: MockIterativeRefiner() : IterativeRefiner(1) {} MOCK_METHOD4(Refine, - Summary(const SparseMatrix& lhs, - const double* rhs, - SparseCholesky* sparse_cholesky, - double* solution)); + void (const SparseMatrix& lhs, + const double* rhs, + SparseCholesky* sparse_cholesky, + double* solution)); };