Main MRPT website > C++ reference
MRPT logo

GeneralizedSelfAdjointEigenSolver.h

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



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:16:28 UTC 2011