IT++ Logo
fastica.cpp
Go to the documentation of this file.
1 
62 #include <itpp/signal/fastica.h>
63 #include <itpp/signal/sigfun.h>
64 #include <itpp/signal/resampling.h>
66 #include <itpp/base/algebra/svd.h>
68 #include <itpp/base/matfunc.h>
69 #include <itpp/base/random.h>
70 #include <itpp/base/sort.h>
71 #include <itpp/base/specmat.h>
72 #include <itpp/base/svec.h>
73 #include <itpp/base/math/min_max.h>
74 #include <itpp/stat/misc_stat.h>
75 
76 
77 using namespace itpp;
78 
79 
84 static void selcol(const mat oldMatrix, const vec maskVector, mat & newMatrix);
85 static int pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat & Es, vec & Ds);
86 static void remmean(mat inVectors, mat & outVectors, vec & meanValue);
87 static void whitenv(const mat vectors, const mat E, const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix);
88 static mat orth(const mat A);
89 static mat mpower(const mat A, const double y);
90 static ivec getSamples(const int max, const double percentage);
91 static vec sumcol(const mat A);
92 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);
95 namespace itpp
96 {
97 
98 // Constructor, init default values
99 Fast_ICA::Fast_ICA(mat ma_mixedSig)
100 {
101 
102  // Init default values
103  approach = FICA_APPROACH_SYMM;
104  g = FICA_NONLIN_POW3;
105  finetune = true;
106  a1 = 1.0;
107  a2 = 1.0;
108  mu = 1.0;
109  epsilon = 0.0001;
110  sampleSize = 1.0;
111  stabilization = false;
112  maxNumIterations = 100000;
113  maxFineTune = 100;
114  firstEig = 1;
115 
116  mixedSig = ma_mixedSig;
117 
118  lastEig = mixedSig.rows();
119  numOfIC = mixedSig.rows();
120  PCAonly = false;
121  initState = FICA_INIT_RAND;
122 
123 }
124 
125 // Call main function
127 {
128 
129  int Dim = numOfIC;
130 
131  mat mixedSigC;
132  vec mixedMean;
133 
134  mat guess;
135  if (initState == FICA_INIT_RAND)
136  guess = zeros(Dim, Dim);
137  else
138  guess = mat(initGuess);
139 
140  VecPr = zeros(mixedSig.rows(), numOfIC);
141 
142  icasig = zeros(numOfIC, mixedSig.cols());
143 
144  remmean(mixedSig, mixedSigC, mixedMean);
145 
146  if (pcamat(mixedSigC, numOfIC, firstEig, lastEig, E, D) < 1) {
147  // no principal components could be found (e.g. all-zero data): return the unchanged input
148  icasig = mixedSig;
149  return false;
150  }
151 
152  whitenv(mixedSigC, E, diag(D), whitesig, whiteningMatrix, dewhiteningMatrix);
153 
154  Dim = whitesig.rows();
155  if (numOfIC > Dim) numOfIC = Dim;
156 
157  ivec NcFirst = to_ivec(zeros(numOfIC));
158  vec NcVp = D;
159  for (int i = 0; i < NcFirst.size(); i++) {
160 
161  NcFirst(i) = max_index(NcVp);
162  NcVp(NcFirst(i)) = 0.0;
163  VecPr.set_col(i, dewhiteningMatrix.get_col(i));
164 
165  }
166 
167  bool result = true;
168  if (PCAonly == false) {
169 
170  result = fpica(whitesig, whiteningMatrix, dewhiteningMatrix, approach, numOfIC, g, finetune, a1, a2, mu, stabilization, epsilon, maxNumIterations, maxFineTune, initState, guess, sampleSize, A, W);
171 
172  icasig = W * mixedSig;
173 
174  }
175 
176  else { // PCA only : returns E as IcaSig
177  icasig = VecPr;
178  }
179  return result;
180 }
181 
182 void Fast_ICA::set_approach(int in_approach) { approach = in_approach; if (approach == FICA_APPROACH_DEFL) finetune = true; }
183 
184 void Fast_ICA::set_nrof_independent_components(int in_nrIC) { numOfIC = in_nrIC; }
185 
186 void Fast_ICA::set_non_linearity(int in_g) { g = in_g; }
187 
188 void Fast_ICA::set_fine_tune(bool in_finetune) { finetune = in_finetune; }
189 
190 void Fast_ICA::set_a1(double fl_a1) { a1 = fl_a1; }
191 
192 void Fast_ICA::set_a2(double fl_a2) { a2 = fl_a2; }
193 
194 void Fast_ICA::set_mu(double fl_mu) { mu = fl_mu; }
195 
196 void Fast_ICA::set_epsilon(double fl_epsilon) { epsilon = fl_epsilon; }
197 
198 void Fast_ICA::set_sample_size(double fl_sampleSize) { sampleSize = fl_sampleSize; }
199 
200 void Fast_ICA::set_stabilization(bool in_stabilization) { stabilization = in_stabilization; }
201 
202 void Fast_ICA::set_max_num_iterations(int in_maxNumIterations) { maxNumIterations = in_maxNumIterations; }
203 
204 void Fast_ICA::set_max_fine_tune(int in_maxFineTune) { maxFineTune = in_maxFineTune; }
205 
206 void Fast_ICA::set_first_eig(int in_firstEig) { firstEig = in_firstEig; }
207 
208 void Fast_ICA::set_last_eig(int in_lastEig) { lastEig = in_lastEig; }
209 
210 void Fast_ICA::set_pca_only(bool in_PCAonly) { PCAonly = in_PCAonly; }
211 
212 void Fast_ICA::set_init_guess(mat ma_initGuess)
213 {
214  initGuess = ma_initGuess;
215  initState = FICA_INIT_GUESS;
216 }
217 
218 mat Fast_ICA::get_mixing_matrix() { if (PCAonly) { it_warning("No ICA performed."); return (zeros(1, 1));} else return A; }
219 
220 mat Fast_ICA::get_separating_matrix() { if (PCAonly) { it_warning("No ICA performed."); return(zeros(1, 1)); } else return W; }
221 
222 mat Fast_ICA::get_independent_components() { if (PCAonly) { it_warning("No ICA performed."); return(zeros(1, 1)); } else return icasig; }
223 
225 
227 
228 mat Fast_ICA::get_whitening_matrix() { return whiteningMatrix; }
229 
230 mat Fast_ICA::get_dewhitening_matrix() { return dewhiteningMatrix; }
231 
232 mat Fast_ICA::get_white_sig() { return whitesig; }
233 
234 } // namespace itpp
235 
236 
237 static void selcol(const mat oldMatrix, const vec maskVector, mat & newMatrix)
238 {
239 
240  int numTaken = 0;
241 
242  for (int i = 0; i < size(maskVector); i++) if (maskVector(i) == 1) numTaken++;
243 
244  newMatrix = zeros(oldMatrix.rows(), numTaken);
245 
246  numTaken = 0;
247 
248  for (int i = 0; i < size(maskVector); i++) {
249 
250  if (maskVector(i) == 1) {
251 
252  newMatrix.set_col(numTaken, oldMatrix.get_col(i));
253  numTaken++;
254 
255  }
256  }
257 
258  return;
259 
260 }
261 
262 static int pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat & Es, vec & Ds)
263 {
264 
265  mat Et;
266  vec Dt;
267  cmat Ec;
268  cvec Dc;
269  double lowerLimitValue = 0.0,
270  higherLimitValue = 0.0;
271 
272  int oldDimension = vectors.rows();
273 
274  mat covarianceMatrix = cov(transpose(vectors), 0);
275 
276  eig_sym(covarianceMatrix, Dt, Et);
277 
278  int maxLastEig = 0;
279 
280  // Compute rank
281  for (int i = 0; i < Dt.length(); i++) if (Dt(i) > FICA_TOL) maxLastEig++;
282  if (maxLastEig < 1) return 0;
283 
284  // Force numOfIC components
285  if (maxLastEig > numOfIC) maxLastEig = numOfIC;
286 
287  vec eigenvalues = zeros(size(Dt));
288  vec eigenvalues2 = zeros(size(Dt));
289 
290  eigenvalues2 = Dt;
291 
292  sort(eigenvalues2);
293 
294  vec lowerColumns = zeros(size(Dt));
295 
296  for (int i = 0; i < size(Dt); i++) eigenvalues(i) = eigenvalues2(size(Dt) - i - 1);
297 
298  if (lastEig > maxLastEig) lastEig = maxLastEig;
299 
300  if (lastEig < oldDimension) lowerLimitValue = (eigenvalues(lastEig - 1) + eigenvalues(lastEig)) / 2;
301  else lowerLimitValue = eigenvalues(oldDimension - 1) - 1;
302 
303  for (int i = 0; i < size(Dt); i++) if (Dt(i) > lowerLimitValue) lowerColumns(i) = 1;
304 
305  if (firstEig > 1) higherLimitValue = (eigenvalues(firstEig - 2) + eigenvalues(firstEig - 1)) / 2;
306  else higherLimitValue = eigenvalues(0) + 1;
307 
308  vec higherColumns = zeros(size(Dt));
309  for (int i = 0; i < size(Dt); i++) if (Dt(i) < higherLimitValue) higherColumns(i) = 1;
310 
311  vec selectedColumns = zeros(size(Dt));
312  for (int i = 0; i < size(Dt); i++) selectedColumns(i) = (lowerColumns(i) == 1 && higherColumns(i) == 1) ? 1 : 0;
313 
314  selcol(Et, selectedColumns, Es);
315 
316  int numTaken = 0;
317 
318  for (int i = 0; i < size(selectedColumns); i++) if (selectedColumns(i) == 1) numTaken++;
319 
320  Ds = zeros(numTaken);
321 
322  numTaken = 0;
323 
324  for (int i = 0; i < size(Dt); i++)
325  if (selectedColumns(i) == 1) {
326  Ds(numTaken) = Dt(i);
327  numTaken++;
328  }
329 
330  return lastEig;
331 
332 }
333 
334 
335 static void remmean(mat inVectors, mat & outVectors, vec & meanValue)
336 {
337 
338  outVectors = zeros(inVectors.rows(), inVectors.cols());
339  meanValue = zeros(inVectors.rows());
340 
341  for (int i = 0; i < inVectors.rows(); i++) {
342 
343  meanValue(i) = mean(inVectors.get_row(i));
344 
345  for (int j = 0; j < inVectors.cols(); j++) outVectors(i, j) = inVectors(i, j) - meanValue(i);
346 
347  }
348 
349 }
350 
351 static void whitenv(const mat vectors, const mat E, const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix)
352 {
353 
354  whiteningMatrix = zeros(E.cols(), E.rows());
355  dewhiteningMatrix = zeros(E.rows(), E.cols());
356 
357  for (int i = 0; i < D.cols(); i++) {
358  whiteningMatrix.set_row(i, std::pow(std::sqrt(D(i, i)), -1)*E.get_col(i));
359  dewhiteningMatrix.set_col(i, std::sqrt(D(i, i))*E.get_col(i));
360  }
361 
362  newVectors = whiteningMatrix * vectors;
363 
364  return;
365 
366 }
367 
368 static mat orth(const mat A)
369 {
370 
371  mat Q;
372  mat U, V;
373  vec S;
374  double eps = 2.2e-16;
375  double tol = 0.0;
376  int mmax = 0;
377  int r = 0;
378 
379  svd(A, U, S, V);
380  if (A.rows() > A.cols()) {
381 
382  U = U(0, U.rows() - 1, 0, A.cols() - 1);
383  S = S(0, A.cols() - 1);
384  }
385 
386  mmax = (A.rows() > A.cols()) ? A.rows() : A.cols();
387 
388  tol = mmax * eps * max(S);
389 
390  for (int i = 0; i < size(S); i++) if (S(i) > tol) r++;
391 
392  Q = U(0, U.rows() - 1, 0, r - 1);
393 
394  return (Q);
395 }
396 
397 static mat mpower(const mat A, const double y)
398 {
399 
400  mat T = zeros(A.rows(), A.cols());
401  mat dd = zeros(A.rows(), A.cols());
402  vec d = zeros(A.rows());
403  vec dOut = zeros(A.rows());
404 
405  eig_sym(A, d, T);
406 
407  dOut = pow(d, y);
408 
409  diag(dOut, dd);
410 
411  for (int i = 0; i < T.cols(); i++) T.set_col(i, T.get_col(i) / norm(T.get_col(i)));
412 
413  return (T*dd*transpose(T));
414 
415 }
416 
417 static ivec getSamples(const int max, const double percentage)
418 {
419 
420  vec rd = randu(max);
421  sparse_vec sV;
422  ivec out;
423  int sZ = 0;
424 
425  for (int i = 0; i < max; i++) if (rd(i) < percentage) { sV.add_elem(sZ, i); sZ++; }
426 
427  out = to_ivec(full(sV));
428 
429  return (out);
430 
431 }
432 
433 static vec sumcol(const mat A)
434 {
435 
436  vec out = zeros(A.cols());
437 
438  for (int i = 0; i < A.cols(); i++) { out(i) = sum(A.get_col(i)); }
439 
440  return (out);
441 
442 }
443 
444 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)
445 {
446 
447  int vectorSize = X.rows();
448  int numSamples = X.cols();
449  int gOrig = g;
450  int gFine = finetune + 1;
451  double myyOrig = myy;
452  double myyK = 0.01;
453  int failureLimit = 5;
454  int usedNlinearity = 0;
455  double stroke = 0.0;
456  int notFine = 1;
457  int loong = 0;
458  int initialStateMode = initState;
459  double minAbsCos = 0.0, minAbsCos2 = 0.0;
460 
461  if (sampleSize * numSamples < 1000) sampleSize = (1000 / (double)numSamples < 1.0) ? 1000 / (double)numSamples : 1.0;
462 
463  if (sampleSize != 1.0) gOrig += 2;
464  if (myy != 1.0) gOrig += 1;
465 
466  int fineTuningEnabled = 1;
467 
468  if (!finetune) {
469  if (myy != 1.0) gFine = gOrig;
470  else gFine = gOrig + 1;
471  fineTuningEnabled = 0;
472  }
473 
474  int stabilizationEnabled = stabilization;
475 
476  if (!stabilization && myy != 1.0) stabilizationEnabled = true;
477 
478  usedNlinearity = gOrig;
479 
480  if (initState == FICA_INIT_GUESS && guess.rows() != whiteningMatrix.cols()) {
481  initialStateMode = 0;
482 
483  }
484  else if (guess.cols() < numOfIC) {
485 
486  mat guess2 = randu(guess.rows(), numOfIC - guess.cols()) - 0.5;
487  guess = concat_horizontal(guess, guess2);
488  }
489  else if (guess.cols() > numOfIC) guess = guess(0, guess.rows() - 1, 0, numOfIC - 1);
490 
491  if (approach == FICA_APPROACH_SYMM) {
492 
493  usedNlinearity = gOrig;
494  stroke = 0;
495  notFine = 1;
496  loong = 0;
497 
498  A = zeros(vectorSize, numOfIC);
499  mat B = zeros(vectorSize, numOfIC);
500 
501  if (initialStateMode == 0) B = orth(randu(vectorSize, numOfIC) - 0.5);
502  else B = whiteningMatrix * guess;
503 
504  mat BOld = zeros(B.rows(), B.cols());
505  mat BOld2 = zeros(B.rows(), B.cols());
506 
507  for (int round = 0; round < maxNumIterations; round++) {
508 
509  if (round == maxNumIterations - 1) {
510 
511  // If there is a convergence problem,
512  // we still want ot return something.
513  // This is difference with original
514  // Matlab implementation.
515  A = dewhiteningMatrix * B;
516  W = transpose(B) * whiteningMatrix;
517 
518  return false;
519  }
520 
521  B = B * mpower(transpose(B) * B , -0.5);
522 
523  minAbsCos = min(abs(diag(transpose(B) * BOld)));
524  minAbsCos2 = min(abs(diag(transpose(B) * BOld2)));
525 
526  if (1 - minAbsCos < epsilon) {
527 
528  if (fineTuningEnabled && notFine) {
529 
530  notFine = 0;
531  usedNlinearity = gFine;
532  myy = myyK * myyOrig;
533  BOld = zeros(B.rows(), B.cols());
534  BOld2 = zeros(B.rows(), B.cols());
535 
536  }
537 
538  else {
539 
540  A = dewhiteningMatrix * B;
541  break;
542 
543  }
544 
545  } // IF epsilon
546 
547  else if (stabilizationEnabled) {
548  if (!stroke && (1 - minAbsCos2 < epsilon)) {
549 
550  stroke = myy;
551  myy /= 2;
552  if (mod(usedNlinearity, 2) == 0) usedNlinearity += 1 ;
553 
554  }
555  else if (stroke) {
556 
557  myy = stroke;
558  stroke = 0;
559  if (myy == 1 && mod(usedNlinearity, 2) != 0) usedNlinearity -= 1;
560 
561  }
562  else if (!loong && (round > maxNumIterations / 2)) {
563 
564  loong = 1;
565  myy /= 2;
566  if (mod(usedNlinearity, 2) == 0) usedNlinearity += 1;
567 
568  }
569 
570  } // stabilizationEnabled
571 
572  BOld2 = BOld;
573  BOld = B;
574 
575  switch (usedNlinearity) {
576 
577  // pow3
578  case FICA_NONLIN_POW3 : {
579  B = (X * pow(transpose(X) * B, 3)) / numSamples - 3 * B;
580  break;
581  }
582  case(FICA_NONLIN_POW3+1) : {
583  mat Y = transpose(X) * B;
584  mat Gpow3 = pow(Y, 3);
585  vec Beta = sumcol(pow(Y, 4));
586  mat D = diag(pow(Beta - 3 * numSamples , -1));
587  B = B + myy * B * (transpose(Y) * Gpow3 - diag(Beta)) * D;
588  break;
589  }
590  case(FICA_NONLIN_POW3+2) : {
591  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
592  B = (Xsub * pow(transpose(Xsub) * B, 3)) / Xsub.cols() - 3 * B;
593  break;
594  }
595  case(FICA_NONLIN_POW3+3) : {
596  mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
597  mat Gpow3 = pow(Ysub, 3);
598  vec Beta = sumcol(pow(Ysub, 4));
599  mat D = diag(pow(Beta - 3 * Ysub.rows() , -1));
600  B = B + myy * B * (transpose(Ysub) * Gpow3 - diag(Beta)) * D;
601  break;
602  }
603 
604  // TANH
605  case FICA_NONLIN_TANH : {
606  mat hypTan = tanh(a1 * transpose(X) * B);
607  B = (X * hypTan) / numSamples - elem_mult(reshape(repeat(sumcol(1 - pow(hypTan, 2)), B.rows()), B.rows(), B.cols()), B) / numSamples * a1;
608  break;
609  }
610  case(FICA_NONLIN_TANH+1) : {
611  mat Y = transpose(X) * B;
612  mat hypTan = tanh(a1 * Y);
613  vec Beta = sumcol(elem_mult(Y, hypTan));
614  vec Beta2 = sumcol(1 - pow(hypTan, 2));
615  mat D = diag(pow(Beta - a1 * Beta2 , -1));
616  B = B + myy * B * (transpose(Y) * hypTan - diag(Beta)) * D;
617  break;
618  }
619  case(FICA_NONLIN_TANH+2) : {
620  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
621  mat hypTan = tanh(a1 * transpose(Xsub) * B);
622  B = (Xsub * hypTan) / Xsub.cols() - elem_mult(reshape(repeat(sumcol(1 - pow(hypTan, 2)), B.rows()), B.rows(), B.cols()), B) / Xsub.cols() * a1;
623  break;
624  }
625  case(FICA_NONLIN_TANH+3) : {
626  mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
627  mat hypTan = tanh(a1 * Ysub);
628  vec Beta = sumcol(elem_mult(Ysub, hypTan));
629  vec Beta2 = sumcol(1 - pow(hypTan, 2));
630  mat D = diag(pow(Beta - a1 * Beta2 , -1));
631  B = B + myy * B * (transpose(Ysub) * hypTan - diag(Beta)) * D;
632  break;
633  }
634 
635  // GAUSS
636  case FICA_NONLIN_GAUSS : {
637  mat U = transpose(X) * B;
638  mat Usquared = pow(U, 2);
639  mat ex = exp(-a2 * Usquared / 2);
640  mat gauss = elem_mult(U, ex);
641  mat dGauss = elem_mult(1 - a2 * Usquared, ex);
642  B = (X * gauss) / numSamples - elem_mult(reshape(repeat(sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / numSamples;
643  break;
644  }
645  case(FICA_NONLIN_GAUSS+1) : {
646  mat Y = transpose(X) * B;
647  mat ex = exp(-a2 * pow(Y, 2) / 2);
648  mat gauss = elem_mult(Y, ex);
649  vec Beta = sumcol(elem_mult(Y, gauss));
650  vec Beta2 = sumcol(elem_mult(1 - a2 * pow(Y, 2), ex));
651  mat D = diag(pow(Beta - Beta2 , -1));
652  B = B + myy * B * (transpose(Y) * gauss - diag(Beta)) * D;
653  break;
654  }
655  case(FICA_NONLIN_GAUSS+2) : {
656  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
657  mat U = transpose(Xsub) * B;
658  mat Usquared = pow(U, 2);
659  mat ex = exp(-a2 * Usquared / 2);
660  mat gauss = elem_mult(U, ex);
661  mat dGauss = elem_mult(1 - a2 * Usquared, ex);
662  B = (Xsub * gauss) / Xsub.cols() - elem_mult(reshape(repeat(sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / Xsub.cols();
663  break;
664  }
665  case(FICA_NONLIN_GAUSS+3) : {
666  mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
667  mat ex = exp(-a2 * pow(Ysub, 2) / 2);
668  mat gauss = elem_mult(Ysub, ex);
669  vec Beta = sumcol(elem_mult(Ysub, gauss));
670  vec Beta2 = sumcol(elem_mult(1 - a2 * pow(Ysub, 2), ex));
671  mat D = diag(pow(Beta - Beta2 , -1));
672  B = B + myy * B * (transpose(Ysub) * gauss - diag(Beta)) * D;
673  break;
674  }
675 
676  // SKEW
677  case FICA_NONLIN_SKEW : {
678  B = (X * (pow(transpose(X) * B, 2))) / numSamples;
679  break;
680  }
681  case(FICA_NONLIN_SKEW+1) : {
682  mat Y = transpose(X) * B;
683  mat Gskew = pow(Y, 2);
684  vec Beta = sumcol(elem_mult(Y, Gskew));
685  mat D = diag(pow(Beta , -1));
686  B = B + myy * B * (transpose(Y) * Gskew - diag(Beta)) * D;
687  break;
688  }
689  case(FICA_NONLIN_SKEW+2) : {
690  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
691  B = (Xsub * (pow(transpose(Xsub) * B, 2))) / Xsub.cols();
692  break;
693  }
694  case(FICA_NONLIN_SKEW+3) : {
695  mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
696  mat Gskew = pow(Ysub, 2);
697  vec Beta = sumcol(elem_mult(Ysub, Gskew));
698  mat D = diag(pow(Beta , -1));
699  B = B + myy * B * (transpose(Ysub) * Gskew - diag(Beta)) * D;
700  break;
701  }
702 
703 
704  } // SWITCH usedNlinearity
705 
706  } // FOR maxIterations
707 
708  W = transpose(B) * whiteningMatrix;
709 
710 
711  } // IF FICA_APPROACH_SYMM APPROACH
712 
713  // DEFLATION
714  else {
715 
716  // FC 01/12/05
717  A = zeros(whiteningMatrix.cols(), numOfIC);
718  // A = zeros( vectorSize, numOfIC );
719  mat B = zeros(vectorSize, numOfIC);
720  W = transpose(B) * whiteningMatrix;
721  int round = 1;
722  int numFailures = 0;
723 
724  while (round <= numOfIC) {
725 
726  myy = myyOrig;
727 
728  usedNlinearity = gOrig;
729  stroke = 0;
730 
731  notFine = 1;
732  loong = 0;
733  int endFinetuning = 0;
734 
735  vec w = zeros(vectorSize);
736 
737  if (initialStateMode == 0)
738 
739  w = randu(vectorSize) - 0.5;
740 
741  else w = whiteningMatrix * guess.get_col(round);
742 
743  w = w - B * transpose(B) * w;
744 
745  w /= norm(w);
746 
747  vec wOld = zeros(vectorSize);
748  vec wOld2 = zeros(vectorSize);
749 
750  int i = 1;
751  int gabba = 1;
752 
753  while (i <= maxNumIterations + gabba) {
754 
755  w = w - B * transpose(B) * w;
756 
757  w /= norm(w);
758 
759  if (notFine) {
760 
761  if (i == maxNumIterations + 1) {
762 
763  round--;
764 
765  numFailures++;
766 
767  if (numFailures > failureLimit) {
768 
769  if (round == 0) {
770 
771  A = dewhiteningMatrix * B;
772  W = transpose(B) * whiteningMatrix;
773 
774  } // IF round
775 
776  return false;
777 
778  } // IF numFailures > failureLimit
779 
780  break;
781 
782  } // IF i == maxNumIterations+1
783 
784  } // IF NOTFINE
785 
786  else if (i >= endFinetuning) wOld = w;
787 
788  if (norm(w - wOld) < epsilon || norm(w + wOld) < epsilon) {
789 
790  if (fineTuningEnabled && notFine) {
791 
792  notFine = 0;
793  gabba = maxFinetune;
794  wOld = zeros(vectorSize);
795  wOld2 = zeros(vectorSize);
796  usedNlinearity = gFine;
797  myy = myyK * myyOrig;
798  endFinetuning = maxFinetune + i;
799 
800  } // IF finetuning
801 
802  else {
803 
804  numFailures = 0;
805 
806  B.set_col(round - 1, w);
807 
808  A.set_col(round - 1, dewhiteningMatrix*w);
809 
810  W.set_row(round - 1, transpose(whiteningMatrix)*w);
811 
812  break;
813 
814  } // ELSE finetuning
815 
816  } // IF epsilon
817 
818  else if (stabilizationEnabled) {
819 
820  if (stroke == 0.0 && (norm(w - wOld2) < epsilon || norm(w + wOld2) < epsilon)) {
821 
822  stroke = myy;
823  myy /= 2.0 ;
824 
825  if (mod(usedNlinearity, 2) == 0) {
826 
827  usedNlinearity++;
828 
829  } // IF MOD
830 
831  }// IF !stroke
832 
833  else if (stroke != 0.0) {
834 
835  myy = stroke;
836  stroke = 0.0;
837 
838  if (myy == 1 && mod(usedNlinearity, 2) != 0) {
839  usedNlinearity--;
840  }
841 
842  } // IF Stroke
843 
844  else if (notFine && !loong && i > maxNumIterations / 2) {
845 
846  loong = 1;
847  myy /= 2.0;
848 
849  if (mod(usedNlinearity, 2) == 0) {
850 
851  usedNlinearity++;
852 
853  } // IF MOD
854 
855  } // IF notFine
856 
857  } // IF stabilization
858 
859 
860  wOld2 = wOld;
861  wOld = w;
862 
863  switch (usedNlinearity) {
864 
865  // pow3
866  case FICA_NONLIN_POW3 : {
867  w = (X * pow(transpose(X) * w, 3)) / numSamples - 3 * w;
868  break;
869  }
870  case(FICA_NONLIN_POW3+1) : {
871  vec Y = transpose(X) * w;
872  vec Gpow3 = X * pow(Y, 3) / numSamples;
873  double Beta = dot(w, Gpow3);
874  w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
875  break;
876  }
877  case(FICA_NONLIN_POW3+2) : {
878  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
879  w = (Xsub * pow(transpose(Xsub) * w, 3)) / Xsub.cols() - 3 * w;
880  break;
881  }
882  case(FICA_NONLIN_POW3+3) : {
883  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
884  vec Gpow3 = Xsub * pow(transpose(Xsub) * w, 3) / (Xsub.cols());
885  double Beta = dot(w, Gpow3);
886  w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
887  break;
888  }
889 
890  // TANH
891  case FICA_NONLIN_TANH : {
892  vec hypTan = tanh(a1 * transpose(X) * w);
893  w = (X * hypTan - a1 * sum(1 - pow(hypTan, 2)) * w) / numSamples;
894  break;
895  }
896  case(FICA_NONLIN_TANH+1) : {
897  vec Y = transpose(X) * w;
898  vec hypTan = tanh(a1 * Y);
899  double Beta = dot(w, X * hypTan);
900  w = w - myy * ((X * hypTan - Beta * w) / (a1 * sum(1 - pow(hypTan, 2)) - Beta));
901  break;
902  }
903  case(FICA_NONLIN_TANH+2) : {
904  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
905  vec hypTan = tanh(a1 * transpose(Xsub) * w);
906  w = (Xsub * hypTan - a1 * sum(1 - pow(hypTan, 2)) * w) / Xsub.cols();
907  break;
908  }
909  case(FICA_NONLIN_TANH+3) : {
910  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
911  vec hypTan = tanh(a1 * transpose(Xsub) * w);
912  double Beta = dot(w, Xsub * hypTan);
913  w = w - myy * ((Xsub * hypTan - Beta * w) / (a1 * sum(1 - pow(hypTan, 2)) - Beta));
914  break;
915  }
916 
917  // GAUSS
918  case FICA_NONLIN_GAUSS : {
919  vec u = transpose(X) * w;
920  vec Usquared = pow(u, 2);
921  vec ex = exp(-a2 * Usquared / 2);
922  vec gauss = elem_mult(u, ex);
923  vec dGauss = elem_mult(1 - a2 * Usquared, ex);
924  w = (X * gauss - sum(dGauss) * w) / numSamples;
925  break;
926  }
927  case(FICA_NONLIN_GAUSS+1) : {
928  vec u = transpose(X) * w;
929  vec Usquared = pow(u, 2);
930 
931  vec ex = exp(-a2 * Usquared / 2);
932  vec gauss = elem_mult(u, ex);
933  vec dGauss = elem_mult(1 - a2 * Usquared, ex);
934  double Beta = dot(w, X * gauss);
935  w = w - myy * ((X * gauss - Beta * w) / (sum(dGauss) - Beta));
936  break;
937  }
938  case(FICA_NONLIN_GAUSS+2) : {
939  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
940  vec u = transpose(Xsub) * w;
941  vec Usquared = pow(u, 2);
942  vec ex = exp(-a2 * Usquared / 2);
943  vec gauss = elem_mult(u, ex);
944  vec dGauss = elem_mult(1 - a2 * Usquared, ex);
945  w = (Xsub * gauss - sum(dGauss) * w) / Xsub.cols();
946  break;
947  }
948  case(FICA_NONLIN_GAUSS+3) : {
949  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
950  vec u = transpose(Xsub) * w;
951  vec Usquared = pow(u, 2);
952  vec ex = exp(-a2 * Usquared / 2);
953  vec gauss = elem_mult(u, ex);
954  vec dGauss = elem_mult(1 - a2 * Usquared, ex);
955  double Beta = dot(w, Xsub * gauss);
956  w = w - myy * ((Xsub * gauss - Beta * w) / (sum(dGauss) - Beta));
957  break;
958  }
959 
960  // SKEW
961  case FICA_NONLIN_SKEW : {
962  w = (X * (pow(transpose(X) * w, 2))) / numSamples;
963  break;
964  }
965  case(FICA_NONLIN_SKEW+1) : {
966  vec Y = transpose(X) * w;
967  vec Gskew = X * pow(Y, 2) / numSamples;
968  double Beta = dot(w, Gskew);
969  w = w - myy * (Gskew - Beta * w / (-Beta));
970  break;
971  }
972  case(FICA_NONLIN_SKEW+2) : {
973  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
974  w = (Xsub * (pow(transpose(Xsub) * w, 2))) / Xsub.cols();
975  break;
976  }
977  case(FICA_NONLIN_SKEW+3) : {
978  mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
979  vec Gskew = Xsub * pow(transpose(Xsub) * w, 2) / Xsub.cols();
980  double Beta = dot(w, Gskew);
981  w = w - myy * (Gskew - Beta * w) / (-Beta);
982  break;
983  }
984 
985  } // SWITCH nonLinearity
986 
987  w /= norm(w);
988  i++;
989 
990  } // WHILE i<= maxNumIterations+gabba
991 
992  round++;
993 
994  } // While round <= numOfIC
995 
996  } // ELSE Deflation
997  return true;
998 } // FPICA
itpp::Vec::ivec
Vec< int > ivec
Definition of integer vector type.
Definition: vec.h:541
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
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
eigen.h
Definitions of eigenvalue decomposition functions.
itpp::Fast_ICA::get_separating_matrix
mat get_separating_matrix()
Get separating matrix.
Definition: fastica.cpp:220
FICA_APPROACH_SYMM
#define FICA_APPROACH_SYMM
Use symmetric approach : compute all ICs at a time.
Definition: fastica.h:71
it_warning
#define it_warning(s)
Display a warning message.
Definition: itassert.h:173
itpp::Fast_ICA::set_stabilization
void set_stabilization(bool in_stabilization)
Set stabilization mode true or off.
Definition: fastica.cpp:200
itpp::Fast_ICA::set_non_linearity
void set_non_linearity(int in_g)
Set non-linearity.
Definition: fastica.cpp:186
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
itpp::to_ivec
ivec to_ivec(const Vec< T > &v)
Converts a Vec<T> to ivec.
Definition: converters.h:79
itpp::Fast_ICA::get_dewhitening_matrix
mat get_dewhitening_matrix()
Get the de-whitening matrix.
Definition: fastica.cpp:230
itpp::Fast_ICA::get_nrof_independent_components
int get_nrof_independent_components()
Get number of independent components.
Definition: fastica.cpp:224
itpp::eps
const double eps
Constant eps.
Definition: misc.h:109
itpp::Fast_ICA::separate
bool separate(void)
Explicit launch of main FastICA function.
Definition: fastica.cpp:126
mpower
static mat mpower(const mat A, const double y)
Local functions for FastICA.
Definition: fastica.cpp:397
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::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::Fast_ICA::get_independent_components
mat get_independent_components()
Get separated signals.
Definition: fastica.cpp:222
itpp::Fast_ICA::set_epsilon
void set_epsilon(double fl_epsilon)
Set convergence parameter .
Definition: fastica.cpp:196
svec.h
Sparse Vector Class definitions.
itpp::Vec::cvec
Vec< std::complex< double > > cvec
Definition of complex<double> vector type.
Definition: vec.h:535
itpp::Vec::sort
void sort(Vec< T > &data, SORTING_METHOD method=INTROSORT)
Sort the data vector in increasing order.
Definition: sort.h:187
itpp
itpp namespace
Definition: itmex.h:37
itpp::mod
int mod(int k, int n)
Calculates the modulus, i.e. the signed reminder after division.
Definition: elem_math.h:166
itpp::size
int size(const Vec< T > &v)
Length of vector.
Definition: matfunc.h:55
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::set_pca_only
void set_pca_only(bool in_PCAonly)
If true, only perform Principal Component Analysis (default = false)
Definition: fastica.cpp:210
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
sigfun.h
Definitions of signal processing functions.
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
FICA_NONLIN_TANH
#define FICA_NONLIN_TANH
Use tanh(x) non-linearity.
Definition: fastica.h:76
itpp::Fast_ICA::get_mixing_matrix
mat get_mixing_matrix()
Get mixing matrix.
Definition: fastica.cpp:218
orth
static mat orth(const mat A)
Local functions for FastICA.
Definition: fastica.cpp:368
random.h
Definition of classes for random number generators.
itpp::min
T min(const Vec< T > &in)
Minimum value of vector.
Definition: min_max.h:125
misc_stat.h
Miscellaneous statistics functions and classes - header file.
matfunc.h
Various functions on vectors and matrices - header file.
itpp::exp
vec exp(const vec &x)
Exp of the elements of a vector x.
Definition: log_exp.h:155
getSamples
static ivec getSamples(const int max, const double percentage)
Local functions for FastICA.
Definition: fastica.cpp:417
itpp::Mat::mat
Mat< double > mat
Default Matrix Type.
Definition: mat.h:482
itpp::Vec::vec
Vec< double > vec
Definition of double vector type.
Definition: vec.h:529
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
FICA_NONLIN_GAUSS
#define FICA_NONLIN_GAUSS
Use Gaussian non-linearity.
Definition: fastica.h:78
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::Sparse_Vec::sparse_vec
Sparse_Vec< double > sparse_vec
Type definition of a double sparse vector.
Definition: svec.h:347
FICA_NONLIN_POW3
#define FICA_NONLIN_POW3
Use x^3 non-linearity.
Definition: fastica.h:74
itpp::Fast_ICA::set_a2
void set_a2(double fl_a2)
Set parameter.
Definition: fastica.cpp:192
fastica.h
Definition of FastICA (Independent Component Analysis) for IT++.
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
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.
itpp::bin::abs
bin abs(const bin &inbin)
absolute value of bin
Definition: binary.h:174
itpp::max
T max(const Vec< T > &v)
Maximum value of vector.
Definition: min_max.h:45
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
svd.h
Definitions of Singular Value Decompositions.
itpp::Vec::to_ivec
ivec to_ivec(const Vec< T > &v)
Converts a Vec<T> to ivec.
Definition: converters.h:79
remmean
static void remmean(mat inVectors, mat &outVectors, vec &meanValue)
Local functions for FastICA.
Definition: fastica.cpp:335
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
resampling.h
Resampling functions - header file.
itpp::Fast_ICA::set_mu
void set_mu(double fl_mu)
Set parameter.
Definition: fastica.cpp:194
FICA_INIT_RAND
#define FICA_INIT_RAND
Set random start for Fast_ICA.
Definition: fastica.h:83
itpp::Fast_ICA::get_whitening_matrix
mat get_whitening_matrix()
Get the whitening matrix.
Definition: fastica.cpp:228
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
selcol
static void selcol(const mat oldMatrix, const vec maskVector, mat &newMatrix)
Local functions for FastICA.
Definition: fastica.cpp:237
itpp::round
ITPP_EXPORT double round(double x)
Round to nearest integer, return result in double.
itpp::Fast_ICA::Fast_ICA
Fast_ICA(mat ma_mixed_sig)
Constructor.
Definition: fastica.cpp:99
itpp::Fast_ICA::set_sample_size
void set_sample_size(double fl_sampleSize)
Set sample size.
Definition: fastica.cpp:198
sumcol
static vec sumcol(const mat A)
Local functions for FastICA.
Definition: fastica.cpp:433
itpp::mean
double mean(const vec &v)
The mean value.
Definition: misc_stat.cpp:36
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::get_principal_eigenvectors
mat get_principal_eigenvectors()
Get nrIC first columns of the de-whitening matrix.
Definition: fastica.cpp:226
itpp::Fast_ICA::set_a1
void set_a1(double fl_a1)
Set parameter.
Definition: fastica.cpp:190
itpp::zeros
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
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::norm
double norm(const cvec &v)
Calculate the 2-norm: norm(v)=sqrt(sum(abs(v).^2))
Definition: misc_stat.cpp:77
sort.h
Sorting functions.
itpp::cov
mat cov(const mat &X, bool is_zero_mean)
Covariance matrix calculation.
Definition: sigfun.cpp:226
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::sqrt
vec sqrt(const vec &x)
Square root of the elements.
Definition: elem_math.h:123
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
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
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
itpp::Fast_ICA::get_white_sig
mat get_white_sig()
Get whitened signals.
Definition: fastica.cpp:232
itpp::randu
double randu(void)
Generates a random uniform (0,1) number.
Definition: random.h:804
itpp::full
Mat< T > full(const Sparse_Mat< T > &s)
Convert a sparse matrix s into its dense representation.
Definition: smat.h:998
trig_hyp.h
Trigonometric and hyperbolic functions - header file.
itpp::Mat::cmat
Mat< std::complex< double > > cmat
Default Complex Matrix Type.
Definition: mat.h:488
min_max.h
Minimum and maximum functions on vectors and matrices.
itpp::Fast_ICA::set_fine_tune
void set_fine_tune(bool in_finetune)
Set fine tuning.
Definition: fastica.cpp:188
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::tanh
vec tanh(const vec &x)
Tan hyperbolic function.
Definition: trig_hyp.h:97
FICA_TOL
#define FICA_TOL
Eigenvalues of the covariance matrix lower than FICA_TOL are discarded for analysis.
Definition: fastica.h:88
FICA_NONLIN_SKEW
#define FICA_NONLIN_SKEW
Use skew non-linearity.
Definition: fastica.h:80
itpp::sum
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition: matfunc.h:59
SourceForge Logo

Generated on Thu Apr 11 2019 00:00:00 for IT++ by Doxygen 1.8.18