eigensolver_generic.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,2012 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 <limits>
13 #include <Eigen/Eigenvalues>
14 
15 template<typename EigType,typename MatType>
16 void check_eigensolver_for_given_mat(const EigType &eig, const MatType& a)
17 {
18  typedef typename NumTraits<typename MatType::Scalar>::Real RealScalar;
19  typedef Matrix<RealScalar, MatType::RowsAtCompileTime, 1> RealVectorType;
20  typedef typename std::complex<RealScalar> Complex;
21  Index n = a.rows();
22  VERIFY_IS_EQUAL(eig.info(), Success);
23  VERIFY_IS_APPROX(a * eig.pseudoEigenvectors(), eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix());
24  VERIFY_IS_APPROX(a.template cast<Complex>() * eig.eigenvectors(),
25  eig.eigenvectors() * eig.eigenvalues().asDiagonal());
26  VERIFY_IS_APPROX(eig.eigenvectors().colwise().norm(), RealVectorType::Ones(n).transpose());
27  VERIFY_IS_APPROX(a.eigenvalues(), eig.eigenvalues());
28 }
29 
30 template<typename MatrixType> void eigensolver(const MatrixType& m)
31 {
32  /* this test covers the following files:
33  EigenSolver.h
34  */
35  Index rows = m.rows();
36  Index cols = m.cols();
37 
38  typedef typename MatrixType::Scalar Scalar;
39  typedef typename NumTraits<Scalar>::Real RealScalar;
40  typedef typename std::complex<RealScalar> Complex;
41 
42  MatrixType a = MatrixType::Random(rows,cols);
43  MatrixType a1 = MatrixType::Random(rows,cols);
44  MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
45 
46  EigenSolver<MatrixType> ei0(symmA);
47  VERIFY_IS_EQUAL(ei0.info(), Success);
48  VERIFY_IS_APPROX(symmA * ei0.pseudoEigenvectors(), ei0.pseudoEigenvectors() * ei0.pseudoEigenvalueMatrix());
49  VERIFY_IS_APPROX((symmA.template cast<Complex>()) * (ei0.pseudoEigenvectors().template cast<Complex>()),
50  (ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal()));
51 
52  EigenSolver<MatrixType> ei1(a);
53  CALL_SUBTEST( check_eigensolver_for_given_mat(ei1,a) );
54 
55  EigenSolver<MatrixType> ei2;
56  ei2.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a);
57  VERIFY_IS_EQUAL(ei2.info(), Success);
58  VERIFY_IS_EQUAL(ei2.eigenvectors(), ei1.eigenvectors());
59  VERIFY_IS_EQUAL(ei2.eigenvalues(), ei1.eigenvalues());
60  if (rows > 2) {
61  ei2.setMaxIterations(1).compute(a);
62  VERIFY_IS_EQUAL(ei2.info(), NoConvergence);
63  VERIFY_IS_EQUAL(ei2.getMaxIterations(), 1);
64  }
65 
66  EigenSolver<MatrixType> eiNoEivecs(a, false);
67  VERIFY_IS_EQUAL(eiNoEivecs.info(), Success);
68  VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
69  VERIFY_IS_APPROX(ei1.pseudoEigenvalueMatrix(), eiNoEivecs.pseudoEigenvalueMatrix());
70 
71  MatrixType id = MatrixType::Identity(rows, cols);
72  VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
73 
74  if (rows > 2 && rows < 20)
75  {
76  // Test matrix with NaN
77  a(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
78  EigenSolver<MatrixType> eiNaN(a);
79  VERIFY_IS_NOT_EQUAL(eiNaN.info(), Success);
80  }
81 
82  // regression test for bug 1098
83  {
84  EigenSolver<MatrixType> eig(a.adjoint() * a);
85  eig.compute(a.adjoint() * a);
86  }
87 
88  // regression test for bug 478
89  {
90  a.setZero();
91  EigenSolver<MatrixType> ei3(a);
92  VERIFY_IS_EQUAL(ei3.info(), Success);
93  VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
94  VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
95  }
96 }
97 
98 template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
99 {
100  EigenSolver<MatrixType> eig;
101  VERIFY_RAISES_ASSERT(eig.eigenvectors());
102  VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
103  VERIFY_RAISES_ASSERT(eig.pseudoEigenvalueMatrix());
104  VERIFY_RAISES_ASSERT(eig.eigenvalues());
105 
106  MatrixType a = MatrixType::Random(m.rows(),m.cols());
107  eig.compute(a, false);
108  VERIFY_RAISES_ASSERT(eig.eigenvectors());
109  VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
110 }
111 
112 
113 template<typename CoeffType>
114 Matrix<typename CoeffType::Scalar,Dynamic,Dynamic>
115 make_companion(const CoeffType& coeffs)
116 {
117  Index n = coeffs.size()-1;
118  Matrix<typename CoeffType::Scalar,Dynamic,Dynamic> res(n,n);
119  res.setZero();
120  res.row(0) = -coeffs.tail(n) / coeffs(0);
121  res.diagonal(-1).setOnes();
122  return res;
123 }
124 
125 template<int>
126 void eigensolver_generic_extra()
127 {
128  {
129  // regression test for bug 793
130  MatrixXd a(3,3);
131  a << 0, 0, 1,
132  1, 1, 1,
133  1, 1e+200, 1;
134  Eigen::EigenSolver<MatrixXd> eig(a);
135  double scale = 1e-200; // scale to avoid overflow during the comparisons
136  VERIFY_IS_APPROX(a * eig.pseudoEigenvectors()*scale, eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()*scale);
137  VERIFY_IS_APPROX(a * eig.eigenvectors()*scale, eig.eigenvectors() * eig.eigenvalues().asDiagonal()*scale);
138  }
139  {
140  // check a case where all eigenvalues are null.
141  MatrixXd a(2,2);
142  a << 1, 1,
143  -1, -1;
144  Eigen::EigenSolver<MatrixXd> eig(a);
145  VERIFY_IS_APPROX(eig.pseudoEigenvectors().squaredNorm(), 2.);
146  VERIFY_IS_APPROX((a * eig.pseudoEigenvectors()).norm()+1., 1.);
147  VERIFY_IS_APPROX((eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()).norm()+1., 1.);
148  VERIFY_IS_APPROX((a * eig.eigenvectors()).norm()+1., 1.);
149  VERIFY_IS_APPROX((eig.eigenvectors() * eig.eigenvalues().asDiagonal()).norm()+1., 1.);
150  }
151 
152  // regression test for bug 933
153  {
154  {
155  VectorXd coeffs(5); coeffs << 1, -3, -175, -225, 2250;
156  MatrixXd C = make_companion(coeffs);
157  EigenSolver<MatrixXd> eig(C);
158  CALL_SUBTEST( check_eigensolver_for_given_mat(eig,C) );
159  }
160  {
161  // this test is tricky because it requires high accuracy in smallest eigenvalues
162  VectorXd coeffs(5); coeffs << 6.154671e-15, -1.003870e-10, -9.819570e-01, 3.995715e+03, 2.211511e+08;
163  MatrixXd C = make_companion(coeffs);
164  EigenSolver<MatrixXd> eig(C);
165  CALL_SUBTEST( check_eigensolver_for_given_mat(eig,C) );
166  Index n = C.rows();
167  for(Index i=0;i<n;++i)
168  {
169  typedef std::complex<double> Complex;
170  MatrixXcd ac = C.cast<Complex>();
171  ac.diagonal().array() -= eig.eigenvalues()(i);
172  VectorXd sv = ac.jacobiSvd().singularValues();
173  // comparing to sv(0) is not enough here to catch the "bug",
174  // the hard-coded 1.0 is important!
175  VERIFY_IS_MUCH_SMALLER_THAN(sv(n-1), 1.0);
176  }
177  }
178  }
179  // regression test for bug 1557
180  {
181  // this test is interesting because it contains zeros on the diagonal.
182  MatrixXd A_bug1557(3,3);
183  A_bug1557 << 0, 0, 0, 1, 0, 0.5887907064808635127, 0, 1, 0;
184  EigenSolver<MatrixXd> eig(A_bug1557);
185  CALL_SUBTEST( check_eigensolver_for_given_mat(eig,A_bug1557) );
186  }
187 
188  // regression test for bug 1174
189  {
190  Index n = 12;
191  MatrixXf A_bug1174(n,n);
192  A_bug1174 << 262144, 0, 0, 262144, 786432, 0, 0, 0, 0, 0, 0, 786432,
193  262144, 0, 0, 262144, 786432, 0, 0, 0, 0, 0, 0, 786432,
194  262144, 0, 0, 262144, 786432, 0, 0, 0, 0, 0, 0, 786432,
195  262144, 0, 0, 262144, 786432, 0, 0, 0, 0, 0, 0, 786432,
196  0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
197  0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
198  0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
199  0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
200  0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
201  0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
202  0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0,
203  0, 262144, 262144, 0, 0, 262144, 262144, 262144, 262144, 262144, 262144, 0;
204  EigenSolver<MatrixXf> eig(A_bug1174);
205  CALL_SUBTEST( check_eigensolver_for_given_mat(eig,A_bug1174) );
206  }
207 }
208 
209 EIGEN_DECLARE_TEST(eigensolver_generic)
210 {
211  int s = 0;
212  for(int i = 0; i < g_repeat; i++) {
213  CALL_SUBTEST_1( eigensolver(Matrix4f()) );
214  s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
215  CALL_SUBTEST_2( eigensolver(MatrixXd(s,s)) );
216  TEST_SET_BUT_UNUSED_VARIABLE(s)
217 
218  // some trivial but implementation-wise tricky cases
219  CALL_SUBTEST_2( eigensolver(MatrixXd(1,1)) );
220  CALL_SUBTEST_2( eigensolver(MatrixXd(2,2)) );
221  CALL_SUBTEST_3( eigensolver(Matrix<double,1,1>()) );
222  CALL_SUBTEST_4( eigensolver(Matrix2d()) );
223  }
224 
225  CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4f()) );
226  s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
227  CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXd(s,s)) );
228  CALL_SUBTEST_3( eigensolver_verify_assert(Matrix<double,1,1>()) );
229  CALL_SUBTEST_4( eigensolver_verify_assert(Matrix2d()) );
230 
231  // Test problem size constructors
232  CALL_SUBTEST_5(EigenSolver<MatrixXf> tmp(s));
233 
234  // regression test for bug 410
235  CALL_SUBTEST_2(
236  {
237  MatrixXd A(1,1);
238  A(0,0) = std::sqrt(-1.); // is Not-a-Number
239  Eigen::EigenSolver<MatrixXd> solver(A);
240  VERIFY_IS_EQUAL(solver.info(), NumericalIssue);
241  }
242  );
243 
244  CALL_SUBTEST_2( eigensolver_generic_extra<0>() );
245 
246  TEST_SET_BUT_UNUSED_VARIABLE(s)
247 }
align_3::a1
Point2 a1
Definition: testPose2.cpp:769
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
Eigen::EigenSolver::pseudoEigenvalueMatrix
MatrixType pseudoEigenvalueMatrix() const
Returns the block-diagonal matrix in the pseudo-eigendecomposition.
Definition: EigenSolver.h:324
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
VERIFY_RAISES_ASSERT
#define VERIFY_RAISES_ASSERT(a)
Definition: main.h:340
Eigen::EigenSolver::pseudoEigenvectors
const MatrixType & pseudoEigenvectors() const
Returns the pseudo-eigenvectors of given matrix.
Definition: EigenSolver.h:199
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::setZero
EIGEN_DEVICE_FUNC Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > & setZero(Index size)
Definition: CwiseNullaryOp.h:562
Eigen::EigenSolver::compute
EigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
VERIFY_IS_NOT_EQUAL
#define VERIFY_IS_NOT_EQUAL(a, b)
Definition: main.h:387
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
eigensolver_verify_assert
void eigensolver_verify_assert(const MatrixType &m)
Definition: eigensolver_generic.cpp:98
Eigen::RealSchur
Performs a real Schur decomposition of a square matrix.
Definition: RealSchur.h:54
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
n
int n
Definition: BiCGSTAB_simple.cpp:1
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
solver
BiCGSTAB< SparseMatrix< double > > solver
Definition: BiCGSTAB_simple.cpp:5
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
Eigen::EigenSolver::getMaxIterations
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: EigenSolver.h:295
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::setOnes
EIGEN_DEVICE_FUNC Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > & setOnes(Index size)
Definition: CwiseNullaryOp.h:714
Complex
std::complex< RealScalar > Complex
Definition: gtsam/3rdparty/Eigen/blas/common.h:94
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
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:386
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
Eigen::EigenSolver::info
ComputationInfo info() const
Definition: EigenSolver.h:281
Eigen::NumericalIssue
Definition: Constants.h:444
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(eigensolver_generic)
Definition: eigensolver_generic.cpp:209
eig
SelfAdjointEigenSolver< PlainMatrixType > eig(mat, computeVectors?ComputeEigenvectors:EigenvaluesOnly)
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
s
RealScalar s
Definition: level1_cplx_impl.h:126
Eigen::EigenSolver::setMaxIterations
EigenSolver & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: EigenSolver.h:288
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Eigen::Success
Definition: Constants.h:442
VERIFY_IS_MUCH_SMALLER_THAN
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition: main.h:390
C
Matrix< Scalar, Dynamic, Dynamic > C
Definition: bench_gemm.cpp:50
Eigen::EigenSolver::eigenvectors
EigenvectorsType eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: EigenSolver.h:345
TEST_SET_BUT_UNUSED_VARIABLE
#define TEST_SET_BUT_UNUSED_VARIABLE(X)
Definition: main.h:121
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
CALL_SUBTEST
#define CALL_SUBTEST(FUNC)
Definition: main.h:399
make_companion
Matrix< typename CoeffType::Scalar, Dynamic, Dynamic > make_companion(const CoeffType &coeffs)
Definition: eigensolver_generic.cpp:115
VERIFY
#define VERIFY(a)
Definition: main.h:380
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition: boostmultiprec.cpp:16
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
scale
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics scale
Definition: gnuplot_common_settings.hh:12
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:38
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
Eigen::EigenSolver::eigenvalues
const EigenvalueType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: EigenSolver.h:244
Eigen::EigenSolver
Computes eigenvalues and eigenvectors of general matrices.
Definition: EigenSolver.h:64
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
check_eigensolver_for_given_mat
void check_eigensolver_for_given_mat(const EigType &eig, const MatType &a)
Definition: eigensolver_generic.cpp:16
eigensolver
void eigensolver(const MatrixType &m)
Definition: eigensolver_generic.cpp:30
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::NoConvergence
Definition: Constants.h:446
eigensolver_generic_extra
void eigensolver_generic_extra()
Definition: eigensolver_generic.cpp:126

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