# Compressing Phase Space Detects State Changes in Nonlinear Dynamical Systems

Authors: Valeria d’Andrea; Manlio De Domenico
Journal: Complexity
 (2020)

---

## Abstract

Equations governing the nonlinear dynamics of complex systems are usually unknown, and indirect methods are used to reconstruct their manifolds. In turn, they depend on embedding parameters requiring other methods and long temporal sequences to be accurate. In this paper, we show that an optimal reconstruction can be achieved by lossless compression of system’s time course, providing a 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.


---

## Body

## 1. Introduction

Dynamics of natural systems is often described by nonlinear equations. When those equations are unknown, we can reproduce the system dynamics through the reconstruction of the manifold from the time course of one of its variables [1–3]. Phase space reconstruction has been widely applied in modeling and predictions of several nonlinear systems, such as ecological, climate, and neural ones [4–7]. The embedding theory proposed by Takens [8] allows one to reconstruct a one-to-one map of the attractor of a dynamical process using time-lagged values of a single system variable. The delay-coordinate map is built from the time series Xt by vectors in Rm of the form Xn=Xn,Xn−τ,xn−2τ,…,Xn−m−1τ, where τ is the time delay. To correctly build the embedding of d-dimensional manifold M, it is crucial to choose adequate values for m and τ, i.e., the embedding parameters.According to the Whitney theorem, the diffeomorphism onM is ensured by choosing an embedding dimension m>2d+1 [9] and the result may be generalized also to noninteger (fractal) dimension [10]. Whitney theorem has been relaxed, for example, in [11, 12], but still those studies provide an upper bound for the estimation of m. Several methods were developed to estimate the minimum possible embedding dimension [13], and usually those methods are based on the fact that when evaluating some quantities on a Rm delay-coordinate map, they do not vary for m higher than the proper embedding dimension. Those diffeomorphism invariants could be, for example, the largest Lyapunov exponent or the percentage of false nearest neighbours [14, 15], where the latter option, in its implementation introduced by Cao [16], is currently the most used method to estimate the minimum m.To estimate the embedding dimension, methods that involve the fact that entropies are diffeomorphism invariants have been proposed and include, for example, differential entropy [17] and permutation entropy [18], where the latter has the advantage to take into account the temporal information contained in the time series [19]. Kolmogorov complexity, also known as algorithmic entropy, was proposed in 1968 as a measure of the amount of information of the trajectory of a dynamical process [20] and is defined as the length of the shortest description that produces the trajectory as output. Even if Kolmogorov complexity cannot be computed, for the trajectories of a dynamical system, it is usually approximated using lossless compression algorithms, following the theorem of Brudno who, in 1978, wrote the equality between Kolmogorov complexity and entropy rate [21]. Nevertheless, to date, estimating embedding dimension is still far from being an easy task, although this parameter is critical to gain insights about the physics of the underlying dynamical system.In this paper, we show that optimal embedding dimension can be estimated through a measure of the Kolmogorov complexity, which is here evaluated using the compression algorithm introduced by Lempel and Ziv [22]. Our dimension estimate could represent a more robust measure than other information estimators because it is independent on the system representation [23], so it may be estimated without prior knowledge of the value of optimal time delay τ [24]. The main advantage of our approach is that we explore the geometry of the manifold of the dynamical system with complexity measures that capture the rich information about the underlying dynamics and reveal change in the system state that is otherwise difficult to detect [25–27]. In particular, here we show that exploring how the system approaches its proper embedding dimension can reveal the emergence of chaotic synchronization phenomena in a coupled drive-response system.


## 2. Low-Dimensional Chaotic Systems

To estimate the optimal embedding dimensionm, we built MXτ,m, an ensemble of delay-coordinate maps from Xt as a function of time delay τ and m. Then, at fixed τ and m, we discretize values of delay coordinate map by using a grid with bin size ε: Xdiscrete=Xnε,τ,m=x1,x2,…,xn. We computed the entropy rate of the resulting sequence of symbols through a Lempel–Ziv data compression algorithm [28]:(1)S=1n∑i=2nLiilogi−1,where Lin is the shortest subsequence starting at index i that does not appear in the window xi−ni−1 of length n. We evaluated entropy rate for the entire ensemble of delay-coordinate maps MXm and estimated as the optimal embedding dimension m the one such that(2)SMXm+1SMXm=cost.This means that the optimum embedding dimension is the one at which entropy rate has at least a component that behaves as a nonlinear function ofm, that is, Sm∼c1ec2m. That choice was suggested by the fact [29] that the system with causal interactions among their elements has entropy that grows as a nonextensive function of their size SN=S0N+S1N, where the nonextensive component is described by a power law function S1∼Nm.To estimate the optimal dimension for the embedding, avoiding the evaluation of the optimal time delayτ, we tested our algorithm with a specific set of τ values and found robust results with respect to the choice of this parameter. Figure 1(a) shows an example for a single realization of a Lorenz signal, where estimation of optimal m does not change for different τ values. Figure 1(b) shows the results of our algorithm for a set of chaotic systems. Specifically, we consider Logistic, Hénon and Ikeda maps, Rossler, Lorenz, and Mackey–Glass systems with three different time delays, widely used to model the dynamics of several natural phenomena, from chemical reactions to climate. For each system, we computed our measures across 50 different realizations and compared our estimates with correlation dimension (d2) measures [30–33]. We found that for most of the tested systems, our dimension estimate is close to Whitney’s upper bound 2d2+1, while for Mackey–Glass systems that we tested at three different time delays, we found that our m measures are close to the lower bound delimited by d2.Figure 1
Estimated embedding dimensions for low-dimensional chaotic system. (a) First derivative ofSMXm+1/SMXm as a function of m at different time delays τ for a Lorenz dynamical system. Across all τ values, the derivative reaches zero for m=4. (b) Estimated m for a set of chaotic systems. For comparison purposes, values are plotted as a function of correlation dimension d2. For each dynamical system, we show mean and standard error of the mean evaluated across N=50 realizations. Dashed line corresponds to m=d2; dotted line is m=2d2+1.
(a)(b)


## 3. Unidirectionally Coupled Systems

In coupled chaotic systems with a drive-response configuration, generalized synchronization (GS) may occur if the state of response systemX does not depend on its initial condition but depend only on the state of the driver Y, that is, if there is a functional relation between trajectories in the phase space, Xt=ΦYt. When Φ is the identity, there is identical synchronization, which is easy to detect because the synchronized motion becomes simply a sharp line in Xt vs Yt plane [34]. Otherwise, when Φ differs from the identity, weak GS 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 synchronization occurs when all of the conditional Lyapunov exponents are negative [36], while it is possible to gain insight into the kind of synchronization that is acting by considering the dimension of the global synchronization manifold dG with respect to the dimension of the driver system dD: if dG=dD, then the response system does not have an effect on the global dimension and there is identical synchronization. Otherwise, if dG>dD, the global manifold has a fractal structure and the synchronization is weak [37]. To reveal weak GS in a coupled system, two different classes of measures are needed, namely, conditional Lyapunov exponents and dimension(s) of the global manifold.Here, we show that the analysis of the dimension of the response system through lossless complexity measures can easily detect the emergence of GS. To this aim, we studied synchronization phenomena between two unidirectional chaotic systems, where GS takes place as a function of coupling factorC. We studied the optimal dimension 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∼Nm. That assumption is especially well posed when the system is weakly sensitive to initial conditions, where it was proven [38, 39] that the usual Shannon entropy measures are not appropriate and a new measure of entropy has to be introduced that depends on sensitivity to initial conditions and the multifractal spectrum.


### 3.1. Heterogeneous Systems

As a first example, we considered a unidirectionally coupled system in which the autonomous driverX is a Rossler oscillator:(3)x˙1=−6x2+x3,x˙2=6x1+0.2x2,x˙3=60.2+x3x1−5.7,and the driven one, Y, is a Lorenz oscillator:(4)y˙1=10−y1+y2,y˙2=28y1−y2−y1y3+Cx22,y˙3=y1y2−2.66y3.This type of system was investigated in previous works [40–42]. In Figure 2(a), we show that, similarly to previous studies, GS arises for a threshold coupling strength C=Cw>2.1, where the conditional Lyapunov exponent becomes negative. We computed Lyapunov exponents using the pull-back method [43, 44] which relies on the Gram–Schmidt orthonormalization of Lyapunov vectors while integrating the dynamical system with a fourth-order Runge–Kutta algorithm (integration time step dt=0.01). We computed exponents with 5000 time points, after discarding the first 10000 iterations. The correlation dimension d2 is estimated by using 25000 time points and looking for the plateau in the function d2m,ε [45], indicating a suitable scaling relationship. As shown in Figure 2(b), d2 of the global manifold is higher than d2 of the driver Rossler system, indicating that at the threshold Cw, the whole system undergoes a regime of weak synchronization.Figure 2
Lorenz driven by Rossler system. (a) Estimatedm and conditional Lyapunov exponents, as a function of coupling strength. (b) Estimated m and correlation dimension d2.
(a)(b)For each coupling value, we estimated the optimal embedding dimension as the average across 50 realizations of the system dynamics. Time series with 1000 time points were used for the estimation. When approaching the synchronization thresholdCw, m increases abruptly and assumes values between the two extremes of two independent Lorenz and Rossler systems. Furthermore, it is worth noting that the trend of m estimates is opposite to the trends of both d2 and conditional Lyapunov exponents, suggesting that those measures are referring to different but complementary properties of the dynamical system. Previous studies investigated how measures of entropy and complexity are both needed to describe natural systems, since they capture different properties of the dynamics [46, 47]. In particular, Lyapunov exponents and fractal dimension measures were usually related to the degree of randomness and disorder of the dynamics, while our hypothesis is that m, which is the dimension at which the entropy rate is described by a nonlinear function, is related to the length of the patterns, i.e., to regularities in the dynamics that allow for its compression.


### 3.2. Identical Systems

A second example we considered is the unidirectionally coupled system formed by two identical Hénon maps [41], where the driver is described by the system:(5)x˙1=1.4−x12+0.3x2,x˙2=x1,and the driven one is described by(6)y˙1=1.4−Cx1y1+1−Cy12+0.3y2,y˙2=y1.We computed Lyapunov exponents using the pull-back method with 5000 time points and we found that the conditional exponent takes negative values in two different intervals of couplings: in a window0.44<Cw<0.54 and then for Ci>0.68 (see Figure 3(a)). As shown in Figure 3(b), in the first window Cw, the correlation dimension of the global manifold is higher than the correlation dimension of an independent Hénon system: dG≈2.2>dHenon=1.2, indicating that the synchronization is weak in this interval. Furthermore, for coupling values higher than Ci, we have that dG=dHenon=1.2, showing that for high couplings, identical synchronization takes place. Both the coupling strength intervals and the two different regimes for GS are revealed with a single embedding measure. Here, we computed for each coupling value the optimal embedding dimension m as the average across 50 realizations, 1000 time points each, using lossless compression of the dynamics. We found that in Cw interval, the complexity of the coupled system increases, giving rise to an increased estimated m of the global manifold. For C>Ci, the optimal m has a drop, showing that there is a change in the system state; in particular, m estimates take values typical of an independent Hénon map, revealing an identical synchronization.Figure 3
Coupled Hénon maps. (a) Estimated embeddingm and conditional Lyapunov exponents. (b) Estimated m and correlation dimension d2 as a function of coupling strength.
(a)(b)


## 3.1. Heterogeneous Systems

As a first example, we considered a unidirectionally coupled system in which the autonomous driverX is a Rossler oscillator:(3)x˙1=−6x2+x3,x˙2=6x1+0.2x2,x˙3=60.2+x3x1−5.7,and the driven one, Y, is a Lorenz oscillator:(4)y˙1=10−y1+y2,y˙2=28y1−y2−y1y3+Cx22,y˙3=y1y2−2.66y3.This type of system was investigated in previous works [40–42]. In Figure 2(a), we show that, similarly to previous studies, GS arises for a threshold coupling strength C=Cw>2.1, where the conditional Lyapunov exponent becomes negative. We computed Lyapunov exponents using the pull-back method [43, 44] which relies on the Gram–Schmidt orthonormalization of Lyapunov vectors while integrating the dynamical system with a fourth-order Runge–Kutta algorithm (integration time step dt=0.01). We computed exponents with 5000 time points, after discarding the first 10000 iterations. The correlation dimension d2 is estimated by using 25000 time points and looking for the plateau in the function d2m,ε [45], indicating a suitable scaling relationship. As shown in Figure 2(b), d2 of the global manifold is higher than d2 of the driver Rossler system, indicating that at the threshold Cw, the whole system undergoes a regime of weak synchronization.Figure 2
Lorenz driven by Rossler system. (a) Estimatedm and conditional Lyapunov exponents, as a function of coupling strength. (b) Estimated m and correlation dimension d2.
(a)(b)For each coupling value, we estimated the optimal embedding dimension as the average across 50 realizations of the system dynamics. Time series with 1000 time points were used for the estimation. When approaching the synchronization thresholdCw, m increases abruptly and assumes values between the two extremes of two independent Lorenz and Rossler systems. Furthermore, it is worth noting that the trend of m estimates is opposite to the trends of both d2 and conditional Lyapunov exponents, suggesting that those measures are referring to different but complementary properties of the dynamical system. Previous studies investigated how measures of entropy and complexity are both needed to describe natural systems, since they capture different properties of the dynamics [46, 47]. In particular, Lyapunov exponents and fractal dimension measures were usually related to the degree of randomness and disorder of the dynamics, while our hypothesis is that m, which is the dimension at which the entropy rate is described by a nonlinear function, is related to the length of the patterns, i.e., to regularities in the dynamics that allow for its compression.


## 3.2. Identical Systems

A second example we considered is the unidirectionally coupled system formed by two identical Hénon maps [41], where the driver is described by the system:(5)x˙1=1.4−x12+0.3x2,x˙2=x1,and the driven one is described by(6)y˙1=1.4−Cx1y1+1−Cy12+0.3y2,y˙2=y1.We computed Lyapunov exponents using the pull-back method with 5000 time points and we found that the conditional exponent takes negative values in two different intervals of couplings: in a window0.44<Cw<0.54 and then for Ci>0.68 (see Figure 3(a)). As shown in Figure 3(b), in the first window Cw, the correlation dimension of the global manifold is higher than the correlation dimension of an independent Hénon system: dG≈2.2>dHenon=1.2, indicating that the synchronization is weak in this interval. Furthermore, for coupling values higher than Ci, we have that dG=dHenon=1.2, showing that for high couplings, identical synchronization takes place. Both the coupling strength intervals and the two different regimes for GS are revealed with a single embedding measure. Here, we computed for each coupling value the optimal embedding dimension m as the average across 50 realizations, 1000 time points each, using lossless compression of the dynamics. We found that in Cw interval, the complexity of the coupled system increases, giving rise to an increased estimated m of the global manifold. For C>Ci, the optimal m has a drop, showing that there is a change in the system state; in particular, m estimates take values typical of an independent Hénon map, revealing an identical synchronization.Figure 3
Coupled Hénon maps. (a) Estimated embeddingm and conditional Lyapunov exponents. (b) Estimated m and correlation dimension d2 as a function of coupling strength.
(a)(b)


## 4. Conclusion

In conclusion, we have shown that complexity measures used to reconstruct the geometry of the manifold of a dynamical system can be used to gain many insights about the system itself, even when the underlying governing equations are not known. We observed how the irregularity of the dynamics, expressed by entropy rate estimates, reaches a plateau and remains constant by increasing the dimension of the manifold, providing a robust and parameter-free estimate of the intrinsic optimal dimension. Our measure is quite stable for different values of time delayτ, providing a desirable method for the reconstruction of the manifold that relies only on a single estimate.We choose to relate complexity of the system to the way at which entropy rate measures depart from extensive functions and become nonlinear functions of the number of system dimensions. How to properly evaluate complexity has been a debated topic in last years. One of the most debated issue is the fact that information theoretic estimates like Shannon entropy measure the degree of randomness of the system and do not take into account system’s dynamical organization, whereas ideal complexity measures should treat both random and lower distributions as minimally complex [48]. In our approach, we focused on the entropy component that deviates from extensivity, arguing that it contains the information that has to be related to effective system’s complexity.To detect synchronization, usually quantities related to the randomness of the dynamics [49, 50], such as Lyapunov exponents and fractal dimension, are investigated. However, to be estimated in a reliable way, those quantities require long time series, in particular to compute correlation dimension, which is also potentially biased by user’s choices about proper scale and dimensions. Our method, on the contrary, gives robust results for shorter time series and has the advantage to capture and distinguish, with a single measure, different synchronization regimes. Furthermore, the dimension at which the time series reaches its maximum disorder is informative and gives us insights about the intrinsic structure of the system. The way in which the optimal embedding dimension varies as a function of the parameters ruling the system dynamics highlights state changes, as long as they affect regularities in dynamical patterns. In this paper, we focused more specifically on the detection of generalized synchronization in coupled chaotic systems, a phenomenon that appears in many biological and physiological processes [51–53], as well as in geophysical fluid dynamics [54], but it is notoriously difficult to unravel. Additionally, the detection of synchronization phenomena permits the identification of causal drivers and leads to a better description and prediction of system dynamics. The key role that causal influence among observables has for the forecasting of their time course has been addressed in many studies related, for example, to ecological [55], financial [56], and multiscale human mobility systems [57, 58]. Our method paves the way for applications to more complex dynamics exhibiting phenomena that usually require multiple complexity measures to be detected, showing that lossless compression of system’s dynamics in the phase space can be suitably used for this purpose.



---

*Source: 8650742-2020-07-22.xml*