Computer Science Mathematics Physics
Computer Science-2021.gz:143391
Computer Science
270278ceaf88107689f99d81a60a517bc96dc35d
Kernel Analog Forecasting: Multiscale Test Problems

Data-driven prediction is becoming increasingly widespread as the volume of data available grows and as algorithmic development matches this growth. The nature of the predictions made, and the manner in which they should be interpreted, depends crucially on the extent to which the variables chosen for prediction are Markovian, or approximately Markovian. Multiscale systems provide a framework in which this issue can be analyzed. In this work kernel analog forecasting methods are studied from the perspective of data generated by multiscale dynamical systems. The problems chosen exhibit a variety of different Markovian closures, using both averaging and homogenization; furthermore, settings where scale-separation is not present and the predicted variables are non-Markovian, are also considered. The studies provide guidance for the interpretation of data-driven prediction methods when used in practice.

1. Introduction. Data-driven prediction holds great promise in many areas of science and engineering. Growth in the volume of data available in numerous application areas has been matched by advances in computational methodologies which are designed to utilize this data for prediction. However fundamental questions arise in this field relating to the choice of variables on which to base prediction, and whether or not the system is Markovian in the chosen variables. Whilst delay embedding can be used to enhance the choice of variables in which Markovian structure is present, prediction is often undertaken using variables in which there is not a Markovian closure or in which this closure is only approximate. The objective of the paper is to use multiscale systems to provide a framework in which the fundamental issue of the role of Markovianity in data-driven prediction can be studied. We work within the setting of kernel analog forecasting (KAF), a methodology that has seen success in a number of application domains, and which is backed by a mature theory. Subsection 1.1 provides an overview of relevant literature in data-driven prediction for dynamical systems and the multiscale setting in which we work. We outline our contributions to the understanding of data-driven prediction within multiscale systems in Subsection 1.2.
1.1. Background And Literature Review. In 1969, Lorenz originally introduced the idea of analog forecasting for prediction of dynamical systems using historical data [27]. Given initial data, the method locates its closest analog among the historical points and reports the historical value of the corresponding observable, shifted by the desired lead time. By construction, analog forecasting avoids model error, but the resulting forecast is not continuous with respect to initial data. It is therefore non-physical and this fact, when combined with the paucity of data available at the time the method was proposed, limited the value of the methodology in practice. An exponential growth in data volume has precipitated the development of improved methodologies which build on Lorenz's original idea, leading to algorithms which are backed by large data theories and which depend smoothly on initial condition. This has been achieved through the use of kernel based methods which result in data-driven prediction based on weighting all the historical data according to its similarity to the initial data; this leads to algorithms which enforce continuity of the forecast with respect to initial data [48] and, building on the theory of reproducing kernel Hilbert spaces [3], to algorithms which can be theoretically justified in the large data limit [15,1].
KAF draws upon several fundamental ideas rooted in kernel methods for machine learning. First, the choice of kernel function is guided by the need for dimension reduction of big data and the specific learning task at hand. For clustering, similarity kernels [38] are used to construct graphs over the data, and clusters determined by its graph Laplacian eigenvectors [4]. This graph Laplacian construction is generalized in [10] to characterize diffusion operators on the manifold upon which the data lie via their eigenfunctions, known as diffusion maps, and further generalized in [7] to a class of variable-bandwidth kernels that control for variations in the density of the sampling distribution. Under different choices of kernel and normalization, the resulting eigenmaps can describe slow coordinates in dynamical systems [32], and can also be used as a basis to approximate evolution operators in stochastic differential equations (SDEs) [6]. Algorithmic development has been aided by advances in the theory for pointwise [40] and spectral [45,42] convergence of these coordinates in the large data limit. Secondly, Markov operators constructed from kernels map into reproducing kernel Hilbert spaces (RKHSs) in which kernel evaluation corresponds to function evaluation, and allow evaluating these eigencoordinates on out-of-sample data [11] in a procedure known as Nyström extension. This has found use in semisupervised classification using support vector machines [39] to extend labels to new data, spline interpolation [46], and forecasting [48].
In addition to exploiting ideas from kernel methods for machine learning, KAF may exploit additional structure from time-ordered data arising from dynamical systems. For example delay embedding is frequently used to identify Markovian structure. In diffusion forecasting [6], diffusion maps are time-shifted to approximate the action of a shift operator on observables in SDEs. Alternatively, time-shifted diffusion maps can be used to approximate the reduction coordinates of this shift operator directly [15], or, temporal structure can be directly embedded into specialized cone kernels for analog forecasting [48]. The latter takes the Nyström extension perspective of KAF, while in fact, KAF evaluates a conditional expectation of this shift operator, conditioned on the observations [1]. The aforementioned shift operator known as the Koopman operator acts by composing observables with the dynamical flow map, and is a linear operator on these function spaces. Hence, the data-driven approximation of the Koopman operator is an exciting area of research. Bernhard Koopman introduced the linear operator that carries his name in the 1930s as part of his study of ergodic and Hamiltonian dynamics [24]. In data-driven identification of coherent structures spectral decompositions of the Koopman operator [30] and the related transfer operator [13] play a central role, driven by the fact that in such a basis, forecasting of nonlinear dynamics amounts to scalar multiplication by eigenvalues. Algorithms used in practice compute finite-dimensional regression onto pre-computed libraries consisting of functions of the time-lagged snapshots [37,36,23]. However, convergence guarantees are limited, requiring stringent assumptions on the libraries and spectrum [2]. In particular mixed spectra resulting from chaotic/mixing systems pose a challenge for numerical methods. Recent data-driven methods which leverage infinite-dimensional feature spaces provided by kernels [16], as well as kernel constructions in spectral space [25], are able to tackle the continuous part of the spectrum of the Koopman operator. For forecasting purposes, pointwise evaluation of the Koopman operator acting on observables is the natural setting, rather than spectral approximation, and is the perspective we take.
Multiscale analysis provides a setting in which to understand the role of rapidly varying (in space or time) system components on the slowly varying variables used for predictive models [47]. In this paper we will work in the framework of averaging and homogenization for partial differential equations (PDEs) and SDEs, as developed in [5]; see [33], and chapters 9, 10 and 11 in particular, for an exposition of the subject that is adapted to the chaotic dynamical systems setting that is the focus of this paper. We will use the example developed in [29], which exploits the (proven) chaotic properties of the Lorenz 63 model [26,43], to provide an example of a chaotic ordinary differential equation (ODE) which homogenizes to give an SDE. And we will use the multiscale Lorenz 96 model [28] to provide an example of an ODE to which the averaging principle may be applied to effect dimension reduction, as pioneered and exploited in [14]. The work of Jiang and Harlim [21] studies data-informed model-driven prediction in partially observed ODEs, using ideas from kernel based approximation; and in the paper [20] the idea is generalized to discrete time dynamical systems, and neural networks and LSTM modeling is used in place of kernel methods. In both the papers [21,20] multiscale systems are used to test their methods in certain regimes.
Data-driven analog forecasting, kernel methods, and Koopman methodologies have each individually found widespread use in real-world forecasting and coherent pattern extraction applications. Analog forecasting, albeit without kernels, has been used to predict weather patterns [8,12,44], yet is known to have limitations predicting chaotic behavior. Khodkar et al. recently developed a Koopman-based framework using delay embedded observables to predict chaotic dynamics [22]. Nonlinear Laplacian spectral analysis, which applied kernel and delay embeddings akin to Koopman observables, successfully recovers intermittent phenomena such as the El Niño Southern oscillation [41]. Koopman operator approximation has also been widely adopted to study high-dimensional complex, even turbulent fluid flows [17], see [31] for a study of these applications.
In this paper we use KAF as developed in the papers [1,15,48]. Technical details on the implementation, including pseudocode and further relevant references, are collected in Appendix A.
1.2. Our Contribution. We use multiscale methodology to introduce four classes of ODE test problems which exhibit Markovian dynamics after elimination of fast variables; stochastic, chaotic, quasiperiodic and periodic behavior may be obtained in the slow variable, depending on the setting considered. Using these test problems, our four main contributions to the understanding of KAF are as follows: 1. We apply KAF techniques to data generated by each of these four classes of multiscale test problems, and use the behaviour of the averaged or homogenized slow system to interpret the resulting predictions. In particular, KAF methods converge, in the large data limit, to a conditional expectation defined via the Koopman operator of the multiscale systems; we use this as the basis for our interpretation. Moreover, we demonstrate the use of KAF to estimate the variance of the predictor itself. In each of the four cases the 2σ-interval captures the real trajectory, even when the KAF predictor does not track the trajectory itself. This can be viewed as a separate application of the KAF methodology which will be useful in cases when forecasting of high probability bounding sets suffices even when the trajectory itself is hard to predict. In all cases we also study problems in which the scale-separation is not present, but KAF prediction of mean and variance is attempted on the basis of data from only a subset of the variables. 2. The KAF method is based on data-driven approximation of the eigenvalue problem arising from a kernel integral operator. In the setting in which the multiscale ODE homogenizes to produce an SDE corresponding to a bistable gradient system with additive noise, a limiting analytical expression is available for the eigenfunctions; we demonstrate that these limiting eigenfunctions are well-approximated by the data-driven method. This comparison gives insight into the empirical methods used to tune free parameters within KAF. 3. In the setting in which the multiscale ODE averages to produce an ODE of lower dimension we use alternative data-driven ODE closures as a benchmark against which to compare the purely data-driven KAF methods. This gives insight into the relative merits of purely data-driven prediction, and prediction which combines model-based knowledge with data. 4. We use the insights from these carefully constructed numerical experiments to make recommendations about deployment, and parameter-tuning, of KAF methods to real data. The paper is organized as follows. In Section 2, we outline the data-driven construction of the prediction function using KAF methodology. We explain the sense in which the construction converges to a conditional expectation defined via the Koopman operator associated to a measure-preserving dynamical system assumed to underlie the data. We also describe two kinds of canonical multiscale systems which give rise to homogenization and averaging effects, and which we use to provide interpretation of this conditional expectation. Section 3 introduces a test problem in the form of a double-well gradient flow driven by chaotic Lorenz 63 dynamics which homogenizes to give an SDE in the scale-separated regime; numerical results applied to prediction of the slow variable exhibit contributions 1 and 2. In section 4, we introduce the multiscale Lorenz 96 system which averages to give an ODE in the slow variables; three different parametric regimes give rise to periodic, quasiperiodic, and chaotic responses in the slow variable. The behaviour of KAF-based prediction in these three regimes is studied, to illustrate contribution 1; and a slow-variable closure model, built using Gaussian process regression, is compared with the KAF to illustrate contribution 3. In section 5 we overview the insights obtained by studying KAF methods through the lens of multiscale systems; and we then make concrete recommendations about interpreting the output of KAF techniques when applied to naturally occurring data, contribution 4.

2.
Methodology. In subsection 2.1 we overview the two key ideas which interact to underpin the studies in this paper: KAF and multiscale methods, tailoring the exposition to the use of the latter as a tool to understand the former. We then give more details on KAF. The two primary components of the KAF methodology are: (i) viewing forecasting as evaluation of a conditional expectation of the Koopman operator applied to the desired observable; (ii) approximation of this conditional expectation in a data-driven fashion. Subsections 2.2 and 2.3 describe (i) and (ii) respectively, whilst subsection A.2 is devoted to a key practical component of the data-driven approximation, namely construction of the kernel, and subsection A.3 to the data-driven choice of integer , the number of (approximate) eigenfunctions used in the data-driven forecast.
2.1. Overview Of Methodology. The problem setting for prediction is as follows. We assume that we are given N time-ordered data samples where x : R → X is a continuous time process, x n = x(n∆t) and ∆t is the sampling rate. We assume that the continuous time process x in X is derived from Markovian dynamics for a coupled pair (x, y) evolving in the larger state space X × Y. Assume that the desired prediction lead time τ is an integer multiple of the sampling interval, that is, τ = q∆t. Included with the data are values of the associated prediction observable advanced by τ time units defined by the Markovian dynamics via an unknown map F : X × Y → R; thus f n = F (x n , y n ). The goal of KAF is to predict F (x(τ ), y(τ )) given only partial information, x(0) = x, and the N data samples x n . We view the data-driven predictor as a map Z τ : X → R which takes initial condition x as input. Given initial data x and lead time τ , the KAF predictor averages over the τshifted observable values provided in the training data and weighted by a kernel p : X × X → R constructed from the data; the resulting algorithm has the following form: The weighting kernel p(x, x n ) determines how much weight to attach to a time-series initialized at point x n , according to its proximity to x, the desired initial point. The features φ j are computed from an eigenvalue problem associated with a data-driven approximation of a kernel integral operator, constructed from x n ; in the large data limit this provides an orthonormal basis for the entire space. The function ψ j is an out-of-sample Nyström extension of φ j , orthonormalized with respect to an underlying RKHS structure. The method may be seen as a smoothed version of Lorenz's original proposal for data-driven prediction -analog forecasting [27]. In particular it addresses the issue of continuity of the prediction with respect to the initial condition, and it does so in a framework which is provably statistically consistent in the large data limit [1,15,48]. Further details of the methodology are given in the next two subsections, and the attendant information in Appendix A. An important challenge addressed by this methodology is that, since the y component of the system is not observed, the sequences {x n } and {f n+q } are non-Markovian. As a consequence the standard idea of constructing a Markov chain from the data is not natural. The kernel analog forecasting method evaluates a conditional expectation of the forecast conditioned, using the observed data {x n }, explicitly incorporating information loss resulting from unobserved y; it is hence a natural approach to the problem at hand. Multiscale systems provide a natural setting for the study of KAF methods, and in particular the issue of prediction of non-Markovian or approximately Markovian systems. In this paper we will consider the variable x as the slow component of a Markovian system for pair (x, y) in which y evolves as a fast variable. We consider averaging and homogenization settings in which the dynamics for x is approximately Markovian, and the conditional expectation arising in the KAF method may be understood explicitly. This will enable us to obtain a deeper understanding of how KAF works, and help users of the methodology interpret it. We now outline the averaging and homogenization settings that we will use. Details of the theory underlying them may be found in [33].
We will study multiscale systems which exhibit averaging, in the form where B : X → Y is linear. The average of By under the invariant measure of the y dynamics, with x frozen, provides a closed approximate ODE dynamics for x, when ε is small. If we denote the x−parameterized invariant measure for the y dynamics with x frozen by ν x (dy) then for ε When the variable By averages to zero a different scaling is required to elicit the effect of the fast variable on the slow one. To this end we also consider multiscale systems which exhibit homogenization, in the form Here we assume that Y By ν(dy) = 0 where ν is the invariant measure of the y dynamics. The approximate dynamics for x, when ε is small, is then governed by an SDE in this setting, provided y is chaotic and mixing. If this is the case then, invoking the homogenization principle, x ≈ X where X is governed by an SDE of the form where W denotes the Wiener process, and σ is a uniquely determined positive constant that can be computed numerically from the mixing properties of the y process.

Koopman Formulation Of Prediction.
We let Ω = X × Y and assume that Φ t : Ω → Ω is an ergodic dynamical system with invariant probability measure µ; we assume t ∈ R + but the extension to discrete time is straightforward. Define the continuous observation map H : Ω → X and the prediction observable F : Ω → R; we assume that F is square-integrable with respect to the invariant measure: Comparison of exact Koopman picture and data-driven approximation. The approximate mapping Zτ (·) is constructed from the data streams We define the Koopman operator U t : We seek, in a sense to be made precise, the function Z τ : X → R such that Z τ • H is the best approximation to a perfect prediction of F • Φ τ . We formalize this by introducing the Hilbert subspace V ⊆ L 2 µ (Ω; R) given by This Hilbert space captures the notion of making predictions based only on information in X . Note that the perfect forecast would satisfy Z τ • H = U τ F , but that such a forecast will not be possible in general because X is a proper subset of Ω, and information needed for perfect prediction of F will be missing. Among all elements of V , the minimal prediction error in L 2 µ (Ω; R) is attained by the conditional expectation This formulation of prediction encapsulates the inherent loss of information incurred through observing only a set of functionals of an ergodic dynamical system, and the effect of this loss of information on prediction. In subsequent sections of this paper we will assume that F = F • H for some F : X → R because it is often natural to try to predict only functionals of the slow variables. Note, however, that the methodology is not restricted to such F and in this subsection, the next subsection and in Appendix A we describe the more general setting for completeness.

Data-Driven Approximation.
The formulation in the preceding subsection encapsulates the inherent loss of predictive power incurred through observing only a set of functionals of an ergodic dynamical system. This is formalized by seeking the best approximation of the Koopman evolution from within a Hilbert subspace capturing the notion of depending only on specified functionals on the state space of the dynamical system. We now demonstrate how data may be used to further approximate this best approximation, and to do so in a manner which is refineable as more data is acquired. The approach is summarized in Figure 2.1.
For observation time ∆t > 0 we define We assume that we are given time-ordered pairs and the objective is to construct, from this data, a function Z τ : X → R which predicts F at lead time τ so that Z τ • H ≈ g where g solves the minimization problem in (2.2). Furthermore we wish to carry this out in a manner which ensures that, in an appropriate topology, To this end we introduce a hypothesis space H ,N , of dimension and depending on the N −dependent data set (2.3), and seek to solve the minimization problem The choice of the hypothesis space is constrained by the need to be able to solve the minimization problem (2.4) explicitly, using only the data (2.3), and by the requirement that Z τ • H recovers g in the large data limit N → ∞. Moreover, in order to be practically useful forecast functions, elements of H ,N should allow pointwise evaluation at any x ∈ X , which is not defined in arbitrary subspaces of L 2 µ (Ω; R). With these considerations in mind, we introduce a kernel function k : X × X → R and RKHS K with the properties We then define H ,N as an -dimensional subspace of K, to be described below. We also note that the kernel k is constructed from a data-stream of length N , but we suppress the explicit dependence of k on N in the notation. In Appendix A we discuss our data-driven construction of k, and choice of . For now we proceed on the assumption that we have a kernel, and hence a RKHS, as well as a method for choosing .
n=0 δ ωn be the sampling measure underlying the training data (2.3) and define The eigenvectors of this matrix lead to an orthonormal basis We may also identify each element φ j ∈ R N with element φ j • H ∈ L 2 µ N (Ω; R) via the definition φ j (x n ) as the n th entry of the vector φ j ∈ R N . Using the same symbols for elements of R N and L 2 µ N (Ω; R), as well as for linear transformations on those spaces, is a useful economy of notation. Then the following functions ψ j : X → R form an orthonormal set in K: This is a form of Nyström extension [11].
As hypothesis space we take noting that the basis functions themselves depend on the data set, and hence on N .
We may now solve the optimization problem (2.4) and an explicit computation yields, for τ = qT , Note that this construction of the predictor Z τ is entirely data-driven: the basis functions ψ j and the eigenvalues λ j are found from the eigenvalues and eigenvectors of the data-defined kernel matrix; and the coefficients c j are computed as sums over the data set. Furthermore, Theorem 14 in [1] proves that Z τ • X converges to g , the solution of the minimization problem (2.2), as N → ∞, followed by → ∞, in an L 2 sense with respect to the invariant measure µ on Ω.
More generally, any function of the observable can be predicted in this data-driven manner, which provides a convenient framework for uncertainty quantification. The conditional variance between forecast and ground truth can also be computed in the hypothesis space as in (2.8) using the coefficients For detail on the data-driven kernel construction, the data-driven choice of and the conditional variance estimator see subsections A.2, A.3 and A.4 respectively.

Homogenization: Lorenz 63 Driven
System. This section is devoted to the setting in which a chaotic ODE of form (H) is approximated by an SDE of form (H0). The goal is to make predictions of the x variable, using data concerning only the x variable from (H); the role of (H0) is simply to help us interpret those predictions. This setting presents unique challenges for forecasting as one cannot expect the outcome of any method to predict a sample path of a stochastic process without knowledge of the driving noise. This fact has direct bearing on prediction in (H) using x-data alone, since (H0) demonstrates that the time series of the x−data is approximately that of an SDE; without knowledge of the noise, which is governed by the unobserved y variable, prediction of the trajectory of x is not possible. In subsection 3.2 we examine instead the long-term statistics predicted by KAF from data generated by (H) -the conditional expectation and variance of the stochastic process -and compare them with estimates computed from (H0) using Monte-Carlo simulation of the SDE. This illustrates our main contribution 1 from the list in subsection 1.2. Then, in subsection 3.3, exploiting the fact that the limiting process is one-dimensional, we find explicit expressions for the kernel eigenfunctions in the limit problem (H0) and compare these with the eigenfunctions obtained from data-driven techniques applied to (H), our main contribution 2 from subsection 1.2. Subsection 3.4 is concerned with non-Markovian prediction, in which there is no scale-separation between observed and unobserved variables. We start, however, in subsection 3.1, introducing the concrete model around which our experiments are organized.
3.1. The Model. The first test problem arises from driving a double-well gradient flow with a chaotic signal obtained from the Lorenz 63 model [26]: This is of form (H). In [29] it is proved that as ε → 0 this system converges weakly in C([0, T ]; R), when projected onto the x variable alone, to solution of the SDE Thus, this white-noise driven gradient system is of form (H0). The value of the constant σ is identified in [18]. For the current work the key point to appreciate is that for small ε the variable x in (3.1) exhibits (approximately) Markovian behaviour, but this behaviour is stochastic. The SDE is ergodic and has invariant probability density function with respect to Lebesgue measure.

Conditional Expectation and Variance.
We aim to predict the x variable from historical data of a long trajectory of x alone. Thus the observation and observable maps are H(ω) = x, F (ω) = x. We will also estimate second moments, enabling us to compute conditional variance, for which F (ω) = x 2 . Observation data is generated by using an implicit time-stepping scheme with time-step 0.01 in the slow variable and built-in Matlab solvers to integrate the fast variables with ε = 0.001. Source data for the slow variable x is gathered for N = 40000 points sampled at the macroscopic time interval ∆t = 0.05. ThenN = 7500 out-of-sample points from a new trajectory {x n } are gathered at the same resolution, which provide the ground truth for assessing forecast error. The natural error metric is the root mean squared error (RMSE ), the L 2 norm of Z τ − U τ F. To account for differences in scale, we normalize the RMSE by the standard deviation of the trajectory. Figure 3.1 depicts the behavior of the KAF forecast Z τ (x) as a function of lead time τ , for fixed x. This forecast exhibits two interesting properties which can be understood through the small-ε homogenization limit. The first relates to the fact that the trajectory itself is not well-predicted; the second explains what is well-predicted.
Firstly, the predictor tracks the conditional mean initialized at x = −1.10, and not the trajectory itself. This is predicted by the theory, since what is predicted is the   Comparison of uncertainty for different initial data. Blue is predictor (conditional mean), grey is trajectory and pink gives two standard deviations bounds computed from the conditional variance. The forecast uncertainty when initialized at x = 0 (right) displays more rapid growth in uncertainty over short lead times than that when x = −1.10 (left). This sensitivity to initial conditions is a desirable feature of KAF, in light of the fact that both plots share the same training set.
long-term conditional expectation E[U τ F |H]. Indeed this latter quantity necessarily converges for large τ to a constant, under mixing assumptions on the (x, y) system, whilst individual trajectories in x exhibit stochastic dynamics, approximately that of X. This explains the growth in RMSE seen in Figure 3.2. Secondly, exploiting the fact that we expect the system (3.1) to behave like (3.2), when projected onto the x co-ordinate, we can provide an objective evaluation of the KAF forecasts by running Monte-Carlo simulations of the SDE for X; to do this we compute the sample mean and variance over 10000 sample paths initialized at the same initial point x = −1.10. Figure 3.1 compares the KAF forecast mean and variance with that predicted by Monte-Carlo mean and variance for the SDE, and they are seen to agree very well over the entire window of computation.
Finally in Figure 3.3 we use the possibility of varying the initial condition in the KAF to demonstrate that the variability encapsulated in the conditional variance is able to pick-up different sensitivities, depending upon initial condition. The panel on the left shows a trajectory initialized at x = −1.10 and the panel on the left shows a trajectory initialized at x = 0. As can be expected from the limiting SDE, the uncertainty when initialized at x = 0 is greater and this is manifest in the conditional variance.
Examination of the KAF technique in this homogenization setting thus clearly reveals the inability of the method to predict trajectories, but shows that it can accurately approximate statistics of trajectories, averaged over the unobserved component of the system. Furthermore, analysis of the SDE provides a means of characterizing the geometry of the underlying hypothesis space, as seen in the next subsection.

Insights Into The Hypothesis
Space. By studying the large data and small kernel bandwidth limit of the matrix G defined in (2.5) we get insights into the structure of the hypothesis space (2.7). The theory in [15], building on the papers [10,7], demonstrates that, for the choice of kernel described in Appendix A, the vectors of G are approximated in the large data and small bandwidth limit by eigenfunctions of the Laplace-Beltrami operator ∆ h on M , associated to a metric h. For the diffusion process (3.2) in dimension m = 1 the manifold M is simply R and the metric h is given by h = ρ −2 , with invariant density ρ given by (3.2). The fact that the density ρ is approximately available to us through the time series x generated by (3.1) is demonstrated in the left panel of Figure 3.4, where we compare the histogram generated by the data with ρ given by (3.3). The conclusion of these various approximations is that we expect the eigenvectors of G, based on data x from (3.1), to be well-approximated by eigenfunctions of ∆ h on R, with h = ρ −2 . We now demonstrate that this is indeed the case.
The action of the Laplace-Beltrami operator on f ∈ C ∞ (M ) is given by where div µ and grad h are the divergence and gradient operators associated with µ and h, respectively. Using the above, we solve the eigenvalue problem −∆ h ϕ = λ 2 ϕ directly. We make the substitution dY = ρdX, mapping X ∈ R into Y ∈ [0, 1]; we note that Y has interpretation as the cumulative distribution function coordinate of ρ. In terms of Y we have Noting that the natural boundary conditions for the Laplace-Beltrami operator are of no-flux type, it follows that, when viewed as functions of Y , the eigenfunctions of  The solutions are the well-known harmonics cos(kπY ), and corresponding eigenvalues λ k = kπ, k ∈ N. Changing back to variable x we obtain We now verify that, for large data sets and small bandwidth, the eigenfunctions of G are indeed close to those associated with the Laplace-Beltrami operator ∆ h . This is demonstrated in the middle and right panels of Figure 3.4. The middle panel shows the first six eigenfunctions of G, computed from data derived from the x variable in (3.1); the right panel shows the eigenfunctions of ∆ h for diffusion process (3.2) in variable X which governs the limiting behaviour of x in the scale-separated case. The agreement is excellent, demonstrating that the heuristics underlying parameter choices within the kernel based methodology (see Appendix A) work well in this setting.

Non-Markovian Regime.
In the preceding subsections we studied predictors for x, based only on time-series data in the x coordinate, for the equation (3.1). We studied the scale-separated regime where ε 1 and x is approximately Markovian -it is approximately governed by an SDE. Here we study the behavior of identical predictors when ε = 1; the system (3.1) then no longer exhibits homogenization and x is no longer Markovian in view of the lack of scale-separation. As a result, the prediction of x, shown in the top of Figure 3.5, is poor even at very short times, and the two standard deviation confidence bands reflect this rapid initial error growth, and then remain large throughout the time window. Z τ converges rapidly to the conditional mean. To render the problem Markovian we may include more data, and in particular the forcing term. To this end we take H(ω) = [x, 4y 2 /90] T ; see the bottom of Figure 3.5. The prediction of x with these augmented observations yields very accurate predictions over a lead time of τ = 3, considerably larger than in the preceding case where H(ω) = x. After τ = 3, however, pathwise predictability again fails, due to the sensitive dependence of solutions with respect to the forcing function. Once again Z τ converges to the conditional mean. However the uncertainty quantification provided by two standard deviation error bars consistently captures the true trajectory, even in this non-Markovian regime. This non scale-separated pair of examples illustrates that prediction of non-Markovian variables is inherently harder than Markovian variables, but that sensitive dependence of trajectories with respect to perturbations also limits predictability, even in the Markovian setting. This second point will be prominent, too, in the next section. The goal is to make predictions of the x variable using data concerning the x variable alone from (A). The role of (A0) is primarily to help us interpret those predictions; however it also serves to motivate a different prediction methodology, which mixes model-based and data-driven methodologies, against which we will compare KAF.
The averaging setting presents interesting opportunities to understand forecasting. In particular, by tuning a parameter in (A), which is also present in (A0), we are able to create settings in which the variable to be predicted is, in turn, approximately periodic, quasi-periodic and chaotic. These different settings give rise to different types of forecasts and we demonstrate this. In subsection 4.2 we examine the long-term statistics predicted by KAF from data generated by (A) -the conditional expectation and variance of one component of the slow variable -and compare them with estimates computed from (A0) using Monte-Carlo simulation. This illustrates our main contribution 1 from the list in subsection 1.2. Then, in subsection 4.3, we compare the purely data-driven method of prediction with a hybrid data-model predictor. This hybrid is computed by using Gaussian process (GP) regression to compute an approximate closure v GP (x) ≈ v(x), in the notation of (A0); for more details on such approximate closures in the context of the specific model we study in the following four subsections see [14] and references therein. The work in subsection 4.3 illustrates our main contribution 3 from subsection 1.2. Subsection 4.4 is concerned with non-Markovian prediction, in which there is no scale-separation between observed and unobserved variables. We start, however, in subsection 4.1, introducing the concrete model around which our experiments are organized.

The Model.
In this section we focus our attention on another chaotic dynamical system, colloquially known as "Lorenz 96 multiscale" [28], which we will simply abbreviate to L-96. Following the notation established in [14], the L-96 equations model K slow variables {x k } K k=1 coupled to JK fast variables {y j,k } J,K j,k=1,1 with evolution given as follows: x k+K = x k , y j,k+K = y j,k , y j+J,k = y j,k+1 .
This is of the form (A). On the assumption that the y-variables, with x frozen, are ergodic, the averaging principle shows the existence of a function C : R K → R K such that, for small ε, the x variables are approximated by X = (X 1 , . . . , X k ) solving with the periodic boundary conditions X k+K = X k and C k : R K → R denoting the k th component of vector-valued function C. This system is of form (A0). Since system (4.1) is index-shift-invariant, it is clear that the closure C k , if it exists, satisfies C k+1 (X) = C k (πX) where π shifts the vector indices by adding one unit, invoking periodicity at the end points. Furthermore, when J is large, empirical evidence [14] suggests that there is a function c : R → R such the approximation C k (X) = c(X k ) is a good one. For the current work the key point to appreciate is that for small ε the variables x in (4.1) exhibit (approximately) Markovian behavior, and this behavior is deterministic and governed by X. However, by tuning F x , different responses arise in the deterministic variable. In the following we fix parameters ε −1 , K, J, h x , h y throughout all our experiments as follows:

Conditional Expectation and Variance.
We aim to predict the x 1 variable from historical data of a long trajectory of x alone. Thus the observation and observable maps are H(ω) = x, F (ω) = x 1 . We will also use F (ω) = x 2 1 when estimating conditional variance. By tuning the scalar parameter F x (not to be confused with function F ) as outlined in the preceding subsection we can obtain periodic, quasi-periodic and chaotic responses in the averaged variable X. It is intuitive that the ability of the KAF to track the true trajectory of the slow variables decreases with increasing complexity; in other words, predictions in the periodic case should be the most accurate whilst those in the chaotic case present a significant challenge. In the experiments that follow the size and sampling interval of the source (training) data remain fixed at (40000, 0.05) and the out-of-sample (test) data set is fixed at N = 7000.
The space of observables X in the current example is the space of all slow variables. Since, under the small-ε limit, an ODE closure of the slow dynamics is obtained, the variable x behaves (approximately) like a deterministic Markov process, and the expectation in (2.2) disappears; the predictor is expected to track the actual trajectory x 1 (t). To see this another way, note that simply knowing the initial values of the x-variables (recall that X is precisely all x-variables) and the closure C(X) in equation (4.2), we are able to predict x 1 (or indeed, any x k ) exactly, given the initial conditions for all x-variables.
However this picture is greatly affected by the sensitivity of the system to initial conditions and sampling errors due to high dimensionality of the attractor. We now describe how these predictions work in practice, in the three regimes shown in Figure 4.1. We display our results in Figure 4.2, where x 1 and standard deviation bands are predicted and compared with the true signal starting from the same point. The long-term predictability in each regime is constrained by the complexity of the underlying Markovian, deterministic, slow dynamics. In the periodic regime, since chaos is absent in the slow variables, a perfect predictor is obtained via the partially observed dynamics; one interpretation of why this occurs is because observables are finite linear combination of Koopman eigenfunctions [2]. Observe that Z τ remains in phase, and the forecast variance is negligible, for long lead times up to the length of the entire out-of-sample trajectory (τ = 350). The quasiperiodic trajectory is tracked imperfectly, but with significant accuracy over the same range of times; errors are visible mainly around the extrema of x 1 as suggested by the phase portrait; the conditional variance reflects the significant accuracy present. Prediction in the fully chaotic regime only tracks the trajectory, however, until a lead time of approximately 1 time unit, exhibiting behaviour at long lead times which is somewhat similar to that seen in the previous, homogenization, section in which the predicted variable behaved as if drawn from a Markov stochastic process. In particular the long-term predictor in the chaotic regime converges to a constant by construction, assuming mixing, and this is consistent with the inherent unpredictability of chaotic dynamics. It is notable that the size of the conditional variance, and the resulting confidence bands, is a useful guideline as to the pathwise accuracy of the data-driven predictor. The observations about the predictability of the system by KAF methods are also manifest in Figure 4.3 which shows the RMSE in each of the periodic, quasi-periodic and chaotic regimes. In the first, periodic response regime, the trajectory is predicted almost perfectly and this accuracy is reflected in the narrow confidence bands. In the second, quasi-periodic response regime, the trajectory is predicted very well, but with growing error reflected accurately in the slowly growing confidence bands. In the third, chaotic response regime, the predictive capability is lost due to sensitivity to initial conditions and this is reflected in the rapidly growing confidence bands and in the convergence of the predictor to a constant, for large τ.
We mention that in the quasi-periodic case the presence of multiple attractors (or multiple lobes of the same attractor), and resulting intermittent switching between these attractors, leads to a loss of predictability that is significant on time-scales much longer than those shown here. For the figure shown here we have ensured that training points and out-of-sample points are gathered from the same (part of the) attractor to maintain accuracy. We train using two different trajectories to gather ample training data.
Recall that at each lead time τ along the horizontal axis there is a potentially different number of eigenfunctions (τ ) used in the data-driven method. See Appendix A.3 for details on the choice of . In the chaotic regime the optimal (τ ) tends to 1 for large times (see subsubsection A.3.2 for an explanation), whilst fluctuates around 50 in the quasiperiodic regime; we obtain ≈ 9 for all τ in the periodic regime.

Comparison Of Data-Driven And Model-Data-Driven Prediction.
The previous subsection concerned purely data-driven prediction of variable x from (A), using only data in the form of a time-series for x. In this subsection we provide  comparison with a different forecasting technique based on a combination of model and data-driven prediction, using data in the form of a time-series for (x, By). Knowledge of By enables the use of Gaussian process regression (GPR) [35] (or kriging) to approximate v(·) by v GP (·) in (A0). Our approach is motivated by the paper [14] which looked at finding such closures for the L-96 model in form (4.1). When applied to (4.1) the methodology leads to an approximate closure for the slow variable X which takes the form subject to periodic boundary conditions X k+K = X k . This should be compared with (4.2), which arises from application of the averaging principle; note that, in addition, we have invoked the hypothesis that C k (X) can be well-approximated by function of c(X k ), as discussed directly after (4.2); and we will determine an approximation c GP for c by GPR.
Explicit details of the procedure we use to build a GP closure are described in Appendix B; here we observe that for training we use tuples {x k (t n ), (By) k (t n )} N n=1 , over all k = 1, · · · , 9. See Figure 4.4 to see the data used (red random subsamples, without replacement, of the total grey data set), and an approximate GP closure c GP determined from that data.
Once we have the closed model appearing in (4.4) we may use it to predict the variable x appearing in (4.1), and we may compare that prediction with the one made by KAF. Figure 4.5 shows the result of doing so. It shows that the KAF approach is superior in the periodic and quasi-periodic settings, but that for predictions of the trajectory itself the model-data based predictor (4.4) is superior to KAF in the chaotic case. Note that the model-data based predictor has access to more data than does the KAF, and requires model knowledge; the KAF is entirely data-driven.
We now dig a little deeper into the comparison. We do this in a systematic way in the periodic, quasi-periodic and chaotic regimes. In each of these three cases we show four RMSE error curves, labelled as follows: a) the standard KAF based on x data alone; b) an enhanced KAF using (x, By) data, the same data used to train the ODE (4.4); c) a prediction using the ODE (4.4); d) a KAF prediction trained on X data alone, generated by the ODE (4.4).  RMSE comparison for the four cases a)-d) described in the text, in the periodic, quasiperiodic, and chaotic regimes. In the periodic regime, KAF (a) is an ideal predictor with negligible growth in error (note the logarithmic scale). In the quasiperiodic response regime, the growth in RMSE with KAF (a,b) is significantly slower than that of the GP-based ODE prediction (c). In the chaotic response regime, the GP-based ODE prediction (c) is more accurate in the near term, yet KAF error stabilizes as the prediction converges to the conditional mean.
Predicted 2σ Fig. 4.7. Prediction in non-Markovian regime, ε = 1, results in much smaller short-term predictability and rapid convergence of the conditional mean to a constant. study the behavior of identical predictors when ε = 1; the system (4.1) then no longer exhibits averaging and x is no longer Markovian because there is no scaleseparation between x and y. Our experiments are conducted with F x = 10. Because of the lack of Markovian behaviour we expect rapid loss of predictability in time, when H(ω) = x, F (ω) = x 1 . The resulting conditional mean and variance, shown in Figure 4.7, confirms this intuition. Indeed the conditional mean is out of phase with the truth at lead time τ = 1, and this is also reflected in the large growth of the conditional variance. Furthermore, the conditional mean tapers to a constant at τ = 6, twice as quickly as it does in the ε 1 setting in which this tapering occurs at τ ≈ 11 (Figures 4.2,4.5).

Conclusions.
1. We have studied KAF for data-driven prediction: • we use multi-scale systems to create dynamical systems in which a subset of the variables (slow variables) evolve in an approximately Markovian fashion; • we study KAF performance for a range of systems in which the slow variables are governed approximately by stochastic, chaotic, quasi-periodic and periodic behaviour; • in the stochastic case we use the homogenized equations for the slow variables to obtain explicit formulae for the eigenfunctions of the operator underlying KAF, and use these to validate the performance of the KAF method; • in the chaotic, quasi-periodic and periodic cases we use the averaged equations for the slow variables to obtain a GPR-based approximate closure model, and use hybrid data-model predictions, from this closure, in order to evaluate the KAF method. 2. What we illustrate about use of the KAF: • when the variable being predicted is (approximately) a component of a stochastic Markov process then prediction of individual trajectories is not possible, whilst the mean and variance, averaged over possible realizations of the stochastic behaviour, can be accurately predicted by KAF; • when the variable being predicted is (approximately) a component of a deterministic but chaotic Markov process then prediction of individual trajectories is also not possible, except over short time horizons; • in both the chaotic and stochastic cases, whilst prediction of individual trajectories is not to be expected, simply bounding the future trajectory may be useful in applications -to this end we show that in all cases two standard deviation bands around the predicted conditional mean always reliably capture the truth; • in both the stochastic and chaotic settings a signature of the lack of predictability is the convergence of the predictor to a constant, for large τ , accompanied by the data-driven choice of parameter converging to one; • when the variable being predicted is (approximately) a component of a deterministic quasi-periodic or periodic process the prediction of individual trajectories over long time horizons is possible; in this case parameter stays away from one, for significantly large τ ; • in all cases the predicted standard deviations around the predictor provide a reliable indicator of the time-scale on which the predictor is accurate trajectory-wise, and on longer time-scales provide a reliable indicator of the scale of the errors incurred, trajectory-wise; • the choice of kernel, and the interpretation of eigenfunctions as harmonics over the observed submanifold, is well-adapted to capturing periodic or quasi-periodic structure whilst the design of eigenfunctions to include the constant eigenfunction permits convergence to the conditional mean for chaotic, mixing dynamics; • the choice of indicates the number of harmonics contained in the predicted variable and is often approximately constant in τ in periodic and quasi-periodic settings, whilst it converges to one in the chaotic or stochastic mixing cases; • when the variable being observed does not evolve in an approximately Markovian fashion then the KAF conditional mean cannot track x for even short times, as observed for ε = 1 for both the Lorenz 63 and Lorenz 96 examples. 3. The work also suggests a number of directions for future study in the area of KAF: • delay embedding can be used to deal with non-Markovian behaviour and it would be of interest to automate the choice of delay embedding dimension within the KAF framework, to get closer to Markovianity; • combining KAF with data assimilation holds the possibility of greater predictability -for work in this direction see [19] which studies the EnKF with a data-driven model update; • it would be of interest to have algorithms to identify slow subspaces, using data living in larger spaces, and hence to conduct KAF using approximately Markovian variables; • it would be of interest to extend KAF predictor Z τ so that it acts on the joint space of initial conditions and key parameters, enabling prediction at as yet unseen parameters. The matrix G is never explicitly formed, instead we exploit the fact that G = SS T , and hence φ are the left singular vectors of S. Thus we bypass the eigendecomposition step with a reduced singular value decomposition (SVD) in step (4). Note that the SVD approach is natural when working with kernels with an explicit factorization G = SS T , including the bistochastic kernels described in Appendix A.2. For more general kernels one typically performs direct eigendecomposition of G. KAF can also be implemented with non-symmetric kernels satisfying a detailed balance condition making them conjugate to positive-definite kernels; see Section 4.2 in [1] for more details. k(·, x n )φ j (x n ), λ j > 0.
• Inputs -Out-of-sample datax 0 ,x 1 , . . .xN −1 ∈ X -Diffusion eigenvectors Φ ∈ R N ×L -Diffusion eigenvalues Λ ∈ R L×L • Outputs -RKHS basis Ψ = [ψ 0 , . . . , ψ L−1 ] • Steps 1. Form matrixŜ ∈ R L×N with entriesŜ ij = κ N (x i , x j ). Note that this requires another kernel density estimation step to evaluate the sampling density on out-of-sample data. where we assume that X is a subset of R d , x, x ∈ X and | · | denotes the Euclidean norm on R d . Note that this kernel is data-independent. From it we can generalize to a data-dependent kernel defined by [7] κ N (x, x ) = exp −|x − x | 2 /(r N (x)r N (x ) ) , Here, , δ and m are lengthscale and dimension parameters, respectively, estimated from the data. The role of including r n , and hence a variable bandwidth kernel, is to compensate for variations in sampling density across the space. A key conceptual idea underlying the construction of q N is that it approximates the Lebesgue density of the measure µ in the large data limit. Thus the kernel κ N weights the distance of x and x in a manner which reflects the sampling density of the data.
From κ N the kernel k N is constructed as follows, invoking a second principle which is to design a bistochastic Markov kernel [9]. Doing so ensures that the top eigenvalue of G is 1 with corresponding eigenvector a constant; then the hypothesis space also contains constants. This is useful for capturing the mean of the predicted quantity U τ F and, in particular, plays a central role in the large τ asymptotics for mixing systems. To achieve the Markov property we proceed as follows. First we define v N (x) = Ω κ N x, H(ω) µ N (dω), Since the above empirically determined quantities only take values at the N sampled points, they are isomorphic to vectors in R N identified by the same name within Algorithm A.1.1. Finally define The operators constructed from the unnormalized and normalized kernels, κ N and k N , are likewise represented by N ×N matrices K and S, respectively, in Algorithm A.1.1.