From 26d7dad1389ca3eb46127dbe8c76c8359e6ded4b Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 21 Jul 2011 19:07:52 +0200 Subject: [PATCH] add a computeDirect method to SelfAdjointEigenSolver for fast eigen decomposition --- Eigen/Eigenvalues | 1 + .../src/Eigenvalues/SelfAdjointEigenSolver.h | 156 +++++++++++++++++- test/eigensolver_selfadjoint.cpp | 9 +- 3 files changed, 164 insertions(+), 2 deletions(-) diff --git a/Eigen/Eigenvalues b/Eigen/Eigenvalues index 250c0f466..bfcd83dd3 100644 --- a/Eigen/Eigenvalues +++ b/Eigen/Eigenvalues @@ -9,6 +9,7 @@ #include "Jacobi" #include "Householder" #include "LU" +#include "Geometry" namespace Eigen { diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 965dda88b..0af11ba1a 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -32,6 +32,10 @@ template class GeneralizedSelfAdjointEigenSolver; +namespace internal { +template struct direct_selfadjoint_eigenvalues; +} + /** \eigenvalues_module \ingroup Eigenvalues_Module * * @@ -86,7 +90,7 @@ template class SelfAdjointEigenSolver Options = MatrixType::Options, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; - + /** \brief Scalar type for matrices of type \p _MatrixType. */ typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; @@ -98,6 +102,8 @@ template class SelfAdjointEigenSolver * complex. */ typedef typename NumTraits::Real RealScalar; + + friend struct internal::direct_selfadjoint_eigenvalues::IsComplex>; /** \brief Type for vector of eigenvalues as returned by eigenvalues(). * @@ -198,6 +204,22 @@ template class SelfAdjointEigenSolver * \sa SelfAdjointEigenSolver(const MatrixType&, int) */ SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors); + + /** \brief Computes eigendecomposition of given matrix using a direct algorithm + * + * This is a variant of compute(const MatrixType&, int options) which + * directly solves the underlying polynomial equation. + * + * Currently only 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d). + * + * This method is usually significantly faster than the QR algorithm + * but it might also be less accurate. It is also worth noting that + * for 3x3 matrices it involves trigonometric operations which are + * not necessarily available for all scalar types. + * + * \sa compute(const MatrixType&, int options) + */ + SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors); /** \brief Returns the eigenvectors of given matrix. * @@ -466,6 +488,138 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver return *this; } + +namespace internal { + +template struct direct_selfadjoint_eigenvalues +{ + inline static void run(SolverType& eig, const typename SolverType::MatrixType& A, int options) + { eig.compute(A,options); } +}; + +template struct direct_selfadjoint_eigenvalues +{ + typedef typename SolverType::MatrixType MatrixType; + typedef typename SolverType::Scalar Scalar; + + template + inline static void computeRoots(const MatrixType& m, Roots& roots) + { + using std::sqrt; + using std::atan2; + using std::cos; + using std::sin; + const Scalar s_inv3 = 1.0/3.0; + const Scalar s_sqrt3 = sqrt(Scalar(3.0)); + + // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The + // eigenvalues are the roots to this equation, all guaranteed to be + // real-valued, because the matrix is symmetric. + Scalar c0 = m(0,0)*m(1,1)*m(2,2) + Scalar(2)*m(1,0)*m(2,0)*m(2,1) - m(0,0)*m(2,1)*m(2,1) - m(1,1)*m(2,0)*m(2,0) - m(2,2)*m(1,0)*m(1,0); + Scalar c1 = m(0,0)*m(1,1) - m(1,0)*m(1,0) + m(0,0)*m(2,2) - m(2,0)*m(2,0) + m(1,1)*m(2,2) - m(2,1)*m(2,1); + Scalar c2 = m(0,0) + m(1,1) + m(2,2); + + // Construct the parameters used in classifying the roots of the equation + // and in solving the equation for the roots in closed form. + Scalar c2_over_3 = c2*s_inv3; + Scalar a_over_3 = (c1 - c2*c2_over_3)*s_inv3; + if (a_over_3 > Scalar(0)) + a_over_3 = Scalar(0); + + Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1)); + + Scalar q = half_b*half_b + a_over_3*a_over_3*a_over_3; + if (q > Scalar(0)) + q = Scalar(0); + + // Compute the eigenvalues by solving for the roots of the polynomial. + Scalar rho = sqrt(-a_over_3); + Scalar theta = atan2(sqrt(-q),half_b)*s_inv3; + Scalar cos_theta = cos(theta); + Scalar sin_theta = sin(theta); + roots(0) = c2_over_3 + Scalar(2)*rho*cos_theta; + roots(1) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta); + roots(2) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta); + + // Sort in increasing order. + if (roots(0) >= roots(1)) + std::swap(roots(0),roots(1)); + if (roots(1) >= roots(2)) + { + std::swap(roots(1),roots(2)); + if (roots(0) >= roots(1)) + std::swap(roots(0),roots(1)); + } + } + + inline static void run(SolverType& solver, const MatrixType& mat, int options) + { + eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows()); + eigen_assert((options&~(EigVecMask|GenEigMask))==0 + && (options&EigVecMask)!=EigVecMask + && "invalid option parameter"); + bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; + + MatrixType& eivecs = solver.m_eivec; + typename SolverType::RealVectorType& eivals = solver.m_eivalues; + + // map the matrix coefficients to [-1:1] to avoid over- and underflow. + Scalar scale = mat.cwiseAbs().maxCoeff(); + scale = std::max(scale,Scalar(1)); + MatrixType scaledMat = mat / scale; + + // Compute the eigenvalues + computeRoots(scaledMat,eivals); + + // compute the eigen vectors + if(computeEigenvectors) + { + if((eivals(2)-eivals(0))<=Eigen::NumTraits::epsilon()) + { + eivecs.setIdentity(); + } + else + { + scaledMat = scaledMat.template selfadjointView(); + MatrixType tmp; + tmp = scaledMat; + tmp.diagonal().array () -= eivals (2); + eivecs.col (2) = tmp.row (0).cross (tmp.row (1)).normalized (); + + tmp = scaledMat; + tmp.diagonal ().array () -= eivals (1); + eivecs.col(1) = tmp.row (0).cross(tmp.row (1)); + Scalar n1 = eivecs.col(1).norm(); + if(n1<=Eigen::NumTraits::epsilon()) + eivecs.col(1) = eivecs.col(2).unitOrthogonal(); + else + eivecs.col(1) /= n1; + + // make sure that evecs[1] is orthogonal to evecs[2] + eivecs.col(1) = eivecs.col(2).cross(eivecs.col(1).cross(eivecs.col(2))).normalized(); + eivecs.col(0) = eivecs.col(2).cross(eivecs.col(1)); + } + } + + // Rescale back to the original size. + eivals *= scale; + + solver.m_info = Success; + solver.m_isInitialized = true; + solver.m_eigenvectorsOk = computeEigenvectors; + } +}; + +} + +template +SelfAdjointEigenSolver& SelfAdjointEigenSolver +::computeDirect(const MatrixType& matrix, int options) +{ + internal::direct_selfadjoint_eigenvalues::IsComplex>::run(*this,matrix,options); + return *this; +} + namespace internal { template static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n) diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 2078b93bc..a98f1dbf8 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -59,6 +59,8 @@ template void selfadjointeigensolver(const MatrixType& m) symmB.template triangularView().setZero(); SelfAdjointEigenSolver eiSymm(symmA); + SelfAdjointEigenSolver eiDirect; + eiDirect.computeDirect(symmA); // generalized eigen pb GeneralizedSelfAdjointEigenSolver eiSymmGen(symmA, symmB); @@ -112,11 +114,16 @@ template void selfadjointeigensolver(const MatrixType& m) VERIFY((symmA.template selfadjointView() * eiSymm.eigenvectors()).isApprox( eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps)); VERIFY_IS_APPROX(symmA.template selfadjointView().eigenvalues(), eiSymm.eigenvalues()); + + VERIFY_IS_EQUAL(eiDirect.info(), Success); + VERIFY((symmA.template selfadjointView() * eiDirect.eigenvectors()).isApprox( + eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal(), largerEps)); + VERIFY_IS_APPROX(symmA.template selfadjointView().eigenvalues(), eiDirect.eigenvalues()); SelfAdjointEigenSolver eiSymmNoEivecs(symmA, false); VERIFY_IS_EQUAL(eiSymmNoEivecs.info(), Success); VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues()); - + // generalized eigen problem Ax = lBx eiSymmGen.compute(symmA, symmB,Ax_lBx); VERIFY_IS_EQUAL(eiSymmGen.info(), Success);