EDMD

For people already familiar with the method and interface, the API docs.

Similarly to DMD, extended DMD (EDMD) [1] is a data-driven method to approximate the Koopman operator. The Koopman operator $\mathcal{K}$ describes how an observable of an (potentially stochastic) dynamical system evolves in time. In a discrete-time setting with $x_{t+1} = T(x_t)$ and an observable $g$, the Koopman operator is defined as

$$(\mathcal{K} g)(x) = g(T(x)).$$

Given observation data $X$ and $Y$, where $y_i = T(x_i)$ and a basis set of functions or observables

$$B = \{ \psi_1,\psi_2,\ldots,\psi_k \},$$

one can define the transformed data $\Psi_X$ and $\Psi_Y$ with

$$(\Psi_*)_i = (\psi_1(*_i),\ldots, \psi_k(*_i))^\top,$$

yielding

$$\Psi_Y^\top = \Psi_X^\top K.$$

Here, $K^\top\in\mathbb{R}^{k\times k}$ acting from the left is the projection of the true Koopman operator $\mathcal{K}$ into the basis $B$.

Usually one performs rank-truncation or has over-/underdetermined systems. In such cases, $K$ is determined via the pseudo-inverse $K=\Psi_Y\Psi_X^\dagger$.

In case of a basis $B$ which just contains the identity, EDMD reduces back to DMD.

To demonstrate EDMD and deeptime’s API, we consider a two-dimensional triple-well potential [2]. It is defined on $[-2,2]\times [-1, 2]$ and looks like the following:

import numpy as np
import matplotlib.pyplot as plt

from deeptime.data import triple_well_2d
system = triple_well_2d()

x = np.linspace(-2, 2, num=100)
y = np.linspace(-1, 2, num=100)
XX, YY = np.meshgrid(x, y)
coords = np.dstack((XX, YY)).reshape(-1, 2)
V = system.potential(coords).reshape(XX.shape)

plt.contourf(x, y, V, levels=np.linspace(-4.5, 4.5, 20), cmap='coolwarm');

Inside this potential landscape we integrate particles subject to the stochastic differential equation

$$\mathrm{d}X_t = \nabla V(X_t) \mathrm{d}t + \sigma(t, X_t)\mathrm{d}W_t$$

with $\sigma=1.09$, using an Euler-Maruyama integrator.

traj = system.trajectory(x0=[[-1, 0]], length=200, seed=66)

plt.contourf(x, y, V, levels=np.linspace(-4.5, 4.5, 20), cmap='coolwarm')
plt.plot(*traj.T, c='black', lw=.5);

To approximate the Koopman operator, we first select 25000 random points inside the domain and integrate them for 10000 steps under an integration step of $h=10^{-5}$:

N = 25000
state = np.random.RandomState(seed=42)
X = np.stack([state.uniform(-2, 2, size=(N,)), state.uniform(-1, 2, size=(N,))]).T
Y = system(X, n_jobs=8)

Now we pick a basis $B$ to be all monomials up to and including degree $\mathrm{deg}=10$.

from deeptime.basis import Monomials

basis = Monomials(p=10, d=2)

Using this basis, we can fit an EDMD estimator and obtain the corresponding model.

from deeptime.decomposition import EDMD

edmd_estimator = EDMD(basis, n_eigs=8)
edmd_model = edmd_estimator.fit((X, Y)).fetch_model()

We can obtain the dominant eigenvalues associated to processes in the system.

plt.plot(np.real(edmd_model.eigenvalues), 'x')
plt.title('Dominant eigenvalues');

We project our input data $X$ using the eigenfunctions:

phi = np.real(edmd_model.transform(X, propagate=False))

# normalize
for i in range(len(edmd_model.eigenvalues)):
    phi[:, i] = phi[:, i] / np.max(np.abs(phi[:, i]))

And visualize the first four eigenfunctions. They contain information about metastable regions and slow transitions.

fig = plt.figure(figsize=(8, 8))
gs = fig.add_gridspec(ncols=2, nrows=2)

ax = fig.add_subplot(gs[0, 0])
ax.scatter(*X.T, c=phi[:, 0], cmap='coolwarm')
ax.set_title('1st eigenfunction')

ax = fig.add_subplot(gs[0, 1])
ax.scatter(*X.T, c=phi[:, 1], cmap='coolwarm')
ax.set_title('2nd eigenfunction')

ax = fig.add_subplot(gs[1, 0])
ax.scatter(*X.T, c=phi[:, 2], cmap='coolwarm')
ax.set_title('3rd eigenfunction')

ax = fig.add_subplot(gs[1, 1])
ax.scatter(*X.T, c=phi[:, 3], cmap='coolwarm')
ax.set_title('4th eigenfunction');

Using a clustering can reveal temporally coherent structures.

from deeptime.clustering import KMeans

kmeans = KMeans(n_clusters=6, n_jobs=4).fit(np.real(phi)).fetch_model()
c = kmeans.transform(np.real(phi))

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
ax1.scatter(*X.T, c=c)
ax1.set_title(r"$t = 0$")
ax1.set_aspect('equal')
ax1.set_xlim([-2.5, 2.5])
ax1.set_ylim([-1.5, 2.5])
ax1.contour(x, y, V, levels=np.linspace(-4.5, 4.5, 20), colors='black')

ax2.scatter(*Y.T, c=c)
ax2.set_title(r"$t = 0.1$")
ax2.set_aspect('equal')
ax2.set_xlim([-2.5, 2.5])
ax2.set_ylim([-1.5, 2.5])
ax2.contour(x, y, V, levels=np.linspace(-4.5, 4.5, 20), colors='black');

References

[1]

Matthew O Williams, Ioannis G Kevrekidis, and Clarence W Rowley. A data-driven approximation of the Koopman operator: extending dynamic mode decomposition. Journal of Nonlinear Science, 25(6):1307–1346, 2015.

[2]

Christof Schütte and Marco Sarich. Metastability and Markov State Models in Molecular Dynamics. Volume 24. American Mathematical Soc., 2013.