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
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:25 for IT++ by Doxygen 1.8.2