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