Note
Go to the end to download the full example code.
Phase-Space Embedding with PhaseSPDNet#
This tutorial demonstrates how to use PhaseSPDNet for EEG classification. PhaseSPDNet applies phase-space embedding (time-delay coordinates) to capture nonlinear dynamical structure before SPDNet processing.
Introduction#
PhaseSPDNet [Carrara* et al., 2024] leverages phase-space embedding from dynamical systems theory. The key idea is that a single time series can be “unfolded” into a higher-dimensional space that reveals the underlying dynamics of the system.
Takens’ Embedding Theorem: For a dynamical system, a time-delayed embedding can reconstruct the topology of the original state space:
\[\mathbf{x}(t) \rightarrow [\mathbf{x}(t), \mathbf{x}(t-\tau), \mathbf{x}(t-2\tau), \ldots]\]
This is particularly useful for EEG, where signals reflect complex brain dynamics that may not be fully captured by linear methods.
Setup and Imports#
import warnings
import matplotlib.pyplot as plt
import numpy as np
import torch
from braindecode import EEGClassifier
from moabb.datasets import BNCI2014_001
from moabb.paradigms import MotorImagery
from sklearn.metrics import ConfusionMatrixDisplay, accuracy_score, confusion_matrix
from sklearn.preprocessing import LabelEncoder
from skorch.callbacks import EpochScoring, GradientNormClipping
from skorch.dataset import ValidSplit
from spd_learn.models import PhaseSPDNet
warnings.filterwarnings("ignore")
/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/braindecode/models/eegpt.py:497: FutureWarning: Montage name 'standard_1020' is deprecated and will be removed in MNE 1.14. Use 'colin27_1020' instead.
montage = mne.channels.make_standard_montage("standard_1020")
/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/braindecode/models/eegpt.py:1452: FutureWarning: Montage name 'standard_1020' is deprecated and will be removed in MNE 1.14. Use 'colin27_1020' instead.
montage = make_standard_montage("standard_1020")
Loading the Dataset#
dataset = BNCI2014_001()
paradigm = MotorImagery(n_classes=4)
print(f"Dataset: {dataset.code}")
print("Sampling rate: 250 Hz")
Choosing from all possible events
Dataset: BNCI2014-001
Sampling rate: 250 Hz
Understanding Phase-Space Embedding#
For a signal with n_chans channels, phase-space embedding with
order=m and lag=τ creates:
\[X_{embedded}(t) = [X(t), X(t-\tau), X(t-2\tau), \ldots, X(t-(m-1)\tau)]\]
This increases the channel dimension by a factor of m:
n_chans → n_chans * order
Choosing parameters:
order: Embedding dimension (typically 2-5)lag: Time delay in samples (often chosen via autocorrelation)
# Visualize embedding concept
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# Original 1D signal
t = np.linspace(0, 4 * np.pi, 200)
x = np.sin(t) + 0.5 * np.sin(2 * t)
ax1 = axes[0]
ax1.plot(t, x, "b-", linewidth=2)
ax1.set_xlabel("Time", fontsize=12)
ax1.set_ylabel("Amplitude", fontsize=12)
ax1.set_title("Original Signal", fontsize=14)
ax1.grid(True, alpha=0.3)
# Phase-space embedding (2D)
lag = 15 # samples
x_delayed = x[:-lag]
x_original = x[lag:]
ax2 = axes[1]
ax2.plot(x_original, x_delayed, "b-", linewidth=1, alpha=0.7)
ax2.scatter(x_original[::10], x_delayed[::10], c=t[lag::10], cmap="viridis", s=30)
ax2.set_xlabel("x(t)", fontsize=12)
ax2.set_ylabel("x(t - τ)", fontsize=12)
ax2.set_title("Phase-Space Embedding (2D)", fontsize=14)
ax2.grid(True, alpha=0.3)
ax2.set_aspect("equal")
plt.tight_layout()
plt.show()
print("Phase-space embedding reveals the underlying attractor structure!")
Phase-space embedding reveals the underlying attractor structure!
Creating the PhaseSPDNet Model#
PhaseSPDNet architecture:
PhaseDelay: Applies time-delay embedding
SPDNet: Processes the embedded signals
The embedding expands channels: 22 channels x order 3 = 66 channels
n_chans = 22
n_outputs = 4
# Phase-space parameters
order = 2 # Embedding dimension (lower for stability)
lag = 10 # Time delay (~40ms at 250Hz)
model = PhaseSPDNet(
n_chans=n_chans,
n_outputs=n_outputs,
order=order,
lag=lag,
subspacedim=22, # BiMap output dimension (half of embedded channels)
threshold=1e-4,
)
print("PhaseSPDNet Configuration:")
print(f" Original channels: {n_chans}")
print(f" Embedding order: {order}")
print(f" Time lag: {lag} samples ({lag / 250 * 1000:.1f} ms)")
print(f" Embedded channels: {n_chans * order}")
print(" Subspace dimension: 22")
print("\nModel Architecture:")
print(model)
PhaseSPDNet Configuration:
Original channels: 22
Embedding order: 2
Time lag: 10 samples (40.0 ms)
Embedded channels: 44
Subspace dimension: 22
Model Architecture:
PhaseSPDNet(
(phase): PhaseDelay()
(spdnet): SPDNet(
(cov): CovLayer()
(bimap): ParametrizedBiMap(
(parametrizations): ModuleDict(
(weight): ParametrizationList(
(0): _Orthogonal()
)
)
)
(reeig): ReEig()
(logeig): LogEig()
(classifier): Linear(in_features=253, out_features=4, bias=True)
)
)
Training the Model#
subject_id = 1
batch_size = 32
max_epochs = 100
learning_rate = 1e-4 # Low learning rate for stable SPD learning
device = "cuda" if torch.cuda.is_available() else "cpu"
print(f"\nUsing device: {device}")
# Cache configuration
cache_config = dict(
save_raw=True,
save_epochs=True,
save_array=True,
use=True,
overwrite_raw=False,
overwrite_epochs=False,
overwrite_array=False,
)
# Load data
X, labels, meta = paradigm.get_data(
dataset=dataset, subjects=[subject_id], cache_config=cache_config
)
# Encode labels
le = LabelEncoder()
y = le.fit_transform(labels)
# Split by session
train_idx = meta.query("session == '0train'").index.to_numpy()
test_idx = meta.query("session == '1test'").index.to_numpy()
print(f"\nData shape: {X.shape}")
print(
f"After embedding: ({X.shape[0]}, {n_chans * order}, {X.shape[2] - (order - 1) * lag})"
)
print(f"Training samples: {len(train_idx)}")
print(f"Test samples: {len(test_idx)}")
# Create classifier
# Note: SPD networks benefit from gradient clipping to prevent
# divergence during training on the Riemannian manifold.
clf = EEGClassifier(
model,
criterion=torch.nn.CrossEntropyLoss,
optimizer=torch.optim.Adam,
optimizer__lr=learning_rate,
train_split=ValidSplit(0.1, stratified=True, random_state=42),
batch_size=batch_size,
max_epochs=max_epochs,
callbacks=[
(
"train_acc",
EpochScoring(
"accuracy", lower_is_better=False, on_train=True, name="train_acc"
),
),
("gradient_clip", GradientNormClipping(gradient_clip_value=1.0)),
],
device=device,
verbose=1,
)
# Train
clf.fit(X[train_idx], y[train_idx])
# Evaluate
y_pred_train = clf.predict(X[train_idx])
y_pred_test = clf.predict(X[test_idx])
train_acc = accuracy_score(y[train_idx], y_pred_train)
test_acc = accuracy_score(y[test_idx], y_pred_test)
print(f"\n{'=' * 50}")
print(f"Results for Subject {subject_id}")
print(f"{'=' * 50}")
print(f"Train Accuracy: {train_acc * 100:.2f}%")
print(f"Test Accuracy: {test_acc * 100:.2f}%")
Using device: cpu
Data shape: (576, 22, 1001)
After embedding: (576, 44, 991)
Training samples: 288
Test samples: 288
epoch train_acc train_loss valid_acc valid_loss dur
------- ----------- ------------ ----------- ------------ ------
1 0.2500 1.4208 0.3448 1.3904 0.1555
2 0.2578 1.3955 0.3448 1.3786 0.1291
3 0.3203 1.3838 0.3448 1.3716 0.1346
4 0.3906 1.3748 0.3103 1.3663 0.1358
5 0.4141 1.3687 0.3448 1.3617 0.1302
6 0.3984 1.3627 0.3793 1.3566 0.1317
7 0.3867 1.3551 0.4138 1.3518 0.1317
8 0.3945 1.3503 0.3793 1.3467 0.1358
9 0.4219 1.3436 0.3103 1.3418 0.1288
10 0.4180 1.3387 0.3448 1.3358 0.1384
11 0.4336 1.3336 0.4483 1.3304 0.1315
12 0.4453 1.3286 0.4483 1.3251 0.0890
13 0.4219 1.3240 0.5517 1.3200 0.0952
14 0.4258 1.3193 0.5172 1.3157 0.1141
15 0.4453 1.3137 0.4828 1.3112 0.0994
16 0.4766 1.3108 0.4828 1.3053 0.0935
17 0.5195 1.3048 0.5172 1.3010 0.1032
18 0.5117 1.2999 0.5517 1.2969 0.1299
19 0.5000 1.2952 0.5172 1.2925 0.1284
20 0.4961 1.2914 0.5172 1.2877 0.1340
21 0.5195 1.2877 0.5517 1.2822 0.1272
22 0.5273 1.2829 0.5172 1.2774 0.1436
23 0.5547 1.2770 0.5517 1.2728 0.1436
24 0.5586 1.2729 0.5517 1.2676 0.1373
25 0.5703 1.2667 0.5862 1.2626 0.1294
26 0.5703 1.2632 0.5517 1.2577 0.1315
27 0.5703 1.2579 0.5862 1.2531 0.1309
28 0.5742 1.2533 0.5517 1.2486 0.1222
29 0.5938 1.2485 0.5862 1.2437 0.1337
30 0.6055 1.2442 0.5172 1.2388 0.1284
31 0.5977 1.2385 0.5517 1.2345 0.1225
32 0.6016 1.2355 0.5862 1.2302 0.1303
33 0.5938 1.2311 0.6207 1.2254 0.1298
34 0.5859 1.2251 0.6552 1.2205 0.1278
35 0.6016 1.2208 0.6207 1.2157 0.1361
36 0.6328 1.2154 0.6552 1.2112 0.1218
37 0.6367 1.2121 0.6897 1.2070 0.1263
38 0.6172 1.2089 0.6207 1.2034 0.1324
39 0.6406 1.2044 0.6897 1.2004 0.1225
40 0.6016 1.2008 0.5862 1.1978 0.0838
41 0.5977 1.1986 0.6207 1.1940 0.0869
42 0.6250 1.1946 0.6897 1.1886 0.1036
43 0.6523 1.1898 0.7241 1.1846 0.1104
44 0.6562 1.1863 0.7241 1.1807 0.0902
45 0.6602 1.1843 0.7586 1.1777 0.0868
46 0.6641 1.1798 0.7241 1.1741 0.1071
47 0.6641 1.1747 0.7931 1.1713 0.1014
48 0.6523 1.1733 0.7241 1.1690 0.0942
49 0.6484 1.1703 0.7586 1.1655 0.0919
50 0.6562 1.1657 0.7241 1.1613 0.0912
51 0.6875 1.1613 0.6552 1.1575 0.1095
52 0.6992 1.1579 0.6897 1.1523 0.0941
53 0.6992 1.1527 0.6897 1.1484 0.0998
54 0.6875 1.1500 0.7241 1.1452 0.1044
55 0.7031 1.1452 0.6897 1.1413 0.0917
56 0.7148 1.1427 0.6897 1.1372 0.0924
57 0.7109 1.1390 0.6207 1.1326 0.0935
58 0.7148 1.1365 0.6207 1.1297 0.1146
59 0.7188 1.1332 0.7241 1.1282 0.0963
60 0.7109 1.1303 0.7931 1.1248 0.0926
61 0.7031 1.1273 0.8276 1.1212 0.1009
62 0.7070 1.1214 0.8621 1.1185 0.1107
63 0.7188 1.1177 0.8621 1.1157 0.0967
64 0.7266 1.1152 0.7241 1.1127 0.0880
65 0.6992 1.1129 0.7586 1.1094 0.0966
66 0.7109 1.1102 0.8276 1.1061 0.1063
67 0.7188 1.1069 0.8276 1.1023 0.0959
68 0.7227 1.1033 0.7586 1.0984 0.0940
69 0.7344 1.0978 0.7241 1.0944 0.0966
70 0.7461 1.0954 0.6897 1.0906 0.1044
71 0.7461 1.0932 0.6552 1.0880 0.1031
72 0.7305 1.0895 0.6897 1.0847 0.0944
73 0.7266 1.0863 0.6897 1.0809 0.1167
74 0.7500 1.0836 0.7241 1.0784 0.0896
75 0.7422 1.0794 0.6552 1.0758 0.0927
76 0.7266 1.0754 0.6207 1.0745 0.0908
77 0.7109 1.0747 0.6207 1.0717 0.1100
78 0.7227 1.0704 0.6207 1.0679 0.1054
79 0.7266 1.0687 0.6552 1.0642 0.0913
80 0.7500 1.0671 0.7241 1.0615 0.0904
81 0.7500 1.0636 0.7241 1.0586 0.1023
82 0.7578 1.0601 0.7241 1.0574 0.1116
83 0.7461 1.0569 0.7586 1.0546 0.0965
84 0.7656 1.0545 0.6897 1.0527 0.0914
85 0.7578 1.0508 0.6897 1.0482 0.0957
86 0.7578 1.0495 0.7241 1.0451 0.1123
87 0.7500 1.0469 0.6897 1.0428 0.0935
88 0.7383 1.0434 0.6897 1.0410 0.0973
89 0.7695 1.0397 0.7586 1.0380 0.0974
90 0.7617 1.0373 0.7931 1.0366 0.1069
91 0.7578 1.0352 0.7241 1.0326 0.0943
92 0.7344 1.0311 0.6897 1.0290 0.0894
93 0.7500 1.0297 0.7241 1.0266 0.1076
94 0.7695 1.0262 0.7586 1.0240 0.1274
95 0.7656 1.0238 0.7931 1.0205 0.0910
96 0.7617 1.0208 0.7931 1.0183 0.0969
97 0.7617 1.0188 0.7241 1.0157 0.1097
98 0.7617 1.0163 0.6897 1.0127 0.1042
99 0.7539 1.0123 0.6897 1.0116 0.0896
100 0.7266 1.0085 0.6897 1.0099 0.0936
==================================================
Results for Subject 1
==================================================
Train Accuracy: 70.49%
Test Accuracy: 68.40%
Visualizing Results#
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
# Training curves
history = clf.history
epochs = range(1, len(history) + 1)
ax1 = axes[0]
ax1.plot(epochs, history[:, "train_loss"], "b-", label="Train Loss", linewidth=2)
ax1.plot(epochs, history[:, "valid_loss"], "r--", label="Valid Loss", linewidth=2)
ax1.set_xlabel("Epoch", fontsize=12)
ax1.set_ylabel("Loss", fontsize=12)
ax1.set_title("Training and Validation Loss", fontsize=14)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax2 = axes[1]
ax2.plot(epochs, history[:, "train_acc"], "b-", label="Train Acc", linewidth=2)
ax2.plot(epochs, history[:, "valid_acc"], "r--", label="Valid Acc", linewidth=2)
ax2.set_xlabel("Epoch", fontsize=12)
ax2.set_ylabel("Accuracy", fontsize=12)
ax2.set_title("Training and Validation Accuracy", fontsize=14)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_ylim([0, 1])
# Confusion matrix
ax3 = axes[2]
cm = confusion_matrix(y[test_idx], y_pred_test)
disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=le.classes_)
disp.plot(ax=ax3, cmap="Blues", values_format="d")
ax3.set_title(f"Confusion Matrix\nAccuracy: {test_acc * 100:.1f}%", fontsize=14)
plt.tight_layout()
plt.show()
Comparing Different Embedding Parameters#
The choice of order and lag affects performance.
Let’s compare different configurations.
print("\nComparing embedding parameters:")
print("-" * 50)
configs = [
{"order": 2, "lag": 5, "name": "order=2, lag=5"},
{"order": 2, "lag": 10, "name": "order=2, lag=10"},
{"order": 3, "lag": 5, "name": "order=3, lag=5"},
{"order": 3, "lag": 10, "name": "order=3, lag=10"},
]
for config in configs:
embedded_chans = n_chans * config["order"]
reduced_time = X.shape[2] - (config["order"] - 1) * config["lag"]
cov_size = embedded_chans * (embedded_chans + 1) // 2
print(
f"{config['name']:20s}: {embedded_chans} channels, "
f"{reduced_time} time points, {cov_size} features"
)
Comparing embedding parameters:
--------------------------------------------------
order=2, lag=5 : 44 channels, 996 time points, 990 features
order=2, lag=10 : 44 channels, 991 time points, 990 features
order=3, lag=5 : 66 channels, 991 time points, 2211 features
order=3, lag=10 : 66 channels, 981 time points, 2211 features
When to Use Phase-Space Embedding#
PhaseSPDNet is particularly effective when:
Nonlinear dynamics: The underlying system has complex, nonlinear behavior (e.g., neural oscillations, chaos)
Limited channels: Embedding can extract more information from fewer channels
Temporal structure: Important features span across time (captured by delay coordinates)
Considerations:
Higher
orderincreases model capacity but also parameterslagshould be chosen based on the signal’s autocorrelationReduces effective time dimension:
T_new = T - (order-1) * lag
Summary#
In this tutorial, we demonstrated:
Phase-space embedding theory and visualization
Creating PhaseSPDNet with different embedding parameters
Training and evaluating on motor imagery data
PhaseSPDNet offers a principled way to incorporate dynamical systems perspectives into EEG classification, potentially capturing nonlinear brain dynamics that linear methods miss.
Total running time of the script: (0 minutes 17.246 seconds)
