Update eigenvalues() and operatorNorm() methods in MatrixBase. * use SelfAdjointView instead of Eigen2's SelfAdjoint flag. * add tests and documentation. * allow eigenvalues() for non-selfadjoint matrices. * they no longer depend only on SelfAdjointEigenSolver, so move them to a separate file
diff --git a/Eigen/Eigenvalues b/Eigen/Eigenvalues index 986f311..f22a3bc 100644 --- a/Eigen/Eigenvalues +++ b/Eigen/Eigenvalues
@@ -42,6 +42,7 @@ #include "src/Eigenvalues/HessenbergDecomposition.h" #include "src/Eigenvalues/ComplexSchur.h" #include "src/Eigenvalues/ComplexEigenSolver.h" +#include "src/Eigenvalues/MatrixBaseEigenvalues.h" // declare all classes for a given matrix type #define EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index c41bbef..9e2afe7 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h
@@ -136,8 +136,8 @@ CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, Eigen::Transpose<Derived> >, Transpose<Derived> >::ret AdjointReturnType; - /** \internal the return type of MatrixBase::eigenvalues() */ - typedef Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1> EigenvaluesReturnType; + /** \internal Return type of eigenvalues() */ + typedef Matrix<std::complex<RealScalar>, ei_traits<Derived>::ColsAtCompileTime, 1> EigenvaluesReturnType; /** \internal the return type of identity */ typedef CwiseNullaryOp<ei_scalar_identity_op<Scalar>,Derived> IdentityReturnType; /** \internal the return type of unit vectors */
diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index f6ae055..277108d 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h
@@ -153,6 +153,16 @@ const LLT<PlainObject, UpLo> llt() const; const LDLT<PlainObject> ldlt() const; +/////////// Eigenvalue module /////////// + + /** Real part of #Scalar */ + typedef typename NumTraits<Scalar>::Real RealScalar; + /** Return type of eigenvalues() */ + typedef Matrix<RealScalar, ei_traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType; + + EigenvaluesReturnType eigenvalues() const; + RealScalar operatorNorm() const; + protected: const typename MatrixType::Nested m_matrix; };
diff --git a/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h b/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h new file mode 100644 index 0000000..7b04e6b --- /dev/null +++ b/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h
@@ -0,0 +1,168 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_MATRIXBASEEIGENVALUES_H +#define EIGEN_MATRIXBASEEIGENVALUES_H + + + +template<typename Derived, bool IsComplex> +struct ei_eigenvalues_selector +{ + // this is the implementation for the case IsComplex = true + static inline typename MatrixBase<Derived>::EigenvaluesReturnType const + run(const MatrixBase<Derived>& m) + { + typedef typename Derived::PlainObject PlainObject; + PlainObject m_eval(m); + return ComplexEigenSolver<PlainObject>(m_eval).eigenvalues(); + } +}; + +template<typename Derived> +struct ei_eigenvalues_selector<Derived, false> +{ + static inline typename MatrixBase<Derived>::EigenvaluesReturnType const + run(const MatrixBase<Derived>& m) + { + typedef typename Derived::PlainObject PlainObject; + PlainObject m_eval(m); + return EigenSolver<PlainObject>(m_eval).eigenvalues(); + } +}; + +/** \brief Computes the eigenvalues of a matrix + * \returns Column vector containing the eigenvalues. + * + * \eigenvalues_module + * This function computes the eigenvalues with the help of the EigenSolver + * class (for real matrices) or the ComplexEigenSolver class (for complex + * matrices). + * + * The eigenvalues are repeated according to their algebraic multiplicity, + * so there are as many eigenvalues as rows in the matrix. + * + * The SelfAdjointView class provides a better algorithm for selfadjoint + * matrices. + * + * Example: \include MatrixBase_eigenvalues.cpp + * Output: \verbinclude MatrixBase_eigenvalues.out + * + * \sa EigenSolver::eigenvalues(), ComplexEigenSolver::eigenvalues(), + * SelfAdjointView::eigenvalues() + */ +template<typename Derived> +inline typename MatrixBase<Derived>::EigenvaluesReturnType +MatrixBase<Derived>::eigenvalues() const +{ + typedef typename ei_traits<Derived>::Scalar Scalar; + return ei_eigenvalues_selector<Derived, NumTraits<Scalar>::IsComplex>::run(derived()); +} + +/** \brief Computes the eigenvalues of a matrix + * \returns Column vector containing the eigenvalues. + * + * \eigenvalues_module + * This function computes the eigenvalues with the help of the + * SelfAdjointEigenSolver class. The eigenvalues are repeated according to + * their algebraic multiplicity, so there are as many eigenvalues as rows in + * the matrix. + * + * Example: \include SelfAdjointView_eigenvalues.cpp + * Output: \verbinclude SelfAdjointView_eigenvalues.out + * + * \sa SelfAdjointEigenSolver::eigenvalues(), MatrixBase::eigenvalues() + */ +template<typename MatrixType, unsigned int UpLo> +inline typename SelfAdjointView<MatrixType, UpLo>::EigenvaluesReturnType +SelfAdjointView<MatrixType, UpLo>::eigenvalues() const +{ + typedef typename SelfAdjointView<MatrixType, UpLo>::PlainObject PlainObject; + PlainObject thisAsMatrix(*this); + return SelfAdjointEigenSolver<PlainObject>(thisAsMatrix).eigenvalues(); +} + + + +/** \brief Computes the L2 operator norm + * \returns Operator norm of the matrix. + * + * \eigenvalues_module + * This function computes the L2 operator norm of a matrix, which is also + * known as the spectral norm. The norm of a matrix \f$ A \f$ is defined to be + * \f[ \|A\|_2 = \max_x \frac{\|Ax\|_2}{\|x\|_2} \f] + * where the maximum is over all vectors and the norm on the right is the + * Euclidean vector norm. The norm equals the largest singular value, which is + * the square root of the largest eigenvalue of the positive semi-definite + * matrix \f$ A^*A \f$. + * + * The current implementation uses the eigenvalues of \f$ A^*A \f$, as computed + * by SelfAdjointView::eigenvalues(), to compute the operator norm of a + * matrix. The SelfAdjointView class provides a better algorithm for + * selfadjoint matrices. + * + * Example: \include MatrixBase_operatorNorm.cpp + * Output: \verbinclude MatrixBase_operatorNorm.out + * + * \sa SelfAdjointView::eigenvalues(), SelfAdjointView::operatorNorm() + */ +template<typename Derived> +inline typename MatrixBase<Derived>::RealScalar +MatrixBase<Derived>::operatorNorm() const +{ + typename Derived::PlainObject m_eval(derived()); + // FIXME if it is really guaranteed that the eigenvalues are already sorted, + // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough. + return ei_sqrt((m_eval*m_eval.adjoint()) + .eval() + .template selfadjointView<Lower>() + .eigenvalues() + .maxCoeff() + ); +} + +/** \brief Computes the L2 operator norm + * \returns Operator norm of the matrix. + * + * \eigenvalues_module + * This function computes the L2 operator norm of a self-adjoint matrix. For a + * self-adjoint matrix, the operator norm is the largest eigenvalue. + * + * The current implementation uses the eigenvalues of the matrix, as computed + * by eigenvalues(), to compute the operator norm of the matrix. + * + * Example: \include SelfAdjointView_operatorNorm.cpp + * Output: \verbinclude SelfAdjointView_operatorNorm.out + * + * \sa eigenvalues(), MatrixBase::operatorNorm() + */ +template<typename MatrixType, unsigned int UpLo> +inline typename SelfAdjointView<MatrixType, UpLo>::RealScalar +SelfAdjointView<MatrixType, UpLo>::operatorNorm() const +{ + return eigenvalues().cwiseAbs().maxCoeff(); +} + +#endif
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 1abbed9..7634364 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -481,59 +481,6 @@ #endif // EIGEN_HIDE_HEAVY_CODE -/** \eigenvalues_module - * - * \returns a vector listing the eigenvalues of this matrix. - */ -template<typename Derived> -inline Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1> -MatrixBase<Derived>::eigenvalues() const -{ - ei_assert(Flags&SelfAdjoint); - return SelfAdjointEigenSolver<typename Derived::PlainObject>(eval(),false).eigenvalues(); -} - -template<typename Derived, bool IsSelfAdjoint> -struct ei_operatorNorm_selector -{ - static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real - operatorNorm(const MatrixBase<Derived>& m) - { - // FIXME if it is really guaranteed that the eigenvalues are already sorted, - // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough. - return m.eigenvalues().cwiseAbs().maxCoeff(); - } -}; - -template<typename Derived> struct ei_operatorNorm_selector<Derived, false> -{ - static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real - operatorNorm(const MatrixBase<Derived>& m) - { - typename Derived::PlainObject m_eval(m); - // FIXME if it is really guaranteed that the eigenvalues are already sorted, - // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough. - return ei_sqrt( - (m_eval*m_eval.adjoint()) - .template marked<SelfAdjoint>() - .eigenvalues() - .maxCoeff() - ); - } -}; - -/** \eigenvalues_module - * - * \returns the matrix norm of this matrix. - */ -template<typename Derived> -inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real -MatrixBase<Derived>::operatorNorm() const -{ - return ei_operatorNorm_selector<Derived, Flags&SelfAdjoint> - ::operatorNorm(derived()); -} - #ifndef EIGEN_EXTERN_INSTANTIATIONS template<typename RealScalar, typename Scalar> static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n)
diff --git a/doc/snippets/MatrixBase_eigenvalues.cpp b/doc/snippets/MatrixBase_eigenvalues.cpp new file mode 100644 index 0000000..039f887 --- /dev/null +++ b/doc/snippets/MatrixBase_eigenvalues.cpp
@@ -0,0 +1,3 @@ +MatrixXd ones = MatrixXd::Ones(3,3); +VectorXcd eivals = ones.eigenvalues(); +cout << "The eigenvalues of the 3x3 matrix of ones are:" << endl << eivals << endl;
diff --git a/doc/snippets/MatrixBase_operatorNorm.cpp b/doc/snippets/MatrixBase_operatorNorm.cpp new file mode 100644 index 0000000..355246f --- /dev/null +++ b/doc/snippets/MatrixBase_operatorNorm.cpp
@@ -0,0 +1,3 @@ +MatrixXd ones = MatrixXd::Ones(3,3); +cout << "The operator norm of the 3x3 matrix of ones is " + << ones.operatorNorm() << endl;
diff --git a/doc/snippets/SelfAdjointView_eigenvalues.cpp b/doc/snippets/SelfAdjointView_eigenvalues.cpp new file mode 100644 index 0000000..be19867 --- /dev/null +++ b/doc/snippets/SelfAdjointView_eigenvalues.cpp
@@ -0,0 +1,3 @@ +MatrixXd ones = MatrixXd::Ones(3,3); +VectorXd eivals = ones.selfadjointView<Lower>().eigenvalues(); +cout << "The eigenvalues of the 3x3 matrix of ones are:" << endl << eivals << endl;
diff --git a/doc/snippets/SelfAdjointView_operatorNorm.cpp b/doc/snippets/SelfAdjointView_operatorNorm.cpp new file mode 100644 index 0000000..f380f55 --- /dev/null +++ b/doc/snippets/SelfAdjointView_operatorNorm.cpp
@@ -0,0 +1,3 @@ +MatrixXd ones = MatrixXd::Ones(3,3); +cout << "The operator norm of the 3x3 matrix of ones is " + << ones.selfadjointView<Lower>().operatorNorm() << endl;
diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp index b3d9ac2..5c5d7b3 100644 --- a/test/eigensolver_complex.cpp +++ b/test/eigensolver_complex.cpp
@@ -26,6 +26,21 @@ #include <Eigen/Eigenvalues> #include <Eigen/LU> +/* Check that two column vectors are approximately equal upto permutations, + by checking that the k-th power sums are equal for k = 1, ..., vec1.rows() */ +template<typename VectorType> +void verify_is_approx_upto_permutation(const VectorType& vec1, const VectorType& vec2) +{ + VERIFY(vec1.cols() == 1); + VERIFY(vec2.cols() == 1); + VERIFY(vec1.rows() == vec2.rows()); + for (int k = 1; k <= vec1.rows(); ++k) + { + VERIFY_IS_APPROX(vec1.array().pow(k).sum(), vec2.array().pow(k).sum()); + } +} + + template<typename MatrixType> void eigensolver(const MatrixType& m) { /* this test covers the following files: @@ -48,11 +63,17 @@ ComplexEigenSolver<MatrixType> ei1(a); VERIFY_IS_APPROX(a * ei1.eigenvectors(), ei1.eigenvectors() * ei1.eigenvalues().asDiagonal()); - + // Note: If MatrixType is real then a.eigenvalues() uses EigenSolver and thus + // another algorithm so results may differ slightly + verify_is_approx_upto_permutation(a.eigenvalues(), ei1.eigenvalues()); + // Regression test for issue #66 MatrixType z = MatrixType::Zero(rows,cols); ComplexEigenSolver<MatrixType> eiz(z); VERIFY((eiz.eigenvalues().cwiseEqual(0)).all()); + + MatrixType id = MatrixType::Identity(rows, cols); + VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1)); } void test_eigensolver_complex()
diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp index f24c3b4..d70f37e 100644 --- a/test/eigensolver_generic.cpp +++ b/test/eigensolver_generic.cpp
@@ -58,7 +58,10 @@ VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix()); VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(), ei1.eigenvectors() * ei1.eigenvalues().asDiagonal()); + VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues()); + MatrixType id = MatrixType::Identity(rows, cols); + VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1)); } template<typename MatrixType> void eigensolver_verify_assert()
diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 70b3e67..25ef280 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp
@@ -103,6 +103,7 @@ VERIFY((symmA * eiSymm.eigenvectors()).isApprox( eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps)); + VERIFY_IS_APPROX(symmA.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues()); // generalized eigen problem Ax = lBx VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox( @@ -111,6 +112,9 @@ MatrixType sqrtSymmA = eiSymm.operatorSqrt(); VERIFY_IS_APPROX(symmA, sqrtSymmA*sqrtSymmA); VERIFY_IS_APPROX(sqrtSymmA, symmA*eiSymm.operatorInverseSqrt()); + + MatrixType id = MatrixType::Identity(rows, cols); + VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1)); } void test_eigensolver_selfadjoint()