source: trunk/lib/utility/PCA.h @ 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: 3.5 KB
Line 
1// $Id: PCA.h 420 2005-12-02 00:15:50Z jari $
2
3#ifndef _theplu_utility_pca_
4#define _theplu_utility_pca_
5
6#include <c++_tools/gslapi/matrix.h>
7#include <c++_tools/gslapi/vector.h>
8
9
10namespace theplu {
11namespace utility {
12
13  /**
14     Class performing PCA using SVD. This class assumes that
15     the columns corresponds to the dimenension of the problem.
16     That means if data has dimension NxM (M=columns) the number
17     of principal-axes will equal M-1. When projecting data into
18     this space, all Nx1 vectors will have dimension Mx1. Hence
19     the projection will have dimension MxM where each column is
20     a point in the new space. Also, it assumes that M>N. The opposite
21     problem is added in the functions: process_transposed_problem and
22     projection_transposed()...
23  */
24  class PCA
25  {
26  public:
27    /**
28       Default constructor (not implemented)
29    */
30    PCA(void); 
31
32    /**
33       Constructor taking the data-matrix as input. No row-centering
34       should have been performed and no products.
35     */
36    inline explicit PCA(const gslapi::matrix& A)
37      : A_(A), process_(false), explained_calc_(false) {}
38
39    /**
40       Will perform PCA according to the following scheme: \n
41       1: Rowcenter A  \n
42       2: SVD(A)  --> USV' \n
43       3: Calculate eigenvalues according to \n
44          \f$ \lambda_{ii} = s_{ii}/N_{rows} \f$ \n
45       4: Sort eigenvectors (from matrix V) according to descending eigenvalues\n
46    */
47    void process(void);
48
49    /**
50       If M<N use this method instead. Using the same format as before
51       where rows in the matrix corresponds to the dimensional coordinate.
52       The only difference is in the SVD step where the matrix V is used
53       after running the transposed matrix. For projections, see
54       projection_transposed() method.
55     */
56    void process_transposed_problem(void);
57
58    /**
59       @return Eigenvector \a i.
60    */
61    inline gslapi::vector get_eigenvector(const size_t& i) const
62      { return gslapi::vector(eigenvectors_,i); }
63
64    /**
65       Returns eigenvalues to covariance matrix
66       \f$ C = \frac{1}{N^2}A^TA \f$
67    */
68    inline double
69    get_eigenvalue(const size_t& i) const { return eigenvalues_[i]; }
70
71    /**
72       Returns the explained intensity of component \a K \n
73       \f$I = \frac{ \sum^{K}_{i=1} \lambda_i }{ \sum^{N}_{j=1} \lambda_j }\f$ \n
74       where \f$N\f$ is the dimension
75    */
76    double get_explained_intensity( const size_t& k );
77
78    /**
79       This function will project data onto the new coordinate-system
80       where the axes are the calculated eigenvectors. This means that
81       PCA must have been run before this function can be used!
82       Output is presented as coordinates in the N-dimensional room
83       spanned by the eigenvectors.
84    */
85    gslapi::matrix projection( const gslapi::matrix& ) const;
86
87    /**
88       Same as projection() but works when used
89       process_transposed_problem().
90    */
91    gslapi::matrix projection_transposed( const gslapi::matrix& ) const;
92
93
94  private:
95    gslapi::matrix  A_; 
96    gslapi::matrix  eigenvectors_;
97    gslapi::vector  eigenvalues_;
98    gslapi::vector  explained_intensity_;
99    gslapi::vector  meanvalues_;
100    bool process_, explained_calc_;
101   
102    /**
103       Private function that will row-center the matrix A,
104       that is, A = A - M, where M is a matrix
105       with the meanvalues of each row
106    */
107    void row_center( gslapi::matrix& A_center );
108
109    /**
110       Private function that will calculate the explained
111       intensity
112    */
113    void calculate_explained_intensity();
114  }; // class PCA 
115
116
117}} // of namespace utility and namespace theplu
118
119#endif
Note: See TracBrowser for help on using the repository browser.