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.

# 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:

\[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 \(\vec e(\theta)\in\mathbb{R}^2\) is a fixed real vector. Thus, \(x(t)\) evolves along the one-dimensional subspace spanned by \(\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 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]
    )
    t_ = t_.sum(axis=1)
    traj_0 = np.zeros((t_.shape[0], 2))
    traj_0[:, 0] = 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_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()