truncatedconjugategradient.hh

Go to the documentation of this file.

// Original File: https://gitlab.com/libeigen/eigen/-/blob/master/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
// SPDX-FileCopyrightText: 2011-2014 Copyright (C) Gael Guennebaud <gael.guennebaud@inria.fr>
// SPDX-License-Identifier: MPL-2.0
// Modifications:
// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Sparse>
 
namespace Eigen {
  enum class TCGStopReason : int {
    negativeCurvature,
    exceededTrustRegion,
    reachedTargetResidualKappaLinear,
    reachedTargetResidualThetaSuperLinear,
    maximumInnerIterations,
    modelIncreased
  };
  template <typename Scalar>
  struct TCGInfo {
    TCGStopReason stop_tCG    = TCGStopReason::maximumInnerIterations;
    Scalar Delta              = 100000;
    Scalar kappa              = 0.1;
    Scalar theta              = 1.0;
    Eigen::Index mininner     = 1;
    Eigen::Index maxinner     = 1000;
    Eigen::Index numInnerIter = 0;
 
    void initRuntimeOptions(int _num_dof_solve) {
      maxinner = _num_dof_solve * 2;
      Delta    = 100000;  // typical distance of manifold
      //            Delta0 = Delta_bar/8.0;
    }
  };
  namespace internal {
 
    template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
    void truncated_conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, const Preconditioner& precond,
                                      Eigen::Index& iters, typename Dest::RealScalar& tol_error,
                                      TCGInfo<typename Dest::RealScalar>& _info) {
      using std::abs;
      using std::sqrt;
      typedef typename Dest::RealScalar RealScalar;
      typedef typename Dest::Scalar Scalar;
      typedef Matrix<Scalar, Dynamic, 1> VectorType;
 
      RealScalar tol = tol_error;
      Index maxIters = iters;
 
      Index n = mat.cols();
 
      VectorType residual = rhs - mat * x;  // initial residual
 
      RealScalar rhsNorm2             = rhs.norm();
      const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
 
      if (rhsNorm2 <= considerAsZero) {
        x.setZero();
        iters     = 0;
        tol_error = 0;
        return;
      }
      RealScalar threshold     = numext::maxi(tol * tol * rhsNorm2 * rhsNorm2, considerAsZero);
      RealScalar residualNorm2 = residual.norm();
      if (residualNorm2 * residualNorm2 < threshold) {
        iters     = 0;
        tol_error = (residualNorm2 / rhsNorm2);
        return;
      }
 
      double e_Pd     = 0.0;
      double e_Pe_new = 0.0;
      double e_Pe     = x.squaredNorm();
      double d_Pd;
      double d_Hd;
      VectorType p(n);
      p = precond.solve(residual);  // initial search direction
                                    // bool coutflag=true;
      _info.stop_tCG = TCGStopReason::maximumInnerIterations;
      VectorType z(n), tmp(n);
      RealScalar absNew = numext::real(residual.dot(p));  // the square of the absolute value of r scaled by invM
      d_Pd              = absNew;
      Index i           = 1;
      while (i < maxIters) {
        tmp.noalias() = mat * p;  // the bottleneck of the algorithm
        d_Hd          = p.dot(tmp);
        Scalar alpha  = absNew / d_Hd;  // the amount we travel on dir
 
        e_Pe_new = e_Pe + 2.0 * alpha * e_Pd + alpha * alpha * d_Pd;
 
        if (d_Hd <= 0 || e_Pe_new >= _info.Delta * _info.Delta)  // negative curvature or execdet trustregion
        {
          double tau = (-e_Pd + sqrt(e_Pd * e_Pd + d_Pd * (_info.Delta * _info.Delta - e_Pe))) / d_Pd;
 
          x += tau * p;
          if (d_Hd <= 0)
            _info.stop_tCG = TCGStopReason::negativeCurvature;
          else
            _info.stop_tCG = TCGStopReason::exceededTrustRegion;
 
          break;
        }
        e_Pe = e_Pe_new;
        x += alpha * p;           // update solution
        residual -= alpha * tmp;  // update residual
 
        residualNorm2 = residual.norm();
 
        if (i >= _info.mininner
            && residualNorm2 <= rhsNorm2 * std::min(rhsNorm2, _info.kappa))  // missing pow(rhsNorm2,_info.theta
        {
          // Residual is small enough to quit
          if (_info.kappa < rhsNorm2)
            _info.stop_tCG = TCGStopReason::reachedTargetResidualKappaLinear;
          else
            _info.stop_tCG = TCGStopReason::reachedTargetResidualThetaSuperLinear;
          break;
        }
        if (residualNorm2 < threshold) break;
 
        z = precond.solve(residual);  // approximately solve for "A z = residual"
 
        RealScalar absOld = absNew;
        absNew            = numext::real(residual.dot(z));  // update the absolute value of r
        RealScalar beta = absNew / absOld;  // calculate the Gram-Schmidt value used to create the new search direction
 
        e_Pd = beta * (e_Pd + alpha * d_Pd);
        d_Pd = absNew + beta * beta * d_Pd;
 
        p = z + beta * p;  // update search direction
        i++;
      }
      tol_error          = (residualNorm2 / rhsNorm2);
      iters              = i;
      _info.numInnerIter = i;
    }
 
  }  // namespace internal
 
  template <typename MatrixType, int UpLo = Lower,
            typename Preconditioner = DiagonalPreconditioner<typename MatrixType::Scalar> >
  class TruncatedConjugateGradient;
 
  namespace internal {
 
    template <typename MatrixType_, int UpLo, typename Preconditioner_>
    struct traits<TruncatedConjugateGradient<MatrixType_, UpLo, Preconditioner_> > {
      typedef MatrixType_ MatrixType;
      typedef Preconditioner_ Preconditioner;
    };
 
  }  // namespace internal
 
  template <typename MatrixType_, int UpLo_, typename Preconditioner_>
  class TruncatedConjugateGradient
      : public IterativeSolverBase<TruncatedConjugateGradient<MatrixType_, UpLo_, Preconditioner_> > {
  public:
    typedef IterativeSolverBase<TruncatedConjugateGradient> Base;
    TruncatedConjugateGradient(TruncatedConjugateGradient&& other) noexcept
        : Base(std::move(other)), algInfo{other.algInfo} {}
 
  private:
    using Base::m_error;
    using Base::m_info;
    using Base::m_isInitialized;
    using Base::m_iterations;
    using Base::matrix;
    mutable TCGInfo<typename MatrixType_::RealScalar> algInfo;
 
  public:
    typedef MatrixType_ MatrixType;
    typedef typename MatrixType::Scalar Scalar;
    typedef typename MatrixType::RealScalar RealScalar;
    typedef Preconditioner_ Preconditioner;
 
    enum { UpLo = UpLo_ };
 
  public:
    TCGInfo<typename MatrixType::RealScalar> getInfo() { return algInfo; }
 
    void setInfo(TCGInfo<typename MatrixType::RealScalar> _alginfo) { this->algInfo = _alginfo; }
    TruncatedConjugateGradient() : Base() {}
 
    template <typename MatrixDerived>
    explicit TruncatedConjugateGradient(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
 
    ~TruncatedConjugateGradient() {}
 
    template <typename Rhs, typename Dest>
    void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const {
      typedef typename Base::MatrixWrapper MatrixWrapper;
      typedef typename Base::ActualMatrixType ActualMatrixType;
      enum {
        TransposeInput = (!MatrixWrapper::MatrixFree) && (UpLo == (Lower | Upper)) && (!MatrixType::IsRowMajor)
                         && (!NumTraits<Scalar>::IsComplex)
      };
      typedef std::conditional_t<TransposeInput, Transpose<const ActualMatrixType>, ActualMatrixType const&>
          RowMajorWrapper;
      EIGEN_STATIC_ASSERT(internal::check_implication(MatrixWrapper::MatrixFree, UpLo == (Lower | Upper)),
                          MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
      typedef std::conditional_t<UpLo == (Lower | Upper), RowMajorWrapper,
                                 typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type>
          SelfAdjointWrapper;
 
      m_iterations = Base::maxIterations();
      m_error      = Base::m_tolerance;
 
      RowMajorWrapper row_mat(matrix());
      internal::truncated_conjugate_gradient(SelfAdjointWrapper(row_mat), b, x, Base::m_preconditioner, m_iterations,
                                             m_error, algInfo);
      m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
    }
 
  protected:
  };
 
}  // end namespace Eigen