1 /*
2 * VMMLib - Tensor Classes
3 *
4 * @author Rafael Ballester
5 *
6 */
7
8 #ifndef __VMML__T4_HOSVD__HPP__
9 #define __VMML__T4_HOSVD__HPP__
10
11 #include <vmmlib/tensor4.hpp>
12 #include <vmmlib/lapack_svd.hpp>
13 #include <vmmlib/lapack_sym_eigs.hpp>
14 #include <vmmlib/blas_dgemm.hpp>
15 #include <vmmlib/blas_daxpy.hpp>
16
17 namespace vmml
18 {
19
20 template< size_t R1, size_t R2, size_t R3, size_t R4, size_t I1, size_t I2, size_t I3, size_t I4, typename T = float >
21 class t4_hosvd
22 {
23 public:
24
25 typedef tensor4< I1, I2, I3, I4, T > t4_type;
26
27 typedef matrix< I1, R1, T > u1_type;
28 typedef matrix< I2, R2, T > u2_type;
29 typedef matrix< I3, R3, T > u3_type;
30 typedef matrix< I4, R4, T > u4_type;
31
32 typedef matrix< I1, I2*I3*I4, T > u1_unfolded_type;
33 typedef matrix< I2, I1*I3*I4, T > u2_unfolded_type;
34 typedef matrix< I3, I1*I2*I4, T > u3_unfolded_type;
35 typedef matrix< I4, I1*I2*I3, T > u4_unfolded_type;
36
37 typedef matrix< I1, I1, T > u1_cov_type;
38 typedef matrix< I2, I2, T > u2_cov_type;
39 typedef matrix< I3, I3, T > u3_cov_type;
40 typedef matrix< I4, I4, T > u4_cov_type;
41
42 static void apply_mode1(const t4_type& data_, u1_type& u1_);
43 static void apply_mode2(const t4_type& data_, u2_type& u2_);
44 static void apply_mode3(const t4_type& data_, u3_type& u3_);
45 static void apply_mode4(const t4_type& data_, u4_type& u4_);
46
47 protected:
48
49 //hosvd on eigenvalue decomposition = hoeigs
50 template< size_t N, size_t R >
51 static void get_eigs_u_red( const matrix< N, N, T >& data_, matrix< N, R, T >& u_ );
52
53 static void eigs_mode1( const t4_type& data_, u1_type& u1_ );
54 static void eigs_mode2( const t4_type& data_, u2_type& u2_ );
55 static void eigs_mode3( const t4_type& data_, u3_type& u3_ );
56 static void eigs_mode4( const t4_type& data_, u4_type& u4_ );
57
58 }; //end hosvd class
59
60 #define VMML_TEMPLATE_STRING template< size_t R1, size_t R2, size_t R3, size_t R4, size_t I1, size_t I2, size_t I3, size_t I4, typename T >
61 #define VMML_TEMPLATE_CLASSNAME t4_hosvd< R1, R2, R3, R4, I1, I2, I3, I4, T >
62
63 VMML_TEMPLATE_STRING
64 void
apply_mode1(const t4_type & data_,u1_type & u1_)65 VMML_TEMPLATE_CLASSNAME::apply_mode1( const t4_type& data_, u1_type& u1_)
66 {
67 eigs_mode1( data_, u1_ );
68 }
69
70 VMML_TEMPLATE_STRING
71 void
apply_mode2(const t4_type & data_,u2_type & u2_)72 VMML_TEMPLATE_CLASSNAME::apply_mode2( const t4_type& data_, u2_type& u2_)
73 {
74 eigs_mode2( data_, u2_ );
75 }
76
77 VMML_TEMPLATE_STRING
78 void
apply_mode3(const t4_type & data_,u3_type & u3_)79 VMML_TEMPLATE_CLASSNAME::apply_mode3( const t4_type& data_, u3_type& u3_)
80 {
81 eigs_mode3( data_, u3_ );
82 }
83
84 VMML_TEMPLATE_STRING
85 void
apply_mode4(const t4_type & data_,u4_type & u4_)86 VMML_TEMPLATE_CLASSNAME::apply_mode4( const t4_type& data_, u4_type& u4_)
87 {
88 eigs_mode4( data_, u4_ );
89 }
90
91 /* EIGS for mode 1, 2, 3 and 4*/
92
93 VMML_TEMPLATE_STRING
94 void
eigs_mode1(const t4_type & data_,u1_type & u1_)95 VMML_TEMPLATE_CLASSNAME::eigs_mode1( const t4_type& data_, u1_type& u1_ )
96 {
97 //unfolding / matricization
98 u1_unfolded_type* unfolding = new u1_unfolded_type;
99 data_.mode1_unfolding_fwd( *unfolding );
100
101 //covariance matrix of unfolded data
102 u1_cov_type* cov = new u1_cov_type;
103 blas_dgemm< I1, I2*I3*I4, I1, T>* blas_cov = new blas_dgemm< I1, I2*I3*I4, I1, T>;
104 blas_cov->compute( *unfolding, *cov );
105 delete blas_cov;
106 delete unfolding;
107
108 //compute x largest magnitude eigenvalues; x = R
109 get_eigs_u_red( *cov, u1_ );
110
111 delete cov;
112 }
113
114 VMML_TEMPLATE_STRING
115 void
eigs_mode2(const t4_type & data_,u2_type & u2_)116 VMML_TEMPLATE_CLASSNAME::eigs_mode2( const t4_type& data_, u2_type& u2_ )
117 {
118 //unfolding / matricization
119 u2_unfolded_type* unfolding = new u2_unfolded_type;
120 data_.mode2_unfolding_fwd( *unfolding );
121
122 //covariance matrix of unfolded data
123 u2_cov_type* cov = new u2_cov_type;
124 blas_dgemm< I2, I1*I3*I4, I2, T>* blas_cov = new blas_dgemm< I2, I1*I3*I4, I2, T>;
125 blas_cov->compute( *unfolding, *cov );
126 delete blas_cov;
127 delete unfolding;
128
129 //compute x largest magnitude eigenvalues; x = R
130 get_eigs_u_red( *cov, u2_ );
131
132 delete cov;
133 }
134
135 VMML_TEMPLATE_STRING
136 void
eigs_mode3(const t4_type & data_,u3_type & u3_)137 VMML_TEMPLATE_CLASSNAME::eigs_mode3( const t4_type& data_, u3_type& u3_)
138 {
139 //unfolding / matricization
140 u3_unfolded_type* unfolding = new u3_unfolded_type;
141 data_.mode3_unfolding_fwd( *unfolding );
142
143 //covariance matrix of unfolded data
144 u3_cov_type* cov = new u3_cov_type;
145 blas_dgemm< I3, I1*I2*I4, I3, T>* blas_cov = new blas_dgemm< I3, I1*I2*I4, I3, T>;
146 blas_cov->compute( *unfolding, *cov );
147 delete blas_cov;
148 delete unfolding;
149
150 //compute x largest magnitude eigenvalues; x = R
151 get_eigs_u_red( *cov, u3_ );
152
153 delete cov;
154 }
155
156 VMML_TEMPLATE_STRING
157 void
eigs_mode4(const t4_type & data_,u4_type & u4_)158 VMML_TEMPLATE_CLASSNAME::eigs_mode4( const t4_type& data_, u4_type& u4_)
159 {
160 //unfolding / matricization
161 u4_unfolded_type* unfolding = new u4_unfolded_type;
162 data_.mode4_unfolding_fwd( *unfolding );
163
164 //covariance matrix of unfolded data
165 u4_cov_type* cov = new u4_cov_type;
166 blas_dgemm< I4, I1*I2*I3, I4, T>* blas_cov = new blas_dgemm< I4, I1*I2*I3, I4, T>;
167 blas_cov->compute( *unfolding, *cov );
168 delete blas_cov;
169 delete unfolding;
170
171 //compute x largest magnitude eigenvalues; x = R
172 get_eigs_u_red( *cov, u4_ );
173
174 delete cov;
175 }
176
177 /* helper methods for SVD and EIGS*/
178
179 VMML_TEMPLATE_STRING
180 template< size_t N, size_t R >
181 void
get_eigs_u_red(const matrix<N,N,T> & data_,matrix<N,R,T> & u_)182 VMML_TEMPLATE_CLASSNAME::get_eigs_u_red( const matrix< N, N, T >& data_, matrix< N, R, T >& u_ )
183 {
184 typedef matrix< N, N, T > cov_matrix_type;
185 typedef vector< R, T > eigval_type;
186 typedef matrix< N, R, T > eigvec_type;
187 //typedef matrix< N, R, T_coeff > coeff_type;
188
189 //compute x largest magnitude eigenvalues; x = R
190 eigval_type* eigxvalues = new eigval_type;
191 eigvec_type* eigxvectors = new eigvec_type;
192
193 lapack_sym_eigs< N, T > eigs;
194 cov_matrix_type* data = new cov_matrix_type;
195 data->cast_from( data_ );
196 if( !eigs.compute_x( *data, *eigxvectors, *eigxvalues) ) {
197
198 /*if( _is_quantify_coeff ){
199 coeff_type* evec_quant = new coeff_type;
200 T min_value = 0; T max_value = 0;
201 u_.cast_from( *eigxvectors );
202 u_.quantize( *evec_quant, min_value, max_value );
203 evec_quant->dequantize( u_, min_value, max_value );
204 delete evec_quant;
205 } else */
206 std::cerr << "Warning: lapack eigenvector computation returned error code" << std::endl;
207 }
208 u_ = *eigxvectors;
209
210 delete eigxvalues;
211 delete eigxvectors;
212 delete data;
213 }
214
215 #undef VMML_TEMPLATE_STRING
216 #undef VMML_TEMPLATE_CLASSNAME
217
218 }//end vmml namespace
219
220 #endif
221