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