eigensolver_selfadjoint.cpp File Reference

#include "main.h"
#include "svd_fill.h"
#include <limits>
#include <Eigen/Eigenvalues>
#include <Eigen/SparseCore>

Functions
template<typename MatrixType >
void	selfadjointeigensolver_essential_check (const MatrixType &m)
template<typename MatrixType >
void	selfadjointeigensolver (const MatrixType &m)
template<int >
void	bug_854 ()
template<int >
void	bug_1014 ()
template<int >
void	bug_1225 ()
template<int >
void	bug_1204 ()
	EIGEN_DECLARE_TEST (eigensolver_selfadjoint)

Function Documentation

bug_1014()

template<int >

void bug_1014

(

)

                {
  Matrix3d m;
  m << 0.11111111111111114658, 0, 0, 0, 0.11111111111111109107, 0, 0, 0, 0.11111111111111107719;
  selfadjointeigensolver_essential_check(m);
}

References m, and selfadjointeigensolver_essential_check().

bug_1204()

template<int >

void bug_1204

(

)

                {
  SparseMatrix<double> A(2, 2);
  A.setIdentity();
  SelfAdjointEigenSolver<Eigen::SparseMatrix<double> > eig(A);
}

bug_1225()

template<int >

void bug_1225

(

)

                {
  Matrix3d m1, m2;
  m1.setRandom();
  m1 = m1 * m1.transpose();
  m2 = m1.triangularView<Upper>();
  SelfAdjointEigenSolver<Matrix3d> eig1(m1);
  SelfAdjointEigenSolver<Matrix3d> eig2(m2.selfadjointView<Upper>());
  VERIFY_IS_APPROX(eig1.eigenvalues(), eig2.eigenvalues());
}

References Eigen::SelfAdjointEigenSolver< MatrixType_ >::eigenvalues(), m1, m2(), Eigen::Upper, and VERIFY_IS_APPROX.

bug_854()

template<int >

void bug_854

(

)

               {
  Matrix3d m;
  m << 850.961, 51.966, 0, 51.966, 254.841, 0, 0, 0, 0;
  selfadjointeigensolver_essential_check(m);
}

References m, and selfadjointeigensolver_essential_check().

EIGEN_DECLARE_TEST()

EIGEN_DECLARE_TEST

(

eigensolver_selfadjoint

)

                                            {
  int s = 0;
  for (int i = 0; i < g_repeat; i++) {
    // trivial test for 1x1 matrices:
    CALL_SUBTEST_1(selfadjointeigensolver(Matrix<float, 1, 1>()));
    CALL_SUBTEST_1(selfadjointeigensolver(Matrix<double, 1, 1>()));
    CALL_SUBTEST_1(selfadjointeigensolver(Matrix<std::complex<double>, 1, 1>()));
 
    // very important to test 3x3 and 2x2 matrices since we provide special paths for them
    CALL_SUBTEST_12(selfadjointeigensolver(Matrix2f()));
    CALL_SUBTEST_12(selfadjointeigensolver(Matrix2d()));
    CALL_SUBTEST_12(selfadjointeigensolver(Matrix2cd()));
    CALL_SUBTEST_13(selfadjointeigensolver(Matrix3f()));
    CALL_SUBTEST_13(selfadjointeigensolver(Matrix3d()));
    CALL_SUBTEST_13(selfadjointeigensolver(Matrix3cd()));
    CALL_SUBTEST_2(selfadjointeigensolver(Matrix4d()));
    CALL_SUBTEST_2(selfadjointeigensolver(Matrix4cd()));
 
    s = internal::random<int>(1, EIGEN_TEST_MAX_SIZE / 4);
    CALL_SUBTEST_3(selfadjointeigensolver(MatrixXf(s, s)));
    CALL_SUBTEST_4(selfadjointeigensolver(MatrixXd(s, s)));
    CALL_SUBTEST_5(selfadjointeigensolver(MatrixXcd(s, s)));
    CALL_SUBTEST_9(selfadjointeigensolver(Matrix<std::complex<double>, Dynamic, Dynamic, RowMajor>(s, s)));
    TEST_SET_BUT_UNUSED_VARIABLE(s)
 
    // some trivial but implementation-wise tricky cases
    CALL_SUBTEST_4(selfadjointeigensolver(MatrixXd(1, 1)));
    CALL_SUBTEST_4(selfadjointeigensolver(MatrixXd(2, 2)));
    CALL_SUBTEST_5(selfadjointeigensolver(MatrixXcd(1, 1)));
    CALL_SUBTEST_5(selfadjointeigensolver(MatrixXcd(2, 2)));
    CALL_SUBTEST_6(selfadjointeigensolver(Matrix<double, 1, 1>()));
    CALL_SUBTEST_7(selfadjointeigensolver(Matrix<double, 2, 2>()));
  }
 
  CALL_SUBTEST_13(bug_854<0>());
  CALL_SUBTEST_13(bug_1014<0>());
  CALL_SUBTEST_13(bug_1204<0>());
  CALL_SUBTEST_13(bug_1225<0>());
 
  // Test problem size constructors
  s = internal::random<int>(1, EIGEN_TEST_MAX_SIZE / 4);
  CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf> tmp1(s));
  CALL_SUBTEST_8(Tridiagonalization<MatrixXf> tmp2(s));
 
  TEST_SET_BUT_UNUSED_VARIABLE(s)
}

References CALL_SUBTEST_1, CALL_SUBTEST_12, CALL_SUBTEST_13, CALL_SUBTEST_2, CALL_SUBTEST_3, CALL_SUBTEST_4, CALL_SUBTEST_5, CALL_SUBTEST_6, CALL_SUBTEST_7, CALL_SUBTEST_8, CALL_SUBTEST_9, Eigen::Dynamic, EIGEN_TEST_MAX_SIZE, Eigen::g_repeat, i, Eigen::RowMajor, s, selfadjointeigensolver(), and TEST_SET_BUT_UNUSED_VARIABLE.

selfadjointeigensolver()

template<typename MatrixType >

void selfadjointeigensolver

(

const MatrixType &

m

)

                                                 {
  /* this test covers the following files:
     EigenSolver.h, SelfAdjointEigenSolver.h (and indirectly: Tridiagonalization.h)
  */
  Index rows = m.rows();
  Index cols = m.cols();
 
  typedef typename MatrixType::Scalar Scalar;
  typedef typename NumTraits<Scalar>::Real RealScalar;
 
  RealScalar largerEps = 10 * test_precision<RealScalar>();
 
  MatrixType a = MatrixType::Random(rows, cols);
  MatrixType a1 = MatrixType::Random(rows, cols);
  MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
  MatrixType symmC = symmA;
 
  svd_fill_random(symmA, Symmetric);
 
  symmA.template triangularView<StrictlyUpper>().setZero();
  symmC.template triangularView<StrictlyUpper>().setZero();
 
  MatrixType b = MatrixType::Random(rows, cols);
  MatrixType b1 = MatrixType::Random(rows, cols);
  MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1;
  symmB.template triangularView<StrictlyUpper>().setZero();
 
  CALL_SUBTEST(selfadjointeigensolver_essential_check(symmA));
 
  SelfAdjointEigenSolver<MatrixType> eiSymm(symmA);
  // generalized eigen pb
  GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmC, symmB);
 
  SelfAdjointEigenSolver<MatrixType> eiSymmNoEivecs(symmA, false);
  VERIFY_IS_EQUAL(eiSymmNoEivecs.info(), Success);
  VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues());
 
  // generalized eigen problem Ax = lBx
  eiSymmGen.compute(symmC, symmB, Ax_lBx);
  VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
  VERIFY((symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())
             .isApprox(symmB.template selfadjointView<Lower>() *
                           (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()),
                       largerEps));
 
  // generalized eigen problem BAx = lx
  eiSymmGen.compute(symmC, symmB, BAx_lx);
  VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
  VERIFY(
      (symmB.template selfadjointView<Lower>() * (symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()))
          .isApprox((eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 
  // generalized eigen problem ABx = lx
  eiSymmGen.compute(symmC, symmB, ABx_lx);
  VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
  VERIFY(
      (symmC.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()))
          .isApprox((eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 
  eiSymm.compute(symmC);
  MatrixType sqrtSymmA = eiSymm.operatorSqrt();
  VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), sqrtSymmA * sqrtSymmA);
  VERIFY_IS_APPROX(sqrtSymmA, symmC.template selfadjointView<Lower>() * eiSymm.operatorInverseSqrt());
 
  MatrixType id = MatrixType::Identity(rows, cols);
  VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1));
 
  SelfAdjointEigenSolver<MatrixType> eiSymmUninitialized;
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.info());
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvalues());
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
 
  eiSymmUninitialized.compute(symmA, false);
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
 
  // test Tridiagonalization's methods
  Tridiagonalization<MatrixType> tridiag(symmC);
  VERIFY_IS_APPROX(tridiag.diagonal(), tridiag.matrixT().diagonal());
  VERIFY_IS_APPROX(tridiag.subDiagonal(), tridiag.matrixT().template diagonal<-1>());
  Matrix<RealScalar, Dynamic, Dynamic> T = tridiag.matrixT();
  if (rows > 1 && cols > 1) {
    // FIXME check that upper and lower part are 0:
    // VERIFY(T.topRightCorner(rows-2, cols-2).template triangularView<Upper>().isZero());
  }
  VERIFY_IS_APPROX(tridiag.diagonal(), T.diagonal());
  VERIFY_IS_APPROX(tridiag.subDiagonal(), T.template diagonal<1>());
  VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()),
                   tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint());
  VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()),
                   tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint());
 
  // Test computation of eigenvalues from tridiagonal matrix
  if (rows > 1) {
    SelfAdjointEigenSolver<MatrixType> eiSymmTridiag;
    eiSymmTridiag.computeFromTridiagonal(tridiag.matrixT().diagonal(), tridiag.matrixT().diagonal(-1),
                                         ComputeEigenvectors);
    VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmTridiag.eigenvalues());
    VERIFY_IS_APPROX(tridiag.matrixT(), eiSymmTridiag.eigenvectors().real() * eiSymmTridiag.eigenvalues().asDiagonal() *
                                            eiSymmTridiag.eigenvectors().real().transpose());
  }
 
  if (rows > 1 && rows < 20) {
    // Test matrix with NaN
    symmC(0, 0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
    SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmC);
    VERIFY_IS_EQUAL(eiSymmNaN.info(), NoConvergence);
  }
 
  // regression test for bug 1098
  {
    SelfAdjointEigenSolver<MatrixType> eig(a.adjoint() * a);
    eig.compute(a.adjoint() * a);
  }
 
  // regression test for bug 478
  {
    a.setZero();
    SelfAdjointEigenSolver<MatrixType> ei3(a);
    VERIFY_IS_EQUAL(ei3.info(), Success);
    VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(), RealScalar(1));
    VERIFY((ei3.eigenvectors().transpose() * ei3.eigenvectors().transpose()).eval().isIdentity());
  }
}

References a, Eigen::ABx_lx, Eigen::HouseholderSequence< VectorsType, CoeffsType, Side >::adjoint(), Eigen::Ax_lBx, b, Eigen::BAx_lx, CALL_SUBTEST, cols, Eigen::SelfAdjointEigenSolver< MatrixType_ >::compute(), Eigen::GeneralizedSelfAdjointEigenSolver< MatrixType_ >::compute(), Eigen::ComputeEigenvectors, Eigen::SelfAdjointEigenSolver< MatrixType_ >::computeFromTridiagonal(), Eigen::Tridiagonalization< MatrixType_ >::diagonal(), Eigen::SelfAdjointEigenSolver< MatrixType_ >::eigenvalues(), Eigen::SelfAdjointEigenSolver< MatrixType_ >::eigenvectors(), Eigen::SelfAdjointEigenSolver< MatrixType_ >::info(), m, Eigen::Tridiagonalization< MatrixType_ >::matrixQ(), Eigen::Tridiagonalization< MatrixType_ >::matrixT(), Eigen::NoConvergence, Eigen::SelfAdjointEigenSolver< MatrixType_ >::operatorInverseSqrt(), Eigen::SelfAdjointEigenSolver< MatrixType_ >::operatorSqrt(), rows, selfadjointeigensolver_essential_check(), Eigen::Tridiagonalization< MatrixType_ >::subDiagonal(), Eigen::Success, svd_fill_random(), Eigen::Symmetric, VERIFY, VERIFY_IS_APPROX, VERIFY_IS_EQUAL, VERIFY_IS_MUCH_SMALLER_THAN, and VERIFY_RAISES_ASSERT.

Referenced by EIGEN_DECLARE_TEST().

selfadjointeigensolver_essential_check()

template<typename MatrixType >

void selfadjointeigensolver_essential_check

(

const MatrixType &

m

)

                                                                 {
  typedef typename MatrixType::Scalar Scalar;
  typedef typename NumTraits<Scalar>::Real RealScalar;
  RealScalar eival_eps =
      numext::mini<RealScalar>(test_precision<RealScalar>(), NumTraits<Scalar>::dummy_precision() * 20000);
 
  SelfAdjointEigenSolver<MatrixType> eiSymm(m);
  VERIFY_IS_EQUAL(eiSymm.info(), Success);
 
  RealScalar scaling = m.cwiseAbs().maxCoeff();
  RealScalar unitary_error_factor = RealScalar(16);
 
  if (scaling < (std::numeric_limits<RealScalar>::min)()) {
    VERIFY(eiSymm.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
  } else {
    VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * eiSymm.eigenvectors()) / scaling,
                     (eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal()) / scaling);
  }
  VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues());
  VERIFY(eiSymm.eigenvectors().isUnitary(test_precision<RealScalar>() * unitary_error_factor));
 
  if (m.cols() <= 4) {
    SelfAdjointEigenSolver<MatrixType> eiDirect;
    eiDirect.computeDirect(m);
    VERIFY_IS_EQUAL(eiDirect.info(), Success);
    if (!eiSymm.eigenvalues().isApprox(eiDirect.eigenvalues(), eival_eps)) {
      std::cerr << "reference eigenvalues: " << eiSymm.eigenvalues().transpose() << "\n"
                << "obtained eigenvalues:  " << eiDirect.eigenvalues().transpose() << "\n"
                << "diff:                  " << (eiSymm.eigenvalues() - eiDirect.eigenvalues()).transpose() << "\n"
                << "error (eps):           "
                << (eiSymm.eigenvalues() - eiDirect.eigenvalues()).norm() / eiSymm.eigenvalues().norm() << "  ("
                << eival_eps << ")\n";
    }
    if (scaling < (std::numeric_limits<RealScalar>::min)()) {
      VERIFY(eiDirect.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
    } else {
      VERIFY_IS_APPROX(eiSymm.eigenvalues() / scaling, eiDirect.eigenvalues() / scaling);
      VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * eiDirect.eigenvectors()) / scaling,
                       (eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal()) / scaling);
      VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues() / scaling, eiDirect.eigenvalues() / scaling);
    }
 
    VERIFY(eiDirect.eigenvectors().isUnitary(test_precision<RealScalar>() * unitary_error_factor));
  }
}

References Eigen::SelfAdjointEigenSolver< MatrixType_ >::computeDirect(), Eigen::SelfAdjointEigenSolver< MatrixType_ >::eigenvalues(), Eigen::SelfAdjointEigenSolver< MatrixType_ >::eigenvectors(), Eigen::SelfAdjointEigenSolver< MatrixType_ >::info(), m, min, Eigen::Success, anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose(), VERIFY, VERIFY_IS_APPROX, and VERIFY_IS_EQUAL.

Referenced by bug_1014(), bug_854(), and selfadjointeigensolver().