Representing turbulent statistics with partitions of state space. Part 1. Theory and methodology

Abstract

This is the first of a two-part paper. We formulate a data-driven method for constructing finite-volume discretizations of an arbitrary dynamical system's underlying Liouville/Fokker–Planck equation. A method is employed that allows for flexibility in partitioning state space, generalizes to function spaces, applies to arbitrarily long sequences of time-series data, is robust to noise and quantifies uncertainty with respect to finite sample effects. After applying the method, one is left with Markov states (cell centres) and a random matrix approximation to the generator. When used in tandem, they emulate the statistics of the underlying system. We illustrate the method on the Lorenz equations (a three-dimensional ordinary differential equation) saving a fluid dynamical application for Part 2 (Souza, J. Fluid Mech., vol. 997, 2024, A2).

1. Introduction

Often, the goal of modelling a complex system is not to determine the dynamical equations but to construct models that converge in distribution to relevant statistics. In the context of turbulence modelling, this can be viewed as one of the goals of a large-eddy simulation, where subsets of statistics (often the kinetic energy spectra) are compared with that of direct numerical simulation. Similarly, in Earth systems modelling, the unpredictability of weather patterns over long time scales necessitates the development of nonlinear models that are queried for relevant statistics. Thus, the models are not meant to converge to dynamical trajectories but rather converge in distribution to target observables.

The present work is motivated by the need to construct simplified statistical models of complex physical phenomena such as turbulence. We take on a dynamical systems view of turbulence original to Hopf (1948), complemented by Lorenz (1963) and found in its modern form in Cvitanović et al. (2016). Thus, the approach is to develop a direct discretization of the statistics associated with chaotic or turbulent dynamics, which we assume to be mixing and associated with a fractal in state space.

There exist many types of discretizations that directly target the statistics, which here means a discretization of the underlying Liouville equation (deterministic dynamics), Fokker–Planck equation (stochastic dynamics), Perron–Frobenius/transfer operator (discrete-time dynamics) or Koopman operator (adjoint of the Perron–Frobenius/transfer operator). Discretizations methods include that of Ulam (1964), Dellnitz, Froyland & Junge (2001), Dellnitz et al. (2005) or, for the stochastic Lorenz equations, Allawala & Marston (2016). Furthermore, there exist efficient extensions of Ulam's method for multidimensional systems such as the box-refinement methods of Dellnitz & Junge (1999), Dellnitz et al. (2001) or the sparse Haar tensor basis of Junge & Koltai (2009). Modern methods take on an operator theoretic plus data-driven approach, leading to Koopman operators that are measure preserving Colbrook (2023), which build off of earlier work such as Schmid (2010) and Rowley et al. (2009). Data-driven construction of the Perron–Frobenius operator is reviewed by Klus et al. (2016), Giannakis (2019) and Fernex, Noack & Semaan (2021). Convergence guarantees under various assumptions are found in, for example, Froyland (1997), Das, Giannakis & Slawinska (2021), Colbrook & Townsend (2023) and Schütte, Klus & Hartmann (2022). Furthermore, modern methods are beginning to use deep learning in order to learn optimal nonlinear dictionaries for observables, see Constante-Amores, Linot & Graham (2023) and Bittracher et al. (2023).

The method presented here can be viewed as a subset of the general method for constructing Koopman operators using the extended dynamic mode decomposition method (Williams, Kevrekidis & Rowley 2015); however, we make a particular choice of a nonlinear dictionary that allows for computational expedience (we do not need to explicitly calculate pseudoinverses), geometric interpretability (which allows for a natural method to increase the size of the nonlinear dictionary) and the quantification of uncertainty due to finite sampling effects. The first two changes allow for much larger-scale computations than what has been previously achieved. Furthermore, we take an approach that closely mirrors a combination of Froyland et al. (2013) and Fernex et al. (2021). When addressing statistics of a partial differential equation, we take a field-theoretic perspective such as that of Hopf (1952). The goal is to construct a discretization of the generator (continuity/Fokker–Planck operator), similar to Rosenfeld et al. (2021). The method of constructing the generator herein is motivated by finite-volume discretizations of advection–diffusion equations and yields an approximation that can be interpreted as a continuous-time Markov process with finite state space.

The rest of the paper is organized as follows. In § 2, we discuss the underlying theory and approximations. The approach is heavily inspired by Hopf (1952) and Cvitanović et al. (2016). The primary idea is to discretize the equations for the statistics (an ‘Eulerian’ quantity) by using the equation for the dynamics (a ‘Lagrangian’ quantity). Further approximations are then made to calculate observables of interest.

In § 3, we introduce a data-driven method with quantified uncertainties for calculating the approximate generator. The technique can be applied to arbitrarily long time-series data, dynamical systems with a large state space, and further provides uncertainty estimates on the entries of the discretized generator.

Section 4 gives an example of utilizing the method on Lorenz (1963). One sees that even a coarse discretization of statistics captures statistical features of the original system.

For those simply interested in what the calculations enable, §§ 2–3 are safely skipped in favour of § 4 or Part 2 of this series (Souza 2024). Furthermore, Appendices A, B, C and D expand the text by discussing symmetries, matrix properties of the generator, an algebraic interpretation of constructing the generator alongside the connection to dynamic mode decomposition and convergence of the data-driven method with respect to the simple harmonic oscillator, respectively.

2. Theory

In the following sections, we outline the framework for calculations involving the Liouville equation of a dynamical system. The focus is not on the ‘Lagrangian’ view (as given by the dynamics) but rather an ‘Eulerian’ one (given by statistics). Consequently, we introduce notation to help distinguish the perspectives inherent in the two viewpoints. We show how to go from the continuous to the discrete and provide formulae.

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 1963). It is, therefore, more natural to study the statistical evolution of probability densities as in Hopf (1952). 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 (1952). 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 (2023b) and Giorgini et al. (2024). 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 (2021) 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 (2014) and Corwin & Shen (2020)), 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 (1919), DeWitt (1972) and Albeverio & Mazzucchi (2016)).

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 (2002) and Cowieson & Young (2005). 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 (2019). Whether or not such a limit generally agrees with a Sinai–Ruelle–Bowen measure (see for example Young (2002), Blank (2017) and Araujo (2023)) 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 2011).

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 (1989), where an approximation to an Ornstein–Uhlenbeck process is constructed. See appendices C.2 and C.3 of Souza et al. (2023b) 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 2005). 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. (2015) 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 2020). 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ć (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 (2020). 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. (2015). Part 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. Methodology

In this section, we outline the general approach to constructing the approximate generator in terms of trajectory data. The most critical component of a discretization comes from defining a classifier $\mathcal {C}: \mathbb {R}^d \rightarrow \{1, 2,\ldots, N \}$ which maps an arbitrary state $\boldsymbol{\mathscr {s}} \in \mathbb {R}^d$ to an integer $n \in \{1, 2,\ldots, N \}$. This function implicitly defines a cell through the relation

(3.1)\begin{equation} \mathcal{M}_j = \{ \boldsymbol{\mathscr{s}} : \mathcal{C}(\boldsymbol{\mathscr{s}}) = j \text{ for each } \boldsymbol{\mathscr{s}} \in \mathbb{R}^d \} \cap \mathcal{M}, \end{equation}

and thus we identify an integer $j$ with a cell $\mathcal {M}_j$. The intersection with the manifold $\mathcal {M}$ is critical to the methodology's success.

Furthermore, the Markov states (cell centres) are chosen to satisfy $\mathcal {C}(\boldsymbol{\sigma }^{[n]}) = n$ for each $n \in \{1,\ldots, N \}$. There is an extraordinary amount of freedom in defining the classifier, and we will go through examples in § 4. We also comment on practical considerations and generalizations in § 5. One can simultaneously solve for a classifier and Markov states using a $K$-means algorithm (see Lloyd (1982)), but we do not wish to restrict ourselves to that choice here. The classifier ‘classifies’ (in the machine learning sense) different flow states with integers as the category labels. The classifier serves as a particular choice of the nonlinear dictionary (basis) for observables of the system; see Appendix C and Klus et al. (2016) for a discussion on this point. Choosing a good basis for observables is difficult and may not even exist, but deep learning methods combined with novel loss functions hold promise for their discovery (Bittracher et al. 2023; Constante-Amores et al. 2023). For now, we will assume that such a function is given and focus on constructing the generator $Q$.

The classifier $\mathcal {C}$ transforms dynamical trajectories into sequences of integers, which we interpret as the realization of a Markov process with finite state space. At this stage, traditional methods can be employed to construct a transfer/Perron–Frobenius operator from data, see Klus et al. (2016) and Fernex et al. (2021). Given that our interest is in constructing a continuous time Markov process and quantifying the uncertainty of matrix entries due to finite sampling effects, the algorithm will modify the traditional approaches.

To construct $Q$, two quantities must be calculated for each cell:

1. (i) the holding times – the amount of time a dynamical trajectory stays in cell $\mathcal {M}_n$ before exiting;

2. (ii) the exit probabilities – the probability of moving from cell $\mathcal {M}_j$ to $\mathcal {M}_i$ upon exiting the cell $\mathcal {M}_j$.

Let $T_j$ be the distribution of holding times associated with cell $j$ and $E_{ij}$ denote entries of the exit probability matrix for $i\neq j$. By convention we take $E_{ii} = -1$ so that $\sum _i E_{ij} = 0$ and for each $j$. The entries of the matrix $Q_{ij}$ are constructed as follows:

(3.2)\begin{equation} Q_{ij} = E_{ij} / \langle T_j \rangle, \end{equation}

where $\langle T_j \rangle$ denotes the expected value of the holding time distribution of cell $j$. The matrix $Q$ is decomposed into the product of the exit probability matrix $E$ and inverse holding time matrix $R \equiv \text {Diagonal}({\boldsymbol {T}})^{-1}$, $Q = E R$ when numerically constructing the matrix entries. The data-driven construction of the Perron–Frobenius operator $\mathcal {F}^{[\Delta t]}$ is similar and one simply keeps track of the transition probabilities from cell $\mathcal {M}_j$ to $\mathcal {M}_i$. A Markovian assumption implies that the holding time is exponentially distributed in the generator case and geometrically distributed for the discrete in time case. We re-examine the Markovian assumption when applying the method in § 4.

In the subsections, we outline an empirical construction of the matrix from finite data and a Bayesian approach incorporating uncertainty due to finite sampling effects. With the latter approach, we do not treat the entries of the $Q_{ij}$ matrix as deterministic numbers but as distributions. The result is a random matrix representation of the generator that incorporates uncertainty. We emphasize that our focus on using the generator of the process is critical to the incorporation uncertainty since we assume that the data comes from a continuous-time dynamical system.

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. 2013). 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. (2013) 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 (2005) and Trendelkamp-Schroer et al. (2015) for other Bayesian examples of calculating uncertainties and Scherer et al. (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 (2023), 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. Illustration of the methodology with the Lorenz equations

We apply the methodology from the previous section to the Lorenz equations. The dynamics are given by

(4.1)$$\begin{gather} \dot{x} ={-} \sigma ( x - y), \end{gather}$$

(4.2)$$\begin{gather}\dot{y} ={-} y + (r-z) x, \end{gather}$$

(4.3)$$\begin{gather}\dot{z} ={-} b z + xy, \end{gather}$$

where we identify $x = s_1$, $y = s_2$, $z = s_3$. The corresponding Liouville equation is given by

(4.4)\begin{align} \partial_t \mathcal{P} + \partial_{\mathscr{x}} \left( \left[ - \sigma ( \mathscr{x} - \mathscr{y})\right] \mathcal{P} \right) + \partial_{\mathscr{y}} \left( \left[ - \mathscr{y} + (r-\mathscr{z}) \mathscr{x} \right] \mathcal{P} \right) + \partial_{\mathscr{z}}\left( \left[ - b \mathscr{z} + \mathscr{x}\mathscr{y} \right] \mathcal{P} \right) = 0 \end{align}

where we use the notation $\mathscr {x}$ for $\mathscr {s}_1$, $\mathscr {y}$ for $\mathscr {s}_2$ and $\mathscr {z}$ for $\mathscr {s}_3$. We shall examine three different methods of partitioning the chaotic attractor.

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ć 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. (2023); 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 (2003), 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 (2023) 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).

5. Conclusion

In summary, we have done three things:

1. (i) § 2 – reviewed a theoretical formulation for transforming a dynamical system into a continuous time Markov process with finite state space;

2. (ii) § 3 – developed a Bayesian stream-based data-driven algorithm for constructing the generator of a continuous-time Markov process with finite state space;

3. (iii) § 4 – applied the methodology to the Lorenz equations.

We have seen that statistics are captured, even with a coarse discretization. In the Lorenz case, we used the fixed points of the dynamical system as the Markov states and the anchors for the partitioning strategy. Three states sufficed to capture the first and second moments. Furthermore, even some autocorrelations and residency times were well-captured with the coarse discretization, depending on the time scale of interest. Furthermore, we explored two further partitioning strategies (using a periodic orbit and random points of the attractor) and reported their resulting emulated statistics.

Future extensions of the present work include a detailed examination of convergence properties by varying the number of states and choosing different classifiers. When the number of cells becomes large, calculating the minimum distance of a state to a cell centre is computationally expensive. In such a case, using a tree structure for the classifier (e.g. a hierarchical k-means) is one option. Similar computationally expedient extensions include the use of a tensor product basis such as Junge & Koltai (2009) or box-refinement strategies used in Dellnitz & Junge (1999) and Dellnitz et al. (2001). Another option is to borrow ideas from thermodynamics and divide state space into ‘macrostates’ and ‘microstates’. For example, we first divide a cell into energy shells and then only compute distances within each energy shell. A noteworthy example of a physics-based partition is found in Jiménez (2023), where the author constructed a Perron–Frobenius operator for wall-bounded flows.

The present work suggests the feasibility of extending this approach to more complex and high-dimensional systems. The only necessary step is to define a method for classifying states. Part 2 of this series focuses on applying the methodology to the compressible Euler equations with gravity and rotation (the model in Souza et al. (2023a)), a high-dimensional dynamical system exhibiting geophysical turbulence. We will detail the particular choices for the classifier and Markov states as applied to the Euler equations. The investigation will include an analysis of different strategies for partitioning the high-dimensional state space, a feature of the methodology that lends itself to adaptation to various dynamical systems. Ultimately, the hope is that the present method enables insights into high-dimensional spaces.