SHOGUN  4.1.0
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
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
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>();
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
67  std::cout << "cov" << std::endl;
68  std::cout << cov << std::endl;
69  #endif
70
71  // Whitening & Projection onto signal subspace
73  eig.compute(cov);
74
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
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
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
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
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
int32_t index_t
Definition: common.h:62
static SGMatrix< float64_t > diagonalize(SGNDArray< float64_t > C, SGMatrix< float64_t > V0=SGMatrix< float64_t >(NULL, 0, 0, false), double eps=CMath::MACHINE_EPSILON, int itermax=200)
Definition: SGMatrix.h:20
index_t num_cols
Definition: SGMatrix.h:378
void init()
#define SG_REF(x)
Definition: SGObject.h:51
index_t num_rows
Definition: SGMatrix.h:376
T * get_matrix(index_t matIdx) const
Definition: SGNDArray.h:72
#define ASSERT(x)
Definition: SGIO.h:201
double float64_t
Definition: common.h:50
virtual CFeatures * apply(CFeatures *features)
SGMatrix< float64_t > get_cumulant_matrix() const
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
class ICAConverter Base class for ICA algorithms
Definition: ICAConverter.h:27
The class Features is the base class of all feature objects.
Definition: Features.h:68
SGMatrix< float64_t > m_mixing_matrix
Definition: ICAConverter.h:83
static void inverse(SGMatrix< float64_t > matrix)
inverses square matrix in-place
Definition: SGMatrix.cpp:885