Lead-field-constrained Block-Term Decomposition for evoked EEG source imaging

A self-contained, reproducible pipeline for the method in the manuscript "Jointly Localized Spatiotemporal Components of Evoked EEG via a Lead-Field-Constrained Block-Term Decomposition" (A. Faes). It fits a rank-(Lr, Lr, 1) Block-Term Decomposition whose spatial factor is constrained to the lead-field column space, so decomposing multi-trial evoked EEG returns, per block, an anatomically grounded source map + time course + trial loading, estimated jointly. Benchmarked against MNE, sLORETA, eLORETA, an LCMV beamformer, the sparse-Bayesian Champagne solver, and the decompose-then-localize CP+eLORETA tensor baseline.

Every number in the manuscript is regenerated by these scripts and archived as a CSV under results/. The table near the end maps each script to its CSV to the manuscript claim it backs.

What the contribution is (and is not)

The contribution is a novel capability, not a leaderboard win. Folding the lead field INTO the decomposition yields components that are cortical by construction - a different object from what a per-sample inverse returns (a map per time sample) and from decompose-then-localize (abstract scalp components localized afterward, a separate ill-posed step). The method is characterized honestly:

• It is the most focal of the classical distributed solvers under the realistic-noise model on the 64-channel SPHERE (dispersion 13.4 mm vs MNE 14.5, eLORETA 17.4, sLORETA 21.8), but the margin over eLORETA is not significant (Wilcoxon p=0.055) and reverses on denser/realistic forwards (BEM: eLORETA 22.1 < BTD 26.8; New York head: BTD 23-33 > eLORETA 17-22).
• It is not the sharpest point localizer: Champagne reaches near-delta maps on idealized point sources, though under a matched sparse prior / inverse-crime configuration.
• It has no multi-source-recovery advantage: under a fair, scale-invariant top-K metric (applied identically to all methods) it is the worst recoverer, not the best. The pipeline reports this metric; no recovery claim is made.
• It is correlation-robust where an LCMV beamformer collapses - but so are the linear distributed inverses, so this is a shared property, not a unique win.
• It is phase-locked only: the rank-one trial factor cannot represent a per-trial-varying waveform, so time-course correlation falls to 0.16 on non-phase-locked (MVAR) sources versus 0.85 for a per-trial eLORETA.

The model

Evoked / multi-trial EEG is a 3-way tensor M of shape (C channels x T time x K trials). We fit M ~ sum_r [(L Gr) Br^T] o cr, where o is the outer product on the trial mode and:

symbol	shape	meaning
L	C x N	lead field (fixed/radial orientation)
Gr	N x Lr	row-sparse source distribution of block r - the source map
Br	T x Lr	temporal subspace (time course) of block r
cr	K	trial loading of block r

The spatial factor is constrained to L Gr, so each block's spatial signature IS a source map; row-sparsity (group-L2 over the Lr columns of each source) gives a focal map. Br and cr act as a consistency prior - only patterns reproducible across time AND trials earn a component. Estimation is greedy/deflationary like BTTR (Faes, Camarrone & Van Hulle, IEEE TNNLS 2022): extract the dominant block, deflate, repeat, then a joint block-coordinate refit. Within a block, constrained ALS alternates cr (mode-3 least squares), Br (least squares), Gr (ISTA with row soft-threshold). The default per-block rank is block_rank=2 (Lr=2); it improves recovered-time-course fidelity over Lr=1 (0.88 to 0.93 on the sphere) at no cost in focality.

Files

core/ is the library, experiments/ the drivers:

file	role
core/leadfield.py	lead field: analytic MNE sphere, fsaverage 3-layer BEM (surface or volume)
core/nyhead.py	realistic New York head lead field (171 elec x 2004 dipoles; Huang 2016)
core/simulate.py	simple ground-truth simulation (white sensor noise) to a tensor M (C,T,K)
core/paper_sim.py	realistic brain-noise simulation (500 noise dipoles; pink/coherent-pink/coherent-AR, Frobenius mixing)
core/realdata.py	real evoked-EEG tensor + matching BEM forward (MNE sample)
core/btd_inverse.py	the solver (LeadfieldBTD)
core/baselines.py	MNE / sLORETA / eLORETA / LCMV / Champagne / CP+eLORETA (numpy)
core/metrics.py	DLE, peak error, detection AUC, spatial dispersion, sparsity, fair top-K recovery, signal correlation
core/paths.py	path anchors (RESULTS_DIR/FIG_DIR -> results/, DATA_DIR -> data/)
experiments/run_pipeline.py	simulation orchestration -> metrics_main.csv, metrics_correlated.csv, metrics_recovery.csv, figures
experiments/run_paper_sim.py	realistic-noise fidelity/focality table -> metrics_paper_sim.csv, stats_paper_sim.csv
experiments/run_nonphaselocked.py	phase-locked-only boundary test -> metrics_nonphaselocked.csv
experiments/run_nyhead.py	New York head depth/proximity focality comparison (prints)
experiments/run_real.py	real-data orchestration (retinotopic contralateral check) -> results/real/

This is a Python package: run drivers as modules from this directory. Outputs are anchored by core/paths.py to results/ (+ results/figs/, results/bem/, results/real/); the committed input is data/ny_leadfield.npz.

Install and run

pip install -r requirements.txt   # numpy scipy scikit-learn tensorly matplotlib mne nibabel
python paper.py                   # rebuild every result + figure, in dependency order

Or run a single driver:

python -m experiments.run_pipeline --forward sphere     # main table + correlated + fair recovery + figures
python -m experiments.run_pipeline --forward bem        # realistic fsaverage 3-layer BEM
python -m experiments.run_paper_sim --forward sphere --seeds 8   # the headline numbers + Wilcoxon
python -m experiments.run_nonphaselocked --seeds 8      # the phase-locked-only limit (0.16 vs 0.85)
python -m experiments.run_nyhead                        # New York head setup-dependence (from the committed .npz)
python -m experiments.run_real                          # real MNE sample EEG: retinotopic contralateral check

External data: the BEM/real-data paths download the open fsaverage (~0.4 GB) and sample (~1.6 GB) datasets via MNE, cached under ~/mne_data (not in this repo). The New York head raw model (~711 MB) is NOT shipped; nyhead.py uses the committed extracted data/ny_leadfield.npz (2.7 MB), so run_nyhead runs out of the box.

Results (realistic brain-noise model, sphere, 64 ch, 8 seeds, Lr=2)

From results/metrics_paper_sim.csv (incoherent-pink rows). sig_corr = mean |corr| of the reconstructed source time course at the true location; disp = energy-weighted spatial spread (mm, lower = more focal); dle = dipole localization error (mm).

method	sig_corr	dispersion (mm)	DLE (mm)
BTD-B	0.93	13.4	29.6
Champagne	1.00	0.1	1.7
MNE	0.98	14.5	31.5
sLORETA	0.98	21.8	28.3
eLORETA	0.99	17.4	27.9
LCMV	0.96	18.4	6.8

Among the classical distributed solvers BTD-B is the most focal, but the margin over eLORETA is not significant (stats_paper_sim.csv: p=0.055; vs sLORETA p=0.008). Champagne and LCMV are sharper point localizers (low DLE); Champagne under a matched sparse prior (inverse-crime caveat). On a correlated-source stress test (metrics_correlated.csv) the LCMV beamformer collapses (dispersion 17 to 53 mm, AUC 1.00 to 0.73) while BTD-B stays focal (17 to 18 mm, AUC 0.93 to 0.90); the linear inverses are correlation-robust too. Denser/realistic forwards reverse the focality edge (BEM eLORETA 22.1 < BTD 26.8; New York head BTD 23-33 > eLORETA 17-22).

Real evoked EEG (run_real.py)

realdata.py builds the same M (C,T,K) tensor from real single-trial epochs of the open MNE sample visual paradigm, with the subject's own BEM forward. With no ground-truth dipole the check is anatomical: each method's peak source is mapped to its Desikan-Killiany label. BTD-B passes the retinotopic contralateral check - a right-visual-field stimulus localizes to left lateral-occipital cortex (and left-field to right). Auditory localization is intrinsically harder in scalp EEG (weak N1) and is not claimed.

Reproducibility map: script -> CSV -> manuscript claim

manuscript claim	script	committed CSV
time-course corr 0.93 +/- 0.09; spread 13.4 vs MNE 14.5 / eLOR 17.4 / sLOR 21.8	run_paper_sim.py	results/metrics_paper_sim.csv
Wilcoxon p=0.055 (BTD vs eLORETA)	run_paper_sim.py	results/stats_paper_sim.csv
Champagne ~0.1 mm, ~1.00 corr	run_paper_sim.py / run_pipeline.py	metrics_paper_sim.csv / metrics_main.csv
LCMV collapse 17 to 53 mm, AUC 1.00 to 0.73; BTD stays focal	run_pipeline.py	results/metrics_correlated.csv
no recovery advantage (fair top-K)	run_pipeline.py	results/metrics_recovery.csv
BEM: eLORETA more focal (22.1 < 26.8)	run_paper_sim.py --forward bem	results/bem/metrics_paper_sim.csv
New York head: BTD 23-33 > eLORETA 17-22	run_nyhead.py	(printed; lead field in data/ny_leadfield.npz)
phase-locked only: 0.16 vs 0.85	run_nonphaselocked.py	results/metrics_nonphaselocked.csv
retinotopic contralateral check passes	run_real.py	results/real/peak_labels.csv

Honest limitations

• The deflationary fit is a greedy approximation to the joint BTD; a joint ALS refit is the natural next step.
• The block count is set by hand (benchmarks use R=3); LeadfieldBTD(L, n_blocks="auto") selects it from a variance-explained elbow, and run_real.py reports the auto R*. Per-block rank Lr is fixed at 2 here; the ACE/BIC machinery in the BTTR repo is the intended selector.
• The BEM uses a template (fsaverage) head with a standard montage; a subject-specific MRI + digitized electrodes would sharpen real-data localization.
• Deep sources degrade for all methods.

License: MIT (LICENSE).