HarmonICA

Code companion to Quasi-Linear ICA for Motor Unit Decomposition during Dynamic Contractions (anonymous NeurIPS 2026 submission).

HarmonICA recovers motor-neuron spike trains from non-stationary surface EMG by combining a static linear separator with a learned, low-rank, time-varying invertible compensation. The two are trained against decoupled losses — a power-law independence contrast on the uncompensated projection (separator only) and a two-component GMM stationarity contrast on the compensated source (compensator only) — so the source-extraction step inherits the identifiability of classical linear ICA while the compensator absorbs all non-stationary distortion. The closed-form Woodbury inverse of the compensator enables per-spike template subtraction during sequential peel-off.

Install

pip install -e .

Requires Python ≥ 3.10 and a working PyTorch install (CUDA recommended; the bundled examples also run on CPU). Optional extras:

pip install -e .[adapt]   # h5py, required by the vendored adapt_decomp baseline
pip install -e .[plot]    # matplotlib, only needed if the compensator's live-plot helper is invoked
pip install -e .[dev]     # pytest
pip install -e .[all]     # everything

Quickstart

python examples/01_quickstart.py

Loads the bundled synthetic recording (examples/data/synthetic_demo.npz, 32 channels × 5 s @ 2048 Hz, 5 motor units, sinusoidal angle profile), runs HarmonICA end-to-end, and prints per-source F1 against ground truth. Expected mean F1 ≥ 0.99, runtime ~60 s on a single A40 GPU (~3 min on CPU).

For a comparison against three adaptive baselines on the same data:

python examples/02_compare_baselines.py

See examples/README.md for the generative model behind the demo and how to regenerate it via make_synthetic.py.

Programmatic use

import numpy as np
from nonstationary_cbss.baselines import run_signal_space
from nonstationary_cbss.evaluate import evaluate_decomposition

# emg: (n_channels, n_samples) float32, fs in Hz, gt: {mu_id: spike_indices}
spikes = run_signal_space(
    emg, fs,
    epochs=200,                # alternating-optimisation steps per source
    sep_only_warmup=50,        # delay compensation by 50 separator-only steps
    R=4,                       # Toeplitz extension factor
    d_eff=128,                 # PCA target dim post-extension
    ica_seed=1909, seed=42,
)
ev = evaluate_decomposition(spikes, gt, fs)
matched = sum(1 for ps in ev["per_source"] if ps["gt_id"] >= 0)
print(f"matched motor units: {matched}")

run_signal_space is the paper's headline method (Algorithm 1): warm-started rank-3 Woodbury compensator with decoupled GMM/independence losses and neural peel-off. The other entry points in baselines/ implement the adaptive baselines compared in the paper: Adapt (Mendez Guerra 2024), Chen 2020, Glaser 2018 (segmented + the paper-faithful cyclostationary variant), Kramberger 2021, and Yeung 2024.

Layout

src/nonstationary_cbss/
├── evaluate.py                       # Spike matching, F1 / precision / recall, shift alignment
├── ica/                              # ICA foundation: bandpass, Toeplitz extension, ZCA, FastICA loop
├── baselines/
│   ├── signal_space.py               # ★ HarmonICA (this paper)
│   ├── adapt.py                      # Mendez Guerra 2024 (vendored — see vendored_adapt_decomp/)
│   ├── chen.py                       # Chen et al. 2020 adaptive CKC
│   ├── glaser.py                     # Glaser & Holobar 2018 segmented cBSS (simplified)
│   ├── glaser_cyclostationary.py     # Glaser & Holobar 2018/2019 paper-faithful Pseudocode 1
│   ├── kramberger.py                 # Kramberger & Holobar 2021 filter-angle Kalman
│   ├── yeung.py                      # Yeung et al. 2024 directional-forgetting RLS
│   └── vendored_adapt_decomp/        # Mendez Guerra 2024 reference implementation (MIT)
└── experimental/_compensator/        # Internal compensator package: rank-r Woodbury layer,
                                      #   joint-step alternating optimisation, peel-off (MIT).
                                      #   Required by run_signal_space.

The _compensator/ package provides the ConditionalLowRank Woodbury layer, the State container, and the DecompositionModule.joint_step alternating-optimisation loop that underlies run_signal_space.

Reproducing the paper

The paper evaluates on muniverse-benchmark-dynamic, a public benchmark of dynamic high-density EMG with ground-truth spike trains. Construction details are in the paper's supplementary §B.

• Dataset: https://huggingface.co/datasets/emgbenchanon/muniverse-dynamic-v2
• Hardware: NVIDIA A40, single GPU. ~44–47 s per recording for HarmonICA; baseline runtimes range 3–42 s per recording (paper supplementary §D).
• Seed mapping: 5-seed test runs use ica_seed = 1909 + idx and nn_seed = 42 + idx for idx ∈ {0..4}. Kramberger is deterministic given calibration data.
• Hyperparameters: paper-headline values are the defaults in run_signal_space (epochs=200, sep_only_warmup=50, R=4, d_eff=128, compensation_rank=3). Hparam sweeps (rank, extension factor, loss family, peel-off variant) are described in the paper's supplementary §E.

After downloading the benchmark, the per-recording API is the same as in examples/01_quickstart.py:

spikes = run_signal_space(emg, fs, epochs=200, sep_only_warmup=50, R=4, d_eff=128)

Tests

pip install -e .[dev]
pytest tests/                       # smoke test: imports + small end-to-end run
# or:
python tests/test_smoke.py

The smoke test runs in ~1 minute on CPU.

License

MIT — see LICENSE. Vendored components retain their original licenses (also MIT).