Blind Source Separation and ICA Equivalence.#

This example demonstrates how Nonlinear DSS can perform Blind Source Separation (BSS), effectively recovering independent sources from mixed signals. It explicitly shows the equivalence between DSS with specific nonlinearities and ICA.

It covers synthetic blind source separation, the link between DSS nonlinearities and FastICA, and a real MEG decomposition that exposes both artifact and brain-like components.

Authors: Sina Esmaeili (sina.esmaeili@umontreal.ca): Hamza Abdelhedi (hamza.abdelhedi@umontreal.ca)

import matplotlib.pyplot as plt
import mne
import numpy as np
from mne.datasets import sample
from scipy import stats

from mne_denoise.dss import IterativeDSS, KurtosisDenoiser, TanhMaskDenoiser, beta_tanh
from mne_denoise.viz import (
    plot_component_summary,
    plot_component_time_series,
    plot_signal_overlay,
)

print(__doc__)

Run DSS with Tanh Nonlinearity (Robust ICA)#

The TanhMaskDenoiser implements the tanh nonlinearity, which is robust to outliers.

This section compares the Newton-style update (beta=beta_tanh) with plain gradient ascent (beta=None). The faster convergence is the same idea that makes FastICA efficient in practice.

print("\nRunning DSS with Tanh Nonlinearity (Robust)...")

# 1. Gradient Ascent (Slow)
print("  Fitting with Gradient Ascent (beta=None)...")
dss_grad = IterativeDSS(
    denoiser=TanhMaskDenoiser(),
    method="deflation",
    n_components=n_sources,
    beta=None,  # Gradient ascent
    random_state=42,
    verbose=False,
)
dss_grad.fit(X)

# 2. Newton Method (Fast - FastICA style)
print("  Fitting with Newton Method (beta=beta_tanh)...")
dss_tanh = IterativeDSS(
    denoiser=TanhMaskDenoiser(),
    method="deflation",
    n_components=n_sources,
    beta=beta_tanh,  # Newton step
    random_state=42,
    verbose=False,
)
dss_tanh.fit(X)
S_dss_tanh = dss_tanh.transform(X)

# Compare iterations
iters_grad = dss_grad.convergence_info_[:, 0].sum()
iters_newton = dss_tanh.convergence_info_[:, 0].sum()
print(f"  Gradient Iterations: {iters_grad:.0f}")
print(
    f"  Newton Iterations:   {iters_newton:.0f} "
    f"(Speedup: {iters_grad / iters_newton:.1f}x)"
)

Running DSS with Tanh Nonlinearity (Robust)...
  Fitting with Gradient Ascent (beta=None)...
  Fitting with Newton Method (beta=beta_tanh)...
  Gradient Iterations: 116
  Newton Iterations:   20 (Speedup: 5.8x)

Run DSS with Kurtosis Nonlinearity (Standard FastICA)#

KurtosisDenoiser with nonlinearity=’cube’ maximizes kurtosis ($s^3$), which is the classic definition of FastICA.

print("Running DSS with Kurtosis Nonlinearity (FastICA standard)...")
dss_kurt = IterativeDSS(
    denoiser=KurtosisDenoiser(nonlinearity="cube"),
    method="deflation",
    n_components=n_sources,
    beta=-3.0,  # Newton step for kurtosis
    random_state=42,
    verbose=False,
)
dss_kurt.fit(X)
S_dss_kurt = dss_kurt.transform(X)

Running DSS with Kurtosis Nonlinearity (FastICA standard)...

Comparison with sklearn FastICA#

We run sklearn.decomposition.FastICA to serve as a ground truth benchmark.

from sklearn.decomposition import FastICA

print("Running sklearn FastICA (Benchmark)...")
ica = FastICA(
    n_components=n_sources, algorithm="deflation", fun="logcosh", random_state=42
)
S_fastica = ica.fit_transform(X.T).T

Running sklearn FastICA (Benchmark)...

Evaluate Performance (Correlation with True Sources)#

We compute the absolute correlation matrix between recovered components and true sources. A perfect recovery would have one 1.0 per row/column (permutation matrix).

def match_sources(S_est, S_true):
    """Calculate best correlation match for each source."""
    n_est = S_est.shape[0]
    n_true = S_true.shape[0]
    corr = np.zeros((n_est, n_true))
    for i in range(n_est):
        for j in range(n_true):
            corr[i, j] = np.abs(np.corrcoef(S_est[i], S_true[j])[0, 1])
    return corr


print("\n--- Evaluation ---")

corr_tanh = match_sources(S_dss_tanh, S_true)
corr_kurt = match_sources(S_dss_kurt, S_true)

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

im0 = axes[0].imshow(corr_tanh, vmin=0, vmax=1, cmap="Greens")
axes[0].set_title(
    f"DSS (Tanh) Match\nMean Max Corr: {np.mean(np.max(corr_tanh, axis=1)):.3f}"
)
axes[0].set_ylabel("Recovered")
axes[0].set_xlabel("True")
plt.colorbar(im0, ax=axes[0])

im1 = axes[1].imshow(corr_kurt, vmin=0, vmax=1, cmap="Blues")
axes[1].set_title(
    f"DSS (Kurtosis) Match\nMean Max Corr: {np.mean(np.max(corr_kurt, axis=1)):.3f}"
)
axes[1].set_xlabel("True")
plt.colorbar(im1, ax=axes[1])

corr_ica = match_sources(S_fastica, S_true)
im2 = axes[2].imshow(corr_ica, vmin=0, vmax=1, cmap="Oranges")
axes[2].set_title(
    f"sklearn FastICA Match\nMean Max Corr: {np.mean(np.max(corr_ica, axis=1)):.3f}"
)
axes[2].set_xlabel("True")
plt.colorbar(im2, ax=axes[2])

plt.suptitle("Source Recovery Quality (Abs Correlation)")
plt.tight_layout()
plt.show()

# Plot recovered time series using viz module
print("Visualizing Recovered Sources (Stacked)...")
# We treat the sources as "components" of the estimator

Source Recovery Quality (Abs Correlation), DSS (Tanh) Match Mean Max Corr: 0.961, DSS (Kurtosis) Match Mean Max Corr: 0.936, sklearn FastICA Match Mean Max Corr: 0.999

--- Evaluation ---
Visualizing Recovered Sources (Stacked)...

Recovered Source Time Series#

plot_component_time_series(dss_tanh, data=X, show=False)
plt.gcf().suptitle("Recovered Sources (DSS Tanh) - Newton Optimization")
plt.show()

Recovered Sources (DSS Tanh) - Newton Optimization, Component Time Series

Candidate Brain Component#

plot_component_summary(
    dss_meg,
    data=raw_meg,
    info=raw_meg.info,
    picks=mne.pick_types(raw_meg.info, meg="grad", eeg=False, eog=False, stim=False),
    n_components=[brain_idx],
    show=False,
)
plt.gcf().suptitle(f"Component #{brain_idx}: Candidate Brain Source")
plt.show()

# Overlay comparison for EOG
# Show how the extracted component matches the EOG channel
eog_raw = mne.io.RawArray(
    eog_ch[None, :], mne.create_info(["EOG"], raw.info["sfreq"], "eog")
)
comp_raw = mne.io.RawArray(
    sources[blink_idx : blink_idx + 1], mne.create_info(1, raw.info["sfreq"], "misc")
)

Component #2: Candidate Brain Source, Comp 2 Pattern, Comp 2 Time Course, PSD

Creating RawArray with float64 data, n_channels=1, n_times=36038
    Range : 0 ... 36037 =      0.000 ...    60.000 secs
Ready.
Creating RawArray with float64 data, n_channels=1, n_times=36038
    Range : 0 ... 36037 =      0.000 ...    60.000 secs
Ready.

EOG Overlay#

plot_signal_overlay(
    eog_raw,
    comp_raw,
    times=eog_raw.times,
    start=10,
    stop=20,
    title="EOG Channel vs Extracted Component (Time Domain)",
    show=False,
)
plt.show()

print("\nBlind Source Separation complete!")

EOG Channel vs Extracted Component (Time Domain)

Blind Source Separation complete!

Total running time of the script: (0 minutes 5.438 seconds)