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

Last change on this file since 789 was 789, checked in by Jari Häkkinen, 16 years ago

Addresses #193. vector now works as outlined here. Added some
functionality. Added a clone function that facilitates resizing of
vectors. clone is needed since assignement operator functionality is
changed.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.0 KB
Line 
1// $Id: PCA.cc 789 2007-03-10 20:07:13Z jari $
2
3/*
4  Copyright (C) The authors contributing to this file.
5
6  This file is part of the yat library, http://lev.thep.lu.se/trac/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 2 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 this program; if not, write to the Free Software
20  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
21  02111-1307, USA.
22*/
23
24#include "PCA.h"
25#include "SVD.h"
26#include "utility.h"
27
28#include <iostream>
29#include <cmath>
30#include <memory>
31
32namespace theplu {
33namespace yat {
34namespace utility {
35
36
37  PCA::PCA(const utility::matrix& A)
38    : A_(A), process_(false), explained_calc_(false)
39  {
40  }
41
42
43  utility::vector PCA::get_eigenvector(size_t i) const
44  {
45    return utility::vector(eigenvectors_,i);
46  }
47
48
49  double PCA::get_eigenvalue(size_t i) const
50  {
51    return eigenvalues_[i];
52  }
53
54
55  void PCA::process(void)
56  {
57    process_ = true;
58    // Row-center the data matrix
59    utility::matrix A_center( A_.rows(), A_.columns() );
60    this->row_center( A_center );
61
62    // Single value decompose the data matrix
63    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
64    pSVD->decompose();
65    utility::matrix U = pSVD->U(); 
66    utility::matrix V = pSVD->V();
67
68    // Read the eigenvectors and eigenvalues
69    eigenvectors_ = U;
70    eigenvectors_ .transpose();
71    eigenvalues_.clone(pSVD->s());
72
73    // T
74    for( size_t i = 0; i < eigenvalues_.size(); ++i )
75      eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
76    eigenvalues_ *= 1.0/A_center.rows();
77
78    // Sort the eigenvectors in order of eigenvalues
79    // Simple (not efficient) algorithm that always
80    // make sure that the i:th element is in its correct
81    // position (N element --> Ordo( N*N ))
82    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
83      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
84        if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
85          std::swap( eigenvalues_[ i ], eigenvalues_[ j ] ); 
86          eigenvectors_.swap_rows(i,j);
87        }
88  }
89
90
91  void PCA::process_transposed_problem(void)
92  {
93    process_ = true;
94    // Row-center the data matrix
95    utility::matrix A_center( A_.rows(), A_.columns() );
96    this->row_center( A_center );
97
98    // Transform into SVD friendly dimension
99    A_.transpose();
100    A_center.transpose();
101
102    // Single value decompose the data matrix
103    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
104    pSVD->decompose();
105    utility::matrix U = pSVD->U(); 
106    utility::matrix V = pSVD->V();
107
108    // Read the eigenvectors and eigenvalues
109    eigenvectors_=V;
110    eigenvectors_.transpose();
111    eigenvalues_.clone(pSVD->s());
112
113    // Transform back when done with SVD!
114    // (used V insted of U now for eigenvectors)
115    A_.transpose();
116    A_center.transpose();
117
118    // T
119    for( size_t i = 0; i < eigenvalues_.size(); ++i )
120      eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
121    eigenvalues_ *= 1.0/A_center.rows();
122
123    // Sort the eigenvectors in order of eigenvalues
124    // Simple (not efficient) algorithm that always
125    // make sure that the i:th element is in its correct
126    // position (N element --> Ordo( N*N ))
127    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
128      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
129        if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
130          std::swap( eigenvalues_[ i ], eigenvalues_[ j ] ); 
131          eigenvectors_.swap_rows(i,j);
132        }
133  }
134
135
136  // This function will row-center the matrix A_,
137  // that is, A_ = A_ - M, where M is a matrix
138  // with the meanvalues of each row
139  void PCA::row_center(utility::matrix& A_center)
140  {
141    meanvalues_.clone(utility::vector(A_.rows()));
142    utility::vector A_row_sum(A_.rows());
143    for (size_t i=0; i<A_row_sum.size(); ++i)
144      A_row_sum(i)=utility::sum(utility::vector(A_,i));
145    for( size_t i = 0; i < A_center.rows(); ++i ) {
146      meanvalues_[i] = A_row_sum(i) / static_cast<double>(A_.columns());
147      for( size_t j = 0; j < A_center.columns(); ++j )
148        A_center(i,j) = A_(i,j) - meanvalues_(i);
149    }
150  }
151
152
153  utility::matrix PCA::projection(const utility::matrix& samples ) const
154  {
155    const size_t Ncol = samples.columns();
156    const size_t Nrow = samples.rows();
157    utility::matrix projs( Ncol, Ncol );
158
159    utility::vector temp(samples.rows());
160    for( size_t j = 0; j < Ncol; ++j ) {
161      for (size_t i=0; i<Ncol; ++i )
162        temp(i) = samples(i,j); 
163      utility::vector centered( Nrow );
164      for( size_t i = 0; i < Nrow; ++i ) 
165        centered(i) = temp(i) - meanvalues_(i);
166      utility::vector proj( eigenvectors_ * centered );
167      for( size_t i = 0; i < Ncol; ++i ) 
168        projs(i,j)=proj(i);
169    }
170    return projs;
171  }
172
173
174  utility::matrix
175  PCA::projection_transposed(const utility::matrix& samples) const
176  {
177    const size_t Ncol = samples.columns();
178    const size_t Nrow = samples.rows();
179    utility::matrix projs( Nrow, Ncol );
180
181    utility::vector temp(samples.rows());
182    for( size_t j = 0; j < Ncol; ++j ) {
183      for (size_t i=0; i<Ncol; ++i )
184        temp(i) = samples(i,j); 
185      utility::vector centered( Nrow );
186      for( size_t i = 0; i < Nrow; ++i ) 
187        centered(i)=temp(i)-meanvalues_(i);
188      utility::vector proj( eigenvectors_ * centered );
189      for( size_t i = 0; i < Nrow; ++i ) 
190        projs(i,j)=proj(i);
191    }
192    return projs;
193  }
194
195
196  void PCA::calculate_explained_intensity(void)
197  {
198    size_t DIM = eigenvalues_.size();
199    explained_intensity_ = utility::vector( DIM );
200    double sum = 0;
201    for( size_t i = 0; i < DIM; ++i ) 
202      sum += eigenvalues_[ i ];
203    double exp_sum = 0;
204    for( size_t i = 0; i < DIM; ++i ) {
205      exp_sum += eigenvalues_[ i ];
206      explained_intensity_[ i ] = exp_sum / sum ;
207    }
208  }
209
210
211  double PCA::get_explained_intensity( const size_t& k )
212  {
213    if( !explained_calc_ ) {
214      this->calculate_explained_intensity();
215      explained_calc_ = true;
216    }
217    return explained_intensity_[ k ];
218  }
219
220
221}}} // of namespace utility, yat, and theplu
Note: See TracBrowser for help on using the repository browser.