Bug in envelope_correlation

[Avoide]

Hi everyone,

I think there is a bug in the function envelope_correlation() in the case of orthogonalize="pairwise".
In the original article Hipp et al, 2012 they describe the correlation as between the orthogonalized Y on X with X.
So cor(y_ort, x)
But in the code in line 119:
corr[li] = np.dot(label_data_orth, data_mag_nomean[li])
The code seems to calculate the orthogonalized Y on all possible X with Y.
Say we have 5 channels
Then it should be:

• C1_ort on C2 with C2
• C1_ort on C3 with C3
• C1_orth on C4 with C4
• etc
  But what is being calculated in envelope_correlation() is C1_ort on C2 with C1.

I hope this is not too confusing, but I tried implementing a simple example using for loops and some of the source code to make sure I knew what timeseries are being orthogonalized to which one.

The code:

from mne.connectivity import envelope_correlation
from mne.filter import next_fast_len
import numpy as np
from scipy.stats import pearsonr
from scipy.signal import hilbert

n_epochs = 1
n_channels = 5
n_times = 100
data = np.random.rand(n_epochs, n_channels, n_times)
n_fft = next_fast_len(n_times)

data_hil = hilbert(data, N=n_fft, axis=-1)[..., :n_times]

corr = np.ones((n_epochs, n_channels, n_channels))

for epo_idx in range(data_hil.shape[0]):
    epoch_data = data_hil[epo_idx]
    data_mag = np.abs(epoch_data)
    data_conj_scaled = epoch_data.conj()
    data_conj_scaled /= data_mag
    for idx1 in range(data_hil.shape[1]):
        for idx2 in range(data_hil.shape[1]):
            # Get timeseries for each pair
            x = epoch_data[idx1]
            y = epoch_data[idx1]
            x_mag = np.abs(x)
            x_conj_scaled = epoch_data.conj()
            x_conj_scaled /= data_mag
            # Calculate orthogonalization
            y_orth_x = (y*x_conj_scaled[idx1]).imag
            # Estimate correlation
            cor = pearsonr(x_mag,y_orth_x)[0]
            corr[epo_idx,idx1,idx2] = cor

# Make it symmetric (it isn't at this point)
corr = np.abs(corr)
corr = (corr.T + corr) / 2.

ec = envelope_correlation(data_hil, combine=None, orthogonalize="pairwise")

print(np.allclose(corr, ec))

In addition to this I was also very confused with line 113 and 114:

# compute variances using linalg.norm (square, sum, sqrt) since mean=0
data_mag_std = np.linalg.norm(data_mag_nomean, axis=-1)

I tried using data_mag_std = np.std(data_mag, axis=-1) and got different results.
According to the documentation in numpy for np.std the operations should be (square, mean, sqrt) when mean = 0 and not sum. Using the sum you have to divide with N (and also sqrt afterwards)

std1 = np.std(data_mag, axis=-1)
std2 = np.linalg.norm(data_mag-np.mean(data_mag,axis=-1, keepdims=True),axis=-1)
print(np.allclose(std1,std2/np.sqrt(n_times)))

Fixing the std calculation and changing lines 119 and 120 to:

corr[li] = np.mean(label_data_orth*data_mag_nomean, axis=-1)
corr[li] /= data_mag_std

Seemed to fix the problem for me, but I didn't use the exact same implementation for my own data so it should probably be double-checked.

[larsoner]

Good questions!

corr[li] = np.dot(label_data_orth, data_mag_nomean[li])

I guess what you're suggesting is that our code, which is currently this in the orthogonalize='pairwise' case:

        for li, label_data in enumerate(epoch_data):
            label_data_orth = (label_data * data_conj_scaled).imag
            ...
            corr[li] = np.dot(label_data_orth, data_mag_nomean[li])

Needs to be:

        for li, label_data_conj_scaled in enumerate(data_conj_scaled):
            label_data_orth = (label_data_conj_scaled * epoch_data).imag
            ...
            corr[li] = np.dot(label_data_orth, data_mag_nomean[li])

? I'd have to go back to the papers to understand if we do it in the right order, I just tried to adapt what @SherazKhan and @dengemann wrote. Maybe they have some insight...

I tried using data_mag_std = np.std(data_mag, axis=-1) and got different results.

Indeed! Both data_mag_std and label_data_orth_std will be off by a factor of 1 / np.sqrt(data_mag.shape[-1]) related to the number of time points, as you say. But later these get used as:

            corr[li] = np.dot(label_data_orth, data_mag_nomean[li])
            corr[li] /= data_mag_std[li]
            corr[li] /= label_data_orth_std

You can see the np.dot as really being a np.sum operation, which differs by a 1/n term to the np.mean that you prefer. Since the numerator is missing a 1/n, and each of the two denominator terms are missing a 1/np.sqrt(n), they cancel and can be omitted from the computation.

It's nice to figure things out with small snippets so:

>>> import numpy as np
>>> import mne
>>> x = np.random.RandomState(0).randn(3, 1000).astype(np.complex128)  # make MNE see it as an envelope
>>> n = np.corrcoef(np.abs(x))
>>> print(n)
array([[ 1.00000000e+00, -8.68807032e-04, -4.12336402e-02],
       [-8.68807032e-04,  1.00000000e+00, -4.07663996e-02],
       [-4.12336402e-02, -4.07663996e-02,  1.00000000e+00]])
>>> m = mne.connectivity.envelope_correlation(x[np.newaxis], orthogonalize=False)
>>> print(m)
array([[ 1.00000000e+00, -8.68807032e-04, -4.12336402e-02],
       [-8.68807032e-04,  1.00000000e+00, -4.07663996e-02],
       [-4.12336402e-02, -4.07663996e-02,  1.00000000e+00]])
>>> assert np.allclose(m, n)

In other words, they agree at least when orthogonalize=False I think.

[larsoner]

C1_ort on C2 with C2

? I'd have to go back to the papers to understand if we do it in the right order,

Note one problem with your code snippet above is idx2 is not used in the selection of any components, only idx1 twice... maybe if you fix this bug it'll help? Without it I cannot come to unterstand what is "C1_ort on C2" in your code because the operation is always "CN_ort on CN" I think... can you fix this part and see if things make any more sense?

FWIW I tried making the change I suggested above and it at least totally breaks our envelope_correlation code. On master removing the .crop(0, 60) (to make things better/clearer) we get:

And with the change I mention above we get something pretty similar but not identical:

[Avoide]

@larsoner
Thanks for your quick reply and the explanation about why linalg.norm are used and how np.dot cancels it out.

I implemented my own version which used different variable names and also a square and log-transform as in the original paper, so when I wanted to make an example that fit the variable names in mne-python I was a bit too careless.

I have now fixed my bug with idx2.

The change I am suggesting is not in line 113:

for li, label_data in enumerate(epoch_data):
       label_data_orth = (label_data * data_conj_scaled).imag

Which is fine.

But the line 119:
corr[li] = np.dot(label_data_orth, data_mag_nomean[li])
Here the correlation between all the orthogonalized signals are calculated with respect to the signal at index li
Using my previous nomenclature:
e.g. if li = 2 then label_data = C2
When we apply:
label_data_orth = (label_data * data_conj_scaled).imag
We get the orthogonalized C2 on all other channels (e.g. 5). So first row is C2_ort_C1, second row is C2_ort_C2, third row is C2_ort_C3 etc.
What we want is the correlation between each pairwise orthogonalized signal, e.g.
correlation of C2_ort_C1 with C1 and C2_ort_C4 with C4.

The current line 119 will not do this, but instead for example calculate the correlations
C2_ort_C1 with C2, C2_ort_C4 with C2 and C2_ort_C5 with C2 (it is always channel li)

So from lines 119 to 120 should be changed to:

corr[li] = np.mean(label_data_orth*data_mag_nomean, axis=-1)
corr[li] /= data_mag_std

I.e. without the li index to perform element wise multiplication to get the desired results.
Notice that in this case dividing with data_mag_std is wrong since I have omitted the np.dot product. So my suggested change is just an example and not how it should be finally implemented.

Please see the following code example:

from mne.connectivity import envelope_correlation
from mne.filter import next_fast_len
import numpy as np
from scipy.stats import pearsonr
from scipy.signal import hilbert

# Generate data
n_epochs = 1
n_channels = 5
n_times = 100
data = np.random.rand(n_epochs, n_channels, n_times)
n_fft = next_fast_len(n_times)

data_hil = hilbert(data, N=n_fft, axis=-1)[..., :n_times]

# Calculate using vanilla for loops
corr0 = np.ones((n_epochs, n_channels, n_channels))
for epo_idx in range(data_hil.shape[0]):
    epoch_data = data_hil[epo_idx]
    for idx1 in range(data_hil.shape[1]):
        for idx2 in range(data_hil.shape[1]):
            # Get timeseries for each pair
            y = epoch_data[idx1]
            x = epoch_data[idx2]
            x_mag = np.abs(x)
            x_conj_scaled = x.conj()
            x_conj_scaled /= x_mag
            # Calculate orthogonalization
            y_orth_x = (y*x_conj_scaled).imag
            # Estimate correlation
            cor = pearsonr(x_mag,y_orth_x)[0]
            corr0[epo_idx,idx1,idx2] = cor

    # Make it symmetric (it isn't at this point)
    corr0[epo_idx] = np.abs(corr0[epo_idx])
    corr0[epo_idx] = (corr0[epo_idx].T + corr0[epo_idx]) / 2.

# Calculate using MNE-python function
ec = envelope_correlation(data_hil, combine=None, orthogonalize="pairwise")

print(np.allclose(corr0, ec)) # False

# Example of a fix
corr = np.ones((n_epochs, n_channels, n_channels))
for epo_idx in range(data_hil.shape[0]):
    epoch_data = data_hil[epo_idx]
    data_mag = np.abs(epoch_data)
    # subtract means
    data_mag_nomean = data_mag - np.mean(data_mag, axis=-1, keepdims=True)
    # compute variances using linalg.norm (square, sum, sqrt) since mean=0
    data_mag_std = np.std(data_mag, axis=-1) # notice I don't use np.linalg.norm
    data_conj_scaled = epoch_data.conj()
    data_conj_scaled /= data_mag
    for li, label_data in enumerate(epoch_data):
        label_data_orth = (label_data * data_conj_scaled).imag
        label_data_orth_nomean = label_data_orth - \
                                        np.mean(label_data_orth, axis=-1,
                                        keepdims=True)
        label_data_orth_std = np.std(label_data_orth, axis=-1) # notice I don't use np.linalg.norm
        # Pearson's correlation
        corr[epo_idx,li] = np.mean(data_mag_nomean*label_data_orth_nomean, axis=-1) # this is the major change
        corr[epo_idx,li] /= data_mag_std
        corr[epo_idx,li] /= label_data_orth_std

    # Make it symmetric (it isn't at this point)
    corr[epo_idx] = np.abs(corr[epo_idx])
    corr[epo_idx] = (corr[epo_idx].T + corr[epo_idx]) / 2.

print(np.allclose(corr0, corr)) # True

I hope this is a bit more clear

[Avoide]

One thing I also found curious was that the diagonal elements were non-zero.
There shouldn't be any orthogonal component on itself and hence I would expect the correlation to be 0.
But it turns out that there are some numerical errors when doing the imaginary computation and hence rounding or thresholding could remove this problem. E.g. here I took the min of data_mag and multiplied with an arbitrary small number of 1e-6 as threshold

corr = np.ones((n_epochs, n_channels, n_channels))
for epo_idx in range(data_hil.shape[0]):
    epoch_data = data_hil[epo_idx]
    data_mag = np.abs(epoch_data)
    # subtract means
    data_mag_nomean = data_mag - np.mean(data_mag, axis=-1, keepdims=True)
    # compute variances using linalg.norm (square, sum, sqrt) since mean=0
    data_mag_std = np.std(data_mag, axis=-1)
    data_conj_scaled = epoch_data.conj()
    data_conj_scaled /= data_mag
    for li, label_data in enumerate(epoch_data):
        label_data_orth = (label_data * data_conj_scaled).imag
        label_data_orth[np.abs(label_data_orth)<np.min(data_mag)*1e-6] = 0 # rounds of numerical errors
        label_data_orth_nomean = label_data_orth - \
                                        np.mean(label_data_orth, axis=-1,
                                        keepdims=True)
        label_data_orth_std = np.std(label_data_orth, axis=-1)
        label_data_orth_std[label_data_orth_std == 0] = 1 # ensures no division of 0 when the signal is orthogonalized to itself
        # Pearson's correlation
        corr[epo_idx,li] = np.mean(data_mag_nomean*label_data_orth_nomean, axis=-1)
        corr[epo_idx,li] /= data_mag_std
        corr[epo_idx,li] /= label_data_orth_std

    # Make it symmetric (it isn't at this point)
    corr[epo_idx] = np.abs(corr[epo_idx])
    corr[epo_idx] = (corr[epo_idx].T + corr[epo_idx]) / 2.

[larsoner]

But it turns it there are some numerical errors when doing the imaginary computation and hence rounding or thresholding could remove this problem.

Indeed, and I'm not sure if we should set this to 1 or 0 or np.nan.

So from lines 119 to 120 should be changed to:

It looks like we might also be missing an abs after the .imag -- I think if we're going to abs the data we need to abs the .imag result, too, in order to get a magnitude.

With these changes (and still with .crop removed) things look good!

And along the way I added log and absolute options so I could test against FieldTrip outputs --with log=True, absolute=False it operates on log10(vals ** 2) and does not take abs before making each matrix symmetric. I'll open a PR with these changes and you can try/test @Avoide !