Compressing phase space detects state changes in nonlinear dynamical systems

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.


I. 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][2][3]. Phase space reconstruction has been widely applied in modeling and predictions of several nonlinear systems, such as ecological, climate and neural ones [4][5][6][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 X(t) by vectors in R m of the form X n = [X(n), X(n − τ ), x(n − 2τ ), ..., X(n − (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 on M is ensured by choosing an embedding dimension m > 2d + 1 [9] and the result may be generalized also to non-integer (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 to the fact that, when evaluating some quantities on a R m delay-coordinate map, they do not vary for m higher than the proper embedding dimension. Those diffeomorphism invariants could be, for examples, 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, that 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 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 a rich information about the underlying dynamics and reveal change in the system state that are otherwise difficult to detect [25][26][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.

II. LOW-DIMENSIONAL CHAOTIC SYSTEMS
To estimate the optimal embedding dimension m, we built M X (τ, m), an ensemble of delay-coordinate maps from X(t) as a function of time delay τ and m. Then, at fixed τ and m, we symbolized each degree of freedom so generating a new discrete variable with a scale that depends on the bin size ǫ used to make discrete the delay-coordinate map: X discrete = X n (ǫ, τ, m) = (x 1 , x 2 , ..., x n ). We computed the entropy rate of the resulting sequence of symbols through a Lempel-Ziv data compression algorithm [28]: i is the shortest sub-sequence starting at index i that does not appear in the window x i−1 i−n of length n. We evaluated entropy rate for the entire ensemble of delay-coordinate maps M X (m) and estimated as the optimal embedding dimension m the one such that S[MX (m+1)] S[MX (m)] = cost. This means that the optimum embedding dimension is the one at which entropy rate has at least a component that behaves as a non linear function of m, that is S(m) ∼ c 1 e c2m . That choice was suggested by the fact [29] that system with causal interactions among their elements have entropy that grows as a non extensive function of their size S(N ) = S 0 N + S 1 (N ), where the non extensive component is described by a power law function S 1 ∼ N m .
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. Fig. 1a shows an example for a single realization of a Lorenz signal, where estimation of optimal m does not change for different τ values. Fig 1b 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 compare our estimates with correlation dimension (d 2 ) measures [30][31][32][33]. We found that for most of the tested systems, our dimension estimate is close to the Whitney's upper bound 2d 2 + 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 d 2 .

III. UNDIRECTIONALLY COUPLED SYSTEMS
In coupled chaotic systems with a drive-response configuration, generalized synchronization (GS) may occur if the state of response system X does not depends on its initial condition but depends only on the state of the driver Y , that is, if there is a functional relation between trajectories in the phase-space: X(t) = Φ(Y (t)). When Φ is the identity, there is identical synchronization, that is easy to detect because the synchronized motion becomes simply a sharp line in X(t) vs Y (t) 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 the kind of synchronization that is acting by considering 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 an effect on the global dimension and there is identical synchronization. Otherwise, if d G > d D , 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 factor C. We studied the optimal dimension m of the systems assuming, as we did for non coupled systems, that entropy is well described by a non extensive function of number of elements S ∼ N m . That assumption is especially well posed when the system is weakly sensitive to initial conditions, where it was proven [38,39] that he 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 from the multifractal spectrum.

A. Heterogeneous systems
As a first example we considered an unidirectionally coupled system in which the autonomous driver X is a Rossler oscillator: and the driven one, Y, is a Lorenz oscillator:     ẏ 1 = 10(−y 1 + y 2 ) y 2 = 28y 1 − y 2 − y1y 3 + Cx 2 2 y 3 = y 1 y 2 − 2.66y 3 (2) This type of system was investigated in previous works [40][41][42]. In Fig. 2a we show that, similarly to previous studies, GS arises for a threshold coupling strength C = C w > 2.1, where the conditional Lyapunov exponent becomes negative. We computed Lyapunov exponents using the pull-back method [43,44], that 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 d 2 is estimated by using 25000 time points and looking for the plateau in the function d 2 (m, ǫ) [45], indicating a suitable scaling relationship. As show in Fig. 2b, d 2 of the global manifold is higher than d 2 of the driver Rossler system, indicating that at the threshold C w the whole system undergoes a regime of weak synchronization.
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. As approaching the synchronization threshold C w , m increases abruptly and assumes values between the two extremes of two independent Lorenz and Rossler systems. Furthermore, is worth noting that the trend of m estimates is opposite to the trends of both d 2 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, that is the dimension at which the entropy rate is described by a non linear function, is related to the length of the patterns, i.e., to regularities in the dynamics that allow for its compression.

B. Identical systems
A second example we considered is the undirectionally coupled system formed by two identical Hénon maps [41], where the driver is described by the system: and the driven one by: We computed Lyapunov exponents using the pullback method with 5000 time points and we found that the conditional exponent takes negative values in two different intervals of couplings: in a window 0.44 < C w < 0.54 and then for C i > 0.68 (see Fig. 3a). As shown in Fig. 3b, in the first window C w , the correlation dimension of the global manifold is higher than the correlation dimension of an independent Hénon system: d G ≈ 2.2 > d Henon = 1.2, indicating that the synchronization is weak in this interval. Furthermore, for coupling values higher than C i we have that d G = d Henon = 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 C w interval the complexity of the coupled system increases, giving rise to an increase estimated m of the global manifold. For C > C i 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.

IV. CONCLUSION
In conclusions, 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 departs from extensive functions and become non linear functions of the number of system dimensions. How to proper 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 don't 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 are investigated quantities related to the randomness of the dynamics [49,50], such as Lyapunov exponents and fractal dimension. However, to be estimated in a reliable way, those quantities require long time series, in particular to compute correlation dimension, which is also poten-tially biased by user's choices about proper scale and dimensions. Our method, on the contrary, give 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 phenomena that appear in many biological and physiological processes [51][52][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 multi-scale 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.