MatrixBaseEigenvalues.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_MATRIXBASEEIGENVALUES_H
00012 #define EIGEN_MATRIXBASEEIGENVALUES_H
00013 
00014 namespace Eigen { 
00015 
00016 namespace internal {
00017 
00018 template<typename Derived, bool IsComplex>
00019 struct eigenvalues_selector
00020 {
00021   // this is the implementation for the case IsComplex = true
00022   static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
00023   run(const MatrixBase<Derived>& m)
00024   {
00025     typedef typename Derived::PlainObject PlainObject;
00026     PlainObject m_eval(m);
00027     return ComplexEigenSolver<PlainObject>(m_eval, false).eigenvalues();
00028   }
00029 };
00030 
00031 template<typename Derived>
00032 struct eigenvalues_selector<Derived, false>
00033 {
00034   static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
00035   run(const MatrixBase<Derived>& m)
00036   {
00037     typedef typename Derived::PlainObject PlainObject;
00038     PlainObject m_eval(m);
00039     return EigenSolver<PlainObject>(m_eval, false).eigenvalues();
00040   }
00041 };
00042 
00043 } // end namespace internal
00044 
00045 /** \brief Computes the eigenvalues of a matrix 
00046   * \returns Column vector containing the eigenvalues.
00047   *
00048   * \eigenvalues_module
00049   * This function computes the eigenvalues with the help of the EigenSolver
00050   * class (for real matrices) or the ComplexEigenSolver class (for complex
00051   * matrices). 
00052   *
00053   * The eigenvalues are repeated according to their algebraic multiplicity,
00054   * so there are as many eigenvalues as rows in the matrix.
00055   *
00056   * The SelfAdjointView class provides a better algorithm for selfadjoint
00057   * matrices.
00058   *
00059   * Example: \include MatrixBase_eigenvalues.cpp
00060   * Output: \verbinclude MatrixBase_eigenvalues.out
00061   *
00062   * \sa EigenSolver::eigenvalues(), ComplexEigenSolver::eigenvalues(),
00063   *     SelfAdjointView::eigenvalues()
00064   */
00065 template<typename Derived>
00066 inline typename MatrixBase<Derived>::EigenvaluesReturnType
00067 MatrixBase<Derived>::eigenvalues() const
00068 {
00069   typedef typename internal::traits<Derived>::Scalar Scalar;
00070   return internal::eigenvalues_selector<Derived, NumTraits<Scalar>::IsComplex>::run(derived());
00071 }
00072 
00073 /** \brief Computes the eigenvalues of a matrix
00074   * \returns Column vector containing the eigenvalues.
00075   *
00076   * \eigenvalues_module
00077   * This function computes the eigenvalues with the help of the
00078   * SelfAdjointEigenSolver class.  The eigenvalues are repeated according to
00079   * their algebraic multiplicity, so there are as many eigenvalues as rows in
00080   * the matrix.
00081   *
00082   * Example: \include SelfAdjointView_eigenvalues.cpp
00083   * Output: \verbinclude SelfAdjointView_eigenvalues.out
00084   *
00085   * \sa SelfAdjointEigenSolver::eigenvalues(), MatrixBase::eigenvalues()
00086   */
00087 template<typename MatrixType, unsigned int UpLo> 
00088 inline typename SelfAdjointView<MatrixType, UpLo>::EigenvaluesReturnType
00089 SelfAdjointView<MatrixType, UpLo>::eigenvalues() const
00090 {
00091   typedef typename SelfAdjointView<MatrixType, UpLo>::PlainObject PlainObject;
00092   PlainObject thisAsMatrix(*this);
00093   return SelfAdjointEigenSolver<PlainObject> (thisAsMatrix, false).eigenvalues();
00094 }
00095 
00096 
00097 
00098 /** \brief Computes the L2 operator norm
00099   * \returns Operator norm of the matrix.
00100   *
00101   * \eigenvalues_module
00102   * This function computes the L2 operator norm of a matrix, which is also
00103   * known as the spectral norm. The norm of a matrix \f$ A \f$ is defined to be
00104   * \f[ \|A\|_2 = \max_x \frac{\|Ax\|_2}{\|x\|_2} \f]
00105   * where the maximum is over all vectors and the norm on the right is the
00106   * Euclidean vector norm. The norm equals the largest singular value, which is
00107   * the square root of the largest eigenvalue of the positive semi-definite
00108   * matrix \f$ A^*A \f$.
00109   *
00110   * The current implementation uses the eigenvalues of \f$ A^*A \f$, as computed
00111   * by SelfAdjointView::eigenvalues(), to compute the operator norm of a
00112   * matrix.  The SelfAdjointView class provides a better algorithm for
00113   * selfadjoint matrices.
00114   *
00115   * Example: \include MatrixBase_operatorNorm.cpp
00116   * Output: \verbinclude MatrixBase_operatorNorm.out
00117   *
00118   * \sa SelfAdjointView::eigenvalues(), SelfAdjointView::operatorNorm()
00119   */
00120 template<typename Derived>
00121 inline typename MatrixBase<Derived>::RealScalar
00122 MatrixBase<Derived>::operatorNorm() const
00123 {
00124   using std::sqrt;
00125   typename Derived::PlainObject m_eval(derived());
00126   // FIXME if it is really guaranteed that the eigenvalues are already sorted,
00127   // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
00128   return sqrt((m_eval*m_eval.adjoint())
00129                  .eval()
00130          .template selfadjointView<Lower>()
00131          .eigenvalues()
00132          .maxCoeff()
00133          );
00134 }
00135 
00136 /** \brief Computes the L2 operator norm
00137   * \returns Operator norm of the matrix.
00138   *
00139   * \eigenvalues_module
00140   * This function computes the L2 operator norm of a self-adjoint matrix. For a
00141   * self-adjoint matrix, the operator norm is the largest eigenvalue.
00142   *
00143   * The current implementation uses the eigenvalues of the matrix, as computed
00144   * by eigenvalues(), to compute the operator norm of the matrix.
00145   *
00146   * Example: \include SelfAdjointView_operatorNorm.cpp
00147   * Output: \verbinclude SelfAdjointView_operatorNorm.out
00148   *
00149   * \sa eigenvalues(), MatrixBase::operatorNorm()
00150   */
00151 template<typename MatrixType, unsigned int UpLo>
00152 inline typename SelfAdjointView<MatrixType, UpLo>::RealScalar
00153 SelfAdjointView<MatrixType, UpLo>::operatorNorm() const
00154 {
00155   return eigenvalues().cwiseAbs().maxCoeff();
00156 }
00157 
00158 } // end namespace Eigen
00159 
00160 #endif