ground/gcs/src/libs/eigen/test/eigensolver_selfadjoint.cpp

   1 // This file is part of Eigen, a lightweight C++ template library
   2 // for linear algebra.
   3 //
   4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
   5 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
   6 //
   7 // This Source Code Form is subject to the terms of the Mozilla
   8 // Public License v. 2.0. If a copy of the MPL was not distributed
   9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
  10
  11 #include "main.h"
  12 #include "svd_fill.h"
  13 #include <limits>
  14 #include <Eigen/Eigenvalues>
  15 #include <Eigen/SparseCore>
  16
  17
  18 template<typename MatrixType> void selfadjointeigensolver_essential_check(const MatrixType& m)
  19 {
  20   typedef typename MatrixType::Scalar Scalar;
  21   typedef typename NumTraits<Scalar>::Real RealScalar;
  22   RealScalar eival_eps = numext::mini<RealScalar>(test_precision<RealScalar>(),  NumTraits<Scalar>::dummy_precision()*20000);
  23
  24   SelfAdjointEigenSolver<MatrixType> eiSymm(m);
  25   VERIFY_IS_EQUAL(eiSymm.info(), Success);
  26
  27   RealScalar scaling = m.cwiseAbs().maxCoeff();
  28
  29   if(scaling<(std::numeric_limits<RealScalar>::min)())
  30   {
  31     VERIFY(eiSymm.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
  32   }
  33   else
  34   {
  35     VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * eiSymm.eigenvectors())/scaling,
  36                      (eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal())/scaling);
  37   }
  38   VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues());
  39   VERIFY_IS_UNITARY(eiSymm.eigenvectors());
  40
  41   if(m.cols()<=4)
  42   {
  43     SelfAdjointEigenSolver<MatrixType> eiDirect;
  44     eiDirect.computeDirect(m);
  45     VERIFY_IS_EQUAL(eiDirect.info(), Success);
  46     if(! eiSymm.eigenvalues().isApprox(eiDirect.eigenvalues(), eival_eps) )
  47     {
  48       std::cerr << "reference eigenvalues: " << eiSymm.eigenvalues().transpose() << "\n"
  49                 << "obtained eigenvalues:  " << eiDirect.eigenvalues().transpose() << "\n"
  50                 << "diff:                  " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).transpose() << "\n"
  51                 << "error (eps):           " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).norm() / eiSymm.eigenvalues().norm() << "  (" << eival_eps << ")\n";
  52     }
  53     if(scaling<(std::numeric_limits<RealScalar>::min)())
  54     {
  55       VERIFY(eiDirect.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
  56     }
  57     else
  58     {
  59       VERIFY_IS_APPROX(eiSymm.eigenvalues()/scaling, eiDirect.eigenvalues()/scaling);
  60       VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * eiDirect.eigenvectors())/scaling,
  61                        (eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal())/scaling);
  62       VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues()/scaling, eiDirect.eigenvalues()/scaling);
  63     }
  64
  65     VERIFY_IS_UNITARY(eiDirect.eigenvectors());
  66   }
  67 }
  68
  69 template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
  70 {
  71   typedef typename MatrixType::Index Index;
  72   /* this test covers the following files:
  73      EigenSolver.h, SelfAdjointEigenSolver.h (and indirectly: Tridiagonalization.h)
  74   */
  75   Index rows = m.rows();
  76   Index cols = m.cols();
  77
  78   typedef typename MatrixType::Scalar Scalar;
  79   typedef typename NumTraits<Scalar>::Real RealScalar;
  80
  81   RealScalar largerEps = 10*test_precision<RealScalar>();
  82
  83   MatrixType a = MatrixType::Random(rows,cols);
  84   MatrixType a1 = MatrixType::Random(rows,cols);
  85   MatrixType symmA =  a.adjoint() * a + a1.adjoint() * a1;
  86   MatrixType symmC = symmA;
  87
  88   svd_fill_random(symmA,Symmetric);
  89
  90   symmA.template triangularView<StrictlyUpper>().setZero();
  91   symmC.template triangularView<StrictlyUpper>().setZero();
  92
  93   MatrixType b = MatrixType::Random(rows,cols);
  94   MatrixType b1 = MatrixType::Random(rows,cols);
  95   MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1;
  96   symmB.template triangularView<StrictlyUpper>().setZero();
  97
  98   CALL_SUBTEST( selfadjointeigensolver_essential_check(symmA) );
  99
 100   SelfAdjointEigenSolver<MatrixType> eiSymm(symmA);
 101   // generalized eigen pb
 102   GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmC, symmB);
 103
 104   SelfAdjointEigenSolver<MatrixType> eiSymmNoEivecs(symmA, false);
 105   VERIFY_IS_EQUAL(eiSymmNoEivecs.info(), Success);
 106   VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues());
 107
 108   // generalized eigen problem Ax = lBx
 109   eiSymmGen.compute(symmC, symmB,Ax_lBx);
 110   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
 111   VERIFY((symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()).isApprox(
 112           symmB.template selfadjointView<Lower>() * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 113
 114   // generalized eigen problem BAx = lx
 115   eiSymmGen.compute(symmC, symmB,BAx_lx);
 116   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
 117   VERIFY((symmB.template selfadjointView<Lower>() * (symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
 118          (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 119
 120   // generalized eigen problem ABx = lx
 121   eiSymmGen.compute(symmC, symmB,ABx_lx);
 122   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
 123   VERIFY((symmC.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
 124          (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 125
 126
 127   eiSymm.compute(symmC);
 128   MatrixType sqrtSymmA = eiSymm.operatorSqrt();
 129   VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), sqrtSymmA*sqrtSymmA);
 130   VERIFY_IS_APPROX(sqrtSymmA, symmC.template selfadjointView<Lower>()*eiSymm.operatorInverseSqrt());
 131
 132   MatrixType id = MatrixType::Identity(rows, cols);
 133   VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1));
 134
 135   SelfAdjointEigenSolver<MatrixType> eiSymmUninitialized;
 136   VERIFY_RAISES_ASSERT(eiSymmUninitialized.info());
 137   VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvalues());
 138   VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
 139   VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
 140   VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
 141
 142   eiSymmUninitialized.compute(symmA, false);
 143   VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
 144   VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
 145   VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
 146
 147   // test Tridiagonalization's methods
 148   Tridiagonalization<MatrixType> tridiag(symmC);
 149   VERIFY_IS_APPROX(tridiag.diagonal(), tridiag.matrixT().diagonal());
 150   VERIFY_IS_APPROX(tridiag.subDiagonal(), tridiag.matrixT().template diagonal<-1>());
 151   Matrix<RealScalar,Dynamic,Dynamic> T = tridiag.matrixT();
 152   if(rows>1 && cols>1) {
 153     // FIXME check that upper and lower part are 0:
 154     //VERIFY(T.topRightCorner(rows-2, cols-2).template triangularView<Upper>().isZero());
 155   }
 156   VERIFY_IS_APPROX(tridiag.diagonal(), T.diagonal());
 157   VERIFY_IS_APPROX(tridiag.subDiagonal(), T.template diagonal<1>());
 158   VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint());
 159   VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint());
 160
 161   // Test computation of eigenvalues from tridiagonal matrix
 162   if(rows > 1)
 163   {
 164     SelfAdjointEigenSolver<MatrixType> eiSymmTridiag;
 165     eiSymmTridiag.computeFromTridiagonal(tridiag.matrixT().diagonal(), tridiag.matrixT().diagonal(-1), ComputeEigenvectors);
 166     VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmTridiag.eigenvalues());
 167     VERIFY_IS_APPROX(tridiag.matrixT(), eiSymmTridiag.eigenvectors().real() * eiSymmTridiag.eigenvalues().asDiagonal() * eiSymmTridiag.eigenvectors().real().transpose());
 168   }
 169
 170   if (rows > 1 && rows < 20)
 171   {
 172     // Test matrix with NaN
 173     symmC(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
 174     SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmC);
 175     VERIFY_IS_EQUAL(eiSymmNaN.info(), NoConvergence);
 176   }
 177
 178   // regression test for bug 1098
 179   {
 180     SelfAdjointEigenSolver<MatrixType> eig(a.adjoint() * a);
 181     eig.compute(a.adjoint() * a);
 182   }
 183
 184   // regression test for bug 478
 185   {
 186     a.setZero();
 187     SelfAdjointEigenSolver<MatrixType> ei3(a);
 188     VERIFY_IS_EQUAL(ei3.info(), Success);
 189     VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
 190     VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
 191   }
 192 }
 193
 194 template<int>
 195 void bug_854()
 196 {
 197   Matrix3d m;
 198   m << 850.961, 51.966, 0,
 199        51.966, 254.841, 0,
 200             0,       0, 0;
 201   selfadjointeigensolver_essential_check(m);
 202 }
 203
 204 template<int>
 205 void bug_1014()
 206 {
 207   Matrix3d m;
 208   m <<        0.11111111111111114658, 0, 0,
 209        0,     0.11111111111111109107, 0,
 210        0, 0,  0.11111111111111107719;
 211   selfadjointeigensolver_essential_check(m);
 212 }
 213
 214 template<int>
 215 void bug_1225()
 216 {
 217   Matrix3d m1, m2;
 218   m1.setRandom();
 219   m1 = m1*m1.transpose();
 220   m2 = m1.triangularView<Upper>();
 221   SelfAdjointEigenSolver<Matrix3d> eig1(m1);
 222   SelfAdjointEigenSolver<Matrix3d> eig2(m2.selfadjointView<Upper>());
 223   VERIFY_IS_APPROX(eig1.eigenvalues(), eig2.eigenvalues());
 224 }
 225
 226 template<int>
 227 void bug_1204()
 228 {
 229   SparseMatrix<double> A(2,2);
 230   A.setIdentity();
 231   SelfAdjointEigenSolver<Eigen::SparseMatrix<double> > eig(A);
 232 }
 233
 234 void test_eigensolver_selfadjoint()
 235 {
 236   int s = 0;
 237   for(int i = 0; i < g_repeat; i++) {
 238     // trivial test for 1x1 matrices:
 239     CALL_SUBTEST_1( selfadjointeigensolver(Matrix<float, 1, 1>()));
 240     CALL_SUBTEST_1( selfadjointeigensolver(Matrix<double, 1, 1>()));
 241     // very important to test 3x3 and 2x2 matrices since we provide special paths for them
 242     CALL_SUBTEST_12( selfadjointeigensolver(Matrix2f()) );
 243     CALL_SUBTEST_12( selfadjointeigensolver(Matrix2d()) );
 244     CALL_SUBTEST_13( selfadjointeigensolver(Matrix3f()) );
 245     CALL_SUBTEST_13( selfadjointeigensolver(Matrix3d()) );
 246     CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) );
 247
 248     s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
 249     CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(s,s)) );
 250     CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(s,s)) );
 251     CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(s,s)) );
 252     CALL_SUBTEST_9( selfadjointeigensolver(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(s,s)) );
 253     TEST_SET_BUT_UNUSED_VARIABLE(s)
 254
 255     // some trivial but implementation-wise tricky cases
 256     CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(1,1)) );
 257     CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(2,2)) );
 258     CALL_SUBTEST_6( selfadjointeigensolver(Matrix<double,1,1>()) );
 259     CALL_SUBTEST_7( selfadjointeigensolver(Matrix<double,2,2>()) );
 260   }
 261
 262   CALL_SUBTEST_13( bug_854<0>() );
 263   CALL_SUBTEST_13( bug_1014<0>() );
 264   CALL_SUBTEST_13( bug_1204<0>() );
 265   CALL_SUBTEST_13( bug_1225<0>() );
 266
 267   // Test problem size constructors
 268   s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
 269   CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf> tmp1(s));
 270   CALL_SUBTEST_8(Tridiagonalization<MatrixXf> tmp2(s));
 271
 272   TEST_SET_BUT_UNUSED_VARIABLE(s)
 273 }
 274