eigensolver_selfadjoint.cpp 11 KB

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