matrix_power.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) 2012, 2013 Chen-Pang He <jdh8@ms63.hinet.net>
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 "matrix_functions.h"
11 
12 template<typename T>
13 void test2dRotation(const T& tol)
14 {
15  Matrix<T,2,2> A, B, C;
16  T angle, c, s;
17 
18  A << 0, 1, -1, 0;
20 
21  for (int i=0; i<=20; ++i) {
22  angle = std::pow(T(10), T(i-10) / T(5.));
23  c = std::cos(angle);
24  s = std::sin(angle);
25  B << c, s, -s, c;
26 
27  C = Apow(std::ldexp(angle,1) / T(EIGEN_PI));
28  std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n';
29  VERIFY(C.isApprox(B, tol));
30  }
31 }
32 
33 template<typename T>
35 {
36  Matrix<std::complex<T>,2,2> A, B, C;
37  T angle, ch = std::cosh((T)1);
38  std::complex<T> ish(0, std::sinh((T)1));
39 
40  A << ch, ish, -ish, ch;
42 
43  for (int i=0; i<=20; ++i) {
44  angle = std::ldexp(static_cast<T>(i-10), -1);
45  ch = std::cosh(angle);
46  ish = std::complex<T>(0, std::sinh(angle));
47  B << ch, ish, -ish, ch;
48 
49  C = Apow(angle);
50  std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n';
51  VERIFY(C.isApprox(B, tol));
52  }
53 }
54 
55 template<typename T>
56 void test3dRotation(const T& tol)
57 {
59  T angle;
60 
61  for (int i=0; i<=20; ++i) {
63  v.normalize();
64  angle = std::pow(T(10), T(i-10) / T(5.));
65  VERIFY(AngleAxis<T>(angle, v).matrix().isApprox(AngleAxis<T>(1,v).matrix().pow(angle), tol));
66  }
67 }
68 
69 template<typename MatrixType>
70 void testGeneral(const MatrixType& m, const typename MatrixType::RealScalar& tol)
71 {
72  typedef typename MatrixType::RealScalar RealScalar;
73  MatrixType m1, m2, m3, m4, m5;
74  RealScalar x, y;
75 
76  for (int i=0; i < g_repeat; ++i) {
79 
80  x = internal::random<RealScalar>();
81  y = internal::random<RealScalar>();
82  m2 = mpow(x);
83  m3 = mpow(y);
84 
85  m4 = mpow(x+y);
86  m5.noalias() = m2 * m3;
87  VERIFY(m4.isApprox(m5, tol));
88 
89  m4 = mpow(x*y);
90  m5 = m2.pow(y);
91  VERIFY(m4.isApprox(m5, tol));
92 
93  m4 = (std::abs(x) * m1).pow(y);
94  m5 = std::pow(std::abs(x), y) * m3;
95  VERIFY(m4.isApprox(m5, tol));
96  }
97 }
98 
99 template<typename MatrixType>
100 void testSingular(const MatrixType& m_const, const typename MatrixType::RealScalar& tol)
101 {
102  // we need to pass by reference in order to prevent errors with
103  // MSVC for aligned data types ...
104  MatrixType& m = const_cast<MatrixType&>(m_const);
105 
107  typedef typename internal::conditional<IsComplex, TriangularView<MatrixType,Upper>, const MatrixType&>::type TriangularType;
108  typename internal::conditional< IsComplex, ComplexSchur<MatrixType>, RealSchur<MatrixType> >::type schur;
109  MatrixType T;
110 
111  for (int i=0; i < g_repeat; ++i) {
112  m.setRandom();
113  m.col(0).fill(0);
114 
115  schur.compute(m);
116  T = schur.matrixT();
117  const MatrixType& U = schur.matrixU();
120 
121  T = T.sqrt();
122  VERIFY(mpow(0.5L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
123 
124  T = T.sqrt();
125  VERIFY(mpow(0.25L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
126 
127  T = T.sqrt();
128  VERIFY(mpow(0.125L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
129  }
130 }
131 
132 template<typename MatrixType>
133 void testLogThenExp(const MatrixType& m_const, const typename MatrixType::RealScalar& tol)
134 {
135  // we need to pass by reference in order to prevent errors with
136  // MSVC for aligned data types ...
137  MatrixType& m = const_cast<MatrixType&>(m_const);
138 
139  typedef typename MatrixType::Scalar Scalar;
140  Scalar x;
141 
142  for (int i=0; i < g_repeat; ++i) {
144  x = internal::random<Scalar>();
145  VERIFY(m.pow(x).isApprox((x * m.log()).exp(), tol));
146  }
147 }
148 
152 
153 EIGEN_DECLARE_TEST(matrix_power)
154 {
155  CALL_SUBTEST_2(test2dRotation<double>(1e-13));
156  CALL_SUBTEST_1(test2dRotation<float>(2e-5f)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64
157  CALL_SUBTEST_9(test2dRotation<long double>(1e-13L));
158  CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14));
159  CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5f));
160  CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14L));
161 
162  CALL_SUBTEST_10(test3dRotation<double>(1e-13));
163  CALL_SUBTEST_11(test3dRotation<float>(1e-5f));
164  CALL_SUBTEST_12(test3dRotation<long double>(1e-13L));
165 
166  CALL_SUBTEST_2(testGeneral(Matrix2d(), 1e-13));
168  CALL_SUBTEST_3(testGeneral(Matrix4cd(), 1e-13));
169  CALL_SUBTEST_4(testGeneral(MatrixXd(8,8), 2e-12));
170  CALL_SUBTEST_1(testGeneral(Matrix2f(), 1e-4f));
171  CALL_SUBTEST_5(testGeneral(Matrix3cf(), 1e-4f));
172  CALL_SUBTEST_8(testGeneral(Matrix4f(), 1e-4f));
173  CALL_SUBTEST_6(testGeneral(MatrixXf(2,2), 1e-3f)); // see bug 614
175  CALL_SUBTEST_10(testGeneral(Matrix3d(), 1e-13));
176  CALL_SUBTEST_11(testGeneral(Matrix3f(), 1e-4f));
178 
179  CALL_SUBTEST_2(testSingular(Matrix2d(), 1e-13));
181  CALL_SUBTEST_3(testSingular(Matrix4cd(), 1e-13));
182  CALL_SUBTEST_4(testSingular(MatrixXd(8,8), 2e-12));
183  CALL_SUBTEST_1(testSingular(Matrix2f(), 1e-4f));
184  CALL_SUBTEST_5(testSingular(Matrix3cf(), 1e-4f));
185  CALL_SUBTEST_8(testSingular(Matrix4f(), 1e-4f));
186  CALL_SUBTEST_6(testSingular(MatrixXf(2,2), 1e-3f));
188  CALL_SUBTEST_10(testSingular(Matrix3d(), 1e-13));
189  CALL_SUBTEST_11(testSingular(Matrix3f(), 1e-4f));
191 
192  CALL_SUBTEST_2(testLogThenExp(Matrix2d(), 1e-13));
194  CALL_SUBTEST_3(testLogThenExp(Matrix4cd(), 1e-13));
195  CALL_SUBTEST_4(testLogThenExp(MatrixXd(8,8), 2e-12));
196  CALL_SUBTEST_1(testLogThenExp(Matrix2f(), 1e-4f));
197  CALL_SUBTEST_5(testLogThenExp(Matrix3cf(), 1e-4f));
198  CALL_SUBTEST_8(testLogThenExp(Matrix4f(), 1e-4f));
199  CALL_SUBTEST_6(testLogThenExp(MatrixXf(2,2), 1e-3f));
201  CALL_SUBTEST_10(testLogThenExp(Matrix3d(), 1e-13));
202  CALL_SUBTEST_11(testLogThenExp(Matrix3f(), 1e-4f));
204 }
B
Matrix< SCALARB, Dynamic, Dynamic, opt_B > B
Definition: bench_gemm.cpp:49
test2dRotation
void test2dRotation(const T &tol)
Definition: matrix_power.cpp:13
gtsam.examples.DogLegOptimizerExample.type
type
Definition: DogLegOptimizerExample.py:111
EIGEN_PI
#define EIGEN_PI
Definition: Eigen/src/Core/MathFunctions.h:16
s
RealScalar s
Definition: level1_cplx_impl.h:126
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
ceres::sin
Jet< T, N > sin(const Jet< T, N > &f)
Definition: jet.h:439
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
x
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 x
Definition: gnuplot_common_settings.hh:12
Eigen::AngleAxis
Represents a 3D rotation as a rotation angle around an arbitrary 3D axis.
Definition: ForwardDeclarations.h:290
MatrixXe
Matrix< long double, Dynamic, Dynamic > MatrixXe
Definition: matrix_power.cpp:151
m1
Matrix3d m1
Definition: IOFormat.cpp:2
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(matrix_power)
Definition: matrix_power.cpp:153
T
Eigen::Triplet< double > T
Definition: Tutorial_sparse_example.cpp:6
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:1947
CALL_SUBTEST_11
#define CALL_SUBTEST_11(FUNC)
Definition: split_test_helper.h:64
CALL_SUBTEST_9
#define CALL_SUBTEST_9(FUNC)
Definition: split_test_helper.h:52
IsComplex
@ IsComplex
Definition: gtsam/3rdparty/Eigen/blas/common.h:98
ceres::cos
Jet< T, N > cos(const Jet< T, N > &f)
Definition: jet.h:426
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
cosh
const EIGEN_DEVICE_FUNC CoshReturnType cosh() const
Definition: ArrayCwiseUnaryOps.h:353
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
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::RealSchur< MatrixType >
relerr
Derived::RealScalar relerr(const MatrixBase< Derived > &A, const MatrixBase< OtherDerived > &B)
Definition: matrix_functions.h:64
CALL_SUBTEST_10
#define CALL_SUBTEST_10(FUNC)
Definition: split_test_helper.h:58
testGeneral
void testGeneral(const MatrixType &m, const typename MatrixType::RealScalar &tol)
Definition: matrix_power.cpp:70
Eigen::MatrixPower
Class for computing matrix powers.
Definition: MatrixPower.h:15
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
L
MatrixXd L
Definition: LLT_example.cpp:6
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
gtsam.examples.DogLegOptimizerExample.run
def run(args)
Definition: DogLegOptimizerExample.py:21
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
Eigen::Triplet< double >
ceres::pow
Jet< T, N > pow(const Jet< T, N > &f, double g)
Definition: jet.h:570
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
schur
ComplexSchur< MatrixXcf > schur(4)
y
Scalar * y
Definition: level1_cplx_impl.h:124
matrix
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Definition: gtsam/3rdparty/Eigen/blas/common.h:110
matrix_functions.h
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Matrix3dRowMajor
Matrix< double, 3, 3, RowMajor > Matrix3dRowMajor
Definition: matrix_power.cpp:149
processTriangularMatrix::run
static void run(MatrixType &, MatrixType &, const MatrixType &)
Definition: matrix_functions.h:17
Eigen::ArrayBase::pow
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar, typename ExponentDerived::Scalar >, const Derived, const ExponentDerived > pow(const Eigen::ArrayBase< Derived > &x, const Eigen::ArrayBase< ExponentDerived > &exponents)
Definition: GlobalFunctions.h:143
m3
static const DiscreteKey m3(M(3), 2)
C
Matrix< Scalar, Dynamic, Dynamic > C
Definition: bench_gemm.cpp:50
test3dRotation
void test3dRotation(const T &tol)
Definition: matrix_power.cpp:56
CALL_SUBTEST_12
#define CALL_SUBTEST_12(FUNC)
Definition: split_test_helper.h:70
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
gtsam::tol
const G double tol
Definition: Group.h:79
U
@ U
Definition: testDecisionTree.cpp:349
test2dHyperbolicRotation
void test2dHyperbolicRotation(const T &tol)
Definition: matrix_power.cpp:34
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
abs
#define abs(x)
Definition: datatypes.h:17
testLogThenExp
void testLogThenExp(const MatrixType &m_const, const typename MatrixType::RealScalar &tol)
Definition: matrix_power.cpp:133
testSingular
void testSingular(const MatrixType &m_const, const typename MatrixType::RealScalar &tol)
Definition: matrix_power.cpp:100
sinh
const EIGEN_DEVICE_FUNC SinhReturnType sinh() const
Definition: ArrayCwiseUnaryOps.h:339
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
Matrix3e
Matrix< long double, 3, 3 > Matrix3e
Definition: matrix_power.cpp:150
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
VERIFY
#define VERIFY(a)
Definition: main.h:380


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:03:03