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

SHOGUN Machine Learning Toolbox - Documentation