SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Jade.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Kevin Hughes
8  */
9 
11 
13 
14 #ifdef HAVE_EIGEN3
15 
19 
20 #ifdef DEBUG_JADE
21 #include <iostream>
22 #endif
23 
24 using namespace shogun;
25 using namespace Eigen;
26 
28 {
29  init();
30 }
31 
33 {
34  m_cumulant_matrix = SGMatrix<float64_t>();
35  SG_ADD(&m_cumulant_matrix, "cumulant_matrix", "m_cumulant_matrix", MS_NOT_AVAILABLE);
36 }
37 
39 {
40 }
41 
43 {
44  return m_cumulant_matrix;
45 }
46 
48 {
49  ASSERT(features);
50  SG_REF(features);
51 
52  SGMatrix<float64_t> X = ((CDenseFeatures<float64_t>*)features)->get_feature_matrix();
53 
54  int n = X.num_rows;
55  int T = X.num_cols;
56  int m = n;
57 
58  Map<MatrixXd> EX(X.matrix,n,T);
59 
60  // Mean center X
61  VectorXd mean = (EX.rowwise().sum() / (float64_t)T);
62  MatrixXd SPX = EX.colwise() - mean;
63 
64  MatrixXd cov = (SPX * SPX.transpose()) / (float64_t)T;
65 
66  #ifdef DEBUG_JADE
67  std::cout << "cov" << std::endl;
68  std::cout << cov << std::endl;
69  #endif
70 
71  // Whitening & Projection onto signal subspace
72  SelfAdjointEigenSolver<MatrixXd> eig;
73  eig.compute(cov);
74 
75  #ifdef DEBUG_JADE
76  std::cout << "eigenvectors" << std::endl;
77  std::cout << eig.eigenvectors() << std::endl;
78 
79  std::cout << "eigenvalues" << std::endl;
80  std::cout << eig.eigenvalues().asDiagonal() << std::endl;
81  #endif
82 
83  // Scaling
84  VectorXd scales = eig.eigenvalues().cwiseSqrt();
85  MatrixXd B = scales.cwiseInverse().asDiagonal() * eig.eigenvectors().transpose();
86 
87  #ifdef DEBUG_JADE
88  std::cout << "whitener" << std::endl;
89  std::cout << B << std::endl;
90  #endif
91 
92  // Sphering
93  SPX = B * SPX;
94 
95  // Estimation of the cumulant matrices
96  int dimsymm = (m * ( m + 1)) / 2; // Dim. of the space of real symm matrices
97  int nbcm = dimsymm; // number of cumulant matrices
98  m_cumulant_matrix = SGMatrix<float64_t>(m,m*nbcm); // Storage for cumulant matrices
99  Map<MatrixXd> CM(m_cumulant_matrix.matrix,m,m*nbcm);
100  MatrixXd R(m,m); R.setIdentity();
101  MatrixXd Qij = MatrixXd::Zero(m,m); // Temp for a cum. matrix
102  VectorXd Xim = VectorXd::Zero(m); // Temp
103  VectorXd Xjm = VectorXd::Zero(m); // Temp
104  VectorXd Xijm = VectorXd::Zero(m); // Temp
105  int Range = 0;
106 
107  for (int im = 0; im < m; im++)
108  {
109  Xim = SPX.row(im);
110  Xijm = Xim.cwiseProduct(Xim);
111  Qij = SPX.cwiseProduct(Xijm.replicate(1,m).transpose()) * SPX.transpose() / (float)T - R - 2*R.col(im)*R.col(im).transpose();
112  CM.block(0,Range,m,m) = Qij;
113  Range = Range + m;
114  for (int jm = 0; jm < im; jm++)
115  {
116  Xjm = SPX.row(jm);
117  Xijm = Xim.cwiseProduct(Xjm);
118  Qij = SPX.cwiseProduct(Xijm.replicate(1,m).transpose()) * SPX.transpose() / (float)T - R.col(im)*R.col(jm).transpose() - R.col(jm)*R.col(im).transpose();
119  CM.block(0,Range,m,m) = sqrt(2)*Qij;
120  Range = Range + m;
121  }
122  }
123 
124  #ifdef DEBUG_JADE
125  std::cout << "cumulatant matrices" << std::endl;
126  std::cout << CM << std::endl;
127  #endif
128 
129  // Stack cumulant matrix into ND Array
130  index_t * M_dims = SG_MALLOC(index_t, 3);
131  M_dims[0] = m;
132  M_dims[1] = m;
133  M_dims[2] = nbcm;
134  SGNDArray< float64_t > M(M_dims, 3);
135 
136  for (int i = 0; i < nbcm; i++)
137  {
138  Map<MatrixXd> EM(M.get_matrix(i),m,m);
139  EM = CM.block(0,i*m,m,m);
140  }
141 
142  // Diagonalize
144  Map<MatrixXd> EQ(Q.matrix,m,m);
145  EQ = -1 * EQ.inverse();
146 
147  #ifdef DEBUG_JADE
148  std::cout << "diagonalizer" << std::endl;
149  std::cout << EQ << std::endl;
150  #endif
151 
152  // Separating matrix
153  SGMatrix<float64_t> sep_matrix = SGMatrix<float64_t>(m,m);
154  Map<MatrixXd> C(sep_matrix.matrix,m,m);
155  C = EQ.transpose() * B;
156 
157  // Sort
158  VectorXd A = (B.inverse()*EQ).cwiseAbs2().colwise().sum();
159  bool swap = false;
160  do
161  {
162  swap = false;
163  for (int j = 1; j < n; j++)
164  {
165  if ( A(j) < A(j-1) )
166  {
167  std::swap(A(j),A(j-1));
168  C.col(j).swap(C.col(j-1));
169  swap = true;
170  }
171  }
172 
173  } while(swap);
174 
175  for (int j = 0; j < m/2; j++)
176  C.row(j).swap( C.row(m-1-j) );
177 
178  // Fix Signs
179  VectorXd signs = VectorXd::Zero(m);
180  for (int i = 0; i < m; i++)
181  {
182  if( C(i,0) < 0 )
183  signs(i) = -1;
184  else
185  signs(i) = 1;
186  }
187  C = signs.asDiagonal() * C;
188 
189  #ifdef DEBUG_JADE
190  std::cout << "un mixing matrix" << std::endl;
191  std::cout << C << std::endl;
192  #endif
193 
194  // Unmix
195  EX = C * EX;
196 
198  Map<MatrixXd> Emixing_matrix(m_mixing_matrix.matrix,m,m);
199  Emixing_matrix = C.inverse();
200 
201  return features;
202 }
203 
204 #endif // HAVE_EIGEN3

SHOGUN Machine Learning Toolbox - Documentation