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

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

Fixes #299. Memory leak in matrix was found and removed.

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