1 // Copyright (C) 2009 Davis E. King (davis@dlib.net) 2 // License: Boost Software License See LICENSE.txt for the full license. 3 4 5 #include <dlib/matrix.h> 6 #include <sstream> 7 #include <string> 8 #include <cstdlib> 9 #include <ctime> 10 #include <vector> 11 #include "../stl_checked.h" 12 #include "../array.h" 13 #include "../rand.h" 14 #include <dlib/string.h> 15 16 #include "tester.h" 17 18 namespace 19 { 20 21 using namespace test; 22 using namespace dlib; 23 using namespace std; 24 25 logger dlog("test.matrix_eig"); 26 27 dlib::rand rnd; 28 29 // ---------------------------------------------------------------------------------------- 30 31 template <typename type> randm(long r,long c)32 const matrix<type> randm(long r, long c) 33 { 34 matrix<type> m(r,c); 35 for (long row = 0; row < m.nr(); ++row) 36 { 37 for (long col = 0; col < m.nc(); ++col) 38 { 39 m(row,col) = static_cast<type>(rnd.get_random_double()); 40 } 41 } 42 43 return m; 44 } 45 46 template <typename type, long NR, long NC> randm()47 const matrix<type,NR,NC> randm() 48 { 49 matrix<type,NR,NC> m; 50 for (long row = 0; row < m.nr(); ++row) 51 { 52 for (long col = 0; col < m.nc(); ++col) 53 { 54 m(row,col) = static_cast<type>(rnd.get_random_double()); 55 } 56 } 57 58 return m; 59 } 60 61 // ---------------------------------------------------------------------------------------- 62 63 template <typename matrix_type, typename U> test_eigenvalue_impl(const matrix_type & m,const eigenvalue_decomposition<U> & test)64 void test_eigenvalue_impl ( const matrix_type& m, const eigenvalue_decomposition<U>& test ) 65 { 66 typedef typename matrix_type::type type; 67 const type eps = 10*max(abs(m))*sqrt(std::numeric_limits<type>::epsilon()); 68 dlog << LDEBUG << "test_eigenvalue(): " << m.nr() << " x " << m.nc() << " eps: " << eps; 69 print_spinner(); 70 71 72 DLIB_TEST(test.dim() == m.nr()); 73 74 // make sure all the various ways of asking for the eigenvalues are actually returning a 75 // consistent set of eigenvalues. 76 DLIB_TEST(equal(real(test.get_eigenvalues()), test.get_real_eigenvalues(), eps)); 77 DLIB_TEST(equal(imag(test.get_eigenvalues()), test.get_imag_eigenvalues(), eps)); 78 DLIB_TEST(equal(real(diag(test.get_d())), test.get_real_eigenvalues(), eps)); 79 DLIB_TEST(equal(imag(diag(test.get_d())), test.get_imag_eigenvalues(), eps)); 80 81 matrix<type> eig1 ( real_eigenvalues(m)); 82 matrix<type> eig2 ( test.get_real_eigenvalues()); 83 sort(&eig1(0), &eig1(0) + eig1.size()); 84 sort(&eig2(0), &eig2(0) + eig2.size()); 85 DLIB_TEST(max(abs(eig1 - eig2)) < eps); 86 87 const matrix<type> V = test.get_pseudo_v(); 88 const matrix<type> D = test.get_pseudo_d(); 89 const matrix<complex<type> > CV = test.get_v(); 90 const matrix<complex<type> > CD = test.get_d(); 91 const matrix<complex<type> > CM = complex_matrix(m, uniform_matrix<type>(m.nr(),m.nc(),0)); 92 93 DLIB_TEST(V.nr() == test.dim()); 94 DLIB_TEST(V.nc() == test.dim()); 95 DLIB_TEST(D.nr() == test.dim()); 96 DLIB_TEST(D.nc() == test.dim()); 97 98 // CD is a diagonal matrix 99 DLIB_TEST(diagm(diag(CD)) == CD); 100 101 // verify that these things are actually eigenvalues and eigenvectors of m 102 DLIB_TEST_MSG(max(abs(m*V - V*D)) < eps, max(abs(m*V - V*D)) << " " << eps); 103 DLIB_TEST(max(norm(CM*CV - CV*CD)) < eps); 104 105 // if m is a symmetric matrix 106 if (max(abs(m-trans(m))) < 1e-5) 107 { 108 dlog << LTRACE << "m is symmetric"; 109 // there aren't any imaginary eigenvalues 110 DLIB_TEST(max(abs(test.get_imag_eigenvalues())) < eps); 111 DLIB_TEST(diagm(diag(D)) == D); 112 113 // only check the determinant against the eigenvalues for small matrices 114 // because for huge ones the determinant might be so big it overflows a floating point number. 115 if (m.nr() < 50) 116 { 117 const type mdet = det(m); 118 DLIB_TEST_MSG(std::abs(prod(test.get_real_eigenvalues()) - mdet) < std::abs(mdet)*sqrt(std::numeric_limits<type>::epsilon()), 119 std::abs(prod(test.get_real_eigenvalues()) - mdet) <<" eps: " << std::abs(mdet)*sqrt(std::numeric_limits<type>::epsilon()) 120 << " mdet: "<< mdet << " prod(eig): " << prod(test.get_real_eigenvalues()) 121 ); 122 } 123 124 // V is orthogonal 125 DLIB_TEST(equal(V*trans(V), identity_matrix<type>(test.dim()), eps)); 126 DLIB_TEST(equal(m , V*D*trans(V), eps)); 127 } 128 else 129 { 130 dlog << LTRACE << "m is NOT symmetric"; 131 DLIB_TEST_MSG(equal(m , V*D*inv(V), eps), max(abs(m - V*D*inv(V)))); 132 } 133 } 134 135 // ---------------------------------------------------------------------------------------- 136 137 template <typename matrix_type> test_eigenvalue(const matrix_type & m)138 void test_eigenvalue ( const matrix_type& m ) 139 { 140 typedef typename matrix_type::type type; 141 typedef typename matrix_type::mem_manager_type MM; 142 matrix<type,matrix_type::NR, matrix_type::NC, MM, row_major_layout> mr(m); 143 matrix<type,matrix_type::NR, matrix_type::NC, MM, column_major_layout> mc(m); 144 145 { 146 eigenvalue_decomposition<matrix_type> test(mr); 147 test_eigenvalue_impl(mr, test); 148 149 eigenvalue_decomposition<matrix_type> test_symm(make_symmetric(mr)); 150 test_eigenvalue_impl(make_symmetric(mr), test_symm); 151 } 152 153 { 154 eigenvalue_decomposition<matrix_type> test(mc); 155 test_eigenvalue_impl(mc, test); 156 157 eigenvalue_decomposition<matrix_type> test_symm(make_symmetric(mc)); 158 test_eigenvalue_impl(make_symmetric(mc), test_symm); 159 } 160 } 161 162 // ---------------------------------------------------------------------------------------- 163 matrix_test_double()164 void matrix_test_double() 165 { 166 167 test_eigenvalue(10*randm<double>(1,1)); 168 test_eigenvalue(10*randm<double>(2,2)); 169 test_eigenvalue(10*randm<double>(3,3)); 170 test_eigenvalue(10*randm<double>(4,4)); 171 test_eigenvalue(10*randm<double>(15,15)); 172 test_eigenvalue(10*randm<double>(150,150)); 173 174 test_eigenvalue(10*randm<double,1,1>()); 175 test_eigenvalue(10*randm<double,2,2>()); 176 test_eigenvalue(10*randm<double,3,3>()); 177 } 178 179 // ---------------------------------------------------------------------------------------- 180 matrix_test_float()181 void matrix_test_float() 182 { 183 184 test_eigenvalue(10*randm<float>(1,1)); 185 test_eigenvalue(10*randm<float>(2,2)); 186 test_eigenvalue(10*randm<float>(3,3)); 187 test_eigenvalue(10*randm<float>(4,4)); 188 test_eigenvalue(10*randm<float>(15,15)); 189 test_eigenvalue(10*randm<float>(50,50)); 190 191 test_eigenvalue(10*randm<float,1,1>()); 192 test_eigenvalue(10*randm<float,2,2>()); 193 test_eigenvalue(10*randm<float,3,3>()); 194 } 195 196 template <int dims> test_eigenvalue2()197 void test_eigenvalue2() 198 { 199 for (int seed = 0; seed < 10; ++seed) 200 { 201 print_spinner(); 202 matrix<double> H = gaussian_randm(dims,dims,seed); 203 H = H*trans(H); 204 205 eigenvalue_decomposition<matrix<double> > eig(H); 206 matrix<double> HH = eig.get_pseudo_v()*diagm(eig.get_real_eigenvalues())*trans(eig.get_pseudo_v()); 207 DLIB_TEST_MSG(max(abs(H - HH))<1e-12, "dims: " << dims << " error: " << max(abs(H - HH))); 208 } 209 } 210 211 // ---------------------------------------------------------------------------------------- 212 213 class matrix_tester : public tester 214 { 215 public: matrix_tester()216 matrix_tester ( 217 ) : 218 tester ("test_matrix_eig", 219 "Runs tests on the matrix eigen decomp component.") 220 { 221 //rnd.set_seed(cast_to_string(time(0))); 222 } 223 perform_test()224 void perform_test ( 225 ) 226 { 227 dlog << LINFO << "seed string: " << rnd.get_seed(); 228 229 dlog << LINFO << "begin testing with double"; 230 matrix_test_double(); 231 dlog << LINFO << "begin testing with float"; 232 matrix_test_float(); 233 234 test_eigenvalue2<10>(); 235 test_eigenvalue2<11>(); 236 test_eigenvalue2<3>(); 237 test_eigenvalue2<2>(); 238 test_eigenvalue2<1>(); 239 } 240 } a; 241 242 } 243 244 245 246