From 450719dc5a1c21aa1307736531b73a3bc730b2cc Mon Sep 17 00:00:00 2001 From: thibaut-germain Date: Fri, 5 Jun 2026 13:39:39 +0200 Subject: [PATCH 1/5] correction SGOT cost matrix, added to contributor list, move sgot example to other --- CITATION.cff | 6 ++++++ CONTRIBUTORS.md | 2 ++ examples/{ => others}/plot_sgot.py | 0 ot/sgot.py | 2 +- 4 files changed, 9 insertions(+), 1 deletion(-) rename examples/{ => others}/plot_sgot.py (100%) diff --git a/CITATION.cff b/CITATION.cff index 144ec668c..f4bdd166a 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -81,6 +81,12 @@ authors: - given-names: David family-names: Coeurjolly affiliation: CNRS, LIRIS + - given-names: Thibaut + family-names: Germain + affiliation: Ecole Polytechnique + - given-names: Sienna + family-names: O'Shea + affiliation: Ecole Polytechnique identifiers: - type: url diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index aa28dd7ab..08f481b33 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -59,6 +59,8 @@ The contributors to this library are: * [Julie Delon](https://judelo.github.io/) (GMM OT) * [Samuel Boïté](https://samuelbx.github.io/) (GMM OT) * [Nathan Neike](https://github.com/nathanneike) (Sparse EMD solver) +* [Thibaut Germain](https://thibaut-germain.github.io) (SGOT) +* Sienna O'Shea (SGOT) ## Acknowledgments diff --git a/examples/plot_sgot.py b/examples/others/plot_sgot.py similarity index 100% rename from examples/plot_sgot.py rename to examples/others/plot_sgot.py diff --git a/ot/sgot.py b/ot/sgot.py index ff67192c5..22ce5798d 100644 --- a/ot/sgot.py +++ b/ot/sgot.py @@ -124,7 +124,7 @@ def _delta_matrix_1d(Rs, Ls, Rt, Lt, nx=None, eps=1e-12): Ltn = _normalize_columns(Lt, nx=nx, eps=eps) Cr = nx.dot(nx.conj(Rsn).T, Rtn) - Cl = nx.dot(nx.conj(Lsn).T, Ltn) + Cl = nx.dot(Lsn.T, nx.conj(Ltn)) delta = nx.abs(Cr * Cl) delta = nx.clip(delta, 0.0, 1.0) From 6e0bc73cc3adbac2c0007aa5e12f7a303add34e7 Mon Sep 17 00:00:00 2001 From: Sienna O'Shea Date: Fri, 19 Jun 2026 17:00:24 +0200 Subject: [PATCH 2/5] updated graphs --- examples/plot_sgot.py | 576 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 576 insertions(+) create mode 100644 examples/plot_sgot.py diff --git a/examples/plot_sgot.py b/examples/plot_sgot.py new file mode 100644 index 000000000..3d66d5bbf --- /dev/null +++ b/examples/plot_sgot.py @@ -0,0 +1,576 @@ +#!/usr/bin/env python +# coding: utf-8 + +""" +========================================= +Spectral-Grassmann OT on dynamical systems operators +========================================= + +This example presents a synthetic example of Spectral Grassmannian-Wasserstein +Optimal Transport (SGOT) on linear dynamical systems. + +We consider a signal formed by the sum of two damped oscillatory modes evolving +along a rotated direction in the plane. The signal is then associated with an +underlying continuous linear dynamical system, and we study how its spectral +representation varies under rotation. The SGOT cost and metric are used to +compare the reference and rotated systems. + +.. [83] T. Germain; R. Flamary; V. R. Kostic; K. Lounici, A Spectral-Grassmann + Wasserstein Metric for Operator Representations of Dynamical Systems, + arXiv preprint arXiv:2509.24920, 2025. +""" + +# Authors: Sienna O'Shea +# Thibaut Germain +# +# License: MIT License + +import numpy as np +import matplotlib.pyplot as plt + +from ot.sgot import sgot_metric, sgot_cost_matrix + +from scipy.linalg import eig + + +# sampling parameters and time grid +fs = 50 +max_t = 5 +time = np.linspace(0, max_t, fs * max_t) +dt = 1 / fs + + +# %% +# Example: rotating a linear dynamical system in 3D +# ------------------------------------------------- +# +# 1. Build a simple observed signal +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# We begin by assuming that the observed signal is made of two oscillatory +# components: +# +# .. math:: +# +# x_{\text{ref}}(t)=e^{-\tau_1 t}\cos(2\pi\omega_1 t)\,\vec e(\theta) +# \;+\; +# e^{-\tau_2 t}\cos(2\pi\omega_2 t)\,\vec e(\theta), +# +# where :math:`\vec e(\theta)\in\mathbb{R}^2` is a fixed real vector. Thus, +# :math:`x(t)` evolves along the one-dimensional subspace spanned by +# :math:`\vec e(\theta)`, while its time dependence exhibits oscillatory and +# dissipative behaviour. + +tau_0 = np.array([0.08, 0.18]) +freq_0 = np.array([1.0, 2.0]) +theta_0 = np.pi / 4 + + +def generate_data(time, tau, freq, theta): + t_ = np.sin(2 * np.pi * freq[None, :] * time[:, None]) * np.exp( + -tau[None, :] * time[:, None] + ) + t_ = t_.sum(axis=1) + traj_0 = np.zeros((t_.shape[0], 2)) + traj_0[:, 0] = t_ + rotation_matrix = np.array( + [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]] + ) + traj_0 = traj_0 @ rotation_matrix.T + return traj_0 + + +traj_0 = generate_data(time, tau_0, freq_0, theta_0) + + +# plot the observed signal components and their sum +plt.figure(figsize=(10, 4)) +plt.plot(time, traj_0, label="base trajectory", linewidth=2) +plt.xlabel("time") +plt.ylabel("amplitude") +plt.legend() +plt.title(r"Observed scalar signal along $\vec{e}(\theta)$") +plt.show() + + +# %% +# 2. Interpret the signal as coming from a continuous linear dynamical system +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# We assume that :math:`x(t)` is generated by an underlying continuous linear +# dynamical system. Since the observed signal is a superposition of two +# sinusoidal modes, the corresponding linear dynamics are naturally described +# by a fourth-order model. We therefore introduce the state vector +# +# .. math:: +# +# z(t)= +# \begin{pmatrix} +# x_1(t)\\ +# x_2(t)\\ +# \vdots\\ +# x_1^{(3)}(t)\\ +# x_2^{(3)}(t) +# \end{pmatrix} +# \in\mathbb{R}^8. +# +# where :math:`x^{(n)}(t)` denotes the n-th derivative of :math:`x(t)`. +# +# This allows us to rewrite the dynamics as a first-order linear system: +# +# .. math:: +# +# \dot{z}(t)=Az(t), +# +# where :math:`A\in\mathbb{R}^{8\times 8}`. Its solution is then given by +# +# .. math:: +# +# z(t)=e^{tA}z_0. + +fig = plt.figure(figsize=(9, 6)) +ax = fig.add_subplot(projection="3d") + +ax.plot(time, traj_0[:, 0], traj_0[:, 1]) +ax.set_xlabel("time") +ax.set_ylabel("x₁(t)") +ax.set_title("Observed trajectory in time") + +ax.text2D(1.08, 0.5, "x₂(t)", transform=ax.transAxes, rotation=90, va="center") + +plt.show() + + +# %% +# 3. Sampling and preprocessing discrete trajectories of the dynamical system +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# We now have a bridge between the continuous system and the operator we later +# aim to infer from sampled data. Since in practice we do not observe the full +# continuous trajectory, we work instead with discrete samples of the signal. +# We take snapshots at uniform time intervals :math:`\Delta t`, and write the +# sampled signal as +# +# .. math:: +# +# S= +# \begin{pmatrix} +# x_1(0) & x_2(0)\\ +# x_1(\Delta t) & x_2(\Delta t)\\ +# \vdots\\ +# x_1(N\Delta t) & x_2(N\Delta t)\\ +# \end{pmatrix} +# +# The goal is now to use these observations to recover the operator governing +# the evolution. To do this, we augment the signal :math:`s` using a sliding +# window of length :math:`w`. For each :math:`k`, define +# +# .. math:: +# +# z_k = +# \begin{pmatrix} +# s_k\\ +# s_{k+1}\\ +# \vdots\\ +# s_{k+w-1} +# \end{pmatrix} +# +# We then form the data matrices +# +# .. math:: +# +# X= +# \begin{pmatrix} +# z_1\\ +# z_2\\ +# \vdots\\ +# z_{N-w} +# \end{pmatrix}, +# \qquad +# Y= +# \begin{pmatrix} +# z_2\\ +# z_3\\ +# \vdots\\ +# z_{N-w+1} +# \end{pmatrix}, +# +# so that :math:`X` contains the present windowed states and :math:`Y` the +# corresponding shifted future states. + + +# build a 4-dimensional state using delay embedding +def augment(traj, window_length=2): + Z = np.lib.stride_tricks.sliding_window_view(traj, (window_length, 1)) + Z = Z.reshape(Z.shape[0], -1) + return Z + + +# create the embedded state matrix Z +Z = augment(traj_0, 4) +Z.shape + +# inspect one embedded state vector +Z[0] + +# create X and Y for the SGOT metric +X = Z[:-1] +Y = Z[1:] + +# inspect shapes of X and Y +print("X shape:", X.shape) +print("Y shape:", Y.shape) + + +# %% +# 4. Estimate the discrete-time operator +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# We now identify the operator that maps :math:`X` to :math:`Y`. From +# +# .. math:: +# +# \dot{z}=Az, +# +# we have +# +# .. math:: +# +# z(t+\Delta t)=e^{\Delta tA}z(t). +# +# Setting +# +# .. math:: +# +# B=e^{\Delta tA}, +# +# the corresponding discrete-time evolution is governed by :math:`B`, and we +# seek the best linear map satisfying +# +# .. math:: +# +# Y\approx X B^T. +# +# Equivalently, we solve the optimisation problem +# +# .. math:: +# +# \min_B \|Y-XB\|^2. +# +# We want to recover the best rank-:math:`r` operator, whose estimator is +# defined as follows: +# +# .. math:: +# +# B = C_{xx}^{-\frac{1}{2}}[C_{xx}^{-\frac{1}{2}}C_{xy}]_r +# \quad \text{s.t} \quad C_{xx} = X^T X \quad \text{and} \quad C_{xy} = X^TY. +# +# Here :math:`[\cdot]_r` denotes the best rank-:math:`r` estimator obtained via +# SVD decomposition. [2] +# +# [2] Kostic, V., Novelli, P., Maurer, A., Ciliberto, C., Rosasco, L. and +# Pontil, M., 2022. Learning dynamical systems via Koopman operator regression +# in reproducing kernel Hilbert spaces. Advances in Neural Information +# Processing Systems, 35, pp.4017-4031. + + +def estimator(X, Y, rank=4): + # X: (n_samples, n_features) + # Y: (n_samples, n_features) + + # estimate operator + cxx = X.T @ X + U, S, Vt = np.linalg.svd(cxx) + S_inv = np.divide(1, S, out=np.zeros_like(S), where=S != 0) + cxx_inv_half = Vt.T @ np.diag(np.sqrt(S_inv)) @ U.T + cxy = X.T @ Y + T = cxx_inv_half @ cxy + U, S, Vt = np.linalg.svd(T) + S[rank:] = 0 + T_rank = U @ np.diag(S) @ Vt + T = cxx_inv_half @ T_rank + + # estimate spectral decomposition + val, vl, vr = eig(T, left=True, right=True) + sort_idx = np.argsort(np.abs(val))[::-1] + val = val[sort_idx][:rank] + vl = vl[:, sort_idx][:, :rank] + vr = vr[:, sort_idx][:, :rank] + + return T, {"eig_val": val, "eig_vec_left": vl, "eig_vec_right": vr} + + +B_0, B_0_spec = estimator(X, Y, rank=4) +Y_pred = X @ B_0 + +plt.figure(figsize=(10, 4)) +plt.plot(Y[:, 0], label="true") +plt.plot(Y_pred[:, 0], "--", label="predicted") +plt.xlabel("sample index") +plt.ylabel("first state coordinate") +plt.title("True Signal vs Predicted Signal") +plt.legend() +plt.show() + + +# %% +# The predicted signal is nearly indistinguishable from the true signal, +# indicating that the estimated operator accurately captures the observed +# dynamics. + +# %% +# 6. Recover continuous-time spectral information from the discrete operator +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# To recover the continuous generator :math:`A`, we study the spectral +# structure of :math:`B`. We diagonalise :math:`B` as +# +# .. math:: +# +# B=PDP^{-1}, +# +# where +# +# .. math:: +# +# D=\operatorname{diag}(\mu_1,\dots,\mu_n). +# +# The continuous-time eigenvalues are of the form +# +# .. math:: +# +# \lambda_k=-\tau_k+2\pi i\,\omega_k, +# \qquad k\in\{1,2\}, +# +# and the corresponding eigenvalues of :math:`B` are +# +# .. math:: +# +# \mu_k=e^{\Delta t\lambda_k} +# =e^{\Delta t(-\tau_k+2\pi i\omega_k)}. +# +# Since :math:`B=e^{\Delta tA}`, we recover :math:`A` by taking the logarithm: +# +# .. math:: +# +# A=P\,\frac{\log(D)}{\Delta t}\,P^{-1}. + +D_0 = np.log(B_0_spec["eig_val"]) * fs +L_0 = B_0_spec["eig_vec_left"] +R_0 = B_0_spec["eig_vec_right"] + +recovered_freqs = D_0.imag / (2 * np.pi) +mask = recovered_freqs > 0 +recovered_freqs = recovered_freqs[mask] +decay = -D_0.real[mask] +print(f"First mode: frequency: {recovered_freqs[0]:.2f} Hz -- decay: {decay[0]:.2f}") +print(f"Second mode: frequency: {recovered_freqs[1]:.2f} Hz -- decay: {decay[1]:.2f}") + + +# %% +# Introduction to SGOT for linear operators +# ----------------------------------------- +# +# To compare two linear operators through their spectral structure, we use the +# SGOT framework introduced in Theorem 1 of [1]. For a non-defective +# finite-rank operator :math:`T \in S_r(\mathcal H)`, the theorem associates a +# discrete spectral measure +# +# .. math:: +# +# \mu(T) \triangleq \sum_{j\in[\ell]} +# \frac{m_j}{m_{\mathrm{tot}}}\,\delta_{(\lambda_j,\mathcal V_j)}, +# +# where :math:`\lambda_j` are the eigenvalues of :math:`T`, :math:`m_j` their +# algebraic multiplicities, and :math:`\mathcal V_j` the corresponding +# eigenspaces. Thus, each spectral component of the operator is represented by +# an atom of the form +# +# .. math:: +# +# (\lambda_j,\mathcal V_j), +# +# combining one eigenvalue with its associated invariant subspace. +# +# Theorem 1 then defines a ground cost between two such atoms by combining a +# spectral discrepancy and a geometric discrepancy: +# +# .. math:: +# +# d_\eta\big((\lambda,\mathcal V),(\lambda',\mathcal V')\big) +# \triangleq +# \eta\,|\lambda-\lambda'| + (1-\eta)\, d_{\mathcal G}(\mathcal V,\mathcal V'), +# +# where :math:`d_{\mathcal G}` denotes the grassmann distance between +# eigenspaces and :math:`\eta\in(0,1)` balances the contribution of eigenvalues +# and eigenspaces. +# +# The SGOT distance between two operators :math:`T` and :math:`T'` is then the +# Wasserstein distance between their associated spectral measures: +# +# .. math:: +# +# d_S(T,T') = W_{d_\eta,p}\big(\mu(T),\mu(T')\big). +# +# In this way, SGOT compares linear operators by optimally matching their +# spectral atoms, taking into account both the location of eigenvalues and the +# relative geometry of their eigenspaces. + +# %% +# SGOT distance versus rotation angle +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# We compare the reference signal with a rotated version obtained by changing +# only the observation direction. The shifted signal is +# +# .. math:: +# +# x_{\mathrm{shift}}^{\mathrm{rot}}(t;\theta) +# = +# \sum_{i=1}^{2} +# e^{-\tau_i t}\cos(2\pi\omega_i t)\,\vec e({\color{red}\theta}), +# +# while the reference signal is recovered at :math:`\theta=\theta_0`. Thus, +# this experiment isolates the effect of rotating the underlying one-dimensional +# subspace in the observation plane. + +thetas = np.linspace(0, np.pi / 2, 51) +rotation_scores = [] + +for theta in thetas: + Z = augment(generate_data(time, tau_0, freq_0, theta), 4) + B, B_spec = estimator(Z[:-1], Z[1:]) + D = np.log(B_spec["eig_val"]) * fs + L = B_spec["eig_vec_left"] + R = B_spec["eig_vec_right"] + rotation_scores.append( + sgot_metric(D_0, R_0, L_0, D, R, L, eta=0.9, grassmann_metric="chordal") + ) + +fig, ax = plt.subplots(figsize=(7, 4)) +ax.plot(thetas, rotation_scores, linewidth=1.8) +ax.axvline(theta_0, color="gray", linestyle="--", linewidth=0.8) +ax.set_xlabel(r"Rotation angle $\theta$ (rad)") +ax.set_ylabel(r"$d_S$") +ax.set_title("SGOT distance vs. rotation angle") +fig.tight_layout() +plt.show() + +# %% +# Comparison across Grassmannian metrics for SGOT distance versus rotation angle +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +metrics = ["chordal", "geodesic", "procrustes", "martin"] +styles = {"chordal": "-", "geodesic": "--", "procrustes": "-.", "martin": ":"} +rotation_results = {m: [] for m in metrics} + +for theta in thetas: + Z = augment(generate_data(time, tau_0, freq_0, theta), 4) + B, B_spec = estimator(Z[:-1], Z[1:]) + D = np.log(B_spec["eig_val"]) * fs + L = B_spec["eig_vec_left"] + R = B_spec["eig_vec_right"] + for m in metrics: + rotation_results[m].append( + sgot_metric(D_0, R_0, L_0, D, R, L, eta=0.9, grassmann_metric=m) + ) + +fig, ax = plt.subplots(figsize=(7, 4)) +for m in metrics: + ax.plot(thetas, rotation_results[m], styles[m], label=m, linewidth=1.8) +ax.axvline( + theta_0, color="gray", linestyle="--", linewidth=0.8, label=r"$\theta_0 = \pi/4$" +) +ax.set_xlabel(r"Rotation angle $\theta$ (rad)") +ax.set_ylabel(r"$d_S$") +ax.set_title("SGOT distance vs. rotation angle across Grassmannian metrics") +ax.legend() +fig.tight_layout() +plt.show() + +# %% +# SGOT distance versus frequency +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# In this experiment, we keep the reference direction fixed and perturb one of +# the oscillatory modes. The shifted signal is +# +# .. math:: +# +# x_{\mathrm{shift}}^{\omega}(t) +# = +# e^{-\tau_1 t}\cos(2\pi\omega_1 t)\,\vec e(\theta_0) +# \;+\; +# e^{-\tau_2 t}\cos(2\pi{\color{red}\omega_2'} t)\,\vec e(\theta_0), +# +# where only the second frequency is modified. We then study how the SGOT +# distance changes as a function of the perturbed frequency :math:`\omega_2'`. + +omegas = np.linspace(0.5, 3.0, 21) +frequency_scores = {m: [] for m in metrics} + +for omega in omegas: + Z = augment(generate_data(time, tau_0, np.array([freq_0[0], omega]), theta_0), 4) + B, B_spec = estimator(Z[:-1], Z[1:]) + D = np.log(B_spec["eig_val"]) * fs + L = B_spec["eig_vec_left"] + R = B_spec["eig_vec_right"] + for m in metrics: + frequency_scores[m].append( + sgot_metric(D_0, R_0, L_0, D, R, L, eta=0.9, grassmann_metric=m) + ) + +fig, ax = plt.subplots(figsize=(7, 4)) +for m in metrics: + ax.plot(omegas, frequency_scores[m], styles[m], label=m, linewidth=1.8) +ax.axvline( + freq_0[1], color="gray", linestyle="--", linewidth=0.8, label=r"$\omega_2 = 2.0$ Hz" +) +ax.set_xlabel(r"Frequency $\omega_2'$ (Hz)") +ax.set_ylabel(r"$d_S$") +ax.set_title("SGOT distance vs. frequency across Grassmannian metrics") +ax.legend() +fig.tight_layout() +plt.show() + +# %% +# SGOT distance versus decay +# ~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# We now study the effect of changing the decay rate while keeping the +# observation direction fixed. The shifted signal is +# +# .. math:: +# +# x_{\mathrm{shift}}^{\mathrm{decay}}(t;\tau) +# = +# e^{-{\color{red}\tau} t}\cos(2\pi\omega_1 t)\,\vec e(\theta_0) +# \;+\; +# e^{-{\color{red}\tau} t}\cos(2\pi\omega_2' t)\,\vec e(\theta_0). +# +# In this way, both modes share the same modified decay parameter +# :math:`\tau`, allowing us to isolate the influence of dissipation on the SGOT +# distance. +taus = np.linspace(0.1, 3.0, 21) +decay_scores = {m: [] for m in metrics} + +for tau in taus: + Z = augment(generate_data(time, np.array([tau, tau]), freq_0, theta_0), 4) + B, B_spec = estimator(Z[:-1], Z[1:]) + D = np.log(B_spec["eig_val"]) * fs + L = B_spec["eig_vec_left"] + R = B_spec["eig_vec_right"] + for m in metrics: + decay_scores[m].append( + sgot_metric(D_0, R_0, L_0, D, R, L, eta=0.9, grassmann_metric=m) + ) + +fig, ax = plt.subplots(figsize=(7, 4)) +for m in metrics: + ax.plot(taus, decay_scores[m], styles[m], label=m, linewidth=1.8) +ax.set_xlabel(r"Decay rate $\tau$") +ax.set_ylabel(r"$d_S$") +ax.set_title("SGOT distance vs. decay across Grassmannian metrics") +ax.legend() +fig.tight_layout() +plt.show() From 8d0ddbdd73e5b2e6a4853b696ade8bac3dd781d7 Mon Sep 17 00:00:00 2001 From: Sienna O'Shea Date: Thu, 2 Jul 2026 09:32:42 +0200 Subject: [PATCH 3/5] fix plots --- examples/plot_sgot.py | 64 ++++++++++++++++++++++++++++++++----------- 1 file changed, 48 insertions(+), 16 deletions(-) diff --git a/examples/plot_sgot.py b/examples/plot_sgot.py index 3d66d5bbf..87f7ff3c8 100644 --- a/examples/plot_sgot.py +++ b/examples/plot_sgot.py @@ -274,12 +274,8 @@ def augment(traj, window_length=2): # Processing Systems, 35, pp.4017-4031. -def estimator(X, Y, rank=4): - # X: (n_samples, n_features) - # Y: (n_samples, n_features) - - # estimate operator - cxx = X.T @ X +def estimator(X, Y, rank=4, eps=1e-8): + cxx = X.T @ X + eps * np.eye(X.shape[1]) U, S, Vt = np.linalg.svd(cxx) S_inv = np.divide(1, S, out=np.zeros_like(S), where=S != 0) cxx_inv_half = Vt.T @ np.diag(np.sqrt(S_inv)) @ U.T @@ -416,6 +412,24 @@ def estimator(X, Y, rank=4): # spectral atoms, taking into account both the location of eigenvalues and the # relative geometry of their eigenspaces. +# %% +# A wider delay window for the SGOT experiments below +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The window of length 4 used above is enough to identify a single reference +# operator, but the experiments below also probe signals whose two modes +# nearly coincide in frequency (e.g. :math:`\omega_2'\to\omega_1`). Telling +# such near-degenerate modes apart requires the delay embedding to span +# enough time to "see" their differing decay, so we re-embed the reference +# signal with a longer window before running the sweeps. + +sgot_window = 10 +Z = augment(traj_0, sgot_window) +_, B_0_spec_sgot = estimator(Z[:-1], Z[1:]) +D_0_sgot = np.log(B_0_spec_sgot["eig_val"]) * fs +L_0_sgot = B_0_spec_sgot["eig_vec_left"] +R_0_sgot = B_0_spec_sgot["eig_vec_right"] + # %% # SGOT distance versus rotation angle # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -438,13 +452,15 @@ def estimator(X, Y, rank=4): rotation_scores = [] for theta in thetas: - Z = augment(generate_data(time, tau_0, freq_0, theta), 4) + Z = augment(generate_data(time, tau_0, freq_0, theta), sgot_window) B, B_spec = estimator(Z[:-1], Z[1:]) D = np.log(B_spec["eig_val"]) * fs L = B_spec["eig_vec_left"] R = B_spec["eig_vec_right"] rotation_scores.append( - sgot_metric(D_0, R_0, L_0, D, R, L, eta=0.9, grassmann_metric="chordal") + sgot_metric( + D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric="chordal" + ) ) fig, ax = plt.subplots(figsize=(7, 4)) @@ -465,21 +481,27 @@ def estimator(X, Y, rank=4): rotation_results = {m: [] for m in metrics} for theta in thetas: - Z = augment(generate_data(time, tau_0, freq_0, theta), 4) + Z = augment(generate_data(time, tau_0, freq_0, theta), sgot_window) B, B_spec = estimator(Z[:-1], Z[1:]) D = np.log(B_spec["eig_val"]) * fs L = B_spec["eig_vec_left"] R = B_spec["eig_vec_right"] for m in metrics: rotation_results[m].append( - sgot_metric(D_0, R_0, L_0, D, R, L, eta=0.9, grassmann_metric=m) + sgot_metric( + D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric=m + ) ) fig, ax = plt.subplots(figsize=(7, 4)) for m in metrics: ax.plot(thetas, rotation_results[m], styles[m], label=m, linewidth=1.8) ax.axvline( - theta_0, color="gray", linestyle="--", linewidth=0.8, label=r"$\theta_0 = \pi/4$" + theta_0, + color="gray", + linestyle="--", + linewidth=0.8, + label=r"$\theta_0 = \pi/4$ (reference)", ) ax.set_xlabel(r"Rotation angle $\theta$ (rad)") ax.set_ylabel(r"$d_S$") @@ -510,21 +532,29 @@ def estimator(X, Y, rank=4): frequency_scores = {m: [] for m in metrics} for omega in omegas: - Z = augment(generate_data(time, tau_0, np.array([freq_0[0], omega]), theta_0), 4) + Z = augment( + generate_data(time, tau_0, np.array([freq_0[0], omega]), theta_0), sgot_window + ) B, B_spec = estimator(Z[:-1], Z[1:]) D = np.log(B_spec["eig_val"]) * fs L = B_spec["eig_vec_left"] R = B_spec["eig_vec_right"] for m in metrics: frequency_scores[m].append( - sgot_metric(D_0, R_0, L_0, D, R, L, eta=0.9, grassmann_metric=m) + sgot_metric( + D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric=m + ) ) fig, ax = plt.subplots(figsize=(7, 4)) for m in metrics: ax.plot(omegas, frequency_scores[m], styles[m], label=m, linewidth=1.8) ax.axvline( - freq_0[1], color="gray", linestyle="--", linewidth=0.8, label=r"$\omega_2 = 2.0$ Hz" + freq_0[1], + color="gray", + linestyle="--", + linewidth=0.8, + label=r"$\omega_2 = 2.0$ Hz (reference)", ) ax.set_xlabel(r"Frequency $\omega_2'$ (Hz)") ax.set_ylabel(r"$d_S$") @@ -555,14 +585,16 @@ def estimator(X, Y, rank=4): decay_scores = {m: [] for m in metrics} for tau in taus: - Z = augment(generate_data(time, np.array([tau, tau]), freq_0, theta_0), 4) + Z = augment(generate_data(time, np.array([tau, tau]), freq_0, theta_0), sgot_window) B, B_spec = estimator(Z[:-1], Z[1:]) D = np.log(B_spec["eig_val"]) * fs L = B_spec["eig_vec_left"] R = B_spec["eig_vec_right"] for m in metrics: decay_scores[m].append( - sgot_metric(D_0, R_0, L_0, D, R, L, eta=0.9, grassmann_metric=m) + sgot_metric( + D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric=m + ) ) fig, ax = plt.subplots(figsize=(7, 4)) From f48d6cea11f357c2c12d0a28430b541025882854 Mon Sep 17 00:00:00 2001 From: thibaut-germain Date: Thu, 2 Jul 2026 13:21:01 +0200 Subject: [PATCH 4/5] update example --- examples/others/plot_sgot.py | 249 +++++++------- examples/plot_sgot.py | 608 ----------------------------------- 2 files changed, 130 insertions(+), 727 deletions(-) delete mode 100644 examples/plot_sgot.py diff --git a/examples/others/plot_sgot.py b/examples/others/plot_sgot.py index 0f4ee1dee..4dc3b6a6e 100644 --- a/examples/others/plot_sgot.py +++ b/examples/others/plot_sgot.py @@ -66,6 +66,10 @@ theta_0 = np.pi / 4 +def rotation_matrix(theta): + return np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]) + + def generate_data(time, tau, freq, theta): t_ = np.sin(2 * np.pi * freq[None, :] * time[:, None]) * np.exp( -tau[None, :] * time[:, None] @@ -73,26 +77,24 @@ def generate_data(time, tau, freq, theta): t_ = t_.sum(axis=1) traj_0 = np.zeros((t_.shape[0], 2)) traj_0[:, 0] = t_ - rotation_matrix = np.array( - [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]] - ) - traj_0 = traj_0 @ rotation_matrix.T + R_ = rotation_matrix(theta) + traj_0 = traj_0 @ R_.T return traj_0 traj_0 = generate_data(time, tau_0, freq_0, theta_0) +traj_0_proj = traj_0 @ rotation_matrix(theta_0)[:, 0] # plot the observed signal components and their sum plt.figure(figsize=(10, 4)) -plt.plot(time, traj_0, label="base trajectory", linewidth=2) +plt.plot(time, traj_0_proj, label="projected trajectory", linewidth=2) plt.xlabel("time") plt.ylabel("amplitude") plt.legend() plt.title(r"Observed scalar signal along $\vec{e}(\theta)$") plt.show() - # %% # 2. Interpret the signal as coming from a continuous linear dynamical system # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -274,12 +276,8 @@ def augment(traj, window_length=2): # Processing Systems, 35, pp.4017-4031. -def estimator(X, Y, rank=4): - # X: (n_samples, n_features) - # Y: (n_samples, n_features) - - # estimate operator - cxx = X.T @ X +def estimator(X, Y, rank=4, eps=1e-8): + cxx = X.T @ X + eps * np.eye(X.shape[1]) U, S, Vt = np.linalg.svd(cxx) S_inv = np.divide(1, S, out=np.zeros_like(S), where=S != 0) cxx_inv_half = Vt.T @ np.diag(np.sqrt(S_inv)) @ U.T @@ -416,6 +414,24 @@ def estimator(X, Y, rank=4): # spectral atoms, taking into account both the location of eigenvalues and the # relative geometry of their eigenspaces. +# %% +# A wider delay window for the SGOT experiments below +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The window of length 4 used above is enough to identify a single reference +# operator, but the experiments below also probe signals whose two modes +# nearly coincide in frequency (e.g. :math:`\omega_2'\to\omega_1`). Telling +# such near-degenerate modes apart requires the delay embedding to span +# enough time to "see" their differing decay, so we re-embed the reference +# signal with a longer window before running the sweeps. + +sgot_window = 10 +Z = augment(traj_0, sgot_window) +_, B_0_spec_sgot = estimator(Z[:-1], Z[1:]) +D_0_sgot = np.log(B_0_spec_sgot["eig_val"]) * fs +L_0_sgot = B_0_spec_sgot["eig_vec_left"] +R_0_sgot = B_0_spec_sgot["eig_vec_right"] + # %% # SGOT distance versus rotation angle # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -434,52 +450,66 @@ def estimator(X, Y, rank=4): # this experiment isolates the effect of rotating the underlying one-dimensional # subspace in the observation plane. -thetas = np.linspace(0, np.pi / 2, 50) -lst = [] -for i, theta in enumerate(thetas): - traj = generate_data(time, tau_0, freq_0, theta) - Z = augment(traj, 4) - X = Z[:-1] - Y = Z[1:] - B, B_spec = estimator(X, Y, rank=4) - D, R, L = B_spec["eig_val"], B_spec["eig_vec_right"], B_spec["eig_vec_left"] - D = np.log(D) * fs - lst.append(sgot_metric(D_0, R_0, L_0, D, R, L, eta=0.01)) - -plt.figure(figsize=(8, 5)) -plt.plot(thetas, lst) -plt.xlabel("theta") -plt.ylabel("SGOT distance") -plt.title("SGOT distance vs rotation angle") +thetas = np.linspace(0, np.pi / 2, 51) +rotation_scores = [] + +for theta in thetas: + Z = augment(generate_data(time, tau_0, freq_0, theta), sgot_window) + B, B_spec = estimator(Z[:-1], Z[1:]) + D = np.log(B_spec["eig_val"]) * fs + L = B_spec["eig_vec_left"] + R = B_spec["eig_vec_right"] + rotation_scores.append( + sgot_metric( + D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric="chordal" + ) + ) + +fig, ax = plt.subplots(figsize=(7, 4)) +ax.plot(thetas, rotation_scores, linewidth=1.8) +ax.axvline(theta_0, color="gray", linestyle="--", linewidth=0.8) +ax.set_xlabel(r"Rotation angle $\theta$ (rad)") +ax.set_ylabel(r"$d_S$") +ax.set_title("SGOT distance vs. rotation angle") +fig.tight_layout() plt.show() # %% # Comparison across Grassmannian metrics for SGOT distance versus rotation angle # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -thetas = np.linspace(0, np.pi / 2, 50) -lst = [] -for i, theta in enumerate(thetas): - traj = generate_data(time, tau_0, freq_0, theta) - Z = augment(traj, 4) - X = Z[:-1] - Y = Z[1:] - B, B_spec = estimator(X, Y, rank=4) - D, R, L = B_spec["eig_val"], B_spec["eig_vec_right"], B_spec["eig_vec_left"] - D = np.log(D) * fs - lst1 = [] - for name in ["chordal", "martin", "geodesic", "procrustes"]: - lst1.append(sgot_metric(D_0, R_0, L_0, D, R, L, eta=0.9, grassmann_metric=name)) - lst.append(lst1) -lst2 = np.array(lst) -plt.figure(figsize=(8, 5)) -for i, name in enumerate(["chordal", "martin", "geodesic", "procrustes"]): - plt.plot(thetas, lst2[:, i], label=name) - -plt.xlabel("theta") -plt.ylabel("SGOT distance") -plt.title("SGOT distance vs rotation angle") -plt.legend() +metrics = ["chordal", "geodesic", "procrustes", "martin"] +styles = {"chordal": "-", "geodesic": "--", "procrustes": "-.", "martin": ":"} +rotation_results = {m: [] for m in metrics} + +for theta in thetas: + Z = augment(generate_data(time, tau_0, freq_0, theta), sgot_window) + B, B_spec = estimator(Z[:-1], Z[1:]) + D = np.log(B_spec["eig_val"]) * fs + L = B_spec["eig_vec_left"] + R = B_spec["eig_vec_right"] + for m in metrics: + rotation_results[m].append( + sgot_metric( + D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric=m + ) + ) + +fig, ax = plt.subplots(figsize=(7, 4)) +for m in metrics: + ax.plot(thetas, rotation_results[m], styles[m], label=m, linewidth=1.8) +ax.axvline( + theta_0, + color="gray", + linestyle="--", + linewidth=0.8, + label=r"$\theta_0 = \pi/4$ (reference)", +) +ax.set_xlabel(r"Rotation angle $\theta$ (rad)") +ax.set_ylabel(r"$d_S$") +ax.set_title("SGOT distance vs. rotation angle across Grassmannian metrics") +ax.legend() +fig.tight_layout() plt.show() # %% @@ -501,38 +531,38 @@ def estimator(X, Y, rank=4): # distance changes as a function of the perturbed frequency :math:`\omega_2'`. omegas = np.linspace(0.5, 3.0, 21) -methods = ["chordal", "martin", "geodesic", "procrustes"] -scores_omega = [] -theta = theta_0 +frequency_scores = {m: [] for m in metrics} -eta_fixed = 0.9 for omega in omegas: - freq_1 = np.array([freq_0[0], omega]) - traj = generate_data(time, tau_0, freq_1, theta) - Z = augment(traj, 4) - X = Z[:-1] - Y = Z[1:] - - B, B_spec = estimator(X, Y, rank=4) - D, R, L = B_spec["eig_val"], B_spec["eig_vec_right"], B_spec["eig_vec_left"] - D = np.log(D) * fs - - row = [] - for name in methods: - row.append( - sgot_metric(D_0, R_0, L_0, D, R, L, eta=eta_fixed, grassmann_metric=name) + Z = augment( + generate_data(time, tau_0, np.array([freq_0[0], omega]), theta_0), sgot_window + ) + B, B_spec = estimator(Z[:-1], Z[1:]) + D = np.log(B_spec["eig_val"]) * fs + L = B_spec["eig_vec_left"] + R = B_spec["eig_vec_right"] + for m in metrics: + frequency_scores[m].append( + sgot_metric( + D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric=m + ) ) - scores_omega.append(row) - -scores_omega = np.array(scores_omega) -plt.figure(figsize=(8, 5)) -for i, name in enumerate(methods): - plt.plot(omegas, scores_omega[:, i], label=name) -plt.xlabel("omega") -plt.ylabel("SGOT distance") -plt.title("SGOT distance vs omega") -plt.legend() +fig, ax = plt.subplots(figsize=(7, 4)) +for m in metrics: + ax.plot(omegas, frequency_scores[m], styles[m], label=m, linewidth=1.8) +ax.axvline( + freq_0[1], + color="gray", + linestyle="--", + linewidth=0.8, + label=r"$\omega_2 = 2.0$ Hz (reference)", +) +ax.set_xlabel(r"Frequency $\omega_2'$ (Hz)") +ax.set_ylabel(r"$d_S$") +ax.set_title("SGOT distance vs. frequency across Grassmannian metrics") +ax.legend() +fig.tight_layout() plt.show() # %% @@ -553,47 +583,28 @@ def estimator(X, Y, rank=4): # In this way, both modes share the same modified decay parameter # :math:`\tau`, allowing us to isolate the influence of dissipation on the SGOT # distance. -decays = np.linspace(0.1, 3.0, 20) # adjust range as needed -methods = ["chordal", "martin", "geodesic", "procrustes"] -scores_decay = [] -theta = theta_0 - -for tau in decays: - freq_1 = np.array([freq_0[0], recovered_freqs[1]]) - tau_1 = np.array([tau, tau]) # or whatever structure your generator expects - - traj = generate_data(time, tau_1, freq_1, theta) - Z = augment(traj, 4) - X = Z[:-1] - Y = Z[1:] - - B, B_spec = estimator(X, Y, rank=4) - D, R, L = B_spec["eig_val"], B_spec["eig_vec_right"], B_spec["eig_vec_left"] - D = np.log(D) * fs - - row = [] - for name in methods: - row.append( +taus = np.linspace(0.1, 3.0, 21) +decay_scores = {m: [] for m in metrics} + +for tau in taus: + Z = augment(generate_data(time, np.array([tau, tau]), freq_0, theta_0), sgot_window) + B, B_spec = estimator(Z[:-1], Z[1:]) + D = np.log(B_spec["eig_val"]) * fs + L = B_spec["eig_vec_left"] + R = B_spec["eig_vec_right"] + for m in metrics: + decay_scores[m].append( sgot_metric( - D_0, - R_0, - L_0, - D, - R, - L, - eta=0.9, # keep eta fixed here - grassmann_metric=name, + D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric=m ) ) - scores_decay.append(row) -scores_decay = np.array(scores_decay) -plt.figure(figsize=(8, 5)) -for i, name in enumerate(methods): - plt.plot(decays, scores_decay[:, i], label=name) - -plt.xlabel("decay") -plt.ylabel("SGOT distance") -plt.title("SGOT distance vs decay") -plt.legend() +fig, ax = plt.subplots(figsize=(7, 4)) +for m in metrics: + ax.plot(taus, decay_scores[m], styles[m], label=m, linewidth=1.8) +ax.set_xlabel(r"Decay rate $\tau$") +ax.set_ylabel(r"$d_S$") +ax.set_title("SGOT distance vs. decay across Grassmannian metrics") +ax.legend() +fig.tight_layout() plt.show() diff --git a/examples/plot_sgot.py b/examples/plot_sgot.py deleted file mode 100644 index 87f7ff3c8..000000000 --- a/examples/plot_sgot.py +++ /dev/null @@ -1,608 +0,0 @@ -#!/usr/bin/env python -# coding: utf-8 - -""" -========================================= -Spectral-Grassmann OT on dynamical systems operators -========================================= - -This example presents a synthetic example of Spectral Grassmannian-Wasserstein -Optimal Transport (SGOT) on linear dynamical systems. - -We consider a signal formed by the sum of two damped oscillatory modes evolving -along a rotated direction in the plane. The signal is then associated with an -underlying continuous linear dynamical system, and we study how its spectral -representation varies under rotation. The SGOT cost and metric are used to -compare the reference and rotated systems. - -.. [83] T. Germain; R. Flamary; V. R. Kostic; K. Lounici, A Spectral-Grassmann - Wasserstein Metric for Operator Representations of Dynamical Systems, - arXiv preprint arXiv:2509.24920, 2025. -""" - -# Authors: Sienna O'Shea -# Thibaut Germain -# -# License: MIT License - -import numpy as np -import matplotlib.pyplot as plt - -from ot.sgot import sgot_metric, sgot_cost_matrix - -from scipy.linalg import eig - - -# sampling parameters and time grid -fs = 50 -max_t = 5 -time = np.linspace(0, max_t, fs * max_t) -dt = 1 / fs - - -# %% -# Example: rotating a linear dynamical system in 3D -# ------------------------------------------------- -# -# 1. Build a simple observed signal -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# We begin by assuming that the observed signal is made of two oscillatory -# components: -# -# .. math:: -# -# x_{\text{ref}}(t)=e^{-\tau_1 t}\cos(2\pi\omega_1 t)\,\vec e(\theta) -# \;+\; -# e^{-\tau_2 t}\cos(2\pi\omega_2 t)\,\vec e(\theta), -# -# where :math:`\vec e(\theta)\in\mathbb{R}^2` is a fixed real vector. Thus, -# :math:`x(t)` evolves along the one-dimensional subspace spanned by -# :math:`\vec e(\theta)`, while its time dependence exhibits oscillatory and -# dissipative behaviour. - -tau_0 = np.array([0.08, 0.18]) -freq_0 = np.array([1.0, 2.0]) -theta_0 = np.pi / 4 - - -def generate_data(time, tau, freq, theta): - t_ = np.sin(2 * np.pi * freq[None, :] * time[:, None]) * np.exp( - -tau[None, :] * time[:, None] - ) - t_ = t_.sum(axis=1) - traj_0 = np.zeros((t_.shape[0], 2)) - traj_0[:, 0] = t_ - rotation_matrix = np.array( - [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]] - ) - traj_0 = traj_0 @ rotation_matrix.T - return traj_0 - - -traj_0 = generate_data(time, tau_0, freq_0, theta_0) - - -# plot the observed signal components and their sum -plt.figure(figsize=(10, 4)) -plt.plot(time, traj_0, label="base trajectory", linewidth=2) -plt.xlabel("time") -plt.ylabel("amplitude") -plt.legend() -plt.title(r"Observed scalar signal along $\vec{e}(\theta)$") -plt.show() - - -# %% -# 2. Interpret the signal as coming from a continuous linear dynamical system -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# We assume that :math:`x(t)` is generated by an underlying continuous linear -# dynamical system. Since the observed signal is a superposition of two -# sinusoidal modes, the corresponding linear dynamics are naturally described -# by a fourth-order model. We therefore introduce the state vector -# -# .. math:: -# -# z(t)= -# \begin{pmatrix} -# x_1(t)\\ -# x_2(t)\\ -# \vdots\\ -# x_1^{(3)}(t)\\ -# x_2^{(3)}(t) -# \end{pmatrix} -# \in\mathbb{R}^8. -# -# where :math:`x^{(n)}(t)` denotes the n-th derivative of :math:`x(t)`. -# -# This allows us to rewrite the dynamics as a first-order linear system: -# -# .. math:: -# -# \dot{z}(t)=Az(t), -# -# where :math:`A\in\mathbb{R}^{8\times 8}`. Its solution is then given by -# -# .. math:: -# -# z(t)=e^{tA}z_0. - -fig = plt.figure(figsize=(9, 6)) -ax = fig.add_subplot(projection="3d") - -ax.plot(time, traj_0[:, 0], traj_0[:, 1]) -ax.set_xlabel("time") -ax.set_ylabel("x₁(t)") -ax.set_title("Observed trajectory in time") - -ax.text2D(1.08, 0.5, "x₂(t)", transform=ax.transAxes, rotation=90, va="center") - -plt.show() - - -# %% -# 3. Sampling and preprocessing discrete trajectories of the dynamical system -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# We now have a bridge between the continuous system and the operator we later -# aim to infer from sampled data. Since in practice we do not observe the full -# continuous trajectory, we work instead with discrete samples of the signal. -# We take snapshots at uniform time intervals :math:`\Delta t`, and write the -# sampled signal as -# -# .. math:: -# -# S= -# \begin{pmatrix} -# x_1(0) & x_2(0)\\ -# x_1(\Delta t) & x_2(\Delta t)\\ -# \vdots\\ -# x_1(N\Delta t) & x_2(N\Delta t)\\ -# \end{pmatrix} -# -# The goal is now to use these observations to recover the operator governing -# the evolution. To do this, we augment the signal :math:`s` using a sliding -# window of length :math:`w`. For each :math:`k`, define -# -# .. math:: -# -# z_k = -# \begin{pmatrix} -# s_k\\ -# s_{k+1}\\ -# \vdots\\ -# s_{k+w-1} -# \end{pmatrix} -# -# We then form the data matrices -# -# .. math:: -# -# X= -# \begin{pmatrix} -# z_1\\ -# z_2\\ -# \vdots\\ -# z_{N-w} -# \end{pmatrix}, -# \qquad -# Y= -# \begin{pmatrix} -# z_2\\ -# z_3\\ -# \vdots\\ -# z_{N-w+1} -# \end{pmatrix}, -# -# so that :math:`X` contains the present windowed states and :math:`Y` the -# corresponding shifted future states. - - -# build a 4-dimensional state using delay embedding -def augment(traj, window_length=2): - Z = np.lib.stride_tricks.sliding_window_view(traj, (window_length, 1)) - Z = Z.reshape(Z.shape[0], -1) - return Z - - -# create the embedded state matrix Z -Z = augment(traj_0, 4) -Z.shape - -# inspect one embedded state vector -Z[0] - -# create X and Y for the SGOT metric -X = Z[:-1] -Y = Z[1:] - -# inspect shapes of X and Y -print("X shape:", X.shape) -print("Y shape:", Y.shape) - - -# %% -# 4. Estimate the discrete-time operator -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# We now identify the operator that maps :math:`X` to :math:`Y`. From -# -# .. math:: -# -# \dot{z}=Az, -# -# we have -# -# .. math:: -# -# z(t+\Delta t)=e^{\Delta tA}z(t). -# -# Setting -# -# .. math:: -# -# B=e^{\Delta tA}, -# -# the corresponding discrete-time evolution is governed by :math:`B`, and we -# seek the best linear map satisfying -# -# .. math:: -# -# Y\approx X B^T. -# -# Equivalently, we solve the optimisation problem -# -# .. math:: -# -# \min_B \|Y-XB\|^2. -# -# We want to recover the best rank-:math:`r` operator, whose estimator is -# defined as follows: -# -# .. math:: -# -# B = C_{xx}^{-\frac{1}{2}}[C_{xx}^{-\frac{1}{2}}C_{xy}]_r -# \quad \text{s.t} \quad C_{xx} = X^T X \quad \text{and} \quad C_{xy} = X^TY. -# -# Here :math:`[\cdot]_r` denotes the best rank-:math:`r` estimator obtained via -# SVD decomposition. [2] -# -# [2] Kostic, V., Novelli, P., Maurer, A., Ciliberto, C., Rosasco, L. and -# Pontil, M., 2022. Learning dynamical systems via Koopman operator regression -# in reproducing kernel Hilbert spaces. Advances in Neural Information -# Processing Systems, 35, pp.4017-4031. - - -def estimator(X, Y, rank=4, eps=1e-8): - cxx = X.T @ X + eps * np.eye(X.shape[1]) - U, S, Vt = np.linalg.svd(cxx) - S_inv = np.divide(1, S, out=np.zeros_like(S), where=S != 0) - cxx_inv_half = Vt.T @ np.diag(np.sqrt(S_inv)) @ U.T - cxy = X.T @ Y - T = cxx_inv_half @ cxy - U, S, Vt = np.linalg.svd(T) - S[rank:] = 0 - T_rank = U @ np.diag(S) @ Vt - T = cxx_inv_half @ T_rank - - # estimate spectral decomposition - val, vl, vr = eig(T, left=True, right=True) - sort_idx = np.argsort(np.abs(val))[::-1] - val = val[sort_idx][:rank] - vl = vl[:, sort_idx][:, :rank] - vr = vr[:, sort_idx][:, :rank] - - return T, {"eig_val": val, "eig_vec_left": vl, "eig_vec_right": vr} - - -B_0, B_0_spec = estimator(X, Y, rank=4) -Y_pred = X @ B_0 - -plt.figure(figsize=(10, 4)) -plt.plot(Y[:, 0], label="true") -plt.plot(Y_pred[:, 0], "--", label="predicted") -plt.xlabel("sample index") -plt.ylabel("first state coordinate") -plt.title("True Signal vs Predicted Signal") -plt.legend() -plt.show() - - -# %% -# The predicted signal is nearly indistinguishable from the true signal, -# indicating that the estimated operator accurately captures the observed -# dynamics. - -# %% -# 6. Recover continuous-time spectral information from the discrete operator -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# To recover the continuous generator :math:`A`, we study the spectral -# structure of :math:`B`. We diagonalise :math:`B` as -# -# .. math:: -# -# B=PDP^{-1}, -# -# where -# -# .. math:: -# -# D=\operatorname{diag}(\mu_1,\dots,\mu_n). -# -# The continuous-time eigenvalues are of the form -# -# .. math:: -# -# \lambda_k=-\tau_k+2\pi i\,\omega_k, -# \qquad k\in\{1,2\}, -# -# and the corresponding eigenvalues of :math:`B` are -# -# .. math:: -# -# \mu_k=e^{\Delta t\lambda_k} -# =e^{\Delta t(-\tau_k+2\pi i\omega_k)}. -# -# Since :math:`B=e^{\Delta tA}`, we recover :math:`A` by taking the logarithm: -# -# .. math:: -# -# A=P\,\frac{\log(D)}{\Delta t}\,P^{-1}. - -D_0 = np.log(B_0_spec["eig_val"]) * fs -L_0 = B_0_spec["eig_vec_left"] -R_0 = B_0_spec["eig_vec_right"] - -recovered_freqs = D_0.imag / (2 * np.pi) -mask = recovered_freqs > 0 -recovered_freqs = recovered_freqs[mask] -decay = -D_0.real[mask] -print(f"First mode: frequency: {recovered_freqs[0]:.2f} Hz -- decay: {decay[0]:.2f}") -print(f"Second mode: frequency: {recovered_freqs[1]:.2f} Hz -- decay: {decay[1]:.2f}") - - -# %% -# Introduction to SGOT for linear operators -# ----------------------------------------- -# -# To compare two linear operators through their spectral structure, we use the -# SGOT framework introduced in Theorem 1 of [1]. For a non-defective -# finite-rank operator :math:`T \in S_r(\mathcal H)`, the theorem associates a -# discrete spectral measure -# -# .. math:: -# -# \mu(T) \triangleq \sum_{j\in[\ell]} -# \frac{m_j}{m_{\mathrm{tot}}}\,\delta_{(\lambda_j,\mathcal V_j)}, -# -# where :math:`\lambda_j` are the eigenvalues of :math:`T`, :math:`m_j` their -# algebraic multiplicities, and :math:`\mathcal V_j` the corresponding -# eigenspaces. Thus, each spectral component of the operator is represented by -# an atom of the form -# -# .. math:: -# -# (\lambda_j,\mathcal V_j), -# -# combining one eigenvalue with its associated invariant subspace. -# -# Theorem 1 then defines a ground cost between two such atoms by combining a -# spectral discrepancy and a geometric discrepancy: -# -# .. math:: -# -# d_\eta\big((\lambda,\mathcal V),(\lambda',\mathcal V')\big) -# \triangleq -# \eta\,|\lambda-\lambda'| + (1-\eta)\, d_{\mathcal G}(\mathcal V,\mathcal V'), -# -# where :math:`d_{\mathcal G}` denotes the grassmann distance between -# eigenspaces and :math:`\eta\in(0,1)` balances the contribution of eigenvalues -# and eigenspaces. -# -# The SGOT distance between two operators :math:`T` and :math:`T'` is then the -# Wasserstein distance between their associated spectral measures: -# -# .. math:: -# -# d_S(T,T') = W_{d_\eta,p}\big(\mu(T),\mu(T')\big). -# -# In this way, SGOT compares linear operators by optimally matching their -# spectral atoms, taking into account both the location of eigenvalues and the -# relative geometry of their eigenspaces. - -# %% -# A wider delay window for the SGOT experiments below -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# The window of length 4 used above is enough to identify a single reference -# operator, but the experiments below also probe signals whose two modes -# nearly coincide in frequency (e.g. :math:`\omega_2'\to\omega_1`). Telling -# such near-degenerate modes apart requires the delay embedding to span -# enough time to "see" their differing decay, so we re-embed the reference -# signal with a longer window before running the sweeps. - -sgot_window = 10 -Z = augment(traj_0, sgot_window) -_, B_0_spec_sgot = estimator(Z[:-1], Z[1:]) -D_0_sgot = np.log(B_0_spec_sgot["eig_val"]) * fs -L_0_sgot = B_0_spec_sgot["eig_vec_left"] -R_0_sgot = B_0_spec_sgot["eig_vec_right"] - -# %% -# SGOT distance versus rotation angle -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# We compare the reference signal with a rotated version obtained by changing -# only the observation direction. The shifted signal is -# -# .. math:: -# -# x_{\mathrm{shift}}^{\mathrm{rot}}(t;\theta) -# = -# \sum_{i=1}^{2} -# e^{-\tau_i t}\cos(2\pi\omega_i t)\,\vec e({\color{red}\theta}), -# -# while the reference signal is recovered at :math:`\theta=\theta_0`. Thus, -# this experiment isolates the effect of rotating the underlying one-dimensional -# subspace in the observation plane. - -thetas = np.linspace(0, np.pi / 2, 51) -rotation_scores = [] - -for theta in thetas: - Z = augment(generate_data(time, tau_0, freq_0, theta), sgot_window) - B, B_spec = estimator(Z[:-1], Z[1:]) - D = np.log(B_spec["eig_val"]) * fs - L = B_spec["eig_vec_left"] - R = B_spec["eig_vec_right"] - rotation_scores.append( - sgot_metric( - D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric="chordal" - ) - ) - -fig, ax = plt.subplots(figsize=(7, 4)) -ax.plot(thetas, rotation_scores, linewidth=1.8) -ax.axvline(theta_0, color="gray", linestyle="--", linewidth=0.8) -ax.set_xlabel(r"Rotation angle $\theta$ (rad)") -ax.set_ylabel(r"$d_S$") -ax.set_title("SGOT distance vs. rotation angle") -fig.tight_layout() -plt.show() - -# %% -# Comparison across Grassmannian metrics for SGOT distance versus rotation angle -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -metrics = ["chordal", "geodesic", "procrustes", "martin"] -styles = {"chordal": "-", "geodesic": "--", "procrustes": "-.", "martin": ":"} -rotation_results = {m: [] for m in metrics} - -for theta in thetas: - Z = augment(generate_data(time, tau_0, freq_0, theta), sgot_window) - B, B_spec = estimator(Z[:-1], Z[1:]) - D = np.log(B_spec["eig_val"]) * fs - L = B_spec["eig_vec_left"] - R = B_spec["eig_vec_right"] - for m in metrics: - rotation_results[m].append( - sgot_metric( - D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric=m - ) - ) - -fig, ax = plt.subplots(figsize=(7, 4)) -for m in metrics: - ax.plot(thetas, rotation_results[m], styles[m], label=m, linewidth=1.8) -ax.axvline( - theta_0, - color="gray", - linestyle="--", - linewidth=0.8, - label=r"$\theta_0 = \pi/4$ (reference)", -) -ax.set_xlabel(r"Rotation angle $\theta$ (rad)") -ax.set_ylabel(r"$d_S$") -ax.set_title("SGOT distance vs. rotation angle across Grassmannian metrics") -ax.legend() -fig.tight_layout() -plt.show() - -# %% -# SGOT distance versus frequency -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# In this experiment, we keep the reference direction fixed and perturb one of -# the oscillatory modes. The shifted signal is -# -# .. math:: -# -# x_{\mathrm{shift}}^{\omega}(t) -# = -# e^{-\tau_1 t}\cos(2\pi\omega_1 t)\,\vec e(\theta_0) -# \;+\; -# e^{-\tau_2 t}\cos(2\pi{\color{red}\omega_2'} t)\,\vec e(\theta_0), -# -# where only the second frequency is modified. We then study how the SGOT -# distance changes as a function of the perturbed frequency :math:`\omega_2'`. - -omegas = np.linspace(0.5, 3.0, 21) -frequency_scores = {m: [] for m in metrics} - -for omega in omegas: - Z = augment( - generate_data(time, tau_0, np.array([freq_0[0], omega]), theta_0), sgot_window - ) - B, B_spec = estimator(Z[:-1], Z[1:]) - D = np.log(B_spec["eig_val"]) * fs - L = B_spec["eig_vec_left"] - R = B_spec["eig_vec_right"] - for m in metrics: - frequency_scores[m].append( - sgot_metric( - D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric=m - ) - ) - -fig, ax = plt.subplots(figsize=(7, 4)) -for m in metrics: - ax.plot(omegas, frequency_scores[m], styles[m], label=m, linewidth=1.8) -ax.axvline( - freq_0[1], - color="gray", - linestyle="--", - linewidth=0.8, - label=r"$\omega_2 = 2.0$ Hz (reference)", -) -ax.set_xlabel(r"Frequency $\omega_2'$ (Hz)") -ax.set_ylabel(r"$d_S$") -ax.set_title("SGOT distance vs. frequency across Grassmannian metrics") -ax.legend() -fig.tight_layout() -plt.show() - -# %% -# SGOT distance versus decay -# ~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# We now study the effect of changing the decay rate while keeping the -# observation direction fixed. The shifted signal is -# -# .. math:: -# -# x_{\mathrm{shift}}^{\mathrm{decay}}(t;\tau) -# = -# e^{-{\color{red}\tau} t}\cos(2\pi\omega_1 t)\,\vec e(\theta_0) -# \;+\; -# e^{-{\color{red}\tau} t}\cos(2\pi\omega_2' t)\,\vec e(\theta_0). -# -# In this way, both modes share the same modified decay parameter -# :math:`\tau`, allowing us to isolate the influence of dissipation on the SGOT -# distance. -taus = np.linspace(0.1, 3.0, 21) -decay_scores = {m: [] for m in metrics} - -for tau in taus: - Z = augment(generate_data(time, np.array([tau, tau]), freq_0, theta_0), sgot_window) - B, B_spec = estimator(Z[:-1], Z[1:]) - D = np.log(B_spec["eig_val"]) * fs - L = B_spec["eig_vec_left"] - R = B_spec["eig_vec_right"] - for m in metrics: - decay_scores[m].append( - sgot_metric( - D_0_sgot, R_0_sgot, L_0_sgot, D, R, L, eta=0.9, grassmann_metric=m - ) - ) - -fig, ax = plt.subplots(figsize=(7, 4)) -for m in metrics: - ax.plot(taus, decay_scores[m], styles[m], label=m, linewidth=1.8) -ax.set_xlabel(r"Decay rate $\tau$") -ax.set_ylabel(r"$d_S$") -ax.set_title("SGOT distance vs. decay across Grassmannian metrics") -ax.legend() -fig.tight_layout() -plt.show() From 267674ddf64d09abf3b886e640158bf3bc006ff7 Mon Sep 17 00:00:00 2001 From: thibaut-germain Date: Thu, 2 Jul 2026 13:34:37 +0200 Subject: [PATCH 5/5] fix releases.md --- RELEASES.md | 1 + 1 file changed, 1 insertion(+) diff --git a/RELEASES.md b/RELEASES.md index f10d7b62b..eea94981d 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -38,6 +38,7 @@ This new release adds support for sparse cost matrices and a new lazy exact OT s - Update the geomloss wrapper to the new version and API (PR #826) - Fix docstrings for `lowrank_gromov_wasserstein_samples` and `lowrank_sinkhorn` (PR #823) - Reorganize all tests per backend (PR #828) +- Update sgot cost function and example (PR #830) #### Closed issues