reduced rank to 1 in mne.minimum_norm.compute_source_psd()

šŸ’¬ Support & Discussions
1

Iā€™m using mne.minimum_norm.compute_source_psd() to compute spectra for 60s-segments of maxfiltered raw data (EEG+MEG, MNE version 1.4.0, Centos 7). Iā€™m currently just looking at the output for the sensor space results.
For some (6 out of 29) datasets I get the message
Reducing data rank 370 ā†’ 1
and indeed the sensor space data have the same time courses across all sensors (EEG and MEG).
When I check ranks using mne.compute_rank(inst, rank=ā€˜infoā€™) everything looks ok.
I included the relevant code snippet and output below.
Any idea what could be wrong with those 6 datasets?

Code snippet:

stc_psd, evo_psd = mne.minimum_norm.compute_source_psd(
raw=raw, inverse_operator=invop, lambda2=1 / 9., method='MNE',
fmin=0., fmax=140., n_fft=60000, overlap=0.5, pca=True,
nave=1, bandwidth='hann', low_bias=True, return_sensor=True)

Output:
Not setting metadata
2 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 1)
1 projection items activated
Considering frequencies 0 ā€¦ 140 Hz
Preparing the inverse operator for useā€¦
Scaled noise and source covariance from nave = 1 to nave = 1
Created the regularized inverter
Created an SSP operator (subspace dimension = 1)
Created the whitener using a noise covariance matrix with rank 134 (236 small eigenvalues omitted)
Picked 370 channels from the data
Computing inverseā€¦
Eigenleads need to be weighted ā€¦
Reducing data rank 370 ā†’ 1
Using hann windowing on at most 2 epochs

2

do you pass rank=ā€œinfoā€ to mne.minimum_norm.make_inverse_operator ?
what do you get with mne.compute_rank(inst, rank=None)?

Alex

3

Yes, I use rank=ā€˜infoā€™ for inverse operator and covariance matrix.
With compute_rank(inst, rank=None) I get {ā€˜megā€™: 306, ā€˜eegā€™: 63}.
I pasted some outputs for my data below.
Thank you,
Olaf

Reading Inverse Operator

Reading inverse operator decomposition from /imaging/hauk/users/olaf/FPVS2/MEG/meg23_113/CBU230249_EEGMEG-inv.fifā€¦
Reading inverse operator infoā€¦
[done]
Reading inverse operator decompositionā€¦
[done]
370 x 370 full covariance (kind = 1) found.
Read a total of 1 projection items:
Average EEG reference (1 x 64) idle
Noise covariance matrix read.
32616 x 32616 diagonal covariance (kind = 2) found.
Source covariance matrix read.
32616 x 32616 diagonal covariance (kind = 6) found.
Orientation priors read.
Did not find the desired covariance matrix (kind = 5)
Did not find the desired covariance matrix (kind = 3)
Reading a source spaceā€¦
Computing patch statisticsā€¦
Patch information addedā€¦
Distance information addedā€¦
[done]
Reading a source spaceā€¦
Computing patch statisticsā€¦
Patch information addedā€¦
Distance information addedā€¦
[done]
2 source spaces read
Read a total of 1 projection items:
Average EEG reference (1 x 64) idle
Source spaces transformed to the inverse solution coordinate frame

mne.compute_rank(raw, info=None)

Computing rank from data with rank=None
Using tolerance 3.9e-11 (2.2e-16 eps * 306 dim * 5.8e+02 max singular value)
Estimated rank (mag + grad): 306
MEG: rank 306 computed from 306 data channels with 0 projectors
:1: RuntimeWarning: Something went wrong in the data-driven estimation of the data rank as it exceeds the theoretical rank from the info (306 > 73). Consider setting rank to ā€œautoā€ or setting it explicitly as an integer.
mne.compute_rank(raw, info=None)
Using tolerance 3e-10 (2.2e-16 eps * 64 dim * 2.1e+04 max singular value)
Estimated rank (eeg): 63
EEG: rank 63 computed from 64 data channels with 1 projector
Out[42]: {ā€˜megā€™: 306, ā€˜eegā€™: 63}

mne.compute_rank(rawface, rank=ā€˜infoā€™)

Computing rank from data with rank=ā€˜infoā€™
MEG: rank 72 after 0 projectors applied to 306 channels
EEG: rank 63 after 1 projector applied to 64 channels
Out[61]: {ā€˜megā€™: 72, ā€˜eegā€™: 63}

4

Iā€™m assuming the crucial code is the following (https://github.com/mne-tools/mne-python/blob/maint/1.4/mne/minimum_norm/time_frequency.py#L742).
This seems to suggest that my inverse operators have one large eigenvalue. I checked coregistration and sensitivity maps, which look ok.

if pca:
        U, s, Vh = linalg.svd(K, full_matrices=False)
        rank = np.sum(s > 1e-8 * s[0])
        K = s[:rank] * U[:, :rank]
        Vh = Vh[:rank]
        logger.info("Reducing data rank %d -> %d" % (len(s), rank))
    else:
        Vh = None
    is_free_ori = inverse_operator["source_ori"] == FIFF.FIFFV_MNE_FREE_ORI
5

mne.minimum_norm.compute_source_psd is one of the oldest functions in
mne. It predates all the rank estimation
cleanum and new similar functions such as
https://mne.tools/stable/generated/mne.minimum_norm.apply_inverse_cov.html
that
allows you to get source power in a frequency by band by computing the
cov on band pass filtered data.

I donā€™t have enough time to look now but I suggest you open an issue
in github so we can actually use correct rank estimation
also in this function.

Alex

6

Thank you. I need a function that provides PSD in sensor and source space. Ideally compute_source_psd() is what Iā€™m looking for.
Iā€™ll open an issue. But does that mean that I canā€™t use compute_source_psd() for maxfiltered data at the moment? Interesting, for 23 out of 29 datasets it seems to work.

Olaf

7

you need to dig into this maybe the 1e-8 constant is not a good guess for your data

Alex

8

Is there a fundamental problem with specifying the rank explicitly, e.g. 60 for MEG, at least as a temporary solution?