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