ENH: compute CSD directly for upper-triangle channel pairs in Fourier/multitaper

[PragnyaKhandelwal]

Reference issue

Fixes #13717

What does this implement/fix?

This PR addresses the FIXME in mne/time_frequency/csd.py by avoiding full channel × channel CSD matrix construction in:

• _csd_fourier
• _csd_multitaper

Instead, it computes CSD directly for upper-triangle channel pairs (np.triu_indices), which matches the internal vectorized storage format and avoids unnecessary intermediate allocations.

Additional information

• Motivation is the existing FIXME around full-matrix computation followed by _sym_mat_to_vector conversion.
• Benchmarks were run on main vs this branch with a self-contained MWE:
  • channel counts tested: 128 and 192
  • both Fourier and multitaper paths show lower peak memory on this branch
  • multitaper also shows clear runtime improvement; Fourier also improved in larger-channel runs
• Output shapes remain identical across branches for equivalent configs.
• Local checks run for this change:
  • pytest mne/time_frequency/tests/test_csd.py -k "csd_fourier or csd_multitaper" -q
  • pre-commit run --files mne/time_frequency/csd.py

[PragnyaKhandelwal]

I looked into the failed Azure job.
From the logs, it seems the failure is happening in the Codecov upload step (bash <(curl -s https://codecov.io/bash))
So this looks like an intermittent CI/network/script-fetch issue....
Would you mind re-running the failed check once?

[larsoner]

Changelog failure is legitimate, you should add a doc/changes/devel/13719.newfeature.rst using the :newcontrib: role and an entry in doc/changes/names.inc

[scott-huberty]

I am not super familiar with csd but assuming that all the unique information in the matrix can be extracted from the upper triangle, then it makes sense to me to do the computation on the upper triangle directly. And the aforementioned FIXME comment seems to acknowledge that this is the case.

Tests are passing and a local micro-benchmark suggests we do get a speedup:

code

from functools import partial
from timeit import repeat

import numpy as np

import mne
from mne.time_frequency import csd_array_fourier, csd_array_multitaper
from mne.time_frequency.tests.test_csd import test_csd_fourier, test_csd_multitaper

mne.set_log_level("WARNING")

rng = np.random.default_rng(0)
n_epochs = 10
n_times = 1000
n_channels = 64
X = rng.standard_normal((n_epochs, n_channels, n_times))
kwargs = dict(sfreq=100, fmin=8, fmax=30, n_jobs=1)

def benchmark(func, n=100, n_repeats=5):
    return repeat(partial(func, X, **kwargs), number=n, repeat=n_repeats)

n = 5
repeats = 3
funcs = {
    "fourier": csd_array_fourier,
    "multitaper": csd_array_multitaper,
}

for name, func in funcs.items():
    times = benchmark(func, n=n, n_repeats=repeats)
    best = min(times)
    per_loop = best / n
    print(f"{name.upper()}: {n} loops per run, best of {repeats} runs: {per_loop:.3f} s per loop")

Main	PR
FOURIER: 5 loops per run, best of 3 runs: 0.438 s per loop	FOURIER: 5 loops per run, best of 3 runs: 0.254 s per loop
MULTITAPER: 5 loops per run, best of 3 runs: 1.511 s per loop	MULTITAPER: 5 loops per run, best of 3 runs: 1.044 s per loop

@wmvanvliet GitLens seems to suggest that you last refactored this code (albeit a long time ago!). Do you have the time to take a quick pass at the diff?

[PragnyaKhandelwal]

Done @larsoner! I have added the changelog entry at doc/changes/devel/13719.newfeature.rst using the :newcontrib: role as requested. I’ve also added my name to doc/changes/names.inc.

All checks are now passing, including the Towncrier verification. This PR is ready for a final look and merge!

[PragnyaKhandelwal]

Hi @tsbinns , just a quick heads-up that I’ve updated the towncrier entry to include the requested files. All CI checks are green now—ready for a final look and merge whenever you have a moment! Thanks!