IT++ Logo Newcom Logo

fastica.cpp

Go to the documentation of this file.
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
SourceForge Logo

Generated on Fri Jun 8 01:07:07 2007 for IT++ by Doxygen 1.5.2