diff options
Diffstat (limited to 'ml/dlib/dlib/test/matrix_eig.cpp')
| -rw-r--r-- | ml/dlib/dlib/test/matrix_eig.cpp | 245 |
1 files changed, 245 insertions, 0 deletions
diff --git a/ml/dlib/dlib/test/matrix_eig.cpp b/ml/dlib/dlib/test/matrix_eig.cpp new file mode 100644 index 000000000..9fbce6598 --- /dev/null +++ b/ml/dlib/dlib/test/matrix_eig.cpp @@ -0,0 +1,245 @@ +// Copyright (C) 2009 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. + + +#include <dlib/matrix.h> +#include <sstream> +#include <string> +#include <cstdlib> +#include <ctime> +#include <vector> +#include "../stl_checked.h" +#include "../array.h" +#include "../rand.h" +#include <dlib/string.h> + +#include "tester.h" + +namespace +{ + + using namespace test; + using namespace dlib; + using namespace std; + + logger dlog("test.matrix_eig"); + + dlib::rand rnd; + +// ---------------------------------------------------------------------------------------- + + template <typename type> + const matrix<type> randm(long r, long c) + { + matrix<type> m(r,c); + for (long row = 0; row < m.nr(); ++row) + { + for (long col = 0; col < m.nc(); ++col) + { + m(row,col) = static_cast<type>(rnd.get_random_double()); + } + } + + return m; + } + + template <typename type, long NR, long NC> + const matrix<type,NR,NC> randm() + { + matrix<type,NR,NC> m; + for (long row = 0; row < m.nr(); ++row) + { + for (long col = 0; col < m.nc(); ++col) + { + m(row,col) = static_cast<type>(rnd.get_random_double()); + } + } + + return m; + } + +// ---------------------------------------------------------------------------------------- + + template <typename matrix_type, typename U> + void test_eigenvalue_impl ( const matrix_type& m, const eigenvalue_decomposition<U>& test ) + { + typedef typename matrix_type::type type; + const type eps = 10*max(abs(m))*sqrt(std::numeric_limits<type>::epsilon()); + dlog << LDEBUG << "test_eigenvalue(): " << m.nr() << " x " << m.nc() << " eps: " << eps; + print_spinner(); + + + DLIB_TEST(test.dim() == m.nr()); + + // make sure all the various ways of asking for the eigenvalues are actually returning a + // consistent set of eigenvalues. + DLIB_TEST(equal(real(test.get_eigenvalues()), test.get_real_eigenvalues(), eps)); + DLIB_TEST(equal(imag(test.get_eigenvalues()), test.get_imag_eigenvalues(), eps)); + DLIB_TEST(equal(real(diag(test.get_d())), test.get_real_eigenvalues(), eps)); + DLIB_TEST(equal(imag(diag(test.get_d())), test.get_imag_eigenvalues(), eps)); + + matrix<type> eig1 ( real_eigenvalues(m)); + matrix<type> eig2 ( test.get_real_eigenvalues()); + sort(&eig1(0), &eig1(0) + eig1.size()); + sort(&eig2(0), &eig2(0) + eig2.size()); + DLIB_TEST(max(abs(eig1 - eig2)) < eps); + + const matrix<type> V = test.get_pseudo_v(); + const matrix<type> D = test.get_pseudo_d(); + const matrix<complex<type> > CV = test.get_v(); + const matrix<complex<type> > CD = test.get_d(); + const matrix<complex<type> > CM = complex_matrix(m, uniform_matrix<type>(m.nr(),m.nc(),0)); + + DLIB_TEST(V.nr() == test.dim()); + DLIB_TEST(V.nc() == test.dim()); + DLIB_TEST(D.nr() == test.dim()); + DLIB_TEST(D.nc() == test.dim()); + + // CD is a diagonal matrix + DLIB_TEST(diagm(diag(CD)) == CD); + + // verify that these things are actually eigenvalues and eigenvectors of m + DLIB_TEST_MSG(max(abs(m*V - V*D)) < eps, max(abs(m*V - V*D)) << " " << eps); + DLIB_TEST(max(norm(CM*CV - CV*CD)) < eps); + + // if m is a symmetric matrix + if (max(abs(m-trans(m))) < 1e-5) + { + dlog << LTRACE << "m is symmetric"; + // there aren't any imaginary eigenvalues + DLIB_TEST(max(abs(test.get_imag_eigenvalues())) < eps); + DLIB_TEST(diagm(diag(D)) == D); + + // only check the determinant against the eigenvalues for small matrices + // because for huge ones the determinant might be so big it overflows a floating point number. + if (m.nr() < 50) + { + const type mdet = det(m); + DLIB_TEST_MSG(std::abs(prod(test.get_real_eigenvalues()) - mdet) < std::abs(mdet)*sqrt(std::numeric_limits<type>::epsilon()), + std::abs(prod(test.get_real_eigenvalues()) - mdet) <<" eps: " << std::abs(mdet)*sqrt(std::numeric_limits<type>::epsilon()) + << " mdet: "<< mdet << " prod(eig): " << prod(test.get_real_eigenvalues()) + ); + } + + // V is orthogonal + DLIB_TEST(equal(V*trans(V), identity_matrix<type>(test.dim()), eps)); + DLIB_TEST(equal(m , V*D*trans(V), eps)); + } + else + { + dlog << LTRACE << "m is NOT symmetric"; + DLIB_TEST_MSG(equal(m , V*D*inv(V), eps), max(abs(m - V*D*inv(V)))); + } + } + +// ---------------------------------------------------------------------------------------- + + template <typename matrix_type> + void test_eigenvalue ( const matrix_type& m ) + { + typedef typename matrix_type::type type; + typedef typename matrix_type::mem_manager_type MM; + matrix<type,matrix_type::NR, matrix_type::NC, MM, row_major_layout> mr(m); + matrix<type,matrix_type::NR, matrix_type::NC, MM, column_major_layout> mc(m); + + { + eigenvalue_decomposition<matrix_type> test(mr); + test_eigenvalue_impl(mr, test); + + eigenvalue_decomposition<matrix_type> test_symm(make_symmetric(mr)); + test_eigenvalue_impl(make_symmetric(mr), test_symm); + } + + { + eigenvalue_decomposition<matrix_type> test(mc); + test_eigenvalue_impl(mc, test); + + eigenvalue_decomposition<matrix_type> test_symm(make_symmetric(mc)); + test_eigenvalue_impl(make_symmetric(mc), test_symm); + } + } + +// ---------------------------------------------------------------------------------------- + + void matrix_test_double() + { + + test_eigenvalue(10*randm<double>(1,1)); + test_eigenvalue(10*randm<double>(2,2)); + test_eigenvalue(10*randm<double>(3,3)); + test_eigenvalue(10*randm<double>(4,4)); + test_eigenvalue(10*randm<double>(15,15)); + test_eigenvalue(10*randm<double>(150,150)); + + test_eigenvalue(10*randm<double,1,1>()); + test_eigenvalue(10*randm<double,2,2>()); + test_eigenvalue(10*randm<double,3,3>()); + } + +// ---------------------------------------------------------------------------------------- + + void matrix_test_float() + { + + test_eigenvalue(10*randm<float>(1,1)); + test_eigenvalue(10*randm<float>(2,2)); + test_eigenvalue(10*randm<float>(3,3)); + test_eigenvalue(10*randm<float>(4,4)); + test_eigenvalue(10*randm<float>(15,15)); + test_eigenvalue(10*randm<float>(50,50)); + + test_eigenvalue(10*randm<float,1,1>()); + test_eigenvalue(10*randm<float,2,2>()); + test_eigenvalue(10*randm<float,3,3>()); + } + + template <int dims> + void test_eigenvalue2() + { + for (int seed = 0; seed < 10; ++seed) + { + print_spinner(); + matrix<double> H = gaussian_randm(dims,dims,seed); + H = H*trans(H); + + eigenvalue_decomposition<matrix<double> > eig(H); + matrix<double> HH = eig.get_pseudo_v()*diagm(eig.get_real_eigenvalues())*trans(eig.get_pseudo_v()); + DLIB_TEST_MSG(max(abs(H - HH))<1e-12, "dims: " << dims << " error: " << max(abs(H - HH))); + } + } + +// ---------------------------------------------------------------------------------------- + + class matrix_tester : public tester + { + public: + matrix_tester ( + ) : + tester ("test_matrix_eig", + "Runs tests on the matrix eigen decomp component.") + { + //rnd.set_seed(cast_to_string(time(0))); + } + + void perform_test ( + ) + { + dlog << LINFO << "seed string: " << rnd.get_seed(); + + dlog << LINFO << "begin testing with double"; + matrix_test_double(); + dlog << LINFO << "begin testing with float"; + matrix_test_float(); + + test_eigenvalue2<10>(); + test_eigenvalue2<11>(); + test_eigenvalue2<3>(); + test_eigenvalue2<2>(); + test_eigenvalue2<1>(); + } + } a; + +} + + + |