Sparse model from optimal nonuniform embedding of time series
Abstract
An approach to obtaining a parsimonious polynomial model from time series is proposed. An optimal minimal nonuniform time series embedding schema is used to obtain a time delay kernel. This scheme recursively optimizes an objective functional that eliminates maximum number of false nearest neighbors between successive state space reconstruction cycles. A polynomial basis is then constructed from this time delay kernel. A sparse model from this polynomial basis is obtained by solving a regularized least squares problem. The constraint satisfaction problem is made computationally tractable by keeping the ratio between the number of constraints to the number of variables small by using fewer samples spanning all regions of the reconstructed state space. This helps the structure selection process from an exponentially large combinatorial search space. A forward stagewise algorithm is then used for fast discovery of the optimization path. Results are presented for the MackeyGlass system.
pacs:
05.45.Tp, 89.70.Eg, 89.75.k, 07.05.Tp, 05.45.Pq, 89.65.Gh, 05.10.a, 05.10.Gg, 05.40.aModeling and forecasting has been a part of the human quest for knowledge since time immemorial. At the heart of all time series modeling methods is the notion of causality. Forecasting relies on the ability to infer the future by observing a system’s past. Modeling involves finding the rules that govern the evolution of the states of the system. Discovering and quantifying these rules lie at the heart of all physical sciences. Over the past decade a new kind of scienceWolfram () has emerged that seeks to identify hidden patterns in apparently complex systems. Complexity is no longer a paranoia with the seminal paper by Watts and StrogatzWatts () elaborating how collective dynamics may emerge from a system with a high number of participating entities. A complex time series has an inherent geometry and Packard et alPackard () were first to show that a representative geometry of a dynamical system can be obtained by using a time series of one of its observables. TakensTakens () later showed that the geometrical reconstruction methodology proposed by Packard et alPackard () is actually an embedding of the original state space of the system. The embedding theorem says that the dynamical attractor obtained using time delayed version of an observable is diffeomorphic to the original state space.
This paper builds a global model for the system from an observed time series on a polynomial basis derived using optimal time delay coordinates. Consider a time series of observation of a system . Let the dynamics of the system be described by the following model:
(1) 
is the embedding dimension of the system and are time delay coordinates. The delays are determined using a procedure that was shown to be both optimal and minimalNichkawde (). The function in general is a nonlinear function of its inputs. Polynomials provide good approximation to with proper structure selectionletellier2009frequently (). Polynomials have been successfully used to model sun spot time seriesaguirre2008forecasting (), Canadian lynx cyclesmaquet2007global (), dynamics of mixing reactorletellier1997recovering (), sleep apnea dynamicsaguirre2004stability () and experimental data from Chua’s circuitaguirre1997nonlinear (). The model accuracy may be improved by considering polynomials of higher degree. However, parsimony principle states that the best model is not just that which minimizes error on the observed time series, but should also do so by having a minimum number of terms. Thus, for the same amount of modeling error on the observed data a sparse model with fewer terms is preferred. Overparameterized models lead to overfitting. The notion of sparsity in parameter space was also shown to be the basis of emergent theories in a recent paperMachta01112013 (). There are two opposing requirements of increasing accuracy: (1) including polynomial terms of higher degree and (2) at the same time ensuring sparsity. The number of terms in the polynomial expansion of grows with degree of the polynomial. One must consider terms for a polynomial of degree with variables. Thus, 20 terms have to be considered for a polynomial of degree 2 with 5 variables, where as for a polynomial of degree 10 with 5 variables 3002 terms will have to be considered. Building a good polynomial model is about structure selection which is a combinatorial optimization problem with search space exponentially large. Modeling a system with high embedding dimension imposes great challenge when a polynomial basis of higher degree is considered, especially when the sample size for modeling is small.
Consider a function of variables .
(2) 
would equal to the future value for an autoregressive function. The example taken here considers a general function approximation. Let be linear and sparse
(3) 
If is sparse most of the coefficients would be zero. Let us collect samples of where is order . The objective is to model using these observations. Let be the predicted value of the model and let be the true value for the data point. and would be the one step ahead predicted and true values respectively for an autoregressive model.
(4) 
The error , as given by Eq. (4), in the modeling ought to minimized. The usual methodology of linear regression is orthogonal least squares (OLS) which uses the GramSchmidt procedure to determine the coefficients s such that is minimized. As noted previously, the best model is not just the one that minimizes error but must do so with minimal complexity. OLS has been used in Ref. aguirre2008forecasting (); maquet2007global (); letellier1997recovering (); aguirre2004stability (); aguirre1997nonlinear () wherein structure selection has been done using error reduction ratio criterion. Another method for structure selection involves use of Akaike information criterionbarahona1996detection (). In this procedure, the following information criterion is minimized in accordance with the parsimony principle:
where is the number of the chosen term in the polynomial model and is representative of a particular combination of polynomial terms, and,
(5) 
is the mean value of or
(6) 
These procedures are obviously time consuming and involve an exhaustive search over exponentially large combinatorial space.
In this paper, learning and structure selection is combined into an efficient single step procedure. A model complexity term is added to the objective functional in Eq. (4) that penalizes the norm of the model parameters. Such a form of linear regression is termed as least absolute shrinkage and selection operator(LASSO)tibshirani1996regression () with the objective functional being of the following form
(7) 
where the hyperparamter . LASSO is known to recover a sparse modeltibshirani1996regression (). The power of based regularization in sparse model recovery is shown in Fig. 1. Let in Eq. (3) be independent and identically distributed random variables. Let and only let 10 of the s be nonzero. 50 samples of and corresponding s are collected. 1% of noise is introduced in . based regularization is able to recover all 10 nonzero coefficients almost perfectly compared to OLS where all coefficients are erroneously found to be nonzero.
In this paper, the power of based regularization is exploited in sparse model recovery to construct parsimonious polynomial model for a chaotic system using measured time series. In the next Section, a probabilistic interpretation for the proposed objective function is provided. In Section II, a constraint satisfaction perspective to the problem is given. In Section III, the solution procedure is explained. Section IV explores the model for a highdimensional version of the MackeyGlass system. Section V concludes the paper.
I Probabilistic interpretation
From a Bayesian perspective, the term in Eq. (7) imposes a Laplacian prior on the parameter space. The Laplace distribution has a sharp probability peak at zero as shown in Fig. 2 for a Laplace distribution with . The prior guess about most parameters being zero is thus imposed in this form.
Bayes rule states that the posterior probability of a model is proportional to the product of likelihood and prior
(8) 
The above is the data, which for our case are values in dimensional space of regressors . is the vector of coefficients in Eq. (3). The residual can be assumed to be Gaussian distributed with zero mean. The likelihood can thus be defined as
(9) 
The Laplacian prior imposed on takes the form
(10) 
Thus, the negative log of the posterior in Eq. (8) can be written as
(11) 
The right hand side of Eq. (11) is exactly equal to the right hand side of Eq. (7). Since logarithm is a monotonic function, the solution to LASSO maximizes the posterior probability of a model with a Laplacian prior. One possible solution strategy could be to use Markov chain Monte Carlo. Such a random walk is designed to be ergodic on the posterior distribution of . Due to ergodicity, it draws more samples from near the hills of the landscape, rather than the ravines. Consequently, the samples drawn are representative of the probability distribution. The values of can be drawn from gamma distribution which is often used in Bayesian inference for width of a distribution. Expectation can be taken on this sample in order to determine . This approach frees us from making a choice of . However, since this is sampling based, it is a slow solution procedure. Moreover, a unique solution is not guaranteed on multiple runs on the same data. A much faster deterministic algorithm is used which will be described later.
Ii Constraint satisfaction interpretation
This can also be seen from a constraint satisfaction problem perspective. constraints are imposed via Eq. (3) on variables . The problem can be viewed as minimization of subject to where is the value of regressor at the point and is a small positive number. Computational tractability of constraint satisfaction problems has been studied in various other contexts such as coloring of a random graphzdeborova2007phase (). Notionally, it is harder to satisfy a larger number of constraints simultaneously and there occurs a phase transition in computational landscape beyond which an efficient solution does not exist. Such problems belong to NP complete classercsey2011optimization (). Thus, the ratio of the number of constraints and the number of variables , , should be kept small so that the optimization task is computationally tractable. Large values of result in NP complete optimization problemsmezard2002analytic (); ercsey2011optimization (). This means that an optimal solution can only be found in exponential time or . The computational landscape analysis performed in the thermodynamic limit of and is borrowed in devising an efficient solution procedure.
Iii Solution
iii.1 Optimal Minimal Embedding
Takens embedding theorem allows us to reconstruct the representative state space of a dynamical systemTakens (). The hallmark of this result is that a state space can be reconstructed using only one or a few of the observables. The first time delay is chosen as the first minimum of mutual information between observed variable and its timedelayed counterpart in order to determine the optimal time delay Mutual (). This function is given by
(12) 
where and . The subsequent delays are chosen in an ad hoc manner as multiples of (). This approach is termed as uniform as the chosen time lags are uniformly separated. The uniform time delay approach is however not always optimal for modelingaguirre2009modeling (). Investigations have also shown that nonuniform embedding performs better than uniform embedding for causality and coupling detectionVlachos (); Luca (); kugiumtzis2013direct (). The embedding dimension above is determined by false nearest neighbors approachKennel (). In this approach, a nearest neighbor is flagged as false, if the same two points become distant in a higher embedding dimension. is thus chosen when there are no false nearest neighbors. Recent investigations have revealed two major caveats with this approachNichkawde ().

That the inability to eliminate false nearest neighbors could be due to a poor choice of the last delay . Please refer to Section III in Ref. Nichkawde ().

The uniform delay approach is not a minimal approach. It is possible to eliminate larger number of false nearest neighbors between successive reconstruction cycles and thus achieve an embedding with a smaller number of delays. Minimality is important in Occam’s razor principle. The minimum description length principlerissanen1978modeling () also stresses minimality by having size or complexity of model a part of the description length.
Another very important concept is irrelevancyCasdagli (). The last few longer delays in the uniform embedding approach may be irrelevant. Causality is lost beyond a certain time in the past.
Addressing the above highlighted issues, an optimal nonuniform embedding of the system is used for modelingNichkawde (). The methodology imposes an upper bound in the search for optimal time delays by using an objective functional that quantifies irrelevancy. The reconstruction methodology is very fast and performs in where is the length of the time series. This methodology recursively chooses delays that maximize derivatives on the project manifold. The objective functional is of the following form:
(13) 
In the above equation, is the absolute value of the directional derivative evaluated in the direction from the to the point of the projected attractor manifold which happens to be the nearest neighbor and stands for the mean value. may also be expressed as follows
(14) 
where be the Euclidean distance of the the point in the reconstruction cycle to its nearest neighbor and is the same distance in the reconstruction cycle with as the additional time delay coordinate. The nearest neighbor in the cycle is deemed as false if:
(15) 
where is a thresholdKennel (). Now using Eq. (14) above can be written as:
(16) 
The right hand side of Eq. (13) is the geometric mean of . It is thus obvious from the above relationship that the recursive maximization of the functional given by Eq. (13) eliminates the largest number of false nearest neighbors between successive reconstruction cycles and thus helps achieve an optimal minimal embeddingNichkawde ().
Fig. 3 shows the loglikelihood of derivatives on projected manifold given by Eq. (13) for the MackeyGlass systemmackey1977oscillation (). MackeyGlass equation is the nonlinear time delay differential equation
(17) 
where represents the value of at time . The equation shows chaos with increasing complexity with an increasing value of . The time series was generated with parameter values and . 2000 points were sampled with of 0.5 after removing the transient. First 1000 points are used for model building. Delays of 33, 15 and 45 shown by round dots are found to be most optimal as shown in Fig. 2(a).
iii.2 Polynomial model from optimal nonuniform delays
Continuing to develop a polynomial model for the MackeyGlass example, a polynomial model of 4 variables and in Eq. (1) is considered. s in Eq. (3) are replaced by a polynomial basis of order of variables and .
(18) 
Terms up to degree 2 have been shown in the equation above. Let there be model terms with corresponding coefficients . Let us collect measurements and plug those in Eq. (18). This system of equations is overdetermined if . Since there are unknowns, only about number of measurements is required. This fact is used and an effective sample of number of measurements is drawn to solve the inverse problem. This is achieved by using means clustering algorithmLloyd () which samples points approximately from all regions of the reconstructed attractor. Given a set of observations , , means clustering aims to partition the observations into sets so that within each cluster, the Euclidean distance of the points from the cluster mean of
is minimized. There are heuristic algorithms which achieve nearly the same result in much faster timesculley2010web (). This accomplishes two simultaneous objectives: 1) using only few representative samples greatly speeds up the computation. 2) Fewer samples imply fewer constraints and therefore this renders the problem more computationally tractable. As noted previously in Section II, from a constraint satisfaction perspective large values in lead to a hard optimization problem. This means that the only possible solution that will work is a sequential search through exponentially large combinatorial search space. For example, a degree polynomial model of MackeyGlass system with 4 variables would have 329 terms and approximately possible structures. Appropriately sampling the reconstructed state space using kmeans clustering keeps the ratio small and brings down the computational complexity to polynomial time or . In this paper, the parameters which minimizes the objective function given by Eq. (7) are learnt by least angle regressionEfron () algorithm which discovers the optimization path in a forward stagewise manner. In this approach, the origin of the residual space (also elsewhere referred to as data spacetranstrum2010nonlinear (); transtrum2011geometry () with manifold formed by various values of parameters as model manifold) is approached by including the regressors one at a time which makes the least angle with the current residue. The optimal value for the parameter is determined using 5fold cross validation. The least angle regression algorithm achieves this in linear time or . The values are rescaled between 0 and 1 before building the model. The proposed approach will work even for small time series as long as the entire range of state space has been sampled.
This is demonstrated with the MackeyGlass time series. 1000 data points are used. These points, however, span almost the entire space of the attractor as shown in Fig. 2(b). The clusters for the reconstructed MackeyGlass attractor and their centroids are also shown in Fig. 2(b). The model building is done by randomly sampling a point from each of the clusters. This model was then used to predict the next 1000 values.
The attractor from the signal generated using the learnt model is shown in Fig. 3(a). Fig. 3(b) shows the predictions for the next 1000 values for models of the MackeyGlass system of with polynomial degrees 3, 5 and 7. The thick solid lines shows the actual values. The model accuracy increases with increased polynomial degree as shown in Fig. 3(c) which shows the normalized mean square error (NMSE) for models of increasing polynomial degree. Model with polynomial degree 7 has the best NMSE which was shown to be 0.07 after 1000 time steps. This model has 64 nonzero terms out of a total of 329 terms.
Iv Highdimensional MackeyGlass system
A highdimensional version of MackeyGlass model is probed. The delay parameter in Eq. (17) is set to 80 and 5000 points are sampled at step size of 0.5. The time series is shown in Fig. 4(a). The KaplanYorke dimension for this system is about 8poon2001titration (). State space reconstruction is performed for this time series. The reconstruction cycles have been shown in Fig. 4(b). Delays of 63, 141, 106, 39, 85, 19, 125 and 157 were obtained with an embedding dimension of 9. A degree polynomial model is built for this series. The model obtained with their terms has been given in Appendix B. The model has a total 89 terms selected out of 11439 terms.
iv.1 Model evaluation using anticipated synchronization
Evaluation of a dynamical model can be done in several different ways. Chaotic attractors have certain geometrical properties such as fractal dimension and Lyapunov exponent which can be compared with that of the original system. In cases, where the available time series is very short and highdimensional, it is not possible to determine these properties from the time series. A different approach is used in this case in order to evaluate the model. An interesting property of matched systems is that they can easily synchronizebrown1994modeling (). This property can be used to evaluate the learnt dynamical model from time series. Let be the true dynamics of the system:
(19) 
Let denote the dynamics of the learnt model. The learnt model can be coupled to the true model as follows
(20) 
If is adequately chosen and is sufficiently close to then . This means that the model will synchronize with the observed values. If and differ, the error will not go to zero but will stay around the origin of the residual space. The average distance to the origin of such a space will depend on the difference between and . Therefore such a distance is a measure of how far the estimated model is from the true dynamics aguirre2006evaluation (). In our case, the true dynamics is unknown. However, the scalar measurement is available. is taken as an appropriate embedding of the scalar observation and is the prediction of the model. A coupling of form where is used. The model learnt for this system is synchronized with another set of observations for the same system. These observations were not used in building the model. Evolution of highdimensional MackeyGlass model given by Eq. (25) coupled to the second set of observations is shown in Fig. 5(a). Since, the model give by Eq. (25) is matched to the time series dynamics, it synchronizes with the external driver.
A particular form of synchronization known as anticipated synchronization can be used to anticipate the future states in advancevoss2001dynamic (). In this form of synchronization, a time delay is introduced in Eq. (20). This introduced memory gives the system
(21) 
The driving system is unaffected and only the past values of driven system are needed. A matched system would synchronize on the manifold given by
(22) 
on which the driven system anticipates driving system’s state . This sort of scheme can also be chained where the output of the first slave in the chain can be made to synchronize with a delayed version of the same model
(23) 
Figure 5(b) shows the anticipation with for two slaves in the chain. Figure 5(c) zooms onto the last 200 values in Fig. 5(b). Slave 1 anticipates 1 time step ahead and slave 2 anticipates 2 time steps ahead. Future states can easily be anticipated in advance using this synchronization scheme. This demonstrates building of a sparse datadriven model for a very high dimensional system.
V Discussion and Conclusions
A new method for developing a polynomial model from complex time series data has been described. A generalized linear model setting has been used with polynomial terms obtained from optimal nonuniform delays as regressors. Penalization of norm of model coefficients leads to recovery of a parsimonious model. The approach combines structure selection and model fitting into a single step procedure by discovering the optimization path in a forward stagewise manner. The solution procedure is guided in terms of both speed and constraint satisfiability by using fewer samples taken from all regions of reconstructed state space. The procedure has been demonstrated on the MackeyGlass system. A highly accurate model of a low dimensional MackeyGlass system was obtained which is able to do long term prediction. The model for highdimensional MackeyGlass system synchronizes very well with the simulated observations. It is interesting to note that the inverse problem for highdimensional MackeyGlass system is underdetermined when a degree polynomial model is considered with an embedding dimension of 9. There are 11439 variables whereas the number of equations is less than this number for the size of the data. The present approach is still able to recover a sparse model. This is because the signal has a sparse representation in a polynomial basis. Thus, for this case, the solution exactly follows the fundamentals of compressed sensingcandes2008introduction () which has emerged as the state of the art technique in signal acquisition. Measurement matrix is chosen to be comprised of random numbers in compressed sensingcandes2008introduction () whereas the ‘measurement matrix’ in this case happens to be values of polynomial regressors. The constraint satisfaction perspective given in Section II is valid only for an overdetermined system (). However, when the system is underdetermined () such as the case for highdimensional MackeyGlass system, all the available data points should be used in building the model. The proposed approach can also be used to detect nonlinear dynamics in time series with the noise titration approachpoon2001titration (). In this approach, nonlinear dynamics is ascertained by presence of nonlinear terms in an autoregressive model. If a linear model performs as well as a nonlinear model, then it is concluded that nonlinear dynamics is not present in the time series. The performance is measured by mean square error over the entire data. As noted previously, it is always possible to reduce such a metric almost indefinitely by overparametrizing. Thus, the robust framework proposed in this paper to build a parsimonious model with small size data would be a great choice for the noise titration approach. Extending this argument, this approach can also be used to detect the presence of colored nonlinearly correlated noisefreitas2009failure (). The dynamics in economics and financial settings are often driven by additive white noise. This noise component is deemed to be conditionally heteroskedastic. This means that they have varying variance whose value depends on few of the past values. This form of modeling is termed as autoregressive conditional heteroskedasticity. The framework proposed in this paper allows us to explore possibly nonlinear conditional heteroskedasticity with perforations. A polynomial approach to nonlinear Granger causalitygranger1969investigating () directly follows from this work. One just needs to check if the predictive power of the model is enhanced by including polynomial terms from the second variable. The methodology can thus detect the causal effect of linear and nonlinear terms composed of the embedding of the second variable.
Appendix A Model for MackeyGlass system
The degree polynomial model for MackeyGlass system with parameters value and :
(24) 
where and .
Appendix B Highdimensional MackeyGlass model
The degree polynomial model for MackeyGlass system with parameters value and :
(25) 
where and .
Acknowledgements.
This research was supported by the Australian Research Council and Sirca Technology Pty Ltd under Linkage Project LP100100312. The author is also supported by International Macquarie University Research Excellence Scholarship (iMQRES). Deb Kane is thanked for her editorial support during preparation of this manuscript.References
 (1) Stephen Wolfram, A New Kind of Science, 1st ed. (Wolfram Media, 2002) ISBN 1579550088
 (2) D. J. Watts and S. H. Strogatz, “Collective dynamics of “smallworld” networks,” Nature 393, 440–442 (1998)
 (3) N. H. Packard, J. P. Crutchfield, J. D. Farmer, and R. S. Shaw, “Geometry from a time series,” Phys. Rev. Lett. 45, 712–716 (Sep 1980), http://link.aps.org/doi/10.1103/PhysRevLett.45.712
 (4) Floris Takens, “Detecting strange attractors in turbulence,” Dynamical Systems and Turbulence, Warwick 1980, Dynamical Systems and Turbulence, Lecture Notes in Mathematics 898, 366–381 (1981), http://dx.doi.org/10.1007/BFb0091924
 (5) Chetan Nichkawde, “Optimal statespace reconstruction using derivatives on projected manifold,” Phys. Rev. E 87, 022905 (Feb 2013), http://link.aps.org/doi/10.1103/PhysRevE.87.022905
 (6) Christophe Letellier, Luis A. Aguirre, and U. S. Freitas, “Frequently asked questions about global modeling,” Chaos: An Interdisciplinary Journal of Nonlinear Science 19, 023103 (2009), http://link.aip.org/link/?CHA/19/023103/1
 (7) LA Aguirre, C Letellier, and J Maquet, “Forecasting the time series of sunspot numbers,” Solar Physics 249, 103–120 (2008)
 (8) J Maquet, C Letellier, and Luis A Aguirre, “Global models from the canadian lynx cycles as a direct evidence for chaos in real ecosystems,” Journal of Mathematical Biology 55, 21–39 (2007)
 (9) C Letellier, L Le Sceller, G Gouesbet, F Lusseyran, A Kemoun, and B Izrar, “Recovering deterministic behavior from experimental time series in mixing reactor,” AIChE Journal 43, 2194–2202 (1997)
 (10) Luis Antonio Aguirre and Álvaro VP Souza, “Stability analysis of sleep apnea time series using identified models: a case study,” Computers in Biology and Medicine 34, 241–257 (2004)
 (11) Luis A Aguirre, Giovani G Rodrigues, and EMAM Mendes, “Nonlinear identification and cluster analysis of chaotic attractors from a real implementation of chua’s circuit,” International Journal of Bifurcation and Chaos 7, 1411–1424 (1997)
 (12) Benjamin B. Machta, Ricky Chachra, Mark K. Transtrum, and James P. Sethna, “Parameter space compression underlies emergent theories and predictive models,” Science 342, 604–607 (2013), http://www.sciencemag.org/content/342/6158/604.full.pdf, http://www.sciencemag.org/content/342/6158/604.abstract
 (13) Mauricio Barahona and ChiSang Poon, “Detection of nonlinear dynamics in short, noisy time series,” Nature 381, 215–217 (1996)
 (14) Robert Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society. Series B (Methodological), 267–288(1996)
 (15) Lenka Zdeborová and Florent Krzakala, “Phase transitions in the coloring of random graphs,” Physical Review E 76, 031131 (2007)
 (16) Mária ErcseyRavasz and Zoltán Toroczkai, “Optimization hardness as transient chaos in an analog approach to constraint satisfaction,” Nature Physics 7, 966–970 (2011)
 (17) Marc Mézard, Giorgio Parisi, and Riccardo Zecchina, ‘‘Analytic and algorithmic solution of random satisfiability problems,” Science 297, 812–815 (2002)
 (18) Andrew M. Fraser and Harry L. Swinney, “Independent coordinates for strange attractors from mutual information,” Phys. Rev. A 33, 1134–1140 (Feb 1986)
 (19) Luis A Aguirre and Christophe Letellier, ‘‘Modeling nonlinear dynamics and chaos: a review,” Mathematical Problems in Engineering 2009 (2009)
 (20) Ioannis Vlachos and Dimitris Kugiumtzis, “Nonuniform statespace reconstruction and coupling detection,” Phys. Rev. E 82, 016207 (Jul 2010)
 (21) Luca Faes, Giandomenico Nollo, and Alberto Porta, “Informationbased detection of nonlinear granger causality in multivariate processes via a nonuniform embedding technique,” Phys. Rev. E 83, 051112 (May 2011)
 (22) D. Kugiumtzis, “Directcoupling information measure from nonuniform embedding,” Phys. Rev. E 87, 062918 (Jun 2013), http://link.aps.org/doi/10.1103/PhysRevE.87.062918
 (23) Matthew B. Kennel, Reggie Brown, and Henry D. I. Abarbanel, “Determining embedding dimension for phasespace reconstruction using a geometrical construction,” Phys. Rev. A 45, 3403–3411 (Mar 1992)
 (24) Jorma Rissanen, “Modeling by shortest data description,” Automatica 14, 465–471 (1978)
 (25) Martin Casdagli, Stephen Eubank, J.Doyne Farmer, and John Gibson, “State space reconstruction in the presence of noise,” Physica D: Nonlinear Phenomena 51, 52 – 98 (1991), ISSN 01672789
 (26) Michael C Mackey, Leon Glass, et al., “Oscillation and chaos in physiological control systems,” Science 197, 287–289 (1977)
 (27) S. Lloyd, “Least squares quantization in pcm,” Information Theory, IEEE Transactions on 28, 129 – 137 (mar 1982), ISSN 00189448
 (28) D Sculley, “Webscale kmeans clustering,” in Proceedings of the 19th international conference on World wide web (ACM, 2010) pp. 1177–1178
 (29) Bradley Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani, “Least angle regression,” Annals of Statistics 32, 407–499 (2004)
 (30) Mark K Transtrum, Benjamin B Machta, and James P Sethna, “Why are nonlinear fits to data so challenging?.” Physical review letters 104, 060201 (2010)
 (31) Mark K Transtrum, Benjamin B Machta, and James P Sethna, “Geometry of nonlinear least squares with applications to sloppy models and optimization,” Physical Review E 83, 036701 (2011)
 (32) ChiSang Poon and Mauricio Barahona, “Titration of chaos with added noise,” Proceedings of the National Academy of Sciences 98, 7107–7112 (2001)
 (33) Reggie Brown, Nikolai F. Rulkov, and Eugene R. Tracy, “Modeling and synchronizing chaotic systems from timeseries data,” Phys. Rev. E 49, 3784–3800 (May 1994), http://link.aps.org/doi/10.1103/PhysRevE.49.3784
 (34) Luis Antonio Aguirre, Edgar Campos Furtado, and Leonardo A. B. Tôrres, ‘‘Evaluation of dynamical models: Dissipative synchronization and other techniques,” Phys. Rev. E 74, 066203 (Dec 2006), http://link.aps.org/doi/10.1103/PhysRevE.74.066203
 (35) Henning U Voss, “Dynamic longterm anticipation of chaotic states,” Physical Review Letters 87, 014102 (2001)
 (36) Emmanuel J Candès and Michael B Wakin, ‘‘An introduction to compressive sampling,” Signal Processing Magazine, IEEE 25, 21–30 (2008)
 (37) Ubiratan S Freitas, Christophe Letellier, and Luis A Aguirre, “Failure in distinguishing colored noise from chaos using the ânoise titrationâ technique,” Physical Review E 79, 035201 (2009)
 (38) C. W. J. Granger, ‘‘Investigating causal relations by econometric models and crossspectral methods,” Econometrica 37, pp. 424–438 (1969), ISSN 00129682, http://www.jstor.org/stable/1912791