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 <jdh8@ms63.hinet.net>
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
13namespace Eigen {
14
15template<typename MatrixType> class MatrixPower;
16
17template<typename MatrixType>
18class 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
40template<typename MatrixType>
41class 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
71template<typename MatrixType>
72MatrixPowerAtomic<MatrixType>::MatrixPowerAtomic(const MatrixType& T, RealScalar p) :
73  m_A(T), m_p(p)
74{ eigen_assert(T.rows() == T.cols()); }
75
76template<typename MatrixType>
77void 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
94template<typename MatrixType>
95void 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)
107template<typename MatrixType>
108void MatrixPowerAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
109{
110  using std::abs;
111  using std::pow;
112
113  ArrayType logTdiag = m_A.diagonal().array().log();
114  res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
115
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));
122    else
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);
125  }
126}
127
128template<typename MatrixType>
129void MatrixPowerAtomic<MatrixType>::computeBig(MatrixType& res) const
130{
131  const int digits = std::numeric_limits<RealScalar>::digits;
132  const RealScalar maxNormForPade = digits <=  24? 4.3386528e-1f:                           // sigle precision
133				    digits <=  53? 2.789358995219730e-1:                    // double precision
134				    digits <=  64? 2.4471944416607995472e-1L:               // extended precision
135				    digits <= 106? 1.1016843812851143391275867258512e-1L:   // double-double
136						   9.134603732914548552537150753385375e-2L; // quadruple precision
137  MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
138  RealScalar normIminusT;
139  int degree, degree2, numberOfSquareRoots = 0;
140  bool hasExtraSquareRoot = false;
141
142  /* FIXME
143   * For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite
144   * loop.  We should move 0 eigenvalues to bottom right corner.  We need not
145   * worry about tiny values (e.g. 1e-300) because they will reach 1 if
146   * repetitively sqrt'ed.
147   *
148   * If the 0 eigenvalues are semisimple, they can form a 0 matrix at the
149   * bottom right corner.
150   *
151   * [ T  A ]^p   [ T^p  (T^-1 T^p A) ]
152   * [      ]   = [                   ]
153   * [ 0  0 ]     [  0         0      ]
154   */
155  for (Index i=0; i < m_A.cols(); ++i)
156    eigen_assert(m_A(i,i) != RealScalar(0));
157
158  while (true) {
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)
165	break;
166      hasExtraSquareRoot = true;
167    }
168    MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
169    T = sqrtT.template triangularView<Upper>();
170    ++numberOfSquareRoots;
171  }
172  computePade(degree, IminusT, res);
173
174  for (; numberOfSquareRoots; --numberOfSquareRoots) {
175    compute2x2(res, std::ldexp(m_p, -numberOfSquareRoots));
176    res = res.template triangularView<Upper>() * res;
177  }
178  compute2x2(res, m_p);
179}
180
181template<typename MatrixType>
182inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(float normIminusT)
183{
184  const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f };
185  int degree = 3;
186  for (; degree <= 4; ++degree)
187    if (normIminusT <= maxNormForPade[degree - 3])
188      break;
189  return degree;
190}
191
192template<typename MatrixType>
193inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(double normIminusT)
194{
195  const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1,
196      1.999045567181744e-1, 2.789358995219730e-1 };
197  int degree = 3;
198  for (; degree <= 7; ++degree)
199    if (normIminusT <= maxNormForPade[degree - 3])
200      break;
201  return degree;
202}
203
204template<typename MatrixType>
205inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(long double normIminusT)
206{
207#if   LDBL_MANT_DIG == 53
208  const int maxPadeDegree = 7;
209  const double maxNormForPade[] = { 1.884160592658218e-2L /* degree = 3 */ , 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 /* degree = 3 */ , 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 /* degree = 3 */ ,
218      1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L,
219      2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L,
220      1.1016843812851143391275867258512e-1L };
221#else
222  const int maxPadeDegree = 10;
223  const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ ,
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 };
228#endif
229  int degree = 3;
230  for (; degree <= maxPadeDegree; ++degree)
231    if (normIminusT <= maxNormForPade[degree - 3])
232      break;
233  return degree;
234}
235
236template<typename MatrixType>
237inline typename MatrixPowerAtomic<MatrixType>::ComplexScalar
238MatrixPowerAtomic<MatrixType>::computeSuperDiag(const ComplexScalar& curr, const ComplexScalar& prev, RealScalar p)
239{
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);
245}
246
247template<typename MatrixType>
248inline typename MatrixPowerAtomic<MatrixType>::RealScalar
249MatrixPowerAtomic<MatrixType>::computeSuperDiag(RealScalar curr, RealScalar prev, RealScalar p)
250{
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);
253}
254
255/**
256 * \ingroup MatrixFunctions_Module
257 *
258 * \brief Class for computing matrix powers.
259 *
260 * \tparam MatrixType  type of the base, expected to be an instantiation
261 * of the Matrix class template.
262 *
263 * This class is capable of computing real/complex matrices raised to
264 * an arbitrary real power. Meanwhile, it saves the result of Schur
265 * decomposition if an non-integral power has even been calculated.
266 * Therefore, if you want to compute multiple (>= 2) matrix powers
267 * for the same matrix, using the class directly is more efficient than
268 * calling MatrixBase::pow().
269 *
270 * Example:
271 * \include MatrixPower_optimal.cpp
272 * Output: \verbinclude MatrixPower_optimal.out
273 */
274template<typename MatrixType>
275class MatrixPower
276{
277  private:
278    enum {
279      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
280      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
281      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
282      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
283    };
284    typedef typename MatrixType::Scalar Scalar;
285    typedef typename MatrixType::RealScalar RealScalar;
286    typedef typename MatrixType::Index Index;
287
288  public:
289    /**
290     * \brief Constructor.
291     *
292     * \param[in] A  the base of the matrix power.
293     *
294     * The class stores a reference to A, so it should not be changed
295     * (or destroyed) before evaluation.
296     */
297    explicit MatrixPower(const MatrixType& A) : m_A(A), m_conditionNumber(0)
298    { eigen_assert(A.rows() == A.cols()); }
299
300    /**
301     * \brief Returns the matrix power.
302     *
303     * \param[in] p  exponent, a real scalar.
304     * \return The expression \f$ A^p \f$, where A is specified in the
305     * constructor.
306     */
307    const MatrixPowerRetval<MatrixType> operator()(RealScalar p)
308    { return MatrixPowerRetval<MatrixType>(*this, p); }
309
310    /**
311     * \brief Compute the matrix power.
312     *
313     * \param[in]  p    exponent, a real scalar.
314     * \param[out] res  \f$ A^p \f$ where A is specified in the
315     * constructor.
316     */
317    template<typename ResultType>
318    void compute(ResultType& res, RealScalar p);
319
320    Index rows() const { return m_A.rows(); }
321    Index cols() const { return m_A.cols(); }
322
323  private:
324    typedef std::complex<RealScalar> ComplexScalar;
325    typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, MatrixType::Options,
326              MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrix;
327
328    typename MatrixType::Nested m_A;
329    MatrixType m_tmp;
330    ComplexMatrix m_T, m_U, m_fT;
331    RealScalar m_conditionNumber;
332
333    RealScalar modfAndInit(RealScalar, RealScalar*);
334
335    template<typename ResultType>
336    void computeIntPower(ResultType&, RealScalar);
337
338    template<typename ResultType>
339    void computeFracPower(ResultType&, RealScalar);
340
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);
346
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);
352};
353
354template<typename MatrixType>
355template<typename ResultType>
356void MatrixPower<MatrixType>::compute(ResultType& res, RealScalar p)
357{
358  switch (cols()) {
359    case 0:
360      break;
361    case 1:
362      res(0,0) = std::pow(m_A.coeff(0,0), p);
363      break;
364    default:
365      RealScalar intpart, x = modfAndInit(p, &intpart);
366      computeIntPower(res, intpart);
367      computeFracPower(res, x);
368  }
369}
370
371template<typename MatrixType>
372typename MatrixPower<MatrixType>::RealScalar
373MatrixPower<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
374{
375  typedef Array<RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime> RealArray;
376
377  *intpart = std::floor(x);
378  RealScalar res = x - *intpart;
379
380  if (!m_conditionNumber && res) {
381    const ComplexSchur<MatrixType> schurOfA(m_A);
382    m_T = schurOfA.matrixT();
383    m_U = schurOfA.matrixU();
384
385    const RealArray absTdiag = m_T.diagonal().array().abs();
386    m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff();
387  }
388
389  if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber, res)) {
390    --res;
391    ++*intpart;
392  }
393  return res;
394}
395
396template<typename MatrixType>
397template<typename ResultType>
398void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p)
399{
400  RealScalar pp = std::abs(p);
401
402  if (p<0)  m_tmp = m_A.inverse();
403  else      m_tmp = m_A;
404
405  res = MatrixType::Identity(rows(), cols());
406  while (pp >= 1) {
407    if (std::fmod(pp, 2) >= 1)
408      res = m_tmp * res;
409    m_tmp *= m_tmp;
410    pp /= 2;
411  }
412}
413
414template<typename MatrixType>
415template<typename ResultType>
416void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
417{
418  if (p) {
419    eigen_assert(m_conditionNumber);
420    MatrixPowerAtomic<ComplexMatrix>(m_T, p).compute(m_fT);
421    revertSchur(m_tmp, m_fT, m_U);
422    res = m_tmp * res;
423  }
424}
425
426template<typename MatrixType>
427template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
428inline 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()); }
433
434template<typename MatrixType>
435template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
436inline 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(); }
441
442/**
443 * \ingroup MatrixFunctions_Module
444 *
445 * \brief Proxy for the matrix power of some matrix (expression).
446 *
447 * \tparam Derived  type of the base, a matrix (expression).
448 *
449 * This class holds the arguments to the matrix power until it is
450 * assigned or evaluated for some other reason (so the argument
451 * should not be changed in the meantime). It is the return type of
452 * MatrixBase::pow() and related functions and most of the
453 * time this is the only way it is used.
454 */
455template<typename Derived>
456class MatrixPowerReturnValue : public ReturnByValue< MatrixPowerReturnValue<Derived> >
457{
458  public:
459    typedef typename Derived::PlainObject PlainObject;
460    typedef typename Derived::RealScalar RealScalar;
461    typedef typename Derived::Index Index;
462
463    /**
464     * \brief Constructor.
465     *
466     * \param[in] A  %Matrix (expression), the base of the matrix power.
467     * \param[in] p  scalar, the exponent of the matrix power.
468     */
469    MatrixPowerReturnValue(const Derived& A, RealScalar p) : m_A(A), m_p(p)
470    { }
471
472    /**
473     * \brief Compute the matrix power.
474     *
475     * \param[out] result  \f$ A^p \f$ where \p A and \p p are as in the
476     * constructor.
477     */
478    template<typename ResultType>
479    inline void evalTo(ResultType& res) const
480    { MatrixPower<PlainObject>(m_A.eval()).compute(res, m_p); }
481
482    Index rows() const { return m_A.rows(); }
483    Index cols() const { return m_A.cols(); }
484
485  private:
486    const Derived& m_A;
487    const RealScalar m_p;
488    MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&);
489};
490
491namespace internal {
492
493template<typename MatrixPowerType>
494struct traits< MatrixPowerRetval<MatrixPowerType> >
495{ typedef typename MatrixPowerType::PlainObject ReturnType; };
496
497template<typename Derived>
498struct traits< MatrixPowerReturnValue<Derived> >
499{ typedef typename Derived::PlainObject ReturnType; };
500
501}
502
503template<typename Derived>
504const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(const RealScalar& p) const
505{ return MatrixPowerReturnValue<Derived>(derived(), p); }
506
507} // namespace Eigen
508
509#endif // EIGEN_MATRIX_POWER
510