SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
FastICA.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation