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 
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 void test_eigensolver_selfadjoint()
234 {
235  int s = 0;
236  for(int i = 0; i < g_repeat; i++) {
237  // trivial test for 1x1 matrices:
238  CALL_SUBTEST_1( selfadjointeigensolver(Matrix<float, 1, 1>()));
239  CALL_SUBTEST_1( selfadjointeigensolver(Matrix<double, 1, 1>()));
240  // very important to test 3x3 and 2x2 matrices since we provide special paths for them
241  CALL_SUBTEST_12( selfadjointeigensolver(Matrix2f()) );
242  CALL_SUBTEST_12( selfadjointeigensolver(Matrix2d()) );
243  CALL_SUBTEST_13( selfadjointeigensolver(Matrix3f()) );
244  CALL_SUBTEST_13( selfadjointeigensolver(Matrix3d()) );
245  CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) );
246 
247  s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
248  CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(s,s)) );
249  CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(s,s)) );
250  CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(s,s)) );
251  CALL_SUBTEST_9( selfadjointeigensolver(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(s,s)) );
252  TEST_SET_BUT_UNUSED_VARIABLE(s)
253 
254  // some trivial but implementation-wise tricky cases
255  CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(1,1)) );
256  CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(2,2)) );
257  CALL_SUBTEST_6( selfadjointeigensolver(Matrix<double,1,1>()) );
258  CALL_SUBTEST_7( selfadjointeigensolver(Matrix<double,2,2>()) );
259  }
260 
261  CALL_SUBTEST_13( bug_854<0>() );
262  CALL_SUBTEST_13( bug_1014<0>() );
263  CALL_SUBTEST_13( bug_1204<0>() );
264  CALL_SUBTEST_13( bug_1225<0>() );
265 
266  // Test problem size constructors
267  s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
268  CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf> tmp1(s));
269  CALL_SUBTEST_8(Tridiagonalization<MatrixXf> tmp2(s));
270 
271  TEST_SET_BUT_UNUSED_VARIABLE(s)
272 }
273 
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
selfadjointeigensolver
void selfadjointeigensolver(const MatrixType &m)
Definition: eigensolver_selfadjoint.cpp:69
Eigen::SelfAdjointEigenSolver::computeDirect
EIGEN_DEVICE_FUNC SelfAdjointEigenSolver & computeDirect(const MatrixType &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix using a closed-form algorithm.
Definition: SelfAdjointEigenSolver.h:799
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:33
svd_fill_random
void svd_fill_random(MatrixType &m, int Option=0)
Definition: svd_fill.h:21
VERIFY_RAISES_ASSERT
#define VERIFY_RAISES_ASSERT(a)
Definition: main.h:285
Eigen::SelfAdjointEigenSolver::operatorSqrt
EIGEN_DEVICE_FUNC MatrixType operatorSqrt() const
Computes the positive-definite square root of the matrix.
Definition: SelfAdjointEigenSolver.h:306
b
Scalar * b
Definition: benchVecAdd.cpp:17
adjoint
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
test_eigensolver_selfadjoint
void test_eigensolver_selfadjoint()
Definition: eigensolver_selfadjoint.cpp:233
min
#define min(a, b)
Definition: datatypes.h:19
m2
MatrixType m2(n_dims)
Eigen::SelfAdjointEigenSolver::compute
EIGEN_DEVICE_FUNC SelfAdjointEigenSolver & compute(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix.
diagonal
void diagonal(const MatrixType &m)
Definition: diagonal.cpp:12
Eigen::SelfAdjointEigenSolver
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Definition: SelfAdjointEigenSolver.h:70
Eigen::ABx_lx
Definition: Constants.h:403
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
Eigen::Tridiagonalization
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:63
Eigen::BAx_lx
Definition: Constants.h:406
A
Matrix< SCALARA, Dynamic, Dynamic > A
Definition: bench_gemm.cpp:35
bug_1204
void bug_1204()
Definition: eigensolver_selfadjoint.cpp:226
bug_854
void bug_854()
Definition: eigensolver_selfadjoint.cpp:194
svd_fill.h
a
Array33i a
Definition: Cwise_product.cpp:1
Eigen::Ax_lBx
Definition: Constants.h:400
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:331
Eigen::ComputeEigenvectors
Definition: Constants.h:395
Eigen::SelfAdjointEigenSolver::info
EIGEN_DEVICE_FUNC ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SelfAdjointEigenSolver.h:343
m1
Matrix3d m1
Definition: IOFormat.cpp:2
Eigen::Symmetric
Definition: Constants.h:222
Eigen::g_repeat
static int g_repeat
Definition: main.h:144
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Eigen::GeneralizedSelfAdjointEigenSolver
Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem.
Definition: GeneralizedSelfAdjointEigenSolver.h:48
eig
SelfAdjointEigenSolver< PlainMatrixType > eig(mat, computeVectors?ComputeEigenvectors:EigenvaluesOnly)
Eigen::SelfAdjointEigenSolver::computeFromTridiagonal
SelfAdjointEigenSolver & computeFromTridiagonal(const RealVectorType &diag, const SubDiagonalType &subdiag, int options=ComputeEigenvectors)
Computes the eigen decomposition from a tridiagonal symmetric matrix.
Definition: SelfAdjointEigenSolver.h:451
s
RealScalar s
Definition: level1_cplx_impl.h:104
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
Eigen::Upper
Definition: Constants.h:206
Eigen::Success
Definition: Constants.h:432
VERIFY_IS_MUCH_SMALLER_THAN
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition: main.h:335
Eigen::SelfAdjointEigenSolver::eigenvalues
EIGEN_DEVICE_FUNC const RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: SelfAdjointEigenSolver.h:282
simple_graph::b1
Vector2 b1(2,-1)
TEST_SET_BUT_UNUSED_VARIABLE
#define TEST_SET_BUT_UNUSED_VARIABLE(X)
Definition: main.h:91
CALL_SUBTEST
#define CALL_SUBTEST(FUNC)
Definition: main.h:342
VERIFY
#define VERIFY(a)
Definition: main.h:325
bug_1225
void bug_1225()
Definition: eigensolver_selfadjoint.cpp:214
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition: boostmultiprec.cpp:16
Eigen::SelfAdjointEigenSolver::operatorInverseSqrt
EIGEN_DEVICE_FUNC MatrixType operatorInverseSqrt() const
Computes the inverse square root of the matrix.
Definition: SelfAdjointEigenSolver.h:331
VERIFY_IS_UNITARY
#define VERIFY_IS_UNITARY(a)
Definition: main.h:340
Eigen::RowMajor
Definition: Constants.h:322
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
Eigen::SelfAdjointEigenSolver::eigenvectors
EIGEN_DEVICE_FUNC const EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: SelfAdjointEigenSolver.h:259
Eigen::GeneralizedSelfAdjointEigenSolver::compute
GeneralizedSelfAdjointEigenSolver & compute(const MatrixType &matA, const MatrixType &matB, int options=ComputeEigenvectors|Ax_lBx)
Computes generalized eigendecomposition of given matrix pencil.
Definition: GeneralizedSelfAdjointEigenSolver.h:163
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:38
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:21
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
Eigen::internal::isApprox
EIGEN_DEVICE_FUNC bool isApprox(const Scalar &x, const Scalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
Definition: Eigen/src/Core/MathFunctions.h:1361
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
bug_1014
void bug_1014()
Definition: eigensolver_selfadjoint.cpp:204
Eigen::NoConvergence
Definition: Constants.h:436
selfadjointeigensolver_essential_check
void selfadjointeigensolver_essential_check(const MatrixType &m)
Definition: eigensolver_selfadjoint.cpp:18
setZero
v setZero(3)

gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:01