diff --git a/Eigen/Eigenvalues b/Eigen/Eigenvalues index 986f31196..f22a3bc30 100644 --- a/Eigen/Eigenvalues +++ b/Eigen/Eigenvalues @@ -42,6 +42,7 @@ namespace Eigen { #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 c41bbefaa..9e2afe7e4 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -136,8 +136,8 @@ template class MatrixBase CwiseUnaryOp, Eigen::Transpose >, Transpose >::ret AdjointReturnType; - /** \internal the return type of MatrixBase::eigenvalues() */ - typedef Matrix::Scalar>::Real, ei_traits::ColsAtCompileTime, 1> EigenvaluesReturnType; + /** \internal Return type of eigenvalues() */ + typedef Matrix, ei_traits::ColsAtCompileTime, 1> EigenvaluesReturnType; /** \internal the return type of identity */ typedef CwiseNullaryOp,Derived> IdentityReturnType; /** \internal the return type of unit vectors */ diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index f6ae05511..277108dd4 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -153,6 +153,16 @@ template class SelfAdjointView const LLT llt() const; const LDLT ldlt() const; +/////////// Eigenvalue module /////////// + + /** Real part of #Scalar */ + typedef typename NumTraits::Real RealScalar; + /** Return type of eigenvalues() */ + typedef Matrix::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 000000000..7b04e6ba7 --- /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 +// Copyright (C) 2010 Jitse Niesen +// +// 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 . + +#ifndef EIGEN_MATRIXBASEEIGENVALUES_H +#define EIGEN_MATRIXBASEEIGENVALUES_H + + + +template +struct ei_eigenvalues_selector +{ + // this is the implementation for the case IsComplex = true + static inline typename MatrixBase::EigenvaluesReturnType const + run(const MatrixBase& m) + { + typedef typename Derived::PlainObject PlainObject; + PlainObject m_eval(m); + return ComplexEigenSolver(m_eval).eigenvalues(); + } +}; + +template +struct ei_eigenvalues_selector +{ + static inline typename MatrixBase::EigenvaluesReturnType const + run(const MatrixBase& m) + { + typedef typename Derived::PlainObject PlainObject; + PlainObject m_eval(m); + return EigenSolver(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 +inline typename MatrixBase::EigenvaluesReturnType +MatrixBase::eigenvalues() const +{ + typedef typename ei_traits::Scalar Scalar; + return ei_eigenvalues_selector::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 +inline typename SelfAdjointView::EigenvaluesReturnType +SelfAdjointView::eigenvalues() const +{ + typedef typename SelfAdjointView::PlainObject PlainObject; + PlainObject thisAsMatrix(*this); + return SelfAdjointEigenSolver(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 +inline typename MatrixBase::RealScalar +MatrixBase::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() + .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 +inline typename SelfAdjointView::RealScalar +SelfAdjointView::operatorNorm() const +{ + return eigenvalues().cwiseAbs().maxCoeff(); +} + +#endif diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 1abbed97b..76343640d 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -481,59 +481,6 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors #endif // EIGEN_HIDE_HEAVY_CODE -/** \eigenvalues_module - * - * \returns a vector listing the eigenvalues of this matrix. - */ -template -inline Matrix::Scalar>::Real, ei_traits::ColsAtCompileTime, 1> -MatrixBase::eigenvalues() const -{ - ei_assert(Flags&SelfAdjoint); - return SelfAdjointEigenSolver(eval(),false).eigenvalues(); -} - -template -struct ei_operatorNorm_selector -{ - static inline typename NumTraits::Scalar>::Real - operatorNorm(const MatrixBase& 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 struct ei_operatorNorm_selector -{ - static inline typename NumTraits::Scalar>::Real - operatorNorm(const MatrixBase& 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() - .eigenvalues() - .maxCoeff() - ); - } -}; - -/** \eigenvalues_module - * - * \returns the matrix norm of this matrix. - */ -template -inline typename NumTraits::Scalar>::Real -MatrixBase::operatorNorm() const -{ - return ei_operatorNorm_selector - ::operatorNorm(derived()); -} - #ifndef EIGEN_EXTERN_INSTANTIATIONS template 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 000000000..039f88701 --- /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 000000000..355246f0d --- /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 000000000..be1986778 --- /dev/null +++ b/doc/snippets/SelfAdjointView_eigenvalues.cpp @@ -0,0 +1,3 @@ +MatrixXd ones = MatrixXd::Ones(3,3); +VectorXd eivals = ones.selfadjointView().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 000000000..f380f5594 --- /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().operatorNorm() << endl; diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp index b3d9ac24b..5c5d7b38f 100644 --- a/test/eigensolver_complex.cpp +++ b/test/eigensolver_complex.cpp @@ -26,6 +26,21 @@ #include #include +/* 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 +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 void eigensolver(const MatrixType& m) { /* this test covers the following files: @@ -48,11 +63,17 @@ template void eigensolver(const MatrixType& m) ComplexEigenSolver 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 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 f24c3b4ed..d70f37ea4 100644 --- a/test/eigensolver_generic.cpp +++ b/test/eigensolver_generic.cpp @@ -58,7 +58,10 @@ template void eigensolver(const MatrixType& m) VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix()); VERIFY_IS_APPROX(a.template cast() * 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 void eigensolver_verify_assert() diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 70b3e6791..25ef280a1 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -103,6 +103,7 @@ template void selfadjointeigensolver(const MatrixType& m) VERIFY((symmA * eiSymm.eigenvectors()).isApprox( eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps)); + VERIFY_IS_APPROX(symmA.template selfadjointView().eigenvalues(), eiSymm.eigenvalues()); // generalized eigen problem Ax = lBx VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox( @@ -111,6 +112,9 @@ template void selfadjointeigensolver(const MatrixType& m) 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().operatorNorm(), RealScalar(1)); } void test_eigensolver_selfadjoint()