Public Member Functions | Private Types | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
Eigen::MatrixPower< MatrixType > Class Template Reference
Matrix functions module

Class for computing matrix powers. More...

Inherits internal::noncopyable.

Public Member Functions

Index cols () const
 
template<typename ResultType >
void compute (ResultType &res, RealScalar p)
 Compute the matrix power. More...
 
 MatrixPower (const MatrixType &A)
 Constructor. More...
 
const MatrixPowerParenthesesReturnValue< MatrixType > operator() (RealScalar p)
 Returns the matrix power. More...
 
Index rows () const
 

Private Types

typedef Matrix< ComplexScalar, Dynamic, Dynamic, 0, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime > ComplexMatrix
 
typedef std::complex< RealScalar > ComplexScalar
 
typedef MatrixType::RealScalar RealScalar
 
typedef MatrixType::Scalar Scalar
 

Private Member Functions

template<typename ResultType >
void computeFracPower (ResultType &res, RealScalar p)
 
template<typename ResultType >
void computeIntPower (ResultType &res, RealScalar p)
 
void initialize ()
 Perform Schur decomposition for fractional power. More...
 
void split (RealScalar &p, RealScalar &intpart)
 Split p into integral part and fractional part. More...
 

Static Private Member Functions

template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
static void revertSchur (Matrix< ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols > &res, const ComplexMatrix &T, const ComplexMatrix &U)
 
template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
static void revertSchur (Matrix< RealScalar, Rows, Cols, Options, MaxRows, MaxCols > &res, const ComplexMatrix &T, const ComplexMatrix &U)
 

Private Attributes

MatrixType::Nested m_A
 Reference to the base of matrix power. More...
 
RealScalar m_conditionNumber
 Condition number of m_A. More...
 
ComplexMatrix m_fT
 Store fractional power of m_T. More...
 
Index m_nulls
 Rank deficiency of m_A. More...
 
Index m_rank
 Rank of m_A. More...
 
ComplexMatrix m_T
 Store the result of Schur decomposition. More...
 
MatrixType m_tmp
 Temporary storage. More...
 
ComplexMatrix m_U
 

Detailed Description

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

Class for computing matrix powers.

Template Parameters
MatrixTypetype of the base, expected to be an instantiation of the Matrix class template.

This class is capable of computing real/complex matrices raised to an arbitrary real power. Meanwhile, it saves the result of Schur decomposition if an non-integral power has even been calculated. Therefore, if you want to compute multiple (>= 2) matrix powers for the same matrix, using the class directly is more efficient than calling MatrixBase::pow().

Example:

#include <unsupported/Eigen/MatrixFunctions>
#include <iostream>
using namespace Eigen;
int main()
{
Matrix4cd A = Matrix4cd::Random();
MatrixPower<Matrix4cd> Apow(A);
std::cout << "The matrix A is:\n" << A << "\n\n"
"A^3.1 is:\n" << Apow(3.1) << "\n\n"
"A^3.3 is:\n" << Apow(3.3) << "\n\n"
"A^3.7 is:\n" << Apow(3.7) << "\n\n"
"A^3.9 is:\n" << Apow(3.9) << std::endl;
return 0;
}
main
int main()
Definition: BVH_Example.cpp:25
A
SparseMatrix< double > A(n, n)
Eigen::DenseBase::Random
static const RandomReturnType Random()
Eigen
: TensorContractionSycl.h, provides various tensor contraction kernel for SYCL backend

Output: 

The matrix A is:
 (-0.211234,0.680375)   (0.10794,-0.444451)   (0.434594,0.271423) (-0.198111,-0.686642)
   (0.59688,0.566198) (0.257742,-0.0452059)  (0.213938,-0.716795) (-0.782382,-0.740419)
 (-0.604897,0.823295) (0.0268018,-0.270431) (-0.514226,-0.967399)  (-0.563486,0.997849)
 (0.536459,-0.329554)    (0.83239,0.904459)  (0.608354,-0.725537)  (0.678224,0.0258648)

A^3.1 is:
   (2.80575,-0.607662) (-1.16847,-0.00660555)    (-0.760385,1.01461)   (-0.38073,-0.106512)
     (1.4041,-3.61891)     (1.00481,0.186263)   (-0.163888,0.449419)   (-0.388981,-1.22629)
   (-2.07957,-1.58136)     (0.825866,2.25962)     (5.09383,0.155736)    (0.394308,-1.63034)
  (-0.818997,0.671026)  (2.11069,-0.00768024)    (-1.37876,0.140165)    (2.50512,-0.854429)

A^3.3 is:
  (2.83571,-0.238717) (-1.48174,-0.0615217)  (-0.0544396,1.68092) (-0.292699,-0.621726)
    (2.0521,-3.58316)    (0.87894,0.400548)  (0.738072,-0.121242)   (-1.07957,-1.63492)
  (-3.00106,-1.10558)     (1.52205,1.92407)    (5.29759,-1.83562)  (-0.532038,-1.50253)
  (-0.491353,-0.4145)     (2.5761,0.481286)  (-1.21994,0.0367069)    (2.67112,-1.06331)

A^3.7 is:
     (1.42126,0.33362)   (-1.39486,-0.560486)      (1.44968,2.47066)   (-0.324079,-1.75879)
    (2.65301,-1.82427)   (0.357333,-0.192429)      (2.01017,-1.4791)    (-2.71518,-2.35892)
   (-3.98544,0.964861)     (2.26033,0.554254)     (3.18211,-5.94352)    (-2.22888,0.128951)
   (0.944969,-2.14683)      (3.31345,1.66075) (-0.0623743,-0.848324)        (2.3897,-1.863)

A^3.9 is:
 (0.0720766,0.378685) (-0.931961,-0.978624)      (1.9855,2.34105)  (-0.530547,-2.17664)
  (2.40934,-0.265286)  (0.0299975,-1.08827)    (1.98974,-2.05886)   (-3.45767,-2.50235)
    (-3.71666,2.3874)        (2.054,-0.303)   (0.844348,-7.29588)    (-2.59136,1.57689)
   (1.87645,-2.38798)     (3.52111,2.10508)    (0.799055,-1.6122)    (1.93452,-2.44408)

Definition at line 339 of file MatrixPower.h.

Member Typedef Documentation

◆ ComplexMatrix

template<typename MatrixType >
typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime> Eigen::MatrixPower< MatrixType >::ComplexMatrix
private

Definition at line 387 of file MatrixPower.h.

◆ ComplexScalar

template<typename MatrixType >
typedef std::complex<RealScalar> Eigen::MatrixPower< MatrixType >::ComplexScalar
private

Definition at line 385 of file MatrixPower.h.

◆ RealScalar

template<typename MatrixType >
typedef MatrixType::RealScalar Eigen::MatrixPower< MatrixType >::RealScalar
private

Definition at line 343 of file MatrixPower.h.

◆ Scalar

template<typename MatrixType >
typedef MatrixType::Scalar Eigen::MatrixPower< MatrixType >::Scalar
private

Definition at line 342 of file MatrixPower.h.

Constructor & Destructor Documentation

◆ MatrixPower()

template<typename MatrixType >
Eigen::MatrixPower< MatrixType >::MatrixPower ( const MatrixType &  A)
inlineexplicit

Constructor.

Parameters
[in]Athe base of the matrix power.

The class stores a reference to A, so it should not be changed (or destroyed) before evaluation.

Definition at line 354 of file MatrixPower.h.

354  :
355  m_A(A),
356  m_conditionNumber(0),
357  m_rank(A.cols()),
358  m_nulls(0)
359  { eigen_assert(A.rows() == A.cols()); }
eigen_assert
#define eigen_assert(x)
Eigen::MatrixPower::m_rank
Index m_rank
Rank of m_A.
Definition: MatrixPower.h:410
Eigen::MatrixPower::m_nulls
Index m_nulls
Rank deficiency of m_A.
Definition: MatrixPower.h:413
Eigen::MatrixPower::m_conditionNumber
RealScalar m_conditionNumber
Condition number of m_A.
Definition: MatrixPower.h:407
Eigen::MatrixPower::m_A
MatrixType::Nested m_A
Reference to the base of matrix power.
Definition: MatrixPower.h:390

Member Function Documentation

◆ cols()

template<typename MatrixType >
Index Eigen::MatrixPower< MatrixType >::cols ( void  ) const
inline

Definition at line 382 of file MatrixPower.h.

382 { return m_A.cols(); }

◆ compute()

template<typename MatrixType >
template<typename ResultType >
void Eigen::MatrixPower< MatrixType >::compute ( ResultType &  res,
RealScalar  p 
)

Compute the matrix power.

Parameters
[in]pexponent, a real scalar.
[out]res\( A^p \) where A is specified in the constructor.

Definition at line 450 of file MatrixPower.h.

451 {
452  using std::pow;
453  switch (cols()) {
454  case 0:
455  break;
456  case 1:
457  res(0,0) = pow(m_A.coeff(0,0), p);
458  break;
459  default:
460  RealScalar intpart;
461  split(p, intpart);
462 
463  res = MatrixType::Identity(rows(), cols());
464  computeIntPower(res, intpart);
465  if (p) computeFracPower(res, p);
466  }
467 }
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Eigen::MatrixPower::cols
Index cols() const
Definition: MatrixPower.h:382
Eigen::MatrixPower::split
void split(RealScalar &p, RealScalar &intpart)
Split p into integral part and fractional part.
Definition: MatrixPower.h:470
Eigen::MatrixPower::RealScalar
MatrixType::RealScalar RealScalar
Definition: MatrixPower.h:343
Eigen::MatrixPower::computeFracPower
void computeFracPower(ResultType &res, RealScalar p)
Definition: MatrixPower.h:553
Eigen::MatrixPower::computeIntPower
void computeIntPower(ResultType &res, RealScalar p)
Definition: MatrixPower.h:530
Eigen::MatrixPower::rows
Index rows() const
Definition: MatrixPower.h:381
Eigen::pow
Eigen::AutoDiffScalar< EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(internal::remove_all_t< DerType >, typename internal::traits< internal::remove_all_t< DerType >>::Scalar, product) > pow(const Eigen::AutoDiffScalar< DerType > &x, const typename internal::traits< internal::remove_all_t< DerType >>::Scalar &y)
Definition: AutoDiffScalar.h:646

◆ computeFracPower()

template<typename MatrixType >
template<typename ResultType >
void Eigen::MatrixPower< MatrixType >::computeFracPower ( ResultType &  res,
RealScalar  p 
)
private

Definition at line 553 of file MatrixPower.h.

554 {
555  Block<ComplexMatrix,Dynamic,Dynamic> blockTp(m_fT, 0, 0, m_rank, m_rank);
556  eigen_assert(m_conditionNumber);
557  eigen_assert(m_rank + m_nulls == rows());
558 
559  MatrixPowerAtomic<ComplexMatrix>(m_T.topLeftCorner(m_rank, m_rank), p).compute(blockTp);
560  if (m_nulls) {
561  m_fT.topRightCorner(m_rank, m_nulls) = m_T.topLeftCorner(m_rank, m_rank).template triangularView<Upper>()
562  .solve(blockTp * m_T.topRightCorner(m_rank, m_nulls));
563  }
564  revertSchur(m_tmp, m_fT, m_U);
565  res = m_tmp * res;
566 }
Eigen::MatrixPower::m_fT
ComplexMatrix m_fT
Store fractional power of m_T.
Definition: MatrixPower.h:399
Eigen::MatrixPower::m_U
ComplexMatrix m_U
Definition: MatrixPower.h:396
Eigen::MatrixPower::m_tmp
MatrixType m_tmp
Temporary storage.
Definition: MatrixPower.h:393
Eigen::MatrixPower::revertSchur
static void revertSchur(Matrix< ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols > &res, const ComplexMatrix &T, const ComplexMatrix &U)
Definition: MatrixPower.h:570
Eigen::MatrixPower::m_T
ComplexMatrix m_T
Store the result of Schur decomposition.
Definition: MatrixPower.h:396

◆ computeIntPower()

template<typename MatrixType >
template<typename ResultType >
void Eigen::MatrixPower< MatrixType >::computeIntPower ( ResultType &  res,
RealScalar  p 
)
private

Definition at line 530 of file MatrixPower.h.

531 {
532  using std::abs;
533  using std::fmod;
534  RealScalar pp = abs(p);
535 
536  if (p<0)
537  m_tmp = m_A.inverse();
538  else
539  m_tmp = m_A;
540 
541  while (true) {
542  if (fmod(pp, 2) >= 1)
543  res = m_tmp * res;
544  pp /= 2;
545  if (pp < 1)
546  break;
547  m_tmp *= m_tmp;
548  }
549 }
fmod
bfloat16 fmod(const bfloat16 &a, const bfloat16 &b)
Eigen::abs
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
adtl::abs
adouble abs(const adouble &x)
Definition: AdolcForward:74

◆ initialize()

template<typename MatrixType >
void Eigen::MatrixPower< MatrixType >::initialize
private

Perform Schur decomposition for fractional power.

Definition at line 491 of file MatrixPower.h.

492 {
493  const ComplexSchur<MatrixType> schurOfA(m_A);
494  JacobiRotation<ComplexScalar> rot;
495  ComplexScalar eigenvalue;
496 
497  m_fT.resizeLike(m_A);
498  m_T = schurOfA.matrixT();
499  m_U = schurOfA.matrixU();
500  m_conditionNumber = m_T.diagonal().array().abs().maxCoeff() / m_T.diagonal().array().abs().minCoeff();
501 
502  // Move zero eigenvalues to the bottom right corner.
503  for (Index i = cols()-1; i>=0; --i) {
504  if (m_rank <= 2)
505  return;
506  if (m_T.coeff(i,i) == RealScalar(0)) {
507  for (Index j=i+1; j < m_rank; ++j) {
508  eigenvalue = m_T.coeff(j,j);
509  rot.makeGivens(m_T.coeff(j-1,j), eigenvalue);
510  m_T.applyOnTheRight(j-1, j, rot);
511  m_T.applyOnTheLeft(j-1, j, rot.adjoint());
512  m_T.coeffRef(j-1,j-1) = eigenvalue;
513  m_T.coeffRef(j,j) = RealScalar(0);
514  m_U.applyOnTheRight(j-1, j, rot);
515  }
516  --m_rank;
517  }
518  }
519 
520  m_nulls = rows() - m_rank;
521  if (m_nulls) {
522  eigen_assert(m_T.bottomRightCorner(m_nulls, m_nulls).isZero()
523  && "Base of matrix power should be invertible or with a semisimple zero eigenvalue.");
524  m_fT.bottomRows(m_nulls).fill(RealScalar(0));
525  }
526 }
i
int i
schurOfA
cout<< "Here is a random 4x4 matrix, A:"<< endl<< A<< endl<< endl;ComplexSchur< MatrixXcf > schurOfA(A, false)
Eigen::MatrixPower::ComplexScalar
std::complex< RealScalar > ComplexScalar
Definition: MatrixPower.h:385
Eigen::PlainObjectBase::coeffRef
constexpr Scalar & coeffRef(Index index)
Eigen::PlainObjectBase::resizeLike
void resizeLike(const EigenBase< OtherDerived > &_other)
Eigen::PlainObjectBase::coeff
const Scalar & coeff(Index index) const
Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
j
std::ptrdiff_t j

◆ operator()()

template<typename MatrixType >
const MatrixPowerParenthesesReturnValue<MatrixType> Eigen::MatrixPower< MatrixType >::operator() ( RealScalar  p)
inline

Returns the matrix power.

Parameters
[in]pexponent, a real scalar.
Returns
The expression \( A^p \), where A is specified in the constructor.

Definition at line 368 of file MatrixPower.h.

369  { return MatrixPowerParenthesesReturnValue<MatrixType>(*this, p); }

◆ revertSchur() [1/2]

template<typename MatrixType >
template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
void Eigen::MatrixPower< MatrixType >::revertSchur ( Matrix< ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols > &  res,
const ComplexMatrix &  T,
const ComplexMatrix &  U 
)
inlinestaticprivate

Definition at line 570 of file MatrixPower.h.

574 { res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint()); }

◆ revertSchur() [2/2]

template<typename MatrixType >
template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
void Eigen::MatrixPower< MatrixType >::revertSchur ( Matrix< RealScalar, Rows, Cols, Options, MaxRows, MaxCols > &  res,
const ComplexMatrix &  T,
const ComplexMatrix &  U 
)
inlinestaticprivate

Definition at line 578 of file MatrixPower.h.

582 { res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }

◆ rows()

template<typename MatrixType >
Index Eigen::MatrixPower< MatrixType >::rows ( void  ) const
inline

Definition at line 381 of file MatrixPower.h.

381 { return m_A.rows(); }

◆ split()

template<typename MatrixType >
void Eigen::MatrixPower< MatrixType >::split ( RealScalar &  p,
RealScalar &  intpart 
)
private

Split p into integral part and fractional part.

Parameters
[in]pThe exponent.
[out]pThe fractional part ranging in \( (-1, 1) \).
[out]intpartThe integral part.

Only if the fractional part is nonzero, it calls initialize().

Definition at line 470 of file MatrixPower.h.

471 {
472  using std::floor;
473  using std::pow;
474 
475  intpart = floor(p);
476  p -= intpart;
477 
478  // Perform Schur decomposition if it is not yet performed and the power is
479  // not an integer.
480  if (!m_conditionNumber && p)
481  initialize();
482 
483  // Choose the more stable of intpart = floor(p) and intpart = ceil(p).
484  if (p > RealScalar(0.5) && p > (1-p) * pow(m_conditionNumber, p)) {
485  --p;
486  ++intpart;
487  }
488 }
p
float * p
Eigen::MatrixPower::initialize
void initialize()
Perform Schur decomposition for fractional power.
Definition: MatrixPower.h:491
Eigen::floor
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_floor_op< typename Derived::Scalar >, const Derived > floor(const Eigen::ArrayBase< Derived > &x)

Member Data Documentation

◆ m_A

template<typename MatrixType >
MatrixType::Nested Eigen::MatrixPower< MatrixType >::m_A
private

Reference to the base of matrix power.

Definition at line 390 of file MatrixPower.h.

◆ m_conditionNumber

template<typename MatrixType >
RealScalar Eigen::MatrixPower< MatrixType >::m_conditionNumber
private

Condition number of m_A.

It is initialized as 0 to avoid performing unnecessary Schur decomposition, which is the bottleneck.

Definition at line 407 of file MatrixPower.h.

◆ m_fT

template<typename MatrixType >
ComplexMatrix Eigen::MatrixPower< MatrixType >::m_fT
private

Store fractional power of m_T.

Definition at line 399 of file MatrixPower.h.

◆ m_nulls

template<typename MatrixType >
Index Eigen::MatrixPower< MatrixType >::m_nulls
private

Rank deficiency of m_A.

Definition at line 413 of file MatrixPower.h.

◆ m_rank

template<typename MatrixType >
Index Eigen::MatrixPower< MatrixType >::m_rank
private

Rank of m_A.

Definition at line 410 of file MatrixPower.h.

◆ m_T

template<typename MatrixType >
ComplexMatrix Eigen::MatrixPower< MatrixType >::m_T
private

Store the result of Schur decomposition.

Definition at line 396 of file MatrixPower.h.

◆ m_tmp

template<typename MatrixType >
MatrixType Eigen::MatrixPower< MatrixType >::m_tmp
private

Temporary storage.

Definition at line 393 of file MatrixPower.h.

◆ m_U

template<typename MatrixType >
ComplexMatrix Eigen::MatrixPower< MatrixType >::m_U
private

Definition at line 396 of file MatrixPower.h.

The documentation for this class was generated from the following file:
Generated on Tue Jul 18 2023 15:28:51 for Eigen-unsupported by Doxygen 1.9.1 and the doxyYoda theme. This is an unofficial themed variant maintained by Rohit Goswami, based on the official documentation.