An MCMC Method for Uncertainty Set Generation
via Operator-Theoretic Metrics

A transfer operator is a linear map $\mathcal{L}$ acting on functions $\phi:\mathcal{M}\to\mathbb{C}$ on phase space $\mathcal{M}$ of a dynamical system $F:\mathcal{M}\to\mathcal{M}$. $\mathcal{L}$ is a linear operator over the function space $\{\phi\}$ and thus is a convenient alternative to the nonlinear state-space function $F$.

In computational fluid dynamics, statistical physics, control theory, and many other fields, the use of transfer operators in describing dynamical systems has seen a recent surge of popularity since it facilitates the usage of linear techniques (e.g. spectral methods) in analyzing nonlinear systems [1].

When $\phi$ are phase-space densities $p_{\mathcal{M}}(x)$, $\mathcal{L}$ is referred to as the Perron–Frobenius operator $\mathcal{P}$ that acts as the pushforward operator for the Markov process $p_{\mathcal{M}}^{+}(x)=\mathcal{P}p_{\mathcal{M}}(x)$. When $\phi$ are arbitrary observables $\phi:\mathcal{M}\to\mathbb{C}$, $\mathcal{L}$ is known as the Koopman (or composition) operator $\mathcal{K}$ [21]. Seminal work by Mezić [1] showed the utility of transfer operator theory in data-driven system identification and model reduction. Using Koopman eigenfunctions to obtain linear descriptions of nonlinear systems facilitated the study of dynamics via the spectrum of the Koopman operator (e.g., characterizations of chaos [22]):

where $\Phi$ and $\Lambda$ are the eigenfunctions and the eigenvalues of the Koopman operator, respectively. This led to the least-squares solutions for $\mathcal{K}$ by Schmid et al. [23], called DMD.

DMD has since been extended with dictionaries of nonlinear basis functions (extended DMD [16]), observables in reproducing kernel Hilbert spaces [24, 25], neural networks (e.g., [26]), and the incorporation of control [4]. In this paper, we build upon this prior work, exploring uncertainty quantification in the context of transfer operators which are learned from observation data.

To discuss transfer operators in the context of uncertain data, let us move to a fundamentally probabilistic setting where observations $X_{t}\in\mathcal{M}$ represent a stochastic process. While the system may not be intrinsically noisy, the probabilistic view facilitates representing uncertainty about future states. We adopt the framework of stochastic Koopman operators as introduced by Klus et al. [25] and Song et al. [27], and provide a brief background below.

Let a dynamical system $F:\mathcal{M}\to\mathcal{M}$ have invariant probability measure $\mu$, compactly supported over a measurable subspace $\Omega\subseteq\mathcal{M}$. Let the sequence $\{X_{t}|t\geq 0\}$ be a Markov process with transition density given by $p(y|x)=\Pr\{F(X)=y|X=x\}$. Given an observable $g\in L^{2}(\mathcal{M},\mu)$, the action of the stochastic Koopman operator $\mathcal{K}$ on $g$ is defined by the conditional mean:

Its right-adjoint, the Perron–Frobenius operator $\mathcal{P}$, directly maps marginal distributions as $p^{+}(x)=\mathcal{P}p(x)=\int_{\Omega}p(y|x)p(x)\mathrm{d}\mu(x)$.

Let observables $\phi$ lie in an inner product space $\mathcal{H}$ ($\mathcal{H}$ is taken to be an RKHS in [25], but for our purposes it can be the completion of a finite observable basis $|\Phi|=d$, and for trajectory visualization one may define these functions such that the preimage can be computed). Now, we may define the Gramian $\mathcal{C}_{XX}$ (i.e., a cross-covariance operator [27]):

Using the relation $\mathcal{C}_{YX}f=\mathbb{E}_{Y|X}[f(Y)|X]\mathcal{C}_{XX}$ [27], we can express the stochastic Koopman and Perron–Frobenius operators (1) in terms of the cross-covariance operators as:

that push forward observables $g\in\mathcal{H}$ and densities $p(x)$ on $\Omega$, respectively. $\epsilon$ is for ensuring invertibility.

Considering perturbations to these operators, which are estimated from cross-covariances matrices using finite trajectories, forms the population of our uncertainty set.

To compare the behavior of two dynamical systems, various metrics exist in the literature. For Markov processes, total-variation distance between density functions is commonly used, and more general classes of linear metrics over Markov chains have been proposed [28]. The standard operator norm can be used for bounded linear operators; when approximated by matrices, one can use an induced matrix norm, such as $\operatorname{Tr}(A^{T}B)$ where $A$ and $B$ are linear finite dynamical systems. However, both of these standard metrics fail to capture the iterated behavior of the respective systems. To address this, Viswanathan et al. leverage the Binet–Cauchy theorem to define a kernel between trajectories or iterated maps [29].

In [10], Ishikawa et al. generalize the Binet-Cauchy kernel to general nonlinear dynamical systems defined by their respective Perron–Frobenius operators as follows. For two dynamical systems $(D_{1},D_{2})$ specified by their initial values and maps $((X_{1,0},f_{1}),(X_{2,0},f_{2}))$,

where $\mathcal{I}:\mathbb{C}^{n}\to\mathcal{H}$ is an initial value operator embedding initial data into a Hilbert space, $\mathcal{P}^{t}$ is the $t$th iterate of the Perron–Frobenius operator on $\mathcal{H}$, and $L_{h}:\mathcal{H}\to\mathcal{H}_{\text{ob}}$ embeds states into an observable space. There is a relation between (4) and the determinant kernel in [29]; please refer to [10] for details. In the proposed sampling method, we make use of this kernel (4) to generate perturbations to Koopman operators.

Whereas the kernel proposed by Ishikawa et al. [10] is defined for Perron–Frobenius operators on RKHS, we adapt it for the Koopman operator as follows. Let the initial value and observable operators $\mathcal{I},\mathcal{L}_{h}$ be already applied, and the kernel (4) evaluated on a dynamical system $\mathcal{K}$ which acts in this (observable) space. Then we simplify (4) to:

As noted in [10], this kernel is convergent in the limit of $T$ for semi-stable $\mathcal{K}$, that is, those with spectral radius $\rho(\mathcal{K})\leq 1$ only. To use $k_{\mathcal{K}}^{m,T}$ for general Koopman operators $\mathcal{K}$, we use an exponential discounting scheme with factor $\lambda\geq 0$:

Note that the discounting factor has been adopted also in previous studies on dynamical system kernels [32, 30].

Let us show the convergence of (6) as $T\to\infty$ informally, for a finite-dimensional approximation $K$ of $\mathcal{K}$. We first observe that $\sum_{t=0}^{\infty}A^{t}$ converges if $\lim_{t\to\infty}||A^{t}||_{F}=0$ for any matrix $A$. As the product of two convergent series is convergent, it suffices to show that $\lim_{t\to\infty}||e^{-\lambda t/2}K^{t}||_{F}=0$ for all $K$. Using Gelfand’s formula,

Even if $\rho(K)>1$, the series is convergent if $\lambda>2\log\rho(K)$. Thus $k^{m,T,\lambda}$ converges for all $K$ and appropriate $\lambda$.

We propose a sampling procedure for dynamical systems models given an inner product defined over transfer operators. We define a (pseudo-)metric bounded in $[0,1]$ for operators $K_{1},K_{2}$, using a cosine similarity, as

where $\langle\cdot,\cdot\rangle_{k}$ denotes the inner product induced by a positive-definite kernel $k$. As a baseline, we may also consider the standard linear kernel $k(A,B)=\operatorname{Tr}(A^{\textsf{T}}B)$, which is not necessarily a good option for transfer operators.

Let $D\in[0,1]$ be a random variable with density $p_{D}$. Furthermore, let $K_{0}$ be a nominal transfer operator (such as ones estimated by DMD). Then we define a likelihood of any dynamical system $K$ as:

For example, we may assume $p_{D}$ is a Beta distribution with $\beta\gg\alpha$, or an exponential distribution with vanishing density past $1$. Using this we may easily construct an uncertainty set of radius $r$ as $\Delta=\{K\sim\mathcal{L}(K|K_{0})\ |\ d_{k}(K,K_{0})\leq r\}$.

Sampling $K$ using $\mathcal{L}$ can be done in a number of ways. The conceptually simplest algorithm is rejection sampling: for any uniformly perturbed $K$, accept if $u<c\mathcal{L}(K|K_{0})$ for $u\sim[0,1]$ and convergence parameter $c$. Unfortunately, even generating the initial uniform perturbations may be computationally intractable due to the high dimensionality of the samples. This results in low acceptance rates for rejection sampling as well as random-walk MCMC methods such as Metropolis–Hastings. We turn our focus to gradient-based MCMC methods which are able to generate distant proposals and achieve dimensionality-independent acceptance rates.

In his seminal work [31], Neal introduced Hamiltonian Monte Carlo (HMC), which uses the gradient of the likelihood to simulate stochastic Hamiltonian dynamics whose stationary distribution is the posterior (7). It is well suited to our case since many kernels over dynamical systems are basically differentiable.

We adapt HMC for transfer operator sampling by introducing an auxiliary momentum variable $R$ which is of the same dimension as $K$, whose relation to $K$ is given via Hamiltonian dynamics. Let us define the potential and the Hamiltonian of an operator $K$ as:

Noting that the potential $U$ must be defined (and continuous) everywhere in order to achieve the correct stationary distribution, $d_{k}$ must be similarly well-behaved for all $K$. Using the discounted kernel (6), we generate samples $\{K\}$ of Koopman operators about a nominal $K_{0}$ via Hamiltonian dynamics in the potential defined by (8) using the leapfrog integrator (we refer the reader to [31] for details).

The mixing time of HMC is highly sensitive to the choice of discretization parameters (in particular, n_leapfrog, step_size). In practice, Neal recommends $\epsilon\sim O(d^{1/4})$ [31], and there also exist adaptive step-determining algorithms such as the No-U-Turn Sampler [33]; however, we find a tradeoff between computational efficiency (samples/step) and sufficient exploration of the model space when using HMC for transfer operators. To compensate, we use a pre-run of HMC in a zero-potential to generate a uniform prior of samples, which are then used as initial conditions for parallel HMC sampling from the distribution (7). We find that this accelerates the mixing time and wall-clock time for sampling significantly.

In the interest of producing meaningful samples, one may wish to constrain the space of models using prior knowledge of the underlying system. Constraints can be readily expressed as boundary conditions on HMC without changing its stationary solution or reversibility (under some assumptions) – known as Reflective HMC [34]. As an example, suppose that we have knowledge that the system is structurally stable under any realistic perturbation; then, we can (roughly) encode this as a constraint on the spectral radius $\rho(K_{0})-\epsilon\leq\rho(K)\leq\rho(K_{0})+\epsilon$.

Extending the reflective HMC procedure from [34], we describe a leapfrog integrator (Algorithm 1) which ensures $f(K)\in[a,b]$ for any differentiable, scalar-valued $f$. In numerical experiments, we use $f(K)=\rho(K)$ and $[a,b]=[\rho(K_{0})-.01,\,\rho(K_{0})+.01]$.

We will give two formulations of HMC (8) for finite arguments, one explicitly over transfer operators, and one implicitly over observed trajectories. Either may be used.

We first consider linear systems of the form $\dot{x}=Ax$ via discretization as $A_{d}=e^{A\Delta t}$. We use simple 2x2 systems in order to clearly characterize the dynamics in a trace-determinant plot. Using the discounted kernel (10), we are able to generate perturbations of both source- and saddle-types in addition to the semistable regimes. We use spread parameter $\beta=5$, HMC step $\epsilon=10^{-4}$, HMC leapfrog $L=100$, and generate $N=1000$ samples with $k=50$ initial conditions for every test shown. We compare the following two kernels:

termed the “trace kernel” and “Koopman kernel”, respectively.

The strength of the proposed method, using $k^{m,T,\lambda}_{K}$, can be seen when the nominal system is within a region of structural instability (see the two center systems, Figure 1). The trace kernel perturbations venture easily into spiral sink or spiral sources, which are distant in dynamics terms but very close in absolute norm. In these cases, the proposed method using the Koopman kernel retains a much tighter spread in the trace-determinant plane. Furthermore, it can be seen from the posterior distributions that the proposed method is able to explore distant dynamics while staying bounded within structurally similar regions.

Next, we consider perturbations of a nonlinear system. We use the unforced Duffing equation:

whose basins of attraction are like in Figure 2.

We use simulated trajectories of length $t=400$ seconds with $8000$ samples per trajectory across $144$ initial conditions in the range $[-2,-2],[2,2]$. Using $15$ polynomial observables with maximum degree $5$, we compute the Koopman operator $K$ as $K=YX^{\textsf{T}}(XX^{\textsf{T}})^{\dagger}$. For all experiments we use spread $\beta=5$, HMC step $\epsilon=5\times 10^{-5}$, HMC leapfrog $L=200$, $N=2000$ samples, and $k=200$ initial conditions. In the shown perturbations, trajectories in the left and right basins of attraction are highlighted in red and blue, respectively (Figure 3).

We immediately observe some qualitative differences between the two perturbation sets. First, it is apparent that perturbations via the Koopman kernel preserve attractor structure in most samples, versus almost none in the baseline setting. In the case of the Duffing oscillator, this is a defining feature, and such preservations are important consideration for any robust prediction or control procedure over dynamics models. Second, a large proportion of samples produced by the baseline method are diverging; these would need to be manually filtered out if used in a robust optimization setting. We observe in experiments that this can be mitigated by restricting the norm of perturbations (i.e., increasing $\beta$), but this comes at the cost of decreased exploration of dynamics space and changes the robustness of an RO solution using the perturbation set (moreover, manual filtering changes the posterior, altering the RO problem). We also find that the spectral radius constraint alleviates many of these concerns with the baseline method, however, non-convex reflection is not a trivial procedure to implement in HMC and has not been a typically used method in generating uncertainty sets.

Finally, while attractors are mostly preserved in our method, the attractor basins seem to undergo some geometry warping. This suggests an interpretation of our perturbation method as warping the underlying potential wells, which may have meaningful physical interpretations.