Enhancing Dynamic-Mode Decomposition via Multi-Scale Analysis

While finite-dimensional dynamical systems are certainly of interest in their own right, we often wish to study nonlinear partial-differential equations (PDEs) of the general form

where $\mathcal{H}$ is some reasonably chosen Hilbert space. Frequently, the only recourse to solving the above initial-value problem is through numerical simulation, which is often facilitated by expanding the solution $u(x,t)$ along an orthonormal basis of $\mathcal{H}$ so that

This generates a system of equations describing the evolution of the coefficients $\hat{u}_{j}(t)$.

We are therefore motivated to study nonlinear dynamical systems of the generic form

where we suppose $N_{s}\gg 1$. We denote the affiliated flow of this equation as ${\bf a}(t)=\varphi_{t}({\bf a}_{0})$. We define the associated Hilbert space of observables, say $L_{2}\left(\mathbb{C}^{N_{s}},\mathbb{C},\mu\right)$, so that $g\in L_{2}\left(\mathbb{C}^{N_{s}},\mathbb{C},\mu\right)$ if

where $\mu$ is some appropriately chosen measure. A great deal of insight can be gained from looking at the affiliated linear representation of the problem given via the infinite-dimensional Koopman operator $\mathcal{K}^{t}$, where

so that if $g\in L_{2}\left(\mathbb{C}^{N_{s}},\mathbb{C},\mu\right)$, then

Note, throughout the remainder of the paper, we refer to the space of observables as $L_{2}\left(O\right)$.

The power in this approach is that by moving to a linear-operator framework, the affiliated dynamics of the nonlinear system as measured via observables is captured via the eigenvalues of $\mathcal{K}^{t}$. Thus, assuming for the moment that $\mathcal{K}^{t}$ has only discrete spectra, then if we can find a basis of $L_{2}\left(O\right)$ via the Koopman modes $h_{j}$ where

then for any other observable $g$, we should have

However, as one would imagine, determining the modes of the Koopman operator is in general impossible in an analytic way. The Dynamic-Mode Decomposition (DMD) method [2, 1, 7, 3] has been developed which allows for the determination of a finite number of the Koopman modes.

To wit, if we sample the flow ${\bf a}(t;{\bf a}_{0})$ at a discrete set of times $t_{j}$, where $j=0,\cdots,N_{T}$, thereby generating the data set ${\bf a}_{j}=\varphi_{t_{j}}({\bf a}_{0})$, and if we select some set of observables ${\bf g}=\{g_{l}\}_{l=1}^{M}$, then the DMD method approximates $\mathcal{K}^{t}$ by computing the spectra of the finite-dimensional operator $\tilde{\mathcal{K}}_{a}$, which solves the the finite-dimensional least-squares problem

where $G_{-}$ and $G_{+}$ are the $M\times N_{T}$ matrices of the past and future states respectively, i.e.

For example, if we chose the observable $g_{l}({\bf a})$ where

which just finds the approximation of the solution to the PDE of interest at the $l^{th}$ point of the computational mesh, then the Koopman operator provides a linear framework to describe the dynamics of a coordinate which evolves along an otherwise nonlinear flow. Note, this particular choice of observation function is called the canonical observation throughout the remainder of this paper.

The efficacy of this approach though is contingent on the choice of observables one makes, and to date no general strategy for observable choosing strategies has been made. Some success has been had by trying to use kernel based means allowing for the efficient exploration of more complicated, nonlinear observables; see [10]. However, as shown in [8], for infinite dimensional problems, or at least the relatively high-dimensional dynamical systems which approximate them, kernel based techniques can be of very limited utility. By using various reproducing kernel techniques from machine learning, [11], more robust methods for determining effective observables can be found. While effective, the methodology in [11] is computationally intensive requiring the manipulation of large matrices.

In this paper we propose a relatively stream-lined, and in some sense, more physically motivated approach to constructing nonlinear observables $g$. We do this by making use of wavelet based, multi-scale approximations of various function space norms. By doing so, we are able to create a more robust implementation of the DMD algorithm without creating a significantly larger computational burden. We likewise show that the greater degrees of freedom introduced by our approach allow for more optimal decomposition strategies than would be available otherwise. Thus our work should be of use in helping to enhance current DMD based approaches as well as the next way of methodologies under current development [12, 13].

The outline of the paper is as follows. In Section 2 we explain our wavelet based multi-scale generation of observables. In Section 3 we present our results on a model test case. Section 4 collects the conclusions of our work.

We take as our prototypical example of a nonlinear PDE the Nonlinear Schrödinger Equation (NLS) which is

While we cannot simulate this equation on all of $\mathbb{R}$, by taking $x\in[-L,L]$, $L\gg 1$, we can generate a highly accurate numerical approximation to the above equation by using a standard periodic pseudo-spectral approximation in space coupled to a 4th-order Runge-Kutta scheme with integrating factor in time. Thus, we decompose $u(x,t)$ so that

where $\delta_{x}=\frac{2L}{K_{T}}$, $\delta_{k}=\frac{\pi}{L}$, thereby generating a $K_{T}$-dimensional dynamical system describing the evolution of $\hat{u}_{j}(t)$. In general, we expect this system to be Hamiltonian, which implies the Koopman operator is unitary, and thus has strictly imaginary spectra $\lambda_{j}$; see [6, 9, 16].

As an initial condition, we take

where $\tilde{\theta}(k)$ is a continuous random-walk generated from the uniform distribution $U(-\delta_{k}/2,\delta_{k}/2)$, and $\epsilon\geq 0$. Thus for nonzero $\epsilon$, we look at an inital condition in which we have a central peak at $x=0$ and a secondary one at $x=x_{s}$ all on a background of normally distributed perturbations. In the limit as $\epsilon\rightarrow 0^{+}$, we should have

Using the canonical observation in the $\epsilon\rightarrow 0^{+}$ limit, we see that if the time step in the numerical integrator is $\delta t$ then we have

where $D^{(s)}$ and $D^{(t)}$ are $K_{T}\times K_{T}$ and $N_{T}\times N_{T}$ diagonal matrices respectively such that

and $O$ is a $K_{T}\times N_{T}$ matrix of ones, i.e.

Thus, in the vanishing $\epsilon$ limit, we see that the approximation to the Koopman operator found via DMD given by $\tilde{K}_{a}$ is

where $\left(D^{(s)}O\right)^{-P}$ denotes the Moore–Penrose pseudoinverse of the matrix $D^{(s)}O$. If the Singular Value Decomposition (SVD) of $O$ is

we readily see that

where $O_{K_{T}}$ is a $K_{T}\times K_{T}$ matrix of ones. Thus the eigenvalues of $\tilde{K}_{a}$ are $e^{i\delta t}$ and $0$, repeated $K_{T}-1$ times, so ignoring the null-space of $\tilde{K}_{a}$, we find the affiliated Koopman eigenvalue $\lambda=i$. Thus in all future computations we anticipate that the computed spectra should reflect some degree of perturbation away from $e^{i\delta t}$.

Following the approach in [2], the practical implementation of the DMD starts with a SVD of the past matrix $G_{-}$ so that Equation (1) becomes

Motivated by the previous example though, we should anticipate the possibility that $\Sigma_{-}$ could be of relatively small rank, which is to say that the conditioning inherent in trying to compute the Moore–Penrose pseudoinverse of $G_{-}$ could compromise the accuracy of any results one might try to compute. Thus, letting ${\bf d}_{\Sigma_{-}}$ be the descending ordered vector of the non-negative primary diagonal entries of $\Sigma_{-}$, we define the index $D_{T_{l}}$ to be the largest index of ${\bf d}_{\Sigma_{-}}$ such that

Thus the non-negative integer tol sets the relative accuracy of any computations and allows for the omission of spurious modes computed by the DMD. This in turn defines a reduced SVD of $G_{-}$ so that

and thus $\tilde{\mathcal{K}}_{a}$ becomes the $M\times M$ matrix

Finding the affiliated eigenvalues and eigenvectors of $\tilde{\mathcal{K}}_{a}$, i.e.

then gives us the affiliated approximation formula at $t_{n}=n\delta t$

where the coefficients reflect the initial values of the observables ${\bf g}({\bf a}(0))$.

With regards to our choice of Besov norms and corresponding multiscale decompositions, we choose to add to the list of canonical observables the decompositions of the quantities

This corresponds to the first conserved quantity of the NLSE, i.e. the energy $\mathcal{E}(u)=\left|\left|u\right|\right|^{2}_{2}$, while the next two quantities reflect the nonlinearity and dispersion of the NLSE as reflected in the Hamiltonian $\mathcal{H}(u)$ where

This adds an overall $3(N_{lvl}+1)$ observables to the $K_{T}$ canonical observables corresponding to finding $u(x_{l},t)$, which we collectively denote as $\tilde{{\bf g}}$. Then the extended observables ${\bf g}$ are given by

We denote the Besov-norm enhanced, multiscale DMD as the MDMD. Note, while we have relied on knowing the model to make our decisions, one would very naturally choose these two Besov norms in so far as their affiliated variational derivatives naturally lead one to cubic nonlinearities and Laplacian based dispersion terms, which from a Ginzburg–Landau analysis are very natural leading order effects which emerge across a wide range of models.

At this point, with the choice of $T_{l}$ and the number of levels across which to perform a wavelet decomposition of the data, it is a nontrivial question to determine parameter choices for optimal performance of the DMD decomposition. To measure this performance, we use two criteria. The first is accuracy of the reconstruction of the given data. The second is that we know the underlying model is a Hamiltonian dynamical system, and thus the affiliated spectra $\mu_{j}$ should lie entirely on the imaginary unit circle. We thus define an affiliated error $E(T_{l},N_{v})$ where

where the reconstruction error $E_{rc}$ is given by

and the spectral error $E_{sp}$ is given by

where $\tilde{{\bf g}}({\bf a}(t_{f}))$ are the canonical observables not including those computed by Besov Norms measured at the final sampled time point, $\tilde{{\bf g}}_{DMD}(t_{f})$ is the affiliated DMD generated approximation of those same observables computed at the final sampled time point, and $\tilde{w}\in[0,2]$ is a relative weighting parameter between the two criteria. Thus for $0\leq\tilde{w}<1$, we favor reconstruction accuracy, while for $1<\tilde{w}\leq 2$, we favor adherance of the spectra $\mu_{j}$ to the unit circle.

Throughout the simulations below then, for given choice of $\tilde{w}$, we use an optimization strategy of sweeping through different choices of $N_{lvl}$ and $T_{l}$ until a minimal error is found. The parameters are swept out over the intervals

Since we examine the role of noise, we also do ensemble runs of where the number of members $N_{e}=16$, which was the largest number of runs which allowed for practical run times of the simulations. Each simulation is run with $L=32$ and $K_{T}=256$ for a time step of $\delta t=.1$ and up to a time of $t_{f}=30$ so that $N_{T}>K_{T}+3(N_{lvl}+1)$ since in this case we have that $1\leq N_{lvl}\leq 7$. The secondary peak in the initial condition is centered at $x_{s}=5$.

We first look at weighting the error norm in favor of approximation accuracy by choosing $\tilde{w}=.01$ for relatively weak perturbations, $\epsilon=.05$, and stronger perturbations in which $\epsilon=.25$. These results are seen respectively in Figures 1 and 2. For the relatively weak perturbation case, we see the reconstruction error is at times almost an order of magnitude better for the MDMD, while only making a modest sacrifice in the spectral error; see Figures 1 (a) and (b). Moreover, looking at Figures 1 (a) and (b), we see that the MDMD method is able to create a more uniform error in the approximation across the 16 different simulations. As we see in Figure 1 (c), this is done in part by allowing for a higher value of $T_{l}$, thereby allowing for the use of more modes emerging from the DMD decomposition, and in part by using $N_{lvl}$ to adjust to varying noise profiles in an adaptive way as seen in Figure 1 (d).

For the larger perturbation when $\epsilon=.25$, as seen in Figure 2 (a), the classic DMD approach often is incapable of maintaining reconstruction accuracy. While the multiscale approach also can suffer, it is far more reliable than its traditional counterpart. As seen in Figure 2 (b), this better reconstruction accuracy comes at a modest loss of spectral accuracy relative to the traditional DMD, though it should be noted that the reconstruction errors of the traditional DMD are $\mathcal{O}(1)$, thus bringing into question the value of it maintaing spectral accuracy. Again, as seen in Figures 2 (b)-(d), the error profile using Besov-norms is slightly better and more uniform, the Besov approach uses more DMD modes, and $N_{lvl}$ is chosen in a highly adaptive manner.

By instead looking at prioritizing the spectral accuracy of the DMD decomposition by setting $\tilde{w}=2$, we see for both levels of perturbations that the Besov based approach provides strong improvements over standard DMD. While for both levels of the perturbation the MDMD approach provides both better and more uniform levels of reconstruction error, again for the largest perturbation, the standard DMD approach completely breaks down while the multiscale approach is markedly more robust; see Figures 3 (a) and 4 (a). Again this comes at a modest loss of spectral accuracy; see Figures 3 (b) and 4 (b). Thus, the MDMD based approach does not compromise approximation fidelity to the extent that the traditional DMD method does through the enhanced adaptability and better capturing of the dynamics facilitated by the Besov norms which help capture the key nonlinear features of the flow.