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 Sun Dec 9 17:38:49 2007 for IT++ by Doxygen 1.5.4