IT++ Logo

fastica.cpp

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

Generated on Sat Apr 19 10:41:15 2008 for IT++ by Doxygen 1.5.5