source: trunk/lib/utility/PCA.cc @ 420

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

Merged better_matrix_class changes r402:419 into the trunk.

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