Add Iterative Refinement Add a class IterativeRefiner which implements iterative refinement for SPD linear systems. Change-Id: I705d4e96cb7de9226ee35e2a9c11d98ffc0ee239
diff --git a/BUILD b/BUILD index 15de209..2d0a5f0 100644 --- a/BUILD +++ b/BUILD
@@ -112,6 +112,7 @@ "inner_product_computer", "invert_psd_matrix", "is_close", + "iterative_refiner", "iterative_schur_complement_solver", "jet", "levenberg_marquardt_strategy",
diff --git a/bazel/ceres.bzl b/bazel/ceres.bzl index d90e5a3..6ba0137 100644 --- a/bazel/ceres.bzl +++ b/bazel/ceres.bzl
@@ -74,6 +74,7 @@ "is_close.cc", "implicit_schur_complement.cc", "inner_product_computer.cc", + "iterative_refiner.cc", "iterative_schur_complement_solver.cc", "lapack.cc", "levenberg_marquardt_strategy.cc",
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index a70f973..8924173 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt
@@ -75,6 +75,7 @@ implicit_schur_complement.cc inner_product_computer.cc is_close.cc + iterative_refiner.cc iterative_schur_complement_solver.cc levenberg_marquardt_strategy.cc lapack.cc @@ -338,6 +339,7 @@ ceres_test(inner_product_computer) ceres_test(invert_psd_matrix) ceres_test(is_close) + ceres_test(iterative_refiner) ceres_test(iterative_schur_complement_solver) ceres_test(jet) ceres_test(levenberg_marquardt_strategy)
diff --git a/internal/ceres/iterative_refiner.cc b/internal/ceres/iterative_refiner.cc new file mode 100644 index 0000000..6a5a0c7 --- /dev/null +++ b/internal/ceres/iterative_refiner.cc
@@ -0,0 +1,112 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2018 Google Inc. All rights reserved. +// http://ceres-solver.org/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: sameeragarwal@google.com (Sameer Agarwal) + +#include <string> +#include "ceres/iterative_refiner.h" + +#include "Eigen/Core" +#include "ceres/sparse_cholesky.h" +#include "ceres/sparse_matrix.h" + +namespace ceres { +namespace internal { + +IterativeRefiner::IterativeRefiner(const int num_cols, + const int max_num_iterations) + : num_cols_(num_cols), + max_num_iterations_(max_num_iterations), + residual_(num_cols), + correction_(num_cols), + lhs_x_solution_(num_cols) {} + +IterativeRefiner::Summary IterativeRefiner::Refine( + const SparseMatrix& lhs, + const double* rhs_ptr, + SparseCholesky* sparse_cholesky, + double* solution_ptr) { + Summary summary; + + 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>(); + + // 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>(); + } + + return summary; +}; + +} // namespace internal +} // namespace ceres
diff --git a/internal/ceres/iterative_refiner.h b/internal/ceres/iterative_refiner.h new file mode 100644 index 0000000..471116c --- /dev/null +++ b/internal/ceres/iterative_refiner.h
@@ -0,0 +1,111 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2018 Google Inc. All rights reserved. +// http://ceres-solver.org/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: sameeragarwal@google.com (Sameer Agarwal) + +#ifndef CERES_INTERNAL_ITERATIVE_REFINER_H_ +#define CERES_INTERNAL_ITERATIVE_REFINER_H_ + +// This include must come before any #ifndef check on Ceres compile options. +#include "ceres/internal/port.h" +#include "ceres/internal/eigen.h" + +namespace ceres { +namespace internal { + +class SparseMatrix; +class SparseCholesky; + +// Iterative refinement +// (https://en.wikipedia.org/wiki/Iterative_refinement) is the process +// of improving the solution to a linear system, by using the +// following iteration. +// +// r_i = b - Ax_i +// Ad_i = r_i +// x_{i+1} = x_i + d_i +// +// 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 +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; + }; + + // num_cols is the number of rows & columns in the linear system + // being solved. + // + // max_num_iterations is the maximum number of refinement iterations + // to perform. + IterativeRefiner(int num_cols, int max_num_iterations); + + // Given an initial estimate of the solution of lhs * x = rhs, use + // iterative refinement to improve it. + // + // sparse_cholesky is assumed to contain an already computed + // factorization (or approximation thereof) of lhs. + // + // solution is expected to contain a approximation to the solution + // to lhs * x = rhs. It can be zero. + Summary Refine(const SparseMatrix& lhs, + const double* rhs, + SparseCholesky* sparse_cholesky, + double* solution); + + private: + int num_cols_; + int max_num_iterations_; + Vector residual_; + Vector correction_; + Vector lhs_x_solution_; +}; + +} // namespace internal +} // namespace ceres + +#endif // CERES_INTERNAL_ITERATIVE_REFINER_H_
diff --git a/internal/ceres/iterative_refiner_test.cc b/internal/ceres/iterative_refiner_test.cc new file mode 100644 index 0000000..4aac253 --- /dev/null +++ b/internal/ceres/iterative_refiner_test.cc
@@ -0,0 +1,192 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2018 Google Inc. All rights reserved. +// http://ceres-solver.org/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: sameeragarwal@google.com (Sameer Agarwal) + +#include "Eigen/Dense" +#include "ceres/iterative_refiner.h" +#include "ceres/internal/eigen.h" +#include "ceres/sparse_cholesky.h" +#include "ceres/sparse_matrix.h" +#include "glog/logging.h" +#include "gtest/gtest.h" + +namespace ceres { +namespace internal { + +// Macros to help us define virtual methods which we do not expect to +// use/call in this test. +#define DO_NOT_CALL \ + { LOG(FATAL) << "DO NOT CALL"; } +#define DO_NOT_CALL_WITH_RETURN(x) \ + { \ + LOG(FATAL) << "DO NOT CALL"; \ + return x; \ + } + +// A fake SparseMatrix, which uses an Eigen matrix to do the real work. +class FakeSparseMatrix : public SparseMatrix { + public: + FakeSparseMatrix(const Matrix& m) : m_(m) {} + virtual ~FakeSparseMatrix() {} + + // 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 { + // We will assume that this is a symmetric matrix. + RightMultiply(x, y); + } + + virtual double* mutable_values() { return m_.data(); } + 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(); } + + // The following methods are not needed for tests in this file. + virtual void SquaredColumnNorm(double* x) const DO_NOT_CALL; + virtual void ScaleColumns(const double* scale) DO_NOT_CALL; + virtual void SetZero() DO_NOT_CALL; + virtual void ToDenseMatrix(Matrix* dense_matrix) const DO_NOT_CALL; + virtual void ToTextFile(FILE* file) const DO_NOT_CALL; + + private: + Matrix m_; +}; + +// A fake SparseCholesky which uses Eigen's Cholesky factorization to +// do the real work. The template parameter allows us to work in +// doubles or floats, even though the source matrix is double. +template <typename Scalar> +class FakeSparseCholesky : public SparseCholesky { + public: + FakeSparseCholesky(const Matrix& lhs) { lhs_ = lhs.cast<Scalar>(); } + virtual ~FakeSparseCholesky() {} + + virtual LinearSolverTerminationType Solve(const double* rhs_ptr, + double* solution_ptr, + std::string* message) { + const int num_cols = lhs_.cols(); + VectorRef solution(solution_ptr, num_cols); + ConstVectorRef rhs(rhs_ptr, num_cols); + solution = lhs_.llt().solve(rhs.cast<Scalar>()).template cast<double>(); + return LINEAR_SOLVER_SUCCESS; + } + + // The following methods are not needed for tests in this file. + virtual CompressedRowSparseMatrix::StorageType StorageType() const + DO_NOT_CALL_WITH_RETURN(CompressedRowSparseMatrix::UPPER_TRIANGULAR); + virtual LinearSolverTerminationType Factorize(CompressedRowSparseMatrix* lhs, + std::string* message) + DO_NOT_CALL_WITH_RETURN(LINEAR_SOLVER_FAILURE); + + virtual LinearSolverTerminationType FactorAndSolve( + CompressedRowSparseMatrix* lhs, + const double* rhs, + double* solution, + std::string* message) DO_NOT_CALL_WITH_RETURN(LINEAR_SOLVER_FAILURE); + + private: + Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> lhs_; +}; + +#undef DO_NOT_CALL +#undef DO_NOT_CALL_WITH_RETURN + +class IterativeRefinerTest : public ::testing::Test { + public: + void SetUp() { + num_cols_ = 5; + max_num_iterations_ = 30; + Matrix m(num_cols_, num_cols_); + m.setRandom(); + lhs_ = m * m.transpose(); + solution_.resize(num_cols_); + solution_.setRandom(); + rhs_ = lhs_ * solution_; + }; + + protected: + int num_cols_; + int max_num_iterations_; + Matrix lhs_; + Vector rhs_; + Vector solution_; +}; + +TEST_F(IterativeRefinerTest, + ExactSolutionWithExactFactorizationReturnsInZeroIterations) { + FakeSparseMatrix lhs(lhs_); + FakeSparseCholesky<double> sparse_cholesky(lhs_); + IterativeRefiner refiner(num_cols_, 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) { + FakeSparseMatrix lhs(lhs_); + FakeSparseCholesky<double> sparse_cholesky(lhs_); + IterativeRefiner refiner(num_cols_, 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); +} + +TEST_F(IterativeRefinerTest, + RandomSolutionWithApproximationFactorizationConverges) { + FakeSparseMatrix lhs(lhs_); + // Use a single precision Cholesky factorization of the double + // precision matrix. This will give us an approximate factorization. + FakeSparseCholesky<float> sparse_cholesky(lhs_); + IterativeRefiner refiner(num_cols_, 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-15); +} + +} // namespace internal +} // namespace ceres