test/kernel_pca.cc

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

Generated on Sat Oct 14 02:21:05 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)