vmmlib  1.7.0
 All Classes Namespaces Functions Pages
t3_hosvd.hpp
1 /*
2  * Copyright (c) 2006-2014, Visualization and Multimedia Lab,
3  * University of Zurich <http://vmml.ifi.uzh.ch>,
4  * Eyescale Software GmbH,
5  * Blue Brain Project, EPFL
6  *
7  * This file is part of VMMLib <https://github.com/VMML/vmmlib/>
8  *
9  * Redistribution and use in source and binary forms, with or without
10  * modification, are permitted provided that the following conditions are met:
11  *
12  * Redistributions of source code must retain the above copyright notice, this
13  * list of conditions and the following disclaimer. Redistributions in binary
14  * form must reproduce the above copyright notice, this list of conditions and
15  * the following disclaimer in the documentation and/or other materials provided
16  * with the distribution. Neither the name of the Visualization and Multimedia
17  * Lab, University of Zurich nor the names of its contributors may be used to
18  * endorse or promote products derived from this software without specific prior
19  * written permission.
20  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
21  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
22  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
23  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
24  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
25  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
26  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
27  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
28  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
29  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30  * POSSIBILITY OF SUCH DAMAGE.
31  */
32 
33 /* @author Susanne Suter
34  *
35  * The Tucker3 tensor class is consists of the same components (core tensor, basis matrices u1-u3) as the tucker3 model described in:
36  * - Tucker, 1966: Some mathematical notes on three-mode factor analysis, Psychometrika.
37  * - Kroonenberg & De Leeuw, 1980: Principal component analysis of three-mode data by means of alternating least squares algorithms. Psychometrika. (TUCKALS)
38  * - De Lathauwer, De Moor, Vandewalle, 2000a: A multilinear singular value decomposition, SIAM J. Matrix Anal. Appl.
39  * - Kolda & Bader, 2009: Tensor Decompositions and Applications, SIAM Review.
40  * - Bader & Kolda, 2006: Algorithm 862: Matlab tensor classes for fast algorithm prototyping. ACM Transactions on Mathematical Software.
41  *
42  */
43 
44 #ifndef __VMML__T3_HOSVD__HPP__
45 #define __VMML__T3_HOSVD__HPP__
46 
47 #include <vmmlib/tensor3.hpp>
48 #include <vmmlib/lapack_svd.hpp>
49 #include <vmmlib/lapack_sym_eigs.hpp>
50 #include <vmmlib/blas_dgemm.hpp>
51 #include <vmmlib/blas_daxpy.hpp>
52 
53 enum hosvd_method {
54  eigs_e,
55  svd_e
56 };
57 
58 
59 namespace vmml
60 {
61 
62  template< size_t R1, size_t R2, size_t R3, size_t I1, size_t I2, size_t I3, typename T = float >
63  class t3_hosvd
64  {
65  public:
66 
67  typedef double T_svd;
68 
69  typedef tensor3< I1, I2, I3, T > t3_type;
70 
71  typedef matrix< I1, R1, T > u1_type;
72  typedef matrix< I2, R2, T > u2_type;
73  typedef matrix< I3, R3, T > u3_type;
74 
75  typedef matrix< I1, I2*I3, T > u1_unfolded_type;
76  typedef matrix< I2, I1*I3, T > u2_unfolded_type;
77  typedef matrix< I3, I1*I2, T > u3_unfolded_type;
78 
79  typedef matrix< I1, I1, T > u1_cov_type;
80  typedef matrix< I2, I2, T > u2_cov_type;
81  typedef matrix< I3, I3, T > u3_cov_type;
82 
83  /* higher-order singular value decomposition (HOSVD) with full rank decomposition (also known as Tucker decomposition).
84  see: De Lathauer et al, 2000a: A multilinear singular value decomposition.
85  the hosvd can be computed (a) with n-mode PCA, i.e., an eigenvalue decomposition on the covariance matrix of every mode's matricization, and
86  (b) by performing a 2D SVD on the matricization of every mode. Matrix matricization means that a tensor I1xI2xI3 is unfolded/sliced into one matrix
87  with the dimensions I1xI2I3, which corresponds to a matrizitation alonge mode I1.
88  other known names for HOSVD: n-mode SVD, 3-mode factor analysis (3MFA, tucker3), 3M-PCA, n-mode PCA, higher-order SVD
89  */
90 
91  static void apply_mode1( const t3_type& data_, u1_type& u1_, hosvd_method method_ = eigs_e );
92  static void apply_mode2( const t3_type& data_, u2_type& u2_, hosvd_method method_ = eigs_e );
93  static void apply_mode3( const t3_type& data_, u3_type& u3_, hosvd_method method_ = eigs_e );
94 
95  static void apply_all( const t3_type& data_, u1_type& u1_, u2_type& u2_, u3_type& u3_, hosvd_method method_ = eigs_e );
96  static void hosvd( const t3_type& data_, u1_type& u1_, u2_type& u2_, u3_type& u3_ );
97  static void hoeigs( const t3_type& data_, u1_type& u1_, u2_type& u2_, u3_type& u3_ );
98 
99  protected:
100 
101  //hosvd
102  template< size_t M, size_t N, size_t R >
103  static void get_svd_u_red( const matrix< M, N, T >& data_, matrix< M, R, T >& u_ );
104 
105  static void svd_mode1( const t3_type& data_, u1_type& u1_ );
106  static void svd_mode2( const t3_type& data_, u2_type& u2_ );
107  static void svd_mode3( const t3_type& data_, u3_type& u3_ );
108 
109  //hosvd on eigenvalue decomposition = hoeigs
110  template< size_t N, size_t R >
111  static void get_eigs_u_red( const matrix< N, N, T >& data_, matrix< N, R, T >& u_ );
112 
113  static void eigs_mode1( const t3_type& data_, u1_type& u1_ );
114  static void eigs_mode2( const t3_type& data_, u2_type& u2_ );
115  static void eigs_mode3( const t3_type& data_, u3_type& u3_ );
116 
117  }; //end hosvd class
118 
119 #define VMML_TEMPLATE_STRING template< size_t R1, size_t R2, size_t R3, size_t I1, size_t I2, size_t I3, typename T >
120 #define VMML_TEMPLATE_CLASSNAME t3_hosvd< R1, R2, R3, I1, I2, I3, T >
121 
122 
123 
124 VMML_TEMPLATE_STRING
125 void
126 VMML_TEMPLATE_CLASSNAME::hosvd( const t3_type& data_, u1_type& u1_, u2_type& u2_, u3_type& u3_ )
127 {
128  svd_mode1( data_, u1_ );
129  svd_mode2( data_, u2_ );
130  svd_mode3( data_, u3_ );
131 }
132 
133 VMML_TEMPLATE_STRING
134 void
135 VMML_TEMPLATE_CLASSNAME::hoeigs( const t3_type& data_, u1_type& u1_, u2_type& u2_, u3_type& u3_ )
136 {
137  eigs_mode1( data_, u1_ );
138  eigs_mode2( data_, u2_ );
139  eigs_mode3( data_, u3_ );
140 }
141 
142 VMML_TEMPLATE_STRING
143 void
144 VMML_TEMPLATE_CLASSNAME::apply_all( const t3_type& data_, u1_type& u1_, u2_type& u2_, u3_type& u3_, hosvd_method method_ )
145 {
146  apply_mode1( data_, u1_, method_ );
147  apply_mode2( data_, u2_, method_ );
148  apply_mode3( data_, u3_, method_ );
149 }
150 
151 VMML_TEMPLATE_STRING
152 void
153 VMML_TEMPLATE_CLASSNAME::apply_mode1( const t3_type& data_, u1_type& u1_, hosvd_method method_ )
154 {
155  switch ( method_ )
156  {
157  case 0:
158  eigs_mode1( data_, u1_ );
159  break;
160  case 1:
161  svd_mode1( data_, u1_ );
162  break;
163  default:
164  eigs_mode1( data_, u1_ );
165  }
166 }
167 
168 
169 VMML_TEMPLATE_STRING
170 void
171 VMML_TEMPLATE_CLASSNAME::apply_mode2( const t3_type& data_, u2_type& u2_, hosvd_method method_ )
172 {
173  switch ( method_ )
174  {
175  case 0:
176  eigs_mode2( data_, u2_ );
177  break;
178  case 1:
179  svd_mode2( data_, u2_ );
180  break;
181  default:
182  eigs_mode2( data_, u2_ );
183  }
184 }
185 
186 
187 
188 VMML_TEMPLATE_STRING
189 void
190 VMML_TEMPLATE_CLASSNAME::apply_mode3( const t3_type& data_, u3_type& u3_, hosvd_method method_ )
191 {
192  switch ( method_ )
193  {
194  case 0:
195  eigs_mode3( data_, u3_ );
196  break;
197  case 1:
198  svd_mode3( data_, u3_ );
199  break;
200  default:
201  eigs_mode3( data_, u3_ );
202  }
203 }
204 
205 /* SVD along mode 1, 2, and 3*/
206 
207 VMML_TEMPLATE_STRING
208 void
209 VMML_TEMPLATE_CLASSNAME::svd_mode1( const t3_type& data_, u1_type& u1_ )
210 {
211  u1_unfolded_type* u = new u1_unfolded_type; // -> u1
212  data_.lateral_unfolding_bwd( *u );
213 
214  get_svd_u_red( *u, u1_ );
215 
216  delete u;
217 }
218 
219 VMML_TEMPLATE_STRING
220 void
221 VMML_TEMPLATE_CLASSNAME::svd_mode2( const t3_type& data_, u2_type& u2_ )
222 {
223  u2_unfolded_type* u = new u2_unfolded_type; // -> u2
224  data_.frontal_unfolding_bwd( *u );
225 
226  get_svd_u_red( *u, u2_ );
227 
228  delete u;
229 }
230 
231 VMML_TEMPLATE_STRING
232 void
233 VMML_TEMPLATE_CLASSNAME::svd_mode3( const t3_type& data_, u3_type& u3_ )
234 {
235  u3_unfolded_type* u = new u3_unfolded_type; // -> u3
236  data_.horizontal_unfolding_bwd( *u );
237 
238  get_svd_u_red( *u, u3_ );
239 
240  delete u;
241 }
242 
243 /* EIGS for mode 1, 2 and 3*/
244 
245 VMML_TEMPLATE_STRING
246 void
247 VMML_TEMPLATE_CLASSNAME::eigs_mode1( const t3_type& data_, u1_type& u1_ )
248 {
249  //unfolding / matricization
250  u1_unfolded_type* m_lateral = new u1_unfolded_type; // -> u1
251  data_.lateral_unfolding_bwd( *m_lateral );
252 
253  //covariance matrix of unfolded data
254  u1_cov_type* cov = new u1_cov_type;
255 #if 1
256  blas_dgemm< I1, I2*I3, I1, T>* blas_cov = new blas_dgemm< I1, I2*I3, I1, T>;
257  blas_cov->compute( *m_lateral, *cov );
258 #else
259  blas_daxpy< I1, T>* blas_cov = new blas_daxpy< I1, T>;
260  blas_cov->compute_mmm( *m_lateral, *cov );
261 #endif
262  delete blas_cov;
263  delete m_lateral;
264 
265  //compute x largest magnitude eigenvalues; x = R
266  get_eigs_u_red( *cov, u1_ );
267 
268  delete cov;
269 }
270 
271 VMML_TEMPLATE_STRING
272 void
273 VMML_TEMPLATE_CLASSNAME::eigs_mode2( const t3_type& data_, u2_type& u2_ )
274 {
275  //unfolding / matricization
276  u2_unfolded_type* m_frontal = new u2_unfolded_type; // -> u2
277  data_.frontal_unfolding_bwd( *m_frontal );
278 
279  //covariance matrix of unfolded data
280  u2_cov_type* cov = new u2_cov_type;
281 #if 1
282  blas_dgemm< I2, I1*I3, I2, T>* blas_cov = new blas_dgemm< I2, I1*I3, I2, T>;
283  blas_cov->compute( *m_frontal, *cov );
284 #else
285  blas_daxpy< I2, T>* blas_cov = new blas_daxpy< I2, T>;
286  blas_cov->compute_mmm( *m_frontal, *cov );
287 #endif
288 
289  delete blas_cov;
290  delete m_frontal;
291 
292  //compute x largest magnitude eigenvalues; x = R
293  get_eigs_u_red( *cov, u2_ );
294 
295  delete cov;
296 }
297 
298 VMML_TEMPLATE_STRING
299 void
300 VMML_TEMPLATE_CLASSNAME::eigs_mode3( const t3_type& data_, u3_type& u3_)
301 {
302  //unfolding / matricization
303  u3_unfolded_type* m_horizontal = new u3_unfolded_type; // -> u3
304  data_.horizontal_unfolding_bwd( *m_horizontal );
305 
306  //covariance matrix of unfolded data
307  u3_cov_type* cov = new u3_cov_type;
308 #if 1
309  blas_dgemm< I3, I1*I2, I3, T>* blas_cov = new blas_dgemm< I3, I1*I2, I3, T>;
310  blas_cov->compute( *m_horizontal, *cov );
311 #else
312  blas_daxpy< I3, T>* blas_cov = new blas_daxpy< I3, T>;
313  blas_cov->compute_mmm( *m_horizontal, *cov );
314 #endif
315 
316  delete blas_cov;
317  delete m_horizontal;
318 
319  //compute x largest magnitude eigenvalues; x = R
320  get_eigs_u_red( *cov, u3_ );
321 
322  delete cov;
323 }
324 
325 
326 
327 /* helper methods for SVD and EIGS*/
328 
329 VMML_TEMPLATE_STRING
330 template< size_t N, size_t R >
331 void
332 VMML_TEMPLATE_CLASSNAME::get_eigs_u_red( const matrix< N, N, T >& data_, matrix< N, R, T >& u_ )
333 {
334  typedef matrix< N, N, T_svd > cov_matrix_type;
335  typedef vector< R, T_svd > eigval_type;
336  typedef matrix< N, R, T_svd > eigvec_type;
337  //typedef matrix< N, R, T_coeff > coeff_type;
338 
339  //compute x largest magnitude eigenvalues; x = R
340  eigval_type* eigxvalues = new eigval_type;
341  eigvec_type* eigxvectors = new eigvec_type;
342 
343  lapack_sym_eigs< N, T_svd > eigs;
344  cov_matrix_type* data = new cov_matrix_type;
345  data->cast_from( data_ );
346  if( eigs.compute_x( *data, *eigxvectors, *eigxvalues) ) {
347 
348  /*if( _is_quantify_coeff ){
349  coeff_type* evec_quant = new coeff_type;
350  T min_value = 0; T max_value = 0;
351  u_.cast_from( *eigxvectors );
352  u_.quantize( *evec_quant, min_value, max_value );
353  evec_quant->dequantize( u_, min_value, max_value );
354  delete evec_quant;
355  } else */ if ( sizeof( T ) != 4 ){
356  u_.cast_from( *eigxvectors );
357  } else {
358  u_ = *eigxvectors;
359  }
360 
361  } else {
362  u_.zero();
363  }
364 
365  delete eigxvalues;
366  delete eigxvectors;
367  delete data;
368 
369 }
370 
371 VMML_TEMPLATE_STRING
372 template< size_t M, size_t N, size_t R >
373 void
374 VMML_TEMPLATE_CLASSNAME::get_svd_u_red( const matrix< M, N, T >& data_, matrix< M, R, T >& u_ )
375 {
376  typedef matrix< M, N, T_svd > svd_m_type;
377  //FIXME: typedef matrix< M, N, T_coeff > coeff_type;
378  typedef matrix< M, N, T > m_type;
379  typedef vector< N, T_svd > lambdas_type;
380 
381  svd_m_type* u_double = new svd_m_type;
382  u_double->cast_from( data_ );
383 
384  //FIXME:: coeff_type* u_quant = new coeff_type;
385  m_type* u_out = new m_type;
386 
387  lambdas_type* lambdas = new lambdas_type;
388  lapack_svd< M, N, T_svd >* svd = new lapack_svd< M, N, T_svd >();
389  if( svd->compute_and_overwrite_input( *u_double, *lambdas )) {
390 
391  /* if( _is_quantify_coeff ){
392  T min_value = 0; T max_value = 0;
393  u_comp->cast_from( *u_double );
394  //FIXME: u_comp->quantize( *u_quant, min_value, max_value );
395  //FIXME: u_quant->dequantize( *u_internal, min_value, max_value );
396  } else */ if ( sizeof( T ) != 4 ){
397  u_out->cast_from( *u_double );
398  } else {
399  *u_out = *u_double;
400  }
401 
402  u_out->get_sub_matrix( u_ );
403 
404  } else {
405  u_.zero();
406  }
407 
408  delete lambdas;
409  delete svd;
410  delete u_double;
411  //delete u_quant;
412  delete u_out;
413 }
414 
415 #undef VMML_TEMPLATE_STRING
416 #undef VMML_TEMPLATE_CLASSNAME
417 
418 }//end vmml namespace
419 
420 #endif
421 
vmml::t3_hosvd
Definition: t3_hosvd.hpp:63