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
void set_sample_size(double fl_sampleSize)
Set sample size.
Definition: fastica.cpp:198
Definition of FastICA (Independent Component Analysis) for IT++.
Various functions on vectors and matrices - header file.
#define FICA_APPROACH_SYMM
Use symmetric approach : compute all ICs at a time.
Definition: fastica.h:71
void set_a1(double fl_a1)
Set parameter.
Definition: fastica.cpp:190
Mat< Num_T > concat_horizontal(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Horizontal concatenation of two matrices.
Definition: mat.h:1194
double randu(void)
Generates a random uniform (0,1) number.
Definition: random.h:804
void set_mu(double fl_mu)
Set parameter.
Definition: fastica.cpp:194
static void remmean(mat inVectors, mat &outVectors, vec &meanValue)
Local functions for FastICA.
Definition: fastica.cpp:335
Definitions of eigenvalue decomposition functions.
bool separate(void)
Explicit launch of main FastICA function.
Definition: fastica.cpp:126
#define FICA_NONLIN_SKEW
Use skew non-linearity.
Definition: fastica.h:80
mat get_whitening_matrix()
Get the whitening matrix.
Definition: fastica.cpp:228
void set_init_guess(mat ma_initGuess)
Set initial guess matrix instead of random (default)
Definition: fastica.cpp:212
double norm(const cvec &v)
Calculate the 2-norm: norm(v)=sqrt(sum(abs(v).^2))
Definition: misc_stat.cpp:77
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
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
int mod(int k, int n)
Calculates the modulus, i.e. the signed reminder after division.
Definition: elem_math.h:166
void set_max_num_iterations(int in_maxNumIterations)
Set maximum number of iterations.
Definition: fastica.cpp:202
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition: matfunc.h:59
void set_epsilon(double fl_epsilon)
Set convergence parameter .
Definition: fastica.cpp:196
#define FICA_NONLIN_TANH
Use tanh(x) non-linearity.
Definition: fastica.h:76
void set_nrof_independent_components(int in_nrIC)
Set number of independent components to separate.
Definition: fastica.cpp:184
static ivec getSamples(const int max, const double percentage)
Local functions for FastICA.
Definition: fastica.cpp:417
#define FICA_NONLIN_GAUSS
Use Gaussian non-linearity.
Definition: fastica.h:78
Sorting functions.
Minimum and maximum functions on vectors and matrices.
Definitions of signal processing functions.
double mean(const vec &v)
The mean value.
Definition: misc_stat.cpp:36
ITPP_EXPORT double round(double x)
Round to nearest integer, return result in double.
#define FICA_APPROACH_DEFL
Use deflation approach : compute IC one-by-one in a Gram-Schmidt-like fashion.
Definition: fastica.h:69
int get_nrof_independent_components()
Get number of independent components.
Definition: fastica.cpp:224
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
const double eps
Constant eps.
Definition: misc.h:109
T min(const Vec< T > &in)
Minimum value of vector.
Definition: min_max.h:125
void set_last_eig(int in_lastEig)
Set last eigenvalue index to take into account.
Definition: fastica.cpp:208
int max_index(const Vec< T > &in)
Return the postion of the maximum element in the vector.
Definition: min_max.h:208
void set_pca_only(bool in_PCAonly)
If true, only perform Principal Component Analysis (default = false)
Definition: fastica.cpp:210
bool svd(const mat &A, vec &S)
Get singular values s of a real matrix A using SVD.
Definition: svd.cpp:185
vec exp(const vec &x)
Exp of the elements of a vector x.
Definition: log_exp.h:155
T max(const Vec< T > &v)
Maximum value of vector.
Definition: min_max.h:45
Trigonometric and hyperbolic functions - header file.
vec pow(const double x, const vec &y)
Calculates x to the power of y (x^y)
Definition: log_exp.h:176
mat cov(const mat &X, bool is_zero_mean)
Covariance matrix calculation.
Definition: sigfun.cpp:226
#define FICA_NONLIN_POW3
Use x^3 non-linearity.
Definition: fastica.h:74
int size(const Vec< T > &v)
Length of vector.
Definition: matfunc.h:55
mat get_separating_matrix()
Get separating matrix.
Definition: fastica.cpp:220
static vec sumcol(const mat A)
Local functions for FastICA.
Definition: fastica.cpp:433
Mat< T > full(const Sparse_Mat< T > &s)
Convert a sparse matrix s into its dense representation.
Definition: smat.h:998
void set_stabilization(bool in_stabilization)
Set stabilization mode true or off.
Definition: fastica.cpp:200
itpp namespace
Definition: itmex.h:36
static mat mpower(const mat A, const double y)
Local functions for FastICA.
Definition: fastica.cpp:397
void set_max_fine_tune(int in_maxFineTune)
Set maximum number of iterations for fine tuning.
Definition: fastica.cpp:204
void set_non_linearity(int in_g)
Set non-linearity.
Definition: fastica.cpp:186
Resampling functions - header file.
Miscellaneous statistics functions and classes - header file.
ivec to_ivec(const Vec< T > &v)
Converts a Vec<T> to ivec.
Definition: converters.h:79
#define FICA_INIT_GUESS
Set predefined start for Fast_ICA.
Definition: fastica.h:85
Definitions of special vectors and matrices.
#define it_warning(s)
Display a warning message.
Definition: itassert.h:173
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
mat get_principal_eigenvectors()
Get nrIC first columns of the de-whitening matrix.
Definition: fastica.cpp:226
mat get_independent_components()
Get separated signals.
Definition: fastica.cpp:222
bool eig_sym(const mat &A, vec &d, mat &V)
Calculates the eigenvalues and eigenvectors of a symmetric real matrix.
Definition: eigen.cpp:252
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
void transpose(const Mat< T > &m, Mat< T > &out)
Transposition of the matrix m returning the transposed matrix in out.
Definition: matfunc.h:308
#define FICA_INIT_RAND
Set random start for Fast_ICA.
Definition: fastica.h:83
Mat< T > reshape(const Mat< T > &m, int rows, int cols)
Reshape the matrix into an rows*cols matrix.
Definition: matfunc.h:822
mat get_dewhitening_matrix()
Get the de-whitening matrix.
Definition: fastica.cpp:230
void set_a2(double fl_a2)
Set parameter.
Definition: fastica.cpp:192
vec tanh(const vec &x)
Tan hyperbolic function.
Definition: trig_hyp.h:97
vec sqrt(const vec &x)
Square root of the elements.
Definition: elem_math.h:123
Definitions of Singular Value Decompositions.
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
static mat orth(const mat A)
Local functions for FastICA.
Definition: fastica.cpp:368
mat get_white_sig()
Get whitened signals.
Definition: fastica.cpp:232
Definition of classes for random number generators.
bin abs(const bin &inbin)
absolute value of bin
Definition: binary.h:174
Fast_ICA(mat ma_mixed_sig)
Constructor.
Definition: fastica.cpp:99
void set_first_eig(int in_firstEig)
Set first eigenvalue index to take into account.
Definition: fastica.cpp:206
Sparse Vector Class definitions.
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
#define FICA_TOL
Eigenvalues of the covariance matrix lower than FICA_TOL are discarded for analysis.
Definition: fastica.h:88
mat get_mixing_matrix()
Get mixing matrix.
Definition: fastica.cpp:218
void set_fine_tune(bool in_finetune)
Set fine tuning.
Definition: fastica.cpp:188
void set_approach(int in_approach)
Set approach : FICA_APPROACH_DEFL or FICA_APPROACH_SYMM (default)
Definition: fastica.cpp:182
static void selcol(const mat oldMatrix, const vec maskVector, mat &newMatrix)
Local functions for FastICA.
Definition: fastica.cpp:237

Generated on Thu Jun 21 2018 16:06:18 for IT++ by Doxygen 1.8.13