DataDriven Reduction for Multiscale Stochastic Dynamical Systems ^{†}^{†}thanks: C.J.D. was supported by the Department of Energy Computational Science Graduate Fellowship (CSGF), grant number DEFG0297ER25308, and the National Science Foundation Graduate Research Fellowship, Grant No. DGE 1148900. R.T. was supported by the European Union’s Seventh Framework Programme (FP7) under Marie Curie Grant 630657 and by the Horev Fellowship. R.T and R.R.C. were supported by the National Science Foundation, Award No. 1309858. I.G.K. was supported by the NSF and the AFOSR.
Abstract
Multiple time scale stochastic dynamical systems are ubiquitous in science and engineering, and the reduction of such systems and their models to only their slow components is often essential for scientific computation and further analysis. Rather than being available in the form of an explicit analytical model, often such systems can only be observed as a data set which exhibits dynamics on several time scales. We will focus on applying and adapting data mining and manifold learning techniques to detect the slow components in such multiscale data. Traditional data mining methods are based on metrics (and thus, geometries) which are not informed of the multiscale nature of the underlying system dynamics; such methods cannot successfully recover the slow variables. Here, we present an approach which utilizes both the local geometry and the local dynamics within the data set through a metric which is both insensitive to the fast variables and more general than simple statistical averaging. Our analysis of the approach provides conditions for successfully recovering the underlying slow variables, as well as an empirical protocol guiding the selection of the method parameters.
Key words. multiscale dynamical systems, Mahalanobis distance, diffusion maps
AMS subject classifications. 37M10, 6207
1 Introduction
Dynamical systems of engineering interest often contain several disparate time scales. When the evolving variables are strongly coupled, resolving the dynamics at all relevant scales can be computationally challenging and pose problems for analysis. Often, the goal is to write a reduced system of equations which accurately captures the dynamics on the slow time scales. These reduced models can greatly accelerate simulation efforts, and are more appropriate for integration into larger modeling frameworks.
Following the methods of Mori [28], Zwanzig [44], and others [5, 8, 20], one can reduce the number of variables needed to describe a system of differential equations. However, in general, this reduction introduces memory terms. It transforms a system of differential equations into a system of (lowerdimensional) integrodifferential equations, so that the reduction of the number of variables is counterpoised by the increased complexity of the reduced model. Here, we will study the special case of evolution equations which contain an inherent time scale separation; in this case, it is possible, in principle, to obtain a reduced system of differential equations in only the slow variables without memory terms. Such an analysis crucially hinges on knowing in which variables (or, more generally, functions of variables) one can write such a reduced system of slow evolution equations.
Moving averages and subsampling have often been used in simple cases as appropriate functions of variables in which to formulate slow lowerdimensional models [29]. However, if the underlying dynamics are sufficiently nonlinear, such statistics may fail to capture the relevant structures and time scales within the data (see Figure LABEL:fig:schematic_fastslow for a schematic illustration). For wellstudied systems, one often has some a priori knowledge of the appropriate observables (such as phase field variables) with which to formulate the reduced dynamics [7, 42]. However, such observables may not be immediately obvious upon inspection for new complex systems, and so we require an automated approach to construct such slow variables.
Given an explicit system of ordinary differential equations, one can make numerical approximations, such as the quasisteady state approximation [32] or the partial equilibrium approximation [17], to reduce the system dimensionality without introducing memory terms. There has been some recent analytical work on extending and generalizing such ideas to more complex systems of equations [1, 6, 11, 12, 19, 29, 36]. However, in many instances, closed form, analytical models are not given explicitly, but can only be inferred from simulation and/or experimental data. We therefore turn to datadriven techniques to analyze such systems and uncover the relevant dynamical modes. In particular, we will use a manifoldlearning based approach, as such methods can accommodate nonlinear structures in highdimensional data.
The core of most manifold learning methods is having a notion of similarity between data points, usually through a distance metric [2, 9, 10, 30, 40]. The distances are then integrated into a global parametrization of the data, typically through the solution of an eigenproblem. In this paper, we will analyze multiple time scale stochastic dynamical systems using datadriven methods. Standard “offtheshelf” manifold learning techniques which utilize the Euclidean distance are not appropriate for analyzing data from such multiscale systems, since this metric does not account for the disparate time scales. Research efforts have addressed the construction of more informative distance metrics, which are less sensitive to noise and can better recover the true underlying structure in the data by suppressing unimportant sources of variability [4, 18, 31, 33, 43]. The Mahalanobis distance is one such metric. It was shown that the Mahalanobis distance can remove the effect of observing the underlying system variables through a complex, nonlinear function [13, 34, 37]. Here, we will show the analogy between removing the effects of such nonlinear observation functions (in the context of data analysis), and reducing a dynamical system to remove the effects of the fast variables. Our approach will build a parametrization of the data which is consistent with the underlying slow variables. Because our approach is datadriven, we require no explicit description of the model, and can extract the underlying slow variables from either simulation or experimental data. Furthermore, the approach implicitly identifies the slow variables within the data and does not require any a priori knowledge of the fast or slow variability sources. Even when the underlying dynamical system is complex with nonlinear coupling between the fast and slow variables, we will show that our approach has the potential to isolate the underlying slow modes.
We will present detailed analysis for our method, and provide conditions under which it will successfully recover the slow variables. Furthermore, based on this analysis, we will present datadriven protocols to tune the parameters of the method appropriately. Our presentation and discussion will address twotimescale stochastic systems; however, we claim that our framework and analysis readily extends to systems with multiple time scale separations.
2 Multiscale Stochastic Systems
Consider the following twotimescale system of stochastic differential equations (SDEs),
\hb@xt@.01(2.1)  
where are independent standard Brownian motions, , and . In the simple case of a linear drift function, i.e., when with , the probability density function of approaches a Gaussian with (finite) variance . The time constant of the approach of the variance to equilibrium is for and for [25]. Thus, the last variables of (LABEL:eq:general_SDE) rapidly approach a local equilibrium measure and exhibit fast dynamics, while the first variables exhibit slow dynamics. A short burst of simulation will yield a cloud of points which is broadly distributed in the fast directions but narrowly distributed in the slow ones. With the appropriate conditions on , the same can be said for more general drift functions, where are the eigenvalues of the Jacobian of [41]. Therefore, (LABEL:eq:general_SDE) defines an dimensional stochastic system with slow variables and fast variables, and defines the time scale separation. The ratio of the powers of in the drift and diffusion terms in (LABEL:eq:general_SDE) is essential, as we require the square of the diffusivity to be of the same order as the drift as [3]. If the diffusivity is larger, then, as , the equilibrium measure will be unbounded. Conversely, if the diffusivity is smaller, the equilibrium measure will go to as .
Assuming the sample average of converges to a distribution which is only a function of the slow variables, then by the averaging principle [16], we can write a reduced SDE in only the slow variables . The aim of our work is to show how we can detect such slow variables automatically from data, in order to help inform modeling efforts and aid in the writing of such reduced stochastic models. In general, we are not given the variables from the original SDE system, but instead, we are given some observations in the form . We assume that , , is a deterministic (possibly nonlinear) function whose image is an dimensional manifold in . For our analysis, we require to be welldefined on , and both and to be continuously differentiable to fourth order. Given data on we would like to recover a parametrization of the data that is onetoone with the slow variables .
3 Local Invariant Metrics
In order to recover the slow variables from data, we will utilize a local metric that collapses the fast directions. Typically, such a metric averages out the fast variables. However, simple averages are inadequate to describe data which is observed through a complicated nonlinear function. Instead, we propose to use the Mahalanobis distance, which measures distances normalized by the respective variances in each local principal direction. Using this metric, we still retain information about both the fast and slow directions and can more clearly observe complex dynamic behavior within the data set.
If two points and are drawn from an dimensional Gaussian distribution with covariance , the Mahalanobis distance between the points is defined as [26]
\hb@xt@.01(3.1) 
In particular, for (LABEL:eq:general_SDE), whose covariance does not depend on , is a constant matrix where
\hb@xt@.01(3.2)  
and the Mahalanobis distance between samples is
\hb@xt@.01(3.3) 
Note that in (LABEL:eq:rescale_x_dist), the fast variables are collapsed and become small, and so this metric is implicitly insensitive to variations in the fast variables. The metric (LABEL:eq:rescale_x_dist) can be rewritten as
\hb@xt@.01(3.4) 
where
\hb@xt@.01(3.5) 
is a stochastic process of the same dimension as , rescaled so that each variable has unit diffusivity. This rescaling transforms our problem from one of detecting the slow variables within dynamic data to one of traditional data mining. The Mahalanobis distance incorporates information about the dynamics and relevant time scales, so that using traditional data mining techniques with this metric will allow us to detect the slow variables in our data [35]. It is important to note that, in practice, we never construct explicitly. It was shown in [34] that, assuming is bilipschitz, the Mahalanobis distance can be extended to approximate (to fourth order) the Euclidean distance between the rescaled samples from accessible ,
\hb@xt@.01(3.6) 
This approximation is accurate when is small. Because we will integrate these distances into a manifold learning algorithm which only considers local distances, we can recover a parametrization of the data which is consistent with the underlying system variables , even when the data are obscured by a function . In Section LABEL:sec:analysis, we will show how we can approximate this distance directly from data .
4 Diffusion Maps for Global Parametrization
From pairwise distances, we want to extract a global parametrization of the data that represents the slow variables. We will use diffusion maps [9, 10], a kernelbased manifold learning technique, to extract a global parametrization using the local distances that we described in the previous section. Given data , we first construct the kernel matrix , where
\hb@xt@.01(4.1) 
Here, denotes the appropriate norm (in our case, the Mahalanobis distance), and is the kernel scale and denotes a characteristic distance within the data set. Note that induces a notion of locality: if , then is negligible. Therefore, we only need our metric to be informative within a ball of radius . We then construct the diagonal matrix , with
\hb@xt@.01(4.2) 
We compute the eigenvalues and eigenvectors of the matrix , and order them such that . is the trivial eigenvector; the next few eigenvectors provide embedding coordinates for the data, so that , the entry of , provides the embedding coordinate for (modulo higher harmonics which characterize the same direction in the data [15]).
5 Estimation of the Mahalanobis Distance
As previously mentioned, we do not have access to the original variables from the underlying original SDE system. Instead, we only have measurements , and we want to estimate the Mahalanobis distance between the variables from observations . The traditional Mahalanobis distance is defined for a fixed distribution, whereas here we are dealing with a distribution that possibly changes as a function of position due to nonlinearities in the observation function and in the drift . Consequently, we use the following modified definition for the Mahalanobis distance between two points,
\hb@xt@.01(5.1) 
where is the covariance of the observed stochastic process at the point , and denotes the MoorePenrose pseudoinverse (since may exceed ).
To motivate this definition of the Mahalanobis distance, we first consider the simple linear case where , with . The covariance of the observed stochastic process is given by . Let be the singular value decomposition (SVD) of , where , , and . The pseudoinverse of the covariance matrix is . Consequently, the Mahalanobis distance (LABEL:eq:mahalanobis_distance) is reduced to
\hb@xt@.01(5.2)  
Hence evaluating the Mahalanobis distances of the observations using (LABEL:eq:mahalanobis_distance) allows us to estimate the Euclidean distances of the rescaled variables (in which the fast coordinates are collapsed).
Following [34], we will show via Taylor expansion that the Mahalanobis distance between the observations (LABEL:eq:mahalanobis_distance) approximates the Euclidean distance in the rescaled variables for general nonlinear observation functions (provided is bilipschitz and both and are differentiable to fourth order). (LABEL:eq:mahalanobis_distance) cannot be evaluated directly since we do not have access to the covariance matrices, so we will instead estimate the covariances directly from data. We can estimate the covariance empirically from a set of values drawn from the local distribution at . One way to obtain such a set of points is to run simulations for a short time, , each starting from . Alternatively, we can consider a single time series of length starting from , and then estimate the covariance from the increments . Although we will present analysis and results for the first type of estimation, the second case is often more practical in practice.
Errors in our estimation of the Mahalanobis distance arise from three sources. One source of error is approximating the function locally as a linear function by truncating the Taylor expansion of at first order. An additional source of error arises from disregarding the drift in the stochastic process, and assuming that samples are drawn from a Gaussian distribution. The third source comes from finite sampling effects. In this work, we will address and discuss the first two sources of error (the finite sampling effects are the subject of future research). We can control the effects of the errors due to truncation of the Taylor expansion by adjusting ; the higherorder terms in this expansion will be small for points which are close, such that adjusting will allow us to only consider distances which are sufficiently accurate in our overall computation scheme. Furthermore, we can control the errors incurred by disregarding the drift by adjusting the time scale of our simulation bursts . Figure LABEL:fig:schematic illustrates some of the issues in choosing the sizes (or if the alternate method is used) and the parameter . We will present both analytical results for the error bounds, as well as an empirical methodology to set the parameters and for our method to accurately recover the slow variable(s).
5.1 Error due to the observation function
We want to relate the distance in the rescaled space, , to the estimated Mahalanobis distance between the observations . We define the error incurred by using the Mahalanobis distance to approximate the true distance as
\hb@xt@.01(5.3) 
By Taylor expansion of around and and averaging the two expansions, we obtain
\hb@xt@.01(5.4)  
where
\hb@xt@.01(5.5)  
In [34], it was shown that the error incurred by using the Mahalanobis distance to approximate the distance between points is (see the Supplementary Materials for details). We now see from (LABEL:eq:mahanaobis_error) that the error is an explicit function of the second and higherorder derivatives of and the distance between samples . We would like to note that this error does not depend on the dynamics of the underlying stochastic process (as we assume the covariances at each point on the manifold are known), but is only a function of the measurement function . The parameter in the diffusion maps calculation determines how much contributes to the overall analysis. From (LABEL:eq:dmaps_kernel), distances which are much greater than are negligible in the diffusion maps computation because of the exponential kernel. Therefore, we want to choose on the order of in a regime where . This is illustrated in Figure LABEL:fig:schematic, where we want to choose small enough so that the curvature and other nonlinear effects (captured in the error term ) are negligible. This will ensure that the errors in the Mahalanobis distance approximation do not greatly effect our overall analysis.
On first inspection, it would appear that our analysis indicates that should be chosen arbitrarily small. However, to obtain a meaningful parametrization of the data set, there must be a nonnegligible number of data points within a ball of radius around each sample. Therefore, the sampling density on the underlying manifold provides a lower bound for .
5.2 Error due to the dynamics
To compute the Mahalanobis distance in (LABEL:eq:mahalanobis_distance), we require , the covariance of the observed stochastic process . We will use simulation bursts to locally explore the dynamics on the manifold of observations in order to estimate the covariance at a point from data [39, 38]. We write the elements of the estimated covariance as
\hb@xt@.01(5.6) 
where is the length of the simulation burst.
Due to the drift in the stochastic process and the (perhaps nonlinear) measurement function , we incur some error by approximating the covariance at a point using simulations of length . Define the error in this approximation as
\hb@xt@.01(5.7) 
By ItôTaylor expansion of and [24],
\hb@xt@.01(5.8)  
where
\hb@xt@.01(5.9)  
From (LABEL:eq:cov_error), the error in the covariance is (as the integrals are each and the integrals are each ) and a function of the derivatives of the observation function and the drift . We want to set such that (this is illustrated in Figure LABEL:fig:schematic), so that the estimated covariances are accurate. Note that in practice, we compute by running many simulations of length starting from , and use the sample average to approximate the expected values in (LABEL:eq:estimated_cov_expected_value). We therefore incur additional error due to finite sampling; this error is ignored for the purposes of this analysis, and quantifying this error is the subject of future research.
Our analysis reveals that the errors decrease with decreasing ; at first inspection, one would want to set arbitrarily small to obtain the highest accuracy possible. However, often in practice, one cannot obtain an arbitrarily refined sampling rate, such that a smaller results in fewer samples with which to approximate the local covariance.When also accounting for these finite sampling errors, and one should take as long as possible while still maintaining negligable errors from the observation function and the drift .
6 Illustrative Examples
For illustrative purposes, we consider the following twodimensional SDE
\hb@xt@.01(6.1)  
where is an constant, as a specific example of (LABEL:eq:general_SDE). is the slow variable, and is a fast noise whose equilibrium measure is bounded and . Figure LABEL:fig:initial_data shows data simulated from this SDE colored by time. We would like to recover a parametrization of this data which is onetoone with the slow variable .
6.1 Linear function
In the first example, our observation function will be the identity function,
\hb@xt@.01(6.2)  
where the fast and slow variables remain uncoupled. In this case, there is no error incurred due to the measurement function (), as the second and higherorder derivatives of are identically 0.
6.1.1 Importance of using the Mahalanobis distance
We want to demonstrate the utility of using the Mahalanobis distance compared to the typical Euclidean distance. We compute the diffusion map embedding for the data in Figure LABEL:fig:initial_data, using both the standard Euclidean distance and the Mahalanobis distance for the computation of the kernel in (LABEL:eq:dmaps_kernel). The data, colored by using the two different metrics, are shown in Figure LABEL:fig:NIV_versus_DMAPS. When using the standard Euclidean distance which does not account for the underlying dynamics, the first diffusion maps recovers the fast variable , suggesting the fast modes is the dominant scale purely in terms of data analysis (Figure LABEL:fig:NIV_versus_DMAPS(a)). In contrast, the slow variable is recovered when using the Mahalanobis distance, as the coloring in Figure LABEL:fig:NIV_versus_DMAPS(b) (where the data are colored by the first diffusion maps variable) is consistent with the coloring in Figure LABEL:fig:initial_data (where the data are colored by time).
6.1.2 Errors in covariance estimation
For the example in (LABEL:eq:linear_transform), the analytical covariance is
\hb@xt@.01(6.3) 
From (LABEL:eq:cov_error), we find
\hb@xt@.01(6.4) 
Therefore, and (provided ; this will be discussed further in Section LABEL:subsec:fastvar). These terms are shown in Figure LABEL:fig:cov_error(a) as a function of . We want to choose in a regime where (the yellow shaded region in Figure LABEL:fig:cov_error indicates where ), so that the errors in the estimated covariance are small with respect to the covariance.
When we do not analytically know the functions or , we can find such a regime empirically by estimating the covariance for several values of . This provides an estimate of as a function of . From Figure LABEL:fig:cov_error(a), we expect a “knee” in the plot of versus when becomes larger than . Figure LABEL:fig:cov_error(b) shows the empirical as a function of for the data in Figure LABEL:fig:initial_data, and the knee in this curve is consistent with the intersection in Figure LABEL:fig:cov_error(a).
6.1.3 Recovery of the fast variable
Note that, for the example in (LABEL:eq:linear_transform), is a constant diagonal matrix. Therefore, taking too large will not lead to nonlinear effects or mixing of the fast and slow variables. Rather, changing will only affect the perceived ratio of the fast and slow timescales.
To see this behavior in our diffusion maps results, we must first discuss the interpretation of the diffusion maps eigenspectrum. The diffusion maps eigenvectors provide embedding coordinates for the data, and the corresponding eigenvalues provide a measure of the importance of each coordinate. However, some eigenvectors can be harmonics of previous eigenvectors; for example, for a data set parameterized by a variable , both and will appear as diffusion maps eigenvectors (see [15] for a more detailed discussion). These harmonics do not capture any new direction within the data set, but do appear as additional eigenvector/eigenvalue pairs. Therefore, for the twodimensional data considered here, the fast variable will not necessarily appear as the second (nontrivial) eigenvector. As the time scale separation increases, the relative importance of the slow and fast directions will also increase. This implies that the eigenvalue corresponding to the eigenvector which parameterizes the fast direction will decrease, and the number of harmonics of the slow mode which appear before the fast mode will increase.
Figure LABEL:fig:recover_fast shows results for three different values of (the corresponding values are indicated by the dashed lines in Figure LABEL:fig:cov_error). When the time scale of the simulation burst used to estimate the local covariance (indicated by the red clouds in the top row of figures), is sufficiently shorter than that of the equilibration time of the fast variable, the estimated local covariance is accurate and the fast variable is collapsed significantly relative to the slow variable. This means that the fast variable is recovered very far down in the diffusion maps eigenvectors. The left two columns of Figure LABEL:fig:recover_fast show that, for this example, when the simulation burst is shorter than the equilibration time, the fast variable is recovered as . However, if the time scale of the burst is longer than the saturation time of the fast variable, the estimated covariance changes: the variance in the slow direction continues to grow, while the variance in the fast direction is fixed. This means that the apparent time scale separation is smaller, the collapse of the fast variable is less pronounced relative to the slow variable, and the fast variable is recovered in an earlier eigenvector (in our ordering of the spectrum). The right column of Figure LABEL:fig:recover_fast shows that, when the burst is now longer than the equilibration time, the fast variable appears earlier in the eigenvalue spectrum and is recovered as .
6.2 Nonlinear observation function
In the second example, our data will be warped into “halfmoon” shapes via the function
\hb@xt@.01(6.5)  
Figure LABEL:fig:initial_data_nonlinear shows the data from Figure LABEL:fig:initial_data transformed by the function in (LABEL:eq:nonlinear_function) and colored by time. It is important to note that this is a difficult class of problem in practice, as none of the observed variables are purely fast or slow, and the observed system appears, at first inspection, to possess no separation of time scales. For this example, the analytical covariance and inverse covariance are
\hb@xt@.01(6.6)  
The fast and slow variables are now coupled through the function , and the Euclidean distance is not informative about the fast or the slow variables. We need to use the Mahalanobis distance to obtain a parametrization that is consistent with the underlying fastslow dynamics.
6.2.1 Errors in Mahalanobis distance
We can bound the Mahalanobis distance by the eigenvalues of ,
\hb@xt@.01(6.7) 
where are the two eigenvalues of . Therefore, for the example in (LABEL:eq:nonlinear_function), we have
\hb@xt@.01(6.8) 
Figure LABEL:fig:cov_error_nonlinear(a) shows and as a function of . The Mahalanobis distance is an accurate approximation to the true intrinsic distance when (the shaded yellow region in the plot indicates where ). We want to choose in a regime where , so that the distances we utilize in the diffusion maps calculation are accurate. We can find such a regime empirically by plotting as a function of , and assessing when the relationship deviates from quadratic. This is shown in Figure LABEL:fig:cov_error_nonlinear(b), and the deviation from quadratic behavior is consistent with the intersection of the analytical expressions plotted in Figure LABEL:fig:cov_error_nonlinear(a).
Figures LABEL:fig:colored_data_nonlinear_cases(a) and (b) show the data from Figure LABEL:fig:initial_data_nonlinear, colored by for two different values of . The corresponding values of are indicated by the dashed lines. When corresponds to a region where , is well correlated with the slow variable. However, when corresponds to a region where , the slow variable is no longer recovered.
6.2.2 Errors in covariance estimation
From (LABEL:eq:cov_error), we find that, for the example in (LABEL:eq:nonlinear_function),
\hb@xt@.01(6.9)  
The error in the covariance is . As expected, the error grows with increasing . We can also see the explicit dependence of the covariance error on the time scale separation ; larger time scale separation results in a larger covariance error, as a more refined simulation burst is required to estimate the covariance of the fast directions. and are plotted as a function of in Figure LABEL:fig:cov_error_nonlinear(c); the shaded yellow portion denotes the region where . As in the previous example, we can empirically find where by plotting as a function of and looking for a knee in the plot. These results are shown in Figure LABEL:fig:cov_error_nonlinear(d).
Figures LABEL:fig:colored_data_nonlinear_cases(a) and (c) show the data from Figure LABEL:fig:initial_data_nonlinear, colored by for two different values of . The corresponding values of are indicated by the dashed lines in Figure LABEL:fig:cov_error_nonlinear(c) and (d). When corresponds to a region where , the slow variable is recovered by the first diffusion maps coordinate. However, when corresponds to a region where , the slow variable is no longer recovered.



7 Conclusions
We have presented a methodology to compute a parametrization of a data set which respects the slow variables in the underlying dynamical system. The approach utilizes diffusion maps, a kernelbased manifold learning technique, with the Mahalanobis distance as the metric. We showed that the Mahalanobis distance collapses the fast directions within a data set, allowing for successful recovery of the slow variables. Furthermore, we showed how to estimate the covariances (required for the Mahalanobis distance) directly from data. A key point in our approach is that the embedding coordinates we compute are not only insensitive to the fast variables, but are also invariant to nonlinear observation functions. Therefore, the approach can be used for data fusion: data collected from the same system via different measurement functions can be combined and merged into a single coordinate system.
In the examples presented, the initial data came from a single trajectory of a dynamical system, and the local covariance at each point in the trajectory was estimated using brief simulation bursts. However, the initial data need not be collected from a single trajectory, and other sampling schemes could be employed. Brief time series are required to estimate the local covariances, but given a simulator, one could reinitialize brief simulation bursts which are sufficiently short and refined from each sample point.
In our examples, we controlled the time scale of sampling and could therefore set the time scale over which to estimate the covariance and the simulation time step arbitrarily small. However, in some settings, such as previously collected historical data, it is not uncommon to have a fixed sampling rate and be unable to reinitialize simulations. In such cases, it is possible that we cannot find an appropriate kernel scale given the fixed such that we can accurately recover the slow variables. For these cases, the data cannot be processed as given, and it is necessary to construct intermediate observers, such as histograms, Fourier coefficients, or scattering transform coefficients [27, 38, 39]. Such intermediates are more complex statistical functions than simple averages and can capture additional structure within the data. They also reduce the effects of noise and permit a larger time step. However, constructing such intermediates often requires additional a priori knowledge about the system dynamics and noise structure.
Clearly, in our analysis, we have ignored the finite sampling effects in our estimation. In reality, both the number of samples used to estimate the covariances, as well as the density of sampled points on the manifold, affect the recovered parametrization and provide additional constraints on and . Future work involves extending our analysis to the finite sample case, and providing guidelines for the amount of data required to apply our methodology.
The methods presented here provide a bridge between traditional data mining and multiple time scale dynamical systems. With this interface established, one can now consider using such datadriven methodologies to extract reduced models (either explicitly, or implicitly via an equationfree framework [14, 21, 22, 23]) which also respect the underlying slow dynamics and geometry of the data. Such reduced models hold the promise of accelerated analysis and reduced simulation of dynamical systems whose effective dynamics are obscure upon simple inspection.
References
 [1] Yacine AïtSahalia, Closedform likelihood expansions for multivariate diffusions, The Annals of Statistics, (2008), pp. 906–937.
 [2] Mikhail Belkin and Partha Niyogi, Laplacian eigenmaps for dimensionality reduction and data representation, Neural Computation, 15 (2003), pp. 1373–1396.
 [3] Nils Berglund and Barbara Gentz, Geometric singular perturbation theory for stochastic differential equations, Journal of Differential Equations, 191 (2003), pp. 1–54.
 [4] Tyrus Berry, John Robert Cressman, Zrinka GreguricÌFerencÌek, and Timothy Sauer, Timescale separation from diffusionmapped delay coordinates, SIAM Journal on Applied Dynamical Systems, 12 (2013), pp. 618–649.
 [5] JJ Brey, R Zwanzig, and JR Dorfman, Nonlinear transport equations in statistical mechanics, Physica A: Statistical Mechanics and its Applications, 109 (1981), pp. 425–444.
 [6] Christopher P Calderon, Fitting effective diffusion models to data associated with a “glassy” potential: Estimation, classical inference procedures, and some heuristics, Multiscale Modeling & Simulation, 6 (2007), pp. 656–687.
 [7] LongQing Chen, Phasefield models for microstructure evolution, Annual Review of Materials Research, 32 (2002), pp. 113–140.
 [8] Alexandre J Chorin, Ole H Hald, and Raz Kupferman, Optimal prediction and the Mori–Zwanzig representation of irreversible processes, Proceedings of the National Academy of Sciences, 97 (2000), pp. 2968–2973.
 [9] Ronald R. Coifman and Stephane Lafon, Diffusion maps, Applied and Computational Harmonic Analysis, 21 (2006), pp. 5–30.
 [10] Ronald R Coifman, Stephane Lafon, Ann B Lee, Mauro Maggioni, Boaz Nadler, Frederick Warner, and Steven W Zucker, Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps, Proceedings of the National Academy of Sciences, 102 (2005), pp. 7426–7431.
 [11] MarieNathalie ContouCarrere, Vassilios Sotiropoulos, Yiannis N Kaznessis, and Prodromos Daoutidis, Model reduction of multiscale chemical Langevin equations, Systems & Control Letters, 60 (2011), pp. 75–86.
 [12] Guang Qiang Dong, Luke Jakobowski, Marco AJ Iafolla, and David R McMillen, Simplification of stochastic chemical reaction models with fast and slow dynamics, Journal of Biological Physics, 33 (2007), pp. 67–95.
 [13] Carmeline J Dsilva, Ronen Talmon, Neta Rabin, Ronald R Coifman, and Ioannis G Kevrekidis, Nonlinear intrinsic variables and state reconstruction in multiscale simulations, The Journal of Chemical Physics, 139 (2013), p. 184109.
 [14] Radek Erban, Ioannis G Kevrekidis, David Adalsteinsson, and Timothy C Elston, Gene regulatory networks: A coarsegrained, equationfree approach to multiscale computation, The Journal of Chemical Physics, 124 (2006), p. 084106.
 [15] Andrew L Ferguson, Athanassios Z Panagiotopoulos, Pablo G Debenedetti, and Ioannis G Kevrekidis, Systematic determination of order parameters for chain dynamics using diffusion maps, Proceedings of the National Academy of Sciences, 107 (2010), pp. 13597–13602.
 [16] Mark I Freidlin, Joseph Szücs, and Alexander D Wentzell, Random perturbations of dynamical systems, vol. 260, Springer, 2012.
 [17] Paul M Gallagher, Amulya L Athayde, and Cornelius F Ivory, The combined flux technique for diffusion–reaction problems in partial equilibrium: Application to the facilitated transport of carbon dioxide in aqueous bicarbonate solutions, Chemical Engineering Science, 41 (1986), pp. 567–578.
 [18] Shai Gepshtein and Yosi Keller, Image completion by diffusion maps and spectral relaxation, IEEE Transactions on Image Processing, 22 (2013), pp. 2983–2994.
 [19] Dror Givon, Raz Kupferman, and Andrew Stuart, Extracting macroscopic dynamics: Model problems and algorithms, Nonlinearity, 17 (2004), p. R55.
 [20] Carmen Hijón, Pep Español, Eric VandenEijnden, and Rafael DelgadoBuscalioni, Mori–Zwanzig formalism as a practical computational tool, Faraday Discussions, 144 (2010), pp. 301–322.
 [21] Ioannis G Kevrekidis, C William Gear, and Gerhard Hummer, Equationfree: The computeraided analysis of complex multiscale systems, AIChE Journal, 50 (2004), pp. 1346–1355.
 [22] Ioannis G Kevrekidis, C William Gear, James M Hyman, Panagiotis G Kevrekidis, Olof Runborg, Constantinos Theodoropoulos, et al., Equationfree, coarsegrained multiscale computation: Enabling microscopic simulators to perform systemlevel analysis, Communications in Mathematical Sciences, 1 (2003), pp. 715–762.
 [23] Ioannis G Kevrekidis and Giovanni Samaey, Equationfree multiscale computation: Algorithms and applications, Annual Review of Physical Chemistry, 60 (2009), pp. 321–344.
 [24] Peter E Kloeden and Eckhard Platen, Numerical Solution of Stochastic Differential Equations, vol. 23, Springer, 1992.
 [25] Tony Lelièvre, Francis Nier, and Grigorios A Pavliotis, Optimal nonreversible linear drift for the convergence to equilibrium of a diffusion, Journal of Statistical Physics, 152 (2013), pp. 237–274.
 [26] Prasanta Chandra Mahalanobis, On the generalized distance in statistics, Proceedings of the National Institute of Sciences (Calcutta), 2 (1936), pp. 49–55.
 [27] Stéphane Mallat, Group invariant scattering, Communications on Pure and Applied Mathematics, 65 (2012), pp. 1331–1398.
 [28] Hazime Mori, Transport, collective motion, and Brownian motion, Progress of Theoretical Physics, 33 (1965), pp. 423–455.
 [29] GA Pavliotis and AM Stuart, Parameter estimation for multiscale diffusions, Journal of Statistical Physics, 127 (2007), pp. 741–781.
 [30] Sam T Roweis and Lawrence K Saul, Nonlinear dimensionality reduction by locally linear embedding, Science, 290 (2000), pp. 2323–2326.
 [31] Yossi Rubner, Carlo Tomasi, and Leonidas J Guibas, The earth mover’s distance as a metric for image retrieval, International Journal of Computer Vision, 40 (2000), pp. 99–121.
 [32] Lee A Segel and Marshall Slemrod, The quasisteadystate assumption: A case study in perturbation, SIAM review, 31 (1989), pp. 446–477.
 [33] Karen Simonyan, Omkar M Parkhi, Andrea Vedaldi, and Andrew Zisserman, Fisher vector faces in the wild, in Proc. BMVC, vol. 1, 2013, p. 7.
 [34] Amit Singer and Ronald R Coifman, Nonlinear independent component analysis with diffusion maps, Applied and Computational Harmonic Analysis, 25 (2008), pp. 226–239.
 [35] Amit Singer, Radek Erban, Ioannis G Kevrekidis, and Ronald R Coifman, Detecting intrinsic slow variables in stochastic dynamical systems by anisotropic diffusion maps, Proceedings of the National Academy of Sciences, 106 (2009), pp. 16090–16095.
 [36] Vassilios Sotiropoulos, MarrieNathalie ContouCarrere, Prodromos Daoutidis, and Yiannis N Kaznessis, Model reduction of multiscale chemical Langevin equations: A numerical case study, IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB), 6 (2009), pp. 470–482.
 [37] Ronen Talmon and Ronald R Coifman, Empirical intrinsic geometry for nonlinear modeling and time series filtering, Proceedings of the National Academy of Sciences, 110 (2013), pp. 12535–12540.
 [38] , Intrinsic modeling of stochastic dynamical systems using empirical geometry, Applied and Computational Harmonic Analysis, (2014).
 [39] Ronen Talmon, Stéphane Mallat, Hitten Zaveri, and Ronald R Coifman, Manifold learning for latent variable inference in dynamical systems, tech. report, submitted, 2014, Tech. Report YALEU/DCS/TR1491, 2014.
 [40] Joshua B Tenenbaum, Vin De Silva, and John C Langford, A global geometric framework for nonlinear dimensionality reduction, Science, 290 (2000), pp. 2319–2323.
 [41] Cédric Villani, Hypocoercivity, American Mathematical Soc., 2009.
 [42] AA Wheeler, WJ Boettinger, and GB McFadden, Phasefield model for isothermal phase transitions in binary alloys, Physical Review A, 45 (1992), p. 7424.
 [43] Eric P Xing, Michael I Jordan, Stuart Russell, and Andrew Y Ng, Distance metric learning with application to clustering with sideinformation, in Advances in Neural Information Processing Systems, 2002, pp. 505–512.
 [44] Robert Zwanzig, Memory effects in irreversible thermodynamics, Physical Review, 124 (1961), p. 983.