diff options
author | Jitse Niesen <jitse@maths.leeds.ac.uk> | 2010-05-24 17:43:50 +0100 |
---|---|---|
committer | Jitse Niesen <jitse@maths.leeds.ac.uk> | 2010-05-24 17:43:50 +0100 |
commit | e7d809d4349fd4048777be71f1c803d0b13f8fe8 (patch) | |
tree | 11c1ef9908d0756958fde6c29de7f1f16d5dc639 /Eigen | |
parent | 8a3f552e39d3fee3ada1cfc1eb75b179c77f2a78 (diff) |
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
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/Eigenvalues | 1 | ||||
-rw-r--r-- | Eigen/src/Core/MatrixBase.h | 4 | ||||
-rw-r--r-- | Eigen/src/Core/SelfAdjointView.h | 10 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h | 168 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 53 |
5 files changed, 181 insertions, 55 deletions
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<typename Derived> class MatrixBase 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 f6ae05511..277108dd4 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -153,6 +153,16 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView 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 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 <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 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<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) |