source: trunk/test/kernel_pca_test.cc @ 2324

Last change on this file since 2324 was 2324, checked in by Peter, 12 years ago

new class KernelPCA. closes #639

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 4.3 KB
Line 
1// $Id: kernel_pca_test.cc 2324 2010-09-21 21:19:25Z peter $
2
3/*
4  Copyright (C) 2010 Peter Johansson
5
6  This file is part of the yat library, http://dev.thep.lu.se/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 3 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with yat. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#include "Suite.h"
23
24#include "yat/normalizer/Centralizer.h"
25#include "yat/normalizer/RowNormalizer.h"
26
27#include "yat/utility/VectorConstView.h"
28#include "yat/utility/Matrix.h"
29#include "yat/utility/PCA.h"
30#include "yat/utility/KernelPCA.h"
31#include "yat/utility/SVD.h"
32
33#include <gsl/gsl_linalg.h>
34
35#include <fstream>
36#include <string>
37
38int main(int argc, char* argv[])
39{
40  using namespace theplu::yat;
41  test::Suite suite(argc, argv);
42
43  std::ifstream is(test::filename("data/nm_data_centralized.txt").c_str());
44  const utility::Matrix data2(is);
45  is.close();
46
47  utility::Matrix data(3,3);
48  for (size_t i=0; i<data.rows(); ++i)
49    for (size_t j=0; j<data.columns(); ++j)
50      data(i,j) = data2(i,j);
51
52  utility::PCA pca(data);
53
54  normalizer::RowNormalizer<normalizer::Centralizer<> > normalizer;
55  utility::Matrix data_centered(data);
56  normalizer(data, data_centered);
57
58
59  // calculate kernel = data' * data
60  utility::Matrix kernel = data_centered;
61  kernel.transpose();
62  kernel *= data_centered;
63  //  kernel *= 1.0/data.rows();
64
65  suite.out() << "check that kernel is symmetric\n";
66  for (size_t i=0; i<kernel.rows(); ++i)
67    for (size_t j=i; j<kernel.columns(); ++j)
68      if (!suite.equal(kernel(i,j), kernel(j,i)))
69        suite.out() << i << " x " << j << "\n";
70
71  //suite.out() << "data:\n" << data_centered << "\n";
72
73  utility::Matrix Z(kernel);
74  // avoid problem with zero eigenvalue
75  Z += 1.0e-10;
76  // kernel = Z Z', where Z is lower-triangle matrix
77  gsl_linalg_cholesky_decomp(Z.gsl_matrix_p());
78  // mask out lower part
79  for (size_t row=0; row<Z.rows(); ++row)
80    for (size_t column=0; column<row; ++column)
81      Z(row, column) = 0.0;
82
83  normalizer(Z, Z);
84  utility::KernelPCA pca2(kernel);
85
86  suite.out() << "check eigenvalues\n";
87  if (pca.eigenvalues().size() != pca2.eigenvalues().size()) {
88    suite.out() << "eigenvalues size incorrect\n";
89    suite.add(false);
90    // skip remaining tests if we are this fundamentally wrong
91    return suite.return_value();
92  }
93
94  // don't compare the last (close to zero) eigenvalues
95  for (size_t i=0; i<pca.eigenvalues().size()-1; ++i) {
96    suite.equal(pca2.eigenvalues()(i), pca.eigenvalues()(i), 1000);
97  }
98  // check that last eigenvalue is zero
99  suite.equal_fix(pca2.eigenvalues()(pca.eigenvalues().size()-1), 0, 1e-15);
100
101  suite.out() << "compare projection\n";
102  utility::Matrix project = pca.projection(data);
103  utility::Matrix project2 = pca2.projection();
104  if (project.rows() == project2.rows() && 
105      project.columns() == project2.columns()) {
106
107    // do not check last eigenvector
108    for (size_t i=0; i<project.rows()-1; ++i) {
109      utility::VectorConstView v1 = project.row_const_view(i);
110      utility::VectorConstView v2 = project2.row_const_view(i);
111      double r = (v1*v2) * (v1*v2) / ((v1*v1) * (v2*v2));
112      suite.add(suite.equal(r, 1.0, 10));
113    }
114    if (!suite.ok()) {
115      suite.out() << "project:\n" << project << "\n"; 
116      suite.out() << "project2:\n" << project2 << "\n"; 
117    }
118  }
119  else {
120    suite.add(false);
121    suite.out() << "projection dimension disagree\n"
122                << "pca: " << project.rows() << " x " 
123                << project.columns() << "\n"
124                << "pca2: " << project2.rows() << " x " 
125                << project2.columns() << "\n";
126  }
127
128  // some test for visual inspection!
129  utility::Matrix x(3,3);
130  x(0,0) = 4;
131  x(0,1) = 6;
132  x(0,2) = -10;
133  x(1,0) = 6;
134  x(1,1) = 4;
135  x(1,2) = -10;
136  x(2,0) = 5;
137  x(2,1) = 5;
138  x(2,2) = -10;
139
140  utility::Matrix k(x);
141  k.transpose();
142  k *= x;
143  suite.out() << k << "\n";
144  utility::KernelPCA pca3(k);
145
146  suite.out() << pca3.projection() << "\n";
147 
148  utility::Matrix new_k = pca3.projection();
149  new_k.transpose();
150  new_k *= pca3.projection();
151  suite.out() << new_k << "\n";
152
153  return suite.return_value();
154}
155
Note: See TracBrowser for help on using the repository browser.