[FieldTrip] More memory-efficient dwPLI connectivity computation algorithm

Schoffelen, J.M. (Jan Mathijs) janmathijs.schoffelen at donders.ru.nl
Wed Oct 6 17:26:01 CEST 2021


Dear Stefan,

Indeed it can be challenging if you’re pressed for RAM to do some of the larger computations. Specifically, as you describe, the wpli computation is inefficient since the current code requires the single observation cross-spectra. Historically, the modular design of ft_connectivityanalysis kind of required for the lower-level connectivity metric computation machinery to receive the bivariate data (crossspctra) in their input. In the early development of the function (starting off from coherence), this restriction was not problematic per se.  For the coherence computation one can easily (and early) drop the observation dimension, so the memory requirements are relatively low. I am fine with adjusting ft_connectivity_wpli in such a way that it also can take fourier data in the input, and that - after a check - the conjugate multiplication will be done only within the function. If that is in place, we can adjust ft_connectivityanalysis to allow the fourierspectra in the input. I’d be happy to guide you through the process of making the changes. Can I suggest you to start this as a pull request on github?

Best wishes and keep up the good work,

Jan-Mathijs



On 5 Oct 2021, at 10:34, Stefan Dvoretskii via fieldtrip <fieldtrip at science.ru.nl<mailto:fieldtrip at science.ru.nl>> wrote:

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<https://urldefense.com/v3/__https://www.dropbox.com/s/cmrv6zc4rsawfvh/virtFreq_example_one_foi.mat?dl=0__;!!HJOPV4FYYWzcc1jazlU!toh8VSFl1ftbd3jDBT-ghnM1lTJvdAF3td8rm6geQk97F7Wi0ueMVta-BQr3afitDCAvlxMH3Lh1neo$>
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
_______________________________________________
fieldtrip mailing list
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
https://urldefense.com/v3/__https://doi.org/10.1371/journal.pcbi.1002202__;!!HJOPV4FYYWzcc1jazlU!toh8VSFl1ftbd3jDBT-ghnM1lTJvdAF3td8rm6geQk97F7Wi0ueMVta-BQr3afitDCAvlxMHMce6rH0$

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20211006/4a0ada77/attachment.htm>


More information about the fieldtrip mailing list