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));
};
