source: trunk/yat/utility/PCA.cc @ 1009

Last change on this file since 1009 was 1009, checked in by Peter, 13 years ago

merging branch peter-dev into trunk delta 1008:994

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.5 KB
Line 
1// $Id: PCA.cc 1009 2008-02-01 04:15:44Z peter $
2
3/*
4  Copyright (C) 2003 Daniel Dalevi
5  Copyright (C) 2004 Jari Häkkinen
6  Copyright (C) 2005 Jari Häkkinen, Peter Johansson
7  Copyright (C) 2006 Jari Häkkinen
8  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
9
10  This file is part of the yat library, http://trac.thep.lu.se/yat
11
12  The yat library is free software; you can redistribute it and/or
13  modify it under the terms of the GNU General Public License as
14  published by the Free Software Foundation; either version 2 of the
15  License, or (at your option) any later version.
16
17  The yat library is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20  General Public License for more details.
21
22  You should have received a copy of the GNU General Public License
23  along with this program; if not, write to the Free Software
24  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
25  02111-1307, USA.
26*/
27
28#include "PCA.h"
29#include "SVD.h"
30#include "utility.h"
31#include "vectorView.h"
32
33#include <iostream>
34#include <cmath>
35#include <memory>
36
37namespace theplu {
38namespace yat {
39namespace utility {
40
41
42  PCA::PCA(const utility::matrix& A)
43    : A_(A)
44  {
45    process();
46  }
47
48
49  const utility::vector& PCA::eigenvalues(void) const
50  {
51    return eigenvalues_;
52  }
53
54
55  const utility::matrix& PCA::eigenvectors(void) const
56  {
57    return eigenvectors_;
58  }
59
60
61  void PCA::process(void)
62  {
63    // Row-center the data matrix
64    utility::matrix A_center( A_.rows(), A_.columns() );
65    this->row_center( A_center );
66
67    // Single value decompose the data matrix
68    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
69    pSVD->decompose();
70    utility::matrix U(pSVD->U());
71    utility::matrix V(pSVD->V());
72
73    // Read the eigenvectors and eigenvalues
74    eigenvectors_.clone(U);
75    eigenvectors_ .transpose();
76    eigenvalues_ = pSVD->s();
77
78    // T
79    for( size_t i = 0; i < eigenvalues_.size(); ++i )
80      eigenvalues_(i) = eigenvalues_(i)*eigenvalues_(i);
81    eigenvalues_ *= 1.0/A_center.rows();
82
83    // Sort the eigenvectors in order of eigenvalues
84    // Simple (not efficient) algorithm that always
85    // make sure that the i:th element is in its correct
86    // position (N element --> Ordo( N*N ))
87    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
88      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
89        if ( eigenvalues_(j) > eigenvalues_(i) ) {
90          std::swap( eigenvalues_(i), eigenvalues_(j) ); 
91          eigenvectors_.swap_rows(i,j);
92        }
93  }
94
95
96  /*
97  void PCA::process_transposed_problem(void)
98  {
99    // Row-center the data matrix
100    utility::matrix A_center( A_.rows(), A_.columns() );
101    this->row_center( A_center );
102
103    // Transform into SVD friendly dimension
104    A_.transpose();
105    A_center.transpose();
106
107    // Single value decompose the data matrix
108    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
109    pSVD->decompose();
110    utility::matrix U(pSVD->U());
111    utility::matrix V(pSVD->V());
112
113    // Read the eigenvectors and eigenvalues
114    eigenvectors_.clone(V);
115    eigenvectors_.transpose();
116    eigenvalues_.clone(pSVD->s());
117
118    // Transform back when done with SVD!
119    // (used V insted of U now for eigenvectors)
120    A_.transpose();
121    A_center.transpose();
122
123    // T
124    for( size_t i = 0; i < eigenvalues_.size(); ++i )
125      eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
126    eigenvalues_ *= 1.0/A_center.rows();
127
128    // Sort the eigenvectors in order of eigenvalues
129    // Simple (not efficient) algorithm that always
130    // make sure that the i:th element is in its correct
131    // position (N element --> Ordo( N*N ))
132    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
133      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
134        if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
135          std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
136          eigenvectors_.swap_rows(i,j);
137        }
138  }
139  */
140
141
142  utility::matrix PCA::projection(const utility::matrix& samples ) const
143  {
144    const size_t Ncol = samples.columns();
145    const size_t Nrow = samples.rows();
146    utility::matrix projs( Ncol, Ncol );
147
148    utility::vector temp(samples.rows());
149    for( size_t j = 0; j < Ncol; ++j ) {
150      for (size_t i=0; i<Ncol; ++i )
151        temp(i) = samples(i,j); 
152      utility::vector centered( Nrow );
153      for( size_t i = 0; i < Nrow; ++i ) 
154        centered(i) = temp(i) - meanvalues_(i);
155      utility::vector proj( eigenvectors_ * centered );
156      for( size_t i = 0; i < Ncol; ++i ) 
157        projs(i,j)=proj(i);
158    }
159    return projs;
160  }
161
162
163  /*
164  utility::matrix
165  PCA::projection_transposed(const utility::matrix& samples) const
166  {
167    const size_t Ncol = samples.columns();
168    const size_t Nrow = samples.rows();
169    utility::matrix projs( Nrow, Ncol );
170
171    utility::vector temp(samples.rows());
172    for( size_t j = 0; j < Ncol; ++j ) {
173      for (size_t i=0; i<Ncol; ++i )
174        temp(i) = samples(i,j);
175      utility::vector centered( Nrow );
176      for( size_t i = 0; i < Nrow; ++i )
177        centered(i)=temp(i)-meanvalues_(i);
178      utility::vector proj( eigenvectors_ * centered );
179      for( size_t i = 0; i < Nrow; ++i )
180        projs(i,j)=proj(i);
181    }
182    return projs;
183  }
184  */
185
186
187  // This function will row-center the matrix A_,
188  // that is, A_ = A_ - M, where M is a matrix
189  // with the meanvalues of each row
190  void PCA::row_center(utility::matrix& A_center)
191  {
192    meanvalues_ = vector(A_.rows());
193    utility::vector A_row_sum(A_.rows());
194    for (size_t i=0; i<A_row_sum.size(); ++i){
195      const vectorView tmp(A_.row_vec(i));
196      A_row_sum(i) = sum(tmp);
197    }
198    for( size_t i = 0; i < A_center.rows(); ++i ) {
199      meanvalues_(i) = A_row_sum(i) / static_cast<double>(A_.columns());
200      for( size_t j = 0; j < A_center.columns(); ++j )
201        A_center(i,j) = A_(i,j) - meanvalues_(i);
202    }
203  }
204
205}}} // of namespace utility, yat, and theplu
Note: See TracBrowser for help on using the repository browser.