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