Delay-Embedding Spatio-Temporal Dynamic Mode Decomposition

Abstract

Spatio-temporal dynamic mode decomposition (STDMD) is an extension of dynamic mode decomposition (DMD) designed to handle spatio-temporal datasets. It extends the framework so that it can analyze data that have both spatial and temporal variations. This facilitates the extraction of spatial structures along with their temporal evolution. The STDMD method extracts temporal and spatial development information simultaneously, including wavenumber, frequencies, and growth rates, which are essential in complex dynamic systems. We provide a comprehensive mathematical framework for sequential and parallel STDMD approaches. To increase the range of applications of the presented techniques, we also introduce a generalization of delay coordinates. The extension, labeled delay-embedding STDMD allows the use of delayed data, which can be both time-delayed and space-delayed. An explicit expression of the presented algorithms in matrix form is also provided, making theoretical analysis easier and providing a solid foundation for further research and development. The novel approach is demonstrated using some illustrative model dynamics.

1. Introduction

Dynamical systems are prevalent in science and engineering, yet analyzing and predicting them remains challenging. While linear systems are well characterized, nonlinear systems are difficult to characterize. They can exhibit an extremely wide range of behaviors, including chaos, and generally do not yield analytical solutions. Koopman operator theory plays an important role in the analysis of such systems [1,2]. The idea is based on transforming the finite-dimensional dynamics of the nonlinear state space into an infinite-dimensional linear dynamical system of functions on the state, represented by the Koopman operator. Through the eigendecomposition of the Koopman operator, we can understand the behavior, stability and long-term dynamics of complex systems. One of the leading algorithms for Koopman spectral analysis is dynamic mode decomposition (DMD), introduced by Schmid in [3]. The method comprises a mathematical technique for identifying spatio-temporal coherent structures from high-dimensional data. After its introduction, the method is now used in a variety of fields, including various jets  [4,5], epidemiology [6], video processing [7], neuroscience  [8], financial trading [9,10,11], robotics [12] and cavity flows [13,14]. For a review of the DMD literature, we refer the reader to [15,16,17,18,19]. For some recent results on DMD extensions, we recommend [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42] to the reader.

While standard DMD is a powerful technique for analyzing dynamic systems, it has limitations related to its assumptions, sensitivity to noise, ability to capture long-term dynamics, computational complexity, parameter sensitivity and others. Researchers continue to develop and refine variations in the DMD method to address these shortcomings and improve its applicability to a wide range of data analysis tasks. Over the last few years, several variants of DMD have been proposed. Chen [20] proposed an optimized DMD method that can reduce numerical sensitivity and calculate the modal growth rate and frequency accurately. Williams et al. [21] suggested extended DMD (EDMD), which can produce improved approximations of the leading Koopman eigenfunctions and eigenvalues. Moreover, Le Clainche et al. [22] developed higher-order DMD (HODMD), which extends DMD to resolve delayed snapshots. In [23], Le Clainche and Vega introduce spatio-temporal Koopman decomposition (STKD), which incorporates higher order DMD and a spatio-temporal approach for the Koopman operator.

One of the modifications of the DMD method, which will play a key role in the exposition of the present work, is the delay-embedding DMD (or Hankel DMD) [43,44]. Delay-embedding methods have also been employed for system identification, most notably by the eigensystem realization algorithm (ERA)[45] and in climate science with singular spectrum analysis (SSA) [46]. Brunton et al. [47] developed a variant of this technique called the Hankel alternative view of Koopman (HAVOK) analysis.

In the present work, we consider the spatio-temporal DMD (STDMD), a generalization of the DMD method designed to handle spatio-temporal datasets. It extends the framework so that it can analyze data that have both spatial and temporal variations, by extracting spatial structures and their temporal evolution. The STDMD method extracts temporal and spatial development information simultaneously, including wavenumber and spatial growth rate. This can be crucial in complex dynamic systems. The “spatio-temporal” aspect refers to the fact that DMD is applied to data that vary both in space and time, such as sequences of images or sensor measurements collected over time and across multiple spatial locations. In such data, patterns and structures can evolve both spatially and temporally, and the approach aims to capture these spatio-temporal dynamics. Applications of spatio-temporal DMD span various fields, including fluid dynamics, neuroscience, climate science, and engineering, where understanding and predicting complex spatio-temporal behaviors is essential. Some recent publications related to the topic suggest applications in the fields of unsteady shear layer flow [48], wake of a circular cylinder [49], urban flow [50], aerodynamic modeling [51], turbulent flow [52] and binary fluid convection [53].

We provide a comprehensive mathematical framework for sequential and parallel STDMD approaches. A clear expression of the presented algorithms in matrix form is also provided. This facilitates theoretical analysis and provides a solid foundation for further research and development. Furthermore, we introduce a delay coordinate generalization of STDMD, enabling the use of both time-delayed and space-delayed snapshots. This extension, labeled delay-embedding STDMD, can be considered as an alternative approach to the STKD method proposed in [23]. The proposed STDMD approach is compared with the results obtained from STKD.

The following is an outline of the paper: in the rest of Section 1, we describe the DMD and some basic concepts related to it; spatio-tempral DMD approaches are in Section 2; in Section 3, we introduce and discuss the framework for delay-embedding STDM; in Section 4, we present the numerical results; in Section 5, we provide the conclusion.

1.1. Dynamic Mode Decomposition

In this paragraph, a brief introduction to the classical dynamic mode decomposition (DMD) framework is provided. For details, we refer the reader to [16,17,19] and the references therein. Consider the system of time-invariant ordinary differential equations of the form

$$\dot{\mathbf{x}}\left(t\right)=f\left(\mathbf{x}\left(t\right)\right),$$

(1)

where $\mathbf{x}\in {\mathcal{R}}^n$ is the state vector and $f:{\mathcal{R}}^n\to {\mathcal{R}}^n$ is a nonlinear map ($n\gg 1$). Let the discrete-time representation of (1) be

$${\mathbf{x}}_{k+1}=\mathbf{F}\left({\mathbf{x}}_k\right),$$

(2)

where ${\mathbf{x}}_k\in {\mathcal{R}}^n$ is a high-dimensional state vector sampled at $t_k=k▵t$ for $k=0,\dots ,m$, and $\mathbf{F}$ is an unknown map that describes the evolution of the state vector between two subsequent sampling times. The initial condition is defined by $\mathbf{x}\left(0\right)={\mathbf{x}}_0$.

Suppose that the evolution of the high-dimensional state $\mathbf{x}$ is governed by some underlying low-dimensional dynamics. Then, the DMD computes a data-driven linear approximation to the system (2) as follows: the sequential set of data

$$\mathcal{D}=[{\mathbf{x}}_0,\dots ,{\mathbf{x}}_m]$$

(3)

is arranged into the following two large data matrices

$$X=[{\mathbf{x}}_0,\dots ,{\mathbf{x}}_{m-1}]\mathrm{and}Y=[{\mathbf{x}}_1,\dots ,{\mathbf{x}}_m].$$

(4)

The goal of the DMD approach is to find a relationship between the future state ${\mathbf{x}}_{k+1}$ and the current state ${\mathbf{x}}_k$, given by

$${\mathbf{x}}_{k+1}=A{\mathbf{x}}_k,$$

(5)

where $A\in {\mathcal{R}}^{n\times n}$ is called the DMD operator. The solution of (5) may be expressed simply in terms of the eigenvalues ${\lambda }_j$ and eigenvectors ${\varphi }_j$ of A:

$${\mathbf{x}}_k=\sum _{j=1}^r{\varphi }_j{b}_j{\lambda }_j^k=\Phi {\Lambda }^k\mathbf{b},$$

(6)

where $\Phi$ is the eigenvector matrix of A, $\Lambda$ is the diagonal matrix of eigenvalues $\Lambda =\mathrm{diag}\left\{{\lambda }_i\right\}$, $\mathbf{b}={\Phi }^{†}{\mathbf{x}}_0$, and ${\Phi }^{†}$ is the Moore–Penrose pseudoinverse of $\Phi$. The parameter r is determined by the low-rank eigendecomposition of matrix A.

Therefore, the corresponding continuous-time approximation of the system (1) can be written as

$$\dot{\mathbf{x}}=\mathcal{A}\mathbf{x},\mathrm{with}A=exp\left(\mathcal{A}\right)$$

(7)

and the initial condition $\mathbf{x}\left(0\right)$. Then, the state-variable evolution in time can be approximated by the following modal expansion

$$\mathbf{x}\left(t\right)=\sum _{j=1}^r{\varphi }_j{b}_{j}exp\left({\omega }_{j}t\right)=\Phi exp(\Omega t)\mathbf{b},$$

(8)

where ${\varphi }_j$ are also the eigenvectors of the approximated matrix $\mathcal{A}$ and matrix $\Omega =\mathrm{diag}\left({\omega }_j\right)$ is a diagonal matrix whose entries are

$${\omega }_j=ln\left({\lambda }_j\right)/▵t$$

(9)

the eigenvalues of $\mathcal{A}$, with ${\lambda }_j$ the eigenvalues of A. The real part of ${\omega }_j$ regulates the growth or decay of the DMD modes, while the imaginary part of ${\omega }_j$ drives oscillations in the DMD modes. In this sense, while the discrete-time eigenvalues ${\lambda }_i$ imply stability when they are inside the unit disc in $\mathbf{C}$, the continuous-time eigenvalues ${\omega }_i$ imply stability when they are in the left half-plane of $\mathbf{C}$. Each component $b_j$ of vector $\mathbf{b}$, in (6) and (8), is a complex scalar that represents the i-th modal contribution of initial vector ${\mathbf{x}}_0$ and can be interpreted as the amplitude of the corresponding DMD mode ${\varphi }_j$.

1.2. Reduced-Order DMD Operator

The relation (5) can be rewritten in terms of snapshot matrices

$$Y=AX.$$

(10)

Then, the dynamic mode decomposition of data matrix $\mathcal{D}$ is given by the eigendecomposition of A. The DMD finds the best-fit solution A, one that minimizes the least-squares distance in the Frobenius norm

$$arg\underset{A}{min}{\parallel Y-AX\parallel }_F,$$

(11)

where ${\parallel .\parallel }_F$ is the Frobenius norm. The solution A to this optimization problem is given by

$$A\approx Y{X}^{†},$$

(12)

where $X^{†}$ denotes the Moore–Penrose pseudo-inverse of X. This is the same as saying that A minimizes $\parallel {\mathbf{x}}_{k+1}-A{\mathbf{x}}_k{\parallel }_2$ across all time steps. The DMD modes and eigenvalues are intended to approximate the eigenvectors and eigenvalues of A.

In practice, the A matrix can be too large and it is computationally inefficient to explicitly compute $A\approx Y{X}^{†}$. It should be noted that calculating the eigendecomposition of the $n\times n$ matrix A can be prohibitively expensive if n is large, i.e., $n\gg 1$. In such cases, DMD aims at finding a reduced representation of A by $\tilde{A}\in {\mathcal{R}}^{r\times r}$ with $r\ll n$. Matrix $\tilde{A}$ can be used to construct DMD modes associated with specific temporal frequencies. Thus, we can use the dynamics of low-rank approximation to represent the full state dynamics. This basis transformation takes the form

$$\mathbf{x}=Q\tilde{\mathbf{x}},$$

(13)

where Q is usually a unitary matrix or such that $Q^{*}Q=I$. The reduced-order model, corresponding to (5), can be derived as follows:

$${\tilde{\mathbf{x}}}_{k+1}=A{\tilde{\mathbf{x}}}_k,$$

(14)

where the corresponding reduced-order matrix is

$$\tilde{A}=Q^{*}AQ,$$

(15)

such that $\tilde{A}\in {\mathcal{R}}^{r\times r}$. The eigenvalues of $\tilde{A}$ and A are equivalent, because of similarity transformation and the eigenvectors are related via a linear transformation.

Let the eigendecomposition of $\tilde{A}$ be

$$\tilde{A}W=W\Lambda$$

(16)

where W is the eigenvector matrix and $\Lambda$ is the diagonal matrix of the associated eigenvalues $\Lambda =\mathrm{diag}\left\{{\lambda }_i\right\}$. Then, the matrix of DMD modes is

$$\Phi =QW$$

(17)

which approximates the eigenvector matrix of A.

Some possible choices for the projection matrix Q in (13) are:

(i). The left singular vector matrix of X. A common approach to choosing the transformation matrix Q is

$$Q=U,$$

(18)

from the truncated SVD of X:

$$X=U\Sigma V^{*},$$

(19)

where $U\in {\mathcal{R}}^{n\times r}$, $\Sigma \in {\mathcal{R}}^{r\times r}$ and $V\in {\mathcal{R}}^{m\times r}$. In this case the reduced order matrix $\tilde{A}$ in (15), can be expresses as

$$\tilde{A}=U^{*}YV{\Sigma }^{†}.$$

(20)

The DMD modes have the following presentation

$$\Phi =YV{\Sigma }^{†}W.$$

(21)

This approach to implementing the DMD method is called exact DMD, since Tu et al. [16] proves that DMD modes computed by (21) are the exact eigenvectors of A. DMD modes computed by (17) are known as projected eigenvectors of A. See [38,39] for some other results.

In this case, the projected matrix of $\mathcal{D}$, in (3), has the following presentation:

$$\tilde{\mathcal{D}}=U^{*}\mathcal{D}=[{\tilde{\mathbf{x}}}_0,\dots ,{\tilde{\mathbf{x}}}_m]$$

(22)

or in equivalent block-matrix form

$$\tilde{\mathcal{D}}=\left[UX|U{\mathbf{x}}_m\right].$$

(23)

If $\mathcal{D}$ is a full-rank matrix, then (23) has the form

$$\tilde{\mathcal{D}}=\left[\Sigma V^{*}|U{\mathbf{x}}_m\right].$$

(24)

(ii). The left singular vector matrix of $\mathcal{D}$. We can choose the transformation matrix Q, in (13), to be

$$Q=U_{\mathcal{D}},$$

(25)

where $U_{\mathcal{D}}$ is from the truncated SVD of the full data matrix $\mathcal{D}$:

$$\mathcal{D}=U_{\mathcal{D}}{\Sigma }_{\mathcal{D}}{V}_{\mathcal{D}}^{*},$$

(26)

where $U_{\mathcal{D}}\in {\mathcal{R}}^{n\times r}$, ${\Sigma }_{\mathcal{D}}\in {\mathcal{R}}^{r\times r}$ and $V_{\mathcal{D}}\in {\mathcal{R}}^{m\times r}$, see  [22].

Then, the projected matrix of $\mathcal{D}$, in (3), has the following presentation:

$$\tilde{\mathcal{D}}=U_{\mathcal{D}}^{*}\mathcal{D}=[{\tilde{\mathbf{x}}}_0,\dots ,{\tilde{\mathbf{x}}}_m]$$

(27)

and if $\mathcal{D}$ is a full-rank matrix, then

$$\tilde{\mathcal{D}}={\Sigma }_{\mathcal{D}}{V}_{\mathcal{D}}^{*}.$$

(28)

The matrix of DMD modes In this case, is

$$\Phi =U_{\mathcal{D}}W,$$

(29)

where W is the eigenvector matrix of $\tilde{A}=U_{\mathcal{D}}^{*}A{U}_{\mathcal{D}}$.

1.3. Optimal Amplitudes of DMD Modes

Finding the DMD mode amplitudes that best fit the DMD modes of a collection of data is referred to as the reconstruction problem. In the context of DMD, reduced-order modeling seeks to identify a subset of DMD modes that perform well in data reconstruction for a data set or a variety of data sets.

Let us consider again Equation (6), which represents the DMD reconstruction of data snapshots $\mathcal{D}$. In the standard DMD approach the vector of amplitudes is computed by

$$\mathbf{b}={\Phi }^{†}{\mathbf{x}}_0$$

(30)

as shown in (6). It is possible to improve this estimate with optimization over all snapshots.

It is straightforward to show that (6) has the following equivalent expression:

$$\mathcal{D}=\Phi \mathrm{diag}\left\{b_i\right\}{V}_{and}\left(\lambda \right),$$

(31)

where $V_{and}\left(\lambda \right)$ is a Vandermonde matrix

$$V_{and}\left(\lambda \right)=\left(\begin{array}{cccc}1& {\lambda }_1& \dots & {\lambda }_1^m\\ 1& {\lambda }_2& \dots & {\lambda }_2^m\\ ⋮& ⋮& \ddots & ⋮\\ 1& {\lambda }_r& \dots & {\lambda }_r^m\end{array}\right).$$

(32)

This demonstrates that the temporal evolution of the dynamic modes is governed by the Vandermonde matrix, which is determined by the r complex eigenvalues ${\lambda }_i$ of $\tilde{A}$ which contain information about the underlying temporal frequencies and growth/decay rates.

Therefore, determination of the unknown vector of amplitudes $\mathbf{b}$ can be considered as the following optimization problem:

$$\underset{\mathbf{b}}{min}{∥\mathcal{D}-\Phi \mathrm{diag}\left\{b_i\right\}{V}_{and}\left(\lambda \right)∥}_F^2.$$

(33)

Using the truncated SVD of $\mathcal{D}=U_{\mathcal{D}}{\Sigma }_{\mathcal{D}}{V}_{\mathcal{D}}^{*}$, and the definition of the matrix $\Phi =U_{\mathcal{D}}W$ in (17), we bring this problem into the following form:

$$\underset{\mathbf{b}}{min}{∥{\Sigma }_{\mathcal{D}}{V}_{\mathcal{D}}^{*}-W\mathrm{diag}\left\{b_i\right\}{V}_{and}\left(\lambda \right)∥}_F^2$$

(34)

or in equivalent form, by using (27) and (28)

$$\underset{\mathbf{b}}{min}{∥\tilde{\mathcal{D}}-W\mathrm{diag}\left\{b_i\right\}{V}_{and}\left(\lambda \right)∥}_F^2$$

(35)

where W is the eigenvector matrix and $\lambda$ is the eigenvalues vector of the reduced-order operator (15). This is a convex optimization problem that can be solved using standard methods. For instance, we can represent (35) in matrix form as

$$M\mathbf{b}=\mathbf{h},$$

(36)

where $M\in {\mathbf{C}}^{r.(m+1)\times r}$ is the coefficient matrix, $\mathbf{h}$ is the forcing term and the unknown amplitude vector $\mathbf{b}$ as given by

$$M=\left[\begin{array}{c}W\\ W\Lambda \\ \dots \\ W{\Lambda }^m\end{array}\right],\mathbf{h}=\left[\begin{array}{c}{\tilde{\mathbf{x}}}_0\\ {\tilde{\mathbf{x}}}_1\\ \dots \\ {\tilde{\mathbf{x}}}_m\end{array}\right],\mathrm{and}\mathbf{b}=\left[\begin{array}{c}{b}_1\\ b_2\\ \dots \\ b_r\end{array}\right],$$

(37)

where $\Lambda =\mathrm{diag}\left\{{\lambda }_i\right\}\in {\mathbf{C}}^{r\times r}$ is a diagonal matrix formed by the eigenvalues of $\tilde{A}$ in (16). Therefore, we can solve the Equation (36) by least-squares approach

$$\mathbf{b}=M^{†}\mathbf{h},$$

(38)

where the pseudoinverse $M^{†}$ may be computed through SVD of M.

1.4. Delay-Embedding Dynamic Mode Decomposition

Delay-embedding is also an important technique when the temporal or spectral complexity of a dynamical system exceeds the spatial complexity, for example, in systems characterized by a broadband spectrum or spatially undersampled. In this case, we arrive at a “short-and-wide”, rather than a “tall-and-skinny”, data matrix $\mathcal{D}$, and the standard algorithm fails at extracting all relevant spectral features.

Delay-Embedding DMD (or Hankel DMD) overcomes several shortcomings of the standard DMD method by extending its capabilities to handle nonlinear dynamics, non-uniformly sampled data, long-term temporal behavior, high-dimensional datasets, and noisy data. This makes it a more versatile and robust technique for dynamic mode decomposition in various applications. The Takens embedding theorem [54] provides a rigorous framework for analyzing the information content of measurements of a nonlinear dynamical system.

To implement delay-embedding DMD, given the data sequence $\mathcal{D}$ in (3), we stack $s\le m$ time-shifted copies of the data to form the augmented input matrix. The following Hankel matrix H is formed:

$${\mathcal{D}}_{aug}=\left(\begin{array}{cccc}{\mathbf{x}}_1& {\mathbf{x}}_2& \dots & {\mathbf{x}}_{m-s+1}\\ {\mathbf{x}}_2& {\mathbf{x}}_3& \dots & {\mathbf{x}}_{m-s+2}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{x}}_s& {\mathbf{x}}_{s+1}& \dots & {\mathbf{x}}_m\end{array}\right),$$

(39)

where the applied embedding dimension is s. The augmented data matrix ${\mathcal{D}}_{aug}$ is then used in place of $\mathcal{D}$ and processed by the standard DMD algorithm. The DMD algorithm prescribed in Equations (3)–(8) is applied to the augmented matrices $X_{aug},Y_{aug}\in {\mathbf{R}}^{(n.s)\times (m-s)}$ in place of X and Y, giving eigenvalues ${\Phi }_{aug}$ and modes ${\Lambda }_{aug}$. The first n rows of ${\Phi }_{aug}$ correspond to the current (not shifted) time and are used to forecast $\mathbf{x}\left(t\right)$.

Arbabi and Mezić [43] have shown the convergence of this time-shifted approach to the eigenfunctions of the Koopman operator. They also illustrated remarkable improvements in the prediction of simple and complex fluid systems. Further examples and theoretical results on delay-embedding and the Hankel viewpoint of Koopman analysis are given by Brunton et al.  [47] and Kamb et al. [44]. They demonstrated that linear time-delayed models are an effective and efficient tool to capture nonlinear and chaotic dynamics.

2. Spatio-Temporal DMD

The idea behind the spatio-temporal extension of the DMD method is to extend the application range of DMD by implementing the simultaneous capture of both spatial and temporal dynamics. This approach is particularly useful for analyzing complex systems where dynamics evolve both in space and time, such as fluid flows, biological systems, and climate phenomena. To our knowledge, the first paper in the literature in which this idea has been attempted is Sharma et al. [55], and later, it was realized by Clainche et al. [23]; see also  [56]. In [23], Le Clainche and Vega introduce spatio-temporal Koopman decomposition (STKD), which incorporates higher order DMD (HODMD) and a spatio-temporal approach for the Koopman operator. For some applications, see [48,49].

In principle, this expansion can be obtained in two ways:

$\left(i\right).$ Sequential method. A temporal DMD algorithm is first applied to the snapshot matrix and a spatial DMD algorithm is applied to the spatial modes. Obviously, the order in which temporal and spatial DMDs are applied can be reversed, and the result of the direct and reverse methods is not identical.

$\left(ii\right).$ Parallel method. Reduced SVD is first applied to the snapshot matrix $\mathcal{D}$, and then, spatial and temporal DMD algorithms are applied to the rescaled left and right singular vector matrices.

In the following, we provide a detailed mathematical description of the parallel STDMD and sequential STDMD approaches.

2.1. Parallel STDMD

The parallel spatio-temporal DMD method simultaneously decomposes spatio-temporal data across both spatial and temporal dimensions, providing insights into the interplay between spatial and temporal dynamics.

Let us recall that the DMD algorithm presented in Section 1.1 uses a low-rank approximation of the linear mapping that best approximates the dynamics of the data $\mathcal{D}$, in (3), collected for the system. Moreover, if we choose the projection matrix to be the matrix $U_{\mathcal{D}}$ from the truncated SVD of the full data matrix $\mathcal{D}$, as shown in (26)

$$\mathcal{D}=U_{\mathcal{D}}{\Sigma }_{\mathcal{D}}{V}_{\mathcal{D}}^{*},$$

we obtain the reduced-order model given by the following data matrix

$$\tilde{\mathcal{D}}=[{\tilde{\mathbf{x}}}_0,\dots ,{\tilde{\mathbf{x}}}_m],$$

which coincides with the scaled right singular vector matrix of $\mathcal{D}$, i.e.,

$$\tilde{\mathcal{D}}={\Sigma }_{\mathcal{D}}{V}_{\mathcal{D}}^{*},$$

according to (27) and (28). Applying the standard DMD approach to reduced model data $\tilde{\mathcal{D}}$, we obtain the following expansion according to (6):

$${\tilde{\mathbf{x}}}_k=W{\Lambda }^k\mathbf{b},$$

(40)

where W is the eigenvector matrix, $\Lambda =\mathrm{diag}\left\{{\lambda }_i\right\}$ is the diagonal matrix of associated eigenvalues of the corresponding DMD operator, and $\mathbf{b}=W^{-1}{\tilde{\mathbf{x}}}_0$. For our purposes, we will call W the matrix of temporal DMD modes and $\Lambda$ the matrix of temporal DMD eigenvalues.

Using equality (13), by multiplying the left side of Equation (40) by matrix $U_{\mathcal{D}}$, we obtain the temporal DMD expansion, in (8):

$${\mathbf{x}}_k=\Phi {\Lambda }^k\mathbf{b}.$$

Following the same idea, we can use the row vectors of the scaled left singular vector matrix $U{\Sigma }_{\mathcal{D}}$ of $\mathcal{D}$ to obtain a spatial expansion similar to (40). Let us denote

$$\overline{\mathcal{D}}={\Sigma }_{\mathcal{D}}{U}^T=[{\overline{\mathbf{y}}}_0,\dots ,{\overline{\mathbf{y}}}_n],$$

(41)

where ${\overline{\mathbf{y}}}_i$ is the i-th column vector of $\overline{\mathcal{D}}$. Applying the standard DMD approach to data $\overline{\mathcal{D}}$, we obtain the following expansion, according to (6):

$${\overline{\mathbf{y}}}_k=\overline{W}{\overline{\Lambda }}^k\overline{\mathbf{b}},$$

(42)

where $\overline{W}$ is the eigenvector matrix, $\overline{\Lambda }=\mathrm{diag}\left\{{\overline{\lambda }}_i\right\}$ is the diagonal matrix of associated eigenvalues of the corresponding DMD operator, and $\overline{\mathbf{b}}={\overline{W}}^{-1}{\overline{\mathbf{y}}}_0$. We will call $\overline{W}$ the matrix of spatial DMD modes and $\overline{\Lambda }$ the matrix of spatial DMD eigenvalues.

From expressions (40) and (42), using (31), obtain

$$\tilde{\mathcal{D}}=W\mathrm{diag}\left\{b_i\right\}{V}_{and}\left(\lambda \right)\mathrm{and}\overline{\mathcal{D}}=\overline{W}\mathrm{diag}\left\{{\overline{b}}_i\right\}{V}_{and}\left(\overline{\lambda }\right).$$

(43)

Then, for the full-data matrix $\mathcal{D}$, using equality

$$\mathcal{D}=\left(U_{\mathcal{D}}{\Sigma }_{\mathcal{D}}\right){\Sigma }_{\mathcal{D}}^{-1}\left({\Sigma }_{\mathcal{D}}{V}_{\mathcal{D}}^{*}\right),$$

we obtain the matrix form presentation

$$\mathcal{D}=V_{and}^T\left(\overline{\lambda }\right)\mathrm{diag}\left\{{\overline{b}}_i\right\}\Psi \mathrm{diag}\left\{b_i\right\}{V}_{and}\left(\lambda \right),$$

(44)

where $r\times r$ the matrix

$$\Psi ={\overline{W}}^T{\Sigma }_{\mathcal{D}}^{-1}W$$

(45)

is the matrix of spatio-temporal DMD modes.

The following algorithm (Algorithm 1) summarizes the steps for parallel STDMD:

Algorithm 1: Parallel STDMD algorithm

1. Compute the (reduced) SVD of $\mathcal{D}$, writing $\mathcal{D}=U_{\mathcal{D}}{\Sigma }_{\mathcal{D}}{V}_{\mathcal{D}}^{*}$.     2. Define spatial and temporal data matrices: $\tilde{\mathcal{D}}={\Sigma }_{\mathcal{D}}{V}_{\mathcal{D}}^{*}\mathrm{and}\overline{\mathcal{D}}={\Sigma }_{\mathcal{D}}{U}^T.$     3. Perform the standard DMD approach to data set $\tilde{\mathcal{D}}$ and compute temporal DMD modes, eigenvalues and amplitudes: $W,\Lambda \mathrm{and}\mathbf{b}.$     4. Perform the standard DMD approach to data set $\overline{\mathcal{D}}$ and compute spatial DMD modes, eigenvalues and amplitudes: $\overline{W},\overline{\Lambda }\mathrm{and}\overline{\mathbf{b}}.$     5. Compute the matrix of spatio-temporal DMD modes $\Psi ={\overline{W}}^T{\Sigma }_{\mathcal{D}}^{-1}W.$

The eigenvalues and DMD modes can then be used to reconstruct the full data ${\mathbf{x}}_k$ in $\mathcal{D}$. Let us denote the elements of snapshot ${\mathbf{x}}_k$ and matrix $\Psi$ as follows:

$${\mathbf{x}}_k={[x_1^{\left(k\right)},x_2^{\left(k\right)},\dots ,x_n^{\left(k\right)}]}^T\mathrm{and}\Psi ={\left[{\psi }_{ij}\right]}_{r\times r}.$$

Then, from (44), for the s-th coordinate of ${\mathbf{x}}_k$ it follows:

$$x_s^{\left(k\right)}=\sum _{i,j=1}^r{\psi }_{ij}{\overline{\lambda }}_i^s{\overline{b}}_i{\lambda }_j^k{b}_j,$$

(46)

where ${\overline{\lambda }}_i$ and ${\lambda }_j$ are the spatial and temporal DMD eigenvalues, respectively.

2.2. Sequential STDMD

In contrast to parallel STDMD, sequential involves decomposing spatio-temporal data sequentially along the temporal axis, capturing both spatial and temporal dynamics separately. This approach enables the identification of spatial structures evolving over time and their corresponding temporal dynamics.

For conventional DMD, the temporal information (temporal growth rate and angular frequency) is explicitly included in the eigenvalue matrix $\Lambda$, whereas the spatial information (spatial growth rate and wavenumber) is implicitly hidden in the dynamic mode matrix $\Phi$. Therefore, this study aims to decompose dynamic modes in a certain way to obtain spatial information.

Let us apply the standard DMD method described in Section 1.1 to the input data $\mathcal{D}$ specified in (3), which results in temporal DMD expansion (8):

$${\mathbf{x}}_k=\Phi {\Lambda }^k\mathbf{b},$$

where $\Phi =YV{\Sigma }^{†}W$ is the matrix of (exact) DMD modes, $\Lambda$ is the matrix of DMD eigenvalues and $\mathbf{b}$ is the vector of amplitudes; see (15)–(21). As we mentioned, this expression is equivalent to (31):

$$\mathcal{D}=\Phi \mathrm{diag}\left\{b_i\right\}{V}_{and}\left(\lambda \right).$$

Note that the spatial information, such as spatial growth rate and wavenumber, of the dynamic in consideration is implicitly hidden in the dynamic mode matrix $\Phi$. We can use the row vectors of the DMD mode matrix $\Phi$ to obtain spatial expansion similar to (40). Let us denote

$$\overline{\mathcal{D}}={\Phi }^T=[{\overline{\mathbf{y}}}_0,\dots ,{\overline{\mathbf{y}}}_n],$$

(47)

where ${\overline{\mathbf{y}}}_i$ is the i-th column vector of $\overline{\mathcal{D}}$. Applying the standard DMD approach to data $\overline{\mathcal{D}}$, we obtain the following expansion according to (6):

$${\overline{\mathbf{y}}}_k=\overline{\Phi }{\overline{\Lambda }}^k\overline{\mathbf{b}},$$

(48)

where $\overline{\Phi }$ is the eigenvector matrix, $\overline{\Lambda }=\mathrm{diag}\left\{{\overline{\lambda }}_i\right\}$ is the diagonal matrix of associated eigenvalues of the corresponding DMD operator, and $\overline{\mathbf{b}}={\overline{\Phi }}^{-1}{\overline{\mathbf{y}}}_0$.

Then, for the full-data matrix $\mathcal{D}$, we obtain the following matrix form presentation:

$$\mathcal{D}=V_{and}^T\left(\overline{\lambda }\right)\mathrm{diag}\left\{{\overline{b}}_i\right\}\Psi \mathrm{diag}\left\{b_i\right\}{V}_{and}\left(\lambda \right),$$

(49)

where $r\times r$ matrix

$$\Psi ={\overline{W}}^T$$

(50)

is the matrix of spatio-temporal DMD modes. The following algorithm (Algorithm 2) summarizes the steps for sequential STDMD:

Algorithm 2: Sequential STDMD algorithm

1. Perform the standard DMD approach to data set $\mathcal{D}$: 1.1. Define the data matrices: X and Y; 1.2. Compute the reduced SVD of $X:X=U\Sigma V^{*};$ 1.3. Construct the reduced-order operator: $\tilde{A}=U^{*}YV{\Sigma }^{†};$ and compute the eigendecomposition of $\tilde{A}:\tilde{A}W=W\Lambda ;$ 1.4. Compute the DMD modes, eigenvalues and amplitudes: $\Phi =YV{\Sigma }^{†}W,\Lambda \mathrm{and}\mathbf{b}.$     2. Define the spatial data matrix as transposed DMD modes: $\overline{\mathcal{D}}={\Phi }^T.$     3. Perform the standard DMD approach to data set $\overline{\mathcal{D}}$ and compute DMD modes, eigenvalues and amplitudes: $\overline{\Phi },\overline{\Lambda }\mathrm{and}\overline{\mathbf{b}}.$     4. Compute the matrix of spatio-temporal DMD modes $\Psi ={\overline{\Phi }}^T.$

For the reconstruction of snapshots in $\mathcal{D}$, we obtain similar to (46) expression

$$x_s^{\left(k\right)}=\sum _{i,j=1}^r{\psi }_{ij}{\overline{\lambda }}_i^s{\overline{b}}_i{\lambda }_j^k{b}_j,$$

(51)

where $x_s^{\left(k\right)}$ is the s-th coordinate of state ${\mathbf{x}}_k$. Note that although the notations of parameters in (51) and (46) are the same, their values are different.

For both cases, in (46) and (51), it is straightforward to obtain the expression for the continuous case, in the form

$$x(s,t)=\sum _{i,j=1}^r{\psi }_{ij}{e}^{{\overline{\omega }}_{i}s}{\overline{b}}_i{e}^{{\omega }_{j}t}{b}_j=\sum _{i,j=1}^r{\psi }_{ij}{\overline{b}}_i{b}_j{e}^{{\overline{\omega }}_{i}s+{\omega }_{j}t},$$

(52)

where s denotes the spatial variable. The spatial DMD eigenvalues ${\overline{\omega }}_i$ give the information about spatial wavenumbers and growth rates, while the temporal DMD eigenvalues ${\omega }_j$ give information about the underlying temporal frequencies and growth rates.

3. Delay-Embedding STDMD

As already mentioned, traditional DMD approaches are limited in their ability to capture the full complexity of nonlinear and non-stationary systems, particularly when dealing with high-dimensional and noisy datasets. Due to the fact that in Algorithms 1 and 2 the standard DMD method is applied sequentially or in parallel, they inherit the disadvantages of the DMD method. To address these limitations, we will propose an extension of STDMD algorithms using the delay-embedding approach described in Section 1.

3.1. Parallel Delay-Embedding STDMD

This approach redesigns the input data of the system, creating new state variables. However, the introduction of the new variables is made at the expense of reducing the number of samples in the training data set. Hence, the number of these new variables (number of rows in the Hankel matrix), in (39), has to be a balance between the ability to detect dominant modes and the accuracy of the estimated model. The following algorithm (Algorithm 3) provides a step-by-step implementation of parallel delay-embedding DMD:

Algorithm 3: Parallel delay-embedding STDMD

1. Compute the (reduced) SVD of $\mathcal{D}$, writing $\mathcal{D}=U_{\mathcal{D}}{\Sigma }_{\mathcal{D}}{V}_{\mathcal{D}}^{*}.$     2. Define spatial and temporal data matrices: $\tilde{\mathcal{D}}={\Sigma }_{\mathcal{D}}{V}_{\mathcal{D}}^{*}\mathrm{and}\overline{\mathcal{D}}={\Sigma }_{\mathcal{D}}{U}^T.$     3. Perform delay-embedding DMD approach to data set $\tilde{\mathcal{D}}$ and compute temporal DMD modes, eigenvalues and amplitudes: $W,\Lambda \mathrm{and}\mathbf{b}.$     4. Perform delay-embedding DMD approach to data set $\overline{\mathcal{D}}$ and compute spatial DMD modes, eigenvalues and amplitudes: $\overline{W},\overline{\Lambda }\mathrm{and}\overline{\mathbf{b}}.$     5. Compute the matrix of spatio-temporal DMD modes $\Psi ={\overline{W}}^T{\Sigma }_{\mathcal{D}}^{-1}W.$

The implementation of the corresponding algorithm for the sequential STDMD approach is similar, so we will omit it here. We note that although delay-embedding is only applied to the reduced input matrices $\tilde{\mathcal{D}}$ and $\tilde{\mathcal{D}}$, the embedding can be applied to the full input data matrix $\mathcal{D}$ as well.

3.2. Delay-Embedding STDMD vs. STKD

Le Clainche et al. [23] introduced spatio-temporal Koopman decomposition (STKD), which incorporates higher order DMD [56] and a spatio-temporal approach for the Koopman operator. Our goal in this section is to present an alternative approach to STKD. Below, we outline some similarities and differences between the two approaches.

• The STKD scheme is similar to delay-embedding STDMD, with the difference that STKD uses higher order DMD instead of augmented DMD. This implies greater computational complexity in STKD than in STDMD, but on the other hand, it allows easily extending the STKD to higher spatial dimensions when the snapshot matrix becomes a tensor.
• In [23], the spatio-temporal expansion, by STKD, corresponding to (52) has the form

  $$x(s,t)=\sum _{i,j=1}^{m,n}{a}_{ij}{q}_{ij}{e}^{{\widehat{\omega }}_{i}s+{\tilde{\omega }}_{j}t},$$

  where $q_{ij}$ are the normalized spatial modes and $a_{ij}$ are the mode amplitudes. In addition, for STKD, the amplitudes $a_{ij}$ are determined by the optimal amplitudes computation scheme described in (34)–(38). We should note that although the calculation of the amplitudes in STDMD is through the standard Formula (30), which is more cost-efficient, a better approximation of the data is achieved when using the schemes (34)–(38).

Among the main advantages of the schemes proposed in this article are the following:

• The matrix form presentations (44) and (49) of the snapshot matrix $\mathcal{D}$, which offer a structured framework that facilitates easier understanding and implementation.
• Additionally, the structured nature of the matrix representation allows for straightforward generalization to parallel computation architectures, enabling seamless scalability and improved computational efficiency in analyzing large-scale spatio-temporal datasets.
• The STDMD approach, augmented with delay-embedding, offers enhanced computational efficiency compared to STKD.
• By augmenting the dataset with delayed observations, the analysis captures underlying dynamics more effectively, reducing the impact of noise on mode identification and reconstruction.

Overall, the delay-embedding STDMD enhances the accessibility and usability of the proposed approaches, making them more practical and widely applicable to researchers and practitioners in various fields.

4. Numerical Examples

In this section, we will illustrate the introduced approach to delay-embedding spatio-temporal DMD. The considered examples are well known in the literature, and through them, we illustrate the ability of the proposed scheme to accurately calculate spatio-temporal DMD modes and eigenvalues, including spatial wavenumbers and growth rates and temporal frequencies and growth rates. We mainly present the results of the application of parallel delay-embedding STDMD (Algorithm 3). Since both methods use extended data matrices and are computationally comparable, we collate the results obtained by Algorithm 3 with those of the STKDM method presented by Le Clainche in [23]. All numerical experiments and simulations were performed on Windows 7 with MATLAB release R2013a on an Acer Aspire 571G laptop with an Intel(R) Core(TM) i3-2328M CPU at 2.2 GHz and 4 GB of RAM.

Example 1.

Combination of travelling wavetrains.

We begin by demonstrating the feature extraction technique for delay embedding STDMD for a spatio-temporal signal:

$$x(s,t)=\left[0.5+sin\left(s\right)\right]\left[2cos(k_{1}s-{\omega }_{1}t)+0.5cos(k_{2}s-{\omega }_{2}t)\right],$$

(53)

defined in a 1D periodic domain, $s\in [0,2\pi )$. This example was taken from [23] and the same example was also discussed in [57]. It represents a simplified model of the signal proposed in [58] with three basic features in the convective variability of the tropical atmosphere as a function of longitude $\left(s\right)$:

(i) 

A time-independent profile, $0.5+sin\left(s\right)$, representing enhanced convective activity over warm oceans over cold oceans such and continental land;

(ii) 

A long-wavelength eastward-propagating wave, $cos(k_{1}s-{\omega }_{1}t)$, representing a large-scale mode of organized convection called Madden-Julian oscillation (MJO);

(iii) 

A short-wavelength westward-propagating wave representing the building blocks of the MJO (so-called convectively coupled equatorial waves).

The natural time units in (53) are days, so the long wave has a period of 45 days and the period of the short wave is approximately 14 days. These periods are comparable to the timescales observed in nature.

In (53), $k_1$ and $k_2$ are integer-valued wavenumbers set to $k_1=2$ and $k_2=10$, and ${\omega }_1$ and ${\omega }_2$ are time-dependent phases for the rationally independent frequencies ${\omega }_1=2\pi /45$ and ${\omega }_2=\sqrt{10}{\omega }_1$. The color map of this pattern is depicted in Figure 1 (left). This pattern is obtained with a spectral spatial and temporal complexity of 12 and 4, respectively. This is because it involves 12 wavenumbers: $\pm k_1,\pm k_2,\pm (k_1\pm 1)$ and $\pm (k_2\pm 1)$, and four frequencies: $\pm {\omega }_1$ and $\pm {\omega }_2$. This pattern is spatially periodic, with a period equal to $2\pi$, but temporally quasi-periodic.

In order to apply the delay-embedding STDMD method, we discretize s and t in the sampled intervals $0\le s\le 10$ and $0\le t\le 5$, using 50 and 25 points, respectively. Generated data are $50\times 25$, but its rank is 4 (see Figure 1 (right)), which yields unsatisfactory results with the pure temporal DMD method.

Performing delay-embedding STDMD (Algorithm 1), with time-delaying index 2 and spatial-delaying index 3, we identify the correct 12 wavenumbers and 4 frequencies. See the dynamic reconstruction with delayed STDMD (Algorithm 1) in Figure 2. Figure 3 depicts the amplitude–frequency and growth rate–frequency diagrams. Figure 3 shows the combinations of spatial modes and temporal modes used in the reconstruction of the data in (52). They are grouped along straight lines in the plane, which may be either horizontal or oblique, and correspond to either standing or travelling patterns, respectively. The results are identical to those in [23], where the STKD method is applied to the same example and input data.

Example 2.

Dynamics of two counter-propagating waves

In this example, we consider the dynamics of two counter-propagating waves

$$x(s,t)=v(s,t)+v(-s,t),$$

(54)

where v is defined as

$$v(s,t)=\frac{1}{2}\sum _{-6}^6{3}^{-\left|m\right|}{e}^{im(10\pi s+30t)}.$$

(55)

The color map of this pattern is depicted in Figure 4. The two counter-propagating waves are visible on the chart, but it is seen that the pattern can also be considered as a modulated standing waves, in which the positions of the nodes and crests do not remain constant, but oscillate left and right. The generated data have a low-rank structure, which can be seen from the singular values depicted in Figure 4.

If we apply the standard DMD approach, we obtain only seven modes and it gives poor reconstruction of the input data. Instead, if we use delay-embedding STDMD (Algorithm 1), with time delay of 2 and spatial delay also of 2, then we obtain 13 modes and reconstruct the input data with greater accuracy. See the dynamic reconstruction with delayed STDMD (Algorithm 1) in Figure 5.

Note that, if we use optimal amplitude computation, in Algorithm 1, as shown in (34)–(38), we obtain a better approximation of the dynamics data and reconstruct the snapshots with a relative RMS error ∼$1.4\times {10}^{-12}$. Figure 6 depicts the amplitude–frequency and growth rate–frequency diagrams. It shows that the relevant points are aligned in two straight lines, which, according to (54), is consistent with the fact that the pattern is the superposition of two counter-propagating travelling waves.

Note that the counterpart to Equation (55) is given by

$$v(s,t)=\frac{1}{2}\sum _{-6}^6{3}^{-\left|m\right|}{e}^{30imt}cos\left(10m\pi s\right),$$

(56)

which implies (from the equality $cos\left(10m\pi s\right)=cos(-10m\pi s)$) that the spatial complexity is 7, while the spectral complexity is 13. It also follows from (56) that the pattern can be seen as modulated standing waves.

5. Conclusions

In this paper, we have provided a detailed exposition of two variants of spatio-temporal dynamic mode decomposition (STDMD), namely the parallel methods STDMD and sequential STDMD. We have introduced the matrix representations underlying these techniques, highlighting their respective computational frameworks for analyzing spatio-temporal data. To address some shortcomings of the presented algorithms, which are inherited from the classic DMD algorithm, we have introduced extensions to these approaches incorporating delay-embedding techniques. Furthermore, we have conducted numerical experiments to validate the efficacy of the proposed extensions in overcoming the identified limitations of traditional DMD methods. Through these experiments, we have illustrated the enhanced performance of delay-embedded STDMD, showcasing its utility in analyzing complex spatio-temporal datasets.

For future work, there are several promising directions that can build upon the methodologies and findings presented in this paper. Firstly, further exploration and refinement of the delay-embedding techniques introduced in our study could lead to more effective approaches for capturing nonlinear dynamics and improving robustness against noise in spatio-temporal DMD analyses. Additionally, investigating the application of our sequential and parallel approaches with delay-embedding across a wide range of spatio-temporal datasets and real-world applications would provide valuable insights into their generalizability and practical utility. Furthermore, exploring hybrid methodologies that combine elements of different spatio-temporal decomposition techniques, such as incorporating machine learning algorithms or Bayesian approaches, could offer new avenues for enhancing the accuracy and interpretability of spatio-temporal analysis. Additionally, future research could focus on the parallel implementation of the approaches introduced in this paper to enhance computational efficiency. Investigating strategies for parallelizing the computation of spatio-temporal DMD algorithms across multiple processing units or distributed computing architectures could significantly reduce computational time and enable the analysis of large-scale datasets. By investigating the parallel implementation of these techniques, future research can enhance their computational efficiency and facilitate their widespread adoption in scientific and engineering domains where timely analysis of spatio-temporal data is critical. Overall, these future directions hold great potential for advancing the state-of-the-art in spatio-temporal DMD methodologies and their applications in diverse fields.

In conclusion, our research contributes to the advancement of spatio-temporal DMD methodologies by introducing extensions that enhance the robustness and accuracy of the analysis. The proposed approaches offer valuable tools for researchers and practitioners in diverse fields, enabling deeper insights into the dynamics of complex spatio-temporal systems. We anticipate that our findings will stimulate further research and development in this area, leading to continued advancements in the analysis and understanding of spatio-temporal phenomena.