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

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

fixes #300 - I simply removed the sorting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.6 KB
Line 
1// $Id: PCA.cc 1476 2008-09-04 13:51:35Z 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, 2008 Jari Häkkinen, Peter Johansson
9
10  This file is part of the yat library, http://dev.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_=U;
75
76    eigenvectors_ .transpose();
77    eigenvalues_ = pSVD->s();
78
79    // T
80    for( size_t i = 0; i < eigenvalues_.size(); ++i )
81      eigenvalues_(i) = eigenvalues_(i)*eigenvalues_(i);
82    eigenvalues_ *= 1.0/A_center.rows();
83
84    // Sort the eigenvectors in order of eigenvalues
85    // Simple (not efficient) algorithm that always
86    // make sure that the i:th element is in its correct
87    // position (N element --> Ordo( N*N ))
88    /*
89      // should not be needed since SVD gives single values ordered
90    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
91      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
92        if ( eigenvalues_(j) > eigenvalues_(i) ) {
93          std::swap( eigenvalues_(i), eigenvalues_(j) );
94          eigenvectors_.swap_rows(i,j);
95        }
96    */
97  }
98
99
100  /*
101  void PCA::process_transposed_problem(void)
102  {
103    // Row-center the data matrix
104    utility::matrix A_center( A_.rows(), A_.columns() );
105    this->row_center( A_center );
106
107    // Transform into SVD friendly dimension
108    A_.transpose();
109    A_center.transpose();
110
111    // Single value decompose the data matrix
112    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
113    pSVD->decompose();
114    utility::matrix U(pSVD->U());
115    utility::matrix V(pSVD->V());
116
117    // Read the eigenvectors and eigenvalues
118    eigenvectors_.clone(V);
119    eigenvectors_.transpose();
120    eigenvalues_.clone(pSVD->s());
121
122    // Transform back when done with SVD!
123    // (used V insted of U now for eigenvectors)
124    A_.transpose();
125    A_center.transpose();
126
127    // T
128    for( size_t i = 0; i < eigenvalues_.size(); ++i )
129      eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
130    eigenvalues_ *= 1.0/A_center.rows();
131
132    // Sort the eigenvectors in order of eigenvalues
133    // Simple (not efficient) algorithm that always
134    // make sure that the i:th element is in its correct
135    // position (N element --> Ordo( N*N ))
136    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
137      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
138        if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
139          std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
140          eigenvectors_.swap_rows(i,j);
141        }
142  }
143  */
144
145
146  utility::Matrix PCA::projection(const utility::Matrix& samples ) const
147  {
148    const size_t Ncol = samples.columns();
149    const size_t Nrow = samples.rows();
150    utility::Matrix projs( Ncol, Ncol );
151
152    utility::Vector temp(samples.rows());
153    for( size_t j = 0; j < Ncol; ++j ) {
154      for (size_t i=0; i<Ncol; ++i )
155        temp(i) = samples(i,j); 
156      utility::Vector centered( Nrow );
157      for( size_t i = 0; i < Nrow; ++i ) 
158        centered(i) = temp(i) - meanvalues_(i);
159      utility::Vector proj( eigenvectors_ * centered );
160      for( size_t i = 0; i < Ncol; ++i ) 
161        projs(i,j)=proj(i);
162    }
163    return projs;
164  }
165
166
167  /*
168  utility::matrix
169  PCA::projection_transposed(const utility::matrix& samples) const
170  {
171    const size_t Ncol = samples.columns();
172    const size_t Nrow = samples.rows();
173    utility::matrix projs( Nrow, Ncol );
174
175    utility::vector temp(samples.rows());
176    for( size_t j = 0; j < Ncol; ++j ) {
177      for (size_t i=0; i<Ncol; ++i )
178        temp(i) = samples(i,j);
179      utility::vector centered( Nrow );
180      for( size_t i = 0; i < Nrow; ++i )
181        centered(i)=temp(i)-meanvalues_(i);
182      utility::vector proj( eigenvectors_ * centered );
183      for( size_t i = 0; i < Nrow; ++i )
184        projs(i,j)=proj(i);
185    }
186    return projs;
187  }
188  */
189
190
191  // This function will row-center the matrix A_,
192  // that is, A_ = A_ - M, where M is a matrix
193  // with the meanvalues of each row
194  void PCA::row_center(utility::Matrix& A_center)
195  {
196    meanvalues_ = Vector(A_.rows());
197    utility::Vector A_row_sum(A_.rows());
198    for (size_t i=0; i<A_row_sum.size(); ++i){
199      A_row_sum(i) = sum(A_.row_const_view(i));
200    }
201    for( size_t i = 0; i < A_center.rows(); ++i ) {
202      meanvalues_(i) = A_row_sum(i) / static_cast<double>(A_.columns());
203      for( size_t j = 0; j < A_center.columns(); ++j )
204        A_center(i,j) = A_(i,j) - meanvalues_(i);
205    }
206  }
207
208}}} // of namespace utility, yat, and theplu
Note: See TracBrowser for help on using the repository browser.