SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FastICA.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  * ported from scikit-learn
9  */
10 
12 
14 
15 #ifdef HAVE_EIGEN3
16 
19 
20 using namespace shogun;
21 using namespace Eigen;
22 
23 namespace {
24 
25  MatrixXd sym_decorrelation(MatrixXd W)
26  {
27  MatrixXd K = W * W.transpose();
28 
29  SelfAdjointEigenSolver<MatrixXd> eig;
30  eig.compute(K);
31 
32  return ((eig.eigenvectors() * eig.eigenvalues().cwiseSqrt().asDiagonal().inverse()) * eig.eigenvectors().transpose()) * W;
33  }
34 
35  float64_t alpha = 1.0; // alpha must be in range [1.0 - 2.0]
36 
37  float64_t gx(float64_t x)
38  {
39  return CMath::tanh(x*alpha);
40  }
41 
42  float64_t g_x(float64_t x)
43  {
44  return alpha * (1.0 - pow(gx(x),2));
45  }
46 
47 };
48 
50 {
51  init();
52 }
53 
55 {
56  whiten = true;
57  SG_ADD(&whiten, "whiten", "flag indicating whether to whiten the data", MS_NOT_AVAILABLE);
58 }
59 
61 {
62 }
63 
64 void CFastICA::set_whiten(bool _whiten)
65 {
66  whiten = _whiten;
67 }
68 
70 {
71  return whiten;
72 }
73 
75 {
76  ASSERT(features);
77  SG_REF(features);
78 
79  SGMatrix<float64_t> X = ((CDenseFeatures<float64_t>*)features)->get_feature_matrix();
80 
81  int n = X.num_rows;
82  int p = X.num_cols;
83  int m = n;
84 
85  Map<MatrixXd> EX(X.matrix,n,p);
86 
87  // Whiten
88  MatrixXd K;
89  MatrixXd WX;
90  if (whiten)
91  {
92  VectorXd mean = (EX.rowwise().sum() / (float64_t)p);
93  MatrixXd SPX = EX.colwise() - mean;
94 
95  Eigen::JacobiSVD<MatrixXd> svd;
96  svd.compute(SPX, Eigen::ComputeThinU);
97 
98  MatrixXd u = svd.matrixU();
99  MatrixXd d = svd.singularValues();
100 
101  // for matching numpy/scikit-learn
102  //u.rightCols(u.cols() - 1) *= -1;
103 
104  // see Hyvarinen (6.33) p.140
105  K = u.transpose();
106  for (int r = 0; r < K.rows(); r++)
107  K.row(r) /= d(r);
108 
109  // see Hyvarinen (13.6) p.267 Here WX is white and data
110  // in X has been projected onto a subspace by PCA
111  WX = K * SPX;
112  WX *= CMath::sqrt((float64_t)p);
113  }
114  else
115  {
116  WX = EX;
117  }
118 
119  // Initial mixing matrix estimate
121  {
123 
124  for (int i = 0; i < m; i++)
125  {
126  for (int j = 0; j < m; j++)
128  }
129  }
130 
131  Map<MatrixXd> W(m_mixing_matrix.matrix, m, m);
132 
133  W = sym_decorrelation(W);
134 
135  int iter = 0;
136  float64_t lim = tol+1;
137  while (lim > tol && iter < max_iter)
138  {
139  MatrixXd wtx = W * WX;
140 
141  MatrixXd gwtx = wtx.unaryExpr(std::ptr_fun(&gx));
142  MatrixXd g_wtx = wtx.unaryExpr(std::ptr_fun(&g_x));
143 
144  MatrixXd W1 = (gwtx * WX.transpose()) / (float64_t)p - (g_wtx.rowwise().sum()/(float64_t)p).asDiagonal() * W;
145 
146  W1 = sym_decorrelation(W1);
147 
148  lim = ((W1 * W.transpose()).diagonal().cwiseAbs().array() - 1).abs().maxCoeff();
149 
150  W = W1;
151 
152  iter++;
153  }
154 
155  // Unmix
156  if (whiten)
157  W = (W*K);
158 
159  EX = W * EX;
160 
161  // set m_mixing_matrix
162  W = W.inverse();
163 
164  return features;
165 }
166 
167 #endif // HAVE_EIGEN3

SHOGUN Machine Learning Toolbox - Documentation