Decomposition (Estimators)

Decomposition (Estimators)#

The main scikit-learn compatible estimators for sparse hemodynamic deconvolution.

`SparseDeconvolution`(*, tr[, te, hrf_model, ...])	Sparse hemodynamic deconvolution using LARS or FISTA.
`LowRankPlusSparse`(*, tr[, te, hrf_model, ...])	Low-rank plus sparse deconvolution (SPLORA).
`StabilitySelection`(*, tr[, te, hrf_model, ...])	Stability selection for robust sparse deconvolution.

Overview#

pySPFM provides three main estimators that follow the scikit-learn API:

SparseDeconvolution: The core estimator for sparse hemodynamic deconvolution. Supports both univariate (voxel-wise) and multivariate (joint spatial) solving modes.
LowRankPlusSparse: Decomposes fMRI data into low-rank (structured noise/drift) and sparse (neural activity) components using robust PCA combined with sparse deconvolution.
StabilitySelection: Provides robust estimation through stability selection, solving the problem across multiple subsampled datasets and computing area-under-curve (AUC) stability scores.

Univariate vs Multivariate Mode#

The SparseDeconvolution estimator’s group parameter controls whether voxels are solved independently or jointly:

Mode	`group` value	Regularization	Parallelization	Criteria
Univariate	`0.0`	L1 (LASSO)	✅ via `n_jobs`	All
Multivariate	`> 0`	L1 + L2,1	❌ Joint solve	FISTA only

See the Usage page for detailed examples and mathematical background.

SparseDeconvolution#

class SparseDeconvolution(*, tr, te=None, hrf_model='spm', block_model=False, criterion='bic', debias=True, group=0.0, pcg=0.8, factor=1.0, lambda_echo=-1, max_iter=400, min_iter=50, tol=1e-06, n_jobs=1, positive=False, weights=None)[source]#

Sparse hemodynamic deconvolution using LARS or FISTA.

This estimator implements Sparse Paradigm Free Mapping (SPFM) for deconvolution of fMRI data. It estimates activity-inducing signals (spike model) or innovation signals (block model) from BOLD timeseries.

The estimator supports two solving modes:

Univariate mode (group=0.0, default):: Each voxel is solved independently using pure L1 (LASSO) regularization. This is efficient and suitable when spatial structure is not important.
Multivariate mode (group > 0.0):: All voxels are solved jointly using L2,1 mixed-norm regularization, which encourages spatial grouping of activity across voxels at the same timepoint. This leverages spatial structure in 4D fMRI data. Requires a FISTA-compatible criterion (‘mad’, ‘factor’, etc.).

Parameters:

tr (float) – Repetition time (TR) of the fMRI acquisition in seconds.
te (list of float, default=None) – Echo times in milliseconds for multi-echo data. If None or [0], assumes single-echo data.
hrf_model (str, default=’spm’) – HRF model to use. Options are ‘spm’, ‘glover’, or a path to a custom HRF file (.1D or .txt).
block_model (bool, default=False) – If True, estimate innovation signals (block model). If False, estimate activity-inducing signals (spike model).
criterion (str, default=’bic’) – Criterion for lambda selection:
- ‘bic’, ‘aic’: Information criteria (LARS solver, univariate only).
- ‘mad’, ‘mad_update’, ‘ut’, ‘lut’, ‘factor’, ‘pcg’, ‘eigval’: Noise-based criteria (FISTA solver, supports multivariate).
debias (bool, default=True) – If True, perform debiasing step to recover true amplitude.
group (float, default=0.0) – Weight for spatial grouping using L2,1 mixed-norm regularization. Range [0, 1] where:
- group=0.0: Pure L1/LASSO (univariate, voxel-wise independent).
- group=1.0: Pure L2,1 (multivariate, maximum spatial grouping).
- 0 < group < 1: Elastic net-like mix of L1 and L2,1.
When group > 0, the solver operates on all voxels jointly (multivariate mode) and requires a FISTA-compatible criterion.
pcg (float, default=0.8) – Percentage of maximum lambda to use (for criterion=’pcg’).
factor (float, default=1.0) – Factor to multiply noise estimate for lambda selection.
max_iter (int, default=400) – Maximum number of iterations for the FISTA solver.
tol (float, default=1e-6) – Convergence tolerance.
n_jobs (int, default=1) – Number of parallel jobs for the per-voxel LARS criteria (‘bic’, ‘aic’). The FISTA criteria solve all voxels in a single batched call, so n_jobs has no effect on them. In multivariate mode, computation is inherently joint.
positive (bool, default=False) – If True, enforce non-negative coefficients.
weights (ndarray of shape (n_voxels,), default=None) – Per-voxel weight map modulating the multivariate (L2,1 + L1 mixed-norm) penalty in multivariate mode (group > 0), adaptive-LASSO style: the effective threshold for voxel j is lambda / w_j, scaling both the sparsity and grouping terms. Higher weights penalize less (retain more activity); lower weights penalize more. Must be finite and strictly positive (w_j = 1 is neutral). Only valid in multivariate mode; passing weights with group == 0 (or a LARS criterion) raises ValueError.

Variables:

coef (ndarray of shape (n_timepoints, n_voxels)) – Estimated activity-inducing (spike) or innovation (block) signals.
lambda (ndarray of shape (n_voxels,) or float) – Selected regularization parameter. In multivariate mode, a single lambda is used for all voxels.
hrf_matrix (ndarray of shape (n_timepoints * n_echoes, n_timepoints)) – The HRF convolution matrix used for deconvolution.
n_features_in (int) – Number of voxels seen during fit.
n_iter (int) – Number of iterations run.

Examples

Univariate deconvolution (each voxel independent):

>>> from pySPFM import SparseDeconvolution
>>> import numpy as np
>>> X = np.random.randn(100, 50)  # 100 timepoints, 50 voxels
>>> model = SparseDeconvolution(tr=2.0, criterion='bic', group=0.0)
>>> model.fit(X)
SparseDeconvolution(criterion='bic', group=0.0, tr=2.0)

Multivariate deconvolution (joint spatial regularization):

>>> model = SparseDeconvolution(tr=2.0, criterion='factor', group=0.5)
>>> model.fit(X)  # All voxels solved jointly
SparseDeconvolution(criterion='factor', group=0.5, tr=2.0)

LowRankPlusSparse#

class LowRankPlusSparse(*, tr, te=None, hrf_model='spm', block_model=False, criterion='mad_update', eigval_threshold=0.1, debias=True, group=0.0, factor=1.0, max_iter=100, min_iter=10, tol=1e-06, n_jobs=1)[source]#

Low-rank plus sparse deconvolution (SPLORA).

This estimator implements the Sparse and Low-Rank Paradigm Free Mapping (SPLORA) algorithm, which separates fMRI data into low-rank (global/ structured) and sparse (transient neuronal) components.

Parameters:

tr (float) – Repetition time (TR) of the fMRI acquisition in seconds.
te (list of float, default=None) – Echo times in seconds for multi-echo data. If None or [0], assumes single-echo data.
hrf_model (str, default=’spm’) – HRF model to use. Options are ‘spm’, ‘glover’, or a path to a custom HRF file (.1D or .txt).
block_model (bool, default=False) – If True, estimate innovation signals (block model). If False, estimate activity-inducing signals (spike model).
criterion (str, default=’mad_update’) – Criterion for lambda selection. Options: ‘mad’, ‘mad_update’, ‘ut’, ‘lut’, ‘factor’, ‘pcg’, ‘eigval’.
eigval_threshold (float, default=0.1) – Eigenvalue threshold for low-rank estimation. Eigenvalues below this fraction of the maximum are set to zero.
debias (bool, default=True) – If True, perform debiasing step to recover true amplitude.
group (float, default=0.0) – Weight for spatial grouping (L2,1-norm). Range [0, 1].
max_iter (int, default=100) – Maximum number of FISTA iterations.
n_jobs (int, default=1) – Number of parallel jobs for voxel-wise computation.

Variables:

coef (ndarray of shape (n_timepoints, n_voxels)) – Estimated sparse (neuronal) activity-inducing signals.
low_rank (ndarray of shape (n_timepoints, n_voxels)) – Estimated low-rank (global/structured) component.
lambda (ndarray of shape (n_voxels,)) – Selected regularization parameter for each voxel.
hrf_matrix (ndarray of shape (n_timepoints * n_echoes, n_timepoints)) – The HRF convolution matrix used for deconvolution.
n_features_in (int) – Number of voxels seen during fit.

Examples

>>> from pySPFM import LowRankPlusSparse
>>> import numpy as np
>>> X = np.random.randn(100, 50)  # 100 timepoints, 50 voxels
>>> model = LowRankPlusSparse(tr=2.0)
>>> model.fit(X)
LowRankPlusSparse(tr=2.0)
>>> sparse_signal = model.coef_
>>> global_signal = model.low_rank_

StabilitySelection#

class StabilitySelection(*, tr, te=None, hrf_model='spm', block_model=False, n_surrogates=50, n_lambdas=None, threshold=0.6, n_jobs=1)[source]#

Stability selection for robust sparse deconvolution.

This estimator uses stability selection to identify robust features in the sparse deconvolution problem. It runs multiple deconvolutions on subsampled data and returns the selection frequency (AUC) for each timepoint.

Parameters:

tr (float) – Repetition time (TR) of the fMRI acquisition in seconds.
te (list of float, default=None) – Echo times in seconds for multi-echo data.
hrf_model (str, default=’spm’) – HRF model to use.
block_model (bool, default=False) – If True, estimate innovation signals (block model).
n_surrogates (int, default=50) – Number of bootstrap surrogates.
n_lambdas (int, default=None) – Number of lambda values in the regularization path. If None, uses n_scans.
threshold (float, default=0.6) – Selection threshold for considering a feature as selected.
n_jobs (int, default=1) – Number of parallel jobs.

Variables:

selection_frequency (ndarray of shape (n_timepoints, n_voxels)) – Selection frequency (AUC) for each timepoint and voxel.
coef (ndarray of shape (n_timepoints, n_voxels)) – Binary selection indicators based on threshold.
hrf_matrix (ndarray) – The HRF convolution matrix.
n_features_in (int) – Number of voxels.

Examples

>>> from pySPFM import StabilitySelection
>>> import numpy as np
>>> X = np.random.randn(100, 50)
>>> model = StabilitySelection(tr=2.0, n_surrogates=20)
>>> model.fit(X)
StabilitySelection(n_surrogates=20, tr=2.0)
>>> selection_freq = model.selection_frequency_

Decomposition (Estimators)

Contents

Decomposition (Estimators)#

Overview#

Univariate vs Multivariate Mode#

SparseDeconvolution#

LowRankPlusSparse#

StabilitySelection#