Extraction and Prediction of Coherent Patterns in Incompressible Flows through Space-Time Koopman Analysis
Dimitrios Giannakis, Suddhasattwa Das
TL;DR
This paper introduces a data-driven method using Koopman operator theory and diffusion maps to detect, analyze, and predict coherent patterns in complex incompressible fluid flows, with applications to vortex and Lorenz systems.
Contribution
It develops a novel space-time Koopman analysis framework that captures global flow patterns and enables model-free predictions of flow evolution.
Findings
Successfully identified coherent structures in vortex flows.
Predicted flow evolution in Lorenz 96 systems.
Demonstrated effectiveness of the method on complex fluid flows.
Abstract
We develop methods for detecting and predicting the evolution of coherent spatiotemporal patterns in incompressible time-dependent fluid flows driven by ergodic dynamical systems. Our approach is based on representations of the generators of the Koopman and Perron-Frobenius groups of operators governing the evolution of observables and probability measures on Lagrangian tracers, respectively, in a smooth orthonormal basis learned from velocity field snapshots through the diffusion maps algorithm. These operators are defined on the product space between the state space of the fluid flow and the spatial domain in which the flow takes place, and as a result their eigenfunctions correspond to global space-time coherent patterns under a skew-product dynamical system. Moreover, using this data-driven representation of the generators in conjunction with Leja interpolation for matrix…
| , | , | ||||
|---|---|---|---|---|---|
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Extraction and Prediction of Coherent Patterns in Incompressible Flows through Space-Time Koopman Analysis
Dimitrios Giannakis
Suddhasattwa Das
Courant Institute of Mathematical Sciences, New York University, New York, NY 10012, USA
Abstract
We develop methods for detecting and predicting the evolution of coherent spatiotemporal patterns in incompressible time-dependent fluid flows driven by ergodic dynamical systems. Our approach is based on representations of the generators of the Koopman and Perron-Frobenius groups of operators governing the evolution of observables and probability measures on Lagrangian tracers, respectively, in a smooth orthonormal basis learned from velocity field snapshots through the diffusion maps algorithm. These operators are defined on the product space between the state space of the fluid flow and the spatial domain in which the flow takes place, and as a result their eigenfunctions correspond to global space-time coherent patterns under a skew-product dynamical system. Moreover, using this data-driven representation of the generators in conjunction with Leja interpolation for matrix exponentiation, we construct model-free prediction schemes for the evolution of observables and probability densities defined on the tracers. We present applications to periodic Gaussian vortex flows and aperiodic flows generated by Lorenz 96 systems.
keywords:
Koopman operators, Perron-Frobenius operators, Lagrangian coherent structures, kernel methods, diffusion maps, nonparametric prediction
††journal: Physica D
1 Introduction
The formation of coherent structures is a ubiquitous feature of time-dependent flows in both natural and engineered systems, including coherent jets and vortices in planetary atmospheres, oceanic currents and eddies, and chemical mixing, among many other examples [1, 2, 3]. Objectively identifying and predicting these patterns has received considerable attention in the mathematical, physical, and engineering disciplines, resulting in a diverse range of techniques to achieve these goals. While virtually all such techniques utilize dynamical systems theory, they generally emphasize fairly distinct aspects of that theory, namely the geometric/state-space perspective [4, 5, 6, 7], or the operator-theoretic/probabilistic perspective [8, 9, 10, 11, 12] (though this dichotomy is not rigid as there are methods that employ aspects of both approaches; e.g., [13, 14, 15, 16, 17]). Among the many references in the literature, here we mention explicitly the paper of Liu and Haller [18] on coherent patterns of diffusive tracers in time-dependent flows. They show that the associated advection-diffusion equation for the tracer concentration field admits a finite-dimensional invariant manifold embedded in its solution space and spanned by a finite set of time-dependent modes, previously identified by Pierrehumbert [19] as strange eigenmodes. We also mention the recent work of Froyland and Koltai [20] who develop a method for recovering coherent patterns in time-periodic flows through a transfer operator defined on an extended state space where the dynamics are autonomous. Both [18, 20] have common aspects with the techniques presented below.
In this work, we study the problem of coherent pattern identification and prediction in time-dependent, incompressible flows having a skew-product structure. In particular, we consider that the time-dependence of the velocity field is itself the outcome of a dynamical system with ergodic properties—this allows us to take advantage of geometrical and operator-theoretic properties of that dynamical system that would otherwise not be present in flows with an arbitrary time dependence. In particular, we consider the product space between the state space for the velocity field dynamics (hereafter, ) and the physical domain (hereafter, ) where tracer motion takes place. Intuitively, that product space, , can be through of as a “space-time manifold” with the component playing the role of time (as it governs the non-autonomous aspect of tracer dynamics), and the role of space. On , the dynamics is autonomous and measure-preserving, though not necessarily ergodic, and has a natural skew-product structure owing to the fact that the state of the velocity field influences the dynamics of the tracers but not vice-versa. Moreover, associated with these dynamics are Koopman and Perron-Frobenius operators [21] governing the evolution of observables and probability measures defined on the tracers.
Our approach for identifying and predicting coherent spatiotemporal patterns is based on data-driven approximations of these operators and their generators constructed from time-ordered velocity field snapshots. In particular, for the purpose of coherent pattern identification, we solve the eigenvalue problem for the Koopman generator of this skew-product system with a small amount of diffusion added for regularization. This extends recently developed approximation techniques for Koopman operators of ergodic dynamical systems [22, 23, 24] to the setting of skew-product systems governing the evolution of Lagrangian tracers. In our approach, the regularized operator is elliptic on the product state space (cf. [18, 20], who consider parabolic operators), and has only discrete spectrum [25]; its eigenfunctions correspond to near-invariant or quasiperiodic observables of the skew-product system, and can be visualized as spatiotemporal patterns for a given trajectory on . If the Koopman generator has discrete spectrum, then the eigenfunctions of the regularized generator behave as small-viscosity perturbations of these eigenfunctions. On the other hand, if the Koopman generator has continuous spectrum, the regularized eigenfunctions acquire a singular behavior on the diffusion regularization parameter [26], and generally exhibit intricate spatial structures characteristic of strange eigenmodes.
As in [22, 23, 24], a key ingredient of our approach is to employ kernel algorithms developed for machine learning [27, 28, 29, 30, 31] to build an orthonormal basis for the Hilbert space of the dynamical system without requiring a priori knowledge of the geometry of the state space or the governing equations. Moreover, in the skew-product setting of interest here, we do not require availability of explicit tracer trajectories. This data-driven basis is chosen such that its elements are optimally smooth, in sense of extremizing a Dirichlet energy functional measuring roughness of functions. We also employ this functional to order Koopman eigenfunctions; thus, we select coherent patterns not on the basis of timescale (though at least some of the low-energy eigenfunctions usually capture slow timescales), but on the basis of small roughness, and thus amenability to robust approximation from data.
For the purpose of prediction of observables and probability measures, we employ our data-driven representation of the generator in conjunction with matrix exponentiation algorithms based on Leja polynomial interpolation [32, 33] to approximate the semigroup action. This allows us to perform model-free statistical prediction in skew-product systems with high numerical stability over long stepsizes. We demonstrate with numerical experiments in periodic and non-periodic flows that the skill of these forecasts compares well against Monte Carlo forecasts with the perfect model, despite our models utilizing no prior information about the equations of motion, and without having access to explicit tracer trajectories.
The plan of this paper is as follows. In Section 2, we present our operator-theoretic framework for identification and prediction of coherent patterns. In Section 3, we illustrate this framework in periodic flows with Gaussian streamfunctions where the matrix elements of the generator can be evaluated analytically (i.e., without the need for a data-driven basis). In Section 4, we describe the construction of our data-driven basis and the implementation of the techniques from Section 2 in that basis. In Section 5, we present applications to a class of aperiodic tracer flows driven by Lorenz 96 (L96) systems [34]. The paper ends in Section 6 with concluding remarks and perspectives on future work. Proofs of some of the Lemmas in the main text, details on numerical implementation, and a discussion on the spectral properties of one of the flows studied in Section 3 are included in Appendices. Movies illustrating the spatiotemporal evolution of Koopman eigenfunctions, observables, and probability densities are provided as supplementary online material.
2 Operator-theoretic framework for skew-product systems
2.1 Notation and preliminaries
Consider a continuous-time ergodic dynamical system operating in a closed (i.e., smooth, compact, orientable, and boundaryless) manifold equipped with its Borel -algebra under the smooth map , , preserving a smooth probability measure . This map is generated by a complete vector field such that , and the divergence with respect to is identically zero. Consider also a closed Riemannian manifold equipped with a smooth metric and a smooth probability measure on its Borel -algebra . In what follows, will be the dynamical system governing the evolution of a time-dependent fluid flow and the physical space where this flow takes place. In particular, we consider that there exists a mapping sending to the vector space of vector fields on , divergence-free with respect to , such that is the velocity field corresponding to the state , and vanishes for all . We equip this space with the Hodge inner product, , . We also assume that is an embedding of into (i.e., it possesses a smooth inverse on ), so that an observation provides complete information about the underlying state .
Since is an embedding, inherits a smooth Riemannian metric such that , where is the pushforward map on tangent vectors associated with . Moreover, because is smooth and and are compact, there exists a unique map such that for every and , is a curve on defined for all , and is a tangent vector to that curve at the point . The curve corresponds to the path of a passive tracer released from and advected by the time-dependent velocity field . In what follows, we use the shorthand notation . Note that satisfies the cocycle property, , for all , , and .
Let now be a fixed sampling interval such that the discrete-time dynamical system with , , is also ergodic. We consider that we have available a dataset consisting of time-ordered observations of the velocity field corresponding to the states with . Given such a dataset, and assuming no availability of explicit tracer trajectories, or prior knowledge of the evolution law and the structure of the underlying state space , in what follows we present techniques to (i) identify coherent spatiotemporal patterns associated with the motion of passive tracers in (Section 2.3); (ii) predict the evolution of observables and probability densities defined on the tracers (Section 2.4).
Example 1**.**
Let and be a rotation with frequency , i.e., (note that we abuse notation by using the symbol to represent both a point in and its canonical angle coordinate). Let also be a periodic two-dimensional domain equipped with the canonical (flat) metric , where and are canonical angle coordinates on . We also consider that is equipped with the normalized Haar measure such that . In this setting, any smooth function gives rise to a time-periodic incompressible velocity field with period , given in the coordinates by
[TABLE]
where . That is, is a time-periodic streamfunction, and is identically zero. Note that (1) can be expressed in coordinate-free form using exterior calculus [35, Section 8.2], namely , where is the Hodge-star operator mapping functions to 2-forms, the codifferential mapping 2-forms to 1-forms, and the Riemannian dual of 1-forms. It should also be noted that not every time-periodic, divergence-free velocity field in is of the class in (1); in particular, (1) does not include harmonic vector fields with vanishing vector Laplacian, . Harmonic vector fields cannot be expressed in terms of a continuous streamfunction, but in the case of the 2-torus domain driven by the periodic dynamics on they take the form of a state-dependent (hence, time-periodic) free stream , where and are real-valued functions in . For the class of vector fields in (1), the Riemannian metric induced on by the embedding is given in the canonical angle coordinate of of by , where
[TABLE]
We will return to this class of time-periodic flows in Section 3, where harmonic vector fields will also appear after a Galilean transformation. This example can be generalized to nonperiodic incompressible flows on nonperiodic domains by appropriate modifications of and .
Remark 2**.**
A number of the assumptions stated above can be potentially relaxed. In particular, it should be possible to extend our framework from smooth manifolds to more general topological spaces by replacing Laplace-Beltrami operators and other operators that depend on a Riemannian geometry by appropriate Hilbert-Schmidt kernel integral operators [24]. Such an extension is beyond the scope of the present work, although in what follows we provide justification (Remark 23) and numerical evidence (Section 5) that our framework is applicable when is not a smooth manifold. Moreover, the embedding assumption on can be generically relaxed by performing delay-coordinate maps [36] on the velocity field data, and using that data as representatives of the underlying dynamical states in .
2.2 Extended state space and the associated Koopman and Perron-Frobenius operators
Our approach for identifying and predicting coherent patterns relies on data-driven approximations of Koopman and Perron-Frobenius operators characterizing the evolution of observables and probability measures, respectively, on the product manifold . Intuitively, we think of as a “space-time manifold” with playing the role of physical space where the motion of tracers takes place and the role of “time” where the velocity field evolves. On , the dynamics is autonomous (though not necessarily ergodic), and is governed by the flow , , having the skew-product form . This flow preserves the product measure defined on the Borel -algebra on . In what follows, we will employ the usual notation for the inner product between complex-valued functions , and use whenever convenient the abbreviation . We also abbreviate and , and use and to denote the inner products of these spaces (defined analogously to ), respectively.
Given an arbitrary smooth Riemannian metric on , we use the notation to represent the weighted Hodge inner product associated with and on vector fields . We also denote the canonical projection maps from into its factors by and . Due to the product structure of , we have the decomposition , , , where and are the pushforward maps on tangent vectors associated with and , respectively. In particular, every tangent vector can be uniquely decomposed as with and .
Given a point , the curve describes the joint evolution of the base dynamics initialized at state and the passive tracer released at the point and advected by the time-dependent velocity field . Moreover, tangent to the family of curves is a vector field such that . Due to the skew-product structure of , we have where and are the canonical lifts of and on , respectively; that is, and . Because preserves , has vanishing divergence with respect to that measure, . Intuitively, one can think of as an analog of the material derivative in fluid dynamics in measuring the rate of change of observables on fluid parcels. Noting that (where is the directional derivative associated with a time-dependent vector field in ), we can also identify with and with .
Next, we introduce Koopman and Perron-Frobenius operators governing the evolution of observables and probability densities, respectively, on and . In the case of , we consider complex-valued observables in , and define as usual the group of Koopman operators , , acting on observables by composition with , i.e., [21, 37]. Since is an invariant measure, is unitary, , and is a Perron-Frobenius (transfer) operator governing the evolution of complex measures on with densities in ; i.e., if is such a measure with density , is the density of the measure with respect to . By Stone’s theorem, the unitary group is generated by a skew-adjoint operator with dense domain , such that , and if (i.e., is a skew-adjoint extension of viewed as a map from to ). Moreover, there exists a unique spectral measure , where is the Borel -algebra on and the set of orthogonal projections on , such that and . Similarly, the Perron-Frobenius group is generated by , and we have .
We employ an analogous construction in the case of the dynamics on the product space . In this case, we consider complex-valued observables in the Hilbert space associated with the invariant product measure . When studying individual dynamical trajectories on , we can think of functions (and thus equivalence classes of functions in ) as spatiotemporal patterns. Specifically, given an initial state , the spatiotemporal pattern associated with is given by . In this setting, the group of Koopman operators , , acts on via composition with the skew-product map, , so that corresponds to the value of at the point reached at time by a passive tracer released at the point and advected by the time-dependent velocity field . Similarly, the group of Perron-Frobenius operators governs the evolution of densities of measures on under , so that if is a measure with density , is the density of . In the case that is a probability measure, corresponds to a time-dependent probability density characterizing uncertainty with respect to both tracer position in and the underlying state of the velocity field in . Moreover, the function is the density of a time-dependent probability measure on which corresponds to a marginalization of over . Note that unlike , does not evolve autonomously.
Consider now the generator of the unitary group . In direct analogy with , the generator of this group is a skew-adjoint extension of the vector field defined on a dense domain , and the skew-adjointness of is a consequence of the fact that preserves . By Stone’s theorem, we have
[TABLE]
where is a unique spectral measure associated with taking values on the set of orthogonal projections on . Thus, if is known, it is possible to compute and , and hence predict the evolution of arbitrary observables and measures for the skew-product tracer system.
2.3 Identification of coherent spatiotemporal patterns
2.3.1 Coherent patterns as eigenfunctions of the generator
We identify coherent spatiotemporal patterns in time-dependent fluid flows through approximate eigenfunctions of the Koopman generator at small corresponding eigenvalues. To motivate our approach, consider the eigenvalue problem , and suppose that this problem has a nonconstant solution at eigenvalue (note that any -a.e. constant function is an eigenfunction of at eigenvalue zero). Because (by (3)), such an eigenfunction is preserved on Lagrangian tracers. Correspondingly, the level sets of create an ergodic quotient [37, 38] of ; that is, a partition into codimension 1 invariant sets on which tracers on are trapped. More specifically, if , then for all since .
In the case that has eigenfunctions with nonzero corresponding eigenvalues, then these eigenfunctions also provide a useful notion of coherent spatiotemporal patterns that vary periodically on the tracers. In particular, because is skew-adjoint, its eigenvalues are purely imaginary, and measures an intrinsic oscillatory frequency associated with eigenfunction . It therefore follows from (3) that , or, equivalently, , where is the rotation map on the complex plane such that . In this case, the sets are not invariant, but the collection (which forms a partition of ) is invariant in the sense that . Such a partition is known as a harmonic partition [38].
Remark 3** (Group structure of Koopman eigenvalues and eigenfunctions).**
Given any two eigenfunctions of with corresponding eigenvalues and , it follows from the Leibniz rule that , so that is an eigenfunction at eigenvalue . Thus, the sets of such eigenvalues and eigenfunctions form commutative groups under addition of complex numbers and multiplication of smooth functions, respectively. In particular, one can generate these sets from any maximal subset of rationally independent eigenvalues and their corresponding eigenfunctions. The number of elements of such sets is at most equal to the dimension of .
The preceding discussion suggests that a useful definition for coherent spatiotemporal patterns on is through eigenfunctions of at zero or small eigenvalue. However, despite its simplicity and attractive theoretical properties, this definition is likely to suffer from at least two major shortcomings, namely that (i) for systems of sufficient complexity (i.e., weak mixing systems) will have no nonconstant eigenfunctions, and (ii) even if has nonconstant eigenfunctions, the set of these eigenfunctions will generally include functions of high roughness which will adversely affect the conditioning of numerical schemes for the eigenvalue problem for .
Example 4**.**
To illustrate the latter point, suppose that the dynamics on and were both circle rotations with , , and an irrational number. Then, is a 2-torus, and the dynamics is ergodic for the Haar measure. Moreover, the eigenvalues and eigenfunctions of are given by and , respectively, with . Because is irrational, the set is dense in the real line; that is, given an eigenfrequency (including ) and a positive number , there exist infinitely many eigenfrequencies in the interval . Given, however, any smooth Riemannian metric on , most of these frequencies will correspond to highly oscillatory eigenfunctions with large gradient. For instance, if is the flat metric on , then , and while can be made arbitrarily close to zero using arbitrarily large (and opposite signed) , the gradient of the corresponding eigenfunctions will be arbitrarily large.
The behavior in Example 4 is typical of other systems too, and renders the approximation of the eigenvalues and eigenfunctions of the raw generator by numerical algorithms problematic—in particular, the classical spectral Galerkin methods for eigenvalue problems on Hilbert spaces [39] require certain coercivity (ellipticity) properties to hold which are not satisfied by skew-adjoint operators such as . Thus, as done elsewhere in the Koopman and Perron-Frobenius operator literature [40, 11, 22, 23, 24, 20], we will regularize the generator by adding a small amount of diffusion that renders the spectrum discrete and its approximation by spectral Galerkin methods well-posed.
2.3.2 Laplace-Beltrami operators on and the associated eigenfunction bases
We employ the framework developed in [22, 23, 24] which involves regularizing the generator by means of the Laplace-Beltrami operator associated with a Riemannian metric on chosen such that it has compatible volume form to the invariant measure, . This operator is defined as , where and are the gradient and divergence operators associated with , respectively. Note that since , we have and for any function and any vector field , . As a result, is a symmetric operator, and it is a standard result in analysis on Riemannian manifolds [41] that there exists an orthonormal basis of consisting of eigenfunctions of (an appropriate self-adjoint extension of) corresponding to the eigenvalues , where . An important property of these eigenfunctions is that they are stationary points of the Rayleigh quotient associated with the Dirichlet energy functional , and . Intuitively, the Dirichlet energy provides a measure of roughness of functions on associated with the Riemannian metric and the invariant measure . Fixing a parameter , the set has minimal roughness with respect to among all -element orthonormal sets in . In what follows, we will use the eigenfunctions , or (in Section 4) approximate eigenfunctions computed from data via kernel algorithms, to construct Galerkin approximation spaces for a regularized generator (to be defined in Section 2.3.3).
Let and be the densities of the probability measures and relative to the Riemannian measures associated and , respectively. In particular, represents the sampling density on with respect to the invariant measure of the dynamics and the observation map . Let also and be the dimensions of and , respectively. In the setting of interest involving the -dimensional product space with , it is natural to construct the metric as the sum , where and are pullbacks of the conformally transformed metrics and on and , respectively. By construction, and , and moreover for any we have , where and are unique vector fields in and , respectively. With this choice of metric, the Laplace-Beltrami operator on becomes , where and . Note that and are lifts of the Laplace-Beltrami operators and associated with and , respectively, in the sense that and for any and . Moreover, and commute, and the eigenfunction basis associated with has the tensor product form , where and are eigenfunctions of and , respectively. Here, the index ordering , is chosen such that and the Dirichlet energies form a non-decreasing sequence with in accordance with our convention. In particular, denoting the eigenvalues corresponding to and by and , respectively, we have , where is the eigenvalue corresponding to . For later convenience, we quote here formulas relating the gradient, divergence, and Laplace-Beltrami operators associated with to the corresponding operators with respect to the ambient space metric , viz.,
[TABLE]
where is the “inverse metric” associated with , and and are an arbitrary smooth function and a vector field on , respectively.
Example 5**.**
In the case of the periodic flows discussed in Example 1, the sampling density and the conformally transformed metric on are given by and , where is a canonical angle coordinate on . Moreover, since is equipped with the canonical flat metric, the density is uniform, , and the metric is also flat; that is, , where and are canonical angle coordinates on . With these Riemannian metrics, the Laplace-Beltrami eigenfunctions are the canonical Fourier functions, and , where , , and are integers, and . The corresponding eigenvalues are and so that .
Besides , in what follows we will make use of Sobolev spaces [41] on of higher regularity. In particular, we consider the space defined the completion of with respect to the norm induced by the inner product with and . Note that the eigenfunctions of have unbounded norm with , and therefore is not a basis of . On the other hand, the rescaled eigenfunctions defined as and have bounded norm, , , and therefore is a basis of . Moreover, for a function expanded in this basis, the Dirichlet energy with respect to can be conveniently computed from the norm of the expansion coefficient sequence excluding ; i.e., . In Section 2.3.3, we will employ the basis for spectral Galerkin approximation of the eigenvalues and eigenfunctions of our regularized Koopman generator. We will also make occasional use of the Sobolev space defined as the completion of with respect to the norm induced by the inner product .
2.3.3 Coherent spatiotemporal patterns with regularization
Following [22, 23, 24], we employ the Laplace-Beltrami operator introduced in Section 2.3.2 to construct regularized generators and defined as
[TABLE]
Proposition 6**.**
The operators have the following properties:
(i) They are dissipative, , for all , and therefore have dissipative closures, .
(ii) They form a dual pair, in the sense that extends ; i.e., and .
(iii) There exists a maximal, closed, dissipative operator with a maximal, dissipative adjoint such that
[TABLE]
(iv) Given a function , we have , where is a constant independent of ; thus, and .
(v) The spectra and are purely discrete; i.e., if and only if is an eigenvalue of , and if and only if is an eigenvalue of . Moreover, if and only if , and if is an eigenfunction of at eigenvalue , then is an eigenfunction of at eigenvalue .
A proof of Proposition 6 can be found in A.1. It follows from Proposition 6(iii) and a classical result due to Phillips [42] that is the generator of a strongly continuous contraction semigroup on ; that is, , and for all , , , and the map is continuous in the strong operator topology. Similarly, is the generator of the adjoint semigroup, , which is also contractive and strongly continuous. Note that, in general, and are nonnormal operators, and therefore are not guaranteed to be spectrally decomposable as and are. In Section 2.4 we will employ these semigroups for prediction of observables and probability densities on Lagrangian tracers, respectively, but for now we focus on using the spectral properties of to recover coherent spatiotemporal patterns. In particular, the fact that consists only of eigenvalues, in conjunction with the discussion in Section 2.3.1, motivates the following definition of coherent spatiotemporal patterns for Lagrangian tracer systems driven by ergodic flows (and for more general skew-product systems):
Definition 7** (Coherent spatiotemporal patterns).**
For fixed , identify coherent spatiotemporal patterns in as the eigenfunctions of with minimal Dirichlet energy with respect to the metric introduced in Section 2.3.2. That is, obtain from the solutions of the problem , , ordered such that , where the are Dirichlet energies given by
[TABLE]
The quantities and correspond to the growth rate (guaranteed non-positive by Proposition 6(iii)) and oscillatory frequency of eigenfunction .
Intuitively, for small we think of the from Definition 7 as approximate eigenfunctions of the raw generator . In the special case that has a pure point spectrum, i.e., there exists a smooth orthonormal basis of consisting of eigenfunctions of , that intuition can be justified using regular perturbation expansions [23]. In particular, it can be verified through heuristic asymptotics that for each with corresponding eigenvalue , there exists an eigenfunction of with corresponding eigenvalue and a smooth function orthogonal to such that , , and . Thus, in systems with pure point spectrum, the main effect of diffusion in is to suppress the pathological eigenfunctions of with large Dirichlet energy (see Example 4), and as the eigenvalues of converge to the eigenvalues of . In fact, in systems with pure point spectrum one can additionally construct a metric (using Takens delay embeddings) such that and commute [23, 24]. In such cases, and have common eigenfunctions, and we have , and without approximation.
On the other hand, in systems with mixed or continuous spectra (or rough Koopman eigenfunctions if the vector field is non-smooth), the relationship between the spectra of and becomes significantly more complicated. At the very least, the fact that always has eigenfunctions whereas may have no eigenfunctions (if the system is weak-mixing) raises questions about the physical significance of the coherent patterns defined according to Definition 7. In fact, Constantin et al. [26] show that even if has eigenfunctions, but these eigenfunctions are not in , a phenomenon known as dissipation enhancement occurs whereby adding diffusion to causes a greater than change in the real part of the eigenvalues of controlling the decay of the evolution semigroup . More specifically, expressed in our notation, the results of [26] imply that if has no eigenfunctions in , then for any there exists such that for all . In other words, phenomena such as dissipation enhancement suggest that in systems of sufficiently high complexity the choice of the diffusion operator in (5) may influence significantly the properties of the coherent patterns recovered via eigenfunctions of . Here, besides this note for caution, we do not pursue the topic of the choice of (which exists more broadly than our study of Lagrangian tracer systems) in the presence of mixed or continuous spectrum further. However, the experiments in Sections 3 and 5 suggest our choice of as the Laplace-Beltrami operator associated with the metric leads to reasonable coherent patterns even in systems where the spectrum is not purely discrete.
2.3.4 Galerkin approximation
We compute weak solutions to the eigenvalue problem for using a spectral Galerkin method. To pass to a weak (variational) formulation, we take the inner product of both sides of the equation , , by an arbitrary test function , integrate by parts the diffusion term, , and introduce the sesquilinear forms
[TABLE]
Note that is the Dirichlet form associated with the heat semigroup generated by , and that by skew-symmetry of .
By the Cauchy-Schwartz inequality and the fact that , these forms satisfy the bounds
[TABLE]
where . Therefore, and can be uniquely extended to bounded sesquilinear forms on . Moreover, one can verify that has the coercivity property , where is any function orthogonal to constant functions and . Together, the boundedness and coercivity of and imply that the following variational eigenvalue problem is well posed [39]:
Definition 8** (variational formulation of eigenvalue problem for coherent patterns).**
Find and such that for all the equality holds, where and are the sesquilinear forms on defined in (6).
We compute approximate solutions of the eigenvalue problem in Definition 8 through a spectral Galerkin method expressed in finite-dimensional subspaces of spanned by the basis introduced in Section 2.3.2. In particular, introducing the -dimensional subspaces , we seek solutions of the equation for and . This is equivalent to solving the matrix generalized eigenvalue problem
[TABLE]
where is an -dimensional complex column vector, and and are matrices with elements
[TABLE]
Given a generalized eigenvector from (8), the corresponding approximate eigenfunction of is given by . Moreover, as stated in Section 2.3.2, the Dirichlet energy of this eigenfunction is given by . In accordance with Definition 7, we order all such solutions in order of increasing Dirichlet energy , and normalize each solution to unit norm, . Convergence of the approximate solutions to the solutions of the full eigenvalue problem in Definition 8 follows from spectral approximation theory for variational eigenvalue problems [39]. A discussion on the numerical solution of (8) can be found in B.4.
Remark 9**.**
Since for every , we can determine the Dirichlet energies of the eigenfunctions directly from the real parts of the corresponding eigenvalues, i.e., . This is useful when solving (8) via iterative solvers with an option to compute subsets of eigenvalues with maximal real parts (e.g., the ARPACK library [43], which is also utilized by our Matlab code). As discussed in B.4, in practice we find that a more rapid numerical convergence occurs if we select eigenvalues with minimal moduli as opposed to maximal real parts. So long that a sufficiently large number of eigenfunctions is computed, these two approaches should yield consistent results for the leading solutions of (8) ordered in order of increasing . Given the computational cost of solving (8) at high (which is typical for the tensor product basis used in this work), hereafter we always compute eigenvalues and eigenfunctions using the minimal option. This approach carries a risk of failing to detect certain eigenvalues of high oscillatory frequency but low Dirichlet energy, and indeed in Section 3 this behavior will occur (although the missed eigenfunctions will turn out to be spatially constant on , and therefore trivial from the point of view of coherent spatiotemporal patterns).
Next, we comment on the numerical conditioning of our scheme. By construction of the basis, is an diagonal matrix with and for . In particular, the operator norm of this matrix is independent of (it is equal to 1 for ). Moreover, it follows from the left-hand equation in (7) and the fact that the have bounded norm that the norm of the matrix can also be bounded by an -independent constant. These facts imply that the norm of is bounded by an -independent constant. In the case of , we can bound its condition number by where the constant is independent of , and ; this follows from the facts that , , and Weyl’s law for the asymptotic growth of the Laplace-Beltrami eigenvalues on a closed -dimensional manifold as . Together, these properties of and ensure the good numerical conditioning of the eigenvalue problem in (8).
2.4 Prediction of observables and probability densities
As stated in Section 2.1, the Koopman and Perron-Frobenius groups of operators and respectively govern the evolution of observables and densities in the space associated with the skew-product tracer system. In this Section, we present a technique for approximating the action of these operators on observables and densities in an analytic or class.
2.4.1 Computing the semigroup action for bounded analytic and observables
Consider first the problem of computing for a general observable . As in the case of the eigenvalue problem for coherent patterns (Section 2.3.3), we relax this problem to the problem of computing the action of the contraction semigroup generated by the regularized operator in proposition 6. An advantage of using this semigroup is numerical stability due to the fact that is dissipative; a disadvantage is that is nonnormal (unlike the raw generator ), and therefore not necessarily decomposable via a spectral integral analogous to (3). A difficulty common to both and is that they are unbounded operators, which means that we cannot evaluate the corresponding and operators via series expansions for arbitrary observables in . To address this issue we restrict attention to a class of analytic observables with respect to [44, Section VIII.5]. Specifically, introducing the space , we consider the space , , consisting of all observables satisfying the bound for all and some constant . We also define . Observables in are known as bounded analytic vectors [45, 46]. For such observables, the sum converges for all , and therefore the action of the semigroup can be evaluated by means of the series expansion
[TABLE]
We now further restrict attention to bandlimited observables in which can be represented through finite expansions of the eigenfunctions of . In particular, let be the -dimensional subspace of spanned by the leading eigenfunctions of (ordered in order of increasing Dirichlet energy in accordance with Definition 7). Let also . We then have:
Lemma 10**.**
(i) The subspace lies in .
(ii) If the eigenfunctions , form a basis of , and the corresponding eigenvalues have no accumulation points, then for every there exists such that lies in .
A proof of this Lemma can be found in A.2. Lemma 10 implies that for observables we can evaluate the semigroup action via the series expansion in (9) for all ; explicitly, we have
[TABLE]
Moreover, if the assumptions in part (ii) of the Lemma hold, then (10) holds for any observable in (which is in turn a dense subspace of ), again for for all . While we do not have results allowing us to determine when these assumptions hold for general , there are specific cases where they are satisfied. In particular, if the Koopman generator has pure point spectrum, then it can be arranged (e.g., using delay-coordinate maps [23, 24]) that and commute; in that case, has an orthonormal eigenbasis with no accumulation points in its eigenvalues (in fact, the eigenfunctions of are also eigenfunctions of ).
Next, we consider how to compute or approximate (9) when for some . First, it follows immediately from (10) that if some eigenvalues and eigenfunctions of are available, say the first , then one can approximate , , by . This is equivalent to introducing the orthogonal projectors and the spectrally truncated generator , and computing
[TABLE]
where it is obvious that for all if .
An expansion analogous to (11) was used in [23] for prediction of observables of ergodic dynamical systems. For the skew-product systems studied here, however, computing enough eigenfunctions of so that (11) is sufficiently accurate may not be feasible (for instance, due to the large number of required basis elements in the tensor product basis ). In such situations, one may opt instead to use orthogonal projections to the subspaces spanned by eigenfunctions of , which may be “easier” to compute. In particular, the eigenfunctions of are also eigenfunctions of the corresponding heat kernel which is a compact operator [41]; a fact which will be useful when building data-driven bases in Section 4 ahead. With this in mind, an alternative to (11) is to define for some , and compute
[TABLE]
where can now be any observable in (as opposed to a bandlimited observable). As with (11), at fixed , the sum over in (12) converges for all since has finite rank. However, unless rather restrictive conditions hold (e.g., that the eigenspaces of are invariant under ), the convergence of to is not guaranteed to be uniform in . Instead, classical results on semidiscrete approximation schemes for linear parabolic partial differential equations [47] lead to the weaker convergence result
[TABLE]
where is a constant and converges to zero as . Note that the weakness of the bound at short times stems from the fact that an arbitrary may not be in the domain of , but lies in for all .
For the remainder of the paper, we restrict attention to approximations of via (12) as opposed to (11). This approximation can be evaluated by computing the matrix exponential of the generator matrix with elements , . Specifically, given , we have , and expanding , , it follows that , where the expansion coefficients can be computed in vector form using
[TABLE]
It should be noted that due to the tensor product nature of the Laplace-Beltrami eigenfunction basis for , the spectral truncation parameter can be very large (e.g., in the experiments in Section 5, will be as large as ); therefore, efficient evaluation of becomes crucial.
2.4.2 Approximation via Leja interpolation
The problem of computing the action of matrix exponentials on vectors has been studied extensively [48, 49, 50]. Here, we approximate the semigroup action in (14) using the method of Caliari et al. [32] and Kandolf et al. [33] which approximates the action of the matrix exponential on vectors via polynomial interpolation at Leja nodes. Specifically, the method in [32, 33] approximates by where is a Newton interpolating polynomial of degree associated with a sequence , , of Leja nodes in the complex plane. By employing an interpolating polynomial, the method does not require explicit evaluation of the matrix exponential; this is important from a computational efficiency standpoint since, for instance, will generally be a full matrix even if is sparse. Moreover, the properties of Leja interpolation nodes ensure that the approximation remains well-conditioned at large . The properties of Leja sequences also allow for efficient iterative approximation refinement by increasing until an error tolerance is met, using a stopping criterion derived from an a posteriori error estimate to terminate the iterations [33].
The method in [32, 33] utilizes Leja nodes on the line segment in the complex plane, where , , and are real numbers that determine the corners of a rectangle bounding the spectrum of . Caliari et al. [32] identify this rectangle by applying Gershgorin’s disk theorem to the symmetric and antisymmetric parts of the matrix ; we also follow a similar approach here. In particular, in the case of the generator matrix associated with the semigroup , the symmetric part is a diagonal matrix with diagonal elements proportional to the Laplace-Beltrami eigenvalues (Dirichlet energies). Thus, we can set and . To place a bound on the imaginary part of the spectrum, we apply Gershgorin’s theorem to the antisymmetric part , where the matrix elements are determined by the matrix representation of the Koopman generator in the Laplace-Beltrami eigenfunction basis of . Applied to this matrix, Gershgorin’s theorem gives .
Once the interval has been established, the first Leja point is selected arbitrarily in that interval, and the remaining points are computed iteratively, setting for . Associated with the Leja points are the divided differences , which can also be computed iteratively (and stably [51]), using
[TABLE]
The Newton interpolating polynomial associated with the Leja nodes is then given by
[TABLE]
Besides the fact that the -th order approximation can be easily computed from , other attractive properties of the Leja method include good numerical conditioning and superlinear convergence to ; additional details on these properties can be found in [32, 33]. Details on the numerical implementation used in the experiments in Sections 3 and 5 ahead are included in B.1.
In applications, we are often interested in evaluating at the times with and a fixed sampling interval (not necessarily equal to the sampling interval in the training data). In such cases, we compute the quantities iteratively using
[TABLE]
In this procedure, the normalization by ensures that the vectors have constant norm, as would be expected for unitary evolution. In effect, this normalization step can be thought of as an inflation to compensate for contraction due to diffusion regularization and spectral truncation associated with the operator in (12). In practice, we have observed that this step is not particularly important, but can sometimes lead to a modest increase of prediction skill.
2.4.3 Computing the semigroup action on probability densities
We now turn to the approximation of the Perron-Frobenius group action on probability measures with densities relative to the invariant measure in . As discussed in Section 2.1, the evolution of such a density is governed by the adjoint of the Koopman group, (i.e., the Perron-Frobenius group), and in order to approximate the action on this group on functions we proceed in an analogous way to the case of in Sections 2.4.1 and 2.4.2. That is, we approximate by , where is the adjoint semigroup generated by , and further approximate this quantity by projection onto the finite-dimensional subspaces , viz., , where . The error of this approximation can be estimated by an analogous expression to (13). Expanding and with and , we obtain an analogous expression to (14), namely, , where and are -dimensional column vectors containing the expansion coefficients of and , respectively, and is the adjoint of the generator matrix. We compute the action of the matrix exponential on using Leja interpolation as described in Section 2.4.2.
3 Semi-analytical examples
In this Section, we apply the framework for coherent pattern identification and prediction presented in Section 2 to time-periodic incompressible flows on a doubly-periodic domain, . According to Examples 1 and 5, the dynamical state space for this class of flows is the circle, , and as a result the product space is the 3-torus. Moreover, the Laplace-Beltrami eigenfunction basis of associated with the metric consists of the Fourier functions , where , , , and and are canonical angle coordinates for and , respectively. The Laplace-Beltrami eigenvalue corresponding to is . Note that in this Section we label the basis elements and the eigenvalues using three integer indices (as opposed to the single index used in Section 2.3.2), as this notation is more convenient for expressing the matrix elements of the generator. Moreover, we use the symbols , , and to represent the spectral truncation parameters for and each of the factors of ; that is, our finite-dimensional approximation spaces have dimension , and are spanned by with , , and , or, as appropriate, their normalized counterparts for the subspaces of . Following Definition 7, whenever we discuss numerical Koopman eigenvalues and eigenfunctions , , we order these eigenvalues and eigenfunctions in order of increasing Dirichlet energy , though as stated in Remark 9, in practice we compute subsets of the spectrum consisting of eigenvalues of minimal modulus using iterative solvers.
In the examples that follow, we consider the same periodic dynamics on as in Example 1 (i.e., ), and set the streamfunction to various linear combinations of circular Gaussian functions. With this choice, the matrix elements of the generator can evaluated in closed form using integral properties of modified Bessel functions. Thus, in this Section, all errors in the computation of eigenfunctions and evolution of observables and densities are caused by projection to finite-dimensional subspaces of , as opposed to ergodic averaging errors due to finite numbers of samples and/or errors associated with the construction of the Laplace-Beltrami eigenfunction basis. We will take up the problem of constructing data-driven approximations to the Laplace-Beltrami eigenfunction basis in Section 4.
3.1 Moving Gaussian vortex in a periodic domain
In our first example, we consider the streamfunction defined as
[TABLE]
where is a positive parameter. Thus, for each , corresponds (up to normalization) to a von Mises density on , which is the analog of a Gaussian density in a periodic domain. The corresponding incompressible velocity field from (1) generates a vortical flow centered at , illustrated in Fig. 1. Since and is zero when and , the velocity field vanishes at the center of the vortex as well as at three points away from the vortex center, leading to a ring-like region around the center where the velocity field is concentrated. Note that the radial extent of the ring decreases with increasing . Moreover, because the dynamics on evolve according to , the vortex moves without deformation at fixed speed along the direction.
For our purposes, a particularly useful property of Gaussian streamfunctions such as (16) is that their integrals against Fourier functions can be evaluated in closed form with the help of the identity
[TABLE]
where is the modified Bessel function of the first kind and order . In particular, recall the decomposition in Section 2.1 of the vector field of the skew-product system in terms of the vector fields and associated with the dynamics on and , respectively. We have
[TABLE]
and using (17) we arrive, after some algebra, at the result
[TABLE]
In addition, we have
[TABLE]
Together, (19) and (20), in conjunction with the fact that and the definitions of the normalized eigenfunction basis for (see Section 2.3.2), are sufficient to evaluate all sesquilinear forms and operator matrix elements appearing in the schemes of Sections 2.3 and 2.4. In what follows, we present results for the flow with and ; with this choice of flow parameters the maximal velocity field norm is , and occurs at a radial distance of from the vortex center (all with respect to the canonical flat metric on ).
3.1.1 Koopman eigenvalues and eigenfunctions
Before presenting our results, it is useful to comment on the general spectral properties of the generator for our choice of flow parameters. As discussed in C, the spectrum of always contains eigenvalues which can be grouped into (1) integer multiples of , i.e., with , corresponding to eigenfunctions that are constant on the spatial domain (); (2) zero eigenvalues corresponding to eigenfunctions constant on the streamlines of the velocity field in a comoving frame with the vortex center; (3) linear combinations of the eigenvalues in class 1 and 2 (i.e., integer multiples of ) corresponding to products of the respective eigenfunctions. In addition, the spectrum of has a continuous component associated with functions in that are nonconstant on the streamlines in class 2.
Eigenfunctions of class 1 exist in all time-periodic flows with the associated skew-product dynamics of Section 2.1; this is because is trivially projectible under . In particular, is equal to the vector field of the periodic dynamics on , and (a skew-adjoint extension of) the latter has pure point spectrum. Eigenfunctions of this class can therefore be expressed as pullbacks of eigenfunctions of the generator dynamics on ; that is, we have , where and . We denote the closed subspace of spanned by these eigenfunctions by (here, the symbol stands for discrete spectrum). Eigenfunctions of class 2 exist because admits a submersion onto a circle where the image of the vector field is identically zero. This submersion can be intuitively thought of as a Galilean transformation to a frame of reference where the vortex is stationary, followed by projection transverse to the streamlines in that frame. As with the eigenfunctions of class 1, the eigenfunctions of class 2 lie in a closed subspace and can be expressed as pullbacks , where are orthonormal functions on equipped with the normalized Lebesgue measure. Eigenfunctions of class 3 are a direct consequence of the group structure of the Koopman eigenvalues and eigenfunctions stated in Remark 3.
Together, eigenfunctions in classes 1–3 span the tensor product space , which is closed and contains all eigenfunctions of . Note that , , and are all invariant under the Koopman group . Thus, we have an invariant orthogonal decomposition , where the is associated with the continuous spectrum of the generator. The existence of the latter can be justified from the fact that the flow in the frame comoving with the center of the vortex is conjugate to a skew rotation on , which is known to have continuous spectrum [52, 53]. Additional details on the properties of the spectrum of can be found in C.
Figure 2 shows snapshots of representative numerical Koopman eigenfunctions for the moving vortex Gaussian flow, computed using the spectral truncation parameters (amounting to a total of basis functions) and the diffusion regularization parameter . To obtain these results, we computed a total of 51 eigenfunctions using Matlab’s eigs solver as described in B.4, targeting eigenvalues of minimal modulus. The eigenfunctions shown in Fig. 2, , , , , , , , and , are also visualized as spatiotemporal patterns in Movie 1. The eigenvalues and Dirichlet energies corresponding to these eigenfunctions are listed in Table 1. Note that the come in complex-conjugate pairs, so we do not show them for consecutive .
The numerical eigenfunctions , , , , , and , are in good agreement with what expected from theory for the Koopman eigenfunctions in the subspace . In particular, comparing Figs. 1(b) and 2, it is evident that the numerical eigenfunctions are to a good approximation constant on the streamlines in the frame comoving with the vortex center, and according to Table 1 the imaginary parts of their corresponding eigenvalues are numerically close to zero. Moreover, the real parts of the eigenvalues are numerically close to ; this is consistent with what expected from heuristic perturbation expansions of the eigenvalues and eigenfunctions of the operator with respect to [23].
Qualitatively, it appears that the numerical eigenfunctions in can be assigned a positive integer “wavenumber” equal to half of the number of zeros of their real and imaginary parts along the line parallel to the coordinate line of passing through the vortex center—as expected, the Dirichlet energy of these eigenfunctions is an increasing function of . This behavior (in conjunction with the structure of and the fact that the eigenfunctions in question are constant on the streamlines) suggests that the numerical in this class could be approximating pullbacks of eigenfunctions of a Laplace-Beltrami operator on the base space of the submersion . That is, a reasonable hypothesis could be that for every we have
[TABLE]
While such a scenario would be desirable, in the sense that the effect of diffusion would lead to no perturbation of the eigenfunctions of compared to the eigenfunctions of , fundamental results from the theory of harmonic maps on manifolds suggest that it is unlikely to hold exactly in practice. In particular, it is known that a necessary and sufficient condition for the relationship
[TABLE]
to hold for arbitrary is that the fibers of the submersion are minimal, i.e., they are surfaces of vanishing mean curvature [54, 55, 56]. Clearly, the latter is not the case given the nonconstant curvature of the streamlines in Fig. 1(b), so we cannot expect (21) to hold for general . Nevertheless, despite that our diffusion operator is not tailored to the flow under study (in the sense of (21) failing to hold), the numerical results in Fig. 2 and Table 1 demonstrate that it is still possible to obtain high-quality numerical Koopman eigenfunctions via the scheme of Section 2.3.
Next, we examine the effects of the continuous part of the spectrum of to the numerically computed eigenfunctions of . As shown in Fig. 2 and Movie 1, the numerical eigenfunctions of include functions (e.g., and shown here) which are highly localized around individual streamlines in the comoving frame with the vortex center, and also exhibit abrupt changes along those streamlines. As a result, the numerical eigenfunctions of this class also have large Dirichlet energy (see Table 1). We believe that such eigenfunctions are remnants of the continuous spectrum of , which appear numerically due to finite spectral truncation in the Galerkin method for the eigenvalues and eigenvectors of . In particular, while does not have nonconstant eigenfunctions on the streamlines (i.e., it does not have eigenfunctions in the Hilbert subspace ), it is nevertheless possible to construct a Koopman eigenvalue problem on a suitable space of distributions [57], such that points in the continuous spectrum of are eigenvalues corresponding to eigendistributions. Such eigendistributions are supported on sets of zero -measure (e.g., finite collections of streamlines), and therefore cannot be represented by functions, but they still can be approximated by functions by molification. Such molified distributions are not expected to solve the eigenvalue problem exactly, but for a given spectral truncation parameter it should be possible to arrange that the residual is -orthogonal to the finite-dimensional approximation space , making a solution of the numerical Koopman eigenvalue problem. With increasing , would have to be concentrated on a set of increasingly small -measure (otherwise, the residual would fail to be orthogonal to ), leading to an unbounded increase of the Dirichlet energy with . The presence of numerical eigenfunctions such as and , which are highly concentrated near individual streamlines (or finite collections of streamlines) is consistent with this hypothesis. In separate calculations, we have also confirmed that increasing indeed causes these eigenfunctions to become increasingly concentrated and appear deeper in the spectrum of (ordered with respect to Dirichlet energy).
Before closing this section, we note that the results presented above do not include eigenvalues and eigenfunctions of class 1. This is likely due to the fact that we have numerically computed only a small (51-element) subset of the spectrum targeting eigenvalues of minimal absolute value—as described in Remark 9, this can cause us to miss eigenfunctions with low Dirichlet energy () but high oscillatory frequency (), such as the class 1 eigenfunctions associated with the periodic dynamics on . In separate calculations, we have confirmed that these eigenfunctions can be found by directly targeting that maximal portion of the spectrum in eigs, at the expense of slower numerical convergence (see B.4). Since the spatial patterns captured by these eigenfunctions are trivial we do not discuss them here.
3.1.2 Prediction of observables and probability densities
We now apply the techniques of Section 2.4 to perform prediction of observables and probability densities for the moving vortex flow, working with the same flow, spectral truncation, and diffusion regularization parameters as in Section 3.1.1. In what follows, we focus on the observables
[TABLE]
where . Our interest in these observables stems from the fact that knowledge of and is sufficient to determine the position of Lagrangian tracers advected by the flow at forecast time . Specifically, using the notation for the canonical angle coordinates at time of a tracer released at time at the point when the state of the time-periodic velocity field is , we have
[TABLE]
and therefore we can recover and from and , respectively.
Figure 3 shows snapshots of the evolution of obtained via the operator-theoretic approach of Section 2.4.1 and explicit solution of the ordinary differential equations (ODEs) governing the evolution of Lagrangian tracers under the time-periodic streamfunction in (16). This evolution is also visualized more directly in Movie 2. To obtain the operator-theoretic results we used Leja interpolation for matrix exponentiation as described in Section 2.4.2 and B.1 with a relative error tolerance of and forecast timestep . We integrated the tracer ODE system using Matlab’s ode45 solver with a relative error tolerance, outputting the solution every time units.
As is evident from both the operator-theoretic and ODE results the evolution of the and coordinates (hence, the and observables) is qualitatively different, with the former exhibiting significantly more mixing than the latter. This can be understood from the facts that varies predominantly transverse to the streamlines in Fig. 1(b), and thus projects more strongly to the discrete spectrum subspace , whereas projects more strongly to the continuous spectrum subspace as it varies predominantly along the streamlines. More qualitatively, after each revolution of the vortex center around the periodic domain, a tracer will tend to experience comparable amounts of displacement in the positive and negative directions, but may accrue a significant net displacement along the direction. The operator-theoretic model successfully captures this behavior for lead times up to (i.e., revolutions of the vortex in ), but since it operates at a finite resolution (determined by the spectral truncation parameters , , and ) and also has diffusion, it eventually develops biases as it fails to capture the increasingly small-scale variations developed by observables such as . For such observables, the effect of diffusion in the semigroup eventually dominates; the signature of this effect in Movie 2 is spurious high-velocity tracer motions leading to the formation of voids in where tracers should be present. Clearly, the temporal extent of useful forecasts for an observable depends on both on how strongly it projects on as opposed to , as well as the spectral truncation and diffusion regularization () used.
Next, we consider prediction of probability densities on tracers via the Perron-Frobenius operator as described in Section 2.4.3. Here, we perform prediction experiments initialized with a density given by a product of circular Gaussians,
[TABLE]
where and . This density represents uncertainty in both the velocity field state in and the tracer positions in . Note that attains its maximum value at the point , and for increasing it becomes increasingly concentrated around that point. In what follows, we visualize the temporal evolution of through the marginal density on ; the quantity is equal to the probability density to find a tracer at the point at lead time marginalized over all flow states given the initial density . We also examine the one-dimensional (1D) marginal densities and . Clearly, the evolution of , , and is neither Markovian nor unitary; that is, the norms of , , and with respect to the appropriate normalized Haar measures are generally not constant (though the corresponding norms are, of course, always equal to 1).
Figure 4 and Movie 3 show the evolution of , and for the initial probability density from (23) with and . With these parameters, at is concentrated upstream of the vortex center (for the state maximizing the initial density ) along the direction, and has a positive offset along the direction (see top-left panel in Fig. 4). As a result, at the initially isotropic is swept and sheared by the moving vortex, producing an anisotropic density. This process is repeated with each revolution of the vortex center around the periodic domain, leading to a filamentary density with highly non-Gaussian features. As noted earlier, this flow produces significantly more mixing along the direction compared to the direction, and as a result, even after periods () most of the probability mass remains concentrated in the portion of the domain.
Overall, the results in Fig. 4 and Movie 3 illustrate that, at least as far as the marginal densities are concerned, the data-driven operator-theoretic approach in Section 2.4.3 agrees well with the results from the Monte Carlo simulation based on the full model for Lagrangian tracer advection. In fact, the range of accurate forecasts from the data-driven model is longer than in the case of the observables and characterizing the tracer positions (Fig. 3 and Movie 2). This is likely due to the fact that the marginal density is given by projection of the full density to the subspace of spanned by functions that do not depend on state (i.e., ; see Section 2.1), and moreover and involve additional projections onto subspaces of spanned by - and -independent functions, respectively. Such projections are expected to cancel at least some of the errors that may be present in the full density , contributing to an increase of forecast skill for the marginal densities. In other words, these examples are a demonstration of the (perhaps, obvious) fact that for a fixed dynamical system and spectral truncation the prediction skill is observable dependent.
3.2 Switching Gaussian vortices
Our second example uses the same torus domain as the moving-vortex example in Section 3.1, but in this case we consider the streamfunction
[TABLE]
where and are positive vortex strength and concentration parameters, respectively. Coupled with the periodic dynamics on , the streamfunction in (24) describes a pair of vortices centered at and in (as usual, and represent canonical angle coordinates), and modulated in time by 90∘ out of phase sinusoids. Representative snapshots over half a period of this flow are shown in Fig. 5. There, it can be seen that an anticklockwise vortex centered at is present at , and this vortex is replaced by a vortex of the same sense centered at at , which is in turn replaced by an anticklockwise vortex centered at at . In the intervening states and , the flow has the structure of shear layers parallel to the and coordinate, respectively. This process is repeated over the interval with a change of vortex sign.
Loosely speaking, this “switching vortex” flow can be thought of as a continuous analog of a piecewise-constant stirring flow by point vortices [58] (also known as blinking vortex flow), which is known to exhibit chaotic advection. At the very least, one would expect the streamfunction in (24) to produce more complex spectral behavior than the moving-vortex example in Section 3.1, for the former is not decomposable via submersions associated with Galilean transformations as the latter is (see C). In particular, as is evident from Fig. 5, the streamlines associated with (24) undergo changes in topology as varies, meaning that it is not possible to construct a submersion analogous to in C mapping the state-dependent velocity field on to a steady velocity field with closed streamlines (the latter would imply the existence of non-trivial eigenfunctions of the generator at eigenvalue zero).
Qualitatively, the larger the parameter in (24) is, the more overturning we expect the tracers to undergo during one period of the dynamics on , leading to stronger mixing. In fact, while we do not have rigorous results to justify this assertion, it appears plausible that for sufficiently large the only eigenfunctions of for the class of streamfunctions in (24) are constant on . If this is indeed the case, then all numerical eigenfunctions we compute are “genuinely approximate” Koopman eigenfunctions, in the sense that they converge to eigenfunctions of the regularized generator as the spectral order parameter tends to infinity (such eigenfunctions always exist by Proposition 6(v)), but the eigenfunctions of do not converge to eigenfunctions of in the limit of vanishing diffusion regularization parameter . Thus, in such situations, diffusion regularization is expected to play an essential role for the recovery of coherent patterns via the approach of Section 2.3.3.
Following a similar approach as in Section 3.1, we compute the matrix elements of the generator in the Fourier basis of using the integral identity for circular Gaussian integrals in (17). Specifically, we have
[TABLE]
and the matrix elements are given by (20) as in the moving-vortex flow. Using these results, we apply the operator-theoretic schemes of Section 2 using the same spectral truncation parameters as in Section 3.1. Throughout, we work with the flow frequency parameter .
3.2.1 Koopman eigenvalues and eigenfunctions
Our objectives in this section are to study the properties of coherent spatiotemporal patterns recovered by numerical Koopman eigenfunctions for different values of the vortex strength parameter , as well as the dependence of these patterns to the diffusion regularization parameter . Figure 6 shows snapshots of representative eigenfunctions computed for the flow with the streamfunction parameters and , and the diffusion regularization parameters and . Videos showing these eigenfunctions are provided in Movies 4 and 5, respectively; Table 2 lists the corresponding eigenvalues and Dirichlet energies. As in Section 3.1, we computed these results using Matlab’s eigs solver, requesting 51 eigenvalues of minimal modulus.
First, consider the results for . The numerical eigenfunctions in this case includes a family, which includes eigenfunctions , , , , , , , , and shown in Fig. 6 and Movie 5 which are characterized by four globular patterns undergoing a sloshing motion as a result of advection by the time-dependent velocity field. According to the results in Table 2, the eigenvalues corresponding to these patterns have essentially vanishing imaginary part (that is, twelve orders of magnitude smaller than ), indicating that, at least to the resolution afforded by our spectral truncation, these patterns are conserved on Lagrangian tracers; that is, the level sets of these eigenfunctions induce an ergodic partition on the state space . As expected, the eigenfunctions in this set with higher Dirichlet energy exhibit increasingly smaller-scale oscillations, which are arranged in a concentric manner relative to the center of each cluster. In other words, the eigenfunctions in this family provide an increasingly fine partition of into quasi-invariant sets. Hereafter, we refer to eigenfunctions of this class as class 1 eigenfunctions.
In addition to the patterns described above, the numerical Koopman eigenfunctions include a family (hereafter, class 2), members of which are , , and displayed in Fig. 6 and Movie 5, which are concentrated in the separatrix regions between the four globular clusters of the previous family. The eigenvalues corresponding to these eigenfunctions have nonzero imaginary parts (e.g., for , , and ), indicating that these eigenfunctions vary periodically in a Lagrangian frame following the tracers. The concentration of these eigenfunctions to subsets of of small Lebesgue measure and the sharp spatial gradients that they exhibit (see, e.g., in Fig. 6 along the diagonal near ) are reminiscent of the behavior of numerical eigenfunctions and in Fig. 2 associated with the continuous spectrum of the moving vortex flow.
Next, we examine the dependence of the patterns described above on the diffusion regularization parameter . Inspecting the results in Fig. 6 and Movie 4, it is evident that the leading class 1 eigenfunctions identified for are relatively robust under an increase of from to ; for example, eigenfunctions have clear counterparts, . More quantitatively, the Dirichlet energies of those eigenfunctions (see Table 2) do not change by more than (and in some cases this change is as little as . Together, these results suggest that class 1 eigenfunctions behave smoothly as and, correspondingly, that they approximate true Koopman eigenfunctions. In contrast, class 2 eigenfunctions such as , , and at have a significantly more sensitive dependence on . That is, while analogs of some of these eigenfunctions can be identified in the results (e.g., and at appear to be related to and at , respectively), the spatial patterns of these eigenfunctions are visibly affected by diffusion (compare, e.g., at with at ). This also reflected in the corresponding Dirichlet energies which change by approximately a factor of five between the two cases.
We now increase the vortex strength parameter to and examine the resulting numerical eigenfunctions, using again the diffusion regularization parameter values and . The resulting eigenfunctions and eigenvalues are shown in Fig. 7, Movie 6 (), Movie 7 (), and Table 3. There, it is evident that in this flow with stronger stirring and overturning, the geometrically smooth patterns recovered in the experiments are replaced by significantly more complex filamentary patterns. Nevertheless, the eigenfunctions still can be grouped into class 1 and class 2 families, the former characterized by approximately zero imaginary part of the corresponding eigenvalue and the latter consisting of eigenfunctions concentrated on subsets of the spatial domain of small Lebesgue measure while having nonzero . Examples of class 1 and class 2 eigenfunctions at shown in Fig. 7 and Movie 7 are and , respectively.
Despite these similarities, the and results have a fundamental difference in that the former include a class of eigenfunctions (class 1) that are only weakly affected by the examined changes in , whereas in the latter case these changes in impart a significant change to all of the eigenfunctions. For instance, when is decreased from to , the Dirichlet energies of class 1 eigenfunctions and increases by a factor of 8 and 6, respectively (see Table 3). This implies that the eigenfunctions acquire increasingly small-scale features with decreasing , as is also evident in Fig. 7 and Movies 6 and 7. In that regard, these eigenfunctions behave similarly to strange eigenmodes for diffusive Lagrangian tracers in time-periodic flows [19], obtained via Floquet analysis of an advection-diffusion operator in [18].
3.2.2 Prediction of observables and probability densities
In this Section, we discuss prediction of observables and probability densities for the switching-vortex flow. Our experiments are analogous to those for the moving-vortex flow in Section 3.1.2; that is, we consider prediction of the observables and in (22) characterizing the position of Lagrangian tracers, and also consider probability densities with the Gaussian initial conditions in (23). Hereafter, we restrict attention to the flow with frequency and vortex concentration and strength parameters and , respectively. Moreover, we work throughout with the diffusion regularization parameter .
Figure 8 and Movie 8 show the evolution of the positions for an ensemble of tracers (initially arranged on a uniform square grid in ) obtained via the operator-theoretic model from Section 2.4.1, compared against the evolution obtained by explicit integration of the ODEs governing the advection of tracers in the switching-vortex flow. We used the same methods and error tolerance parameters to compute these results as in Section 3.1.2; here, the forecast timestep is . In Fig. 8 and Movie 8 it can be seen that, at least over short to moderate lead times (), the operator-theoretic model is able to predict the evolution of the tracer positions with comparable accuracy to the explicit ODE model, but eventually it develops biases due to the combined effects of diffusion and spectral truncation. Thus, the operator-theoretic model performs comparably to the moving-vortex experiment in Section 3.1.2, but note that in this case the dynamics is more complex since mixing takes place with respect to both the and coordinates of the tracers.
Figure 9 and Movie 9 show the evolution of the marginal densities , , and (see Section 3.1.2) over the forecast interval obtained via the operator-theoretic model of Section 2.4.3 and a Monte Carlo ensemble of 300,000 particles. The initial density on is given by (23) with location and concentration parameters and ; as shown in Fig. 9, the corresponding initial marginal density is concentrated to the right of the vortex center associated with the state (i.e., the mode of the density ). At time , the initially radially symmetric density is sheared by the flow, and probability mass is gradually expelled from the vicinity of where is concentrated. During this process, and the 1D densities and become highly non-Gaussian. As is evident from Fig. 9 and Movie 9, the operator-theoretic model agrees well with the evolution of the marginal densities over the full forecast interval.
4 Approximation in a data-driven basis
Thus far, the development of our methods for coherent pattern extraction and nonparametric prediction, as well as the associated numerical experiments, were made under the strong assumption that orthonormal bases for the Hilbert spaces and associated with the state space and the physical domain , respectively, are available. In this Section, we partially relax this assumption using kernel algorithms [28, 30] to formulate the techniques of Section 2 in a basis of functions on learned from the observed velocity field snapshots . In particular, our objectives are that (1) this basis approximates (in a sense that will be made precise below) the Laplace-Beltrami eigenfunction basis of associated with the Riemannian metric from Section 2.3.2, and (2) the Dirichlet energies of the basis functions can also be approximated.
Our approach follows closely [22, 59], where an analogous data-driven basis was used for mode decomposition and prediction in ergodic dynamical systems. Here, we extend this framework to the setting of skew-product systems under the additional assumption that observations of the velocity field driving the flow in are available. Moreover, we introduce a more general framework for computing Dirichlet energies for the data-driven basis functions. It should be noted that in this work we do not address the question of constructing a data-driven basis of , and we will continue to assume that the basis is available. Nevertheless, building data-driven bases of from an arbitrary collection of points in should be possible through analogous kernel techniques to those used to build the basis for , or via the variety of meshing techniques developed in finite element methods. It is also worthwhile noting that despite being formulated for finitely sampled data, the schemes for coherent pattern extraction and prediction introduced below are structurally very similar to those already presented in Section 2. Thus, it should be possible for readers mainly interested in numerical results to skip ahead to Section 5 where we discuss experiments with aperiodic tracer flows driven by L96 systems.
4.1 Kernel function
The starting point for the construction of our data-driven basis is a specification of a kernel function providing a measure of similarity of an arbitrary pair of states in which is computable using observed velocity field data in . Mathematically, we require that
is symmetric, i.e., for any ; 2. 2.
is smooth (hence, bounded above by compactness) on ; 3. 3.
is bounded below, i.e., there exists a constant such that .
In what follows, we work with the variable-bandwidth kernel introduced in [30], viz.
[TABLE]
where is a positive bandwidth parameter, and a positive function in , bounded away from zero. In particular, the bandwidth function is an estimate of the function , where is the density of the invariant measure relative to the volume of the ambient-space metric introduced in Section 2.3.2, and the dimension of . This bandwidth function is constructed such that converges to almost surely with respect to the starting state as and uniformly with respect to as . Details on this construction can be found in B.2. With the definition in (26), the conditions on kernels stated above are satisfied by construction and the fact that is compact. Moreover, by the pointwise ergodic theorem, as , converges almost surely to the value of the kernel
[TABLE]
which also meets conditions 1–3 stated above.
Remark 11**.**
The role of the bandwidth functions and in (26) and (27), respectively, is to implement the conformal transformation leading to the Riemannian metric (which has uniform volume form relative to the invariant measure; see Section 2.3.2) from the ambient space metric .
Remark 12**.**
The analysis that follows relies on the assumption made in Section 2.3.2 that is a smooth embedding of into the space of velocity fields . As stated in remark 2, if that condition is not met then it is possible to replace in (26) (as well the definition of ) by the delay-coordinate map given by
[TABLE]
where is a positive integer parameter. According to the theory of delay-coordinate maps [36], will generically be an embedding for sufficiently large and under mild assumptions on and . Moreover, as discussed in [23], performing delays can be beneficial even if is already an embedding as it improves the noise robustness of the data-driven basis functions. Our scheme is also applicable when takes values in other spaces than (e.g., if the fluid dynamical driving system has additional degrees of freedom such as temperature), so long as and/or or are embeddings.
4.2 Kernel integral operators
Having specified an appropriate kernel, we now apply the normalization procedure introduced in the diffusion maps algorithm [28] and further developed in [31] to construct a kernel integral operator which can be used to approximate the eigenfunctions of the Laplace-Beltrami operator associated with the Riemannian metric and the associated Dirichlet energies. In what follows, we describe the construction of this integral operator and discuss spectral convergence results justifying its utility in providing the data-driven basis used in our mode decomposition and prediction schemes. The discussion below is mainly based on results established in a number of works in the literature, including [28, 29, 30, 60, 31]. We believe that presenting this material here in a fairly self-contained manner is useful given that each of these works addresses a fairly distinct aspect of the problem of interest here, but readers familiar with these topics may wish to skip to Section 4.3.
We begin by establishing the relevant function spaces for our data-driven scheme. In particular, note that in practical applications with finitely many samples we cannot work with the Hilbert space directly. Instead, we consider the Hilbert space of complex-valued functions associated with the sampling measure . That is, consists of equivalence classes of functions which have common values and are square-summable on the sampled states . Its associated inner product and norm are
[TABLE]
Of course, is isomorphic to , but we prefer to work with the former space to emphasize that its elements are equivalence classes of functions defined on the same underlying space as in the case of . Nevertheless, given , it is convenient to represent the corresponding equivalence class in by the column vector storing the values at the sampled states ; we will use this notation whenever we wish to distinguish between functions on and equivalence classes of functions in . Similarly, we will use the notation to represent the equivalence class associated with . Clearly, we have , where is the standard dot product for vectors in . Moreover, by the pointwise ergodic theorem, given such that and are both in , then for -a.e. starting state , and are also in and as , converges to .
Besides these Hilbert spaces, we will make use of the Banach space of continuous functions on equipped with the uniform norm (denoted ). This space will provide a common space to define kernel integral operators associated with the finite dataset and the whole state space , allowing us to take advantage of existing spectral convergence results in the machine learning and statistics literature [29]. It is then straightforward to lift those results to the Hilbert space setting of and .
With these definitions in place, consider the bounded integral operators and defined as
[TABLE]
Clearly, we can also define analogs and of by and , where is the canonical restriction map, . Although related by trivial restrictions, the distinction between these operators will be useful below. Similarly, we introduce and , where is the canonical inclusion map, . By the assumptions on and stated in Section 4.1, the range of both and is included in . Note that can be evaluated with no approximation at arbitrary given only the values at the sampled states . Moreover, by the pointwise ergodic theorem, converges to in uniform norm for -a.e. starting state . In the case of and , the symmetry of and implies that these operators are self-adjoint. In addition, we have:
Lemma 13**.**
The operators and are compact.
A proof of Lemma 13 can be found in A.3.
Next, using and , respectively, we construct normalized operators and by applying the normalization procedure in [28, 60]. Specifically, we set
[TABLE]
where and are defined as
[TABLE]
and is the function equal to 1 at every point in . Both and can be alternatively be expressed as kernel integral operators, viz.
[TABLE]
where the kernels and satisfy
[TABLE]
for all . By our conditions on kernels and compactness of , the functions , , , and introduced above are all positive functions bounded away from zero, and , for every by the pointwise ergodic theorem. As a result, and are smooth with respect to both of their arguments, positive, and bounded away from zero. In [60], the normalization steps in (28) and (29) are referred to as left and right normalization, respectively.
Remark 14**.**
Intuitively, the quantity can be thought of as an estimate (up to a proportionality constant) of the density of the invariant measure relative to the volume form of a Riemannian metric on that depends on . For the choice of kernel in (26), that metric is in fact the desired metric , whose volume form is equal to the invariant measure, and one can show that for every , where is a positive constant. As a result, we have . Note that despite that and are constants in the limit, their behavior at is also important (e.g., for Lemma 15(ii) ahead), and at that order they are non-constant on .
Next, through an appropriate use of restriction/inclusion operators, we define operators and having the same kernel as , and similarly we define and having the same kernel as . Note that by smoothness of and , both and are Hilbert-Schmidt (hence compact) kernel integral operators. In particular, (31) and the fact that and are bounded away from zero imply imply that and generate ergodic Markov semigroups, respectively and , and thus each has a simple eigenvalue at 1 corresponding to a constant eigenfunction. The following Lemma lists four key properties of the averaging operators introduced thus far.
Lemma 15**.**
(i) For -a.e. starting point in the training data, converges as to in norm.
(ii) For every and , the mapping is ; in particular, for every we have
[TABLE]
where the estimate holds uniformly on .
(iii) For every , converges as to the heat operator , , in norm.
(iv) and are similar to the self-adjoint Hilbert-Schmidt operators and , respectively, where
[TABLE]
Moreover, we have as .
A proof of this Lemma can be found in A.4.
Let now and be orthonormal bases of and , respectively, consisting of eigenfunctions of and . Let also and , respectively, be the corresponding (real) eigenvalues. Since is similar to (resp. is similar to ) and (resp. ) is an ergodic Markov operator, it follows that the eigenvalues admit the ordering (resp. ). Moreover, one can verify that eigenfunctions of and corresponding to these eigenvalues are given by
[TABLE]
respectively, where and are the densities of the stationary measures of the Markov semigroups generated by (resp. ) relative to (resp. ). Explicitly, we have , , and , where is a positive smooth function, bounded away from zero. Analogous formulas hold for with . Note that while and are bases of and , respectively, they are generally not orthonormal (although, as we will see in Corollary 17 ahead, the lack of orthonormality vanishes in a suitable limit and ). For later convenience, we introduce the kernel integral operators and which are defined using the same kernels as in the definitions of and , respectively.
By smoothness of and compactness of , and are Hilbert-Schmidt integral operators, and therefore compact. Moreover, by Lemma 13, , and are also compact. (In fact, and have even finite rank.) As a result, the eigenvalues of , , , and can only accumulate at zero, and their nonzero eigenvalues have finite multiplicities. For such eigenvalues, convergence in operator norm implies convergence of eigenvalues and eigenfunctions corresponding to simple eigenvalues. For non-simple eigenvalues, convergence in operator norm implies convergence of projection operators to the corresponding eigenspaces [29]. In particular, we have:
Lemma 16**.**
(i) Let be an eigenfunction of (resp. ) at eigenvalue . Then, (resp. ) is an eigenfunction of (resp. ), also at eigenvalue .
(ii) Let be an eigenfunction of at nonzero eigenvalue . Then is an eigenfunction of , also at eigenvalue . The analogous result holds for eigenfunctions of and .
(iii) As , the eigenvalue of converges -a.s. to the eigenvalue of , assuming that the latter is nonzero. Moreover, for each eigenfunction of corresponding to , there exist eigenfunctions of such that , where the limit is taken in uniform norm.
(iv) The eigenvalues of satisfy
[TABLE]
where is the -th eigenvalue of . Moreover, for each eigenfunction of at eigenvalue there exist eigenfunctions of such that in norm.
(v) The analogous results to (i)–(iv) hold for the self-adjoint operators and .
A proof of Lemma 16 is included in A.5. Together, the results in Lemma 16 imply that we can approximate the eigenvalues and eigenfunctions of using the eigenvalues and eigenfunctions of either or (both of which are accessible from data), in the following sense.
Corollary 17**.**
For -a.e. starting state in the training data,
[TABLE]
Moreover, there exist eigenfunctions of (resp. eigenfunctions of ) that approximate the eigenfunctions of , in the sense that
[TABLE]
in uniform norm, where is the unique smooth representative of . In particular, for any fixed state , the values and converge -a.s. to .
Details on the numerical solution of the eigenvalue problems for and are provided in B.3.
4.3 Data-driven bases and the associated Dirichlet energies
Having established the spectral convergence results of Section 4.2 (in particular, Corollary 17), we are ready to specify the data-driven bases and their associated Dirichlet energies to be used in the mode decomposition and forecasting techniques of Sections 2.3 and 2.4, respectively. In essence, given a fixed dataset of samples, this procedure entails three steps, namely (1) selection of an appropriate value of the kernel bandwidth parameter ; (2) specification of basis functions for the finite-dimensional space ; (3) an assignment of Dirichlet energies for the basis functions and construction of an associated Sobolev space (needed for the variational formulation of the Koopman eigenvalue problem in Section 2.3.3). We now describe these steps.
4.3.1 Kernel bandwidth selection
We select using the method developed in [61, 30], which was also employed in [62, 22, 23, 24]. The method tunes by examining the behavior of the function , which is approximated in the case of finite data by . As and , that function tends to the constant values 0 and 1, respectively. Moreover, as shown in [61, 30], there exists a regime of intermediate values where exhibits the scaling , where is a constant independent of the intrinsic dimension of . In particular, in that regime, is approximately equal to . In [30], the quantity is interpreted as a “resolution” of the kernel, and they suggest choosing as the maximizer of that function, yielding also an estimate of . In what follows, we select using that criterion, approximating through finite differences of on a logarithmic grid. Further details on this procedure can be found in [62, 23].
4.3.2 Basis functions
In light of the results of Section 4.2, obvious candidates for data-driven basis functions of are the eigenfunctions of and , i.e., and , respectively. According to Corollary 17, both of these choices lead to consistent approximations of the Laplace-Beltrami eigenfunctions in the limit of infinitely many samples and vanishing bandwidth parameter, but away from that limit they each have their advantages and disadvantages.
More specifically, the main advantage of the basis is that it is orthonormal on . Its main disadvantage is that it is not an eigenbasis of a Markov operator, and as a result, estimating the associated Dirichlet energies using the eigenvalues from diffusion maps would lead to certain biases (e.g., corresponding to eigenvalue would be assigned zero Dirichlet energy while being non-constant; the latter follows from (34) and the fact that is constant whereas is non-constant). Working with the basis avoids this issue, but only at the expense of losing orthogonality of the basis functions.
One way of alleviating the latter issue is to observe that the are orthonormal on the Hilbert space associated with the stationary distribution of the associated Markov semigroup (see Section 4.2). That is, the are orthonormal with respect to the inner product
[TABLE]
where is the stationary density of introduced in Section 4.2. This suggests that we can formulate our data-driven techniques in this space instead of . However, in doing so one must remember that taking the limit at fixed will not lead to convergence to the inner product of the Hilbert space associated with the invariant measure , but rather to the inner product of a different Hilbert space, , associated with the stationary density of (in particular, the Koopman group is non-unitary in this Hilbert space). However, the limit can be taken in conjunction with a decreasing function (e.g., via the tuning procedure in Section 4.3.1), so that converges to . In particular, in this limit, the functions and both converge to (see also Remark 14).
Due to the important role that Dirichlet energies play in our schemes, in the numerical experiments presented in Section 5 ahead, we always work with the basis and the weighted inner product space . However, in practice we find that the range of values of is usually not too large (e.g., in the applications of Section 5 it does not vary by more than ), and therefore comparable results can also be obtained using the orthonormal basis of .
For later convenience, we note that since is finite-dimensional, the pairwise products of the basis elements are also in . Thus, we can write , where the coefficients
[TABLE]
can be thought of as “structure constants” for the algebra of equivalence classes of functions . Since the are all real, the are invariant under arbitrary permutations of . Moreover, because is a constant, we have . Analogous coefficients can also be defined for the basis of by continuity of the , but note that is not closed under multiplication.
4.3.3 Dirichlet energies
We now discuss how to construct data-driven analogs of the Dirichlet form and the associated Laplace-Beltrami operator needed to implement the variatonal formulation of the eigenvalue problem for coherent patterns in Section 2.3.4. First, before stating definitions for these quantities suppose that we have at our disposal quantities and approximating the Laplace-Beltrami eigenvalues . In particular, we assume that the and are ordered in a non-decreasing sequence, that , and that all other and are positive and finite. We also assume that for fixed and by ergodicity. Then, for , we define the Sobolev spaces
[TABLE]
Clearly, since these spaces are finite-dimensional (though we prefer to maintain their distinction for consistency with the infinite-dimensional case), but are strict subspaces of . However, we have:
Lemma 18**.**
The spaces are dense in .
Proof.
Let be an arbitrary function in , and let be any positive number. Then, there exists an integer such that for all the function satisfies , and moreover is finite so . ∎
Let now and be arbitrary functions in and , respectively. As usual, we equip these spaces with the inner products
[TABLE]
and the norms and . In addition, for , we define the sesquilinear forms and ,
[TABLE]
and for we define the operators and with
[TABLE]
These forms and operators have the following properties:
Proposition 19**.**
(i) and obey the bounds
[TABLE]
respectively.
(ii) and are self-adjoint. Moreover, and are dissipative; i.e.,
[TABLE]
A proof of this Proposition can be found in A.6. In Section 4.4 ahead, and will be used to construct data-driven analogs of the Dirichlet form employed in the Koopman eigenvalue problem for coherent patterns in Section 2.3.4. It is worthwhile noting that Proposition 19(ii) in conjunction with [42] implies that and (or their maximal self-adjoint extensions) are generators of strongly-continuous contraction semigroups, and , respectively, with and . For appropriate choices of and (e.g., (40) below), these semigroups approximate the heat semigroup , , associated with the Laplace-Beltrami operator .
We now return to the question of defining and . As with the definition of the data-driven basis functions in Section 4.3.2, here too there are options. Below, we list three possibilities for ,
[TABLE]
which we refer to as option 1, 2, and 3 respectively. By Lemma 16(iii), as , these quantities converge to
[TABLE]
respectively.
Remark 20**.**
In our definitions of and via options 2 and 3, we have tacitly assumed that the and are all positive, or, equivalently, that and are positive operators. In the event that these operators have negative eigenvalues, one can modify the definitions of and for (with the index of the smallest positive eigenvalue) to a positive, increasing sequence with no accumulation points. For example, one could set , which is consistent with Weyl’s law for the asymptotic growth of Laplace-Beltrami eigenvalues on compact Riemannian manifolds. It should be noted that this modification is mainly formal since for any finite spectral truncation parameter used in our schemes there exist and such that, for all and , and are positive. For the rest of the paper, we will treat the and as positive eigenvalues.
With option 1, the are equal to eigenvalues of the operator , i.e., the normalized graph Laplacian widely used in machine learning [27, 28, 29, 60]. It is a well-known fact that for an appropriate decreasing sequence , the from this option converge to ; that is, are consistent estimators of the Laplace-Beltrami eigenvalues. At fixed , the limit quantities , and hence the corresponding operator , are bounded (in particular, since is Markov). The from option 2 also converge to the true Laplace-Beltrami eigenvalues for a decreasing sequence in accordance with Corollary 17, but notice that since , in this case and become unbounded as . Finally, the from option 3 converge to (again, for an appropriate decreasing sequence ) under the additional assumption that the remainder in the Taylor expansion is asymptotically smaller than as . At finite and nonzero we have , and therefore the Dirichlet energies from this option dominate those from option 2. Moreover the ratio increases with .
Note now that each of the definitions of in (40) will lead to a different regularized Koopman eigenvalue problem for coherent spatiotemporal patterns (after employing these quantities to define Dirichlet forms for functions on , as described in Section 4.4 below). At a minimum, whatever the specific choice of (and hence ) is, the associated Sobolev space should have the property that the sesquilinear form associated with the generator of the Koopman group on is bounded on (see Section 2.3.4), at least in the limit . While all of the definitions for in (40) meet that asymptotic requirement, they will of course behave differently at finite and nonzero . In particular, in light of the above arguments, we expect the from option 3 to provide the strongest and most selective regularization among the three options, followed by option 2, and then option 1. Indeed, in experiments we found that the Dirichlet forms from option 3 generally led to higher-quality numerical coherent patterns, so in what follows we nominally work with that option.
4.4 Data-driven techniques for skew-product systems
We now construct analogs of the schemes of coherent pattern extraction and nonparametric prediction of Section 2 using the basis functions and Dirichlet forms from Section 4.3. Throughout this Section, we use the symbols and in place of and , respectively, in order to avoid notational ambiguities when building tensor product bases on the product space .
4.4.1 Approximate Koopman generator on the Hilbert space associated with the sampling measure
As stated in Section 4.2, in practical applications involving finite numbers of samples we do not have access to the Hilbert space where the Koopman group is naturally defined, and we work instead with Hilbert spaces such as and associated with the sampling measure . While the latter two spaces form natural settings for defining kernel integral operators and their associated eigenfunctions and Dirichlet forms, they are unfortunately not suitable for defining a Koopman operators analogous to via composition with the flow map. In essence, this is because the dynamical flow does not preserve null sets with respect to ; that is, for some measurable set does not imply that vanishes too. More specifically, Koopman operators on and would have to act on equivalence classes of functions that have common values on the sampled states , but given two representatives and in the functions and could lie in different equivalence classes, resulting in “” being ill-defined.
Despite this issue, it is nevertheless possible to define operators on and that approximate and the generator without explicit use of the flow map. Here, we focus on approximations of the generator that appears in the schemes of Section 2. As a concrete example of such an approximation, consider the operator (here, is the sampling interval; see Section 2.1), whose action on is given by with
[TABLE]
Intuitively, , can be thought of as a (second-order, central) finite-difference approximation of for . Such finite-difference approximations were previously used in [59] to construct -dependent kernels and in [22, 23, 24] for approximating Koopman generators as done here.
The following Lemma (proved in A.7) establishes a consistency property of the sesquilinear form associated with the approximate generator, evaluated on the data-driven basis elements:
Lemma 21**.**
Let be any two eigenfunctions of from Section 4.2 converging, in the sense of Corollary 17, to eigenfunctions of . Then,
[TABLE]
Moreover, the order of the limits and can be interchanged.
As we discuss below, the consistency results in Corollary 17 and Lemma 21 are sufficient to establish consistency of the techniques presented in Section 2 in the presence of errors due to finite training datasets and discrete-time sampling.
Remark 22**.**
Of course, besides the second-order, central finite difference scheme in (42), other approaches could be employed to construct the approximate generator , including higher-order finite differences or interpolation techniques. For instance, a straightforward improvement of (42) would be to replace and with values obtained from linear extrapolation of and , respectively, though this modification would have only minor consequences at large .
4.4.2 Koopman eigenvalue problem for coherent spatiotemporal patterns
Recall that we assume that the Laplace-Beltrami eigenfunction basis for the Hilbert space of functions on the spatial domain and the associated eigenvalues are given. Thus, in order to formulate an analog of the Galerkin method in Section 2.3.4 utilizing the data-driven constructions of Sections 4.1–4.4.1, it suffices to consider the Hilbert space associated with the product measure . Our orthonormal basis for this space has a tensor product form, consisting of the functions , where and . We also define , the associated Sobolev spaces
[TABLE]
with , , the Dirichlet form with , and the Laplace-Beltrami operator with . As in Section 2.3.2, we order the basis elements in order of increasing Dirichlet energy , and when convenient we use single-index notation, and , , consistent with that ordering. We also introduce the normalized eigenfunctions and , which form an orthogonal basis of with and . Note that as with the corresponding basis of in Section 2.3.2, we have if and .
It follows by construction of the basis that a consistency result analogous to Corollary 17 holds; that is, as and , for every eigenfunction of the Laplace-Beltrami operator there exist converging to in the sense of Corollary 17. In the case of the , if these quantities are computed using options 1 and 2 in (40), then Corollary 17 implies that they converge to the Laplace-Beltrami eigenvalues . If they are computed using option 3 (i.e., our nominal choice), then the may either (1) converge to , (2) converge to finite quantities , (3) diverge. In addition to the first of these cases, case 2 is also sufficient for the well-posedness of the Koopman eigenvalue problem introduced below. In particular, associated with is a Sobolev space on (defined in the obvious way), which is dense in and satisfies . This space can be used as an alternative to in the variational eigenvalue problem in Definition 8. In what follows, we will use the symbols and to represent the normalized eigenfunctions and Dirichlet form, respectively, for . Then, for any and converging to and , respectively,
[TABLE]
Next, using the results of Section 4.4.1, we introduce an operator on approximating the generator of the skew-product dynamical system on . Specifically, we define , , where is the natural lift of on (that is, where and ), and is the lift of the state-dependent vector field on as in Section 2.1. This operator has the analogous consistency property to that established for in Lemma 21, i.e.,
[TABLE]
Moreover, for all Dirichlet energy options in (40) we have
[TABLE]
In the case of options 1 and 2 we additionally have , and thus we obtain a consistent approximation of , where is the sesquilinear form on associated with (see (6)). In the case of option 3, that approximation may not be consistent, but so long as the are finite, we nevertheless have
[TABLE]
Thus, can also be defined as a sesquilinear form on .
Remark 23**.**
One of the properties that must necessarily hold in order for the variational eigenvalue problems employed in this work to be well posed is that is bounded on the domain of the Dirichlet form used for regularization (e.g., and with domains and , respectively). Depending on how the are defined, it may be possible to satisfy this condition without taking the limit . For instance, the from options 2 and 3 in (41) are unbounded for all , so it is possible that is bounded on ; the order-1 Sobolev space associated with the . If it indeed holds, this property would suggest that our schemes are applicable for more general state spaces than smooth manifolds (since limits may be problematic on such spaces), including non-smooth compact topological attractors and/or compact supports of physical measures. Addressing these issues (which arise more broadly than the skew-product systems of interest here) is beyond the scope of this work, but in Section 5 we provide numerical evidence that our approach may indeed be applicable if is not smooth.
Fixing now the diffusion regularization parameter and the spectral truncation parameters , , and , we introduce the -dimensional Galerkin approximation spaces and consider the following variational eigenvalue problem, which is analogous to the eigenvalue problem in Definition 8:
Definition 24** (data-driven eigenvalue problem for coherent patterns).**
Find and such that for all the equality holds, where and are the sesquilinear forms on , given by
[TABLE]
Solving the eigenvalue problem in Definition 24 is equivalent to solving an matrix generalized eigenvalue problem with an entirely analogous structure to (8). Moreover, it follows from (44) and (45) that for any , the Galerkin method in Definition 24 converges as and to a Galerkin method for the variational eigenvalue problem on given by , , restricted to the -dimensional approximation space. By the properties of (in particular, (46)) and , this problem is well-posed and the Galerkin approximation converges as . If we have in addition (which is necessarily the case for options 1 and 2 in (40) but not for option 3), then it follows that , and the Galerkin method in Definition 24 converges to the one in Definition 8. It should be noted that and essentially differ only with respect to diffusion regularization, and are both expected to yield adequate approximate Koopman eigenfunctions at small values of the regularization parameter (at least when has discrete or mixed spectrum). As stated in Section 4.3.3, in practical applications we find that the diffusion regularization from option 3 leads to a modest improvement in the quality of numerical Koopman eigenfunctions. Explicit formulas for the matrices appearing in the generalized eigenvalue problem from Definition 24 in the case of a doubly-periodic spatial domain, , can be found in A.8.
4.4.3 Prediction of observables and probability densities
In this section, we introduce prediction schemes for observables and probability densities analogous to those in Section 2.4, but formulated in the data-driven basis of . In what follows, we consider a prediction observable such that lies in ; this ensures that for -a.e. starting state in the training data, lies in . It then follows from the convergence of the data-driven eigenfunctions to the continuous functions established in Corollary 17, that the expansion coefficients converge to for all ; specifically,
[TABLE]
Let now be the data-driven regularized generator defined via , where , and and are the data-driven generator and Laplace-Beltrami operators from Section 4.4.2, respectively. Let also be orthogonal projectors from to , and define (cf. (12))
[TABLE]
Then, assuming that the eigenvalues of converge to the eigenvalues of , we have
[TABLE]
Together, (47) and (49) lead to the following consistency result:
Lemma 25**.**
If the eigenvalues converge to , then for any observable such that and -a.e. starting state in the training data,
[TABLE]
Now, as with the eigenvalue problem in Section 4.4.2, if the do not converge to as and , but instead converge to finite quantities , then these quantities are bounded below by , and the matrix elements converge to the matrix elements of , where is a symmetric, dissipative operator such that . An analogous consistency result to that in Lemma 25 then holds for the contraction semigroup generated by , and this semigroup can also be used as an approximation of the Koopman group .
Following the approach of Section 2.4.3, we approximate the action of the Perron-Frobenius group on probability densities with by evaluating
[TABLE]
This approximation has an analogous consistency property to that in Lemma 25 if the converge to ; otherwise, if the converge to , approximates the Perron-Frobenius group through a different semigroup than .
As with the schemes of Sections 2.4, we evaluate (48) and (50) by forming the generator matrices and , and computing and , respectively, via Leja interpolation, where , , and , , are the expansion coefficients of and in the basis of , respectively. Explicit formulas for the elements of and in the case are provided in A.8.
5 Demonstration in flows driven by Lorenz 96 systems
5.1 Model description
In this Section, we apply the methods presented in Section 4 to a class of time-dependent incompressible flows introduced by Qi and Majda [34]. In this class of flows the spatial domain is doubly periodic, i.e., as in the examples of Section 3. Moreover, in canonical angle coordinates the velocity field has the general form
[TABLE]
where is the state in (to be specified below), and are continuous functions on taking values in and satisfying for every . Thus, by construction, has vanishing divergence with respect to the normalized Haar measure on for all . Based on earlier work by Majda et al. [63, 64], in [34] they choose the functions and so as to produce a flow involving a large-scale (spatially homogeneous) cross-sweep flow along the direction and a spatially inhomogeneous shear along the direction. Specifically,
[TABLE]
where is a positive integer, are complex-valued functions on , and , , are the Fourier functions on . Note that (51) determines the Fourier expansion coefficients and in the generator formula for arbitrary state-dependent velocity fields in (57); that is,
[TABLE]
The model in (51) was introduced as a low-order model for passive tracers in a jet; correspondingly, the dynamics driving the evolution of the Fourier coefficients was chosen so as to mimic some of the key features of turbulent tracer advection such as intermittency and non-Gaussian statistics. In particular, in [34] the evolution of the Fourier coefficients is driven by a L96 system with degrees of freedom, which we take here to be , with , , and
[TABLE]
This system was introduced by Lorenz in 1996 [65] as a low-order model for atmospheric flow at a constant latitude circle, with the variables representing the zonal (west–east) component of the velocity field at zonally equispaced gridpoints.
Formally, we assume that for each parameter choice used below, the L96 system admits an invariant ergodic measure supported on a compact (invariant) set , which is physical with respect to the Lebesgue measure of . That is, given solutions of the L96 system sampled discretely in time at the sampling interval , the set
[TABLE]
called the basin of , is assumed to have positive Lebesgue measure in . Note that unlike the Lorenz 63 model, where this property is known to hold [66], to our knowledge the existence of physical measures for L96 systems is an open problem.
Under these assumptions, we can construct a skew-product flow on governing the evolution of Lagrangian tracers as described in Section 2.1, and approximate the associated Koopman generator from long-time integrations of the L96 driving system. In particular, the assumed compactness of ensures that the kernel integral operators and from Section 4 used to build our data-driven basis are well defined and compact. However, is not expected to be smooth, meaning that the various limits taken in Section 4 to establish convergence of our data-driven basis functions and Dirichlet energies to corresponding quantities associated with heat kernels on smooth manifolds may not exist. Nevertheless, according to Remark 23, our schemes for coherent pattern identification and prediction may still converge away from that limit.
Under the assumption that such a set exists, we choose it as the state space of the ergodic system driving the non-autonomous flow on . Moreover, following [34], we define such that for ,
[TABLE]
Thus, the functions used to define the state-dependent velocity field in (51) are obtained by taking the discrete Fourier transform of the L96 variables . Taking also the generator of the dynamics on to be the restriction of the L96 vector field from (52) completes the specification of the skew-product system on with vector field and invariant measure . Note that the specification of via (51) and (52) also defines the observation map from Section 2.1; i.e., . By standard properties of Fourier transforms, we have so that computing pairwise distances between velocity field snapshots is equivalent to computing pairwise distances between L96 state vectors.
As discussed in A.8, given a data-driven basis of (computed as described in Section 5.2 ahead) and the corresponding triple product coefficients from (37), the matrix elements of the generator appearing in the schemes of Section 4 can be computed through (59) and (60). In the case of the L96-driven flow in (51), the latter equation becomes
[TABLE]
where
[TABLE]
and is a spectral truncation parameter for the state-dependent velocity field . In what follows, we set .
5.2 Dynamical regimes and data-driven basis functions
The L96 model in (52) includes a number of important features of atmospheric dynamics, namely energy-preserving nonlinear advection (), damping (), and external forcing (). It is well known that as increases, the system transitions through different dynamical regimes ranging from steady, periodic, quasiperiodic, weakly chaotic, strongly chaotic, and fully turbulent [67]. Here, we work with the parameter values (corresponding to 41 degrees of freedom) and or ; representative solutions for these regimes, computed “on the attractor”, and the root mean square (RMS) amplitudes of the Fourier modes are shown in Figs. 10 and 11, respectively. Published numerical results on the L96 model [68] indicate that the case , corresponds to a weakly chaotic regime with Lyapunov attractor dimension . We were not able to find analogous results for the , case, but based on available numerical results for nearby regimes [69] and qualitative inspection of numerical integrations, it is likely that this regime exhibits quasiperiodic or low-dimensional chaotic behavior (referred to hereafter as quasiperiodic behavior). In particular, it can be seen from Figs. 10 and 11 that the regime is dominated by wavenumber 7 waves with negative phase velocity, whose duration varies aperiodically. Wavenumber 7 waves are also clearly visible in the regime, but now the system also exhibits pronounced intermittent excursions from coherent wavenumber 7 activity.
Numerically, we generate training data from our chosen dynamical regimes by integrating the L96 model (using Matlab’s ode45 solver) with initial conditions . Waiting for transients to decay over a long spinup time (approximately 5,000 natural time units), we sample the solution times every time units after the spinup period, where and 512,000 for and , respectively (Figs. 10 and 11 were produced using those samples). Thus, taking into account the assumption on physical measures made in Section 5.1, we are effectively assuming that for both of the above choices of the point , lies in the basin of some ergodic probability measure supported on an invariant compact set . In what follows, we use the shorthand notation , for the L96 training states for notational compatibility with the rest of the paper, but note that we do not require that the lie exactly on .
Using this training data, we compute data-driven basis functions and their associated eigenvalues as described in Sections 4.1 and 4.2. Representative eigenfunctions form these calculations, plotted as time series for the same portions of the training datasets as those shown in Fig. 10, are displayed in Fig. 12. There, it can be seen that in both the and cases the eigenfunctions capture a range of timescales generated by the L96 system. As discussed in Section 4.3.1, our procedure for tuning the kernel bandwidth also yields an effective intrinsic dimension for ; this was found to be 4.1 and 6.4, respectively for and 5, respectively. Note that the effective dimension for is smaller than the corresponding Lyapunov attractor dimension reported in [68] by a factor of 2.5. This is not too surprising given that if is indeed -dimensional our 512,000 sample dataset only provides a coarse sampling of that set.
Next, we examine how many basis functions are needed to accurately reconstruct the Fourier modes in (53) (and hence the state-dependent velocity field ). Figure 10 shows the absolute and relative RMS reconstruction error for the Fourier modes as a function of the spectral truncation parameter , computed as and , respectively, where . Two features that stand out in these results are that (1) for both and , the high-energy Fourier modes require fewer basis functions for the same amount of reconstruction error than the low-energy modes, and (2) as expected from the more chaotic nature of the dynamics in that regime, the Fourier modes for require a significantly larger number of basis functions for accurate reconstruction. The large number of basis functions required for accurate reconstruction of becomes an especially pertinent issue associated with computational cost since eventually we will be working with the tensor product basis on , where, e.g., for and Fourier resolution the number of basis functions becomes . As discussed in B, this issue can be mitigated in numerical implementations through the use of code that computes the action of the various operator matrices in our schemes on vectors without explicit formation of the matrices themselves and/or code that exploits any sparsity that these matrices might possess. Indeed, it is evident from (54) and (59) that the generator matrices associated with the L96-driven system are sparse, and the code used to obtain the results presented below takes advantage of that sparsity (see, in particular, B.1),
5.3 Coherent spatiotemporal patterns
We begin by considering the system. Figure 13 shows snapshots of the real parts of approximate Koopman eigenfunctions computed for this system through the variatonal eigenvalue problem in Definition 24. In this case, we used the spectral truncation and diffusion regularization parameters , and . For this choice of spectral truncation parameters the total number of basis functions in the tensor product basis is . Note that the number of Fourier basis functions in each dimension on , , is greater than the number of L96 modes; this is important since the generator of the dynamics on can generate smaller scales on observables than those present in the Fourier coefficients on the velocity field . The dynamical evolution of the numerical Koopman eigenfunctions in Fig. 13 can also be visualized more directly as a video in Movie 10. Table 4 lists the corresponding eigenvalues and Dirichlet energies.
In general, the numerical eigenfunctions for have the structure of propagating patterns, advected along the positive direction by the spatially uniform cross-sweep flow in (51). Qualitatively, we can identify four distinct classes of eigenfunctions, namely (1) eigenfunctions such as and that do not depend on the coordinate; (2) eigenfunctions such as and that depend on but are approximately constant on the instantaneous streamlines of the flow in ; (3) eigenfunctions such as , , , and that depend on both and and whose level sets form approximately straight lines in the plane; (4) eigenfunctions such as , , , and that have the structure of propagating wavepackets. Note that eigenfunctions of class 2 have small values of ( and for and , respectively) indicating that the level sets of these functions are nearly invariant under the skew-product flow. On the other hand, the rest of the eigenfunctions shown here have , indicating that they vary on the tracers more strongly, albeit still at a slow timescale compared to the natural timescale of the flow. Unlike the moving-vortex example in Section 3.1 (see also C), in the case of the L96-driven system we do not have a detailed justification of the properties of the numerical Koopman eigenfunctions. Nevertheless, at least in some cases, we can provide partial justification provided that certain simplifying assumptions are made.
First, consider the eigenfunctions of class 1. Since the velocity field component is uniform on for all , the vector field is projectible under the map such that . In particular, we have . Assume now that (1) is constant, and (2) the generator on has eigenfunctions. The validity of assumption 1 can be tested from the time series of in the training data; in the regime we find the mean and standard deviation of to be 1.46 and 0.038, respectively, indicating that can indeed be treated as a constant to a good approximation. Assumption 2 is supported from the quasiperiodic nature of the dynamics in the regime (see Fig. 10). Under these assumptions, the product of any eigenfunction of at frequency with any Fourier function at wavenumber will be an eigenfunction of at eigenvalue . Moreover, the pullback of that eigenfunction onto will be an eigenfunction of at the same eigenvalue having constant values at the fibers of (i.e., constant hypersurfaces), as observed in Fig. 13 and Movie 10.
Next, to interpret the eigenfunctions of class 2, observe that if is nonzero (as is the case for all in the regime), the instantaneous streamlines of time-dependent flows of the class (51) are closed. In particular, since is non-vanishing and independent of , the streamlines at fixed are solutions of the ODE
[TABLE]
and since has no wavenumber zero Fourier component, integration of the right-hand side shows that , which implies that the streamlines are closed. Thus, as done in the moving-vortex example (see C), we can consider a projection onto a circle transverse to the streamlines. In particular, we can define such that where is the intercept of the streamline at state passing through the point with the circle . Since is everywhere tangent to the streamlines by definition, the component of the generator vanishes under , and we have . Thus, every function on which is constant on is an eigenfunction of at eigenvalue zero, and that function pulls back to an eigenfunction of , also at eigenvalue zero, which is constant on the streamlines of the L96-driven flow. The numerical eigenfunctions of class 2 described above are consistent with this behavior.
While we do not have a justification of the properties of class 3 and 4 eigenfunctions, it is worthwhile noting that the spatially localized character of the latter suggests that they are associated with the continuous spectrum of the generator (cf. Section 3.1 and C).
We now consider the regime . For the reasons discussed in B.4, our current code and available computational resources did not allow us to work with as high values of the spectral truncation parameter as the results in Fig. 11 indicate are required for accurate reconstruction of the state-dependent velocity field . Thus, aiming for a compromise between resolution on and , we decreased and to 20, which allowed us to increase to 500. With these parameter choices, the number of Fourier modes per dimension of is equal to the number of degrees of freedom of the L96 system; that is, we are fully resolving the velocity field but not any smaller scales that could develop as a result of the action of that velocity field on tracers. Moreover, the total number of basis functions is . Note that even though in this case is smaller than in the experiments, the memory requirements are actually higher due to the presence of -sized blocks of nonzero elements in the generator matrix ; see (54).
Figure 13, Movie 11 and Table 4 show representative Koopman eigenfunctions and eigenvalues computed with these parameters for the regime. The compromises we had to make on spectral resolution somewhat diminish our confidence in these results, but nevertheless some of he qualitative features identified in the experiments are also present here. In particular, the leading eigenfunction, , exhibits the characteristics of class 2 eigenfunctions discussed above; i.e., its level sets tend to follow the instantaneous streamlines of the flow (though with discrepancies at times, likely due to spectral truncation) and the corresponding eigenfrequency is small, .
Eigenfunctions broadly resembling the class 1 eigenfunctions at are also present at (e.g., , , and ), although they differ notably from the previous regime in that they display significant amplitude modulations. This is actually not too surprising, for our interpretation of class 1 eigenfunctions relied on the simplifying assumptions that the cross-sweep flow component is constant and that the generator of the dynamics on has eigenfunctions. At , the mean and standard deviation of are 1.67 and 0.17, respectively; i.e., the strength of fluctuations relative to the mean is about an order of magnitude stronger than in the case. Moreover, the assumption that has eigenfunctions becomes more questionable at . Thus, while we cannot exclude the possibility that the amplitude fluctuations in , , and are due to spectral truncation, our results are consistent with these functions being associated with the continuous spectrum of the system, exhibiting amplitude modulations as a result of localization on .
5.4 Prediction of observables and probability densities
We now examine prediction of observables and probability densities using the methods of Section 4.4.3. In what follows, we discuss a suite of numerical experiments that are analogous to those performed for the vortex flow models in Section 3; namely, we consider prediction of observables representing the position of tracers in and probability densities with a circular Gaussian structure at forecast initialization.
Since the data-driven basis function is constant by construction (see Section 4.2), we can define and in direct analogy with the corresponding observables in (22), and use the semigroup action from (48) to approximate the Koopman group action (recall that the latter gives the position at lead time of a tracer released from point when the L96 state is ). Figures 14 and 15 show the evolution of tracer positions for the systems with and , respectively, obtained via the operator-theoretic approach of Section 4.4.3 and explicit calculation by numerical ODE integration of the time-dependent flow in (51). These results were computed for initial L96 states in the training data; that is, in the operator-theoretic and ODE-based experiments we compute and for starting states in the training data, . This avoids out-of-sample extension errors that would occur if , so that the model error generated by the operator theoretic model is entirely due to diffusion regularization and spectral truncation. The evolution depicted in Figs. 14 and 15 is also more directly visualized as videos in Movies 12 and 13, respectively. In the case, we show results for 1,024 timesteps, i.e., forecast times up to 10.24 natural time units. In the case, we were only able to compute for 256 time steps with the computational resources available to us.
As expected from the structure of the velocity field in (51), in both the and cases the system generates a spatially uniform (but time-dependent) tracer flow along the direction and spatially nonuniform displacements along the direction, creating traveling-wave-like disturbances in observables and . As expected, the spatiotemporal structure of these disturbances is more complex in the case, but even at the system generates patterns at a broader range of spatial scales than that suggested by the dominant Fourier modes of the velocity field (see Fig. 11). In general, over the examined time interval, the results from the operator-theoretic model and the the Monte Carlo simulation are in very good agreement for , with only minor discrepancies becoming noticeable towards the end of the interval (). The results from the two methods are also in reasonably good agreement in the experiments, although there the operator theoretic model fails to capture some of the small-scale spatial patterns developed by the true model, even at times as early as . These discrepancies are likely due to both spectral truncation (notice in Fig. 11 that even at there are significant reconstruction errors in the velocity field) and chaotic dynamics of the L96 system. Nevertheless, the results in Fig. 15 and Movie 13 indicate that, at least over short to moderate lead times, the operator-theoretic model is able to perform useful forecasts of tracer positions.
Next, we consider prediction of probability densities. As stated above, we set the initial density on to a circular Gaussian, given by the same formula as (23). In the experiments that follow, we work throughout with the location and concentration parameters and , respectively. To construct an initial density on we fix an arbitrary state , and define , where is the Markov kernel from (30). In what follows, we choose to be a state “on the attractor” but outside of the training dataset . In particular, we compute the solution of the L96 system beyond the latest time in the training dataset, and set . Taking the product then leads to our initial density in . Note that being an element of , is unspecified for , but its expansion coefficients in the basis depend only on the values . Given these expansion coefficients, we approximate the evolution of the initial density through the semigroup action in (50). As in Section 3, we visualize the evolution of through the corresponding marginal densities , , and .
Figures 16 and 17 show the evolution of these marginal densities under the L96-driven flows with and , respectively, computed using the operator-theoretic approach of Section 4.4.3 and explicit Monte Carlo simulation for ensembles of 300,000 tracers drawn independently from . The spectral truncation parameters and simulation time intervals are identical to the corresponding values used for prediction of the and observables for and discussed above. A more direct visualization of the evolution in Figs. 16 and 17 is provided in Movies 14 and 15, respectively.
In both the and cases, the initially globular density is distorted by the velocity field fluctuations in the direction while simultaneously being advected along the positive direction by the spatially uniform cross-sweep flow. The spatial distortions of the initial density are more pronounced at , likely due to the more pronounced spectral peaks of the velocity field in that regime. In both cases, the 1D marginal density is advected along the direction with little change of shape, whereas undergoes oscillations in its variance due to the velocity field fluctuations in the direction. In general, the agreement between the operator-theoretic and Monte Carlo results is very good for . At , the agreement is also reasonably good (at least over the examined lead times), but biases in start becoming noticeable at ; as with prediction of , this is at least partly caused by spectral truncation.
6 Concluding remarks
In this paper, we have developed a family of techniques for coherent pattern identification and statistical prediction in measure-preserving, skew-product dynamical systems, with a particular emphasis on systems of Lagrangian tracers evolving under time-dependent fluid flows driven by ergodic dynamical systems. Central to these techniques is a data-driven representation of the generator of the Koopman and Perron-Frobenius groups of operators on an extended state space (as recently done in [20] in the case of periodic flows), given by the product of the state space of the dynamical system driving the flow and the physical domain in which the flow takes place. In particular, our approach is to represent the generator of the associated skew-product system (governing the joint evolution of the velocity field and the tracers) in a smooth basis of the space of , constructed as a tensor product of bases for the spaces on and . Importantly, the basis for is built from Laplace-Beltrami eigenfunctions learned from time-ordered velocity field snapshots using kernel algorithms for machine learning [28, 30], requiring no a priori knowledge of the geometry or the dynamics on , or availability of explicit tracer trajectories.
With the availability of this basis, we constructed a Galerkin method for the eigenvalue problem of a regularized generator for the skew-product system having only point spectrum [25], and used the eigenfunctions of this operator to define coherent spatiotemporal patterns. Following earlier work [22, 23], our approach has been to select eigenfunctions with minimal roughness; the latter measured through a Dirichlet energy functional accessible to data-driven approximation through the same kernel algorithms as those used to build the basis. Taking the exponential of the generator and its adjoint, we also constructed representations of the Koopman and Perron-Frobenius operators governing the evolution of observables and probability densities of the skew-product system, and used these representations in model-free statistical prediction schemes. Here, we computed the action of the generator exponential on functions by means of Leja polynomial interpolation algorithms [32, 33], used in exponential integrators for ODEs.
We demonstrated the utility of these techniques in periodic flows with Gaussian streamfunctions on two-dimensional periodic domains and in aperiodic flows driven by Lorenz 96 systems [34] in weakly chaotic regimes. The periodic-flow examples where chosen so as to have mixed spectra and generator matrices computable in closed form, allowing us to assess the properties of these schemes in the absence of sampling errors. We saw that in these flows the numerically computed eigenfunctions of the generator provide useful notions of coherent spatiotemporal patterns, with properties consistent with those expected from theory (when such results are available), and in some cases providing numerical access to the continuous spectrum through patterns behaving as “strange eigenmodes” [18, 19]. We also saw that our methods for prediction of observables and probability densities perform comparably to explicit integration of the full model, including in their ability to reproduce highly non-Gaussian probability density functions. In the case of L96-driven flows, our numerical results provide evidence that our framework is also applicable when the state space of the driving system is not a smooth manifold, but rather a more general compact invariant set of a dynamical system. Here, a major challenge was that the L96 driving system can require large numbers of basis functions in order to accurately represent the generator, especially in chaotic regimes with many positive Lyapunov exponents. Nevertheless, we found that our approach yields reasonably good results, both in terms of coherent pattern detection and prediction of observables and densities.
Two aspects of this work that warrant future study and potential improvement are (1) the way that diffusion regularization interacts with the skew-adjoint Koopman generator in the presence of continuous spectrum, and (2) the high computational cost associated with the tensor product basis for the Hilbert spaces encountered in skew-product systems. The first issue was addressed in a specific case in [23, 24] where it was shown that through an appropriate use of delay-coordinate maps one can construct a data-driven diffusion operator the commutes with the Koopman generator on the discrete subspace of . However, that work did not address how objects such as strange eigenmodes behave when the diffusion operators and the Koopman generator do not commute (which, to our knowledge, is an open problem in PDE theory). Possible modifications of our approach to address the second issue would be to replace the tensor product basis with a basis constructed directly on the extended state space of the skew-product system, possibly by using explicit tracer trajectories. Alternatively, one could seek efficiency improvements by employing multi-resolution bases [70, 71], or through sparse representations of the generator [72]. We plan to pursue these topics in future work.
Acknowledgments
Dimitrios Giannakis acknowledges support by 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 wish to thank Pierre Germain and Igor Mezić for stimulating conversations.
Appendix A Technical results
A.1 Proof of Proposition 6
(i) That are dissipative follows immediately from the skew-symmetry of and positivity of . The latter imply that for any , and , leading to . It is a classical result [42] that every densely-defined dissipative operator is closable and and its closure is dissipative.
(ii) Let and be arbitrary functions in . It follows by integration by parts in conjunction with the skew-symmetry of and the symmetry of that . Since this relationship holds for any , we can conclude that and . It is also straightforward to check that is not closed. Together, these results in conjunction with the fact that is closed imply that . The result follows similarly.
(iii) The existence of follows from [73, Chapter IV, Proposition 4.2] in conjunction with the fact that is dissipative (in their terminology, accretive). That is maximally dissipative is equivalent to also being maximally dissipative [42, 74], leading to the second claim.
(iv) First, observe that for any function we have
[TABLE]
where , , and we have used the fact that . Next, let be any function in , and consider a Cauchy sequence converging to in norm. Equation (55), in conjunction with the fact that is Cauchy in , implies that is Cauchy in . Therefore, since is closed, converges to a point in the graph of , where . Hence, we can conclude that
[TABLE]
as claimed. The result for follows similarly.
(v) A proof of the claim can be found in [25, Proposition 1].
A.2 Proof of Lemma 10
(i) Given , there exist unique coefficients such that . For later convenience, let , and define the norm with
[TABLE]
so that . By equivalence of norms on finite-dimensional vector spaces, there exist constants and such that , where is the canonical 2-norm on . Thus, for we have
[TABLE]
from which it follows that .
(ii) Since is assumed to be basis of , for any there exist unique coefficients such that . Moreover, there exist constants and such that for any . This implies that for any , but because is unbounded and is assumed to have no accumulation points, this bound cannot hold unless there exists such that for all . Thus, it follows that , as desired.
A.3 Proof of Lemma 13
That is compact follows immediately from the fact that it has finite rank. Showing that is compact is equivalent to showing that for any sequence with bounded uniform norm , the sequence has a limit point in the uniform-norm topology. Since the underlying space is compact, this will follow from the Arzela-Ascoli theorem, provided that it can be shown that is equicontinuous and uniformly bounded. Indeed, uniform boundedness and equicontinuity follow from the inequalities
[TABLE]
respectively, proving the Lemma.
A.4 Proof of Lemma 15
(i) The claim follows from Proposition 11 in [29].
(ii) The smoothness of the mapping follows from the smoothness of . In particular, using the results in Appendices A.4 and A.5 in [30] and the relation between the Laplace-Beltrami operator associated with and the ambient space metric in (4), it is possible to derive the uniform asymptotic approximation
[TABLE]
where ,
[TABLE]
are Gaussian integrals, and is a smooth function. Using this expansion in the normalization steps in (28) and (29) (as done in [28, 60]) leads to the desired result. In particular, note that the term cancels as a result of normalization.
(iii) The claim was proved in [28] in the case of radial (fixed-bandwidth) Gaussian kernels. The corresponding result for the variable-bandwidth kernels in (26) can be obtained following the same approach, employing the asymptotic expansion in (32) in place of the corresponding expansion for radial Gaussian kernels as necessary.
(iv) Introducing the multiplication operators and with , , and , , it follows immediately that and . Note that we can also relate and to the corresponding unnormalized operators and via
[TABLE]
To establish convergence of to , note first that, as , converges in operator norm to the multiplication operator associated with the positive, constant function . Correspondingly, and converge to and , respectively. We thus have , and
[TABLE]
where denotes the commutator of operators. As , both terms in the right-hand side of the last inequality vanish.
A.5 Proof of Lemma 16
(i, ii) The claims are obvious.
(iii) Since is nonzero, it follows from part (ii) that it is an eigenvalue of , and is an eigenfunction. As shown in [29], because converges (in this case, -a.s.) to (Lemma 15(i)), there exist eigenvalues of and corresponding eigenfunctions converging as to and , respectively. It then follows from part (i) that and are also eigenvalues and eigenfunctions of , respectively. Thus, converges to , as desired. Moreover, since , there exists a stage such that for all we have , and therefore (by part (ii)) . The fact that in uniform norm completes the proof of the claim.
(iv) Since converges to the heat operator (Lemma 15(iv)), we have , which proves (35). Similarly, the existence of a family of eigenfunctions converging to follows from convergence of to .
(v) The results follow by obvious modifications of the arguments in the proof of parts (i)–(iv).
A.6 Proof of Proposition 19
(i) As in the case of the Dirichlet form in Section 2.3.4, the claim follows from the Cauchy-Schwartz inequality and the facts that and .
(ii) In the case of , the claim follows immediately from the fact that the operator is symmetric and is finite-dimensional. In the case of , since the operator is symmetric and densely defined (by Lemma 18), is an extension of . It therefore suffices to show that is also an extension of , i.e., . In particular, it suffices to show that for every there exists such that . The last equation is satisfied for , proving that is self-adjoint. The fact that and are dissipative is obvious from the definition of these operators in (39).
A.7 Proof of Lemma 21
First, note that and are equivalence classes of functions, so they lie in the domain of and the right-hand side of (43) is finite. Next, consider the bounded operator defined via , and the functions and (the existence of follows from Lemma 16). Then, using Corollary 17, we obtain
[TABLE]
Since and are both in the domain of , and the inner product weight function for satisfies (see Section 4.2), we have
[TABLE]
leading to the desired result.
A.8 Explicit formulas for the generator for periodic spatial domains
Since all the numerical experiments in this paper (Sections 3 and 5) were performed on a doubly-periodic spatial domain, , equipped with the normalized Haar measure , in this Appendix we provide explicit formulas for the matrix representation of the generator in this domain for arbitrary state-dependent incompressible velocity fields . Throughout, we work in the tensor product basis , , with , , set to the data-driven basis functions from Section 4, and set to the Fourier basis, , .
A special property of (with obvious generalizations to tori of other dimensions) is that the canonical basis vector fields and span the tangent space at every . Thus, we can expand , , as
[TABLE]
where and are smooth, complex-valued functions on . The corresponding equivalence classes and in can be further expanded as
[TABLE]
with and .
Next, we make use of an important property of the Fourier basis, namely that is an eigenfunction of both and with and . Using this property in conjunction with the expansion in (57), we obtain (cf. (18))
[TABLE]
where the coefficients are given in (37). Moreover, we have (cf. (20))
[TABLE]
where the quantities are evaluated by finite differences using (42).
Equations (58) and (59), together with the definition of the basis of and the facts that and are in principle sufficient to form the matrices appearing in the eigenvalue problem in Definition 24. Equations (19) and (20) in conjunction with are also in principle sufficient to compute the generator matrices and employed in the data-driven prediction schemes for observables and densities in Section 4.4.3. In practice, however, evaluation of (58) is likely to be unwieldy since it requires computation of eigenfunctions , which may not be feasible at large . To address this issue, we take advantage of the expansion of in (57) in terms of the globally defined vector fields and , defining
[TABLE]
where is a spectral truncation parameter (usually, ). Using this smoothed version of the velocity field, we compute
[TABLE]
and implement the schemes of Section 4 using these quantities instead of (58). Note that this approximation does not affect the asymptotic consistency of our schemes, so long as an appropriate increasing sequence of is employed; e.g., . All of the numerical experiments of Section 5 were performed with the latter choice for .
Appendix B Numerical implementation
All numerical experiments presented in this paper were carried out in Matlab running on a medium-size Linux cluster with “big memory” capabilities (up to 1.5 TiB of memory per compute node). In the case of the L96-driven experiments of Section 5, we also interface with a Fortran 77 code (using Matlab’s MEX framework) to carry out matrix-vector products involving the Koopman and Perron-Frobenius generators used in Leja interpolation. It should be noted that while big-memory capabilities are highly convenient for numerical implementation of our schemes, in many cases they are not essential. In particular, it is frequently the case that matrix-vector products involving the generator can be evaluated with tolerable computational overhead without explicit formation of the generator matrix itself; our Fortran-based implementation for the L96-driven system is an example of this approach.
B.1 Leja interpolation
We have implemented the data-driven prediction schemes of Sections 2.4 and 4.4.3 using a publicly available implementation of the Leja method in Matlab authored by Caliari and Kandolf.111Available at https://www.mathworks.com/matlabcentral/fileexchange/44039-matrix-exponential-times-a-vector. In addition to the Leja interpolation procedure, this Matlab code also performs a subdivision step, expressing as for a positive integer parameter , and repeatedly performing Leja interpolation to approximate each of the substeps for , , and . This subdivision procedure is introduced in order to overcome the so-called “hump problem” [48], which is a transient error growth with the polynomial order of approximation sometimes observed at large . The code uses an automatic procedure to determine the number of substeps , as well as in each substep, needed to achieve a desired accuracy.
As stated above, in the case of the L96-driven flow, we interface the code of Caliari and Kandolf with a Fortran code to evaluate matrix-vector products involving the generator matrix (computed via (54) and (59)). Due to the sparse structure of that matrix, it is possible to evaluate matrix-vector products without explicit formation of , at the expense of a moderate computation cost overhead. This is particularly advantageous for reducing memory use at large values of the spectral truncation parameter , but typically requires coding the evaluation of using nested for loops which are inefficient in an array-oriented language such as Matlab. On the other hand, these loops can be efficiently coded in a compiled language such as Fortran using multithreading and other high-performance features, allowing scalability of the method to large . A similar approach can be applied in more general settings to take advantage of any special structure that the generator matrix might possess.
B.2 Density estimation
Following [30, 62, 22, 23], we compute the kernel bandwidth functions in (26) using kernel density estimation. First, we introduce an integer-valued function such that is equal to the time index of the -th nearest point in the training dataset to an arbitrary data point measured with respect to the norm . Fixing a positive integer parameter , we then define the bandwidth function such that and the kernel given by
[TABLE]
Using this kernel, we select a value for the bandwidth parameter and compute an estimate of the manifold dimension using the approach outlined in Section 4.3.1. We then compute the function given by
[TABLE]
and define the bandwidth function used in the variable-bandwidth kernel in (26) as . Note that the optimal bandwidth parameter and the dimension estimate from are generally not equal to the corresponding values computed using , though the dimension estimates from the two approaches are generally in reasonably good agreement. Throughout this paper we work with the neighborhood size which was also used in [62, 22, 23].
By compactness of and smoothness of , and are both smooth functions for all . Moreover, as , and converge pointwise and -a.s. to smooth functions and , respectively, and by properties of Gaussian integrals on compact manifolds, converges as to the sampling density , as desired. Furthermore, the error is uniformly ; the latter is a necessary condition for Lemma 15(ii) to hold [30]. Note that besides the density estimation procedure outlined above any other technique achieving this level of accuracy can be employed in the definition of the bandwidth functions in (26).
B.3 Kernel eigenvalue problems
The eigenvalue problems for the finite-rank operators and are equivalent to the eigenvalue problems for the matrices and ; that is, in the notation of Section 4.2, we have
[TABLE]
Numerically, it is generally preferable to solve the eigenvalue problem for to take advantage of the symmetry of that matrix, and compute the eigenvectors of via , where is the diagonal matrix with diagonal elements . The numerical results presented in this paper were obtained via this approach, using Matlab’s eigs solver for sparse arrays. For symmetric matrices such as , this solver utilizes an implicitly restarted Lanczos method implemented in the ARPACK library [43].
As is customary, we take advantage of the exponential decay of the kernel and sparsify . Specifically, we start by computing the kernel values on the training dataset, and for each state retain the largest values , where is a neighborhood size parameter (different from the one used for density estimation in B.2). Using these values, we form an , sparse, symmetric kernel matrix whose nonzero entries are if is in the neighborhood of or is in the neighborhood of . We then perform the diffusion maps normalization procedure as described in Section 4.2 to construct based on . A pseudocode for this procedure can be found in [23]. In the L96 experiments of Section 5 we use the value which amounts to approximately and of the training datasets for the and cases, respectively. A Matlab code for computing and its eigenvalues and eigenvectors is available for download at the first author’s personal website (http://cims.nyu.edu/~dimitris).
B.4 Koopman eigenvalue problem
As discussed in the main text, the Koopman eigenvalue problems in Definitions 8 and 24 are equivalent to matrix generalized eigenvalue problems of the form (8), involving the stiffness and mass matrices and , respectively. An obvious issue with these eigenvalue problems is that the tensor product bases of and employed in our schemes generally require large values of the spectral truncation parameter for accurate numerical solutions, leading to high computation and memory cost. Arguably the most straightforward way of overcoming this issue (at least in the context of a Matlab implementation) is to take advantage of any sparsity that might possess (note that is trivially sparse since our bases are orthogonal). Such sparsity is present in both the vortical flow examples of Section 3 (see (19), (20), and (25)) and the L96-driven flows of Section 5 (see (54) and (59)). In such cases, so long as can be stored in memory as a sparse array, (8) can be solved efficiently using Matlab’s eigs ARPACK-based solver. This solver has an option LR which allows one to compute a subset of the eigenvalues with the largest real part, which, according to Remark 9 is consistent with our Dirichlet energy ordering criterion in Definition 7. However, in practice we find that this option can be slow to converge, and faster convergence is generally achieved using option SM targeting eigenvalues with the smallest modulus. All of the eigenvalues and eigenfunction results reported in this paper were computed using the latter approach, which carries a risk of missing certain eigenvalues with small real part but large modulus.
As the problem size increases, the approach described above requiring explicit formation of , will eventually become infeasible. For example, in the case of the L96-driven flow in Section 5 with , we were not able to increase the spectral truncation parameters and beyond 500 and 20, respectively, despite the sparsity of the generator matrix. This limitation can in principle be overcome using iterative solvers such as ARPACK, since what is required by such solvers are functions for computing matrix-vector products or solving linear systems involving and , and these functions can be coded so as not to require explicit formation of (see B.1). In practice, however, we found that attempting to compute via this approach in Matlab’s eigs solver led to unstable behavior, with ARPACK iterations frequently failing to converge. Addressing this issue is beyond the scope of this work, but with appropriate technical modifications (perhaps avoiding use of high-level functions such as eigs) it should be possible to overcome the problem size limitations our current code experiences.
Appendix C The spectrum of the generator for the moving Gaussian vortex flow
In this Appendix, we discuss the properties of the eigenfunctions of class 2 for the moving Gaussian vortex flow in Section 3.1, and also establish the existence of a continuous part of the spectrum for this system. First, note that the streamfunction in (16) is invariant under transformations of that leave unchanged. In particular, we can define the transformation such that with and , and consider that is the pullback of the streamfunction on , where and are canonical angle coordinates. That is, we have . This transformation can be thought of as a Galilean change of frame of reference to a frame comoving with the vortex center. In particular, the generating vector field of the dynamics on is projectible under that transformation; i.e., is a well-defined vector field on , and its components in the coordinate basis are given by
[TABLE]
respectively. Denoting the flow on generated by by , one can check that . Therefore, preserves the pushforward measure , and the spectrum of includes the spectrum of (an appropriate skew-adjoint extension of) . In particular, if is an eigenfunction of at eigenvalue , then is an eigenfunction of , also at eigenvalue . Note that the pushforward measure associated with the submersion is the normalized Haar measure on the two-torus, i.e., . (Note that in this section we abuse notation and do not employ a distinct symbol to represent the skew-adjoint extension of generating the Koopman group on , but the distinction is important since the spectrum of a non-closed operator such as is equal to ; see, e.g., [75, Section 2.2].)
Inspecting (61), we see that the flow on consists of a superposition of a stationary Gaussian vortex flow and a free stream (i.e., a harmonic velocity field; see Example 1) along the negative direction associated with the Galilean transformation . In this domain, the equation of motion for Lagrangian tracers is , where . For our choice of the and parameters, is everywhere negative so that tracers move along the negative direction without stagnating or turning around. This, in conjunction with the fact that possesses the reflection symmetry , implies that the streamlines of the flow are closed. In particular, the coordinate of the streamline passing through the point can be expressed as a function by solving the ODE
[TABLE]
The collection of all such streamlines partitions into a one-dimensional foliation such that the vector field is tangent to the leaves; that is, is a vertical vector field with respect to this foliation [76]. In particular, picking the closed line , we can construct a submersion sending every point in streamline passing through to the point . The preimages of this map correspond to the leaves of the foliation, and are streamlines by construction.
Since is vertical, we have , so that is trivially projectible and the projected “dynamics” on reduce to the identity map (i.e., ). As a result, the elements of any orthonormal basis of (note that is the normalized Haar measure on ) are orthonormal eigenfunctions of at eigenvalue 0, and correspondingly the functions are orthonormal eigenfunctions of in , also at eigenvalue 0. Composing with , we obtain the submersion , , and the corresponding subspace discussed in Section 3.1. Note that this submersion is associated with a two-dimensional foliation of .
Next, to establish that , and hence , also have continuous spectrum, it suffices to demonstrate that is conjugate by diffeomorphism to a skew rotation on [52]. In canonical angle coordinates , the vector field generating this skew rotation takes the form
[TABLE]
where is a smooth function on . Note that the flow generated by (62) takes place along lines of constant coordinate with an -dependent frequency. To map the original flow taking place along curved streamlines to this skew rotation, fix canonical angle coordinates on with , and define to be the time taken for a Lagrangian tracer released at the point (i.e., the intercept of the streamline passing through with the line) to reach the point . Define also to be the time taken for the tracer released at to first return to this point; i.e., . Note that in the domain of definition of the coordinates we have . With these definitions, we introduce the diffeomorphism such that , where
[TABLE]
Because and are constant along the streamlines in , this map rectifies these streamlines into “straight” lines in at constant coordinate. Moreover, is projectible under this map, and because , it maps to the vector field such that with (note that the components of in the coordinate basis are given by , and ). It is also possible to verify that , which implies (in conjunction with the fact that is a diffeomorphism) that and have the same spectra.
The skew rotation generated by the vector field in (62) is a special case of the flows studied by Broer and Takens in [52]. There, it is shown that for nonconstant the generator has continuous spectrum associated with nonconstant functions along the integral curves of (i.e., the lines ). Intuitively, if were an eigenfunction of then there should exist a time such that is mapped back into itself, but this is not possible if is nonconstant and , for the the values of at different coordinates map back to themselves at different times .
The subspace of functions in (where is the normalized Haar measure) associated with the continuous spectrum is the orthogonal complement of . Correspondingly, the subspace of associated with the continuous spectrum of is , and pulling back this subspace to we obtain a subspace associated with the continuous spectrum of the generator of the full skew-product system. This subspace lies in the orthogonal complement of the discrete spectrum subspace identified in Section 3.1; in fact, we also have the decomposition (note that for any and the product lies in [24]). In summary, the full space of the skew-product system on associated with the moving vortex flow admits the decomposition
[TABLE]
References
- [1]
J. M. Ottino, The Kinematics of Mixing: Stretching, Chaos, and Transport, no. 3 in Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, 1989.
- [2]
R. M. Samelson, S. Wiggins, Lagrangian Transport in Geophysical Jets and Waves: The Dynamical Systems Approach, Vol. 31 of Interdisciplinary Applied Mathematics, Springer, New York, 2006.
- [3]
G. Haller, Lagrangian coherent structures, Annu. Rev. Fluid Mech. 47 (2015) 137–162.
doi:10.1146/annurev-fluid-010313-141322.
- [4]
G. Haller, G. Yuan, Lagrangian coherent structures and mixing in two-dimensional turbulence, Phys. D 147 (2000) 352–370.
doi:10.1016/S0167-2789(00)00142-1.
- [5]
S. C. Shadden, F. Lekien, J. E. Marsden, Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows, Phys. D 212 (2005) 271–304.
doi:10.1016/j.physd.2005.10.007.
- [6]
G. Haller, A variational theory of hyperbolic Lagrangian coherent structures, Phys. D 240 (2011) 574–598.
doi:10.1016/j.physd.2010.11.010.
- [7]
M. Serra, G. Haller, Objective Eulerian coherent structures, Chaos 26 (2016) 053110.
- [8]
G. Froyland, K. Padberg, M. H. England, A. M. Treguier, Detection of coherent oceanic structures via transfer operators, Phys. Rev. Lett. 98 (22) (2007) 224503.
doi:10.1103/PhysRevLett.98.224503.
- [9]
G. Froyland, S. Lloyd, N. Santitissadeekorn, Coherent sets for nonautonomous dynamical systems, Phys. D 239 (2010) 1527–1541.
doi:doi:10.1016/j.physd.2010.03.009.
- [10]
G. Froyland, N. Santitissadeekorn, A. Monahan, Transport in time-dependent dynamical systems: Finite-time coherent sets, Chaos 20 (2010) 0431116.
- [11]
G. Froyland, An analytic framework for identifying finite-time coherent sets in time-dependent dynamical systems, Phys. D 250 (2013) 1–19.
doi:10.1016/j.physd.2013.01.013.
- [12]
M. Budisić, I. Mezić, Geometry of the ergodic quotient reveals coherent structures in flows, Phys. D 241 (2012) 1255–1269.
doi:10.1016/j.physd.2012.04.006.
- [13]
G. Froyland, K. Padberg, Almost-invariant sets and invariant manifolds – connecting probabilistic and geometric descriptions of coherent structures in flows, Phys. D 238 (2009) 1507–1523.
doi:10.1016/j.physd.2009.03.002.
- [14]
G. Froyland, Dynamic isoperimetry and the geometry of Lagrangian coherent structures, Nonlinearity (2015) 3587–3622doi:10.1088/0951-7715/28/10/3587.
- [15]
D. Karrasch, J. Keller, A geometric heat-flow theory of Lagrangian coherent structures.
- [16]
R. Banisch, P. Koltai, Understanding the geometry of transport: Diffusion maps for Lagrangian trajectory data unravel coherent sets, Chaos 27 (2017) 035804.
- [17]
G. Froyland, O. Junge, Robust FEM-based extraction of finite-time coherent sets using scattered, sparse, and incomplete trajectories (2017).
- [18]
W. Liu, H. Haller, Strange eigenmodes and decay of variance in the mixing of diffusive tracers, Phys. D 188 (2004) 1–39.
doi:10.1016/S0167-2789(03)00287-2.
- [19]
R. T. Pierrehumbert, Tracer microstructure in the large-eddy dominated regime, Chaos Soliton. Fract. 4 (6) (1994) 1091–1110.
doi:10.1016/0960-0779(94)90139-2.
- [20]
G. Froyland, P. Koltai, Estimating long-term behavior of periodically driven flows without trajectory integration, Nonlinearity 30 (2017) 1948–1986.
- [21]
T. Eisner, B. Farkas, M. Haase, R. Nagel, Operator Theoretic Aspects of Ergodic Theory, Vol. 272 of Graduate Texts in Mathematics, Springer, 2015.
- [22]
D. Gannakis, J. Slawinska, Z. Zhao, Spatiotemporal feature extraction with data-driven Koopman operators, J. Mach. Learn. Res. Proceedings 44 (2015) 103–115.
- [23]
D. Giannakis, Data-driven spectral decomposition and forecasting of ergodic dynamical systems, Appl. Comput. Harmon. Anal.In minor revision.
- [24]
S. Das, D. Giannakis, Delay-coordinate maps and the spectra of Koopman operators, NonlinearityIn preparation.
- [25]
B. Franke, C.-R. Hwang, H.-M. Pai, S. Sheu, The behavior of the spectral gap under growing drift, Trans. Amer. Math. Soc. 362 (3) (2010) 1325–1350.
doi:10.1090/s0002-9947-09-04939-3.
- [26]
P. Constantin, A. Kiselev, L. Ryzhik, A. Zlatoš, Diffusion and mixing in fluid flow, Ann. Math. 168 (2008) 643–674.
URL http://www.jstor.org/stable/40345422
- [27]
M. Belkin, P. Niyogi, Laplacian eigenmaps for dimensionality reduction and data representation, Neural Comput. 15 (2003) 1373–1396.
doi:10.1162/089976603321780317.
- [28]
R. R. Coifman, S. Lafon, Diffusion maps, Appl. Comput. Harmon. Anal. 21 (2006) 5–30.
doi:10.1016/j.acha.2006.04.006.
- [29]
U. von Luxburg, M. Belkin, O. Bousquet, Consitency of spectral clustering, Ann. Stat. 26 (2) (2008) 555–586.
doi:10.1214/009053607000000640.
- [30]
T. Berry, J. Harlim, Variable bandwidth diffusion kernels, J. Appl. Comput. Harmon. Anal. 40 (1) (2015) 68–96.
doi:10.1016/j.acha.2015.01.001.
- [31]
T. Berry, T. Sauer, Local kernels and the geometric structure of data, J. Appl. Comput. Harmon. Anal. 40 (3) (2016) 439–469.
doi:10.1016/j.acha.2015.03.002.
- [32]
M. Caliari, M. Vianello, L. Bergamaschi, Interpolating discrete advection–-diffusion propagators at Leja sequences, J. Comput. Appl. Math. 172 (79–-99).
doi:10.1016/j.cam.2003.11.015.
- [33]
P. Kandolf, A. Ostermann, S. Rainer, A residual based error estimate for Leja interpolation of matrix functions, Linear Alg. Appl. 456 (2014) 157–173.
- [34]
D. Qi, A. J. Majda, Predicting fat-tailed intermittent probability distributions in passive scalar turbulence with imperfect models through empirical information theory, Commun. Math. Sci. 14 (6) (2016) 1687–1722.
doi:10.4310/CMS.2016.v14.n6.a11.
- [35]
R. Abraham, J. E. Marsden, T. Ratiu, Manifolds, Tensor Analysis, and Applications, Vol. 75 of Applied Mathematical Sciences, Springer, New York, 1988.
- [36]
T. Sauer, J. A. Yorke, M. Casdagli, Embedology, J. Stat. Phys. 65 (3–4) (1991) 579–616.
- [37]
I. Mezić, Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dyn. 41 (2005) 309–325.
doi:10.1007/s11071-005-2824-x.
- [38]
M. Budisić, R. Mohr, I. Mezić, Applied Koopmanism, Chaos 22 (2012) 047510.
- [39]
I. Babuška, J. Osborn, Eigenvalue problems, in: P. G. Ciarlet, J. L. Lions (Eds.), Finite Element Methods (Part 1), Vol. II of Handbook of Numerical Analysis, North-Holland, Amsterdam, 1991, pp. 641–787.
- [40]
M. Dellnitz, O. Junge, On the approximation of complicated dynamical behavior, SIAM J. Numer. Anal. 36 (1999) 491.
doi:10.1137/S0036142996313002.
- [41]
S. Rosenberg, The Laplacian on a Riemannian Manifold, Vol. 31 of London Mathematical Society Student Texts, Cambridge University Press, Cambridge, 1997.
- [42]
R. S. Phillips, Dissipative operators and hyperbolic systems of partial differential equations, Trans. Amer. Math. Soc. 90 (2) (1959) 193–254.
- [43]
R. B. Lehoucq, D. C. Sorensen, C. Yang, ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems With Implictly Restarted Arnoldi Methods, SIAM, Philadelphia, 1998.
- [44]
M. Reed, B. Simon, Methods of Modern Mathematical Physics. I: Functional Analysis, Academic Press, San Diego, 1980.
- [45]
E. Nelson, Analytic vectors, Ann. Math. 70 (3) (1959) 572–615.
- [46]
F. H. Stochel, J. Szafraniec, -vectors and boundedness, Ann. Polon. Math. (1997) 223–238.
- [47]
H. Fujita, T. Suzuki, Evolution problems, in: P. G. Ciarlet, J. L. Lions (Eds.), Finite Element Methods (Part 1), Vol. II of Handbook of Numerical Analysis, North-Holland, Amsterdam, 1991, pp. 789–923.
- [48]
C. Moler, C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Rev. 45 (1) (2003) 3–49.
- [49]
N. J. Higham (Ed.), Functions of Matrices, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 2008.
- [50]
M. Caliari, P. Kandolf, A. Ostermann, S. Rainer, Comparison of methods for computing the action of the matrix exponential, BIT Numer. Math. 52 (1) (2014) 113–128.
- [51]
M. Caliari, Accurate evaluation of divided differences for polynomial interpolation of exponential propagators, Computing 80 (2007) 189–201.
- [52]
H. Broer, F. Takens, Mixed spectra and rotational symmetry, Arch. Rational Mech. Anal. 124 (1993) 13–42.
- [53]
I. Mezić, Koopman operator spectrum and data analysis (2017).
- [54]
B. Watson, Manifold maps commuting with the Laplacian, J. Diff. Geom. 8 (1973) 85–94.
- [55]
S. I. Goldberg, T. Isihara, Riemannian submersions commuting with the Laplacian, J. Diff. Geom. 13 (1978) 139–144.
- [56]
P. B. Gilkey, J. V. Leahy, J. H. Park, Eigenvalues of the form valued Laplacian for Riemannian submersions, Proc. Amer. Math. Soc. 126 (6) (1998) 1845–1850.
- [57]
M. Putinar, Generalized eigenfunction expansions and spectral decompositions, Banach Center Publications 38 (1) (1997) 265–286.
- [58]
H. Aref, Stirring by chaotic advection, J. Fluid. Mech. 143 (1984) 1–21.
doi:dx.doi.org/10.1017/S0022112084001233.
- [59]
D. Giannakis, Dynamics-adapted cone kernels, SIAM J. Appl. Dyn. Sys. 14 (2) (2015) 556–608.
- [60]
T. Berry, T. Sauer, Consistent manifold representation for topological data analysis (2016).
- [61]
R. R. Coifman, Y. Shkolnisky, F. J. Sigworth, A. Singer, Graph Laplacian tomography from unknown random projections, IEEE Trans. Image Process. 17 (10) (2008) 1891–1899.
- [62]
T. Berry, D. Giannakis, J. Harlim, Nonparametric forecasting of low-dimensional dynamical systems, Phys. Rev. E. 91 (2015) 032915.
doi:10.1103/PhysRevE.91.032915.
- [63]
A. J. Majda, P. R. Kramer, Simplified models for turbulent diffusion: Theory, numerical modeling, and physical phenomena, Phys. Rep. (1999) 237–574doi:10.1016/S0370-1573(98)00083-0.
- [64]
A. J. Majda, B. Gershgorin, Quantifying uncertainty in climage change science through empirical information theory, Proc. Natl. Acad. Sci. 107 (34) (2010) 14958–14963.
- [65]
E. N. Lorenz, Predictability of weather and climate, in: T. Palmer, R. Hagedorn (Eds.), Predictability of Weather and Climate, Cambridge University Press, Cambridge, 1996, Ch. 3, pp. 40–58.
- [66]
S. Luzzatto, I. Melbourne, F. Paccaut, The Lorenz attractor is mixing, Comm. Math. Phys. 260 (2) (2005) 393–401.
- [67]
A. J. Majda, R. V. Abramov, M. J. Grote, Information Theory and Stochastics for Multiscale Nonlinear Systems, Vol. 25 of CRM Monograph Series, American Mathematical Society, Providence, 2005.
- [68]
A. Karimi, M. R. Paul, Extensive chaos in the Lorenz-96 model, Chaos (2010) 043105doi:10.1063/1.3496397.
- [69]
D. L. van Kekem, A. E. Sterk, Travelling waves and their bifurcations in the Lorenz-96 model (2017).
- [70]
R. R. Coifman, M. Maggioni, Diffusion wavelets, Appl. Comput. Harmon. Anal. 21 (1) (2006) 53–94.
doi:10.1016/j.acha.2006.04.004.
- [71]
W. K. Allard, G. Chen, M. Maggioni, Multi-scale geometric methods for data sets II: Geometric multi-resolution analysis, Appl. Comput. Harmon. Anal. 32 (3) (2012) 435–462.
doi:10.1016/j.acha.2011.08.001.
- [72]
S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proc. Natl. Acad. Sci. 113 (15) (2016) 3932–3937.
- [73]
B. Sz.-Nagy, C. Foias, H. Bercovici, L. Kérchy, Harmonic Analysis of Operators on Hilbert Space, 2nd Edition, Springer-Verlag, New York, 2010.
doi:10.1007/978-1-4419-6094-8.
- [74]
C. Firschbacher, S. Naboko, I. Wood, The proper dissipative extensions of a dual pair (2016).
- [75]
K. Schmüdgen, Unbounded and Self-Adjoint Operators on Hilbert Space, Vol. 265 of Graduate Texts in Mathematics, Springer, Dordrecht, 2012.
doi:10.1007/978-94-007-4753-1_2.
- [76]
B. O’Neill, The fundamental equations of a submersion, Michigan Math. J. 13 (4) (1966) 459–469.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] J. M. Ottino, The Kinematics of Mixing: Stretching, Chaos, and Transport, no. 3 in Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, 1989.
- 2[2] R. M. Samelson, S. Wiggins, Lagrangian Transport in Geophysical Jets and Waves: The Dynamical Systems Approach, Vol. 31 of Interdisciplinary Applied Mathematics, Springer, New York, 2006.
- 3[3] G. Haller, Lagrangian coherent structures, Annu. Rev. Fluid Mech. 47 (2015) 137–162. doi:10.1146/annurev-fluid-010313-141322 . · doi ↗
- 4[4] G. Haller, G. Yuan, Lagrangian coherent structures and mixing in two-dimensional turbulence, Phys. D 147 (2000) 352–370. doi:10.1016/S 0167-2789(00)00142-1 . · doi ↗
- 5[5] S. C. Shadden, F. Lekien, J. E. Marsden, Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows, Phys. D 212 (2005) 271–304. doi:10.1016/j.physd.2005.10.007 . · doi ↗
- 6[6] G. Haller, A variational theory of hyperbolic Lagrangian coherent structures, Phys. D 240 (2011) 574–598. doi:10.1016/j.physd.2010.11.010 . · doi ↗
- 7[7] M. Serra, G. Haller, Objective Eulerian coherent structures, Chaos 26 (2016) 053110. doi:10.1063/1.4951720 . · doi ↗
- 8[8] G. Froyland, K. Padberg, M. H. England, A. M. Treguier, Detection of coherent oceanic structures via transfer operators, Phys. Rev. Lett. 98 (22) (2007) 224503. doi:10.1103/Phys Rev Lett.98.224503 . · doi ↗