Deep Learning Enhanced Dynamic Mode Decomposition

One particular collection of methods that has garnered significant attention and development is those built around the approximation of Koopman operators [1], which are broadly described as dynamic mode decomposition (DMD). These methods in some sense fit within the larger context of modal decompositions [2, 3], but in several ways, they go further insofar as they also readily generate proxies for the flow maps connected to the otherwise unknown but underlying dynamical systems generating the data in question. The first papers in this area [4, 5] showed across a handful of problems that with a surprisingly straightforward approach, one could readily generate accurate models with nothing but measured time series. Further extensions soon appeared by way of the extended DMD (EDMD) and the kernel DMD (KDMD) [6, 7] which generalized the prior approaches in such a way to better exploit the approximation methodologies used in Reproducing Kernel Hilbert Spaces. These methods were used to develop an all purpose approach to model discovery [8]. See [9] for a textbook presentation and summary of what one might describe as the first wave of DMD literature.

However, in the above methods, user choices for how to represent data are necessary, and as shown in [10], even what appear to be natural choices can produce misleading or even outright incorrect results. In response, researchers have brought several supervised learning, or what one might broadly call machine learning (ML), approaches to bear on the problem of finding optimal representations of data for the purpose of approximating Koopman operators. These approaches have included dictionary learning approaches [11, 12, 13] and related kernel learning methods [14]. Likewise, building on theoretical results in [15, 16], more general methods seeking the discovery of optimal topological conjugacies of the flow have been developed [17, 18]. Related work on discovering optimal embeddings of chaotic time series [19] and ML driven model discovery [20, 21, 22] have also recently appeared.

Between the dictionary learning and topological conjugacy approaches then, in this work we present a hybrid methodology that uses ML to discover optimal dictionaries of observables which also act as embeddings of the dynamics that enforce both one-time step accuracy as well as global stability in the approximation of the time series. A similar philosophy is used in [23]. However, in that work autoencoders are used to obtain low-dimensional structure from high-dimensional data; see also the results in [20, 24]. In contrast, in this paper we look at low-dimensional dynamical systems and use autoencoders to find optimal embeddings of the phase-space coordinates. Running EDMD on the embedded coordinates then generates global linearizations of the flow which are characterized by a single set of eigenvalues that govern the dynamics of a given trajectory. Our spectral results may be contrasted with that of [17], in which their autoencoder is paired with an auxiliary network that parameterizes new eigenvalues at each point along the trajectory. We forego this additional parameterization and obtain a more clearly global spectral representation for any given path in the flow.

To assess the validity and robustness of our approach, we explore a number of examples motivated by several different phase-space geometries. We examine how our method deals with the classic harmonic oscillator, the Duffing system, the multiscale Van der Pol oscillator, and the chaotic Lorenz-63 equations. This presents a range of problems in which one has a global center, two centers, slow/fast dynamics on an attractive limit cycle, and finally a strange attractor. As we show, our method is able to learn across these different geometries, thereby establishing confidence that it should be adaptable to more complex problems. While accuracy for test data is high for the planar problems, as expected, the chaotic dynamics in Lorenz-63 cause prediction accuracy to degrade relatively quickly. Nevertheless, with a minimal number of assumptions and guidance in our learning algorithm, we are able to generate a reasonable approximation to the strange attractor, which shows our method should be applicable in a wider range of use cases than those explored here.

In Section 2, we introduce the Koopman operator and outline the DMD for finding a finite-dimensional approximation to the Koopman operator. The role of deep autoencoders in our method is presented in Section 3, and our algorithm and implementation details are explained in Section 4. Results are presented in Section 5. Discussion of results and explorations of future directions are presented in Section 6.

In this work, we seek to generate approximate predictive models for time series, say $\left\{{\bf y}_{j}\right\}_{j=1}^{N_{T}+1}$, which are generated by an unknown dynamical system of the form

where $\mathcal{M}$ is some connected, compact subset of $\mathbb{R}^{N_{s}}$. We denote the affiliated flow of this equation as ${\bf y}(t)=\varphi(t;{\bf x})$, and observables of the state as $g({\bf y}(t))$, where $g:\mathcal{M}\mapsto\mathbb{C}$. To build these approximate predictive models, we build on the seminal work in [1] which shows that there exists a linear representation of the flow map given by the Koopman operator $\mathcal{K}^{t}$ where

Linearity is gained at the cost of turning the problem of predicting a finite-dimensional flow into one of examining the infinite-dimensional structure of the linear operator $\mathcal{K}^{t}$. To wit, if we consider the observables $g$ to be members of the associated Hilbert space of observables, say $L_{2}\left(\mathcal{M},\mathbb{C}\right)$, $g\in L_{2}\left(\mathcal{M},\mathbb{C}\right)$ if $\left|\left|g\right|\right|_{2}<\infty$ where $\left|\left|g\right|\right|_{2}$ is given by

For concreteness, we have chosen to integrate with respect to Lebesgue measure, though of course others could be more convenient or appropriate. We see then that this makes the Koopman operator $\mathcal{K}^{t}$ an infinite-dimensional map such that

In this case, one has

where

Based on the above computations, we observe that the Koopman operator is an isometry for volume preserving flows.

If we further suppose that $\mathcal{M}$ is invariant with respect to the flow, or $\varphi_{t}(\mathcal{M})\subset\mathcal{M}$ for $t>0$, we can simplify the above statement so that

We will assume this requirement going forward since it is necessary to make much in the way of concrete statements about the analytic properties of the Koopman operator. In particular, using Equation (2), if we further suppose that we have an observable $g$ such that $g\in C_{1}(\mathcal{M})$ then we see that

where $\left|\left|\cdot\right|\right|_{\infty}$ denotes the supremum norm over the compact subset $\mathcal{M}$. From this computation, we find that the infinitesimal generator, say $\mathcal{L}$, affiliated with $\mathcal{K}^{t}$ is given by

and we likewise can use the Koopman operator to solve the initial-value problem

so that $u({\bf x},t)=\mathcal{K}^{t}g({\bf x})$. From here, the central question in Koopman analysis is to determine the spectrum and affiliated modes of $\mathcal{L}$ since these then completely determine the behavior of $\mathcal{K}^{t}$. This can, of course, involve solving the classic eigenvalue problem

However, as the reader may have noticed, there has been no discussion of boundary conditions. Therefore, while one can get many useful results focusing on this problem, one must also allow that continuous spectra are a natural feature, and eigenfunctions can have pathological behavior as well as being difficult to categorize completely and in detail; see [16] and [25]. As this issue represents an evolving research area, we do not attempt to make any further claim and simply follow in the existing trend of the literature and focus on trying to numerically solve the classic eigenvalue problem.

To this end, if we can find the Koopman eigenfunctions $\{\phi_{l}\}_{l=1}^{\infty}$ with affiliated eigenvalues $\{\lambda_{l}\}_{l=1}^{\infty}$, where

then for any observable $g$ one, in principle, has a modal decomposition such that

as well as an analytic representation of the associated dynamics

Note, that the analytic details we provide regarding the Koopman operator in this section are those we find most essential to understanding the present work and for providing a reasonably complete reading experience. Essentially all of it exists, and more detail and depth can be found, in [26, 27, 28], among many other sources.

Now, the challenge of determining the modes and eigenvalues of the infinite-dimensional operator, $\mathcal{K}^{t}$, remains. In general, this is impossible to obtain in an analytic way, however, the DMD and its extensions, the EDMD and the KDMD [4, 5, 9, 6, 7], allow for the numerical determination of a finite number of the Koopman modes and eigenvalues. In this work, we focus on the EDMD since it most readily aligns with the methodologies of ML.

The EDMD begins with the data set $\left\{{\bf y}_{j}\right\}_{j=1}^{N_{T}+1}$ where

where $\delta t$ is the time step at which data are sampled. Note, per this definition, we have that ${\bf y}_{1}={\bf x}$. Following [6, 7], given our time snapshots $\left\{{\bf y}_{j}\right\}_{j=1}^{N_{T}+1}$, we suppose that any observable $g({\bf x})$ of interest lives in a finite-dimensional subspace $\mathcal{F}_{D}\subset L_{2}\left(\mathcal{O}\right)$ described by a given basis of observables $\left\{\psi_{l}\right\}_{l=1}^{N_{o}}$ so that

Given this ansatz, we then suppose that

where $r({\bf x};{\bf K}_{o})$ is the associated error which results from the introduction of the finite-dimensional approximation of the Koopman operator represented by the $N_{o}\times N_{o}$ matrix ${\bf K}_{o}$.

While we ultimately find a matrix ${\bf K}_{o}$ to minimize the error $r({\bf x};{\bf K}_{o})$ relative to the choice of observables, for this approach to make any real computational sense, we tacitly make the following assumption when using the EDMD

One can then demonstrate, in line with the larger DMD literature, that finding the error minimizing matrix ${\bf K}_{o}$ is equivalent to solving the optimization problem

where the $N_{o}\times N_{T}$ matrices $\bm{\Psi}_{\pm}$ are given by

where each column in the above matrices is a $N_{o}\times 1$ vector of observables of the form

where $\left|\left|\cdot\right|\right|_{F}$ is the Frobenius norm. In practice, we solve (6) using the Singular Value Decomposition (SVD) of $\bm{\Psi}_{-}$ so that

This then gives us

where $-P$ denotes the Moore–Penrose pseudoinverse. The corresponding error in the Frobenius norm $E_{r}(\bm{K}_{o})$ is given by

We see that $E_{r}({\bf K}_{o})$ serves as a proxy for the error function $r({\bf x};{\bf K}_{o})$.

Following existing methodology [6], if we wished to find Koopman eigenfunctions and eigenvalues, then after diagonalizing ${\bf K}_{o}$ so that

then one can show that the quantity $\lambda_{l}=\ln(\tilde{t}_{l})/\delta t$ should be an approximation to an actual Koopman eigenvalue and

From here, in the traditional EDMD algorithm, one approximates the dynamics via the reconstruction formula

where the Koopman modes ${\bf k}_{l}\in\mathbb{C}^{N_{s}}$ in principle solve the initial-value problem

In the original coordinates of the chosen observables, using Equation (8) we obtain the equivalent formula

where ${\bf K}_{m}$ is the $N_{s}\times N_{o}$ matrix whose columns are the Koopman modes ${\bf k}_{l}$ and $\bm{\Lambda}$ is the $N_{o}\times N_{o}$ diagonal matrix with diagonal entries $\lambda_{l}$. In practice, we find the Koopman modes via the fitting formula

where ${\bf H}$ is any complex $N_{s}\times N_{o}$ matrix and the norm $\left|\left|\cdot\right|\right|_{2}$ refers to the standard Euclidean norm. Of course, the appropriateness of this fit is completely contingent on the degree to which Ansatz 1 holds. We address this issue in Section 3.

Finally, we note that the standard DMD is given by letting $N_{o}=N_{s}$ and $\psi_{l}({\bf x})=x_{l}$. In this case, we have that ${\bf K}_{m}={\bf V}$. Due to its popularity and ease of implementation, we use the DMD for reference in Section 5.

The central dilemma when using the EDMD is finding suitable collections of observables relative to the data stream $\left\{{\bf y}_{j}\right\}_{j=1}^{N_{T}+1}$. To accomplish this in an algorithmic way, we suppose that an optimal set of observables can be found using a deep neural network $\mathcal{E}$ where

and where across dimensions we define

We call this the encoder network and the transformed coordinates, $\tilde{{\bf x}}$, the latent variables. Given that the neural network representing $\mathcal{E}$ consists of almost everywhere smooth transformations, we note that by Sard’s Theorem [29], we can generically assume that $\mathcal{E}$ has a Jacobian of full rank at almost all points in the domain of the encoder. We generically assume that the latent dimension $N_{o}\geq N_{s}$, making $\mathcal{E}$ an immersion [29] from the flow space into the space of observables. In this case, the matrices $\bm{\Psi}_{\pm}$ from equation (7) are now given by

so that we compute ${\bf K}_{o}$ in the latent-space coordinates . Beyond the lower bound given above, we do not enforce any constraint on the number of dimensions of the latent space. Instead, $N_{o}$ is treated as a hyperparameter. This is a key feature of our encoder networks as they are allowed to lift the data into higher-dimensional latent spaces, giving the EDMD a rich, and flexibly defined, space of observables with which to work.

Corresponding to the immersion $\mathcal{E}$, we introduce the decoder network $\mathcal{D}$ which is meant to act as the inverse of $\mathcal{E}$, i.e. we seek a decoder mapping acting as a submersion

so that

Upon finding such a mapping, we can show the following lemma.

As shown in [27, 16], flows which are diffeomorphic to one another share Koopman eigenvalues and have eigenfunctions that are identical up to diffeomorphism. In our case, $\mathcal{E}$ and $\mathcal{D}$ are typically not invertible since the reverse composition, $\mathcal{E}\circ\mathcal{D}$, does not necessarily yield the identity. However, we can define the affiliated flow $\tilde{\varphi}(t;\tilde{{\bf x}})$ such that

Immediately, we find that there must of course be an affiliated Koopman operator, $\tilde{\mathcal{K}}^{t}$, corresponding to this encoded, or lifted, flow. This then allows us to show the following theorem.

Thus we see that every Koopman mode and eigenvalue in the original flow space is lifted up by the embedding $\mathcal{E}$. That said, it is of course possible then that new spectral information affiliated with $\tilde{\mathcal{K}}^{t}$ can appear, and if we perform the EDMD in the lifted variables, there is no immediate guarantee we are, in fact, computing spectral information affiliated with the primary Koopman operator $\mathcal{K}^{t}$. That said, if in fact Ansatz 1 holds for the given choice of observables, which is to say we have made the right choice of embedding $\mathcal{E}$, then we can prove the following theorem.

So we see on the one hand that the EDMD left to its own devices is prone to introducing perhaps spurious spectral information, and of course, without recourse to a known reference, we have no way in advance of knowing how to tell which results generated via the EDMD produce relevant spectra. We note that this issue was numerically illustrated in [10]. On the other hand, if we can somehow ensure Ansatz 1 holds, then the EDMD is guaranteed to produce meaningful results with essentially zero error. Of course, what remains in either case is the fundamental dilemma of how to choose observables such that Ansatz 1 is enforced, or at least such that the error is guaranteed to be controlled in some uniform way.

To address this dilemma then, we propose the deep learning dynamic mode decomposition (DLDMD), in which we determine the encoder/decoder pair $\mathcal{E}$ and $\mathcal{D}$ by minimizing the following loss function $\mathcal{L}$, where

such that

and ${\bf W}_{g}$ contains the weights of the $\mathcal{E}$ and $\mathcal{D}$ networks making the final term in $\mathcal{L}$ a regularization term. The hyperparameters $\alpha_{1},\alpha_{2},\alpha_{3},$ and $\alpha_{4}$ provide appropriate scalings for each loss component. The autoencoder reconstruction loss, $\mathcal{L}_{\text{recon}}$, demands that $\mathcal{D}$ approximates the inverse of $\mathcal{E}$. We see then that $\mathcal{L}_{\text{recon}}$ in effect replaces Equation (9). The DMD loss, $\mathcal{L}_{\text{dmd}}$, keeps $r({\bf x};{\bf K}_{o})$ as small as possible relative to advancing one timestep. In contrast, we see $\mathcal{L}_{\text{pred}}$ makes the overall method remain stable under repeated applications of the finite-dimensional approximation to $\mathcal{K}^{\delta t}$. Thus $\mathcal{L}_{\text{dmd}}$ plays the role of ensuring we have a consistent time-stepping scheme, while $\mathcal{L}_{\text{pred}}$ ensures we have a globally stable method, and so by combining the two, we address the fundamental dilemma presented by Ansatz 1. It is this ensured global stability that motivates us to call the diagonalization of the matrix ${\bf K}_{o}={\bf V}{\bf T}{\bf V}^{-1}$ from our algorithm’s EDMD the global linearization of the flow.

We further see that this loss function implements the diagram in Figure 1. Through this diagram we impose a one-to-one mapping between trajectories in the original coordinates, $\bm{y}$, and trajectories in the latent space, $\tilde{\bm{y}}$. Note that this diagram allows us to circumvent the unknown flow map, $\varphi(\delta t;\bm{y}_{j})$, which is equivalent to the Koopman operator of interest, $\mathcal{K}^{\delta t}$, from Equations (4) and (5).

We build the autoencoder in the Python programming language using Tensorflow version 2. The deep neural networks we construct that act as $\mathcal{E}$ and $\mathcal{D}$ from Equation (10) must transform each vector of coordinates along a sample trajectory to and from the latent-space coordinates , respectively. We chose to use dense layers in each network, though other layer types should suffice so long as they encode each point along the trajectory separately, are densely connected, and output the correct dimensions.

As is the case when training any type of neural network, there are a number of hyperparameters that the researcher must take care in selecting. However, we found that the encoder and decoder networks did not require significantly different hyperparameters from dataset to dataset. Notably, the architecture of the neural networks remained identical across all examples. We found that 3 hidden layers each with 128 neurons were sufficient for all of the test problems presented in Section 5. The primary tunable parameter was the dimension of the latent space, $N_{o}$, which was tuned manually. We used Rectified Linear Units (ReLU) for the activation functions and chose the Adam optimizer with a learning rate of $10^{-3}$ for the harmonic oscillator, and $10^{-4}$ for Duffing, Van der Pol, and Lorenz 63. For the loss function hyperparameters in Equation (11), $\alpha_{1},\alpha_{2}$, and $\alpha_{3}$ were all set to $1$, while $\alpha_{4}=10^{-9}$. The harmonic oscillator was trained with a batch size of 512 while all other systems had batch sizes of 256. All systems were trained for 1,000 epochs. The hardware used was an Nvidia Tesla V100.

See Algorithm 1 for the complete pseudocode for the DLDMD training method. The trained DLDMD model is applied by sending a trajectory through the encoder network, performing the EDMD using the encoded coordinates as observables, then using the modes, eigenvalues, and eigenfunctions to reconstruct the full length of the trajectory and beyond in the latent space. The decoder network then allows us to map the entire EDMD reconstruction back into the original coordinate system.

We test the DLDMD method on several datasets generated from dynamical systems that each exhibit some unique flow feature. In Sections 5.1 - 5.3, we examine much of the range of planar dynamics by way of studying the harmonic oscillator, Duffing, and Van der Pol equations. For example, when limited to the first separatrix, the harmonic oscillator gives closed orbits about a single center. Proceeding in complexity, the Duffing equation is comprised of not only closed orbits but also a homoclinic connection that separates two distinct regions of the phase space, requiring now that the DLDMD compute trajectories on either side of a separatrix. The Van der Pol oscillator has trajectories both growing and decaying onto a limit cycle and exhibits multi-scale, slow/fast dynamics. Finally, in Section 5.4, we demonstrate the DLDMD on chaotic trajectories from the Lorenz-63 system, which extends our approach to three-dimensional results evolving over a strange attractor.

The training, validation, and test data for each example were generated using a fourth-order Runge-Kutta scheme. To generate test-data, for the harmonic oscillator and the Duffing system, we chose a time step of $\delta t=0.05$ and ran simulations out to $t_{f}=20$ in order to get closed orbits for all initial conditions. The Van der Pol system used $\delta_{t}=0.02$ and $t_{f}=15$. This system required a shorter integration step size to sufficiently sample the slow and fast parts of each trajectory. Finally, the Lorenz-63 system used $\delta_{t}=0.01$ and $t_{f}=3$. For testing, we applied the trained encoder and decoder with an EDMD decomposition on the latent trajectories and generated reconstructions up to the simulation times used for training. Then, in each case we ran the trajectories out further in time in their latent-space coordinates using only the spectral information from their respective EDMD to evaluate stability. The prediction time for the harmonic oscillator and Duffing system was from $t_{f}=20$ to $t_{f}=40$, for the Van der Pol equation it was $t_{f}=15$ to $t_{f}=30$, and finally for the Lorenz-63 equations, it was from $t_{f}=3$ to $t_{f}=6$.

We generate 15,000 trajectories for each system, using 10,000 for training, 3,000 for validation, and 2,000 for testing. Each trajectory is generated by uniform random sampling of initial conditions from some pre-defined region of phase space. For the harmonic oscillator, we used $x_{1}\in\{-3.1,3.1\}$, $x_{2}\in\{-2,2\}$, and we limited our choice of trajectories to those within the first separatrix using the potential function $0.5x_{2}^{2}-\cos{x_{1}}<0.99$. The Duffing system was generated over initial conditions sampled from $x_{1}\in\{-1,1\}$ and $x_{2}\in\{-1,1\}$. The Van der Pol system used $x_{1}\in\{-2,2\}$ and $x_{2}\in\{-2,2\}$ to generate trajectories and was then scaled by its standard deviation. The Lorenz-63 system used $x\in\{-15,15\}$, $y\in\{-20,20\}$, and $z\in\{0,40\}$.

In this paper, we have developed a deep learning extension to the EDMD, dubbed DLDMD, which ensures, through the inclusion in our loss function of the terms $\mathcal{L}_{\text{dmd}}$ and $\mathcal{L}_{\text{pred}}$, both the local one step accuracy and global stability, respectively, in the generated approximations. This keeps us numerically close to satisfying the requirements stated in Ansatz (1), which is necessary for the success of the EDMD in the latent variables. Likewise, by constructing a loss function to train the autoencoder using the diagram in Figure 1, the DLDMD learns mappings to and from the latent space which are one-to-one. Thus we ensure that all Koopman modes and eigenvalues are captured in this latent space as well. These results taken together ensure that the DLDMD finds a global linearization of the flow, facilitating a straightforward spectral characterization of any analyzed trajectory and an accurate prediction model for future dynamics.

Moreover, we have developed a ML approach which requires a relatively minimal number of assumptions to ensure successful training. We have demonstrated this across a relatively wide array of dynamical phenomena including chaotic dynamics; see Table 1 for a performance summary of the DLDMD across all presented examples. We are able to cover so much ground by way of implementing relatively straightforward higher-dimensional embeddings via the encoder network, $\mathcal{E}$. The latent-space coordinates of each trajectory generally show that the encoder produces modes with near monochromatic Fourier spectra; again see Figures 4, 8, 11, and 15. This is an especially compelling result of the DLDMD, and it speaks to the power of deep learning methods that such elegant representations can be discovered with so little user guidance. However, by lifting into higher-dimensional latent spaces we do lose the strict topological equivalence in [16] and cannot guarantee that all spectra are invariant. This issue does not seem to manifest in any measurable way in our results, but it should be kept in mind when pursuing future work.

The DLDMD is not without other drawbacks. In particular, the latent dimension parameter $N_{o}$ is critical to the success or failure of the DLDMD. Unfortunately, there is no readily apparent method for choosing this embedding dimension before training. Therefore, the optimal $N_{o}$ had to be determined by simply training the model at successively larger values for $N_{o}$ and stopping once the error became too large or grew unstable. This approach is time consuming and an obvious disadvantage to the method. Determining more optimal ways to approach the latent dimensionality will be the subject of future research.

Another drawback of the DLDMD is the number of trajectories used during training. For each example we generated 10,000 initial conditions sampled uniformly over the phase space of interest for training. Rarely do real-world datasets provide such a uniform sampling of the space. Methods to cope with sparse observations could potentially add far more utility to the method. Very high-dimensional systems that exhibit low-rank structure also present problems for DLDMD and the more conventional use of autoencoder networks for compression rather than lifting [20, 23] may be more applicable. The chaotic trajectories of the Lorenz-63 system were clearly the most challenging for DLDMD to reproduce. This is evident in the low average MSE loss and in the relatively poor agreement between the test and predicted trajectory. For this reason, chaotic systems will likely require additional innovations, such as incorporating delay embeddings into the EDMD step [30] and is the topic of future work.

Lastly, systems with noisy observations will pose significant challenges to DLDMD as it is reliant on the one-step prediction to enforce consistency in the latent spectra. In this regard, the method would seem to benefit from some of the alternative DMD approaches as found in [31], an issue we will also pursue in future work.

The data that support the findings of this study are openly available at
https://github.com/JayLago/DLDMD.git

This work was supported in part by the Naval Information Warfare Center Pacific (NIWC Pacific) and the Office of Naval Research (ONR).

The authors would like to thank Dr. Erik M. Bollt from the Department of Electrical and Computer Engineering at Clarkson University for his feedback and insights, as well as Joseph Diaz and Robert Simpson for their insightful comments.