diff options
Diffstat (limited to 'Eigen/src/EigenSolver/SelfAdjointEigenSolver.h')
| -rw-r--r-- | Eigen/src/EigenSolver/SelfAdjointEigenSolver.h | 366 |
1 files changed, 366 insertions, 0 deletions
diff --git a/Eigen/src/EigenSolver/SelfAdjointEigenSolver.h b/Eigen/src/EigenSolver/SelfAdjointEigenSolver.h new file mode 100644 index 000000000..40b06df2c --- /dev/null +++ b/Eigen/src/EigenSolver/SelfAdjointEigenSolver.h @@ -0,0 +1,366 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// 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_SELFADJOINTEIGENSOLVER_H +#define EIGEN_SELFADJOINTEIGENSOLVER_H + +/** \eigensolver_module \ingroup EigenSolver_Module + * \nonstableyet + * + * \class SelfAdjointEigenSolver + * + * \brief Eigen values/vectors solver for selfadjoint matrix + * + * \param MatrixType the type of the matrix of which we are computing the eigen decomposition + * + * \note MatrixType must be an actual Matrix type, it can't be an expression type. + * + * \sa MatrixBase::eigenvalues(), class EigenSolver + */ +template<typename _MatrixType> class SelfAdjointEigenSolver +{ + public: + + enum {Size = _MatrixType::RowsAtCompileTime }; + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef std::complex<RealScalar> Complex; + typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType; + typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX; + typedef Tridiagonalization<MatrixType> TridiagonalizationType; +// typedef typename TridiagonalizationType::TridiagonalMatrixType TridiagonalMatrixType; + + SelfAdjointEigenSolver() + : m_eivec(int(Size), int(Size)), + m_eivalues(int(Size)) + { + ei_assert(Size!=Dynamic); + } + + SelfAdjointEigenSolver(int size) + : m_eivec(size, size), + m_eivalues(size) + {} + + /** Constructors computing the eigenvalues of the selfadjoint matrix \a matrix, + * as well as the eigenvectors if \a computeEigenvectors is true. + * + * \sa compute(MatrixType,bool), SelfAdjointEigenSolver(MatrixType,MatrixType,bool) + */ + SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) + : m_eivec(matrix.rows(), matrix.cols()), + m_eivalues(matrix.cols()) + { + compute(matrix, computeEigenvectors); + } + + /** Constructors computing the eigenvalues of the generalized eigen problem + * \f$ Ax = lambda B x \f$ with \a matA the selfadjoint matrix \f$ A \f$ + * and \a matB the positive definite matrix \f$ B \f$ . The eigenvectors + * are computed if \a computeEigenvectors is true. + * + * \sa compute(MatrixType,MatrixType,bool), SelfAdjointEigenSolver(MatrixType,bool) + */ + SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) + : m_eivec(matA.rows(), matA.cols()), + m_eivalues(matA.cols()) + { + compute(matA, matB, computeEigenvectors); + } + + SelfAdjointEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); + + SelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true); + + /** \returns the computed eigen vectors as a matrix of column vectors */ + MatrixType eigenvectors(void) const + { + #ifndef NDEBUG + ei_assert(m_eigenvectorsOk); + #endif + return m_eivec; + } + + /** \returns the computed eigen values */ + RealVectorType eigenvalues(void) const { return m_eivalues; } + + /** \returns the positive square root of the matrix + * + * \note the matrix itself must be positive in order for this to make sense. + */ + MatrixType operatorSqrt() const + { + return m_eivec * m_eivalues.cwise().sqrt().asDiagonal() * m_eivec.adjoint(); + } + + /** \returns the positive inverse square root of the matrix + * + * \note the matrix itself must be positive definite in order for this to make sense. + */ + MatrixType operatorInverseSqrt() const + { + return m_eivec * m_eivalues.cwise().inverse().cwise().sqrt().asDiagonal() * m_eivec.adjoint(); + } + + + protected: + MatrixType m_eivec; + RealVectorType m_eivalues; + #ifndef NDEBUG + bool m_eigenvectorsOk; + #endif +}; + +#ifndef EIGEN_HIDE_HEAVY_CODE + +/** \internal + * + * \eigensolver_module + * + * Performs a QR step on a tridiagonal symmetric matrix represented as a + * pair of two vectors \a diag and \a subdiag. + * + * \param matA the input selfadjoint matrix + * \param hCoeffs returned Householder coefficients + * + * For compilation efficiency reasons, this procedure does not use eigen expression + * for its arguments. + * + * Implemented from Golub's "Matrix Computations", algorithm 8.3.2: + * "implicit symmetric QR step with Wilkinson shift" + */ +template<typename RealScalar, typename Scalar> +static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n); + +/** Computes the eigenvalues of the selfadjoint matrix \a matrix, + * as well as the eigenvectors if \a computeEigenvectors is true. + * + * \sa SelfAdjointEigenSolver(MatrixType,bool), compute(MatrixType,MatrixType,bool) + */ +template<typename MatrixType> +SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) +{ + #ifndef NDEBUG + m_eigenvectorsOk = computeEigenvectors; + #endif + assert(matrix.cols() == matrix.rows()); + int n = matrix.cols(); + m_eivalues.resize(n,1); + + if(n==1) + { + m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0)); + m_eivec.setOnes(); + return *this; + } + + m_eivec = matrix; + + // FIXME, should tridiag be a local variable of this function or an attribute of SelfAdjointEigenSolver ? + // the latter avoids multiple memory allocation when the same SelfAdjointEigenSolver is used multiple times... + // (same for diag and subdiag) + RealVectorType& diag = m_eivalues; + typename TridiagonalizationType::SubDiagonalType subdiag(n-1); + TridiagonalizationType::decomposeInPlace(m_eivec, diag, subdiag, computeEigenvectors); + + int end = n-1; + int start = 0; + while (end>0) + { + for (int i = start; i<end; ++i) + if (ei_isMuchSmallerThan(ei_abs(subdiag[i]),(ei_abs(diag[i])+ei_abs(diag[i+1])))) + subdiag[i] = 0; + + // find the largest unreduced block + while (end>0 && subdiag[end-1]==0) + end--; + if (end<=0) + break; + start = end - 1; + while (start>0 && subdiag[start-1]!=0) + start--; + + ei_tridiagonal_qr_step(diag.data(), subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n); + } + + // Sort eigenvalues and corresponding vectors. + // TODO make the sort optional ? + // TODO use a better sort algorithm !! + for (int i = 0; i < n-1; ++i) + { + int k; + m_eivalues.segment(i,n-i).minCoeff(&k); + if (k > 0) + { + std::swap(m_eivalues[i], m_eivalues[k+i]); + m_eivec.col(i).swap(m_eivec.col(k+i)); + } + } + return *this; +} + +/** Computes the eigenvalues of the generalized eigen problem + * \f$ Ax = lambda B x \f$ with \a matA the selfadjoint matrix \f$ A \f$ + * and \a matB the positive definite matrix \f$ B \f$ . The eigenvectors + * are computed if \a computeEigenvectors is true. + * + * \sa SelfAdjointEigenSolver(MatrixType,MatrixType,bool), compute(MatrixType,bool) + */ +template<typename MatrixType> +SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>:: +compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors) +{ + ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); + + // Compute the cholesky decomposition of matB = L L' + LLT<MatrixType> cholB(matB); + + // compute C = inv(L) A inv(L') + MatrixType matC = matA; + cholB.matrixL().solveInPlace(matC); + // FIXME since we currently do not support A * inv(L'), let's do (inv(L) A')' : + matC.adjointInPlace(); + cholB.matrixL().solveInPlace(matC); + matC.adjointInPlace(); + // this version works too: +// matC = matC.transpose(); +// cholB.matrixL().conjugate().template marked<LowerTriangular>().solveTriangularInPlace(matC); +// matC = matC.transpose(); + // FIXME: this should work: (currently it only does for small matrices) +// Transpose<MatrixType> trMatC(matC); +// cholB.matrixL().conjugate().eval().template marked<LowerTriangular>().solveTriangularInPlace(trMatC); + + compute(matC, computeEigenvectors); + + if (computeEigenvectors) + { + // transform back the eigen vectors: evecs = inv(U) * evecs + cholB.matrixU().solveInPlace(m_eivec); + for (int i=0; i<m_eivec.cols(); ++i) + m_eivec.col(i) = m_eivec.col(i).normalized(); + } + return *this; +} + +#endif // EIGEN_HIDE_HEAVY_CODE + +/** \eigensolver_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&SelfAdjointBit); + return SelfAdjointEigenSolver<typename Derived::PlainMatrixType>(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().cwise().abs().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::PlainMatrixType 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() + ); + } +}; + +/** \eigensolver_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&SelfAdjointBit> + ::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) +{ + RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5); + RealScalar e2 = ei_abs2(subdiag[end-1]); + RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * ei_sqrt(td*td + e2)); + RealScalar x = diag[start] - mu; + RealScalar z = subdiag[start]; + + for (int k = start; k < end; ++k) + { + PlanarRotation<RealScalar> rot; + rot.makeGivens(x, z); + + // do T = G' T G + RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k]; + RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1]; + + diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]); + diag[k+1] = rot.s() * sdk + rot.c() * dkp1; + subdiag[k] = rot.c() * sdk - rot.s() * dkp1; + + if (k > start) + subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z; + + x = subdiag[k]; + + if (k < end - 1) + { + z = -rot.s() * subdiag[k+1]; + subdiag[k + 1] = rot.c() * subdiag[k+1]; + } + + // apply the givens rotation to the unit matrix Q = Q * G + if (matrixQ) + { + Map<Matrix<Scalar,Dynamic,Dynamic> > q(matrixQ,n,n); + q.applyOnTheRight(k,k+1,rot); + } + } +} +#endif + +#endif // EIGEN_SELFADJOINTEIGENSOLVER_H |