Delay-coordinate maps and the spectra of Koopman operators

TL;DR

This paper introduces a method to extract Koopman eigenfunctions from high-dimensional data using delay-coordinate maps and kernel operators, facilitating analysis of complex dynamical systems with mixed spectra.

Contribution

It presents a novel approach combining delay coordinates and kernel integral operators to approximate Koopman eigenfunctions, especially in systems with mixed spectra.

Findings

Efficient approximation of Koopman eigenfunctions from high-dimensional data.

Delay-coordinate maps enable convergence to the discrete spectrum of the Koopman operator.

Method applicable to systems with pure point or mixed spectra.

Abstract

The Koopman operator induced by a dynamical system is inherently linear and provides an alternate method of studying many properties of the system, including attractor reconstruction and forecasting. Koopman eigenfunctions represent the non-mixing component of the dynamics. They factor the dynamics, which can be chaotic, into quasiperiodic rotations on tori. Here, we describe a method through which these eigenfunctions can be obtained from a kernel integral operator, which also annihilates the continuous spectrum. We show that incorporating a large number of delay coordinates in constructing the kernel of that operator results, in the limit of infinitely many delays, in the creation of a map into the discrete spectrum subspace of the Koopman operator. This enables efficient approximation of Koopman eigenfunctions from high-dimensional data in systems with pure point or mixed spectra.

Figures

Equations

$$(U^{t}f):x\mapsto f\left({\Phi }^t(x)\right),\text{for μ-a.e. x∈X}.$$

$$U^{t}z=\exp (i\omega t)z.$$

$$L^2(X,\mu )=\mathcal{D}\oplus {\mathcal{D}}^{\perp }.$$

$$Kf(x):={\int }_{X}k(x,y)f(y)d\mu (y).$$

$$d_Q^2(x,y)=\frac{1}{Q}\sum _{q=0}^{Q-1}{\Vert F({\Phi }^{q\Delta t}(x))-F({\Phi }^{q\Delta t}(y))\Vert }^2.$$

$$k_Q(x,y)=e^{-d_Q^2(x,y)/ϵ}.$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{n=0}^{N-1}f({\Phi }^{n\Delta t}(x))={\int }_{X}f(y)d\mu (y),\forall f\in C^0(M).$$

$$Vf:=\underset{t\to 0}{\lim }\frac{1}{t}\left(U^{t}f-f\right),f\in D(V).$$

$$Vz=i\omega z.$$

$$U^t{z}_i=\exp (it{\omega }_i)z_i,i\in \{1,2\},$$

$$⟹U^t(z_1{z}_2)=(U^t{z}_1)(U^t{z}_2)=\exp (it({\omega }_1+{\omega }_2))z_1{z}_2.$$

$$z_{\vec{a}}=\prod _{j=1}^m{z}_1^{a_1}\cdots z_m^{a_m},$$

$$U^{t}f=\sum _{\vec{a}\in {\mathbb{Z}}^m}{\hat{f}}_{\vec{a}}{e}^{i{\omega }_{\vec{a}}t}{z}_{\vec{a}},{\hat{f}}_{\vec{a}}=⟨z_{\vec{a}},f⟩,$$

$$\pi (x)=(z_{{\vec{a}}_1}(x),\dots ,z_{{\vec{a}}_l}(x)),\text{for μ-a.e. x∈X},$$

$${\pi }_{\omega }f=\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{n=0}^{N-1}{e}^{-i\omega n\Delta t}{U}^{n\Delta t}f.$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{n=0}^{N-1}{g}^{\ast }({\Phi }^{n\Delta t}(x))f({\Phi }^{n\Delta t}(y))=0.$$

$$\begin{gathered} \underset{N\to \infty }{\lim }\frac{1}{N}\sum _{n=0}^{N-1}{g}^{\ast }({\Phi }^{n\Delta t}(x))f({\Phi }^{n\Delta t}(y)) \\ =g^{\ast }(x)\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{n=0}^{N-1}{e}^{-in\omega \Delta t}f({\Phi }^{n\Delta t}(y)), \end{gathered}$$

$$k_Q(x,y)=h(d_Q(x,y)),$$

$${\tilde{k}}_Q(x,y)=\frac{1}{Q}\sum _{q=0}^{Q-1}h(d_1({\Phi }^{q\Delta t}(x),{\Phi }^{q\Delta t}(y))),$$

$$\begin{array}{c}\Psi (F)(x,y):=\underset{Q\to \infty }{\lim }{\Psi }_Q(F)(x,y),\\ {\Psi }_Q(F)(x,y):=\frac{1}{Q}\sum _{q=0}^{Q-1}{\Vert F({\Phi }^{q\Delta t}(x))-F({\Phi }^{q\Delta t}(y))\Vert }^2.\end{array}$$

$$\Psi (F)=\Psi (F_{\mathcal{D}})+2\Vert F_{{\mathcal{D}}^{\perp }}{\Vert }_{L^2}^2.$$

$$\begin{array}{rl}\Psi (F)(x,y)& =\underset{Q\to \infty }{\lim }\frac{1}{Q}\sum _{q=0}^{Q-1}\left({\Vert G_{\mathcal{D}}(x_q,y_q)\Vert }^2+{\Vert G_{{\mathcal{D}}^{\perp }}(x_q,y_q)\Vert }^2\right)\\ & -2\underset{Q\to \infty }{\lim }\frac{1}{Q}\sum _{q=0}^{Q-1}{G}_{\mathcal{D}}(x_q,y_q)\cdot G_{{\mathcal{D}}^{\perp }}(x_q,y_q),\end{array}$$

$$\underset{Q\to \infty }{\lim }\frac{1}{Q}\sum _{q=0}^{Q-1}\left(F_{{\mathcal{D}}^{\perp }}(x_q)-F_{{\mathcal{D}}^{\perp }}(y_q)\right)\cdot \left(F_{\mathcal{D}}(x_q)-F_{\mathcal{D}}(y_q)\right)=0,$$

$$\begin{array}{rl}(\Psi F_{{\mathcal{D}}^{\perp }})(x,y)& =\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{n=0}^{N-1}|F_{{\mathcal{D}}^{\perp }}(x_n){|}^2+\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{n=0}^{N-1}|F_{{\mathcal{D}}^{\perp }}(y_n){|}^2\\ & +2\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{n=0}^{N-1}{F}_{{\mathcal{D}}^{\perp }}(x_n)F_{{\mathcal{D}}^{\perp }}(y_n).\end{array}$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{n=0}^{N-1}{F}_{{\mathcal{D}}^{\perp }}(x_n)F_{{\mathcal{D}}^{\perp }}(y_n)=⟨J,1_{X\times X}⟩=0.$$

$$d_{\infty }(x,y):=\underset{Q\to \infty }{\lim }{d}_Q(x,y);(x,y)\in D(d_{\infty })\subseteq X\times X$$

$$k_{\infty }(x,y)=h(d_{\infty }(x,y)),(x,y)\in D(d_{\infty }),$$

$$\left|Kf(x_1)-Kf(x_2)\right|$$

$$\le \Vert k_{\infty }(x_1,\cdot )-k_{\infty }(x_2,\cdot ){\Vert }_{L^2}\Vert f{\Vert }_{L^2}$$

$$\le \Vert k_{\infty }(x_1,\cdot )-k_{\infty }(x_2,\cdot ){\Vert }_{L^{\infty }}.$$

$$Kf(x)={\int }_X{k}_{\infty }(x,y)f(y)d\mu (y)={\int }_X{\bar{g}}_{\infty }(x,y)f(y)d\mu (y).$$

$$\begin{array}{rl}\left|Kf(x)\right|& \le \left|{\int }_{X^{\prime }}{k}_{\infty }(x,y)f(y)d\mu (y)\right|\\ & \le \Vert k_{\infty }(x,\cdot ){\Vert }_{L^2}\Vert f{\Vert }_{L^2}\le \Vert k_{\infty }(x,\cdot ){\Vert }_{L^{\infty }}\Vert f{\Vert }_{L^2}\\ & \le \Vert k_{\infty }{\Vert }_{L^{\infty }(X\times X)}\Vert f{\Vert }_{L^2},\end{array}$$

Full text

1 Introduction

The tasks of dimension reduction and forecasting of time series are very common in physical and engineering sciences, where the time-series studied are often partial observations of a nonlinear dynamical system. A classical example of such time series is data collected from the Earth’s climate system, where many of the active degrees of freedom are difficult to access via direct observations (e.g., subsurface ocean circulation). Moreover, the available observations typically mix together different physical processes operating on a wide range of spatial and temporal scales. For instance, in the climate system, the seasonal cycle and the El Niño Southern Oscillation (the latter evolving on interannual timescales) both have strong associated signals in sea surface temperature [66]. In such applications, identifying dynamically important, coherent patterns of variability from the data can enhance our scientific understanding and predictive capabilities of complex phenomena.

Ergodic theory, and in particular its operator-theoretic formulation [15, 23], provides a natural framework to address these objectives. In this framework, the focus is on the action of the dynamical system on spaces of observables (functions of the state), as opposed to the dynamical flow itself. The advantage of this approach, first realized in the seminal work of Koopman [39], is that the action of a general dynamical system on spaces of observables is always linear. As a result, with appropriate regularity assumptions, the problem of identification and prediction of dynamically intrinsic coherent patterns can be formulated as an estimation problem for the spectrum of a linear evolution operator. In addition, for systems exhibiting ergodic behavior, spectral quantities such as eigenvalues and eigenfunctions can be statistically estimated from time-ordered data without prior knowledge of the state space geometry or the equations of motion. At the same time, spaces of observables are also infinite dimensional, so the issue of finite-dimensional approximation of (potentially unbounded) operators becomes relevant.

Starting from the techniques proposed in [22, 52, 50], the operator-theoretic approach to ergodic theory has stimulated the development of a broad range of techniques for data-driven modeling of dynamical systems. These methods employ either the Koopman [52, 50, 55, 57, 15, 62, 67, 30, 14, 2, 33, 34, 40, 41] or the Perron-Frobenius (transfer) operators [22, 21, 28, 29, 27], which are duals to one another in appropriate function spaces. The goal common to these techniques is to approximate spectral quantities for the operator in question, such as eigenvalues, eigenfunctions, and spectral projections, from measured values of observables along orbits of the dynamics. To that end, a diverse range of approaches has been employed, including state space partitions [22, 22, 21, 28, 29, 27], harmonic averaging [52, 50, 20], iterative methods [57, 55], dictionary/basis representations [62, 67, 30, 33, 34, 40], delay-coordinate embeddings [14, 2, 30, 33], and spectral-moment estimation [41].

Compared to observables identified by eigendecomposition techniques based on kernel integral operators that do not depend on the dynamics (e.g., covariance [3, 38] or heat operators [5, 16, 11], the latter of which have been popular in manifold learning applications), eigenfunctions of evolution operators are likely to offer higher physical interpretability and predictability, as they are determined from an operator intrinsic to the dynamical system. In particular, one of the key properties of Koopman or Perron-Frobenius eigenfunctions for ergodic dynamical systems is that they evolve periodically and with a single frequency (even if the underlying dynamical system is aperiodic), and thus have high predictability. This and a number of other attractive properties motivate the identification of such eigenfunctions of data.

Yet, for systems of sufficient complexity, Koopman and Perron-Frobenius operators have significantly more complicated spectral behavior than kernel integral operators, generally exhibiting a continuous spectral component and/or non-isolated eigenvalues, which presents challenges to the construction of data-driven approximation techniques with spectral convergence guarantees. Indeed, to our knowledge, spectral convergence results for the data-driven approximation of Koopman eigenvalues and eigenfunctions have been limited to special cases such as quasiperiodic rotations on tori [33], or systems observed through measurement functions lying in finite-dimensional invariant subspaces [2].

The main contribution of our work is the construction of a data-driven approximation scheme for Koopman eigenvalues and eigenfunctions that provably converges for a broad class of ergodic dynamical systems and observation maps, encompassing many of the applications encountered in the physical and engineering sciences. Our approach will be based on a combination of ideas from delay-coordinate maps of dynamical systems [56], kernel integral operators for machine learning [5, 16, 11, 9, 47], and Galerkin approximation techniques for variational eigenvalue problems [4]. Using these tools, we will construct a compact kernel integral operator that commutes with the Koopman operator in an asymptotic limit of infinitely many delays, and employ the finite-dimensional common eigenspaces of these operators as Galerkin approximation spaces for the Koopman eigenvalue problem. We will show that orthonormal bases of these spaces can be stably and efficiently approximated from finitely many measurements taken near the attractor, and the resulting data-driven Galerkin schemes converge in the asymptotic limit of large data.

2 Assumptions and statement of main results

A common underlying assumption in the statistical modeling of dynamical systems is ergodicity. This assumption encapsulates the working principle that the global properties (with respect to $\mu$) of an observable $F$ can be obtained from a time series for $F$, namely, $F(x_{0}),\ldots,F(x_{N-1})$, where $x_{0},\ldots,x_{N-1}$ is an unobserved trajectory on the state space of the dynamical system. Moreover, ergodicity implies that $L^{2}$ inner products between observables can be approximated by time-correlations. Also, our methods rely on integral operators, and these can be approximated as matrices under the ergodic hypothesis. We now make our assumptions more precise.

Assumption 1**.**

Let $M$ be a topological manifold, equipped with its Borel $\sigma$-algebra. $\Phi^{t}:M\to M$, $t\in\mathbb{R}$, is a continuous flow on $M$ with an ergodic, Borel probability measure $\mu$ with a compact support $X$. $F:M\to\mathbb{R}^{d}$ is a continuous measurement function through which we collect a time-ordered data set consisting of $N$ samples $F(x_{0}),F(x_{1}),\ldots,F(x_{N-1})$, each $F(x_{n})$ lying in $d$-dimensional data space. Here, $x_{n}=\Phi^{n\,\Delta t}(x_{0})$, and $\Delta t$ is a fixed sampling interval such that the map $\Phi^{\Delta t}$ is ergodic for the invariant measure $\mu$.

The Koopman operator.

Central to all our following discussions will be the concept of the Koopman operator. Koopman operators [15, 23, 53] act on observables by composition with the flow map, i.e., by time shifts. The space $L^{2}(X,\mu)$ of square-integrable, complex-valued functions on $X$ will be our space of observables. Given an observable $f\in L^{2}(X,\mu)$ and time $t\in\mathbb{R}$, $U^{t}:L^{2}(X,\mu)\to L^{2}(X,\mu)$ is the operator defined as

[TABLE]

$U^{t}$ is called the Koopman operator at time $t$ associated with the flow. For measure-preserving systems, $U^{t}$ is unitary, and has a well-defined spectral expansion consisting in general of both point and continuous parts lying in the unit circle [50]. The problems of mode decomposition and non-parametric prediction can both be stated in terms of the Koopman operator [33]. We will now describe an important tool for studying Koopman operators, namely their eigenfunctions.

Koopman eigenfunctions.

Every eigenfunction $z$ of $U^{t}$ satisfies the following equation for some $\omega\in\mathbb{R}$:

[TABLE]

Koopman eigenfunctions are particularly useful for prediction and dimension reduction in dynamical systems. This is because, as seen in (1), the knowledge of an eigenfunction $z$ at time $t=0$ enables accurate predictions of $z$ up to any time $t$, since $U^{t}$ operates on $z$ as a multiplication operator by a time-periodic, single-frequency multiplication factor. Moreover, it is possible to construct a dimension reduction map, sending the high-dimensional data $F(x)\in\mathbb{R}^{d}$ to the vector $(z_{1}(x),\ldots,z_{l}(x))\in\mathbb{C}^{l}$, where $l\ll d$, and the $z_{1},\ldots,z_{l}$ are Koopman eigenfunctions corresponding to rationally independent frequencies $\omega_{1},\ldots,\omega_{l}$ [50, 30, 33]. In this representation, the $z_{j}$ can be thought of as “coordinates” corresponding to distinct periodic processes operating at the timescales $2\pi/\omega_{j}$. Also of interest (and in some cases easier to compute) are the projections of the observation map $F$ onto the Koopman eigenfunctions, called Koopman modes [50]. Data-driven techniques for computing Koopman eigenvalues, eigenfunctions, and modes that have been explored in the past include methods based on generalized Laplace analysis [52, 50], dynamic mode decomposition (DMD) [58, 55, 57, 62], extended DMD (EDMD) [67, 40], Hankel matrix analysis [62, 2, 14], and data-driven Galerkin methods [30, 33, 34]. The latter approach, as well as the related work in [8], additionally address the problem of nonparametric prediction of observables and probability densities.

Let $\mathcal{D}$ be the closed subspace of $L^{2}(X,\mu)$ spanned by the eigenfunctions of $U^{t}$, and $\mathcal{D}^{\bot}$ its orthogonal complement. As is well known [37], and will be discussed in more detail in Section 3, the subspaces $\mathcal{D}$ and $\mathcal{D}^{\perp}$ represent the quasiperiodic and mixing (chaotic) components of the dynamics, respectively. Moreover, they are both invariant under $U^{t}$ for every time $t\in\mathbb{R}$, thus inducing an invariant splitting [50]

[TABLE]

Systems for which $\mathcal{D}$ contains non-constant functions and $\mathcal{D}^{\bot}$ is non-empty are called mixed-spectrum systems.

Kernel integral operators.

The method that we will describe in this paper relies heavily on kernel integral operators. A kernel is a function $k:M\times M\to\mathbb{R}$, measuring the similarity between pairs of points on $M$. Kernel functions can be of various designs, and are meant to capture the nonlinear geometric structures of data; see for example [5, 59, 16]. One advantage of using kernels is that they can be defined so as to operate directly on the data space, e.g., $k(x,y)=\kappa(F(x),F(y))$ for some function $\kappa:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ of appropriate regularity. Defined in this manner, $k$ can be evaluated using measured quantities $F(x)$ without explicit knowledge of the underlying state $x$. Associated with a square-integrable kernel $k\in L^{2}(X\times X,\mu\times\mu)$ is a compact integral operator $K:L^{2}(X,\mu)\to L^{2}(X,\mu)$ such that

[TABLE]

In some cases, we will make the following assumptions on kernels.

Assumption 2**.**

The kernel $k:M\times M\to\mathbb{R}$ is (i) symmetric and continuous; (ii) strictly positive-valued.

Overview of approach.

We will address the eigenvalue problem for $U^{t}$ by solving an eigenvalue problem for a kernel integral operator $P_{Q}$, which is accessible from data, and in the limit of $Q\to\infty$, commutes with $U^{t}$. Since commuting operators have common eigenspaces, this will allow us to compute eigenfunctions of $U^{t}$ through expansions in eigenbases obtained from $P_{Q}$. These operators have Markovian kernels $p_{Q}:M\times M\to\mathbb{R}$ (i.e., $p_{Q}\geq 0$ and $\int_{X}p_{Q}(x,\cdot)\,d\mu=1$, for $\mu$-a.e. $x\in M$), whose construction begins from a family of distance-like functions $d_{Q}:M\times M\to\mathbb{R}$, defined by

[TABLE]

Here, $Q$ is a positive integer parameter, and $\lVert\cdot\rVert$ the canonical 2-norm on $\mathbb{R}^{d}$. Intuitively, $d_{Q}(x,y)$ assigns a distance-like quantity between points $x$ and $y$ equal to the root-mean square distance between $Q$ consecutive “snapshots” of the observable $F$, measured along dynamical trajectories starting from $x$ and $y$. In other words, $d_{Q}$ corresponds to a distance between data in delay-coordinate space with $Q$ delays. Several of our results will depend on the asymptotic behavior of $d_{Q}$ as $Q\to\infty$, which we will study in detail.

Composing $d_{Q}$ with a continuous shape function $h:\mathbb{R}\to\mathbb{R}$, leads to a kernel $k_{Q}:M\times M\to\mathbb{R}$, $k_{Q}=h\circ d_{Q}$, assigning a pairwise measure of similarity between points in $M$. In this paper, we will nominally work with Gaussian shape functions, $h(s)=e^{-s^{2}/\epsilon}$, parameterized by a bandwidth parameter $\epsilon>0$, so that

[TABLE]

Such kernels satisfy Assumption 2(i), (ii). They are popular in manifold learning applications [5, 16, 11] due to their localizing behavior as $\epsilon\to 0$ and their ability to approximate heat kernels, but our results also hold for many other kernel choices; e.g., [31]. Having constructed $k_{Q}$, the kernel $p_{Q}$ associated with the integral operator $P_{Q}$ is obtained via a Markov normalization procedure [16, 11], described in Section 4.3. With these definitions, we are ready to state our main results.

Theorem 1**.**

Under Assumption 1, there exists a real, self-adjoint, ergodic, compact Markov operator $P:L^{2}(X,\mu)\to L^{2}(X,\mu)$, which commutes with $U^{t}$, and is a limit of operators $P_{1},P_{2},\ldots$ (also real, self-adjoint, ergodic, compact, and Markov) in the $L^{2}(X,\mu)$ operator-norm topology. The operators $P_{Q}$ have Markov kernels $p_{Q}:M\times M\to\mathbb{R}$ satisfying the conditions in Assumption 2, and determined from delay-coordinate mapped observations $F(x),F(\Phi^{\Delta t}(x)),\ldots,F(\Phi^{(Q-1)\Delta t}(x))$ with $Q$ delays. Moreover, the kernel $p:M\times M\to\mathbb{R}$ of $P$ lies in $L^{\infty}(X\times X,\mu\times\mu)$, and $p_{Q}$ converges to $p$ in $L^{p}(X\times X,\mu\times\mu)$ norm with $1\leq p<\infty$.

The strong convergence of the compact operators $P_{Q}$ to $P$ leads to the following spectral convergence result (e.g., Section 7 in [4] and [1]).

Corollary 2** (spectral convergence).**

Under the assumptions of Theorem 1, the following hold:

1. (i)

For every nonzero eigenvalue $\lambda$ of $P$ with multiplicity $\alpha$ and every neighborhood $S\subset\mathbb{R}$ of $\lambda$ such that $\operatorname{spec}(P)\cap S=\{\lambda\}$, there exists $Q_{0}\in\mathbb{N}_{0}$ such that for all $Q>Q_{0}$, $\operatorname{spec}(P_{Q})\cap S$ contains $\alpha$ elements converging as $Q\to\infty$ to $\lambda$. 2. (ii)

*Let $\Pi$ be any projector to the eigenspace $W_{\lambda}$ of $P$ at eigenvalue $\lambda$. Let also $\Pi_{Q}$ be any projector to the union of the eigenspaces of $P_{Q}$ corresponding to the eigenvalues in $\operatorname{spec}(P_{Q})\cap S$. Then, as $Q\to\infty$, $\Pi_{Q}$ converges strongly to $\Pi$. Moreover, the *gap (distance) between $W_{\lambda}$ and $\operatorname{ran}\Pi_{Q}$, defined as in [4], converges to zero.

Theorem 3 below is a continuation of Theorem 1, and can be used to conclude some useful properties of the operator $P$.

Theorem 3**.**

Let $\Phi^{t}$ be a measurable flow on a compact set $X$ supporting an invariant ergodic probability measure $\mu$, and $T$ be a kernel integral operator with a real-valued, symmetric kernel $\tau\in L^{2}(X\times X,\mu\times\mu)$ such that $T$ commutes with $U^{t}$ (e.g., $T=P$ from Theorem 1). Then:

1. (i)

$\tau$* lies in the tensor product subspace $\mathcal{D}\otimes\mathcal{D}$, and is invariant under the flow $\Phi^{t}\times\Phi^{t}$.* 2. (ii)

$\mathcal{D}$* and $\mathcal{D}^{\perp}$ are invariant under $T$. Moreover, $\overline{\operatorname{ran}T}$ is a subspace of $\mathcal{D}$, $\mathcal{D}^{\perp}$ is a subspace of $\ker T$, and both $\operatorname{ran}T$ and $\ker T$ are invariant under $U^{t}$.*

Moreover, if $\operatorname{ran}T$ contains non-constant functions:

1. (iii)

There exists a measurable map $\pi:X\to\mathbb{T}^{D}$ for some $D\in\mathbb{N}$, whose components consist of joint eigenfunctions of $T$ and $U^{t}$, such that $\pi$ factors $\Phi^{t}$ into a rotation on the torus by a vector $\vec{\omega}\in\mathbb{R}^{D}$, i.e., $\pi(\Phi^{t}(x))=\pi(x)+\vec{\omega}t\bmod{1}$ for $\mu$-a.e. $x\in X$. 2. (iv)

There exists a choice of dimension $D$ from (iii) and a symmetric kernel $\hat{\tau}\in L^{2}(\mathbb{T}^{D}\times\mathbb{T}^{D},\operatorname{Leb})$ on the $D$-torus, such that $\tau(x,y)=\hat{\tau}(\pi(x),\pi(y))$ for $\mu\times\mu$-a.e. $(x,y)\in X\times X$.

Note that Theorems 1 and 3 hold for operators acting on $L^{2}$ spaces only. To be able to say more about the behavior of these operators on spaces of continuous functions, an additional assumption on the Koopman eigenfunctions and the observation map will be needed. In what follows, $F_{D}$ will denote the orthogonal projection of $F$ onto the quasiperiodic subspace $\mathcal{D}$ from (2).

Assumption 3**.**

All Koopman eigenfunctions, as well as the quasiperiodic component of the observation map $F_{\mathcal{D}}$, are continuous.

Although we will explicitly assume that $F_{\mathcal{D}}$ is continuous, we are not aware of a counter-example where the observation map $F$ is continuous (in accordance with Assumption 1), the Koopman eigenfunctions $z$ are continuous, but $F_{\mathcal{D}}$ is not continuous. On the other hand, smooth dynamical systems on smooth manifolds with discontinuous Koopman eigenfunctions (and in fact, pure point spectra) are known to exist, in both discrete- [13] and continuous-time settings [19]. This indicates that the continuity requirement on Koopman eigenfunctions in Assumption 3 is complementary to the assumed continuity of the dynamical flow in Assumption 1. The following theorem establishes a number of properties of $P$ under these additional continuity assumptions.

Theorem 4**.**

Let Assumptions 1 and 3 hold. Then, the kernel $p$ of the operator $P$ from Theorem 1 is uniformly continuous on a full-measure, dense subset of $X\times X$. As a result:

1. (i)

$P$* maps $L^{2}(X,\mu)$ into the space of $\mu$-a.e. continuous functions on $X$.* 2. (ii)

$P$* compactly maps $C^{0}(X)$ into itself.* 3. (iii)

The norms of the operators $P$ in (i) and (ii) are bounded above by $\lVert p\rVert_{L^{\infty}(X\times X)}$. 4. (iv)

For every $f\in C^{0}(X)$, $P_{Q}f$ is a sequence of continuous functions converging $\mu$-a.e. to $Pf$.

Remark**.**

The class of integral operators $P_{Q}$ studied in this work has previously been used for dimension reduction and mode decomposition of high-dimensional time series (e.g., [35, 36, 7, 60]). In these works, a phenomenon called in [7] “timescale separation” was observed; namely, it was observed that at increasingly large $Q$ the eigenfunctions of $P_{Q}$ capture increasingly distinct timescales of a multiscale input signal. Theorems 1 and 3 provide an interpretation of this observation from the point of view of spectral properties of Koopman operators; in particular, from the fact that $P_{Q}$ has, in the limit $Q\to\infty$, common eigenfunctions with $U^{t}$ and the latter capture distinct timescales associated with the eigenfrequencies $\omega$. Even though in this work we focus on the class of Markov operators $P_{Q}$, analogous results also hold for other classes of integral operators for data analysis that employ delays, including the covariance operators used in singular spectrum analysis (SSA) [54, 12, 65] and the related Hankel matrix analysis [62, 14, 2]. Collectively, these results establish a connection between two major branches of data analysis techniques for dynamical systems, namely those based on Koopman operators, and those based on kernel integral operators.

Theorems 1–4 are proved in Section 5. A result analogous to Theorem 1, but restricted to smooth manifolds, smooth observation maps, and Koopman operators with pure point spectrum and smooth eigenfunctions, was presented in [33]. Theorem 1 generalizes this result to non-smooth state spaces and Koopman operators with mixed spectra. With this result, the eigenvalues and eigenfunctions of $P_{Q}$ consistently approximate those of $P$, and the latter can be used in turn to construct orthonormal bases of Koopman eigenspaces. The availability of such bases is useful in many applications, including approximation techniques for the eigenvalues and eigenfunctions of $U^{t}$ or its generator (defined in Section 3 ahead). One such technique will be presented in Section 6, utilizing the eigenvalues and eigenfunctions of $P$ to perform diffusion regularization of the generator, and then solving the eigenvalue problem for the generator via a Petrov-Galerkin method. Note that the Markov property of $P$ is not trivial; for instance, it does not hold for covariance kernels. The commutativity between $U^{t}$ and $P$, in conjunction with the Markov property, lead to well posedness of these schemes despite the presence of a continuous spectrum of $V$.

Physical measures.

A point $x\in M$ is said to be in the basin of the measure $\mu$ with respect to the discrete-time map $\Phi^{\Delta t}$ if

[TABLE]

The basin ${\mathcal{B}_{\mu}}$ of an invariant ergodic measure $\mu$ always includes $\mu$-a.e. point in the support of $\mu$ (in this case, $X$), and is a forward-invariant set. An important property that we need the invariant measure $\mu$ to have is that it is physical [68]. Moreover, we will require that the dynamics has a suitable absorbing ball property. These assumptions can be summarized as follows:

Assumption 4**.**

The set ${\mathcal{B}_{\mu}}$ of points satisfying (6) has positive Lebesgue measure, i.e., the measure $\mu$ is physical. Moreover, there exists a subset $\mathcal{V}\subseteq\mathcal{B}_{\mu}$, also of positive Lebesgue measure, such that for every $x_{0}\in\mathcal{V}$ there exists a compact set $\mathcal{U}$ (which may depend on $x_{0}$, and necessarily includes $X$), such that the orbit $x_{n}=\Phi^{n\,\Delta t}(x_{0})$ enters $\mathcal{U}$ and never leaves it.

Examples where Assumption 4 is satisfied include: (i) ergodic flows on compact manifolds with Lebesgue absolutely continuous, fully supported, invariant measures, in which case $\mathcal{U}=\mathcal{V}=\overline{\mathcal{B}_{\mu}}=M=X$; (ii) certain classes of dissipative flows on potentially noncompact manifolds (e.g., the Lorenz 63 (L63) system on $M=\mathbb{R}^{3}$ [45] studied in Section 8 ahead); and (iii) certain classes of dissipative partial differential equations possessing inertial manifolds and physical measures [46, 44].

The following result shows that under Assumptions 1–4, the nonzero eigenvalues of $P_{Q}$ and the corresponding (continuous) eigenfunctions can be approximated to any degree of accuracy by data-driven operators $P_{Q,N}$, acting on the finite-dimensional Hilbert space $L^{2}(\mathcal{U},\mu_{N})$ associated with the sampling probability measure $\mu_{N}=\sum_{n=0}^{N-1}\delta_{x_{n}}/N$. These operators are constructed from time-ordered measurements $F(x_{0}),\ldots,F(x_{N-1})$ of the observable $F$ analogously to (3)–(5), replacing throughout integrals with respect to the invariant measure $\mu$ by integrals with respect to the sampling measure $\mu_{N}$. Moreover, because $P_{Q}$ and $P_{Q,N}$ act on different Hilbert spaces, we will approach the problem of comparing their eigenvalues and eigenfunctions through integral operators $P^{\prime\prime}_{Q}:C^{0}(\mathcal{U})\to C^{0}(\mathcal{U})$ and $P^{\prime\prime}_{Q,N}:C^{0}(\mathcal{U})\to C^{0}(\mathcal{U})$, defined analogously to $P_{Q}$ and $P_{Q,N}$, respectively, but acting on the same Banach space of continuous functions on $\mathcal{U}$. A complete description of these constructions will be made in Section 7.

Theorem 5**.**

Let Assumptions 1–4 hold. Then, for any initial point $x_{0}\in\mathcal{V}$,

1. (i)

Every eigenfunction of $P_{Q}$ ($P_{Q,N}$) at nonzero eigenvalue extends to a continuous eigenfunction of $P^{\prime\prime}_{Q}$ ($P^{\prime\prime}_{Q,N}$), corresponding to the same eigenvalue. 2. (ii)

As $N\to\infty$, $P^{\prime\prime}_{Q,N}$ converges in spectrum to $P^{\prime\prime}_{Q}$ in the sense of Corollary 2.

Theorem 5 will be proved in Section 7. Figure 1 shows numerical eigenfunctions of $P_{Q,N}$ obtained from data generated by two mixed-spectrum dynamical systems, described in (40) and (41), respectively. In both examples, we start with a $C^{\infty}$ vector field $\vec{V}$ on a smooth manifold $M$. In the first example, $M=X=\mathbb{T}^{4}$, so $\mathcal{U}=X=M$; in the second example, $M=\mathbb{R}^{3}\times S^{1}$ and $X=X_{\text{Lor}}\times S^{1}\subset M$, where $X_{\text{Lor}}$ is the Lorenz 63 attractor embedded in $\mathbb{R}^{3}$. Eigenfunctions of the operator $P_{Q,N}$ are then computed using a large number of delays, $Q=2000$.

Using the eigenvalues and eigenfunctions of $P_{Q,N}$, we will also construct data-driven Galerkin schemes for the eigenvalue problem of the generator, which are structurally identical to its counterparts formulated in terms of the eigenvalues and eigenfunctions of $P$. Because we do not assume a priori knowledge of the vector field of the dynamics and/or closed-form expressions for the eigenfunctions of $P_{Q,N}$, these schemes will estimate the action of the generator on eigenfunctions through finite-difference approximations at the sampling interval $\Delta t$. In effect, $\Delta t$ will play the role of an additional asymptotic approximation parameter, such that the data-driven solutions converge in a suitable joint limit of vanishing sampling interval ($\Delta t\to 0$), large data ($N\to\infty$), infinitely many delays ($Q\to\infty$), and infinite Galerkin approximation space dimension. This convergence result, along with minimal regularity requirements on the dynamical flow and the kernel, will be stated in a precise manner in Proposition 26 and Assumption 6, respectively. Note that, intuitively, our data-driven Galerkin framework for the generator $V$ requires $\Delta t$ as an additional approximation parameter over methods that approximate the Koopman subgroup generated by $U^{\Delta t}$ at a fixed time step $\Delta t$, since $V$ encodes the information of the entire Koopman group, parameterized by the real time parameter $t$.

Outline of the paper.

In Section 3, we review some important concepts from the spectral theory of dynamical systems. In Section 4, we construct the integral operator $P_{Q}$, which is the key tool of our methods and is also the operator described in Theorems 1–4. Next, we prove these theorems and Corollary 2 in Section 5. In Section 6, we present a Galerkin method for the eigenvalue problem for the Koopman generator, with a small amount of diffusion added for regularization, formulated in the eigenbasis of $P$. In Section 7, we introduce the data-driven realization of $P_{Q}$, and establish the spectral convergence properties stated in Theorem 5, along with the convergence properties of the associated data-driven Galerkin scheme for the generator. In Section 8, the methods are applied to two mixed-spectrum flows, followed by a discussion of the results.

3 Overview of spectral methods for dynamical systems

In this section, we review some concepts from the spectral theory of dynamical systems and establish some facts about Koopman eigenfunctions. Henceforth, we use the notations $\langle f,g\rangle=\int_{X}f^{*}g\,d\mu$ and $\lVert f\rVert=\langle f,f\rangle^{1/2}$ to represent the inner product and norm of $L^{2}(X,\mu)$, respectively.

Generator of a flow.

By continuity of the flow $\Phi^{t}$, the family of operators $U^{t}$ is a strongly continuous, 1-parameter group of unitary transformations of the Hilbert space $L^{2}(X,\mu)$. By Stone’s theorem [61], any such family has a generator $V$, which is a skew-adjoint operator with a dense domain $D(V)\subset L^{2}(X,\mu)$, defined as

[TABLE]

The operators $U^{t}$ and $V$ share the same eigenfunctions; in particular, $z\in D(V)$ with $U^{t}z=e^{i\omega t}z$ satisfies

[TABLE]

In light of (7) and the above relation, we can interpret the quantity $\omega\in\mathbb{R}$ as a frequency intrinsic to the dynamical system (which we sometimes refer to as an “eigenfrequency”).

Vector fields as generators.

If we start with a vector field $\vec{V}$ on a $C^{1}$ manifold $M$, then under appropriate regularity conditions (for example, $\vec{V}$ is locally Lipschitz continuous and satisfies suitable growth bounds at infinity), this vector field induces a $C^{1}$ flow $\Phi^{t}:M\to M$ defined for all $t\in\mathbb{R}$. Suppose that there is a compact invariant set $X\subseteq M$ with an ergodic invariant measure $\mu$. This set $X$ is not necessarily a submanifold, and may not even have any differentiability properties. Nevertheless, $(X,\Phi^{t},\mu)$ is an ergodic dynamical system with an associated strongly-continuous, unitary group of Koopman operators $U^{t}$. Acting on $C^{1}(M)$ functions restricted to $X$, the generator $V$ of this group coincides with the vector field $\vec{V}$, the latter viewed as an operator $\vec{V}:C^{1}(M)\to C^{0}(M)$. For example, in quasiperiodic systems, $X=M=\mathbb{T}^{m}$, $\vec{V}$ generates a rotation, and $\mu$ is equivalent to the Lebesgue volume measure. On the other hand, for the Lorenz attractor (see (39)), $M=\mathbb{R}^{3}$, $\vec{V}$ is smooth and dissipative, $X$ is a compact subset with non-integer fractal dimension [49], and $\mu$ is supported on $X$.

Eigenfunctions as factor maps.

We state the following properties of a Koopman eigenfunction $z$ of an ergodic dynamical system.

If $z$ corresponds to a nonzero eigenfrequency $\omega$, then it has zero mean with respect to the invariant measure $\mu$. This can be concisely expressed as $\langle 1,z\rangle=0$. 2. 2.

The flow $\Phi^{t}$ is semi-conjugate to the irrational rotation by $\omega t$ on the unit circle, with $z$ acting as a semiconjugacy map. This follows directly from (1). Since the eigenfunctions are $L^{2}$ equivalence classes, the semiconjugacy is measure-theoretic (holds $\mu$-a.e.), but would be $C^{r}$ if the eigenfunctions have a $C^{r}$ representation. 3. 3.

Normalized eigenfunctions with $\lVert z\rVert=1$ have $\lvert z(x)\rvert=1$ for $\mu$-a.e. $x\in X$, by (1). As a result, the map $z$ can now be viewed as a projection onto a circle in a measure-theoretic sense, $z(x)\in S^{1}$ for $\mu$-a.e. $x\in X$.

Eigenfunctions form a group.

Another important property of Koopman eigenfunctions for ergodic dynamical systems is that they form a group under multiplication. That is, the product of two eigenfunctions of $U^{t}$ is again an eigenfunction, because of the following relation:

[TABLE]

Moreover, an analogous relation holds for the eigenfunctions and eigenvalues of $V$. The fact that products of Koopman eigenfunctions are Koopman eigenfunctions leads to the following result about products of elements of $\mathcal{D}$ with elements of $\mathcal{D}^{\bot}$.

Lemma 6**.**

Let $\Phi^{t}$ be an ergodic flow on a probability space $(X,\mu)$ such that $U^{t}$ has a mixed spectrum. Then, for every $f\in\mathcal{D}$ and $g\in\mathcal{D}^{\perp}$ for which $fg\in L^{2}(X,\mu)$, $fg$ lies in $\mathcal{D}^{\bot}$.

The eigenvalues of $V$ are closed under integer linear combinations. Moreover, if all corresponding eigenfunctions are continuous, they are generated by a finite set of rationally independent eigenvalues $i\omega_{1},\ldots,i\omega_{m}$. That is, every eigenvalue of $V$ is simple, and has the form $\omega_{\vec{a}}=\sum_{j=1}^{m}a_{j}\omega_{j}$ for some $\vec{a}=(a_{1},\ldots,a_{m})\in\mathbb{Z}^{m}$. Moreover, the corresponding eigenfunction is given by

[TABLE]

where $z_{j}$ is the eigenfunction at eigenvalue $i\omega_{j}$. By virtue of (8) the evolution of every observable $f\in\mathcal{D}$ under $U^{t}$ has the closed-form expression

[TABLE]

which can be evaluated given knowledge of finitely many generating eigenfunctions and eigenfrequencies. The following is a generalization of Property 2 of Koopman eigenfunctions listed above.

Proposition 7**.**

Given an arbitrary collection $\{z_{{\vec{a}}_{1}},z_{{\vec{a}}_{2}},\ldots,z_{{\vec{a}}_{l}}\}$ of $l$ Koopman eigenfunctions, there exists a map $\pi:X\to\mathbb{C}^{l}$ with

[TABLE]

such that:

1. (i)

The image $\pi(X)$ is a torus of dimension $D\leq\min\{m,l\}$, with $D=l$ if $\omega_{{\vec{a}}_{1}},\ldots,\omega_{{\vec{a}}_{l}}$ are rationally independent. 2. (ii)

The flow $(\Phi^{t},\mu)$ on $X$ is semi-conjugate to an ergodic rotation $(\Omega^{t},\operatorname{Leb})$ on $\mathbb{T}^{D}$ (i.e., $\pi\circ\Phi^{t}=\Omega^{t}\circ\pi$, $\mu$-a.e.) associated with a frequency vector whose components are a subset of $\{\omega_{{\vec{a}_{1}}},\ldots,\omega_{{\vec{a}_{l}}}\}$. 3. (iii)

Every Koopman eigenfunction $z$ whose corresponding eigenfrequency is a linear combination of the $\omega_{{\vec{a}}_{1}},\ldots,\omega_{{\vec{a}}_{l}}$ satisfies $z(x)=\zeta(\pi(x))$ for $\mu$-a.e. $x\in X$, where $\zeta\in C^{\infty}(\mathbb{T}^{D})$ is a smooth Koopman eigenfunction of the ergodic rotation on the $D$-torus corresponding to the same eigenfrequency.

Remark**.**

If $m>1$, the set of eigenvalues $\{i\omega_{\vec{a}}\}_{\vec{a}\in\mathbb{Z}^{m}}$ is dense on the imaginary axis. This property adversely affects the stability of numerical approximations of Koopman eigenvalues and eigenfunctions even in systems with pure point spectrum, necessitating the use of regularization [33]. We will return to this point in Section 6.

Lemma 8** ([51], Section 2.3).**

Let $\Delta t>0$ be as in Assumption 1. Then, the orthogonal projection $\pi_{\omega}f$ of an observable $f\in L^{2}(X,\mu)$ onto the eigenspace of $U^{\Delta t}$ corresponding to the eigenvalue $e^{i\omega\,\Delta t}$ of $U^{\Delta t}$ is given by

[TABLE]

Moreover, $\pi_{\omega}\equiv 0$ if $i\omega$ is not an eigenvalue of the generator. Otherwise, $U^{\Delta t}\pi_{\omega}f=e^{i\omega\,\Delta t}\pi_{\omega}f$.

Mixing and weak mixing.

An observable $f\in L^{2}(X,\mu)$ is said to be mixing if for all $g\in L^{2}(X,\mu)$, $\lim_{t\to\infty}\langle g,U^{t}f\rangle=0$; it is said to be weak-mixing if $\lim_{t\to\infty}t^{-1}\int_{0}^{T}\lvert\langle g,U^{s}f\rangle\rvert\,ds=0$. The latter, is equivalent to the requirement that for Lebesgue almost every $\Delta t\in\mathbb{R}$, $\lim_{N\to\infty}N^{-1}\sum_{n=0}^{N-1}\lvert\langle g,U^{n\,\Delta t}f\rangle\rvert=0$. The flow $\Phi^{t}$ is said to be (weak-) mixing if $f$ is (weak-) mixing for all $f\in L^{2}(X,\mu)$. It is known that every $f\in\mathcal{D}^{\bot}$ is weak mixing (see, e.g., Mixing Theorem, p. 45 in [37]), whereas no observable in $\mathcal{D}$ is weak-mixing. Thus, the component $\mathcal{D}$, often called the quasiperiodic subspace, shows no decay of correlation, unlike its complement $\mathcal{D}^{\bot}$, which represents the chaotic component of the dynamics. In addition, weak-mixing observables in $\mathcal{D}^{\perp}$ and observables in $\mathcal{D}$ have a useful pointwise decorrelation property:

Lemma 9**.**

Let $f\in\mathcal{D}^{\perp}$ and $g\in\mathcal{D}$. Then, for $\mu$-a.e. $x,y\in X$,

[TABLE]

Proof.

Without loss of generality, we may assume that $g$ is an eigenfunction of $U^{\Delta t}$ with eigenvalue $e^{i\omega\,\Delta t}$. Then,

[TABLE]

which is equal to $g^{*}(x)\pi_{\omega}f(y)$ by Lemma 8. The latter is equal to zero since $f\in\mathcal{D}^{\bot}$. ∎∎

4 Kernel integral operators from delay-coordinate mapped data

4.1 Choice of kernel

Consider a kernel integral operator of the class (3) associated with an $L^{2}$ kernel $k:M\times M\to\mathbb{R}$. Then, the following properties hold (e.g., [26, 25]):

$K$ is a Hilbert-Schmidt, and therefore compact, operator on $L^{2}(X,\mu)$, with operator norm bounded by $\lVert k\rVert_{L^{2}(X\times X)}$. 2. 2.

If $k$ is symmetric, then $K$ is self-adjoint. 3. 3.

If $k$ is $C^{0}$, then $Kf$ is also $C^{0}$ for every $f\in L^{2}(X,\mu)$. 4. 4.

If $M$ is a $C^{r}$ manifold and $k$ is $C^{r}$, then $Kf$ is also $C^{r}$ for every $f\in L^{2}(X,\mu)$.

As stated in Section 2, we will work with kernels of the form

[TABLE]

where $h$ is a continuous shape function on $\mathbb{R}$, and $d_{Q}:M\times M\to\mathbb{R}_{0}$ is the distance-like function on $M$ from (4), parameterized by the number of delays $Q$. Kernels of this class are sometimes referred to as stationary kernels [31], as they only depend on distances between datapoints. For example, in (5), we used a Gaussian shape function, which is popular in manifold learning and other related geometrical data analysis techniques. Note that $d_{Q}$ is symmetric, non-negative, and satisfies the triangle inequality, but depending on the properties of $F$ and the number of delays it may vanish on distinct points. That is, $d_{Q}$ is a pseudo-distance on $M$, induced from delay-coordinate mapped data with $Q$ delays.

The kernels in (10) satisfy Assumption 2(i), and the associated kernel integral operators $K_{Q}$ have all four properties listed above. In addition, if $h$ is strictly positive, $k_{Q}$ satisfies Assumption 2(ii). The behavior of integral operators associated with other classes of kernels, e.g., the covariance operators employed in SSA and Hankel matrix analysis induced by inner products in data space, can be studied via similar techniques to those presented below. However, it should be kept in mind that the Markov normalization procedure described in Section 4.3 (which will be important for the well-posedness of the Galerkin schemes in Sections 6 and 7) requires that the kernel be sign-definite. Another consideration to keep in mind is that the ability to approximate Koopman eigenfunctions with our techniques depends on the “richness” of the range of $K_{Q}$. As can be readily verified, the operator $K_{Q}$ constructed from covariance kernels in $d$-dimensional data space (as in Assumption 1) has at most a $dQ$-dimensional range, whereas the corresponding operators associated with Gaussian kernels, as well as other non-polynomial kernels, have typically infinite-dimensional range for any $Q$. Our approach should also be applicable with little modification to families of kernels of the form

[TABLE]

where averaging takes place after application of the shape function. Lemma 10 below states some useful properties of $K_{Q}$ associated with strictly positive kernels. In what follows, $1_{S}$ will denote the constant function equal to 1 on a set $S$.

Lemma 10**.**

Under Assumptions 1 and 2(ii), for any $Q\in\mathbb{N}$, the functions $\rho_{Q}=K_{Q}1_{X}$, and $\sigma_{Q}=K_{Q}\left(1/\rho_{Q}\right)$ are continuous and positive. Moreover, restricted on $X$, they are bounded away from zero.

Proof.

The claims follow directly by compactness of $X$ and the fact that $k_{Q}\rvert_{X\times X}$ is a continuous function, bounded away from zero. ∎∎

Intuitively, $\rho_{Q}$ can be thought of as a “sampling density” on $X$. For instance, if $X$ were a manifold embedded in $\mathbb{R}^{Qd}$ by a delay-coordinate map constructed from $F$, then up to an $\epsilon$-dependent scaling, $\rho_{Q}$ would approximate the density of the invariant measure $\mu$ relative to the volume measure associated with that embedding.

Remark**.**

In a number of applications, such as statistical learning on manifolds [5, 16, 11, 9], one-parameter families of integral operators such as $K_{Q}$ and $P_{Q}$ are studied in the limit $\epsilon\to 0$, where under certain conditions they can be used to approximate generators of Markov semigroups; one of the primary examples being the Laplace-Beltrami operator on Riemannian manifolds. Here, the fact that the state space $X$ may not (and in general, will not) be smooth precludes us from taking such limits unconditionally. However, according to Theorem 3(ii), passing first to the limit $Q\to\infty$ allows one to view $K$ and $P$ as operators on functions on a smooth manifold, namely a $D$-dimensional torus, and study the small-$\epsilon$ behavior of these operators in that setting.

4.2 Asymptotic behavior in the infinite-delay limit

To study the behavior of $K_{Q}$ in the limit of infinitely many delays, $Q\to\infty$, we first consider the properties of the pseudometric $d_{Q}$ in the same limit. The latter can be studied in turn through a useful (nonlinear) map $\Psi:C^{0}(X)\to L^{\infty}(X\times X,\mu\times\mu)$, which maps a given observation function $F$ into a (pseudo)metric on $X$, namely,

[TABLE]

In what follows, $d_{X}:X\times X\to\mathbb{R}$ will denote the metric $X$ inherits from $M$.

Theorem 11**.**

Let Assumption 1 hold, and $F=F_{\mathcal{D}}+F_{\mathcal{D}^{\bot}}$ be the $L^{2}$ decomposition of $F$ from (2). Then, $\Psi(F)$ in (11) is well-defined as a function in $L^{\infty}(X\times X,\mu\times\mu)$, and $\Psi_{Q}(F)$ converges to $\Psi(F)$ in $L^{p}(X\times X,\mu\times\mu)$ norm for $1\leq p<\infty$. Moreover:

1. (i)

For every $t\in\mathbb{R}$ and $\mu$-a.e. $x,y\in X$, $\Psi(F)(\Phi^{t}(x),\Phi^{t}(y))=\Psi(F)(x,y)$. 2. (ii)

For $\mu$-a.e. $x,y\in X$, $\Psi(F)(x,y)=\Psi(F_{\mathcal{D}^{\bot}})(x,y)+\Psi(F_{\mathcal{D}})(x,y)$. 3. (iii)

$\Psi(F_{\mathcal{D}^{\bot}})$* is a constant almost everywhere and equals $2\|F_{\mathcal{D}^{\bot}}\|_{L^{2}}^{2}$. Therefore,*

[TABLE]

In particular, $\Psi(F)\in\mathcal{D}\times\mathcal{D}$.

If, moreover, Assumption 3 holds:

1. (iv)

$\Psi(F_{\mathcal{D}})\in C^{0}(X\times X)$* and $\Psi_{Q}(F_{\mathcal{D}})$ converges to $\Psi(F_{\mathcal{D}})$ uniformly on $X\times X$.* 2. (v)

$\Psi(F)$* is uniformly continuous on a full-measure, dense subset of $X\times X$.* 3. (vi)

$\Psi(F)$* has a unique continuous extension $\bar{\Psi}(F)\in C^{0}(X\times X)$, and $\Psi_{Q}(F)$ converges to $\bar{\Psi}(F)$ $\mu$-almost uniformly.*

Proof.

To prove well-definition of $\Psi$, note that $\Psi(F)$ exists $\mu$-a.e. since it is the pointwise limit of the Birkhoff averages $\Psi_{Q}(F)$ of the continuous function $d_{1}:(x,y)\to\lVert F(x)-F(y)\rVert$ with respect to the product flow $\Phi^{t}\times\Phi^{t}$ on $X\times X$. By compactness of $X\times X$, each of the functions $\Psi_{Q}(F)$ is bounded above by $\lVert d_{1}\rVert_{C^{0}(X\times X)}$. Therefore, $\Psi(F)$ lies in $L^{\infty}(X\times X,\mu\times\mu)$, and thus in $L^{p}(X\times X,\mu\times\mu)$, $1\leq p<\infty$, since $\mu\times\mu$ is a probability measure. The $\Psi_{Q}(F)\to\Psi(F)$ convergence in $L^{p}(X\times X,\mu\times\mu)$, 1 $\leq p<\infty$, then follows from the $L^{p}$ Von Neumann ergodic theorem.

By the invariance of the infinite Birkhoff averages, $\Psi(F)$ is invariant under the flow $\Phi^{\Delta t}\times\Phi^{\Delta t}$. Thus $\Psi(F)$ must lie in the kernel of $V\otimes V$ and thus is invariant under the flow $\Phi^{t}\times\Phi^{t}$ for all $t\in\mathbb{R}$, proving Claim (i).

To prove Claim (ii), let $x_{q}$ and $y_{q}$ denote $\Phi^{q\,\Delta t}(x)$ and $\Phi^{q\,\Delta t}(y)$ respectively. Let $G_{\mathcal{D}}:X\times X\to\mathbb{R}^{d}$:= $(x,y)\to F_{\mathcal{D}}(x_{q})-F_{\mathcal{D}}(y_{q})$, and similarly define $G_{\mathcal{D}^{\bot}}:X\times X\to\mathbb{R}^{d}$. Expanding the right-hand side of (11) gives,

[TABLE]

and the first two terms in the equation above are $\Psi(F_{\mathcal{D}})(x,y)$ and $\Psi(F_{\mathcal{D}^{\bot}})(x,y)$ respectively. Therefore, to prove Claim (ii), it suffices to prove that the third term vanishes. This is equivalent to showing that for $\mu$-a.e. $x,y\in X$,

[TABLE]

which follows from Lemma 9. This completes the proof of Claim (ii).

To prove Claim (iii), let $x_{n}$ and $y_{n}$ denote $\Phi^{n\,\Delta t}(x)$ and $\Phi^{n\,\Delta t}(y)$, respectively. Then, (11) can be rewritten for $F_{\mathcal{D}^{\bot}}$ as

[TABLE]

The first two terms converge to the constant $\|F_{\mathcal{D}^{\bot}}\|_{L^{2}}^{2}$. It is therefore sufficient to show that the last term vanishes. Indeed, since the function $J:(x,y)\to F_{\mathcal{D}^{\bot}}(x)F_{\mathcal{D}^{\bot}}(y)$ lies in the continuous spectrum subspace of the product-system $(X\times X,\Phi^{t}\times\Phi^{t},\mu\times\mu)$, we have

[TABLE]

Since $F_{\mathcal{D}}$ is continuous, $\Psi(F_{\mathcal{D}})$ is continuous by a classic result of Krengel ([42], Theorem 1.2.7). This proves Claim (iv).

Turning to Claim (v), it follows directly from Claims (iii) and (iv) that there exists a full-measure subset $S\subseteq X\times X$ on which $k_{\infty}$ is uniformly continuous. Suppose that $S$ were not dense in $X\times X$. Then, there would exist an open set $B\subset X\times X$ disjoint from $S$, and with positive measure (since $X\times X$ is the support of $\mu\times\mu$, and every open subset of the support of a Borel measure has positive measure), which would in turn imply that $(\mu\times\mu)(S)<1$, leading to a contradiction. Therefore, $S$ is a full-measure, dense subset of $X\times X$, completing the proof of the claim.

Finally, the existence of $\bar{\Psi}(F)$ in Claim (vi) follows from the fact that $\Psi(F)$ is uniformly continuous on the dense subset $S$ of the compact metric space $X\times X$, and the almost uniform convergence of $\Psi_{Q}(F)$ to $\bar{\Psi}(F)$ is a consequence of Egorov’s theorem. ∎∎

Remark**.**

*Although the measure $\mu\times\mu$ is invariant under $\Phi^{t}\times\Phi^{t}$, it is not ergodic. In fact, it is ergodic iff $(\Phi^{t},\mu)$ is mixing (equivalently, $U^{t}$ has purely continuous spectrum and a simple eigenvalue at 1), in which case the metric $d_{\infty}$ would be constant almost everywhere, in accordance with (12). *

Theorem 11 establishes that the function $d_{\infty}:D(d_{\infty})\to\mathbb{R}$, such that

[TABLE]

is well-defined as a function in $L^{p}(X\times X,\mu\times\mu)$, $1\leq p\leq\infty$, with $\sup d_{\infty}\leq\lVert d_{1}\rVert_{C^{0}(X\times X)}$. It can also be verified that $d_{\infty}$ satisfies the triangle inequality and is non-negative. However, depending on the properties of the dynamical system and observation map, it may be a degenerate metric as $d_{\infty}(x,y)$ may vanish for some $x\neq y$, even if $d_{Q}(x,y)$ is non-vanishing. In fact, it is easy to check that if $y$ lies in the stable manifold of $x$, then $d_{\infty}(x,y)=0$. Analogously to the finite-delay case in (10), we employ $d_{\infty}$ and the shape function $h$ to define a corresponding kernel $k_{\infty}:M\times M\to\mathbb{R}$, where

[TABLE]

and $k_{\infty}(x,y)=0$ otherwise. We also let $K$ be the kernel integral operator from (3) associated with $k_{\infty}$.

Proposition 12 shows that the operator $K$ depends only on the quasiperiodic component of $F$, and is a direct consequence of Theorem 11 and (12).

Proposition 12**.**

Let $(X,\Phi^{t},\mu)$ and $F$ be as in Theorem 1. Then, the integral operator $K$ is a constant scaling operator iff its kernel $k_{\infty}$ is a constant $\mu$-a.e., which occurs iff $F_{\mathcal{D}}$ is a constant.

In general, $k_{\infty}$ may not be continuous. Nevertheless, it has a number of other useful properties, which follow directly from Theorem 11 in conjunction with the boundedness and continuity of the Gaussian shape function.

Lemma 13**.**

Under Assumption 1, the following hold:

1. (i)

$k_{\infty}$* is the $L^{p}(X,\mu)$-norm limit, $1\leq p<\infty$, of the sequence of continuous kernels $k_{1},k_{2},\ldots$.* 2. (ii)

$k_{\infty}$* is invariant under $U^{t}\times U^{t}$ for all $t\in\mathbb{R}$.* 3. (iii)

$k_{\infty}$* lies in $L^{\infty}(X\times X,\mu\times\mu)$, and under Assumption 2(ii), $1/k_{\infty}$ also lies in that space.

Moreover, if Assumption 3 additionally holds:

1. (iv)

$k_{\infty}$* is uniformly continuous on a dense, full-measure subset of $X\times X$.* 2. (v)

$k_{\infty}$* has a unique continuous representative $\bar{k}_{\infty}\in C^{0}(X\times X)$, and as $Q\to\infty$, $k_{Q}$ converges to $k_{\infty}$ almost uniformly.*

The stronger regularity properties of $k_{\infty}$ under Assumption 3 have the following important implications on the behavior of the corresponding integral operator.

Lemma 14**.**

Under Assumptions 1 and 3, the kernel integral operator $K$ associated with $k_{\infty}$ has the following properties:

1. (i)

For every $f\in L^{2}(X,\mu)$, $Kf$ has a unique continuous representative. 2. (ii)

For every $f\in C^{0}(X)$, $Kf$ is continuous. 3. (iii)

$\lVert K\rVert\leq\lVert k_{\infty}\rVert_{L^{\infty}(X\times X)}$* in either $L^{2}$ or $C^{0}$ operator norm.* 4. (iv)

As an operator on $C^{0}(X)$, $K$ is compact. 5. (v)

For every $f\in C^{0}(X)$, $K_{Q}f$ is a sequence of continuous functions converging $\mu$-a.e. to $Kf$.

Proof.

(i) Since $k_{\infty}$ is uniformly continuous on a set $S\subseteq X\times X$ of full $\mu\times\mu$ measure, there exists a full $\mu$-measure set $X^{\prime}\subseteq X$, such that for every $x\in X^{\prime}$, $k_{\infty}(x,\cdot)$ is continuous $\mu$-a.e. on $X$. Moreover, proceeding analogously to the proof of Theorem 11(v), it can be shown that $X^{\prime}$ is dense in $X$. Let now $f\in L^{2}(X,\mu)$, $\|f\|_{L^{2}}=1$. Then, for every $x_{1},x_{2}\in X^{\prime}$,

[TABLE]

Since $k_{\infty}$ is uniformly continuous on $S$, for every $\epsilon>0$, there exists $\delta>0$ such that if $d_{X}(x_{1},x_{2})<\delta$, $\|k_{\infty}(x_{1},\cdot)-k_{\infty}(x_{2},\cdot)\|_{L^{\infty}}<\epsilon$. Thus, for all such $x_{1}$ and $x_{2}$, we have $\left|Kf(x_{1})-Kf(x_{2})\right|<\epsilon$, which implies that $Kf$, restricted to $X^{\prime}$, is uniformly continuous. As a result, since $X^{\prime}$ is dense in the compact metric space $X$, $Kf\rvert_{X^{\prime}}$ has a unique continuous extension $g\in C^{0}(X)$. Moreover, since $X^{\prime}$ has full measure, $g$ lies in the same $L^{2}$ equivalence class as $Kf$, proving the claim.

(ii) Since $k_{\infty}$ is uniformly continuous on a dense set of full measure, for any $f\in C^{0}(X)$, the function $g:X\times X\to\mathbb{C}$ with $g(x,y)=k_{\infty}(x,y)f(y)$ has a unique continuous representative $\bar{g}\in C^{0}(X\times X)$. Therefore, for every $x\in X$, the function $k_{\infty}(x,\cdot)f$ is $\mu$-a.e. equal to $\bar{g}(x,\cdot)$ by $\mu$-a.e. continuity of $k_{\infty}(x,\cdot)$, and

[TABLE]

It then follows that $Kf$ is continuous by continuity of integrals of $X$-sections of continuous functions on $X\times X$.

(iii) To verify the claim on the $L^{2}$ and $C^{0}$ operator norms, observe that for every $f\in L^{2}(X,\mu)$ and $x\in X^{\prime}$, where $X^{\prime}$ is as in the proof of Claim (i),

[TABLE]

and therefore

[TABLE]

The bound on the $L^{2}$ operator norm follows by setting $\|f\rVert_{L^{2}}=1$ in (15), together with the fact that $\lVert Kf\rVert_{L^{2}}\leq\lVert Kf\rVert_{L^{\infty}}$. The bound on the $C^{0}$ operator norm follows from (15) with $f\in C^{0}(X)$, in conjunction with the facts that $\lVert f\rVert_{L^{2}}\leq\lVert f\rVert_{C^{0}}$ and $\lVert Kf\rVert_{L^{\infty}}=\lVert Kf\rVert_{C^{0}}$.

(iv) Since, by the Arzelà-Ascoli theorem, every equicontinuous sequence of functions on a compact metric space has a limit point, it suffices to show that for every sequence $f_{n}\in C^{0}(X)$ with $\lVert f_{n}\rVert_{C^{0}}\leq 1$, the sequence $g_{n}=Kf_{n}$ has a limit point with respect to $C^{0}$ norm. Let $\bar{k}_{\infty}\in C^{0}(X\times X)$ be the unique continuous representative of $k_{\infty}$. For every $x_{1},x_{2}\in X$, we have

[TABLE]

and by uniform continuity of $\bar{k}_{\infty}$, for any $\epsilon>0$, there exists $\delta>0$, independent of $n$, such that, for every $x_{1},x_{2}\in X$ with $d(x_{1},x_{2})<\delta$, $\lvert g_{n}(x_{1})-g_{n}(x_{2})\rvert<\epsilon$. This establishes equicontinuity of $g_{n}$, and thus compactness of $K$ on $C^{0}(X)$.

(v) The continuity of $K_{Q}f$ and $Kf$ follows from Claim (ii). The $\mu$-a.e. convergence follows from Lemma 13(v). ∎

We end this section with two important corollaries of Theorem 11 and Lemmas 13, 14, which are central to both Theorems 1 and 3.

Corollary 15**.**

The operators $U^{t}$ and $K$ commute.

Proof.

Since $\mu$ is an invariant measure, for every $x$ in $X$ and $t\in\mathbb{R}$ we have

[TABLE]

It therefore follows from Lemma 13(ii) that

[TABLE]

and the claim of the corollary follows. ∎∎

Corollary 16**.**

Under Assumptions 1 and 2(ii), the function $\rho=K1_{X}$ is $\mu$-a.e. equal to a constant bounded away from zero (i.e., $1/\rho$ lies in $L^{\infty}(X,\mu)$). Further, if Assumption 3 holds, then $\rho\rvert_{X}$ and $1/\rho\rvert_{X}$ are continuous.

Proof.

Corollary 15 and the fact that $U^{t}1_{X}=1_{X}$ imply that $U^{t}\rho=\rho$, and it then follows by ergodicity that $\rho$ is constant $\mu$-a.e. That $\lVert 1/\rho\rVert_{L^{\infty}}$ is finite follows from Lemma 13(iii). Finally, the continuity of $\rho$ under Assumption 3 is a direct consequence of Lemma 14. ∎∎

4.3 Markov normalization

Next, we construct the Markov operators $P_{Q}$ and $P$ appearing in Theorems 1 and 3 by normalization of $K_{Q}$ and $K$. Throughout this section, we consider that Assumptions 1 and 2 hold. Under these assumptions, we employ a normalization procedure introduced in the diffusion maps algorithm [16] and further developed in [11], although there are also other approaches with the same asymptotic behavior. Specifically, using the normalizing functions $\rho_{Q}$ and $\sigma_{Q}$ from Lemma 10 and $\rho$ from Corollary 16, we introduce the kernels $p_{Q}:M\times M\to\mathbb{R}$ and $p:M\times M\to\mathbb{R}$, given by

[TABLE]

respectively. By Lemma 10, $p_{Q}$ satisfies the boundedness and continuity properties in Assumption 2. On the other hand, $p$ is neither guaranteed to be continuous nor bounded on arbitrary compact sets, but it nevertheless follows from Lemma 13 and Corollary 16 that both $p$ and $1/p$ lie in $L^{\infty}(X\times X)$. Based on these facts, we can therefore define the kernel integral operators $P_{Q}:L^{2}(X,\mu)\to L^{2}(X,\mu)$ and $P:L^{2}(X,\mu)\to L^{2}(X,\mu)$ from (3) associated with the kernels $p_{Q}$ and $p$, respectively, and these operators are both Hilbert-Schmidt (see Section 4.1). Note that $p$ and $P$ have analogous properties to those stated for $k_{\infty}$ and $K$ in Lemmas 13, 14 and Corollary 15. In particular, $p$ is invariant under $U^{t}\times U^{t}$, and $P$ commutes with $U^{t}$.

The operators $P_{Q}$ and $P$ can also be obtained directly from $K_{Q}$ and $K$, respectively, through the sequence of operations

[TABLE]

In [11], the steps leading to $\tilde{K}_{Q}$ from $K_{Q}$ and to $P_{Q}$ from $\tilde{K}_{Q}$ are called right and left normalization, respectively. In the case of $P$, the effects of right normalization cancel since $K1_{X}$ is $\mu$-a.e. constant by Corollary 16, so it is sufficient to construct this operator directly via left normalization of $K$.

As is evident from (17), $P_{Q}$ and $P$ are both Markov operators preserving constant functions. Moreover, for all $x\in M$ we have $\int_{X}p_{Q}(x,\cdot)\,d\mu=1$, and for $\mu$-a.e. $x\in M$, $\int_{X}p(x,\cdot)\,d\mu=1$, i.e., both $p_{Q}$ and $p$ are transition probability kernels. In particular, since $X$ is compact and $p_{Q}$ and $p$ are essentially bounded below, $P_{Q}$ and $P$ are both ergodic Markov operators; that is, their eigenspaces at eigenvalue 1 are one-dimensional.

The Markov kernel $p$ is $\mu$-a.e. symmetric by symmetry of $k_{\infty}$ and the fact that $\rho$ is $\mu$-a.e. constant. As a result, $P$ is self-adjoint, its eigenvalues admit the ordering $1=\lambda_{0}>\lambda_{1}\geq\lambda_{2}\geq\cdots$, and there exists a real orthonormal basis $L^{2}(X,\mu)$ consisting of corresponding eigenfunctions, $\phi_{j}$, with $\phi_{0}$ being constant. On the other hand, because $p_{Q}$ is not symmetric, the operator $P_{Q}$ is not self-adjoint, but is nevertheless related to a self-adjoint operator via a similarity transformation by a bounded multiplication operator with a bounded inverse. To verify this, define

[TABLE]

where $\rho_{Q}$ and $\sigma_{Q}$ are as in Lemma 10. Let also $D_{Q}$ be the multiplication operator which multiplies by $\tilde{\sigma}_{Q}$, and $\hat{P}_{Q}$ the kernel integral operator with kernel $\hat{p}_{Q}:M\times M\to\mathbb{R}$

[TABLE]

Observe now that $\hat{P}_{Q}$ is a symmetric operator, and $P_{Q}$ is related to $\hat{P}_{Q}$ via the similarity transformation

[TABLE]

that is, for every $f\in L^{2}(X,\mu)$,

[TABLE]

The following are useful properties of $\hat{P}_{Q}$ that follow from its relation to $P_{Q}$.

$\hat{P}_{Q}$ has the same discrete spectrum as $P_{Q}$, consisting of eigenvalues $\lambda_{j,Q}$ with $1=\lambda_{0,Q}>\lambda_{1,Q}\geq\lambda_{2,Q}\geq\cdots$. 2. 2.

Let $\phi_{j,Q}$ denote the eigenfunctions of $\hat{P}_{Q}$ corresponding to the nonzero eigenvalues $\lambda_{j,Q}$. These form an orthonormal basis for the closed subspace $\overline{\operatorname{ran}\hat{P}_{Q}}=(\ker\hat{P}_{Q})^{\bot}$. Moreover, the $\phi_{j,Q}$ can be chosen to be real-valued. 3. 3.

The eigenfunction $\phi_{0,Q}$ of $\hat{P}_{Q}$ is equal up to proportionality constant to $\rho_{Q}\sigma^{1/2}_{Q}$.

Remark**.**

In applications, it may be the case that $\rho_{Q}$ and $1/\rho_{Q}$ take a large range of values. In such situations, it may be warranted to replace (4) by a variable-bandwidth kernel of the form $k_{Q}(x,y)=\exp\left(-\frac{d^{2}_{Q}(x,y)}{\epsilon r_{Q}(x)r_{Q}(y)}\right)$, with a bandwidth function $r_{Q}$ introduced so as to control the decay of the kernel away from the diagonal, $x=y$. Various types of bandwidth functions have been proposed in the literature, including functions based on neighborhood distances [69, 10], state space velocities [36, 32], and local density estimates [9]. While we do not study variable bandwidth techniques in this work, our approach should be applicable in that setting too, so long as Corollary 16 holds.

5 Proof of Theorems 1–4 and Corollary 2

Proof of Theorem 1.

That $P$ and $U^{t}$ commute follows from the invariance of $p$ under $U^{t}\times U^{t}$ and an analogous calculation to that in the proof of Corollary 15. Next, as $Q\to\infty$, $p_{Q}$ converges to $p$ in any $L^{p}(X\times X,\mu\times\mu)$ norm with $1\leq p<\infty$ by the analogous result to Lemma 13(i) that holds for these kernels (see Section 4.3). In particular, that $p_{Q}$ converges to $p$ in $L^{2}(X\times X,\mu\times\mu)$ norm implies that $P_{Q}$ converges to $P$ in $L^{2}(X,\mu)$ operator norm, since $P_{Q}-P$ is Hilbert-Schmidt and thus bounded in operator norm by $\lVert p_{Q}-p\rVert_{L^{2}(X\times X)}$. ∎

Proof of Theorem 3.

We first establish that $\tau$ is a.e. invariant under $\Phi^{t}\times\Phi^{t}$. Since the integral operator $T$ commutes with $U^{t}$, for $\mu$-a.e. $x\in X$,

[TABLE]

where the second equality was obtained by the change of variables $y=\Phi^{t}(y^{\prime})$, and utilizes the invariance of the measure $\mu$ under $\Phi^{t}$. The only way the terms at the two ends of the equation can be equal for $\mu$-a.e. $x\in X$ is if $\tau(\Phi^{t}(x),\Phi^{t}(y^{\prime}))$ = $\tau(x,y^{\prime})$ $\mu$-a.e.

Next, observe that, by (2), the space $L^{2}(X\times X,\mu\times\mu)$ splits as the $U^{t}\times U^{t}$-invariant orthogonal sum of $\mathcal{D}\otimes\mathcal{D}$, $\mathcal{D}^{\bot}\otimes\mathcal{D}^{\bot}$, $\mathcal{D}^{\bot}\otimes\mathcal{D}$, and $\mathcal{D}\otimes\mathcal{D}^{\bot}$. Since $\tau$ is an $L^{2}$ kernel, it has orthogonal projections onto each of these subspaces, all of which are $U^{t}\otimes U^{t}$-invariant by the invariance of $\tau$ just established. By symmetry of $\tau$, the projections onto $\mathcal{D}^{\bot}\otimes\mathcal{D}$ and $\mathcal{D}\otimes\mathcal{D}^{\bot}$ vanish. Moreover, the projection $\tau_{\mathcal{D}^{\perp}\otimes\mathcal{D}^{\perp}}\in\mathcal{D}^{\perp}\otimes\mathcal{D}^{\perp}$ is orthogonal to constant functions, and it follows by the Birkhoff ergodic theorem that for $\mu\times\mu$-a.e. $x,y\in X\times X$,

[TABLE]

This completes the proof of Claim (i). The statements in Claim (ii) that $\mathcal{D}^{\bot}\subset\ker(T)$ and that $\mathcal{D}$ and $\mathcal{D}^{\perp}$ are invariant under $T$ are direct consequences of Claim (i).

The remaining two claims in the theorem, which requires that both $\mathcal{D}$ and $\operatorname{ran}T$ contain non-constant functions, can be proved by means of the following, slightly stronger, result.

Proposition 17**.**

For any nonzero eigenvalue $\lambda$ of $T$, the corresponding eigenspace $W_{\lambda}$ is invariant under the action of the Koopman generator $V$, and $V\rvert_{W_{\lambda}}$ is diagonalizable. Moreover, the constant function $1_{X}$ is an eigenfunction of $T$. If $W_{\lambda}$ does not contain $1_{X}$, its dimension is an even number.

Proof.

Since $T$ is compact, every nonzero eigenvalue $\lambda$ has finite multiplicity and its corresponding eigenspace $W_{\lambda}$ has finite dimension, $l=\dim W_{\lambda}$. Since $U^{t}$ commutes with $T$, $U^{t}$ and hence $V$ leave $W_{\lambda}$ invariant. Similarly, since the constant function is an eigenfunction of $V$, it is an eigenfunction of $T$.

Let $\lambda_{0}$ be the eigenvalue of $T$ corresponding to the constant eigenfunction, and $\lambda\neq\lambda_{0}$ be any other eigenvalue of $T$. Then, $V\rvert_{W_{\lambda}}$ is a skew-symmetric operator on a finite-dimensional space, and thus can be diagonalized with respect to a basis of simultaneous eigenfunctions of $T$ and $V$. Fix any element $\zeta$ of this basis. By our choice of $\lambda$, $\zeta$ is a non-constant eigenfunction of $V$, hence $\langle\zeta,1\rangle=0$. Therefore, by ergodicity of $(\Phi^{t},\mu)$, $V\zeta=i\omega\zeta$ for some $\omega\neq 0$. This implies that $\zeta$ has non-zero real and imaginary parts. Hence, the conjugate $\zeta^{*}$ is linearly independent from $\zeta$ and corresponds to eigenvalue $-i\omega$ of $V$. However, since $T$ is a real operator, $\zeta^{*}$ lies in $W_{\lambda}$. We therefore conclude that $W_{\lambda}$ can be split into disjoint 2-dimensional spaces spanned by the conjugate pair of eigenfunctions $\zeta$ and $\zeta^{*}$. Therefore $\dim W_{\lambda}$ is an even number. ∎∎

Corollary 18**.**

The representation of $V\rvert_{\overline{\operatorname{ran}P}}$ in the basis $\{\phi_{0},\phi_{1},\ldots\}$ has a block-diagonal structure, consisting of even-sized blocks associated with the eigenspaces $W_{\lambda\neq 1}$, and a $1\times 1$ block with the element 0, associated with $W_{1}$. Moreover, the range of $P$ lies in the domain of $V$, and $V\rvert_{\overline{\operatorname{ran}P}}$ and $P\rvert_{\overline{\operatorname{ran}P}}$ are simultaneously diagonalizable.

By Proposition 17, $U^{t}$ and $T$ have joint eigenfunctions, each of which factors the dynamics into a rotation on the circle in accordance with (1). According to Proposition 7, any collection of $D$ such eigenfunctions factors the dynamics into a rotation on $\mathbb{T}^{D}$. This proves Claim (iii).

To prove Claim (iv), we use (8) to expand the kernel as

[TABLE]

In this expansion, there is a minimal number $D\leq m$ of generating eigenfunctions $z_{j}$ from (8), arranged without loss of generality as $z_{1},\ldots,z_{D}$, such that the expansion coefficients $\tilde{\tau}_{\vec{a}\vec{b}}$ corresponding to $\vec{a}=(a_{1},\ldots,a_{m})$ and $\vec{b}=(b_{1},\ldots,b_{m})$ with nonzero $a_{D+1},\ldots,a_{m}$ and $b_{D+1},\ldots,b_{m}$, respectively, vanish (in other words, the kernel $\tau$ does does not project onto the subspaces generated by $z_{D+1},\ldots,z_{m}$ and their powers). By Proposition 7, the Koopman eigenfunctions corresponding to non-vanishing $\tilde{\tau}_{\vec{a}\vec{b}}$ can be expressed as $z_{\vec{a}}=\zeta_{\vec{a}}\circ\pi$, where the $\zeta_{\vec{a}}$ are smooth Koopman eigenfunctions on $\mathbb{T}^{D}$ associated with an ergodic rotation. Thus, denoting the index set for the nonzero $\tilde{\tau}_{\vec{a}\vec{b}}$ by coefficients by $I\in\mathbb{Z}^{m}\times\mathbb{Z}^{m}$, we have $\tau(x,y)=\hat{\tau}(\pi(x),\pi(y))$ for $\mu\times\mu$-a.e., $(x,y)\in X\times X$, where $\hat{\tau}$ is the $L^{2}$ kernel on $\mathbb{T}^{D}$ given by

[TABLE]

This completes the proof of Claim (v) and of Theorem 3. ∎

Proof of Theorem 4.

That $p$ is uniformly continuous on a full-measure, dense subset of $X\times X$ follows from the analogous result to Lemma 13(iv), which holds for $p$ (see Section 4.3). Claims (i)–(iv) of the theorem follow analogously to Lemma 14. ∎

Rates of convergence in the continuous case.

As an auxiliary result, we state a lemma that establishes rates of convergence with respect to the number of delays $Q$ of the kernel integral operators studied in this work.

Lemma 19** (Convergence of commutators).**

Let the assumptions of Theorem 4 hold, and the shape function $h$ from (10) be continuously differentiable. Then, the following operators converge in $C^{0}(X)$ operator norm to [math] as $Q\to\infty$, with rates given below:

1. (i)

$\left\|U^{\Delta t}K_{Q}-K_{Q}U^{\Delta t}\right\|=O\left(Q^{-1}\right)$, 2. (ii)

$\left\|U^{\Delta t}\tilde{K}_{Q}-\tilde{K}_{Q}U^{\Delta t}\right\|=O\left(Q^{-1}\right)$, 3. (iii)

$\left\|U^{\Delta t}P_{Q}-P_{Q}U^{\Delta t}\right\|=O\left(Q^{-1}\right)$.

Proof.

Let $\tilde{F}_{Q,\Delta t}(x,y):=\left\|F(x)-F(y)\right\|-\left\|F(\Phi^{Q\,\Delta t}x)-F(\Phi^{Q\,\Delta t}(y))\right\|$, and notice that by continuity of $F$ and compactness of $X$ this quantity is bounded on $X\times X$. Note that (i) $d_{Q}(\Phi^{\Delta t}(x),\Phi^{\Delta t}(y))=d_{Q}(x,y)+Q^{-1}\tilde{F}_{Q,\Delta t}(x,y)$; and (ii) $h(u+\Delta u)=h(u)+\Delta u\,h^{\prime}(u)+o(\Delta u)$, as $\Delta u\to 0$. Thus,

[TABLE]

where the estimate holds uniformly with respect to $x,y\in X$. Therefore, for every $f\in L^{2}(X,\mu)$ and $x\in X$ we have

[TABLE]

Note that we have used the fact that $\mu$ is an invariant measure in the second-to-last line. Since $k_{Q}$ is continuous, it follows from the Cauchy-Schwarz inequality that $\|K_{Q}f\|_{C^{0}}\leq\|K_{Q}\|_{C^{0}}\|f\|_{L^{2}}$. Substituting this result in the right-hand side and taking the supremum over $x\in X$ yields

[TABLE]

Claim (i) then follows from the fact that $\lVert\cdot\rVert_{L^{2}}\leq\lVert\cdot\rVert_{C^{0}}$. Claims (ii) and (iii) can be proved in a similar manner. ∎∎

6 Galerkin approximation of Koopman eigenvalue problems

In this section, we formulate a Galerkin method for the eigenvalue problem of the Koopman generator $V$ in the eigenbasis of $P$, under the implicit assumption that the latter operator is available to us from $P_{Q}$ after having taken a large number of delays $Q$. The task of finding the eigenvalues of $V$ has two challenges, namely, (i) $V$ is an unbounded operator defined on a proper subspace $D(V)\subset L^{2}(X,\mu)$ which is not known a priori; (ii) the spectrum of $V$ could be dense in $i\mathbb{R}$ (even for a pure point spectrum system such an ergodic rotation on $\mathbb{T}^{D}$ with $D\geq 2$; e.g., [33], Remark 8), in which case, solving for its eigenvalues is a numerically ill-posed problem. Following [30, 33], we will address these issues by employing a Galerkin scheme for the eigenvalue problem of $V$, with a small amount of judiciously constructed diffusion added for regularization. Throughout this section, we consider that Assumptions 1, 2, and 3 hold. Further, we assume the following.

Assumption 5**.**

The kernels $k_{Q}$, and thus $k_{\infty}$, are symmetric positive-definite. That is, (i) $k_{Q}(x,y)$, for every $x,y\in M$; (ii) for every $x_{0},x_{1},\ldots,x_{n}\in M$ and $c_{0},c_{1},\ldots,c_{n}\in\mathbb{C}$, $\sum_{i,j=0}^{n-1}c_{i}^{*}k_{Q}(x_{i},x_{j})c_{j}\geq 0$; and (iii) the analogous conditions hold for $k_{\infty}$.

Our approach has the following steps.

Step 1. Sobolev spaces.

We first construct subspaces of $L^{2}$ in which we search for eigenfunctions. These spaces will be shown to be dense in $H$, defined as the closed subspace of $\overline{\operatorname{ran}P}$ orthogonal to constant functions (that is, $H$ only consists of zero-mean functions). Note that $\{\phi_{j}\}_{j\in J}$, where $J$ is an index set for the nonzero eigenvalues $\lambda_{j}$ of $P$, strictly less than 1, is an orthonormal basis of $H$. For any $p\geq 0$, we define

[TABLE]

The spaces $H^{p}$ are analogous to the usual Sobolev spaces associated with self-adjoint, positive semidefinite, unbounded operators with compact resolvents and discrete spectra (here, $\{\eta_{j}\}_{j\in J}$). In particular, when $(X,g)$ is a smooth Riemannian manifold with a metric tensor $g$ satisfying $\operatorname{vol}_{g}=\mu$, and $(\eta_{j},\phi_{j})$ are the eigenvalues and orthonormal eigenfunctions of the corresponding Laplace-Beltrami operator, then $H^{p}$ becomes the canonical Sobolev space $H^{p}(X,g)$, restricted to be orthogonal to constant functions. $H^{p}$ from (20) is a Hilbert space with the inner product

[TABLE]

where $f=\sum_{j\in J}c_{j}\phi_{j}$ and $g=\sum_{j\in J}d_{j}\phi_{j}$. Moreover, $\{\phi_{j}^{(p)}\}_{j\in J}$ with $\phi_{j}^{(p)}=\phi_{j}/\lVert\phi_{j}\rVert_{p}$, $\lVert\phi_{j}\rVert_{H^{p}}^{2}=\sum_{q=0}^{p}\lambda_{j}^{q}$, forms an orthonormal basis of $H^{p}$.

Proposition 20**.**

For every $p>0$, the space $H^{p}$ is dense in $H$ and moreover, the inclusion map $H^{p}\to H$, and thus $H^{p}\to L^{2}(X,\mu)$, is compact.

Proof.

To see that $H^{p}$ is dense, note that $H^{p}$ includes all finite linear combinations of the $\phi_{j}$. Since the $\phi_{j}$ are an orthonormal basis of $H$, these finite linear combinations are dense in $H$. Next, the embedding of $H^{p}$ in $\mathcal{D}$ can be represented by a diagonal operator $G:H^{p}\to H$ such that $G_{jj}:=\langle\phi_{j},G\phi_{j}^{(p)}\rangle=\eta_{j}^{-p/2}$. This operator is compact iff $G_{jj}$ converges to [math] as $j\to\infty$. This is true by (22) below. The compactness of the inclusion $H^{p}\to L^{2}(X,\mu)$ follows immediately. ∎∎

Step 2. Regularized generator.

For every $\theta>0$, we define the unbounded operators $\Delta:D(\Delta)\to H$ and $L_{\theta}:D(L_{\theta})\to H$, where $D(\Delta)=D(L_{\theta})\subset D(V)$, and

[TABLE]

As we will see in Step 3 below, the role of the diffusion term $\theta\Delta$ is to penalize the eigenfunctions of $V$ with large eigenvalues of a Dirichlet energy functional. Theorem 21 below identifies a domain in which the operators in (21) are continuous, and establishes that the eigensolutions of $L_{\theta}$ converge to eigensolutions of $V$ as $\theta\to 0$.

Theorem 21**.**

Viewed as operators from $H^{2}$ to $H$, the generator $V$, as well as the operators $L_{\theta}$ and $\Delta$ from (21), are bounded. In particular, we can set $D(\Delta)=D(L_{\theta})=H^{2}$. Finally, for every eigenvalue $i\omega$ of $V$, whose corresponding eigenspace lies in $H^{2}$, there exists an eigenvalue $\eta$ of $\Delta$ such that the smooth curve $\theta\mapsto\gamma_{\theta}:=i\omega-\theta\eta$ consists of eigenvalues $\gamma_{\theta}$ of $L_{\theta}$, converging to $i\omega$ as $\theta\to 0^{+}$.

Proof.

First, by Corollary 18, we can consider that the basis $\{\phi_{j}\}_{j\in J}$ of $H$ consists of simultaneous eigenfunctions of $V$ and $P$ (and thus $\Delta$), without loss of generality. Then, to verify that $V$ is a bounded operator on $H^{2}$, first observe that $\omega_{j}=O(j)$, which follows from the fact that the eigenvalues of $V$ are integer linear combinations of finitely many rationally independent frequencies (by Assumption 3; see also Section 3). By Theorem 4, the kernel $p$ associated with $P$ is $L^{2}$ integrable, and thus by a result of Ferreira and Menegatto on integral operators ([25], Corollary 2.5), $\lambda_{j}=o(j^{-1})$. Combining these estimates, we obtain

[TABLE]

and therefore deduce that there exists a constant $C>0$ such that

[TABLE]

Hence, for $f=\sum_{j\in J}c_{j}\phi_{j}\in H^{2}$,

[TABLE]

proving that $V$ is a bounded operator on $H^{2}$. The same reasoning applies for $L_{\theta}$ and $\Delta$. Finally, to establish convergence of the eigenvalues of $L_{\theta}$ to those of $V\rvert_{H^{2}}$, let $i\omega_{j}$ be the eigenvalue of $V$ corresponding to $\phi_{j}$. Then, by definition of $L_{\theta}$ and the basis $\{\phi_{j}\}_{j\in J}$,

[TABLE]

and the claim follows immediately. This completes the proof of Theorem 21. ∎

Remark**.**

Theorem 21 establishes that $H^{2}$ is a domain on which $V$ is a bounded operator, but if $X$ had a smooth manifold structure, it is possible to show that the standard $H^{1}$ Sobolev space associated with a Riemannian metric on $X$ is also a suitable domain. In this work, $X$ has no smooth structure, and we can state Theorem 21 above only for $V\rvert_{H^{2}}$. In separate calculations, we have observed that an analog of the weak eigenvalue problem for $L_{\theta}$ formulated in $H^{1}\times H^{1}$ actually performs well numerically.

Step 3. Galerkin method.

By virtue of Theorem 21, the eigenvalues of $L_{\theta}$ can be considered to be approximations of the eigenvalues of $V$. We will take the Galerkin approach in finding the eigenvalues of $L_{\theta}$ by solving for $z\in H^{2}$ and $\gamma\in\mathbb{C}$ in the following variational (weak) eigenvalue problem:

Definition 22** (Regularized Koopman eigenvalue problem).**

Find $\gamma\in\mathbb{C}$ and $z\in H^{2}$ such that for all $f\in H$,

[TABLE]

where $A:H\times H^{2}\to\mathbb{C}$ is the sesquilinear form defined by

[TABLE]

In the above, the form $E:H\times H^{2}\to\mathbb{C}$ induces a Dirichlet energy functional $E(f)=E(f,f)$, $f\in H^{2}$, providing a measure of roughness of functions in $H^{2}$. In particular, if $X$ were a smooth Riemannian manifold, and the $(\eta_{j},\phi_{j})$ were set to Laplace-Beltrami eigenvalues and eigenfunctions, respectively, we would have $E(f)=\int_{X}\lVert\operatorname{grad}f\rVert^{2}\,d\mu$. While the lack of smoothness of $X$ in our setting precludes us from defining $E$ by means of a gradient operator, its definition in terms of the $\eta_{j}$ from (20) still provides a meaningful measure of roughness of functions. For instance, it follows from results in spectral graph theory that the variance of estimates $\eta_{j}^{(N)}$ of the $\eta_{j}$ computed from finite data sets (e.g., as described in Section 7 ahead) increases with $k$ [47, 10], which is consistent with the intuitive expectation that rough (highly oscillatory) functions require larger numbers of samples for accurate approximations.

Following [30, 33], we will order all solutions $(\gamma_{j},z_{j})$ of the problem in Definition 22 in order of increasing Dirichlet energy $E(z_{j})$. Since $A(f,f)=-\theta E(f,f)$ by skew-symmetry of $V$, we can compute the Dirichlet energy of eigenfunction $z_{j}$ directly from the corresponding eigenvalue, viz. $E(z_{j})=-\operatorname{Re}\gamma_{j}/\theta$. Similarly, we have $\omega_{j}=\operatorname{Im}\gamma_{j}$. By (23), there exist constants $C_{1},C_{2}>0$ such that

[TABLE]

To justify the well-posedness of the eigenvalue problem in Definition 22, we will state three important properties of $A$, namely,

[TABLE]

[TABLE]

[TABLE]

We now give brief proofs of these results. In the following, $v=\sum_{j\in J}d_{j}\phi_{j}$ and $u=\sum_{j\in J}c_{j}\phi_{j}$ will be arbitrary functions in $H^{2}$ and $H$, respectively. Moreover, as in the proof of Theorem 21, we will assume that the basis $\{\phi_{j}\}_{j\in J}$ consists of simultaneous eigenfunctions of $V$ and $\Delta$. First, note that,

[TABLE]

By the Cauchy-Schwartz inequality on $\ell^{2}$ and (24),

[TABLE]

proving (25). To prove (26), let $f=\sum_{j\in J}a_{j}\phi_{j}\in H$. Then, the left-hand side of that equation becomes $\sum_{j\in J}(i\omega_{j}/\eta_{j}-\theta)\eta_{j}a_{j}^{*}d_{j}$. Let $R_{j}:=i\omega_{j}/\eta_{j}-\theta$, where $|R_{j}|\geq C_{2}$ by (24), By the Cauchy-Schwarz inequality, under the constraint $\sum_{j\in J}|a_{j}|^{2}=1$, the sum $\left|\sum_{j\in J}a_{j}^{*}\eta_{j}d_{j}\right|$ attains the maximum value of $\sum_{j\in J}|\eta_{j}^{2}d_{j}|^{2}$. Therefore,

[TABLE]

This proves (26). The proof of (27) is similar to that of (26), with $f$ replaced by a trial function $g=\sum_{j\in J}b_{j}\phi_{j}\in H^{2}$ and the constraint $\|g\|_{H^{2}}^{2}=\sum_{j\in J}|b_{j}|^{2}=1$. A direct consequence of (26) and (27) is,

[TABLE]

Equations (25), (26), (28), and the compact embedding of $H^{2}$ in $H$ by Proposition 20 together guarantee that the eigenvalues of $A$ restricted to the finite-dimensional subspaces of $H\times H^{2}$ spanned by the leading $m$ eigenfunctions $\phi_{1},\ldots,\phi_{m}$ converge, as $m\to\infty$, to the weak eigenvalues of $L_{\theta}$. See [4], Section 8, for an exposition on this classic result. The resulting finite-dimensional Galerkin approximations of the weak eigenvalue problem for $L_{\theta}$ can be summarized as follows:

Definition 23** (Koopman eigenvalue problem, Galerkin approximation).**

Set $\tilde{H}_{m}=\operatorname{span}\{\phi_{1},\ldots,\phi_{m}\}$ and $\tilde{H}^{2}_{m}=\operatorname{span}\{\phi^{(2)}_{1},\ldots,\phi^{(2)}_{m}\}$, $m\geq 1$. Then, find $\gamma\in\mathbb{C}$ and $z\in\tilde{H}^{2}_{m}$ such that for all $f\in\tilde{H}_{m}$,

[TABLE]

where the sesquilinear form $A:H\times H^{2}\to\mathbb{C}$ is as in Definition 22.

This problem is equivalent to solving a matrix generalized eigenvalue problem

[TABLE]

where $\boldsymbol{A}$ and $\boldsymbol{B}$ are $m\times m$ matrices with elements

[TABLE]

respectively, and $\vec{c}=(c_{1},\ldots,c_{m})^{\top}$ is a column vector in $\mathbb{C}^{m}$ containing the expansion coefficients of the solution $z$ in the $\{\phi_{j}^{(2)}\}$ basis of $\tilde{H}^{2}_{m}$, viz. $z=\sum_{k=1}^{m}c_{j}\phi_{j}^{(2)}$. It is important to note that, unlike the proofs of Theorem 21 and (25)–(27), in (30) we do not require that the $\phi_{j}$ be simultaneous eigenfunctions of $V$ and $P$. This concludes the description of our Galerkin approximation of the eigenvalue problem for $L_{\theta}$ and therefore for $V$.

7 Data-driven approximation

In this section, we discuss the numeric procedures used to approximate the integral operators described in Sections 4, 5, and implement the Galerkin method of Section 6 using a finite, time-ordered dataset of observations $\left(F(x_{n})\right)_{n=0}^{N-1}$. In addition, we will prove Theorem 5. Throughout this section, we will assume that Assumptions 1–3 hold. In particular, by Assumption 4, we can assume without loss of generality that the underlying trajectory $\left(x_{n}\right)_{n=0}^{N-1}$ starts at a point $x_{0}$ in the compact set $\mathcal{U}$ (for, if $x_{0}$ were to lie in $\mathcal{V}\setminus\mathcal{U}$, the trajectory would enter $\mathcal{U}$ after finitely many steps, and its portion lying in $\mathcal{V}\setminus\mathcal{U}$ would not affect the asymptotic behavior of our schemes as $N\to\infty$). Besides this assumption, the trajectory $\left(x_{n}\right)_{n=0}^{N-1}$ is assumed to be unknown, and note that it need not lie on $X$.

For the purposes of the analysis that follows, it will be important to distinguish between operators that act on $L^{2}$ and $C^{0}$ spaces. Specifically, to every kernel $k:M\times M\to\mathbb{R}$ satisfying Assumption 2, we will assign a bounded operator $K^{\prime}:L^{2}(X,\mu)\to C^{0}(\mathcal{U})$, acting on $f\in L^{2}(X,\mu)$ via the same integral formula as in (3), but with the image $K^{\prime}f$ understood as an everywhere-defined, continuous function on $\mathcal{U}$. With this definition, the operator $K:L^{2}(X,\mu)\to L^{2}(X,\mu)$ acting on $L^{2}$ equivalence classes can be expressed as as $K^{\prime\prime}=\iota\circ K^{\prime}$, where $\iota:C^{0}(\mathcal{U})\to L^{2}(X,\mu)$ is the canonical $L^{2}$ inclusion map on $C^{0}(\mathcal{U})$, and we can also define an analog $K^{\prime\prime}:C^{0}(\mathcal{U})\to C^{0}(\mathcal{U})$ acting on continuous functions via $K^{\prime\prime}=K^{\prime}\circ\iota$. It can be verified using the Arzelà-Ascoli theorem that $K^{\prime\prime}$ is compact.

Data-driven Hilbert spaces.

Let $\mu_{N}:=N^{-1}\sum_{n=0}^{N-1}\delta_{x_{n}}$ be the sampling probability measure associated with the finite trajectory $(x_{n})_{n=0}^{N-1}$. The compact set $\mathcal{U}$ from Assumption 4 always contains the support of $\mu_{N}$. Moreover, since $x_{0}$ lies in the basin of the physical measure $\mu$, as $N\to\infty$, $\mu_{N}$ converges weakly to $\mu$, in the sense that

[TABLE]

Our data-driven analog of the space $L^{2}(X,\mu)$ will be $L^{2}(\mathcal{U},\mu_{N})$; the set of equivalence classes of complex-valued functions on $M$ which are square-summable and have common values at the sampled states $x_{n}$. Note that $L^{2}(\mathcal{U},\mu_{N})\cong\mathbb{C}^{N}$, and therefore every element $f\in L^{2}(\mathcal{U},\mu_{N})$ can be represented in the canonical basis of $\mathbb{C}^{N}$ as an $N$-vector $\vec{f}=(f(x_{0}),\ldots,f(x_{N-1}))$. In fact, $L^{2}(\mathcal{U},\mu_{N})$ is the image of $C^{0}(\mathcal{U})$ under the restriction map $\pi_{N}:C^{0}(\mathcal{U})\to L^{2}(\mathcal{U},\mu_{N})$, where $\pi_{N}f=\left(f(x_{0}),\ldots,f(x_{N-1})\right)$. Moreover, given any $f,g\in L^{2}(\mathcal{U},\mu_{N})$, we have $\langle f,g\rangle_{L^{2}(\mathcal{U},\mu_{N})}=\vec{f}\cdot\vec{g}/N$, where $\cdot$ denotes the canonical inner product on $\mathbb{C}^{N}$.

Kernel integral operators.

In the data-driven setting, given a continuous kernel $k:M\times M\to\mathbb{R}$, we define a kernel integral operator $K_{N}:L^{2}(\mathcal{U},\mu_{N})\to C^{0}(\mathcal{U})$ by (cf. (3))

[TABLE]

and we also set $K_{N}:L^{2}(\mathcal{U},\mu_{N})\to L^{2}(\mathcal{U},\mu_{N})$ and $K^{\prime\prime}_{N}:C^{0}(\mathcal{U})\to C^{0}(\mathcal{U})$ with $K_{N}=\pi_{N}\circ K^{\prime}_{N}$ and $K^{\prime\prime}_{N}=K^{\prime}_{N}\circ\pi_{N}$. Note that $K_{N}$ can be represented by an $N\times N$ matrix $\boldsymbol{K}$ with elements $K_{ij}=k(x_{i},x_{j})$. In this representation, the function $g=K_{N}f$, $f\in L^{2}(\mathcal{U},\mu_{N})$, is represented by $\vec{g}=\boldsymbol{K}\vec{f}$.

When $k=k_{Q}$ from (10), one can similarly define operators $K^{\prime}_{Q,N}:L^{2}(\mathcal{U},\mu_{N})\to C^{0}(\mathcal{U})$, $K_{Q,N}:L^{2}(\mathcal{U},\mu_{N})\to L^{2}(\mathcal{U},\mu_{N})$, and $K^{\prime\prime}_{Q,N}:C^{0}(\mathcal{U})\to C^{0}(\mathcal{U})$. This family of operators has the analogous properties to those stated for $K_{Q}$ in Lemma 10; namely, the functions $\rho_{Q,N}=K^{\prime\prime}_{Q,N}1_{\mathcal{U}}$ and $\sigma_{Q,N}=K^{\prime\prime}_{Q,N}(1/\rho_{Q,N})$ are both continuous, positive, and bounded away from zero on $\mathcal{U}$. Therefore, one can define a kernel $p_{Q,N}:M\times M\to\mathbb{R}$ by

[TABLE]

The kernel $p_{Q,N}$ has the Markov property, i.e., $\int_{\mathcal{U}}p_{Q,N}(x,\cdot)\,d\mu_{N}=1$ for every $x\in M$. Associated to $p_{Q,N}$ are the Markov operators $P^{\prime}_{Q,N}:L^{2}(\mathcal{U},\mu_{N})\to C^{0}(\mathcal{U})$ , $P_{Q,N}:L^{2}(\mathcal{U},\mu_{N})\to L^{2}(\mathcal{U},\mu_{N})$ and $P^{\prime\prime}_{Q,N}:C^{0}(\mathcal{U})\to C^{0}(\mathcal{U})$. Moreover, $P_{Q,N}$ is related to the self-adjoint operator $\hat{P}_{Q,N}:L^{2}(\mathcal{U,}\mu_{N})\to L^{2}(\mathcal{U},\mu_{N})$ with kernel $\hat{p}_{Q,N}:M\times M\to\mathbb{R}$,

[TABLE]

via a similarity transformation analogous to (19). From the kernel $\hat{p}_{Q,N}$ one can construct the operators $\hat{P}_{Q,N}$, $\hat{P}^{\prime}_{Q,N}$, and $\hat{P}^{\prime\prime}_{Q,N}$ as above.

Data-driven basis.

We will use the eigenvectors $\phi_{j,Q,N}$ of $\hat{P}_{Q,N}$ as an orthonormal basis of $L^{2}(\mathcal{U},\mu_{N})$, and employ the corresponding eigenvalues, $1=\lambda_{0,Q,N}>\lambda_{1,Q,N}\geq\cdots\geq\lambda_{N-1,Q,N}\geq 0$, to define data-driven analogs

[TABLE]

of the $\eta_{j}$ in (20), where $J_{N}=\{j:\lambda_{j,Q,N}>0\}$. The eigenvalue problem for $\hat{P}_{Q,N}$ is equivalent to a matrix eigenvalue problem for the $N\times N$ symmetric matrix $\hat{\boldsymbol{P}}=[\hat{p}_{Q,N}(x_{i},x_{j})]$ representing $\hat{P}_{Q,N}$. Details on the numerical solution of this problem can be found in [33, 34]. Note that for kernels $k_{Q}$ with exponential decay, such as the Gaussian kernels in (4), $\hat{\boldsymbol{P}}$ can be well approximated by a sparse matrix, allowing scalability of our techniques to large $N$.

To establish convergence of our schemes in the limit of large data, $N\to\infty$, we would like to establish a correspondence between the eigenvalues and eigenvectors of $\hat{P}_{Q,N}$ accessible from data and those of $\hat{P}_{Q}$, but because these operators act on the different spaces, a direct comparison of their eigenvectors is not possible. Therefore, as stated in Section 2, we will first establish a correspondence between the eigenvalues and eigenvectors of $\hat{P}_{Q,N}$ ($\hat{P}_{Q})$ and those of $\hat{P}^{\prime\prime}_{Q,N}$ ($\hat{P}^{\prime\prime}_{Q}$), and show that $\hat{P}^{\prime\prime}_{Q,N}$ spectrally converges to $\hat{P}^{\prime\prime}_{Q}$. The latter problem is meaningful since both $\hat{P}^{\prime\prime}_{Q,N}$ and $\hat{P}^{\prime\prime}_{Q}$ act on $C^{0}(\mathcal{U})$.

Lemma 24**.**

The following correspondence between the spectra of operators holds:

1. (i)

${\lambda}_{j,Q,N}$* is a nonzero eigenvalue of $\hat{P}_{Q,N}$ iff it is a nonzero eigenvalue of $\hat{P}^{\prime\prime}_{Q,N}$. Moreover, if ${\phi}_{j,Q,N}\in L^{2}(\mathcal{U},\mu_{N})$ is an eigenfunction of $\hat{P}_{Q,N}$ corresponding to $\lambda_{j,Q,N}$, then ${\varphi}_{j,Q,N}={\lambda}_{j,Q,N}^{-1}\hat{P}^{\prime}_{Q,N}{\phi}_{j,Q,N}\in C^{0}(\mathcal{U})$ is an eigenfunction of $\hat{P}^{\prime\prime}_{Q,N}$ corresponding to the same eigenvalue.* 2. (ii)

${\lambda}_{j,Q}$* is a nonzero eigenvalue of $\hat{P}_{Q}$ iff it is a nonzero eigenvalue of $\hat{P}^{\prime\prime}_{Q,N}$. Moreover, if ${\phi}_{j,Q}\in L^{2}(\mathcal{U},\mu)$ is an eigenfunction of $\hat{P}_{Q}$ corresponding to $\lambda_{j,Q}$, then ${\varphi}_{j,Q}={\lambda}_{j,Q}^{-1}\hat{P}^{\prime}_{Q}{\phi}_{j,Q}\in C^{0}(\mathcal{U})$ is an eigenfunction of $\hat{P}^{\prime\prime}_{Q}$ corresponding to the same eigenvalue.*

Lemma 24 is a direct consequence of the definitions of $\hat{P}_{Q,N}$ and $P_{Q,N}^{\prime\prime}$. Next, we establish spectral convergence of $\hat{P}^{\prime\prime}_{Q,N}$ to $\hat{P}^{\prime\prime}_{Q}$. For that, we will need the following notion of convergence of operators.

Compact convergence.

A sequence of operators $A_{n}$ on a Banach space $B$ is said to be compactly convergent to an operator $A$ if $A_{n}\to A$ pointwise, and for every bounded sequence of vectors $(f_{n})_{n\in\mathbb{N}}$, $f_{n}\in B$, the sequence $((A-A_{n})f_{n})_{n\in\mathbb{N}}$ has compact closure. The following proposition states that the data-driven operators $\hat{P}_{Q,N}$ converge compactly, and as result in spectrum; for a proof, see [47], Proposition 11, and [64], Theorem 2.4.1.

Proposition 25**.**

Let Assumptions 1–5 hold. Given a trajectory $(x_{n})_{n\in\mathbb{N}}$ starting in ${\mathcal{B}_{\mu}}$, the corresponding sequence of operators $\hat{P}^{\prime\prime}_{Q,N}$ constructed from the observations $F(x_{0}),\ldots,F(x_{N-1})$ converges compactly as $N\to\infty$ to $\hat{P}^{\prime\prime}_{Q}$. As a result, the sequence $\hat{P}_{Q,N}$ converges spectrally, in the sense of Corollary 2, to $\hat{P}_{Q}$. In particular, since the nonzero spectrum of a compact operator only consists of isolated eigenvalues, the convergence holds for all nonzero eigenvalues of $\hat{P}^{\prime\prime}_{Q,N}$ and the corresponding eigenspaces.

The spectral convergence results above follow from Proposition 13 in [47]. We will now prove Theorem 5. Note that there is some similarity between our methods and papers on spectral convergence of kernel algorithms, e.g., [6, 47], but our assumptions distinguishes Theorem 5 from previously studied cases. In particular, we do not assume an i.i.d. sequence of observed quantities, or that the sampled sequence $(x_{n})_{n=0}^{N-1}$ lies on the support $X$ of the invariant measure (as assumed in [6, 47]). Finally, $X$ need not have a manifold structure (as assumed in [6] and other manifold learning algorithms).

Proof of Theorem 5.

The claims of the theorem follow from analogous results to Lemma 24 and Proposition 25 for the operators $P_{Q,N}$, $P^{\prime}_{Q,N}$, $P^{\prime\prime}_{Q,N}$ and $P_{Q}$, $P^{\prime}_{Q}$, $P^{\prime\prime}_{Q}$. ∎

Together, Lemma 24 and Proposition 25 imply that every eigenpair $(\lambda_{j,Q},\phi_{j,Q})$ of $\hat{P}_{Q}$ can be consistently approximated by a sequence of eigenpairs $({\lambda}_{j,Q,N},{\phi}_{j,Q,N})$ of $\hat{P}_{Q,N}$. Moreover, by Corollary 2, as $Q\to\infty$, $(\lambda_{j,Q},\phi_{j,Q})$ approximates in turn the eigenpair $(\lambda_{j},\phi_{j})$ of $P$; that is,

[TABLE]

where the second limit is taken with respect to the $L^{2}(X,\mu)$ norm. Since, as can be seen in (30), the Galerkin scheme in Section 6 can be entirely formulated using the $\lambda_{j}$ and the matrix elements $\langle\phi_{i},V\phi_{j}^{(2)}\rangle$ of the generator, (34) indicates in turn that we can construct a consistent data-driven Galerkin scheme if we can consistently compute approximate generator matrix elements using the data-driven eigenfunctions ${\phi}_{j,Q,N}$. To that end, we will employ finite-difference approximations, as described below.

Finite-difference approximation.

The action $Vf$ of the generator on a function $f\in D(V)$ is defined via the limit in (7). This suggests that for data sampled discretely at sampling interval $\Delta t$, we can approximate $Vf$ by a finite-difference approximation [32, 30, 33]. For example, the following are first- and second-order, approximation schemes for $V$, respectively:

[TABLE]

In the finite-sample case, we approximate $V_{\Delta t}$ by a corresponding $r$-th order finite-difference operator $V_{\Delta t,N}:L^{2}(\mathcal{U},\mu_{N})\to L^{2}(\mathcal{U},\mu_{N})$. For example, in the case of the first-order scheme in (35), $V_{\Delta t,N}$ becomes

[TABLE]

and $V_{\Delta t,N}f(x_{N-1})=0$. To ensure that the approximations $V_{\Delta t,N}f$ converge to the true function $Vf$ for a class of functions of sufficient regularity, the following smoothness conditions are sufficient:

Assumption 6**.**

$\mathcal{U}$* is a $C^{1+\alpha}$ compact manifold for some $\alpha>0$, and $\Phi^{t}\rvert_{\mathcal{U}}$ is generated by a $C^{\alpha}$ vector field $\vec{V}$. Moreover, $F\rvert_{\mathcal{U}}\in C^{1+\alpha}(\mathcal{U};\mathbb{R}^{d})$, and the kernel shape function $h:\mathbb{R}\to\mathbb{R}$ is $C^{1+\alpha}$. $V_{\Delta t}$ and $V_{\Delta t,N}$ are first-order finite difference schemes, as in (35) and (36), respectively.*

Under Assumption 6, the flow $\Phi^{t}$ is generated by a $C^{\alpha}$ vector field $\vec{V}:C^{1}(\mathcal{U})\to C^{0}(\mathcal{U})$, and the generator $V$ of the Koopman group is an extension of the latter differential operator. Moreover, we can approximate $\vec{V}$ by finite-difference schemes $\vec{V}_{\Delta t}:C^{0}(\mathcal{U})\to C^{0}(\mathcal{U})$, defined analogously to (35) with $U^{\Delta t}$ replaced by $\Phi^{t}$. We then have:

Proposition 26**.**

Let Assumptions 1, 2, and 6 hold. Then for every $i,j\in\mathbb{N}$:

1. (i)

The eigenfunctions $\varphi_{j,Q,N}$ and $\varphi_{j,Q}$ from Lemma 24 lie in $C^{1+\alpha}(\mathcal{U})$. Moreover, as $\Delta t\to 0$,

[TABLE]

where the estimate holds uniformly on $\mathcal{U}$. 2. (ii)

$\lim_{\Delta t\to 0}\lim_{N\to\infty}\langle\hat{\phi}_{i},V_{\Delta t,N}\hat{\phi}_{j}\rangle_{L^{2}(\mathcal{U},\mu_{N})}=\langle\phi_{i,Q},V\phi_{j,Q}\rangle$.

Proof.

To prove Claim (i), note that under Assumption 6, for a finite number of delays $Q$, by (4), $\hat{p}_{Q}$ is a $C^{1+\alpha}$-smooth kernel. Hence, according to [25], the ranges of the integral operators $\hat{P}^{\prime}_{Q}$ and $\hat{P}^{\prime}_{Q,N}$, and thus $\varphi_{j,Q}$ and $\varphi_{j,Q,N}$ lie in $C^{1+\alpha}(\mathcal{U})$. Since the vector field $\vec{V}$ is $C^{\alpha}$, the trajectories are $C^{1+\alpha}$-smooth, and therefore, $\vec{V}\varphi_{j,Q}$, which is the time derivative along the orbit, has a first-order Taylor expansion. The first-order finite-difference scheme gives the $\|\varphi_{j,Q}\|_{C^{1+\alpha}(\mathcal{U})}O(\Delta t)^{\alpha}$ error. Claim (ii) is a consequence of Claim (i), in conjunction with the weak convergence of measures in (31) and Lemma 24. ∎∎

Remark**.**

In many cases, such as flows induced on inertial manifolds in dissipative PDEs [18], the $C^{1+\alpha}$ regularity in Assumption 6 cannot be strengthened. Proposition 26 provides the basis for numerically approximating $V$ for these cases. If $M$, $\mathcal{U}$, $\vec{V}$, $F$ and $h$ have a higher degree of smoothness, say $C^{r}$ for some $r\geq 2$, then taking $V_{\Delta t}$ to be an $r$-th order finite-difference scheme would lead to an improved, $O(\Delta t)^{r}$, convergence.

Data-driven Galerkin method.

Using the ${\eta}_{j,Q,N}$ from (33), we define the data-driven normalized basis vectors $\hat{\phi}_{j}^{(p)}=\phi_{j,Q,N}/\hat{\eta}^{p/2}_{j,Q,N}$, $j\in J_{N}$ (cf. the $\phi_{j}^{(p)}$ from Step 1 in Section 6), and the associated Galerkin approximation spaces $H^{p}_{N,Q,m}=\operatorname{span}\{\hat{\phi}^{(p)}_{j}\}_{j=1}^{m}\subseteq L^{2}(\mathcal{U},\mu_{N})$, $m\leq J_{N}$, where we abbreviate $H^{p}_{Q,N,J_{N}}=:H^{p}_{Q,N}$ and $H^{0}_{Q,N}=:H_{Q,N}$. We also define the positive semidefinite, self-adjoint operator $\Delta_{Q,N}:H_{Q,N}\to H_{Q,N}$, where

[TABLE]

This operator is a data-driven analog of $\Delta$ in (21). With these definitions and the finite-difference approximation of $V$ described above, we pose the following data-driven analog of the Galerkin approximation in Definition 23:

Definition 27** (Koopman eigenvalue problem, data-driven form).**

Find $\gamma\in\mathbb{C}$ and $z\in H^{2}_{Q,N,m}$ such that for all $f\in H_{Q,N,m}$,

[TABLE]

where $A_{\Delta t,Q,N}:H_{Q,N}\times H_{Q,N}\to\mathbb{C}$ is the sesquilinear form defined as

[TABLE]

Numerically, this is equivalent to solving a matrix generalized eigenvalue problem analogous to that in (29), viz.

[TABLE]

where $\boldsymbol{A}$ and $\boldsymbol{B}$ are $m\times m$ matrices with elements

[TABLE]

respectively, and $\vec{c}=(c_{1},\ldots,c_{m})^{\top}$ is a column vector in $\mathbb{C}^{m}$ containing the expansion coefficients of the solution $z=\sum_{j=1}^{m}c_{j}\hat{\phi}_{j,Q,N}^{(2)}$ in the $\{\hat{\phi}_{j}^{(2)}\}$ basis of $H^{2}_{Q,N,m}$. Analogously to the continuous case, we define a data-driven Dirichlet energy functional $E_{Q,N}$ on $H^{2}_{Q,N}$, given by $E_{Q,N}(f)=\langle f,\Delta_{Q,N}f\rangle_{L^{2}(\mathcal{U},\mu_{N})}$, and use that functional to order the computed eigenfunctions in order of increasing Dirichlet energy. Note that, unless an antisymmetrization is explicitly performed, in the data driven setting, $V_{ij}$ will generally not be equal to $-V_{ji}$, and thus $\operatorname{Re}\gamma$ will not be equal to $-\theta E_{Q,N}(z)$ (cf. Section 6). Nevertheless, in practice we observe that $\operatorname{Re}\gamma\approx-\theta E_{Q,N}(z)$, at least for the leading eigenfunctions.

For any fixed $m$, and up to similarity transformations, the matrices $\boldsymbol{A}$ and $\boldsymbol{B}$ converge in the limits $Q\to\infty$, after $\Delta t\to 0$, after $N\to\infty$ (in that order) to the corresponding matrices in the variational eigenvalue problem in (29). We therefore conclude that the data-driven Galerkin method in Definition 27 is consistent (as $\Delta t\to 0$ and $Q,N\to\infty$) with the Galerkin method in Definition 23, which is in turn consistent (as $m\to\infty$) with the weak eigenvalue problem for the regularized generator $L_{\theta}$ in Definition 22.

8 Results and discussion

In this section, we apply the methods described in Sections 4–7 to two ergodic dynamical systems with mixed spectrum, constructed as products of either a mixing flow on the 3-torus, or the L63 system, with circle rotations. Our objectives are to demonstrate that (i) the results of Theorem 1 and Corollaries 2 and 3 hold, that is, the eigenfunctions $\hat{\phi}_{j}$ of $P_{Q,N}$ from (17) are eigenfunctions of $U^{t}$; (ii) the eigenvalues obtained using the Galerkin scheme in Definition 27 are consistent with those expected theoretically.

8.1 Two systems with mixed spectrum

The first system studied below is based on a strongly mixing flow on the 3-torus introduced by Fayad [24]. The flow, denoted by $\Phi_{\mathbb{T}^{3}}^{t}$, is given by the solution of the ordinary differential equation (ODE) $d(x,y,z)/dt=\vec{V}(x,y,z)$, where $(x,y,z)\in\mathbb{T}^{3}$, and $\vec{V}$ is the smooth vector field

[TABLE]

parameterized by the constant frequency vector $\vec{\nu}$. Hereafter, we set $\vec{\nu}=(\sqrt{2},\sqrt{10},1)^{\top}$. Note that the orbits under $\Phi_{\mathbb{T}^{3}}^{t}$ are the same as that of the ergodic, non-mixing linear flow with constant vector $\vec{\nu}$. $\Phi_{\mathbb{T}^{3}}^{t}$ has a unique Borel, invariant, ergodic probability measure $\mu$ with density $\varphi/\int_{M}\varphi\,d\operatorname{Leb}$ relative to Lebesgue measure. Such flows are also called reparameterized flows as $\vec{\nu}$ is scaled by the function $\varphi$ at each point $(x,y,z)\in\mathbb{T}^{3}$.

This system is strongly mixing with respect to its invariant measure, i.e., its generator has purely continuous spectrum [24]. To construct an associated mixed-spectrum system, we take the product $\Phi_{\mathbb{T}^{3}}^{t}\times\Phi_{\omega}^{t}$ with a periodic flow $\Phi_{\omega}^{t}$ on $S^{1}$, defined as

[TABLE]

Thus, the state space of the product system is $M=\mathbb{T}^{3}\times S^{1}=\mathbb{T}^{4}$. Note that in this example the attracting set $X$ is smooth and coincides with the state space, $M=X$; in particular, all states sampled experimentally lie exactly on $X$. Moreover, the Koopman generator $V:D(V)\to L^{2}(X,\mu)$ is a skew-adjoint extension of the differential operator $\vec{V}\oplus\vec{\omega}:C^{\infty}(X)\to L^{2}(X,\mu)$, where $\vec{\omega}:C^{\infty}(S^{1})\to C^{\infty}(S^{1})$ is the differential operator $f\mapsto\vec{\omega}(f):=\omega f^{\prime}$. Since $\Phi_{\omega}^{t}$ has a pure point spectrum consisting of integer multiples of $i\omega$ and $\Phi_{\mathbb{T}^{3}}^{t}$ has no eigenvalues, the discrete spectrum of the product system is $\{ik\omega,k\in\mathbb{Z}\}$.

The second system that we study is based on the L63 system [45]. This system is known to have a chaotic attractor $X_{\text{Lor}}\subset\mathbb{R}^{3}$ with fractal dimension $2.0627160$ [49], supporting a physical invariant measure [63] with a compact absorbing ball [43] and mixing dynamics [48]. The latter, implies that the generator $V$ of the system has only constant eigenfunctions, corresponding to eigenvalue [math]. The flow, denoted by $\Phi^{t}_{\text{Lor}}$, is generated by a smooth vector field $\vec{V}\in C^{\infty}(\mathbb{R}^{3};\mathbb{R}^{3})$, whose components at $(x,y,z)\in\mathbb{R}^{3}$ are

[TABLE]

Throughout, we use the standard parameter values $\beta=8/3$, $\rho=28$, $\sigma=10$. As in the torus case, we form the product $\Phi^{t}_{\text{Lor}}\times\Phi_{\omega}^{t}$ with the rotation $\Phi_{\omega}^{t}$ in (38), leading to a mixed spectrum system with the same discrete spectrum $\{ik\omega,k\in\mathbb{Z}\}$. Note that unlike the torus-based system, the attracting set $X=X_{\text{Lor}}\times S^{1}$ is a strict subset of the state space $M=\mathbb{R}^{3}\times S^{1}$.

For each product system, we define a continuous map $F:M\to\mathbb{R}^{3}$ coupling the degrees of freedom of the continuous-spectrum subsystem with the rotation. In the case of the torus-based system, we define $F(x,y,z,\alpha)=(F_{1},F_{2},F_{3})$, $(x,y,z)\in\mathbb{T}^{3}$, $\alpha\in S^{1}$, via additive coupling, viz.

[TABLE]

In the case of the L63-based system, the coupling is nonlinear with $F(x,y,z,\alpha)=(F_{1},F_{2},F_{3})$, $(x,y,z)\in\mathbb{R}^{3}$, $\alpha\in S^{1}$, and

[TABLE]

8.2 Experimental results

We generated numerical trajectories $x_{0},x_{1},\ldots,x_{N-1}$ of the torus- and L63-based systems described above starting in each case from an arbitrary initial condition $y\in M$. In the torus experiments, the system is always on the attractor, so the starting state $x_{0}$ in the training data was set to $y$. In the L63 experiments, we let the system relax towards the attractor, and set $x_{0}$ to a state sampled after a long spinup time (4000 time units); that is, we formally assume that $y$ (and therefore $x_{0}$) lie in the basin ${\mathcal{B}_{\mu}}$ of the physical measure associated with $X$. In both cases, the number of samples was $N=\text{50,000}$, the integration time-step was $0.01$, and the number of delays was $Q=2000$. Gaussian kernels $k_{Q}$ from (5) were used throughout. We employed the ode45 solver of Matlab to compute the trajectories and generated time series $F(x_{0}),F(x_{1}),\ldots,F(x_{N-1})$ by applying the observation maps in (40) and (41) to the respective states $x_{n}$. Portions of the observable time series from each system are displayed in Fig. 2. Note that the $x_{n}$ were not presented to our kernel algorithm.

We computed data-driven eigenpairs $(\lambda_{j,Q,N},\phi_{j,Q,N})$ by solving the eigenvalue problem for the operator $\hat{P}_{Q,N}$ from (32) , using Matlab’s eigs iterative solver. Henceforth, for ease of notation, we abbreviate $\lambda_{j,Q,N}$ and $\phi_{j,Q,N}$ by $\hat{\lambda}_{j}$ and $\hat{\phi}_{j}$, respectively. The bandwidth parameter $\epsilon$ of the Gaussian kernels was selected using the tuning procedure described in [9, 17, 8], which yielded $\epsilon\approx 3.6$ and $\approx 2.053$ for the torus and L63 systems, respectively. Representative eigenfunctions $\hat{\phi}_{j}$, plotted as time series $n\to\hat{\phi}_{j}(x_{n})$, and the corresponding eigenvalues are displayed in Figs. 1 and 3, respectively. We now describe these results in more detail.

According to Theorem 1 and Proposition 17, at large numbers of delays (here, $Q=2000$), the eigenfunctions $\hat{\phi}_{j}$ of $\hat{P}_{Q,N}$ should form doubly degenerate pairs, and each pair should exhibit a single frequency associated with an eigenvalue of $V$. More precisely, $\hat{\phi}_{j}\pm i\hat{\phi}_{j+1}$ with $j\in\{1,3,\ldots\}$ should approximate an eigenfunction of $V$. Both systems studied here have exactly one rationally independent eigenvalue $i\omega=i$, so the eigenfunctions of $P_{Q,N}$ are expected to evolve at frequencies $j\omega$, $j\in\mathbb{N}$. This is evidently the case in the time series plots in Fig. 1. Also, each of the $\hat{\phi}_{j}$ has multiplicity $2$ (note that only one eigenfunction from each eigenspace is shown in Fig. 1). The left-hand panels of Fig. 1 show a matrix representation of the generator $V$ (approximated via the finite-difference scheme in (36)) in the 51-dimensional data-driven subspace spanned by $\hat{\phi}_{0},\ldots,\hat{\phi}_{50}$. The matrix is, to a good approximation, skew-symmetric, consistent with the fact that $V$ is a skew-symmetric operator, and exhibits prominent $2\times 2$ diagonal blocks associated with the eigenspaces of $V$ approximated by $(\hat{\phi}_{1},\hat{\phi}_{2}),(\hat{\phi}_{3},\hat{\phi}_{4}),\ldots$, in agreement with Corollary 18.

Figure 4 shows the approximated eigenvalues $\gamma_{j}$ of the regularized generator $L_{\theta}$ obtained from this basis using the Galerkin scheme in Definition 27 with the diffusion regularization and spectral order parameters $\theta=10^{-4}$ and $m=50$, respectively. Each plot in Fig. 4 shows the first 20 eigenvalues corresponding to eigenfunctions of increasing Dirichlet energy $E(z_{j})$ of the corresponding eigenfunction $z_{j}$ (recall that $\operatorname{Re}\gamma_{j}\approx-\theta E_{Q,N}(z_{j})$). According to Section 6, the imaginary parts of the $\gamma_{j}$ should approximate the Koopman eigenfrequencies $j(k)\omega$, where $j$ is an integer-valued function giving the frequency of the Koopman eigenfunction with the $k$-th smallest Dirichlet energy. In Fig. 4, the $\operatorname{Im}\gamma_{j}$ are indeed equal to integer multiples of $\omega=1$ to a good approximation for the first $m$ eigenvalues (ordered in order of increasing Dirichlet energy). For indices $k$ close to $m$, the accuracy of the eigenvalues begins to deteriorate. This is due to the facts that (i) even with a “perfect” basis $\{\hat{\phi}_{j}\}$, eigenfunctions of higher Dirichlet energy (and stronger oscillatory behavior) require increasingly higher-order Galerkin approximation spaces; (ii) at finite sample numbers $N$, the quality of the data-driven elements $\hat{\phi}_{j}$ degrades at large $j$.

8.3 Discussion

The examples presented in Sections 8.1 and 8.2 are Cartesian products of weak mixing and quasiperiodic flows, with their phase variables combined through some observation map. We begin with some observations about our kernel method applied to Cartesian products.

Cartesian products.

Let $(X,\Phi_{X}^{t},\mu_{X})$ and $(Y,\Phi_{Y}^{t},\mu_{Y})$ be two ergodic flows on compact metric spaces with purely continuous and pure-point spectra, respectively. We are interested in the measure-preserving mixed spectrum dynamical system $(X\times Y,\Phi_{X}^{t}\times\Phi_{Y}^{t},\mu_{X}\times\mu_{Y})$. It is well known that the space $L^{2}(X\times Y,\mu_{X}\times\mu_{Y})$ is densely spanned by products of the form {$f\otimes g$ : $f\in L^{2}(X,\mu_{X})$, $g\in L^{2}(Y,\mu_{Y})$}. Recall that the observation map $F$ is the basis of our construction of all our data-driven operators. Corollary 28 below is a direct consequence of Proposition 12, and gives an “observability” condition that must be fulfilled by the observation map $F$ in order for the methods presented here to yield non-trivial results.

Corollary 28**.**

Let $(X,\Phi_{X}^{t},\mu_{X})$ and $(Y,\Phi_{Y}^{t},\mu_{Y})$ be as described above, and $F\in L^{2}(X\times Y,\mu_{X}\times\mu_{Y})$ be the sum $F=\sum_{n=1}^{\infty}f_{n}\otimes g_{n}$. Then, $F_{\mathcal{D}}=\sum_{n=1}^{\infty}E(f_{n})g_{n}$, where $E(f_{n})=\int_{X}f_{n}\,d\mu_{X}$. Hence, a necessary and sufficient condition that $P$ is not trivial is that $E(f_{n})\neq 0$ for at least one $n\in\mathbb{N}$.

Kernels with a small number of delays.

An implicit assumption in the approximation of the operator $P$ in (17) by the operator $P_{Q}$ in (17) with finitely many delays $Q$, is that $Q$ is large-enough for the asymptotic analysis of Lemma 19 to hold. When $Q$ is small, $d_{Q}$ is closer to a proper metric and therefore, the entries $K_{ij}=\exp(-d_{Q}(x_{i},x_{j})^{2}/\epsilon)$ of the kernel matrix $\boldsymbol{K}$ decay rapidly away from the diagonal $i=j$. Then $K_{ij}$ is close to a diagonal matrix, and $P_{i,j}$ is close to the identity matrix. On the other hand, for $Q$ large, $d_{Q}$ becomes a pseudo-metric and $P_{i,j}$ is not necessarily close to a diagonal matrix. Figure 3 shows how the Koopman eigenvalues computed for the two examples from (40) and (41) cluster near $1$ for $Q=1$ and decay more rapidly for $Q=2000$.

Systems with purely continuous spectra.

An important assumption of our kernel-based method is that the dynamics has Koopman eigenvalues, i.e., $\mathcal{D}$ contains non-constant functions. This underlies the ability of our regularized operator $L_{\theta}$ in (21) to be a suitable substitute of $V$ (Theorem 21). In fact, by Theorem 3, in the limit of infinitely many delays $Q\to\infty$, if $\mathcal{D}$ only contains constant functions, then the kernels $k_{Q}$, $p_{Q}$ converge to [math] (in the $L^{2}$ sense). However, when using finitely many delays, $k_{Q}\neq 0$, and correspondingly $p$ obtained by normalization of $k_{\infty}$ is not close to [math]. It is not currently understood how the operator $P_{Q}$ should behave for purely continuous spectrum systems (i.e., $\mathcal{D}=\operatorname{span}\{1_{X}\}$) and $Q<\infty$.

One of the consequences of Theorem 3(ii) is that in the limit of $Q\to\infty$, the continuous spectrum subspace $\mathcal{D}^{\perp}$ is annihilated by the integral operator $P_{\infty}$, thus rendering this operator ineffective for studying or reconstructing the mixing component of the dynamics. In particular, for weak-mixing systems, $P_{\infty}$ should have all but one of its eigenvalues to equal to zero. Numerical results shown in Fig. 5 indicate that the finite-rank, data-driven operator $P_{Q,N}$ for the L63 system still has nonzero eigenvalues strictly less than 1, but these eigenvalues are clustered around a small value ($\hat{\lambda}_{j}\approx 0.1$). This behavior is in agreement with Theorem 3, according to which all the eigenvalues of $P_{Q}$ other than $1$ should converge to zero as $Q\to\infty$. Note that the matrix representation of $V$ (also shown in Fig. 5) is still skew-symmetric to a good approximation, since $V$ is a skew-symmetric operator. Intriguingly, the matrix has a $2\times 2$ block-diagonal form, despite $V$ having no eigenfunctions. This form of the generator matrix has some aspects in common with the recent results of Brunton et al. [14], who obtained a bi-diagonal matrix representation of the L63 generator in a data-driven basis from Hankel matrix analysis. In Fig. 5, the lack of Koopman eigenfunctions is evident from the time-series plots of the numerical eigenfunctions $\hat{\phi}_{j}$, which are clearly non-periodic. Moreover, a phase space plot of $\hat{\phi}_{2}$ illustrates that it is a highly rough function on the Lorenz attractor.

In light of the above, the results established in this work have implications for delay-embedding techniques, as they point to a tradeoff between reconstruction of the system’s state space topology in delay embedding space (favored by large numbers of delays) and the ability of operators for data analysis, such as $P_{Q}$, to adequately represent the mixing component of the dynamics. Nevertheless, the ability to consistently approximate the quasiperiodic dynamics through Koopman eigenfunctions is still useful, as it allows identification and efficient modeling (e.g., via (9) of observables with high predictability). At the very least, the “negative” result described above provides a reference point that may aid the design of delay-embedding methodologies aiming to reconstruct the full structure of the dynamics. One of the goals of our future work is to investigate the behavior of the techniques presented here away from the asymptotic limit $Q\to\infty$ in the presence of a continuous spectrum.

Acknowledgements

Dimitrios Giannakis received support from ONR YIP grant N00014-16-1-2649, NSF grant DMS-1521775, and DARPA grant HR0011-16-C-0116. Suddhasattwa Das is supported as a postdoctoral research fellow from the first grant. The authors are grateful to L S Young for her suggestions.