fastica.cpp
Go to the documentation of this file.
#include <itpp/signal/fastica.h>
#include <itpp/signal/sigfun.h>
#include <itpp/signal/resampling.h>
#include <itpp/base/algebra/eigen.h>
#include <itpp/base/algebra/svd.h>
#include <itpp/base/math/trig_hyp.h>
#include <itpp/base/matfunc.h>
#include <itpp/base/random.h>
#include <itpp/base/sort.h>
#include <itpp/base/specmat.h>
#include <itpp/base/svec.h>
#include <itpp/base/math/min_max.h>
#include <itpp/stat/misc_stat.h>
using namespace itpp;
static void selcol(const mat oldMatrix, const vec maskVector, mat & newMatrix);
static int pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat & Es, vec & Ds);
static void remmean(mat inVectors, mat & outVectors, vec & meanValue);
static void whitenv(const mat vectors, const mat E, const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix);
static mat orth(const mat A);
static mat mpower(const mat A, const double y);
static ivec getSamples(const int max, const double percentage);
static vec sumcol(const mat A);
static bool fpica(const mat X, const mat whiteningMatrix, const mat dewhiteningMatrix, const int approach, const int numOfIC, const int g, const int finetune, const double a1, const double a2, double myy, const int stabilization, const double epsilon, const int maxNumIterations, const int maxFinetune, const int initState, mat guess, double sampleSize, mat & A, mat & W);
namespace itpp
{
// Constructor, init default values
Fast_ICA::Fast_ICA(mat ma_mixedSig)
{
// Init default values
approach = FICA_APPROACH_SYMM;
g = FICA_NONLIN_POW3;
finetune = true;
a1 = 1.0;
a2 = 1.0;
mu = 1.0;
epsilon = 0.0001;
sampleSize = 1.0;
stabilization = false;
maxNumIterations = 100000;
maxFineTune = 100;
firstEig = 1;
mixedSig = ma_mixedSig;
lastEig = mixedSig.rows();
numOfIC = mixedSig.rows();
PCAonly = false;
initState = FICA_INIT_RAND;
}
// Call main function
bool Fast_ICA::separate(void)
{
int Dim = numOfIC;
mat mixedSigC;
vec mixedMean;
mat guess;
if (initState == FICA_INIT_RAND)
guess = zeros(Dim, Dim);
else
guess = mat(initGuess);
VecPr = zeros(mixedSig.rows(), numOfIC);
icasig = zeros(numOfIC, mixedSig.cols());
remmean(mixedSig, mixedSigC, mixedMean);
if (pcamat(mixedSigC, numOfIC, firstEig, lastEig, E, D) < 1) {
// no principal components could be found (e.g. all-zero data): return the unchanged input
icasig = mixedSig;
return false;
}
whitenv(mixedSigC, E, diag(D), whitesig, whiteningMatrix, dewhiteningMatrix);
Dim = whitesig.rows();
if (numOfIC > Dim) numOfIC = Dim;
ivec NcFirst = to_ivec(zeros(numOfIC));
vec NcVp = D;
for (int i = 0; i < NcFirst.size(); i++) {
NcFirst(i) = max_index(NcVp);
NcVp(NcFirst(i)) = 0.0;
VecPr.set_col(i, dewhiteningMatrix.get_col(i));
}
bool result = true;
if (PCAonly == false) {
result = fpica(whitesig, whiteningMatrix, dewhiteningMatrix, approach, numOfIC, g, finetune, a1, a2, mu, stabilization, epsilon, maxNumIterations, maxFineTune, initState, guess, sampleSize, A, W);
icasig = W * mixedSig;
}
else { // PCA only : returns E as IcaSig
icasig = VecPr;
}
return result;
}
void Fast_ICA::set_approach(int in_approach) { approach = in_approach; if (approach == FICA_APPROACH_DEFL) finetune = true; }
void Fast_ICA::set_nrof_independent_components(int in_nrIC) { numOfIC = in_nrIC; }
void Fast_ICA::set_non_linearity(int in_g) { g = in_g; }
void Fast_ICA::set_fine_tune(bool in_finetune) { finetune = in_finetune; }
void Fast_ICA::set_a1(double fl_a1) { a1 = fl_a1; }
void Fast_ICA::set_a2(double fl_a2) { a2 = fl_a2; }
void Fast_ICA::set_mu(double fl_mu) { mu = fl_mu; }
void Fast_ICA::set_epsilon(double fl_epsilon) { epsilon = fl_epsilon; }
void Fast_ICA::set_sample_size(double fl_sampleSize) { sampleSize = fl_sampleSize; }
void Fast_ICA::set_stabilization(bool in_stabilization) { stabilization = in_stabilization; }
void Fast_ICA::set_max_num_iterations(int in_maxNumIterations) { maxNumIterations = in_maxNumIterations; }
void Fast_ICA::set_max_fine_tune(int in_maxFineTune) { maxFineTune = in_maxFineTune; }
void Fast_ICA::set_first_eig(int in_firstEig) { firstEig = in_firstEig; }
void Fast_ICA::set_last_eig(int in_lastEig) { lastEig = in_lastEig; }
void Fast_ICA::set_pca_only(bool in_PCAonly) { PCAonly = in_PCAonly; }
void Fast_ICA::set_init_guess(mat ma_initGuess)
{
initGuess = ma_initGuess;
initState = FICA_INIT_GUESS;
}
mat Fast_ICA::get_mixing_matrix() { if (PCAonly) { it_warning("No ICA performed."); return (zeros(1, 1));} else return A; }
mat Fast_ICA::get_separating_matrix() { if (PCAonly) { it_warning("No ICA performed."); return(zeros(1, 1)); } else return W; }
mat Fast_ICA::get_independent_components() { if (PCAonly) { it_warning("No ICA performed."); return(zeros(1, 1)); } else return icasig; }
int Fast_ICA::get_nrof_independent_components() { return numOfIC; }
mat Fast_ICA::get_principal_eigenvectors() { return VecPr; }
mat Fast_ICA::get_whitening_matrix() { return whiteningMatrix; }
mat Fast_ICA::get_dewhitening_matrix() { return dewhiteningMatrix; }
mat Fast_ICA::get_white_sig() { return whitesig; }
} // namespace itpp
static void selcol(const mat oldMatrix, const vec maskVector, mat & newMatrix)
{
int numTaken = 0;
for (int i = 0; i < size(maskVector); i++) if (maskVector(i) == 1) numTaken++;
newMatrix = zeros(oldMatrix.rows(), numTaken);
numTaken = 0;
for (int i = 0; i < size(maskVector); i++) {
if (maskVector(i) == 1) {
newMatrix.set_col(numTaken, oldMatrix.get_col(i));
numTaken++;
}
}
return;
}
static int pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat & Es, vec & Ds)
{
mat Et;
vec Dt;
cmat Ec;
cvec Dc;
double lowerLimitValue = 0.0,
higherLimitValue = 0.0;
int oldDimension = vectors.rows();
mat covarianceMatrix = cov(transpose(vectors), 0);
eig_sym(covarianceMatrix, Dt, Et);
int maxLastEig = 0;
// Compute rank
for (int i = 0; i < Dt.length(); i++) if (Dt(i) > FICA_TOL) maxLastEig++;
if (maxLastEig < 1) return 0;
// Force numOfIC components
if (maxLastEig > numOfIC) maxLastEig = numOfIC;
vec eigenvalues = zeros(size(Dt));
vec eigenvalues2 = zeros(size(Dt));
eigenvalues2 = Dt;
sort(eigenvalues2);
vec lowerColumns = zeros(size(Dt));
for (int i = 0; i < size(Dt); i++) eigenvalues(i) = eigenvalues2(size(Dt) - i - 1);
if (lastEig > maxLastEig) lastEig = maxLastEig;
if (lastEig < oldDimension) lowerLimitValue = (eigenvalues(lastEig - 1) + eigenvalues(lastEig)) / 2;
else lowerLimitValue = eigenvalues(oldDimension - 1) - 1;
for (int i = 0; i < size(Dt); i++) if (Dt(i) > lowerLimitValue) lowerColumns(i) = 1;
if (firstEig > 1) higherLimitValue = (eigenvalues(firstEig - 2) + eigenvalues(firstEig - 1)) / 2;
else higherLimitValue = eigenvalues(0) + 1;
vec higherColumns = zeros(size(Dt));
for (int i = 0; i < size(Dt); i++) if (Dt(i) < higherLimitValue) higherColumns(i) = 1;
vec selectedColumns = zeros(size(Dt));
for (int i = 0; i < size(Dt); i++) selectedColumns(i) = (lowerColumns(i) == 1 && higherColumns(i) == 1) ? 1 : 0;
selcol(Et, selectedColumns, Es);
int numTaken = 0;
for (int i = 0; i < size(selectedColumns); i++) if (selectedColumns(i) == 1) numTaken++;
Ds = zeros(numTaken);
numTaken = 0;
for (int i = 0; i < size(Dt); i++)
if (selectedColumns(i) == 1) {
Ds(numTaken) = Dt(i);
numTaken++;
}
return lastEig;
}
static void remmean(mat inVectors, mat & outVectors, vec & meanValue)
{
outVectors = zeros(inVectors.rows(), inVectors.cols());
meanValue = zeros(inVectors.rows());
for (int i = 0; i < inVectors.rows(); i++) {
meanValue(i) = mean(inVectors.get_row(i));
for (int j = 0; j < inVectors.cols(); j++) outVectors(i, j) = inVectors(i, j) - meanValue(i);
}
}
static void whitenv(const mat vectors, const mat E, const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix)
{
whiteningMatrix = zeros(E.cols(), E.rows());
dewhiteningMatrix = zeros(E.rows(), E.cols());
for (int i = 0; i < D.cols(); i++) {
whiteningMatrix.set_row(i, std::pow(std::sqrt(D(i, i)), -1)*E.get_col(i));
dewhiteningMatrix.set_col(i, std::sqrt(D(i, i))*E.get_col(i));
}
newVectors = whiteningMatrix * vectors;
return;
}
static mat orth(const mat A)
{
mat Q;
mat U, V;
vec S;
double eps = 2.2e-16;
double tol = 0.0;
int mmax = 0;
int r = 0;
svd(A, U, S, V);
if (A.rows() > A.cols()) {
U = U(0, U.rows() - 1, 0, A.cols() - 1);
S = S(0, A.cols() - 1);
}
mmax = (A.rows() > A.cols()) ? A.rows() : A.cols();
tol = mmax * eps * max(S);
for (int i = 0; i < size(S); i++) if (S(i) > tol) r++;
Q = U(0, U.rows() - 1, 0, r - 1);
return (Q);
}
static mat mpower(const mat A, const double y)
{
mat T = zeros(A.rows(), A.cols());
mat dd = zeros(A.rows(), A.cols());
vec d = zeros(A.rows());
vec dOut = zeros(A.rows());
eig_sym(A, d, T);
dOut = pow(d, y);
diag(dOut, dd);
for (int i = 0; i < T.cols(); i++) T.set_col(i, T.get_col(i) / norm(T.get_col(i)));
return (T*dd*transpose(T));
}
static ivec getSamples(const int max, const double percentage)
{
vec rd = randu(max);
sparse_vec sV;
ivec out;
int sZ = 0;
for (int i = 0; i < max; i++) if (rd(i) < percentage) { sV.add_elem(sZ, i); sZ++; }
out = to_ivec(full(sV));
return (out);
}
static vec sumcol(const mat A)
{
vec out = zeros(A.cols());
for (int i = 0; i < A.cols(); i++) { out(i) = sum(A.get_col(i)); }
return (out);
}
static bool fpica(const mat X, const mat whiteningMatrix, const mat dewhiteningMatrix, const int approach, const int numOfIC, const int g, const int finetune, const double a1, const double a2, double myy, const int stabilization, const double epsilon, const int maxNumIterations, const int maxFinetune, const int initState, mat guess, double sampleSize, mat & A, mat & W)
{
int vectorSize = X.rows();
int numSamples = X.cols();
int gOrig = g;
int gFine = finetune + 1;
double myyOrig = myy;
double myyK = 0.01;
int failureLimit = 5;
int usedNlinearity = 0;
double stroke = 0.0;
int notFine = 1;
int loong = 0;
int initialStateMode = initState;
double minAbsCos = 0.0, minAbsCos2 = 0.0;
if (sampleSize * numSamples < 1000) sampleSize = (1000 / (double)numSamples < 1.0) ? 1000 / (double)numSamples : 1.0;
if (sampleSize != 1.0) gOrig += 2;
if (myy != 1.0) gOrig += 1;
int fineTuningEnabled = 1;
if (!finetune) {
if (myy != 1.0) gFine = gOrig;
else gFine = gOrig + 1;
fineTuningEnabled = 0;
}
int stabilizationEnabled = stabilization;
if (!stabilization && myy != 1.0) stabilizationEnabled = true;
usedNlinearity = gOrig;
if (initState == FICA_INIT_GUESS && guess.rows() != whiteningMatrix.cols()) {
initialStateMode = 0;
}
else if (guess.cols() < numOfIC) {
mat guess2 = randu(guess.rows(), numOfIC - guess.cols()) - 0.5;
guess = concat_horizontal(guess, guess2);
}
else if (guess.cols() > numOfIC) guess = guess(0, guess.rows() - 1, 0, numOfIC - 1);
if (approach == FICA_APPROACH_SYMM) {
usedNlinearity = gOrig;
stroke = 0;
notFine = 1;
loong = 0;
A = zeros(vectorSize, numOfIC);
mat B = zeros(vectorSize, numOfIC);
if (initialStateMode == 0) B = orth(randu(vectorSize, numOfIC) - 0.5);
else B = whiteningMatrix * guess;
mat BOld = zeros(B.rows(), B.cols());
mat BOld2 = zeros(B.rows(), B.cols());
for (int round = 0; round < maxNumIterations; round++) {
if (round == maxNumIterations - 1) {
// If there is a convergence problem,
// we still want ot return something.
// This is difference with original
// Matlab implementation.
A = dewhiteningMatrix * B;
W = transpose(B) * whiteningMatrix;
return false;
}
B = B * mpower(transpose(B) * B , -0.5);
minAbsCos = min(abs(diag(transpose(B) * BOld)));
minAbsCos2 = min(abs(diag(transpose(B) * BOld2)));
if (1 - minAbsCos < epsilon) {
if (fineTuningEnabled && notFine) {
notFine = 0;
usedNlinearity = gFine;
myy = myyK * myyOrig;
BOld = zeros(B.rows(), B.cols());
BOld2 = zeros(B.rows(), B.cols());
}
else {
A = dewhiteningMatrix * B;
break;
}
} // IF epsilon
else if (stabilizationEnabled) {
if (!stroke && (1 - minAbsCos2 < epsilon)) {
stroke = myy;
myy /= 2;
if (mod(usedNlinearity, 2) == 0) usedNlinearity += 1 ;
}
else if (stroke) {
myy = stroke;
stroke = 0;
if (myy == 1 && mod(usedNlinearity, 2) != 0) usedNlinearity -= 1;
}
else if (!loong && (round > maxNumIterations / 2)) {
loong = 1;
myy /= 2;
if (mod(usedNlinearity, 2) == 0) usedNlinearity += 1;
}
} // stabilizationEnabled
BOld2 = BOld;
BOld = B;
switch (usedNlinearity) {
// pow3
case FICA_NONLIN_POW3 : {
B = (X * pow(transpose(X) * B, 3)) / numSamples - 3 * B;
break;
}
case(FICA_NONLIN_POW3+1) : {
mat Y = transpose(X) * B;
mat Gpow3 = pow(Y, 3);
vec Beta = sumcol(pow(Y, 4));
mat D = diag(pow(Beta - 3 * numSamples , -1));
B = B + myy * B * (transpose(Y) * Gpow3 - diag(Beta)) * D;
break;
}
case(FICA_NONLIN_POW3+2) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
B = (Xsub * pow(transpose(Xsub) * B, 3)) / Xsub.cols() - 3 * B;
break;
}
case(FICA_NONLIN_POW3+3) : {
mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
mat Gpow3 = pow(Ysub, 3);
vec Beta = sumcol(pow(Ysub, 4));
mat D = diag(pow(Beta - 3 * Ysub.rows() , -1));
B = B + myy * B * (transpose(Ysub) * Gpow3 - diag(Beta)) * D;
break;
}
// TANH
case FICA_NONLIN_TANH : {
mat hypTan = tanh(a1 * transpose(X) * B);
B = (X * hypTan) / numSamples - elem_mult(reshape(repeat(sumcol(1 - pow(hypTan, 2)), B.rows()), B.rows(), B.cols()), B) / numSamples * a1;
break;
}
case(FICA_NONLIN_TANH+1) : {
mat Y = transpose(X) * B;
mat hypTan = tanh(a1 * Y);
vec Beta = sumcol(elem_mult(Y, hypTan));
vec Beta2 = sumcol(1 - pow(hypTan, 2));
mat D = diag(pow(Beta - a1 * Beta2 , -1));
B = B + myy * B * (transpose(Y) * hypTan - diag(Beta)) * D;
break;
}
case(FICA_NONLIN_TANH+2) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
mat hypTan = tanh(a1 * transpose(Xsub) * B);
B = (Xsub * hypTan) / Xsub.cols() - elem_mult(reshape(repeat(sumcol(1 - pow(hypTan, 2)), B.rows()), B.rows(), B.cols()), B) / Xsub.cols() * a1;
break;
}
case(FICA_NONLIN_TANH+3) : {
mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
mat hypTan = tanh(a1 * Ysub);
vec Beta = sumcol(elem_mult(Ysub, hypTan));
vec Beta2 = sumcol(1 - pow(hypTan, 2));
mat D = diag(pow(Beta - a1 * Beta2 , -1));
B = B + myy * B * (transpose(Ysub) * hypTan - diag(Beta)) * D;
break;
}
// GAUSS
case FICA_NONLIN_GAUSS : {
mat U = transpose(X) * B;
mat Usquared = pow(U, 2);
mat ex = exp(-a2 * Usquared / 2);
mat gauss = elem_mult(U, ex);
mat dGauss = elem_mult(1 - a2 * Usquared, ex);
B = (X * gauss) / numSamples - elem_mult(reshape(repeat(sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / numSamples;
break;
}
case(FICA_NONLIN_GAUSS+1) : {
mat Y = transpose(X) * B;
mat ex = exp(-a2 * pow(Y, 2) / 2);
mat gauss = elem_mult(Y, ex);
vec Beta = sumcol(elem_mult(Y, gauss));
vec Beta2 = sumcol(elem_mult(1 - a2 * pow(Y, 2), ex));
mat D = diag(pow(Beta - Beta2 , -1));
B = B + myy * B * (transpose(Y) * gauss - diag(Beta)) * D;
break;
}
case(FICA_NONLIN_GAUSS+2) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
mat U = transpose(Xsub) * B;
mat Usquared = pow(U, 2);
mat ex = exp(-a2 * Usquared / 2);
mat gauss = elem_mult(U, ex);
mat dGauss = elem_mult(1 - a2 * Usquared, ex);
B = (Xsub * gauss) / Xsub.cols() - elem_mult(reshape(repeat(sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / Xsub.cols();
break;
}
case(FICA_NONLIN_GAUSS+3) : {
mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
mat ex = exp(-a2 * pow(Ysub, 2) / 2);
mat gauss = elem_mult(Ysub, ex);
vec Beta = sumcol(elem_mult(Ysub, gauss));
vec Beta2 = sumcol(elem_mult(1 - a2 * pow(Ysub, 2), ex));
mat D = diag(pow(Beta - Beta2 , -1));
B = B + myy * B * (transpose(Ysub) * gauss - diag(Beta)) * D;
break;
}
// SKEW
case FICA_NONLIN_SKEW : {
B = (X * (pow(transpose(X) * B, 2))) / numSamples;
break;
}
case(FICA_NONLIN_SKEW+1) : {
mat Y = transpose(X) * B;
mat Gskew = pow(Y, 2);
vec Beta = sumcol(elem_mult(Y, Gskew));
mat D = diag(pow(Beta , -1));
B = B + myy * B * (transpose(Y) * Gskew - diag(Beta)) * D;
break;
}
case(FICA_NONLIN_SKEW+2) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
B = (Xsub * (pow(transpose(Xsub) * B, 2))) / Xsub.cols();
break;
}
case(FICA_NONLIN_SKEW+3) : {
mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
mat Gskew = pow(Ysub, 2);
vec Beta = sumcol(elem_mult(Ysub, Gskew));
mat D = diag(pow(Beta , -1));
B = B + myy * B * (transpose(Ysub) * Gskew - diag(Beta)) * D;
break;
}
} // SWITCH usedNlinearity
} // FOR maxIterations
W = transpose(B) * whiteningMatrix;
} // IF FICA_APPROACH_SYMM APPROACH
// DEFLATION
else {
// FC 01/12/05
A = zeros(whiteningMatrix.cols(), numOfIC);
// A = zeros( vectorSize, numOfIC );
mat B = zeros(vectorSize, numOfIC);
W = transpose(B) * whiteningMatrix;
int round = 1;
int numFailures = 0;
while (round <= numOfIC) {
myy = myyOrig;
usedNlinearity = gOrig;
stroke = 0;
notFine = 1;
loong = 0;
int endFinetuning = 0;
vec w = zeros(vectorSize);
if (initialStateMode == 0)
w = randu(vectorSize) - 0.5;
else w = whiteningMatrix * guess.get_col(round);
w = w - B * transpose(B) * w;
w /= norm(w);
vec wOld = zeros(vectorSize);
vec wOld2 = zeros(vectorSize);
int i = 1;
int gabba = 1;
while (i <= maxNumIterations + gabba) {
w = w - B * transpose(B) * w;
w /= norm(w);
if (notFine) {
if (i == maxNumIterations + 1) {
round--;
numFailures++;
if (numFailures > failureLimit) {
if (round == 0) {
A = dewhiteningMatrix * B;
W = transpose(B) * whiteningMatrix;
} // IF round
return false;
} // IF numFailures > failureLimit
break;
} // IF i == maxNumIterations+1
} // IF NOTFINE
else if (i >= endFinetuning) wOld = w;
if (norm(w - wOld) < epsilon || norm(w + wOld) < epsilon) {
if (fineTuningEnabled && notFine) {
notFine = 0;
gabba = maxFinetune;
wOld = zeros(vectorSize);
wOld2 = zeros(vectorSize);
usedNlinearity = gFine;
myy = myyK * myyOrig;
endFinetuning = maxFinetune + i;
} // IF finetuning
else {
numFailures = 0;
B.set_col(round - 1, w);
A.set_col(round - 1, dewhiteningMatrix*w);
W.set_row(round - 1, transpose(whiteningMatrix)*w);
break;
} // ELSE finetuning
} // IF epsilon
else if (stabilizationEnabled) {
if (stroke == 0.0 && (norm(w - wOld2) < epsilon || norm(w + wOld2) < epsilon)) {
stroke = myy;
myy /= 2.0 ;
if (mod(usedNlinearity, 2) == 0) {
usedNlinearity++;
} // IF MOD
}// IF !stroke
else if (stroke != 0.0) {
myy = stroke;
stroke = 0.0;
if (myy == 1 && mod(usedNlinearity, 2) != 0) {
usedNlinearity--;
}
} // IF Stroke
else if (notFine && !loong && i > maxNumIterations / 2) {
loong = 1;
myy /= 2.0;
if (mod(usedNlinearity, 2) == 0) {
usedNlinearity++;
} // IF MOD
} // IF notFine
} // IF stabilization
wOld2 = wOld;
wOld = w;
switch (usedNlinearity) {
// pow3
case FICA_NONLIN_POW3 : {
w = (X * pow(transpose(X) * w, 3)) / numSamples - 3 * w;
break;
}
case(FICA_NONLIN_POW3+1) : {
vec Y = transpose(X) * w;
vec Gpow3 = X * pow(Y, 3) / numSamples;
double Beta = dot(w, Gpow3);
w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
break;
}
case(FICA_NONLIN_POW3+2) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
w = (Xsub * pow(transpose(Xsub) * w, 3)) / Xsub.cols() - 3 * w;
break;
}
case(FICA_NONLIN_POW3+3) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
vec Gpow3 = Xsub * pow(transpose(Xsub) * w, 3) / (Xsub.cols());
double Beta = dot(w, Gpow3);
w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
break;
}
// TANH
case FICA_NONLIN_TANH : {
vec hypTan = tanh(a1 * transpose(X) * w);
w = (X * hypTan - a1 * sum(1 - pow(hypTan, 2)) * w) / numSamples;
break;
}
case(FICA_NONLIN_TANH+1) : {
vec Y = transpose(X) * w;
vec hypTan = tanh(a1 * Y);
double Beta = dot(w, X * hypTan);
w = w - myy * ((X * hypTan - Beta * w) / (a1 * sum(1 - pow(hypTan, 2)) - Beta));
break;
}
case(FICA_NONLIN_TANH+2) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
vec hypTan = tanh(a1 * transpose(Xsub) * w);
w = (Xsub * hypTan - a1 * sum(1 - pow(hypTan, 2)) * w) / Xsub.cols();
break;
}
case(FICA_NONLIN_TANH+3) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
vec hypTan = tanh(a1 * transpose(Xsub) * w);
double Beta = dot(w, Xsub * hypTan);
w = w - myy * ((Xsub * hypTan - Beta * w) / (a1 * sum(1 - pow(hypTan, 2)) - Beta));
break;
}
// GAUSS
case FICA_NONLIN_GAUSS : {
vec u = transpose(X) * w;
vec Usquared = pow(u, 2);
vec ex = exp(-a2 * Usquared / 2);
vec gauss = elem_mult(u, ex);
vec dGauss = elem_mult(1 - a2 * Usquared, ex);
w = (X * gauss - sum(dGauss) * w) / numSamples;
break;
}
case(FICA_NONLIN_GAUSS+1) : {
vec u = transpose(X) * w;
vec Usquared = pow(u, 2);
vec ex = exp(-a2 * Usquared / 2);
vec gauss = elem_mult(u, ex);
vec dGauss = elem_mult(1 - a2 * Usquared, ex);
double Beta = dot(w, X * gauss);
w = w - myy * ((X * gauss - Beta * w) / (sum(dGauss) - Beta));
break;
}
case(FICA_NONLIN_GAUSS+2) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
vec u = transpose(Xsub) * w;
vec Usquared = pow(u, 2);
vec ex = exp(-a2 * Usquared / 2);
vec gauss = elem_mult(u, ex);
vec dGauss = elem_mult(1 - a2 * Usquared, ex);
w = (Xsub * gauss - sum(dGauss) * w) / Xsub.cols();
break;
}
case(FICA_NONLIN_GAUSS+3) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
vec u = transpose(Xsub) * w;
vec Usquared = pow(u, 2);
vec ex = exp(-a2 * Usquared / 2);
vec gauss = elem_mult(u, ex);
vec dGauss = elem_mult(1 - a2 * Usquared, ex);
double Beta = dot(w, Xsub * gauss);
w = w - myy * ((Xsub * gauss - Beta * w) / (sum(dGauss) - Beta));
break;
}
// SKEW
case FICA_NONLIN_SKEW : {
w = (X * (pow(transpose(X) * w, 2))) / numSamples;
break;
}
case(FICA_NONLIN_SKEW+1) : {
vec Y = transpose(X) * w;
vec Gskew = X * pow(Y, 2) / numSamples;
double Beta = dot(w, Gskew);
w = w - myy * (Gskew - Beta * w / (-Beta));
break;
}
case(FICA_NONLIN_SKEW+2) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
w = (Xsub * (pow(transpose(Xsub) * w, 2))) / Xsub.cols();
break;
}
case(FICA_NONLIN_SKEW+3) : {
mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
vec Gskew = Xsub * pow(transpose(Xsub) * w, 2) / Xsub.cols();
double Beta = dot(w, Gskew);
w = w - myy * (Gskew - Beta * w) / (-Beta);
break;
}
} // SWITCH nonLinearity
w /= norm(w);
i++;
} // WHILE i<= maxNumIterations+gabba
round++;
} // While round <= numOfIC
} // ELSE Deflation
return true;
} // FPICA
itpp::Fast_ICA::set_sample_size
void set_sample_size(double fl_sampleSize)
Set sample size.
Definition: fastica.cpp:198
fastica.h
Definition of FastICA (Independent Component Analysis) for IT++.
matfunc.h
Various functions on vectors and matrices - header file.
FICA_APPROACH_SYMM
#define FICA_APPROACH_SYMM
Use symmetric approach : compute all ICs at a time.
Definition: fastica.h:71
itpp::Fast_ICA::set_a1
void set_a1(double fl_a1)
Set parameter.
Definition: fastica.cpp:190
itpp::Mat::mat
Mat< double > mat
Default Matrix Type.
Definition: mat.h:482
itpp::concat_horizontal
Mat< Num_T > concat_horizontal(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Horizontal concatenation of two matrices.
Definition: mat.h:1194
itpp::randu
double randu(void)
Generates a random uniform (0,1) number.
Definition: random.h:804
itpp::Fast_ICA::set_mu
void set_mu(double fl_mu)
Set parameter.
Definition: fastica.cpp:194
remmean
static void remmean(mat inVectors, mat &outVectors, vec &meanValue)
Local functions for FastICA.
Definition: fastica.cpp:335
eigen.h
Definitions of eigenvalue decomposition functions.
itpp::Fast_ICA::separate
bool separate(void)
Explicit launch of main FastICA function.
Definition: fastica.cpp:126
FICA_NONLIN_SKEW
#define FICA_NONLIN_SKEW
Use skew non-linearity.
Definition: fastica.h:80
itpp::Fast_ICA::get_whitening_matrix
mat get_whitening_matrix()
Get the whitening matrix.
Definition: fastica.cpp:228
itpp::Fast_ICA::set_init_guess
void set_init_guess(mat ma_initGuess)
Set initial guess matrix instead of random (default)
Definition: fastica.cpp:212
itpp::Vec::cvec
Vec< std::complex< double > > cvec
Definition of complex<double> vector type.
Definition: vec.h:535
itpp::norm
double norm(const cvec &v)
Calculate the 2-norm: norm(v)=sqrt(sum(abs(v).^2))
Definition: misc_stat.cpp:77
itpp::Sparse_Vec::sparse_vec
Sparse_Vec< double > sparse_vec
Type definition of a double sparse vector.
Definition: svec.h:347
itpp::elem_mult
Mat< Num_T > elem_mult(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise multiplication of two matrices.
Definition: mat.h:1582
pcamat
static int pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat &Es, vec &Ds)
Local functions for FastICA.
Definition: fastica.cpp:262
itpp::mod
int mod(int k, int n)
Calculates the modulus, i.e. the signed reminder after division.
Definition: elem_math.h:166
itpp::Fast_ICA::set_max_num_iterations
void set_max_num_iterations(int in_maxNumIterations)
Set maximum number of iterations.
Definition: fastica.cpp:202
itpp::sum
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition: matfunc.h:59
itpp::Fast_ICA::set_epsilon
void set_epsilon(double fl_epsilon)
Set convergence parameter .
Definition: fastica.cpp:196
FICA_NONLIN_TANH
#define FICA_NONLIN_TANH
Use tanh(x) non-linearity.
Definition: fastica.h:76
itpp::Fast_ICA::set_nrof_independent_components
void set_nrof_independent_components(int in_nrIC)
Set number of independent components to separate.
Definition: fastica.cpp:184
getSamples
static ivec getSamples(const int max, const double percentage)
Local functions for FastICA.
Definition: fastica.cpp:417
FICA_NONLIN_GAUSS
#define FICA_NONLIN_GAUSS
Use Gaussian non-linearity.
Definition: fastica.h:78
sort.h
Sorting functions.
min_max.h
Minimum and maximum functions on vectors and matrices.
sigfun.h
Definitions of signal processing functions.
itpp::mean
double mean(const vec &v)
The mean value.
Definition: misc_stat.cpp:36
itpp::round
ITPP_EXPORT double round(double x)
Round to nearest integer, return result in double.
itpp::Vec::sort
void sort(Vec< T > &data, SORTING_METHOD method=INTROSORT)
Sort the data vector in increasing order.
Definition: sort.h:187
FICA_APPROACH_DEFL
#define FICA_APPROACH_DEFL
Use deflation approach : compute IC one-by-one in a Gram-Schmidt-like fashion.
Definition: fastica.h:69
itpp::Fast_ICA::get_nrof_independent_components
int get_nrof_independent_components()
Get number of independent components.
Definition: fastica.cpp:224
itpp::zeros
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
itpp::eps
const double eps
Constant eps.
Definition: misc.h:109
itpp::min
T min(const Vec< T > &in)
Minimum value of vector.
Definition: min_max.h:125
itpp::Fast_ICA::set_last_eig
void set_last_eig(int in_lastEig)
Set last eigenvalue index to take into account.
Definition: fastica.cpp:208
itpp::max_index
int max_index(const Vec< T > &in)
Return the postion of the maximum element in the vector.
Definition: min_max.h:208
itpp::Fast_ICA::set_pca_only
void set_pca_only(bool in_PCAonly)
If true, only perform Principal Component Analysis (default = false)
Definition: fastica.cpp:210
itpp::svd
bool svd(const mat &A, vec &S)
Get singular values s of a real matrix A using SVD.
Definition: svd.cpp:185
itpp::exp
vec exp(const vec &x)
Exp of the elements of a vector x.
Definition: log_exp.h:155
itpp::Mat::cmat
Mat< std::complex< double > > cmat
Default Complex Matrix Type.
Definition: mat.h:488
itpp::max
T max(const Vec< T > &v)
Maximum value of vector.
Definition: min_max.h:45
trig_hyp.h
Trigonometric and hyperbolic functions - header file.
itpp::pow
vec pow(const double x, const vec &y)
Calculates x to the power of y (x^y)
Definition: log_exp.h:176
itpp::cov
mat cov(const mat &X, bool is_zero_mean)
Covariance matrix calculation.
Definition: sigfun.cpp:226
itpp::bin::abs
bin abs(const bin &inbin)
absolute value of bin
Definition: binary.h:174
FICA_NONLIN_POW3
#define FICA_NONLIN_POW3
Use x^3 non-linearity.
Definition: fastica.h:74
itpp::size
int size(const Vec< T > &v)
Length of vector.
Definition: matfunc.h:55
itpp::Fast_ICA::get_separating_matrix
mat get_separating_matrix()
Get separating matrix.
Definition: fastica.cpp:220
sumcol
static vec sumcol(const mat A)
Local functions for FastICA.
Definition: fastica.cpp:433
itpp::full
Mat< T > full(const Sparse_Mat< T > &s)
Convert a sparse matrix s into its dense representation.
Definition: smat.h:998
itpp::Fast_ICA::set_stabilization
void set_stabilization(bool in_stabilization)
Set stabilization mode true or off.
Definition: fastica.cpp:200
mpower
static mat mpower(const mat A, const double y)
Local functions for FastICA.
Definition: fastica.cpp:397
itpp::Fast_ICA::set_max_fine_tune
void set_max_fine_tune(int in_maxFineTune)
Set maximum number of iterations for fine tuning.
Definition: fastica.cpp:204
itpp::Fast_ICA::set_non_linearity
void set_non_linearity(int in_g)
Set non-linearity.
Definition: fastica.cpp:186
resampling.h
Resampling functions - header file.
misc_stat.h
Miscellaneous statistics functions and classes - header file.
itpp::to_ivec
ivec to_ivec(const Vec< T > &v)
Converts a Vec<T> to ivec.
Definition: converters.h:79
FICA_INIT_GUESS
#define FICA_INIT_GUESS
Set predefined start for Fast_ICA.
Definition: fastica.h:85
specmat.h
Definitions of special vectors and matrices.
it_warning
#define it_warning(s)
Display a warning message.
Definition: itassert.h:173
fpica
static bool fpica(const mat X, const mat whiteningMatrix, const mat dewhiteningMatrix, const int approach, const int numOfIC, const int g, const int finetune, const double a1, const double a2, double myy, const int stabilization, const double epsilon, const int maxNumIterations, const int maxFinetune, const int initState, mat guess, double sampleSize, mat &A, mat &W)
Local functions for FastICA.
Definition: fastica.cpp:444
itpp::Fast_ICA::get_principal_eigenvectors
mat get_principal_eigenvectors()
Get nrIC first columns of the de-whitening matrix.
Definition: fastica.cpp:226
itpp::Fast_ICA::get_independent_components
mat get_independent_components()
Get separated signals.
Definition: fastica.cpp:222
itpp::eig_sym
bool eig_sym(const mat &A, vec &d, mat &V)
Calculates the eigenvalues and eigenvectors of a symmetric real matrix.
Definition: eigen.cpp:252
whitenv
static void whitenv(const mat vectors, const mat E, const mat D, mat &newVectors, mat &whiteningMatrix, mat &dewhiteningMatrix)
Local functions for FastICA.
Definition: fastica.cpp:351
itpp::transpose
void transpose(const Mat< T > &m, Mat< T > &out)
Transposition of the matrix m returning the transposed matrix in out.
Definition: matfunc.h:308
FICA_INIT_RAND
#define FICA_INIT_RAND
Set random start for Fast_ICA.
Definition: fastica.h:83
itpp::reshape
Mat< T > reshape(const Mat< T > &m, int rows, int cols)
Reshape the matrix into an rows*cols matrix.
Definition: matfunc.h:822
itpp::Vec::ivec
Vec< int > ivec
Definition of integer vector type.
Definition: vec.h:541
itpp::Fast_ICA::get_dewhitening_matrix
mat get_dewhitening_matrix()
Get the de-whitening matrix.
Definition: fastica.cpp:230
itpp::Fast_ICA::set_a2
void set_a2(double fl_a2)
Set parameter.
Definition: fastica.cpp:192
itpp::tanh
vec tanh(const vec &x)
Tan hyperbolic function.
Definition: trig_hyp.h:97
itpp::sqrt
vec sqrt(const vec &x)
Square root of the elements.
Definition: elem_math.h:123
svd.h
Definitions of Singular Value Decompositions.
itpp::diag
Mat< T > diag(const Vec< T > &v, const int K=0)
Create a diagonal matrix using vector v as its diagonal.
Definition: matfunc.h:557
orth
static mat orth(const mat A)
Local functions for FastICA.
Definition: fastica.cpp:368
itpp::Fast_ICA::get_white_sig
mat get_white_sig()
Get whitened signals.
Definition: fastica.cpp:232
random.h
Definition of classes for random number generators.
itpp::Fast_ICA::Fast_ICA
Fast_ICA(mat ma_mixed_sig)
Constructor.
Definition: fastica.cpp:99
itpp::Fast_ICA::set_first_eig
void set_first_eig(int in_firstEig)
Set first eigenvalue index to take into account.
Definition: fastica.cpp:206
svec.h
Sparse Vector Class definitions.
itpp::dot
Num_T dot(const Vec< Num_T > &v1, const Vec< Num_T > &v2)
Inner (dot) product of two vectors v1 and v2.
Definition: vec.h:1005
FICA_TOL
#define FICA_TOL
Eigenvalues of the covariance matrix lower than FICA_TOL are discarded for analysis.
Definition: fastica.h:88
itpp::Fast_ICA::get_mixing_matrix
mat get_mixing_matrix()
Get mixing matrix.
Definition: fastica.cpp:218
itpp::Fast_ICA::set_fine_tune
void set_fine_tune(bool in_finetune)
Set fine tuning.
Definition: fastica.cpp:188
itpp::Fast_ICA::set_approach
void set_approach(int in_approach)
Set approach : FICA_APPROACH_DEFL or FICA_APPROACH_SYMM (default)
Definition: fastica.cpp:182
itpp::Vec::vec
Vec< double > vec
Definition of double vector type.
Definition: vec.h:529
selcol
static void selcol(const mat oldMatrix, const vec maskVector, mat &newMatrix)
Local functions for FastICA.
Definition: fastica.cpp:237

Generated on Thu Mar 20 2014 12:26:43 for IT++ by Doxygen 1.8.6