eigensolver_selfadjoint.cpp (11419B)
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 /* this test covers the following files: 72 EigenSolver.h, SelfAdjointEigenSolver.h (and indirectly: Tridiagonalization.h) 73 */ 74 Index rows = m.rows(); 75 Index cols = m.cols(); 76 77 typedef typename MatrixType::Scalar Scalar; 78 typedef typename NumTraits<Scalar>::Real RealScalar; 79 80 RealScalar largerEps = 10*test_precision<RealScalar>(); 81 82 MatrixType a = MatrixType::Random(rows,cols); 83 MatrixType a1 = MatrixType::Random(rows,cols); 84 MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1; 85 MatrixType symmC = symmA; 86 87 svd_fill_random(symmA,Symmetric); 88 89 symmA.template triangularView<StrictlyUpper>().setZero(); 90 symmC.template triangularView<StrictlyUpper>().setZero(); 91 92 MatrixType b = MatrixType::Random(rows,cols); 93 MatrixType b1 = MatrixType::Random(rows,cols); 94 MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1; 95 symmB.template triangularView<StrictlyUpper>().setZero(); 96 97 CALL_SUBTEST( selfadjointeigensolver_essential_check(symmA) ); 98 99 SelfAdjointEigenSolver<MatrixType> eiSymm(symmA); 100 // generalized eigen pb 101 GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmC, symmB); 102 103 SelfAdjointEigenSolver<MatrixType> eiSymmNoEivecs(symmA, false); 104 VERIFY_IS_EQUAL(eiSymmNoEivecs.info(), Success); 105 VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues()); 106 107 // generalized eigen problem Ax = lBx 108 eiSymmGen.compute(symmC, symmB,Ax_lBx); 109 VERIFY_IS_EQUAL(eiSymmGen.info(), Success); 110 VERIFY((symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()).isApprox( 111 symmB.template selfadjointView<Lower>() * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps)); 112 113 // generalized eigen problem BAx = lx 114 eiSymmGen.compute(symmC, symmB,BAx_lx); 115 VERIFY_IS_EQUAL(eiSymmGen.info(), Success); 116 VERIFY((symmB.template selfadjointView<Lower>() * (symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox( 117 (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps)); 118 119 // generalized eigen problem ABx = lx 120 eiSymmGen.compute(symmC, symmB,ABx_lx); 121 VERIFY_IS_EQUAL(eiSymmGen.info(), Success); 122 VERIFY((symmC.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox( 123 (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps)); 124 125 126 eiSymm.compute(symmC); 127 MatrixType sqrtSymmA = eiSymm.operatorSqrt(); 128 VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), sqrtSymmA*sqrtSymmA); 129 VERIFY_IS_APPROX(sqrtSymmA, symmC.template selfadjointView<Lower>()*eiSymm.operatorInverseSqrt()); 130 131 MatrixType id = MatrixType::Identity(rows, cols); 132 VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1)); 133 134 SelfAdjointEigenSolver<MatrixType> eiSymmUninitialized; 135 VERIFY_RAISES_ASSERT(eiSymmUninitialized.info()); 136 VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvalues()); 137 VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors()); 138 VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt()); 139 VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt()); 140 141 eiSymmUninitialized.compute(symmA, false); 142 VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors()); 143 VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt()); 144 VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt()); 145 146 // test Tridiagonalization's methods 147 Tridiagonalization<MatrixType> tridiag(symmC); 148 VERIFY_IS_APPROX(tridiag.diagonal(), tridiag.matrixT().diagonal()); 149 VERIFY_IS_APPROX(tridiag.subDiagonal(), tridiag.matrixT().template diagonal<-1>()); 150 Matrix<RealScalar,Dynamic,Dynamic> T = tridiag.matrixT(); 151 if(rows>1 && cols>1) { 152 // FIXME check that upper and lower part are 0: 153 //VERIFY(T.topRightCorner(rows-2, cols-2).template triangularView<Upper>().isZero()); 154 } 155 VERIFY_IS_APPROX(tridiag.diagonal(), T.diagonal()); 156 VERIFY_IS_APPROX(tridiag.subDiagonal(), T.template diagonal<1>()); 157 VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint()); 158 VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint()); 159 160 // Test computation of eigenvalues from tridiagonal matrix 161 if(rows > 1) 162 { 163 SelfAdjointEigenSolver<MatrixType> eiSymmTridiag; 164 eiSymmTridiag.computeFromTridiagonal(tridiag.matrixT().diagonal(), tridiag.matrixT().diagonal(-1), ComputeEigenvectors); 165 VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmTridiag.eigenvalues()); 166 VERIFY_IS_APPROX(tridiag.matrixT(), eiSymmTridiag.eigenvectors().real() * eiSymmTridiag.eigenvalues().asDiagonal() * eiSymmTridiag.eigenvectors().real().transpose()); 167 } 168 169 if (rows > 1 && rows < 20) 170 { 171 // Test matrix with NaN 172 symmC(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN(); 173 SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmC); 174 VERIFY_IS_EQUAL(eiSymmNaN.info(), NoConvergence); 175 } 176 177 // regression test for bug 1098 178 { 179 SelfAdjointEigenSolver<MatrixType> eig(a.adjoint() * a); 180 eig.compute(a.adjoint() * a); 181 } 182 183 // regression test for bug 478 184 { 185 a.setZero(); 186 SelfAdjointEigenSolver<MatrixType> ei3(a); 187 VERIFY_IS_EQUAL(ei3.info(), Success); 188 VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1)); 189 VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity()); 190 } 191 } 192 193 template<int> 194 void bug_854() 195 { 196 Matrix3d m; 197 m << 850.961, 51.966, 0, 198 51.966, 254.841, 0, 199 0, 0, 0; 200 selfadjointeigensolver_essential_check(m); 201 } 202 203 template<int> 204 void bug_1014() 205 { 206 Matrix3d m; 207 m << 0.11111111111111114658, 0, 0, 208 0, 0.11111111111111109107, 0, 209 0, 0, 0.11111111111111107719; 210 selfadjointeigensolver_essential_check(m); 211 } 212 213 template<int> 214 void bug_1225() 215 { 216 Matrix3d m1, m2; 217 m1.setRandom(); 218 m1 = m1*m1.transpose(); 219 m2 = m1.triangularView<Upper>(); 220 SelfAdjointEigenSolver<Matrix3d> eig1(m1); 221 SelfAdjointEigenSolver<Matrix3d> eig2(m2.selfadjointView<Upper>()); 222 VERIFY_IS_APPROX(eig1.eigenvalues(), eig2.eigenvalues()); 223 } 224 225 template<int> 226 void bug_1204() 227 { 228 SparseMatrix<double> A(2,2); 229 A.setIdentity(); 230 SelfAdjointEigenSolver<Eigen::SparseMatrix<double> > eig(A); 231 } 232 233 EIGEN_DECLARE_TEST(eigensolver_selfadjoint) 234 { 235 int s = 0; 236 for(int i = 0; i < g_repeat; i++) { 237 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 CALL_SUBTEST_1( selfadjointeigensolver(Matrix<std::complex<double>, 1, 1>())); 242 243 // very important to test 3x3 and 2x2 matrices since we provide special paths for them 244 CALL_SUBTEST_12( selfadjointeigensolver(Matrix2f()) ); 245 CALL_SUBTEST_12( selfadjointeigensolver(Matrix2d()) ); 246 CALL_SUBTEST_12( selfadjointeigensolver(Matrix2cd()) ); 247 CALL_SUBTEST_13( selfadjointeigensolver(Matrix3f()) ); 248 CALL_SUBTEST_13( selfadjointeigensolver(Matrix3d()) ); 249 CALL_SUBTEST_13( selfadjointeigensolver(Matrix3cd()) ); 250 CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) ); 251 CALL_SUBTEST_2( selfadjointeigensolver(Matrix4cd()) ); 252 253 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); 254 CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(s,s)) ); 255 CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(s,s)) ); 256 CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(s,s)) ); 257 CALL_SUBTEST_9( selfadjointeigensolver(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(s,s)) ); 258 TEST_SET_BUT_UNUSED_VARIABLE(s) 259 260 // some trivial but implementation-wise tricky cases 261 CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(1,1)) ); 262 CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(2,2)) ); 263 CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(1,1)) ); 264 CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(2,2)) ); 265 CALL_SUBTEST_6( selfadjointeigensolver(Matrix<double,1,1>()) ); 266 CALL_SUBTEST_7( selfadjointeigensolver(Matrix<double,2,2>()) ); 267 } 268 269 CALL_SUBTEST_13( bug_854<0>() ); 270 CALL_SUBTEST_13( bug_1014<0>() ); 271 CALL_SUBTEST_13( bug_1204<0>() ); 272 CALL_SUBTEST_13( bug_1225<0>() ); 273 274 // Test problem size constructors 275 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); 276 CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf> tmp1(s)); 277 CALL_SUBTEST_8(Tridiagonalization<MatrixXf> tmp2(s)); 278 279 TEST_SET_BUT_UNUSED_VARIABLE(s) 280 } 281