SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2013 Kevin Hughes 00008 * ported from scikit-learn 00009 */ 00010 00011 #include <shogun/converter/ica/FastICA.h> 00012 00013 #include <shogun/features/DenseFeatures.h> 00014 00015 #ifdef HAVE_EIGEN3 00016 00017 #include <shogun/mathematics/Math.h> 00018 #include <shogun/mathematics/eigen3.h> 00019 00020 using namespace shogun; 00021 using namespace Eigen; 00022 00023 namespace { 00024 00025 MatrixXd sym_decorrelation(MatrixXd W) 00026 { 00027 MatrixXd K = W * W.transpose(); 00028 00029 SelfAdjointEigenSolver<MatrixXd> eig; 00030 eig.compute(K); 00031 00032 return ((eig.eigenvectors() * eig.eigenvalues().cwiseSqrt().asDiagonal().inverse()) * eig.eigenvectors().transpose()) * W; 00033 } 00034 00035 float64_t alpha = 1.0; // alpha must be in range [1.0 - 2.0] 00036 00037 float64_t gx(float64_t x) 00038 { 00039 return CMath::tanh(x*alpha); 00040 } 00041 00042 float64_t g_x(float64_t x) 00043 { 00044 return alpha * (1.0 - pow(gx(x),2)); 00045 } 00046 00047 }; 00048 00049 CFastICA::CFastICA() : CICAConverter() 00050 { 00051 init(); 00052 } 00053 00054 void CFastICA::init() 00055 { 00056 whiten = true; 00057 SG_ADD(&whiten, "whiten", "flag indicating whether to whiten the data", MS_NOT_AVAILABLE); 00058 } 00059 00060 CFastICA::~CFastICA() 00061 { 00062 } 00063 00064 void CFastICA::set_whiten(bool _whiten) 00065 { 00066 whiten = _whiten; 00067 } 00068 00069 bool CFastICA::get_whiten() const 00070 { 00071 return whiten; 00072 } 00073 00074 CFeatures* CFastICA::apply(CFeatures* features) 00075 { 00076 ASSERT(features); 00077 SG_REF(features); 00078 00079 SGMatrix<float64_t> X = ((CDenseFeatures<float64_t>*)features)->get_feature_matrix(); 00080 00081 int n = X.num_rows; 00082 int p = X.num_cols; 00083 int m = n; 00084 00085 Map<MatrixXd> EX(X.matrix,n,p); 00086 00087 // Whiten 00088 MatrixXd K; 00089 MatrixXd WX; 00090 if (whiten) 00091 { 00092 VectorXd mean = (EX.rowwise().sum() / (float64_t)p); 00093 MatrixXd SPX = EX.colwise() - mean; 00094 00095 Eigen::JacobiSVD<MatrixXd> svd; 00096 svd.compute(SPX, Eigen::ComputeThinU); 00097 00098 MatrixXd u = svd.matrixU(); 00099 MatrixXd d = svd.singularValues(); 00100 00101 // for matching numpy/scikit-learn 00102 //u.rightCols(u.cols() - 1) *= -1; 00103 00104 // see Hyvarinen (6.33) p.140 00105 K = u.transpose(); 00106 for (int r = 0; r < K.rows(); r++) 00107 K.row(r) /= d(r); 00108 00109 // see Hyvarinen (13.6) p.267 Here WX is white and data 00110 // in X has been projected onto a subspace by PCA 00111 WX = K * SPX; 00112 WX *= CMath::sqrt((float64_t)p); 00113 } 00114 else 00115 { 00116 WX = EX; 00117 } 00118 00119 // Initial mixing matrix estimate 00120 if (m_mixing_matrix.num_rows != m || m_mixing_matrix.num_cols != m) 00121 { 00122 m_mixing_matrix = SGMatrix<float64_t>(m,m); 00123 00124 for (int i = 0; i < m; i++) 00125 { 00126 for (int j = 0; j < m; j++) 00127 m_mixing_matrix(i,j) = CMath::randn_double(); 00128 } 00129 } 00130 00131 Map<MatrixXd> W(m_mixing_matrix.matrix, m, m); 00132 00133 W = sym_decorrelation(W); 00134 00135 int iter = 0; 00136 float64_t lim = tol+1; 00137 while (lim > tol && iter < max_iter) 00138 { 00139 MatrixXd wtx = W * WX; 00140 00141 MatrixXd gwtx = wtx.unaryExpr(std::ptr_fun(&gx)); 00142 MatrixXd g_wtx = wtx.unaryExpr(std::ptr_fun(&g_x)); 00143 00144 MatrixXd W1 = (gwtx * WX.transpose()) / (float64_t)p - (g_wtx.rowwise().sum()/(float64_t)p).asDiagonal() * W; 00145 00146 W1 = sym_decorrelation(W1); 00147 00148 lim = ((W1 * W.transpose()).diagonal().cwiseAbs().array() - 1).abs().maxCoeff(); 00149 00150 W = W1; 00151 00152 iter++; 00153 } 00154 00155 // Unmix 00156 if (whiten) 00157 W = (W*K); 00158 00159 EX = W * EX; 00160 00161 // set m_mixing_matrix 00162 W = W.inverse(); 00163 00164 return features; 00165 } 00166 00167 #endif // HAVE_EIGEN3