GeneralizedSelfAdjointEigenSolver.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2010 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_GENERALIZEDSELFADJOINTEIGENSOLVER_H
00012 #define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
00013 
00014 #include "./Tridiagonalization.h"
00015 
00016 namespace Eigen { 
00017 
00018 /** \eigenvalues_module \ingroup Eigenvalues_Module
00019   *
00020   *
00021   * \class GeneralizedSelfAdjointEigenSolver
00022   *
00023   * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem
00024   *
00025   * \tparam _MatrixType the type of the matrix of which we are computing the
00026   * eigendecomposition; this is expected to be an instantiation of the Matrix
00027   * class template.
00028   *
00029   * This class solves the generalized eigenvalue problem
00030   * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be
00031   * selfadjoint and the matrix \f$ B \f$ should be positive definite.
00032   *
00033   * Only the \b lower \b triangular \b part of the input matrix is referenced.
00034   *
00035   * Call the function compute() to compute the eigenvalues and eigenvectors of
00036   * a given matrix. Alternatively, you can use the
00037   * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
00038   * constructor which computes the eigenvalues and eigenvectors at construction time.
00039   * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues()
00040   * and eigenvectors() functions.
00041   *
00042   * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
00043   * contains an example of the typical use of this class.
00044   *
00045   * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver
00046   */
00047 template<typename _MatrixType>
00048 class GeneralizedSelfAdjointEigenSolver  : public SelfAdjointEigenSolver <_MatrixType>
00049 {
00050     typedef SelfAdjointEigenSolver<_MatrixType>  Base ;
00051   public:
00052 
00053     typedef typename Base::Index Index;
00054     typedef _MatrixType MatrixType;
00055 
00056     /** \brief Default constructor for fixed-size matrices.
00057       *
00058       * The default constructor is useful in cases in which the user intends to
00059       * perform decompositions via compute(). This constructor
00060       * can only be used if \p _MatrixType is a fixed-size matrix; use
00061       * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices.
00062       */
00063     GeneralizedSelfAdjointEigenSolver() : Base () {}
00064 
00065     /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
00066       *
00067       * \param [in]  size  Positive integer, size of the matrix whose
00068       * eigenvalues and eigenvectors will be computed.
00069       *
00070       * This constructor is useful for dynamic-size matrices, when the user
00071       * intends to perform decompositions via compute(). The \p size
00072       * parameter is only used as a hint. It is not an error to give a wrong
00073       * \p size, but it may impair performance.
00074       *
00075       * \sa compute() for an example
00076       */
00077     GeneralizedSelfAdjointEigenSolver(Index size)
00078         : Base (size)
00079     {}
00080 
00081     /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil.
00082       *
00083       * \param[in]  matA  Selfadjoint matrix in matrix pencil.
00084       *                   Only the lower triangular part of the matrix is referenced.
00085       * \param[in]  matB  Positive-definite matrix in matrix pencil.
00086       *                   Only the lower triangular part of the matrix is referenced.
00087       * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
00088       *                     Default is #ComputeEigenvectors|#Ax_lBx.
00089       *
00090       * This constructor calls compute(const MatrixType&, const MatrixType&, int)
00091       * to compute the eigenvalues and (if requested) the eigenvectors of the
00092       * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the
00093       * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix
00094       * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property
00095       * \f$ x^* B x = 1 \f$. The eigenvectors are computed if
00096       * \a options contains ComputeEigenvectors.
00097       *
00098       * In addition, the two following variants can be solved via \p options:
00099       * - \c ABx_lx: \f$ ABx = \lambda x \f$
00100       * - \c BAx_lx: \f$ BAx = \lambda x \f$
00101       *
00102       * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp
00103       * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out
00104       *
00105       * \sa compute(const MatrixType&, const MatrixType&, int)
00106       */
00107     GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
00108                                       int options = ComputeEigenvectors|Ax_lBx)
00109       : Base (matA.cols())
00110     {
00111       compute(matA, matB, options);
00112     }
00113 
00114     /** \brief Computes generalized eigendecomposition of given matrix pencil.
00115       *
00116       * \param[in]  matA  Selfadjoint matrix in matrix pencil.
00117       *                   Only the lower triangular part of the matrix is referenced.
00118       * \param[in]  matB  Positive-definite matrix in matrix pencil.
00119       *                   Only the lower triangular part of the matrix is referenced.
00120       * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
00121       *                     Default is #ComputeEigenvectors|#Ax_lBx.
00122       *
00123       * \returns    Reference to \c *this
00124       *
00125       * Accoring to \p options, this function computes eigenvalues and (if requested)
00126       * the eigenvectors of one of the following three generalized eigenproblems:
00127       * - \c Ax_lBx: \f$ Ax = \lambda B x \f$
00128       * - \c ABx_lx: \f$ ABx = \lambda x \f$
00129       * - \c BAx_lx: \f$ BAx = \lambda x \f$
00130       * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite
00131       * matrix \f$ B \f$.
00132       * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$.
00133       *
00134       * The eigenvalues() function can be used to retrieve
00135       * the eigenvalues. If \p options contains ComputeEigenvectors, then the
00136       * eigenvectors are also computed and can be retrieved by calling
00137       * eigenvectors().
00138       *
00139       * The implementation uses LLT to compute the Cholesky decomposition
00140       * \f$ B = LL^* \f$ and computes the classical eigendecomposition
00141       * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx
00142       * and of \f$ L^{*} A L \f$ otherwise. This solves the
00143       * generalized eigenproblem, because any solution of the generalized
00144       * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution
00145       * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the
00146       * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements
00147       * can be made for the two other variants.
00148       *
00149       * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp
00150       * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out
00151       *
00152       * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
00153       */
00154     GeneralizedSelfAdjointEigenSolver & compute(const MatrixType& matA, const MatrixType& matB,
00155                                                int options = ComputeEigenvectors|Ax_lBx);
00156 
00157   protected:
00158 
00159 };
00160 
00161 
00162 template<typename MatrixType>
00163 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
00164 compute(const MatrixType& matA, const MatrixType& matB, int options)
00165 {
00166   eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
00167   eigen_assert((options&~(EigVecMask|GenEigMask))==0
00168           && (options&EigVecMask)!=EigVecMask
00169           && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
00170            || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
00171           && "invalid option parameter");
00172 
00173   bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
00174 
00175   // Compute the cholesky decomposition of matB = L L' = U'U
00176   LLT<MatrixType> cholB(matB);
00177 
00178   int type = (options&GenEigMask);
00179   if(type==0)
00180     type = Ax_lBx;
00181 
00182   if(type==Ax_lBx)
00183   {
00184     // compute C = inv(L) A inv(L')
00185     MatrixType matC = matA.template selfadjointView<Lower>();
00186     cholB.matrixL ().template solveInPlace<OnTheLeft>(matC);
00187     cholB.matrixU ().template solveInPlace<OnTheRight>(matC);
00188 
00189     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
00190 
00191     // transform back the eigen vectors: evecs = inv(U) * evecs
00192     if(computeEigVecs)
00193       cholB.matrixU ().solveInPlace(Base::m_eivec);
00194   }
00195   else if(type==ABx_lx)
00196   {
00197     // compute C = L' A L
00198     MatrixType matC = matA.template selfadjointView<Lower>();
00199     matC = matC * cholB.matrixL ();
00200     matC = cholB.matrixU () * matC;
00201 
00202     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
00203 
00204     // transform back the eigen vectors: evecs = inv(U) * evecs
00205     if(computeEigVecs)
00206       cholB.matrixU ().solveInPlace(Base::m_eivec);
00207   }
00208   else if(type==BAx_lx)
00209   {
00210     // compute C = L' A L
00211     MatrixType matC = matA.template selfadjointView<Lower>();
00212     matC = matC * cholB.matrixL ();
00213     matC = cholB.matrixU () * matC;
00214 
00215     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
00216 
00217     // transform back the eigen vectors: evecs = L * evecs
00218     if(computeEigVecs)
00219       Base::m_eivec = cholB.matrixL () * Base::m_eivec;
00220   }
00221 
00222   return *this;
00223 }
00224 
00225 } // end namespace Eigen
00226 
00227 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H