10 #ifndef EIGEN_MATRIX_POWER
11 #define EIGEN_MATRIX_POWER
17 template<
typename MatrixType>
18 class MatrixPowerRetval :
public ReturnByValue< MatrixPowerRetval<MatrixType> >
21 typedef typename MatrixType::RealScalar RealScalar;
22 typedef typename MatrixType::Index Index;
27 template<
typename ResultType>
28 inline void evalTo(ResultType& res)
const
29 { m_pow.compute(res, m_p); }
31 Index rows()
const {
return m_pow.rows(); }
32 Index cols()
const {
return m_pow.cols(); }
35 MatrixPower<MatrixType>& m_pow;
37 MatrixPowerRetval& operator=(
const MatrixPowerRetval&);
40 template<
typename MatrixType>
41 class MatrixPowerAtomic
45 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
46 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
48 typedef typename MatrixType::Scalar Scalar;
49 typedef typename MatrixType::RealScalar RealScalar;
50 typedef std::complex<RealScalar> ComplexScalar;
51 typedef typename MatrixType::Index Index;
52 typedef Array<Scalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime> ArrayType;
54 const MatrixType& m_A;
57 void computePade(
int degree,
const MatrixType& IminusT, MatrixType& res)
const;
58 void compute2x2(MatrixType& res, RealScalar p)
const;
59 void computeBig(MatrixType& res)
const;
60 static int getPadeDegree(
float normIminusT);
61 static int getPadeDegree(
double normIminusT);
62 static int getPadeDegree(
long double normIminusT);
63 static ComplexScalar computeSuperDiag(
const ComplexScalar&,
const ComplexScalar&, RealScalar p);
64 static RealScalar computeSuperDiag(RealScalar, RealScalar, RealScalar p);
67 MatrixPowerAtomic(
const MatrixType& T, RealScalar p);
68 void compute(MatrixType& res)
const;
71 template<
typename MatrixType>
72 MatrixPowerAtomic<MatrixType>::MatrixPowerAtomic(
const MatrixType& T, RealScalar p) :
74 { eigen_assert(T.rows() == T.cols()); }
76 template<
typename MatrixType>
77 void MatrixPowerAtomic<MatrixType>::compute(MatrixType& res)
const
84 res(0,0) = std::pow(m_A(0,0), m_p);
94 template<
typename MatrixType>
95 void MatrixPowerAtomic<MatrixType>::computePade(
int degree,
const MatrixType& IminusT, MatrixType& res)
const
98 res = (m_p-degree) / ((i-1)<<1) * IminusT;
100 res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView<Upper>()
101 .solve((i==1 ? -m_p : i&1 ? (-m_p-(i>>1))/(i<<1) : (m_p-(i>>1))/((i-1)<<1)) * IminusT).eval();
103 res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
107 template<
typename MatrixType>
108 void MatrixPowerAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p)
const
113 ArrayType logTdiag = m_A.diagonal().array().log();
114 res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
116 for (Index i=1; i < m_A.cols(); ++i) {
117 res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
118 if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i))
119 res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
120 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)))
121 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));
123 res.coeffRef(i-1,i) = computeSuperDiag(m_A.coeff(i,i), m_A.coeff(i-1,i-1), p);
124 res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
128 template<
typename MatrixType>
129 void MatrixPowerAtomic<MatrixType>::computeBig(MatrixType& res)
const
131 const int digits = std::numeric_limits<RealScalar>::digits;
132 const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f:
133 digits <= 53? 2.789358995219730e-1:
134 digits <= 64? 2.4471944416607995472e-1L:
135 digits <= 106? 1.1016843812851143391275867258512e-1L:
136 9.134603732914548552537150753385375e-2L;
137 MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
138 RealScalar normIminusT;
139 int degree, degree2, numberOfSquareRoots = 0;
140 bool hasExtraSquareRoot =
false;
155 for (Index i=0; i < m_A.cols(); ++i)
156 eigen_assert(m_A(i,i) != RealScalar(0));
159 IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
160 normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
161 if (normIminusT < maxNormForPade) {
162 degree = getPadeDegree(normIminusT);
163 degree2 = getPadeDegree(normIminusT/2);
164 if (degree - degree2 <= 1 || hasExtraSquareRoot)
166 hasExtraSquareRoot =
true;
168 MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
169 T = sqrtT.template triangularView<Upper>();
170 ++numberOfSquareRoots;
172 computePade(degree, IminusT, res);
174 for (; numberOfSquareRoots; --numberOfSquareRoots) {
175 compute2x2(res, std::ldexp(m_p, -numberOfSquareRoots));
176 res = res.template triangularView<Upper>() * res;
178 compute2x2(res, m_p);
181 template<
typename MatrixType>
182 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(
float normIminusT)
184 const float maxNormForPade[] = { 2.8064004e-1f , 4.3386528e-1f };
186 for (; degree <= 4; ++degree)
187 if (normIminusT <= maxNormForPade[degree - 3])
192 template<
typename MatrixType>
193 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(
double normIminusT)
195 const double maxNormForPade[] = { 1.884160592658218e-2 , 6.038881904059573e-2, 1.239917516308172e-1,
196 1.999045567181744e-1, 2.789358995219730e-1 };
198 for (; degree <= 7; ++degree)
199 if (normIminusT <= maxNormForPade[degree - 3])
204 template<
typename MatrixType>
205 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(
long double normIminusT)
207 #if LDBL_MANT_DIG == 53
208 const int maxPadeDegree = 7;
209 const double maxNormForPade[] = { 1.884160592658218e-2L , 6.038881904059573e-2L, 1.239917516308172e-1L,
210 1.999045567181744e-1L, 2.789358995219730e-1L };
211 #elif LDBL_MANT_DIG <= 64
212 const int maxPadeDegree = 8;
213 const double maxNormForPade[] = { 6.3854693117491799460e-3L , 2.6394893435456973676e-2L,
214 6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L };
215 #elif LDBL_MANT_DIG <= 106
216 const int maxPadeDegree = 10;
217 const double maxNormForPade[] = { 1.0007161601787493236741409687186e-4L ,
218 1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L,
219 2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L,
220 1.1016843812851143391275867258512e-1L };
222 const int maxPadeDegree = 10;
223 const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L ,
224 6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L,
225 9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L,
226 3.908166513900489428442993794761185e-2L, 6.266780814639442865832535460550138e-2L,
227 9.134603732914548552537150753385375e-2L };
230 for (; degree <= maxPadeDegree; ++degree)
231 if (normIminusT <= maxNormForPade[degree - 3])
236 template<
typename MatrixType>
237 inline typename MatrixPowerAtomic<MatrixType>::ComplexScalar
238 MatrixPowerAtomic<MatrixType>::computeSuperDiag(
const ComplexScalar& curr,
const ComplexScalar& prev, RealScalar p)
240 ComplexScalar logCurr = std::log(curr);
241 ComplexScalar logPrev = std::log(prev);
242 int unwindingNumber = std::ceil((numext::imag(logCurr - logPrev) - M_PI) / (2*M_PI));
243 ComplexScalar w = numext::atanh2(curr - prev, curr + prev) + ComplexScalar(0, M_PI*unwindingNumber);
244 return RealScalar(2) * std::exp(RealScalar(0.5) * p * (logCurr + logPrev)) * std::sinh(p * w) / (curr - prev);
247 template<
typename MatrixType>
248 inline typename MatrixPowerAtomic<MatrixType>::RealScalar
249 MatrixPowerAtomic<MatrixType>::computeSuperDiag(RealScalar curr, RealScalar prev, RealScalar p)
251 RealScalar w = numext::atanh2(curr - prev, curr + prev);
252 return 2 * std::exp(p * (std::log(curr) + std::log(prev)) / 2) * std::sinh(p * w) / (curr - prev);
274 template<
typename MatrixType>
279 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
280 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
281 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
282 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
284 typedef typename MatrixType::Scalar Scalar;
285 typedef typename MatrixType::RealScalar RealScalar;
286 typedef typename MatrixType::Index Index;
297 explicit MatrixPower(
const MatrixType& A) : m_A(A), m_conditionNumber(0)
298 { eigen_assert(A.rows() == A.cols()); }
308 {
return MatrixPowerRetval<MatrixType>(*
this, p); }
317 template<
typename ResultType>
318 void compute(ResultType& res, RealScalar p);
320 Index rows()
const {
return m_A.rows(); }
321 Index cols()
const {
return m_A.cols(); }
324 typedef std::complex<RealScalar> ComplexScalar;
325 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, MatrixType::Options,
326 MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrix;
328 typename MatrixType::Nested m_A;
330 ComplexMatrix m_T, m_U, m_fT;
331 RealScalar m_conditionNumber;
333 RealScalar modfAndInit(RealScalar, RealScalar*);
335 template<
typename ResultType>
336 void computeIntPower(ResultType&, RealScalar);
338 template<
typename ResultType>
339 void computeFracPower(ResultType&, RealScalar);
341 template<
int Rows,
int Cols,
int Options,
int MaxRows,
int MaxCols>
342 static void revertSchur(
343 Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
344 const ComplexMatrix& T,
345 const ComplexMatrix& U);
347 template<
int Rows,
int Cols,
int Options,
int MaxRows,
int MaxCols>
348 static void revertSchur(
349 Matrix<RealScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
350 const ComplexMatrix& T,
351 const ComplexMatrix& U);
354 template<
typename MatrixType>
355 template<
typename ResultType>
362 res(0,0) = std::pow(m_A.coeff(0,0), p);
365 RealScalar intpart, x = modfAndInit(p, &intpart);
366 computeIntPower(res, intpart);
367 computeFracPower(res, x);
371 template<
typename MatrixType>
372 typename MatrixPower<MatrixType>::RealScalar
375 typedef Array<RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime> RealArray;
377 *intpart = std::floor(x);
378 RealScalar res = x - *intpart;
380 if (!m_conditionNumber && res) {
381 const ComplexSchur<MatrixType> schurOfA(m_A);
382 m_T = schurOfA.matrixT();
383 m_U = schurOfA.matrixU();
385 const RealArray absTdiag = m_T.diagonal().array().abs();
386 m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff();
389 if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber, res)) {
396 template<
typename MatrixType>
397 template<
typename ResultType>
398 void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p)
400 RealScalar pp = std::abs(p);
402 if (p<0) m_tmp = m_A.inverse();
405 res = MatrixType::Identity(rows(), cols());
407 if (std::fmod(pp, 2) >= 1)
414 template<
typename MatrixType>
415 template<
typename ResultType>
416 void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
419 eigen_assert(m_conditionNumber);
420 MatrixPowerAtomic<ComplexMatrix>(m_T, p).compute(m_fT);
421 revertSchur(m_tmp, m_fT, m_U);
426 template<
typename MatrixType>
427 template<
int Rows,
int Cols,
int Options,
int MaxRows,
int MaxCols>
428 inline void MatrixPower<MatrixType>::revertSchur(
429 Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
430 const ComplexMatrix& T,
431 const ComplexMatrix& U)
432 { res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint()); }
434 template<
typename MatrixType>
435 template<
int Rows,
int Cols,
int Options,
int MaxRows,
int MaxCols>
436 inline void MatrixPower<MatrixType>::revertSchur(
437 Matrix<RealScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
438 const ComplexMatrix& T,
439 const ComplexMatrix& U)
440 { res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }
455 template<
typename Derived>
459 typedef typename Derived::PlainObject PlainObject;
460 typedef typename Derived::RealScalar RealScalar;
461 typedef typename Derived::Index Index;
478 template<
typename ResultType>
479 inline void evalTo(ResultType& res)
const
482 Index rows()
const {
return m_A.rows(); }
483 Index cols()
const {
return m_A.cols(); }
487 const RealScalar m_p;
493 template<
typename MatrixPowerType>
494 struct traits< MatrixPowerRetval<MatrixPowerType> >
495 {
typedef typename MatrixPowerType::PlainObject ReturnType; };
497 template<
typename Derived>
498 struct traits< MatrixPowerReturnValue<Derived> >
499 {
typedef typename Derived::PlainObject ReturnType; };
503 template<
typename Derived>
504 const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(
const RealScalar& p)
const
505 {
return MatrixPowerReturnValue<Derived>(derived(), p); }
509 #endif // EIGEN_MATRIX_POWER
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: AutoDiffJacobian.h:13
MatrixPowerReturnValue(const Derived &A, RealScalar p)
Constructor.
Definition: MatrixPower.h:469
Class for computing matrix powers.
Definition: MatrixPower.h:15
void evalTo(ResultType &res) const
Compute the matrix power.
Definition: MatrixPower.h:479
const MatrixPowerRetval< MatrixType > operator()(RealScalar p)
Returns the matrix power.
Definition: MatrixPower.h:307
void compute(ResultType &res, RealScalar p)
Compute the matrix power.
Definition: MatrixPower.h:356
MatrixPower(const MatrixType &A)
Constructor.
Definition: MatrixPower.h:297
Proxy for the matrix power of some matrix (expression).
Definition: MatrixPower.h:456