[FieldTrip] More memory-efficient dwPLI connectivity computation algorithm

[Stefan Dvoretskii]

Dear Fieldtrip community,

my name is Stefan Dvoretskii and I am working in the PainLabMunich on a
pain-related data preprocessing pipeline for EEG. Currently, I am analyzing
connectivity in resting-state EEG data. In this e-mail, I would like to
suggest an alternative, more memory-efficient algorithm for the dwPLI
computation.

During the analysis, we have faced memory problems when trying to compute
the debiased weight phase lag index (dwPLI) on source-level data. We tried
to compute a connectivity matrix from band passed ‘virtual channel data’
(2020 voxels x 500 samples x 299 trials) after spatial filtering with this
code:

[Snippet 1]

fois = 8:0.5:13; % Alpha band with a freq. resolution of 0.5

connMatrix = zeros(2020,2020,length(fois));



for idx = 1:length(fois)

    % Fourier components

    cfg = [];

    cfg.method = 'mtmfft';

    Cfg.output = 'fourier';

    cfg.keeptrials = 'yes';

    cfg.foi = fois(idx);

    cfg.tapsmofrq = 1;

    virtFreq = ft_reqanalysis(cfg,virtChan_data);



    % Connectivity

    cfg = [];

    cfg.method = 'wpli_debiased';

    source_conn = ft_connectiviyanalysis(cfg,virtFreq);

    connMatrix(:,:,idx) = source_conn.(['wpli_debiased','spctrm']);

end

Matlab was just being killed by the system attempting to take too much
memory, without any error message. By debugging, I figured that in our case
the problem was a check of the data in the function ft_connectivityanalysis
(line 251):

data = ft_checkdata(data,'datatype',{'freqmvar' 'freq'});

The cross-spectral density was computed resulting in an array of
2020*2020*300 complex positions for one recording, totaling to 10GB of
working memory. It was too much for the current working machine.

I have addressed this memory problem by creating a function combining the
ft_connectivityanalysis() and ft_connectivity_wpli() that avoids the
memory-costly precomputation of the CSD and computes the cross-spectrum
from the Fourier spectrum. It processes one frequency of interest at a go
[Snippet 2]:

[Snippet 2]

%% Compute dwPLI using only Fourier spectrum to save space

% input:

%       - data: virtual channels frequency decomposition for one frequency
of interest (e.g. 8.5 Hz). Output of ft_frequencyanalysis; includes Tapers
x Channels Fourier spectrum as a field.

function connMatrix = conn_dwpli_fourier(data)

    nchan = size(data.label,1);

    nrpt = size(data.cumtapcnt,1);

    sumtapcnt = [0;cumsum(data.cumtapcnt(:,1))];



    % see ft_connectivity_wpli.m

    outsum = complex(zeros(nchan, nchan));

    outsumW = zeros(nchan, nchan);

    outssq = complex(zeros(nchan, nchan));

    % alternatively: loop along channel combinations. reduces the csdTrial
size

    % (converts to the float), but increases the amount of operations. good
for

    % cluster/peer network jobs.

    for p=1:nrpt

        indx=(sumtapcnt(p)+1):sumtapcnt(p+1);

        fourierTrial = transpose(data.fourierspctrm(indx,:));

        csdTrial = (fourierTrial*fourierTrial')./length(indx);

        % imaginary part (phase info) only

        input = imag(csdTrial);

        outsum = outsum + input;

        outsumW = outsumW + abs(input);

        outssq = outssq + (input.^2);

    end

    connMatrix = (outsum.^2 - outssq)./(outsumW.^2 - outssq);

end

This function [Snippet 2] is to be used in place of ft_connectivityanalysis
in [Snippet 1]. We attach an example of virtChan_data for reproducibility:

https://www.dropbox.com/s/cmrv6zc4rsawfvh/virtFreq_example_one_foi.mat?dl=0

We hope to help people facing similar problems and would be looking forward
to eventual suggestions of integrating this algorithm into fieldtrip
codebase.

Best regards,

Stefan Dvoretskii
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20211005/55a7fccf/attachment.htm>