2.1. Finite-dimensional dynamical system

We start with a generic continuous time dynamical system in $d$-dimensions given by

(2.1)\begin{equation} \dot{\boldsymbol{s}} = \boldsymbol{U}(\boldsymbol{s}), \end{equation}

where $\boldsymbol{s}(t): \mathbb {R} \rightarrow \mathbb {R}^d$ is the state of the system and $\boldsymbol{U}: \mathbb {R}^d \rightarrow \mathbb {R}^d$ is the evolution equation. Equation (2.1) provides a succinct rule for determining the evolution of a dynamical system; however, uncertain initial conditions stymie future predictions in the presence of chaos (Lorenz Reference Lorenz1963). It is, therefore, more natural to study the statistical evolution of probability densities as in Hopf (Reference Hopf1952). Thus, we focus not on the $d$-dimensional ordinary differential equation given by (2.1) but rather the $d$-dimensional partial differential equation that governs the evolution of probability densities in state space.

We denote fixed vector in state space by $\boldsymbol{\mathscr {s}} \in \mathbb {R}^d$, the components of the state $\boldsymbol{\mathscr {s}}$ by $\boldsymbol{\mathscr {s}} = (\mathscr {s}_1, \mathscr {s}_2,\ldots,\mathscr {s}_d)$ and the components of the evolution rule by $\boldsymbol{U} = (U_1, U_2,\ldots,U_d)$. The evolution equation for the statistics of (2.1), as characterized by a probability density function,

(2.2)\begin{equation} \mathcal{P} = \mathcal{P}(\mathscr{s}_1, \mathscr{s}_2,\ldots, \mathscr{s}_d, t) = \mathcal{P}(\boldsymbol{\mathscr{s}}, t), \end{equation}

is given by the Liouville equation

(2.3)\begin{equation} \partial_t \mathcal{P} + \sum_{i=1}^d \frac{\partial}{\partial \mathscr{s}_i} \left( U_i(\boldsymbol{\mathscr{s}}) \mathcal{P} \right) = 0. \end{equation}

The above equation is a statement of probability conservation. It is precisely analogous to the mass conservation equation from the compressible Navier–Stokes equations. However, the ‘mass’ is being interpreted as a probability density. The distribution, $\mathcal {P}$, is guided by the flow dynamics $\boldsymbol{U}$ to likely regions of state space.

The presence of stochastic white noise, $\boldsymbol {\xi }$ in a dynamical system,

(2.4)\begin{equation} \dot{\boldsymbol{s}} = \boldsymbol{U}(\boldsymbol{s}) + \boldsymbol{\mathcal{D}^{1/2}} \boldsymbol{\xi} \end{equation}

where the ensemble average of $\boldsymbol{\xi }$ satisfies $\langle \xi _i(t) \xi _j(t') \rangle = \delta _{ij} \delta (t - t')$ and $\boldsymbol{\mathcal {D}^{1/2}}$ is the matrix square root of the covariance matrix $\boldsymbol{\mathcal {D}}$, e.g. $\boldsymbol{\mathcal {D}} = \boldsymbol{\mathcal {D}^{1/2}} \boldsymbol{\mathcal {D}^{1/2}}$, modifies (2.3) through the inclusion of a diffusion term

(2.5)\begin{equation} \partial_t \mathcal{P} + \sum_{i=1}^d \frac{\partial}{\partial \mathscr{s}_i} \left( U_i(\boldsymbol{\mathscr{s}}) \mathcal{P} + \frac{1}{2} \sum_{i=j}^d \mathcal{D}_{ij} \frac{\partial}{\partial \mathscr{s}_j} \mathcal{P} \right) = 0, \end{equation}

i.e. the Fokker–Planck equation. The presence of the diffusion tensor $\boldsymbol{\mathcal {D}}$ regularizes (2.3).

2.2. Infinite-dimensional dynamical system

When the underlying dynamical system is a partial differential equation, we assume a suitably well-defined discretization exists to reduce it to a formally $\mathbb {R}^d$-dimensional dynamical system. We contend ourselves to the study of the statistics of the $\mathbb {R}^d$ approximation. One hopes that different discretizations lead to similar statistical statements of the underlying partial differential equation; thus, it is worth introducing notation for the analogous Liouville equation for a partial differential equation, as was done by Hopf (Reference Hopf1952). The calculations and notation that follow are formal, and there are mathematical difficulties in assigning them rigorous meaning; however, we find value in the approach, as it often allows for expedient calculations as was done in Souza, Lutz & Flierl (Reference Souza, Lutz and Flierl2023b) and Giorgini et al. (Reference Giorgini, Deck, Bischoff and Souza2024). The overall (heuristic) recipe for the transition to function spaces is to replace sums with integrals, derivatives with variational derivatives, discrete indices with continuous indices and Kronecker deltas with Dirac deltas.

The $d$-dimensional vector from before now becomes a vector in function space whose components are labelled by a continuous index, $\boldsymbol{x}$, a position in a domain $\varOmega$, and discrete index $j$, the index for the field of interest. Thus, the component choice $s_i(t)$ for a fixed index $i$ is analogous to $s_{(\boldsymbol{x}, j)}(t)$ for a fixed position $\boldsymbol{x}$ and field index $j \in \{1,\ldots, d_s \}$, e.g. $d_s = 3$ for the three velocity components of the incompressible Navier–Stokes. In the discrete case, the single index $i$ loops over all velocity components and all simulation grid points.

Specifically, we consider a partial differential equation for a state $\boldsymbol{s}$ defined over a domain $\varOmega$, with suitable boundary conditions,

(2.6)\begin{equation} \partial_t \boldsymbol{s} = \boldsymbol{\mathcal{U}}[\boldsymbol{s}], \end{equation}

where the operator $\boldsymbol{\mathcal {U}}: \mathcal {X} \rightarrow \mathcal {X}$ characterizes the evolution of system and $\mathcal {X}$ is a function space. The component of $\boldsymbol{\mathcal {U}}$ at position $\boldsymbol{x}$ and field index $j$ is denoted by $\mathcal {U}_{(\boldsymbol{x}, j)}$.

The analogous evolution for the probability density functional,

(2.7)\begin{equation} \mathcal{P} =\mathcal{P}[\mathscr{s}_{(\boldsymbol{x},1)}, \mathscr{s}_{(\boldsymbol{x},2)},\ldots, \mathscr{s}_{(\boldsymbol{x},d_s)}, t] = \mathcal{P}[\boldsymbol{\mathscr{s}}, t], \end{equation}

is denoted by

(2.8)\begin{equation} \partial_t \mathcal{P} + \sum_{j=1}^{d_s} \int_{\varOmega} {\rm d}\kern0.07em\boldsymbol{x}\, \frac{\delta }{\delta \mathscr{s}_{(\boldsymbol{x},j)} } \left( \mathcal{U}_{(\boldsymbol{x},j)}[\boldsymbol{\mathscr{s}}] \mathcal{P} \right) = 0. \end{equation}

The sum in (2.3) is replaced by an integral over position indices and a sum over field indices in (2.8). Furthermore, the partial derivatives are replaced by variational derivatives. The variational derivative is being used in the physicist's sense, that is to say,

(2.9)\begin{equation} \frac{\delta \mathscr{s}_{(\boldsymbol{x }, i)} }{\delta \mathscr{s}_{ (\boldsymbol{y}, j)} } = \delta (\boldsymbol{x} - \boldsymbol{y}) \delta_{ij} \Leftrightarrow \frac{\partial \mathscr{s}_{i'}}{\partial \mathscr{s}_{j'}} = \delta_{i' j'}, \end{equation}

where $\delta (\boldsymbol{x} - \boldsymbol{y})$ is the Dirac delta function and $\delta _{ij}$ is the Kronecker delta function, $\boldsymbol{x}, \boldsymbol{y} \in \varOmega$, $i,j \in \{1,\ldots, d_s \}$, and $i',j' \in \{1,\ldots, d\}$. In the typical physics notation, it is common to drop the dependence on the position $\boldsymbol{x}$ and explicitly write out the field variable in terms of its components (as opposed to the indexing that we do here), e.g.

(2.10)\begin{equation} \frac{\delta }{\delta \mathscr{s}_{ (\boldsymbol{y },1)} } \sum_{i=1}^{d_s} \int_{\varOmega} {\rm d}\kern0.07em\boldsymbol{x} (\mathscr{s}_{(\boldsymbol{x}, i)})^2 = 2 \mathscr{s}_{(\boldsymbol{y},1)} \Rightarrow \frac{\delta }{\delta u} \int_{\varOmega} (u^2 + v^2 + w^2)= 2 u \end{equation}

for the prognostic variables $u,v,w$ of the incompressible Navier–Stokes equations. See § 35 of Zinn-Justin (Reference Zinn-Justin2021) for examples of the Fokker–Planck equation in the field-theoretic context.

To derive (2.8), we suppose that (2.1) is a discretization of (2.6). Starting from (2.3), first introduce a control volume at index i as $\Delta \boldsymbol{x}_i$ to rewrite the equation as

(2.11)\begin{equation} \partial_t \mathcal{P} + \sum_{i=1}^d \Delta\kern0.07em \boldsymbol{x}_i \frac{1}{\Delta\kern0.07em \boldsymbol{x}_i}\frac{\partial}{\partial \mathscr{s}_i} \left( U_i(\boldsymbol{\mathscr{s}}) \mathcal{P} \right) = 0. \end{equation}

In the ‘limit’ as the grid is refined and $|\Delta \boldsymbol{x}_i| \rightarrow 0$, we have

(2.12a,b)\begin{equation} \sum_{i=1}^d \Delta\kern0.07em \boldsymbol{x}_i \frac{1}{\Delta\kern0.07em \boldsymbol{x}_i}\frac{\partial}{\partial \mathscr{s}_i} \rightarrow \sum_{j=1}^{d_s} \int_{\varOmega} {\rm d}\kern0.07em\boldsymbol{x}\,\frac{\delta }{\delta \mathscr{s}_{(\boldsymbol{x},j)}}\quad \text{and}\quad U_i \rightarrow \mathcal{U}_{(\boldsymbol{x},j)}. \end{equation}

A stochastic partial differential equation is given by

(2.13)\begin{equation} \partial_t \boldsymbol{s} = \boldsymbol{\mathcal{U}}[\boldsymbol{s}] + \boldsymbol{\mathcal{D}^{1/2}}[\boldsymbol{\xi}], \end{equation}

where $\boldsymbol{\xi }$ is space–time (and state) noise with covariance

(2.14)\begin{equation} \langle \xi_{(\boldsymbol{x}, i)}(t) \xi_{(\boldsymbol{y}, j)}(t') \rangle = \delta_{ij} \delta (\boldsymbol{x} - \boldsymbol{y}) \delta (t - t'). \end{equation}

The noise $\boldsymbol{\mathcal {D}^{1/2}}[\boldsymbol{\xi }]$ is interpreted as a Gaussian process with a spatial and field covariance given by properties of $\boldsymbol{\mathcal {D}}$. Furthermore, the action of the symmetric positive definite linear operator $\boldsymbol{\mathcal {D}}$, is expressed as

(2.15)\begin{equation} \mathcal{D}_{(\boldsymbol{x}, i)}[\boldsymbol{\xi}] = \sum_{j=1}^{d_s} \int_{\varOmega} {\rm d} \boldsymbol{y}\, \mathcal{K}_{ij}(\boldsymbol{x}, \boldsymbol{y}) \xi_{(\boldsymbol{y},j)}, \end{equation}

where $\mathcal {K}_{ij}$ is the kernel of the integral operator. Formally, (2.13) has the corresponding Fokker–Planck equation

(2.16)\begin{equation} \partial_t \mathcal{P} + \sum_{j=1}^{d_s} \int_{\varOmega} {\rm d}\kern0.07em \boldsymbol{x} \, \frac{\delta }{\delta \mathscr{s}_{(\boldsymbol{x},j)} } \left( \mathcal{U}_{(\boldsymbol{x},j)}[\boldsymbol{\mathscr{s}}] \mathcal{P} + \frac{1}{2}\sum_{k=1}^{d_s} \int_{\varOmega} {\rm d} \boldsymbol{y} \, \mathcal{K}_{jk}(\boldsymbol{x}, \boldsymbol{y}) \frac{\delta \mathcal{P} }{\delta \mathscr{s}_{(\boldsymbol{y},k)} } \right) = 0. \end{equation}

Making sense of (2.13) is an area of active research (see Hairer (Reference Hairer2014) and Corwin & Shen (Reference Corwin and Shen2020)), where one must confront defining probability distributions over function spaces. Defining integrals over function spaces has met with considerable difficulties, although progress has been made (see Daniell (Reference Daniell1919), DeWitt (Reference DeWitt1972) and Albeverio & Mazzucchi (Reference Albeverio and Mazzucchi2016)).

With the formalism now set, the focus of this work is on methods for discretizing (2.3) and (2.8) on subsets of state space that are typically thought of as chaotic or turbulent given only trajectory information from (2.1) and (2.6), respectively. The data-driven methods developed herein apply without change to the stochastic analogues.

2.3. Finite-volume discretization

To focus our discussion, we use the finite-dimensional and deterministic setting. We first assume that the underlying dynamics are on a chaotic attractor associated with a compact subset of state space $\mathcal {M} \subset \mathbb {R}^d$. We further assume that the dynamical system is robust to noise in order to regularize the attractor and probability densities as in Young (Reference Young2002) and Cowieson & Young (Reference Cowieson and Young2005). Statements about deterministic chaos in this manuscript implicitly use a noise regularization where the zero noise limit is always the last limit taken in any computation, similar to Giannakis (Reference Giannakis2019). Whether or not such a limit generally agrees with a Sinai–Ruelle–Bowen measure (see for example Young (Reference Young2002), Blank (Reference Blank2017) and Araujo (Reference Araujo2023)) is unknown. Furthermore, in the presence of correlated or state-space-dependent noise, one must be careful how a zero-noise limit is taken.

We introduce a partition of $\mathcal {M}$ into $N$ cells which we denote by $\mathcal {M}_n$ for $n = 1,\ldots, N$. The coarse-grained discretization variables $P_n$ are

(2.17)\begin{equation} \int_{\mathcal{M}_n} {\rm d} \boldsymbol{\mathscr{s}} \,\mathcal{P} = \int_{\mathcal{M}_n} \mathcal{P} = P_n \end{equation}

as is common in finite-volume methods. When unambiguous, we drop the infinitesimal state space volume ${\rm d}\boldsymbol{\mathscr {s}}$. Here $P_n(t)$ is the probability in time of being found in the subset of state space $\mathcal {M}_n$ at time $t$. Integrating (2.3) with respect to the cells yields

(2.18)\begin{equation} \frac{{\rm d}}{{\rm d} t} P_n = \int_{\mathcal{M}_n} \left[ \sum_{i=1}^d \frac{\partial}{\partial \mathscr{s}_i} \left( U_i(\boldsymbol{\mathscr{s}}) \mathcal{P} \right) \right] = \int_{\partial \mathcal{M}_n} \boldsymbol{U} \boldsymbol{\cdot} \hat{\boldsymbol{n}} \mathcal{P}, \end{equation}

where $\partial \mathcal {M}_n$ is the boundary of the cell and $\hat {n}$ is a normal vector. The art of finite-volume methods comes from expressing the right-hand side of (2.18) in terms of the coarse-grained variables $P_n$ through a suitable choice of numerical flux.

We go about calculating the numerical flux in a roundabout way. We list some desiderata for a numerical discretization.

1. (i) The discrete equation is expressed in terms of the instantaneous coarse-grained variables $P_n$.

2. (ii) The discrete equation is linear, in analogy to the infinite-dimensional one.

3. (iii) The equation must conserve probability.

4. (iv) Probability must be positive at all times.

The first two requirements state

(2.19)\begin{equation} \int_{\partial \mathcal{M}_n} \boldsymbol{U} \boldsymbol{\cdot} \hat{\boldsymbol{n}} \mathcal{P} \approx \sum_{m} Q_{nm} P_m \end{equation}

for some matrix $Q$. Thus, we want an equation of the form

(2.20)\begin{equation} \frac{{\rm d}}{{\rm d}t} \hat{P}_n = \sum_m Q_{nm} \hat{P}_m. \end{equation}

We introduced a ‘hat’ to distinguish the numerical approximation, $\hat {P}_n$, with the exact solution $P_n$. The third requirement states

(2.21)\begin{equation} \sum_n \hat{P}_n = 1 \Rightarrow \frac{{\rm d}}{{\rm d}t} \sum_n \hat{P}_n = 0 = \sum_{mn} Q_{nm} \hat{P}_m \end{equation}

for each $\hat {P}_m$, thus

(2.22)\begin{equation} \sum_{n} Q_{nm} = 0 \end{equation}

for each $m$, i.e. the columns of the matrix must add to zero. Moreover, the last requirement states that the off-diagonal terms of $Q_{nm}$ must all be positive. To see this last point, we do a proof by contradiction. Suppose there is a negative off-diagonal entry, without loss of generality, component $Q_{21}$. Then if at time zero our probability vector starts at $\hat {P}_1(0) = 1$ and $\hat {P}_n(0) = 0$ for $n > 1$, an infinitesimal time step ${\rm d}t$ later we have

(2.23)\begin{equation} \hat{P}_2({\rm d}t ) = {\rm d}t \sum_{m} Q_{2m} \hat{P}_m(0) = Q_{21} < 0, \end{equation}

a contradiction since probabilities must remain positive at all times. Of all the requirements, the combination of the second and the fourth ones restricts us to a lower-order discretization since it is not possible to have a higher-order discretization that is both positivity preserving and linear (Zhang & Shu Reference Zhang and Shu2011).

Motivation for (2.20) yielding an approximation of the underlying partial differential equation comes directly from the ability of finite-volume methods to converge (under suitable assumptions) to an underlying linear partial differential equation. An example of a matrix, $Q$, with such a structure, can be seen in the appendix of Hagan, Doering & Levermore (Reference Hagan, Doering and Levermore1989), where an approximation to an Ornstein–Uhlenbeck process is constructed. See appendices C.2 and C.3 of Souza et al. (Reference Souza, Lutz and Flierl2023b) for examples on analytically constructing coarse-grained operators of stochastic systems with these properties and Appendix D for an example with the simple harmonic oscillator.

These four requirements, taken together, are enough to identify the matrix $Q$ as the generator of a continuous-time Markov process with finite state space. This observation forms the backbone of the data-driven approach towards discretizing (2.3) and (2.8). The diagonal entries of the matrix are related to the average amount of time spent in a cell, and the off-diagonal entries within a column are proportional to the probabilities of entering a given cell upon exit of a cell. The implication is that we construct the numerical fluxes on a boundary through Monte Carlo integration of the equations of motion.

Intuitively, as a state trajectory enters passes through a cell, it becomes associated with the cell. The time a trajectory spends within a cell is called the holding time from whence it will eventually exit to some other cell in state space. A sufficiently long integration of the equations of motion constructs the holding time distributions of a cell and exit probabilities to different cells in state space. Furthermore, to perform calculations, we associate each region of state space with a ‘cell centre’, which, in this paper, we call a Markov state. The Markov state will serve as a centre of a delta function approximation to the steady distribution in that region of state space. With the transition matrix $Q$ and the Markov states associated with a partition, we can perform calculations of moments, steady-state distributions and autocorrelations of any variable of interest.

2.4. Time versus ensemble calculations

At this point, we have discussed the equation for statistics of a dynamical system, the notation for the infinite-dimensional case and how to associate a continuous time Markov process with a finite-volume discretization of the Liouville equation. We now discuss how to perform statistical calculations from the discretization and how we will confirm that the discretization captures statistics of the underlying Liouville equation. In short, we compare temporal averages with ensemble averages and analogous calculations for autocorrelations.

We must introduce additional notation. As stated before, we assume that a chaotic attractor exists and that there is an appropriately regularized invariant measure, which we denote by $\mathcal {I}(\boldsymbol{\mathscr {s}})$, associated with it (Cowieson & Young Reference Cowieson and Young2005). We further assume that the zero noise limit is the last limit taken in all computations in order to relate computations to deterministic systems. With regards to a turbulent flow, we are assuming that a statistically steady-state exists and denote its probability density (formally) by $\mathcal {I}(\boldsymbol{\mathscr {s}})$. The conditional invariant measure with respect to a cell $\mathcal {M}_n$ is denoted by $\mathcal {I}_n(\boldsymbol{\mathscr {s}} )$ and the probability of a state being found in a cell $\mathcal {M}_n$ is $\mathbb {P}(\mathcal {M}_n) = \int _{\mathcal {M}_n} \mathcal {I}$ so that the invariant measure is decomposed as

(2.24)\begin{equation} \mathcal{I}(\boldsymbol{\mathscr{s}}) = \sum_n \mathcal{I}_n(\boldsymbol{\mathscr{s}} ) \mathbb{P}(\mathcal{M}_n). \end{equation}

In addition, we introduce notation for the transfer operator $\mathscr {T}^\tau$, which is defined through the relation

(2.25)\begin{equation} \mathcal{P}(\boldsymbol{\mathscr{s}}, t + \tau ) = \mathscr{T}^\tau \mathcal{P}(\boldsymbol{\mathscr{s}}, t ), \end{equation}

where $\mathcal {P}$ is a solution to (2.3). Thus, the transfer operator is an instruction to evolve the density, $\mathcal {P}$, via (2.3) to a time $\tau$ in the future. Furthermore,

(2.26)\begin{equation} \lim_{\tau \rightarrow \infty} \mathscr{T}^\tau \mathcal{P}(\boldsymbol{\mathscr{s}}, t ) = \mathcal{I}( \boldsymbol{\mathscr{s}}) \end{equation}

for arbitrary densities $\mathcal {P}(\boldsymbol{\mathscr {s}}, t )$, including $\delta ( \boldsymbol{\mathscr {s}} - \boldsymbol{\mathscr {s}}')$ for an initial state $\boldsymbol{\mathscr {s}}' \in \mathcal {M}$ from our assumption of ergodicity. (We are being sloppy with limits here, but this should be understood where the limit to a delta function density is the last limit taken.)

For an observable $g: \mathcal {M} \rightarrow \mathbb {R}$, we calculate long time averages

(2.27)\begin{equation} \langle g \rangle_T = \lim_{T \rightarrow \infty} \frac{1}{T} \int_0^T g( \boldsymbol{s} (t)) \,{\rm d}t \end{equation}

and compare with ensemble averages

(2.28)\begin{equation} \langle g \rangle_E = \int_{\mathcal{M}} g(\boldsymbol{\mathscr{s}} ) \mathcal{I}(\boldsymbol{\mathscr{s}}). \end{equation}

Furthermore, we compare time-correlated observables. The time-series calculation is

(2.29)\begin{equation} R_T(g, \tau) = \lim_{T \rightarrow \infty} \frac{1}{T} \int_0^T g(\boldsymbol{s}(t+\tau))g(\boldsymbol{s}(t)) \,{\rm d}t, \end{equation}

whence we obtain the autocovariance, $C_T$, and autocorrelation, $\tilde {C}_T$,

(2.30a,b)\begin{equation} C_T(g, \tau) \equiv R_T(g, \tau) - \langle g \rangle_T^2\quad\text{and}\quad \tilde{C}_T(g, \tau) \equiv C_T(g, \tau) / C_T(g, 0). \end{equation}

The ensemble average version requires more explanation. We correlate a variable $g(\boldsymbol {s}(t))$ with $g(\boldsymbol{s}(t+\tau ))$, which involves the joint distribution of two variables. We first review a fact about random variables $X, Y$ with joint density $\rho (x,y)$, conditional density $\rho (x |y)$ and marginal density $\rho _y(y)$. The correlation of two observables is calculated as

(2.31)\begin{align} \langle g(X) g(Y) \rangle &= \iint {{\rm d}\kern0.07em x}\, {{\rm d}y}\, g(x) g(y) \rho(x,y) = \iint {{\rm d}\kern0.07em x} \,{{\rm d}y}\, g(x) g(y) \rho(x|y)\rho_y(y) \end{align}

(2.32)\begin{align} &= \int {{\rm d}y}\, g(y) \rho_y(y) \left[ \int {{\rm d}\kern0.07em x}\, g(x) \rho(x|y) \right]. \end{align}

To translate the above calculation to the present case, we consider $\rho _y$ as the invariant measure, $\mathcal {I}$. The conditional distribution $\rho (x|y)$ is thought of as the probability density at a time $\tau$ in the future, given that we know that it is initially at state $\boldsymbol{\mathscr {s}}$ at $\tau = 0$.

Thus, in our present case, $\rho (x|y)$ becomes $\mathscr {T}^\tau \delta ( \boldsymbol{\mathscr {s}} - \boldsymbol{\mathscr {s}}')$ where the $\delta$ function density is a statement of the exact knowledge of the state at time $\tau = 0$. In total, the ensemble time-autocorrelation is calculated as

(2.33)\begin{equation} R_E(g, \tau)= \int_{\mathcal{M}} {\rm d}\boldsymbol{\mathscr{s}}' \,g(\boldsymbol{\mathscr{s}}' ) \mathcal{I}( \boldsymbol{\mathscr{s}}' )\left[ \int_{\mathcal{M}} {\rm d}\boldsymbol{\mathscr{s}}\, g(\boldsymbol{\mathscr{s}}) \mathscr{T}^\tau \delta ( \boldsymbol{\mathscr{s}} - \boldsymbol{\mathscr{s}}') \right]. \end{equation}

The autocovariance, $C_E$, and autocorrelation, $\tilde {C}_E$, are

(2.34a,b)\begin{equation} C_E(g, \tau) \equiv R_E(g, \tau) - \langle g \rangle_E^2 \quad\text{and}\quad \tilde{C}_E(g, \tau) \equiv C_E(g, \tau) / C_E(g, 0). \end{equation}

These calculations summarize the exact relations we wish to compare. However, first, we will approximate the temporal averages via long-time finite trajectories and the ensemble averages via the finite-volume discretization from § 2.3.

2.5. Approximations to time versus ensemble calculations

The prior section represents the mathematical ideal with which we would like to perform calculations; however, given that we use a data-driven construction, we are faced with performing calculations in finite-dimensional spaces and over finite-dimensional time.

Given the time series of a state at evenly spaced times at times $t_n$ for $n = 1$ to $N_t$ with time spacing $\Delta t$, we approximate the mean and long time averages of an observable $g$ as

(2.35)$$\begin{gather} \langle g \rangle_T \approx \frac{1}{N_t} \sum_{n = 1}^{N_t} g(\boldsymbol{s}(t_n)), \end{gather}$$

(2.36)$$\begin{gather}R_T(g, \tau)\approx \frac{1}{N_t'} \sum_{n = 1}^{N_t'} g(\boldsymbol{s}(t_n + \text{round}(\tau / \Delta t) \Delta t)) g(\boldsymbol{s}(t_n)), \end{gather}$$

where the $\text {round}$ function computes the closest integer and $N_t' = N_t - \text {round}(\tau / \Delta t)$.

We use the construction from § 2.3 to calculate ensemble averages. Recall that in the end, we had approximated the generator of the process with a matrix $Q$, which described the evolution of probabilities associated with cells of state space. In addition, we select a state $\boldsymbol{\sigma }^{[n]}$ associated with a cell $\mathcal {M}_n$ as the ‘cell centre’ to perform calculations. We use superscripts to denote different states since subscripts are reserved for the evaluation of the component of a state. Furthermore, we do not require the Markov state $\boldsymbol{\sigma }^{[n]}$ to be a member of the cell $\mathcal {M}_n$. For example, we could choose the $\boldsymbol{\sigma }^{[n]}$ as fixed points of the dynamical system or a few points along a periodic orbit within the chaotic attractor $\mathcal {M}$.

The ensemble average of an observable is calculated by making use of the decomposition of the invariant measure (2.24), but then approximating

(2.37)\begin{equation} \mathcal{I}_n(\boldsymbol{\mathscr{s}} ) \approx \delta(\boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[n]} ) \end{equation}

which is a simple but crude approximation. Thus, the ensemble averages are calculated as

(2.38)\begin{align} \langle g \rangle_E &= \int_{\mathcal{M}} g( \boldsymbol{\mathscr{s}} ) \left[ \sum_n \mathcal{I}_n( \boldsymbol{\mathscr{s}} ) \mathbb{P}(\mathcal{M}_n) \right] = \sum_n \left[ \int_{\mathcal{M}} g( \boldsymbol{\mathscr{s}} ) \mathcal{I}_n(\boldsymbol{\mathscr{s}} ) \right] \mathbb{P}(\mathcal{M}_n), \end{align}

(2.39)\begin{align} &\approx \sum_n \left[ \int_{\mathcal{M}} g( \boldsymbol{\mathscr{s}} ) \delta \left( \boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[n]} \right)\right] \mathbb{P}(\mathcal{M}_n) = \sum_n g( \boldsymbol{\sigma}^{[n]} ) \mathbb{P}(\mathcal{M}_n). \end{align}

For the ensemble average version of time autocorrelations, we must, in addition to approximating the invariant measure, approximate the transfer operator acting delta function density of the state, $\mathscr {T}^\tau \delta (\boldsymbol{\mathscr {s}} - \boldsymbol{\mathscr {s}}')$. We calculate

(2.40)\begin{align} R_E(g, \tau)&= \int_{\mathcal{M} } \int_{\mathcal{M} } {\rm d}\boldsymbol{\mathscr{s}}'\,{\rm d}\boldsymbol{\mathscr{s}} \,g(\boldsymbol{\mathscr{s}}' ) \mathcal{I}( \boldsymbol{\mathscr{s}}' ) g(\boldsymbol{\mathscr{s}}) \mathscr{T}^\tau \delta ( \boldsymbol{\mathscr{s}} - \boldsymbol{\mathscr{s}}') \end{align}

(2.41)\begin{align} &\approx \sum_{n} \int_{\mathcal{M} } \int_{\mathcal{M} } {\rm d}\boldsymbol{\mathscr{s}}'\,{\rm d}\boldsymbol{\mathscr{s}} \,g(\boldsymbol{\mathscr{s}}' ) \delta(\boldsymbol{\mathscr{s}}' - \boldsymbol{\sigma}^{[n]} ) \mathbb{P}(\mathcal{M}_n) g(\boldsymbol{\mathscr{s}}) \mathscr{T}^\tau \delta ( \boldsymbol{\mathscr{s}} - \boldsymbol{\mathscr{s}}') \end{align}

(2.42)\begin{align} &= \sum_{n} g(\boldsymbol{\sigma}^{[n]} ) \mathbb{P}(\mathcal{M}_n) \int_{\mathcal{M} } {\rm d}\boldsymbol{\mathscr{s}} \, g(\boldsymbol{\mathscr{s}}) \mathscr{T}^\tau \delta ( \boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[n]}) \end{align}

then additionally approximate

(2.43)\begin{equation} \mathscr{T}^\tau \delta(\boldsymbol{\mathscr{s}}- \boldsymbol{\sigma}^{[n]}) \approx \sum_{m=1}^{N} \delta(\boldsymbol{\sigma}^{[m]} - \boldsymbol{\mathscr{s}}) [\exp(Q \tau)]_{mn} . \end{equation}

The matrix exponential $\exp (Q \tau )$ is the discrete Perron–Frobenius/transfer operator. We also directly consider the Perron–Frobenius operator at time scale $\tau$ and denote this construction by $\mathcal {F}^{[\tau ]}$. The matrix is a (column) stochastic matrix whose entries sum to one. The intuition behind the approximation is to treat the forward evolution of the transfer operator for delta distribution centred at state $\boldsymbol{\sigma }^{[n]}$ as a weighted sum of delta functions centred at state $\boldsymbol{\sigma }^{[m]}$. The weights are given by the probability of being in cell $\mathcal {M}_m$ a time $\tau$ in the future, given that we started in cell $\mathcal {M}_n$. Putting together the pieces results in

(2.44)\begin{equation} R_E(g, \tau)\approx \sum_{n = 1}^{N} g(\boldsymbol{\sigma}^{[n]}) \mathbb{P}(\mathcal{M}_n) \left[ \sum_{m=1}^{N} g( \boldsymbol{\sigma}^{[m]}) [\exp(Q \tau)]_{mn} \right]. \end{equation}

For finite state space ergodic Markov processes in statistical equilibrium, continuous or discrete in time, the above equation is exact. In practice, we use (2.44) as an a posteriori check on the fidelity of a partition by computing the autocorrelation of a few select observables. When eigenvectors and eigenvalues of the underlying continuous operator exist, the ability of (2.44) to properly represent autocorrelations relies on the data-driven method's ability to represent the underlying operator's eigenvectors and eigenvalues.

In addition, all covariances and correlations are calculated using the above approximations. This review completes the discussion of how to approximate ensemble averages and covariances from the finite-volume discretization of the generator. However, it remains to be shown how to construct the matrix, and Markov states, from data. The construction of the generator is the subject of § 3.

2.6. Eigenvalues and eigenvectors of the generator

The generator has eigenvalues as well as (left and right) eigenvectors. The simplest eigenvectors are the ones associated with the eigenvalue $0$. Recall the requirement that the column sum of matrix entries should add up to zero, for each column. This requirement is a statement about the left eigenvector of $Q$. To see this denote $\boldsymbol {1}$ as the eigenvector of all $1$s; then

(2.45)\begin{equation} \boldsymbol{1}^T Q = \boldsymbol{0}^T = 0 \boldsymbol{1}^T. \end{equation}

The first equality is exactly the statement that the column sum of the matrix should add up to zero for each column and the second equality rewrites the zero vector as the zero scalar times the original vector. The latter inequality shows that it is an eigenvector of the system. The assumption that the $Q$ matrix is ergodic then implies that there is only one eigenvector corresponding to the eigenvalue 0. We shall now discuss the other eigenpairs.

2.7. Global Koopman eigenvectors and modes

We do not ask ‘can we predict an observable of interest?’, but rather, ‘what can we predict?’. The latter question is an emergent property of the system and captured by the Koopman eigenvectors of the underlying system. Those Koopman eigenvectors whose decorrelation time scales are long-lived constitute the most predictable features of the system.

We use the same terminology introduced in Williams et al. (Reference Williams, Kevrekidis and Rowley2015) to discuss the Koopman operator and its eigenvectors. Koopman eigenvectors are observables as well as left eigenvectors of the transition probability operator $\mathscr {T}^\tau$. For example, if $g_\lambda$ is a left eigenvector of $\mathscr {T}^\tau$ with eigenvalue ${\rm e}^{\lambda \tau }$ then we have the following:

(2.46)\begin{align} R_E(g_\lambda, \tau) &= \int_{\mathcal{M}} {\rm d}\boldsymbol{\mathscr{s}}'\,g_\lambda(\boldsymbol{\mathscr{s}}' ) \mathcal{I}( \boldsymbol{\mathscr{s}}' )\left[ \int_{\mathcal{M}} {\rm d}\boldsymbol{\mathscr{s}} \, g_\lambda(\boldsymbol{\mathscr{s}}) \mathscr{T}^\tau \delta ( \boldsymbol{\mathscr{s}} - \boldsymbol{\mathscr{s}}') \right] \end{align}

(2.47)\begin{align} &= \int_{\mathcal{M}} {\rm d}\boldsymbol{\mathscr{s}}' \,g_\lambda(\boldsymbol{\mathscr{s}}' ) \mathcal{I}( \boldsymbol{\mathscr{s}}' )\left[ \int_{\mathcal{M}} {\rm d}\boldsymbol{\mathscr{s}} \, g_\lambda(\boldsymbol{\mathscr{s}}) {\rm e}^{\lambda \tau} \delta ( \boldsymbol{\mathscr{s}} - \boldsymbol{\mathscr{s}}') \right] \end{align}

(2.48)\begin{align} &= {\rm e}^{\lambda \tau} \int_{\mathcal{M}} {\rm d}\boldsymbol{\mathscr{s}}' \, g_\lambda(\boldsymbol{\mathscr{s}}' )^2 \mathcal{I}( \boldsymbol{\mathscr{s}}' ) \end{align}

(2.49)\begin{align} &= {\rm e}^{\lambda \tau} \langle g_\lambda ^2 \rangle_E. \end{align}

Thus, the most useful Koopman eigenvectors, from a predictability standpoint, are those such that they decorrelate slowly in time, i.e. $\text {real}(\lambda ) \approx 0$, but additionally have an oscillatory component so that the ratio $\text {real}(\lambda )/\text {imaginary}(\lambda ) \approx 0$ holds.

If $\text {real}(\lambda ) = 0$ on a chaotic attractor, then we expect this to be the ‘trivial’ observable $g_\lambda (\boldsymbol{\mathscr {s}}) = c$ for a constant $c$. (The presence of pure-imaginary eigenvalues would imply the existence of observables that are predictable for arbitrary times in the future on a chaotic attractor.) Otherwise, we expect that $\text {real}(\lambda ) < 0$ for all eigenvalues of $Q$, i.e. we expect that all non-trivial observables will eventually decorrelate. This implies $\langle g_\lambda \rangle _E = 0$ since

(2.50)\begin{equation} \langle g_\lambda \rangle_E = \int_{\mathcal{M}} {\rm d} \boldsymbol{\mathscr{s}} \,g_\lambda( \boldsymbol{\mathscr{s}} ) \mathcal{I}(\boldsymbol{\mathscr{s}} ) = \lim_{\tau \rightarrow \infty} \int_{\mathcal{M}} {\rm d} \boldsymbol{\mathscr{s}} \,g_\lambda( \boldsymbol{\mathscr{s}}) \mathscr{T}^\tau \delta ( \boldsymbol{\mathscr{s}} - \boldsymbol{\mathscr{s}}') = \lim_{\tau \rightarrow \infty} {\rm e}^{\lambda \tau} g_\lambda(\boldsymbol{\mathscr{s}}') = 0, \end{equation}

where $\boldsymbol{\mathscr {s}}'$ is an arbitrary state on the attractor $\mathcal {M}$.

The following four statements about a Koopman eigenvectors $g_{\lambda }$ cannot hold simultaneously:

1. (i) the Koopman eigenvector satisfies the relation $g_\lambda (\boldsymbol{s}(t+\tau )) = {\rm e}^{\lambda \tau } g_\lambda (\boldsymbol{s}(t))$;

2. (ii) the Koopman eigenvector $g_\lambda$ is a continuous function of state space;

3. (iii) there exist an arbitrary number of near recurrences on the dynamical trajectory;

4. (iv) the eigenvalue associated with the Koopman eigenvector satisfies $\text {real}(\lambda ) < 0$.

The proof is as follows. Suppose that all four criteria are satisfied. Let $\boldsymbol{\mathscr {s}}'_n$ be near-recurrences of $\boldsymbol{\mathscr {s}}$ some sequence of times $\{\tau _n \}$, where $\tau _n \rightarrow \infty$, in the future so that $\| \boldsymbol{\mathscr {s}} - \boldsymbol{\mathscr {s}}'_n \| < \epsilon$ for some norm and all $n$ uniformly. Continuity of $g_\lambda$ with respect to the norm implies

(2.51)\begin{equation} | g_\lambda( \boldsymbol{\mathscr{s}}) - g_\lambda( \boldsymbol{\mathscr{s}}'_n) | < \delta \end{equation}

but $g_\lambda ( \boldsymbol{\mathscr {s}}'_n) = {\rm e}^{\lambda \tau _n} g_\lambda ( \boldsymbol{\mathscr {s}})$ by assumption, hence

(2.52)\begin{equation} | g_\lambda( \boldsymbol{\mathscr{s}})| |1 - {\rm e}^{\lambda \tau_n}| < \delta \end{equation}

which is a contradiction since $\tau _n$ can be made arbitrarily large and $\delta$ arbitrarily small. The non-existence of Koopman eigenvectors satisfying $g_\lambda (\boldsymbol{s}(t+\tau )) = {\rm e}^{\lambda \tau } g_\lambda (\boldsymbol{s}(t))$ over all of state space is corroborated by numerical evidence (Parker & Page Reference Parker and Page2020). In so far as a turbulent attractor is mixing one does not expect a finite-dimensional linear subspace for the Koopman eigenvectors (except for the constant observable). See Arbabi & Mezić (Reference Arbabi and Mezić2017) for a similar statement with regards to the Lorenz attractor. We take the above proof as a plausible argument for the use of a piecewise constant basis in the representation of Koopman eigenvectors.

For stochastic dynamical systems, one expects that the Koopman eigenvectors (the Koopman operator is defined in terms of the adjoint of the Fokker–Planck operator in that context) are more likely to be continuous functions of the state but no longer obey the relation $g_\lambda (\boldsymbol{s}(t+\tau )) = {\rm e}^{\lambda \tau } g_\lambda (\boldsymbol{s}(t))$. Heuristically, we expect better regularity because the stochastic noise in the dynamics acts as a diffusion in probability space, which smooths out non-smooth fields. For a stochastic differential equation

(2.53)\begin{equation} \dot{\boldsymbol{s}} = \boldsymbol{U}(\boldsymbol{s}) + \epsilon \boldsymbol{\xi}, \end{equation}

where $\epsilon$ is the noise variance and $\boldsymbol{\xi }$ is a $d$-dimensional white noise, the Koopman eigenvector evolves according to

(2.54)\begin{equation} \dot{g}_{\lambda}(\boldsymbol{s}(t)) = \lambda g_{\lambda}(\boldsymbol{s}(t)) + \epsilon \boldsymbol{\nabla} g_{\lambda}(\boldsymbol{s}(t)) \boldsymbol{\cdot} \boldsymbol{\xi}, \end{equation}

see (6.55) of Maćešić & Črnjarić-Žic (Reference Maćešić and Črnjarić-Žic2020). The formula above follows from Ito's lemma and using that $g_\lambda$ as an eigenvector with eigenvalue $\lambda$. The composition with the state variable is necessary because $g_\lambda : \mathbb {R}^d \rightarrow \mathbb {R}$, but one cannot consider the evolution of $g_\lambda$ independently from where it is being evaluated in state space. In the limit that the noise goes to zero, $\epsilon \rightarrow 0$, the gradient term, $\boldsymbol {\nabla } g_{\lambda }$, can go to infinity at particular points in state space, as would be expected in a one-dimensional stochastic dynamical system with a two-well potential. Consequently, pathologies are unexpected in linear systems.

As another point, although they are often called Koopman eigenfunctions when the dynamical system is a partial differential equation, they should perhaps more appropriately be called Koopman eigenoperators (in analogy to eigenvectors and eigenfunctions in lower-dimensional contexts), i.e. functionals that act on a state. The right eigenvectors of the transfer operator act as projection operators of an observable to Koopman eigenvectors. When a continuum of observables describe a field the result of the projection of this continuum is called a Koopman mode, see Williams et al. (Reference Williams, Kevrekidis and Rowley2015). Reference SouzaPart 2 goes through a concrete example, but we remain abstract here.

A continuum of observables indexed by spatial index $\boldsymbol{x}$ define statistical modes as

(2.55)\begin{equation} G_{\lambda}(\boldsymbol{x}) \equiv \int_{\mathcal{M}} {\rm d} \boldsymbol{\mathscr{s}}\, g^{\boldsymbol{x}}(\boldsymbol{\mathscr{s}}) v_{\lambda}( \boldsymbol{\mathscr{s}} ), \end{equation}

where $v_{\lambda }$ is a right eigenvector of the transfer operator, i.e.

(2.56)\begin{equation} \mathscr{T}^{\tau}v_{\lambda} = {\rm e}^{\lambda \tau} v_{\lambda}. \end{equation}

Equation (2.55) projects the observable $g^{\boldsymbol{x}}$ onto the appropriate Koopman eigenvector, due to biorthogonality of left and right eigenvectors. The object $G_{\lambda }(\boldsymbol{x})$, for a fixed $\lambda$ is a Koopman mode. Whether or not the set of Koopman eigenvectors form a complete basis so that $g^{\boldsymbol{x}} = \sum _{\lambda } G_\lambda (\boldsymbol{x}) g_{\lambda }$ is unclear. The implication is that an arbitrary observable $g^{\boldsymbol{x}}$ could be fundamentally unpredictable if it cannot be expressed as a sum of Koopman eigenvectors.

In the next section, we discuss the numerical approximation to Koopman eigenvectors.

2.8. Numerical approximation

The numerical Koopman eigenvectors are the left eigenvectors of the matrix $Q$, denoted by $\boldsymbol{g}_{\lambda }$ and their approximation as functionals acting on the state $\boldsymbol{\mathscr {s}} \in \mathcal {M}_n$ is given by

(2.57)\begin{equation} g_\lambda(\boldsymbol{\mathscr{s} }) \approx [\boldsymbol{g}_\lambda]_{n}, \end{equation}

where $[\boldsymbol{g}_\lambda ]_n$ is the $n$th component of the eigenvector $\boldsymbol{g}_{\lambda }$. Thus, we first determine which cell the state $\boldsymbol{\mathscr {s}}$ belongs to and then use the integer label to pick out the component of the eigenvector $\boldsymbol{g}_\lambda$. Thus, we use a piecewise constant approximation to the Koopman eigenvector.

3.1. Empirical construction

We start with an empirical construction of the generator. It suffices to focus on cell $j$ associated with the $j$th column of the matrix $Q_{ij}$. To calculate the empirical holding time distribution and empirical mean, we count how often we see cell $j$ before transitioning to cell $i \neq j$. For example, suppose that we have three cells, $j=1$, and consider the following sequence of integers given by the classifier applied to a time series with $\Delta t$ spacing in time:

(3.3)\begin{equation} 1, 1, 1, 2, 2, 1, 1, 3, 1, 2, 1, 1. \end{equation}

We group the sequence as follows:

(3.4)\begin{equation} (1, 1, 1), 2, 2, (1, 1), 3, (1), 2, (1, 1) \end{equation}

to determine the holding times. Thus, the holding times for cell 1 would be

(3.5)\begin{equation} 3 \Delta t, 2\Delta t, \Delta t, 2\Delta t \end{equation}

whose empirical expected value is $2 \Delta t$ implying a transition rate $1 / (2 \Delta t)$.

To calculate exit probabilities for cell $j$, we count how often we see transitions to cells $i$ and divide by the total number of transitions. In the example, to calculate the exit probabilities for cell $1$ into cell $2$ or $3$, we group them together as follows:

(3.6)\begin{equation} 1, 1, (1, 2), 2, 1, (1, 3), (1, 2), 1, 1. \end{equation}

Thus, we observed three exits, two of which went to state $2$ and one to state $3$; hence, the exit probabilities are $E_{21} = 2/3$ and $E_{31} = 1/3$.

The rest are constructed analogously to produce the matrix

(3.7) \begin{align} Q = E R = \begin{bmatrix} -1 & 1 & 1 \\ { }2/3 & -1 & 0 \\ { }1/3 & 0 & -1 \end{bmatrix} \begin{bmatrix} \dfrac{1}{2 \Delta t} & 0 & 0 \\ 0 & \dfrac{2}{3 \Delta t} & 0 \\ 0 & 0 & \dfrac{1}{\Delta t} \end{bmatrix} = \dfrac{1}{\Delta t}\begin{bmatrix} -1/2 & { }2/3 & { } 1 \\ { }1/3 & -2/3 & 0 \\ { }1/6 & 0 & -1 \end{bmatrix} . \end{align}

We give an alternative method of constructing $Q$ in Appendix C where we interpret operations from an algebraic perspective. As described, the generator is only accurate to order $\Delta t$ since we do not interpolate in time to find the ‘exact’ holding time. We do not preoccupy ourselves with improving this since we believe that the primary source of error comes from finite sampling effects.

In certain cases, the data-driven construction given above can be obtained directly from a data-driven construction of the Perron–Frobenius operator. For example, when a time series is sampled with a time step of $\Delta t$ and a trajectory only spends one time step $\Delta t$ within each cell, then the $Q$ matrix is discretely related to the data-driven construction Perron–Frobenius operator, $\mathcal {F}^{[\Delta t]}$, through the relation $\mathcal {F}^{[\Delta t]} = \mathbb {I} + \Delta t Q$ where $\mathbb {I}$ is the identity matrix. This happens because the only transition probabilities observed are exit probabilities. Furthermore, the relation $\mathcal {F}^{[\Delta t]} \approx \mathbb {I} + \Delta t Q$ also holds in the limit that a trajectory spends sufficient time within each cell $j$, i.e. the average holding times are much greater than the time step size. See Appendix D for concrete examples using the simple harmonic oscillator. The advantage of the continuous-time formulation comes from uncertainty-quantification which yields sensible uncertainty estimates in both limits.

In the following section, we augment the empirical construction with uncertainty estimates based on finite sampling and a Bayesian framework.

3.2. Bayesian construction

Our goal is to quantify the uncertainty of the discrete generator's matrix entries due to finite sampling effects from data-driven methods. To this end, we use a Bayesian method, which requires choices for likelihood functions and prior distributions. In the case of a generator matrix, we require four ingredients:

1. (i) a likelihood distribution for the holding times of each cell;

2. (ii) a prior distribution for the transition rates associated with the holding times;

3. (iii) a likelihood distribution for the exit probabilities of each cell;

4. (iv) a prior distribution for the probability values associated with the exit probabilities.

The likelihood distribution is the distribution we believe we are sampling from, and the prior distribution encapsulates our uncertainty with respect to the parameters of the likelihood function. Bayesian methods are computationally efficient under certain likelihood and prior distribution assumptions (Gelman et al. Reference Gelman, Carlin, Stern, Dunson, Vehtari and Rubin2013). In particular, if the prior distribution and the posterior distribution are from the same family of distributions, then it is only necessary to update the parameters of the prior/posterior distribution. For a fixed likelihood function, in special cases, it is possible to determine a prior distribution such that the posterior distribution comes from the same family as the prior. Such distributions are called conjugate priors (conjugate with respect to the likelihood), and we shall use them.

We make choices for our likelihood function compatible with that of a continuous-time Markov process:

1. (i) the likelihood distribution for the holding times is exponentially distributed with rate parameter $\lambda _i$ for each cell independently;

2. (ii) the likelihood distribution for exit probabilities is a multinomial distribution with parameters $\boldsymbol{p} \in [0, 1]^{N-1}$ satisfying the relation $\sum _{i=1}^{N-1} p_i = 1$ for each cell independently.

Both the exponential distribution and the multinomial distribution have known conjugate priors:

1. (i) the conjugate prior distribution for the rate parameter of the exponential distribution is distributed according to the Gamma distribution with parameters $(\alpha, \beta )$, denoted by $\varGamma (\alpha, \beta )$;

2. (ii) the conjugate prior distribution for the exit probabilities of the multinomial distribution comes from a Dirichlet distribution with parameter vector $\boldsymbol {\alpha }$ of length $N-1$, which we denote by Dirichlet ($\boldsymbol {\alpha }$).

Thus, the posterior distributions then come from the same family as the prior distributions, e.g. posterior distributions are $\varGamma (\alpha ', \beta ')$ and Dirichlet ($\boldsymbol{\alpha }'$) upon the observation of data; see Gelman et al. (Reference Gelman, Carlin, Stern, Dunson, Vehtari and Rubin2013) and the section below for an example. Under the construction of this section, a $3 \times 3$ matrix will always be of the form

(3.8)\begin{equation} Q = \begin{bmatrix} - 1 & [{\boldsymbol{D}}_{2}]_1 & [{\boldsymbol{D}}_{3}]_1 \\ [{\boldsymbol{D}}_{1}]_1 & -1 & [{\boldsymbol{D}}_{3}]_2 \\ [{\boldsymbol{D}}_{1}]_2 & [{\boldsymbol{D}}_{2}]_2 & - 1 \end{bmatrix} \begin{bmatrix} G_1 & 0 & 0 \\ 0 & G_2 & 0 \\ 0 & 0 & G_3 \end{bmatrix} \end{equation}

where $G_i \sim \varGamma (\alpha _i, \beta _i)$, $\boldsymbol {D}_{i} \sim \text {Dirichlet}(\boldsymbol {\alpha }_i)$ and $[\boldsymbol {D}_{i}]_j$ denotes the $j$th component of the random vector $\boldsymbol {D}_{i}$.

As mentioned in the previous paragraphs, conjugate priors greatly expedite the Bayesian update procedure; only the parameters of the Gamma and Dirichlet distributions need to be updated to construct a posterior distribution. The parameters $(\alpha _i, \beta _i)$ and $\boldsymbol {\alpha }_i$ are updated according to Bayes’ rule for each column upon data acquisition. For example, suppose that we have observed the following empirical counts associated with cell $i$:

1. (i) $M$ exits from cell $i$;

2. (ii) $[\boldsymbol {M}]_j$ exits from cell $i$ to cell $j$;

3. (iii) $\hat {T}_1, \hat {T}_2,\ldots, \hat {T}_M$ empirically observed holding times;

and that we start with $(\alpha ^0, \beta ^0)$ and $\boldsymbol {\alpha }^0$ as the parameters for our prior distributions. The relation $\sum _j [\boldsymbol {M}]_j = M$ holds. (There can be an ‘off-by-one’ error here which we ignore for presentation purposes.) The posterior distribution parameters $(\alpha ^1, \beta ^1)$ and $\boldsymbol {\alpha }^1$ are

(3.9a–c)\begin{equation} \alpha^1 = \alpha^0 + M { , }\quad \beta^1 = \beta^0 + \sum_{i}^M \hat{T}_i \quad {\rm and}\quad \boldsymbol{\alpha}^1 = \boldsymbol{\alpha}^0 + \boldsymbol{M} . \end{equation}

In the limit that $\alpha ^0 , \beta ^0$ and $|\boldsymbol {\alpha }^0|$ go to zero, then the empirical approach from the previous section agrees with the expected value from the Bayesian approach.

The current approach is one of many approaches to constructing matrices with quantified uncertainties. See Singhal & Pande (Reference Singhal and Pande2005) and Trendelkamp-Schroer et al. (Reference Trendelkamp-Schroer, Wu, Paul and Noé2015) for other Bayesian examples of calculating uncertainties and Scherer et al. (Reference Scherer, Trendelkamp-Schroer, Paul, Pérez-Hernández, Hoffmann, Plattner, Wehmeyer, Prinz and Noé2015) for numerical implementations. We comment that the present methodology also allows for uncertainty quantification of generators that satisfy detailed balance, see Appendix B. Other probabilistically inspired methods of estimating generator entries consist of methods such as the expectation maximization method in Otto, Peitz & Rowley (Reference Otto, Peitz and Rowley2023), which yield maximum likelihood estimates.

The current uncertainty quantification is imperfect (e.g. when holding times do not follow an exponential distribution or the system is not Markovian over infinitesimal steps). Still, we hold the position that some quantification of uncertainty is better than none. One of the benefits of the construction of this section is that it is robust to infinite temporal resolution over a fixed period and, hence, is consistent with data that comes from a continuous time process. We use uncertainty quantification to dismiss spurious results rather than increase confidence in the correctness of an inference. Furthermore, the assumptions that we made for the posterior and prior distribution still yield the same empirical construction from (3.1) upon using uninformative priors. In the large data limit, the Bayesian update procedure eventually yields a sharply peaked distribution that converges to a Gaussian.

An example now follows.

4.1. Fixed point partition

We choose the classic parameter values $r=28$, $\sigma = 10$ and $b= 8/3$ for the Lorenz system, which is known to exhibit chaotic solutions. Construction of the generator is automated through the methodology of § 3 upon choosing the Markov states $\boldsymbol{\sigma }^{[n]}$ and a classifier $\mathcal {C}$. We use the following fiction to guide our choices.

It is said that the coherent structures of a flow organize and guide the dynamics of chaos (Cvitanović Reference Cvitanović2013). As a trajectory wanders through state space, it spends a disproportionate time near coherent structures and inherits their properties. The coherent structures then imprint their behaviour on the chaotic trajectory, manifesting in ensemble averages. Thus, chaotic averages are understood in terms of transitions between simpler structures.

This picturesque story motivates the use of fixed points, as Markov states,

(4.5)$$\begin{gather} \boldsymbol{\sigma}^{[1]} = [-\sqrt{72} ,-\sqrt{72}, 27], \end{gather}$$

(4.6)$$\begin{gather}\boldsymbol{\sigma}^{[2]} = [0, 0, 0], \end{gather}$$

(4.7)$$\begin{gather}\boldsymbol{\sigma}^{[3]} =[\sqrt{72} ,\sqrt{72}, 27], \end{gather}$$

and partitioning state space according to the closest fixed point,

(4.8)\begin{equation} \mathcal{C}(\boldsymbol{\mathscr{s}}) = \begin{cases} 1 & \text{if } \| \boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[1]} \| < \| \boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[2]} \| \text{ and } \| \boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[3]} \|\\ 2 & \text{if } \| \boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[2]} \| < \| \boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[3]} \| \text{ and } \| \boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[1]} \|\\ 3 & \text{if } \| \boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[3]} \| < \| \boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[1]} \| \text{ and } \| \boldsymbol{\mathscr{s}} - \boldsymbol{\sigma}^{[2]} \| \end{cases} \end{equation}

where $\| \cdot \|$ denotes the standard Euclidean norm. The classifier determines the partition by associating a trajectory with the closest fixed point. This partitioning strategy is the intersection of the chaotic attractor, $\mathcal {M}$, with a Voronoi tessellation over the full state space, $\mathbb {R}^3$. Figure 1 shows the partition induced by this choice. The regions are colour-coded according to the closest fixed points.

Figure 1. Lorenz fixed point partition. Here, we show the emerging partition from several angles. The colours correspond to the different partitions associated with trajectories that are ‘closest’ to a given fixed point.

We construct a time series from the Lorenz equations using a fourth-order Runge–Kutta time stepping scheme with time step $\Delta t = 5 \times 10^{-3}$. We take the initial condition to be $(x(0), y(0), z(0)) = (14, 20, 27)$ and integrate to time $T = 10^5$, leading to $2 \times 10^7$ time snapshots. At each moment in time, we apply the classifier to create a sequence of integers representing the partition dynamics. (One can think of this as defining a symbol sequence.) Figure 2 visualizes this process.

Figure 2. Lorenz fixed point partition Markov chain. The dynamics of the $x, y, z$ variables are shown in (a), and the associated partition dynamics is shown in (b). As a dynamical trajectory moves through state space, it is labelled according to its proximity to the closest fixed point.

From the sequence of integers, we apply the method from § 3.2 to construct the data-driven approximation to the generator with quantified uncertainty. For our prior distribution, we use an uninformative prior – i.e. initial parameters $\alpha = \beta = 0$ for the Gamma distribution and $\boldsymbol {\alpha } = \boldsymbol {0}$ for the Dirichlet distribution – so that the mean of the random matrix agrees with the empirical construction from 3.1. The mean for each entry of the generator (reported to two decimal places) is

(4.9)\begin{align} \langle Q \rangle = \begin{bmatrix} -1.17 & 1.93 & 0.65 \\ 0.52 & -3.86 & 0.52 \\ 0.65 & 1.93 & -1.17 \end{bmatrix} \approx \begin{bmatrix} -1.0 & 0.5 & 0.55 \\ 0.45 & -1.0 & 0.45 \\ 0.55 & 0.5 & -1.0 \end{bmatrix} \begin{bmatrix} 1.17 & 0.0 & 0.0 \\ 0.0 & 3.86 & 0.0 \\ 0.0 & 0.0 & 1.17 \end{bmatrix}, \end{align}

where we have decomposed the matrix into the exit probability matrix and the diagonal rate matrix on the right-hand side.

In the matrix decomposition, the off-diagonals of the left matrix correspond to the exit probabilities, and the right matrix is the rate matrix, whose entries are the inverse of the time spent within a partition. From the latter matrix, we see that trajectories spend less time in the partition associated with the zero fixed point (blue) since $1/3.86 < 1/1.17$. The apparent symmetry in the matrix results from the truncation to two decimal places and the abundance of data. In Appendix A, we show how to incorporate symmetries of the Lorenz equation and report ensemble mean statistics.

The utility of using a random matrix to represent uncertainty is summarized in figure 3. The distribution of each matrix entry for various subsets of time is displayed. Using fewer data (represented by a shorter gathering time, $T$) results in significant uncertainty for the matrix entries. Additionally, using unconnected subsets of time demonstrates an apparent convergence of matrix entries.

Figure 3. Lorenz fixed point partition distributions of the generator. The uncertainty estimates for the entries of the $3 \times 3$ generator are shown in the above figure. A one-to-one correspondence exists between the distributions in the panel and the matrix entries. The different coloured distributions within a panel represent different estimates of the entries based on the amount of available data, here presented in terms of the simulation time of the Lorenz system. We see that as we increase the time interval of the simulation and thus have more data, we become more confident about the matrix entries. Furthermore, the distributional spreads overlap with one another.

We are now in a position to calculate statistical quantities. For simplicity, we only report first-, second- and third-order moments calculated from the mean value of the generator, $\langle Q \rangle$. The steady-state distribution of $\langle Q \rangle$, corresponding to eigenvalue $\lambda = 0$, is reported to two decimal places as

(4.10)\begin{equation} [\mathbb{P}(\mathcal{M}_1) , \mathbb{P}(\mathcal{M}_2), \mathbb{P}(\mathcal{M}_3) ] \approx [0.44, 0.12, 0.44] \end{equation}

from whence we calculate the steady state statistics for any observable using the approximations in § 2.4 and the Markov states $\boldsymbol{\sigma }^{[n]}$ for $n = 1, 2,3$. Explicitly, the ensemble average of the observables,

(4.11)\begin{equation} g^{[1]}(\boldsymbol{\mathscr{s}})=\mathscr{s}_{3} = \mathscr{z}, \quad g^{[2]}(\boldsymbol{\mathscr{s}})= (\mathscr{s}_{3})^2 = \mathscr{z}^2,\quad \text{or}\quad g^{[3]}(\boldsymbol{\mathscr{s}})=(\mathscr{s}_{1})^2\mathscr{s}_{3} = \mathscr{x}^2 \mathscr{z} \end{equation}

is approximated via (2.39), repeated here for convenience,

(4.12)\begin{equation} \langle g^{[j]} \rangle_E = g^{[j]}(\boldsymbol{\sigma}^{[1]}) \mathbb{P}(\mathcal{M}_1) + g^{[j]}(\boldsymbol{\sigma}^{[2]}) \mathbb{P}(\mathcal{M}_2) + g^{[j]}(\boldsymbol{\sigma}^{[3]}) \mathbb{P}(\mathcal{M}_3)\quad \text{for each } j \end{equation}

to yield

(4.13)$$\begin{gather} \langle \mathscr{z} \rangle_E \approx 27 \times 0.44 + 0 \times 0.12 + 27 \times 0.44 \approx 24, \end{gather}$$

(4.14)$$\begin{gather}\langle \mathscr{z}^2 \rangle_E \approx 27^2 \times 0.44 + 0^2 \times 0.12 + 27^2 \times 0.44 \approx 642, \end{gather}$$

(4.15)$$\begin{gather}\langle \mathscr{x}^2 \mathscr{z} \rangle_E \approx \left(-\sqrt{72} \right)^2 \times 27 \times 0.44 + 0^3 \times 0.12 + \left(\sqrt{72}\right)^2 \times 27 \times 0.44 \approx 1711 . \end{gather}$$

Table 1 shows the result from both the temporal and ensemble average (using full machine precision for computations). There is a correspondence for all averages, with the most significant discrepancy being those involving $y^2$ terms, for which the relative error is within $25\,\%$. The fixed points of a dynamical system are unique in that they satisfy all the same dynamical balances of a statistically steady state, e.g. chaotic trajectories, periodic orbits and fixed points of the Lorenz equation satisfy the relation

(4.16)\begin{equation} \langle \mathscr{x} \mathscr{y} \rangle = b \langle \mathscr{z} \rangle, \end{equation}

and more generally

(4.17)\begin{equation} 0 = \langle \boldsymbol{U}(\boldsymbol{\mathscr{s}}) \boldsymbol{\cdot} \boldsymbol{\nabla} g \rangle , \end{equation}

for any bounded and differentiable observable $g$, where $\boldsymbol {U}$ is vector field defined by the right-hand side of the Lorenz equations, and the averaging brackets $\langle \cdot \rangle$ are defined over the trajectory. Thus, an accurate estimation of $\langle \mathscr {x} \mathscr {y} \mathscr {z}^{n-1}\rangle$ only depends on an accurate estimate for $\langle \mathscr {z}^n \rangle$ for the fixed points of the Lorenz equation (using $g(\boldsymbol{\mathscr {s}}) = \mathscr {z}^{n}$). In this case the ‘closure problem’ in fluid mechanics is a boon rather than a curse. A good representation of a lower-order moment automatically yields a good representation of a higher-order moment through the closure ‘problem’.

Table 1. Empirical moments of the Lorenz attractor. A comparison between ensemble averaging and time averaging.

Although we focused on moments, one can compare the statistics of any observable, e.g.

(4.18a,b)\begin{equation} \langle \mathscr{z} \log(\mathscr{z}) \rangle_E \approx 78.4\quad\text{and}\quad \langle z \log(z) \rangle_T \approx 76.0, \end{equation}

where we used $z \log (z) \rightarrow 0$ as $z \rightarrow 0$. By symmetry, one expects

(4.19)\begin{equation} \langle x \rangle = \langle y \rangle = \langle x z \rangle = \langle y z \rangle = \langle y y y \rangle = \langle xx y\rangle = \langle x yy \rangle = \langle x zz \rangle = \langle y zz \rangle = 0 \end{equation}

but finite sampling effects prevent this from happening. As done in Appendix A, incorporating the symmetries allows ensemble calculations to achieve this to machine precision.

In addition to containing information about steady-state distributions, the generator $Q$ provides temporal information: autocorrelations and the average holding time within a given cell. We show the autocorrelation of six observables,

(4.20a–d)\begin{gather} g^{[1]}(\boldsymbol{\mathscr{s}}) = \mathscr{x} {,}\quad g^{[2]}(\boldsymbol{\mathscr{s}}) = \mathscr{y} {,}\quad g^{[3]}(\boldsymbol{\mathscr{s}}) = \mathscr{z} {,}\quad g^{[4]}(\boldsymbol{\mathscr{s}}) = \begin{cases} 1 & \text{if } \mathcal{C}(\boldsymbol{\mathscr{s}}) = 1\\ 0 & \text{otherwise } \end{cases}, \end{gather}

(4.21a,b)\begin{gather} g^{[5]}(\boldsymbol{\mathscr{s}}) = \begin{cases} 1 & \text{if } \mathscr{x} > 0\\ -1 & \text{if } \mathscr{x} < 0\\ 0 & \text{otherwise} \end{cases}\quad\text{and}\quad g^{[6]}(\boldsymbol{\mathscr{s}}) = \begin{cases} 1 & \text{if } \mathcal{C}(\boldsymbol{\mathscr{s}}) = 2\\ 0 & \text{otherwise } \end{cases} \end{gather}

in figure 4, which are calculated via (2.36) and (2.44), with appropriate modifications accounting for means and normalizing the height to one. Here, we see both the success and limitations of the method at capturing autocorrelations. In general, the decorrelation of an observable is captured by the Markov model if it is approximately constant within a given cell, e.g. the observables $g^{[4]}$ and $g^{[5]}$. However, sometimes it is possible to do ‘well’, such as $g^{[1]}$ or $g^{[2]}$, despite not being constant within a region.

Figure 4. Lorenz autocorrelations generator versus time series. Six autocorrelations of observables are shown. The transparent purple line is calculated from the generator, and the black line is calculated from the time series. Even a coarse partition can capture observables $g^{[4]}$ and $g^{[5]}$ but struggles with oscillatory correlations.

The inability to capture the autocorrelation of $g^{[6]}$, which is constant within $\mathcal {M}_2$, is partially due to the holding time distribution being far from exponentially distributed. To see this mode of failure, we plot the holding time distribution of the cells in figure 5. We show several binning strategies of the distribution to demonstrate the ability of an exponential distribution to capture quantiles of the empirical holding time distribution.

Figure 5. Lorenz fixed point partition holding times. An underlying assumption of using a generator for a given partitioning strategy is that the time spent in a cell is exponentially distributed. Here, we examine quantiles of the holding time distribution for a cell as given by the different binning numbers. The black dots correspond to the equivalent exponential distribution quantile, where the generator gives the rate parameter.

Depending on the time scale of interest, the $\mathcal {M}_1$ and $\mathcal {M}_3$ cells are approximately exponentially distributed, although they become fractal-like in terms of the distribution of holding times. In contrast, the holding time distribution of cell $\mathcal {M}_2$ is far from exponentially distributed upon refining the bins of the histogram. This calls into question using a Markovian, i.e. ‘memoryless’ model. There is an inherent assumption in the construction of the generator that transition probabilities are independent of the amount of time spent in a particular subset of state space. A better statistical model would incorporate exit probabilities conditioned on the time spent in a cell. Stated differently, memory is necessary to correctly reproduce the autocorrelation of a coarse discretization, see Lin et al. (Reference Lin, Tian, Perez and Livescu2023); however, the approach taken here is to view a given discrete representation as inherently imperfect and subject to improvement upon refinement of a partition.

Figure 6 summarizes the resulting statistical dynamics, where the generator and transition probabilities define a graph structure. The graph structure contains information about the topological connectivity between different regions of state space and the ‘strength’ of connectivity over different time scales as encapsulated by the transition probabilities.

Figure 6. Lorenz fixed point partition graph. The generator (a) and transition probabilities over several time scales are visualized as a graph. The transition probabilities change depending on the time scale.

In the next section we consider two different partitioning strategies of the Lorenz equations, over the same dataset.

4.2. The AB and sampling partitions

We now show how the partition choice affects steady-state statistics and temporal autocorrelations. We choose two partitions, both using 128 cells. The first partition uses 128 evenly spaced (in time) points on the AB periodic orbit of the Lorenz equation, see Viswanath (Reference Viswanath2003), to define the cells of the partition, and our second partition uses 128 points sampled on the attractor. We call the former the ‘AB partition’ and the latter the ‘sampling partition’. Both methods define cells by first labelling the cell centres with integers within 1 to 128, computing the distance to all 128 points, and then assigning a cell by picking out the cell centre with the smallest distance, similar to what was done for the three fixed points in (4.8). More succinctly, the points to define a Voronoi tessellation of the domain using the same time series from § 4.1. In both cases, we use an uninformative prior to construct the generator. See figure 7 for a visualization of the two partitioning strategies.

Figure 7. Lorenz AB and sampling partitions. In (a–c), we show the partition of the Lorenz attractor according to 128 points on the AB periodic orbit of Lorenz, and in (d–f), we show the sampling partition associated with 128 randomly chosen points. The black dots are the Markov states (cell centres) associated with the cells of each partition.

The AB periodic yields partitions that are thin wedges. In this case, further refinement by using more points on the periodic orbit does not yield a refinement strategy that converges to the statistics of the Lorenz equations. For example, the maximum value of $\mathscr {z}$ can never be approximated by points on the AB periodic orbit. On the other hand, the sampling partition yields a partition with disproportionate cell sizes and clusters in regions of high probability on the attractor. We expect that further refinement by choosing more random points yields a higher fidelity partition, but have no proof on this matter. This latter strategy is readily employable on any data set and can be thought of as a ‘go-to’ strategy in the absence of system knowledge.

We show two statistical measures to assess the quality of the partition. The first is in table 2. We see that the statistics of the AB partition are farther from the attractor for many of the variables than the sampling partition. This discrepancy is because the AB partition uses cell centres that are all on the periodic orbit and thus will represent a compromise between the periodic orbit statistics and the attractor statistics. In figure 8, we show the autocorrelations using the two partitions where the observables are defined in (4.20a–d). We see that the AB periodic orbit partition better captures oscillatory behaviour at the expense of the decay behaviour in the $g^{[1]}$ and $g^{[2]}$ autocorrelations. The $g^{[3]}$ autocorrelations, exhibiting less decay than its $g^{[1]}$ and $g^{[2]}$ counterparts are better captured by the periodic orbit partition. The oscillatory behaviour in $g^{[3]}$ is also captured by the sampling partition, albeit with more dissipation. We expect that using more points on the attractor for the sampling partition would yield better correspondence with the attractor. See Giorgini, Souza & Schmid (Reference Giorgini, Souza and Schmid2023) for a data-driven method that uses the methods developed here as a baseline for improving convergence of autocorrelations in coarse settings.

Table 2. Empirical moments of the Lorenz attractor. A comparison between ensemble averaging with two different partitions, ‘AB’ and ‘sampling’, and time averaging.

Figure 8. Autocorrelation for observables of the Lorenz system. The AB partition (red) captures oscillatory behaviour, while the sampling partition (blue) better captures the decay of the temporally obtained autocorrelation (black).