Machine Learning Enhanced Hankel Dynamic-Mode Decomposition

While the acquisition of time series has become more straightforward, developing dynamical models from time series is still a challenging and evolving problem domain. Within the last several years, to address this problem, there has been a merging of machine learning tools with what is called the dynamic mode decomposition (DMD). This general approach has been shown to be an especially promising avenue for accurate model development. Building on this prior body of work, we develop a deep learning DMD based method which makes use of the fundamental insight of Takens’ Embedding Theorem to build an adaptive learning scheme that better approximates higher dimensional and chaotic dynamics. We call this method the Deep Learning Hankel DMD (DLHDMD). We likewise explore how our method learns mappings which tend, after successful training, to significantly change the mutual information between dimensions in the dynamics. This appears to be a key feature in enhancing the DMD overall, and it should help provide further insight for developing other deep learning methods for time series analysis and model generation.

This work uses machine learning to develop an accurate method for generating models of chaotic dynamical systems using measurements alone. A number of challenging examples are examined which show the broad utility of the method and point towards its potential impacts in advancing data analysis and modeling in the physical sciences. Finally, we present quantitative studies of the information theoretic behavior of the machine learning tools used in our work, thereby allowing for a more detailed understanding of what can otherwise be an inscrutable method.

The incorporation of modern machine learning methodology into dynamical systems is creating an ever expanding array of techniques pushing the boundaries of what is possible with regards to describing and predicting nonlinear multi-dimensional time series. Longstanding problems such as finding optimal Takens’ embeddings [1, 2] now have powerful and novel deep learning based algorithmic approaches [3] which would not have been feasible even ten years ago. Likewise, the field of equation free modeling using Koopman operator methods, broadly described by Dynamic Mode Decomposition (DMD), has seen several innovative deep learning based methods emerge over the last several years [4, 5, 6] which have been shown to greatly expand the accuracy and flexibility of DMD based approaches. There have also been related and significant advances in model identification and solving nonlinear partial differential equations via deep learning techniques [7, 8, 9, 10].

With this background in mind, in this work we focus on extending the methods in [6] which were called Deep Learning DMD (DLDMD). In that work, a relatively straightforward method merging auto-encoders with the extended DMD (EDMD) was developed. This was done by using an encoder to embed dynamics in a sufficiently high enough dimensional space which then generated a sufficiently large enough space of observables for the EDMD to generate accurate linear models of the embedded dynamics. Decoding then returned the embedded time series to the original variables in such a way as to guarantee the global stability of iterating the linear model to generate both reconstructions and forecasts of the dynamics. The DLDMD was shown to be very effective in finding equation-free models which were able to both reconstruct and then forecast from data coming from planar dynamical systems.

However, when chaotic time series from the Lorenz-63 system were examined, the performance of the DLDMD was found to degrade. While this clearly makes the DLDMD approach limited in its scope, we note that the successful use of DMD based approaches to accurately reconstruct or forecast chaotic dynamics are not readily available. Other methods such as HAVOK [11] or SINDy [12] are more focused on the analysis of chaotic time series or the discovery of model equations which generate chaotic dynamics, though of course if one has an accurate model, then one should be able to generate accurate forecasts. In this vein, there are also methods using reservoir computing (RC) [13, 14], though again, nonlinear models are essentially first learned and then used to generate forecasts. However, both SINDy and RC rely on proposing libraries of terms to build models which are then fit (or learned from) data.

While effective, such approaches do not allow for the spectral or modal analysis which has proven to be such an attractive and useful feature of DMD based methods. Likewise, they require a number of user decisions about how to construct the analytic models used in later regressive fitting that amount to a guess and check approach to generating accurate reconstructions and forecasts. Therefore in this work, using insights coming from the Takens’ Embedding Theorem (TET) [15, 3], we expand the DLDMD framework so as to make it accurate in both generating reconstructions and forecasts of chaotic time series. This is done by first making the EDMD over embedded coordinates global as opposed to the local approach of [6]; see also [4, 5]. Second, we develop an adaptive Hankel matrix based ordering of the embedded coordinates which adds more expressive power for approximating dynamics to the deep learning framework. To study our method, we use data generated by the Lorenz-63 and Rossler systems as well as twelve-dimensional projections of data from the Kuramoto–Sivashinksky (KS) equation. In all of these cases, we show that by combining our proposed modifications to the DLDMD that we are able to generate far more accurate reconstructions and forecasts for chaotic systems than with DLDMD alone. Moreover, we have built a method which still allows for the straightforward modal analysis which DMD affords and keeps user choices to a handful of real-valued hyperparamters while still producing results competitive with other approaches in the literature.

Further, motivated by the classic information theory (IT) studies of the TET [1], as well as modern insights into the role that information plays in deep learning [16, 17], we study how the fully trained encoder changes the information content of the dynamics coming from the Lorenz-63 and Rossler systems. For the Lorenz-63 system, the encoder tends to either slightly decrease the mutual information or cause strong phase shifts which decrease the coupling times across dimensions. However, the characteristic timescales corresponding to lobe switching in the Lorenz ‘butterfly’ are clearly seen to be preserved in the dynamics of the information for the Lorenz-63 system. In contrast then, for the Rossler system, the slow/fast dichotomy in the dynamics seen in the original coordinates is made more uniform so that rapid transients in the information coupling are removed by the encoder. Thus in either case, we see that the encoder generates significant differences in the information content between dimensions in the latent coordinates relative to the original ones, and that this strong change in information content is a critical feature in successful training.

Of course, the present work is ultimately preliminary, and there are a number of important questions left to be resolved. First, while we are able to easily display computed spectra, the affiliated global Koopman modes we find are not as straightforward to show. We generate our results from random initial conditions, so the most effective means of constructing the global Koopman modes would be via radial-basis functions, but the implementation would be nontrivial due to the infamous ill-conditioning issues which can plague the approach [18]. Second, there is a clear need for a comparison across SINDy, RC, and our DLHDMD methods. In particular, the present work generates excellent reconstructions, but so far the predictive horizon is relatively short and difficult to increase. How well other methods address this issue relative to their reconstruction and other diagnostic properties, and then how all of these methods compare in these several different ways is as yet unclear. While acknowledging then the limitations of the present work, we defer addressing the above issues till later works where each of the above issues can be dealt with in the detail that is needed.

The structure of this paper is as follows. In Section 2, we provide an introduction to the Extended DMD and then explain the extensions we develop which are critical to the success of the present work. In Section 3, we introduce the Hankel DMD and, incorporating the extensions introduced in Section 2, we show how well it does and does not perform on several examples. Then in Section 4 we introduce the Deep Learning Hankel DMD and provide results on its performance as well as an analysis of how the mutual information changes in the latent variables. Section 5 presents our results on mutual information. Section 6 provides conclusion and discussion.

To begin, we suppose that we have the data set $\left\{{\bf y}_{j}\right\}_{j=1}^{N_{T}+1}$ where

where $\delta t$ is the time step at which data is sampled and $\varphi(t;{\bf x})$ is a flow map such that $\varphi(t_{1},{\bf x})={\bf x}$. From the flow map, we define the affiliated Koopman operator $\mathcal{K}^{t}$ such that for a given scalar observable $g({\bf x})$, one has

so that the Koopman operator linearly tracks the evolution of the observable along the flow. We likewise define the associated Hilbert space of observables, say $L_{2}\left(\mathbb{R}^{N_{s}},\mathbb{R},\mu\right)$, or more tersely as $L_{2}\left(\mathcal{O}\right)$, so that $g\in L_{2}\left(\mathcal{O}\right)$ if

where $\mu$ is some appropriately chosen measure. This makes the infinite-dimensional Koopman operator $\mathcal{K}^{t}$ a map such that $\mathcal{K}^{t}:L_{2}\left(\mathcal{O}\right)\rightarrow L_{2}\left(\mathcal{O}\right).$

Following [19, 20], 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_{ob}}$ so that

Given this ansatz, we then suppose that

where $r({\bf x};\mathbf{K}_{a})$ is the associated error which results from the introduction of the finite-dimensional approximation of the Koopman operator represented by $\mathbf{K}_{a}$. We can then find ${\bf K}_{a}$ by solving the following minimization problem

where ${\bf a}=(a_{1}\cdots a_{N_{ob}})^{T}$, $\mathbf{\mathbf{\Psi}}_{j}=\left(\psi_{1}({\bf y}_{j})\cdots\psi_{N_{ob}}({\bf y}_{j})\right)^{T}$, the inner product $\left<,\right>$ is the standard one over $\mathbb{C}^{N_{ob}}$, and the $\ast$ symbol denotes complex conjugation. It is straightforward to show that an equivalent and easier to solve form of this optimization problem is given by

where $\left|\left|\cdot\right|\right|_{F}$ is the Frobenius norm, and the $N_{ob}\times N_{T}$ matrices $\mathbf{\Psi}_{\pm}$ are given by

In practice, we solve this equation using the Singular-Value Decomposition (SVD) of $\mathbf{\Psi}_{-}$ so that

This then gives us

To complete the algorithm, after diagonalizing $\mathbf{K}_{a}$ so that

then one can show that the Koopman eigenfunctions $\phi_{l}({\bf y}_{j})$ are found via the equations

From here, one can, starting from the initial conditions, approximate the dynamics via the reconstruction formula

where $\lambda_{l}=\ln(\ell_{l})/\delta t$ and the Koopman modes ${\bf k}_{l}\in\mathbb{C}^{N_{s}}$ solve the initial-value problem

Again, in matrix/vector notation, keeping in mind that ${\bf x}\in\mathbb{R}^{N_{s}}$ and that in general $N_{s}\neq N_{ob}$, we have

where ${\bf K}_{M}$ is the $N_{s}\times N_{ob}$ matrix whose columns are the Koopman modes ${\bf k}_{j}$. As can be seen then, generically, one can only find the Koopman modes through least-squares solutions of the non-square problem. In this regard, one would do well to have information from as many initial conditions as possible to over-determine the problem.

When implementing EDMD, the most natural observables are the projections along the canonical Cartesian axes, i.e.

If we stick to this space of observables, the EDMD method reduces to the standard DMD method. Thus the idea with EDMD is to include more nonlinear observables to hopefully represent a richer subspace of dynamics and thereby make the approximation of corresponding Koopman operator more accurate and sophisticated.

With this in mind, [15] built upon the classic idea of Takens embeddings [21] and explored using affiliated Hankel matrices to generate natural spaces of observables for EDMD, and approach we describe as Hankel DMD (HDMD). Also of note in this direction is the HAVOK method developed in [11], though in some ways HAVOK is more akin to the embedology methods explored in such classic works as [22, 2].

HDMD thus begins with an affiliated scalar measurement of our time series, say $\left\{g({\bf y}_{j})\right\}_{j=1}^{N_{T}+1}$. From this, by introducing a window size $N_{w}$ one builds the affiliated Hankel matrix ${\bf H}_{g}\left({\bf x}\right)$ where

where the number of observables $N_{ob}=N_{T}+1-(N_{w}-1)$.

What one sees then is that each row of ${\bf H}_{g}({\bf x})$ is some iteration of the Koopman operator $\mathcal{K}^{\delta t}$. From here then, each row of $N_{w}$ time steps is defined to be its own separate observable $\psi_{l}({\bf x})$, i.e.

One then proceeds as above with the EDMD algorithm, where we emphasize that $N_{T}$ is replaced by $N_{w}-1$. This is an interesting feature, or arguably limitation, of the HDMD method in which we generate matrices $\mathbf{\Phi}_{\pm}$ (see Equation (4)) up to the time index $N_{w}-1\leq N_{T}$. Thus later times are used to build approximations at prior times. This makes the issue of forecasting data more difficult since one must iterate the EDMD results, as is done via Equation (6), from time index $N_{w}-1$ up to $N_{T}$ to reconstruct the original data that was used in the first place. Throughout the remainder of the paper then, we take care to distinguish between iterated reconstructions and actual forecasts which make novel predictions beyond the given data.

Finishing our explanation of HDMD, if one has data along multiple initial conditions, say $\left\{{\bf x}_{k}\right\}_{k=1}^{N_{C}}$, we can extend the above algorithm by concatenating Hankel matrices so that we perform EDMD on the combined matrix ${\bf H}_{C}$ so that

The inclusion of other observables can be done in a similar fashion.

To improve the HDMD such that it is able to deal with chaotic systems such as the Lorenz equation, we now turn to and adapt the framework of the deep learning DMD (DLDMD) developed in [6]. Our deep learning enhanced HDMD begins with an autoencoder composed of neural networks $\mathcal{E}$ (the encoder) and $\mathcal{D}$ (the decoder) such that

and such that our auto-encoder is a near identity, i.e.

Note, we call the encoded coordinates latent variables or latent dimensions in line with the larger literature on machine learning. Motivated by our results in [6], we take $N_{e}\geq N_{s}$. We then have in the later EDMD step that $N_{ob}=N_{e}\bar{N}_{ob}$, so that larger values of $N_{e}$ correspond to more observables to represent the dynamics through an approximation to the Koopman operator.

The encoded coordinates should represent a set of observables which should enhance the overall accuracy of HDMD approximations of the dynamics. To train for this, after making reasonable choices for how to initialize the weights of the auto-encoder, and fixing a choice for $N_{st}$, given the training data $\left\{{\bf y}_{j,k}\right\}_{j,k=1}^{N_{T}+1,N_{C}}$, after shuffling into $\lfloor N_{C}/N_{B}\rfloor$ batches of batch size $N_{B}$, over a given batch ${\bf Y}=\left\{{\bf y}_{j,k}\right\}_{j,k=1}^{N_{T}+1,N_{B}}$, we use the following loss function

where

with $[\cdot]_{k=1}^{N_{B}}$ denoting averaging over a given batch, and $[\cdot]_{j=j_{1}}^{j_{2}}$ denoting averaging over the $j_{1}$ to $j_{2}$ timesteps. We also emphasize here that the matrices $\bar{{\bf K}}_{M}$, ${\bf K}_{a}$, and $\mathbf{\Psi}_{-}$ are defined globally over a given batch of trajectories. This is in contrast to the methods in [6] or [4] which computes equivalent matrices over each trajectory. Arguably then, our approach is learning more global representations of the Koopman operator, the consequences of which we will explore in future work, though see the Appendix for a brief exploration on this point.

Given the success of the DLHDMD in reconstructing and forecasting dynamics along a strange attractor, especially when compared to the relative failure of trying to do the same using just the HDMD or DLDMD alone, it is of further interest to try to assess exactly what role the auto-encoder plays in improving the outcome of the HDMD. While we can certainly point to the performance of the components of the loss function $\mathcal{L}_{tot}$ to explain the impact of the encoder, this does not provide us with any more explanatory power. In [6], it was empirically shown that the role of the encoder was to generally transform time series to nearly monochromatic periodic signals, which is to say, the effect of encoding was to generate far more localized Fourier spectral representations of the original time series. This does not turn out to be the case though for the DLHDMD. Instead, inspired both by the evolving understanding of how mutual information better explains results in dynamical systems [1, 31] and machine learning [16, 17], we assess the impact of the encoder on the DLHDMD by tracking how the information across dimensions and time lags changes in the original and latent variables.

For two random variables ${\bf X}$ and ${\bf Y}$ with joint density $p({\bf X},{\bf Y})$, the mutual information (MI) between them $I({\bf X},{\bf Y})$ is defined to be

where $p({\bf X})$ and $p({\bf Y})$ are the affiliated marginals. One can readily show that $I({\bf X},{\bf Y})\geq 0$ and $I({\bf X},{\bf Y})=0$ if and only if ${\bf X}$ and ${\bf Y}$ are independent. Thus information gives us a stronger metric of statistical coupling between random variables than more traditional tools in time series analysis such as correlation measurements. We also should note here that $I({\bf X},{\bf Y})=I({\bf Y},{\bf X})$, which is to say it is symmetric. We also note that MI is invariant under the action of diffeomorphisms of the variables. Thus we cannot expect to get much use from computing the multidimensional MI of the original and latent variables, thereby allowing for meaningful differences to appear between original and latent variable computations.

Instead, using the $N_{C}=1000$ trajectories in the test data, we define the $m$-step averaged lagged self-information (ALSI) between the $n^{th}$ and $v^{th}$ dimensions $I_{nv}(m)$ to be

We refer to the parameter $m$ as a lag. In words then, after averaging over the ensemble of initial conditions in the test data, we compute the degree to which the signal becomes statistically independent from itself across all of the dimensions along which the dynamics evolve. We emphasize that due to the strong nonlinearities in our dynamics, we compute the lagged information as opposed to the more traditional auto-correlation so as to get a more accurate understanding of the degree of self-dependence across dimension in our dynamics. Further, by measuring the lagged MI across isolated dimensions, we break the invariance of MI with respect to diffeomorphisms.

In this work, we have developed a machine learning enhanced version of the HDMD and DLDMD which we call the DLHDMD. We have shown that its performance is significantly better than just the HDMD, and when comparing against existing results in [6] we see radical improvement over the DLDMD method for the Lorenz-63 system. Likewise, we find that our method is successful across several challenging chaotic dynamical systems varying in dynamical features and size. Thus, we have an accurate parallel approach fitting within the larger framework of Koopman operator based methods. Moreover, we have a method which computes Koopman modes globally. Finally, our analysis of the relative information dynamics across physical dimensions in the original and latent variables provides us a means of understanding the impact of the encoder network on the dynamics in line with modern thinking in machine learning as well as better pointing towards an understanding that the HDMD is enhanced by decreasing the relative statistical dependence across physical dimensions. As explained in detail in the Introduction, there are of course a number of questions that remain to be addressed, and they will certainly be the subject of future research.

C.W. Curtis would like to acknowledge the generous support of the Office of Naval Research and their Summer Research Faculty Program for providing the support of this project. D.J. Alford-Lago acknowledges the support of the Naval Information Warfare Center. E. Bollt was funded in part by the U.S. Army Research Office grant W911NF- 16-1-0081, by the U.S. Naval Research Office, the Defense Advanced Research Projects Agency, the U.S. Air Force Research Office STTR program, and the National Institutes of Health through the CRCNS. We also wish to thank both anonymous reviewers for a number of insightful comments and questions which helped dramatically improve the paper.