matrix_function.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) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #include "main.h"
11 #include <unsupported/Eigen/MatrixFunctions>
12 
13 // Variant of VERIFY_IS_APPROX which uses absolute error instead of
14 // relative error.
15 #define VERIFY_IS_APPROX_ABS(a, b) VERIFY(test_isApprox_abs(a, b))
16 
17 template<typename Type1, typename Type2>
18 inline bool test_isApprox_abs(const Type1& a, const Type2& b)
19 {
20  return ((a-b).array().abs() < test_precision<typename Type1::RealScalar>()).all();
21 }
22 
23 
24 // Returns a matrix with eigenvalues clustered around 0, 1 and 2.
25 template<typename MatrixType>
27 {
28  typedef typename MatrixType::Scalar Scalar;
29  typedef typename MatrixType::RealScalar RealScalar;
30  MatrixType diag = MatrixType::Zero(size, size);
31  for (Index i = 0; i < size; ++i) {
32  diag(i, i) = Scalar(RealScalar(internal::random<int>(0,2)))
33  + internal::random<Scalar>() * Scalar(RealScalar(0.01));
34  }
35  MatrixType A = MatrixType::Random(size, size);
37  return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
38 }
39 
42 {
43  // Returns a matrix with eigenvalues clustered around 0 and +/- i.
44  static MatrixType run(const Index size);
45 };
46 
47 // Partial specialization for real matrices
48 template<typename MatrixType>
50 {
51  static MatrixType run(const Index size)
52  {
53  typedef typename MatrixType::Scalar Scalar;
54  MatrixType diag = MatrixType::Zero(size, size);
55  Index i = 0;
56  while (i < size) {
57  Index randomInt = internal::random<Index>(-1, 1);
58  if (randomInt == 0 || i == size-1) {
59  diag(i, i) = internal::random<Scalar>() * Scalar(0.01);
60  ++i;
61  } else {
62  Scalar alpha = Scalar(randomInt) + internal::random<Scalar>() * Scalar(0.01);
63  diag(i, i+1) = alpha;
64  diag(i+1, i) = -alpha;
65  i += 2;
66  }
67  }
68  MatrixType A = MatrixType::Random(size, size);
70  return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
71  }
72 };
73 
74 // Partial specialization for complex matrices
75 template<typename MatrixType>
77 {
78  static MatrixType run(const Index size)
79  {
80  typedef typename MatrixType::Scalar Scalar;
81  typedef typename MatrixType::RealScalar RealScalar;
82  const Scalar imagUnit(0, 1);
83  MatrixType diag = MatrixType::Zero(size, size);
84  for (Index i = 0; i < size; ++i) {
85  diag(i, i) = Scalar(RealScalar(internal::random<Index>(-1, 1))) * imagUnit
86  + internal::random<Scalar>() * Scalar(RealScalar(0.01));
87  }
88  MatrixType A = MatrixType::Random(size, size);
90  return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
91  }
92 };
93 
94 
95 template<typename MatrixType>
97 {
99  typedef typename NumTraits<Scalar>::Real RealScalar;
100  typedef std::complex<RealScalar> ComplexScalar;
101 
102  VERIFY_IS_APPROX(A.exp(), A.matrixFunction(internal::stem_function_exp<ComplexScalar>));
103 }
104 
105 template<typename MatrixType>
107 {
109  typedef typename NumTraits<Scalar>::Real RealScalar;
110 
111  MatrixType scaledA;
112  RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff();
113  if (maxImagPartOfSpectrum >= RealScalar(0.9L * EIGEN_PI))
114  scaledA = A * RealScalar(0.9L * EIGEN_PI) / maxImagPartOfSpectrum;
115  else
116  scaledA = A;
117 
118  // identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X
119  MatrixType expA = scaledA.exp();
120  MatrixType logExpA = expA.log();
121  VERIFY_IS_APPROX(logExpA, scaledA);
122 }
123 
124 template<typename MatrixType>
126 {
127  // Need to use absolute error because of possible cancellation when
128  // adding/subtracting expA and expmA.
129  VERIFY_IS_APPROX_ABS(A.sinh(), (A.exp() - (-A).exp()) / 2);
130  VERIFY_IS_APPROX_ABS(A.cosh(), (A.exp() + (-A).exp()) / 2);
131 }
132 
133 template<typename MatrixType>
135 {
136  typedef typename MatrixType::Scalar Scalar;
137  typedef typename NumTraits<Scalar>::Real RealScalar;
138  typedef std::complex<RealScalar> ComplexScalar;
139  typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime,
140  MatrixType::ColsAtCompileTime, MatrixType::Options> ComplexMatrix;
141 
142  ComplexScalar imagUnit(0,1);
143  ComplexScalar two(2,0);
144 
145  ComplexMatrix Ac = A.template cast<ComplexScalar>();
146 
147  ComplexMatrix exp_iA = (imagUnit * Ac).exp();
148  ComplexMatrix exp_miA = (-imagUnit * Ac).exp();
149 
150  ComplexMatrix sinAc = A.sin().template cast<ComplexScalar>();
151  VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit));
152 
153  ComplexMatrix cosAc = A.cos().template cast<ComplexScalar>();
154  VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2);
155 }
156 
157 template<typename MatrixType>
158 void testMatrix(const MatrixType& A)
159 {
164 }
165 
166 template<typename MatrixType>
168 {
169  // Matrices with clustered eigenvalue lead to different code paths
170  // in MatrixFunction.h and are thus useful for testing.
171 
172  const Index size = m.rows();
173  for (int i = 0; i < g_repeat; i++) {
174  testMatrix(MatrixType::Random(size, size).eval());
175  testMatrix(randomMatrixWithRealEivals<MatrixType>(size));
177  }
178 }
179 
180 template<typename MatrixType>
181 void testMapRef(const MatrixType& A)
182 {
183  // Test if passing Ref and Map objects is possible
184  // (Regression test for Bug #1796)
185  Index size = A.rows();
186  MatrixType X; X.setRandom(size, size);
187  MatrixType Y(size,size);
188  Ref< MatrixType> R(Y);
189  Ref<const MatrixType> Rc(X);
190  Map< MatrixType> M(Y.data(), size, size);
191  Map<const MatrixType> Mc(X.data(), size, size);
192 
193  X = X*X; // make sure sqrt is possible
194  Y = X.sqrt();
195  R = Rc.sqrt();
196  M = Mc.sqrt();
197  Y = X.exp();
198  R = Rc.exp();
199  M = Mc.exp();
200  X = Y; // make sure log is possible
201  Y = X.log();
202  R = Rc.log();
203  M = Mc.log();
204 
205  Y = X.cos() + Rc.cos() + Mc.cos();
206  Y = X.sin() + Rc.sin() + Mc.sin();
207 
208  Y = X.cosh() + Rc.cosh() + Mc.cosh();
209  Y = X.sinh() + Rc.sinh() + Mc.sinh();
210 }
211 
212 
213 EIGEN_DECLARE_TEST(matrix_function)
214 {
216  CALL_SUBTEST_2(testMatrixType(Matrix3cf()));
217  CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8)));
218  CALL_SUBTEST_4(testMatrixType(Matrix2d()));
220  CALL_SUBTEST_6(testMatrixType(Matrix4cd()));
221  CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13)));
222 
224  CALL_SUBTEST_2(testMapRef(Matrix3cf()));
225  CALL_SUBTEST_3(testMapRef(MatrixXf(8,8)));
226  CALL_SUBTEST_7(testMapRef(MatrixXd(13,13)));
227 }
Matrix3f m
int array[24]
void testMatrix(const MatrixType &A)
void testMatrixLogarithm(const MatrixType &A)
SCALAR Scalar
Definition: bench_gemm.cpp:46
const char Y
#define CALL_SUBTEST_6(FUNC)
#define CALL_SUBTEST_4(FUNC)
Matrix< RealScalar, Dynamic, Dynamic > M
Definition: bench_gemm.cpp:51
void testMatrixExponential(const MatrixType &A)
Scalar * b
Definition: benchVecAdd.cpp:17
bool test_isApprox_abs(const Type1 &a, const Type2 &b)
#define EIGEN_PI
Matrix diag(const std::vector< Matrix > &Hs)
Definition: Matrix.cpp:206
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
#define CALL_SUBTEST_3(FUNC)
#define CALL_SUBTEST_7(FUNC)
void testMapRef(const MatrixType &A)
Rot2 R(Rot2::fromAngle(0.1))
void testMatrixType(const MatrixType &m)
MatrixXd L
Definition: LLT_example.cpp:6
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
MatrixType randomMatrixWithRealEivals(const Index size)
EIGEN_DECLARE_TEST(matrix_function)
void testHyperbolicFunctions(const MatrixType &A)
EIGEN_DEVICE_FUNC const ExpReturnType exp() const
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
AdjointReturnType inverse() const
Inverse of the Householder sequence (equals the adjoint).
#define VERIFY_IS_APPROX(a, b)
#define CALL_SUBTEST_1(FUNC)
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
RealScalar alpha
HouseholderSequenceType householderQ() const
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
static MatrixType run(const Index size)
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
#define CALL_SUBTEST_5(FUNC)
static MatrixType run(const Index size)
static MatrixType run(const Index size)
#define VERIFY_IS_APPROX_ABS(a, b)
#define CALL_SUBTEST_2(FUNC)
internal::nested_eval< T, 1 >::type eval(const T &xpr)
The matrix class, also used for vectors and row-vectors.
#define X
Definition: icosphere.cpp:20
#define abs(x)
Definition: datatypes.h:17
void testGonioFunctions(const MatrixType &A)


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