All Classes Namespaces Functions Variables Typedefs Enumerator Groups Pages
MatrixPower.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012, 2013 Chen-Pang He <[email protected]>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_MATRIX_POWER
11 #define EIGEN_MATRIX_POWER
12 
13 namespace Eigen {
14 
15 template<typename MatrixType> class MatrixPower;
16 
17 template<typename MatrixType>
18 class MatrixPowerRetval : public ReturnByValue< MatrixPowerRetval<MatrixType> >
19 {
20  public:
21  typedef typename MatrixType::RealScalar RealScalar;
22  typedef typename MatrixType::Index Index;
23 
24  MatrixPowerRetval(MatrixPower<MatrixType>& pow, RealScalar p) : m_pow(pow), m_p(p)
25  { }
26 
27  template<typename ResultType>
28  inline void evalTo(ResultType& res) const
29  { m_pow.compute(res, m_p); }
30 
31  Index rows() const { return m_pow.rows(); }
32  Index cols() const { return m_pow.cols(); }
33 
34  private:
35  MatrixPower<MatrixType>& m_pow;
36  const RealScalar m_p;
37  MatrixPowerRetval& operator=(const MatrixPowerRetval&);
38 };
39 
40 template<typename MatrixType>
41 class MatrixPowerAtomic
42 {
43  private:
44  enum {
45  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
46  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
47  };
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;
53 
54  const MatrixType& m_A;
55  RealScalar m_p;
56 
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);
65 
66  public:
67  MatrixPowerAtomic(const MatrixType& T, RealScalar p);
68  void compute(MatrixType& res) const;
69 };
70 
71 template<typename MatrixType>
72 MatrixPowerAtomic<MatrixType>::MatrixPowerAtomic(const MatrixType& T, RealScalar p) :
73  m_A(T), m_p(p)
74 { eigen_assert(T.rows() == T.cols()); }
75 
76 template<typename MatrixType>
77 void MatrixPowerAtomic<MatrixType>::compute(MatrixType& res) const
78 {
79  res.resizeLike(m_A);
80  switch (m_A.rows()) {
81  case 0:
82  break;
83  case 1:
84  res(0,0) = std::pow(m_A(0,0), m_p);
85  break;
86  case 2:
87  compute2x2(res, m_p);
88  break;
89  default:
90  computeBig(res);
91  }
92 }
93 
94 template<typename MatrixType>
95 void MatrixPowerAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res) const
96 {
97  int i = degree<<1;
98  res = (m_p-degree) / ((i-1)<<1) * IminusT;
99  for (--i; i; --i) {
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();
102  }
103  res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
104 }
105 
106 // This function assumes that res has the correct size (see bug 614)
107 template<typename MatrixType>
108 void MatrixPowerAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
109 {
110  using std::abs;
111  using std::pow;
112 
113  res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
114 
115  for (Index i=1; i < m_A.cols(); ++i) {
116  res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
117  if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i))
118  res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
119  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)))
120  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));
121  else
122  res.coeffRef(i-1,i) = computeSuperDiag(m_A.coeff(i,i), m_A.coeff(i-1,i-1), p);
123  res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
124  }
125 }
126 
127 template<typename MatrixType>
128 void MatrixPowerAtomic<MatrixType>::computeBig(MatrixType& res) const
129 {
130  const int digits = std::numeric_limits<RealScalar>::digits;
131  const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision
132  digits <= 53? 2.789358995219730e-1: // double precision
133  digits <= 64? 2.4471944416607995472e-1L: // extended precision
134  digits <= 106? 1.1016843812851143391275867258512e-1L: // double-double
135  9.134603732914548552537150753385375e-2L; // quadruple precision
136  MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
137  RealScalar normIminusT;
138  int degree, degree2, numberOfSquareRoots = 0;
139  bool hasExtraSquareRoot = false;
140 
141  /* FIXME
142  * For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite
143  * loop. We should move 0 eigenvalues to bottom right corner. We need not
144  * worry about tiny values (e.g. 1e-300) because they will reach 1 if
145  * repetitively sqrt'ed.
146  *
147  * If the 0 eigenvalues are semisimple, they can form a 0 matrix at the
148  * bottom right corner.
149  *
150  * [ T A ]^p [ T^p (T^-1 T^p A) ]
151  * [ ] = [ ]
152  * [ 0 0 ] [ 0 0 ]
153  */
154  for (Index i=0; i < m_A.cols(); ++i)
155  eigen_assert(m_A(i,i) != RealScalar(0));
156 
157  while (true) {
158  IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
159  normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
160  if (normIminusT < maxNormForPade) {
161  degree = getPadeDegree(normIminusT);
162  degree2 = getPadeDegree(normIminusT/2);
163  if (degree - degree2 <= 1 || hasExtraSquareRoot)
164  break;
165  hasExtraSquareRoot = true;
166  }
167  MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
168  T = sqrtT.template triangularView<Upper>();
169  ++numberOfSquareRoots;
170  }
171  computePade(degree, IminusT, res);
172 
173  for (; numberOfSquareRoots; --numberOfSquareRoots) {
174  compute2x2(res, std::ldexp(m_p, -numberOfSquareRoots));
175  res = res.template triangularView<Upper>() * res;
176  }
177  compute2x2(res, m_p);
178 }
179 
180 template<typename MatrixType>
181 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(float normIminusT)
182 {
183  const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f };
184  int degree = 3;
185  for (; degree <= 4; ++degree)
186  if (normIminusT <= maxNormForPade[degree - 3])
187  break;
188  return degree;
189 }
190 
191 template<typename MatrixType>
192 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(double normIminusT)
193 {
194  const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1,
195  1.999045567181744e-1, 2.789358995219730e-1 };
196  int degree = 3;
197  for (; degree <= 7; ++degree)
198  if (normIminusT <= maxNormForPade[degree - 3])
199  break;
200  return degree;
201 }
202 
203 template<typename MatrixType>
204 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(long double normIminusT)
205 {
206 #if LDBL_MANT_DIG == 53
207  const int maxPadeDegree = 7;
208  const double maxNormForPade[] = { 1.884160592658218e-2L /* degree = 3 */ , 6.038881904059573e-2L, 1.239917516308172e-1L,
209  1.999045567181744e-1L, 2.789358995219730e-1L };
210 #elif LDBL_MANT_DIG <= 64
211  const int maxPadeDegree = 8;
212  const double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L,
213  6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L };
214 #elif LDBL_MANT_DIG <= 106
215  const int maxPadeDegree = 10;
216  const double maxNormForPade[] = { 1.0007161601787493236741409687186e-4L /* degree = 3 */ ,
217  1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L,
218  2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L,
219  1.1016843812851143391275867258512e-1L };
220 #else
221  const int maxPadeDegree = 10;
222  const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ ,
223  6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L,
224  9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L,
225  3.908166513900489428442993794761185e-2L, 6.266780814639442865832535460550138e-2L,
226  9.134603732914548552537150753385375e-2L };
227 #endif
228  int degree = 3;
229  for (; degree <= maxPadeDegree; ++degree)
230  if (normIminusT <= maxNormForPade[degree - 3])
231  break;
232  return degree;
233 }
234 
235 template<typename MatrixType>
236 inline typename MatrixPowerAtomic<MatrixType>::ComplexScalar
237 MatrixPowerAtomic<MatrixType>::computeSuperDiag(const ComplexScalar& curr, const ComplexScalar& prev, RealScalar p)
238 {
239  ComplexScalar logCurr = std::log(curr);
240  ComplexScalar logPrev = std::log(prev);
241  int unwindingNumber = std::ceil((numext::imag(logCurr - logPrev) - M_PI) / (2*M_PI));
242  ComplexScalar w = numext::atanh2(curr - prev, curr + prev) + ComplexScalar(0, M_PI*unwindingNumber);
243  return RealScalar(2) * std::exp(RealScalar(0.5) * p * (logCurr + logPrev)) * std::sinh(p * w) / (curr - prev);
244 }
245 
246 template<typename MatrixType>
247 inline typename MatrixPowerAtomic<MatrixType>::RealScalar
248 MatrixPowerAtomic<MatrixType>::computeSuperDiag(RealScalar curr, RealScalar prev, RealScalar p)
249 {
250  RealScalar w = numext::atanh2(curr - prev, curr + prev);
251  return 2 * std::exp(p * (std::log(curr) + std::log(prev)) / 2) * std::sinh(p * w) / (curr - prev);
252 }
253 
273 template<typename MatrixType>
274 class MatrixPower
275 {
276  private:
277  enum {
278  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
279  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
280  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
281  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
282  };
283  typedef typename MatrixType::Scalar Scalar;
284  typedef typename MatrixType::RealScalar RealScalar;
285  typedef typename MatrixType::Index Index;
286 
287  public:
296  explicit MatrixPower(const MatrixType& A) : m_A(A), m_conditionNumber(0)
297  { eigen_assert(A.rows() == A.cols()); }
298 
306  const MatrixPowerRetval<MatrixType> operator()(RealScalar p)
307  { return MatrixPowerRetval<MatrixType>(*this, p); }
308 
316  template<typename ResultType>
317  void compute(ResultType& res, RealScalar p);
318 
319  Index rows() const { return m_A.rows(); }
320  Index cols() const { return m_A.cols(); }
321 
322  private:
323  typedef std::complex<RealScalar> ComplexScalar;
324  typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, MatrixType::Options,
325  MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrix;
326 
327  typename MatrixType::Nested m_A;
328  MatrixType m_tmp;
329  ComplexMatrix m_T, m_U, m_fT;
330  RealScalar m_conditionNumber;
331 
332  RealScalar modfAndInit(RealScalar, RealScalar*);
333 
334  template<typename ResultType>
335  void computeIntPower(ResultType&, RealScalar);
336 
337  template<typename ResultType>
338  void computeFracPower(ResultType&, RealScalar);
339 
340  template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
341  static void revertSchur(
342  Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
343  const ComplexMatrix& T,
344  const ComplexMatrix& U);
345 
346  template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
347  static void revertSchur(
348  Matrix<RealScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
349  const ComplexMatrix& T,
350  const ComplexMatrix& U);
351 };
352 
353 template<typename MatrixType>
354 template<typename ResultType>
355 void MatrixPower<MatrixType>::compute(ResultType& res, RealScalar p)
356 {
357  switch (cols()) {
358  case 0:
359  break;
360  case 1:
361  res(0,0) = std::pow(m_A.coeff(0,0), p);
362  break;
363  default:
364  RealScalar intpart, x = modfAndInit(p, &intpart);
365  computeIntPower(res, intpart);
366  computeFracPower(res, x);
367  }
368 }
369 
370 template<typename MatrixType>
371 typename MatrixPower<MatrixType>::RealScalar
372 MatrixPower<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
373 {
374  typedef Array<RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime> RealArray;
375 
376  *intpart = std::floor(x);
377  RealScalar res = x - *intpart;
378 
379  if (!m_conditionNumber && res) {
380  const ComplexSchur<MatrixType> schurOfA(m_A);
381  m_T = schurOfA.matrixT();
382  m_U = schurOfA.matrixU();
383 
384  const RealArray absTdiag = m_T.diagonal().array().abs();
385  m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff();
386  }
387 
388  if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber, res)) {
389  --res;
390  ++*intpart;
391  }
392  return res;
393 }
394 
395 template<typename MatrixType>
396 template<typename ResultType>
397 void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p)
398 {
399  RealScalar pp = std::abs(p);
400 
401  if (p<0) m_tmp = m_A.inverse();
402  else m_tmp = m_A;
403 
404  res = MatrixType::Identity(rows(), cols());
405  while (pp >= 1) {
406  if (std::fmod(pp, 2) >= 1)
407  res = m_tmp * res;
408  m_tmp *= m_tmp;
409  pp /= 2;
410  }
411 }
412 
413 template<typename MatrixType>
414 template<typename ResultType>
415 void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
416 {
417  if (p) {
418  eigen_assert(m_conditionNumber);
419  MatrixPowerAtomic<ComplexMatrix>(m_T, p).compute(m_fT);
420  revertSchur(m_tmp, m_fT, m_U);
421  res = m_tmp * res;
422  }
423 }
424 
425 template<typename MatrixType>
426 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
427 inline void MatrixPower<MatrixType>::revertSchur(
428  Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
429  const ComplexMatrix& T,
430  const ComplexMatrix& U)
431 { res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint()); }
432 
433 template<typename MatrixType>
434 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
435 inline void MatrixPower<MatrixType>::revertSchur(
436  Matrix<RealScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
437  const ComplexMatrix& T,
438  const ComplexMatrix& U)
439 { res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }
440 
454 template<typename Derived>
455 class MatrixPowerReturnValue : public ReturnByValue< MatrixPowerReturnValue<Derived> >
456 {
457  public:
458  typedef typename Derived::PlainObject PlainObject;
459  typedef typename Derived::RealScalar RealScalar;
460  typedef typename Derived::Index Index;
461 
468  MatrixPowerReturnValue(const Derived& A, RealScalar p) : m_A(A), m_p(p)
469  { }
470 
477  template<typename ResultType>
478  inline void evalTo(ResultType& res) const
479  { MatrixPower<PlainObject>(m_A.eval()).compute(res, m_p); }
480 
481  Index rows() const { return m_A.rows(); }
482  Index cols() const { return m_A.cols(); }
483 
484  private:
485  const Derived& m_A;
486  const RealScalar m_p;
487  MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&);
488 };
489 
490 namespace internal {
491 
492 template<typename MatrixPowerType>
493 struct traits< MatrixPowerRetval<MatrixPowerType> >
494 { typedef typename MatrixPowerType::PlainObject ReturnType; };
495 
496 template<typename Derived>
497 struct traits< MatrixPowerReturnValue<Derived> >
498 { typedef typename Derived::PlainObject ReturnType; };
499 
500 }
501 
502 template<typename Derived>
503 const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(const RealScalar& p) const
504 { return MatrixPowerReturnValue<Derived>(derived(), p); }
505 
506 } // namespace Eigen
507 
508 #endif // EIGEN_MATRIX_POWER
Eigen::MatrixPower
Class for computing matrix powers.
Definition: MatrixPower.h:15
Eigen::MatrixPowerReturnValue::MatrixPowerReturnValue
MatrixPowerReturnValue(const Derived &A, RealScalar p)
Constructor.
Definition: MatrixPower.h:468
Eigen::MatrixPower::MatrixPower
MatrixPower(const MatrixType &A)
Constructor.
Definition: MatrixPower.h:296
Eigen::MatrixPower::compute
void compute(ResultType &res, RealScalar p)
Compute the matrix power.
Definition: MatrixPower.h:355
Eigen::MatrixPower::operator()
const MatrixPowerRetval< MatrixType > operator()(RealScalar p)
Returns the matrix power.
Definition: MatrixPower.h:306
Eigen::MatrixPowerReturnValue::evalTo
void evalTo(ResultType &res) const
Compute the matrix power.
Definition: MatrixPower.h:478
Eigen::MatrixPowerReturnValue
Proxy for the matrix power of some matrix (expression).
Definition: MatrixPower.h:455