diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index 27bdcddd0..041d3b7ec 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -228,15 +228,15 @@ const MatrixPowerReturnValue MatrixBase::pow(con \endcode \param[in] M base of the matrix power, should be a square matrix. -\param[in] p exponent of the matrix power, should be an integer or -the same type as the real scalar in \p M. +\param[in] p exponent of the matrix power, should be real. The matrix power \f$ M^p \f$ is defined as \f$ \exp(p \log(M)) \f$, 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. +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(). diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h index 781d7bf93..7b40c0a43 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h @@ -51,7 +51,7 @@ private: void compute2x2(const MatrixType& A, MatrixType& result); void computeBig(const MatrixType& A, MatrixType& result); - static Scalar atanh(Scalar x); + static Scalar atanh2(Scalar y, Scalar x); int getPadeDegree(float normTminusI); int getPadeDegree(double normTminusI); int getPadeDegree(long double normTminusI); @@ -93,16 +93,18 @@ MatrixType MatrixLogarithmAtomic::compute(const MatrixType& A) return result; } -/** \brief Compute atanh (inverse hyperbolic tangent). */ +/** \brief Compute atanh (inverse hyperbolic tangent) for \f$ y / x \f$. */ template -typename MatrixType::Scalar MatrixLogarithmAtomic::atanh(typename MatrixType::Scalar x) +typename MatrixType::Scalar MatrixLogarithmAtomic::atanh2(Scalar y, Scalar x) { using std::abs; using std::sqrt; - if (abs(x) > sqrt(NumTraits::epsilon())) - return Scalar(0.5) * log((Scalar(1) + x) / (Scalar(1) - x)); + + Scalar z = y / x; + if (abs(z) > sqrt(NumTraits::epsilon())) + return Scalar(0.5) * log((x + y) / (x - y)); else - return x + x*x*x / Scalar(3); + return z + z*z*z / Scalar(3); } /** \brief Compute logarithm of 2x2 triangular matrix. */ @@ -128,8 +130,8 @@ void MatrixLogarithmAtomic::compute2x2(const MatrixType& A, MatrixTy } else { // computation in previous branch is inaccurate if A(1,1) \approx A(0,0) int unwindingNumber = static_cast(ceil((imag(logA11 - logA00) - M_PI) / (2*M_PI))); - Scalar z = (A(1,1) - A(0,0)) / (A(1,1) + A(0,0)); - result(0,1) = A(0,1) * (Scalar(2) * atanh(z) + Scalar(0,2*M_PI*unwindingNumber)) / (A(1,1) - A(0,0)); + Scalar y = A(1,1) - A(0,0), x = A(1,1) + A(0,0); + result(0,1) = A(0,1) * (Scalar(2) * atanh2(y,x) + Scalar(0,2*M_PI*unwindingNumber)) / y; } } diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 86ef24eac..2dff28080 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -21,16 +21,31 @@ namespace Eigen { * * \brief Class for computing matrix powers. * - * \tparam MatrixType type of the base, expected to be an instantiation + * \tparam MatrixType type of the base, expected to be an instantiation * of the Matrix class template. - * \tparam RealScalar type of the exponent, a real scalar. - * \tparam PlainObject type of the multiplier. - * \tparam IsInteger used internally to select correct specialization. + * \tparam ExponentType type of the exponent, a real scalar. + * \tparam PlainObject type of the multiplier. + * \tparam IsInteger used internally to select correct specialization. */ -template ::IsInteger> +template ::IsInteger> class MatrixPower { + private: + typedef internal::traits 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 ComplexScalar; + typedef typename MatrixType::Index Index; + typedef Matrix ComplexMatrix; + typedef Array ComplexArray; + public: /** * \brief Constructor. @@ -39,7 +54,7 @@ class MatrixPower * \param[in] p the exponent of the matrix power. * \param[in] b the multiplier. */ - MatrixPower(const MatrixType& A, const RealScalar& p, const PlainObject& b) : + MatrixPower(const MatrixType& A, RealScalar p, const PlainObject& b) : m_A(A), m_p(p), m_b(b), @@ -55,19 +70,6 @@ class MatrixPower template void compute(ResultType& result); private: - typedef internal::traits 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 std::complex ComplexScalar; - typedef typename MatrixType::Index Index; - typedef Matrix ComplexMatrix; - typedef Array ComplexArray; - /** * \brief Compute the matrix power. * @@ -112,8 +114,8 @@ class MatrixPower */ void getFractionalExponent(); - /** \brief Compute atanh (inverse hyperbolic tangent). */ - ComplexScalar atanh(const ComplexScalar& x); + /** \brief Compute atanh (inverse hyperbolic tangent) for \f$ y / x \f$. */ + ComplexScalar atanh2(const ComplexScalar& y, const ComplexScalar& x); /** \brief Compute power of 2x2 triangular matrix. */ void compute2x2(const RealScalar& p); @@ -237,9 +239,9 @@ class MatrixPower /******* Specialized for real exponents *******/ -template +template template -void MatrixPower::compute(ResultType& result) +void MatrixPower::compute(ResultType& result) { using std::floor; using std::pow; @@ -264,9 +266,9 @@ void MatrixPower::compute(ResultTyp } } -template +template template -void MatrixPower::computeIntPower(ResultType& result) +void MatrixPower::computeIntPower(ResultType& result) { if (m_dimb > m_dimA) { MatrixType tmp = MatrixType::Identity(m_A.rows(), m_A.cols()); @@ -278,9 +280,9 @@ void MatrixPower::computeIntPower(R } } -template +template template -void MatrixPower::computeChainProduct(ResultType& result) +void MatrixPower::computeChainProduct(ResultType& result) { using std::frexp; using std::ldexp; @@ -312,8 +314,8 @@ void MatrixPower::computeChainProdu result = m_tmp * result; } -template -int MatrixPower::computeCost(RealScalar p) +template +int MatrixPower::computeCost(RealScalar p) { using std::frexp; using std::ldexp; @@ -326,25 +328,25 @@ int MatrixPower::computeCost(RealSc return cost; } -template +template template -void MatrixPower::partialPivLuSolve(ResultType& result, RealScalar p) +void MatrixPower::partialPivLuSolve(ResultType& result, RealScalar p) { const PartialPivLU Asolver(m_A); for (; p >= RealScalar(1); p--) result = Asolver.solve(result); } -template -void MatrixPower::computeSchurDecomposition() +template +void MatrixPower::computeSchurDecomposition() { const ComplexSchur schurOfA(m_A); m_T = schurOfA.matrixT(); m_U = schurOfA.matrixU(); } -template -void MatrixPower::getFractionalExponent() +template +void MatrixPower::getFractionalExponent() { using std::pow; @@ -373,21 +375,24 @@ void MatrixPower::getFractionalExpo } } -template -std::complex MatrixPower::atanh(const ComplexScalar& x) +template +std::complex +MatrixPower::atanh2(const ComplexScalar& y, const ComplexScalar& x) { using std::abs; using std::log; using std::sqrt; - if (abs(x) > sqrt(NumTraits::epsilon())) - return RealScalar(0.5) * log((RealScalar(1) + x) / (RealScalar(1) - x)); + const ComplexScalar z = y / x; + + if (abs(z) > sqrt(NumTraits::epsilon())) + return RealScalar(0.5) * log((x + y) / (x - y)); else - return x + x*x*x / RealScalar(3); + return z + z*z*z / RealScalar(3); } -template -void MatrixPower::compute2x2(const RealScalar& p) +template +void MatrixPower::compute2x2(const RealScalar& p) { using std::abs; using std::ceil; @@ -414,15 +419,15 @@ void MatrixPower::compute2x2(const else { // computation in previous branch is inaccurate if abs(m_T(j,j)) \approx abs(m_T(i,i)) unwindingNumber = static_cast(ceil((imag(m_logTdiag[j] - m_logTdiag[i]) - M_PI) / (2 * M_PI))); - w = atanh((m_T(j,j) - m_T(i,i)) / (m_T(j,j) + m_T(i,i))) + ComplexScalar(0, M_PI * unwindingNumber); + w = 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 -void MatrixPower::computeBig() +template +void MatrixPower::computeBig() { using std::ldexp; const int digits = std::numeric_limits::digits; @@ -458,8 +463,8 @@ void MatrixPower::computeBig() compute2x2(m_pfrac); } -template -inline int MatrixPower::getPadeDegree(float normIminusT) +template +inline int MatrixPower::getPadeDegree(float normIminusT) { const float maxNormForPade[] = { 2.7996156e-1f /* degree = 3 */ , 4.3268868e-1f }; int degree = 3; @@ -469,8 +474,8 @@ inline int MatrixPower::getPadeDegr return degree; } -template -inline int MatrixPower::getPadeDegree(double normIminusT) +template +inline int MatrixPower::getPadeDegree(double normIminusT) { const double maxNormForPade[] = { 1.882832775783710e-2 /* degree = 3 */ , 6.036100693089536e-2, 1.239372725584857e-1, 1.998030690604104e-1, 2.787629930861592e-1 }; @@ -481,8 +486,8 @@ inline int MatrixPower::getPadeDegr return degree; } -template -inline int MatrixPower::getPadeDegree(long double normIminusT) +template +inline int MatrixPower::getPadeDegree(long double normIminusT) { #if LDBL_MANT_DIG == 53 const int maxPadeDegree = 7; @@ -514,8 +519,8 @@ inline int MatrixPower::getPadeDegr break; return degree; } -template -void MatrixPower::computePade(const int& degree, const ComplexMatrix& IminusT) +template +void MatrixPower::computePade(const int& degree, const ComplexMatrix& IminusT) { int i = degree << 1; m_fT = coeff(i) * IminusT; @@ -526,8 +531,8 @@ void MatrixPower::computePade(const m_fT += ComplexMatrix::Identity(m_A.rows(), m_A.cols()); } -template -inline RealScalar MatrixPower::coeff(const int& i) +template +inline typename MatrixType::RealScalar MatrixPower::coeff(const int& i) { if (i == 1) return -m_pfrac; @@ -537,13 +542,13 @@ inline RealScalar MatrixPower::coef return (m_pfrac - RealScalar(i >> 1)) / RealScalar(i-1 << 1); } -template -void MatrixPower::computeTmp(RealScalar) +template +void MatrixPower::computeTmp(RealScalar) { m_tmp = (m_U * m_fT * m_U.adjoint()).real(); } -template -void MatrixPower::computeTmp(ComplexScalar) -{ m_tmp = (m_U * m_fT * m_U.adjoint()).eval(); } +template +void MatrixPower::computeTmp(ComplexScalar) +{ m_tmp = m_U * m_fT * m_U.adjoint(); } /******* Specialized for integral exponents *******/ diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index b34b151b1..ff0137ec6 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -33,6 +33,7 @@ endif() ei_add_test(matrix_exponential) ei_add_test(matrix_function) +ei_add_test(matrix_power) ei_add_test(matrix_square_root) ei_add_test(alignedvector3) ei_add_test(FFT) diff --git a/unsupported/test/matrix_exponential.cpp b/unsupported/test/matrix_exponential.cpp index 695472f91..50dec083d 100644 --- a/unsupported/test/matrix_exponential.cpp +++ b/unsupported/test/matrix_exponential.cpp @@ -7,8 +7,7 @@ // 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/. -#include "main.h" -#include +#include "matrix_functions.h" double binom(int n, int k) { @@ -18,12 +17,6 @@ double binom(int n, int k) return res; } -template -double relerr(const MatrixBase& A, const MatrixBase& B) -{ - return std::sqrt((A - B).cwiseAbs2().sum() / (std::min)(A.cwiseAbs2().sum(), B.cwiseAbs2().sum())); -} - template T expfn(T x, int) { @@ -109,8 +102,7 @@ void randomTest(const MatrixType& m, double tol) */ typename MatrixType::Index rows = m.rows(); typename MatrixType::Index cols = m.cols(); - MatrixType m1(rows, cols), m2(rows, cols), m3(rows, cols), - identity = MatrixType::Identity(rows, rows); + MatrixType m1(rows, cols), m2(rows, cols), identity = MatrixType::Identity(rows, cols); typedef typename NumTraits::Scalar>::Real RealScalar; diff --git a/unsupported/test/matrix_functions.h b/unsupported/test/matrix_functions.h new file mode 100644 index 000000000..5817caef6 --- /dev/null +++ b/unsupported/test/matrix_functions.h @@ -0,0 +1,47 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009-2011 Jitse Niesen +// +// 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/. + +#include "main.h" +#include + +template ::Scalar>::IsComplex> +struct generateTestMatrix; + +// for real matrices, make sure none of the eigenvalues are negative +template +struct generateTestMatrix +{ + static void run(MatrixType& result, typename MatrixType::Index size) + { + MatrixType mat = MatrixType::Random(size, size); + EigenSolver es(mat); + typename EigenSolver::EigenvalueType eivals = es.eigenvalues(); + for (typename MatrixType::Index i = 0; i < size; ++i) { + if (eivals(i).imag() == 0 && eivals(i).real() < 0) + eivals(i) = -eivals(i); + } + result = (es.eigenvectors() * eivals.asDiagonal() * es.eigenvectors().inverse()).real(); + } +}; + +// for complex matrices, any matrix is fine +template +struct generateTestMatrix +{ + static void run(MatrixType& result, typename MatrixType::Index size) + { + result = MatrixType::Random(size, size); + } +}; + +template +double relerr(const MatrixBase& A, const MatrixBase& B) +{ + return std::sqrt((A - B).cwiseAbs2().sum() / (std::min)(A.cwiseAbs2().sum(), B.cwiseAbs2().sum())); +} diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp new file mode 100644 index 000000000..f3ef57157 --- /dev/null +++ b/unsupported/test/matrix_power.cpp @@ -0,0 +1,104 @@ +// 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/. + +#include "matrix_functions.h" + +template +void test2dRotation(double tol) +{ + Matrix A, B, C; + T angle, c, s; + + A << 0, 1, -1, 0; + 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))); + } +} + +template +void test2dHyperbolicRotation(double tol) +{ + Matrix,2,2> A, B, C; + T angle, ch = std::cosh(1); + std::complex ish(0, std::sinh(1)); + + A << ch, ish, -ish, ch; + for (int i = 0; i <= 20; i++) { + angle = std::ldexp(T(i-10), -1); + ch = std::cosh(angle); + ish = std::complex(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))); + } +} + +template +void testExponentLaws(const MatrixType& m, double tol) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + + typename MatrixType::Index rows = m.rows(); + typename MatrixType::Index cols = m.cols(); + MatrixType m1, m1x, m1y, m2, m3; + RealScalar x = internal::random(), y = internal::random(); + double err[3]; + + for(int i = 0; i < g_repeat; i++) { + generateTestMatrix::run(m1, m.rows()); + m1x = m1.pow(x); + m1y = m1.pow(y); + + m2 = m1.pow(x + y); + m3 = m1x * m1y; + err[0] = relerr(m2, m3); + VERIFY(m2.isApprox(m3, static_cast(tol))); + + m2 = m1.pow(x * y); + m3 = m1x.pow(y); + err[1] = relerr(m2, m3); + VERIFY(m2.isApprox(m3, static_cast(tol))); + + m2 = (std::abs(x) * m1).pow(y); + m3 = std::pow(std::abs(x), y) * m1y; + err[2] = relerr(m2, m3); + VERIFY(m2.isApprox(m3, static_cast(tol))); + + std::cout << "testExponentLaws: error powerm = " << err[0] << " " << err[1] << " " << err[2] << "\n"; + } +} + +void test_matrix_power() +{ + CALL_SUBTEST_2(test2dRotation(1e-13)); + CALL_SUBTEST_1(test2dRotation(2e-5)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64 + CALL_SUBTEST_8(test2dRotation(1e-13)); + CALL_SUBTEST_2(test2dHyperbolicRotation(1e-14)); + CALL_SUBTEST_1(test2dHyperbolicRotation(1e-5)); + CALL_SUBTEST_8(test2dHyperbolicRotation(1e-14)); + CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13)); + CALL_SUBTEST_7(testExponentLaws(Matrix(), 1e-13)); + CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13)); + CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 1e-13)); + CALL_SUBTEST_1(testExponentLaws(Matrix2f(), 1e-4)); + CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4)); + CALL_SUBTEST_1(testExponentLaws(Matrix4f(), 1e-4)); + CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4)); + CALL_SUBTEST_9(testExponentLaws(Matrix(7,7), 1e-13)); +} diff --git a/unsupported/test/matrix_square_root.cpp b/unsupported/test/matrix_square_root.cpp index 508619a7a..ea541e1ea 100644 --- a/unsupported/test/matrix_square_root.cpp +++ b/unsupported/test/matrix_square_root.cpp @@ -7,38 +7,7 @@ // 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/. -#include "main.h" -#include - -template ::Scalar>::IsComplex> -struct generateTestMatrix; - -// for real matrices, make sure none of the eigenvalues are negative -template -struct generateTestMatrix -{ - static void run(MatrixType& result, typename MatrixType::Index size) - { - MatrixType mat = MatrixType::Random(size, size); - EigenSolver es(mat); - typename EigenSolver::EigenvalueType eivals = es.eigenvalues(); - for (typename MatrixType::Index i = 0; i < size; ++i) { - if (eivals(i).imag() == 0 && eivals(i).real() < 0) - eivals(i) = -eivals(i); - } - result = (es.eigenvectors() * eivals.asDiagonal() * es.eigenvectors().inverse()).real(); - } -}; - -// for complex matrices, any matrix is fine -template -struct generateTestMatrix -{ - static void run(MatrixType& result, typename MatrixType::Index size) - { - result = MatrixType::Random(size, size); - } -}; +#include "matrix_functions.h" template void testMatrixSqrt(const MatrixType& m)