diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index 124175841..0991817d5 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -59,7 +59,6 @@ #include "src/MatrixFunctions/MatrixFunction.h" #include "src/MatrixFunctions/MatrixSquareRoot.h" #include "src/MatrixFunctions/MatrixLogarithm.h" -#include "src/MatrixFunctions/MatrixPowerBase.h" #include "src/MatrixFunctions/MatrixPower.h" diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index e75fc25b4..12628df1c 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2012 Chen-Pang He +// Copyright (C) 2012, 2013 Chen-Pang He // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -12,6 +12,248 @@ namespace Eigen { +namespace MatrixPowerHelper { + +template +class ReturnValue : public ReturnByValue< ReturnValue > +{ + public: + typedef typename MatrixPowerType::PlainObject::RealScalar RealScalar; + typedef typename MatrixPowerType::PlainObject::Index Index; + + ReturnValue(MatrixPowerType& pow, RealScalar p) : m_pow(pow), m_p(p) + { } + + template + inline void evalTo(ResultType& res) const + { m_pow.compute(res, m_p); } + + Index rows() const { return m_pow.rows(); } + Index cols() const { return m_pow.cols(); } + + private: + MatrixPowerType& m_pow; + const RealScalar m_p; + ReturnValue& operator=(const ReturnValue&); +}; + +} // namespace MatrixPowerHelper + +template +class MatrixPowerAtomic +{ + private: + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime + }; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef std::complex ComplexScalar; + typedef typename MatrixType::Index Index; + typedef Array< Scalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > ArrayType; + + const MatrixType& m_A; + RealScalar m_p; + + void computePade(int degree, const MatrixType& IminusT, MatrixType& res) const; + void compute2x2(MatrixType& res, RealScalar p) const; + void computeBig(MatrixType& res) const; + 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); + + public: + MatrixPowerAtomic(const MatrixType& T, RealScalar p); + void compute(MatrixType& res) const; +}; + +template +MatrixPowerAtomic::MatrixPowerAtomic(const MatrixType& T, RealScalar p) : + m_A(T), m_p(p) +{ eigen_assert(T.rows() == T.cols()); } + +template +void MatrixPowerAtomic::compute(MatrixType& res) const +{ + res.resizeLike(m_A); + switch (m_A.rows()) { + case 0: + break; + case 1: + res(0,0) = std::pow(m_A(0,0), m_p); + break; + case 2: + compute2x2(res, m_p); + break; + default: + computeBig(res); + } +} + +template +void MatrixPowerAtomic::computePade(int degree, const MatrixType& IminusT, MatrixType& res) const +{ + int i = degree<<1; + res = (m_p-degree) / ((i-1)<<1) * IminusT; + for (--i; i; --i) { + res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView() + .solve((i==1 ? -m_p : i&1 ? (-m_p-(i>>1))/(i<<1) : (m_p-(i>>1))/((i-1)<<1)) * IminusT).eval(); + } + res += MatrixType::Identity(IminusT.rows(), IminusT.cols()); +} + +// This function assumes that res has the correct size (see bug 614) +template +void MatrixPowerAtomic::compute2x2(MatrixType& res, RealScalar p) const +{ + using std::abs; + using std::pow; + + ArrayType logTdiag = m_A.diagonal().array().log(); + 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); + } +} + +template +void MatrixPowerAtomic::computeBig(MatrixType& res) const +{ + const int digits = std::numeric_limits::digits; + const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision + digits <= 53? 2.789358995219730e-1: // 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(); + RealScalar normIminusT; + int degree, degree2, numberOfSquareRoots = 0; + bool hasExtraSquareRoot = false; + + /* FIXME + * For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite + * loop. We should move 0 eigenvalues to bottom right corner. We need not + * worry about tiny values (e.g. 1e-300) because they will reach 1 if + * repetitively sqrt'ed. + * + * If the 0 eigenvalues are semisimple, they can form a 0 matrix at the + * bottom right corner. + * + * [ T A ]^p [ T^p (T^-1 T^p A) ] + * [ ] = [ ] + * [ 0 0 ] [ 0 0 ] + */ + 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; + } + MatrixSquareRootTriangular(T).compute(sqrtT); + T = sqrtT.template triangularView(); + ++numberOfSquareRoots; + } + computePade(degree, IminusT, res); + + for (; numberOfSquareRoots; --numberOfSquareRoots) { + compute2x2(res, std::ldexp(m_p, -numberOfSquareRoots)); + res = res.template triangularView() * res; + } + compute2x2(res, m_p); +} + +template +inline int MatrixPowerAtomic::getPadeDegree(float normIminusT) +{ + 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; +} + +template +inline int MatrixPowerAtomic::getPadeDegree(double normIminusT) +{ + 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; +} + +template +inline int MatrixPowerAtomic::getPadeDegree(long double normIminusT) +{ +#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 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 <= maxNormForPade[degree - 3]) + break; + return degree; +} + +template +inline typename MatrixPowerAtomic::ComplexScalar +MatrixPowerAtomic::computeSuperDiag(const ComplexScalar& curr, const ComplexScalar& prev, RealScalar p) +{ + ComplexScalar logCurr = std::log(curr); + ComplexScalar logPrev = std::log(prev); + int unwindingNumber = std::ceil((numext::imag(logCurr - logPrev) - M_PI) / (2*M_PI)); + ComplexScalar w = numext::atanh2(curr - prev, curr + prev) + ComplexScalar(0, M_PI*unwindingNumber); + return RealScalar(2) * std::exp(RealScalar(0.5) * p * (logCurr + logPrev)) * std::sinh(p * w) / (curr - prev); +} + +template +inline typename MatrixPowerAtomic::RealScalar +MatrixPowerAtomic::computeSuperDiag(RealScalar curr, RealScalar prev, RealScalar p) +{ + RealScalar w = numext::atanh2(curr - prev, curr + prev); + return 2 * std::exp(p * (std::log(curr) + std::log(prev)) / 2) * std::sinh(p * w) / (curr - prev); +} + /** * \ingroup MatrixFunctions_Module * @@ -24,10 +266,22 @@ namespace Eigen { * to an arbitrary real power. */ template -class MatrixPowerTriangular : public MatrixPowerBase,MatrixType> +class MatrixPowerTriangular { + private: + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + public: - EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(MatrixPowerTriangular) + typedef MatrixType PlainObject; /** * \brief Constructor. @@ -37,10 +291,9 @@ class MatrixPowerTriangular : public MatrixPowerBase,MatrixType> operator()(RealScalar p); - #endif + const MatrixPowerHelper::ReturnValue operator()(RealScalar p) + { return MatrixPowerHelper::ReturnValue(*this, p); } /** * \brief Compute the matrix power. @@ -59,34 +312,19 @@ class MatrixPowerTriangular : public MatrixPowerBase - void compute(const Derived& b, ResultType& res, RealScalar p); + + Index rows() const { return m_A.rows(); } + Index cols() const { return m_A.cols(); } private: - EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(MatrixPowerTriangular) - - const TriangularView m_T; - + typename MatrixType::Nested m_A; + MatrixType m_tmp; + RealScalar m_conditionNumber; RealScalar modfAndInit(RealScalar, RealScalar*); - template - void apply(const Derived&, ResultType&, bool&); - template void computeIntPower(ResultType&, RealScalar); - template - void computeIntPower(const Derived&, ResultType&, RealScalar); - template void computeFracPower(ResultType&, RealScalar); }; @@ -94,41 +332,25 @@ class MatrixPowerTriangular : public MatrixPowerBase void MatrixPowerTriangular::compute(MatrixType& res, RealScalar p) { - switch (m_A.cols()) { + switch (cols()) { case 0: break; case 1: - res(0,0) = std::pow(m_T.coeff(0,0), p); + res(0,0) = std::pow(m_A.coeff(0,0), p); break; default: RealScalar intpart, x = modfAndInit(p, &intpart); - res = MatrixType::Identity(m_A.rows(), m_A.cols()); computeIntPower(res, intpart); computeFracPower(res, x); } } -template -template -void MatrixPowerTriangular::compute(const Derived& b, ResultType& res, RealScalar p) -{ - switch (m_A.cols()) { - case 0: - break; - case 1: - res = std::pow(m_T.coeff(0,0), p) * b; - break; - default: - RealScalar intpart, x = modfAndInit(p, &intpart); - computeIntPower(b, res, intpart); - computeFracPower(res, x); - } -} - template typename MatrixPowerTriangular::RealScalar MatrixPowerTriangular::modfAndInit(RealScalar x, RealScalar* intpart) { + typedef Array< RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > RealArray; + *intpart = std::floor(x); RealScalar res = x - *intpart; @@ -137,95 +359,39 @@ MatrixPowerTriangular::modfAndInit(RealScalar x, RealScalar* intpart m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff(); } - if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber,res)) { + if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber, res)) { --res; ++*intpart; } return res; } -template -template -void MatrixPowerTriangular::apply(const Derived& b, ResultType& res, bool& init) -{ - if (init) - res = m_tmp1.template triangularView() * res; - else { - init = true; - res.noalias() = m_tmp1.template triangularView() * b; - } -} - template template void MatrixPowerTriangular::computeIntPower(ResultType& res, RealScalar p) { RealScalar pp = std::abs(p); - if (p<0) m_tmp1 = m_T.solve(MatrixType::Identity(m_A.rows(), m_A.cols())); - else m_tmp1 = m_T; + if (p<0) m_tmp = m_A.template triangularView().solve(MatrixType::Identity(rows(), cols())); + else m_tmp = m_A.template triangularView(); + res = MatrixType::Identity(rows(), cols()); while (pp >= 1) { if (std::fmod(pp, 2) >= 1) - res = m_tmp1.template triangularView() * res; - m_tmp1 = m_tmp1.template triangularView() * m_tmp1; + res.template triangularView() = m_tmp.template triangularView() * res; + m_tmp.template triangularView() = m_tmp.template triangularView() * m_tmp; pp /= 2; } } -template -template -void MatrixPowerTriangular::computeIntPower(const Derived& b, ResultType& res, RealScalar p) -{ - if (b.cols() >= m_A.cols()) { - m_tmp2 = MatrixType::Identity(m_A.rows(), m_A.cols()); - computeIntPower(m_tmp2, p); - res.noalias() = m_tmp2.template triangularView() * b; - } - else { - RealScalar pp = std::abs(p); - int squarings, applyings = internal::binary_powering_cost(pp, &squarings); - bool init = false; - - if (p==0) { - res = b; - return; - } - else if (p>0) { - m_tmp1 = m_T; - } - else if (b.cols()*(pp-applyings) <= m_A.cols()*squarings) { - res = m_T.solve(b); - for (--pp; pp >= 1; --pp) - res = m_T.solve(res); - return; - } - else { - m_tmp1 = m_T.solve(MatrixType::Identity(m_A.rows(), m_A.cols())); - } - - while (b.cols()*(pp-applyings) > m_A.cols()*squarings) { - if (std::fmod(pp, 2) >= 1) { - apply(b, res, init); - --applyings; - } - m_tmp1 = m_tmp1.template triangularView() * m_tmp1; - --squarings; - pp /= 2; - } - for (; pp >= 1; --pp) - apply(b, res, init); - } -} - template template void MatrixPowerTriangular::computeFracPower(ResultType& res, RealScalar p) { if (p) { eigen_assert(m_conditionNumber); - MatrixPowerTriangularAtomic(m_A).compute(m_tmp1, p); - res = m_tmp1.template triangularView() * res; + MatrixPowerAtomic(m_A, p).compute(m_tmp); + res = m_tmp * res; } } @@ -249,10 +415,22 @@ void MatrixPowerTriangular::computeFracPower(ResultType& res, RealSc * Output: \verbinclude MatrixPower_optimal.out */ template -class MatrixPower : public MatrixPowerBase,MatrixType> +class MatrixPower { + private: + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + public: - EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(MatrixPower) + typedef MatrixType PlainObject; /** * \brief Constructor. @@ -262,10 +440,9 @@ class MatrixPower : public MatrixPowerBase,MatrixType> * The class stores a reference to A, so it should not be changed * (or destroyed) before evaluation. */ - explicit MatrixPower(const MatrixType& A) : Base(A) - { } + explicit MatrixPower(const MatrixType& A) : m_A(A), m_conditionNumber(0) + { eigen_assert(A.rows() == A.cols()); } - #ifdef EIGEN_PARSED_BY_DOXYGEN /** * \brief Returns the matrix power. * @@ -273,8 +450,8 @@ class MatrixPower : public MatrixPowerBase,MatrixType> * \return The expression \f$ A^p \f$, where A is specified in the * constructor. */ - const MatrixPowerBaseReturnValue,MatrixType> operator()(RealScalar p); - #endif + const MatrixPowerHelper::ReturnValue operator()(RealScalar p) + { return MatrixPowerHelper::ReturnValue(*this, p); } /** * \brief Compute the matrix power. @@ -284,45 +461,45 @@ class MatrixPower : public MatrixPowerBase,MatrixType> * constructor. */ void compute(MatrixType& res, RealScalar p); - - /** - * \brief Compute the matrix power multiplied by another matrix. - * - * \param[in] b a matrix with the same rows as A. - * \param[in] p exponent, a real scalar. - * \param[out] res \f$ A^p b \f$, where A is specified in the - * constructor. - */ - template - void compute(const Derived& b, ResultType& res, RealScalar p); + + Index rows() const { return m_A.rows(); } + Index cols() const { return m_A.cols(); } private: - EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(MatrixPower) + typedef std::complex ComplexScalar; + typedef Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, + MaxColsAtCompileTime > ComplexMatrix; - typedef Matrix, RowsAtCompileTime, ColsAtCompileTime, - Options,MaxRowsAtCompileTime,MaxColsAtCompileTime> ComplexMatrix; - static const bool m_OKforLU = RowsAtCompileTime == Dynamic || RowsAtCompileTime > 4; + typename MatrixType::Nested m_A; + MatrixType m_tmp; ComplexMatrix m_T, m_U, m_fT; + RealScalar m_conditionNumber; RealScalar modfAndInit(RealScalar, RealScalar*); - template - void apply(const Derived&, ResultType&, bool&); - template void computeIntPower(ResultType&, RealScalar); - template - void computeIntPower(const Derived&, ResultType&, RealScalar); - template void computeFracPower(ResultType&, RealScalar); + + template + static void revertSchur( + Matrix< ComplexScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res, + const ComplexMatrix& T, + const ComplexMatrix& U); + + template + static void revertSchur( + Matrix< RealScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res, + const ComplexMatrix& T, + const ComplexMatrix& U); }; template void MatrixPower::compute(MatrixType& res, RealScalar p) { - switch (m_A.cols()) { + switch (cols()) { case 0: break; case 1: @@ -330,32 +507,17 @@ void MatrixPower::compute(MatrixType& res, RealScalar p) break; default: RealScalar intpart, x = modfAndInit(p, &intpart); - res = MatrixType::Identity(m_A.rows(), m_A.cols()); computeIntPower(res, intpart); computeFracPower(res, x); } } template -template -void MatrixPower::compute(const Derived& b, ResultType& res, RealScalar p) +typename MatrixPower::RealScalar +MatrixPower::modfAndInit(RealScalar x, RealScalar* intpart) { - switch (m_A.cols()) { - case 0: - break; - case 1: - res = std::pow(m_A.coeff(0,0), p) * b; - break; - default: - RealScalar intpart, x = modfAndInit(p, &intpart); - computeIntPower(b, res, intpart); - computeFracPower(res, x); - } -} + typedef Array< RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > RealArray; -template -typename MatrixPower::RealScalar MatrixPower::modfAndInit(RealScalar x, RealScalar* intpart) -{ *intpart = std::floor(x); RealScalar res = x - *intpart; @@ -375,100 +537,51 @@ typename MatrixPower::RealScalar MatrixPower::modfAndIni return res; } -template -template -void MatrixPower::apply(const Derived& b, ResultType& res, bool& init) -{ - if (init) - res = m_tmp1 * res; - else { - init = true; - res.noalias() = m_tmp1 * b; - } -} - template template void MatrixPower::computeIntPower(ResultType& res, RealScalar p) { RealScalar pp = std::abs(p); - if (p<0) m_tmp1 = m_A.inverse(); - else m_tmp1 = m_A; + if (p<0) m_tmp = m_A.inverse(); + else m_tmp = m_A; + res = MatrixType::Identity(rows(), cols()); while (pp >= 1) { if (std::fmod(pp, 2) >= 1) - res = m_tmp1 * res; - m_tmp1 *= m_tmp1; + res = m_tmp * res; + m_tmp *= m_tmp; pp /= 2; } } -template -template -void MatrixPower::computeIntPower(const Derived& b, ResultType& res, RealScalar p) -{ - if (b.cols() >= m_A.cols()) { - m_tmp2 = MatrixType::Identity(m_A.rows(), m_A.cols()); - computeIntPower(m_tmp2, p); - res.noalias() = m_tmp2 * b; - } - else { - RealScalar pp = std::abs(p); - int squarings, applyings = internal::binary_powering_cost(pp, &squarings); - bool init = false; - - if (p==0) { - res = b; - return; - } - else if (p>0) { - m_tmp1 = m_A; - } - else if (m_OKforLU && b.cols()*(pp-applyings) <= m_A.cols()*squarings) { - PartialPivLU A(m_A); - res = A.solve(b); - for (--pp; pp >= 1; --pp) - res = A.solve(res); - return; - } - else { - m_tmp1 = m_A.inverse(); - } - - while (b.cols()*(pp-applyings) > m_A.cols()*squarings) { - if (std::fmod(pp, 2) >= 1) { - apply(b, res, init); - --applyings; - } - m_tmp1 *= m_tmp1; - --squarings; - pp /= 2; - } - for (; pp >= 1; --pp) - apply(b, res, init); - } -} - template template void MatrixPower::computeFracPower(ResultType& res, RealScalar p) { if (p) { eigen_assert(m_conditionNumber); - MatrixPowerTriangularAtomic(m_T).compute(m_fT, p); - internal::recompose_complex_schur::IsComplex>::run(m_tmp1, m_fT, m_U); - res = m_tmp1 * res; + MatrixPowerAtomic(m_T, p).compute(m_fT); + revertSchur(m_tmp, m_fT, m_U); + res = m_tmp * res; } } -namespace internal { +template +template +inline void MatrixPower::revertSchur( + Matrix< ComplexScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res, + const ComplexMatrix& T, + const ComplexMatrix& U) +{ res.noalias() = U * (T.template triangularView() * U.adjoint()); } -template -struct traits > -{ typedef typename Derived::PlainObject ReturnType; }; - -} // namespace internal +template +template +inline void MatrixPower::revertSchur( + Matrix< RealScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res, + const ComplexMatrix& T, + const ComplexMatrix& U) +{ res.noalias() = (U * (T.template triangularView() * U.adjoint())).real(); } /** * \ingroup MatrixFunctions_Module @@ -484,7 +597,7 @@ struct traits > * time this is the only way it is used. */ template -class MatrixPowerReturnValue : public ReturnByValue > +class MatrixPowerReturnValue : public ReturnByValue< MatrixPowerReturnValue > { public: typedef typename Derived::PlainObject PlainObject; @@ -510,21 +623,6 @@ class MatrixPowerReturnValue : public ReturnByValue(m_A.eval()).compute(res, m_p); } - /** - * \brief Return the expression \f$ A^p b \f$. - * - * \p A and \p p are specified in the constructor. - * - * \param[in] b the matrix (expression) to be applied. - */ - template - const MatrixPowerProduct,PlainObject,OtherDerived> - operator*(const MatrixBase& b) const - { - MatrixPower Apow(m_A.eval()); - return MatrixPowerProduct,PlainObject,OtherDerived>(Apow, b.derived(), m_p); - } - Index rows() const { return m_A.rows(); } Index cols() const { return m_A.cols(); } @@ -534,6 +632,18 @@ class MatrixPowerReturnValue : public ReturnByValue +struct traits< MatrixPowerHelper::ReturnValue > +{ typedef typename MatrixPowerType::PlainObject ReturnType; }; + +template +struct traits< MatrixPowerReturnValue > +{ typedef typename Derived::PlainObject ReturnType; }; + +} + template const MatrixPowerReturnValue MatrixBase::pow(const RealScalar& p) const { return MatrixPowerReturnValue(derived(), p); } diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h deleted file mode 100644 index 25b3f4496..000000000 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h +++ /dev/null @@ -1,404 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2012 Chen-Pang He -// -// This Source Code Form is subject to the terms of the Mozilla -// Public License v. 2.0. If a copy of the MPL was not distributed -// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#ifndef EIGEN_MATRIX_POWER_BASE -#define EIGEN_MATRIX_POWER_BASE - -namespace Eigen { - -#define EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(Derived) \ - typedef MatrixPowerBase Base; \ - using Base::RowsAtCompileTime; \ - using Base::ColsAtCompileTime; \ - using Base::Options; \ - using Base::MaxRowsAtCompileTime; \ - using Base::MaxColsAtCompileTime; \ - typedef typename Base::Scalar Scalar; \ - typedef typename Base::RealScalar RealScalar; \ - typedef typename Base::RealArray RealArray; - -#define EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(Derived) \ - using Base::m_A; \ - using Base::m_tmp1; \ - using Base::m_tmp2; \ - using Base::m_conditionNumber; - -template -class MatrixPowerBaseReturnValue : public ReturnByValue > -{ - public: - typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; - - MatrixPowerBaseReturnValue(Derived& pow, RealScalar p) : m_pow(pow), m_p(p) - { } - - template - inline void evalTo(ResultType& res) const - { m_pow.compute(res, m_p); } - - template - const MatrixPowerProduct operator*(const MatrixBase& b) const - { return MatrixPowerProduct(m_pow, b.derived(), m_p); } - - Index rows() const { return m_pow.rows(); } - Index cols() const { return m_pow.cols(); } - - private: - Derived& m_pow; - const RealScalar m_p; - MatrixPowerBaseReturnValue& operator=(const MatrixPowerBaseReturnValue&); -}; - -template -class MatrixPowerBase -{ - private: - Derived& derived() { return *static_cast(this); } - - public: - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime, - Options = MatrixType::Options, - MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, - MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime - }; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; - - explicit MatrixPowerBase(const MatrixType& A) : - m_A(A), - m_conditionNumber(0) - { eigen_assert(A.rows() == A.cols()); } - - #ifndef EIGEN_PARSED_BY_DOXYGEN - const MatrixPowerBaseReturnValue operator()(RealScalar p) - { return MatrixPowerBaseReturnValue(derived(), p); } - #endif - - void compute(MatrixType& res, RealScalar p) - { derived().compute(res,p); } - - template - void compute(const OtherDerived& b, ResultType& res, RealScalar p) - { derived().compute(b,res,p); } - - Index rows() const { return m_A.rows(); } - Index cols() const { return m_A.cols(); } - - protected: - typedef Array RealArray; - - typename MatrixType::Nested m_A; - MatrixType m_tmp1, m_tmp2; - RealScalar m_conditionNumber; -}; - -template -class MatrixPowerProduct : public MatrixBase > -{ - public: - typedef MatrixBase Base; - EIGEN_DENSE_PUBLIC_INTERFACE(MatrixPowerProduct) - - MatrixPowerProduct(Derived& pow, const Rhs& b, RealScalar p) : - m_pow(pow), - m_b(b), - m_p(p) - { eigen_assert(pow.cols() == b.rows()); } - - template - inline void evalTo(ResultType& res) const - { m_pow.compute(m_b, res, m_p); } - - inline Index rows() const { return m_pow.rows(); } - inline Index cols() const { return m_b.cols(); } - - private: - Derived& m_pow; - typename Rhs::Nested m_b; - const RealScalar m_p; -}; - -template -template -Derived& MatrixBase::lazyAssign(const MatrixPowerProduct& other) -{ - other.evalTo(derived()); - return derived(); -} - -namespace internal { - -template -struct traits > -{ typedef MatrixType ReturnType; }; - -template -struct traits > -{ - typedef MatrixXpr XprKind; - typedef typename remove_all<_Lhs>::type Lhs; - typedef typename remove_all<_Rhs>::type Rhs; - typedef typename scalar_product_traits::ReturnType Scalar; - typedef Dense StorageKind; - typedef typename promote_index_type::type Index; - - enum { - RowsAtCompileTime = traits::RowsAtCompileTime, - ColsAtCompileTime = traits::ColsAtCompileTime, - MaxRowsAtCompileTime = traits::MaxRowsAtCompileTime, - MaxColsAtCompileTime = traits::MaxColsAtCompileTime, - Flags = (MaxRowsAtCompileTime==1 ? RowMajorBit : 0) - | EvalBeforeNestingBit | EvalBeforeAssigningBit | NestByRefBit, - CoeffReadCost = 0 - }; -}; - -template -struct recompose_complex_schur -{ - template - static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U) - { res.noalias() = U * (T.template triangularView() * U.adjoint()); } -}; - -template<> -struct recompose_complex_schur<0> -{ - template - static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U) - { res.noalias() = (U * (T.template triangularView() * U.adjoint())).real(); } -}; - -template::IsComplex> -struct matrix_power_unwinder -{ - static inline Scalar run(const Scalar& eival, const Scalar& eival0, int unwindingNumber) - { return numext::atanh2(eival-eival0, eival+eival0) + Scalar(0, M_PI*unwindingNumber); } -}; - -template -struct matrix_power_unwinder -{ - static inline Scalar run(Scalar eival, Scalar eival0, int) - { return numext::atanh2(eival-eival0, eival+eival0); } -}; - -template -inline int binary_powering_cost(T p, int* squarings) -{ - int applyings=0, tmp; - - frexp(p, squarings); - --*squarings; - - while (std::frexp(p, &tmp), tmp > 0) { - p -= std::ldexp(static_cast(0.5), tmp); - ++applyings; - } - return applyings; -} - -inline int matrix_power_get_pade_degree(float normIminusT) -{ - 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; -} - -inline int matrix_power_get_pade_degree(double normIminusT) -{ - 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; -} - -inline int matrix_power_get_pade_degree(long double normIminusT) -{ -#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 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 <= maxNormForPade[degree - 3]) - break; - return degree; -} - -} // namespace internal - -template -class MatrixPowerTriangularAtomic -{ - private: - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime - }; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; - typedef Array ArrayType; - - const MatrixType& m_A; - - static void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p); - void compute2x2(MatrixType& res, RealScalar p) const; - void computeBig(MatrixType& res, RealScalar p) const; - - public: - explicit MatrixPowerTriangularAtomic(const MatrixType& T); - void compute(MatrixType& res, RealScalar p) const; -}; - -template -MatrixPowerTriangularAtomic::MatrixPowerTriangularAtomic(const MatrixType& T) : - m_A(T) -{ eigen_assert(T.rows() == T.cols()); } - -template -void MatrixPowerTriangularAtomic::compute(MatrixType& res, RealScalar p) const -{ - res.resizeLike(m_A); - switch (m_A.rows()) { - case 0: - break; - case 1: - res(0,0) = std::pow(m_A(0,0), p); - break; - case 2: - compute2x2(res, p); - break; - default: - computeBig(res, p); - } -} - -template -void MatrixPowerTriangularAtomic::computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) -{ - int i = degree<<1; - res = (p-degree) / ((i-1)<<1) * IminusT; - for (--i; i; --i) { - res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView() - .solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : (p-(i>>1))/((i-1)<<1)) * IminusT).eval(); - } - res += MatrixType::Identity(IminusT.rows(), IminusT.cols()); -} - -// this function assumes that res has the correct size (see bug 614) -template -void MatrixPowerTriangularAtomic::compute2x2(MatrixType& res, RealScalar p) const -{ - using std::abs; - using std::pow; - - ArrayType logTdiag = m_A.diagonal().array().log(); - 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 { - int unwindingNumber = std::ceil((numext::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI)); - Scalar w = internal::matrix_power_unwinder::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber); - res.coeffRef(i-1,i) = RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) * std::sinh(p * w) - / (m_A.coeff(i,i) - m_A.coeff(i-1,i-1)); - } - res.coeffRef(i-1,i) *= m_A.coeff(i-1,i); - } -} - -template -void MatrixPowerTriangularAtomic::computeBig(MatrixType& res, RealScalar p) const -{ - const int digits = std::numeric_limits::digits; - const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision - digits <= 53? 2.789358995219730e-1: // 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(); - RealScalar normIminusT; - int degree, degree2, numberOfSquareRoots = 0; - bool hasExtraSquareRoot = false; - - /* FIXME - * For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite - * loop. We should move 0 eigenvalues to bottom right corner. We need not - * worry about tiny values (e.g. 1e-300) because they will reach 1 if - * repetitively sqrt'ed. - * - * [ T A ]^p [ T^p (T^-1 T^p A) ] - * [ ] = [ ] - * [ 0 0 ] [ 0 0 ] - */ - 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 = internal::matrix_power_get_pade_degree(normIminusT); - degree2 = internal::matrix_power_get_pade_degree(normIminusT/2); - if (degree - degree2 <= 1 || hasExtraSquareRoot) - break; - hasExtraSquareRoot = true; - } - MatrixSquareRootTriangular(T).compute(sqrtT); - T = sqrtT.template triangularView(); - ++numberOfSquareRoots; - } - computePade(degree, IminusT, res, p); - - for (; numberOfSquareRoots; --numberOfSquareRoots) { - compute2x2(res, std::ldexp(p,-numberOfSquareRoots)); - res = res.template triangularView() * res; - } - compute2x2(res, p); -} - -} // namespace Eigen - -#endif // EIGEN_MATRIX_POWER diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index 4bb2b4d03..b9d513b45 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -109,62 +109,9 @@ void testExponentLaws(const MatrixType& m, double tol) } } -template -void testProduct(const MatrixType& m, const VectorType& v, double tol) -{ - typedef typename MatrixType::RealScalar RealScalar; - MatrixType m1; - VectorType v1, v2, v3; - RealScalar p; - - for (int i=0; i < g_repeat; ++i) { - generateTestMatrix::run(m1, m.rows()); - MatrixPower mpow(m1); - - v1 = VectorType::Random(v.rows(), v.cols()); - p = internal::random(); - - v2.noalias() = mpow(p) * v1; - v3.noalias() = mpow(p).eval() * v1; - std::cout << "testProduct: error powerm = " << relerr(v2, v3) << '\n'; - VERIFY(v2.isApprox(v3, static_cast(tol))); - } -} - -template -void testTriangularProduct(const MatrixType& m, const VectorType& v, double tol) -{ - typedef typename MatrixType::RealScalar RealScalar; - MatrixType m1; - VectorType v1, v2, v3; - RealScalar p; - - for (int i=0; i < g_repeat; ++i) { - generateTriangularMatrix::run(m1, m.rows()); - MatrixPowerTriangular mpow(m1); - - v1 = VectorType::Random(v.rows(), v.cols()); - p = internal::random(); - - v2.noalias() = mpow(p) * v1; - v3.noalias() = mpow(p).eval() * v1; - std::cout << "testTriangularProduct: error powerm = " << relerr(v2, v3) << '\n'; - VERIFY(v2.isApprox(v3, static_cast(tol))); - } -} - -template -void testMatrixVector(const MatrixType& m, const VectorType& v, double tol) -{ - testExponentLaws(m,tol); - testProduct(m,v,tol); - testTriangularProduct(m,v,tol); -} - typedef Matrix Matrix3dRowMajor; typedef Matrix MatrixXe; -typedef Matrix VectorXe; - + void test_matrix_power() { CALL_SUBTEST_2(test2dRotation(1e-13)); @@ -174,13 +121,13 @@ void test_matrix_power() CALL_SUBTEST_1(test2dHyperbolicRotation(1e-5)); CALL_SUBTEST_9(test2dHyperbolicRotation(1e-14)); - CALL_SUBTEST_2(testMatrixVector(Matrix2d(), Vector2d(), 1e-13)); - CALL_SUBTEST_7(testMatrixVector(Matrix3dRowMajor(), MatrixXd(3,5), 1e-13)); - CALL_SUBTEST_3(testMatrixVector(Matrix4cd(), Vector4cd(), 1e-13)); - CALL_SUBTEST_4(testMatrixVector(MatrixXd(8,8), VectorXd(8), 2e-12)); - CALL_SUBTEST_1(testMatrixVector(Matrix2f(), Vector2f(), 1e-4)); - CALL_SUBTEST_5(testMatrixVector(Matrix3cf(), Vector3cf(), 1e-4)); - CALL_SUBTEST_8(testMatrixVector(Matrix4f(), Vector4f(), 1e-4)); - CALL_SUBTEST_6(testMatrixVector(MatrixXf(2,2), VectorXf(2), 1e-3)); // see bug 614 - CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7), VectorXe(7), 1e-13)); + CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13)); + CALL_SUBTEST_7(testExponentLaws(Matrix3dRowMajor(), 1e-13)); + CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13)); + CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 2e-12)); + CALL_SUBTEST_1(testExponentLaws(Matrix2f(), 1e-4)); + CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4)); + CALL_SUBTEST_8(testExponentLaws(Matrix4f(), 1e-4)); + CALL_SUBTEST_6(testExponentLaws(MatrixXf(2,2), 1e-3)); // see bug 614 + CALL_SUBTEST_9(testExponentLaws(MatrixXe(7,7), 1e-13)); }