Physics-Informed Probabilistic Learning of Linear Embeddings of Non-linear Dynamics With Guaranteed Stability

TL;DR

This paper introduces a physics-informed probabilistic deep learning framework for stable, measure-theoretic learning of Koopman operators, enabling accurate, uncertainty-aware modeling of nonlinear dynamical systems.

Contribution

It formalizes continuous-time Koopman operator learning with neural networks, ensuring stability and incorporating uncertainty quantification via hierarchical Bayesian inference.

Findings

Successfully applied to polynomial, oscillator, and fluid flow systems.

Achieved stable and interpretable linear embeddings of nonlinear dynamics.

Quantified uncertainties in dynamic predictions.

Abstract

The Koopman operator has emerged as a powerful tool for the analysis of nonlinear dynamical systems as it provides coordinate transformations to globally linearize the dynamics. While recent deep learning approaches have been useful in extracting the Koopman operator from a data-driven perspective, several challenges remain. In this work, we formalize the problem of learning the continuous-time Koopman operator with deep neural networks in a measure-theoretic framework. Our approach induces two types of models: differential and recurrent form, the choice of which depends on the availability of the governing equations and data. We then enforce a structural parameterization that renders the realization of the Koopman operator provably stable. A new autoencoder architecture is constructed, such that only the residual of the dynamic mode decomposition is learned. Finally, we employ mean-field variational inference (MFVI) on the aforementioned framework in a hierarchical Bayesian setting to quantify uncertainties in the characterization and prediction of the dynamics of observables. The framework is evaluated on a simple polynomial system, the Duffing oscillator, and an unstable cylinder wake flow with noisy measurements.

Tables

Table 1. Table 1: Comparison with other approaches in the literature

Previous works	continuous /discrete	nonlinear reconstruction	continuous spectrum	uncertainty	stability
Yeung et al. [81]	discrete	✗	✗	✗	✗
Li et al. [43]	discrete	✗	✗	✗	✗
Takeishi et al. [73]	discrete	✓	✗	✗	✗
Otto and Rowley [60]	discrete	✓	✗	✗	✗
Lusch et al. [48]	discrete	✓	✓	✗	✗
Morton et al. [58]	discrete	✓	✗	✓	✗
Erichson et al. [17]	discrete	✓	✗	✗	✓
Our framework	continuous	✓	✗	✓	✓

Table 2. Table 2: Hyperparameters of differential form model for 2D fixed point attractor.

layer structure	optimizer	learning rate	total epoch	batch size
2-8-16-16-8-2-8-16-16-8-2	Adam	1e-4	20000	128

Table 3. Table 3: Hyperparameters of differential form model for unforced Duffing system.

layer structure	optimizer	learning rate	total epoch	batch size
2-16-16-24-16-16-3-16-16-24-16-16-2	Adam	1e-3	20000	128

Table 4. Table 4: Hyperparameters of recurrent form model for flow past cylinder at R ​ e D = 100 𝑅 subscript 𝑒 𝐷 100 Re_{D}=100 .

layer structure	optimizer	learning rate	total epoch	batch size	$T$
50-100-50-20-20-20-50-100-50	Adam	1e-3	20000	64	100

Equations

$$\dot{\mathbf{x}}=\mathbf{F}(\mathbf{x}),\mathbf{x}\in \mathcal{M}\subset {\mathbb{R}}^N.$$

$${\mathcal{F}}_D=\text{span}\{{\varphi }_1,\dots ,{\varphi }_D\}\subset \mathcal{F},$$

$$\mathbf{\Phi }(x)=\left[\begin{array}{cccc}{\varphi }_1(x)& {\varphi }_2(x)& \dots & {\varphi }_D(x)\end{array}\right]\in {\mathcal{F}}^D.$$

$$\dot{\mathbf{\Phi }}\triangleq d\mathbf{\Phi }/dt=\mathbf{\Phi K}.$$

$$\mathbf{F}\cdot {\nabla }_{\mathbf{x}}\mathbf{\Phi }=\mathbf{\Phi K}.$$

$${\eta }_l={\sigma }_l({\eta }_{l-1}{W}_l+b_l),$$

$$g({\eta }_0;{\mathbf{W}}_g)={\eta }_{L-1}{W}_L+b_L,$$

$${\Vert \varphi \Vert }_{\mathcal{F}}\triangleq \sqrt{{\int }_{\mathcal{M}}|\varphi {|}^{2}d\mu }<\infty .$$

$${\Vert \mathbf{\Phi }\Vert }_{{\mathcal{F}}^n}\triangleq \sqrt{{\int }_{\mathcal{M}}\sum _{j=1}^n|{\varphi }_j{|}^{2}d\mu }.$$

$$\mathcal{J}[\mathbf{\Phi }]=\underset{\psi \in \{{\varphi }_1,\dots ,{\varphi }_D\}}{\max }\underset{f\in {\mathcal{F}}_D}{\min }{\Vert f-\mathcal{K}\psi \Vert }_{\mathcal{F}}^2,$$

$${\mathbf{\Phi }}^{\ast }=\underset{\begin{array}{c}\mathbf{\Phi }\in {\mathcal{F}}^D,\exists \mathbf{\Psi }:{\mathbb{R}}^D\mapsto \mathcal{M}\\ \mathbf{\Psi }\circ \mathbf{\Phi }=\mathcal{I}\end{array}}{\mathrm{argmin}}\mathcal{J}[\mathbf{\Phi }].$$

$$0\le \underset{\begin{array}{c}\mathbf{\Phi }\in {\mathcal{F}}^D,\exists \mathbf{\Psi }:{\mathbb{R}}^D\mapsto \mathcal{M}\\ \mathbf{\Psi }\circ \mathbf{\Phi }=\mathcal{I}\end{array}}{\min }\mathcal{J}[\mathbf{\Phi }]\le \underset{\begin{array}{c}{\mathbf{W}}_{\mathbf{\Phi }},\exists \mathbf{\Psi }\in C({\mathbb{R}}^D,{\mathbb{R}}^N)\\ \mathbf{\Psi }\circ \mathbf{\Phi }(\cdot ;{\mathbf{W}}_{\mathbf{\Phi }})=\mathcal{I}\end{array}}{\min }\mathcal{J}[\mathbf{\Phi }(\cdot ;{\mathbf{W}}_{\mathbf{\Phi }})].$$

$${\mathbf{W}}_{\mathbf{\Phi }}^{\ast }$$

$${\hat{\mathbf{W}}}_{\mathbf{\Phi }},{\hat{\mathbf{W}}}_{\mathbf{\Psi }}$$

$$\mathcal{R}[\mathbf{\Phi },\mathbf{\Psi }]={\Vert \mathbf{\Psi }\circ \mathbf{\Phi }-\mathcal{I}\Vert }_{{\mathcal{F}}^N}^2,$$

$$\mathcal{J}[\mathbf{\Phi }]\le \hat{\mathcal{J}}[\mathbf{\Phi }]$$

$$\tilde{\mathcal{J}}[\mathbf{\Phi },\mathbf{K}]$$

$${\hat{\mathbf{W}}}_{\mathbf{\Phi }},{\hat{\mathbf{W}}}_{\mathbf{\Psi }},\hat{\mathbf{K}}=\underset{\begin{array}{c}{\mathbf{W}}_{\mathbf{\Phi }},{\mathbf{W}}_{\mathbf{\Psi }},\mathbf{K}\end{array}}{\mathrm{argmin}}\tilde{\mathcal{J}}[\mathbf{\Phi }(\cdot ;{\mathbf{W}}_{\mathbf{\Phi }}),\mathbf{K}]+\mathcal{R}[\mathbf{\Phi }(\cdot ;{\mathbf{W}}_{\mathbf{\Phi }}),\mathbf{\Psi }(\cdot ;{\mathbf{W}}_{\mathbf{\Psi }})].$$

$${\mathcal{J}}_r[\mathbf{\Phi }]=\underset{\psi \in \{{\varphi }_1,\dots ,{\varphi }_D\}}{\max }\underset{\mathbf{K}\in {\mathbb{R}}^{D\times D}}{\min }{\Vert \mathbf{\Phi }{e}^{t\mathbf{K}}{c}_{\psi }-{\mathcal{K}}_t\psi \Vert }_{\mathcal{G}}^2,$$

$$\tilde{\mathcal{P}}[\mathbf{\Phi },\mathbf{\Psi },\mathbf{K}]={\Vert \mathbf{\Psi }\circ \mathbf{\Phi }{e}^{t\mathbf{K}}-{\mathcal{K}}_t\mathcal{I}\Vert }_{{\mathcal{G}}^N}^2.$$

$${\tilde{\mathcal{J}}}_r[\mathbf{\Phi }(\cdot ;{\mathbf{W}}_{\mathbf{\Phi }}),\mathbf{K}]={\Vert \mathbf{\Phi }{e}^{t\mathbf{K}}-{\mathcal{K}}_t\mathbf{\Phi }\Vert }_{{\mathcal{G}}^D}^2,$$

$${\hat{\mathbf{W}}}_{\mathbf{\Phi }},{\hat{\mathbf{W}}}_{\mathbf{\Psi }},\hat{\mathbf{K}}=\underset{\begin{array}{c}{\mathbf{W}}_{\mathbf{\Phi }},{\mathbf{W}}_{\mathbf{\Psi }},\mathbf{K}\end{array}}{\mathrm{argmin}}{\tilde{\mathcal{J}}}_r[\mathbf{\Phi }(\cdot ;{\mathbf{W}}_{\mathbf{\Phi }}),\mathbf{K}]+\hat{\mathcal{P}}[\mathbf{\Phi }(\cdot ;{\mathbf{W}}_{\mathbf{\Phi }}),\mathbf{\Psi }(\cdot ;{\mathbf{W}}_{\mathbf{\Psi }}),\mathbf{K}].$$

$$\mathcal{R}[\mathbf{\Phi },\mathbf{\Psi }]={\Vert {\mathcal{K}}_t(\mathbf{\Psi }\circ \mathbf{\Phi }-\mathcal{I})\Vert }_{{\mathcal{G}}^N}^2,$$

$${\Vert \mathbf{\Psi }\circ \mathbf{\Phi }{e}^{t\mathbf{K}}-{\mathcal{K}}_t\mathcal{I}\Vert }_{{\mathcal{G}}^N}$$

$$\le L_{\mathbf{\Psi }}{\Vert \mathbf{\Phi }{e}^{t\mathbf{K}}-{\mathcal{K}}_t\mathbf{\Phi }\Vert }_{{\mathcal{G}}^D}+{\Vert {\mathcal{K}}_t(\mathbf{\Psi }\circ \mathbf{\Phi }-\mathcal{I})\Vert }_{{\mathcal{G}}^N},$$

$${\hat{\mu }}_M=\frac{1}{M}\sum _{m=1}^M{\delta }_{{\mathbf{x}}_m},$$

$${\hat{\mathbf{W}}}_{\mathbf{\Phi }},{\hat{\mathbf{W}}}_{\mathbf{\Psi }},\hat{\mathbf{K}}=\underset{M\to \infty }{\lim }{\hat{\mathbf{W}}}_{\mathbf{\Phi },M},{\hat{\mathbf{W}}}_{\mathbf{\Psi },M},{\hat{\mathbf{K}}}_M$$

$$=\underset{M\to \infty }{\lim }\underset{\begin{array}{c}{\mathbf{W}}_{\mathbf{\Phi }},{\mathbf{W}}_{\mathbf{\Psi }},\mathbf{K}\end{array}}{\mathrm{argmin}}{\tilde{\mathcal{J}}}_M[\mathbf{\Phi }(\cdot ;{\mathbf{W}}_{\mathbf{\Phi }}),\mathbf{K}]+{\mathcal{R}}_M[\mathbf{\Phi }(\cdot ;{\mathbf{W}}_{\mathbf{\Phi }}),\mathbf{\Psi }(\cdot ;{\mathbf{W}}_{\mathbf{\Psi }})]$$

$$=\underset{M\to \infty }{\lim }\underset{\begin{array}{c}{\mathbf{W}}_{\mathbf{\Phi }},{\mathbf{W}}_{\mathbf{\Psi }},\mathbf{K}\end{array}}{\mathrm{argmin}}{\Vert \mathbf{\Phi }(\cdot ;{\mathbf{W}}_{\mathbf{\Phi }})\mathbf{K}-\mathcal{K}\mathbf{\Phi }\Vert }_{{\hat{\mathcal{F}}}_M^D}^2+{\Vert \mathbf{\Psi }(\cdot ;{\mathbf{W}}_{\mathbf{\Psi }})\circ \mathbf{\Phi }(\cdot ;{\mathbf{W}}_{\mathbf{\Phi }})-\mathcal{I}\Vert }_{{\hat{\mathcal{F}}}_M^N}^2$$

$$\displaystyle=\lim_{M\rightarrow\infty}\operatorname*{\arg\!\min}_{\begin{subarray}{c}\mathbf{W}_{\mathbf{\Phi}},\mathbf{W}_{\mathbf{\Psi}},\mathbf{K}\end{subarray}}\frac{1}{M}\sum_{m=1}^{M}\bigg{(}\lVert\mathbf{\Phi}(\mathbf{x}_{m};\mathbf{W}_{\mathbf{\Phi}})\mathbf{K}-\mathbf{F}\cdot\nabla_{\mathbf{x}}\mathbf{\Phi}(\mathbf{x}_{m})\rVert^{2}$$

$$\displaystyle\qquad\qquad\qquad\qquad+\left\lVert\mathbf{\Psi}(\mathbf{\Phi}(\mathbf{x}_{m};\mathbf{W}_{\mathbf{\Phi}});\mathbf{W}_{{\mathbf{\Psi}}})-\mathbf{x}_{m}\right\rVert^{2}\bigg{)},$$

$$\nu ({\mathcal{E}}_{\mathcal{M}}\times {\mathcal{E}}_{\mathcal{T}})={\int }_{{\mathcal{E}}_{\mathcal{M}}}\kappa ({\mathcal{E}}_{\mathcal{T}},\mathbf{x})d\mu (\mathbf{x}),$$

$${\hat{\nu }}_{M,\hat{T}}=\frac{1}{M}\sum _{m=1}^M{\delta }_{{\mathbf{x}}_m}\left(\frac{1}{T_m}\sum _{j=1}^{T_m}{\delta }_{t_{m,j}}\right),$$

$${\hat{\mathbf{W}}}_{\mathbf{\Phi }},{\hat{\mathbf{W}}}_{\mathbf{\Psi }},\hat{\mathbf{K}}=\underset{M,T\to \infty }{\lim }{\hat{\mathbf{W}}}_{\mathbf{\Phi },M,\hat{T}},{\hat{\mathbf{W}}}_{\mathbf{\Psi },M,\hat{T}},{\hat{\mathbf{K}}}_{M,\hat{T}}$$

Full text

{AMS}

37M10, 37M25, 47B33, 62F15

1 Introduction

Nonlinear dynamical systems are prevalent in physical, mathematical and engineering problems [78]. Nonlinearity gives rise to rich phenomena, and is much more challenging to characterize in comparison to linear systems for which comprehensive techniques have been developed for decades [32].

Consider the autonomous system $\dot{\mathbf{x}}=\mathbf{F}(\mathbf{x})$, $\mathbf{x}\in\mathcal{M}\subset\mathbb{R}^{N}$, where $\mathcal{M}$ is a smooth manifold in the state space, $\mathbf{F}:\mathcal{M}\mapsto T\mathcal{M}$ is a vector-valued smooth function and $T\mathcal{M}$ is the tangent bundle, i.e., $\forall p\in\mathcal{M},\mathbf{F}(p)\in T_{p}\mathcal{M}$. Instead of a geometric viewpoint [23, 28], Koopman [37] offered an operator-theoretic perspective by describing the “dynamics of observables” on the measure space $(\mathcal{M},\Sigma,\mu)$ via the Koopman operator $\mathcal{K}_{t}:\mathcal{F}\mapsto\mathcal{F}$, such that for an observable on the manifold, $\forall f\in\mathcal{F}$, $t\in\mathbb{R}_{\geq 0}$, $f:\mathcal{M}\mapsto\mathbb{C}$, $\mathcal{K}_{t}f\triangleq f\circ S_{t}$, where $S_{t}(\mathbf{x}_{0})\triangleq S(\mathbf{x}_{0},t)$ is the flow map that takes the initial condition $x_{0}$ and advances it by time $t$ by solving the initial value problem for the aforementioned nonlinear dynamics, with $\mathcal{F}=L^{2}(\mathcal{M},\mu)$, where $\mu$ is some positive measure on $(\mathcal{M},\Sigma)$.

The semigroup of Koopman operators $\{\mathcal{K}_{t}\}_{t\in\mathbb{R}^{+}}$ is generated by $\mathcal{K}:\mathcal{D}(\mathcal{K})\mapsto\mathcal{F}$, $\mathcal{K}f\triangleq\lim_{t\rightarrow 0}(\mathcal{K}_{t}f-f)/t$ where $\mathcal{D}(\mathcal{K})$ is the domain in which the aforementioned limit is well-defined and $\overline{\mathcal{D}(\mathcal{K})}=\mathcal{F}$. This operator is linear, but defined on $\mathcal{F}$, which makes it inherently infinite-dimensional. Physically, $\mathcal{K}_{t}$ governs the temporal evolution of any observable in $\mathcal{F}$.

Specifically, for measure-preserving systems, e.g., Hamiltonian system or the dynamics on an attractor of the Naiver–Stokes equations, one can guarantee a spectral decomposition of any observable $f\in\mathcal{F}$, or the existence of a spectral resolution of the Koopman operator restricted on the attractor, i.e., $\mathcal{K}_{t}|_{\mathcal{A}}$, where $\mathcal{A}\subset\mathcal{M}$ is the attractor. Such a spectral decomposition [54] is a sum of three essentials: temporal averaging of the observable, the contribution from the point spectrum, which is almost periodic in time, and that from the continuous spectrum, which is chaotic [36]. For a more comprehensive discussion, readers are referred to the excellent review by Budišić et al. [11]. It should be stressed that for measure-preserving systems, the Koopman operator is not only well-defined on $L^{2}(\mathcal{M},\mu)$ but also can be shown to be unitary, which ensures properties such as simple eigenvalues and the existence of the aforementioned spectral resolution [11]. Throughout this work, we employ the simple but practical assumption [39] $\mathcal{F}=L^{2}(\mathcal{M},\mu)$, where $\mu$ is some positive measure with support equal to $\mathcal{M}$. Note that this implies that the Koopman operator is well-defined on $\mathcal{F}$.

The appeal of the Koopman operator lies in the possibility of globally linearizing the nonlinear dynamics of certain observables, i.e., the quantities of interest. In practice, we are mostly interested in the subspace $\mathcal{F}_{D}$, which consists of a finite number of Koopman eigenfunctions, namely the finite-dimensional Koopman invariant subspace, where $D\in\mathbb{N}$ is the dimension of the subspace. Naturally, for a given observable $f$, we are interested in the minimal $\mathcal{F}_{D}$, such that $f\in\mathcal{F}_{D}$. Koopman analysis entails the search for the Koopman eigenfunctions that span $\mathcal{F}_{D}$ and its associated eigenvalues. While we assume that the eigenvalues are simple, it is possible to extend to the generalized case [11]. The identity mapping, i.e., $f=\mathbf{x}$ is of central interest since it corresponds to a linear system that is topologically conjugate to the nonlinear dynamical system. This essentially generalizes the Hartman-Grobman theorem to the entire basin of attraction of the equilibrium point or periodic orbit [41]. From the viewpoint of understanding the behavior of the dynamical system in state space, the level sets of Koopman eigenfunctions form the invariant partition in the state space [10], which can help study mixing. Further, the ergodic partition can be analyzed with Koopman eigenfunctions [57].

Koopman analysis has also gained increased attention in different communities. Koopman eigenfunctions have been leveraged in nonlinear optimal control [31, 56, 8, 66, 38]. In modal analysis, the equivalence between the dynamic mode decomposition and the approximation of Koopman modes was recognized by Rowley et al. [69]. In a reduced order modeling context, the operator-theoretic viewpoint of the Koopman operator is critical to the formulation of Mori-Zwanzig (M-Z) formalism [15], which has recently found utility in addressing closure and stabilization [21, 63, 64] of multiscale problems.

Despite its appealing properties, the Koopman decomposition cannot be pursued in its original form described above for practical scientific applications, because the operator is defined on the infinite dimensional Hilbert space. To accommodate practical computation of the Koopman operator, we consider two simplifications:

• •

Only a finite dimensional space invariant to $\mathcal{K}_{t}$ is of interest. This excludes the possibility of dealing with a chaotic system since it is impossible for a finite dimensional linear system to be topologically conjugate to a chaotic system [11]. However, one could assign a state-dependent Koopman operator to work around this limitation [48].

• •

Nonlinear reconstruction of observables from Koopman eigenfunctions. The original concept of the Koopman operator studies the evolution of any observable in the Koopman invariant subspace. Thus, the observables of interest can be linearly reconstructed from the Koopman eigenfunctions. Regarding the latter Linear reconstruction is desirable especially in the content of control [31]. Li et al. [43] considered augmenting the Koopman invariant subspace with neural network-trained observables together with the system state $\mathbf{x}$ to force linear reconstruction. However, for systems with multiple fixed points, there does not exist a finite-dimensional Koopman invariant subspace that also spans the state globally111In this paper, the term “global embedding” is used to imply non-locality in phase space. Note that the Hartman-Grobman theorem [4] establishes a topological conjugacy to a linear system with the same eigenvalues in a small neighborhood of the fixed point. due to the impossibility of establishing such a topological conjugacy [8] with the linear system. Therefore, attention naturally moves to the introduction of nonlinear reconstruction that is more expressive than linear reconstruction, together with extra modes that indicate different basins of each attractor [73, 60]. Unfortunately, nonlinear reconstruction is equivalent to removing Koopman modes for the observables of interest [60]. In other words, this is equivalent to removing the constraint that the system state lies in the finite dimensional Koopman invariant subspace. Although for systems with a single attractor, one could still retain a linear reconstruction representation, which might lead to a decrease in accuracy [60].

To determine the Koopman eigenfunctions, we consider addressing the problem from a data-driven perspective. A common approximation is to assume that the Koopman invariant subspace $\mathcal{F}_{D}$ is spanned by a finite dictionary of functions such as monomials, and then minimize the empirical risk of the residual that comes from the imperfect $\mathcal{F}_{D}$. This is the general idea behind Extended Dynamic Mode Decomposition (EDMD) [79] and Kernel-based Dynamic Mode Decomposition (KDMD) [80] where an implicitly-defined infinite dimensional Koopman invariant subspace can be computationally affordable if the kernel satisfies Mercer’s condition [29]. Note that EDMD is essentially a least-squares regression of the (linear) action of the Koopman operator for the features in the dictionary, with samples drawn from some measure $\mu$ with fixed features, which can be proved as a $L^{2}(\mathcal{M},\mu)$ orthogonal projection in the limit of large independent identically distributed data [35, 79].

By framing the data-derived Koopman operator in the Hilbert space endowed with a proper measure, one can prove optimality in the asymptotic sense for EDMD [39] and Hankel-DMD [2]. As a side note, if the dynamics on the attractor is of central interest, it has been shown that Hankel-DMD with appropriate delay embedding and sampling rate can recover the entire system, even without sampling the full attractor [62]. Further, one can establish a duality between the measure-preserving dynamics on the attractor with a stationary stochastic process, which reflects a close relationship between the estimation of a continuous spectrum from the trajectories, and the extraction of the power spectral density of stochastic signals [3]. Other algorithms include generalized Laplace analysis (GLA) [11, 50], and the Ulam Galerkin method [18]. If the representation of the Koopman eigenfunctions is sparse in the pre-defined dictionary in EDMD, then one can leverage sparse regression techniques [9] to select observables in an iterative fashion [8]. In this work, we do not address continuous spectra. One way of addressing continuous spectra is to use [48] an auxiliary sub-network to characterize the eigenvalues as state-dependent.

However, it has been shown that traditional EDMD/KDMD method is prone to overfitting and is hard to interpret [60]. For instance, the number of approximated Koopman eigenfunctions scales with either the number of features in the dictionary or the number of training snapshots even though only a few of them are genuine Koopman eigenfunctions. Recently, there have been several alternative attempts to leverage deep learning architectures [60, 43, 73, 77, 81], to extract Koopman decompositions. Yeung et al. [81] used feedforward neural networks to learn the dictionary, but the reconstruction loss was found to be non-optimal. Li et al. [43] enforced several non-trainable functions, e.g., components of the system state, in the Koopman observables to ensure an accurate reconstruction but one that could be inefficient in terms of obtaining a finite-dimensional Koopman observable subspace [60]. Further, Takeishi et al. [73] utilized linear time delay embedding in the feedforward neural network framework to construct Koopman observables with nonlinear reconstruction which is critical for partially observed systems [61]. Lusch et al. [48] further extended the deep learning framework to chaotic systems. Otto and Rowley [60] considered a recurrent loss for better performance on long time prediction on trajectories that transit to the attractor. Recently, Morton et al. [58] addressed the uncertainty in a deep learning model with a focus on control. The benefit of formulating the search for the Koopman operator in an optimization setting enables the enforcement of stability. For example, it is also feasible to constrain eigenvalues in optimized DMD [14, 5]. Specifically in the neural network context, Erichson et al. [17] considered a stability promoting loss to encourage Lyapunov stability of the dynamical system.

A unified approach towards uncertainty quantification, stabilization, and incorporation of physics information is lacking, particularly in a continuous setting. This motivates us to establish a probabilistic stabilized deep learning framework specifically to learn the Koopman decomposition for a continuous dynamical system. We employ automatic differentiation variational inference (ADVI) [40] to quantify parametric uncertainty in deep neural networks, and structural parameterization to enforce stability of the Koopman operator extracted from the deep neural networks. A broad comparison of the present work with other approaches is shown in Table 1.

The outline of this paper is as follows: In Section 2, our framework of deep learning of Koopman operators of continuous dynamics, together with structured parameterization that enforces stabilization is presented. In Section 3, Bayesian deep learning, specifically the employment of variational inference is introduced. Results and discussion on several problems are given in Section 4. Conclusions are drawn in Section 5.

2 Data-driven framework to learn continuous-time Koopman decompositions

In contrast to previous approaches [60, 48, 81, 43], our framework pursues continuous-time Koopman decompositions. The continuous formulation is more amenable to posit desired constraints and contend with non-uniform sampling [60, 5], which is frequently encountered in experiments or temporal multi-scale data. To begin with, recall the general form of autonomous continuous nonlinear dynamical systems,

[TABLE]

We seek a finite dimensional Koopman invariant subspace $\mathcal{F}_{D}$ with $D$ linearly independent smooth observation functions, defined as

[TABLE]

where $\phi_{i}\in C^{1}(\mathcal{M},\mathbb{R})$ and $\phi_{i}\in\mathcal{F}$, $i=1,\ldots,D$. Correspondingly, the observation vector $\mathbf{\Phi}$ is defined as

[TABLE]

Based on the aforementioned structure, we require the following two conditions. First, $\mathbf{\Phi}(x)$ evolves linearly in time, i.e., $\exists$ $\mathbf{K}\in\mathbb{R}^{D\times D}$ s.t.

[TABLE]

By the chain rule, the relationship between Eq. 4 and Eq. 1 is

[TABLE]

Second, $\exists$ $\mathbf{\Psi}:\mathbb{R}^{D}\mapsto\mathcal{M}$, s.t. $\mathbf{\Psi}\circ\mathbf{\Phi}=\mathcal{I}$, $\mathcal{I}:\mathcal{M}\mapsto\mathcal{M}$ is the identity map. Therefore, we can recover the state $\mathbf{x}$ from $\mathbf{\Phi}$.

The goal is to find a minimal set of basis functions that spans a Koopman invariant subspace. There are several existing methods for approximation:

• •

Dynamic mode decomposition (DMD): observation functions are linear transformations of the state, usually POD modes for robustness purposes [70].

• •

Extended/Kernel DMD: observation functions are pre-specified, either explicit polynomials [79] or implicitly by the kernel [80].

• •

Deep learning Koopman/Neural Network DMD [60, 48] searching a set of fixed size for the nonlinear observations by artificial neural networks.

In this work, we specifically focus on using neural networks.

2.1 Neural network

The basic building block we use is a feedforward neural network (FNN). A typical FNN with $L$ layers is the mapping $g(\cdot;\mathbf{W}_{g}):\mathbb{R}^{n}\mapsto\mathbb{R}^{m}$ such that

[TABLE]

for $l=1,\dots,L-1$, where $\eta_{0}$ stands for the input of the neural network $x$, $\eta_{l}\in\mathbb{R}^{n_{l}}$ is the number of hidden units in the layer $l$, $\sigma_{l}$ is activation function of layer $l$. Note that the last layer is linear, i.e., $\sigma_{L}(x)=x$:

[TABLE]

where $\mathbf{W}_{g}=\{W_{1},b_{1},\ldots,W_{L},b_{L}\}$, $W_{l}\in\mathbb{R}^{n_{l-1}\times n_{l}}$, $b_{l}\in\mathbb{R}^{n_{l}}$, for $l=1,\ldots,L$. Note that $n_{0}=n\in\mathbb{N}$, $n_{L}=m\in\mathbb{N}$. Such a mapping, given an arbitrary number of hidden units, even with a single hidden layer, has been shown to be a universal approximator in the $L^{p}$ sense $1\leq p<\infty$, as long as the activation function is not polynomial almost everywhere [42]. Throughout this work, we use the Swish [67] activation function, which is continuously differentiable and found to achieve strong empirical performance over many variants of ReLU and ELU on typical deep learning tasks.

2.2 Problem formulation

We first define the ideal problem of learning the Koopman decomposition as a constrained variational problem, and incorporate assumptions to make it tractable step by step. Then, we introduce the effect of data as an empirical measure into the optimization in the function space. Specifically, we propose two slightly different formulations based on varying requirements.

Differential form: for low-dimensional systems when the governing equations, i.e., Eq. 5 are known, they can be leveraged without the trajectory data. 2. 2.

Recurrent form: for high-dimensional systems where only discrete trajectory data is obtained, which implies the access to the action of $\mathcal{K}_{t}$ over discrete $t$.

Recall that we are searching the Koopman operator defined on $\mathcal{F}=L^{2}(\mathcal{M},\mu)$, which is the space of all measurable functions $\phi:\mathcal{M}\mapsto\mathbb{R}$ such that,

[TABLE]

As a natural extension, for any finite $n$, given the vector-valued function $\mathbf{\Phi}=\begin{bmatrix}\phi_{1}&\ldots&\phi_{n}\end{bmatrix}\in\mathcal{F}^{n}$, we define the corresponding norm as,

[TABLE]

2.2.1 Leveraging known physics: differential form

For any observation vector $\mathbf{\Phi}\in\mathcal{F}^{D}$, we can define the following Koopman error functional $\mathcal{J}[\cdot]:\mathcal{F}^{D}\mapsto\mathbb{R}_{\geq 0}$, $\mathbb{R}_{\geq 0}=[0,+\infty)$ as the square of maximal distance for all components in $\Phi$ between the action of $\mathcal{K}$ and its $L^{2}$ projection onto $\mathcal{F}_{D}$ [39],

[TABLE]

which describes the extent to which $\mathcal{F}_{D}$ is invariant to the Koopman operator $\mathcal{K}_{t}$ with $t\rightarrow 0$ with respect to each basis. Ideally, if we can find $\mathbf{\Phi}$ such that the corresponding $\mathcal{F}$ is invariant to $\mathcal{K}$, i.e., $\mathcal{J}[\mathbf{\Phi}]=0$, then $\mathcal{F}_{D}$ is invariant to $\mathcal{K}_{t}$, $\forall t>0$, i.e., a perfect linear embedding in the $L^{2}(\mu)$ sense. Once such an embedding is determined, the realization of $\mathcal{K}$ is simply the matrix $\mathbf{K}$. In this work, we are interested in $\mathbf{\Phi}$ such that one can recover $\mathbf{x}$ and $\mathcal{J}[\mathbf{\Phi}]\geq 0$.

Although the above problem setup only contains minimal necessary assumptions, it is both mathematically and computationally challenging. For practical purposes, following previous studies [60, 80, 79, 70, 48, 73], we consider the following assumptions to make the problem tractable:

Instead of solving the equation $\mathcal{J}[\mathbf{\Phi}]=0$, we search for $\mathbf{\Phi}$ by finding the minimum of the following constrained variational problem,

[TABLE] 2. 2.

Instead of directly solving the variational problem in the infinite dimensional $\mathcal{F}^{D}$, we optimize $\mathbf{\Phi}$ in the finite dimensional parameter space of feedforward neural networks with fixed architecture in which the number of layers and the number of hidden units in each layer is determined heuristically. Note that we are searching in a subset of $\mathcal{F}^{D}$ described by $\mathbf{W}_{\mathbf{\Phi}}$, which might induce a gap due to the choice of the neural network architecture,

[TABLE]

Clearly, the gap is bounded by the right hand side of the second inequality above. In addition, it should be noted that the requirement of linear independence in $\{\phi_{1},\ldots,\phi_{D}\}$ is relaxed, but $\dim\mathcal{F}_{D}$ is bounded by $D$. 3. 3.

As a standard procedure in deep learning [20], we use the penalty method to approximate the constrained optimization problem with an unconstrained optimization with unity penalty coefficient. Since this still entails nonconvex optimization, we define a global minima as follows:

[TABLE]

where the reconstruction error functional $\mathcal{R}[\cdot,\cdot]:\mathcal{F}^{D}\times C(\mathbb{R}^{D},\mathbb{R}^{N})\mapsto\mathbb{R}_{\geq 0}$, is defined as

[TABLE]

for $\mathbf{\Phi}\in\mathcal{F}^{D}$, $\mathbf{\Psi}\in C(\mathbb{R}^{D},\mathbb{R}^{N})$. Then we assume $\mathbf{\Phi}(\cdot;\widehat{\mathbf{W}}_{\mathbf{\Phi}})$ approximates one of the global minima, i.e., $\mathbf{\Phi}(\cdot;\mathbf{W}^{*}_{\mathbf{\Phi}})$. Note that the convergence to a global minimum of the constrained optimization can be guaranteed if one is given the global minima of the sequential unconstrained optimization and by increasing the penalty coefficient to infinity [47]. 4. 4.

We then optimize the sum of square error $\widehat{\mathcal{J}}$ over all components of $\mathbf{\Phi}$, which serves a upper bound of $\mathcal{J}$,

[TABLE]

The above formulation also implies equal importance among all components of $\mathbf{\Phi}$. 5. 5.

Despite the non-convex nature of the problem, we employ gradient-descent optimization to search for a local minimum.

In summary based on the above assumptions, we will solve the following optimization problem by gradient-descent:

[TABLE]

2.2.2 Unknown physics with only the trajectory data: recurrent form

From the viewpoint of approximation, it can be argued that the most natural way to determine the continuous Koopman operator is in the aforementioned differential form. However, higher accuracy can be achieved by taking advantage of trajectory information and minimizing the error over multiple time steps in the trajectory. This is the key idea behind optimized DMD [5]. Recently, Lusch et al. [48], Otto and Rowley [60] extended this idea to determine the discrete-time Koopman operator using deep learning.

Recall that we assume that the space of observation functions $\mathcal{F}=L^{2}(\mathcal{M},\mu)$ is invariant to $\mathcal{K}_{t}$, $\forall t\in\mathbb{R}_{\geq 0}$. Thus, we consider the $t$-time evolution of any function $\phi\in\mathcal{F}$, i.e., $\mathcal{K}_{t}\phi(\mathbf{x})$, as a function on $\mathcal{U}=\mathcal{M}\times\mathcal{T}$, given an initial condition $\mathbf{x}\in\mathcal{M}$ and time evolution $t\in\mathcal{T}$ where $\mathcal{T}$ is a one dimensional smooth manifold, sometimes also referred to as the time horizon. We assume $\mathcal{K}_{t}\phi(\mathbf{x})\in\mathcal{G}=L^{2}(\mathcal{U},\nu)$. Based on the fact that $\mathcal{F}$ is invariant to $\mathcal{K}_{t}$, such an assumption can be shown to be valid for compact $\mathcal{T}$ for a proper measure $\nu$.

Similar to the differential form, for any observation vector $\mathbf{\Phi}\in\mathcal{F}^{D}$, we define the following Koopman error functional $\mathcal{J}_{r}[\cdot]:\mathcal{F}^{D}\mapsto\mathbb{R}_{\geq 0}$ as the square of maximal $L^{2}$ distance for all components in $\mathbf{\Phi}$ between the predicted and ground truth trajectory,

[TABLE]

where $\psi=\mathbf{\Phi}c_{\psi}$, $c_{\psi}\in\mathbb{R}^{D\times 1}$. Following similar assumptions in Section 2.2.1, we need to define a reconstruction error functional to describe the discrepancy between the reconstructed and original states. Indeed, one can directly define the prediction error functional in the recurrent form, $\widetilde{\mathcal{P}}[\cdot,\cdot,\cdot]:\mathcal{F}^{D}\times C(\mathbb{R}^{D},\mathbb{R}^{N})\times\mathbb{R}^{D\times D}\mapsto\mathbb{R}_{\geq 0}$ as

[TABLE]

Similarly, we define

[TABLE]

and solve the following optimization problem via a gradient-based algorithm:

[TABLE]

As a side note, one might also define the following reconstruction functional similar to previous differential form that is independent of $\mathbf{K}$,

[TABLE]

which can bound the prediction error functional together with the Koopman error functional by triangular inequality,

[TABLE]

where $L_{\mathbf{\Psi}}$ is the Lipschitz constant for $\mathbf{\Psi}$. However, in this work, we do not have control over $L_{\mathbf{\Psi}}$, and thus we prefer to directly minimize the prediction error function as in previous studies [60, 48].

2.3 Measures

2.3.1 Measure for data generated by sampling in phase space

We consider the situation where $M$ data points on $\mathcal{M}$, i.e., $\{\mathbf{x}_{m}\}_{m=1}^{M}$ are drawn independently from some measure $\mu$ on $\mathcal{M}$, e.g., uniform distribution. This induces the following empirical measure $\hat{\mu}_{M}$,

[TABLE]

where $\delta_{x}$ is the Dirac measure for $x$. Note that $\hat{\mu}_{M}$ uniformly converges to $\mu$ [76] as $M\rightarrow\infty$. Thus, one can rewrite the differential form in Eq. 18 as an empirical risk minimization [75],

[TABLE]

where $\widehat{\mathcal{F}}_{M}=L^{2}(\mathcal{M},\hat{\mu}_{M})$.

2.3.2 Measure for trajectory data generated by solving the initial value problem

In the general case, information content in trajectory data resulting from the solution of the initial value problem is strongly dependent on the initial condition. For instance, when the initial condition is in a region of phase space with sharp changes, it would be sensible to use a high sampling rate. On the other hand, if the initial condition is near a fixed point attractor, one might prefer to stop collecting the data after the system arrives at the fixed point. Such a sampling pattern for the specific initial state can be summarized as a Markov kernel [34] $\kappa:\mathcal{M}\times\Sigma_{\mathcal{T}}\mapsto[0,1]$. For every fixed initial state $\mathbf{x}\in\mathcal{M}$, the map $\kappa(\mathcal{E}_{\mathcal{T}},\mathbf{x})$ is a measure on $(\mathcal{T},\Sigma_{\mathcal{T}})$, $\mathcal{E}_{\mathcal{T}}\in\Sigma_{\mathcal{T}}$, where $\Sigma_{\mathcal{T}}$ is the $\sigma$-algebra on $\mathcal{T}$. Then, there exists a unique measure $\nu$ on $\mathcal{U}$ [34] such that,

[TABLE]

for $\mathcal{E}_{\mathcal{T}}\in\Sigma_{\mathcal{T}},\mathcal{E}_{\mathcal{M}}\in\Sigma_{\mathcal{M}}$.

Assume we are given $M$ trajectories, $\{\{\mathbf{x}_{m,j}\}_{j=1}^{T_{m}}\}_{m=1}^{M}$, $\{\{t_{m,j}\}_{j=1}^{T_{m}}\}_{m=1}^{M}$. The initial condition for each trajectory is drawn independently from $\mu$. For the $m$-th trajectory with initial condition $\mathbf{x}_{m}$, there are $T_{m}$ samples drawn independently from measure $\kappa(\cdot,\mathbf{x}_{m,1})$ where the time elapse of $j$-th sample away from initial condition $\mathbf{x}_{m,1}$ is $t_{m,j}$ ($t_{m,1}=0$), where $m=1,\ldots,M$. Then, we can define the corresponding empirical measure from Eq. 28,

[TABLE]

where $\displaystyle\widehat{T}=\max_{m=1,\ldots,M}T_{m}$, and $\hat{\nu}_{M,\widehat{T}}$ uniformly converges to $\nu$ as $M,\widehat{T}\rightarrow\infty$.

Similarly, we can rewrite Eq. 22 as the following,

[TABLE]

where $\widehat{G}=L^{2}(\mathcal{U},\widehat{\nu}_{M,\widehat{T}})$.Note that Eq. 30 generalizes the LRAN model with a simpler loss function [60] and the discrete spectrum model in the paper of Lusch et al. [48].

Note that the above data setup also generalizes to cases where data along a single long trajectory is “Hankelized” [60, 2], i.e., dividing a long trajectory into several smaller-sized consecutive trajectories. This would lead to a truncated time horizon in the model, which leads to better computational efficiency compared to the original data, at the cost of some loss in prediction accuracy.

2.4 Guaranteed stabilization of the Koopman operator

Eigenvalues of the Koopman operator are critically important in understanding the temporal behavior of certain modes in the dynamical system. For a measure-preserving system [11] or systems on an attractor, e.g., post-transient flow dynamics [3], even if the system is chaotic, the eigenspectrum of the continuous-time Koopman operator would still be on the imaginary axis. It should be noted that although the Koopman operator still accepts unstable modes, i.e., the real part of the eigenvalues of continuous-time Koopman operator being positive, its absence in systems governed by Navier-Stokes equation in fluid mechanics has been documented [55]. Hence, in this work, we assume that the Koopman eigenvalues corresponding to the finite dimensional Koopman invariant subspace of interest have non-positive real parts. It is important to note that the concept of unstable Koopman modes should not be confused with that of flow instability. Prior models (for instance,  [60]) have not explicitly taken stability into account, and thus resulted in slightly unstable Koopman modes. While this is acceptable for relatively short time predictions, long time predictions will be problematic.

Before presenting the stabilization technique, it is instructive to note the non-uniqueness of ideal observation functions $\mathbf{\Phi}$, i.e., the one corresponding to the exact Koopman invariant subspace. This is because one can simultaneously multiply any $D\times D$ invertible real matrix and its inverse before and after the observation vector $\mathbf{\Phi}$ while keeping the Koopman eigenfunction, eigenvalues, and the output from the reconstruction the same. Thus, observation functions described by neural networks cannot be expected to be uniquely determined. We will leverage this non-uniqueness to enforce stability.

Enlightened by recent studies [24, 13, 12] in the design of stable deep neural network structures where skew-symmetric weights are used inside nonlinear activations, we devised a novel parameterization for the realization of the Koopman operator in the following form:

[TABLE]

where $\zeta_{1},\ldots,\zeta_{D-1}$, $\sigma_{1},\ldots,\sigma_{D}\in\mathbb{R}$.

In Appendix A, we show that the real parts of the eigenvalues of a $n\times n$ real (possibly non-symmetric) negative semi-definite matrix $\mathbf{A}$ are non-positive. We then prove that the above parameterization posits a constraint that the time evolution associated with $\mathbf{K}_{stable}$ in Eq. 31 for any choice of parameters in $\mathbf{R}$ is always stable. Further, the constraint from the parameterization in Eq. 31 is actually rich enough such that any diagonalizable matrix corresponding to a stable Koopman operator can be represented without loss of expressivity.

Theorem 2.1**.**

*For any real square diagonalizable matrix $\mathbf{K}\in\mathbb{R}^{D\times D}$ that only has non-positive real parts of the eigenvalues $D\geq 2$, there exists a set of $\zeta_{1},\ldots,\zeta_{D-1}$, $\sigma_{1},\ldots,\sigma_{D}\in\mathbb{R}$ such that the corresponding $\mathbf{K}_{stable}$ in Eq. 31 is similar to $\mathbf{K}$ over $\mathbb{R}$. Moreover, for any $\zeta_{1},\ldots,\zeta_{D-1}$, $\sigma_{1},\ldots,\sigma_{D}\in\mathbb{R}$, the real part of the eigenvalues of the corresponding $\mathbf{K}_{stable}$ is non-positive. *

Proof 2.2**.**

For any real square diagonalizable matrix $\mathbf{K}\in\mathbb{R}^{D\times D}$ that only has non-positive real parts of the eigenvalues, there exists an eigendecomposition,

[TABLE]

Without loss of generality, the diagonal matrix $\mathbf{J}$ contains $2D_{c}$ complex eigenvalues $\{\lambda^{c}_{j}\}_{j=1}^{2D_{c}}$ and $D_{r}$ real eigenvalues $\{\lambda^{r}_{j}\}_{j=1}^{D_{r}}$ where $2D_{c}+D_{r}=D$.

Consider a $2\times 2$ real matrix, $\mathbf{A}_{j}=\begin{bmatrix}-\sigma_{j}^{2}&\zeta_{j}\\ -\zeta_{j}&-\sigma_{j}^{2}\end{bmatrix}$ where the eigenvalues are $\lambda_{1,2}=-\sigma_{j}^{2}\pm j\zeta_{j}^{2}$. We use this matrix to construct a $2\times 2$ matrix that has eigenvalues $\lambda^{c}_{j}$. For each pair of complex eigenvalues, $j=1,\ldots,D_{c}$, we have $\mathbf{A}_{j}=\begin{bmatrix}\textrm{Re}(\lambda^{c}_{j})&\sqrt{|\textrm{Im}(\lambda^{c}_{j})|}\\ -\sqrt{|\textrm{Im}(\lambda^{c}_{j})|}&\textrm{Re}(\lambda^{c}_{j})\end{bmatrix}$. Next, combining with the $D_{r}$ real eigenvalues, we have the following block diagonal matrix that shares the same eigenvalues as $\mathbf{J}$,

[TABLE]

Since two matrices with the same (complex) eigenvalues are similar [53], $\mathbf{K}$ is similar to $\widetilde{\mathbf{K}}$ over $\mathbb{C}$. However, since $\mathbf{K}$ and $\widetilde{\mathbf{K}}$ are matrices over $\mathbb{R}$, they are already similar over $\mathbb{R}$ (Refer Corollary 1 on p. 312 in Ref. [26]). Notice that the general form of $\widetilde{\mathbf{K}}$ is just a special form of $\mathbf{K}_{stable}$ in Eq. 31. Therefore, one can always find a set of $\zeta_{1},\ldots,\zeta_{D}$ and $\sigma_{1},\ldots,\sigma_{D}$ such that $\mathbf{K}_{stable}$ is similar to $\mathbf{K}$ over $\mathbb{R}$. Next, notice that given $\sigma_{1},\ldots,\sigma_{D}$, $\zeta_{1},\ldots,\zeta_{D}$, $\forall\mathbf{v}\in\mathbb{R}^{D\times 1}$, we have

[TABLE]

*Therefore, $\mathbf{K}_{stable}$ is a negative semi-definite matrix. Following Lemma A.2, the real part of any eigenvalue of $\mathbf{K}_{stable}$ is non-positive. *

Thanks to the conjugation over $\mathbb{R}$, the last layer of encoder, i.e., $\mathbf{\Phi},$ and the first layer of the decoder, i.e., $\mathbf{\Psi}$, will absorb any linear transformation necessary. Hence, the above theorem simply means one can parameterize $\mathbf{K}$ with Eq. 31 that guarantees stability of the linear dynamics of $\mathbf{\Phi}$ without loss of expressibility for cases where unstable modes are absent.

Careful readers might notice that indeed one can further truncate the parameterization with a $2\times 2$ block matrix instead of the tridiagonal form used in Eq. 31. However, since this reduction takes more effort in the implementation while the parameterization in Eq. 31 has already reduced the number of parameters from $O(D^{2})$ to $O(D)$, we prefer the tridiagonal form. For the rest of the work, we will use this parameterization for all of the cases concerned.

2.5 Design of neural network architecture with SVD-DMD embedding

In the previous subsection, we described the deep learning formulation for the Koopman operator as a finite dimensional optimization problem that approximates the constrained variational problem, together with a stable parameterization for $\mathbf{K}$. In this subsection, we will describe the design of our neural network that further embeds the SVD-DMD for differential and recurrent forms, which were presented in Sections 2.2.1 and 2.2.2.

2.5.1 Input normalization

As a standard procedure, we consider normalization on the snapshot matrix of state variable $\mathbf{X}$,

[TABLE]

for better training performance on neural networks [20]. Specifically, we consider Z-normalization shown in Eq. 36, i.e., subtracting the mean of $\mathbf{x}$ then dividing the standard deviation to obtain $\mathbf{z}$.

[TABLE]

where $\overline{\mathbf{x}}=\frac{1}{M}\sum_{m=1}^{M}\mathbf{x}_{m}$, $\mathbf{\Lambda}=\textrm{diag}\{d_{1},\ldots,d_{N}\}$, $d_{j}$ is the uncorrected standard deviation of $j$-th component of $\mathbf{x}$, i.e.,

[TABLE]

where $\mathbf{x}_{m,j}$ is $j$-th component of $m$-th data, $j=1,\ldots,N$.

While such a normalization is helpful for neural network training in most cases, in some cases where the data is a set of POD coefficients, it cannot differentiate between components that could be more significant than others. Therefore, for those cases, we consider a different normalization with the $\mathbf{\Lambda}$ that sets the ratio of standard deviation between components as:

[TABLE]

where $d_{max}=\max_{j=1,\ldots,N}d_{j}$, and $\mathbf{I}$ is the identity matrix.

2.5.2 Embedding with SVD-DMD

Instead of directly using the standard feedforward neural network structure employed in previous works, we embed SVD-DMD [70] into the framework and learn the residual. Recall for the standard SVD-DMD algorithm, given $M$ sequential snapshots in Eq. 35 uniformly sampled at intervals $\Delta t$, one linearly approximates the action of $\mathcal{K}_{\Delta t}$ of the first $r$ dominant SVD modes of centered snapshots, i.e., $\mathbf{q}=(\mathbf{x}-\overline{\mathbf{x}})\mathbf{V}_{r}$. Here, $r$ is empirically chosen as a balance between numerical robustness and reconstruction accuracy, and $\mathbf{V}_{r}$ are the first $r$ columns from $\mathbf{V}$ of the SVD of centered snapshots,

[TABLE]

Then the SVD-DMD operator is simply the matrix $\mathbf{A}$ that minimizes $\sum_{m=1}^{M-1}\lVert\mathbf{q}_{m+1}-\mathbf{q}_{m}\mathbf{A}\rVert$, where $\mathbf{q}_{m}$ is the corresponding orthogonal projection of $\mathbf{x}_{m}$.

To embed the above SVD-DMD structure into the neural network, we introduce two modifications. First, we take $r=D$. Since $D$ is arbitrary, if $N<D$, we simply append zero columns in $\mathbf{V}_{D}$. Second, we would also need to accommodate $\mathbf{z}$ with the input normalization in Section 2.5.1. Thus, we cast $\mathbf{q}=\mathbf{z}\mathbf{\Lambda}\mathbf{V}_{D}$. We then have,

[TABLE]

where \mathbf{\Phi}_{nn}\triangleq\mathbf{\Phi}(\cdot;{\color[rgb]{0,0,0}\mathbf{W}_{\mathbf{\Phi}}\setminus{\{W_{enc,L}\}}}), $\mathbf{W}_{\mathbf{\Phi}}=\{W_{enc,1},b_{enc,1},\ldots,W_{enc,L}\}$, $\mathbf{\Psi}_{nn}\triangleq\mathbf{\Psi}(\cdot;\mathbf{W}_{\mathbf{\Psi}})$, $\mathbf{W}_{\mathbf{\Psi}}=\{W_{dec,1},\ldots,W_{dec,L},b_{dec,L}\}$, $L$ is the number of layers for encoder or decoder neural network. The embedding is illustrated in Fig. 1.

The intuition behind the embedding of the SVD-DMD into the framework is given below:

Schmid’s DMD algorithm [70] relies on the dominant POD modes to reduce the effect of noise, and has been shown to approximate the Koopman invariant subspace [69]. This has been demonstrated even for high-dimensional nonlinear dynamical systems with millions of degrees of freedom [71]. Moreover, for systems with continuous spectra, POD appears to be a robust alternative to Koopman mode decomposition [56]. Thus, we assume that the true Koopman invariant subspace is easier to obtain by only learning the residual with respect to the dominant POD subspace. 2. 2.

Although we did not precisely implement ResNet blocks [25] (empirically, for the problem concerned in this paper, our network does not have to be as deep as common architectures in the deep learning community [25]), we believe that such fixed mappings may have similar benefits in ResNets. 3. 3.

For an ideal linear dynamical system, using the aforementioned neural network model with nonlinear reconstruction can result in an infinite number of unnecessary Koopman modes as global minima. For example, consider the case of a linear dynamical system, i.e., $\mathbf{F}(\mathbf{x})=\mathbf{x}\mathbf{A}$, $\mathbf{A}\in\mathbb{R}^{N\times N}$. The desired Koopman invariant subspace is trivially the span of the projections of $\mathbf{x}$ onto each component and the desired Koopman eigenvalues are simply that of $\mathbf{A}$. Assuming simple and real eigenvalues, one can have $\dot{\mathbf{y}}=\mathbf{y}\mathbf{\Lambda}$, where $\mathbf{y}=\mathbf{x}\mathbf{M}$, $\mathbf{A}=\mathbf{M}\mathbf{J}\mathbf{M}^{-1}$ as the eigen-decomposition with $\mathbf{J}=\textrm{diag}\{\lambda_{j}\}_{j=1}^{N}$. Then for each component $y_{j}$, we have $\dot{y_{j}}=\lambda_{j}y_{j}$. For any $n_{j}\in\mathbb{N}$, consider the observable $\phi_{j}=y_{j}^{2n_{j}+1}$, one can have $\dot{\phi_{j}}=(2n_{j}+1)\lambda_{j}\phi_{j}$, i.e., $\textrm{span}\{\phi_{1},\ldots,\phi_{N}\}$ is invariant to $\mathcal{K}_{t}$ 222Budišić et al. [11] showed that the set of eigenfunctions naturally forms an Abelian semigroup under pointwise products.. Then, consider the nonlinear decoder as one that simply takes the $2n_{j}+1$-th root on $\phi_{j}$, and one can recover $\mathbf{y}$ exactly. Finally, augmenting the decoder with $\mathbf{M}^{-1}$ and the encoder with $\mathbf{M}$, the neural network model can find a spurious linear embedding, with eigenvalues as $\{(2n_{1}+1)\lambda_{1},\ldots,(2n_{N}+1)\lambda_{N}\}$ rather than the desired $\{\lambda_{1},\ldots,\lambda_{N}\}$, which is an over-complicated nonlinear reconstruction. It is trivial to generalize the above thought experiment to cases where eigenvalues are complex for a real-input-real-output neural network to accommodate. On the other hand, since most often the neural network is initialized with small weights near zero, the effect of the nonlinear encoder and decoder can be small initially compared to the DMD part. Thus, if the system can be exactly represented by DMD, the optimization for the embedded architecture is initialized near the desired minimum. We note that this could lead to attenuation of the spurious modes due to the nonlinear reconstruction for essentially linear dynamics.333It has to be mentioned that that such an issue could exist also in cases where DMD is not desired.

The neural network architecture for the differential form is shown in Fig. 2,

and for the recurrent form in Fig. 3.

Note that, if the nonlinear part, i.e., the feedforward neural network is not activated, the above formulation reverts to an over-parameterized SVD-DMD. Specifically, the recurrent form model with neural network deactivated can be viewed as a simplified version of optimized DMD [5, 14].

2.6 Implementation

The framework is built using Tensorflow [1]. Neural network parameters $\mathbf{W}_{\mathbf{\Phi}},\mathbf{W}_{\mathbf{\Psi}}$, are initialized with the standard truncated normal distribution. $\mathbf{K}$ is initialized with the corresponding DMD approximation. The objective function is optimized using Adam [33], which is an adaptive first order stochastic optimization method using gradient updates scaled by square roots of exponential moving averages of previous squared gradients. Note that we also include weight decay regularization, $\mathcal{L}_{reg}=\left\lVert\bm{W}_{\mathbf{\Phi}}\right\rVert^{2}+\left\lVert\bm{W}_{\mathbf{\Psi}}\right\rVert^{2}$, in the objective function, to avoid spurious oscillations in the learned Koopman functions, which helps generalization in an interpolation sense [20].

3 Probabilistic formulation

3.1 Bayesian neural networks

A Bayesian formalism is adopted to quantify the impact of several sources of uncertainty in model construction on the model predictions. Bayes’ rule is

[TABLE]

where $\mathcal{D}$ is data, and $\mathbf{\Theta}$ is the set of parameters. For simplicity, $P$ represents the probability density function (PDF) on the measure space generated by the data and parameters. From a traditional Bayesian standpoint, as the number of parameters in the neural network is large, it is impossible for common inference tools such as Markov Chain Monte Carlo (MCMC) to be practical. To overcome the curse of dimensionality, several general approaches such as Laplacian approximation [16] and variational inference [40, 7] have been proposed. The former is computationally economical but has two major limitations: First, computing the full Hessian is impossible and expensive for a high dimensional problem and most often approximations such as $\mathcal{J}^{\top}\mathcal{J}$, where $\mathcal{J}$ is the Jacobian are employed [49]. Second, it only provides a local approximation, which can often be far removed from the true posterior. Variational inference has become popular in the deep learning community as it offers an informed balance between the computationally expensive MCMC method, and the cheap but less descriptive models such as the Laplacian approximation. Historically, variational inference for neural networks [27] has been difficult [59] to formulate, largely due to the difficulty of deriving analytical solutions to the required integrals over the variational posteriors [22] even for simple network structures. Graves [22] proposed a stochastic method for variational inference with a diagonal Gaussian posterior that can be applied to almost any differentiable log-loss parametric model, including neural networks. However, there is always a trade-off between the complexity of the posterior and scalability and robustness [82]. In this work, we adopt the mean-field variational inference [40].

3.2 Mean-field variational inference

As illustrated in the left figure of Fig. 4, the key idea in variational inference [6] is to recast Bayesian inference as an optimization problem by searching the best parameterized probability density function $q(\mathbf{\Theta};\hat{\xi})$ in a family of approximating densities, namely the variational posterior, $\{q(\mathbf{\Theta};\xi)|\xi\in\Xi\}$, such that it is closest to the true posterior $P(\mathbf{\Theta}|\mathcal{D})$. Most often, the Kullback–Leibler (KL) divergence is employed to measure the distance, which is defined as \textrm{KL}(q(\mathbf{\Theta};{\xi})\|P(\mathbf{\Theta}|\mathcal{D}))=\int_{\Omega}q(\mathbf{\Theta};\xi)\log\frac{q(\mathbf{\Theta};\xi)}{P(\mathbf{\Theta}|\mathcal{D})}d\mathbf{\Theta}=\mathbb{E}_{q(\mathbf{\Theta};\xi)}\big{[}\log{q(\mathbf{\Theta};\xi)}{P(\mathbf{\Theta}|\mathcal{D})}\big{]}, where $\Omega$ is the support of $q(\mathbf{\Theta};\xi)$. This implies $\textrm{supp}(q(\mathbf{\Theta};\xi))\subseteq\textrm{supp}(P(\mathbf{\Theta}|\mathcal{D}))$. Further, we assume $\textrm{supp}(P(\mathbf{\Theta}|\mathcal{D}))=\textrm{supp}(P(\mathbf{\Theta}))$. $\Xi$ is the domain of $\xi$, depending on the parameterization and family of approximating densities.

However, direct computation of the KL divergence is intractable, since we do not have access to $\log P(\mathcal{D})$. Instead, we choose an alternative, the evidence lower bound (ELBO), i.e., the negative KL divergence plus $\log P(\mathcal{D})$, in Eq. 44 to be maximized.

[TABLE]

To maximize the ELBO, we leverage automatic differentiation from Tensorflow to compute the gradients with respect to $\xi$, following the framework of Automatic Differentiation Variational Inference (ADVI) [40] where Gaussian distributions are considered as the variational family. Specifically, we employ the mean-field assumption in Eq. 45 such that,

[TABLE]

where $Z$ is the total number of parameters. For $j=1,\ldots,Z$, $\theta_{j}$ is the $j$-th parameter as random variable, and $\xi_{j}$ is the corresponding variational parameter that describes the distribution. This is particularly convenient for neural network models constructed in Tensorflow since weights and biases are naturally defined on some real coordinate space. If the support of the parameter distribution is restricted, one can simply consider a one-to-one differentiable coordinate transformation $\Upsilon(\mathbf{\Theta})=\mathbf{Z}$, such that $\mathbf{Z}$ is not restricted in some real coordinate space, and posit a Gaussian distribution on $\mathbf{Z}$. Note that this naturally induces non-Gaussian distribution. Here, we employ the ADVI functionality in Edward [74], which is built upon Tensorflow to implement ADVI. Interested readers should refer to the original paper of ADVI [40] for the specific details of implementing mean-field variational inference including the usage of the reparameterization-trick to compute the gradients. In contrast, note that the maximum a posteriori (MAP) estimation of the posterior can be cast as a regularized deterministic model as illustrated in Fig. 4. Since weight decay is employed in previous deep learning models [60, 48] to learn the Koopman decomposition, one can show that the previous model is essentially a MAP estimation of the corresponding posterior.

Note that the mean-field Gaussian assumption is still simplified, yet effective and scalable for deep neural nets [7]. It is interesting to note that several recent works [83, 84] leverage Stein Variational Gradient Descent (SVGD) [44], a non-parametric, particle-based inference method, which is able to capture multi-modal posteriors. Robustness and scalability to high dimensions, e.g., deep neural nets, is still an area of active research [85, 45, 72, 46].

3.3 Bayesian hierarchical model setup

Recall that we have the following parameters for the deep learning models introduced in Section 2.2.1 and Section 2.2.2: 1) weights and biases for the “encoder”, $\mathbf{W}_{\mathbf{\Phi}}$, 2) weights and biases for the “decoder”, $\mathbf{W}_{\mathbf{\Psi}}$, and 3) stabilized $\mathbf{K}$ with $\zeta_{1},\ldots,\zeta_{D}$ and $\sigma_{1},\ldots,\sigma_{D}$. Based on mean-field assumptions, we just need to prescribe the prior and variational posterior for each parameter. For each weight and bias, we posit a Gaussian prior with zero mean, with the scaling associated with each parameter to follow the recommended half–Cauchy distribution [19, 65], which has zero mean and scale as 1 (empirical). The variational posterior for each weight and bias is Gaussian, and Log-normal for scale parameters. For the off-diagonal part of $\mathbf{K}$, $\zeta_{i}$, we posit the same type of Gaussian prior as before with the scale parameter following a hierarchical model for each $i=1,\ldots,D$. The corresponding variational posterior for $\zeta_{i}$ is still Gaussian while log-normal for the scale parameter. For the non-negative diagonal part of $\mathbf{K}$, we posit a Gamma distribution for each $\sigma_{i}^{2}$, $i=1,\ldots,D$, with rate parameter as 0.5 and shape parameter following the previous half-Cauchy distribution. Variational posteriors for both $\sigma^{2}_{i}$ and its shape parameter are log-normal.

For the differential form in Section 2.2.1, we set up the following likelihood function444Note that independence between the data and the structure of the aleatoric uncertainty is assumed. Such an assumption correlates well with existing deterministic models based on mean-square-error. based on Z-normalized data $\mathcal{D}=\{\mathbf{z}_{m},\mathbf{\dot{z}}_{m}\}_{m=1}^{M}$:

[TABLE]

where $\mathbf{\Lambda}_{rec},\mathbf{\Lambda}_{lin}$ are diagonal covariance matrices with the prior of each diagonal element following the previous half–Cauchy distribution. We also posit log-normals to infer the posterior of $\mathbf{\Lambda}_{rec},\mathbf{\Lambda}_{lin}$. We denote $\mathbf{\Theta}$ as $\mathbf{W}_{\mathbf{\Phi}},\mathbf{W}_{\mathbf{\Psi}}$, and associated scale parameters together with $\mathbf{\Lambda}_{rec}$, $\mathbf{\Lambda}_{lin}$.

For the recurrent form in Section 2.2.2, given normalized data,

[TABLE]

we consider the following likelihood, 555Alternative likelihoods can be chosen to account for the variation of aleatoric noise in time, which is well-suited for short-horizon forecasting. However, we are more interested in a free-run situation where only the initial condition is given.:

[TABLE]

3.4 Propagation of uncertainties

Given data $\mathcal{D}$ and the inferred posterior $P(\mathbf{\Theta}|\mathcal{D})$, assuming a noise-free initial condition $\mathbf{z}_{0}$, we are interested in future state predictions with uncertainties. For the differential form in Eq. 46, we have

[TABLE]

However, $P(\mathbf{\Phi}(t)|\mathbf{K},\mathbf{\Phi}(\mathbf{z}_{0};\mathbf{W}_{\mathbf{\Phi}}))$ is unknown because the differential form does not use trajectory information. If we assume multivariate Gaussian white noise with the same covariance in the linear loss $\mathbf{\Lambda}_{lin}$, then one can forward propagate the aleatoric uncertainty associated with the likelihood function of the linear loss. Then one immediately recognizes that the continuous-time random process of $\mathbf{\Phi}(t)$ becomes a multivariate Ornstein–Uhlenbeck process,

[TABLE]

where $\mathbf{B}(t)$ is a $D$-dimensional Gaussian white noise vector with unit variance. Note that [68]

[TABLE]

where $\int_{0}^{t}e^{s\mathbf{K}}\mathbf{\Lambda}_{lin}e^{s\mathbf{K}^{\top}}ds=\mathbf{vec}^{-1}(-(\mathbf{K}\oplus\mathbf{K})^{-1}(\mathbf{I}-e^{t(\mathbf{K}\oplus\mathbf{K})})\mathbf{vec}(\mathbf{\Lambda}_{lin}))$. $\mathbf{vec}(\cdot)$ is the stack operator, and $\oplus$ is the Kronecker sum [52].

It is interesting to note that, since $\mathbf{K}$ is restricted by Eq. 31 and does not contain any eigenvalues with positive real part, the variance in Eq. 51 will not diverge in finite time. One can thus simply draw samples of $\mathbf{\Phi}(t)$ from Eq. 51. Thus, we approximate Eq. 49 with Monte Carlo sampling from the corresponding variational posterior,

[TABLE]

where the superscript with parentheses represents the index of samples, $N_{mc}$, $M_{mc}$ are the number of samples corresponding to variational posteriors and the Ornstein–Uhlenbeck process.

For the recurrent form in Eq. 48, the posterior predictive distribution of $\mathbf{z}(t)$ given the initial condition is straightforward:

[TABLE]

4 Numerical Results & discussion

To demonstrate and analyze the approaches presented herein, three numerical examples are pursued.

4.1 2D fixed point attractor

Consider the two-dimensional nonlinear dynamical system [48] with a fixed point,

[TABLE]

where $\mu=-0.05$, $\lambda=-1$. For this low dimensional system with known governing equations, we consider the differential form in Section 2.2.1, with 1600 states as training samples, sampled from $[-0.5,0.5]$ using the standard Latin-Hypercube-Sampling method [51] for both $x_{1}$ and $x_{2}$. The embedding of SVD-DMD is not employed in this case since it is not very meaningful. The hyperparameters for training are given in  Table 2.

After we obtained the inferred posterior, we consider Monte Carlo sampling described in Section 3.4 with $N_{mc}=100$, $M_{mc}=10$ to approximate the posterior distribution of the Koopman observables and prediction on the dynamical system given an initial condition $\mathbf{x}_{0}$. The mean Koopman eigenvalues are $\lambda_{1}=-0.99656$ and $\lambda_{2}=-0.05049$, and the amplitude and phase angle for the mean eigenfunctions are shown in Fig. 5, which resembles the analytic solution with $\varphi_{1}=x_{2}-\lambda x_{1}^{2}/(\lambda-2\mu)$, $\varphi_{2}=x_{1}$.

Finally, since we have a learned a distribution over the Koopman operator, we can obtain the posterior distribution of the predicted dynamics, given an arbitrary unseen initial condition following Section 3.4, for example, at $\mathbf{x}_{0}=(0.4,-0.4)$. The effect of the number of data samples $M$ on the confidence of the predicted dynamics can be seen in Fig. 6. Clearly, as more data is collected in the region of interest, the propagated uncertainty of the evolution of the dynamics predicted on the testing data decreases as expected. It is interesting to note that, even when the data is halved, the standard deviation of the Koopman eigenvalue is very small compared to the mean.

4.2 2D unforced duffing oscillator

Next, we consider the unforced duffing system:

[TABLE]

where $\delta=0.5$, $\beta=-1$, $\alpha=1$. We use 10000 samples of the state $\mathbf{x}$ distributed on $x_{1},x_{2}\in[-2,2]$ from Latin-Hypercube-Sampling. We infer the posterior using the differential form model in Section 2.2.1, with the following hyperparameters in Table 3.

To assess the uncertainty in the Koopman eigenfunctions, we draw 100 samples from the variational posterior of $\mathbf{W}_{\mathbf{\Phi}}$ and $\mathbf{K}$. First, the mean of the non-trivial Koopman eigenvalues are $\lambda_{1,2}=-0.535\pm 0.750i$ and the mean of the third eigenvalue is $\lambda_{3}=-2\times 10^{-5}$. The mean of the module and phase angle of the Koopman eigenfunctions is shown in Fig. 7. The results are similar to Ref. [60] in which a deterministic model is employed. The Koopman eigenfunction associated with $\lambda_{3}$ acts as an indicator of the basin of attraction. Second, to better visualize the effect of uncertainty on unseen data, we normalize the standard deviation of the module of Koopman eigfunctions by the minimal standard deviation over $[-4,4]$. Fig. 7 shows a uniformly distributed standard deviation within $[-2,2]\times[-2,2]$ where sampling data is distributed, and a drastic increase outside that training region as expected. Note that the area where the normalized standard deviation is larger than ten is cropped.

To obtain the posterior distribution of the evolution of the dynamics predicted by the Koopman operator, we arbitrarily choose an initial condition within the range of training data as $\mathbf{x}(0)=[1.2,1.2]$. Note that we did not input the model with any trajectory data other than the state and the corresponding time derivative obtained from the governing equation. With Monte Carlo sampling, we obtain the distribution of the trajectory in Fig. 7. Note that the uncertainty is quite small since it is within the training data region with enough data.

4.3 Flow wake behind a cylinder

We consider the velocity field of a two-dimensional laminar flow past a cylinder in a transient regime at a Reynolds number $Re_{D}=U_{\infty}D/\nu=100$ , which is above 47, the critical Reynolds number associated with Hopf bifurcation. $U_{\infty}$ is the freestream velocity, $D$ is the cylinder diameter, and $\nu$ is the kinematic viscosity. The transient regime is rather difficult due to the high nonlinearity of post-Hopf-bifurcation dynamics between unstable equilibrium and the limit cycle [14]. Data is generated by solving the two dimensional incompressible Navier–Stokes equations using OpenFOAM [30]. Grid convergence was verified using a sequence of successively refined meshes.

The initial condition is a uniform flow with the freestream velocity superimposed with standard Gaussian random noise, and with pressure initialized with Gaussian random noise. Although this initial condition is a rough approximation to the development of true instabilities from equilibrium [14], the flow is observed to rapidly converge to a quasi-steady solution after a few steps, and then starts to oscillate and form a long separation bubble with two counter-rotating vortices. The first 50 POD snapshots of velocity are considered with the kinetic energy captured upto 99%. We sample 1245 snapshots of data on the trajectory starting from the unstable equilibrium point to the vortex shedding attractor with $\Delta t=0.1t_{ref}=0.1D/U_{\infty}$ where the characteristic advection time scale $t_{ref}=D/U_{\infty}=2\textrm{ sec}$. The first 600 snapshots are considered as training data, i.e., $0\leq t\leq 60t_{ref}$.

To further analyze the robustness of the model to noisy data, Gaussian white noise is added 666Note that noise was not added to the original flow field since taking the dominant POD modes on the flow field would contribute to de-noising. to the temporal data of POD coefficients by considering a fixed signal-to-noise ratio as 5%, 10%, 20%, 30%, for each component at each time instance. The model performance is evaluated by predicting the entire trajectory with the (noisy) initial condition given, including the remaining 645 snapshots. We consider the recurrent form model described in Section 2.2.2 together with SVD-DMD described in Section 2.5.2 with the corresponding hyperparameters given in Table 4. Note that we consider 20 intrinsic modes, i.e., at most 10 distinct frequencies can be captured, which are empirically chosen. Also, the finite-horizon window length corresponding to $T=100$ is $10t_{ref}$, which is much less than the time required for the system to transit from unstable equilibrium to the attractor.

The continuous-time Koopman eigenvalue distribution of 20 modes for the training data with five different signal-to-noise ratios is shown in Fig. 8. First, all Koopman eigenvalues are stable according to the stable parameterization of $\mathbf{K}$ in Eq. 31. Second, when noise is added, the eigenvalues are seen to deviate except the one on the imaginary axis with $\lambda=\pm 0.528j$, which corresponds to the dominant vortex shedding frequency on the limiting cycle with $St=\lambda D/(2\pi U_{\infty})=0.168$. This exercise verifies the robustness of the present approach to Gaussian noise.

Recall that we have a posterior distribution of the predicted trajectory of POD components for with five different noise ratios. Monte Carlo sampling of the posterior distribution is shown for the noisy training data in Fig. 9. Clearly, uncertainty from the data due to the Gaussian white noise is well characterized by the ensemble of Monte Carlo sampling on the distribution. To further analyze the effect of the difference between the mean of the posterior distribution of the prediction and the ground truth clean trajectory on the flow field, we show the (projected) mean and standard deviation of the predicted POD coefficients at $t=100t_{ref}$, in Fig. 10 and Fig. 11. As seen in Fig. 10, there is hardly any difference between the mean of the posterior distribution and the ground truth. The contour of standard deviation projected onto the flow field shows a similar pattern to the vortex shedding, and the standard deviation near the wake region is relatively small compared to other domains in the flow field.

5 Conclusions

A probabilistic, stabilized deep learning framework was presented towards the end of extracting the Koopman decomposition for continuous nonlinear dynamical systems. We formulated the deep learning problem from a measure-theoretic perspective with a clear layout of all the assumptions. Two different forms: differential and recurrent, suitable for different situations were proposed and discussed. A parameterization of Koopman operator was proposed, with guaranteed stability. Further, a novel deep learning architecture was devised such that the SVD-DMD [70] is naturally embedded. Finally, mean-field variational inference was used to quantify uncertainty in the modeling. To evaluate the posterior distribution, Monte Carlo sampling procedures corresponding to different forms were derived. Finally, the model was evaluated on three continuous nonlinear dynamical systems ranging from toy polynomial problems to an unstable wake flow behind a cylinder, from three different aspects: the uncertainty with respect to the density of data in the domain, unseen data, and noise in the data. The results show that the proposed model is able to capture the uncertainty in all the above cases and is robust to noise.

Appendix A Real part of eigenvalues of a negative semi-definite real square matrix

Definition A.1**.**

*An $n\times n$ real matrix (and possibly non-symmetric) $\mathbf{A}$ is called negative semi-definite, if $\mathbf{x}^{\top}\mathbf{A}\mathbf{x}\leq 0$, for all non-zero vectors $\mathbf{x}\in\mathbb{R}^{n\times 1}$. *

Lemma A.2**.**

*For a negative semi-definite real square matrix, the real part of all of its eigenvalues is non-positive. *

Proof A.3**.**

Consider the general negative semi-definite matrix defined in Definition A.1 as $\mathbf{A}\in\mathbb{R}^{n\times n}$, such that $\mathbf{x}^{\top}\mathbf{A}\mathbf{x}\leq 0$ for any non-zero vector $\mathbf{x}\in\mathbb{R}^{n\times 1}$. Without loss of generality, denote its eigenvalues as $\lambda=\alpha+j\beta$, where $\alpha,\beta\in\mathbb{R}$ and corresponding eigenvectors as $\mathbf{v}=\mathbf{v}_{r}+j\mathbf{v}_{i}$ where $\mathbf{v}_{r},\mathbf{v}_{i}\in\mathbb{R}^{n}$.

Then we have $0=(\mathbf{A}-\lambda)\mathbf{v}=(\mathbf{A}-\alpha-j\beta)(\mathbf{v}_{r}+j\mathbf{v}_{i})$, which further leads to $(\mathbf{A}-\alpha)\mathbf{v}_{r}=-\beta\mathbf{v}_{i}$ and $(\mathbf{A}-\alpha)\mathbf{v}_{i}=\beta\mathbf{v}_{r}$. Then we have $\mathbf{v}_{r}^{\top}(\mathbf{A}-\alpha)\mathbf{v}_{r}=-\beta\mathbf{v}_{i}^{\top}\mathbf{v}_{r}$, and $\mathbf{v}_{i}^{\top}(\mathbf{A}-\alpha)\mathbf{v}_{i}=\mathbf{v}_{i}^{\top}\beta\mathbf{v}_{r}$. Thus adding the two equations, we have $\mathbf{v}_{r}^{\top}(\mathbf{A}-\alpha)\mathbf{v}_{r}+\mathbf{v}_{i}^{\top}(\mathbf{A}-\alpha)\mathbf{v}_{i}=0$, which implies $\alpha=(\mathbf{v}_{r}^{\top}\mathbf{A}\mathbf{v}_{r}+\mathbf{v}^{\top}_{i}\mathbf{A}\mathbf{v}_{i})/(\mathbf{v}_{r}^{\top}\mathbf{v}_{r}+\mathbf{v}_{i}^{\top}\mathbf{v}_{i})$.

*Using the definition of negative semi-definite matrices, we have $\mathbf{v}_{r}^{\top}\mathbf{A}\mathbf{v}_{r}\leq 0$ and $\mathbf{v}_{i}^{\top}\mathbf{A}\mathbf{v}_{i}\leq 0$, even if $\mathbf{v}_{r}$ or $\mathbf{v}_{i}=\mathbf{0}$. Since there is at least one non-zero vector between $\mathbf{v}_{r}$ and $\mathbf{v}_{i}$, one can safely arrive at $\alpha\leq 0$ that the real part of any eigenvalue of $\mathbf{A}$ is non-positive. *

Acknowledgements

This work was supported by DARPA under the grant titled Physics Inspired Learning and Learning the Order and Structure Of Physics,. Computing resources were provided by the NSF via grant 1531752 MRI: Acquisition of Conflux, A Novel Platform for Data-Driven Computational Physics.