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());
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>());
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);
256  CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(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 
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:828
VERIFY_IS_MUCH_SMALLER_THAN
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition: main.h:390
Eigen::Symmetric
@ Symmetric
Definition: Constants.h:227
Eigen::BAx_lx
@ BAx_lx
Definition: Constants.h:416
Eigen::SparseMatrix< double >
simple_graph::b1
Vector2 b1(2, -1)
Eigen::Tridiagonalization
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:64
s
RealScalar s
Definition: level1_cplx_impl.h:126
Eigen::Tridiagonalization::subDiagonal
SubDiagonalReturnType subDiagonal() const
Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:316
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
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:472
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
b
Scalar * b
Definition: benchVecAdd.cpp:17
selfadjointeigensolver
void selfadjointeigensolver(const MatrixType &m)
Definition: eigensolver_selfadjoint.cpp:69
bug_1014
void bug_1014()
Definition: eigensolver_selfadjoint.cpp:204
Eigen::Upper
@ Upper
Definition: Constants.h:211
m1
Matrix3d m1
Definition: IOFormat.cpp:2
Eigen::Success
@ Success
Definition: Constants.h:442
Eigen::SelfAdjointEigenSolver::compute
EIGEN_DEVICE_FUNC SelfAdjointEigenSolver & compute(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix.
selfadjointeigensolver_essential_check
void selfadjointeigensolver_essential_check(const MatrixType &m)
Definition: eigensolver_selfadjoint.cpp:18
Eigen::SelfAdjointEigenSolver::eigenvalues
const EIGEN_DEVICE_FUNC RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: SelfAdjointEigenSolver.h:300
Eigen::RowMajor
@ RowMajor
Definition: Constants.h:321
Eigen::Tridiagonalization::matrixQ
HouseholderSequenceType matrixQ() const
Returns the unitary matrix Q in the decomposition.
Definition: Tridiagonalization.h:241
Eigen::ComputeEigenvectors
@ ComputeEigenvectors
Definition: Constants.h:405
CALL_SUBTEST_9
#define CALL_SUBTEST_9(FUNC)
Definition: split_test_helper.h:52
Eigen::Tridiagonalization::matrixT
MatrixTReturnType matrixT() const
Returns an expression of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:266
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
VERIFY_RAISES_ASSERT
#define VERIFY_RAISES_ASSERT(a)
Definition: main.h:340
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
Eigen::SelfAdjointEigenSolver::info
EIGEN_DEVICE_FUNC ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SelfAdjointEigenSolver.h:361
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(eigensolver_selfadjoint)
Definition: eigensolver_selfadjoint.cpp:233
align_3::a1
Point2 a1
Definition: testPose2.cpp:769
m2
MatrixType m2(n_dims)
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
Eigen::SelfAdjointEigenSolver::eigenvectors
const EIGEN_DEVICE_FUNC EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: SelfAdjointEigenSolver.h:277
Eigen::NoConvergence
@ NoConvergence
Definition: Constants.h:446
svd_fill_random
void svd_fill_random(MatrixType &m, int Option=0)
Definition: svd_fill.h:21
CALL_SUBTEST_13
#define CALL_SUBTEST_13(FUNC)
Definition: split_test_helper.h:76
Eigen::SelfAdjointEigenSolver
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Definition: SelfAdjointEigenSolver.h:76
bug_1225
void bug_1225()
Definition: eigensolver_selfadjoint.cpp:214
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:22
Eigen::Ax_lBx
@ Ax_lBx
Definition: Constants.h:410
bug_1204
void bug_1204()
Definition: eigensolver_selfadjoint.cpp:226
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
Eigen::Triplet< double >
Eigen::SelfAdjointEigenSolver::operatorInverseSqrt
EIGEN_DEVICE_FUNC MatrixType operatorInverseSqrt() const
Computes the inverse square root of the matrix.
Definition: SelfAdjointEigenSolver.h:349
eig
SelfAdjointEigenSolver< PlainMatrixType > eig(mat, computeVectors?ComputeEigenvectors:EigenvaluesOnly)
CALL_SUBTEST_6
#define CALL_SUBTEST_6(FUNC)
Definition: split_test_helper.h:34
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
Eigen::HouseholderSequence::adjoint
AdjointReturnType adjoint() const
Adjoint (conjugate transpose) of the Householder sequence.
Definition: HouseholderSequence.h:266
VERIFY_IS_UNITARY
#define VERIFY_IS_UNITARY(a)
Definition: main.h:395
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
Eigen::SelfAdjointEigenSolver::operatorSqrt
EIGEN_DEVICE_FUNC MatrixType operatorSqrt() const
Computes the positive-definite square root of the matrix.
Definition: SelfAdjointEigenSolver.h:324
Eigen::Tridiagonalization::diagonal
DiagonalReturnType diagonal() const
Returns the diagonal of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:308
main.h
Eigen::GeneralizedSelfAdjointEigenSolver
Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem.
Definition: GeneralizedSelfAdjointEigenSolver.h:48
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition: boostmultiprec.cpp:16
CALL_SUBTEST_12
#define CALL_SUBTEST_12(FUNC)
Definition: split_test_helper.h:70
svd_fill.h
min
#define min(a, b)
Definition: datatypes.h:19
Eigen::ABx_lx
@ ABx_lx
Definition: Constants.h:413
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
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
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
CALL_SUBTEST_7
#define CALL_SUBTEST_7(FUNC)
Definition: split_test_helper.h:40
CALL_SUBTEST_8
#define CALL_SUBTEST_8(FUNC)
Definition: split_test_helper.h:46
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
TEST_SET_BUT_UNUSED_VARIABLE
#define TEST_SET_BUT_UNUSED_VARIABLE(X)
Definition: main.h:121
bug_854
void bug_854()
Definition: eigensolver_selfadjoint.cpp:194
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
CALL_SUBTEST
#define CALL_SUBTEST(FUNC)
Definition: main.h:399
VERIFY
#define VERIFY(a)
Definition: main.h:380
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:29