Changeset 15


Ignore:
Timestamp:
Jun 25, 2003, 12:54:05 PM (18 years ago)
Author:
daniel
Message:

Added the transposed problem to PCA!

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/PCA.cc

    r14 r15  
    5959
    6060
     61void PCA::process_transposed_problem()
     62{
     63  process_ = true;
     64  // Row-center the data matrix
     65  thep_gsl_api::matrix A_center( A_.rows(), A_.cols() );
     66  this->row_center( A_center );
     67
     68  // Transform into SVD friendly dimension
     69  ////////////////////////////////////////
     70  A_ = A_.transpose();
     71  A_center = A_center.transpose();
     72
     73  // Single value decompose the data matrix
     74  auto_ptr<thep_cpp_tools::SVD> pSVD( new thep_cpp_tools::SVD( A_center ) );
     75  pSVD->process();
     76  thep_gsl_api::matrix U = pSVD->get_u();
     77  thep_gsl_api::matrix V = pSVD->get_v();
     78
     79  // Read the eigenvectors and eigenvalues
     80  eigenvectors_ = V.transpose();//U.transpose();
     81  eigenvalues_ = pSVD->get_s_vec();
     82
     83  // Transform back when done with SVD!
     84  // (used V insted of U now for eigenvectors)
     85  ////////////////////////////////////////////
     86  A_ = A_.transpose();
     87  A_center = A_center.transpose();
     88
     89
     90  // T
     91  for( size_t i = 0; i < eigenvalues_.size(); ++i )
     92    {
     93     eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
     94    }
     95
     96  eigenvalues_ /= A_center.rows(); 
     97 
     98  // Sort the eigenvectors in order of eigenvalues
     99  // Simple (not efficient) algorithm that always
     100  // make sure that the i:th element is in its correct
     101  // position (N element --> Ordo( N*N ))
     102  for( size_t i = 0; i < eigenvalues_.size(); ++i )
     103    {
     104     for( size_t j = i + 1; j < eigenvalues_.size(); ++j )
     105       {
     106  if( eigenvalues_[ j ] > eigenvalues_[ i ] )
     107    {
     108     swap( eigenvalues_[ i ], eigenvalues_[ j ] );
     109     eigenvectors_.swap_rows( i, j );
     110    }
     111       }
     112    }
     113}
     114
     115
    61116// This function will row-center the matrix A_,
    62117// that is, A_ = A_ - M, where M is a matrix
     
    106161
    107162
     163thep_gsl_api::matrix PCA::projection_transposed( const thep_gsl_api::matrix&
     164             samples ) const
     165{
     166  assert( process_ );
     167  const size_t Ncol = samples.cols();
     168  const size_t Nrow = samples.rows();
     169  thep_gsl_api::matrix projs( Nrow, Ncol );
     170 
     171  for( size_t j = 0; j < Ncol; ++j )
     172    {
     173     thep_gsl_api::matrix temp = samples.col( j );   
     174     thep_gsl_api::vector centered( Nrow );
     175
     176     for( size_t i = 0; i < Nrow; ++i )
     177       {
     178  centered[ i ] = temp.get( i, 0 ) - meanvalues_[ i ];
     179  //cout << temp.get( i, 0 ) << " " << meanvalues_[ i ] << " " << centered[ i ] << endl;
     180       }
     181
     182     thep_gsl_api::vector proj( eigenvectors_ * centered );
     183     for( size_t i = 0; i < Nrow; ++i )
     184       {
     185    projs.set( i, j,  proj[ i ] );
     186       }
     187    }
     188  return projs;
     189}
     190
     191
     192
    108193void PCA::calculate_explained_intensity()
    109194{
     
    137222{
    138223  cout << "1. Print original matrix A = \n";//Testing row-centering of matrix\n";
    139   thep_gsl_api::matrix A( 4, 3 );
    140   for( size_t i = 0; i < 4; ++i )
    141     {
    142      for( size_t j = 0; j < 3; ++j )
     224  thep_gsl_api::matrix A( 3, 4 );
     225  for( size_t i = 0; i < 3; ++i )
     226    {
     227     for( size_t j = 0; j < 4; ++j )
    143228       {
    144229  A.set( i, j, sin( static_cast<double>( i + j + 2 + (i+1)*(j+1) ) ) );
     
    153238
    154239  cout << "2. Testing PCA\n";
    155   this->process();
     240  this->process_transposed_problem();
    156241
    157242  for( size_t i = 0; i < eigenvalues_.size(); ++i )
     
    224309   
    225310  cout << "\nThe following projection is obtained:\n";
    226   cout << this->projection( A_ ) << endl;
     311  cout << this->projection_transposed( A_ ).transpose() << endl;
    227312 
    228313
  • trunk/src/PCA.h

    r14 r15  
    2525     this space, all Nx1 vectors will have dimension Mx1. Hence
    2626     the projection will have dimension MxM where each column is
    27      a point in the new space. Also, it assumes that M>N. The
    28      opposite problem is yet not implemented. However, could easily
    29      be done by using that, if A=input-data=MxN, then SVM(A)->A=USV,
    30      and SVM(A')->A'=VSU. So by changing eigenvectors=U to
    31      eigenvectors=V and rewrite the projection-method in terms of columns
    32      instead of rows ...
     27     a point in the new space. Also, it assumes that M>N. The opposite
     28     problem is added in the functions: process_transposed_problem and
     29     projection_transposed()...
    3330  */
    3431 
     
    5855    void process();
    5956
     57    /**
     58       If M<N
     59     */
     60    void process_transposed_problem();
    6061   
    6162    /**
     
    100101    */
    101102    thep_gsl_api::matrix projection( const thep_gsl_api::matrix& ) const;
     103    thep_gsl_api::matrix projection_transposed( const thep_gsl_api::matrix& ) const;
    102104
     105   
    103106   
    104107  private:
Note: See TracChangeset for help on using the changeset viewer.