1 Introduction↩︎

Modeling and forecasting the behavior of noisy, high-dimensional, and nonlinear systems remains an active area of research. Yet, the governing laws of these systems are often unknown or prohibitively expensive to simulate. With the proliferation of measurement sensors and the advent of efficient computational hardware, our modern world is abundant with multi-fidelity data that enables us to construct models directly from data without relying solely on first principles. Beyond merely generating predictions, it is also possible to extract interpretable spatio-temporal patterns and coherent structures from the data, thereby deepening our understanding and enhancing our analysis of complex dynamical systems.

Transfer operator theory and reproducing kernel Hilbert spaces (RKHSs) are popular approaches that result in analytically informative algorithms [1]–[9], while also contributing widely to practical applications [10]–[17]. Both Koopman operator and RKHS techniques frame their corresponding problems in a functional analytic setting. Therefore, it stands to reason that by combining methods from both of these two fields, we can design robust algorithms that provide compact and interpretable models for (stochastic) time-series.

1.1 Dynamic Mode Decomposition↩︎

Transfer operators, such as the Koopman and Perron–Frobenius operators, describe the theory of linearizing a dynamical system, not merely around its fixed points, but globally without any approximation. The operator-theoretic perspective is particularly appealing because it is integrated with well-understood statistical and geometric perspectives, which ties in with insights into the physical behavior of dynamical systems mezic2019spectrumkoopmanoperatorspectral?, koopman1931hamiltonian?, koopman1932dynamical?.

The Koopman operator, being an infinite-dimensional object, must often be empirically approximated before it can be implemented on a computer. This involves discretizing the operator and identifying a suitable set of basis functions that span a finite-dimensional subspace wherein the spectral properties of the operator can be represented — a task that falls under the heading of dictionary learning [18]–[24].

One such suite of algorithms that discretizes the Koopman operator is dynamic mode decomposition (DMD; [25]–[27]). DMD can be described as a Galerkin method that blends together principal component analysis (PCA) and the discrete-time Fourier transform [28]. In other words, one may think of DMD as rotating the (lower) dimensional PCA space such that each basis has dynamics oscillating according to a single frequency and a single growth (or decay) rate.

In part, DMD remains attractive because it is based on an eigenvalue decomposition of a model obtained through linear regression. After projecting the system to a higher-dimensional space, DMD can provide accurate models of nonlinear [26], [29] and stochastic systems [30]–[32], vcrnjaric2020koopman?. In other words, all the techniques belonging to linear time-invariant ODEs, such as spectral analysis, multi-step forecasting, numerical simulation, and even control [33], [34], peitz2019koopman?, otto2021koopman? can be readily transferred to nonlinear systems. Moreover, should the system contain a large number of variables (e.g., in fluid dynamics), DMD can facilitate dimensionality reduction.

1.2 Gaussian Process Regression↩︎

Following the authors [29], [35]–[40], das2020koopman?, ikeda2022koopman?, bevanda2024koopman?, bevanda2025koopman?, the underlying machinery of our DMD algorithm is a reproducing kernel function (i.e., a covariance function). Kernel methods are non-parametric, enhance interpretability of features, and are generally more flexible compared to parametric models (e.g., multilayer perceptrons), since they encompass all functions sharing the same degree of smoothness. It is the non-parametric nature of kernel algorithms that grant them the intrinsic ability to implicitly construct a functional basis for a model directly from data. Thereby, they provide a natural choice for the hypotheses space: an RKHS. It is these properties that often lead to strong theoretical performance guarantees [9], [41], [42].

One popular kernel regression framework (for modeling dynamical systems), is Gaussian process (GP) regression [4], [43]–[46]. GP regression is a function approximation technique that is capable of providing accurate estimations of an unknown function based on a limited set of noisy training targets and high-level assumptions about the function. Moreover, a GP model is numerically straightforward to estimate, the predictions of the model are fully probabilistic, and the theory provides a well-founded framework for model selection.

1.3 Contributions of this work↩︎

The core idea is to transfer the attractive features of GP regression to DMD. Thereby, we develop a Bayesian interpretation of the embedded Perron–Frobenius operator [39] which opens up avenues to extend kernel-based DMD algorithms with established sparsity and hyperparameter tuning techniques. The result is a modification of ‘Extended DMD’ (EDMD; [26]) that constructs models from noisy measurement data of a nonlinear dynamical system and provides uncertainty estimates for multi-step predictions.

The remainder of this paper is structured as follows. Section 2 introduces a Bayesian formulation of DMD, and Section 3 explores its mathematical consequences and numerical performance. Lastly, in Section 4, we conclude with some final remarks, open questions, and point toward straightforward extensions. In Appendix 5 we review (sparse) GP regression and the notation of [4], and for convenience and completeness we also present an overview of (kernel) transfer operators [39] in Appendix 6.

2 Gaussian Process Dynamic Mode Decomposition↩︎

Because kernel-based DMD algorithms often quickly become impractical for large data sets due to the high training cost, one of our goals is to apply the variational free energy (VFE) method [47] to kernel Koopman operators [39]. For example, in [39], to estimate the Koopman matrix we need to take the inverse of a Gramian which has a complexity of \(\mathcal{O}(N^3)\), where \(N\) is the number of training samples; moreover, the eigenvalue decomposition also has a complexity of \(\mathcal{O}(N^3)\). Therefore, we are motivated to find a better trade-off between the model’s generalization error and its computational complexity. We want a non-parametric model that automatically adapts to the complexity of the dynamics, but not to the point where all the training inputs define the basis dictionary.

We will demonstrate how our algorithm is a special case of EDMD [26]. Which means that it inherits some of EDMD’s limitations [48]–[50]. This should not be surprising, since it is quite common to derive the same expressions within either a frequentist or Bayesian framework. This duality between ridge-regression and Bayesian models is well-known, and has been pointed out in [4], [51]–[53]. One might then assume that the choice of framework merely reflects an attitude, with each attitude corresponding to a preferred mode of reasoning about how to construct a (data-driven) model. However, we could argue that attitudes matter because they suggest different approaches to improving upon the state-of-the-art. Hence, one of our goals is to shed some light on this debate.

For example, EDMD does not intrinsically provide an optimization scheme for selecting either the hyperparameters or the dictionary, and therefore requires another explicit strategy, e.g., [54]. Neither does EDMD provide uncertainty-aware predictions, nor does it take into account measurement noise from the outset. Our amalgamation of DMD and GP regression wraps up the entire model inference pipeline into one coherent theory that simultaneously addresses hyperparameter optimization and sparse dictionary learning. This naturally raises the question of whether such an approach is overly restrictive and sacrifices the flexibility of frequentist methods.

2.1 Operators as Random Variables↩︎

r0.5 4pt

We start with a full Bayesian treatment of the embedded Perron–Frobenius operator (PFO) [39], which we denote by \(\mathcal{P}_\varepsilon : \mathbb{H} \rightarrow \mathbb{H}\). The embedded Perron–Frobenius operator, is the adjoint of the kernel Koopman operator. Unlike the PFO which propagates densities, the embedded PFO can be thought of as a push-forward map that acts on representations of probability densities. In short, the idea of feature maps is extended to the space of probability distributions by thinking of \(\varphi(\boldsymbol{x}) \, \forall \, \boldsymbol{x} \in \mathbb{X}\) as an injective representation of a density function. See Appendix 7 for a more detailed discussion of kernel mean embeddings.

Let us consider the ‘implicitly’ lifted training data set obtained from a dynamical system as \(\left(\Phi_X, \Phi_Y\right) =\{ \left(\varphi\left(\boldsymbol{x}_i\right), \varphi\left(\boldsymbol{y}_i\right) \right)\}_{i = 1}^N\), where the targets \(\{\boldsymbol{y}_i\}_{i=1}^N\) are produced by the discrete-time flow map \(\boldsymbol{F}: \mathbb{X} \rightarrow \mathbb{X}\), and \(\mathbb{X} \subseteq \mathbb{R}^D\). We further assume only isotropic measurement noise: \[\begin{align} {3} Y_i = \boldsymbol{F}(X_i) + \boldsymbol{\epsilon}_i, \quad \text{where } \boldsymbol{\epsilon}_i \sim \mathcal{N}^D(\boldsymbol{0}, \sigma_Y^2 \mathbf{I}_D). \end{align}\]

In Appendix 5, we reviewed Gaussian process regression with outputs in \(\mathbb{R}\). We now extend this setting by taking the target space to be the RKHS \(\mathbb{H}\). This presents a challenge of handling not just finitely many outputs, but an infinite number of output channels. To address this, we employ an operator-valued kernel, and subsequently derive the posterior GP associated with the embedded PFO.

For the sake of clarity, allow a brief oversimplification and consider the \(D\! =\! 1\) case. In this setting, the feature map \(\varphi(\boldsymbol{x}) \! \mathrel{\vcenter{:}}= \kappa_\textrm{pr}(\boldsymbol{x},\cdot) \, \in \, \mathbb{H}\) is defined by the GP prior over a one-dimensional flow map, i.e. a real- and scalar-valued kernel \(\kappa_\textrm{pr} : \mathbb{X} \times \mathbb{X} \rightarrow \mathbb{R}\).

The purpose of the Bayesian-consistency kernel \(\kappa_\textrm{bc}: \mathbb{X} \times \mathbb{X} \rightarrow \mathbb{R}\) is to relate the observation model in the state space to its counterpart in the RKHS. It models the correlations between the ‘components’ of \(\varphi(\cdot)\) such that consistency is maintained with the flow map, i.e., ensuring that that \(\boldsymbol{F}(\cdot)\) is an observable of the Koopman operator. The random variable \(\nu_i \in \mathbb{H}\) is defined by the covariance operator \(\mathcal{C}_\textrm{bc}: \mathbb{H} \rightarrow \mathbb{H}\), which in practice is required to be self-adjoint and invertible, and represents the statistical properties of the noisy full-state observable within the RKHS. For more details refer to Appendix 9.

Similar to scalar-valued GP regression, our prior beliefs about the statistical properties of the operator \(\mu_{Y \mid \cdot}: \mathbb{X} \rightarrow \mathbb{H}\) are represented by a GP, \(\mu_{Y \mid \cdot} \sim \mathcal{GP}^\infty \left(\mathcal{M}_\textrm{pr}, {\xi}_\textrm{pr} \right)\)¹. However, here we are defining the prior covariance \({\xi}_\textrm{pr} \colon \mathbb{X} \times \mathbb{X} \to \mathbb{H} \, \otimes \, \mathbb{H}\) with an operator-valued kernel [55]. Then, without loss of generality, we choose the mean operator \(\mathcal{M}_\textrm{pr}:\mathbb{X} \rightarrow \mathbb{H}\) such that it produces an identically zero function. It follows that the prior distribution over the latent operator evaluations \(\mu_{Y \mid \mathbf{X}} \mathrel{\vcenter{:}}= \mathcal{P}_\varepsilon \; \Phi_X\in \mathbb{H}^{1 \times N}\) is \[\begin{align} {3} p_\textrm{pr}(\mu_{Y \mid \mathbf{X}}) =\mathcal{GP}^N\left(\mu_{Y \mid \mathbf{X}}; \mathcal{M}_\textrm{pr}(\mathbf{X}), \Xi_{XX} \right), \end{align}\] where \(\Xi_{XX}\) is the prior covariance evaluated at \(\left(\mathbf{X}, \, \mathbf{X}\right)\) and hence an element of \((\mathbb{H} \, \otimes \, \mathbb{H})^{N \times N}\).

The Riesz—Fréchet representation theorem [56], [57] says that for each \(\boldsymbol{x}\! \in\! \mathbb{X}\) and \(g \! \in \! \mathbb{H}\), there is an element \(\Xi(\boldsymbol{x})g\! : \mathbb{X} \rightarrow \mathbb{H}\), such that the reproducing property is \[\begin{align} {3} \langle \mu_{Y \mid \cdot}, \Xi(\boldsymbol{x})g \rangle_{\mathbb{H}^{\mathbb{X}}} = \langle \mu_{Y \mid \boldsymbol{x}}, g\rangle_{\mathbb{H}}. \end{align}\] This implies that the feature map \(\Xi(\boldsymbol{x})\) is a bounded linear operator, \(\Xi(\boldsymbol{x}) \in \mathfrak{B}(\mathbb{H}, \mathbb{H}^{\mathbb{X}})\), where \(\mathbb{H}^{\mathbb{X}}\) denotes the \(\mathbb{H}\)-valued RKHS on \(\mathbb{X}\) (i.e., \(\mu_{Y \mid \cdot} \in \mathbb{H}^{\mathbb{X}}\)). Thereby, the relationship between the reproducing kernel \(\xi_\textrm{pr}(\cdot, \cdot)\) and its associated feature map \(\Xi(\boldsymbol{x})\) is \[\begin{align} {3} \langle \Xi(\boldsymbol{x}) \, g, \, \Xi(\boldsymbol{x}') \, g'\rangle_{\mathbb{H}^{\mathbb{X}}} = \langle g , \, \xi_\textrm{pr}(\boldsymbol{x}, \boldsymbol{x}') \, g'\rangle_{\mathbb{H}}. \end{align}\] The entry at \(\left(\xi_\textrm{pr}\left(\boldsymbol{x}, \boldsymbol{x}'\right)\right)_{\boldsymbol{y},\boldsymbol{y}'}\) determines the prior covariance between \(\mu_{Y \mid \boldsymbol{x}}(\boldsymbol{y})\) and \(\mu_{Y \mid \boldsymbol{x}'}(\boldsymbol{y}')\), i.e., it expresses the degree to which the response at \((\boldsymbol{x}, \boldsymbol{y})\) is affected by \((\boldsymbol{x}', \boldsymbol{y}')\) — before incorporating data.

Definition 1 models the lifted targets as scalar-valued GPs, and therefore the likelihood can be expressed as \(p(\varphi(\boldsymbol{y}) \mid \mathcal{P}_\varepsilon, \boldsymbol{x}) = \mathcal{GP}(\mu_{Y \mid \boldsymbol{x}}, \sigma_{Y}^2\kappa_\textrm{bc})\). Thereby, the posterior retains the same structure as in the scalar-valued case [58]–[60]: \[\begin{align} {3} \mathcal{M}_\textrm{pst}(\boldsymbol{x}) &\mathrel{\vcenter{:}}= \mathcal{A} \; \xi_{X}({\boldsymbol{x}})= \sum_{i = 1}^N \mathcal{A}_i \;{\xi}_\textrm{pr} (\boldsymbol{x}_i, \boldsymbol{x}), \label{eq:32GP32mean32operator}\\ {\xi}_\textrm{pst}(\boldsymbol{x}, \boldsymbol{x}') &\mathrel{\vcenter{:}}= {\xi}_\textrm{pr}(\boldsymbol{x}, \boldsymbol{x}') - \mathbf{\xi}_{X}^\top(\boldsymbol{x}) \text{ } \tilde{\Xi}_{XX}^{-1}\text{ }\mathbf{\xi}_{X}(\boldsymbol{x}'). \end{align}\tag{1}\] We unpack these equations by concatenating the noisy targets ‘end-to-end’ in \(\Upsilon^\top \! \mathrel{\vcenter{:}}= \! \operatorname{vec}(\Phi_Y)\! \in \!\mathbb{H}^{\oplus N}\). We denote the weight operator by \(\mathcal{A} \mathrel{\vcenter{:}}= \Upsilon \; \tilde{\Xi}_{XX}^{-1}\), the noisy covariance operator by \(\tilde{\Xi}_{XX} \mathrel{\vcenter{:}}= \Xi_{XX} + \Sigma\), and the feature “vector" by \(\mathbf{\xi}_{X}(\boldsymbol{x}) \mathrel{\vcenter{:}}= [{\xi}_\textrm{pr}(\boldsymbol{x}_1, \boldsymbol{x}), \ldots, {\xi}_\textrm{pr}(\boldsymbol{x}_N, \boldsymbol{x})]^\top\).

We will elaborate on this point in Section 2.4, but the structure the operator-valued kernel is restricted to \(\xi_\textrm{pr}(\boldsymbol{x}, \boldsymbol{x}') \! = \! \kappa_\textrm{pr}(\boldsymbol{x}, \boldsymbol{x}') \, \mathcal{C}_\textrm{bc}\). This form allows us to express the covariance operators in terms of the prior covariances over latent values of the flow map and Kronecker tensor products with \(\mathcal{C}_\textrm{bc}\), i.e., \(\Xi_{XX} \mathrel{\vcenter{:}}= \mathbf{K}_{XX} \otimes \mathcal{C}_\textrm{bc}\) where \(\mathbf{K}_{XX} \mathrel{\vcenter{:}}= \Phi_X^\top \, \Phi_X \in \mathbb{R}^{N \times N}\). Similarly, the “sensor" noise covariance operator is denoted by \(\Sigma \mathrel{\vcenter{:}}= \mathbf{I}_N \, \otimes \, \sigma_{Y}^2\,\mathcal{C}_\textrm{bc}\), and the feature vector is \(\mathbf{\xi}_{X}(\boldsymbol{x}) = \mathbf{k}_{X}(\boldsymbol{x}) \otimes \mathcal{C}_\textrm{bc}\) where \(\mathbf{k}_{X}(\boldsymbol{x}) \mathrel{\vcenter{:}}= \Phi_X^\top \varphi(\boldsymbol{x})\). Thereby, in essence we are implementing an intrinsic model of coregionalization [58].

Then by recalling the mixed matrix-vector product property: \((\mathbf{A}\otimes\mathbf{B})\, \operatorname{vec}(\mathbf{C}) = \operatorname{vec}(\mathbf{BCA}^\top)\), the posterior mean 1 of the embedded PFO reduces to \[\begin{align} {3} \label{eq:32embedded32PFO} \mathcal{M}_\textrm{pst}(\boldsymbol{x})\! &= (\mathbf{k}^\top_{X}(\boldsymbol{x}) \! \otimes \! \mathcal{C}_\textrm{bc}) (\mathbf{K}_{XX} \! \otimes \! \mathcal{C}_\textrm{bc} + \mathbf{I}_N \! \otimes \! \sigma_{Y}^2\mathcal{C}_\textrm{bc})^{-1} \Upsilon^\top \nonumber\\ &= \left(\mathbf{k}_{X}^{\top}(\boldsymbol{x}) \,\tilde{\mathbf{K}}_{XX}^{-1} \otimes \mathcal{I} \right) \operatorname{vec}(\Phi_Y)\\ &= \Phi_Y \tilde{\mathbf{K}}_{XX}^{-1} \Phi_X^\top \; \varphi(\boldsymbol{x}) = \mathbb{E}\left[ {\mathcal{P}}_\varepsilon \, \varphi(X) \mid X = \boldsymbol{x}\right].\nonumber \end{align}\tag{2}\] Hence, we recover the same expression for \({\mathcal{P}}_\varepsilon\) derived in [39], while also mirroring the frequentist approach of [61]. In the expression of the mean, the operator \(\mathcal{C}_\textrm{bc}\) plays no role whatsoever, yet we observe that the posterior covariance over \({\mathcal{P}}_\varepsilon\varphi(\boldsymbol{x})\) retains the same form as the prior: \(\xi_\textrm{pst}(\boldsymbol{x}, \boldsymbol{x}') = \kappa_\textrm{pst}(\boldsymbol{x}, \boldsymbol{x}') \, \mathcal{C}_\textrm{bc}\), i.e., a linear transformation of the covariance over the flow map.

In the case where \(D > 1\), we may want to treat the prior and Bayesian-consistency kernels as matrix-valued and include anisotropic measurement noise [58], [62]. However, for the rest of this manuscript we will continue assuming that the channels in the flow map are unrelated and that the sensor noise is isotropic.

2.2 Sparse Variational Bayesian Transfer Operators↩︎

Compressing the exact posterior to a smaller basis dictionary \(\Phi_Z \mathrel{\vcenter{:}}= \{ \varphi\left(\boldsymbol{z}_i \right)\}_{i = 1}^M\in \mathbb{H}^{1 \times M}\) with \(M < N\), is possibly the easiest step, and demonstrates how straightforward it is to apply Gaussian process techniques to kernel transfer operators. Because, once we have framed the problem within a Bayesian framework, we can simply apply the variational free energy (VFE) method (see Appendix 5.3 and [47]): \[\begin{align} {3} \label{eq:32sparse32embedded32PFO} \begin{aligned} \mathcal{M}_\textrm{pst}(\boldsymbol{x}) &= \xi_{Z}^\top({\boldsymbol{x}}) \left( \Xi_{ZX} \Xi^\top_{ZX} + \Sigma\,\Xi_{ZZ} \right)^{-1} \Xi_{ZX}\Upsilon^\top \\ &= \left(\mathbf{k}^\top_{Z}(\boldsymbol{x}) \,\tilde{\mathbf{C}}_{XX}^{-1} \otimes \mathcal{I} \right) \, \operatorname{vec}\left(\Phi_Y \mathbf{K}_{ZX}^\top \right)\\ &= \Phi_Y \mathbf{K}_{ZX}^\top \tilde{\mathbf{C}}_{XX}^{-1} \Phi_Z^\top \, \varphi(\boldsymbol{x}), \end{aligned} \end{align}\tag{3}\] where we have the compressed covariance operators \(\Xi_{ZX} \, \mathrel{\vcenter{:}}= \, \mathbf{K}_{ZX} \, \otimes \, \mathcal{C}_\textrm{bc}\) and \(\Xi_{ZZ} \, \mathrel{\vcenter{:}}= \, \mathbf{K}_{ZZ} \, \otimes \, \mathcal{C}_\textrm{bc}\). The sparse covariance matrices \(\mathbf{K}_{ZX} \, \mathrel{\vcenter{:}}= \, \Phi_Z^\top \, \Phi_X, \; \mathbf{K}_{ZZ} \, \mathrel{\vcenter{:}}= \, \Phi_Z^\top \, \Phi_Z\) allow us to define the noisy Gramian \(\tilde{\mathbf{C}}_{XX} \, \mathrel{\vcenter{:}}= \,\mathbf{K}_{ZX} \, \mathbf{K}_{ZX}^\top \, + \, \sigma^2_Y \mathbf{K}_{ZZ}\). Lastly the feature vectors are \(\mathbf{k}_{Z}(\boldsymbol{x}) \! \mathrel{\vcenter{:}}= \! \Phi_Z^\top \varphi(\boldsymbol{x})\) and \(\xi_{Z}({\boldsymbol{x}}) \! \mathrel{\vcenter{:}}= \! \mathbf{k}_{Z}(\boldsymbol{x}) \otimes \mathcal{C}_\textrm{bc}\).

In addition to providing a means of learning a sparse dictionary, the VFE method offers convergence-rate guarantees [63] and can also be used for hyperparameter optimization. Hyperparameter tuning is another obstacle for kernel-based Koopman methods, since one should ideally design the reproducing kernel function such that the Koopman operator is represented in an RKHS that is both densely defined and closable ikeda2022koopman?. Consequently, finding such a kernel function tends to be a delicate and time-consuming task.

2.3 Koopman Mode Decomposition↩︎

To obtain an explicit, numerically tractable expression for the kernel Koopman matrix we first assume that any observable function \(g \in \mathbb{H}\) can be approximated with a negligible error (i.e., \(\mathbb{H}\) is an invariant subspace of the Koopman operator). We then proceed by taking the adjoint of the embedded PFO in 3 to find the posterior expression for the kernel Koopman operator. Discretizing this operator yields the Koopman matrix (see 21 in Appendix 6.3): \[\begin{align} {3} \begin{aligned} \label{eq:32CME32matrix} {\mathbf{U}} &= \tilde{\mathbf{C}}_{XX}^{-1} \mathbf{C}_{XY}\in \mathbb{R}^{M \times M}. \end{aligned} \end{align}\tag{4}\]

Notice that since \(\mathbf{C}_{XY} \mathrel{\vcenter{:}}= \mathbf{K}_{ZX} \mathbf{K}_{ZY}^\top\) and \(\mathbf{K}_{ZY} \mathrel{\vcenter{:}}= \Phi_Z^\top \, \Phi_Y\), equation 4 recovers the standard algorithm for EDMD [2]. In other words, the Bayesian approach we have followed produces the same expression as the least-squares approach of EDMD with a Tikhonov parameter, \(\sigma_{Y}^2\). Of course, since EDMD is valid for more general choices of dictionaries, the previous statement is only true for a kernel-based dictionary [26].

Provided that the eigenvalues are non-degenerate and unique, the eigenvalue decomposition of \({\mathbf{U}}\) determines the coefficients of the spectral and spatial components: \[\begin{align} \label{eq:32eigendecomp} {\mathbf{U}} = \mathbf{W}\; \mathbf{\Lambda} \; {\mathbf{V}}_{\kappa}, \end{align}\tag{5}\] where \(\mathbf{\Lambda}\in \mathbb{C}^{M \times M}\) is a diagonal matrix containing the discrete-time eigenvalues. The right eigenvectors \(\mathbf{W}\in \mathbb{C}^{M \times M}\) contain the projection vectors onto the eigenfunction basis, whereas the left eigenvectors \({\mathbf{V}}_{\kappa} = \mathbf{W}^{-1}\) correspond to the modes of the observables \(\mathbf{k}_Z(\cdot) = [\kappa_\textrm{pr}(\boldsymbol{z}_1, \cdot), \, \ldots, \, \kappa_\textrm{pr}(\boldsymbol{z}_M, \cdot)]^\top\).

The projected modes with respect to the full-state observables, \({\mathbf{V}_f} \in \mathbb{C}^{M \times D}\), are found via \[\begin{align} {3} \label{eq:32GP32DMD32modes} \begin{aligned} {\mathbf{V}}_{f} &= \mathbf{\Lambda}^{-1} \, {\mathbf{V}}_{\kappa} \, \tilde{\mathbf{C}}_{XX}^{-1} \, \mathbf{K}_{ZX} \mathbf{Y}^\top. \end{aligned} \end{align}\tag{6}\] The implication of 6 is that for a single-step prediction, the model in 9 is also the posterior mean of a GP prior placed on the flow map (using the dynamically evolved states \(\mathbf{Y}\) as the targets). The point is that we can perform DMD on a random variable associated with the Koopman operator while maintaining consistency with the nonlinear dynamics on the state manifold.

In Appendix 8, we present the numerically stable and efficient algorithm (GP-TCCA) used in the experiments of Section 3 for computing equations 4 –6 .

2.4 Prediction Theory↩︎

To construct a procedure by which we can produce multi-step forecasts and quantify the model’s confidence in those future estimates, we interpret the posterior covariance as heteroskedastic process noise (which naturally decreases as the size of the training set grows). Put differently, the Bayesian Koopman operator is a random variable which represents stochastic knowledge-dependent dynamics. Thereby, since the projected dynamical system is linear in \(\mathbb{R}^M\), the time series \(\{\mathbf{k}_Z(X_i)\}_{i=0}^{k}\) of random vectors can be described by a stochastic difference equation [64]: \[\begin{align} {3} \label{eq:32LSDE} \mathbf{k}_Z(X_{k}) = \mathbf{A}_\kappa^1 \; \mathbf{k}_Z(X_{k-1}) \, + \, {\mathbf{L}}_\textrm{pst}(X_{k-1})\boldsymbol{\omega}_{k}, \end{align}\tag{7}\] where \(\boldsymbol{\omega}_k \sim \mathcal{N}^M(\boldsymbol{0}, \mathbf{I}_M)\) and the multi-step transition matrix is \(\mathbf{A}_\kappa^k \mathrel{\vcenter{:}}= \left( \mathbf{W}\,\mathbf{\Lambda}^k \, {\mathbf{V}}_{\kappa} \right)^\top\). The process \(\{ \boldsymbol{\omega}_{i} \}_{i=1}^{k}\) is an i.i.d. martingale difference sequence adapted to the sequence of information sets \(\{J_i\}_{i=0}^{k}\). Each information set \(J_k\) represents all the information available to an agent at time-step \(k\), and contains all measurable functions of \(\{X_0, \boldsymbol{\omega}_1, \ldots, \boldsymbol{\omega}_k \}\).

The lower triangular matrix \({\mathbf{L}}_\textrm{pst}(X_k)\) is determined by the Cholesky decomposition of \({\mathbf{\Xi}}_{\textrm{pst}}^{1}(X_k)\). This covariance matrix encodes the one-step-ahead forecasting error of \(\mathbf{k}_Z(X_{k})\) conditioned on \(J_{k-1}\). Given the operator-valued kernel defined earlier, in the current setting this matrix is \({\mathbf{\Xi}}_{\textrm{pst}}^{1}(X_k) = {\kappa}_\textrm{pst}(X_k, X_k) \, \mathbf{K}_\textrm{bc}\). Now recall that the expression for the sparse posterior kernel (see 17 in Appendix 5.3) is \[\begin{align} {3} \label{eq:32posterior32kernel32VFE} {\kappa}_\textrm{pst}(\boldsymbol{x}, \boldsymbol{x}') &= \kappa_\textrm{pr}(\boldsymbol{x}, \boldsymbol{x}') - \mathbf{k}_{Z}^\top(\boldsymbol{x}) \text{ }\tilde{\mathbf{B}}\text{ }\mathbf{k}_{Z}(\boldsymbol{x}') \end{align}\tag{8}\] where \(\tilde{\mathbf{B}} \mathrel{\vcenter{:}}= \left(\mathbf{K}_{ZZ}^{-1} -\sigma_Y^2 \tilde{\mathbf{C}}_{XX}^{-1}\right)\). Hence, the single-step prediction error is represented by the discretized version of the posterior kernel \(\xi_\textrm{pst}(\cdot, \cdot)\) and subsequently the operator \(\mathcal{C}_\textrm{bc}\) (see Appendix 9 for details).

Expanding 7 into a recursive expression that is a function of \(\{X_{k-j-1}, \boldsymbol{\omega}_{k-j} \}_{j = 0}^{k-1}\) results in \[\begin{align} {3} \begin{aligned} \mathbf{k}_Z(X_{k}) &= \mathbf{A}_\kappa^k\; \mathbf{k}_Z(X_{0}) + \sum_{j=0}^{k-1} \mathbf{A}_\kappa^j\; {\mathbf{L}}_\textrm{pst}(X_{k-j-1})\boldsymbol{\omega}_{k-j}. \end{aligned} \end{align}\] Therefore, the mean vector of \(X_{k}\) given the realization \(\boldsymbol{x}_{0}\) is \[\begin{align} {3} \label{eq:32GP32DMD} \begin{aligned} \hat{\boldsymbol{x}}_k &= \mathbb{E}\left[X_{k} \mid X_0 =\boldsymbol{x}_0\right] =\mathbf{A}_f^k \; \mathbf{k}_Z(\boldsymbol{x}_{0}), \end{aligned} \end{align}\tag{9}\] where the weight matrix is \(\mathbf{A}_f^k \mathrel{\vcenter{:}}= \left( \mathbf{W}\,\mathbf{\Lambda}^k \, {\mathbf{V}}_{f} \right)^\top \in \mathbb{R}^{D \times M}\). In deriving 9 we used 6 , the relationship \(X_k = (\mathbf{W}{\mathbf{V}}_{f} )^\top \mathbf{k}_Z(X_{k})\), and the fact that a linear transformation of a Gaussian random vector yields another Gaussian random vector [65]. This last property was the key to calculating the Bayesian-consistency matrix \(\mathbf{K}_\textrm{bc} \in \mathbb{R}^{M \times M}\) and enforcing the structure of the operator–valued prior covariance function to \(\xi_\textrm{pr}(\boldsymbol{x}, \boldsymbol{x}') \! = \! \kappa_\textrm{pr}(\boldsymbol{x}, \boldsymbol{x}') \, \mathcal{C}_\textrm{bc}\).

The consistency between the generative models of the lifted space and the state space, rests on the knowledge that the following expression has to hold when linearly projecting the posterior covariance of the Gaussian random vectors: \[\begin{align} {3} \label{eq:32consistency32check} \begin{aligned} \mathbf{K}_\textrm{pst}^k(\boldsymbol{x}) + \sigma_Y^2 \mathbf{I}_D &= \mathbf{A}_f^0 \left({\mathbf{\Xi}}_{\textrm{pst}}^{k}(\boldsymbol{x}) + \sigma_Y^2 \mathbf{K}_\textrm{bc} \right) \left(\mathbf{A}_f^0\right)^\top. \end{aligned} \end{align}\tag{10}\] Given the independence assumption we placed on the components of the flow map, the single-step covariance matrix with respect to the full-state observables is \(\mathbf{K}_\textrm{pst}^1(\boldsymbol{x}) = \kappa_\textrm{pst}(\boldsymbol{x}, \boldsymbol{x}) \, \mathbf{I}_D\). Hence, it follows that discretizing the operator \(\mathcal{C}_\textrm{bc}\) leads to \[\begin{align} {3} \begin{aligned} \mathbf{K}_\textrm{bc} &= \left( \left(\mathbf{A}_f^0\right)^\top \mathbf{A}_f^0 \right)^{\dagger} = {\mathbf{V}}_{\kappa}^* \left(\mathbf{V}_f \mathbf{V}_f^* \right)^{\dagger} {\mathbf{V}}_{\kappa}. \end{aligned} \end{align}\] Where the Hermitian adjoint is indicated by \(^*\), and since the matrix \(\mathbf{K}_\textrm{bc}\) always has a rank less than or equal to \(D\), we employed the Moore–Penrose pseudoinverse \(^\dagger\).

Now we are in a position to delve into the theory that describes how to propagate the covariance matrix \(\mathbf{\Xi}_{\textrm{pst}}^{k}(\boldsymbol{x}_0)\). Thus, we start with the expression for the \(k\)-step prediction error: \[\begin{align} {3} \zeta_{k} \mathrel{\vcenter{:}}= \mathbf{k}_Z(X_{k}) - \mathbf{k}_Z(\hat{\boldsymbol{x}}_k) = \sum_{j=0}^{k-1}\mathbf{A}_\kappa^j \;{\mathbf{L}}_\textrm{pst}(X_{k-j-1})\boldsymbol{\omega}_{k-j}, \end{align}\] which results in the propagated covariance matrix with respect to the observables \(\mathbf{k}_Z(\cdot)\): \[\begin{align} {3} \mathbf{\Xi}_{\textrm{pst}}^{k}(X_0) &= \mathbb{E}\left[ \zeta_{k} \, \zeta_{k}^{\top} \mid J_k\right]\\ &= \mathbb{E}\left[ \sum_{j = 0}^{k-1} \mathbf{A}_\kappa^j \; \mathbf{\Xi}_{\textrm{pst}}^{1}\left( X_{k-j-1}\right) \; \left(\mathbf{A}_\kappa^j\right)^\top \mid J_k \right]. \end{align}\] As this expression is in general not analytically tractable, similar to [44], [66] we can approximate it by performing Monte Carlo sampling or, if efficiency is of greater concern than accuracy, a second-order Taylor series expansion about the predicted mean²: \[\begin{align} {3} \mathbb{E}\left[ \mathbf{\Xi}_{\textrm{pst}}^{1}\left( X_{k}\right) \mid \boldsymbol{x}_0\right] &\approx \mathbf{\Xi}_{\textrm{pst}}^{1}\left( \hat{\boldsymbol{x}}_{k}\right) + {\mathbf{H}}^k \left( \boldsymbol{x}_0 \right), \end{align}\] where we have defined the curvature correction term as \[\begin{align} {3} {\mathbf{H}}^k\left(\boldsymbol{x}_0\right) &\mathrel{\vcenter{:}}= \frac{1}{2} \operatorname{Tr}\left( \nabla^2 \left( \hat{\boldsymbol{x}}_k \right) {\mathbf{K}}_{\textrm{pst}}^{k}\left({\boldsymbol{x}}_{0}\right)\right) \mathbf{K}_\textrm{bc}, \end{align}\] and \(\nabla^2(\boldsymbol{x})\) as the \(D \times D\) Hessian matrix³ of \({\kappa}_\textrm{pst}(\boldsymbol{x}, \boldsymbol{x})\).

Putting it all together we can obtain a recursive expression for the \(k\)-step covariance matrix: \[\begin{align} {3} \label{eq:32GP32DMD32cov} \mathbf{\Xi}_{\textrm{pst}}^{k}(\boldsymbol{x}_0) &\approx \mathbf{\Xi}_{\textrm{pst}}^{1}\left(\hat{\boldsymbol{x}}_{k-1}\right) + {\mathbf{H}}^{k-1}(\boldsymbol{x}_{0}) + \mathbf{U}^\top\mathbf{\Xi}_{\textrm{pst}}^{k-1}\left( \boldsymbol{x}_{0} \right)\mathbf{U}, \end{align}\tag{11}\] which we can interpret as iteratively accumulating the prediction errors arising from random ‘shocks’ at each intermediate step, and where the magnitude of these shocks are dependent on the system’s state.

If all the eigenvalues of \(\mathbf{U}\) have magnitudes less than \(1\), then the effect of the shocks will diminish as the errors are propagated (i.e. the process is covariance-stationary); for eigenvalues with magnitudes equal to \(1\) the effect of the shocks will persist within the process, and for eigenvalues with amplitudes greater than \(1\) the effects are amplified. However, we should be cognizant of the fact that the Koopman operator is a positive Markov operator [1], and will always have one “trivial" mode with a unity eigenvalue and a constant-one eigenfunction as a fixed point.

Additionally, it is straightforward to adapt 11 such that it is valid for arbitrary choices of observables. For example, the expression, \(\mathbf{W}^{*}\mathbf{\Xi}_{\textrm{pst}}^{k}(\boldsymbol{x}_0)\mathbf{W}\), is the \(k\)-step covariance matrix of the eigenfunctions. This ability to compute the prediction error for eigenfunctions are showcased in Figures 5 and 6, while the propagated errors are illustrated in Figure 1.

2.5 Model Selection↩︎

Hyperparameter optimization proceeds as if one were only interested in learning the optimal feature space for the flow map with the VFE method [47]. This is an optimization problem that is solved by maximizing the sum of the VFEs of each of the components of the flow map \(\boldsymbol{F} = \{ f_i \}_{i = 1}^D\). The reasoning is that the theoretical Bayesian-consistent structure of our model also ensures a coherent parameterization of the kernel Koopman and Perron–Frobenius operators.

It is important to recognize that optimizing the pseudo-inputs \(\mathbf{Z}\) even with analytic gradients can be computationally overwhelming. As a compromise, we recommend a layered greedy approach starting with the active learning methods of [69]–[71]. These methods are computationally very efficient, since they only require rank-one updates. While these approaches focus more on finding a good ‘spread’ of dictionary functions, the VFE method, in comparison, learns a dictionary that shapes the feature space to suit the observed targets.

Thereby, we used the approximate linear dependence (ALD) selection criterion [70] as a first pass to determine a subset of candidate locations for the pseudo-inputs from the training inputs. Since the ALD approach is a local ‘whack-a-mole’ active learning scheme, we followed it up by implementing active learning Cohn [72]: a more globally oriented methodology that approximates a full A-optimal design [73]. As a last step, only after a predefined number of new dictionary functions had been identified was the VFE hyperparameter tuning algorithm activated, with the greedy results given as a warm start⁴. This process is interleaved with the adaptation of generative hyperparameters, and repeats until no more pseudo-inputs are added to the basis dictionary or a threshold is passed, e.g., an upper limit on the dictionary size.

3 Numerical Experiments↩︎

With the foundational unification of GP regression and DMD complete, we present the results obtained from a few experiments. In all of the experiments, we used the \(\text{Matérn}^{5/2}\) kernel [76] with automatic relevance determination (ARD; [4]): \[\begin{align} {3} \kappa(\boldsymbol{x}, \boldsymbol{x}') = \sigma_f^2 \, \left(1 + \sqrt{5}\; r + \frac{5}{3}r^2\right) \text{exp}\left(-\sqrt{5}\; r\right). \end{align}\] where \(r^2 = (\boldsymbol{x} - \boldsymbol{x}')^\top \mathbf{T}^{-1} (\boldsymbol{x} - \boldsymbol{x}')\). The matrix \(\mathbf{T} \! \succ \! 0\) is diagonal and contains generative hyperparameters referred to as the characteristic lengthscales. These variables effectively prioritize the different components of \(\boldsymbol{x}\) and \(\boldsymbol{x}'\) and quantify how far it is needed to move along a particular axis for the function values to become uncorrelated. The hyperparameter \(\sigma_f^2\), called the signal variance, represents the average distance squared of the regression function from its mean.

For all numerical experiments, the data was standardized, and we emphasize that the data collection setup was asymmetric with respect to measurement noise. Since our model is not theoretically equipped to account for sensor noise on the inputs, we restricted the corruption to the dynamically evolved states. This simplification is admissible in surrogate modeling, but unrealistic in the majority of physical data acquisition contexts.

3.1 Generalization Errors↩︎

We compared the multi-step errors of our algorithm (referred to as GP-TCCA) using the Van der Pol oscillator as a benchmark: \[\begin{align} {3} \begin{aligned} \dot{x}_{t,1} &= x_{t,2},\\ \dot{x}_{t,2} &= \alpha (1 - x_{t,1}^2) x_{t,2} - x_{t,1}, \end{aligned} \end{align}\] where \(\alpha = 2\) controls the degree of nonlinearity. The dynamical system was simulated by a numerical ODE solver and sampled in time intervals of \(\Delta t = 50.0 \text{ [ms]}\). Figure 1 depicts an example trajectory of the oscillator and its associated forecast.

The metric used to quantify the generalization errors over a test set of \(N_\textrm{tst} = 5000\) noise-free samples was the symmetric mean absolute percentage error (SMAPE; [77]): \[\begin{align} {3} SMAPE = 100\% \frac{3}{N_\textrm{tst}} \sum_{i = 1}^{N_\textrm{tst}} \frac{||\boldsymbol{y}_i - \hat{\boldsymbol{y}}_i||}{||\boldsymbol{y}_i|| + ||\hat{\boldsymbol{y}}_i||}. \end{align}\] The SMAPE equals \(100\%\) if all the model predictions are twice the testing values, i.e., \(\hat{\boldsymbol{y}}_i = 2 \boldsymbol{y}_i\).

As depicted in Figure 2, alongside our model we implemented two more models and compared their accuracies over multi-step predictions. The first algorithm that we compared ours to is ‘Exact (Extended) DMD’ [78], and the second is a sparse GP model of the flow map [47]. When comparing across the three algorithms, all three models shared the same set of hyperparameters.

We also compared the VFE method to the ‘variational approach to Markov processes’ (VAMP) for hyperparameter optimization [79]. Specifically we maximized the VAMP-\(2\) score of the model.

When optimizing the VAMP score, the two DMD algorithms are virtually indistinguishable, suggesting that the Tikhonov regularization has little to no effect. This indicates that the VAMP-\(2\) objective does not capture sensor noise as effectively as maximizing the VFE of the flow map. Figure 2 shows that our GP-DMD variant outperforms Exact EDMD under noisy conditions. This aligns with our expectations given the explicit treatment of sensor noise and the additional structure imposed by the Bayesian model. However, the nonlinear growth did prompt a closer examination of our regularization strategy.

3.2 Decoupling the Lifted Observation Model↩︎

To assess the sensitivity of our algorithm to the approximation introduced in Definition 1, we separately train a regularization parameter for the lifted model using the kernelized targets \(\mathbf{K}_{ZY}\) (after estimating the optimal kernel hyperparameters and pseudo-inputs as discussed in Section 2.5). In effect, we are defining the RKHS noise model with \(\sigma_\mathcal{Y}^2\kappa_\textrm{bc}\left(\cdot, \cdot\right)\) where \(\sigma_\mathcal{Y}^2\) is not directly coupled to the variance of the sensor noise in the state space, \(\sigma_{Y}^2\). Put differently, the regularization parameter in 4 now corresponds to \(\sigma_\mathcal{Y}\), whereas in 6 it is still \(\sigma_{Y}\). This establishes a division of roles between the two regularization terms: \(\sigma_Y\) is predominantly responsible for projecting to the state space, while \(\sigma_\mathcal{Y}\) controls the accuracies of the multi-step forecasts. Consequently, the relationship in 10 does not hold exactly anymore, but our model is more flexible.

The effect of decoupling is illustrated in Figure 3 where we observe that the growth rate of the GP-TCCA model has markedly reduced, leading to more accurate long-term predictions. We therefore increase the credence in our hypothesis: rigorously accounting for the distortion introduced by the nonlinear feature map will yield a more robust model. For example, we are naturally tempted to model the statistical properties of \(\nu_i\) as non-Gaussian and input-dependent (i.e., heteroskedastic). This is apparent when we notice that \(\kappa_{\mathrm{pr}}(\boldsymbol{x}, Y)\) is typically a nonlinear many-to-one transformation of \(Y\).

3.3 Reprojections↩︎

Although the GP model attains the lowest error, the DMD variants are more computationally efficient, because multi-step forecasts follow from spectral propagation (i.e., raising the eigenvalue matrix to the desired power). In contrast, a purely GP-based rollout must repeatedly lift each intermediate state, increasing computational overhead. Consequently, for long horizons, DMD methods often employ a reprojection scheme that combines both approaches; intermittent reprojections onto the state manifold mitigates drift arising from the fact that the learned finite-dimensional subspace may not be perfectly invariant under the true Koopman operator [80].

In Figure 4 we demonstrate how reprojections based on a mechanism defined on 11 can be effectively used to decrease the computational cost while also maintaining accuracy. The idea was simply to reproject when the Euclidean norm of the diagonal entries in \({\mathbf{K}}_{\textrm{pst}}^{k}\left({\boldsymbol{x}}_{0}\right)\) exceeded some predefined tolerance. At such points, the current estimate \(\hat{\boldsymbol{x}}_k\) was treated as a noise-free measurement, and the feature vector \(\mathbf{k}_Z(\hat{\boldsymbol{x}}_k)\) was recomputed⁵.

For the depicted trajectory, on average a reprojection was performed after every \(5\) time-steps. In other words, the length of the forecasting horizons vary over the state space as the model automatically incorporates its epistemic uncertainty into the predictions. Interestingly, the regions where the vector field changes slowly is where the forecasting horizons are shorter, which suggest that errors accumulate faster in these regions.

3.4 Eigenfunction Uncertainty Quantification↩︎

The second dynamical system that we analyzed was the stochastic well problem described by the overdamped Langevin equation: \[\begin{align} {3} \begin{aligned} \mathrm{d}{x}_{t,1} &= -\nabla_{x_1} V(\boldsymbol{x}_t) \, \mathrm{d}t + \sigma_T \; \mathrm{d}W_{t,1},\\ \mathrm{d}{x}_{t,2} &= -\nabla_{x_2} V(\boldsymbol{x}_t) \, \mathrm{d}t + \sigma_T \; \mathrm{d}W_{t,2},\\ \end{aligned} \end{align}\] where \(W_{t,1}\) and \(W_{t,2}\) are two independent standard Wiener processes. In our experiments, we fixed the noise intensity to \(\sigma_T = 0.7\) and employed a sampling interval of \(\Delta t = 10.0\, [\text{s}]\).

When the potential function is \(V(\boldsymbol{x}_t) = (x_{t,1}^2 - 1)^2 + x_{t,2}^2\), the system has two minima (or “wells"), separated by a barrier, representing the two metastable states of the system. Physically, a particle will spend long periods of time near one of the two minima and only rarely jump over the barrier to the other well due to random thermal fluctuations. By slightly altering the potential function we can increase the number of wells and dominant eigenfunctions. For example, the expression of the potential function for the quadruple-well system is \(V(\boldsymbol{x}_t) = (x_{t,1}^2 - 1)^2 + (x_{t,2} - 1)^2\).

In the context of the double-well problem, the second eigenfunction, depicted in Figure 5, is the first and only non-trivial mode. This eigenfunction captures the switching between the two wells, and reveals the two-state structure and the slow timescale of the transitions [2]. In turn, the quadruple-well has three dominant non-trivial eigenfunctions, one of which we are depicting in Figure 6. Since these systems are reversible, the eigenvalues of the associated PFO and Koopman operators are exclusively real-valued [1].

Figures 5 and 6 show that the well locations align with regions of high certainty. This is because the dynamical systems spend most of their time around the wells, leading to dense sampling in these areas. While the kernel lengthscales also influence the geometry of the low-uncertainty regions, the dominant effect in these plots is the sampling distribution of the training data sets: initial conditions were chosen based on an optimized Latin hypercube design [81], [82], and succeeded by rolling out trajectories with lengths of \(20\) time-steps.

In the bottom plot of Figure 6, propagating the distribution over the eigenfunction forward in time reveals both a mixing and rescaling phenomenon: the eigenfunction collapses while the credible regions around the wells expand, ultimately spanning the entire sampling domain.

4 Conclusion↩︎

We have proposed a sparse kernel-based DMD algorithm which we formulated within a Bayesian framework that treats the embedded Perron–Frobenius operator [39] as a random variable. By unifying operator learning with the mature field of Gaussian process regression the method extends EDMD by incorporating hyperparameter optimization, sparse dictionary learning, enables propagating uncertainties in the eigenfunctions of the flow map, and provides a criteria for reprojections.

The primary insight we have gained from the ablation studies is that kernel transfer operators could benefit from a non-Gaussian and a heteroskedastic observation model; one that is closely intertwined with a Bayesian consistent notion of an operator-valued RKHS. However, we have another reason to modify our model with a more complex noise model; thus far, we have only focused on compensating for measurement noise on the target variables. Yet, DMD is known to be susceptible to adverse effects arising from noise on the inputs [83]–[86]. Incidentally, one approach to dealing with input-noise is a heteroskedastic GP model [62], [87], [88]. Alternatively, in the Koopman community, several variants of DMD have been proposed to correct for the bias induced by input-noise [89]–[92].

Earlier we posed the question of which philosophical stance is most appropriate for kernel transfer operators: a Bayesian or a frequentist perspective. We have argued that the Bayesian viewpoint offers ample opportunities to enrich the interpretation of DMD models. At the same time, we have also demonstrated that diverse perspectives can be especially valuable, as uncovering the reasons for discrepancies between the two paradigms led to the insights that solved the problem. Thus, rather than exclusively favoring one approach above the other, we advocate for employing Bayesian and frequentist methods in tandem whenever possible.

References↩︎