Multiscale Proper Orthogonal Decomposition (mPOD)

This document offers a brief reference guide to the use of the Multiscale Proper Orthogonal Decomposition (mPOD).

Theory

The main idea of the mPOD is to perform a POD on non-overlapping portions of the frequency domain of the signal :cite: ninni_modulo_2020. Considering a signal \(\mathbf{u}(t)\) defined in the time domain, sampled with a sampling frequency \(f_s\), one can define to isolate phenomena of interest belonging in different frequency bands identified by :math: M frequency bands. This divides the frequency in blocks, \([0, f_1]\), \([f_1, f_2]\), …, \([f_{M-1}, f_M]\) where \(f_1, f_2, \ldots, f_M\) are the frequencies of interest.

The mPOD performs the MRA on the temporal correlation matrix \(\mathbf{K} = \mathbf{D}^T \mathbf{D} \in \mathbb{R}^{n_t \times n_t}\). Given a list of suitable transfer functions \(\left\{H_m\right\}_{m=1}^M\), the mPOD splits the matrix \(\mathbf{K}\) as:

\[\mathbf{K} = \sum_{m=1}^M \mathbf{K}_m = \sum_{m=1}^M \mathbf{\Psi}_F \big[ \hat{\mathbf{K}} \odot H_m \big] \mathbf{\Psi}_F\]

where \(\mathbf{\Psi}_F\) is the Fourier transform matrix, \(\hat{\mathbf{K}}=\bar{\mathbf{\Psi}_F \mathbf{K} \mathbf{\Psi}\) is the 2D Fourier transform of the correlation matrix, and \(\odot\) is the entry by entry product.

The filtering ensures that the key properties of \(\mathbf{K}\) are preserved, and thus characterised with a set of orthonormal eigenvectors and eigenvalues:

\[\mathbf{K}_m = \sum_{j=1}^{n_m} \lambda_{m,j} \Psi_{P, mj} \Psi_{P, mj}^T\]

and \(n_m\) is the number of nonzero eigenvalues at each scale. This procedure also ensures that the eigenvectors of all scales are orthogonal complements that span \(\mathbb{R}^{n_t}\).

Finally, the mPOD temporal basis is constructed by collecting the POD basis of all scales, such that:

\[\mathbf{\Psi}_M = [\mathbf{\Psi}_{1}, \mathbf{\Psi}_{2}, \ldots, \mathbf{\Psi}_{M}] \mathbf{P}_\mathbf{\Sigma}\]

where \(\mathbf{P}_\mathbf{\Sigma}\) is a permutation matrix that ranks the structures in decreasing order of energy contribution.

In conclusion, we note that the mPOD generalizes POD and DFT: for \(M=1\), the mPOD is equivalent to POD, while for \(M=n_t\), the mPOD is equivalent to DFT.

fast mPOD

The fast mPOD changes the workflow of classical mPOD by performing the analysis in frequency domain and replacing filters with compact-support masks [7]. This allows to reduce the computational cost of the mPOD by orders of magnitude, while preserving the key properties of the decomposition. Since it uses spectral masks instead of filter bank, it sacrifices the advanced boundary treatment of the classical mPOD to gain the aforementioned computational benefits.

Example in MODULO

An example of the usage of the mPOD routine, extracted from the examples folder is reported below.

import numpy as np
from modulo_vki.modulo import ModuloVKI
from modulo_vki.utils import plot_mPOD

FOLDER_MPOD_RESULTS=FOLDER+os.sep+'mPOD_Results_Jet'
if not os.path.exists(FOLDER_MPOD_RESULTS):
    os.mkdir(FOLDER_MPOD_RESULTS)

# We here perform the mPOD as done in the previous tutorials.
# This is mostly a copy paste from those, but we include it for completenetss

# Updated in MODULO 2.1.0: mPOD requires now Keep and Nf of size len(F_V) + 1, to ensure all scales are computed. First entry specifies the low pass filter.

Keep = np.array([1, 1, 1, 1, 1])
Nf = np.array([201, 201, 201, 201, 201])
# --- Test Case Data:
# + Stand off distance nozzle to plate
H = 4 / 1000
# + Mean velocity of the jet at the outlet
U0 = 6.5
# + Input frequency splitting vector in dimensionless form (Strohual Number)
ST_V = np.array([0.1, 0.2, 0.25, 0.4])
# + Frequency Splitting Vector in Hz
F_V = ST_V * U0 / H
# + Size of the extension for the BC (Check Docs)
Ex = 203  # This must be at least as Nf.
dt = 1/2000; boundaries = 'reflective'; MODE = 'reduced'
# Here 's the mPOD
Phi_M, Psi_M, Sigmas_M = m.mPOD(Nf, Ex, F_V, Keep, 20 ,boundaries, MODE, dt, False)

The variable Keep is a vector of size len(F_V) + 1, defines whether the scale is processed or not, while Nf is a vector of size len(F_V) + 1 is the number of points in the frequency domain, and Ex is the size of the extension for the BC. The boundary conditions can be set to reflect, nearest, wrap or extract.

Fast mPOD

Fast mPOD is fully functional without Fortran. The default useFortran=False path is portable and is the supported configuration on Windows. See Installation before opting into the Linux/macOS Fortran accelerator.

An example of the usage of the fast mPOD routine, extracted from the examples folder is reported below.

import numpy as np
from modulo_vki.modulo import ModuloVKI
import os
from math import pi

if not os.path.exists("./ex6_results"):
    os.makedirs("./ex6_results")

# Create artificial spatial structure
toposes = np.zeros((3000,3))
toposes[:1000,0] = 1
toposes[1000:2000,1] = 1
toposes[2000:,2] = 1
toposes = toposes/np.linalg.norm(toposes,axis=0)

# Associate importance with the test modes
sigmas = np.diag([6000,4000,2000])

# Create artificial temporal structure
nt = 6000   # number of time steps
chronoses = np.zeros((nt,3))
for i in range(2000):
    chronoses[i,0] = np.sin(2*pi*i/60)+(i-500)/1200     # f1/fs = 1/60 + slope
for i in range(4000,6000):
    chronoses[i,1] = np.sin(2*pi*(i-4000)/30)-(i-4000)/800   # f2/fs = 1/30 + slope

# single wavelet
for i in range(2200,3000):
    chronoses[i,2] = 1
for i in range(3000,3800):
    chronoses[i,2] = -1

for i in range(3):
    chronoses[:,i] = chronoses[:,i]/np.linalg.norm(chronoses[:,i])

# Assemble the data matrix
D = toposes @ sigmas @ chronoses.T

# Save temporal evolutions
chronTot = np.zeros((chronoses.shape[0],chronoses.shape[1]+1))
chronTot[:,0] = np.arange(nt)
chronTot[:,1:] = chronoses
np.save(f"{target_dir}/chronoses.npy", chronTot)

# Frequency split and selection of scales
nModes = 10     # Number of modes to compute
F_V = [1/120,3/120,5/120,10/120]   # Frequency splitting vector in Hz, assuming fs = 1 (normalized frequencies) for this tutorial
Keep = [1, 1, 1, 1, 0]  # Which scales to keep (0 means not processed, 1 means processed)
winType = "hann" # Type of window for the tapering (hann, hamming, blackman, etc)
taper = [1/300, 1/300, 1/300, 1/300, 1/300] # Taper width in Hz for each scale, assuming fs = 1 (normalized frequencies) for this tutorial
mode = "fullK" # Type of computation (supports "fullK", "bandK", "fullSVD", "randSVD")


# MODULO object initialization when NOT using Fortran.
m = ModuloVKI(data=np.nan_to_num(D), n_Modes = nModes)

# Fast mPOD call
phi, psi, sig = m.fastmPOD(F_V, 1, Keep, winType = winType, taper = taper, mode = mode, GThresh=5, ncpus=9)

An example of the usage of the fast mPOD routine with Fortran support for correlation matrix assembly. Only the parts that differ from the previous example are shown below.

# MODULO object initialization when USING Fortran.
m = ModuloVKI(data=np.nan_to_num(D), n_Modes = nModes, dtype=np.float64) # dtype = np.float64 is MANDATORY for Fortran support.

# Fast mPOD call
phi, psi, sig = m.fastmPOD(F_V, fs, Keep, winType = winType, taper = taper, mode = mode, GThresh=5, ncpus=9, useFortran=True)