diff options
| -rw-r--r-- | Eigen/QR | 1 | ||||
| -rw-r--r-- | Eigen/src/QR/EigenSolver.h | 293 | ||||
| -rw-r--r-- | Eigen/src/QR/SelfAdjointEigenSolver.h | 172 | ||||
| -rwxr-xr-x | Eigen/src/QR/Tridiagonalization.h | 4 |
4 files changed, 194 insertions, 276 deletions
@@ -8,6 +8,7 @@ namespace Eigen { #include "src/QR/QR.h" #include "src/QR/Tridiagonalization.h" #include "src/QR/EigenSolver.h" +#include "src/QR/SelfAdjointEigenSolver.h" } // namespace Eigen diff --git a/Eigen/src/QR/EigenSolver.h b/Eigen/src/QR/EigenSolver.h index 4e875757e..2433781cd 100644 --- a/Eigen/src/QR/EigenSolver.h +++ b/Eigen/src/QR/EigenSolver.h @@ -27,19 +27,17 @@ /** \class EigenSolver * - * \brief Eigen values/vectors solver + * \brief Eigen values/vectors solver for non selfadjoint matrices * * \param MatrixType the type of the matrix of which we are computing the eigen decomposition - * \param IsSelfadjoint tells the input matrix is guaranteed to be selfadjoint (hermitian). In that case the - * return type of eigenvalues() is a real vector. * * Currently it only support real matrices. * * \note this code was adapted from JAMA (public domain) * - * \sa MatrixBase::eigenvalues() + * \sa MatrixBase::eigenvalues(), SelfAdjointEigenSolver */ -template<typename _MatrixType, bool IsSelfadjoint=false> class EigenSolver +template<typename _MatrixType> class EigenSolver { public: @@ -47,7 +45,7 @@ template<typename _MatrixType, bool IsSelfadjoint=false> class EigenSolver typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; typedef std::complex<RealScalar> Complex; - typedef Matrix<typename ei_meta_if<IsSelfadjoint, Scalar, Complex>::ret, MatrixType::ColsAtCompileTime, 1> EigenvalueType; + typedef Matrix<Complex, MatrixType::ColsAtCompileTime, 1> EigenvalueType; typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType; typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX; @@ -55,7 +53,7 @@ template<typename _MatrixType, bool IsSelfadjoint=false> class EigenSolver : m_eivec(matrix.rows(), matrix.cols()), m_eivalues(matrix.cols()) { - _compute(matrix); + compute(matrix); } MatrixType eigenvectors(void) const { return m_eivec; } @@ -64,15 +62,7 @@ template<typename _MatrixType, bool IsSelfadjoint=false> class EigenSolver private: - void _compute(const MatrixType& matrix) - { - computeImpl(matrix, typename ei_meta_if<IsSelfadjoint, ei_meta_true, ei_meta_false>::ret()); - } - void computeImpl(const MatrixType& matrix, ei_meta_true isSelfadjoint); - void computeImpl(const MatrixType& matrix, ei_meta_false isNotSelfadjoint); - - void tridiagonalization(RealVectorType& eivalr, RealVectorType& eivali); - void tql2(RealVectorType& eivalr, RealVectorType& eivali); + void compute(const MatrixType& matrix); void orthes(MatrixType& matH, RealVectorType& ort); void hqr2(MatrixType& matH); @@ -82,274 +72,27 @@ template<typename _MatrixType, bool IsSelfadjoint=false> class EigenSolver EigenvalueType m_eivalues; }; -template<typename MatrixType, bool IsSelfadjoint> -void EigenSolver<MatrixType,IsSelfadjoint>::computeImpl(const MatrixType& matrix, ei_meta_true) -{ - assert(matrix.cols() == matrix.rows()); - int n = matrix.cols(); - m_eivalues.resize(n,1); - - RealVectorType eivali(n); - m_eivec = matrix; - - // Tridiagonalize. - tridiagonalization(m_eivalues, eivali); - - // Diagonalize. - tql2(m_eivalues, eivali); -} -template<typename MatrixType, bool IsSelfadjoint> -void EigenSolver<MatrixType,IsSelfadjoint>::computeImpl(const MatrixType& matrix, ei_meta_false) +template<typename MatrixType> +void EigenSolver<MatrixType>::compute(const MatrixType& matrix) { assert(matrix.cols() == matrix.rows()); int n = matrix.cols(); m_eivalues.resize(n,1); - bool isSelfadjoint = true; - for (int j = 0; (j < n) && isSelfadjoint; j++) - for (int i = 0; (i < j) && isSelfadjoint; i++) - isSelfadjoint = (matrix(i,j) == matrix(j,i)); - - if (isSelfadjoint) - { - RealVectorType eivalr(n); - RealVectorType eivali(n); - m_eivec = matrix; - - // Tridiagonalize. - tridiagonalization(eivalr, eivali); - - // Diagonalize. - tql2(eivalr, eivali); - - m_eivalues = eivalr.template cast<Complex>(); - } - else - { - MatrixType matH = matrix; - RealVectorType ort(n); - - // Reduce to Hessenberg form. - orthes(matH, ort); - - // Reduce Hessenberg to real Schur form. - hqr2(matH); - } -} - - -// Symmetric Householder reduction to tridiagonal form. -template<typename MatrixType, bool IsSelfadjoint> -void EigenSolver<MatrixType,IsSelfadjoint>::tridiagonalization(RealVectorType& eivalr, RealVectorType& eivali) -{ - -// This is derived from the Algol procedures tred2 by -// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for -// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding -// Fortran subroutine in EISPACK. - - int n = m_eivec.cols(); - eivalr = m_eivec.row(eivalr.size()-1); - - // Householder reduction to tridiagonal form. - for (int i = n-1; i > 0; i--) - { - // Scale to avoid under/overflow. - Scalar scale = 0.0; - Scalar h = 0.0; - scale = eivalr.start(i).cwiseAbs().sum(); - - if (scale == 0.0) - { - eivali.coeffRef(i) = eivalr.coeff(i-1); - eivalr.start(i) = m_eivec.row(i-1).start(i); - m_eivec.corner(TopLeft, i, i) = m_eivec.corner(TopLeft, i, i).diagonal().asDiagonal(); - } - else - { - // Generate Householder vector. - eivalr.start(i) /= scale; - h = eivalr.start(i).cwiseAbs2().sum(); - - Scalar f = eivalr.coeff(i-1); - Scalar g = ei_sqrt(h); - if (f > 0) - g = -g; - eivali.coeffRef(i) = scale * g; - h = h - f * g; - eivalr.coeffRef(i-1) = f - g; - eivali.start(i).setZero(); - - // Apply similarity transformation to remaining columns. - for (int j = 0; j < i; j++) - { - f = eivalr.coeff(j); - m_eivec.coeffRef(j,i) = f; - g = eivali.coeff(j) + m_eivec.coeff(j,j) * f; - int bSize = i-j-1; - if (bSize>0) - { - g += (m_eivec.col(j).block(j+1, bSize).transpose() * eivalr.block(j+1, bSize))(0,0); - eivali.block(j+1, bSize) += m_eivec.col(j).block(j+1, bSize) * f; - } - eivali.coeffRef(j) = g; - } - - f = (eivali.start(i).transpose() * eivalr.start(i))(0,0); - eivali.start(i) = (eivali.start(i) - (f / (h + h)) * eivalr.start(i))/h; - - m_eivec.corner(TopLeft, i, i).template part<Lower>() -= - ( (eivali.start(i) * eivalr.start(i).transpose()).lazy() - + (eivalr.start(i) * eivali.start(i).transpose()).lazy()); - - eivalr.start(i) = m_eivec.row(i-1).start(i); - m_eivec.row(i).start(i).setZero(); - } - eivalr.coeffRef(i) = h; - } - - // Accumulate transformations. - for (int i = 0; i < n-1; i++) - { - m_eivec.coeffRef(n-1,i) = m_eivec.coeff(i,i); - m_eivec.coeffRef(i,i) = 1.0; - Scalar h = eivalr.coeff(i+1); - // FIXME this does not looks very stable ;) - if (h != 0.0) - { - eivalr.start(i+1) = m_eivec.col(i+1).start(i+1) / h; - m_eivec.corner(TopLeft, i+1, i+1) -= eivalr.start(i+1) - * ( m_eivec.col(i+1).start(i+1).transpose() * m_eivec.corner(TopLeft, i+1, i+1) ); - } - m_eivec.col(i+1).start(i+1).setZero(); - } - eivalr = m_eivec.row(eivalr.size()-1); - m_eivec.row(eivalr.size()-1).setZero(); - m_eivec.coeffRef(n-1,n-1) = 1.0; - eivali.coeffRef(0) = 0.0; -} - - -// Symmetric tridiagonal QL algorithm. -template<typename MatrixType, bool IsSelfadjoint> -void EigenSolver<MatrixType,IsSelfadjoint>::tql2(RealVectorType& eivalr, RealVectorType& eivali) -{ - // This is derived from the Algol procedures tql2, by - // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for - // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutine in EISPACK. - - int n = eivalr.size(); - - for (int i = 1; i < n; i++) { - eivali.coeffRef(i-1) = eivali.coeff(i); - } - eivali.coeffRef(n-1) = 0.0; - - Scalar f = 0.0; - Scalar tst1 = 0.0; - Scalar eps = std::pow(2.0,-52.0); - for (int l = 0; l < n; l++) - { - // Find small subdiagonal element - tst1 = std::max(tst1,ei_abs(eivalr.coeff(l)) + ei_abs(eivali.coeff(l))); - int m = l; - - while ( (m < n) && (ei_abs(eivali.coeff(m)) > eps*tst1) ) - m++; + MatrixType matH = matrix; + RealVectorType ort(n); - // If m == l, eivalr.coeff(l) is an eigenvalue, - // otherwise, iterate. - if (m > l) - { - int iter = 0; - do - { - iter = iter + 1; + // Reduce to Hessenberg form. + orthes(matH, ort); - // Compute implicit shift - Scalar g = eivalr.coeff(l); - Scalar p = (eivalr.coeff(l+1) - g) / (2.0 * eivali.coeff(l)); - Scalar r = hypot(p,1.0); - if (p < 0) - r = -r; - - eivalr.coeffRef(l) = eivali.coeff(l) / (p + r); - eivalr.coeffRef(l+1) = eivali.coeff(l) * (p + r); - Scalar dl1 = eivalr.coeff(l+1); - Scalar h = g - eivalr.coeff(l); - if (l+2<n) - eivalr.end(n-l-2) -= RealVectorTypeX::constant(n-l-2, h); - f = f + h; - - // Implicit QL transformation. - p = eivalr.coeff(m); - Scalar c = 1.0; - Scalar c2 = c; - Scalar c3 = c; - Scalar el1 = eivali.coeff(l+1); - Scalar s = 0.0; - Scalar s2 = 0.0; - for (int i = m-1; i >= l; i--) - { - c3 = c2; - c2 = c; - s2 = s; - g = c * eivali.coeff(i); - h = c * p; - r = hypot(p,eivali.coeff(i)); - eivali.coeffRef(i+1) = s * r; - s = eivali.coeff(i) / r; - c = p / r; - p = c * eivalr.coeff(i) - s * g; - eivalr.coeffRef(i+1) = h + s * (c * g + s * eivalr.coeff(i)); - - // Accumulate transformation. - for (int k = 0; k < n; k++) - { - h = m_eivec.coeff(k,i+1); - m_eivec.coeffRef(k,i+1) = s * m_eivec.coeff(k,i) + c * h; - m_eivec.coeffRef(k,i) = c * m_eivec.coeff(k,i) - s * h; - } - } - p = -s * s2 * c3 * el1 * eivali.coeff(l) / dl1; - eivali.coeffRef(l) = s * p; - eivalr.coeffRef(l) = c * p; - - // Check for convergence. - } while (ei_abs(eivali.coeff(l)) > eps*tst1); - } - eivalr.coeffRef(l) = eivalr.coeff(l) + f; - eivali.coeffRef(l) = 0.0; - } - - // Sort eigenvalues and corresponding vectors. - // TODO use a better sort algorithm !! - for (int i = 0; i < n-1; i++) - { - int k = i; - Scalar minValue = eivalr.coeff(i); - for (int j = i+1; j < n; j++) - { - if (eivalr.coeff(j) < minValue) - { - k = j; - minValue = eivalr.coeff(j); - } - } - if (k != i) - { - std::swap(eivalr.coeffRef(i), eivalr.coeffRef(k)); - m_eivec.col(i).swap(m_eivec.col(k)); - } - } + // Reduce Hessenberg to real Schur form. + hqr2(matH); } - // Nonsymmetric reduction to Hessenberg form. -template<typename MatrixType, bool IsSelfadjoint> -void EigenSolver<MatrixType,IsSelfadjoint>::orthes(MatrixType& matH, RealVectorType& ort) +template<typename MatrixType> +void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort) { // This is derived from the Algol procedures orthes and ortran, // by Martin and Wilkinson, Handbook for Auto. Comp., @@ -432,8 +175,8 @@ std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi) // Nonsymmetric reduction from Hessenberg to real Schur form. -template<typename MatrixType, bool IsSelfadjoint> -void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH) +template<typename MatrixType> +void EigenSolver<MatrixType>::hqr2(MatrixType& matH) { // This is derived from the Algol procedure hqr2, // by Martin and Wilkinson, Handbook for Auto. Comp., diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h new file mode 100644 index 000000000..2ac12fc6e --- /dev/null +++ b/Eigen/src/QR/SelfAdjointEigenSolver.h @@ -0,0 +1,172 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// 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 + +/** \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 + * + * \sa MatrixBase::eigenvalues(), class EigenSolver + */ +template<typename _MatrixType> class SelfAdjointEigenSolver +{ + public: + + 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> EigenvalueType; + typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType; + typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX; + + SelfAdjointEigenSolver(const MatrixType& matrix) + : m_eivec(matrix.rows(), matrix.cols()), + m_eivalues(matrix.cols()) + { + compute(matrix); + } + + void compute(const MatrixType& matrix); + + MatrixType eigenvectors(void) const { return m_eivec; } + + RealVectorType eigenvalues(void) const { return m_eivalues; } + + + protected: + MatrixType m_eivec; + RealVectorType m_eivalues; +}; + +// from Golub's "Matrix Computations", algorithm 5.1.3 +template<typename Scalar> +static void ei_givens_rotation(Scalar a, Scalar b, Scalar& c, Scalar& s) +{ + if (b==0) + { + c = 1; s = 0; + } + else if (ei_abs(b)>ei_abs(a)) + { + Scalar t = -a/b; + s = Scalar(1)/ei_sqrt(1+t*t); + c = s * t; + } + else + { + Scalar t = -b/a; + c = Scalar(1)/ei_sqrt(1+t*t); + s = c * t; + } +} + +/** \internal + * Performs a QR step on a tridiagonal symmetric matrix represented as a + * pair of two vectors \a diag \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 Scalar> +static void ei_tridiagonal_qr_step(Scalar* diag, Scalar* subdiag, int n) +{ + Scalar td = (diag[n-2] - diag[n-1])*0.5; + Scalar e2 = ei_abs2(subdiag[n-2]); + Scalar mu = diag[n-1] - e2 / (td + (td>0 ? 1 : -1) * ei_sqrt(td*td + e2)); + Scalar x = diag[0] - mu; + Scalar z = subdiag[0]; + + for (int k = 0; k < n-1; ++k) + { + Scalar c, s; + ei_givens_rotation(x, z, c, s); + + // do T = G' T G + Scalar sdk = s * diag[k] + c * subdiag[k]; + Scalar dkp1 = s * subdiag[k] + c * diag[k+1]; + + diag[k] = c * (c * diag[k] - s * subdiag[k]) - s * (c * subdiag[k] - s * diag[k+1]); + diag[k+1] = s * sdk + c * dkp1; + subdiag[k] = c * sdk - s * dkp1; + + if (k > 0) + subdiag[k - 1] = c * subdiag[k-1] - s * z; + + x = subdiag[k]; + z = -s * subdiag[k+1]; + + if (k < n - 2) + subdiag[k + 1] = c * subdiag[k+1]; + } +} + +template<typename MatrixType> +void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix) +{ + assert(matrix.cols() == matrix.rows()); + int n = matrix.cols(); + m_eivalues.resize(n,1); + m_eivec = matrix; + + Tridiagonalization<MatrixType> tridiag(m_eivec); + RealVectorType& diag = m_eivalues; + RealVectorType subdiag(n-1); + diag = tridiag.diagonal(); + subdiag = tridiag.subDiagonal(); + + 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.coeffRef(start), &subdiag.coeffRef(start), end-start+1); + } + + std::cout << "ei values = " << m_eivalues.transpose() << "\n\n"; +} + +#endif // EIGEN_SELFADJOINTEIGENSOLVER_H diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h index a85068e38..96a4292b3 100755 --- a/Eigen/src/QR/Tridiagonalization.h +++ b/Eigen/src/QR/Tridiagonalization.h @@ -126,7 +126,9 @@ template<typename _MatrixType> class Tridiagonalization * \param matA the input selfadjoint matrix * \param hCoeffs returned Householder coefficients * - * The result is written in the lower triangular part of \a matA: + * The result is written in the lower triangular part of \a matA. + * + * Implemented from Golub's "Matrix Computations", algorithm 8.3.1. * * \sa packedMatrix() */ |