SHOGUN  6.1.3
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  REQUIRE(X.data(), "Features have not been provided.\n");
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 
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 
static T sqrt(T x)
Definition: Math.h:428
virtual ~CFastICA()
Definition: FastICA.cpp:59
Definition: SGMatrix.h:25
static float64_t randn_double()
Definition: Math.h:924
#define REQUIRE(x,...)
Definition: SGIO.h:181
std::enable_if<!std::is_same< T, complex128_t >::value, float64_t >::type mean(const Container< T > &a)
#define SG_REF(x)
Definition: SGObject.h:52
#define ASSERT(x)
Definition: SGIO.h:176
double float64_t
Definition: common.h:60
index_t num_rows
Definition: SGMatrix.h:495
index_t num_cols
Definition: SGMatrix.h:497
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:604
The class Features is the base class of all feature objects.
Definition: Features.h:69
bool get_whiten() const
Definition: FastICA.cpp:68
T * data() const
Definition: SGMatrix.h:269
SGMatrix< float64_t > m_mixing_matrix
Definition: ICAConverter.h:82
void svd(const SGMatrix< T > &A, SGVector< T > &s, SGMatrix< T > &U, bool thin_U=true, SVDAlgorithm alg=SVDAlgorithm::BidiagonalDivideConquer)
void set_whiten(bool whiten)
Definition: FastICA.cpp:63
#define SG_ADD(...)
Definition: SGObject.h:93

SHOGUN Machine Learning Toolbox - Documentation