PERF: Default OpenBLAS thread count causes up to 13x regression for decomposition-heavy operations

[sharifhsn]

Summary

Default OpenBLAS thread settings cause severe performance regressions for decomposition-heavy MNE operations on multi-core Linux machines. Shrunk covariance estimation is 13.3x slower, Maxwell filtering (tSSS) is 6.5x slower, and Maxwell SSS is 9.2x slower at the default thread count compared to the optimal setting.

This affects any Linux user with OpenBLAS (the default NumPy BLAS on most systems) and ≥8 cores.

Benchmarks

System: AMD EPYC 7R13, 16 vCPU, OpenBLAS 0.3.31, 60s of sample data.

Operation	BLAS=1	BLAS=2	BLAS=4	BLAS=8	BLAS=16 (default)	Best	Regression
Cov (shrunk)	3.31s	2.67s	2.38s	2.97s	31.58s	4	13.3x
Maxwell SSS	1.02s	0.83s	0.75s	0.72s	6.59s	8	9.2x
Maxwell tSSS	4.42s	2.90s	2.33s	2.40s	15.14s	4	6.5x
apply_inverse_epochs	13.41s	7.99s	5.23s	3.93s	4.16s	8	3.4x
Covariance (emp)	0.58s	0.46s	0.41s	0.43s	0.96s	4	2.3x
make_inverse_op	1.59s	1.20s	1.01s	0.90s	1.11s	8	1.8x
ICA (20 comp)	1.60s	1.18s	1.01s	0.92s	1.14s	8	1.7x
ICA (50 comp)	7.95s	6.80s	6.50s	5.92s	8.23s	8	1.4x
Bandpass filter	0.29s	0.27s	0.29s	0.28s	0.29s	—	1.0x
TFR morlet	6.06s	6.35s	6.19s	6.26s	6.22s	—	1.0x
PSD Welch	0.22s	0.21s	0.22s	0.21s	0.21s	—	1.0x
Source morph	0.67s	0.69s	0.67s	0.66s	0.67s	—	1.0x

Operations fall into three categories:

1. Catastrophic at default threads: Maxwell filter, shrunk covariance — dominated by SVD/QR/eigh on tall-skinny or small matrices where thread synchronization overhead exceeds compute benefit
2. Moderate scaling (peak at 4-8 threads): ICA, inverse operator, empirical covariance
3. BLAS-insensitive: bandpass/notch filter, TFR, PSD, source morph — these are FFT-dominated or sparse and completely unaffected

Root cause

Per-operation benchmarks confirm the issue is specific to decomposition routines on medium-sized matrices:

Operation	BLAS=1	BLAS=4	BLAS=8	BLAS=16
SVD (10000×306)	230ms	120ms	99ms	111ms ↑
QR+pivot (10000×306)	200ms	86ms	71ms	76ms ↑
matmul (306×306)@(306×10000)	39ms	13ms	8ms	7ms
matmul (5000×5000)	4848ms	1229ms	645ms	631ms

SVD and QR on tall-skinny matrices (the kind used in Maxwell tSSS, shrunk covariance, and ICA) peak at 4-8 threads and regress at 16 due to OpenBLAS thread synchronization overhead. Large square matmul scales linearly to 16 threads — but the decomposition-heavy operations don't benefit.

Immediate user workaround

export OPENBLAS_NUM_THREADS=4  # add to .bashrc

This single environment variable provides larger speedups than any code-level optimization for decomposition-heavy workflows.

Possible fixes

1. Documentation: Add a note to the installation/performance docs recommending OPENBLAS_NUM_THREADS=4 for multi-core Linux systems
2. Runtime detection: Use threadpoolctl (already an MNE dependency) to detect and warn when OpenBLAS is running with too many threads
3. Auto-tuning: Wrap decomposition-heavy functions with threadpool_limits(limits=N) to cap BLAS threads — as proposed in Improve control over number of threads used in an mne call #10522 but never implemented

Option 3 is the most robust but also the most invasive. Option 1 is low-effort and immediately helpful.

Reproduction

Benchmark script (requires MNE sample data)

"""Run with: OPENBLAS_NUM_THREADS=N python bench.py"""
import time, warnings, numpy as np
warnings.filterwarnings('ignore')
import mne
from mne.preprocessing import maxwell_filter

sample_path = mne.datasets.sample.data_path()
raw = mne.io.read_raw_fif(
    sample_path / 'MEG/sample/sample_audvis_raw.fif', verbose=False)
raw.crop(tmax=60.).load_data(verbose=False)

# Maxwell tSSS
raw.fix_mag_coil_types()
t0 = time.perf_counter()
maxwell_filter(raw.copy(),
    cross_talk=sample_path / 'SSS/ct_sparse_mgh.fif',
    calibration=sample_path / 'SSS/sss_cal_mgh.dat',
    st_duration=10., verbose=False)
print(f'Maxwell tSSS: {time.perf_counter()-t0:.3f}s')

# Shrunk covariance
raw2 = raw.copy().set_eeg_reference(projection=True, verbose=False)
events = mne.find_events(raw2, verbose=False)
epochs = mne.Epochs(raw2, events, tmin=-0.2, tmax=0.5,
                    baseline=(None, 0), verbose=False, preload=True)
t0 = time.perf_counter()
mne.compute_covariance(epochs, tmax=0., method='shrunk', verbose=False)
print(f'Cov (shrunk): {time.perf_counter()-t0:.3f}s')

Compare: OPENBLAS_NUM_THREADS=4 python bench.py vs OPENBLAS_NUM_THREADS=16 python bench.py

Related

• Improve control over number of threads used in an mne call #10522 — proposed using threadpoolctl to control BLAS threads via n_jobs (closed, never implemented)

[sharifhsn]

End-to-end pipeline impact

Complete MEG preprocessing pipeline (load → Maxwell tSSS → bandpass filter → ICA → epoch → shrunk covariance → inverse operator → source estimation) on 60s of sample data.

System: AMD EPYC 7R13, 16 vCPU, OpenBLAS 0.3.31. Median of 3 runs (warmup excluded).

Full pipeline: fine-grained thread sweep

BLAS= 1   24.4s  ########################
BLAS= 2   21.8s  #####################
BLAS= 3   17.6s  #################
BLAS= 4   16.0s  ################
BLAS= 6   13.7s  #############
BLAS= 7   13.1s  #############
BLAS= 8   12.9s  ############        ← optimum
BLAS= 9   17.3s  #################   ← +35% cliff
BLAS=10   24.2s  ########################
BLAS=12   32.3s  ################################
BLAS=14   41.1s  #########################################
BLAS=16   46.8s  ##############################################  ← default

The pipeline is 3.6x slower at BLAS=16 (default) vs BLAS=8 (optimal). There is a sharp cliff between 8 and 9 threads.

Per-step breakdown at every thread count

BLAS	tSSS	ICA (20)	Cov (shrunk)	Inv operator	Inv epochs
1	3.74s	1.21s	3.47s	1.38s	14.16s
2	5.57s	1.13s	4.59s	2.24s	8.55s
3	4.32s	0.95s	4.13s	1.78s	6.57s
4	3.84s	0.89s	4.15s	1.59s	5.64s
5	3.47s	0.84s	3.81s	1.45s	5.08s
6	3.34s	0.82s	3.86s	1.35s	4.69s
7	3.14s	0.80s	3.53s	1.24s	4.39s
8	3.16s	0.79s	3.51s	1.16s	4.22s
9	5.24s	0.73s	5.85s	1.22s	4.03s
10	7.36s	0.79s	8.77s	1.32s	3.93s
12	10.78s	0.82s	14.15s	1.36s	3.73s
14	14.94s	0.79s	19.09s	1.37s	3.62s
16	14.15s	0.75s	23.67s	1.14s	3.52s

Each operation has a different optimal thread count:

Operation	Optimal BLAS	Behavior past optimum	Worst/Best ratio
Maxwell tSSS	7-8	Cliff at 9, catastrophic past 10	4.8x
Covariance (shrunk)	1	Monotonically worse with more threads	6.8x
ICA fit	9+	Gentle improvement, no cliff	1.7x
Inverse operator	8-16	Flat past 8	2.0x
apply_inverse_epochs	16	Linear scaling, no cliff	4.0x

The shrunk covariance is fastest at BLAS=1 and gets steadily worse — cross-validated eigendecompositions on small matrices have zero parallelism to exploit.

Per-step breakdown (coarse, with totals)

Step	BLAS=1	BLAS=4	BLAS=8	BLAS=16
Load raw + crop	0.25s	0.25s	0.25s	0.25s
Maxwell tSSS	4.27s	2.53s	2.65s	15.90s
Bandpass filter	0.34s	0.34s	0.35s	0.37s
ICA fit (20 comp)	2.02s	1.45s	1.40s	1.59s
ICA apply	0.30s	0.16s	0.14s	0.17s
Epoching + baseline	0.32s	0.20s	0.17s	0.23s
Covariance (shrunk)	2.00s	1.99s	2.21s	31.80s
Forward + inverse op	1.49s	0.99s	1.00s	1.40s
apply_inverse_epochs	13.42s	5.22s	3.92s	4.19s
TOTAL	25.26s	13.13s	12.02s	55.86s

Fix

# Add to .bashrc or cluster job script
export OPENBLAS_NUM_THREADS=8  # optimal for this pipeline on 16-core
# or OPENBLAS_NUM_THREADS=4   # safer for shared machines (within 9% of optimal)

This provides a 3.6x speedup for a typical MEG pipeline with zero code changes.

Why BLAS=8 and not fewer/more?

The pipeline contains two opposing forces:

• Decomposition-heavy steps (tSSS, covariance): optimal at 1-8 threads, catastrophic past 8
• Matmul-heavy steps (apply_inverse_epochs): scales linearly to 16

At BLAS=8, the decomposition steps are still near-optimal while the matmul step gets reasonable parallelism. Past 8, the decomposition steps cliff sharply and dominate the total.

[sharifhsn]

Cross-validation on 48-core machine (EPYC 9R14, no hyperthreading)

To test whether the cliff is caused by hyperthreading, I ran the same benchmark on a c7a.12xlarge: 48 physical cores, SMT disabled (Threads per core: 1).

System: AMD EPYC 9R14, 48 vCPU (= 48 physical cores), 1 socket, 6 CCDs, OpenBLAS 0.3.31.

Individual operations

BLAS	tSSS	ICA	CovShrunk	InvOp	InvEpochs
1	3.74s	1.21s	3.56s	1.38s	14.18s
12	3.09s	0.80s	3.34s	1.08s	3.99s
22	3.09s	0.79s	3.54s	1.04s	3.56s
23	3.02s	0.77s	3.54s	1.03s	3.53s
24	3.01s	0.77s	3.56s	0.99s	3.52s
25	4.62s	0.79s	3.71s	1.06s	3.54s
26	4.85s	0.82s	3.86s	1.29s	3.48s
28	5.10s	0.81s	3.87s	1.30s	3.46s
48	16.15s	0.81s	15.24s	1.65s	3.30s

Pipeline

BLAS= 1   24.4s  ########################
BLAS=12   12.4s  ############
BLAS=22   12.0s  ############
BLAS=23   12.0s  ############
BLAS=24   12.0s  ############         ← optimum
BLAS=25   13.7s  #############        ← +14% cliff
BLAS=26   14.6s  ##############

Hypothesis update: NOT hyperthreading — it's NUMA/CCD topology

The original hypothesis was that the cliff occurs at the physical core count (8 on the 16-vCPU machine). This is disproven — the 48-core machine has no hyperthreading, yet the cliff is at 24 (= N_cores / 2).

The EPYC 9R14 has 6 CCDs (Core Complex Dies) with 8 cores each. Each CCD has its own L3 cache. The cliff at 24 = 3 CCDs suggests OpenBLAS thread scheduling crosses a NUMA boundary when using more than half the cores, causing cache-thrashing for the memory-intensive SVD/QR workloads.

Comparison across machines

Machine	vCPUs	Physical cores	SMT	Cliff at	Cliff / cores
c7a.4xlarge (EPYC 7R13)	16	8	Yes (2x)	8→9	1.0x phys
c7a.12xlarge (EPYC 9R14)	48	48	No	24→25	0.5x phys

The consistent pattern: the cliff occurs when BLAS threads exceed the capacity of the "local" NUMA domain. On the 16-vCPU machine with 8 physical cores on a single CCD, 8 threads fill one CCD. On the 48-core machine with 6 CCDs, 24 threads fill 3 CCDs (half the die).

Updated recommendation

# Safe recommendation for any machine:
export OPENBLAS_NUM_THREADS=$(( $(nproc) / 2 ))
# or simply:
export OPENBLAS_NUM_THREADS=4  # conservative, works everywhere

The N/2 heuristic works on both tested machines and avoids the NUMA cliff while still getting good parallelism for matmul-heavy steps.

[sharifhsn]

Cross-architecture validation: AMD vs Intel vs no-SMT

Tested three machines to identify the root cause of the performance cliff.

Machine specs

Machine	CPU	Phys cores	SMT	NUMA nodes	vCPUs
c7a.4xlarge	AMD EPYC 7R13	8	Yes (2x)	1	16
c6i.4xlarge	Intel Xeon 8375C	8	Yes (2x)	1	16
c7a.12xlarge	AMD EPYC 9R14	48	No	1	48

Pipeline results (median of 3 runs)

BLAS	AMD 8c+HT	Intel 8c+HT	AMD 48c (no HT)
1	26.2s	28.2s	24.4s
4	13.2s	15.9s	16.0s
7	11.5s	14.5s	—
8	12.0s	14.8s	—
9	19.2s (+60%)	21.5s (+45%)	—
10	34.8s	27.4s	—
12	55.4s	36.2s	12.4s
16	56.4s	50.6s	—
24	—	—	12.0s
25	—	—	13.7s (+14%)
48	—	—	~47s (est.)

Per-operation comparison (AMD vs Intel, 16 vCPU)

BLAS		tSSS		CovShrunk		InvEpochs
	AMD	Intel	AMD	Intel	AMD	Intel
7	2.17s	3.07s	2.46s	3.92s	4.51s	5.25s
8	2.50s	3.38s	3.14s	5.13s	4.21s	4.95s
9	4.53s	5.25s	6.30s	6.99s	5.26s	6.04s

Conclusions

1. The cliff is at the physical core count on both AMD and Intel. On machines with 8 physical cores + hyperthreading, the cliff is consistently at 8→9 threads (AMD: +60%, Intel: +45%). This is vendor-independent.

2. It's NOT NUMA or CCD topology (on the 16-vCPU machines). Both the AMD and Intel instances have 1 socket, 1 NUMA node. The cliff is caused by hyperthreads sharing physical execution units, which hurts cache-intensive decomposition workloads (SVD, QR, eigendecomposition).

3. Intel degrades more gracefully past the cliff. At BLAS=12, AMD is at 55.4s (4.8x optimal) while Intel is at 36.2s (2.5x optimal). At BLAS=16, AMD=56.4s vs Intel=50.6s. Intel's shared LLC mesh may handle cross-hyperthread contention better than AMD's CCD architecture.

4. The no-SMT machine (48 core EPYC 9R14) also has a cliff at ~50% of cores, suggesting a secondary effect beyond hyperthreading — possibly L3 cache contention when too many cores compete for the same cache lines during decomposition.

Updated recommendation

# Robust recommendation for any machine:
# Use number of physical cores (not vCPUs)
export OPENBLAS_NUM_THREADS=$(lscpu | awk '/^Core\(s\) per socket:/{cores=$NF} /^Socket\(s\):/{sockets=$NF} END{print cores*sockets}')

# Or simply half of nproc (works for SMT and is conservative for non-SMT):
export OPENBLAS_NUM_THREADS=$(( $(nproc) / 2 ))

# Or just use a safe constant:
export OPENBLAS_NUM_THREADS=4

The physical-cores formula is optimal for SMT machines. The nproc/2 heuristic is safe everywhere. The constant 4 is conservative but never catastrophic.

[larsoner]

I am not sure we should ever set this value ourselves, but documenting this sort of thing could certainly be useful.

We could also consider a code change where the number of threads is queried using threadpoolctl and if it exceeds the number of physical CPU cores we could logger.info something in some of these functions that suggests setting the number of threads might lead to improved performance

[sharifhsn]

I agree on documentation and warning. The threadpoolctl infrastructure is also already there. I'm not a user so I wouldn't really know the best way to surface this to users. But I think that someone with better understanding of MNE should be able to implement quite easily.

We should also keep in mind that this is specific to Linux OpenBLAS. I believe that this is default for many scientists with HPCs, but none of this applies to MacBook M-series architecture which I have been testing performance on, which automatically accelerates BLAS operations regardless of threading. So any automated infrastructure would have to check for operating system, and to be fully robust, check for usage of OpenBLAS.

[JiwaniZakir]

The catastrophic regressions in the decomposition-heavy paths (Maxwell SSS/tSSS, shrunk covariance) are consistent with OpenBLAS's fork-join threading model incurring synchronization overhead that far outweighs the compute benefit when operating on the tall-skinny or small matrices typical of those routines — mne/preprocessing/maxwell.py calls several SVD/QR steps on matrices whose inner dimension is bounded by the number of SSS basis functions (~300–400), which is far below the threshold where 16-thread BLAS parallelism pays off. The empirical sweet spot of 4 threads for shrunk covariance and tSSS matches the known OpenBLAS crossover point for matrix sizes in the low hundreds. A targeted fix would be to use threadpoolctl.threadpool_limits(limits=1, user_api="blas") as a context manager around the innermost decomposition calls in mne/cov.py (_compute_covariance_auto) and the SVD steps in Maxwell filtering, rather than relying on a global thread count — this avoids penalizing users who do benefit from BLAS threading in other contexts like apply_inverse_epochs where the matrices are larger.

[larsoner]

A lot of this stuff is going to be OS, CPU, system load, and/or BLAS backend dependent. That makes it hard to come up with a general-purpose solution. We need to balance the potential user gains against potential user harms (e.g., setting the context ourselves could make things worse for some people sometimes!) and our own maintainability.

One easy, non-controversial thing to do here would be to document this stuff / how to control these things somewhere clearly, and users can tweak themselves as needed.

However, given the fact that Linux + x86_64 + OpenBLAS running stuff locally is going to be a very large contingent of our user base, maybe we can start by just trying to optimize a bit there. I wouldn't mind adding some context manager around some of these sensitive operations like:

def maxwell_filter(...):
    ...
    # more realistically, some private function that this function calls could have this in it:
    with _linalg_context(svd_size=(n_channels, n_channels)):
        ...

with our mne/utils/linalg.py having something like:

@functools.cache
def _is_linux_x86_64_openblas():
    return ...  # some conditional thing that can tell if this system is using this config


@contextlib.contextmanager
def _linalg_context(*, svd_size=None):
    context = nullcontext()  # default is to respect user config
    if _has_threadpoolctl:
        if _is_linux_x86_64_openblas():
            if svd_size is not None and np.prod(svd_size) >= ...:  # some heuristic that is reasonable
                context = threadpoolctl.threadpool_limits(limits=1, user_api="blas")
    with context:
        yield context

At a minimum we'd want to check the overhead of entering and exiting the threadpool_limits, but if we put it at the correct stack level (i.e., more toward the maxwell_filter level than each individual linalg.svd call) the overhead shouldn't really matter too much.