[Eeglablist] pop_runica rank calculation question

[Sven Hoffmann]

Hi Matthew,

very interesting. You should report this issue to the bugzilla. I do not know why in pop_runica
actually the upper number of ICs is chosen, I haven't ever experienced this problem on linux if
using the code of Aapo, which is mathematically absolutely correct.

Which eeglab version do you use? I cannot replicate this bug (at least on my machine) with the
most current version.
However, some quick thoughts about the nature of ICA and average reference (may be I'm wrong, so
please correct me if so):

1.) Consider what the average reference transformation does: It relates each channel to the
average of all channels, thus inducing that all channels become,  at least to some degree, a
linear combination of all channels. Check reref.m and the algebra will become clear:

%nchansin: number of channels

refmatrix = eye(nchansin)-ones(nchansin)*1/nchansin;
chansout = chansin;
data(chansout,:) = refmatrix*data(chansin,:);

2.)  The rank calculation estimates the number of linearly independent vectors of a matrix. As a
result, the dimensionality of a average referenced data set is reduced. Eigenvalue decomposition
(snippet of Aapo) is only one way to calculate this. The matlab rank function utilizes singular
value decomposition (SVD). Compare the rank.m algorithm with the one of Aapo:

s = svd(A);
tol = max(size(A))*eps(max(s));
r = sum(s > tol);

Decreasing the rank tolerance in the rank.m function of matlab  should yield a comparable result
like Aapo's snippet.


3.) Now have a look a look at the "backprojection" (from pop_subcomp.m):


compproj = EEG.icawinv(:, component_keep)*eeg_getdatact(EEG, 'component', component_keep,
'reshape', '2d');

which means:

cleaned data= W(:, comps to keep) * S(comps to keep,:)

4.) ICA yields independent components. Also, prior ICA (at least prior infomax) the data are
sphered, which means that they they are "decorrelated'.

The reason for "loosing track" might be thus  3.) and 4.): The backprojected data are linear
combinations of independent components, which are reference free. As a result the linear
combination of them will produce a matrix with a rank according to the number of components
which were kept (varying with the value for the rank tolerance). A comparable experiment can be
made by removing some channels from an average referenced data set. To my mind, this isn't
really a problem, as long as you do not re-reference again.



Best,
Sven Hoffmann.




Am 18.02.2012 02:30, schrieb Matthew Stief:
> Greetings!
>
> I have a question regarding the way that pop_runica calculates rank.  It does not seem to be
> doing so correctly for us.
>
> I have been experimenting, and when I try to run it on relatively unprocessed data, it seems
> to work properly, throws up no warning messages, and tries to calculate components equal to
> the number of electrodes.  However, if I reduce the rank of the data by referencing to the
> average electrode, it throws up this error:
>
> Warning: fixing rank computation inconsistency (132 vs 131) most likely because running under
> Linux 64-bit Matlab...
> Finding 132 ICA components using extended ICA.
>
> Now I have seen this error discussed on the list, but in this case I am sure that it is
> wrong.  I know what I did to the dataset (i.e. changed to the average reference), and I know
> that the correct rank should be 131.  I'm also not running in Linux, I'm running in Windows 7
> with 64-bit Matlab R2011b 7.13.0.564.
>
> Now I've familiarized myself with the bug reported by Sven Hoffman on the issues with the rank
> function and double versus single precision, and switched to double precision in the memory
> options, and I notice that in pop_runica the correct rank of 131 is being calculated by the
> code he borrowed from fastica.
>
> *   covarianceMatrix = cov(tmpdata', 1);
>    [E, D] = eig (covarianceMatrix);
>    rankTolerance = 1e-7;
>    tmprank2=sum (diag (D) > rankTolerance);
>
> *And that this is then compared with the output of the regular rank function and whichever is
> largest is chosen, with a reference to the Linux issue.  However, of course in this case the
> correct rank is actually the smaller of the two!
>
> Of course it would be simple to just change the code to instead select the smaller of the two
> numbers, but before changing pop_runica I would love to get the advice of the experts.
>
> There is an additional issue as well, which is that if I reduce the rank by changing to
> average reference, and then reduce the rank again by removing artefactual components,
> Hoffman's code seems to lose track of the rank reduction from the average reference and
> returns a rank as if only the components have been removed!  This double-ICA approach isn't
> really the one that I'm going to take with my data, but I'd still like to understand how and
> why incorrect ranks are being calculated.
>
> I'm guessing that the general solution to this issue will be to just keep manual track of how
> much rank has been lost on theoretical grounds, but it would certainly be nice if there were
> an automatic rank detection method that worked.
>
> -- 
> _________________________________________________________________
> Matthew Stief
> Human Development | Sex & Gender Lab | Cornell University
> http://www.human.cornell.edu/HD/sexgender
>
>
> Heterosexuality isn't normal, it's just common.
> -Dorothy Parker
>
>
> _______________________________________________
> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> To unsubscribe, send an empty email to eeglablist-unsubscribe at sccn.ucsd.edu
> For digest mode, send an email with the subject "set digest mime" to eeglablist-request at sccn.ucsd.edu

-- 
*Dr. rer. nat. Sven Hoffmann*
Leibniz-Institut für Arbeitsforschung an der TU Dortmund
/Leibniz Research Centre for Working Environment and Human Factors/
Ardeystr. 67
D-44139 Dortmund
Phone: +49-231-1084-265
Fax: +49-231-1084-401
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20120228/578b588e/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 554 bytes
Desc: OpenPGP digital signature
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20120228/578b588e/attachment.sig>