| 1 | // This file is part of Eigen, a lightweight C++ template library
|
|---|
| 2 | // for linear algebra.
|
|---|
| 3 | //
|
|---|
| 4 | // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
|
|---|
| 5 | // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
|
|---|
| 6 | //
|
|---|
| 7 | // This Source Code Form is subject to the terms of the Mozilla
|
|---|
| 8 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
|---|
| 9 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
|---|
| 10 |
|
|---|
| 11 | #ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
|
|---|
| 12 | #define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
|
|---|
| 13 |
|
|---|
| 14 | #include "./Tridiagonalization.h"
|
|---|
| 15 |
|
|---|
| 16 | namespace Eigen {
|
|---|
| 17 |
|
|---|
| 18 | /** \eigenvalues_module \ingroup Eigenvalues_Module
|
|---|
| 19 | *
|
|---|
| 20 | *
|
|---|
| 21 | * \class GeneralizedSelfAdjointEigenSolver
|
|---|
| 22 | *
|
|---|
| 23 | * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem
|
|---|
| 24 | *
|
|---|
| 25 | * \tparam _MatrixType the type of the matrix of which we are computing the
|
|---|
| 26 | * eigendecomposition; this is expected to be an instantiation of the Matrix
|
|---|
| 27 | * class template.
|
|---|
| 28 | *
|
|---|
| 29 | * This class solves the generalized eigenvalue problem
|
|---|
| 30 | * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be
|
|---|
| 31 | * selfadjoint and the matrix \f$ B \f$ should be positive definite.
|
|---|
| 32 | *
|
|---|
| 33 | * Only the \b lower \b triangular \b part of the input matrix is referenced.
|
|---|
| 34 | *
|
|---|
| 35 | * Call the function compute() to compute the eigenvalues and eigenvectors of
|
|---|
| 36 | * a given matrix. Alternatively, you can use the
|
|---|
| 37 | * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
|
|---|
| 38 | * constructor which computes the eigenvalues and eigenvectors at construction time.
|
|---|
| 39 | * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues()
|
|---|
| 40 | * and eigenvectors() functions.
|
|---|
| 41 | *
|
|---|
| 42 | * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
|
|---|
| 43 | * contains an example of the typical use of this class.
|
|---|
| 44 | *
|
|---|
| 45 | * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver
|
|---|
| 46 | */
|
|---|
| 47 | template<typename _MatrixType>
|
|---|
| 48 | class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
|
|---|
| 49 | {
|
|---|
| 50 | typedef SelfAdjointEigenSolver<_MatrixType> Base;
|
|---|
| 51 | public:
|
|---|
| 52 |
|
|---|
| 53 | typedef typename Base::Index Index;
|
|---|
| 54 | typedef _MatrixType MatrixType;
|
|---|
| 55 |
|
|---|
| 56 | /** \brief Default constructor for fixed-size matrices.
|
|---|
| 57 | *
|
|---|
| 58 | * The default constructor is useful in cases in which the user intends to
|
|---|
| 59 | * perform decompositions via compute(). This constructor
|
|---|
| 60 | * can only be used if \p _MatrixType is a fixed-size matrix; use
|
|---|
| 61 | * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices.
|
|---|
| 62 | */
|
|---|
| 63 | GeneralizedSelfAdjointEigenSolver() : Base() {}
|
|---|
| 64 |
|
|---|
| 65 | /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
|
|---|
| 66 | *
|
|---|
| 67 | * \param [in] size Positive integer, size of the matrix whose
|
|---|
| 68 | * eigenvalues and eigenvectors will be computed.
|
|---|
| 69 | *
|
|---|
| 70 | * This constructor is useful for dynamic-size matrices, when the user
|
|---|
| 71 | * intends to perform decompositions via compute(). The \p size
|
|---|
| 72 | * parameter is only used as a hint. It is not an error to give a wrong
|
|---|
| 73 | * \p size, but it may impair performance.
|
|---|
| 74 | *
|
|---|
| 75 | * \sa compute() for an example
|
|---|
| 76 | */
|
|---|
| 77 | GeneralizedSelfAdjointEigenSolver(Index size)
|
|---|
| 78 | : Base(size)
|
|---|
| 79 | {}
|
|---|
| 80 |
|
|---|
| 81 | /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil.
|
|---|
| 82 | *
|
|---|
| 83 | * \param[in] matA Selfadjoint matrix in matrix pencil.
|
|---|
| 84 | * Only the lower triangular part of the matrix is referenced.
|
|---|
| 85 | * \param[in] matB Positive-definite matrix in matrix pencil.
|
|---|
| 86 | * Only the lower triangular part of the matrix is referenced.
|
|---|
| 87 | * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
|
|---|
| 88 | * Default is #ComputeEigenvectors|#Ax_lBx.
|
|---|
| 89 | *
|
|---|
| 90 | * This constructor calls compute(const MatrixType&, const MatrixType&, int)
|
|---|
| 91 | * to compute the eigenvalues and (if requested) the eigenvectors of the
|
|---|
| 92 | * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the
|
|---|
| 93 | * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix
|
|---|
| 94 | * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property
|
|---|
| 95 | * \f$ x^* B x = 1 \f$. The eigenvectors are computed if
|
|---|
| 96 | * \a options contains ComputeEigenvectors.
|
|---|
| 97 | *
|
|---|
| 98 | * In addition, the two following variants can be solved via \p options:
|
|---|
| 99 | * - \c ABx_lx: \f$ ABx = \lambda x \f$
|
|---|
| 100 | * - \c BAx_lx: \f$ BAx = \lambda x \f$
|
|---|
| 101 | *
|
|---|
| 102 | * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp
|
|---|
| 103 | * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out
|
|---|
| 104 | *
|
|---|
| 105 | * \sa compute(const MatrixType&, const MatrixType&, int)
|
|---|
| 106 | */
|
|---|
| 107 | GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
|
|---|
| 108 | int options = ComputeEigenvectors|Ax_lBx)
|
|---|
| 109 | : Base(matA.cols())
|
|---|
| 110 | {
|
|---|
| 111 | compute(matA, matB, options);
|
|---|
| 112 | }
|
|---|
| 113 |
|
|---|
| 114 | /** \brief Computes generalized eigendecomposition of given matrix pencil.
|
|---|
| 115 | *
|
|---|
| 116 | * \param[in] matA Selfadjoint matrix in matrix pencil.
|
|---|
| 117 | * Only the lower triangular part of the matrix is referenced.
|
|---|
| 118 | * \param[in] matB Positive-definite matrix in matrix pencil.
|
|---|
| 119 | * Only the lower triangular part of the matrix is referenced.
|
|---|
| 120 | * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
|
|---|
| 121 | * Default is #ComputeEigenvectors|#Ax_lBx.
|
|---|
| 122 | *
|
|---|
| 123 | * \returns Reference to \c *this
|
|---|
| 124 | *
|
|---|
| 125 | * Accoring to \p options, this function computes eigenvalues and (if requested)
|
|---|
| 126 | * the eigenvectors of one of the following three generalized eigenproblems:
|
|---|
| 127 | * - \c Ax_lBx: \f$ Ax = \lambda B x \f$
|
|---|
| 128 | * - \c ABx_lx: \f$ ABx = \lambda x \f$
|
|---|
| 129 | * - \c BAx_lx: \f$ BAx = \lambda x \f$
|
|---|
| 130 | * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite
|
|---|
| 131 | * matrix \f$ B \f$.
|
|---|
| 132 | * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$.
|
|---|
| 133 | *
|
|---|
| 134 | * The eigenvalues() function can be used to retrieve
|
|---|
| 135 | * the eigenvalues. If \p options contains ComputeEigenvectors, then the
|
|---|
| 136 | * eigenvectors are also computed and can be retrieved by calling
|
|---|
| 137 | * eigenvectors().
|
|---|
| 138 | *
|
|---|
| 139 | * The implementation uses LLT to compute the Cholesky decomposition
|
|---|
| 140 | * \f$ B = LL^* \f$ and computes the classical eigendecomposition
|
|---|
| 141 | * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx
|
|---|
| 142 | * and of \f$ L^{*} A L \f$ otherwise. This solves the
|
|---|
| 143 | * generalized eigenproblem, because any solution of the generalized
|
|---|
| 144 | * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution
|
|---|
| 145 | * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the
|
|---|
| 146 | * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements
|
|---|
| 147 | * can be made for the two other variants.
|
|---|
| 148 | *
|
|---|
| 149 | * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp
|
|---|
| 150 | * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out
|
|---|
| 151 | *
|
|---|
| 152 | * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
|
|---|
| 153 | */
|
|---|
| 154 | GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
|
|---|
| 155 | int options = ComputeEigenvectors|Ax_lBx);
|
|---|
| 156 |
|
|---|
| 157 | protected:
|
|---|
| 158 |
|
|---|
| 159 | };
|
|---|
| 160 |
|
|---|
| 161 |
|
|---|
| 162 | template<typename MatrixType>
|
|---|
| 163 | GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
|
|---|
| 164 | compute(const MatrixType& matA, const MatrixType& matB, int options)
|
|---|
| 165 | {
|
|---|
| 166 | eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
|
|---|
| 167 | eigen_assert((options&~(EigVecMask|GenEigMask))==0
|
|---|
| 168 | && (options&EigVecMask)!=EigVecMask
|
|---|
| 169 | && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
|
|---|
| 170 | || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
|
|---|
| 171 | && "invalid option parameter");
|
|---|
| 172 |
|
|---|
| 173 | bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
|
|---|
| 174 |
|
|---|
| 175 | // Compute the cholesky decomposition of matB = L L' = U'U
|
|---|
| 176 | LLT<MatrixType> cholB(matB);
|
|---|
| 177 |
|
|---|
| 178 | int type = (options&GenEigMask);
|
|---|
| 179 | if(type==0)
|
|---|
| 180 | type = Ax_lBx;
|
|---|
| 181 |
|
|---|
| 182 | if(type==Ax_lBx)
|
|---|
| 183 | {
|
|---|
| 184 | // compute C = inv(L) A inv(L')
|
|---|
| 185 | MatrixType matC = matA.template selfadjointView<Lower>();
|
|---|
| 186 | cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
|
|---|
| 187 | cholB.matrixU().template solveInPlace<OnTheRight>(matC);
|
|---|
| 188 |
|
|---|
| 189 | Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
|
|---|
| 190 |
|
|---|
| 191 | // transform back the eigen vectors: evecs = inv(U) * evecs
|
|---|
| 192 | if(computeEigVecs)
|
|---|
| 193 | cholB.matrixU().solveInPlace(Base::m_eivec);
|
|---|
| 194 | }
|
|---|
| 195 | else if(type==ABx_lx)
|
|---|
| 196 | {
|
|---|
| 197 | // compute C = L' A L
|
|---|
| 198 | MatrixType matC = matA.template selfadjointView<Lower>();
|
|---|
| 199 | matC = matC * cholB.matrixL();
|
|---|
| 200 | matC = cholB.matrixU() * matC;
|
|---|
| 201 |
|
|---|
| 202 | Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
|
|---|
| 203 |
|
|---|
| 204 | // transform back the eigen vectors: evecs = inv(U) * evecs
|
|---|
| 205 | if(computeEigVecs)
|
|---|
| 206 | cholB.matrixU().solveInPlace(Base::m_eivec);
|
|---|
| 207 | }
|
|---|
| 208 | else if(type==BAx_lx)
|
|---|
| 209 | {
|
|---|
| 210 | // compute C = L' A L
|
|---|
| 211 | MatrixType matC = matA.template selfadjointView<Lower>();
|
|---|
| 212 | matC = matC * cholB.matrixL();
|
|---|
| 213 | matC = cholB.matrixU() * matC;
|
|---|
| 214 |
|
|---|
| 215 | Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
|
|---|
| 216 |
|
|---|
| 217 | // transform back the eigen vectors: evecs = L * evecs
|
|---|
| 218 | if(computeEigVecs)
|
|---|
| 219 | Base::m_eivec = cholB.matrixL() * Base::m_eivec;
|
|---|
| 220 | }
|
|---|
| 221 |
|
|---|
| 222 | return *this;
|
|---|
| 223 | }
|
|---|
| 224 |
|
|---|
| 225 | } // end namespace Eigen
|
|---|
| 226 |
|
|---|
| 227 | #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
|
|---|
