{'fraction_non_alphanumeric': 0.04359633453869501, 'fraction_numerical': 0.012431199988382055, 'mean_word_length': 4.1977657004830915, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An information-geometric approach to feature extraction and moment reconstruction in dynamical systems

We propose a dimension reduction framework for feature extraction and moment reconstruction in dynamical systems that operates on spaces of probability measures induced by observables of the system rather than directly in the original data space of the observables themselves as in more conventional methods. Our approach is based on the fact that orbits of a dynamical system induce probability measures over the measurable space defined by (partial) observations of the system. We equip the space of these probability measures with a divergence, i.e., a distance between probability distributions, and use this divergence to define a kernel integral operator. The eigenfunctions of this operator create an orthonormal basis of functions that capture different timescales of the dynamical system. One of our main results shows that the evolution of the moments of the dynamics-dependent probability measures can be related to a time-averaging operator on the original dynamical system. Using this result, we show that the moments can be expanded in the eigenfunction basis, thus opening up the avenue for nonparametric forecasting of the moments. If the collection of probability measures is itself a manifold, we can in addition equip the statistical manifold with the Riemannian metric and use techniques from information geometry. We present applications to ergodic dynamical systems on the 2-torus and the Lorenz 63 system, and show on a real-world example that a small number of eigenvectors is sufficient to reconstruct the moments (here the first four moments) of an atmospheric time series, i.e., the realtime multivariate Madden-Julian oscillation index.


Introduction
Extracting temporal patterns from data generated by complex dynamical systems is an important problem in modern science with applications in virtually every scientific and engineering domain. The datasets acquired from such systems are increasingly large both in sample size and dimensionality, and it is worthwhile exploring the multitude of data analysis techniques available in the machine learning literature to analyze them. However, many machine learning techniques consider the data points to be independent and identically distributed, and do not take into account the temporal information (i.e., the time ordering of the data), which is a direct outcome of the dynamical evolution taking place in the system's state space. The feature extraction and moment reconstruction method presented here takes into account this information about the dynamics, therefore placing the current work at the intersection of three different fields, namely machine learning, dynamical systems theory, and information geometry.
In this paper, we take a data-driven approach to the study of dynamical systems. Data-driven methods perform feature extraction by computing the eigenfunctions of a kernel integral operator, i.e., a covariance [1][2][3][4] or heat [5,6] operator, acting on functions defined on the state space (observables). These methods often incorporate information about the dynamics by embedding the data in delay-coordinate spaces (Takens time-lagged embedding) [7][8][9][10]. An alternative approach that has been successfully applied to the analysis of nonlinear dynamical systems is to compute the eigenfunctions of groups or semigroups of operators, e.g., the Koopman or Perron-Frobenius operators [11,12], governing the time evolution of observables under the dynamics [13][14][15][16][17][18][19][20]. Recently, it was established that these two families of techniques can yield equivalent results in an asymptotic limit of infinitely many delays [17,19]. While the above methods are applied directly in the ambient data space, in this paper we propose to work in spaces of probability measures induced by observables of the dynamical system rather than directly in the data space of these observables. These probability measures form together a collection or an ensemble of probability measures that we denote by S. As we will discuss in Sect. 4.3, under certain assumptions the collection S forms a statistical manifold, i.e., a differentiable manifold whose points are probability measures. This allows us to make connections between the framework presented here and the field of information geometry that studies statistical manifolds.
In recent years there has been an increased interest in extracting dynamical patterns by working in probability spaces [21][22][23][24] instead of the ambient data space. In particular, [22] and [23] consider that the observations are drawn from conditional time-varying probability density functions (PDFs). Representative (temporal) patterns on the manifold can be extracted by employing feature extraction techniques, such as the Diffusion Maps algorithm [25] with a kernel based on statistical distances, i.e., divergences, between the PDFs. Another approach [21] is to embed multiple dynamical systems into a single low-dimensional space by comparing their invariant measures. In this approach, phase space trajectories are interpreted as probability measures, and a distance-based embedding method, namely multidimensional scaling [MDS;26], is applied to the Wasserstein distance matrix to perform feature extraction and uncover the low-dimensional manifold. Information-geometric techniques have also been employed in data analysis applications such as flow cytometry or document classification [27,28]. In these cases, the PDFs are estimated over subsets of sample populations, and the distances between (similar enough) PDFs on the statistical manifold are computed using approximations of the Fisher information distance, such as the Kullback-Leibler divergence, the Hellinger distance or the cosine distance. Subsequently, the distances between all pairs of PDFs are approximated by geodesic distances on the statistical manifold, and the full distance matrix is then embedded into a low-dimensional manifold using MDS techniques.
In this paper, we consider to have access to (potentially partial) observations of a deterministic ergodic dynamical system through some scalar or vector-valued observable f . Finite-time trajectories along the dynamical system induce probability measures over the measurable space defined by the observation map f . One of our standing assumptions is that the data is drawn from observations on an ergodic trajectory. We will make our assumptions more precise later, but intuitively, ergodicity ensures that time averages of a quantity on a typical trajectory converge to its space average, as one takes longer and longer averaging windows. Thus an ergodic trajectory fills out the ambient space with a distribution that is the same as that of an invariant measure of the dynamics. We equip the probability space defined by the collection of probability measures with a divergence, here the Hellinger distance, and use it to define a symmetric and positive definite kernel over the space of measures. We demonstrate that the eigenfunctions of a kernel integral operator computed using the Diffusion Maps algorithm [22,23,25] capture temporal and spatiotemporal patterns of interest of the dynamical system, and are meaningful for dimension reduction and moment reconstruction. Under suitable assumptions, the probability measures described above lie on a statistical manifold equipped with a Riemannian metric, namely the Fisher information metric. The distance induced by this metric, i.e., the Fisher information distance, can be approximated by the Hellinger distance, thus further justifying our choice of divergence in the probability space.
An important contribution that we propose in this paper is a novel method to reconstruct the n-th moment of the probability measures induced by trajectories along the dynamical system, for every n ∈ N. The moment reconstruction is done in the Hilbert space defined by the basis of eigenfunctions of a kernel integral operator constructed in the space of probability measures. One of our key results is an identity (Theorem 1) that expresses these moments as time-averaging operators of a certain function on the collection of probability measures. We expand this result into a more general identity (Theorem 2) that applies to any continuous function on the data space, not just moments. We demonstrate the reconstruction of moments  on an atmospheric time series index (the realtime multivariate Madden-Jullian oscilation index), where we show that a small number of leading eigenvectors is sufficient to accurately recover the moments. This opens up the avenue for nonparametric forecasting of the moments of the distributions, and therefore of the distributions themselves.

Contributions and outline of the paper
The main goal of the paper is to present a data-driven technique for the reconstruction of finite-time statistics, e.g., moments, of probability measures defined on trajectories of a dynamical system. To enable that, we lay down a rigorous theoretical framework and prove that the data-driven vectors and matrices that we use in our numerical methods converge to functions and operators on associated underlying spaces. In order to do that, we introduce several spaces and mappings over the course of the paper, and we list the most important ones in Table 1.
Assumption 1 in Sect. 3 introduces the basic assumptions that we impose on the system, and Assumption 3 in Sect. 5 introduces the assumption made on the data available. One of our main contributions is Theorem 1 where we show that the moments of probability measures induced by finite-time trajectories of a dynamical system can be obtained by averaging an associated observation map over the respective trajectories. In Sect. 3 we give a more general version (Theorem 2) of Theorem 1. Theorem 2 is stated for integrals of arbitrary continuous functions with respect to a collection of probability measures. The various moments of a probability measure discussed in Theorem 1 are a special case of such integrals. Section 4 describes several geometric aspects, e.g., embeddings of probability measures and invariance under isometries of the data, using a divergence between probability measures, i.e., the Hellinger distance. If the collection S is a manifold, it can be equipped with a Riemannian metric and the natural divergence associated with it can be approximated by the Hellinger distance, thus allowing us to make connections to the field of information geometry.
The data-driven aspect of our framework and numerical approximations of the continuous quantities introduced in the theorems is addressed in Sect. 5. Algorithm 1 shows the feature extraction in the probability space, and Algorithm 2 shows the moment reconstruction using a small number of leading eigenvectors. Theorem 5 establishes the almost sure convergence of our numerical methods. The proof of all the theorems is done in Sect. 6. In Sect. 7 we study three low-dimensional dynamical systems and a real-world atmospheric time series index, namely the realtime multivariate Madden-Julian oscillation (RMM) index. We end with some concluding remarks and future perspectives in Sect. 8.

Dynamics-dependent probability measures on the data space
The following will be a standing assumption in the rest of our discussions. Assumption 1.X is a C 1 -manifold (differentiable manifold) with a C 1 deterministic flow Ψ t :X →X. There exists a Ψ t -invariant ergodic, probability measure µ for the flow, with a compact support X ⊂X. X is equipped with its Borel σ-algebra B(X). f :X → Y is a C 1 observation map taking values on a manifold Y .
In the examples from Sects. 7.1 and 7.2 ahead, X =X = T 2 , the 2-torus, whereas in Sect. 7.3,X = R 3 and X is the Lorenz 63 strange attractor. The space Y plays the role of a data space, and is often the Euclidean space R d . The triple (X, Ψ t , µ) defines a measure-preserving dynamical system. One can associate to each state x ∈ X a trajectory {Ψ t (x)} t∈R . Our focus will be on the collection of probability measures induced by the observation map f on finite-time trajectories at each point, and their associated statistics, such as the moments of these probability measures.

Probability measures on the data space
Let T = [−∆t, 0], ∆t > 0, be a closed time interval, B(T ) be its Borel σ-algebra, and λ be the Lebesgue probability measure on T . Next define a map g which assigns to every state x ∈ X a continuous, Y -valued map g x = g(x) on T . g x is defined as g x (t) = f (Ψ t (x)) for every t ∈ T . Thus g is a mapping from X to the set of continuous functions defined on T with values in Y , g : X → C 0 (T ; Y ). Let Prob(Y ) denote the set of Borel probability measures on Y , and Prob c (Y ) denote its subset of compactly supported measures.
Here h * λ is the push-forward of the Lebesgue probability measure λ on T to a Borel probability measure on Y , defined as h * λ(A) = λ(h −1 (A)) for every Borel set A ⊂ Y . Now define the map p : X → Prob c (Y ) as The image of the observation map f will be denoted as M ⊂ Y , and the image of p, namely the set p(X) = {p x | x ∈ X}, will be denoted as S. Thus S ⊂ Prob(M) ⊂ Prob c (Y ) ⊂ Prob(Y ). Our main focus in this paper will be this collection S of probability measures p x . They represent the set of all possible probability measures obtained from trajectories of length ∆t from initial states in X.
In a typical data-driven setting, X and its dynamics are unknown, and we only observe data points y = f (x), y ∈ Y , as functions of the unobserved states x. In the operator-theoretic framework, e.g., Koopman or Perron-Frobenius operators, one studies the effect of the dynamics on the space of observables on X, i.e., the data space Y . In the framework presented here, we go one step further from a typical data-driven setting and consider observables γ on the data space Y . These indirectly lead to observables on X through composition with f . Namely, any function γ : Y → R m induces the function γ • f : X → R m . Here we consider the set of probability measures Prob(M) or Prob c (Y ). As we will see in the following, the moments of the probability measures p x are a special case of the observables γ, and are of particular interest in this paper.

Moments of the probability measures in S
Consider the case when Y = R d . Every p x in the collection S is a probability measure with a compact (bounded) support in the set M ⊂ Y = R d . Therefore all of its moments exist and are finite. Let γ (n) : R d → R d be the function that raises each component of a vector to the power n, i.e., Then for every p x ∈ S and n ∈ N, the d-dimensional vector of the n-th moment E n of p x is defined as Therefore, E n (p x ) is the d-dimensional vector whose i-th coordinate E (i) n (p x ) = Y (y i ) n dp x (y)dy. Our first main result (Theorem 1 below) expresses these moments as a function of the initial state x ∈ X, in terms of the following time-averaging operator A ∆t : Theorem 1. Let Assumption 1 hold and Y = R d . Then for every n ∈ N and x ∈ X, the n-th moment E n (p x ) is finite and Taking for example n = 1, we have γ (n) = γ (1) = I Y , the identity map on Y , and the first moment (mean) E 1 is a time-averaging operator acting on the observable f since A ∆t (γ (1) Remark. A (probability) measure is completely characterized by its moments, and Theorem 1 expresses these moments in terms of a time-averaging operator. Note that in a data-driven setting, the observation map f is fixed and its values are known at some sampling points x i . The function γ (n) • f can be computed to any desired degree of accuracy at these sampled points. The operator A ∆t involves an integral which can be numerically approximated by averaging along a sampled trajectory. We later show in (4) that Theorem 1 is a special case of an identity that involves more general functions γ : Y → R m , e.g., Fourier functions, spherical harmonics or polynomials, not just γ = γ (n) from (2).

Observables on the data space
Since for every function γ ∈ C 0 (M; R m ), J(γ) is a linear map on Mes(M), it automatically becomes a continuous map on S, which is a collection of probability measures contained in Mes(M). Moreover, for every a 1 , a 2 ∈ C and every γ 1 , γ 2 ∈ C 0 (M; R m ), making J a linear map. Another way to view the action of J is as a dual mapping. The space Mes(M) lies in the dual space to the Banach space C 0 (M; C), and C 0 (M; C) embeds isomorphically and canonically into the double dual C 0 (M; C) * * = Mes(M) * . Then it follows from (3) that, by definition, J maps each component of γ into its double dual. The following important theorem expresses the action of J(γ) on elements of S, p x ∈ S, in terms of the time-averaging operator A ∆t , for every ∆t > 0.
Theorem 2. For every γ ∈ C 0 (M; R m ) and every x ∈ X, The diagram below illustrates the domains and codomains of the various maps defined so far, and how they are connected through (4).
We will next discuss various geometric aspects of the collection S of probability measures.

Geometrical structures on the collection of probability measures
In the following we define a metric on the space Prob c (Y ) of compactly supported Borel probability measures on Y in a way such that the map p : X → Prob c (Y ) is continuous. One initial difficulty in defining such a metric is that for different x, x ′ ∈ X, the measures p x , p x ′ are highly singular measures with respect to (wrt) the Lebesgue measure on Y . Moreover, their supports are one-dimensional non-intersecting curves. Thus the task is to define a concept of distance between measures supported on almost disjoint curves lying in high-dimensional ambient data spaces. This can be achieved by smoothing out these singular measures using a C r kernel function k : Y × Y → R + 0 , by a procedure called kernel density estimation [29]. Let Mes c (Y ) denote the space of finite signed measures on Y with compact support, and C r (Y ) denote the set of functions r times differentiable on Y . Then one has the following map It is easy to check that when restricted to Prob c (Y ), we have a map K : Prob c (Y ) → C r (Y ; R + 0 ) . Now let α be any reference measure on Y . Given a C 0 (continuous) non-negative function h : Y → R + 0 , hα is a new measure which is absolutely continuous wrt α and with a C r density h. Let Prob(Y ; α) denote the set of Borel measures on Y which are α-a.c. (i.e., absolutely continuous wrt α), and with C 0 density functions. Thus we have the canonical embedding We summarize the various maps we have described in the commutative diagram 1 below.
The dashed arrow defines the map which embeds the state (phase) space X of the dynamics into a space of measures which are absolutely continuous wrt α.
Choice of kernel. The kernel k is usually taken to be an isotropic kernel, i.e., of the form where η : R → R + 0 is the so-called shape function, and d(y, y ′ ) is the distance between data points y and y ′ , e.g., the Euclidean distance. We are going to need the following assumption on the kernel : 0 is a C r , strictly positive definite, isotropic kernel as in (6). Further, its shape function η satisfies lim |d|→∞ η(d) = 0, with d the distance measure.

The Hellinger distance
In order to study the geometrical properties of the collection of probability measures S we need to equip the probability space with a divergence, i.e., a notion of dissimilarity between probability measures. We choose the Hellinger distance over other divergences (e.g., Kullback Leibler divergence (relative entropy), Wasserstein distance, total variation distance) due to several reasons that will become clearer in Sects. 4.1-4.3.
Given two probability distributions π and π ′ which are absolutely continuous wrt the probability measure α with densities ρ, ρ ′ , the squared Hellinger distance is defined as This definition is usually accompanied by a multiplicative factor of 1/2, which we have dropped. We will now establish some conditions under which p α is an embedding into Prob c (Y ).
Theorem 3. Let Assumption 1 hold, then for every x ∈ X, p x has a compact support contained in M, and p : X → Prob(M) is a continuous map in the weak topology on Prob c (Y ) (and Prob c (M)). Further, let Assumption 2 hold. Then p α : X → Prob(Y ; α) is a continuous map wrt the Hellinger distance. If in addition, the observation map f : X → Y is one-to-one, then p is a homeomorphism between X and S, and p α is injective.
Remark. If X is a manifold and p is injective, then S can be assigned the same manifold structure as X, and S becomes a statistical manifold of probability measures, parameterized by the manifold X. With this manifold property, techniques from the field of information geometry [30] (Sect. 4.3) can be employed to design data analysis algorithms taking advantage of the geometrical structure of S. In Sect. 7.3, we will demonstrate with numerical experiments on the Lorenz 63 system that these techniques remain useful even if the state space X and/or S are not smooth manifolds. Statistical manifolds carry a lot of information intrinsic to the underlying dynamical system and are easily tractable from a data point of view. Note that the maps g, and therefore p and p α , depend on the observation map f : X → Y . We would like to have invariance of the collection of probability distributions under transformations of the dataset M = f (X) such as translation and rotation. We show in the following that this is achieved by our kernelbased embedding of Prob c (Y ) into the space Prob(Y ; α) of α-a.c. measures.
Diffeomorphisms of the image. Let D : Y → Y ′ be a C r map which maps Y diffeomorphically into its image.
The following commutative diagrams succinctly display the relations between these two sets of maps.
The map D * : S → S ′ in the diagram is the push-forward of probability measures on M under the map D, and the map D• : Theorem 4. Let Assumption 1 hold, and D : Y → Y ′ be the diffeomorphism as above. Suppose that the kernel functions k : Then we have the following commutative diagram.
is an isomorphism wrt the Hellinger distances of the respective spaces.
The theorem is proved in Section 6.3.
Remark. If k is an isotropic kernel (6), then Theorem 4 says that the pull-back metric on X, given by The measure α. Note that we did not place any restriction on the choice of the measure α. A natural choice for α is the push-forward of the invariant measure of the dynamics µ on X, i.e., α = f * µ. Thereby, we can take advantage of the ergodicity of µ and approximate α by the samples y i = f (x i ) on the data space. We use here for α the Lebesgue measure which is consistent with our definition of probability measures in (1), but the above formula does not depend on the type of measure used as reference, i.e., the Hellinger distance will not change if the densities are defined relative to a different equivalent measure. Moreover, we have the following commutative diagram, similar to Theorem 4.

Kernel integral operators on S
Given a choice of a reference measure α on Y , one gets the distance d α on the state space X (8) induced by the Hellinger distance. Using d α we will define a new kernel on Prob c (Y ) (or S), defined similarly to the isotropic kernel (6), with the shape function η chosen to be the Gaussian function where ǫ > 0 is a bandwidth parameter. The resulting Gaussian kernel on Prob c (Y ) is given by In the particular instance when both π, π ′ lie on S and equals p x , p x ′ respectively, Since p maps X into S, it induces a probability measure ν on S, which is the push-forward of the Ψ tinvariant probability measure µ on X, i.e., ν = p * µ . Thus for any Borel measurable set U ⊆ S, ν(U ) = µ{x ∈ X | p x ∈ U }. As a result, one can define a kernel integral operator G : where the second equality above follows from the change of variables formula for integrals and the definition of ν. Each function φ in the range of G, also a member of L 2 (ν), is a pointwise defined function continuous wrt the weak topology on S. Let 1 S denote the constant function equal to 1 on S. Next, we define a sequence of normalizations as in the Diffusion Maps algorithm [25] : to get a Markov integral operator H : L 2 (ν) → L 2 (ν). The function v is a continuous function on S and is positive everywhere, therefore it can be interpreted as a density function. The operator H is not symmetric, but if V and Q denote the multiplication operators by the functions v and q, respectively, theñ is a symmetric operator. Moreover, since v is a continuous function, positive and uniformly bounded away from 0, V and V −1 are both bounded operators. Thus,H is also compact and self-adjoint, and it has a complete eigenbasis consisting of functions {ψ l } ∞ l=1 with eigenvalues 1 − λ l : The eigenfunctions ψ l are orthonormal in L 2 (ν), i.e., wrt the inner product and the measure ν: Since H andH are related by a similarity transformation, H has the same spectrum asH. Thus H has the same eigenvalues asH with the corresponding eigenfunctions φ l satisfying Let ω be the measure whose density wrt ν is v, i.e., dω(π) = v(π)dν(π). Then the φ l s form a basis, which is not orthonormal in L 2 (ν) but orthonormal in L 2 (ω), namely, The eigenfunctions φ l of the Markov matrix H are also the eigenfunctions of the random walk Laplacian operator L = I − H, i.e., with eigenvalues λ l . Similarly, the eigenfunctions ψ l ofH are the eigenfunctions of the symmetric normalized Laplacian operatorL = I −H, where I is the identity operator. The normalized LaplacianL is positive semidefinite and therefore has real non-negative eigenvalues 0 = λ 1 < λ 2 ≤ . . . ≤ 2.
The eigenfunctions φ l correspond to temporal patterns of the dynamical flow and capture different timescales of the system induced by the probability measures in (1) (as will be illustrated in Sect. 7) acting as filters on L 2 (ν). Moreover, since the function J(γ) with γ ∈ C 0 (M; R m ) from Theorem 2 is in L 2 (ν), there exist expansion coefficients c l ∈ R m such that we can expand J(γ) in the {φ l } eigenfunction basis: Using (14), (2) and J(γ) = E n , the moments of the probability distributions p x can be written as where c l ∈ R d are the expansion coefficients used for the moment reconstruction.
In Sect. 5 we will consider the data-driven implementation of these quantities and of Theorems 1 and 2. We end this section by considering a special case when S has an additional geometric structure.

Statistical manifolds
A case of particular interest is when the collection S of probability measures is itself a manifold, thus justifying the name of statistical manifold, i.e., a manifold where each point is a probability measure. Statistical manifolds and their properties are studied in the field of information geometry [30,31] using techniques from differential geometry. If the probability measures are all absolutely continuous with respect to a reference measure, say τ , then S can be locally parameterized with a coordinate system θ = {θ 1 , . . . , θ W } inducing a parametric family of probability densities ρ(·, θ) on R d with respect to τ . To track the statistical distances, i.e., divergences, between probability distributions locally near a point p ∈ S, one needs to define a Riemannian metric g p . There is a canonical Riemannian metric, called the Fisher information metric, which measure the amount of information between two PDFs. It is an inner product on the tangent spaces of S defined by the natural basis of tangent vectors ∂ ∂θ 1 , . . . , ∂ ∂θ W , given by The Fisher information metric as defined above is positive definite and transforms as a type (0, 2) tensor on S under changes of coordinate system θ. Given any two tangent vectors u = , is called the Fisher information matrix (FIM) and is a positive definite matrix. The natural divergence on statistical manifolds associated with FIM is the Fisher information distance. For infinitesimal small changes in the probability distributions (i.e., θ x ′ = θ x + dθ, for p x , p x ′ ∈ S), the Fisher information distance can be expressed using the quadratic differential form, i.e., Approximations of the parametric Fisher information distance can be achieved by well-known nonparametric divergences in the literature, such as the Kullback-Leibler (KL) divergence, the Hellinger distance or the cosine distance [32]. From the family of f -divergences, the Hellinger distance satisfies all metric properties, and it has been shown to outperform the symmetric KL divergence (Jeffreys divergence) and Bhattacharyya distance (at least in the Gaussian case) due to the fact that the latter divergences do not obey the triangle inequality [33]. Its metric properties also allow for isometric embeddings in reproducible kernel Hilbert spaces [33]. These properties of the Hellinger distance provide an additional argument for our choice of divergence in the probability space, and the Hellinger distance remain a good choice of a metric as an approximation of the Fisher information distance even for the case when S has the additional manifold structure.

Numerical approximations
Though we introduced our framework in the continuous case, in practice we work with discrete versions of the continuous-time dynamical systems described above. For the discrete case, we make the following assumption.
Assumption 3. There is a sequence of states x i = Ψ i δt (x 0 ) for some initial state x 0 ∈ X and sampling interval δt, and a time-ordered sequence of measurements y i = f (x i ). The sampling interval δt is such that the discrete-time dynamical system (X, Ψ δt , µ) is also ergodic.
The assumption on ergodicity of µ will be key to ensuring that our results of data analysis converge in the limit of large data (N → ∞ and δt → 0). In most practical situations, the invariant measure µ is unknown to us. Instead, given a dataset consisting of N samples, µ is approximated using the discrete sampling measure the average of Dirac-δ measures supported on a trajectory {x 0 , . . . , x N −1 }. By the ergodicity of µ (Assumption 1), for µ-a.e. (almost everywhere) x 0 ∈ X, the discrete sampling measures µ N converge weakly to µ. This means that for a set of initial points x 0 with µ measure equal to one, Ergodicity is a property of the system that is often implicitly assumed in data-driven studies of dynamical systems. It provides a justification for the principle that the global statistical properties of an observable f : X → R d wrt µ can be obtained from a time series for f , namely, f (x 0 ), . . . , f (x N −1 ). The second discretization required is for the measures p x . These are the push-forward measures under g x of the Lebesgue probability measure λ on T (1). The measure λ will now be discretized by the measure λ R = 1 R R−1 r=0 δ −rδt on T , where ∆t = Rδt, and R is the number of samples within an embedding window. The interval T = [−∆t, 0] has been discretized as {−rδt : 0 ≤ r < R}. As a result we have, . Just as with the measure ν, p andp also push the sampling measure µ N into the discrete measures ν N = p * µ N andν N =p * µ N , respectively. ν N andν N are measures on the collections S andŜ respectively, whereŜ :=p(X) = {p x | x ∈ X}. Note that both S andŜ are contained in Prob c (Y ). In the data-driven approximation scheme, the space L 2 (µ) will be approximated by L 2 (µ N ), and L 2 (ν) by L 2 (ν N ) and L 2 (ν N ). If the set X is not a fixed point, then for µ-a.e. x 0 , these three spaces are isomorphic to C N , and the functions in these spaces are N -dimensional vectors. The inner product on L 2 (ν N ) is given by where the last inner product is the usual inner product on C N . L 2 (µ N ) and L 2 (ν N ) carry analogous inner products. Under the kernel density estimation,p x =p(x) is mapped intô Finally, we choose the reference measure α on Y to be a discrete measure supported on Q points {z q : q = 0, . . . , Q − 1}. The squared Hellinger distance becomes We will overuse notation and use d 2 H to also denote the distance induced on the setŜ :=p(X) by the squared Hellinger distance in (20), namely, Remark. A typical choice for α is any set of Q points Y with an independent and identical distribution wrt some measureᾱ. As Q → ∞, α would converge weakly toᾱ. Another choice is α = f * µ Q . This has the added advantage that it requires no extra points of evaluation other than on the original data points. This relates our work to kernel mean embedding techniques [34,35] that have been used to embed probability measures into a Hilbert spaces of functions, called Reproducing Kernel Hilbert Spaces (RKHS). This Hilbert structure provides additional tools, such as orthogonal projections, and allows to use pointwise evaluations as bounded linear functionals. We however do not use the RKHS aspect of the technique.
The infinite-dimensional space L 2 (ν) will be approximated by the N -dimensional space L 2 (ν N ), the operators defined in Sect. 4.1 become N × N matrices. We will denote the matrix versions of all the operators using boldface notation. The kernel integral operator G will be approximated by a N × N matrix G = {G ij } acting on L 2 (ν N ), as We have similar to the continuous case, Let V be the diagonal matrix V := diag (v 0 , . . . , v N −1 ). Then note that H is a Markov matrix of transition probabilities and is similar to the symmetric normalized matrixH below: SinceH is a normalized symmetric matrix, it has a complete eigenbasis consisting of orthonormal vectors with eigenvalues 1 − λ l . The eigenvectors ψ l are orthonormal with respect to the standard inner product ψ k , ψ l νN = 1 N N −1 i=0 ψ k,i ψ l,i = δ kl . Moreover, since H andH are related by a similarity transformation, H has the same spectrum asH, and we have The asymmetric random walk Laplacian matrix associated with the random walk matrix H is L = I − H = I − V −1G , and it has the same eigenvectors φ l as H: L has the same eigenvalues λ l as the symmetric normalized Laplacian matrixL = I −H which is positive semidefinite and therefore has non-negative eigenvalues 0 = λ 1 < λ 2 ≤ . . . ≤ λ N ≤ 2.L has the same eigenvectors asH,Lψ l = λ l ψ l . The Laplacian eigenvectors associated with the lowest eigenvalues λ l vary slowly on the graph, i.e., two vertices that are connected by an edge with a large weight will have similar values of the leading eigenvectors [36].
The top eigenvector ψ 1 ofH coincides with v 1 2 and corresponds to the eigenvalue 1 − λ 1 = 1. This follows from the fact that since H is a Markov operator, the eigenvector corresponding to eigenvalue 1 is the constant vector φ 1 = 1. The eigenvectors φ l are orthonormal with respect to the L 2 (ω N ) inner product, where the measure ω N is absolutely continuous wrt the sampling measureν N and has density v = ψ 2 1 : We have outlined a summary of the entire numerical procedure to obtain the λ l and φ l in Algorithm 1.

Algorithm 1: Feature extraction (temporal patterns)
Input: ∆t -embedding window R -number of delay-coordinates, so that the sampling interval is δt = ∆t/R

Moment reconstruction
Using Theorems 1 and 2 and (14), we can expand the moments of the distributions in the {φ l } eigenvector basis. The moments can be reconstructed (Algorithm 2) using only the leading M eigenvectors and the expansion coefficients c l ∈ R d as follows:Ê where v i = v(p xi ), and

Algorithm 2: Moments reconstruction
Input: Φ = {φ l } -eigenvectorŝ p xi -the estimated probability distributions for the embedding windows v -the density vector Output:Ê n (p xi ) -reconstructed moments 1: Compute the true (first few) moments of each probability measure E n (p xi ).
3: Reconstruct the momentsÊ n (p xi ) from (26) using only the leading M eigenvectors φ l .

Continuous extensions
In the following, we will need to make more clear the dependence of the data-driven quantities used previously in this section on the parameters N and R, and we will add these parameters as subscripts. For example, the vectors v and φ l from (22) and (24) will be denoted as v N,R and φ N,R,l , respectively. The dependence on R comes from the fact that these quantities are constructed using the mapp from (17) which depends on R.
The eigenvectors φ l from (24), denoted as φ N,R,l , are N -dimensional vectors and are interpreted as functions on the set {p xi : i = 0, . . . , N − 1}. These functions can be continuously extended to the space Prob c (Y ), which contains both S andŜ. This is done in a manner similar to (22). We use a boldface notation for vectors, and have dropped the boldface for their continuous extensions. In the equations below π, π ′ ∈ Prob c (Y ): As a result of these definitions, there exists a continuous function φ N,R,l : Prob c (Y ) → R such that the i-th coordinate of the vector φ N,R,l is given by We now have the following theorem establishing the convergence of our data-driven methods.
Theorem 5. Let Assumptions 1 and 3 hold, and v N,R and φ N,R,l be the data-driven functions as above. Then for every γ ∈ C 0 (M; R m ) and µ-a.e. x ∈ X, we have Remark. The moments E n of the distributions p x can be expanded in the {φ l } eigenvector basis. If we consider only the subspace span{φ 1 , . . . , φ M } spanned by the first leading M eigenvectors, then we can approximate the moments using only this reduced basis. Results on a real-world time series show that a small number of eigenvectors is already sufficient to accurately approximate the moments (see Sect. 7).
Out-of-sample extension methods such as the kernel analog forecasting method [37] can be used to forecast future values of the eigenfunctions, and therefore also to predict the moments of the distributions p x .
The role of dimensionality. While density estimation in high dimensions is known to pose numerical problems because of the curse of dimensionality (i.e., the number of samples needs to grow exponentially with the number of dimensions [38]), nonparametric multivariate density estimation such as kernel estimators have proven to perform well in low dimensions (d ≤ 3). If the sampling density δt → 0, nonparametric density estimation techniques will perform well even in higher dimensions, i.e., for higher-dimensional observables.
In the numerical experiments presented in this paper we mostly consider the observable f to be one-or two-dimensional, f : X → R or f : X → R 2 , and we will use 1D or 2D kernel density estimation techniques.
Strength and limitations. We proposed this framework for ergodic deterministic dynamical systems, however our experiments and additional initial results indicate that in practice the approach is also applicable to non-stationary and stochastic systems (such as RMM (Sect. 7.4)) with promising results. The probability measures capture temporal information on the dynamics, thus going beyond standard analysis in ambient data spaces, however part of this dynamical information is lost within the embedding windows themselves.
Comparison with previous works. Our framework is related to the works of [22] and [23]. In [22], the timevarying PDFs are approximated by histograms, and the pairwise distances between them are computed using the Mahalanobis distance. Computation of these distances requires the estimation of time-window local covariance matrices between histograms, which can prove to be costly, especially if the features (histograms) are high-dimensional, i.e., a high number of histogram bins. When taking the features to be equal to the logarithm of the histograms, the local covariance matrix becomes related to the Fisher information matrix, thus revealing an implicit connection to the field of information geometry. The Fisher information matrix and the associated statistical manifold are explicitly employed in [23], where a parametric model of the underlying PDFs is assumed. More precisely, the observations at each time step are assumed to be drawn from a multivariate Gaussian distribution with time-evolving parameters. The Kullback-Leibler divergence is used to compute the distances between the PDFs on the statistical manifold, and to define a Gaussian kernel used further for dimension reduction with Diffusion Maps. In [22] and [23], the dynamical systems are stochastic, i.e., dynamical noise is present, and the observations are assumed to be drawn from a time-varying probability density function. Empirical PDFs over time windows are subsequently introduced for estimation purposes. We work with deterministic dynamical systems where the probability measures are induced by trajectories of the dynamical system over specified time windows. Thus, even though methodologically our approach has certain aspects in common with the methods of [22] and [23], there are differences in perspective since we consider how a deterministic dynamical system acts on probability densities of observables over time windows, rather than using such windows for estimation purposes in a stochastic setting. Our framework does allow for observational noise which is different from stochastic/dynamical noise. Another difference between their approach and ours is the way we use the empirical histograms to compute pairwise distances between states of the dynamical system. We use nonparametric multivariate density estimation to compute joint densities between the vector components, i.e., the vector-valued observables f , and assign pairwise distances between dynamical states using the Hellinger distance between those densities. On the other hand, [22] assigns pairwise distances through the Mahalanobis distances of concatenated histograms of each vector component, while [23] assumes a parametric model for the generating PDFs, i.e., a multivariate Gaussian model with time-evolving local covariance matrices.

Proofs of Theorems 1 and 2
The following basic result from measure theory will be needed. It is commonly referred to as the "change of variable" formula for integrals. Lemma 6. Let (A, ξ) be a measure space, and T : A → B a measurable one-to-one map. Let χ ∈ L 1 (B, T * ξ) and T * ξ be the push-forward of the measure ξ under T (see the diagram below).
This lemma is a basic result from Analysis [e.g., see 39, Thm 1.6.9] and we will skip the proof. We next prove Theorem 2.
Proof of Theorem 2. Note that dλ = dt ∆t and the right hand side (RHS) of (4) can be rewritten as

Now invoke Lemma 6 with the substitutions
This completes the proof of the theorem.
Proof of Theorem 1. The n-th moment E n of p x can be rewritten in terms of the linear transformation J(γ) as E n (p x ) = Y γ (n) (y)dp x (y) = J(γ (n) )(p x ).
It follows from this equation and by substituting γ = γ (n) in (4) that, This completes the proof of the theorem.

Proof of Theorem 3
Since M is the image of f , for every x ∈ X, the image of g x is contained in M and therefore p x , which is the push-forward of λ under g x , will have a support contained in M. This proves the first part of the claim. To prove the continuity of p, we have to show that for every continuous map ζ : Y → R , the map x → Y ζ dp x is continuous. So let ǫ > 0 and x ∈ X be fixed. ζ restricted to the compact set M is uniformly continuous and therefore there is a δ > 0 such that for z, z ′ ∈ M, if d(z, z ′ ) < δ, then |ζ(z) − ζ(z ′ )| < ǫ. Now since X is a compact set, and Ψ t , f are continuous maps, for x ′ sufficiently close to x, we will have Thus for every such x ′ close to x, by the choice of δ, Invoking Lemma 6, we get for every fixed ζ ∈ C 0 (Y ) and every y close to x, This proves the continuity of p as claimed. Note that if f is injective, the support of each p x , which is the curve {f (ψ t (x)) : t ∈ T } are distinct, and hence the p x are necessarily distinct. If k is strictly positive definite, then K is injective by a result of Lemma 7 (Mercer's theorem [41]). Let β be a Borel measure on a first-countable, topological space, with compact support Z. Let k : Z × Z → R be a continuous, symmetric, strictly positive definite kernel on Z. Then there exists an orthonormal eigenbasis {e i : i ∈ N} of L 2 (Z, β) of continuous functions such that Moreover, the convergence is absolute and uniform on Z.
Now let x ∈ X and x n be a sequence of points converging to x. We will show that for every ǫ > 0, Let δ 1 > 0 be fixed. Then for x ′ sufficiently close to x, the Hausdorff distance of the images of the functions g x and g x ′ is less than δ 1 . Let Z ′ be the closed δ 1 -neighborhood of the image of g x . Thus for x ′ close to x, the support of p(x ′ ) will lie in Z ′ . Further, by Assumption 2, there is a δ 2 > 0 such that if d(y, y ′ ) > δ 2 , then k(y, y ′ ) < ǫ. Let Z be the closed δ 2 neighborhood of Z ′ . Note that for every probability measure p with support in Z ′ and every y ∈ Y \ Z, Therefore, since the supports of the p(x n ) eventually lie in Z ′ , we have lim sup n→∞ sup y∈Y \Z |K(p(x n )) − K(p(x))(y)| < 2ǫ.
Thus it only remains to prove an inequality analogous to the one above, but for y ∈ Z. By Mercer's theorem, there are continuous functions {e i : i ∈ N} and λ i > 0 such that k(y, y ′ ) = i∈N e i (y)e i (y ′ ) uniformly and absolutely over Z. Thus there exists an L ∈ N such that Therefore for any probability measure q with support in Z ′ , we have Now replacing q by p(x) and p(x n ) respectively, gives sup y∈Z |K(p(x n )) − K(p(x))(y)| < L i=1 e i (y) e i (y ′ )dp x (y ′ ) − e i (y ′ )dp xn (y ′ ) + 2ǫ.
where we have used the fact that p(x n ) converges weakly to p(x). This completes the proof of continuity of ι r • K • p. Its injectivity follows from the injectivity of all these three maps.
We will now show that ι α : is a continuous map wrt the Hellinger distance. So let ρ, ρ ′ ∈ C 0 c (Y, R + 0 ). Since ρ has compact support, it has a bounded range, [a, b]. If ρ − ρ ′ C 0 (Y ) < ǫ, then by the continuity of the square-root function on the interval [a, b], Since ǫ was arbitrary, this proves the continuity of ι α at each point ρ ∈ C 0 c (Y, R + 0 ). The proof of injectivity of ι α is trivial and is left to the reader. This completes the proof of Theorem 3.

Proof of Theorem 4
To prove that a commutative diagram is true, it is necessary and sufficient to prove that the commuting relations in the smallest loops hold. We will now break down the commuting diagram in Theorem 4 into its smallest commutative components. First consider the relations The first two relations follow from the definitions of f ′ and g ′ . The next two follow from the definition in (1). Next consider the relations The first relation is already explained in Sect. 4.1. To prove the second relation, take some π ∈ Prob c (Y ). Then, for every y ∈ Y Therefore, (K ′ D * π)•D = Kπ. Finally, to check the last relation, take a ̺ ∈ C r (Y ′ ; R + 0 ) and let β := ι α (̺•D). To show that the two measures D * β and ι D * α (̺) are equal, it is enough to show that their integrals with every continuous function with compact support χ ∈ C 0 c (Y ′ ) are equal. This holds because This completes the proof of Theorem 4.

Proof of Theorem 5
As discussed in Subsect. 5.2, we will use R and N in the subscripts to denote the dependence on the parameters R and N . Similarly to v N,R and q N,R , we denote the data-driven matrices G from (21) and H from (24) as G (N,R) and H (N,R) , respectively, since they are constructed using the mapp which depends on R as in (17). If they are constructed in an analogous manner using the map p, the resulting matrices will be denoted as G (N ) and H (N ) . For the former pair of matrices, the Hilbert space in question is L 2 (ν N ), while in the latter case it is L 2 (ν N ). One similarly gets vectors v N,l and φ N,l and their continuous extensions v N,l and φ N,l , instead of v N,R,l , φ N,R,l , v N,R,l and φ N,R,l . The expansion coefficients c N,R,l will be denoted as c N,l if we use the map p instead ofp.
Given a function/observable γ ∈ C 0 (M; R m ), since J(γ) ∈ L (Mes(M ); R m ), J(γ) is an m-dimensional vector of continuous functions on S. For simplicity and without loss of generality, in the proof we make the assumption that m = 1. Using the φ l basis from (12), J(γ) can be represented as an L 2 (ν) function: Here, the inner product is taken between ψ l = v 1 2 φ l and each component of J(γ) wrt the measure ν. Based on results from [42,43], J(γ) can be approximated from data in the L 2 (ν) norm. The following limit holds in the L 2 (ν) sense for µ-a.e. x: As mentioned above, the functions v N and φ N,l are continuous extensions of the vectors v N and φ N,l when using the map p instead ofp. From the theory of graph Laplacians [e.g., 44], we have Since the measures λ R converge weakly to λ, their push-forwards under g x also converge for every x ∈ X, namely,p where w − − → denotes weak convergence of measures. Thus, for every 1 ≤ l ≤ N , Next, since X is compact, the image M = f (X) is compact too, and thus by (31), for every x ∈ X. As a result, for every 1 ≤ i, j ≤ N , Thus, the N × N matrix G (N,R) converges to G (N ) as R → ∞. It can be similarly shown that H (N,R) converges to H (N ) . Therefore, by the theory of perturbation of compact operators [45,Sect. 7], their eigenvectors and eigenvalues also converge. In particular, the functions and vectors resulting from the continuous extensions in (28) converge uniformly, too. In summary, we have The first consequence of the uniform convergence in (33) is that The parameter R can now be inserted in (29) using (30), (31) and (34)  This completes the proof of Theorem 5.

Experiments
In this section, we illustrate our feature extraction and moment reconstruction framework on applications to three low-dimensional dynamical systems on: 1) an integrable ergodic flow on the 2-torus with two different dynamical regimes, Model I ( Fig. 1(a)) and Model II ( Fig. 1(b)), 2) the Oxtoby system on the 2-torus [46] with a fixed point (Fig. 1(c)), and 3) the Lorenz 63 system (Fig. 1(d)). Each dynamical system is observed only through partial observations in the sense that the map f is not one-to-one. Some of these systems, and their associated eigenvectors, have been previously analyzed using different kernels and operators -Models I and II on the 2-torus using the cone kernel [47], and the Oxtoby system with a fixed point, Models I and II using the Koopman operator [17]. These methods relied on a full observation map.

Integrable flow on the 2-torus
In our first example, we consider an ergodic dynamical system whose state space X is the 2-torus. Denoting the azimuthal and polar angles by ϑ 1 and ϑ 2 , respectively, the vector field of this system is given where ζ is an irrational angular frequency parameter. The speed variations of the flow are controlled through the parameter β ∈ (0, 1], such that the flow along the manifold speeds up or slows down when β < 1 (here β = 0.5). We consider two dynamical models on the 2-torus with two different angular frequency, ζ = 30 1/2 and ζ = 30 −1/2 , denoted Models I and II, respectively. The datasets for Models I and II were generated using N = 64, 000 samples at a sampling interval δt = 2π/S, where S = 500 controls the number of samples in each quasi-period (equivalent to approximately 128 periods). The full observation map is the standard embedding of the 2-torus in R 3 , that is, where both the azimuthal radius r 1 and polar radius r 2 are equal to 0.5. All kernel density estimators use Q = 50 evaluation points per dimension. In both models, we varied the parameter values and the results seem robust within certain ranges. See Appendix A for more details.

Model I
The angular frequency of Model I was set to ζ = 30 1/2 . Here, the partial observation map f is given by f = (f 1 , f 2 ). Figure 2 shows two examples of PDFs over trajectories of multivariate (2D) observables for the 2-torus (35), Model I. The eigenvectors are computed using (24), and some of the most representative ones are shown in Fig. 3. The embedding window in (18) is set to R = 40 timesteps, i.e., ∆t = Rδt = 40 * 2π/500. The parameters of the Diffusion Maps algorithm were set to k = 500 nearest neighbors, the width of the Gaussian kerne ǫ = 1, and the normalization parameter α = 1. The normalization parameter α in the Diffusion Maps algorithm will be henceforth set to α = 1, to fully decouple the geometry of the data from the density of the data [25]. The eigenvectors φ l in Fig. 3 capture different timescales of the dynamical flow. The probability measures introduced in (1) uncover temporal patterns of the underlying dynamical patterns despite having access to only partial observations of the system through f instead of F from (36). This is due to the fact that information on the dynamical evolution of the system is implicitly captured in the individual trajectories. The regions where the system evolves slowly (negative values along the horizontal axis in Fig. 3 corresponding to negative values of the azimuthal angle ϑ 1 ) correspond to regions where the eigenvectors vary to a high degree.

Model II
The angular frequency of Model II was set to ζ = 30 −1/2 . Some of the most representative eigenvectors are shown in Fig. 4. The embedding window in (18) was set to R = 80 timesteps. The parameters of the Diffusion Maps algorithm were set to k = 7, 000 nearest neighbors, the width of the Gaussian kernel ǫ = 0.18, and the normalization parameter α = 1.
The eigenvectors φ l capture different timescales of the dynamical system. The slow timescale on the 2-torus has time period T = 2π/ min{1, ζ}, i.e., T = 2π for Model I, and T = 2π × 30 1/2 for Model II. This is equivalent to saying that the slow timescales happen along ϑ 1 for Model I, and along ϑ 2 for Model II.  are able to capture the characteristic dynamical patterns of each model, with swirl patterns for Model I and azimuthal patterns for Model II. In the R 3 standard embedding representation of the 2-torus, the slow timescales for Model II happen along the third dimension f 3 . Our approach captures these slow timescales despite the fact that the system is observed only partially through f = (f 1 , f 2 ) without any information on f 3 .

Kernels based on Euclidean distances in the ambient data space
For comparison, we also computed the eigenvectors directly in the full ambient data space using k = 10, 000 nearest neighbors, a standard Gaussian kernel with width ǫ = 0.25, and Diffusion Maps normalization with α = 1. We used the Euclidean distance to compute the pairwise distances. The eigenvectors for Models I and II are shown in Figs. 5 and 6, respectively. The distances are computed using the fully observed system in R 3 , F = (f 1 , f 2 , f 3 ). However, despite the fact that the models generating the data are different, the embeddings are identical for the same set of parameter values. Thus, the ambient data space is not capable of capturing through the eigenvectors the difference in the dynamical evolution of the two systems.

Dynamics on the 2-torus with a fixed point
For this second experiment, we consider the dynamical system on the 2-torus proposed by [46] with a fixed point and the vector field v = (v 1 , v 2 ) given by where ζ is an irrational frequency parameter. The flow in (37) has a fixed point at coordinates ϑ 1 = ϑ 2 = 0. Trajectories along this dynamical system pass by the fixed point at arbitrarily small distances, but they circumvent the fixed point by developing "bumps" around the fixed point. The standard (flat) embedding of the 2-torus in R 4 is given by  We generate N = 64, 000 points from this dynamical system for the frequency ζ = 20 1/2 and a timestep δt = 0.01. The probability measures were estimated using KDE with an embedding window of R = 40 timesteps and Q = 50 evaluation points per dimension. The parameters of the Diffusion Maps algorithm were set to k = 3, 000 nearest neighbors, the width of the Gaussian kernel ǫ = 1, and the normalization parameter α = 1.
The leading eigenvectors for the system partially observed through f = (f 1 , f 2 ) instead of the full observational map F from (38) are shown in Fig. 7. Their associated scatterplots, i.e., spatial patterns, follow the orbits of the dynamical flow, despite partial observation of the system. We performed the same analysis for all partially-observed systems through all 1D and 2D combination of observables of the full coordinate space F = (f 1 , f 2 , f 3 , f 4 ). The results are consistent for the different cases, and the eigenvectors capture mainly the slowly-varying timescales of the system, i.e., trajectories along the orbits of the flow.

Lorenz attractor
In the third experiment, we consider the Lorenz 63 mathematical model initially proposed as a simple model for atmospheric convection, which consists of three ordinary differential equations: Here ω 1 , ω 2 , ω 3 are the system states, and σ, ρ, β are the system parameters. We consider here the typical parameter values for the Lorenz system: ρ = 28, σ = 10, β = 8/3. The embedding in R 3 is given by We generated N = 66, 828 points starting at the initial point (0, 1, 1.05) for the time interval ∆T = [0, 500], after having removed the first 150 transient points. The probability measures were estimated using KDE with an embedding window of R = 30 timesteps and Q = 50 evaluation points per dimension. The parameters of the Diffusion Maps algorithm were set to k = 2, 000 nearest neighbors, the width of the Gaussian kernel ǫ = 0.32, and the normalization parameter α = 1. The Lorenz 63 system is highly nonlinear and non-periodic, and feature extraction is therefore a challenging problem.
The leading eigenvectors for the system partially observed through f = (f 1 , f 3 ) and f = (f 1 , f 2 ) are shown in Figs. 8 and 9, respectively. The eigenvectors φ l capture different patterns of the slowlyvarying timescales of the system, while faster-varying timescales emerge as we go deeper in the eigenfunction spectrum. We observe that some of the eigenvectors in Figs. 8 and 9 overlap (with sometimes only a change of sign), indicating that the framework proposed here recovers the underlying dynamics under different partial observables of the system. Figure 10 shows examples of one-dimensional time series and two-dimensional representations of the eigenvectors for the system partially observed through f = (f 1 , f 2 ). In this case, eigenfunction φ 2 represents the two wings of the Lorenz attractor, i.e., positive and negative values correspond to the left and right wing of the attractor, respectively (see also Fig. 9); eigenfunction φ 3 represents the variation within each wing with positive values corresponding to points further apart from the intersection of the wings, while negative values correspond to points closer to the intersection; and eigenfunction φ 4 represents a switch between the wings.

Moment reconstruction for the realtime multivariate MJO (RMM) index
The phenomenon from climate science that we study here is known as the Madden-Julian oscillation (MJO; [48,49]). MJO is the main tropical intraseasonal oscillation (ISO), and it corresponds to a 30-90-day eastward-propagating wave pattern with zonal wavenumber 1-4.
The observation map. Among the multitude of indices that measure the MJO, the most common is the realtime multivariate Madden-Jullian oscillation (RMM) index [50]. RMM is a combined measure of the first two empirical orthogonal functions (EOFs) -or principal components (PCs) -of bandpass-filtered, and equatorially averaged outgoing longwave radiation, and 200hPa and 850hPa zonal wind data. In this experiment we use our framework to extract temporal patterns from the RMM index, and then show that the first four moments (mean, standard deviation, skewness and kurtosis) of the PDFs can be accurately reconstructed using only a small number of the leading eigenvectors.
The dataset covers 23 years from September 1983 to June 2006, sampled once a day δt = 1, that is, a total of N = 8337 samples. We set the number of nearest neighbours to k = 100, and ǫ = 0.02 as the width of the Gaussian kernel. We choose the embedding window to be ∆t = 60 days (R = 60 timesteps) as it represents the average time of an MJO (30-90 days). Figure 11 (Left) shows the absolute values of the expansion coefficients c l of the M = 50 leading eigenvectors computed using (27) for the first four moments of the distributions -the sign of the coefficients depends on the sign of the eigenvectors. We estimate the reconstruction error in Fig. 11 (Right) for the leading M = 50 eigenvectors, using Root Mean Square Error (RMSE) between the observed and the reconstructed time series: where E n is the observed n-th moment of the distributions, andÊ n is the reconstructed n-th moment using only the leading few eigenvectors as in (26). Examples of reconstructions of the first four moments using M ∈ {5, 15, 50} eigenvectors are shown in Figs. 12, 13, 14, respectively. The optimal value of the number of eigenvectors needed to reconstruct each moment can be chosen using the decay of the error in Fig. 11. Thus, the mean and standard deviation are faithfully reconstructed using a very small number of eigenvectors (here approximately M = 5), while the skewness and kurtosis need a slightly higher number of eigenvectors (approximately M = 15 and M = 50, respectively). The correlation between the first nonconstant eigenfunction φ 2 and the mean of RMM is 0.9611 (both have been normalized). Thus, some of the eigenvectors detected using our framework recover intrinsic properties of the statistical manifold, i.e., here φ 2 recovers the mean of the PDFs on the manifold. We also found that some of the eigenvectors are strongly correlated (≈ 0.6 − 0.7) with the standard deviation of the RMM index. Having the mean and the standard deviation could be very useful for example for prediction and uncertainty quantification when the past trajectory of the dynamical system is known. We tested our algorithm for robustness using the following parameter values: ∆t ∈ [30,90], k ∈ [50, 1000], ǫ ∈ [0.005, 2], and the results are very robust within these ranges.

Conclusion
In this paper, we introduced a novel framework for feature extraction and moment reconstruction in dynamical systems that integrates ideas from machine learning, dynamical systems theory, and information geometry. Through dimension reduction, we extract temporal and spatiotemporal patterns of interest that describe the evolution of the dynamical system. While more conventional approaches for dimension reduction act directly in the original data spaces, our method acts on probability spaces. We use divergences between probability distributions on finite-time trajectories of the dynamical system, and thus are able to capture the dynamical evolution of the system. The divergences in the probability space allows us to define a kernel integral operator whose orthonormal eigenfunctions capture different timescales of the dynamical system (with emphasis on the slow timescales). If the collection of probability measures is in addition a manifold, we can equip this statistical manifold with the canonical Riemannian metric, allowing us to make connections to the field of information geometry.
One of our main results shows that linear transformations on general functions over observables of the dynamical system can be written in terms of a time-averaging operator. For particular choices of these functions, these time-averaging operators are used to compute the moments of the collection of probability measures. We next exploited the fact that these transformations are linear combinations of the eigenfunctions of a kernel integral operator, and we showed that we can expand the moments of the distributions in this eigenfunction basis. This property provides a powerful tool that allows us to use nonparametric forecasting techniques based on out-of-sample extensions to predict the moments of the distributions at times further down the trajectory. We applied these techniques to three toy examples and a real-world atmosphere-ocean time series of the Madden-Julian oscillation, and showed that the (first four) moments of the distributions are reconstructed faithfully using only a few leading eigenfunctions.
The framework presented here opens up multiple possibilities for future work. For example, complex phenomena where data is generated by multiple heterogeneous sources, i.e., different ambient data spaces with different units, pose challenging issues for problems such as manifold alignment or multimodal data fusion and integration. Being intrinsically homogeneous, techniques that act directly on probability spaces allow for a coherent analysis of heterogeneous data. In this paper, the probability measures are estimated through kernel density estimation techniques, however this becomes intractable for high-dimensional observables (we worked here with one and two-dimensional observables). Kernel mean embedding of distributions [34,51] are nonparametric techniques that map probability distributions into a reproducing kernel Hilbert space without requiring an explicit estimation of the probability distributions. Integration of these techniques with our framework would be very useful when dealing with high-dimensional observables. We work here with ergodic deterministic dynamical systems, but the framework could also be extended to stochastic dynamical systems [52]. In terms of applications, the framework presented here can be used to extract temporal patterns and forecast moments of the probability distributions of any dynamical system. In the future, we plan to apply the method to study other observables of the climate atmosphere ocean system, such as state-of-the-art El Niño indices.

A. Choice of parameters
In Fig. 15 we display the distance distributions of the first k nearest neighbors for the four models in Fig. 1 for different values of k, and we discuss an empirical way of choosing the optimal number of nearest neighbors in order to obtain good embeddings. As the number of nearest neighbors k increases, the distance distributions shift from a positive skewness towards a negative skewness. In our experiments we observed that the best results were obtained when minimizing the skewness of the distance distributions of the nearest neighbors (here we used a symmetrized distance matrix). This is in part explained by the fact that a large positive skewness corresponds to a neighborhood graph where only very small neighborhoods are connected, while a large negative skewness indicates that there are a lot of edges connecting far away neighborhoods, thus leading to dense graphs. On the other hand, the neutral (no) skewness gives equal weight to both small and large distances creating a balanced neighborhood graph.
For our experiments, we chose the values of k that minimize the skewness, that is k = 500 for Model I on the 2-torus, k = 7, 000 for Model II on the 2-torus, k = 3, 000 for the fixed point torus, and k = 2, 000 for the Lorenz system. For all models in Fig. 1, the results were robust for ranges of values that guaranteed a small skewness. Thus, for Model I embeddings using a number of nearest neighbors in the range k = [500, 1000] were performing the best, as were embeddings using k = [5000, 10000] for Model II.
Concerning the width of the Gaussian kernel, its choice is tightly related to the mean of the distance distribution to the k nearest neighbors, which is around ε = 0.5, with the bandwidth in (9) ǫ = 2ε 2 . In the numerical experiments, we show results using ε that provides the best visual embeddings, but our results have shown to be robust for values in the range ε = [0.2, 1]. Kernel similarities for five random data points for each of the four models are shown in Fig. 16. We see for example that the kernel decays significantly faster for Model I compared to Model II on the 2-torus, thus directly influencing the choice of the optimal parameter values (both k and ǫ). By truncating the kernel at k nearest neighbors, we remove the smoothness, however previous results [53] have shown convergence of the graph Laplacian even for non-smooth kernels. The Hellinger distance being upper bounded by 1, the kernel similarities will be lower bounded by e − 1 2ε 2 . In this paper, the parameters k and ǫ are global to all data points, but as the densities and the neighborhoods on the manifold change (e.g., for Model I the flow along the manifold speeds up or slows down in different regions leading to different local behaviors), one solution would be to use an adaptive nearest neighbor and a variable-bandwidth kernel [53]. We believe adaptive algorithms would improve our results as for example we see that the decay of the kernel similarities behaves differently for different data points (Fig. 16). We leave this analysis for future work.