Eigen::MatrixPowerAtomic< MatrixType > Class Template Reference

Class for computing matrix powers. More...

#include <MatrixPower.h>

Inheritance diagram for Eigen::MatrixPowerAtomic< MatrixType >:

Public Member Functions
	MatrixPowerAtomic (const MatrixType &T, RealScalar p)
	Constructor. More...
void	compute (ResultType &res) const
	Compute the matrix power. More...

Private Types
enum	{ RowsAtCompileTime = MatrixType::RowsAtCompileTime , MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime }
typedef MatrixType::Scalar	Scalar
typedef MatrixType::RealScalar	RealScalar
typedef std::complex< RealScalar >	ComplexScalar
typedef Block< MatrixType, Dynamic, Dynamic >	ResultType

Private Member Functions
void	computePade (int degree, const MatrixType &IminusT, ResultType &res) const
void	compute2x2 (ResultType &res, RealScalar p) const
void	computeBig (ResultType &res) const
Private Member Functions inherited from Eigen::internal::noncopyable
EIGEN_DEVICE_FUNC	noncopyable ()
EIGEN_DEVICE_FUNC	~noncopyable ()

Static Private Member Functions
static int	getPadeDegree (float normIminusT)
static int	getPadeDegree (double normIminusT)
static int	getPadeDegree (long double normIminusT)
static ComplexScalar	computeSuperDiag (const ComplexScalar &, const ComplexScalar &, RealScalar p)
static RealScalar	computeSuperDiag (RealScalar, RealScalar, RealScalar p)

Private Attributes
const MatrixType &	m_A
RealScalar	m_p

Detailed Description

template<typename MatrixType>
class Eigen::MatrixPowerAtomic< MatrixType >

Class for computing matrix powers.

Template Parameters

MatrixType

type of the base, expected to be an instantiation of the Matrix class template.

This class is capable of computing triangular real/complex matrices raised to a power in the interval \( (-1, 1) \).

Note

Currently this class is only used by MatrixPower. One may insist that this be nested into MatrixPower. This class is here to facilitate future development of triangular matrix functions.

Member Typedef Documentation

ComplexScalar

template<typename MatrixType >

typedef std::complex<RealScalar> Eigen::MatrixPowerAtomic< MatrixType >::ComplexScalar

private

RealScalar

template<typename MatrixType >

typedef MatrixType::RealScalar Eigen::MatrixPowerAtomic< MatrixType >::RealScalar

private

ResultType

template<typename MatrixType >

typedef Block<MatrixType, Dynamic, Dynamic> Eigen::MatrixPowerAtomic< MatrixType >::ResultType

private

Scalar

template<typename MatrixType >

typedef MatrixType::Scalar Eigen::MatrixPowerAtomic< MatrixType >::Scalar

private

Member Enumeration Documentation

anonymous enum

template<typename MatrixType >

anonymous enum

private

Enumerator
RowsAtCompileTime
MaxRowsAtCompileTime

{ RowsAtCompileTime = MatrixType::RowsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime };

Constructor & Destructor Documentation

MatrixPowerAtomic()

template<typename MatrixType >

Eigen::MatrixPowerAtomic< MatrixType >::MatrixPowerAtomic	(	const MatrixType &	T,
		RealScalar	p
	)

Constructor.

Parameters

[in]	T	the base of the matrix power.
[in]	p	the exponent of the matrix power, should be in \( (-1, 1) \).

The class stores a reference to T, so it should not be changed (or destroyed) before evaluation. Only the upper triangular part of T is read.

                                                                                  : m_A(T), m_p(p) {
  eigen_assert(T.rows() == T.cols());
  eigen_assert(p > -1 && p < 1);
}

References eigen_assert, and p.

Member Function Documentation

compute()

template<typename MatrixType >

void Eigen::MatrixPowerAtomic< MatrixType >::compute

(

ResultType &

res

)

const

Compute the matrix power.

Parameters

[out]

res

\( A^p \) where A and p are specified in the constructor.

                                                                 {
  using std::pow;
  switch (m_A.rows()) {
    case 0:
      break;
    case 1:
      res(0, 0) = pow(m_A(0, 0), m_p);
      break;
    case 2:
      compute2x2(res, m_p);
      break;
    default:
      computeBig(res);
  }
}

References Eigen::bfloat16_impl::pow(), and res.

Referenced by Eigen::MatrixPower< MatrixType >::computeFracPower().

compute2x2()

template<typename MatrixType >

void Eigen::MatrixPowerAtomic< MatrixType >::compute2x2 ( ResultType &  res, RealScalar  p  ) const

private

                                                                                  {
  using std::abs;
  using std::pow;
  res.coeffRef(0, 0) = pow(m_A.coeff(0, 0), p);
 
  for (Index i = 1; i < m_A.cols(); ++i) {
    res.coeffRef(i, i) = pow(m_A.coeff(i, i), p);
    if (m_A.coeff(i - 1, i - 1) == m_A.coeff(i, i))
      res.coeffRef(i - 1, i) = p * pow(m_A.coeff(i, i), p - 1);
    else if (2 * abs(m_A.coeff(i - 1, i - 1)) < abs(m_A.coeff(i, i)) ||
             2 * abs(m_A.coeff(i, i)) < abs(m_A.coeff(i - 1, i - 1)))
      res.coeffRef(i - 1, i) =
          (res.coeff(i, i) - res.coeff(i - 1, i - 1)) / (m_A.coeff(i, i) - m_A.coeff(i - 1, i - 1));
    else
      res.coeffRef(i - 1, i) = computeSuperDiag(m_A.coeff(i, i), m_A.coeff(i - 1, i - 1), p);
    res.coeffRef(i - 1, i) *= m_A.coeff(i - 1, i);
  }
}

References abs(), i, p, Eigen::bfloat16_impl::pow(), and res.

computeBig()

template<typename MatrixType >

void Eigen::MatrixPowerAtomic< MatrixType >::computeBig ( ResultType &  res ) const

private

                                                                    {
  using std::ldexp;
  const int digits = std::numeric_limits<RealScalar>::digits;
  const RealScalar maxNormForPade =
      RealScalar(digits <= 24    ? 4.3386528e-1L                              // single precision
                 : digits <= 53  ? 2.789358995219730e-1L                      // double precision
                 : digits <= 64  ? 2.4471944416607995472e-1L                  // extended precision
                 : digits <= 106 ? 1.1016843812851143391275867258512e-1L      // double-double
                                 : 9.134603732914548552537150753385375e-2L);  // quadruple precision
  MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
  RealScalar normIminusT;
  int degree, degree2, numberOfSquareRoots = 0;
  bool hasExtraSquareRoot = false;
 
  for (Index i = 0; i < m_A.cols(); ++i) eigen_assert(m_A(i, i) != RealScalar(0));
 
  while (true) {
    IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
    normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
    if (normIminusT < maxNormForPade) {
      degree = getPadeDegree(normIminusT);
      degree2 = getPadeDegree(normIminusT / 2);
      if (degree - degree2 <= 1 || hasExtraSquareRoot) break;
      hasExtraSquareRoot = true;
    }
    matrix_sqrt_triangular(T, sqrtT);
    T = sqrtT.template triangularView<Upper>();
    ++numberOfSquareRoots;
  }
  computePade(degree, IminusT, res);
 
  for (; numberOfSquareRoots; --numberOfSquareRoots) {
    compute2x2(res, ldexp(m_p, -numberOfSquareRoots));
    res = res.template triangularView<Upper>() * res;
  }
  compute2x2(res, m_p);
}

References constants::degree, eigen_assert, i, Eigen::matrix_sqrt_triangular(), and res.

computePade()

template<typename MatrixType >

void Eigen::MatrixPowerAtomic< MatrixType >::computePade ( int  degree, const MatrixType &  IminusT, ResultType &  res  ) const

private

                                                                                                            {
  int i = 2 * degree;
  res = (m_p - RealScalar(degree)) / RealScalar(2 * i - 2) * IminusT;
 
  for (--i; i; --i) {
    res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res)
              .template triangularView<Upper>()
              .solve((i == 1  ? -m_p
                      : i & 1 ? (-m_p - RealScalar(i / 2)) / RealScalar(2 * i)
                              : (m_p - RealScalar(i / 2)) / RealScalar(2 * i - 2)) *
                     IminusT)
              .eval();
  }
  res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
}

References constants::degree, eval(), i, and res.

computeSuperDiag() [1/2]

template<typename MatrixType >

MatrixPowerAtomic< MatrixType >::ComplexScalar Eigen::MatrixPowerAtomic< MatrixType >::computeSuperDiag ( const ComplexScalar &  curr, const ComplexScalar &  prev, RealScalar  p  )

inlinestaticprivate

                                                                        {
  using std::ceil;
  using std::exp;
  using std::log;
  using std::sinh;
 
  ComplexScalar logCurr = log(curr);
  ComplexScalar logPrev = log(prev);
  RealScalar unwindingNumber =
      ceil((numext::imag(logCurr - logPrev) - RealScalar(EIGEN_PI)) / RealScalar(2 * EIGEN_PI));
  ComplexScalar w =
      numext::log1p((curr - prev) / prev) / RealScalar(2) + ComplexScalar(0, RealScalar(EIGEN_PI) * unwindingNumber);
  return RealScalar(2) * exp(RealScalar(0.5) * p * (logCurr + logPrev)) * sinh(p * w) / (curr - prev);
}

References Eigen::bfloat16_impl::ceil(), EIGEN_PI, Eigen::bfloat16_impl::exp(), imag(), Eigen::bfloat16_impl::log(), Eigen::bfloat16_impl::log1p(), p, Eigen::bfloat16_impl::sinh(), and w.

computeSuperDiag() [2/2]

template<typename MatrixType >

MatrixPowerAtomic< MatrixType >::RealScalar Eigen::MatrixPowerAtomic< MatrixType >::computeSuperDiag ( RealScalar  curr, RealScalar  prev, RealScalar  p  )

inlinestaticprivate

                                                    {
  using std::exp;
  using std::log;
  using std::sinh;
 
  RealScalar w = numext::log1p((curr - prev) / prev) / RealScalar(2);
  return 2 * exp(p * (log(curr) + log(prev)) / 2) * sinh(p * w) / (curr - prev);
}

References Eigen::bfloat16_impl::exp(), Eigen::bfloat16_impl::log(), Eigen::bfloat16_impl::log1p(), p, Eigen::bfloat16_impl::sinh(), and w.

getPadeDegree() [1/3]

template<typename MatrixType >

int Eigen::MatrixPowerAtomic< MatrixType >::getPadeDegree ( double  normIminusT )

inlinestaticprivate

                                                                          {
  const double maxNormForPade[] = {1.884160592658218e-2 /* degree = 3 */, 6.038881904059573e-2, 1.239917516308172e-1,
                                   1.999045567181744e-1, 2.789358995219730e-1};
  int degree = 3;
  for (; degree <= 7; ++degree)
    if (normIminusT <= maxNormForPade[degree - 3]) break;
  return degree;
}

References constants::degree.

getPadeDegree() [2/3]

template<typename MatrixType >

int Eigen::MatrixPowerAtomic< MatrixType >::getPadeDegree ( float  normIminusT )

inlinestaticprivate

                                                                         {
  const float maxNormForPade[] = {2.8064004e-1f /* degree = 3 */, 4.3386528e-1f};
  int degree = 3;
  for (; degree <= 4; ++degree)
    if (normIminusT <= maxNormForPade[degree - 3]) break;
  return degree;
}

References constants::degree.

getPadeDegree() [3/3]

template<typename MatrixType >

int Eigen::MatrixPowerAtomic< MatrixType >::getPadeDegree ( long double  normIminusT )

inlinestaticprivate

                                                                               {
#if LDBL_MANT_DIG == 53
  const int maxPadeDegree = 7;
  const double maxNormForPade[] = {1.884160592658218e-2L /* degree = 3 */, 6.038881904059573e-2L, 1.239917516308172e-1L,
                                   1.999045567181744e-1L, 2.789358995219730e-1L};
#elif LDBL_MANT_DIG <= 64
  const int maxPadeDegree = 8;
  const long double maxNormForPade[] = {6.3854693117491799460e-3L /* degree = 3 */,
                                        2.6394893435456973676e-2L,
                                        6.4216043030404063729e-2L,
                                        1.1701165502926694307e-1L,
                                        1.7904284231268670284e-1L,
                                        2.4471944416607995472e-1L};
#elif LDBL_MANT_DIG <= 106
  const int maxPadeDegree = 10;
  const double maxNormForPade[] = {1.0007161601787493236741409687186e-4L /* degree = 3 */,
                                   1.0007161601787493236741409687186e-3L,
                                   4.7069769360887572939882574746264e-3L,
                                   1.3220386624169159689406653101695e-2L,
                                   2.8063482381631737920612944054906e-2L,
                                   4.9625993951953473052385361085058e-2L,
                                   7.7367040706027886224557538328171e-2L,
                                   1.1016843812851143391275867258512e-1L};
#else
  const int maxPadeDegree = 10;
  const double maxNormForPade[] = {5.524506147036624377378713555116378e-5L /* degree = 3 */,
                                   6.640600568157479679823602193345995e-4L,
                                   3.227716520106894279249709728084626e-3L,
                                   9.619593944683432960546978734646284e-3L,
                                   2.134595382433742403911124458161147e-2L,
                                   3.908166513900489428442993794761185e-2L,
                                   6.266780814639442865832535460550138e-2L,
                                   9.134603732914548552537150753385375e-2L};
#endif
  int degree = 3;
  for (; degree <= maxPadeDegree; ++degree)
    if (normIminusT <= static_cast<long double>(maxNormForPade[degree - 3])) break;
  return degree;
}

References constants::degree.

Member Data Documentation

m_A

template<typename MatrixType >

const MatrixType& Eigen::MatrixPowerAtomic< MatrixType >::m_A

private

m_p

template<typename MatrixType >

RealScalar Eigen::MatrixPowerAtomic< MatrixType >::m_p

private

---

The documentation for this class was generated from the following file:

• MatrixPower.h