xref: /aosp_15_r20/external/eigen/unsupported/test/matrix_function.cpp (revision bf2c37156dfe67e5dfebd6d394bad8b2ab5804d4)
1*bf2c3715SXin Li // This file is part of Eigen, a lightweight C++ template library
2*bf2c3715SXin Li // for linear algebra.
3*bf2c3715SXin Li //
4*bf2c3715SXin Li // Copyright (C) 2010 Jitse Niesen <[email protected]>
5*bf2c3715SXin Li //
6*bf2c3715SXin Li // This Source Code Form is subject to the terms of the Mozilla
7*bf2c3715SXin Li // Public License v. 2.0. If a copy of the MPL was not distributed
8*bf2c3715SXin Li // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9*bf2c3715SXin Li 
10*bf2c3715SXin Li #include "main.h"
11*bf2c3715SXin Li #include <unsupported/Eigen/MatrixFunctions>
12*bf2c3715SXin Li 
13*bf2c3715SXin Li // Variant of VERIFY_IS_APPROX which uses absolute error instead of
14*bf2c3715SXin Li // relative error.
15*bf2c3715SXin Li #define VERIFY_IS_APPROX_ABS(a, b) VERIFY(test_isApprox_abs(a, b))
16*bf2c3715SXin Li 
17*bf2c3715SXin Li template<typename Type1, typename Type2>
test_isApprox_abs(const Type1 & a,const Type2 & b)18*bf2c3715SXin Li inline bool test_isApprox_abs(const Type1& a, const Type2& b)
19*bf2c3715SXin Li {
20*bf2c3715SXin Li   return ((a-b).array().abs() < test_precision<typename Type1::RealScalar>()).all();
21*bf2c3715SXin Li }
22*bf2c3715SXin Li 
23*bf2c3715SXin Li 
24*bf2c3715SXin Li // Returns a matrix with eigenvalues clustered around 0, 1 and 2.
25*bf2c3715SXin Li template<typename MatrixType>
randomMatrixWithRealEivals(const Index size)26*bf2c3715SXin Li MatrixType randomMatrixWithRealEivals(const Index size)
27*bf2c3715SXin Li {
28*bf2c3715SXin Li   typedef typename MatrixType::Scalar Scalar;
29*bf2c3715SXin Li   typedef typename MatrixType::RealScalar RealScalar;
30*bf2c3715SXin Li   MatrixType diag = MatrixType::Zero(size, size);
31*bf2c3715SXin Li   for (Index i = 0; i < size; ++i) {
32*bf2c3715SXin Li     diag(i, i) = Scalar(RealScalar(internal::random<int>(0,2)))
33*bf2c3715SXin Li       + internal::random<Scalar>() * Scalar(RealScalar(0.01));
34*bf2c3715SXin Li   }
35*bf2c3715SXin Li   MatrixType A = MatrixType::Random(size, size);
36*bf2c3715SXin Li   HouseholderQR<MatrixType> QRofA(A);
37*bf2c3715SXin Li   return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
38*bf2c3715SXin Li }
39*bf2c3715SXin Li 
40*bf2c3715SXin Li template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
41*bf2c3715SXin Li struct randomMatrixWithImagEivals
42*bf2c3715SXin Li {
43*bf2c3715SXin Li   // Returns a matrix with eigenvalues clustered around 0 and +/- i.
44*bf2c3715SXin Li   static MatrixType run(const Index size);
45*bf2c3715SXin Li };
46*bf2c3715SXin Li 
47*bf2c3715SXin Li // Partial specialization for real matrices
48*bf2c3715SXin Li template<typename MatrixType>
49*bf2c3715SXin Li struct randomMatrixWithImagEivals<MatrixType, 0>
50*bf2c3715SXin Li {
runrandomMatrixWithImagEivals51*bf2c3715SXin Li   static MatrixType run(const Index size)
52*bf2c3715SXin Li   {
53*bf2c3715SXin Li     typedef typename MatrixType::Scalar Scalar;
54*bf2c3715SXin Li     MatrixType diag = MatrixType::Zero(size, size);
55*bf2c3715SXin Li     Index i = 0;
56*bf2c3715SXin Li     while (i < size) {
57*bf2c3715SXin Li       Index randomInt = internal::random<Index>(-1, 1);
58*bf2c3715SXin Li       if (randomInt == 0 || i == size-1) {
59*bf2c3715SXin Li         diag(i, i) = internal::random<Scalar>() * Scalar(0.01);
60*bf2c3715SXin Li         ++i;
61*bf2c3715SXin Li       } else {
62*bf2c3715SXin Li         Scalar alpha = Scalar(randomInt) + internal::random<Scalar>() * Scalar(0.01);
63*bf2c3715SXin Li         diag(i, i+1) = alpha;
64*bf2c3715SXin Li         diag(i+1, i) = -alpha;
65*bf2c3715SXin Li         i += 2;
66*bf2c3715SXin Li       }
67*bf2c3715SXin Li     }
68*bf2c3715SXin Li     MatrixType A = MatrixType::Random(size, size);
69*bf2c3715SXin Li     HouseholderQR<MatrixType> QRofA(A);
70*bf2c3715SXin Li     return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
71*bf2c3715SXin Li   }
72*bf2c3715SXin Li };
73*bf2c3715SXin Li 
74*bf2c3715SXin Li // Partial specialization for complex matrices
75*bf2c3715SXin Li template<typename MatrixType>
76*bf2c3715SXin Li struct randomMatrixWithImagEivals<MatrixType, 1>
77*bf2c3715SXin Li {
runrandomMatrixWithImagEivals78*bf2c3715SXin Li   static MatrixType run(const Index size)
79*bf2c3715SXin Li   {
80*bf2c3715SXin Li     typedef typename MatrixType::Scalar Scalar;
81*bf2c3715SXin Li     typedef typename MatrixType::RealScalar RealScalar;
82*bf2c3715SXin Li     const Scalar imagUnit(0, 1);
83*bf2c3715SXin Li     MatrixType diag = MatrixType::Zero(size, size);
84*bf2c3715SXin Li     for (Index i = 0; i < size; ++i) {
85*bf2c3715SXin Li       diag(i, i) = Scalar(RealScalar(internal::random<Index>(-1, 1))) * imagUnit
86*bf2c3715SXin Li         + internal::random<Scalar>() * Scalar(RealScalar(0.01));
87*bf2c3715SXin Li     }
88*bf2c3715SXin Li     MatrixType A = MatrixType::Random(size, size);
89*bf2c3715SXin Li     HouseholderQR<MatrixType> QRofA(A);
90*bf2c3715SXin Li     return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
91*bf2c3715SXin Li   }
92*bf2c3715SXin Li };
93*bf2c3715SXin Li 
94*bf2c3715SXin Li 
95*bf2c3715SXin Li template<typename MatrixType>
testMatrixExponential(const MatrixType & A)96*bf2c3715SXin Li void testMatrixExponential(const MatrixType& A)
97*bf2c3715SXin Li {
98*bf2c3715SXin Li   typedef typename internal::traits<MatrixType>::Scalar Scalar;
99*bf2c3715SXin Li   typedef typename NumTraits<Scalar>::Real RealScalar;
100*bf2c3715SXin Li   typedef std::complex<RealScalar> ComplexScalar;
101*bf2c3715SXin Li 
102*bf2c3715SXin Li   VERIFY_IS_APPROX(A.exp(), A.matrixFunction(internal::stem_function_exp<ComplexScalar>));
103*bf2c3715SXin Li }
104*bf2c3715SXin Li 
105*bf2c3715SXin Li template<typename MatrixType>
testMatrixLogarithm(const MatrixType & A)106*bf2c3715SXin Li void testMatrixLogarithm(const MatrixType& A)
107*bf2c3715SXin Li {
108*bf2c3715SXin Li   typedef typename internal::traits<MatrixType>::Scalar Scalar;
109*bf2c3715SXin Li   typedef typename NumTraits<Scalar>::Real RealScalar;
110*bf2c3715SXin Li 
111*bf2c3715SXin Li   MatrixType scaledA;
112*bf2c3715SXin Li   RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff();
113*bf2c3715SXin Li   if (maxImagPartOfSpectrum >= RealScalar(0.9L * EIGEN_PI))
114*bf2c3715SXin Li     scaledA = A * RealScalar(0.9L * EIGEN_PI) / maxImagPartOfSpectrum;
115*bf2c3715SXin Li   else
116*bf2c3715SXin Li     scaledA = A;
117*bf2c3715SXin Li 
118*bf2c3715SXin Li   // identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X
119*bf2c3715SXin Li   MatrixType expA = scaledA.exp();
120*bf2c3715SXin Li   MatrixType logExpA = expA.log();
121*bf2c3715SXin Li   VERIFY_IS_APPROX(logExpA, scaledA);
122*bf2c3715SXin Li }
123*bf2c3715SXin Li 
124*bf2c3715SXin Li template<typename MatrixType>
testHyperbolicFunctions(const MatrixType & A)125*bf2c3715SXin Li void testHyperbolicFunctions(const MatrixType& A)
126*bf2c3715SXin Li {
127*bf2c3715SXin Li   // Need to use absolute error because of possible cancellation when
128*bf2c3715SXin Li   // adding/subtracting expA and expmA.
129*bf2c3715SXin Li   VERIFY_IS_APPROX_ABS(A.sinh(), (A.exp() - (-A).exp()) / 2);
130*bf2c3715SXin Li   VERIFY_IS_APPROX_ABS(A.cosh(), (A.exp() + (-A).exp()) / 2);
131*bf2c3715SXin Li }
132*bf2c3715SXin Li 
133*bf2c3715SXin Li template<typename MatrixType>
testGonioFunctions(const MatrixType & A)134*bf2c3715SXin Li void testGonioFunctions(const MatrixType& A)
135*bf2c3715SXin Li {
136*bf2c3715SXin Li   typedef typename MatrixType::Scalar Scalar;
137*bf2c3715SXin Li   typedef typename NumTraits<Scalar>::Real RealScalar;
138*bf2c3715SXin Li   typedef std::complex<RealScalar> ComplexScalar;
139*bf2c3715SXin Li   typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime,
140*bf2c3715SXin Li                  MatrixType::ColsAtCompileTime, MatrixType::Options> ComplexMatrix;
141*bf2c3715SXin Li 
142*bf2c3715SXin Li   ComplexScalar imagUnit(0,1);
143*bf2c3715SXin Li   ComplexScalar two(2,0);
144*bf2c3715SXin Li 
145*bf2c3715SXin Li   ComplexMatrix Ac = A.template cast<ComplexScalar>();
146*bf2c3715SXin Li 
147*bf2c3715SXin Li   ComplexMatrix exp_iA = (imagUnit * Ac).exp();
148*bf2c3715SXin Li   ComplexMatrix exp_miA = (-imagUnit * Ac).exp();
149*bf2c3715SXin Li 
150*bf2c3715SXin Li   ComplexMatrix sinAc = A.sin().template cast<ComplexScalar>();
151*bf2c3715SXin Li   VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit));
152*bf2c3715SXin Li 
153*bf2c3715SXin Li   ComplexMatrix cosAc = A.cos().template cast<ComplexScalar>();
154*bf2c3715SXin Li   VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2);
155*bf2c3715SXin Li }
156*bf2c3715SXin Li 
157*bf2c3715SXin Li template<typename MatrixType>
testMatrix(const MatrixType & A)158*bf2c3715SXin Li void testMatrix(const MatrixType& A)
159*bf2c3715SXin Li {
160*bf2c3715SXin Li   testMatrixExponential(A);
161*bf2c3715SXin Li   testMatrixLogarithm(A);
162*bf2c3715SXin Li   testHyperbolicFunctions(A);
163*bf2c3715SXin Li   testGonioFunctions(A);
164*bf2c3715SXin Li }
165*bf2c3715SXin Li 
166*bf2c3715SXin Li template<typename MatrixType>
testMatrixType(const MatrixType & m)167*bf2c3715SXin Li void testMatrixType(const MatrixType& m)
168*bf2c3715SXin Li {
169*bf2c3715SXin Li   // Matrices with clustered eigenvalue lead to different code paths
170*bf2c3715SXin Li   // in MatrixFunction.h and are thus useful for testing.
171*bf2c3715SXin Li 
172*bf2c3715SXin Li   const Index size = m.rows();
173*bf2c3715SXin Li   for (int i = 0; i < g_repeat; i++) {
174*bf2c3715SXin Li     testMatrix(MatrixType::Random(size, size).eval());
175*bf2c3715SXin Li     testMatrix(randomMatrixWithRealEivals<MatrixType>(size));
176*bf2c3715SXin Li     testMatrix(randomMatrixWithImagEivals<MatrixType>::run(size));
177*bf2c3715SXin Li   }
178*bf2c3715SXin Li }
179*bf2c3715SXin Li 
180*bf2c3715SXin Li template<typename MatrixType>
testMapRef(const MatrixType & A)181*bf2c3715SXin Li void testMapRef(const MatrixType& A)
182*bf2c3715SXin Li {
183*bf2c3715SXin Li   // Test if passing Ref and Map objects is possible
184*bf2c3715SXin Li   // (Regression test for Bug #1796)
185*bf2c3715SXin Li   Index size = A.rows();
186*bf2c3715SXin Li   MatrixType X; X.setRandom(size, size);
187*bf2c3715SXin Li   MatrixType Y(size,size);
188*bf2c3715SXin Li   Ref<      MatrixType> R(Y);
189*bf2c3715SXin Li   Ref<const MatrixType> Rc(X);
190*bf2c3715SXin Li   Map<      MatrixType> M(Y.data(), size, size);
191*bf2c3715SXin Li   Map<const MatrixType> Mc(X.data(), size, size);
192*bf2c3715SXin Li 
193*bf2c3715SXin Li   X = X*X; // make sure sqrt is possible
194*bf2c3715SXin Li   Y = X.sqrt();
195*bf2c3715SXin Li   R = Rc.sqrt();
196*bf2c3715SXin Li   M = Mc.sqrt();
197*bf2c3715SXin Li   Y = X.exp();
198*bf2c3715SXin Li   R = Rc.exp();
199*bf2c3715SXin Li   M = Mc.exp();
200*bf2c3715SXin Li   X = Y; // make sure log is possible
201*bf2c3715SXin Li   Y = X.log();
202*bf2c3715SXin Li   R = Rc.log();
203*bf2c3715SXin Li   M = Mc.log();
204*bf2c3715SXin Li 
205*bf2c3715SXin Li   Y = X.cos() + Rc.cos() + Mc.cos();
206*bf2c3715SXin Li   Y = X.sin() + Rc.sin() + Mc.sin();
207*bf2c3715SXin Li 
208*bf2c3715SXin Li   Y = X.cosh() + Rc.cosh() + Mc.cosh();
209*bf2c3715SXin Li   Y = X.sinh() + Rc.sinh() + Mc.sinh();
210*bf2c3715SXin Li }
211*bf2c3715SXin Li 
212*bf2c3715SXin Li 
EIGEN_DECLARE_TEST(matrix_function)213*bf2c3715SXin Li EIGEN_DECLARE_TEST(matrix_function)
214*bf2c3715SXin Li {
215*bf2c3715SXin Li   CALL_SUBTEST_1(testMatrixType(Matrix<float,1,1>()));
216*bf2c3715SXin Li   CALL_SUBTEST_2(testMatrixType(Matrix3cf()));
217*bf2c3715SXin Li   CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8)));
218*bf2c3715SXin Li   CALL_SUBTEST_4(testMatrixType(Matrix2d()));
219*bf2c3715SXin Li   CALL_SUBTEST_5(testMatrixType(Matrix<double,5,5,RowMajor>()));
220*bf2c3715SXin Li   CALL_SUBTEST_6(testMatrixType(Matrix4cd()));
221*bf2c3715SXin Li   CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13)));
222*bf2c3715SXin Li 
223*bf2c3715SXin Li   CALL_SUBTEST_1(testMapRef(Matrix<float,1,1>()));
224*bf2c3715SXin Li   CALL_SUBTEST_2(testMapRef(Matrix3cf()));
225*bf2c3715SXin Li   CALL_SUBTEST_3(testMapRef(MatrixXf(8,8)));
226*bf2c3715SXin Li   CALL_SUBTEST_7(testMapRef(MatrixXd(13,13)));
227*bf2c3715SXin Li }
228