Compressing phase space detects state changes in nonlinear dynamical systems
aa r X i v : . [ n li n . C D ] J un Compressing phase space detects state changes in nonlinear dynamical systems
Valeria d’Andrea ∗ and Manlio De Domenico † CoMuNe Lab, Fondazione Bruno Kessler, Via Sommarive 18, 38123 Povo (TN), Italy (Dated: June 24, 2020)Equations governing the nonlinear dynamics of complex systems are usually unknown and indirectmethods are used to reconstruct their manifolds. In turn, they depend on embedding parametersrequiring other methods and long temporal sequences to be accurate. In this paper, we show that anoptimal reconstruction can be achieved by lossless compression of system’s time course, providinga self-consistent analysis of its dynamics and a measure of its complexity, even for short sequences.Our measure of complexity detects system’s state changes such as weak synchronization phenomena,characterizing many systems, in one step, integrating results from Lyapunov and fractal analysis.
I. INTRODUCTION
Dynamics of natural systems is often describedby nonlinear equations. When those equations areunknown, we can reproduce the system dynamicsthrough the reconstruction of the manifold from thetime course of one of its variables [1–3]. Phase spacereconstruction has been widely applied in modelingand predictions of several nonlinear systems, suchas ecological, climate and neural ones [4–7]. Theembedding theory proposed by Takens [8] allows oneto reconstruct a one-to-one map of the attractor of adynamical process using time-lagged values of a singlesystem variable. The delay-coordinate map is builtfrom the time series X ( t ) by vectors in R m of the form X n = [ X ( n ) , X ( n − τ ) , x ( n − τ ) , ..., X ( n − ( m − τ )],where τ is the time delay. To correctly build theembedding of d -dimensional manifold M it is crucial tochoose adequate values for m and τ , i.e. the embeddingparameters.According to the Whitney theorem, the diffeomorphismon M is ensured by choosing an embedding dimension m > d + 1 [9] and the result may be generalized also tonon-integer (fractal) dimension [10]. Whitney theoremhas been relaxed, for example in [11, 12], but still thosestudies provide an upper bound for the estimation of m . Several methods were developed to estimate theminimum possible embedding dimension [13] and usuallythose methods are based to the fact that, when evaluat-ing some quantities on a R m delay-coordinate map, theydo not vary for m higher than the proper embeddingdimension. Those diffeomorphism invariants could be,for examples, the largest Lyapunov exponent or thepercentage of false nearest neighbours [14, 15], wherethe latter option, in its implementation introduced byCao [16], is currently the most used method to estimatethe minimum m .To estimate the embedding dimension, methods thatinvolve the fact that entropies are diffeomorphisminvariants have been proposed and include, for example, ∗ Corresponding author: [email protected] † Corresponding author: [email protected] differential entropy [17] and permutation entropy [18],where the latter has the advantage to take into accountthe temporal information contained in the time series[19]. Kolmogorov complexity, also known as algorithmicentropy, was proposed in 1968 as a measure of theamount of information of the trajectory of a dynamicalprocess [20] and is defined as the length of the shortestdescription that produces the trajectory as output.Even if Kolmogorov complexity cannot be computed,for the trajectories of a dynamical system it is usuallyapproximated using lossless compression algorithms,following the theorem of Brudno who, in 1978, wrote theequality between Kolmogorov complexity and entropyrate [21]. Nevertheless, to date, estimating embeddingdimension is still far from being an easy task, althoughthis parameter is critical to gain insights about thephysics of the underlying dynamical system.In this paper we show that optimal embeddingdimension can be estimated through a measure of theKolmogorov complexity, that is here evaluated usingthe compression algorithm introduced by Lempel andZiv [22]. Our dimension estimate could represent amore robust measure than other information estimatorsbecause is independent on the system representation [23]so it may be estimated without prior knowledge of thevalue of optimal time-delay τ [24]. The main advantageof our approach is that we explore the geometry ofthe manifold of the dynamical system with complexitymeasures that capture a rich information about theunderlying dynamics and reveal change in the systemstate that are otherwise difficult to detect [25–27]. Inparticular, here we show that exploring how the systemapproaches its proper embedding dimension can revealthe emergence of chaotic synchronization phenomena ina coupled drive-response system. II. LOW-DIMENSIONAL CHAOTIC SYSTEMS
To estimate the optimal embedding dimension m ,we built M X ( τ, m ), an ensemble of delay-coordinatemaps from X ( t ) as a function of time delay τ and m . Then, at fixed τ and m , we symbolized eachdegree of freedom so generating a new discrete vari-able with a scale that depends on the bin size ǫ used to make discrete the delay-coordinate map: X discrete = X n ( ǫ, τ, m ) = ( x , x , ..., x n ). We computedthe entropy rate of the resulting sequence of symbolsthrough a Lempel-Ziv data compression algorithm[28]: S = ( n n P i =2 L ii log i ) − , where L ni is the shortestsub-sequence starting at index i that does not appearin the window x i − i − n of length n . We evaluated entropyrate for the entire ensemble of delay-coordinate maps M X ( m ) and estimated as the optimal embeddingdimension m the one such that S [ M X ( m +1)] S [ M X ( m )] = cost . Thismeans that the optimum embedding dimension is theone at which entropy rate has at least a componentthat behaves as a non linear function of m , that is S ( m ) ∼ c e c m . That choice was suggested by the fact[29] that system with causal interactions among theirelements have entropy that grows as a non extensivefunction of their size S ( N ) = S N + S ( N ), where thenon extensive component is described by a power lawfunction S ∼ N m .To estimate the optimal dimension for the embedding,avoiding the evaluation of the optimal time delay τ , wetested our algorithm with a specific set of τ values andfound robust results with respect to the choice of thisparameter. Fig. 1a shows an example for a single real-ization of a Lorenz signal, where estimation of optimal m does not change for different τ values. Fig 1b showsthe results of our algorithm for a set of chaotic systems.Specifically, we consider Logistic, H´enon and Ikedamaps, Rossler, Lorenz and Mackey-Glass systems withthree different time delays, widely used to model thedynamics of several natural phenomena, from chemicalreactions to climate. For each system we computed ourmeasures across 50 different realizations and compareour estimates with correlation dimension ( d ) measures[30–33]. We found that for most of the tested systems,our dimension estimate is close to the Whitney’s upperbound 2 d + 1, while, for Mackey- Glass systems, that wetested at three different time delays, we found that our m measures are close to the lower bound delimited by d . III. UNDIRECTIONALLY COUPLED SYSTEMS
In coupled chaotic systems with a drive-response con-figuration, generalized synchronization (GS) may occurif the state of response system X does not depends onits initial condition but depends only on the state of thedriver Y , that is, if there is a functional relation betweentrajectories in the phase-space: X ( t ) = Φ( Y ( t )). WhenΦ is the identity, there is identical synchronization, thatis easy to detect because the synchronized motion be-comes simply a sharp line in X ( t ) vs Y ( t ) plane [34].Otherwise, when Φ differs from the identity, weak GS a)b) FIG. 1. Estimated embedding dimensions for low dimensionalchaotic system. (a) First derivative of S [ M X ( m +1)] S [ M X ( m )] as a func-tion of m at different time delays τ for a Lorenz dynamicalsystem. Across all τ values the derivative reaches zero for m = 4. (b) estimated m for a set of chaotic systems. Forcomparison purposes values are plotted as a function of cor-relation dimension d . For each dynamical system we showmean and standard error of the mean evaluated across N = 50realizations. Dashed line corresponds to m = d , dotted lineis m = 2 d + 1 may emerge and this phenomenon is difficult to detect.Different methods to detect GS have been proposed [35].For instance, it has been proven that synchronizationoccurs when all of the conditional Lyapunov exponentsare negative [36], while it is possible to gain insight intothe the kind of synchronization that is acting by consider-ing the dimension of the global synchronization manifold d G with respect to the dimension of the driver system d D :if d G = d D then the response system does not have aneffect on the global dimension and there is identical syn-chronization. Otherwise, if d G > d D , the global manifoldhas a fractal structure and the synchronization is weak[37]. To reveal weak GS in a coupled system, two differ-ent classes of measures are needed, namely conditionalLyapunov exponents and dimension(s) of the global man-ifold.Here we show that the analysis of the dimension ofthe response system through lossless complexity mea-sures can easily detect the emergence of GS. To this aim,we studied synchronization phenomena between two uni-directional chaotic systems, where GS takes place as afunction of coupling factor C . We studied the optimaldimension m of the systems assuming, as we did for noncoupled systems, that entropy is well described by a nonextensive function of number of elements S ∼ N m . Thatassumption is especially well posed when the system isweakly sensitive to initial conditions, where it was proven[38, 39] that he usual Shannon entropy measures are notappropriate, and a new measure of entropy has to be in-troduced, that depends on sensitivity to initial conditionsand from the multifractal spectrum. A. Heterogeneous systems
As a first example we considered an unidirectionallycoupled system in which the autonomous driver X is aRossler oscillator: ˙ x = − { x + x } ˙ x = 6 { x + 0 . x } ˙ x = 6 { . x ( x − . } (1)and the driven one, Y , is a Lorenz oscillator: ˙ y = 10( − y + y )˙ y = 28 y − y − y y + Cx ˙ y = y y − . y (2)This type of system was investigated in previous works[40–42]. In Fig. 2a we show that, similarly to previ-ous studies, GS arises for a threshold coupling strength C = C w > .
1, where the conditional Lyapunov expo-nent becomes negative. We computed Lyapunov expo-nents using the pull-back method [43, 44], that relieson the Gram-Schmidt orthonormalization of Lyapunovvectors while integrating the dynamical system with afourth order Runge-Kutta algorithm (integration timestep dt = 0 . d is estimated by using 25000time points and looking for the plateau in the function d ( m, ǫ ) [45], indicating a suitable scaling relationship.As show in Fig. 2b, d of the global manifold is higherthan d of the driver Rossler system, indicating that atthe threshold C w the whole system undergoes a regimeof weak synchronization.For each coupling value, we estimated the optimal em-bedding dimension as the average across 50 realizationsof the system dynamics. Time series with 1000 timepoints were used for the estimation. As approachingthe synchronization threshold C w , m increases abruptlyand assumes values between the two extremes of two a)b) FIG. 2. Lorenz driven by Rossler system. ( a ) Estimated m and correlation dimension d as a function of couplingstrength. ( b ) Estimated m and conditional Lyapunov expo-nents . independent Lorenz and Rossler systems. Furthermore,is worth noting that the trend of m estimates is oppositeto the trends of both d and conditional Lyapunov ex-ponents, suggesting that those measures are referring todifferent but complementary properties of the dynamicalsystem. Previous studies investigated how measures ofentropy and complexity are both needed to describenatural systems, since they capture different propertiesof the dynamics [46, 47]. In particular, Lyapunovexponents and fractal dimension measures were usuallyrelated to the degree of randomness and disorder of thedynamics, while our hypothesis is that m , that is thedimension at which the entropy rate is described bya non linear function, is related to the length of thepatterns, i.e., to regularities in the dynamics that allowfor its compression. B. Identical systems
A second example we considered is the undirectionallycoupled system formed by two identical H´enon maps [41],where the driver is described by the system: (cid:26) ˙ x = 1 . − x + 0 . x ˙ x = x (3)and the driven one by: ( ˙ y = 1 . − ( Cx y + (1 − C ) y ) + 0 . y ˙ y = y (4)We computed Lyapunov exponents using the pull-back method with 5000 time points and we foundthat the conditional exponent takes negative valuesin two different intervals of couplings: in a window0 . < C w < .
54 and then for C i > .
68 (see Fig. 3a).As shown in Fig. 3b, in the first window C w , thecorrelation dimension of the global manifold is higherthan the correlation dimension of an independent H´enonsystem: d G ≈ . > d Henon = 1 .
2, indicating thatthe synchronization is weak in this interval. Further-more, for coupling values higher than C i we have that d G = d Henon = 1 .
2, showing that for high couplingsidentical synchronization takes place. Both the couplingstrength intervals and the two different regimes for GSare revealed with a single embedding measure. Here wecomputed for each coupling value the optimal embeddingdimension m as the average across 50 realizations, 1000time points each, using lossless compression of thedynamics. We found that in C w interval the complexityof the coupled system increases, giving rise to an increaseestimated m of the global manifold. For C > C i theoptimal m has a drop, showing that there is a change inthe system state, in particular m estimates take valuestypical of an independent H´enon map, revealing anidentical synchronization. IV. CONCLUSION
In conclusions, we have shown that complexity mea-sures used to reconstruct the geometry of the manifoldof a dynamical system can be used to gain many insightsabout the system itself, even when the underlyinggoverning equations are not known. We observed howthe irregularity of the dynamics, expressed by entropyrate estimates, reaches a plateau and remains constantby increasing the dimension of the manifold, providinga robust and parameter-free estimate of the intrinsicoptimal dimension. Our measure is quite stable fordifferent values of time-delay τ , providing a desirablemethod for the reconstruction of the manifold that relies a)b) FIG. 3. Coupled H´enon maps. ( a ), Estimated embedding m and correlation dimension d as a function of couplingstrength. ( b ), Estimated m and conditional Lyapunov ex-ponents. only on a single estimate.We choose to relate complexity of the system to theway at which entropy rate measures departs from exten-sive functions and become non linear functions of thenumber of system dimensions. How to proper evaluatecomplexity has been a debated topic in last years. Oneof the most debated issue is the fact that informationtheoretic estimates like Shannon entropy measure thedegree of randomness of the system and don’t take intoaccount system’s dynamical organization, whereas idealcomplexity measures should treat both random andlower distributions as minimally complex [48]. In ourapproach we focused on the entropy component thatdeviates from extensivity, arguing that it contains theinformation that has to be related to effective system’scomplexityTo detect synchronization usually are investigatedquantities related to the randomness of the dynamics[49, 50], such as Lyapunov exponents and fractal di-mension. However, to be estimated in a reliable way,those quantities require long time series, in particularto compute correlation dimension, which is also poten-tially biased by user’s choices about proper scale anddimensions. Our method, on the contrary, give robustresults for shorter time series and has the advantageto capture and distinguish, with a single measure,different synchronization regimes. Furthermore, thedimension at which the time series reaches its maximumdisorder is informative and gives us insights about theintrinsic structure of the system. The way in which theoptimal embedding dimension varies as a function of theparameters ruling the system dynamics highlights statechanges, as long as they affect regularities in dynamicalpatterns. In this paper, we focused more specifically onthe detection of generalized synchronization in coupledchaotic systems, a phenomena that appear in manybiological and physiological processes [51–53], as well as in geophysical fluid dynamics [54], but it is notoriouslydifficult to unravel. Additionally, the detection ofsynchronization phenomena permits the identificationof causal drivers and leads to a better description andprediction of system dynamics. The key role that causalinfluence among observables has for the forecasting oftheir time course has been addressed in many studiesrelated, for example, to ecological [55], financial [56]and multi-scale human mobility systems [57, 58]. Ourmethod paves the way for applications to more complexdynamics exhibiting phenomena that usually requiremultiple complexity measures to be detected, showingthat lossless compression of system’s dynamics in thephase space can be suitably used for this purpose. [1] E. Ott, T. D. Sauer, and J. A. Yorke, Coping with chaos:analysis of chaotic data and the exploitation of chaoticsystems (John Wiley and Sons Ltd, New York, NY, USA,1994).[2] N. Packard, J. Crutchfield, J. Farmer, and R. Shaw,Physical Review Letters , 712 (1980).[3] J. D. Farmer and J. J. Sidorowich, Physical Review Let-ters , 845 (1987).[4] M. Casdagli, Physica D: Nonlinear Phenomena , 335(1989).[5] M. Chavez, M. Valencia, V. Navarro, V. Latora, andJ. Martinerie, Physical Review Letters , 1 (2010).[6] G. Sugihara and H. Ye, Science , 922 (2016).[7] S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser,and J. Nathan Kutz, Nature Communications , 1(2017).[8] F. Takens, Springer (1981).[9] H. Whitney, Annals of Mathematics , 645 (1936).[10] T. Sauer, J. A. Yorke, and M. Casdagli, Journal of sta-tistical Physics , 579 (1991).[11] A. Haefliger and M. W. Hirsch, Topology , 129 (1963).[12] F. Fuquan, Topology , 447 (1994).[13] C. J. Cellucci, A. M. Albano, and P. E. Rapp, Physi-cal Review E - Statistical Physics, Plasmas, Fluids, andRelated Interdisciplinary Topics , 13 (2003).[14] M. B. Kennel, R. Brown, and D. I. Abarbanel, PhysicalReview A (1992).[15] H. D. Abarbanel and M. B. Kennel, Physical Review E , 3057 (1993).[16] L. Cao, Physica D: Nonlinear Phenomena , 43 (1997).[17] T. Gautama, D. P. Mandic, and M. M. Van Hulle,ICASSP, IEEE International Conference on Acoustics,Speech and Signal Processing - Proceedings , 29 (2003).[18] M. Riedl, A. M¨uller, and N. Wessel, European PhysicalJournal: Special Topics , 249 (2013).[19] C. Bandt and B. Pompe, Physical Review Letters , 4(2002).[20] A. N. Kolmogorov, International Journal of ComputerMathematics , 157 (1968).[21] Brudno, Russian Math. Surveys , 197 (1978).[22] J. Ziv and A. Lempel, IEEE Transactions on InformationTheory , 530 (1978). [23] F. Argenti, V. Benci, P. Cerrai, A. Cordelli, S. Galatolo,and G. Menconi, Chaos, Solitons and Fractals , 461(2002).[24] A. M. Fraser and H. L. Swinney, Phys. Rev. A , 1134(1986).[25] O. Melchert and A. K. Hartmann, Physical Review E- Statistical, Nonlinear, and Soft Matter Physics , 1(2015).[26] R. Avinery, M. Kornreich, and R. Beck, Physical reviewletters , 178102 (2017).[27] S. Martiniani, P. M. Chaikin, and D. Levine, PhysicalReview X , 11031 (2019).[28] M. Lin, W.-J. Hsu, and Z. Q. Lee, UbiComp , 381 (2012).[29] W. Bialek, I. Nemenman, and N. Tishby, Neural Comput. , 2409 (2001).[30] P. Grassberger and I. Procaccia, Physical Review Letters (1983).[31] P. Grassberger and I. Procaccia, Physica D , 170 (1983).[32] P. Grassberger, Physics Letters A , 224 (1983).[33] P. Grassberger and I. Procaccia, Physical Review A ,2591 (1983).[34] N. Rulkov and D. I. Abarbanel, Physical Review (1995).[35] S. Boccaletti, A. N. Pisarchik, C. I. del Genio, andA. Amann, Synchronization: From Coupled Systems toComplex Networks (Cambridge University Press, 2018).[36] L. M. Pecora and T. L. Carroll, Physical Review Letters (1990).[37] K. Pyragas, Physical Review E - Statistical Physics, Plas-mas, Fluids, and Related Interdisciplinary Topics ,5183 (1997).[38] V. Latora, M. Baranger, A. Rapisarda, and C. Tsallis,Physical Letters A , 97 (2000).[39] M. Ponce-Flores, J. Frausto-Sols, G. Santamara-Bonfil,J. Prez Ortega, and J. J. Gonzlez-Barbosa, Entropy ,1 (2020).[40] K. Pyragas, Physical Review E , 63 (1996).[41] R. Q. Quiroga, J. Arnhold, and P. Grassberger, Physi-cal Review E - Statistical Physics, Plasmas, Fluids, andRelated Interdisciplinary Topics , 5142 (2000).[42] R. Quian Quiroga, A. Kraskov, T. Kreuz, and P. Grass-berger, Physical Review E - Statistical Physics, Plas-mas, Fluids, and Related Interdisciplinary Topics , 14 (2002).[43] G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn,Meccanica , 9 (1980).[44] A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano,Physica D: Nonlinear Phenomena , 285 (1985).[45] R. Hegger, H. Kantz, and T. Schreiber, Chaos , 413(1999).[46] J. P. Crutchfield and D. P. Feldman, Chaos , 25 (2003).[47] D. P. Feldman, C. S. Mctague, J. P. Crutchfield, D. P.Feldman, C. S. Mctague, and J. P. Crutchfield, Chaos (2008).[48] T. Deacon and S. Koutroufinis, Information , 404(2014).[49] H. D. I. Abarbanel, N. F. Rulkov, and M. M. Sushchik,Physical Review E , 4528 (1996).[50] S. Boccaletti, J. Kurths, G. Osipov, D. L. Valladares, andC. S. Zhou, Physics Reports , 1 (2002).[51] L. Glass, Nature , 277 (2001).[52] C. Chen, S. Liu, X.-q. Shi, H. Chat´e, and Y. Wu, NaturePublishing Group (2017). [53] M. Frasca, A. Buscarino, A. Rizzo, L. Fortuna, andS. Boccaletti, Physical Review Letters , 1 (2008).[54] G. S. Duane and J. J. Tribbia, Physical Reviews Letters , 4298 (2001).[55] H. Ye, R. J. Beamish, S. M. Glaser, S. C. Grant, C. H.Hsieh, L. J. Richards, J. T. Schnute, and G. Sugihara,Proceedings of the National Academy of Sciences of theUnited States of America , E1569 (2015).[56] R. M. May, Levin Simon, and G. Sugihara, Nature ,893 (2008).[57] H. Barbosa, M. Barth´elemy, G. Ghoshal, C. R. James,M. Lenormand, T. Louail, R. Menezes, J. J. Ramasco,F. Simini, and M. Tomasini, Physics Reports , 1(2018).[58] M. De Domenico, A. Lima, and M. Musolesi, Pervasiveand Mobile Computing9