This example shows how to use non-negative matrix factorization (NMF) to decompose a fUSI recording into non-negative spatial maps and their associated non-negative time courses. It complements the PCA and FastICA examples in the same gallery.

NMF is unique among the decomposers in ConfUSIus because it requires strictly non-negative inputs¹. Power Doppler fUSI signals are non-negative by construction, so they pass the constraint directly. However, raw power is dominated by each voxel's baseline intensity, which can make bright vessels dominate the factorization.

A practical workaround is to center and scale each voxel across time, then split the standardized signal into separate positive and negative channels, preserving NMF's non-negativity constraint. NMF can then discover additive components in above-baseline and below-baseline fluctuations separately.

Load a fUSI recording¶

We use the same spontaneous activity recording from the Nunez-Elizalde 2022 dataset as in the PCA and FastICA examples. See the Datasets user guide for more details on how to download this dataset using ConfUSIus.

from pathlib import Path

import matplotlib as mpl
import matplotlib.pyplot as plt
import xarray as xr

import confusius as cf

# Adapt background color to the current Matplotlib style.
bg_color = mpl.colors.to_hex(mpl.rcParams["figure.facecolor"])

# Keep notebook output compact for large DataArray displays.
xr.set_options(display_expand_data=False)

bids_root = cf.datasets.fetch_nunez_elizalde_2022(
    subjects="CR022",
    sessions="20201011",
    tasks="spontaneous",
    acqs="slice03",
)

pwd_path = (
    Path(bids_root)
    / "sub-CR022"
    / "ses-20201011"
    / "fusi"
    / "sub-CR022_ses-20201011_task-spontaneous_acq-slice03_pwd.nii.gz"
)
data = cf.load(pwd_path).compute()
data

Correct for brain motion¶

As in the PCA and FastICA examples, we first perform a rigid transformation correction with register_volumewise to mitigate brain motion.

data = cf.registration.register_volumewise(data, learning_rate=1e-2)

/home/runner/work/confusius/confusius/.venv/lib/python3.14/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Standardize for NMF¶

NMF requires non-negative inputs, but raw Power Doppler is dominated by each voxel's baseline intensity rather than the temporal fluctuations we want to group. We therefore:

Z-score each voxel's time series with standardize to remove its mean and put voxels on a comparable scale.
Split the standardized signal into separate positive and negative parts.

This keeps the sign information— above-baseline versus below-baseline fluctuations—while still presenting a non-negative matrix to NMF.

z = cf.signal.standardize(data)
data_nmf = xr.concat(
    [z.clip(min=0), (-z).clip(min=0)],
    dim=xr.IndexVariable("sign", ["pos", "neg"]),
)

Fit temporal NMF¶

NMF wraps the familiar scikit-learn NMF estimator while preserving fUSI DataArray metadata and coordinates. With mode="temporal" (the default), it fits on (time, voxels) and returns:

maps_: non-negative spatial maps. Because we split the input into positive and negative channels, the maps here have shape (component, sign, z, y, x).
fit_transform: non-negative time courses of shape (time, component).

nmf = cf.decomposition.NMF(n_components=10, random_state=0, max_iter=500)
signals = nmf.fit_transform(data_nmf)
signals

/home/runner/work/confusius/confusius/.venv/lib/python3.14/site-packages/sklearn/decomposition/_nmf.py:1720: ConvergenceWarning: Maximum number of iterations 500 reached. Increase it to improve convergence.
  warnings.warn(

xarray.DataArray

time: 751
component: 10

0.1658 0.01306 0.0 0.1819 0.0007141 ... 0.04795 0.0 0.3565 0.0 0.06798

array([[0.1657736 , 0.01305686, 0.        , ..., 0.17466769, 0.55248916,
        0.        ],
       [0.07842413, 0.        , 0.        , ..., 0.0567409 , 0.521537  ,
        0.        ],
       [0.02165581, 0.        , 0.02380166, ..., 0.01371775, 0.36238965,
        0.        ],
       ...,
       [0.10035805, 0.01997683, 0.        , ..., 0.09034111, 0.25914687,
        0.        ],
       [0.20754123, 0.0196385 , 0.        , ..., 0.11031527, 0.12536858,
        0.        ],
       [0.18679054, 0.01474215, 0.        , ..., 0.35648656, 0.        ,
        0.06798259]], shape=(751, 10), dtype=float32)

Coordinates: (2)
- time
  (time)
  float64
  10.61 10.91 11.21 ... 235.4 235.7
  units :
  s
  volume_acquisition_reference :
  start
  volume_acquisition_duration :
  0.3
```
array([ 10.608,  10.908,  11.208, ..., 235.095, 235.395, 235.695], shape=(751,))
```
- component
  (component)
  int64
  0 1 2 3 4 5 6 7 8 9
```
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
```
Attributes: (1)
long_name :
NMF signals

Reconstruction error¶

reconstruction_err_ is the Frobenius norm of X - WH, where W are the spatial maps and H the time courses. It gives a sense of how well the chosen number of components explains the standardized data. A quantitative model-selection procedure is out of scope here, but the trace is useful when sweeping n_components and looking for diminishing returns.

print(f"reconstruction_err_: {nmf.reconstruction_err_:.3f}")
print(f"n_iter_: {nmf.n_iter_}")

reconstruction_err_: 1911.476
n_iter_: 500

Spatial maps and time courses¶

Looking at the spatial maps and the associated time courses side by side is a useful first sanity check. Here each component has two map panels: one for above-baseline fluctuations (pos) and one for below-baseline fluctuations (neg). Localized, anatomically plausible structure paired with a clear transient in the time course tends to reflect a coherent spatiotemporal pattern, while diffuse maps paired with noisy or drift-like fluctuations often indicate residual motion or physiological artefacts.

n_show = 10
fig = plt.figure(figsize=(11.5, 12.0), constrained_layout=True)
fig.patch.set_facecolor(bg_color)
gs = fig.add_gridspec(n_show, 2, width_ratios=[1.4, 3])

axes_tc = [fig.add_subplot(gs[i, 1]) for i in range(n_show)]
for ax in axes_tc[1:]:
    ax.sharex(axes_tc[0])

for i, comp in enumerate(range(n_show)):
    component_map = nmf.maps_.isel(component=comp, drop=True)
    vmax = float(component_map.max())
    map_gs = gs[i, 0].subgridspec(1, 2, wspace=0.02)

    for j, sign in enumerate(["pos", "neg"]):
        cf.plotting.plot_volume(
            component_map.sel(sign=sign, drop=True),
            axes=fig.add_subplot(map_gs[0, j]),
            slice_mode="z",
            cmap="viridis",
            vmin=0,
            vmax=vmax,
            show_axes=False,
            show_colorbar=False,
            show_titles=False,
            bg_color=bg_color,
        )

    signals.sel(component=comp).plot(ax=axes_tc[i], lw=1.1)
    axes_tc[i].set_title(f"Component {comp + 1}")
    axes_tc[i].set_ylabel("Signal")
    axes_tc[i].set_xlabel("")

for ax in axes_tc[:-1]:
    ax.tick_params(labelbottom=False)
axes_tc[-1].set_xlabel("Time (s)")
_ = fig.suptitle(
    "Temporal NMF: positive/negative maps and time courses (first 10 components)",
    fontsize=21,
)

Spatial NMF¶

NMF also accepts mode="spatial", which transposes the data to (voxels, time) before fitting. The output convention is identical to temporal mode — maps_ still holds the non-negative spatial maps (here with pos/neg channels) and fit_transform returns their non-negative time courses — so the choice between the two modes mirrors the temporal/spatial choice offered by PCA and FastICA.

nmf_spatial = cf.decomposition.NMF(
    n_components=10, mode="spatial", random_state=0, max_iter=500
)
signals_s = nmf_spatial.fit_transform(data_nmf)
signals_s

/home/runner/work/confusius/confusius/.venv/lib/python3.14/site-packages/sklearn/decomposition/_nmf.py:1720: ConvergenceWarning: Maximum number of iterations 500 reached. Increase it to improve convergence.
  warnings.warn(

xarray.DataArray

time: 751
component: 10

-49.06 -108.7 23.6 477.0 273.2 ... 31.01 -59.16 840.9 351.1 -65.79

array([[ -49.060493, -108.72904 ,   23.599054, ..., 1360.9792  ,
         301.0807  ,  -79.03104 ],
       [ -46.849873, -102.89783 ,   88.44416 , ..., 1054.8505  ,
         202.54056 , -314.424   ],
       [ -48.20063 ,  -89.69527 ,   50.060436, ...,  992.90845 ,
          60.874683, -348.29483 ],
       ...,
       [ -52.879208, -147.47987 ,    5.829556, ..., 1231.0189  ,
         251.24892 , -261.63306 ],
       [ -48.73986 , -141.4787  ,   35.300552, ..., 1235.6964  ,
         480.693   , -212.94316 ],
       [ -41.16797 , -107.332   ,  -33.120617, ...,  840.94934 ,
         351.05264 ,  -65.78985 ]], shape=(751, 10), dtype=float32)

Coordinates: (2)
- time
  (time)
  float64
  10.61 10.91 11.21 ... 235.4 235.7
  units :
  s
  volume_acquisition_reference :
  start
  volume_acquisition_duration :
  0.3
```
array([ 10.608,  10.908,  11.208, ..., 235.095, 235.395, 235.695], shape=(751,))
```
- component
  (component)
  int64
  0 1 2 3 4 5 6 7 8 9
```
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
```
Attributes: (1)
long_name :
NMF signals

Spatial maps and time courses¶

As in temporal mode, we inspect the positive and negative spatial maps and their corresponding time courses side by side.

n_show = 10
fig = plt.figure(figsize=(11.5, 12.0), constrained_layout=True)
fig.patch.set_facecolor(bg_color)
gs = fig.add_gridspec(n_show, 2, width_ratios=[1.4, 3])

axes_tc = [fig.add_subplot(gs[i, 1]) for i in range(n_show)]
for ax in axes_tc[1:]:
    ax.sharex(axes_tc[0])

for i, comp in enumerate(range(n_show)):
    component_map = nmf_spatial.maps_.isel(component=comp, drop=True)
    vmax = float(component_map.max())
    map_gs = gs[i, 0].subgridspec(1, 2, wspace=0.02)

    for j, sign in enumerate(["pos", "neg"]):
        cf.plotting.plot_volume(
            component_map.sel(sign=sign, drop=True),
            axes=fig.add_subplot(map_gs[0, j]),
            slice_mode="z",
            cmap="viridis",
            vmin=0,
            vmax=vmax,
            show_axes=False,
            show_colorbar=False,
            show_titles=False,
            bg_color=bg_color,
        )

    signals_s.sel(component=comp).plot(ax=axes_tc[i], lw=1.1)
    axes_tc[i].set_title(f"Component {comp + 1}")
    axes_tc[i].set_ylabel("Signal")
    axes_tc[i].set_xlabel("")

for ax in axes_tc[:-1]:
    ax.tick_params(labelbottom=False)
axes_tc[-1].set_xlabel("Time (s)")
_ = fig.suptitle(
    "Spatial NMF: positive/negative maps and time courses (first 10 components)",
    fontsize=21,
)

Total running time: 290.8 s

Launch in Binder Download .py Download .ipynb

Lee, D. D., and Seung, H. S. (1999). "Learning the parts of objects by non-negative matrix factorization". Nature, 401(6755), 788-791. ↩