mirror of
https://github.com/ceres-solver/ceres-solver.git
synced 2026-08-29 08:34:37 +08:00
Add Iterative Refinement
Add a class IterativeRefiner which implements iterative refinement for SPD linear systems. Change-Id: I705d4e96cb7de9226ee35e2a9c11d98ffc0ee239
This commit is contained in:
Sameer Agarwal
@@ -112,6 +112,7 @@ CERES_TESTS = [
|
||||
"inner_product_computer",
|
||||
"invert_psd_matrix",
|
||||
"is_close",
|
||||
"iterative_refiner",
|
||||
"iterative_schur_complement_solver",
|
||||
"jet",
|
||||
"levenberg_marquardt_strategy",
|
||||
|
||||
@@ -74,6 +74,7 @@ CERES_SRCS = ["internal/ceres/" + filename for filename in [
|
||||
"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",
|
||||
|
||||
@@ -75,6 +75,7 @@ set(CERES_INTERNAL_SRC
|
||||
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 @@ if (BUILD_TESTING AND GFLAGS)
|
||||
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)
|
||||
|
||||
@@ -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
|
||||
@@ -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_
|
||||
@@ -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
|
||||
Reference in New Issue
Block a user