dmd.py

# Author : Jean-Christophe Loiseau <[email protected]>
# Date : July 2020

# --> Standard python libraries.
import numpy as np
from scipy.linalg import pinv, eigh, eig

def dmd_analysis(x, y=None, rank=2):
  
  # --> Partition data matrix into X and Y.
  if y is None:
    x, y = x[:, :-1], x[:, 1:]
    
  # --> Pre-compute the covariance and cross-covariance matrices.
  Cxx, Cyx = x @ x.T, y @ x.T
  pinv_Cxx = pinv(Cxx)
  
  # --> Compute the eigenspectrum of the symmetric positive definite matrix.
  #     appearing in the closed-form solution of the DMD problem.
  #     Inspecting its eigenspectrum will help choosing an appropriate rank for the model.
  sigma, P = eigh(np.linalg.multi_dot([Cyx, pinv_Cxx, Cyx.T]))
  
  # --> Sort and normalize the eigenvalues.
  idx = np.argsort(-sigma)
  sigma, P = sigma[idx], P[:, idx]
  sigma /= np.sum(sigma)
  
  # --> Truncate the output basis to the desired rank.
  P = P[:, :rank].reshape(-1, rank)
  
  # --> Compute the basis for the input space.
  Q = (Cyx @ pinv_Cxx).T @ P
  
  # --> Compute the eigenpairs of the low-dimensional DMD operator.
  mu, psi, phi = eig(Q.T @ P, left=True, right=True)
  
  # --> Compute the high-dimensional DMD modes.
  phi, psi = P @ phi, Q @ psi
  
  # --> Normalize using their bi-orthogonal property.
  M = psi.T.conj() @ phi
  
  if rank > 1:
    psi /= np.diag(M).conj()
   else:
    psi /= M.conj()
    
  return P, Q, sigma, mu