1 Introduction↩︎

Dynamical systems are a vital tool to describe deterministic and stochastic processes in science and engineering: the motion of celestial bodies, dynamics of molecules, or the development of the human brain. Even the most complex dynamical systems can often be analysed by studying certain associated linear operators such as the Koopman operator, the Perron–Frobenius operator, and their generators [1]–[4]. These operators have been used in a wide range of fields, such as molecular dynamics [5], [6], fluid dynamics [7], [8], and engineering [9], [10]. Koopman and Perron–Frobenius operators allow us to study the evolution of observables, such as the velocity, acceleration, or energy of the system, and probability densities. The generators of these operators are used to study the temporal rate of change of observables and densities, respectively.

As a result, the investigation of our dynamical system reduces to having access to a linear operator \(\mathcal{A}\), which acts on observables or probability densities of our system. The goal of data-driven methods is to obtain an approximation \(\widehat{\mathcal{A}}_{NM}\) of \(\mathcal{A}\) using the information derived from studying the system using only a finite amount of observables \(\psi_1,\dots,\psi_N\) and training points \(\boldsymbol{x}_1, \dots, \boldsymbol{x}_M\). These methods gained considerable interest in the literature in recent years, and various methods have been developed to address this problem; some of the most notable are dynamic mode decomposition (DMD) [11]–[14], which relies on a least-squares estimate of the Koopman operator using a set of linear basis functions, extended dynamic mode decomposition EDMD [15], [16], which can be regarded as a nonlinear generalisation of DMD, and generator extended dynamic mode decomposition (gEDMD) [17], [18], which approximates the infinitesimal generator of the Koopman operator.

More recently, the convergence of these methods has been analysed. Convergence properties of EDMD were studied in [19], though no error bounds were given. The rate of convergence of gEDMD has been studied in [20]–[22]. However, all these works require strong assumptions that may be impossible to verify in practice. In this work, our goal is to resolve these limitations. Our main contributions are as follows:

1. A common framework: We introduce a unified framework to study methods for the data-driven analysis of dynamical systems (such as, but not limited to, EDMD and gEDMD). Within this framework, a linear operator \(\mathcal{A}\) is estimated using a Monte Carlo approximation, denoted by \(\widehat{\mathcal{A}}_{NM}\). This approximation utilises \(M\) Monte Carlo samples to estimate \(\mathcal{A}_N\), which is the projection of \(\mathcal{A}\) onto the dictionary space \(\mathcal{F}_N = \mathop{\mathrm{span}}(\psi_1,\ldots, \psi_N)\).

2. Convergence: We show that \(\widehat{\mathcal{A}}_{NM}\) is the projection of \(\mathcal{A}\) onto the space of empirical samples \(\widehat{\mathcal{F}}_{NM}\), almost sure convergence of \(\widehat{\mathcal{A}}_{NM}\) to \(\mathcal{A}_N\), convergence of eigenvalues and weak convergence of eigenfunctions of \(\widehat{\mathcal{A}}_{NM}\) along a subsequence.

3. Error bounds: We derive explicit bounds for the approximation errors \(\|\widehat{\mathcal{A}}_{NM}-\mathcal{A}_N\|\) and \(\|\widehat{\mathcal{A}}_{NM}-\mathcal{A}\|\) and then extend these results to the case of noisy observations.

4. Relaxed assumptions: We derive our results without the restrictive conditions prevalent in prior literature. Specifically, we do not assume the data is sampled from an invariant measure or a single trajectory, nor do we require an orthonormal dictionary, bounded operators on \(\mathcal{F}\), or an invertible empirical Gram matrix. The precise conditions are detailed at the end of Section 2.

The outline of the article is as follows: In Section 2, we introduce Koopman and Perron–Frobenius operators and the mathematical setting for our problem. In Sections 3, 4 and 6, we prove Contribution [part322]. In Section 5, we prove Contribution [part324]. In Section 7, we provide numerical simulations to illustrate our results. We conclude with a discussion and some recommendations based on the theoretical and numerical results in Section 8. Our approach in Sections 3 and 4 is heavily influenced by the results in [19], and our work in Section 5 is more closely related to that of [22], [23]. However, we develop these results under a more general framework in which an arbitrary operator is estimated and under weaker assumptions.

2 A general framework for data-based recovery of dynamics↩︎

In this section, we will introduce the required mathematical concepts and the notation used throughout the paper.

2.1 Notation↩︎

Let \(E_1,E_2\) be generic vector spaces and vectors \(\Psi_1=\left\{\psi_n\right\}_{n=1}^{N_1} \subseteq E_1, \Psi_2=\left\{\phi_n\right\}_{n=1}^{N_2} \subseteq E_2\). We define \(\mathop{\mathrm{span}}(\Psi_1),\mathop{\mathrm{span}}(\Psi_2)\) to be the smallest vector space containing \(\Psi_1,\Psi_2\) respectively. An operator \(\mathcal{T}\colon \mathop{\mathrm{span}}(\Psi_1) \to \mathop{\mathrm{span}}(\Psi_2)\) can be represented by any matrix \(\boldsymbol{T}^{\Psi_1 \to \Psi_2} \in \mathbb{R}^{N_2 \times N_1}\) verifying \[\mathcal{T}\psi_j = \sum_{i=1}^{N_2}\boldsymbol{T}^{\Psi_1 \to \Psi _2}_{ij}\phi_i, \quad j=1,\dots, N_1.\] The matrix \(\boldsymbol{T}^{\Psi_1 \to \Psi _2}\) always exists and if \(\Psi_1\) and \(\Psi _2\) are a collection of independent vectors, \(\boldsymbol{T}^{\Psi_1 \to \Psi _2}\) is unique. If \(\Psi _1=\Psi _2\) we write \(\boldsymbol{T}^{\Psi_1}:=\boldsymbol{T}^{\Psi_1 \to \Psi _1}\).

We denote the Euclidean norm of a vector \(\boldsymbol{v}\in\mathbb{C}^N\) by \(\left| \boldsymbol{v} \right|\). Given an operator \(\mathcal{T}\colon\mathcal{D}\to\mathcal{F}\) between normed spaces and a matrix \(\boldsymbol{T} \in \mathbb{C}^{N\times N}\), we denote the induced operator norms of \(\mathcal{T}\) and \(\boldsymbol{T}\) and the Frobenius norm of \(\boldsymbol{T}\) by \[\|\mathcal{T}\| := \sup_{\left\lVert \phi \right\rVert_{\mathcal{D}}=1}\left\lVert \mathcal{T}\phi \right\rVert_{\mathcal{F}}, \qquad \|\boldsymbol{T}\| := \sup_{\left| \boldsymbol{v} \right|=1}\left| \boldsymbol{T} \boldsymbol{v} \right|, \qquad \left\lVert \boldsymbol{T} \right\rVert_F := \left(\sum_{i,j=1}^N\left| \boldsymbol{T}_{ij} \right|^2\right)^{\frac{1}{2}},\] respectively.

We define random objects, such as random variables, stochastic processes, and random flows on a common underlying probability space \((\Omega, \mathfrak{A}, \mathbb{P})\), but usually simplify our presentation by ignoring this dependence. Moreover, we let \((\mathbb{X}, \mathfrak{B}(\mathbb{X}), \mu)\) be an additional probability space with \(\mathfrak{B}(\mathbb{X})\) being the Borel \(\sigma\)-algebra on the topological space \(\mathbb{X}\) and denote the space of \(\mu\) square-integrable functions on \(\mathbb{X}\) by \(\mathcal{F}:= L^2(\mathbb{X}\to\mathbb{C}, \mu)\).

In what follows, we will denote vectors by bold lowercase letters, matrices by bold uppercase letters, and operators by calligraphic letters. Moreover, \(\boldsymbol{x}\) are elements of \(\mathbb{X}\) and \(\psi_1, \dots, \psi_N\) denote dictionary functions. Finally, we will always use \(i,j,n \in \left\{1, \dots,N\right\}\) to index the dictionary and \(m \in \left\{1, \dots,M\right\}\) to index samples where \(M,N \in \mathbb{N}\) . For the sake of convenience, we include an overview of the notation used throughout the paper in Appendix 11.

2.2 Dynamical systems and linear operators↩︎

Linear operators are a powerful tool for studying dynamical systems. Two such operators are the Koopman and Perron–Frobenius operators, along with their generators. Let \(\mathbb{X}\) be the state space and consider a deterministic discrete-time dynamical system \(\Phi \colon \mathbb{X} \to \mathbb{X}\) defined by \[\label{discrete} \boldsymbol{x}_{\ell+1} = \Phi(\boldsymbol{x}_\ell), \quad \ell \in \mathbb{N}_0,\tag{1}\] with an appropriate initial condition \(\boldsymbol{x}_0 \in \mathbb{X}\). We refer to \(\Phi\) as the flow of the dynamical system. In practice, we may not have access to \(\Phi\). In order to recover some information about the dynamical system, we may measure some quantity \(f\colon\mathbb{X}\to\mathbb{C}\) of the system, such as its velocity, acceleration, or energy and investigate how it evolves. That is, we study \[\label{Koopman32operator} \mathcal{K}f(\boldsymbol{x}) := f(\Phi(\boldsymbol{x})),\tag{2}\] where the operator \(\mathcal{K}\colon C(\mathbb{X}) \to C(\mathbb{X})\) is known as the Koopman32operator. Given a measure \(\mu\) on \(\mathbb{X}\), if \(\Phi\) preserves \(\mu\)-null sets, one can extend the Koopman operator to \(\mathcal{K}\colon L^\infty(\mathbb{X}) \to L^\infty(\mathbb{X})\). Its preadjoint operator, acting on \(L^1(\mathbb{X})\) and often denoted \(\mathcal{K}_*\), is the Perron–Frobenius operator.

These ideas translate to continuous-time dynamical systems, such as \[\boldsymbol{x}_t = \Phi^{t-s}(\boldsymbol{x}_s) \in \mathbb{X}, \quad s,t\in\mathbb{R}^+, s < t,\] where now \(\Phi^t \colon \mathbb{X} \to \mathbb{X}\) defines the flow in continuous time, and we study the time evolution of observables \(f\) through the semigroup of Koopman operators \(\mathcal{K}^t\), \(t \geq 0\), defined by \[\mathcal{K}^t f(\boldsymbol{x}) := f(\Phi^t(\boldsymbol{x})).\] The infinitesimal generator of the Koopman semigroup \(\mathcal{K}^t\) is then defined on suitable \(f\) as \[\mathcal{L} f(\boldsymbol{x}) := \lim _{t \downarrow 0} \frac{1}{t}\left(\mathcal{K}^t f(\boldsymbol{x})-f(\boldsymbol{x})\right).\] The above operators can also be extended to stochastic dynamical systems, i.e., to systems in which the flow \(\Phi\) or \(\Phi^t\) is a random object mapping from \(\Omega \times \mathbb{X}\) to \(\mathbb{X}\). In this setting, we still refer to the Koopman operators and semigroups as well as the associated Perron–Frobenius operators and their generators by \(\mathcal{K}\), \(\mathcal{K}^t\), \(\mathcal{K}_*\), \((\mathcal{K}^t)_*\), \(\mathcal{L}\), and \(\mathcal{L}^*\), respectively. In this case, the Koopman operators and semigroups are given by \[\begin{align} \label{eq:stochKK} \mathcal{K}f(\boldsymbol{x}) := \mathbb{E}[f(\Phi(\boldsymbol{x}))], \qquad \mathcal{K}^t f(\boldsymbol{x}) := \mathbb{E}[f(\Phi^t(\boldsymbol{x}))]. \end{align}\tag{3}\] The corresponding Perron–Frobenius operators and the generators are then defined in an analogous fashion. A typical case is when \(\boldsymbol{X}_t\in \mathbb{X}{ \subset \mathbb{R}^d}\) is a continuous time dynamical system evolving according to the stochastic differential equation \[\label{SDE} \mathrm{d} \boldsymbol{X}_t = \boldsymbol{b}(\boldsymbol{X}_t) \mathrm{d} t + \boldsymbol{\sigma}(\boldsymbol{X}_t) \, \mathrm{d} \boldsymbol{W}_t, \quad \boldsymbol{X}_0=\boldsymbol{x}.\tag{4}\]

By applying Itô’s formula, one can show that the generator \(\mathcal{L}\) associated with the stochastic differential equation 4 is defined on \(f \in C^2(\mathbb{X})\) as \[\require{physics} \label{Koopman32generator} \mathcal{L} f= \boldsymbol{b}\cdot \nabla f+\frac{1}{2} {\mathrm{Tr}}\qty(\boldsymbol{\Sigma} \nabla^2 f),\tag{5}\] where \(\boldsymbol{\Sigma} := \boldsymbol{\sigma \sigma}^\top\) [24]. In physical terms, \(\mathcal{L}\) describes the infinitesimal rate of change of observables \(f\) evolving under our dynamical system. That is, writing \(u(t, \boldsymbol{x}) := \mathcal{K}^t f(\boldsymbol{x})\), we have \[\partial _t u = \mathcal{L}u.\] The above is called the Kolmogorov backward equation. Similarly, \(\mathcal{L}^*\) describes the evolution of probability distributions of our system. If \(X_t\) has density \(\nu(t, \boldsymbol{x}) \in C^1_tC^2_{\boldsymbol{x}}([0, T]\times \mathbb{X})\), i.e., is once differentiable in time and twice in space with bounded derivatives, it holds that \[\partial _t \nu = \mathcal{L}^* \nu,\] see [25]. This is known as the forward Kolmogorov equation or Fokker–Planck equation. As a result, knowledge of \(\mathcal{L}\) gives us complete knowledge of the evolution of \(\boldsymbol{X}\).

Through the expressions above, we see that the Koopman and Perron–Frobenius operators, along with their generators, are critical tools to describe how a dynamical system evolves.

2.3 Mathematical framework↩︎

We now move on to the data-driven approximation of an operator of interest, which we denote by \(\mathcal{A}\). In what follows, this operator can be one of the operators we have introduced in the context of dynamical systems, i.e., \(\mathcal{K}\) or \(\mathcal{K}_*\) in discrete time and \(\mathcal{K}^t\), \((\mathcal{K}^t)_*\), \(\mathcal{L}\), or \(\mathcal{L}^*\) in continuous time, in which case it is necessary that the flow \(\Phi\) preserves \(\mu\)-null sets. However, our theory is not limited to these operators. As before, we consider the state space \((\mathbb{X}, \mathcal{B}, \mu)\) and write \(\mathcal{F}:= L^2(\mathbb{X}, \mu)\). The observables on which \(\mathcal{A}\) can be evaluated define the domain of \(\mathcal{A}\), i.e., \[\mathcal{D}:= \left\{ f\in \mathcal{F}: \mathcal{A}f \in \mathcal{F}\right\}.\] We say that \(\mathcal{A}\) is a closed operator if the graph of \(\mathcal{A}\), \[\{(f, \mathcal{A}f) \in \mathcal{F}\times \mathcal{F}: f \in \mathcal{D}\},\] is a closed subspace of \(\mathcal{F}\times\mathcal{F}\). If \(\mathcal{A}\) is a closed operator, \(\mathcal{D}\) is a Hilbert space with the inner product \[\left\langle f,g\right\rangle_{\mathcal{D}} := \left\langle f,g\right\rangle_\mathcal{F}+\left\langle\mathcal{A}f, \mathcal{A}g\right\rangle_\mathcal{F}, \quad \left\lVert f \right\rVert_\mathcal{D}^2 := \left\langle f,f\right\rangle_{\mathcal{D}}, \quad \forall f,g \in \mathcal{D},\] see [26] for more details. A direct consequence of the above definitions is that \(\left\lVert f \right\rVert_{\mathcal{F}}\leq \left\lVert f \right\rVert_{\mathcal{D}}\) and \(\mathcal{A}\colon\mathcal{D}\to\mathcal{F}\) is a continuous operator with \(\left\lVert \mathcal{A} \right\rVert\leq 1\).

In practice, there are two main cases. First, if \(\mathcal{A}\) is a bounded operator on \(\mathcal{F}\) (e.g., the Koopman operator \(\mathcal{K}^t\) provided \(\Phi\) preserves \(\mu\)-null sets), then it is defined on all of \(\mathcal{F}\) and closed, so \(\mathcal{D}=\mathcal{F}\), and the norms \(\left\lVert \cdot \right\rVert_\mathcal{D}\) and \(\left\lVert \cdot \right\rVert_\mathcal{F}\) are equivalent. Second, \(\mathcal{A}\) may be an unbounded operator, such as the infinitesimal generator \(\mathcal{L}\) of a semigroup (such as in 5 ). For \(\mathcal{L}\) to be a closed operator, it is sufficient that the associated semigroup \(\mathcal{K}^t\) is strongly continuous; standard semigroup theory dictates that the generator of any strongly continuous semigroup is densely defined and closed. In this case, \(\mathcal{D}\) is a dense proper subspace of \(\mathcal{F}\), see [28] for details.

A common example is where \(\mathbb{X}\) is an open subset of \(\mathbb{R}^d\) and \(\mathcal{D}\) is a weighted Sobolev space \(H^k(\mathbb{X}, \mu )\) supplemented by appropriate boundary conditions [17]. Specifically, if \(\mu\) is absolutely continuous with respect to the Lebesgue measure with a locally integrable density, the weighted Sobolev space \(H^k(\mathbb{X},\mu)\) is defined as the space of functions in \(L^2(\mathbb{X},\mu)\) whose classical weak derivatives \(D^\alpha\) up to order \(k\) remain in \(L^2(\mathbb{X},\mu)\), equipped with the norm [29], [30] \[\left\|\psi \right\|_{H^k(\mathbb{X}, \mu )}=\sum_{|\alpha|\leq k}\left(\int_{\mathbb{X}}\left|D^\alpha \psi\right|^2 \,\mathrm{d}\mu \right)^{\frac{1}{2}}.\] For example, ignoring boundary conditions, let \(\mu\) be absolutely continuous with respect to the Lebesgue measure and let \(\mathcal{A}\) be the generator in 5 , then \(\mathcal{D}= H^2(\mathbb{X}, \mu)\). If the system is deterministic, on the other hand, then \(\mathcal{D}=H^1(\mathbb{X}, \mu)\).

We emphasize that this is strictly a concrete example for differential operators on Euclidean domains. For systems with singular measures (e.g., supported on attractors), classical weak derivatives are undefined, and the domain \(\mathcal{D}\) is instead defined abstractly via the generator. Provided the Koopman semigroup is strongly continuous on \(L^2(\mathbb{X}, \mu)\), standard semigroup theory guarantees that this abstract domain \(\mathcal{D}\) is dense in \(L^2(\mathbb{X}, \mu)\), ensuring our general framework applies without requiring absolute continuity.

Our goal is to study a given dynamical system through the action of \(\mathcal{A}\) on observables \(f \in \mathcal{D}\). A practical difficulty is that we will typically not know the value of \(\mathcal{A}f(\boldsymbol{x})\) for every observable \(f\) and every point \(\boldsymbol{x} \in \mathbb{X}\), as this would correspond to having an infinite amount of information about our system. The most we can hope for is to know this information for a finite amount of observables \({\psi_1, \dots, \psi_N}\subset \mathcal{D}\), called the dictionary, and a finite amount of training data points \(\boldsymbol{x}_1, \dots, \boldsymbol{x}_M\), which are assumed to be sampled independently from some arbitrary probability measure \(\mu\). That is, our total information is \[\label{info} \big\{\psi_n(\boldsymbol{x}_m), \mathcal{A}\psi_n (\boldsymbol{x}_m)\big\}_{m,n=1}^{M,N}.\tag{6}\] Here, the terms \(\mathcal{A}\psi_n (\boldsymbol{x}_m)\) may be evaluated exactly or approximated using existing techniques (e.g., finite differences, Kramers–Moyal expansions, etc.), see [17], [31]. Since \(\mathcal{A}\) is infinite-dimensional, this finite amount of information cannot fully describe \(\mathcal{A}\) except in simple cases. In general, the best we can hope is to obtain some approximation \(\widehat{\mathcal{A}}\) of \({\mathcal{A}}\).

Many approaches for the data-driven description of dynamics, such as DMD, EDMD, and gEDMD, can be described as trying to obtain the best approximation \(\widehat{\mathcal{A}}\) of \(\mathcal{A}\) given the limited data in 6 . In fact, all of them can be subsumed under the framework that we will introduce below.

There are two aspects that need to be considered: a finite number of dictionary functions \(\psi_1,\ldots,\psi_N\) that may not span the space \(\mathcal{F}\) and a finite number of samples \(\boldsymbol{x}_1,\ldots,\boldsymbol{x}_M\) that are not space-filling in \(\mathbb{X}\). We now discuss the implications of these aspects, starting with the finite dictionary. Using the information in 6 , the operator \(\mathcal{A}\) can be approximated on \[\mathcal{F}_N := \mathop{\mathrm{span}}(\Psi) = \mathop{\mathrm{span}}(\psi_1, \dots, \psi_N)\] by taking sampling-based approximations of the projection of \(\mathcal{A}\) onto \(\mathcal{F}_N\). In all the preceding methods, this is done as follows: We would like our finite-dimensional approximation \(\mathcal{A}_N\colon \mathcal{F}_N \to \mathcal{F}_N\) to \(\mathcal{A}\) to satisfy \[\require{physics} \label{Galerkin32matrices} \big[\boldsymbol{C}_N\big]_{ij} :=\left\langle\mathcal{A}\psi_i, \psi_j\right\rangle_{L^2(\mu )}=\left\langle\mathcal{A}_N\psi_i, \psi_j\right\rangle_{L^2(\mu )}= \qty[\qty(\boldsymbol{A}^\Psi_N)^\top \boldsymbol{G}_N]_{ij},\tag{7}\] for all \(i,j\), where \[\big[\boldsymbol{C}_N\big]_{ij} :=\left\langle\mathcal{A}\psi_i, \psi_j\right\rangle_{L^2(\mu )},\quad \big[\boldsymbol{G}_N\big]_{ij} := \left\langle\psi_i, \psi_j\right\rangle_{L^2(\mu )}\] are the structure matrix and the Gram matrix of \(\mathcal{A}\) and \(\Psi\), respectively. If we could evaluate the above integrals exactly, we would obtain our approximation via \[\label{exact32Galerkin} \big(\boldsymbol{A}_N^\Psi\big)^\top = \boldsymbol{C}_N\boldsymbol{G}_N^{-1},\tag{8}\] where the invertibility of the Gram matrix is equivalent to \(\Psi\) being linearly independent. The operator \(\mathcal{A}_N\colon\mathcal{F}_N\to\mathcal{F}_N\) defined through 8 satisfies 7 and, as a result, is necessarily the projection of \(\mathcal{A}\) onto \(\mathcal{F}_N\). That is, if we write \(\mathcal{P}_{\mathcal{F}_N}\) for the projection onto \(\mathcal{F}_N\), \[\label{Galerkin32operator} \mathcal{A}_N=\left.\mathcal{P}_{\mathcal{F}_N} \mathcal{A}\right|_{\mathcal{F}_N} .\tag{9}\] Due to the connection with the finite element method (FEM), 9 is often called the Galerkin approximation (or projection) of \(\mathcal{A}\). This is also known in some contexts as the finite section of \(\mathcal{A}\).

However, since we only have access to the information in 6 , the best we can do is to use the samples \(\boldsymbol{x}_1,\ldots,\boldsymbol{x}_M \sim \mu\) and define the empirical structure matrix and the empirical Gram matrix as \[\label{approximate32matrices} \big[\widehat{\boldsymbol{C}}_{NM}\big]_{ij} := \frac{1}{M} \sum_{m=1}^M \overline{\psi_j(\boldsymbol{x}_m)}\mathcal{A}\psi_i (\boldsymbol{x}_m), \quad \big[\widehat{\boldsymbol{G}}_{NM}\big]_{ij} := \frac{1}{M}\sum_{m=1}^{M} \psi_i(\boldsymbol{x}_m) \overline{\psi_j(\boldsymbol{x}_m)},\tag{10}\] respectively. Let us denote the empirical measure associated with \(\{\boldsymbol{x}_1,\ldots,\boldsymbol{x}_M\}\) by \[\widehat{\mu}_M := \frac{1}{M} \sum_{m=1}^M \delta _{\boldsymbol{x}_m},\] then the above can also be written as \[\big[\widehat{\boldsymbol{C}}_{NM}\big]_{ij}= \left\langle\mathcal{A}\psi_i, \psi_j\right\rangle_{L^2(\widehat{ \mu }_M)} , \quad \big[\widehat{\boldsymbol{G}}_{NM}\big]_{ij}= \left\langle\psi_i, \psi_j\right\rangle_{L^2(\widehat{ \mu }_M)} .\] We define, analogously to 8 , \[\label{approximate32Galerkin} \widehat{\boldsymbol{A}}_{NM}^\top := \widehat{\boldsymbol{C}}_{NM}\widehat{\boldsymbol{G}}_{NM}^+.\tag{11}\] Here, \(\widehat{\boldsymbol{G}}_{NM}^+\) denotes the Moore–Penrose pseudoinverse [32] of \(\widehat{\boldsymbol{G}}_{NM}\) and by construction is such that 11 minimises the empirical error, \[\left\lVert \widehat{\boldsymbol{A}}_{NM}^\top\widehat{\boldsymbol{G}}_{NM} - \widehat{\boldsymbol{C}}_{NM} \right\rVert^2_F := \sum_{i,j=1}^N \left| \big[\widehat{\boldsymbol{A}}_{NM}^\top\widehat{\boldsymbol{G}}_{NM}\big]_{ij} - \big[\widehat{\boldsymbol{C}}_{NM}\big]_{ij} \right|^2.\] The Gram matrix \(\boldsymbol{G}_N\) is always invertible as the basis functions contained in \(\Psi\) are assumed to be linearly independent. However, its Monte Carlo approximation \(\widehat{\boldsymbol{G}}_{NM}\) may not be invertible. For this reason, the pseudoinverse is required and ensures that we obtain the matrix \(\widehat{\boldsymbol{A}}_{NM} \in \mathbb{C}^{N\times N}\). This can give rise to theoretical issues, such as the discontinuity of the map from a matrix to its pseudoinverse. These problems can be resolved with careful analysis. We will touch upon this in more detail in Theorem 1.

By the strong law of large numbers and the definitions in 7 and 10 , almost surely \[\lim_{M \to \infty}\widehat{\boldsymbol{C}}_{NM} = \boldsymbol{C}_{N}, \qquad \lim_{M \to \infty}\widehat{\boldsymbol{G}}_{NM} = \boldsymbol{G}_N .\] As a result, we expect \(\widehat{\boldsymbol{A}}_{NM}\) to converge to \(\boldsymbol{A}^\Psi_{N}\) as the number of data points goes to infinity. This will be studied in Section 3. Convergence to \(\mathcal{A}\), when the number of basis functions goes to infinity, is studied in Sections 4 and 5, where also error bounds are established. Convergence to the spectrum of \(\mathcal{A}\) is studied in Section 6. In practice, one may not have access to the exact values of \(\psi_n(\boldsymbol{x}_m),\mathcal{A}\psi_n(\boldsymbol{x}_m)\), but only to an approximation or a noisy measurement. To deal with this case, we also study in Subsection [Noise32section] the case where our data has some noise \(\boldsymbol{\eta}, \boldsymbol{\xi}\), i.e., we have access to the mappings \[{\boldsymbol{x}}_m \to \psi_n(\boldsymbol{x}_m)+\eta_N^{m,n}, \quad {\boldsymbol{x}}_m \to \mathcal{A}\psi_n(\boldsymbol{x}_m)+\xi_N^{m,n}, \quad n=1,\dots,N.\] We show how a direct implementation in this case introduces some bias, whereas by splitting the evaluations of the basis functions into two batches, convergence of the approximations to the true dynamics is obtained.

We stress once more that all the previously mentioned methods for data-driven recovery of dynamics (DMD, EDMD, and gEDMD) fall into this framework and are thus covered by the analysis. Here, we use minimal assumptions. More precisely, we require that the observables be linearly independent, sufficiently regular and that the data be i.i.d.(Assumptions 1 and 2). To show convergence and error bounds, it is also necessary that, as \(N\) goes to infinity, the basis functions \(\psi_1,\dots,\psi_N\) cover the whole space (Assumption 3 and 4) and are bounded (Assumption 5). Finally, for the case where the measurements are noisy, we assume that the noise is centred at zero and independent (Assumption 7).

3 Data-driven approximation as a projection↩︎

The following section is inspired by [19]. However, in our analysis we do not require the empirical Gram matrix \(\widehat{\boldsymbol{G}}_{NM}\) to be invertible. We discuss natural situations in which \(\widehat{\boldsymbol{G}}_{NM}\) may not be invertible throughout this section. We begin by imposing the basic assumptions that will be used in what follows.

Point \(\ref{a1}\) of Assumption 1 is necessary to ensure that the Gram matrix \(\boldsymbol{G}_N\) is invertible. Point [continuous] is required so that the pointwise evaluation in the Monte Carlo approximations 10 is well-defined. Lastly, Point [a3] is the basic assumption underlying Monte Carlo approximations.

As we will see in Theorem 1, the approximation \(\widehat{\boldsymbol{A}}_{NM}\) is best defined on the empirical spaces \[\widehat{\mathcal{F}}_{M} := L^2(\mathbb{X}\to\mathbb{C},\widehat{\mu}_M), \quad {\widehat{\mathcal{D}}_M:= \widehat{ \mathcal{F}} _M \times \widehat{\mathcal{F}} _M}.\] Elements in \(\mathcal{F}\) are not in \(\widehat{\mathcal{F}}_M\) as functions defined \(\mu\) almost everywhere are not in general well-defined \(\widehat{\mu }_M\) almost everywhere. However, if \(\phi \in \mathcal{F}\) is continuous almost everywhere, we can view \(\phi\) as an element of \(\widehat{\mathcal{F}}_M\) through the (non-injective) mapping \[\label{identification320} \phi \mapsto \phi^{{\widehat{\mathcal{F}}}} := \sum_{m=1}^M \phi(\boldsymbol{x}_m)\delta_{\boldsymbol{x}_m}\in \widehat{\mathcal{F}}_M.\tag{12}\] Likewise, given \(\psi \in \mathcal{D}\) such that \(\mathcal{A}\psi\) is continuous almost everywhere, we can view \(\psi\) as an element of \({\widehat{\mathcal{D}}}_M\) through \[\require{physics} \label{identification} \psi \mapsto \psi^{{\widehat{\mathcal{D}}}} := \sum_{m=1}^M \qty(\psi(\boldsymbol{x}_m), {\mathcal{A}\psi ({\boldsymbol{x}}_m)})\delta_{\boldsymbol{x}_m}\in {\widehat{\mathcal{D}}}_M.\tag{13}\] We also use the notation \[\begin{align} \Psi^{{\widehat{\mathcal{F}}} } &:= \{\psi_n^{{\widehat{\mathcal{F}}}}\}_{n=1}^N \subset \widehat{\mathcal{F}}_M, \quad {\widehat{\mathcal{F}}}_{NM} := \mathop{\mathrm{span}}(\Psi^{{\widehat{\mathcal{F}}}})\subset \widehat{\mathcal{F}}_{M}, \\ {\Psi^{\widehat{\mathcal{D}}} } &:= \{\psi_n^{\widehat{\mathcal{D}}}\}_{n=1}^N \subset \widehat{\mathcal{D}}_M, \quad \widehat{\mathcal{D}}_{NM} := \mathop{\mathrm{span}}(\Psi^{\widehat{\mathcal{D}}})\subset \widehat{\mathcal{D}}_{M} \end{align}\] to denote the basis \(\Psi\) and subspace \(\mathcal{F}_N\) when viewed as objects in \(\widehat{\mathcal{F}}_M\) and \(\widehat{\mathcal{D}}_M\) using 12 and 13 , respectively. We denote by \(\mathcal{A}^{{\mathrm{emp}}}\) the operator induced by \(\mathcal{A}\), i.e., \[\require{physics} \label{induced32operator} {\mathcal{A}^{{\mathrm{emp}}} \colon \widehat{\mathcal{D}}_M \to \widehat{\mathcal{F}}_M}, \quad \mathcal{A}^{{\mathrm{emp}}}\psi^{{\widehat{\mathcal{D}}}} := \qty(\mathcal{A}\psi)^{{\widehat{\mathcal{F}}}}.\tag{14}\] There may be various functions \(\psi \in \mathcal{D}\) with the same empirical representative \(\psi ^{\widehat{\mathcal{D}} }\). However, by construction of \({\cdot }^{\widehat{\mathcal{D}}}\), all of them map to the same value \(\require{physics} \qty(\mathcal{A}\psi)^{\widehat{\mathcal{F}}}\). As a result, \(\mathcal{A}^{{\mathrm{emp}}}\) is well defined.

The functions \(\psi^{{\widehat{\mathcal{D}}} }\) and \((\mathcal{A}\psi)^{{\widehat{\mathcal{F}}} }\) are only well-defined if \(\psi\) and \(\mathcal{A}\psi\) are continuous almost everywhere. In our case, this is satisfied for all \(\psi\in \mathcal{F}_N\) by Point [continuous] of Assumption 1. Using this notation, the Monte Carlo approximation of the structure and Gram matrices in 10 become \[\label{approximate32matrices32prod} \big[\widehat{\boldsymbol{C}}_{NM}\big]_{ij}=\left\langle\mathcal{A}^{{\mathrm{emp}}} \psi^{{\widehat{\mathcal{D}}} }_i, \psi_j^{{\widehat{\mathcal{F}}} }\right\rangle_{L^2(\widehat{\mu }_M)}, \quad \big[\widehat{\boldsymbol{G}}_{NM}\big]_{ij}=\left\langle\psi_i^{{\widehat{\mathcal{F}}} }, \psi_j^{{\widehat{\mathcal{F}}} }\right\rangle_{L^2(\widehat{\mu }_M)}.\tag{15}\] Here, the Gram matrix \(\widehat{\boldsymbol{G}}_{NM}\) with respect to \(\Psi^{{\widehat{\mathcal{F}}} }\) is in general not invertible as \(\Psi^{{\widehat{\mathcal{F}}} }\) may no longer be linearly independent (consider for example the case \(M=1\) and \(N=2\)). By construction, \(\widehat{\mathcal{F}}_M\) is a Hilbert space, so we can define the projection of \(\widehat{\mathcal{F}}_M\) onto \(\widehat{\mathcal{F}}_{NM}\) by \[\mathcal{P}_{\widehat{\mathcal{F}}_{NM}}\colon\widehat{\mathcal{F}}_M\to \widehat{\mathcal{F}}_{NM}.\] We prove that our data-driven approximation \(\widehat{\boldsymbol{A}}_{NM}\) corresponds to the projection of \(\mathcal{A}^{{\mathrm{emp}}}\) onto \(\widehat{\mathcal{F}}_{NM}\).

We now show that if \(\mathcal{F}_N\) is invariant under the action of \(\mathcal{A}\), then \(\widehat{\mathcal{A}}_{NM}^{{\mathrm{emp}}}\) is an exact approximation of the Galerkin projection \(\mathcal{A}_N\) defined in 9 of \(\mathcal{A}\) with probability \(1\).

As in 14 , we define \[\begin{align} \mathcal{A}_N^{{\mathrm{emp}}} : \widehat{\mathcal{D}}_{NM} \to \widehat{\mathcal{F}}_{NM}, \quad \psi^{\widehat{\mathcal{D}} } \to (\mathcal{A}_N \psi)^{\widehat{\mathcal{F}} }. \end{align}\]

Theorem 1 is a generalisation of Theorem 1 in [19], where it was also required that \(\Psi^{{\widehat{\mathcal{F}} }}\) be linearly independent in \(\widehat{\mathcal{F}}_M\) (or equivalently, that \(\widehat{\boldsymbol{G}}_{NM}\) be invertible). However, this assumption is not altogether benign. The assumption will always fail when \(M< N\) (as \(\widehat{\mathcal{F}}_M\) has dimension \(M\)) and also in many cases of interest, such as the following example.

In the above example, if there is a cell \(T_{l_0}\) which holds no data point, then \(\widehat{\psi}_{l_0}=0\). As a result, \(\widehat{\Psi}\) will not be linearly independent in general. We now study the convergence of \(\widehat{\mathcal{A}}_{NM}\) to the projection of \(\mathcal{A}\) onto \(\mathcal{F}_N\). This convergence uses the law of large numbers and requires that the empirical matrices have finite variances.

For example, the above assumption will hold by Assumption 1 if \(\mathbb{X}\) is a compact domain. Though \(\widehat{\Psi}\) are not linearly independent in general, they will be (almost surely) when we have a large enough training data set. This is shown in the following lemma.

Theorem 3.2 establishes that \(\widehat{\boldsymbol{A}}_{NM}\) represents a projection operator, \(\widehat{\mathcal{A}}_{NM}^{\mathrm{emp}}\), in the discrete space of the data samples. This guarantees our approximation is optimal with respect to the given data.

For convergence and error analysis, however, we must compare our approximation to the true operator \(\mathcal{A}_N\), which acts on the continuous space \(\mathcal{F}_N\). We therefore use the empirically derived matrix to define our primary object of analysis.

This definition is our object of interest for the remainder of this paper and is valid for all \(M\) and \(N\). We now connect it back to the empirical setting.

Using Lemma 1, we now prove that in the infinite sample limit, the data-driven operator \(\widehat{\mathcal{A}}_{NM}\) converges to \(\mathcal{A}_N\), i.e., the Galerkin projection of \(\mathcal{A}\) onto \(\mathcal{F}_N\).

4 Convergence of the projections↩︎

In the previous section, we have shown that the data-driven approximation in 11 defines an operator \(\widehat{\mathcal{A}}_{NM}\) that converges to the Galerkin projection \(\mathcal{A}_N=\left. \mathcal{P}_{\mathcal{F}_N} \mathcal{A}\right|_{\mathcal{F}_N}\). In this section, our goal is to show that \(\mathcal{A}_N\) converges to \(\mathcal{A}\). This was also done in [19], where it was assumed that the (Koopman) operator is bounded on \(\mathcal{F}\) and that \(\left\{\psi_n\right\}_{n=1}^\infty\) form an orthonormal basis of \(\mathcal{F}\). In many cases, however, \(\mathcal{A}\) will not be bounded (for example, if \(\mathcal{A}\) is the generator \(\mathcal{L}\) of the Koopman operator). Furthermore, the requirement that the basis be orthonormal is restrictive, and in practice, it is often preferable to work with a dictionary that does not have to be orthonormalised; a finite element basis, for instance, will often produce sparse operators, but it is not orthonormal. Additionally, the sampling measure \(\mu\) is typically unknown, so it is not possible to orthonormalise with respect to the norm on \(\mathcal{F}\).

We will work in a more general setting. Firstly, we require that \(\mathcal{F}_N\) converges to \(\mathcal{F}\) as \(N \rightarrow \infty\). That is, as we take more basis functions, we fill \(\mathcal{F}\).

We will also need to approximate functions in the domain \(\mathcal{D}\). This necessitates the following assumption.

If \(\mathcal{A}\) is bounded on \(\mathcal{F}\), then \(\mathcal{F}=\mathcal{D}\) and Assumptions 3 and 4 are equivalent. In general, though, this is not the case. One assumption may imply that functions are approximated well in \(L^2(\mathbb{X}, \mu)\) and the other in \(H^r(\mathbb{X}, \mu)\).

To see an example motivated by spectral approximation of the Koopman operator when \(\mu\) is invariant, see [34].

The above assumptions only require that the subspaces \(\mathcal{F}_N\) are good approximations of \(\mathcal{F}\) and \(\mathcal{D}\), with the error vanishing in the limit. We now introduce the notation \[\mathcal{F}_\infty := \bigcup_{n=1}^\infty \mathcal{F}_n.\] Under the above assumptions, we can prove the convergence of the Galerkin approximation.

The first part of Theorem 3 is useful in the case where we have a finite-dimensional space of observables we are interested in. In this case, these can be incorporated directly into \(\mathcal{F}_N\). The second part of the theorem is useful when we want to know the evolution of every possible observable. The proof also shows that the order of convergence depends completely on the rate of convergence of \(\mathcal{P}_{\mathcal{F}_N} f\) and \(\mathcal{P}_{\mathcal{D}_N}f\) to \(f\). We summarise this result in the following corollary.

5 Joint limit in data and dictionary↩︎

In Sections 3 and 4, we studied the iterated limits of \(\widehat{\mathcal{A}}_{NM}\) when the size of the training data set \(M\) and the number of basis functions \(N\) go to infinity. In this section, we study the behaviour of \(\widehat{\mathcal{A}}_{NM}\) when \(M\) and \(N\) increase simultaneously.

Some concentration bounds for EDMD were derived in [35]. The projection error was studied for various approximation spaces, such as reproducing kernel Hilbert spaces and those generated by finite-dimensional bases of wavelets, in [20]. In [21], the authors work with the generator of an ordinary differential equation and derive a projection error in the context of a finite element basis. They also provide a finite-data error bound on the approximation of the generator in the case where the data is sampled from the Lebesgue measure. In [22], the authors derive an error bound for the approximation error of gEDMD under the assumptions that the Koopman semigroup is exponentially stable and the points are sampled from a probability measure invariant under the flow and from a single ergodic trajectory. The error bounds in [21], [22] both require \(M\) to be “sufficiently large” so that the empirical Gram matrix is invertible. However, no bound on how large \(M\) must be is given. In fact, a deterministic bound on \(M\) is impossible, but rather a probabilistic one is needed.

We work under relaxed assumptions, for example, we do not impose that the dictionary functions be orthonormal, we do not impose that \(\boldsymbol{x}_1, \dots, \boldsymbol{x}_M\) be sampled from the Lebesgue measure, a measure invariant under the flow or even from the same trajectory, and we do not impose that the empirical Gram matrix be invertible. We first formulate an existence result of the following type:

This existence result forms a good foundation. In practice, however, we may need an explicit dependence between \(M\) and \(N\) and explicit error bounds. Before we can proceed, we first need to state additional assumptions.

Let us illustrate these assumptions with two examples.

The first ingredient for our error analysis is an error bound for the Gram matrix \(\boldsymbol{G}\).

In an analogous fashion, we can bound the error due to using \(\boldsymbol{\widehat{C}}_{NM}\) with an arbitrarily large probability. In fact, the result is more straightforward as we no longer have to deal with the matrix inversion necessary for \(\boldsymbol{\widehat{G}}_{NM}\).

Having bounded the error due to using \(\widehat{\boldsymbol{G}}_{NM}\) and \(\widehat{\boldsymbol{C}}_{NM}\) instead of \(\boldsymbol{G}_{N}\) and \(\boldsymbol{C}_{N}\), we can now bound the error in approximating \(\mathcal{A}_N\) by \(\widehat{\mathcal{A}}_{NM}\). To do so, we need to relate the bounds in the matrix norms back to the operator norms. As shown in the Lemma 5, this depends on the condition number of \(\boldsymbol{G}_N\), which can be written as \[\kappa (\boldsymbol{G}_N) = \frac{\lambda_{\max}(\boldsymbol{G}_N)}{\lambda_{\min}(\boldsymbol{G}_N)},\] since \(\boldsymbol{G}_N\) is Hermitian. The condition number is \(1\) for an orthonormal basis but can become very large for other bases, such as the basis of monomials. This will be discussed in more detail below.

The bound on the error in approximating the projection \(\mathcal{A}_N\) in the just proved proposition can be combined with an error in the approximation of \(\mathcal{A}\) through \(\mathcal{A}_N\). This gives the following result.

This theorem allows us to estimate the cost of an approximation of \(\mathcal{A}\) with accuracy \(\varepsilon > 0\) that holds with probability \(p\). In terms of the order of convergence of \(\widehat{\mathcal{A}}_{NM}\) to \(\mathcal{A}_N\), we have that \[\require{physics} \label{order32of32convergence} M =\mathcal{O}\qty(\gamma_N^{} \max \left\{\left\lVert \boldsymbol{G}_N \right\rVert,\left\lVert \boldsymbol{T}_N \right\rVert\right\}\kappa \qty(\boldsymbol{G}_N)\left\lVert \boldsymbol{C}_N \right\rVert^2 \left\lVert \boldsymbol{G}_N^{-1} \right\rVert^4 \log \left(\frac{ N}{1-p}\right) \varepsilon ^{-2}).\tag{29}\] An important consequence of Theorem 6 is that it provides an explicit lower bound on \(M\) as a function of \(N\), which allows one to choose \(M\) adaptively depending on \(N\). The ability to prescribe \(M\) as an explicit increasing function of \(N\), rather than merely requiring \(M\) to be “sufficiently large” without quantification, is a rare and highly non-trivial feature in the operator approximation literature. See [36] for a related discussion in the context of Koopman learning.

Let us now consider two examples.

5.1 Accounting for measurement error within the operators↩︎

In this section we show the effect of measurement error in the data-driven approximation of the operator \(\mathcal{A}\). In the first subsection we show how a direct implementation leads to bias. In the second subsection we show how to correct for this bias by dividing the evaluation of the basis functions into two batches.

5.1.1 Data-driven approximation with noisy measurements: The issue of bias↩︎

In practice, it is often not possible to evaluate the dictionary \((\psi_n)_{n=1}^N\) and operator applied to the dictionary functions \((\mathcal{A}\psi_n)_{n=1}^N\) precisely. The evaluations may be perturbed by measurement errors, or we may need to approximate the operator \(\mathcal{A}\), as mentioned in the beginning of Section 2. This error will influence the accuracy of the approximation of \(\mathcal{A}\). We will assume in what follows that instead of 6 , we only have access to \[\label{data22} \{\psi_n(\boldsymbol{x}_m)+\eta_N^{m,n}, \mathcal{A}\psi_n(\boldsymbol{x}_m) +\xi_N^{m,n}\}_{m,n=1}^{M,N},\tag{30}\] where \((\eta_N^{m,1},\dots, \eta_N^{m, N}) =: \boldsymbol{\eta }_N^m\) and \((\xi^{m,1},\dots,\xi^{m, N}) =: \boldsymbol{\xi }_N^m\) represent the measurement or evaluation error in dictionary and operator, respectively. Given this perturbed data, a direct implementation of the data-driven approximation is \[\widetilde{{\boldsymbol{A}}}_{NM}^\top := \widetilde{{\boldsymbol{C}}}_{NM}\widetilde{{\boldsymbol{G}}}_{NM}^+,\] with perturbed structure and Gram matrices given by \[\require{physics} \begin{align} \label{noisy32matrices2} [\widetilde{\boldsymbol{{C}}}_{NM}]_{ij} & :=\frac{1}{M}\sum_{m=1}^M\qty(\mathcal{A}\psi_i(\boldsymbol{x}_m)+\boldsymbol{\xi}_N^{m,i})\qty(\psi_j(\boldsymbol{x}_m)+\boldsymbol{\eta}_N^{m,j})^\dagger, \\ [\widetilde{\boldsymbol{{G}}}_{NM}]_{ij} & :=\frac{1}{M}\sum_{m=1}^M\qty(\psi_i(\boldsymbol{x}_m)+\boldsymbol{\eta}_N^{m,i})\qty(\psi_j(\boldsymbol{x}_m)+\boldsymbol{\eta}_N^{m,j})^\dagger, \end{align}\tag{31}\] respectively. In this section, we show how the error bounds in Lemmas 2 and 3 can be modified to account for this noise. This leads to an analogous error bound on the data-driven operator as in Proposition 5 and Theorem 6 where a bias term appears. We make the following assumptions on the noise.

Given \(\widetilde{\gamma }_N \geq 0\) we will use the notation \[\require{physics} \widetilde{p}_N:=\mathbb{P}\qty[\left| (\boldsymbol{\eta}_N^m,\boldsymbol{\xi}_N^m) \right|^2\leq \widetilde{\gamma}_N^{}].\] Let \(\mathbb{E}[\cdot \mid \cdot]\) denote the conditional expectation and write \[\require{physics} \begin{align} {\boldsymbol{\Sigma}}_{N}^{\boldsymbol{\eta}}&:=\mathbb{E}\qty[\boldsymbol{\eta}_N^m(\boldsymbol{\eta}_N^{m})^\dagger\mid \left\{\left| (\boldsymbol{\eta}_N^m,\boldsymbol{\xi}_N^m) \right|^2\leq \widetilde{\gamma}_N^{}, ~\forall m\right\} ],\\ {\boldsymbol{\Sigma}}_{N}^{\boldsymbol{\xi}}&:=\mathbb{E}\qty[\boldsymbol{\xi}_N^m(\boldsymbol{\xi}_N^{m})^\dagger\mid \left\{\left| (\boldsymbol{\eta}_N^m,\boldsymbol{\xi}_N^m) \right|^2\leq \widetilde{\gamma}_N^{}, ~\forall m\right\}], \\ {\boldsymbol{\Sigma}}_{N}^{{\boldsymbol{\xi}},{\boldsymbol{\eta}}}&:=\mathbb{E}\qty[\boldsymbol{\xi}_N^m(\boldsymbol{\eta}_N^{m})^\dagger\mid\left\{\left| (\boldsymbol{\eta}_N^m,\boldsymbol{\xi}_N^m) \right|^2\leq \widetilde{\gamma}_N^{}, ~\forall m\right\}], \end{align}\] for the covariance matrices of the noise knowing that its norm squared is bounded by \(\widetilde{\gamma }_N\). We also write the perturbed exact matrices as \[\begin{align} \label{noisy32matrices3} \widetilde{{\boldsymbol{G}}}_N & := \boldsymbol{G}_N + {\Sigma_N^{\boldsymbol{\eta }}}, \quad \widetilde{{\boldsymbol{C}}}_N:= \boldsymbol{C}_N+ {\Sigma_N^{{\boldsymbol{\xi}},{\boldsymbol{\eta}}}}, \quad \widetilde{{\boldsymbol{T}}}_N:= \boldsymbol{T}_N + {\Sigma_N^{\boldsymbol{\xi }}} \end{align}\tag{32}\] and define the perturbed matrix of the operator as \[\begin{align} \widetilde{{\boldsymbol{A}}}_N ^T& := \widetilde{{\boldsymbol{C}}}_N \widetilde{{\boldsymbol{G}}}_N^+. \end{align}\]

Proposition 7 shows that, in this setting and stressing the dependence of \(\delta\) on \(M\), the error is of the form \(\lambda_N(\delta_M+\sigma _N^2)\). As a result, the error will not be small unless the variance of \({\boldsymbol{\eta }}\) and covariance of \({\boldsymbol{\eta }}\) and \({\boldsymbol{\xi }}\) are small. Additionally, the mass matrix \(\widetilde{{\boldsymbol{G}}}_{NM}\) may cease to be invertible unless \(\delta +{\left\lVert {\boldsymbol{\Sigma }}_N^{\boldsymbol{\eta}} \right\rVert<\left\lVert \boldsymbol{G}_N^{-1} \right\rVert^{-1}}\). This is in contrast to Proposition 6 and Proposition 8 of the next subsection where the error is of the form \(\lambda_N \delta_M\). In these cases, the error can be made small by taking \(M\) large enough.

5.2 Unbiased estimation with noisy measurements via batching↩︎

In this section we show how to correct for the bias in the data-driven approximation of the operator \(\mathcal{A}\). The idea is to evaluate the dictionary functions twice independently, and then use the average of these evaluations to form the data-driven approximation. Consider the noisy maps \[\begin{align} {\boldsymbol{x}} \to \psi_n(\boldsymbol{x})+\eta_N^{{\boldsymbol{x}},n}, \quad {\boldsymbol{x}} \to \mathcal{A}\psi_n(\boldsymbol{x})+\xi_N^{{\boldsymbol{x}},n}, \quad n=1,\dots,N, \end{align}\] where \(\boldsymbol{\eta}_N^{{\boldsymbol{x}}}= (\eta_N^{{\boldsymbol{x}},1},\ldots,\eta_N^{{\boldsymbol{x}},N})\) and \(\boldsymbol{\xi}_N^{{\boldsymbol{x}}}=(\xi_N^{{\boldsymbol{x}},1},\ldots,\xi_N^{{\boldsymbol{x}},N})\) represent the measurement or evaluation error in dictionary and operator, respectively. Consider now \({\boldsymbol{x}}_1, \cdots, {\boldsymbol{x}}_M \sim \mu\). For each \({\boldsymbol{x}}_m\) we perform these noisy evaluations, evaluating the dictionary twice independently to obtain the data \[\label{data2} \{\psi_n(\boldsymbol{x}_m)+\eta_N^{m,n}, \psi_n(\boldsymbol{x}_m)+\overset{\circ}{\eta}_N^{m,n},\mathcal{A}\psi_n(\boldsymbol{x}_m)+\xi_N^{m,n}\}_{m,n=1}^{M,N}\tag{40}\] where, we wrote \(\boldsymbol{\eta }_N^m:=(\eta_N^{m,1},\dots, \eta_N^{m, N}), \overset{\circ}{\boldsymbol{\eta}}_N^m:=(\overset{\circ}{\eta}_N^{m,1},\dots, \overset{\circ}{\eta}_N^{m, N})\), and \(\boldsymbol{\xi }_N^m:=(\xi_N^{m,1},\dots, \xi_N^{m, N})\). Given the perturbed data 40 , we form the data-driven approximation \[\overset{\circ}{{\boldsymbol{A}}}_{NM}^\top := \overset{\circ}{{\boldsymbol{C}}}_{NM}\overset{\circ}{{\boldsymbol{G}}}_{NM}^+,\] with perturbed structure and Gram matrices defined respectively by \[\require{physics} \begin{align} \label{noisy32matrices} [\overset{\circ}{\boldsymbol{{C}}}_{NM}]_{ij} & :=\frac{1}{M}\sum_{m=1}^M\qty(\mathcal{A}\psi_i(\boldsymbol{x}_m)+{\xi}_N^{m,i})\qty(\psi_j(\boldsymbol{x}_m)+\overset{\circ}{{\eta}}_N^{m,j})^\dagger, \\ [\overset{\circ}{\boldsymbol{{G}}}_{NM}]_{ij} & :=\frac{1}{M}\sum_{m=1}^M\qty(\psi_i(\boldsymbol{x}_m)+{\eta}_N^{m,i})\qty(\psi_j(\boldsymbol{x}_m)+\overset{\circ}{{\eta}}_N^{m,j})^\dagger. \end{align}\tag{41}\] In the same line as Assumption 6 we make the following assumption on the noise

Given \(\overset{\circ}{\gamma }_N \geq 0\) we will use the notation \[\require{physics} \overset{\circ}{p}_N:=\mathbb{P}\qty[\left| (\boldsymbol{\eta}_N^m,\boldsymbol{\xi}_N^m,\overset{\circ}{\boldsymbol{\eta}}_N^m) \right|^2\leq \overset{\circ}{\gamma}_N^{}],\] and define the conditioned noise covariances \[\require{physics} \begin{align} {\boldsymbol{\Sigma}}_{N}^{\boldsymbol{\eta}}&:=\mathbb{E}\qty[\boldsymbol{\eta}_N^m(\boldsymbol{\eta}_N^m)^\dagger\mid \left\{\left| (\boldsymbol{\eta}_N^{_m}, \boldsymbol{\xi}_N^{m},\overset{\circ}{\boldsymbol{\eta }}_N^m) \right|^2\leq \overset{\circ}{\gamma}_N^{}, ~\forall m\right\}],\\ {\boldsymbol{\Sigma}}_{N}^{\boldsymbol{\xi}}&:=\mathbb{E}\qty[\boldsymbol{\xi}_N^m(\boldsymbol{\xi}_N^m)^\dagger\mid \left\{\left| (\boldsymbol{\eta}_N^{_m}, \boldsymbol{\xi}_N^{m},\overset{\circ}{\boldsymbol{\eta }}_N^m) \right|^2\leq \overset{\circ}{\gamma}_N^{}, ~\forall m\right\}]. \end{align}\] Note that since \(\boldsymbol{\eta}_N^m\) and \(\overset{\circ}{\boldsymbol{\eta}}_N^m\) are identically distributed, the conditioned covariance of \(\overset{\circ}{\boldsymbol{\eta}}_N^m\) equals \({\boldsymbol{\Sigma}}_{N}^{\boldsymbol{\eta}}\).

Assumption 7 is satisfied by Gaussian measurement error and, by studying the cumulative distribution function of the \(\chi^2\) distribution, an explicit expression for \(\overset{\circ}{p}_N\) can be given. We do so in Example 11 below, after adapting our previous estimates to the setting of noisy data, starting with the result analogous to Proposition 5.

The generalisation of Theorem 6 to the case with noise can be proved in the same way as said theorem, but this time using Proposition 8 instead of Proposition 5. This gives the following result.

Thus, Theorem 9 shows an order of convergence of \(\overset{\circ}{\mathcal{A}}_{NM}\) to \(\mathcal{A}_N\) of \[\require{physics} \label{order32of32convergence32noise} \begin{align} M &=\mathcal{O}\left((\gamma_N^{}+\overset{\circ}{\gamma}_N^{}) \max \left\{\left\lVert \boldsymbol{G}_N \right\rVert+\left\lVert {\boldsymbol{\Sigma}}_{N}^{{\boldsymbol{\eta}}} \right\rVert,\left\lVert \boldsymbol{T}_N+ {\Sigma_N^{\boldsymbol{\xi }}} \right\rVert\right\}\kappa \qty(\boldsymbol{G}_N)\left\lVert \boldsymbol{C}_N \right\rVert^2 \left\lVert \boldsymbol{G}_N^{-1} \right\rVert^4\right.\\&\cdot\left. \log \left(\frac{ N}{1-p/\overset{\circ}{p}_N^M}\right) \varepsilon ^{-2}\right). \end{align}\tag{47}\]

We now discuss this result in the context of Gaussian noise.

6 Weak spectral convergence↩︎

In Sections 3, 4, and 5, we established the convergence of the data-driven operators \(\widehat{\mathcal{A}}_{NM}\) when \(M\) goes to infinity, the convergence of the Galerkin approximations \(\mathcal{A}_{N}\) when \(N\) goes to infinity, and the joint convergence of \(\hat{\mathcal{A}}_{NM}\) when \(N,M\) go to infinity, respectively. This section establishes convergence results for the eigenvalues and eigenfunctions of the approximations \(\widehat{\boldsymbol{A}}_{NM}\) and \(\widehat{\mathcal{A}}_{NM}\).

Eigenvalues and eigenfunctions of transfer operators play a vital role in the global analysis of complex dynamical systems. The dominant eigenvalues are related to the slowest timescales of the underlying system and the corresponding eigenfunctions contain important information about slowly evolving spatiotemporal patterns and have been, for instance, used to detect stable conformations of molecules or gyres in the ocean, see [6], [7] for more details.

For a general non-normal operator \(\mathcal{A}\), the spectrum can be highly sensitive to the perturbations introduced by the numerical approximation. This is characterized by the operator’s pseudospectrum, which can explain the appearance of spurious eigenvalues (spectral pollution) [38]. Additional challenges arise when the underlying dynamical system has a continuous spectrum, which can cause approximate eigenfunctions to converge weakly to zero [39]. For practical approaches designed to mitigate spectral pollution by discarding spurious eigenpairs, we refer the reader to residual-based methods, see [37] and the approximation of the Koopman generator of continuous time measure preserving flows by operators with discrete spectrum [40].

Conversely, the spectrum of a normal operator is stable and insensitive to small perturbations. However, this stability is only realized computationally if the discretization scheme preserves this normality. Standard Galerkin methods do not generally satisfy this property. As a result, spectral pollution can still appear even for well-behaved normal operators, as demonstrated in [38]. The most reliable cases for spectral approximation are therefore situations where the operator is normal and the chosen numerical method is designed to preserve this structure.

We again discuss large data and large dictionary limits, first separately and then jointly, establishing three convergence results below. We first remind the reader of the definition of weak convergence in a Hilbert space. Let \((f_n)_{n=1}^\infty\) be a sequence taking values in a Hilbert space \((\mathcal{H}, \langle\cdot, \cdot\rangle_\mathcal{H})\) and \(f \in \mathcal{H}\). Then \((f_n)_{n=1}^\infty\) converges weakly to \(f\), denoted by \(f_n \stackrel{w}{\to} f\), if \[\left\langle f_n, g\right\rangle_\mathcal{H}\rightarrow\langle f, g\rangle_\mathcal{H}\text{ for every } g\in \mathcal{H}.\] In what follows it will be important that the weak limit \(f\) of the eigenfunctions is non-zero, see Example 2 in [39].

We now establish a convergence result for eigenpairs in the large data limit.

We now move on to the convergence result in the large dictionary limit. The following theorems include the additional assumption that the limit of the eigenfunctions is nonzero. This does not always hold, as mentioned at the beginning of the section.

Conditions under which such a sequence exists can be found in, for example, [42], which rely on \(\mathcal{A}_n \approx \mathcal{A}\) and \(\mathcal{D}_n \approx \mathcal{D}\).

From Theorem 4, we finally obtain the following result regarding the convergence of eigenvalues and eigenfunctions in the joint large data and dictionary limit.

7 Numerical experiments↩︎

We will now illustrate the derived convergence results and error bounds and test their sharpness.

7.1 Benchmark problems↩︎

We analyse transfer operators associated with deterministic and stochastic systems of the form 4 . In particular, we consider:

1. The ODE defined by \[\label{ODE1} \boldsymbol{b}({\boldsymbol{x}}) = \begin{bmatrix} \gamma x_1 \\ \delta (x_2 -x_1^2) \end{bmatrix},\tag{52}\] with \(\gamma =-0.8, \delta = -0.7\).

2. The overdamped Langevin dynamics corresponding to the two-dimensional double-well potential \(V(x)=(x_1^2-1)^2 + x_2^2\) with anisotropic diffusion, whose drift and diffusion terms are given by \[\label{double32well} \boldsymbol{b}(\boldsymbol{x})=\left[\begin{array}{c} 4 x_1-4 x_1^3 \\ -2 x_2 \end{array}\right], \quad \boldsymbol{\sigma(x)}=\left[\begin{array}{cc} 0.7 & x_1 \\ 0 & 0.5 \end{array}\right],\tag{53}\] respectively.

3. The one-dimensional Ornstein–Uhlenbeck process given by \[\label{OU} b(x)=- \alpha x , \quad \sigma(x)= \sqrt{\frac{1}{2 \beta}},\tag{54}\] This SDE has a unique invariant distribution \(\require{physics} \mathcal{N}\qty(0,\frac{1}{2\beta})\). We choose \(\alpha=1\) and \(\beta=2\).

For each of the dynamical systems, we aim to reconstruct the Koopman generator, the Perron Frobenius generator or the Koopman operator. We will specify the precise setup below. In each case, we compute the normalized error \[\label{normalized32error} \varepsilon_N := \frac{\big\|\boldsymbol{A}_N-\widehat{\boldsymbol{A}}_{NM}\big\|}{\big\|\boldsymbol{A}_N\big\|}\tag{55}\] of the data-driven matrix \(\widehat{\boldsymbol{A}}_{NM}\) with respect to the true Galerkin projection \(\boldsymbol{A}_N\) in 5 , where we consider as a proxy for \(\boldsymbol{A}_N\) the result of gEDMD with a large number of data points. In the two-dimensional case, we use the domain \(\mathbb{X}=[-2,2] \times [-1,1]\) and in the one-dimensional case, we define \(\mathbb{X}=[-2,2]\). In all cases, we consider the Lebesgue measure on \(\mathbb{X}\). For the three dynamical systems introduced above, we now study three different scenarios: the large data limit, the large dictionary limit, and the large data limit when using noisy data.

7.2 Numerical results as the number of data points tends to infinity↩︎

In our first experiment, we choose as basis functions monomials of order up to \(k\), defined by \[\Psi^{{\mathrm{MON}}} = \left\{\boldsymbol{x}^\alpha \right\}_{\left| \alpha \right| \leq k},\] and Gaussians centred at a collection of equidistant grid points \(\{\boldsymbol{p}_{n}\}_{n=1}^N\), i.e., \[\require{physics} \Psi^{{\mathrm{GSN}}} = \left\{\exp\qty(-\frac{\left\lVert \boldsymbol{x}-\boldsymbol{p}_n \right\rVert^2}{2 \theta^2})\right\}_{n=1}^N,\] where for two-dimensional problems we use \[\require{physics} \big\{\boldsymbol{p}_{n}\big\}_{n=1}^N=\left\{\qty(\frac{i}{2}-2,\frac{j}{2}-1), i=0,\dots,8, j=0,\dots,4\right\},\] and for one-dimension problems \[\{\boldsymbol{p}_{n}\}_{n=1}^N=\left\{\frac{i}{2}-2, i=0,\dots,8\right\}.\]

These basis functions are independent and smooth, thus satisfying Assumption 1. Moreover, since \(\mathbb{X}\) is compact, they satisfy the boundedness assumption of Assumption 5. All these basis functions are dense in \(H^k(\mathbb{X})\) for any \(k \in \mathbb{N}_{\geq 0}\) , see [43] and [44]. All the operators considered have \(\mathcal{D}\) as a Sobolev space and as a result, Assumptions 3 and 4 are satisfied.

We also use a finite element method (FEM) basis of piecewise linear functions with \(0\) boundary condition on a uniform mesh. That is, if we write \(\left\{\boldsymbol{v}_j\right\}_{j=1}^N\) for the non-boundary vertices of the mesh, we have that \(\psi^{{\mathrm{FEM}}}\) are the only continuous piecewise linear functions that satisfy \[\psi^{{\mathrm{FEM}}}_i(\boldsymbol{v}_j)= \delta _{ij}, \quad\forall i,j=1,\dots,N,\] and write \(\Psi^{{\mathrm{FEM}}} = \left\{\psi^{{\mathrm{FEM}}}_n\right\}_{n=1}^N\). The basis functions in \(\Psi^{{\mathrm{FEM}}}\) are only once weakly differentiable and as a result do not belong to \(\mathcal{D}\) when the operator is a second-order differential operator. However, the structure matrix \(\boldsymbol{C}\) can still be calculated for second-order operators using Observation 1.

We set the maximum degree of the monomials to \(k=8\) so that \(N = \binom{k+d}{k}\). That is, \(N=45\) when the dimension is \(d=2\) and \(N=9\) when \(d=1\). Additionally, we take a uniform mesh with \(45\) and \(9\) non-boundary nodes in the \(2D\) and \(1D\) cases, respectively. In this way, there are the same number of observables in \(\Psi^{{\mathrm{MON}}},\Psi^{{\mathrm{GSN}}}\) and \(\Psi^{{\mathrm{FEM}}}\) . We then set \(\theta = \frac{1}{2N}\) so that the basis functions are well separated. That is, \(\theta < \frac{1}{2}\min_{i\ne j}\left\lVert \boldsymbol{p}_i-\boldsymbol{p}_j \right\rVert\). The points \(\boldsymbol{x}_1,\ldots,\boldsymbol{x}_M\) are sampled uniformly and independently from \(\mathbb{X}\) with the Lebesgue measure.

Since we know the dynamical system, the operator \(\mathcal{A}\) and its action on the basis functions is known exactly. This allows for the exact computation of \(\widehat{{\boldsymbol{C}}}_{NM}, \widehat{{\boldsymbol{G}}}_{NM}\) and \(\widehat{{\boldsymbol{A}}}_{NM}\) in 10 –11 . We apply this to approximate the Koopman semigroup, and the Koopman and Perron–Frobenius generators associated with 52 , 53 , and 54 for an increasing number of data points using \(M=2^8,2^9,\dots,2^{19}\) and compute 55 , where, since in general we do not have access to \(\boldsymbol{A}_N\), we approximate it by \(\widehat{\boldsymbol{A}}_{NM}\) with \(M=2^{20}\). We repeat this process \(50\) times for each \(M\) to calculate the average normalized operator error \(\varepsilon\). In Figure 1, we plot in log-log scale the relationship between \(M\) and \(\varepsilon\), including a 95% confidence interval for the error. To serve as a reference, we show dashed lines with slope \(-\frac{1}{2}\) and \(-1\), respectively.

As can be seen, for all choices of basis functions and all systems, the error has a slope of approximately \(-\frac{1}{2}\). This is in accordance with Theorem 6 as when \(N\) is fixed we obtain \(\varepsilon = \mathcal{O}(M^{-\frac{1}{2}})\). In Figures 1 (b) and 1 (c), we see that the error of the approximation using FEM basis functions is quite large. This is to be expected as for small \(M\) it is possible for an element of the mesh to have few points. If this is the case, the empirical Gram matrix is close to singular (see the comments around Example 3). This also results in initially large confidence interval for the error, which, when represented on a log-log plot, creates a strong visual effect. The error decays at the expected rate. In Figure 1 (d), the error using monomial basis functions becomes zero. This is because the subspace spanned by the monomial basis functions is invariant under the Koopman generator of the OU system, see Corollary 1.

7.3 Numerical results as the number of dictionary elements tends to infinity↩︎

We now study the effect of the number of dictionary elements on the operator error. To do so, we consider as our basis functions the dictionary comprising Gaussian functions from the previous subsection. We partition the domain into \(N\) equally sized quadrants and define \(\{\boldsymbol{p}_n\}_{n=1}^N\) to be the centres of these quadrants. We apply the data-driven algorithms with \(M=10^4\) data points. We repeat this process \(50\) times and calculate the average normalized operator error where again \(\boldsymbol{A}_N\) is approximated by \(\widehat{\boldsymbol{A}}_{NM}\) with \(M=10^{5}\). We then increase \(N\) from \(4\) to \(1024\). In this case, the monomial basis functions are not chosen as for higher orders, the matrix \(\boldsymbol{G}_N\) becomes ill-conditioned (for example, for monomials of order \(10\) in two dimensions, we have that \(\kappa(\boldsymbol{G}_N)\geq 10^{28}\)). When approximating the Koopman operator for the OU process, we also use the piecewise linear FEM basis functions \(\Psi^{{\mathrm{FEM}}}\). However, these are not used for gEDMD as the theoretical error requires the calculation of \(\boldsymbol{T}_{ij}= \left\langle\mathcal{L}\psi_i,\mathcal{L}\psi_j \right\rangle\). For gEDMD, \(\mathcal{L}\) is a second order operator and \(\psi_i^{{\mathrm{FEM}}}\) are not twice differentiable. Since we plot normalised errors, we cap the theoretical bound in Theorem 6 in the plots at one. As can be seen in Figure 2, the error increases with the number of observables as expected in view of Theorem 6. We observe that the theoretical bound on the error increases faster than the simulation error. This indicates that perhaps there is some underlying structure which could make the bounds tighter. Additionally, we see that the confidence interval for the FEM basis functions becomes larger as the number of basis functions increases. This is because of the increased likelihood that very few particles end up in some interval of the partition. This leads to close to a singular mass matrix and a large condition number. See the comments around Example 3.

7.4 Numerical results with noise↩︎

In this section, we repeat the experiments carried out in Section 7.2 with the addition of a noise term as in 40 , where we take normal i.i.d.noise, i.e., \[\require{physics} \qty(\boldsymbol{\eta}_n,\boldsymbol{\xi}_n) \sim \mathcal{N}\qty(0,\sigma^2 \boldsymbol{I}_{2n \times 2n}).\] We increase the noise using \(\sigma = 10^{-3}, 10^{-2}, 10^{-1}\) and study its effect on the normalised error. To limit the number of plots, we restrict ourselves to gEDMD for the ODE 52 and gEDMD for the Perron–Frobenius operator associated with the Ornstein–Uhlenbeck process.

As we can see, the Gaussian and FEM observables are more resilient to increased noise. The FEM basis functions perform the best in the presence of noise. This is expected as we evaluated the FEM basis functions exactly to \(0\) on all points outside of their support and only added noise to the non-zero ones. The lower resilience of the monomials compared to the other two is due to the fact that the condition number of the Gram matrix of the monomials is larger and thus the inverse of the Gram matrix is very sensitive to noise. In the second and third rows of Figure 3, the monomials are slightly more resilient to the noise. This, however, is due to the fact that, since the dimension of the domain is \(1\) instead of \(2\), there are fewer monomials and thus the condition number of the Gram matrix is smaller. It can be seen that even when the error without noise is exactly zero (see Figure 1 (d)), the monomial basis functions are not resilient to noise.

7.5 Numerical results on spectral convergence↩︎

The results contained in Section 6 do not provide bounds on the speed of convergence of the spectrum. In this section, we numerically compute the difference between the eigenvalues of \(\widehat{{\boldsymbol{A}}}_{NM}\) and those of \(\widehat{{\boldsymbol{A}}}_N\). As an example, we consider the Ornstein–Uhlenbeck process 54 since, when using the monomials as basis functions, the eigenvalues of their projection onto the space of monomials are given by

\[\begin{align} \mathrm{spec}(\mathcal{L}_N )&= -n \alpha \quad \quad n=0,1,\ldots,N-1, \\ \mathrm{spec}(\mathcal{K}_N^t )&= e^{-n \alpha t}, \quad n=0,1,\ldots,N-1. \end{align}\] We use the monomial, Gaussian, and FEM basis functions of Section 7.2 and calculate the spectral error \[\begin{align} \label{spectral32error} \varepsilon_{\mathrm{spec}}:= \left\lVert \mathrm{spec}(\widehat{{\boldsymbol{A}}}_{NM})-\mathrm{spec}(\widehat{{\boldsymbol{A}}}_N) \right\rVert_{\mathbb{R}^N} = \left(\sum_{n=1}^N \left| \lambda^{(n)}_{NM} - \lambda^{(n)}_{N} \right|^2\right)^{1/2}, \end{align}\tag{56}\] where \(\lambda^{(n)}_{NM}\) and \(\lambda^{(n)}_{N}\) are the ordered \(n\)-th eigenvalues of \(\widehat{{\boldsymbol{A}}}_{NM}\) and \(\widehat{{\boldsymbol{A}}}_{N}\), respectively.

As before, we use an increasing number of data points using \(M=2^8,2^9,\dots,2^{19}\), compute 56 and repeat the approximation \(50\) times for each \(M\) to calculate the average spectral error \(\varepsilon_{\mathrm{spec}}\). In Figure 4, we plot in log-log scale the relationship between \(M\) and \(\varepsilon_{\mathrm{spec}}\), including a 95% confidence interval for the error. To serve as a reference, we show dashed lines with slope \(-\frac{1}{2}\) and \(-1\), respectively.

Since the subspace spanned by the monomial basis functions is invariant under the Koopman generator, the error using monomial basis functions becomes zero, see Corollary 1. Otherwise, the error has a slope of approximately \(-\frac{1}{2}\). This suggests that, in some particular cases, error bounds similar to those of Section 4 may hold for the spectrum.

8 Conclusion↩︎

In this article, we have investigated the approximation of an operator \(\mathcal{A}\) from point evaluations of a dictionary to which the operator has been applied. That is, we have assumed that we have training data of the form \[\big\{\psi_n(\boldsymbol{x}_m), \mathcal{A}\psi_n (\boldsymbol{x}_m)\big\}_{m,n=1}^{M,N},\] where the evaluations could additionally be subject to random noise. After describing the estimation procedure for linear operators from data, we have presented a thorough convergence and error analysis in the dictionary limit (\(N \rightarrow \infty\)), the data limit (\(M \rightarrow \infty\)), and their joint limit (\(N, M \rightarrow \infty\)). We have studied the convergence of the full operators as well as their spectra.

Throughout this work, we have usually thought of the approximation of transfer operators in the context of dynamical systems, such as Koopman operators, Perron–Frobenius operators, and their generators. The framework we have presented is clearly not limited to such operators. Indeed, it generalises approximation techniques for transfer operators, such as EDMD and gEDMD, to general operators. EDMD and gEDMD fall naturally into our framework and we have shown significant new convergence results and error analyses for these methods. Additionally, we have verified our analytical results in numerical experiments.

9 When does a matrix determine an operator?↩︎

As discussed in Section 2.1, a matrix may not always define an operator when it is interpreted as acting on a set of vectors that are not linearly independent. The following lemma gives necessary and sufficient conditions for such an operator to be well-defined.

10 Bernstein inequality↩︎

The following results can be found in [45].

11 Table of notation↩︎

The notation used throughout the manuscript is summarised in Table 1.

References↩︎

---

1. L.L.and S.L.contributed equally.↩︎