Spectral dynamic mode decomposition
H. Lange et al.
Spectral dynamic mode decomposition is my term for Algorithm 1 from the paper From Fourier to Koopman: Spectral Methods for Long-term Time Series Prediction by H. Lange, S.L. Brunton, J.N. Kutz (video) (arXiv). Necessary background is a familiarity with the dynamic mode decomposition (video) (Wikipedia entry).
This project walks through a toy example I coded up in Python. The example is similar to one from the paper which gets across the main points. I have packed away the main functions for the algorithm in the utility spectral_help.py which has been linked in this project’s description.
import numpy as np
from numpy.linalg import norm, solve
import matplotlib.pyplot as plt
cmap = plt.get_cmap('tab20')
from spectral_help import *
Example Definition
class toy_data:
"""
Generate toy oscillation data at specific frequencies.
"""
def __init__(self, tf, npts, noise):
"""
Parameters:
tf: final time
npts: number of samples between 0 and tf
noise: standard deviation of additive Gaussian noise
"""
self.tf = tf
self.npts = npts
self.ts = np.linspace(0, self.tf, self.npts)
self.dt = self.ts[1] - self.ts[0]
self.noise = noise
# Manual toy data (a stack of sines)
self.freqs = [0.5, 2, 0.75, 3]
self.X_fn = lambda ts: np.vstack([(np.sin([2 * np.pi * self.freqs[0] * ts]) + np.sin([2 * np.pi * self.freqs[1] * ts])),
(np.sin([2 * np.pi * self.freqs[2] * ts]) + np.sin([2 * np.pi * self.freqs[3] * ts]))])
# Normalize, add noise
X = self.X_fn(self.ts)
self.X_mean = np.mean(X, axis=1).reshape(-1,1)
self.X_std = np.std(X, axis=1).reshape(-1,1)
X = (X - self.X_mean)/self.X_std
self.X_true = np.copy(X)
self.X = X + np.random.randn(*X.shape)*self.noise
# Construct test vars
self.ts_test = None
self.X_test = None
def run_test(self, tf_predict, npts_predict):
"""
Parameters:
tf_predict: final time
npts_predict: number of points between 0 and tf_predict
Updates:
self.ts_test: test time series
self.X_test_true: normalized ground truth for the test
"""
self.ts_test = np.linspace(0, tf_predict, npts_predict)
X_test = self.X_fn(self.ts_test)
# Use training std_dev and mean
self.X_test_true = (X_test - self.X_mean)/self.X_std
Run the Experiment
We run the spectral dynamic mode decomposition algorithm to learn the frequencies of an operator generating our toy time-series data. Because of the nature of the algorithm, we can do this with very noisy data (Here, we set the standard deviation at 0.5 for mean-zero variance-one toy training data).
First, we configure the model. Then we set up the algorithm and run the optimization. The optimization is performed using the Accelerated Proximal Gradient Descent Method, or AGPD. The reason is because we are imposing sparsity using a $\ell_1$ norm–that is, a LASSO-type optimization.
1. Configure the model
# Config
# ======
np.random.seed(1)
# Data params
npts = 400 # number of time points
tf = 8 # max time
std_noise = .5 # set the amount of noise on the mean-0, variance-1 data
predict_factor = 3 # prediction goes out to factor*T
exper = toy_data(tf, npts, std_noise)
# Algorithm params
freq_dim = 24 # freq. for algo to try
learning_rate = 1e-3 # LR = 1/beta from beta-smooth obj. bound (at least ideally--I'm just choosing a number here)
reg_factor = 5 # regularization on sparsity
# SVD parameters
threshold_type = 'percent' # choose 'count' for number or 'percent' for s/max(s) > threshold
threshold = 1e-1
# Plot toggle
print_omega_updates = True
def print_update(omg, title):
if print_omega_updates:
print('{} $\omega$:\t'.format(title), np.sort(np.round(omg[omg.astype(bool)], 3)))
2. Set up the algorithm & 3. Optimize
# Algorithm
# =========
# 1. Initialize
omega = np.zeros(freq_dim)*2
print_update(omega, 'Initial')
A = np.random.rand(exper.X.shape[0], freq_dim*2)
obj_his = []
err_his = []
# 2. FFT to obtain the initial starting point for the optimization
for ifreq in range(len(omega)):
# - Construct the residual via the current frequencies
res = residual_j(ifreq, exper.X, A, omega, exper.ts)
# - Select the maximum fft frequency as the initial value
omega[ifreq] = max_fft_update(res, exper.dt)
# - Update A
A = update_A(exper.X, omega, exper.ts, threshold, threshold_type)
# 3. Perform proximal gradient descent from the initial point
# - Construct optimization functions
lam_cs = reg_factor*norm(A.T.dot(exper.X), np.inf)
def f(w):
return loss(exper.X, A, w.flatten(), exper.ts)
def gradf(w):
return grad_loss(exper.X, A, w.flatten(), exper.ts)
def func_g(w):
return lam_cs*np.linalg.norm(w, ord=1)
def prox_g(w, t):
res = []
r = t*lam_cs
for wi in w.flatten():
if wi > r:
res.append(wi-r)
elif wi < -r:
res.append(wi+r)
else:
res.append(0)
return np.array(res)
# - Optimization algorithm
w, iobj_his, ierr_his, cond = optimizeWithAPGD(omega.reshape(1,-1), f, func_g, gradf, prox_g, (1/learning_rate)*npts, max_iter=5000, verbose=True)
obj_his.append(iobj_his)
err_his.append(ierr_his)
omega = w
print_update(omega, 'Final ')
# 4. Final operator update
A = update_A(exper.X, omega, exper.ts, threshold, threshold_type)
print_update(np.array(exper.freqs), 'Expected')
Initial $\omega$: []
Final $\omega$: [0.494 0.743 1.991 2.99 ]
Expected $\omega$: [0.5 0.75 2. 3. ]
We printed out the frequencies (previous cell). In the next cell, we show the objective value and the gradient of the objective for the iterations of the optimization. The characteristic oscillations of an accelerated gradient descent are observed.
# Plot
# ====
# Inspect convergence results
fig,axes = plt.subplots(2,1,figsize=[10,10])
ax = axes[0]
ax.plot(iobj_his)
ax.set_ylabel('Obj. value')
ax.set_xlabel('Iterations')
ax.set_yscale('log')
ax = axes[1]
ax.plot(ierr_his)
ax.set_ylabel('Gradient magn.')
ax.set_xlabel('Iterations')
ax.set_yscale('log')
Result
We see the training (top) and test (bottom) results for our two-dimensional multi-frequency dynamics.
Notice the order of magnitude increase in the horizontal time axis on the bottom plot. On this plot, we observe that the solutions match the test simulations and are stable because of the model assumptions.
# Make prediction
exper.run_test(tf*predict_factor, 10*npts) # keep sample freq the same
bigOmg_test = BigOmg(omega, exper.ts_test)
X_pred = A@(bigOmg_test)
# Plot data
fig,axes = plt.subplots(1,2,figsize=[20,3])
fig.subplots_adjust(hspace=0.3, wspace=0.2)
leg_params = {'loc': 'upper right', 'shadow': True, 'fancybox': True}
ax = axes[0]
ax.plot(exper.ts, exper.X_true[0], color=cmap(1), label='Simulation')
ax.plot(exper.ts, exper.X[0], ls='', marker='+', color=cmap(0), label='Training Data')
ax.set_xlabel('t')
ax.set_ylabel('x')
ax.legend(**leg_params)
ax = axes[1]
ax.plot(exper.ts, exper.X_true[1], color=cmap(3), label='Simulation')
ax.plot(exper.ts, exper.X[1], ls='', marker='+', color=cmap(2), label='Training Data')
ax.legend(**leg_params)
ax.set_xlabel('t')
ax.set_ylabel('y')
ax.set_ylim([-3.5,3.5])
# Plot model
fig,axes = plt.subplots(2,1,figsize=[20,5])
fig.subplots_adjust(hspace=0.3, wspace=0.2)
leg_params = {'loc': 'upper right', 'shadow': True, 'fancybox': True}
ax = axes[0]
ax.plot(exper.ts_test, exper.X_test_true[0], color=cmap(1), label='Test Simulation')
ax.plot(exper.ts_test, X_pred[0], ls='-', color=cmap(0), label='Model Prediction')
ax.legend(**leg_params)
ax.set_xlabel('t')
ax.set_ylabel('x')
ax = axes[1]
ax.plot(exper.ts_test, exper.X_test_true[1], color=cmap(3), label='Test Simulation')
ax.plot(exper.ts_test, X_pred[1], ls='-', color=cmap(2), label='Model Prediction')
ax.legend(**leg_params)
ax.set_xlabel('t')
ax.set_ylabel('y')
ax.set_ylim([-3.5,3.5])
(Click on the figures to zoom.)
Training data:
Test data:
