conjugate_gradient.h

Go to the documentation of this file.

#ifndef SOPT_CONJUGATE_GRADIENT
#define SOPT_CONJUGATE_GRADIENT
 
#include "sopt/config.h"
#include <limits>
#include <type_traits>
#include "sopt/logging.h"
#include "sopt/maths.h"
#include "sopt/types.h"
#include "sopt/wrapper.h"
 
namespace sopt {
class ConjugateGradient {
  template <typename T>
  class ApplyMatrix;
 
 public:
  struct Diagnostic {
    t_uint niters;
    t_real residual;
    bool good;
  };
  template <typename T>
  struct DiagnosticAndResult : public Diagnostic {
    Vector<T> result;
  };
  ConjugateGradient(t_uint itermax = std::numeric_limits<t_uint>::max(), t_real tolerance = 1e-8)
      : tolerance_(tolerance), itermax_(itermax) {}
  virtual ~ConjugateGradient() {}
 
  template <typename VECTOR, typename T1, typename T2>
  Diagnostic operator()(VECTOR &x, Eigen::MatrixBase<T1> const &A,
                        Eigen::MatrixBase<T2> const &b) const {
    return implementation(x, A, b);
  }
  template <typename VECTOR, typename T1, typename T2>
  typename std::enable_if<not std::is_base_of<Eigen::EigenBase<T1>, T1>::value, Diagnostic>::type
  operator()(VECTOR &x, T1 const &A, Eigen::MatrixBase<T2> const &b) const {
    return implementation(x, details::wrap<typename VECTOR::PlainObject>(A), b);
  }
  template <typename T0, typename A_TYPE>
  DiagnosticAndResult<typename T0::Scalar> operator()(A_TYPE const &A,
                                                      Eigen::MatrixBase<T0> const &b) const {
    DiagnosticAndResult<typename T0::Scalar> result;
    result.result = Vector<typename T0::Scalar>::Zero(b.size());
    *static_cast<Diagnostic *>(&result) = operator()(result.result, A, b);
    return result;
  }
 
  t_uint itermax() const { return itermax_; }
  void itermax(t_uint const &itermax) { itermax_ = itermax; }
  t_real tolerance() const { return tolerance_; }
  void tolerance(t_real const &tolerance) {
    if (tolerance <= 0e0) throw std::domain_error("Incorrect tolerance input");
    tolerance_ = tolerance;
  }
 
 protected:
  t_real tolerance_;
  t_uint itermax_;
  Image<> work_v;
  Image<> work_r;
  Image<> work_p;
 
 private:
  template <typename VECTOR, typename T1, typename MATRIXLIKE>
  Diagnostic implementation(VECTOR &x, MATRIXLIKE const &A, Eigen::MatrixBase<T1> const &b) const;
};
 
template <typename VECTOR, typename T1, typename MATRIXLIKE>
ConjugateGradient::Diagnostic ConjugateGradient::implementation(
    VECTOR &x, MATRIXLIKE const &A, Eigen::MatrixBase<T1> const &b) const {
  using Scalar = typename T1::Scalar;
  using Real = typename real_type<Scalar>::type;
 
  x.resize(b.size());
  if (std::abs((b.transpose().conjugate() * b)(0)) < tolerance()) {
    x.fill(0);
    return {0, 0, true};
  }
 
  Vector<Scalar> Ap(b.size());
 
  x = b;
  Vector<Scalar> residuals = b - A * x;
  Vector<Scalar> p = residuals;
  Real residual = std::abs((residuals.transpose().conjugate() * residuals)(0));
 
  t_uint i(0);
  for (; i < itermax(); ++i) {
    Ap = A * p;
    Scalar const alpha = residual / (p.transpose().conjugate() * Ap)(0);
    x += alpha * p;
    residuals -= alpha * Ap;
 
    Real new_residual = std::abs((residuals.transpose().conjugate() * residuals)(0));
    SOPT_LOW_LOG("CG iteration {} - residuals: {}", i, new_residual);
    if (std::abs(new_residual) < tolerance()) {
      residual = new_residual;
      break;
    }
 
    p = residuals + new_residual / residual * p;
    residual = new_residual;
  }
  return {i, residual, residual < tolerance()};
}
}  // namespace sopt
#endif