eigensolver_selfadjoint.cpp
Go to the documentation of this file.
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 
25  VERIFY_IS_EQUAL(eiSymm.info(), Success);
26 
27  RealScalar scaling = m.cwiseAbs().maxCoeff();
28 
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  {
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  }
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 
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 
98 
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  {
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  {
180  eig.compute(a.adjoint() * a);
181  }
182 
183  // regression test for bug 478
184  {
185  a.setZero();
187  VERIFY_IS_EQUAL(ei3.info(), Success);
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;
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;
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>();
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();
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:
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
246  CALL_SUBTEST_12( selfadjointeigensolver(Matrix2cd()) );
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)) );
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)) );
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);
278 
280 }
281 
Point2 a1
Definition: testPose2.cpp:769
Matrix3f m
void selfadjointeigensolver(const MatrixType &m)
EIGEN_DEVICE_FUNC SelfAdjointEigenSolver & computeDirect(const MatrixType &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix using a closed-form algorithm.
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define CALL_SUBTEST_12(FUNC)
void svd_fill_random(MatrixType &m, int Option=0)
Definition: svd_fill.h:21
#define VERIFY_RAISES_ASSERT(a)
Definition: main.h:340
#define CALL_SUBTEST_9(FUNC)
#define CALL_SUBTEST_6(FUNC)
#define CALL_SUBTEST_4(FUNC)
Scalar * b
Definition: benchVecAdd.cpp:17
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
EIGEN_DEVICE_FUNC ComputationInfo info() const
Reports whether previous computation was successful.
#define min(a, b)
Definition: datatypes.h:19
#define CALL_SUBTEST_13(FUNC)
#define CALL_SUBTEST_3(FUNC)
MatrixType m2(n_dims)
EIGEN_DEVICE_FUNC SelfAdjointEigenSolver & compute(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix.
#define CALL_SUBTEST_7(FUNC)
void diagonal(const MatrixType &m)
Definition: diagonal.cpp:12
Computes eigenvalues and eigenvectors of selfadjoint matrices.
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
Tridiagonal decomposition of a selfadjoint matrix.
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
void bug_1204()
void bug_854()
#define VERIFY_IS_APPROX(a, b)
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
#define CALL_SUBTEST_1(FUNC)
Matrix3d m1
Definition: IOFormat.cpp:2
EIGEN_DEVICE_FUNC const RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
static int g_repeat
Definition: main.h:169
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem.
SelfAdjointEigenSolver< PlainMatrixType > eig(mat, computeVectors?ComputeEigenvectors:EigenvaluesOnly)
SelfAdjointEigenSolver & computeFromTridiagonal(const RealVectorType &diag, const SubDiagonalType &subdiag, int options=ComputeEigenvectors)
Computes the eigen decomposition from a tridiagonal symmetric matrix.
#define CALL_SUBTEST_8(FUNC)
RealScalar s
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition: main.h:390
#define TEST_SET_BUT_UNUSED_VARIABLE(X)
Definition: main.h:121
#define CALL_SUBTEST_5(FUNC)
#define CALL_SUBTEST(FUNC)
Definition: main.h:399
#define VERIFY(a)
Definition: main.h:380
void bug_1225()
#define EIGEN_TEST_MAX_SIZE
#define VERIFY_IS_UNITARY(a)
Definition: main.h:395
EIGEN_DECLARE_TEST(eigensolver_selfadjoint)
EIGEN_DEVICE_FUNC MatrixType operatorInverseSqrt() const
Computes the inverse square root of the matrix.
GeneralizedSelfAdjointEigenSolver & compute(const MatrixType &matA, const MatrixType &matB, int options=ComputeEigenvectors|Ax_lBx)
Computes generalized eigendecomposition of given matrix pencil.
#define CALL_SUBTEST_2(FUNC)
internal::nested_eval< T, 1 >::type eval(const T &xpr)
EIGEN_DEVICE_FUNC const EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
const int Dynamic
Definition: Constants.h:22
The matrix class, also used for vectors and row-vectors.
EIGEN_DEVICE_FUNC bool isApprox(const Scalar &x, const Scalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
EIGEN_DEVICE_FUNC MatrixType operatorSqrt() const
Computes the positive-definite square root of the matrix.
void bug_1014()
void selfadjointeigensolver_essential_check(const MatrixType &m)
v setZero(3)
Vector2 b1(2, -1)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:12