From 2de8da70fd0b35849845dc76b2741bb0689f0643 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 12 Dec 2018 17:30:08 +0100 Subject: [PATCH] bug #1557: fix RealSchur and EigenSolver for matrices with only zeros on the diagonal. --- Eigen/src/Eigenvalues/RealSchur.h | 15 +++++-- test/eigensolver_generic.cpp | 74 +++++++++++++++++++++++++++---- 2 files changed, 77 insertions(+), 12 deletions(-) diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h index aca8a8279..8dbd9e314 100644 --- a/Eigen/src/Eigenvalues/RealSchur.h +++ b/Eigen/src/Eigenvalues/RealSchur.h @@ -236,7 +236,7 @@ template class RealSchur typedef Matrix Vector3s; Scalar computeNormOfT(); - Index findSmallSubdiagEntry(Index iu); + Index findSmallSubdiagEntry(Index iu, const Scalar& considerAsZero); void splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift); void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo); void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector); @@ -307,12 +307,16 @@ RealSchur& RealSchur::computeFromHessenberg(const HessMa Index totalIter = 0; // iteration count for whole matrix Scalar exshift(0); // sum of exceptional shifts Scalar norm = computeNormOfT(); + // sub-diagonal entries smaller than considerAsZero will be treated as zero. + // We use eps^2 to enable more precision in small eigenvalues. + Scalar considerAsZero = numext::maxi( norm * numext::abs2(NumTraits::epsilon()), + (std::numeric_limits::min)() ); if(norm!=Scalar(0)) { while (iu >= 0) { - Index il = findSmallSubdiagEntry(iu); + Index il = findSmallSubdiagEntry(iu,considerAsZero); // Check for convergence if (il == iu) // One root found @@ -369,14 +373,17 @@ inline typename MatrixType::Scalar RealSchur::computeNormOfT() /** \internal Look for single small sub-diagonal element and returns its index */ template -inline Index RealSchur::findSmallSubdiagEntry(Index iu) +inline Index RealSchur::findSmallSubdiagEntry(Index iu, const Scalar& considerAsZero) { using std::abs; Index res = iu; while (res > 0) { Scalar s = abs(m_matT.coeff(res-1,res-1)) + abs(m_matT.coeff(res,res)); - if (abs(m_matT.coeff(res,res-1)) <= NumTraits::epsilon() * s) + + s = numext::maxi(s * NumTraits::epsilon(), considerAsZero); + + if (abs(m_matT.coeff(res,res-1)) <= s) break; res--; } diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp index e0e435151..086ecdf5e 100644 --- a/test/eigensolver_generic.cpp +++ b/test/eigensolver_generic.cpp @@ -12,6 +12,21 @@ #include #include +template +void check_eigensolver_for_given_mat(const EigType &eig, const MatType& a) +{ + typedef typename NumTraits::Real RealScalar; + typedef Matrix RealVectorType; + typedef typename std::complex Complex; + Index n = a.rows(); + VERIFY_IS_EQUAL(eig.info(), Success); + VERIFY_IS_APPROX(a * eig.pseudoEigenvectors(), eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()); + VERIFY_IS_APPROX(a.template cast() * eig.eigenvectors(), + eig.eigenvectors() * eig.eigenvalues().asDiagonal()); + VERIFY_IS_APPROX(eig.eigenvectors().colwise().norm(), RealVectorType::Ones(n).transpose()); + VERIFY_IS_APPROX(a.eigenvalues(), eig.eigenvalues()); +} + template void eigensolver(const MatrixType& m) { /* this test covers the following files: @@ -22,8 +37,7 @@ template void eigensolver(const MatrixType& m) typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Matrix RealVectorType; - typedef typename std::complex::Real> Complex; + typedef typename std::complex Complex; MatrixType a = MatrixType::Random(rows,cols); MatrixType a1 = MatrixType::Random(rows,cols); @@ -36,12 +50,7 @@ template void eigensolver(const MatrixType& m) (ei0.pseudoEigenvectors().template cast()) * (ei0.eigenvalues().asDiagonal())); EigenSolver ei1(a); - VERIFY_IS_EQUAL(ei1.info(), Success); - VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix()); - VERIFY_IS_APPROX(a.template cast() * ei1.eigenvectors(), - ei1.eigenvectors() * ei1.eigenvalues().asDiagonal()); - VERIFY_IS_APPROX(ei1.eigenvectors().colwise().norm(), RealVectorType::Ones(rows).transpose()); - VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues()); + CALL_SUBTEST( check_eigensolver_for_given_mat(ei1,a) ); EigenSolver ei2; ei2.setMaxIterations(RealSchur::m_maxIterationsPerRow * rows).compute(a); @@ -100,6 +109,19 @@ template void eigensolver_verify_assert(const MatrixType& m VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors()); } + +template +Matrix +make_companion(const CoeffType& coeffs) +{ + Index n = coeffs.size()-1; + Matrix res(n,n); + res.setZero(); + res.row(0) = -coeffs.tail(n) / coeffs(0); + res.diagonal(-1).setOnes(); + return res; +} + template void eigensolver_generic_extra() { @@ -126,6 +148,42 @@ void eigensolver_generic_extra() VERIFY_IS_APPROX((a * eig.eigenvectors()).norm()+1., 1.); VERIFY_IS_APPROX((eig.eigenvectors() * eig.eigenvalues().asDiagonal()).norm()+1., 1.); } + + // regression test for bug 933 + { + { + VectorXd coeffs(5); coeffs << 1, -3, -175, -225, 2250; + MatrixXd C = make_companion(coeffs); + EigenSolver eig(C); + CALL_SUBTEST( check_eigensolver_for_given_mat(eig,C) ); + } + { + // this test is tricky because it requires high accuracy in smallest eigenvalues + VectorXd coeffs(5); coeffs << 6.154671e-15, -1.003870e-10, -9.819570e-01, 3.995715e+03, 2.211511e+08; + MatrixXd C = make_companion(coeffs); + EigenSolver eig(C); + CALL_SUBTEST( check_eigensolver_for_given_mat(eig,C) ); + Index n = C.rows(); + for(Index i=0;i Complex; + MatrixXcd ac = C.cast(); + ac.diagonal().array() -= eig.eigenvalues()(i); + VectorXd sv = ac.jacobiSvd().singularValues(); + // comparing to sv(0) is not enough here to catch the "bug", + // the hard-coded 1.0 is important! + VERIFY_IS_MUCH_SMALLER_THAN(sv(n-1), 1.0); + } + } + } + // regression test for bug 1557 + { + // this test is interesting because it contains zeros on the diagonal. + MatrixXd A_bug1557(3,3); + A_bug1557 << 0, 0, 0, 1, 0, 0.5887907064808635127, 0, 1, 0; + EigenSolver eig(A_bug1557); + CALL_SUBTEST( check_eigensolver_for_given_mat(eig,A_bug1557) ); + } } EIGEN_DECLARE_TEST(eigensolver_generic)