# source:trunk/src/PCA.cc@189

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

• Property svn:eol-style set to native
• Property svn:keywords set to Author Date Id Revision
File size: 7.9 KB
Line
1// \$Id: PCA.cc 189 2004-10-21 13:00:06Z jari \$
2
3#include <iostream>
4#include <cmath>
5#include <memory>
6
7#include "PCA.h"
8#include "SVD.h"
9
10
11// using namespace std;
12// using namespace thep_cpp_tools;
13
14namespace theplu {
15namespace cpptools {
16
17PCA::PCA() : process_( false ), explained_calc_( false )
18{
19}
20
21PCA::PCA( const gslapi::matrix& A ) : A_( A ),
22              process_( false ),
23              explained_calc_( false )
24{
25}
26
27void PCA::process()
28{
29  process_ = true;
30  // Row-center the data matrix
31  gslapi::matrix A_center( A_.rows(), A_.columns() );
32  this->row_center( A_center );
33
34
35  // Single value decompose the data matrix
36  std::auto_ptr<cpptools::SVD> pSVD( new cpptools::SVD( A_center ) );
37  pSVD->process();
38  gslapi::matrix U = pSVD->get_u();
39  gslapi::matrix V = pSVD->get_v();
40
41  // Read the eigenvectors and eigenvalues
42  eigenvectors_ = U.transpose();
43  eigenvalues_ = pSVD->get_s_vec();
44
45  // T
46  for( size_t i = 0; i < eigenvalues_.size(); ++i )
47    {
48     eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
49    }
50
51  eigenvalues_ *= 1.0/A_center.rows();
52
53  // Sort the eigenvectors in order of eigenvalues
54  // Simple (not efficient) algorithm that always
55  // make sure that the i:th element is in its correct
56  // position (N element --> Ordo( N*N ))
57  for( size_t i = 0; i < eigenvalues_.size(); ++i )
58    {
59     for( size_t j = i + 1; j < eigenvalues_.size(); ++j )
60       {
61  if( eigenvalues_[ j ] > eigenvalues_[ i ] )
62    {
63      std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
64     eigenvectors_.row_swap(i,j);
65    }
66       }
67    }
68}
69
70
71void PCA::process_transposed_problem()
72{
73  process_ = true;
74  // Row-center the data matrix
75  gslapi::matrix A_center( A_.rows(), A_.columns() );
76  this->row_center( A_center );
77
78  // Transform into SVD friendly dimension
79  ////////////////////////////////////////
80  A_ = A_.transpose();
81  A_center = A_center.transpose();
82
83  // Single value decompose the data matrix
84  std::auto_ptr<cpptools::SVD> pSVD( new cpptools::SVD( A_center ) );
85  pSVD->process();
86  gslapi::matrix U = pSVD->get_u();
87  gslapi::matrix V = pSVD->get_v();
88
89  // Read the eigenvectors and eigenvalues
90  eigenvectors_ = V.transpose();//U.transpose();
91  eigenvalues_ = pSVD->get_s_vec();
92
93  // Transform back when done with SVD!
94  // (used V insted of U now for eigenvectors)
95  ////////////////////////////////////////////
96  A_ = A_.transpose();
97  A_center = A_center.transpose();
98
99
100  // T
101  for( size_t i = 0; i < eigenvalues_.size(); ++i )
102    {
103     eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
104    }
105
106  eigenvalues_ *= 1.0/A_center.rows();
107
108  // Sort the eigenvectors in order of eigenvalues
109  // Simple (not efficient) algorithm that always
110  // make sure that the i:th element is in its correct
111  // position (N element --> Ordo( N*N ))
112  for( size_t i = 0; i < eigenvalues_.size(); ++i )
113    {
114     for( size_t j = i + 1; j < eigenvalues_.size(); ++j )
115       {
116  if( eigenvalues_[ j ] > eigenvalues_[ i ] )
117    {
118     std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
119     eigenvectors_.row_swap(i,j);
120    }
121       }
122    }
123}
124
125
126// This function will row-center the matrix A_,
127// that is, A_ = A_ - M, where M is a matrix
128// with the meanvalues of each row
129void PCA::row_center( gslapi::matrix& A_center )
130{
131  meanvalues_ = gslapi::vector( A_.rows() );
132  gslapi::vector A_row_sum(A_.rows());
133  for (size_t i=0; i<A_row_sum.size(); ++i)
134    A_row_sum(i)=A_[i].sum();
135  for( size_t i = 0; i < A_center.rows(); ++i )
136    {
137     meanvalues_[i] = A_row_sum(i) / static_cast<double>(A_.columns());
138     for( size_t j = 0; j < A_center.columns(); ++j )
139       A_center(i,j) = A_(i,j) - meanvalues_(i);
140    }
141}
142
143
144gslapi::matrix PCA::projection( const gslapi::matrix&
145              samples ) const
146{
147  const size_t Ncol = samples.columns();
148  const size_t Nrow = samples.rows();
149  gslapi::matrix projs( Ncol, Ncol );
150
151  gslapi::vector temp(samples.rows());
152  for( size_t j = 0; j < Ncol; ++j )
153    {
154      for (size_t i=0; i<Ncol; ++i )
155        temp(i) = samples(i,j);
156     gslapi::vector centered( Nrow );
157
158     for( size_t i = 0; i < Nrow; ++i )
159       centered(i) = temp(i) - meanvalues_(i);
160
161     gslapi::vector proj( eigenvectors_ * centered );
162      for( size_t i = 0; i < Ncol; ++i )
163        projs(i,j)=proj(i);
164    }
165  return projs;
166}
167
168
169gslapi::matrix PCA::projection_transposed( const gslapi::matrix&
170             samples ) const
171{
172  const size_t Ncol = samples.columns();
173  const size_t Nrow = samples.rows();
174  gslapi::matrix projs( Nrow, Ncol );
175
176  gslapi::vector temp(samples.rows());
177  for( size_t j = 0; j < Ncol; ++j )
178    {
179      for (size_t i=0; i<Ncol; ++i )
180        temp(i) = samples(i,j);
181     gslapi::vector centered( Nrow );
182
183     for( size_t i = 0; i < Nrow; ++i )
184       centered(i)=temp(i)-meanvalues_(i);
185
186     gslapi::vector proj( eigenvectors_ * centered );
187     for( size_t i = 0; i < Nrow; ++i )
188       projs(i,j)=proj(i);
189    }
190  return projs;
191}
192
193
194
195void PCA::calculate_explained_intensity()
196{
197  size_t DIM = eigenvalues_.size();
198  explained_intensity_ = gslapi::vector( DIM );
199  double sum = 0;
200  for( size_t i = 0; i < DIM; ++i )
201    {
202     sum += eigenvalues_[ i ];
203    }
204  double exp_sum = 0;
205  for( size_t i = 0; i < DIM; ++i )
206    {
207     exp_sum += eigenvalues_[ i ];
208     explained_intensity_[ i ] = exp_sum / sum ;
209    }
210}
211
212double PCA::get_explained_intensity( const size_t& k )
213{
214  if( !explained_calc_ )
215    {
216     this->calculate_explained_intensity();
217     explained_calc_ = true;
218    }
219  return explained_intensity_[ k ];
220}
221
222
223bool PCA::test()
224{
225  std::cout << "1. Print original matrix A = \n";
226  gslapi::matrix A( 3, 4 );
227  for( size_t i = 0; i < 3; ++i )
228    for( size_t j = 0; j < 4; ++j )
229      A(i,j)= sin( static_cast<double>(i+j+2+(i+1)*(j+1)) );
230  A_=A;
231
232// Print only matrix that is row-centered as PCA ( process() )
233// will do that before calculation
234  this->row_center( A );
235  std::cout << A << std::endl;
236
237  std::cout << "2. Testing PCA\n";
238  this->process_transposed_problem();
239
240  for( size_t i = 0; i < eigenvalues_.size(); ++i )
241    {
242     std::cout << "lambda(" << i << ") = " << eigenvalues_[i] << " "
243    << "has eigenvector = " << std::endl;
244     std::cout << eigenvectors_[i] << std::endl;
245    }
246
247  // The eigenvalues and eigenvectors are
248  // from matlab using eig(A'A)/3
249//   double answer_lambda[ 3 ] = { 1.8908, 0.0327, 0 };
250//   double answer_eigenvec[ 3 ][ 3 ] =
251//     {
252//       { 0.4763, -0.6738, 0.5649  },
253//       { 0.8770, 0.3182,  -0.3599 },
254//       { 0.0627, 0.6669 , 0.7425  },
255//     };
256
257//   double ERROR_TOL = 0.0001;
258
259//   // error eigenvalues
260//   cout << "abs( lambda(i) - lambda_matlab(i) ) ) = \n( ";
261//   double max_err = 0;
262//   for( size_t i = 0; i < 3; ++i )
263//     {
264//      double err = fabs( eigenvalues_[ i ] - answer_lambda[ i ] );
265//      if( err > max_err ) max_err = err;
266//      cout << err << " ";
267//     }
268//   cout << ")\n";
269
270//   cerr << "max_error = " << max_err << ", to be compared with ERROR_TOL = "
271//        << ERROR_TOL << endl;
272//   if( max_err < ERROR_TOL )
273//     {
274//      cout << "ok, compared with precision in given matlab-references!\n";
275//     }
276//   else
277//     {
278//      cerr << "Not ok, error to big!\n";
279//      return false;
280//     }
281
282//   // error eigenvectors
283//   max_err = 0;
284//   for( size_t i = 0; i < 3; ++i )
285//     {
286//      cout << "err_eigvec( " << i << " ) = \n";
287//      for( size_t j = 0; j < 3; ++j )
288//        {
289//  double err = fabs( fabs( eigenvectors_.get( i, j ) ) - fabs( answer_eigenvec[ i ][ j ] ) );
290//  if( err > max_err ) max_err = err;
291//  cout << err << endl;
292//        }
293//      cout << ")\n\n";
294//     }
295
296//   cerr << "max_error = " << max_err << ", to be compared with ERROR_TOL = "
297//        << ERROR_TOL << endl;
298//   if( max_err < ERROR_TOL )
299//     {
300//      cout << "ok, compared with precision in given matlab-references!\n";
301//     }
302//   else
303//     {
304//      cerr << "Not ok, error to big!\n";
305//      return false;
306//     }
307
308  std::cout << "\nThe following projection is obtained:\n";
309  std::cout << this->projection_transposed( A_ ).transpose() << std::endl;
310
311
312  return true;
313}
314
315}} // of namespace cpptools and namespace theplu
Note: See TracBrowser for help on using the repository browser.