diff options
-rw-r--r-- | unsupported/Eigen/MatrixFunctions | 48 | ||||
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixPower.h | 849 | ||||
-rw-r--r-- | unsupported/doc/examples/MatrixPower.cpp | 2 | ||||
-rw-r--r-- | unsupported/doc/examples/MatrixPower_optimal.cpp | 17 | ||||
-rw-r--r-- | unsupported/test/matrix_power.cpp | 102 |
5 files changed, 463 insertions, 555 deletions
diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index 002c1f71c..ffebe0324 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -211,8 +211,9 @@ documentation of \ref matrixbase_exp "exp()". \include MatrixLogarithm.cpp Output: \verbinclude MatrixLogarithm.out -\note \p M has to be a matrix of \c float, \c double, \c long double -\c complex<float>, \c complex<double>, or \c complex<long double> . +\note \p M has to be a matrix of \c float, \c double, <tt>long +double</tt>, \c complex<float>, \c complex<double>, or \c complex<long +double> . \sa MatrixBase::exp(), MatrixBase::matrixFunction(), class MatrixLogarithmAtomic, MatrixBase::sqrt(). @@ -234,27 +235,14 @@ where exp denotes the matrix exponential, and log denotes the matrix logarithm. The matrix \f$ M \f$ should meet the conditions to be an argument of -matrix logarithm. If \p p is neither an integer nor the real scalar -type of \p M, it is casted into the real scalar type of \p M. - -This function computes the matrix logarithm using the -Schur-Padé algorithm as implemented by MatrixBase::pow(). -The exponent is split into integral part and fractional part, where -the fractional part is in the interval \f$ (-1, 1) \f$. The main -diagonal and the first super-diagonal is directly computed. - -The actual work is done by the MatrixPower class, which can compute -\f$ M^p v \f$, where \p v is another matrix with the same rows as -\p M. The matrix \p v is set to be the identity matrix by default. -Therefore, the expression <tt>M.pow(p) * v</tt> is specialized for -this. No temporary storage is created for the result. The code below -directly evaluates R-values into L-values without aliasing issue. Do -\b NOT try to \a optimize with noalias(). It won't compile. -\code -v = m.pow(p) * v; -m = m.pow(q); -// v2.noalias() = m.pow(p) * v1; Won't compile! -\endcode +matrix logarithm. If \p p is not of the real scalar type of \p M, it +is casted into the real scalar type of \p M. + +This function computes the matrix power using the Schur-Padé +algorithm as implemented by class MatrixPower. The exponent is split +into integral part and fractional part, where the fractional part is +in the interval \f$ (-1, 1) \f$. The main diagonal and the first +super-diagonal is directly computed. Details of the algorithm can be found in: Nicholas J. Higham and Lijing Lin, "A Schur-Padé algorithm for fractional powers of a @@ -277,8 +265,18 @@ the z-axis. \include MatrixPower.cpp Output: \verbinclude MatrixPower.out -\note \p M has to be a matrix of \c float, \c double, \c long double -\c complex<float>, \c complex<double>, or \c complex<long double> . +MatrixBase::pow() is user-friendly. However, there are some +circumstances under which you should use class MatrixPower directly. +MatrixPower can save the result of Schur decomposition, so it's +better for computing various powers for the same matrix. + +Example: +\include MatrixPower_optimal.cpp +Output: \verbinclude MatrixPower_optimal.out + +\note \p M has to be a matrix of \c float, \c double, <tt>long +double</tt>, \c complex<float>, \c complex<double>, or \c complex<long +double> . \sa MatrixBase::exp(), MatrixBase::log(), class MatrixPower. diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 7238501ed..9b61502ed 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -10,395 +10,69 @@ #ifndef EIGEN_MATRIX_POWER #define EIGEN_MATRIX_POWER -#ifndef M_PI -#define M_PI 3.141592653589793238462643383279503L -#endif - namespace Eigen { -/** - * \ingroup MatrixFunctions_Module - * - * \brief Class for computing matrix powers. - * - * \tparam MatrixType type of the base, expected to be an instantiation - * of the Matrix class template. - * \tparam PlainObject type of the multiplier. - */ -template<typename MatrixType, typename PlainObject = MatrixType> -class MatrixPower -{ - private: - typedef internal::traits<MatrixType> Traits; - static const int Rows = Traits::RowsAtCompileTime; - static const int Cols = Traits::ColsAtCompileTime; - static const int Options = Traits::Options; - static const int MaxRows = Traits::MaxRowsAtCompileTime; - static const int MaxCols = Traits::MaxColsAtCompileTime; - - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef std::complex<RealScalar> ComplexScalar; - typedef typename MatrixType::Index Index; - typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix; - typedef Array<ComplexScalar, Rows, 1, ColMajor, MaxRows> ComplexArray; - - public: - /** - * \brief Constructor. - * - * \param[in] A the base of the matrix power. - * \param[in] p the exponent of the matrix power. - * \param[in] b the multiplier. - */ - MatrixPower(const MatrixType& A, RealScalar p, const PlainObject& b) : - m_A(A), - m_p(p), - m_b(b), - m_dimA(A.cols()), - m_dimb(b.cols()) - { /* empty body */ } - - /** - * \brief Compute the matrix power. - * - * \param[out] result \f$ A^p b \f$, as specified in the constructor. - */ - template<typename ResultType> void compute(ResultType& result); - - private: - /** - * \brief Compute the matrix to integral power. - * - * If \p b is \em fatter than \p A, it computes \f$ A^{p_{\textrm int}} - * \f$ first, and then multiplies it with \p b. Otherwise, - * #computeChainProduct optimizes the expression. - * - * \sa computeChainProduct(ResultType&); - */ - template<typename ResultType> - void computeIntPower(ResultType& result); - - /** - * \brief Convert integral power of the matrix into chain product. - * - * This optimizes the matrix expression. It automatically chooses binary - * powering or matrix chain multiplication or solving the linear system - * repetitively, according to which algorithm costs less. - */ - template<typename ResultType> - void computeChainProduct(ResultType&); - - /** \brief Compute the cost of binary powering. */ - static int computeCost(RealScalar); - - /** \brief Solve the linear system repetitively. */ - template<typename ResultType> - void partialPivLuSolve(ResultType&, RealScalar); - - /** \brief Compute Schur decomposition of #m_A. */ - void computeSchurDecomposition(); - - /** - * \brief Split #m_p into integral part and fractional part. - * - * This method stores the integral part \f$ p_{\textrm int} \f$ into - * #m_pInt and the fractional part \f$ p_{\textrm frac} \f$ into - * #m_pFrac, where #m_pFrac is in the interval \f$ (-1,1) \f$. To - * choose between the possibilities below, it considers the computation - * of \f$ A^{p_1} \f$ and \f$ A^{p_2} \f$ and determines which of these - * computations is the better conditioned. - */ - void getFractionalExponent(); - - /** \brief Compute power of 2x2 triangular matrix. */ - void compute2x2(RealScalar p); - - /** - * \brief Compute power of triangular matrices with size > 2. - * \details This uses a Schur-Padé algorithm. - */ - void computeBig(); - - /** \brief Get suitable degree for Pade approximation. (specialized for RealScalar = double) */ - inline int getPadeDegree(double); - - /** \brief Get suitable degree for Pade approximation. (specialized for RealScalar = float) */ - inline int getPadeDegree(float); - - /** \brief Get suitable degree for Pade approximation. (specialized for RealScalar = long double) */ - inline int getPadeDegree(long double); - - /** \brief Compute Padé approximation to matrix fractional power. */ - void computePade(const int& degree, const ComplexMatrix& IminusT); - - /** \brief Get a certain coefficient of the Padé approximation. */ - inline RealScalar coeff(const int& degree); - - /** - * \brief Store the fractional power into #m_tmp. - * - * This intended for complex matrices. - */ - void computeTmp(ComplexScalar); - - /** - * \brief Store the fractional power into #m_tmp. - * - * This is intended for real matrices. It takes the real part of - * \f$ U T^{p_{\textrm frac}} U^H \f$. - * - * \sa computeTmp(ComplexScalar); - */ - void computeTmp(RealScalar); - - const MatrixType& m_A; ///< \brief Reference to the matrix base. - const RealScalar m_p; ///< \brief The real exponent. - const PlainObject& m_b; ///< \brief Reference to the multiplier. - const Index m_dimA; ///< \brief The dimension of #m_A, equivalent to %m_A.cols(). - const Index m_dimb; ///< \brief The dimension of #m_b, equivalent to %m_b.cols(). - MatrixType m_tmp; ///< \brief Used for temporary storage. - RealScalar m_pInt; ///< \brief Integral part of #m_p. - RealScalar m_pFrac; ///< \brief Fractional part of #m_p. - ComplexMatrix m_T; ///< \brief Triangular part of Schur decomposition. - ComplexMatrix m_U; ///< \brief Unitary part of Schur decomposition. - ComplexMatrix m_fT; ///< \brief #m_T to the power of #m_pFrac. - ComplexArray m_logTdiag; ///< \brief Logarithm of the main diagonal of #m_T. -}; - -template<typename MatrixType, typename PlainObject> -template<typename ResultType> -void MatrixPower<MatrixType,PlainObject>::compute(ResultType& result) -{ - using std::floor; - using std::pow; - - m_pInt = floor(m_p + RealScalar(0.5)); - m_pFrac = m_p - m_pInt; - - if (!m_pFrac) { - computeIntPower(result); - } else if (m_dimA == 1) - result = pow(m_A(0,0), m_p) * m_b; - else { - computeSchurDecomposition(); - getFractionalExponent(); - computeIntPower(result); - if (m_dimA == 2) - compute2x2(m_pFrac); - else - computeBig(); - computeTmp(Scalar()); - result = m_tmp * result; - } -} +namespace internal { -template<typename MatrixType, typename PlainObject> -template<typename ResultType> -void MatrixPower<MatrixType,PlainObject>::computeIntPower(ResultType& result) +template<int IsComplex> +struct recompose_complex_schur { - MatrixType tmp; - if (m_dimb > m_dimA) { - tmp = MatrixType::Identity(m_dimA, m_dimA); - computeChainProduct(tmp); - result.noalias() = tmp * m_b; - } else { - result = m_b; - computeChainProduct(result); - } -} + template<typename ResultType, typename MatrixType> + static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U) + { res = U * (T.template triangularView<Upper>() * U.adjoint()); } +}; -template<typename MatrixType, typename PlainObject> -template<typename ResultType> -void MatrixPower<MatrixType,PlainObject>::computeChainProduct(ResultType& result) +template<> +struct recompose_complex_schur<0> { - using std::abs; - using std::fmod; - using std::ldexp; - - RealScalar p = abs(m_pInt); - int cost = computeCost(p); - - if (m_pInt < RealScalar(0)) { - if (p * m_dimb <= cost * m_dimA && m_dimA > 2) { - partialPivLuSolve(result, p); - return; - } else { - m_tmp = m_A.inverse(); - } - } else { - m_tmp = m_A; - } - while (p * m_dimb > cost * m_dimA) { - if (fmod(p, RealScalar(2)) >= RealScalar(1)) { - result = m_tmp * result; - cost--; - } - m_tmp *= m_tmp; - cost--; - p = ldexp(p, -1); - } - for (; p >= RealScalar(1); p--) - result = m_tmp * result; -} + template<typename ResultType, typename MatrixType> + static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U) + { res = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); } +}; -template<typename MatrixType, typename PlainObject> -int MatrixPower<MatrixType,PlainObject>::computeCost(RealScalar p) +template<typename T> +inline int binary_powering_cost(T p) { - using std::frexp; - using std::ldexp; int cost, tmp; - frexp(p, &cost); - while (frexp(p, &tmp), tmp > 0) { - p -= ldexp(RealScalar(0.5), tmp); - cost++; + while (std::frexp(p, &tmp), tmp > 0) { + p -= std::ldexp(static_cast<T>(0.5), tmp); + ++cost; } return cost; } -template<typename MatrixType, typename PlainObject> -template<typename ResultType> -void MatrixPower<MatrixType,PlainObject>::partialPivLuSolve(ResultType& result, RealScalar p) -{ - const PartialPivLU<MatrixType> Asolver(m_A); - for (; p >= RealScalar(1); p--) - result = Asolver.solve(result); -} - -template<typename MatrixType, typename PlainObject> -void MatrixPower<MatrixType,PlainObject>::computeSchurDecomposition() -{ - const ComplexSchur<MatrixType> schurOfA(m_A); - m_T = schurOfA.matrixT(); - m_U = schurOfA.matrixU(); -} - -template<typename MatrixType, typename PlainObject> -void MatrixPower<MatrixType,PlainObject>::getFractionalExponent() -{ - using std::pow; - typedef Array<RealScalar, Rows, 1, ColMajor, MaxRows> RealArray; - - const ComplexArray Tdiag = m_T.diagonal(); - const RealArray absTdiag = Tdiag.abs(); - const RealScalar maxAbsEival = absTdiag.maxCoeff(); - const RealScalar minAbsEival = absTdiag.minCoeff(); - - m_logTdiag = Tdiag.log(); - if (m_pFrac > RealScalar(0.5) && // This is just a shortcut. - m_pFrac > (RealScalar(1) - m_pFrac) * pow(maxAbsEival/minAbsEival, m_pFrac)) { - m_pFrac--; - m_pInt++; - } -} - -template<typename MatrixType, typename PlainObject> -void MatrixPower<MatrixType,PlainObject>::compute2x2(RealScalar p) -{ - using std::abs; - using std::ceil; - using std::exp; - using std::imag; - using std::ldexp; - using std::pow; - using std::sinh; - - int i, j, unwindingNumber; - ComplexScalar w; - - m_fT(0,0) = pow(m_T(0,0), p); - for (j = 1; j < m_dimA; j++) { - i = j - 1; - m_fT(j,j) = pow(m_T(j,j), p); - - if (m_T(i,i) == m_T(j,j)) { - m_fT(i,j) = p * pow(m_T(i,j), p - RealScalar(1)); - } else if (abs(m_T(i,i)) < ldexp(abs(m_T(j,j)), -1) || abs(m_T(j,j)) < ldexp(abs(m_T(i,i)), -1)) { - m_fT(i,j) = m_T(i,j) * (m_fT(j,j) - m_fT(i,i)) / (m_T(j,j) - m_T(i,i)); - } else { - // computation in previous branch is inaccurate if abs(m_T(j,j)) \approx abs(m_T(i,i)) - unwindingNumber = ceil((imag(m_logTdiag[j] - m_logTdiag[i]) - M_PI) / (2 * M_PI)); - w = internal::atanh2(m_T(j,j) - m_T(i,i), m_T(j,j) + m_T(i,i)) + ComplexScalar(0, M_PI * unwindingNumber); - m_fT(i,j) = m_T(i,j) * RealScalar(2) * exp(RealScalar(0.5) * p * (m_logTdiag[j] + m_logTdiag[i])) * - sinh(p * w) / (m_T(j,j) - m_T(i,i)); - } - } -} - -template<typename MatrixType, typename PlainObject> -void MatrixPower<MatrixType,PlainObject>::computeBig() -{ - using std::ldexp; - const int digits = std::numeric_limits<RealScalar>::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-01: // double-double - 9.134603732914548552537150753385375e-02; // quadruple precision - int degree, degree2, numberOfSquareRoots = 0, numberOfExtraSquareRoots = 0; - ComplexMatrix IminusT, sqrtT, T = m_T; - RealScalar normIminusT; - - while (true) { - IminusT = ComplexMatrix::Identity(m_dimA, m_dimA) - T; - normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff(); - if (normIminusT < maxNormForPade) { - degree = getPadeDegree(normIminusT); - degree2 = getPadeDegree(normIminusT * RealScalar(0.5)); - if (degree - degree2 <= 1 || numberOfExtraSquareRoots) - break; - numberOfExtraSquareRoots++; - } - MatrixSquareRootTriangular<ComplexMatrix>(T).compute(sqrtT); - T = sqrtT; - numberOfSquareRoots++; - } - computePade(degree, IminusT); - - for (; numberOfSquareRoots; numberOfSquareRoots--) { - compute2x2(ldexp(m_pFrac, -numberOfSquareRoots)); - m_fT *= m_fT; - } - compute2x2(m_pFrac); -} - -template<typename MatrixType, typename PlainObject> -inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(float normIminusT) +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++) + for (; degree <= 4; ++degree) if (normIminusT <= maxNormForPade[degree - 3]) break; return degree; } -template<typename MatrixType, typename PlainObject> -inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(double normIminusT) +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 }; + 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++) + for (; degree <= 7; ++degree) if (normIminusT <= maxNormForPade[degree - 3]) break; return degree; } -template<typename MatrixType, typename PlainObject> -inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(long double normIminusT) +inline int matrix_power_get_pade_degree(long double normIminusT) { -#if LDBL_MANT_DIG == 53 +#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 }; - + 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 */ , @@ -414,101 +88,369 @@ inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(long double normIm 9.134603732914548552537150753385375e-2L }; #endif int degree = 3; - for (; degree <= maxPadeDegree; degree++) + for (; degree <= maxPadeDegree; ++degree) if (normIminusT <= maxNormForPade[degree - 3]) break; return degree; } -template<typename MatrixType, typename PlainObject> -void MatrixPower<MatrixType,PlainObject>::computePade(const int& degree, const ComplexMatrix& IminusT) +} // namespace internal + +/* (non-doc) + * \brief Class for computing triangular matrices to fractional power. + * + * \tparam MatrixType type of the base. + */ +template<typename MatrixType, int UpLo = Upper> class MatrixPowerTriangularAtomic { - int i = degree << 1; - m_fT = coeff(i) * IminusT; - for (i--; i; i--) { - m_fT = (ComplexMatrix::Identity(m_dimA, m_dimA) + m_fT).template triangularView<Upper>() - .solve(coeff(i) * IminusT).eval(); + private: + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Array<Scalar, + EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), + 1,ColMajor, + EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::MaxRowsAtCompileTime,MatrixType::MaxColsAtCompileTime)> ArrayType; + const MatrixType& m_T; + + void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) const; + 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<typename MatrixType, int UpLo> +MatrixPowerTriangularAtomic<MatrixType,UpLo>::MatrixPowerTriangularAtomic(const MatrixType& T) : + m_T(T) +{ eigen_assert(T.rows() == T.cols()); } + +template<typename MatrixType, int UpLo> +void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute(MatrixType& res, RealScalar p) const +{ + switch (m_T.rows()) { + case 0: + break; + case 1: + res(0,0) = std::pow(m_T(0,0), p); + break; + case 2: + compute2x2(res, p); + break; + default: + computeBig(res, p); } - m_fT += ComplexMatrix::Identity(m_dimA, m_dimA); } -template<typename MatrixType, typename PlainObject> -inline typename MatrixType::RealScalar MatrixPower<MatrixType,PlainObject>::coeff(const int& i) +template<typename MatrixType, int UpLo> +void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computePade(int degree, const MatrixType& IminusT, MatrixType& res, + RealScalar p) const { - if (i == 1) - return -m_pFrac; - else if (i & 1) - return (-m_pFrac - RealScalar(i >> 1)) / RealScalar(i << 1); - else - return (m_pFrac - RealScalar(i >> 1)) / RealScalar((i - 1) << 1); + int i = degree<<1; + res = (p-(i>>1)) / ((i-1)<<1) * IminusT; + for (--i; i; --i) { + res = (MatrixType::Identity(m_T.rows(), m_T.cols()) + res).template triangularView<UpLo>() + .solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : (p-(i>>1))/((i-1)<<1)) * IminusT).eval(); + } + res += MatrixType::Identity(m_T.rows(), m_T.cols()); } -template<typename MatrixType, typename PlainObject> -void MatrixPower<MatrixType,PlainObject>::computeTmp(RealScalar) -{ m_tmp = (m_U * m_fT * m_U.adjoint()).real(); } +template<typename MatrixType, int UpLo> +void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute2x2(MatrixType& res, RealScalar p) const +{ + using std::abs; + using std::pow; + + ArrayType logTdiag = m_T.diagonal().array().log(); + res(0,0) = pow(m_T(0,0), p); + + for (int i=1; i < m_T.cols(); ++i) { + res(i,i) = pow(m_T(i,i), p); + if (m_T(i-1,i-1) == m_T(i,i)) { + res(i-1,i) = p * pow(m_T(i-1,i), p-1); + } else if (2*abs(m_T(i-1,i-1)) < abs(m_T(i,i)) || 2*abs(m_T(i,i)) < abs(m_T(i-1,i-1))) { + res(i-1,i) = m_T(i-1,i) * (res(i,i)-res(i-1,i-1)) / (m_T(i,i)-m_T(i-1,i-1)); + } else { + // computation in previous branch is inaccurate if abs(m_T(i,i)) \approx abs(m_T(i-1,i-1)) + int unwindingNumber = std::ceil(((logTdiag[i]-logTdiag[i-1]).imag() - M_PI) / (2*M_PI)); + Scalar w = internal::atanh2(m_T(i,i)-m_T(i-1,i-1), m_T(i,i)+m_T(i-1,i-1)) + Scalar(0, M_PI*unwindingNumber); + res(i-1,i) = m_T(i-1,i) * RealScalar(2) * std::exp(RealScalar(0.5) * p * (logTdiag[i]+logTdiag[i-1])) * + std::sinh(p * w) / (m_T(i,i) - m_T(i-1,i-1)); + } + } +} -template<typename MatrixType, typename PlainObject> -void MatrixPower<MatrixType,PlainObject>::computeTmp(ComplexScalar) -{ m_tmp = m_U * m_fT * m_U.adjoint(); } +template<typename MatrixType, int UpLo> +void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computeBig(MatrixType& res, RealScalar p) const +{ + const int digits = std::numeric_limits<RealScalar>::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-01: // double-double + 9.134603732914548552537150753385375e-02; // quadruple precision + int degree, degree2, numberOfSquareRoots=0, numberOfExtraSquareRoots=0; + MatrixType IminusT, sqrtT, T=m_T; + RealScalar normIminusT; + + while (true) { + IminusT = MatrixType::Identity(m_T.rows(), m_T.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 || numberOfExtraSquareRoots) + break; + ++numberOfExtraSquareRoots; + } + MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT); + T = sqrtT; + ++numberOfSquareRoots; + } + computePade(degree, IminusT, res, p); + + for (; numberOfSquareRoots; --numberOfSquareRoots) { + compute2x2(res, std::ldexp(p,-numberOfSquareRoots)); + res *= res; + } + compute2x2(res, p); +} /** * \ingroup MatrixFunctions_Module * - * \brief Proxy for the matrix power multiplied by other matrix. + * \brief Class for computing matrix powers. * - * \tparam Derived type of the base, a matrix (expression). - * \tparam RhsDerived type of the multiplier. + * \tparam MatrixType type of the base, expected to be an instantiation + * of the Matrix class template. * - * This class holds the arguments to the matrix power until it is - * assigned or evaluated for some other reason (so the argument - * should not be changed in the meantime). It is the return type of - * MatrixPowerReturnValue::operator*() and related functions and most - * of the time this is the only way it is used. + * 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 MatrixPower_optimal.cpp + * Output: \verbinclude MatrixPower_optimal.out */ -template<typename Derived, typename RhsDerived> -class MatrixPowerProductReturnValue : public ReturnByValue<MatrixPowerProductReturnValue<Derived,RhsDerived> > +template<typename MatrixType> class MatrixPower { private: - typedef typename Derived::PlainObject MatrixType; - typedef typename RhsDerived::PlainObject PlainObject; - typedef typename RhsDerived::RealScalar RealScalar; - typedef typename RhsDerived::Index Index; + static const int Rows = MatrixType::RowsAtCompileTime; + static const int Cols = MatrixType::ColsAtCompileTime; + static const int Options = MatrixType::Options; + static const int MaxRows = MatrixType::MaxRowsAtCompileTime; + static const int MaxCols = MatrixType::MaxColsAtCompileTime; + + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + typedef Matrix<std::complex<RealScalar>,Rows,Cols,Options,MaxRows,MaxCols> ComplexMatrix; + + const MatrixType& m_A; + MatrixType m_tmp1, m_tmp2; + ComplexMatrix m_T, m_U, m_fT; + bool m_init; + + RealScalar modfAndInit(RealScalar, RealScalar*); + + template<typename PlainObject, typename ResultType> + void apply(const PlainObject&, ResultType&, bool&); + + template<typename ResultType> + void computeIntPower(ResultType&, RealScalar); + + template<typename PlainObject, typename ResultType> + void computeIntPower(const PlainObject&, ResultType&, RealScalar); + + template<typename ResultType> + void computeFracPower(ResultType&, RealScalar); public: /** * \brief Constructor. * - * \param[in] A %Matrix (expression), the base of the matrix power. - * \param[in] p scalar, the exponent of the matrix power. - * \prarm[in] b %Matrix (expression), the multiplier. + * \param[in] A the base of the matrix power. */ - MatrixPowerProductReturnValue(const Derived& A, RealScalar p, const RhsDerived& b) - : m_A(A), m_p(p), m_b(b) { } + explicit MatrixPower(const MatrixType& A); /** - * \brief Compute the expression. + * \brief Return the expression \f$ A^p \f$. * - * \param[out] result \f$ A^p b \f$ where \p A, \p p and \p bare as - * in the constructor. + * \param[in] p exponent, a real scalar. */ - template<typename ResultType> - inline void evalTo(ResultType& result) const - { - const MatrixType A = m_A; - const PlainObject b = m_b; - MatrixPower<MatrixType, PlainObject> mp(A, m_p, b); - mp.compute(result); - } + const MatrixPowerReturnValue<MatrixPower<MatrixType> > operator()(RealScalar p) + { return MatrixPowerReturnValue<MatrixPower<MatrixType> >(*this, p); } - Index rows() const { return m_b.rows(); } - Index cols() const { return m_b.cols(); } + /** + * \brief Compute the matrix power. + * + * \param[in] p exponent, a real scalar. + * \param[out] res \f$ A^p \f$ where A is specified in the + * constructor. + */ + void compute(MatrixType& res, RealScalar p); - private: - const Derived& m_A; - const RealScalar m_p; - const RhsDerived& m_b; - MatrixPowerProductReturnValue& operator=(const MatrixPowerProductReturnValue&); + /** + * \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[in] noalias + * \param[out] res \f$ A^p b \f$, where A is specified in the + * constructor. + */ + template<typename PlainObject, typename ResultType> + void compute(const PlainObject& b, ResultType& res, RealScalar p); + + Index rows() const { return m_A.rows(); } + Index cols() const { return m_A.cols(); } }; +template<typename MatrixType> +MatrixPower<MatrixType>::MatrixPower(const MatrixType& A) : + m_A(A), + m_init(false) +{ /* empty body */ } + +template<typename MatrixType> +void MatrixPower<MatrixType>::compute(MatrixType& res, RealScalar p) +{ + switch (m_A.cols()) { + case 0: + break; + case 1: + res(0,0) = std::pow(m_A(0,0), p); + break; + default: + RealScalar intpart; + RealScalar x = modfAndInit(p, &intpart); + res = MatrixType::Identity(m_A.rows(),m_A.cols()); + computeIntPower(res, intpart); + computeFracPower(res, x); + } +} + +template<typename MatrixType> +template<typename PlainObject, typename ResultType> +void MatrixPower<MatrixType>::compute(const PlainObject& b, ResultType& res, RealScalar p) +{ + switch (m_A.cols()) { + case 0: + break; + case 1: + res = std::pow(m_A(0,0), p) * b; + break; + default: + RealScalar intpart; + RealScalar x = modfAndInit(p, &intpart); + computeIntPower(b, res, intpart); + computeFracPower(res, x); + } +} + +template<typename MatrixType> +typename MatrixType::RealScalar MatrixPower<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart) +{ + static RealScalar maxAbsEival, minAbsEival; + *intpart = std::floor(x); + RealScalar res = x - *intpart; + + if (!m_init && res) { // !init && res + const ComplexSchur<MatrixType> schurOfA(m_A); + m_T = schurOfA.matrixT(); + m_U = schurOfA.matrixU(); + m_init = true; + + const Array<RealScalar,EIGEN_SIZE_MIN_PREFER_FIXED(Rows,Cols),1,ColMajor,EIGEN_SIZE_MIN_PREFER_FIXED(MaxRows,MaxCols)> + absTdiag = m_T.diagonal().array().abs(); + maxAbsEival = absTdiag.maxCoeff(); + minAbsEival = absTdiag.minCoeff(); + } + + if (res > RealScalar(0.5) && res > (1-res) * std::pow(maxAbsEival/minAbsEival, res)) { + --res; + ++*intpart; + } + return res; +} + +template<typename MatrixType> +template<typename PlainObject, typename ResultType> +void MatrixPower<MatrixType>::apply(const PlainObject& b, ResultType& res, bool& init) +{ + if (init) + res = m_tmp1 * res; + else { + init = true; + res.noalias() = m_tmp1 * b; + } +} + +template<typename MatrixType> +template<typename ResultType> +void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p) +{ + RealScalar pp = std::abs(p); + + if (p<0) m_tmp1 = m_A.inverse(); + else m_tmp1 = m_A; + + while (pp >= 1) { + if (std::fmod(pp, 2) >= 1) + res = m_tmp1 * res; + m_tmp1 *= m_tmp1; + pp /= 2; + } +} + +template<typename MatrixType> +template<typename PlainObject, typename ResultType> +void MatrixPower<MatrixType>::computeIntPower(const PlainObject& 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 cost = internal::binary_powering_cost(pp); + bool init = false; + + if (p==0) { + res = b; + return; + } + if (p<0) m_tmp1 = m_A.inverse(); + else m_tmp1 = m_A; + + while (b.cols()*pp > m_A.cols()*cost) { + if (std::fmod(pp, 2) >= 1) { + apply(b, res, init); + --cost; + } + m_tmp1 *= m_tmp1; + --cost; + pp /= 2; + } + for (; pp >= 1; --pp) + apply(b, res, init); + } +} + +template<typename MatrixType> +template<typename ResultType> +void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p) +{ + if (p) { + MatrixPowerTriangularAtomic<ComplexMatrix>(m_T).compute(m_fT, p); + internal::recompose_complex_schur<NumTraits<Scalar>::IsComplex>::run(m_tmp1, m_fT, m_U); + res = m_tmp1 * res; + } +} + /** * \ingroup MatrixFunctions_Module * @@ -525,11 +467,10 @@ class MatrixPowerProductReturnValue : public ReturnByValue<MatrixPowerProductRet template<typename Derived> class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Derived> > { - private: + public: typedef typename Derived::RealScalar RealScalar; typedef typename Derived::Index Index; - public: /** * \brief Constructor. * @@ -540,40 +481,14 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv : m_A(A), m_p(p) { } /** - * \brief Return the matrix power multiplied by %Matrix \p b. - * - * The %MatrixPower class can optimize \f$ A^p b \f$ computing, and - * this method provides an elegant way to call it. - * - * Unlike general matrix-matrix / matrix-vector product, this does - * \b NOT produce a temporary storage for the result. Therefore, - * the code below is \a already optimal: - * \code - * v = A.pow(p) * b; - * // v.noalias() = A.pow(p) * b; Won't compile! - * \endcode - * - * \param[in] b %Matrix (expression), the multiplier. - */ - template<typename RhsDerived> - const MatrixPowerProductReturnValue<Derived,RhsDerived> operator*(const MatrixBase<RhsDerived>& b) const - { return MatrixPowerProductReturnValue<Derived,RhsDerived>(m_A, m_p, b.derived()); } - - /** * \brief Compute the matrix power. * * \param[out] result \f$ A^p \f$ where \p A and \p p are as in the * constructor. */ template<typename ResultType> - inline void evalTo(ResultType& result) const - { - typedef typename Derived::PlainObject PlainObject; - const PlainObject A = m_A; - const PlainObject Identity = PlainObject::Identity(m_A.rows(), m_A.cols()); - MatrixPower<PlainObject> mp(A, m_p, Identity); - mp.compute(result); - } + inline void evalTo(ResultType& res) const + { MatrixPower<typename Derived::PlainObject>(m_A).compute(res, m_p); } Index rows() const { return m_A.rows(); } Index cols() const { return m_A.cols(); } @@ -584,18 +499,42 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&); }; +template<typename MatrixType> +class MatrixPowerReturnValue<MatrixPower<MatrixType> > +: public ReturnByValue<MatrixPowerReturnValue<MatrixPower<MatrixType> > > +{ + public: + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + + MatrixPowerReturnValue(MatrixPower<MatrixType>& ref, RealScalar p) + : m_pow(ref), m_p(p) { } + + template<typename ResultType> + 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: + MatrixPower<MatrixType>& m_pow; + const RealScalar m_p; + MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&); +}; + namespace internal { - template<typename Derived> - struct traits<MatrixPowerReturnValue<Derived> > - { - typedef typename Derived::PlainObject ReturnType; - }; - - template<typename Derived, typename RhsDerived> - struct traits<MatrixPowerProductReturnValue<Derived,RhsDerived> > - { - typedef typename RhsDerived::PlainObject ReturnType; - }; +template<typename Derived> +struct traits<MatrixPowerReturnValue<Derived> > +{ typedef typename Derived::PlainObject ReturnType; }; + +template<typename MatrixType> +struct traits<MatrixPowerReturnValue<MatrixPower<MatrixType> > > +{ typedef MatrixType ReturnType; }; + +template<typename Derived> +struct traits<MatrixPowerProductBase<Derived> > +{ typedef typename traits<Derived>::ReturnType ReturnType; }; } template<typename Derived> diff --git a/unsupported/doc/examples/MatrixPower.cpp b/unsupported/doc/examples/MatrixPower.cpp index 6ade0b8af..222452476 100644 --- a/unsupported/doc/examples/MatrixPower.cpp +++ b/unsupported/doc/examples/MatrixPower.cpp @@ -11,6 +11,6 @@ int main() sin(1), cos(1), 0, 0 , 0 , 1; std::cout << "The matrix A is:\n" << A << "\n\n" - << "The matrix power A^(pi/4) is:\n" << A.pow(pi/4) << std::endl; + "The matrix power A^(pi/4) is:\n" << A.pow(pi/4) << std::endl; return 0; } diff --git a/unsupported/doc/examples/MatrixPower_optimal.cpp b/unsupported/doc/examples/MatrixPower_optimal.cpp new file mode 100644 index 000000000..86470ba0a --- /dev/null +++ b/unsupported/doc/examples/MatrixPower_optimal.cpp @@ -0,0 +1,17 @@ +#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; +} diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index d891641f4..c5967d2eb 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -16,15 +16,17 @@ void test2dRotation(double tol) T angle, c, s; A << 0, 1, -1, 0; - for (int i = 0; i <= 20; i++) { + MatrixPower<Matrix<T,2,2> > Apow(A); + + for (int i=0; i<=20; ++i) { angle = pow(10, (i-10) / 5.); c = std::cos(angle); s = std::sin(angle); B << c, s, -s, c; - C = A.pow(std::ldexp(angle, 1) / M_PI); - std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C, B) << '\n'; - VERIFY(C.isApprox(B, T(tol))); + C = Apow(std::ldexp(angle,1) / M_PI); + std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n'; + VERIFY(C.isApprox(B, static_cast<T>(tol))); } } @@ -36,15 +38,17 @@ void test2dHyperbolicRotation(double tol) std::complex<T> ish(0, std::sinh(1)); A << ch, ish, -ish, ch; - for (int i = 0; i <= 20; i++) { - angle = std::ldexp(T(i-10), -1); + MatrixPower<Matrix<std::complex<T>,2,2> > Apow(A); + + for (int i=0; i<=20; ++i) { + angle = std::ldexp(static_cast<T>(i-10), -1); ch = std::cosh(angle); ish = std::complex<T>(0, std::sinh(angle)); B << ch, ish, -ish, ch; - C = A.pow(angle); - std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C, B) << '\n'; - VERIFY(C.isApprox(B, T(tol))); + C = Apow(angle); + std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n'; + VERIFY(C.isApprox(B, static_cast<T>(tol))); } } @@ -55,71 +59,37 @@ void testExponentLaws(const MatrixType& m, double tol) MatrixType m1, m2, m3, m4, m5; RealScalar x, y; - for (int i = 0; i < g_repeat; i++) { + for (int i=0; i<g_repeat; ++i) { generateTestMatrix<MatrixType>::run(m1, m.rows()); + MatrixPower<MatrixType> mpow(m1); + x = internal::random<RealScalar>(); y = internal::random<RealScalar>(); - m2 = m1.pow(x); - m3 = m1.pow(y); + m2 = mpow(x); + m3 = mpow(y); - m4 = m1.pow(x + y); + m4 = mpow(x+y); m5.noalias() = m2 * m3; + std::cout << "testExponentLaws: error powerm = " << relerr(m4, m5); - VERIFY(m4.isApprox(m5, RealScalar(tol))); + VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol))); - if (!NumTraits<typename MatrixType::Scalar>::IsComplex) { - m4 = m1.pow(x * y); - m5 = m2.pow(y); - std::cout << " " << relerr(m4, m5); - VERIFY(m4.isApprox(m5, RealScalar(tol))); - } + m4 = mpow(x*y); + m5 = m2.pow(y); + std::cout << " " << relerr(m4, m5); + VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol))); m4 = (std::abs(x) * m1).pow(y); m5 = std::pow(std::abs(x), y) * m3; std::cout << " " << relerr(m4, m5) << '\n'; - VERIFY(m4.isApprox(m5, RealScalar(tol))); - } -} - -template<typename MatrixType, typename VectorType> -void testMatrixVectorProduct(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<MatrixType>::run(m1, m.rows()); - v1 = VectorType::Random(v.rows(), v.cols()); - p = internal::random<RealScalar>(); - - v2.noalias() = m1.pow(p).eval() * v1; - v1 = m1.pow(p) * v1; - std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v1) << '\n'; - VERIFY(v2.isApprox(v1, RealScalar(tol))); - } -} - -template<typename MatrixType> -void testAliasing(const MatrixType& m) -{ - typedef typename MatrixType::RealScalar RealScalar; - MatrixType m1, m2; - RealScalar p; - - for (int i = 0; i < g_repeat; i++) { - generateTestMatrix<MatrixType>::run(m1, m.rows()); - p = internal::random<RealScalar>(); - - m2 = m1.pow(p); - m1 = m1.pow(p); - VERIFY(m1 == m2); + VERIFY(m4.isApprox(m5, static_cast<RealScalar>(tol))); } } void test_matrix_power() { + typedef Matrix<long double,Dynamic,Dynamic> MatrixXe; + CALL_SUBTEST_2(test2dRotation<double>(1e-13)); CALL_SUBTEST_1(test2dRotation<float>(2e-5)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64 CALL_SUBTEST_9(test2dRotation<long double>(1e-13)); @@ -135,20 +105,4 @@ void test_matrix_power() CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4)); CALL_SUBTEST_8(testExponentLaws(Matrix4f(), 1e-4)); CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4)); - - CALL_SUBTEST_2(testMatrixVectorProduct(Matrix2d(), Vector2d(), 1e-13)); - CALL_SUBTEST_7(testMatrixVectorProduct(Matrix<double,3,3,RowMajor>(), Vector3d(), 1e-13)); - CALL_SUBTEST_3(testMatrixVectorProduct(Matrix4cd(), Vector4cd(), 1e-13)); - CALL_SUBTEST_4(testMatrixVectorProduct(MatrixXd(8,8), MatrixXd(8,2), 1e-13)); - CALL_SUBTEST_1(testMatrixVectorProduct(Matrix2f(), Vector2f(), 1e-4)); - CALL_SUBTEST_5(testMatrixVectorProduct(Matrix3cf(), Vector3cf(), 1e-4)); - CALL_SUBTEST_8(testMatrixVectorProduct(Matrix4f(), Vector4f(), 1e-4)); - CALL_SUBTEST_6(testMatrixVectorProduct(MatrixXf(8,8), VectorXf(8), 1e-4)); - CALL_SUBTEST_9(testMatrixVectorProduct(Matrix<long double,Dynamic,Dynamic>(7,7), Matrix<long double,7,9>(), 1e-13)); - - CALL_SUBTEST_7(testAliasing(Matrix<double,3,3,RowMajor>())); - CALL_SUBTEST_3(testAliasing(Matrix4cd())); - CALL_SUBTEST_4(testAliasing(MatrixXd(8,8))); - CALL_SUBTEST_5(testAliasing(Matrix3cf())); - CALL_SUBTEST_6(testAliasing(MatrixXf(8,8))); } |