coco
246 lines
8.1 KiB
C++
246 lines
8.1 KiB
C++
// 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;
|
|
|
|
}
|
|
|
|
|
|
|