00001 00062 #include <itpp/signal/fastica.h> 00063 #include <itpp/signal/sigfun.h> 00064 #include <itpp/signal/resampling.h> 00065 #include <itpp/base/algebra/eigen.h> 00066 #include <itpp/base/algebra/svd.h> 00067 #include <itpp/base/math/trig_hyp.h> 00068 #include <itpp/base/matfunc.h> 00069 #include <itpp/base/random.h> 00070 #include <itpp/base/sort.h> 00071 #include <itpp/base/specmat.h> 00072 #include <itpp/base/svec.h> 00073 #include <itpp/base/math/min_max.h> 00074 #include <itpp/stat/misc_stat.h> 00075 00076 00077 using namespace itpp; 00078 00079 00084 static void selcol(const mat oldMatrix, const vec maskVector, mat & newMatrix); 00085 static void pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat & Es, vec & Ds); 00086 static void remmean(mat inVectors, mat & outVectors, vec & meanValue); 00087 static void whitenv(const mat vectors, const mat E, const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix); 00088 static mat orth(const mat A); 00089 static mat mpower(const mat A, const double y); 00090 static ivec getSamples(const int max, const double percentage); 00091 static vec sumcol(const mat A); 00092 static void 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); 00095 namespace itpp 00096 { 00097 00098 // Constructor, init default values 00099 Fast_ICA::Fast_ICA(mat ma_mixedSig) 00100 { 00101 00102 // Init default values 00103 approach = FICA_APPROACH_SYMM; 00104 g = FICA_NONLIN_POW3; 00105 finetune = true; 00106 a1 = 1.0; 00107 a2 = 1.0; 00108 mu = 1.0; 00109 epsilon = 0.0001; 00110 sampleSize = 1.0; 00111 stabilization = false; 00112 maxNumIterations = 100000; 00113 maxFineTune = 100; 00114 firstEig = 1; 00115 00116 mixedSig = ma_mixedSig; 00117 00118 lastEig = mixedSig.rows(); 00119 numOfIC = mixedSig.rows(); 00120 PCAonly = false; 00121 initState = FICA_INIT_RAND; 00122 00123 } 00124 00125 // Call main function 00126 void Fast_ICA::separate(void) 00127 { 00128 00129 int Dim = numOfIC; 00130 00131 mat mixedSigC; 00132 vec mixedMean; 00133 00134 mat guess; 00135 if (initState == FICA_INIT_RAND) 00136 guess = zeros(Dim, Dim); 00137 else 00138 guess = mat(initGuess); 00139 00140 VecPr = zeros(mixedSig.rows(), numOfIC); 00141 00142 icasig = zeros(numOfIC, mixedSig.cols()); 00143 00144 remmean(mixedSig, mixedSigC, mixedMean); 00145 00146 pcamat(mixedSigC, numOfIC, firstEig, lastEig, E, D); 00147 00148 whitenv(mixedSigC, E, diag(D), whitesig, whiteningMatrix, dewhiteningMatrix); 00149 00150 00151 ivec NcFirst = to_ivec(zeros(numOfIC)); 00152 vec NcVp = D; 00153 for (int i = 0; i < NcFirst.size(); i++) { 00154 00155 NcFirst(i) = max_index(NcVp); 00156 NcVp(NcFirst(i)) = 0.0; 00157 VecPr.set_col(i, dewhiteningMatrix.get_col(i)); 00158 00159 } 00160 00161 if (PCAonly == false) { 00162 00163 Dim = whitesig.rows(); 00164 00165 if (numOfIC > Dim) numOfIC = Dim; 00166 00167 fpica(whitesig, whiteningMatrix, dewhiteningMatrix, approach, numOfIC, g, finetune, a1, a2, mu, stabilization, epsilon, maxNumIterations, maxFineTune, initState, guess, sampleSize, A, W); 00168 00169 icasig = W * mixedSig; 00170 00171 } 00172 00173 else { // PCA only : returns E as IcaSig 00174 icasig = VecPr; 00175 } 00176 } 00177 00178 void Fast_ICA::set_approach(int in_approach) { approach = in_approach; if (approach == FICA_APPROACH_DEFL) finetune = true; } 00179 00180 void Fast_ICA::set_nrof_independent_components(int in_nrIC) { numOfIC = in_nrIC; } 00181 00182 void Fast_ICA::set_non_linearity(int in_g) { g = in_g; } 00183 00184 void Fast_ICA::set_fine_tune(bool in_finetune) { finetune = in_finetune; } 00185 00186 void Fast_ICA::set_a1(double fl_a1) { a1 = fl_a1; } 00187 00188 void Fast_ICA::set_a2(double fl_a2) { a2 = fl_a2; } 00189 00190 void Fast_ICA::set_mu(double fl_mu) { mu = fl_mu; } 00191 00192 void Fast_ICA::set_epsilon(double fl_epsilon) { epsilon = fl_epsilon; } 00193 00194 void Fast_ICA::set_sample_size(double fl_sampleSize) { sampleSize = fl_sampleSize; } 00195 00196 void Fast_ICA::set_stabilization(bool in_stabilization) { stabilization = in_stabilization; } 00197 00198 void Fast_ICA::set_max_num_iterations(int in_maxNumIterations) { maxNumIterations = in_maxNumIterations; } 00199 00200 void Fast_ICA::set_max_fine_tune(int in_maxFineTune) { maxFineTune = in_maxFineTune; } 00201 00202 void Fast_ICA::set_first_eig(int in_firstEig) { firstEig = in_firstEig; } 00203 00204 void Fast_ICA::set_last_eig(int in_lastEig) { lastEig = in_lastEig; } 00205 00206 void Fast_ICA::set_pca_only(bool in_PCAonly) { PCAonly = in_PCAonly; } 00207 00208 void Fast_ICA::set_init_guess(mat ma_initGuess) 00209 { 00210 initGuess = ma_initGuess; 00211 initState = FICA_INIT_GUESS; 00212 } 00213 00214 mat Fast_ICA::get_mixing_matrix() { if (PCAonly) { it_warning("No ICA performed."); return (zeros(1, 1));} else return A; } 00215 00216 mat Fast_ICA::get_separating_matrix() { if (PCAonly) { it_warning("No ICA performed."); return(zeros(1, 1)); } else return W; } 00217 00218 mat Fast_ICA::get_independent_components() { if (PCAonly) { it_warning("No ICA performed."); return(zeros(1, 1)); } else return icasig; } 00219 00220 int Fast_ICA::get_nrof_independent_components() { return numOfIC; } 00221 00222 mat Fast_ICA::get_principal_eigenvectors() { return VecPr; } 00223 00224 mat Fast_ICA::get_whitening_matrix() { return whiteningMatrix; } 00225 00226 mat Fast_ICA::get_dewhitening_matrix() { return dewhiteningMatrix; } 00227 00228 mat Fast_ICA::get_white_sig() { return whitesig; } 00229 00230 } // namespace itpp 00231 00232 00233 static void selcol(const mat oldMatrix, const vec maskVector, mat & newMatrix) 00234 { 00235 00236 int numTaken = 0; 00237 00238 for (int i = 0; i < size(maskVector); i++) if (maskVector(i) == 1) numTaken++; 00239 00240 newMatrix = zeros(oldMatrix.rows(), numTaken); 00241 00242 numTaken = 0; 00243 00244 for (int i = 0; i < size(maskVector); i++) { 00245 00246 if (maskVector(i) == 1) { 00247 00248 newMatrix.set_col(numTaken, oldMatrix.get_col(i)); 00249 numTaken++; 00250 00251 } 00252 } 00253 00254 return; 00255 00256 } 00257 00258 static void pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat & Es, vec & Ds) 00259 { 00260 00261 mat Et; 00262 vec Dt; 00263 cmat Ec; 00264 cvec Dc; 00265 double lowerLimitValue = 0.0, 00266 higherLimitValue = 0.0; 00267 00268 int oldDimension = vectors.rows(); 00269 00270 mat covarianceMatrix = cov(transpose(vectors), 0); 00271 00272 eig_sym(covarianceMatrix, Dt, Et); 00273 00274 int maxLastEig = 0; 00275 00276 // Compute rank 00277 for (int i = 0; i < Dt.length(); i++) if (Dt(i) > FICA_TOL) maxLastEig++; 00278 00279 // Force numOfIC components 00280 if (maxLastEig > numOfIC) maxLastEig = numOfIC; 00281 00282 vec eigenvalues = zeros(size(Dt)); 00283 vec eigenvalues2 = zeros(size(Dt)); 00284 00285 eigenvalues2 = Dt; 00286 00287 sort(eigenvalues2); 00288 00289 vec lowerColumns = zeros(size(Dt)); 00290 00291 for (int i = 0; i < size(Dt); i++) eigenvalues(i) = eigenvalues2(size(Dt) - i - 1); 00292 00293 if (lastEig > maxLastEig) lastEig = maxLastEig; 00294 00295 if (lastEig < oldDimension) lowerLimitValue = (eigenvalues(lastEig - 1) + eigenvalues(lastEig)) / 2; 00296 else lowerLimitValue = eigenvalues(oldDimension - 1) - 1; 00297 00298 for (int i = 0; i < size(Dt); i++) if (Dt(i) > lowerLimitValue) lowerColumns(i) = 1; 00299 00300 if (firstEig > 1) higherLimitValue = (eigenvalues(firstEig - 2) + eigenvalues(firstEig - 1)) / 2; 00301 else higherLimitValue = eigenvalues(0) + 1; 00302 00303 vec higherColumns = zeros(size(Dt)); 00304 for (int i = 0; i < size(Dt); i++) if (Dt(i) < higherLimitValue) higherColumns(i) = 1; 00305 00306 vec selectedColumns = zeros(size(Dt)); 00307 for (int i = 0; i < size(Dt); i++) selectedColumns(i) = (lowerColumns(i) == 1 && higherColumns(i) == 1) ? 1 : 0; 00308 00309 selcol(Et, selectedColumns, Es); 00310 00311 int numTaken = 0; 00312 00313 for (int i = 0; i < size(selectedColumns); i++) if (selectedColumns(i) == 1) numTaken++; 00314 00315 Ds = zeros(numTaken); 00316 00317 numTaken = 0; 00318 00319 for (int i = 0; i < size(Dt); i++) 00320 if (selectedColumns(i) == 1) { 00321 Ds(numTaken) = Dt(i); 00322 numTaken++; 00323 } 00324 00325 return; 00326 00327 } 00328 00329 00330 static void remmean(mat inVectors, mat & outVectors, vec & meanValue) 00331 { 00332 00333 outVectors = zeros(inVectors.rows(), inVectors.cols()); 00334 meanValue = zeros(inVectors.rows()); 00335 00336 for (int i = 0; i < inVectors.rows(); i++) { 00337 00338 meanValue(i) = mean(inVectors.get_row(i)); 00339 00340 for (int j = 0; j < inVectors.cols(); j++) outVectors(i, j) = inVectors(i, j) - meanValue(i); 00341 00342 } 00343 00344 } 00345 00346 static void whitenv(const mat vectors, const mat E, const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix) 00347 { 00348 00349 whiteningMatrix = zeros(E.cols(), E.rows()); 00350 dewhiteningMatrix = zeros(E.rows(), E.cols()); 00351 00352 for (int i = 0; i < D.cols(); i++) { 00353 whiteningMatrix.set_row(i, std::pow(std::sqrt(D(i, i)), -1)*E.get_col(i)); 00354 dewhiteningMatrix.set_col(i, std::sqrt(D(i, i))*E.get_col(i)); 00355 } 00356 00357 newVectors = whiteningMatrix * vectors; 00358 00359 return; 00360 00361 } 00362 00363 static mat orth(const mat A) 00364 { 00365 00366 mat Q; 00367 mat U, V; 00368 vec S; 00369 double eps = 2.2e-16; 00370 double tol = 0.0; 00371 int mmax = 0; 00372 int r = 0; 00373 00374 svd(A, U, S, V); 00375 if (A.rows() > A.cols()) { 00376 00377 U = U(0, U.rows() - 1, 0, A.cols() - 1); 00378 S = S(0, A.cols() - 1); 00379 } 00380 00381 mmax = (A.rows() > A.cols()) ? A.rows() : A.cols(); 00382 00383 tol = mmax * eps * max(S); 00384 00385 for (int i = 0; i < size(S); i++) if (S(i) > tol) r++; 00386 00387 Q = U(0, U.rows() - 1, 0, r - 1); 00388 00389 return (Q); 00390 } 00391 00392 static mat mpower(const mat A, const double y) 00393 { 00394 00395 mat T = zeros(A.rows(), A.cols()); 00396 mat dd = zeros(A.rows(), A.cols()); 00397 vec d = zeros(A.rows()); 00398 vec dOut = zeros(A.rows()); 00399 00400 eig_sym(A, d, T); 00401 00402 dOut = pow(d, y); 00403 00404 diag(dOut, dd); 00405 00406 for (int i = 0; i < T.cols(); i++) T.set_col(i, T.get_col(i) / norm(T.get_col(i))); 00407 00408 return (T*dd*transpose(T)); 00409 00410 } 00411 00412 static ivec getSamples(const int max, const double percentage) 00413 { 00414 00415 vec rd = randu(max); 00416 sparse_vec sV; 00417 ivec out; 00418 int sZ = 0; 00419 00420 for (int i = 0; i < max; i++) if (rd(i) < percentage) { sV.add_elem(sZ, i); sZ++; } 00421 00422 out = to_ivec(full(sV)); 00423 00424 return (out); 00425 00426 } 00427 00428 static vec sumcol(const mat A) 00429 { 00430 00431 vec out = zeros(A.cols()); 00432 00433 for (int i = 0; i < A.cols(); i++) { out(i) = sum(A.get_col(i)); } 00434 00435 return (out); 00436 00437 } 00438 00439 static void 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) 00440 { 00441 00442 int vectorSize = X.rows(); 00443 int numSamples = X.cols(); 00444 int gOrig = g; 00445 int gFine = finetune + 1; 00446 double myyOrig = myy; 00447 double myyK = 0.01; 00448 int failureLimit = 5; 00449 int usedNlinearity = 0; 00450 double stroke = 0.0; 00451 int notFine = 1; 00452 int loong = 0; 00453 int initialStateMode = initState; 00454 double minAbsCos = 0.0, minAbsCos2 = 0.0; 00455 00456 if (sampleSize * numSamples < 1000) sampleSize = (1000 / (double)numSamples < 1.0) ? 1000 / (double)numSamples : 1.0; 00457 00458 if (sampleSize != 1.0) gOrig += 2; 00459 if (myy != 1.0) gOrig += 1; 00460 00461 int fineTuningEnabled = 1; 00462 00463 if (!finetune) { 00464 if (myy != 1.0) gFine = gOrig; 00465 else gFine = gOrig + 1; 00466 fineTuningEnabled = 0; 00467 } 00468 00469 int stabilizationEnabled = stabilization; 00470 00471 if (!stabilization && myy != 1.0) stabilizationEnabled = true; 00472 00473 usedNlinearity = gOrig; 00474 00475 if (initState == FICA_INIT_GUESS && guess.rows() != whiteningMatrix.cols()) { 00476 initialStateMode = 0; 00477 00478 } 00479 else if (guess.cols() < numOfIC) { 00480 00481 mat guess2 = randu(guess.rows(), numOfIC - guess.cols()) - 0.5; 00482 guess = concat_horizontal(guess, guess2); 00483 } 00484 else if (guess.cols() > numOfIC) guess = guess(0, guess.rows() - 1, 0, numOfIC - 1); 00485 00486 if (approach == FICA_APPROACH_SYMM) { 00487 00488 usedNlinearity = gOrig; 00489 stroke = 0; 00490 notFine = 1; 00491 loong = 0; 00492 00493 A = zeros(vectorSize, numOfIC); 00494 mat B = zeros(vectorSize, numOfIC); 00495 00496 if (initialStateMode == 0) B = orth(randu(vectorSize, numOfIC) - 0.5); 00497 else B = whiteningMatrix * guess; 00498 00499 mat BOld = zeros(B.rows(), B.cols()); 00500 mat BOld2 = zeros(B.rows(), B.cols()); 00501 00502 for (int round = 0; round < maxNumIterations; round++) { 00503 00504 if (round == maxNumIterations - 1) { 00505 00506 // If there is a convergence problem, 00507 // we still want ot return something. 00508 // This is difference with original 00509 // Matlab implementation. 00510 A = dewhiteningMatrix * B; 00511 W = transpose(B) * whiteningMatrix; 00512 00513 return; 00514 } 00515 00516 B = B * mpower(transpose(B) * B , -0.5); 00517 00518 minAbsCos = min(abs(diag(transpose(B) * BOld))); 00519 minAbsCos2 = min(abs(diag(transpose(B) * BOld2))); 00520 00521 if (1 - minAbsCos < epsilon) { 00522 00523 if (fineTuningEnabled && notFine) { 00524 00525 notFine = 0; 00526 usedNlinearity = gFine; 00527 myy = myyK * myyOrig; 00528 BOld = zeros(B.rows(), B.cols()); 00529 BOld2 = zeros(B.rows(), B.cols()); 00530 00531 } 00532 00533 else { 00534 00535 A = dewhiteningMatrix * B; 00536 break; 00537 00538 } 00539 00540 } // IF epsilon 00541 00542 else if (stabilizationEnabled) { 00543 if (!stroke && (1 - minAbsCos2 < epsilon)) { 00544 00545 stroke = myy; 00546 myy /= 2; 00547 if (mod(usedNlinearity, 2) == 0) usedNlinearity += 1 ; 00548 00549 } 00550 else if (stroke) { 00551 00552 myy = stroke; 00553 stroke = 0; 00554 if (myy == 1 && mod(usedNlinearity, 2) != 0) usedNlinearity -= 1; 00555 00556 } 00557 else if (!loong && (round > maxNumIterations / 2)) { 00558 00559 loong = 1; 00560 myy /= 2; 00561 if (mod(usedNlinearity, 2) == 0) usedNlinearity += 1; 00562 00563 } 00564 00565 } // stabilizationEnabled 00566 00567 BOld2 = BOld; 00568 BOld = B; 00569 00570 switch (usedNlinearity) { 00571 00572 // pow3 00573 case FICA_NONLIN_POW3 : { 00574 B = (X * pow(transpose(X) * B, 3)) / numSamples - 3 * B; 00575 break; 00576 } 00577 case(FICA_NONLIN_POW3+1) : { 00578 mat Y = transpose(X) * B; 00579 mat Gpow3 = pow(Y, 3); 00580 vec Beta = sumcol(pow(Y, 4)); 00581 mat D = diag(pow(Beta - 3 * numSamples , -1)); 00582 B = B + myy * B * (transpose(Y) * Gpow3 - diag(Beta)) * D; 00583 break; 00584 } 00585 case(FICA_NONLIN_POW3+2) : { 00586 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00587 B = (Xsub * pow(transpose(Xsub) * B, 3)) / Xsub.cols() - 3 * B; 00588 break; 00589 } 00590 case(FICA_NONLIN_POW3+3) : { 00591 mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B; 00592 mat Gpow3 = pow(Ysub, 3); 00593 vec Beta = sumcol(pow(Ysub, 4)); 00594 mat D = diag(pow(Beta - 3 * Ysub.rows() , -1)); 00595 B = B + myy * B * (transpose(Ysub) * Gpow3 - diag(Beta)) * D; 00596 break; 00597 } 00598 00599 // TANH 00600 case FICA_NONLIN_TANH : { 00601 mat hypTan = tanh(a1 * transpose(X) * B); 00602 B = (X * hypTan) / numSamples - elem_mult(reshape(repeat(sumcol(1 - pow(hypTan, 2)), B.rows()), B.rows(), B.cols()), B) / numSamples * a1; 00603 break; 00604 } 00605 case(FICA_NONLIN_TANH+1) : { 00606 mat Y = transpose(X) * B; 00607 mat hypTan = tanh(a1 * Y); 00608 vec Beta = sumcol(elem_mult(Y, hypTan)); 00609 vec Beta2 = sumcol(1 - pow(hypTan, 2)); 00610 mat D = diag(pow(Beta - a1 * Beta2 , -1)); 00611 B = B + myy * B * (transpose(Y) * hypTan - diag(Beta)) * D; 00612 break; 00613 } 00614 case(FICA_NONLIN_TANH+2) : { 00615 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00616 mat hypTan = tanh(a1 * transpose(Xsub) * B); 00617 B = (Xsub * hypTan) / Xsub.cols() - elem_mult(reshape(repeat(sumcol(1 - pow(hypTan, 2)), B.rows()), B.rows(), B.cols()), B) / Xsub.cols() * a1; 00618 break; 00619 } 00620 case(FICA_NONLIN_TANH+3) : { 00621 mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B; 00622 mat hypTan = tanh(a1 * Ysub); 00623 vec Beta = sumcol(elem_mult(Ysub, hypTan)); 00624 vec Beta2 = sumcol(1 - pow(hypTan, 2)); 00625 mat D = diag(pow(Beta - a1 * Beta2 , -1)); 00626 B = B + myy * B * (transpose(Ysub) * hypTan - diag(Beta)) * D; 00627 break; 00628 } 00629 00630 // GAUSS 00631 case FICA_NONLIN_GAUSS : { 00632 mat U = transpose(X) * B; 00633 mat Usquared = pow(U, 2); 00634 mat ex = exp(-a2 * Usquared / 2); 00635 mat gauss = elem_mult(U, ex); 00636 mat dGauss = elem_mult(1 - a2 * Usquared, ex); 00637 B = (X * gauss) / numSamples - elem_mult(reshape(repeat(sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / numSamples; 00638 break; 00639 } 00640 case(FICA_NONLIN_GAUSS+1) : { 00641 mat Y = transpose(X) * B; 00642 mat ex = exp(-a2 * pow(Y, 2) / 2); 00643 mat gauss = elem_mult(Y, ex); 00644 vec Beta = sumcol(elem_mult(Y, gauss)); 00645 vec Beta2 = sumcol(elem_mult(1 - a2 * pow(Y, 2), ex)); 00646 mat D = diag(pow(Beta - Beta2 , -1)); 00647 B = B + myy * B * (transpose(Y) * gauss - diag(Beta)) * D; 00648 break; 00649 } 00650 case(FICA_NONLIN_GAUSS+2) : { 00651 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00652 mat U = transpose(Xsub) * B; 00653 mat Usquared = pow(U, 2); 00654 mat ex = exp(-a2 * Usquared / 2); 00655 mat gauss = elem_mult(U, ex); 00656 mat dGauss = elem_mult(1 - a2 * Usquared, ex); 00657 B = (Xsub * gauss) / Xsub.cols() - elem_mult(reshape(repeat(sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / Xsub.cols(); 00658 break; 00659 } 00660 case(FICA_NONLIN_GAUSS+3) : { 00661 mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B; 00662 mat ex = exp(-a2 * pow(Ysub, 2) / 2); 00663 mat gauss = elem_mult(Ysub, ex); 00664 vec Beta = sumcol(elem_mult(Ysub, gauss)); 00665 vec Beta2 = sumcol(elem_mult(1 - a2 * pow(Ysub, 2), ex)); 00666 mat D = diag(pow(Beta - Beta2 , -1)); 00667 B = B + myy * B * (transpose(Ysub) * gauss - diag(Beta)) * D; 00668 break; 00669 } 00670 00671 // SKEW 00672 case FICA_NONLIN_SKEW : { 00673 B = (X * (pow(transpose(X) * B, 2))) / numSamples; 00674 break; 00675 } 00676 case(FICA_NONLIN_SKEW+1) : { 00677 mat Y = transpose(X) * B; 00678 mat Gskew = pow(Y, 2); 00679 vec Beta = sumcol(elem_mult(Y, Gskew)); 00680 mat D = diag(pow(Beta , -1)); 00681 B = B + myy * B * (transpose(Y) * Gskew - diag(Beta)) * D; 00682 break; 00683 } 00684 case(FICA_NONLIN_SKEW+2) : { 00685 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00686 B = (Xsub * (pow(transpose(Xsub) * B, 2))) / Xsub.cols(); 00687 break; 00688 } 00689 case(FICA_NONLIN_SKEW+3) : { 00690 mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B; 00691 mat Gskew = pow(Ysub, 2); 00692 vec Beta = sumcol(elem_mult(Ysub, Gskew)); 00693 mat D = diag(pow(Beta , -1)); 00694 B = B + myy * B * (transpose(Ysub) * Gskew - diag(Beta)) * D; 00695 break; 00696 } 00697 00698 00699 } // SWITCH usedNlinearity 00700 00701 } // FOR maxIterations 00702 00703 W = transpose(B) * whiteningMatrix; 00704 00705 00706 } // IF FICA_APPROACH_SYMM APPROACH 00707 00708 // DEFLATION 00709 else { 00710 00711 // FC 01/12/05 00712 A = zeros(whiteningMatrix.cols(), numOfIC); 00713 // A = zeros( vectorSize, numOfIC ); 00714 mat B = zeros(vectorSize, numOfIC); 00715 W = transpose(B) * whiteningMatrix; 00716 int round = 1; 00717 int numFailures = 0; 00718 00719 while (round <= numOfIC) { 00720 00721 myy = myyOrig; 00722 00723 usedNlinearity = gOrig; 00724 stroke = 0; 00725 00726 notFine = 1; 00727 loong = 0; 00728 int endFinetuning = 0; 00729 00730 vec w = zeros(vectorSize); 00731 00732 if (initialStateMode == 0) 00733 00734 w = randu(vectorSize) - 0.5; 00735 00736 else w = whiteningMatrix * guess.get_col(round); 00737 00738 w = w - B * transpose(B) * w; 00739 00740 w /= norm(w); 00741 00742 vec wOld = zeros(vectorSize); 00743 vec wOld2 = zeros(vectorSize); 00744 00745 int i = 1; 00746 int gabba = 1; 00747 00748 while (i <= maxNumIterations + gabba) { 00749 00750 w = w - B * transpose(B) * w; 00751 00752 w /= norm(w); 00753 00754 if (notFine) { 00755 00756 if (i == maxNumIterations + 1) { 00757 00758 round--; 00759 00760 numFailures++; 00761 00762 if (numFailures > failureLimit) { 00763 00764 if (round == 0) { 00765 00766 A = dewhiteningMatrix * B; 00767 W = transpose(B) * whiteningMatrix; 00768 00769 } // IF round 00770 00771 break; 00772 00773 } // IF numFailures > failureLimit 00774 00775 break; 00776 00777 } // IF i == maxNumIterations+1 00778 00779 } // IF NOTFINE 00780 00781 else if (i >= endFinetuning) wOld = w; 00782 00783 if (norm(w - wOld) < epsilon || norm(w + wOld) < epsilon) { 00784 00785 if (fineTuningEnabled && notFine) { 00786 00787 notFine = 0; 00788 gabba = maxFinetune; 00789 wOld = zeros(vectorSize); 00790 wOld2 = zeros(vectorSize); 00791 usedNlinearity = gFine; 00792 myy = myyK * myyOrig; 00793 endFinetuning = maxFinetune + i; 00794 00795 } // IF finetuning 00796 00797 else { 00798 00799 numFailures = 0; 00800 00801 B.set_col(round - 1, w); 00802 00803 A.set_col(round - 1, dewhiteningMatrix*w); 00804 00805 W.set_row(round - 1, transpose(whiteningMatrix)*w); 00806 00807 break; 00808 00809 } // ELSE finetuning 00810 00811 } // IF epsilon 00812 00813 else if (stabilizationEnabled) { 00814 00815 if (stroke == 0.0 && (norm(w - wOld2) < epsilon || norm(w + wOld2) < epsilon)) { 00816 00817 stroke = myy; 00818 myy /= 2.0 ; 00819 00820 if (mod(usedNlinearity, 2) == 0) { 00821 00822 usedNlinearity++; 00823 00824 } // IF MOD 00825 00826 }// IF !stroke 00827 00828 else if (stroke != 0.0) { 00829 00830 myy = stroke; 00831 stroke = 0.0; 00832 00833 if (myy == 1 && mod(usedNlinearity, 2) != 0) { 00834 usedNlinearity--; 00835 } 00836 00837 } // IF Stroke 00838 00839 else if (notFine && !loong && i > maxNumIterations / 2) { 00840 00841 loong = 1; 00842 myy /= 2.0; 00843 00844 if (mod(usedNlinearity, 2) == 0) { 00845 00846 usedNlinearity++; 00847 00848 } // IF MOD 00849 00850 } // IF notFine 00851 00852 } // IF stabilization 00853 00854 00855 wOld2 = wOld; 00856 wOld = w; 00857 00858 switch (usedNlinearity) { 00859 00860 // pow3 00861 case FICA_NONLIN_POW3 : { 00862 w = (X * pow(transpose(X) * w, 3)) / numSamples - 3 * w; 00863 break; 00864 } 00865 case(FICA_NONLIN_POW3+1) : { 00866 vec Y = transpose(X) * w; 00867 vec Gpow3 = X * pow(Y, 3) / numSamples; 00868 double Beta = dot(w, Gpow3); 00869 w = w - myy * (Gpow3 - Beta * w) / (3 - Beta); 00870 break; 00871 } 00872 case(FICA_NONLIN_POW3+2) : { 00873 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00874 w = (Xsub * pow(transpose(Xsub) * w, 3)) / Xsub.cols() - 3 * w; 00875 break; 00876 } 00877 case(FICA_NONLIN_POW3+3) : { 00878 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00879 vec Gpow3 = Xsub * pow(transpose(Xsub) * w, 3) / (Xsub.cols()); 00880 double Beta = dot(w, Gpow3); 00881 w = w - myy * (Gpow3 - Beta * w) / (3 - Beta); 00882 break; 00883 } 00884 00885 // TANH 00886 case FICA_NONLIN_TANH : { 00887 vec hypTan = tanh(a1 * transpose(X) * w); 00888 w = (X * hypTan - a1 * sum(1 - pow(hypTan, 2)) * w) / numSamples; 00889 break; 00890 } 00891 case(FICA_NONLIN_TANH+1) : { 00892 vec Y = transpose(X) * w; 00893 vec hypTan = tanh(a1 * Y); 00894 double Beta = dot(w, X * hypTan); 00895 w = w - myy * ((X * hypTan - Beta * w) / (a1 * sum(1 - pow(hypTan, 2)) - Beta)); 00896 break; 00897 } 00898 case(FICA_NONLIN_TANH+2) : { 00899 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00900 vec hypTan = tanh(a1 * transpose(Xsub) * w); 00901 w = (Xsub * hypTan - a1 * sum(1 - pow(hypTan, 2)) * w) / Xsub.cols(); 00902 break; 00903 } 00904 case(FICA_NONLIN_TANH+3) : { 00905 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00906 vec hypTan = tanh(a1 * transpose(Xsub) * w); 00907 double Beta = dot(w, Xsub * hypTan); 00908 w = w - myy * ((Xsub * hypTan - Beta * w) / (a1 * sum(1 - pow(hypTan, 2)) - Beta)); 00909 break; 00910 } 00911 00912 // GAUSS 00913 case FICA_NONLIN_GAUSS : { 00914 vec u = transpose(X) * w; 00915 vec Usquared = pow(u, 2); 00916 vec ex = exp(-a2 * Usquared / 2); 00917 vec gauss = elem_mult(u, ex); 00918 vec dGauss = elem_mult(1 - a2 * Usquared, ex); 00919 w = (X * gauss - sum(dGauss) * w) / numSamples; 00920 break; 00921 } 00922 case(FICA_NONLIN_GAUSS+1) : { 00923 vec u = transpose(X) * w; 00924 vec Usquared = pow(u, 2); 00925 00926 vec ex = exp(-a2 * Usquared / 2); 00927 vec gauss = elem_mult(u, ex); 00928 vec dGauss = elem_mult(1 - a2 * Usquared, ex); 00929 double Beta = dot(w, X * gauss); 00930 w = w - myy * ((X * gauss - Beta * w) / (sum(dGauss) - Beta)); 00931 break; 00932 } 00933 case(FICA_NONLIN_GAUSS+2) : { 00934 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00935 vec u = transpose(Xsub) * w; 00936 vec Usquared = pow(u, 2); 00937 vec ex = exp(-a2 * Usquared / 2); 00938 vec gauss = elem_mult(u, ex); 00939 vec dGauss = elem_mult(1 - a2 * Usquared, ex); 00940 w = (Xsub * gauss - sum(dGauss) * w) / Xsub.cols(); 00941 break; 00942 } 00943 case(FICA_NONLIN_GAUSS+3) : { 00944 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00945 vec u = transpose(Xsub) * w; 00946 vec Usquared = pow(u, 2); 00947 vec ex = exp(-a2 * Usquared / 2); 00948 vec gauss = elem_mult(u, ex); 00949 vec dGauss = elem_mult(1 - a2 * Usquared, ex); 00950 double Beta = dot(w, Xsub * gauss); 00951 w = w - myy * ((Xsub * gauss - Beta * w) / (sum(dGauss) - Beta)); 00952 break; 00953 } 00954 00955 // SKEW 00956 case FICA_NONLIN_SKEW : { 00957 w = (X * (pow(transpose(X) * w, 2))) / numSamples; 00958 break; 00959 } 00960 case(FICA_NONLIN_SKEW+1) : { 00961 vec Y = transpose(X) * w; 00962 vec Gskew = X * pow(Y, 2) / numSamples; 00963 double Beta = dot(w, Gskew); 00964 w = w - myy * (Gskew - Beta * w / (-Beta)); 00965 break; 00966 } 00967 case(FICA_NONLIN_SKEW+2) : { 00968 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00969 w = (Xsub * (pow(transpose(Xsub) * w, 2))) / Xsub.cols(); 00970 break; 00971 } 00972 case(FICA_NONLIN_SKEW+3) : { 00973 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize)); 00974 vec Gskew = Xsub * pow(transpose(Xsub) * w, 2) / Xsub.cols(); 00975 double Beta = dot(w, Gskew); 00976 w = w - myy * (Gskew - Beta * w) / (-Beta); 00977 break; 00978 } 00979 00980 } // SWITCH nonLinearity 00981 00982 w /= norm(w); 00983 i++; 00984 00985 } // WHILE i<= maxNumIterations+gabba 00986 00987 round++; 00988 00989 } // While round <= numOfIC 00990 00991 } // ELSE Deflation 00992 00993 } // FPICA
Generated on Sat Jul 9 2011 15:21:32 for IT++ by Doxygen 1.7.4