% generate 3 random super-Gaussian "EEG sources"
S=randn(3,100000).^2;
% mix the sources to produce the "EEG observations"
A = rand(3);
X = A*S;
% the covariance matrix of our data is:
C = cov(X');
% we can check that this covariance is well-conditioned (i.e. our data is not
% rank defitient because rcond(C)>>0
rcond(C); % significantly greater than 0
% 1st ICA run
[weights,sphere] = runica(X);
% the estimated separating and mixing matrix in the first ica run
W1 = weights*sphere;
A1 = inv(W1);
% the estimated sources
S1 = W1*X;
% we reject the third component (assume it is an EOG component)
S1 = W1(1:2,:)*X;
% so that after projecting back to the electrodes, the filtered EEG is:
X1 = A1(:,1:2)*S1;
% now we can check that our "filtered EEG" is rank deficient, i.e. its
% covariance matrix is ill-conditioned:
C1 = cov(X1');
rcond(C1) % very close to zero (eps)
% if we now attempt to run ICA on this rank-deficient data:
[weights,sphere] = runica(X1); 
% you better press CTRL+C because the command above will run forever, make
% MATLAB crash, produce imaginary numbers or give you very incaccurate results

% SOLUTION:
% do PCA before running ICA for the second time
[V,D] = eig(C1);
% observe that the smallest eigenvalue in the diagonal matrix D is almost
% zero so we can reconstruct the original data almost perfectly using only
% two principal components. The PCA transformation matrix will be:
Wpca = (V(:,2:3)*D(2:3,2:3)^(-1/2))';
% so that the two principal components will be
X1pca = Wpca*X1;
% now we can safely run ICA on this principal components
[weights,sphere] = runica(X1pca); 
% so that the actual estimated separating matrix in the second ICA run is:
W2 = weights*sphere*Wpca; 
% the estimated sources in the second ICA run are:
S2 = W2*X1;
% and the estimated mixing matrix is:
A2 = pinv(W2); %note that I use pinv because W2 is not a square matrix
% we can easily check that this second ICA run did not improve our results.
% We just obtained the same two sources that we kept in the first ICA run:
W2*A1(:,1:2)  
% since this matrix is almost a permutation (and/or sign-change) matrix
% the estimated components in the second ICA run are just a permuted (and
% possibly sign-changed) version of S1. Note that: S2=W2*A1(:,1:2)*S1;