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

Last change on this file since 1706 was 1487, checked in by Jari Häkkinen, 13 years ago

Addresses #436. GPL license copy reference should also be updated.

  • 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 1487 2008-09-10 08:41:36Z jari $
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 3 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 yat. If not, see <http://www.gnu.org/licenses/>.
24*/
25
26#include "PCA.h"
27#include "SVD.h"
28#include "utility.h"
29#include "VectorView.h"
30
31#include <iostream>
32#include <cmath>
33#include <memory>
34
35namespace theplu {
36namespace yat {
37namespace utility {
38
39
40  PCA::PCA(const utility::Matrix& A)
41    : A_(A)
42  {
43    process();
44  }
45
46
47  const utility::Vector& PCA::eigenvalues(void) const
48  {
49    return eigenvalues_;
50  }
51
52
53  const utility::Matrix& PCA::eigenvectors(void) const
54  {
55    return eigenvectors_;
56  }
57
58
59  void PCA::process(void)
60  {
61    // Row-center the data matrix
62    utility::Matrix A_center( A_.rows(), A_.columns() );
63    this->row_center( A_center );
64
65    // Single value decompose the data matrix
66    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
67    pSVD->decompose();
68    utility::Matrix U(pSVD->U());
69    utility::Matrix V(pSVD->V());
70
71    // Read the eigenvectors and eigenvalues
72    eigenvectors_=U;
73
74    eigenvectors_ .transpose();
75    eigenvalues_ = pSVD->s();
76
77    // T
78    for( size_t i = 0; i < eigenvalues_.size(); ++i )
79      eigenvalues_(i) = eigenvalues_(i)*eigenvalues_(i);
80    eigenvalues_ *= 1.0/A_center.rows();
81
82    // Sort the eigenvectors in order of eigenvalues
83    // Simple (not efficient) algorithm that always
84    // make sure that the i:th element is in its correct
85    // position (N element --> Ordo( N*N ))
86    /*
87      // should not be needed since SVD gives single values ordered
88    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
89      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
90        if ( eigenvalues_(j) > eigenvalues_(i) ) {
91          std::swap( eigenvalues_(i), eigenvalues_(j) );
92          eigenvectors_.swap_rows(i,j);
93        }
94    */
95  }
96
97
98  /*
99  void PCA::process_transposed_problem(void)
100  {
101    // Row-center the data matrix
102    utility::matrix A_center( A_.rows(), A_.columns() );
103    this->row_center( A_center );
104
105    // Transform into SVD friendly dimension
106    A_.transpose();
107    A_center.transpose();
108
109    // Single value decompose the data matrix
110    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
111    pSVD->decompose();
112    utility::matrix U(pSVD->U());
113    utility::matrix V(pSVD->V());
114
115    // Read the eigenvectors and eigenvalues
116    eigenvectors_.clone(V);
117    eigenvectors_.transpose();
118    eigenvalues_.clone(pSVD->s());
119
120    // Transform back when done with SVD!
121    // (used V insted of U now for eigenvectors)
122    A_.transpose();
123    A_center.transpose();
124
125    // T
126    for( size_t i = 0; i < eigenvalues_.size(); ++i )
127      eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
128    eigenvalues_ *= 1.0/A_center.rows();
129
130    // Sort the eigenvectors in order of eigenvalues
131    // Simple (not efficient) algorithm that always
132    // make sure that the i:th element is in its correct
133    // position (N element --> Ordo( N*N ))
134    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
135      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
136        if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
137          std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
138          eigenvectors_.swap_rows(i,j);
139        }
140  }
141  */
142
143
144  utility::Matrix PCA::projection(const utility::Matrix& samples ) const
145  {
146    const size_t Ncol = samples.columns();
147    const size_t Nrow = samples.rows();
148    utility::Matrix projs( Ncol, Ncol );
149
150    utility::Vector temp(samples.rows());
151    for( size_t j = 0; j < Ncol; ++j ) {
152      for (size_t i=0; i<Ncol; ++i )
153        temp(i) = samples(i,j); 
154      utility::Vector centered( Nrow );
155      for( size_t i = 0; i < Nrow; ++i ) 
156        centered(i) = temp(i) - meanvalues_(i);
157      utility::Vector proj( eigenvectors_ * centered );
158      for( size_t i = 0; i < Ncol; ++i ) 
159        projs(i,j)=proj(i);
160    }
161    return projs;
162  }
163
164
165  /*
166  utility::matrix
167  PCA::projection_transposed(const utility::matrix& samples) const
168  {
169    const size_t Ncol = samples.columns();
170    const size_t Nrow = samples.rows();
171    utility::matrix projs( Nrow, Ncol );
172
173    utility::vector temp(samples.rows());
174    for( size_t j = 0; j < Ncol; ++j ) {
175      for (size_t i=0; i<Ncol; ++i )
176        temp(i) = samples(i,j);
177      utility::vector centered( Nrow );
178      for( size_t i = 0; i < Nrow; ++i )
179        centered(i)=temp(i)-meanvalues_(i);
180      utility::vector proj( eigenvectors_ * centered );
181      for( size_t i = 0; i < Nrow; ++i )
182        projs(i,j)=proj(i);
183    }
184    return projs;
185  }
186  */
187
188
189  // This function will row-center the matrix A_,
190  // that is, A_ = A_ - M, where M is a matrix
191  // with the meanvalues of each row
192  void PCA::row_center(utility::Matrix& A_center)
193  {
194    meanvalues_ = Vector(A_.rows());
195    utility::Vector A_row_sum(A_.rows());
196    for (size_t i=0; i<A_row_sum.size(); ++i){
197      A_row_sum(i) = sum(A_.row_const_view(i));
198    }
199    for( size_t i = 0; i < A_center.rows(); ++i ) {
200      meanvalues_(i) = A_row_sum(i) / static_cast<double>(A_.columns());
201      for( size_t j = 0; j < A_center.columns(); ++j )
202        A_center(i,j) = A_(i,j) - meanvalues_(i);
203    }
204  }
205
206}}} // of namespace utility, yat, and theplu
Note: See TracBrowser for help on using the repository browser.