1 | // $Id: PCA.cc 1476 2008-09-04 13:51:35Z peter $ |
---|
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, 2008 Jari Häkkinen, Peter Johansson |
---|
9 | |
---|
10 | This file is part of the yat library, http://dev.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 | |
---|
37 | namespace theplu { |
---|
38 | namespace yat { |
---|
39 | namespace 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 | /* |
---|
89 | // should not be needed since SVD gives single values ordered |
---|
90 | for ( size_t i = 0; i < eigenvalues_.size(); ++i ) |
---|
91 | for ( size_t j = i + 1; j < eigenvalues_.size(); ++j ) |
---|
92 | if ( eigenvalues_(j) > eigenvalues_(i) ) { |
---|
93 | std::swap( eigenvalues_(i), eigenvalues_(j) ); |
---|
94 | eigenvectors_.swap_rows(i,j); |
---|
95 | } |
---|
96 | */ |
---|
97 | } |
---|
98 | |
---|
99 | |
---|
100 | /* |
---|
101 | void PCA::process_transposed_problem(void) |
---|
102 | { |
---|
103 | // Row-center the data matrix |
---|
104 | utility::matrix A_center( A_.rows(), A_.columns() ); |
---|
105 | this->row_center( A_center ); |
---|
106 | |
---|
107 | // Transform into SVD friendly dimension |
---|
108 | A_.transpose(); |
---|
109 | A_center.transpose(); |
---|
110 | |
---|
111 | // Single value decompose the data matrix |
---|
112 | std::auto_ptr<SVD> pSVD( new SVD( A_center ) ); |
---|
113 | pSVD->decompose(); |
---|
114 | utility::matrix U(pSVD->U()); |
---|
115 | utility::matrix V(pSVD->V()); |
---|
116 | |
---|
117 | // Read the eigenvectors and eigenvalues |
---|
118 | eigenvectors_.clone(V); |
---|
119 | eigenvectors_.transpose(); |
---|
120 | eigenvalues_.clone(pSVD->s()); |
---|
121 | |
---|
122 | // Transform back when done with SVD! |
---|
123 | // (used V insted of U now for eigenvectors) |
---|
124 | A_.transpose(); |
---|
125 | A_center.transpose(); |
---|
126 | |
---|
127 | // T |
---|
128 | for( size_t i = 0; i < eigenvalues_.size(); ++i ) |
---|
129 | eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ]; |
---|
130 | eigenvalues_ *= 1.0/A_center.rows(); |
---|
131 | |
---|
132 | // Sort the eigenvectors in order of eigenvalues |
---|
133 | // Simple (not efficient) algorithm that always |
---|
134 | // make sure that the i:th element is in its correct |
---|
135 | // position (N element --> Ordo( N*N )) |
---|
136 | for ( size_t i = 0; i < eigenvalues_.size(); ++i ) |
---|
137 | for ( size_t j = i + 1; j < eigenvalues_.size(); ++j ) |
---|
138 | if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) { |
---|
139 | std::swap( eigenvalues_[ i ], eigenvalues_[ j ] ); |
---|
140 | eigenvectors_.swap_rows(i,j); |
---|
141 | } |
---|
142 | } |
---|
143 | */ |
---|
144 | |
---|
145 | |
---|
146 | utility::Matrix PCA::projection(const utility::Matrix& samples ) const |
---|
147 | { |
---|
148 | const size_t Ncol = samples.columns(); |
---|
149 | const size_t Nrow = samples.rows(); |
---|
150 | utility::Matrix projs( Ncol, Ncol ); |
---|
151 | |
---|
152 | utility::Vector temp(samples.rows()); |
---|
153 | for( size_t j = 0; j < Ncol; ++j ) { |
---|
154 | for (size_t i=0; i<Ncol; ++i ) |
---|
155 | temp(i) = samples(i,j); |
---|
156 | utility::Vector centered( Nrow ); |
---|
157 | for( size_t i = 0; i < Nrow; ++i ) |
---|
158 | centered(i) = temp(i) - meanvalues_(i); |
---|
159 | utility::Vector proj( eigenvectors_ * centered ); |
---|
160 | for( size_t i = 0; i < Ncol; ++i ) |
---|
161 | projs(i,j)=proj(i); |
---|
162 | } |
---|
163 | return projs; |
---|
164 | } |
---|
165 | |
---|
166 | |
---|
167 | /* |
---|
168 | utility::matrix |
---|
169 | PCA::projection_transposed(const utility::matrix& samples) const |
---|
170 | { |
---|
171 | const size_t Ncol = samples.columns(); |
---|
172 | const size_t Nrow = samples.rows(); |
---|
173 | utility::matrix projs( Nrow, Ncol ); |
---|
174 | |
---|
175 | utility::vector temp(samples.rows()); |
---|
176 | for( size_t j = 0; j < Ncol; ++j ) { |
---|
177 | for (size_t i=0; i<Ncol; ++i ) |
---|
178 | temp(i) = samples(i,j); |
---|
179 | utility::vector centered( Nrow ); |
---|
180 | for( size_t i = 0; i < Nrow; ++i ) |
---|
181 | centered(i)=temp(i)-meanvalues_(i); |
---|
182 | utility::vector proj( eigenvectors_ * centered ); |
---|
183 | for( size_t i = 0; i < Nrow; ++i ) |
---|
184 | projs(i,j)=proj(i); |
---|
185 | } |
---|
186 | return projs; |
---|
187 | } |
---|
188 | */ |
---|
189 | |
---|
190 | |
---|
191 | // This function will row-center the matrix A_, |
---|
192 | // that is, A_ = A_ - M, where M is a matrix |
---|
193 | // with the meanvalues of each row |
---|
194 | void PCA::row_center(utility::Matrix& A_center) |
---|
195 | { |
---|
196 | meanvalues_ = Vector(A_.rows()); |
---|
197 | utility::Vector A_row_sum(A_.rows()); |
---|
198 | for (size_t i=0; i<A_row_sum.size(); ++i){ |
---|
199 | A_row_sum(i) = sum(A_.row_const_view(i)); |
---|
200 | } |
---|
201 | for( size_t i = 0; i < A_center.rows(); ++i ) { |
---|
202 | meanvalues_(i) = A_row_sum(i) / static_cast<double>(A_.columns()); |
---|
203 | for( size_t j = 0; j < A_center.columns(); ++j ) |
---|
204 | A_center(i,j) = A_(i,j) - meanvalues_(i); |
---|
205 | } |
---|
206 | } |
---|
207 | |
---|
208 | }}} // of namespace utility, yat, and theplu |
---|