# Synchronization Measure Based on a Geometric Approach to Attractor Embedding Using Finite Observation Windows

Authors: Inga Timofejeva; Kristina Poskuviene; Maosen Cao; Minvydas Ragulskis
Journal: Complexity
 (2018)

---

## Abstract

A simple and effective algorithm for the identification of optimal time delays based on the geometrical properties of the embedded attractor is presented in this paper. A time series synchronization measure based on optimal time delays is derived. The approach is based on the comparison of optimal time delay sequences that are computed for segments of the considered time series. The proposed technique is validated using coupled chaotic Rössler systems.


---

## Body

## 1. Introduction

The reconstruction of geometry from a time series is one of the paradigmatic algorithms used for the computational analysis of nonlinear dynamical systems. There is a number of algorithms for the analysis of time series that consider phase spaces, that is, the set of all possible system states, of dynamical systems. For a deterministic dynamical system, the phase space is known from the mathematical model (equations of motion). Thus, with the knowledge of the system state at a particular time, it is always possible to determine the future states. However, real-world or experimental dynamical systems are usually too complex, and phase spaces of such systems are unknown due to the chaotic nature of the analyzed data. Therefore, phase space reconstruction methods had to be developed. A branch of applied dynamical system theory, analyzing such algorithms, that is, extracting the information about geometrical and topological properties of the phase space of time series obtained by performing measurements on an evolving system, is called embedology [1]. One of the first methods for reconstructing the phase space system from a time series by using a time delay variable was proposed several decades ago [2]. A complex time series has an inherent geometry, and Packard et al. [2] 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.The much celebrated Takens embedding theorem proved that it is possible to reconstruct the attractor of a dynamical system using only a single time sequence of scalar measurements [3]. The vectors accomplishing this reconstruction have the form x→t=xt,xt−τ,…,xt−m−1τ, where xt are the observations on the system, the integer m is the embedding dimension, and τ is the time difference between consecutive components (also called the time lag or the delay time) and is a multiple of the sampling time step. Takens assumed an infinite sequence of noise-free measurements and proved the existence of a diffeomorphism between the original and the reconstructed attractors for almost any choice of positive delay times τ and a sufficiently big dimension m.Although the application of the Takens time delay technique is straightforward for almost any time series, the selection of optimal reconstruction parametersτ and m is a nontrivial problem due to their dependence on the nature of the analyzed data. Thus, numerous studies have been made in this area. The selection of optimal time delay value is usually based on the optimization of a particular target function that corresponds to some measures of the quality of phase space reconstruction. Classical statistical approaches for the construction of the target function include the autocorrelation function method and the mutual information method [4–6]. Such techniques are based on maximizing the independence between the coordinates of the delay vector. This approach is heavily reliant on the geometry of the phase space reconstruction, which is appropriate for the study of parameters such as fractal dimension and Lyapunov exponents. However, the selection of such target function does not always yield best results for nonlinear time series prediction problems [7].Another class of methods for the selection of an optimalτ is based on the concept of phase space expansion. The “fill factor” as a spatial measure of the phase space was first introduced by Buzug and Pfister in [8, 9]. The fill factor quantifies an attractor’s utilization of the embedding space as a function of τ. While this method does not possess the drawbacks of statistical techniques with respect to nonlinear data, it is computationally intensive (for m>3) and is prone to overfolding of the attractor [10]. A more computationally effective algorithm based on the expansion of the area for all pairwise planar projections of the embedded attractors is presented in [11]. Similar approaches based on the spatial properties of phase spaces include the average displacement method [10], the SVF method [12], and the wavering product [13] method.A different viewpoint considers both phase space reconstruction parametersm,τ simultaneously since some experiments indicate that an irrelevant relationship between m and τ may impact the congruence between the original system and the reconstructed phase space. This class of methods includes the small-window solution [14] and the C-C method [15].Novel approaches to the selection of embedding parameters have been developed in recent years. It is shown in [16] that an improved phase space reconstruction method based on manifold embedding and Laplacian eigenmaps can reveal the structure of the chaotic attractor. A new algorithm for the estimation of the dimension of chaotic dynamical systems using neural networks and robust location estimate is proposed in [17]. Statistical analysis of the nearest neighbors is exploited for the reliable estimation of minimum embedding dimension for noisy time series in [18].It has been demonstrated in [11, 19–21] that nonuniform embedding (when all time delays are not equal) performs better than uniform embedding in a variety of applications—typical examples are causality and coupling detection and time series prediction. However, the procedure for the selection of time delays is usually related to a particular application and is implemented by introducing a specific target function which determines the utility generated by a concrete nonuniform embedding.This paper has two main objectives. The first objective is to introduce a simple and effective algorithm for the identification of optimal time delays which is based on the geometrical properties of the embedded attractor and is applicable for short time series.The second objective is to introduce a time series synchronization measure based on optimal time delays. A short review of commonly applied methods for the detection of synchronization is given below.The first linear synchronization measures based on correlation analysis, namely, cross correlation and coherence functions, are widely applied due to their computational effectiveness [22]. However, such measures can only detect the most straightforward regimes of synchronization.Nonlinear synchronization measures were introduced in order to quantify more complex synchronization effects. The mutual information measure is based on Shannon entropy and takes into account not only linear but also nonlinear dependencies [23]. Phase synchronization measure is used to quantify similarity between cyclic signals and time series. Two approaches—based on the Hilbert or wavelet transform—can be used to implement the phase synchronization detection algorithm [22].A new form of synchronization between coupled chaotic oscillators called amplitude envelope synchronization has been discovered in [24]. Generalized synchronization of dynamical systems occurs if dynamical variables from one subsystem are a function of the variable of another subsystem [25]. A nonlinear interdependence of dynamical systems based on state space reconstruction is a similar approach to generalized synchronization but does not require a functional relationship between the dynamics of the underlying systems [26]. The function projective synchronization in relay coupled systems is studied in [27]. A novel sort of synchronization called complex antilag synchronization is introduced in [28]. A new approach for the investigation of hybrid chaos synchronization in discrete-time hyperchaotic dynamical systems based on stability theory of linear discrete-time systems and Lyapunov stability theory was proposed in [29].This paper presents a novel time series synchronization measure based on the geometrical approach towards optimal time delays. This approach is based on the determination of time lags that maximize the volume of the state space occupied by the embedded segments of the time series. The sequences of obtained time lags reveal the level of synchronization between 2 time series. The presented algorithm for the identification of synchronization between two time series is validated using coupled chaotic Rössler systems.


## 2. Preliminaries

### 2.1. Two-Dimensional Delay Coordinate Space

Let us consider a harmonic functionxt=sinωt+ϕ, where ω is the angular frequency and ϕ is the phase of harmonic oscillations, and the amplitude is set to 1. A pair of function values xt and xt+τ, where τ is the time lag, is mapped into a point xt;xt+τ in the phase plane X1OX2, where X1 and X2 are the coordinates of the delay coordinate plane. Two-dimensional time delay embedding maps the harmonic function into an ellipse; the equation of that ellipse reads
(1)X2=X1cosτ+1−X12sinτ.The geometrical shape of the ellipse can be exploited for the quantification of the quality of the embedding. Such geometric approach (for a two-dimensional phase plane) was firstly proposed by Buzug and Pfister in 1992 [9]. The radiuses of the embedded ellipse r1 and r2 can be directly expressed from (1):
(2)r1=sinωτ1−cosωτ,r2=sinωτ1+cosωτ.Thus, the area of the ellipseE reads
(3)E=π⋅r1r2=πsinωτ.The maximum area of the ellipse isπ (when the ellipse becomes a circle). Now, the function representing the quality of embedding into a two-dimensional delay coordinate space can be defined as a ratio between the area of the ellipse and the area of the circle:
(4)Qτ,ω=sinωτ.


### 2.2.m-Dimensional Delay Coordinate Space

#### 2.2.1. The Embedding Quality Function

The geometrical approach for embedding a harmonic function into a two-dimensional phase plane is generalized for them-dimensional delay coordinate space in [11]. The coordinates of the reconstructed point in the m-dimensional delay coordinate space read
(5)xt,xt+τ1,xt+τ1+τ2,…,xt+τ1+⋯+τm−1,where τ1,τ2,…,τm−1 is the time delay vector; ∑k=1m−1τk is the total embedding window. It can be observed that there exist mm−1/2 planar projections of the embedded ellipse. The function representing the quality of embedding into the m-dimensional delay coordinate space is constructed as an arithmetic average of all quality functions is all possible planar projections [30]:
(6)Qτ1,τ2,…,τm−1,ω=2mm−1∑k=1m−1sinωτk+∑k=1m−1sinωτk+τk+1+⋯+sinω∑k=1m−1τk.In the case of the uniform embedding (whenτ1=τ2=⋯=τm−1=τ), (6) reads
(7)Qτ,ω=2mm−1∑k=1m−1m−ksinkωτ.The equalityQ=0 implies that the harmonic function is compressed into line segments in all possible planar projections. The maximization of Q in respect of time lags results in a maximum average area of ellipses in all possible planar projections.


#### 2.2.2. A Nonharmonic Function

Let us consider a functionxt defined in the time interval t0≤t≤t1, where t0 and t1 are the limits of the observation window. Every harmonic component of xt is affected by the quality function Qτ1,τ2,…,τm−1,ω when the appropriate harmonic signal is embedded into the m-dimensional delay coordinate space. Harmonic components with frequencies where Q is small will be suppressed (in average) in all possible planar projections more than those harmonic components where Q is large. Such motivation did help to construct the optimization problem for the identification of an optimal set of time delays which do result in the richest representation of the attractor in the delay coordinate space:
(8)τ1∗,…,τm−1∗=argmaxSτ1,…,τm−1,where
(9)Sτ1,…,τm−1=π2∫0∞AωQτ1,…,τm−1,ωdω∫0∞Aωdω,τ1∗,…,τm−1∗ is the optimal set of time delays and Aω is the Fourier amplitude spectrum of the signal in the observation window [11]. The integral in the denominator is used in order to normalize the target function in respect of the signal (this integral can be computed once at the beginning of the optimization process). The multiplier π/2 is used to normalize the target function in respect of the white noise (the value of the target function for the white noise signal now becomes equal to 1 for any time delays not equal to 0) [11].


## 2.1. Two-Dimensional Delay Coordinate Space

Let us consider a harmonic functionxt=sinωt+ϕ, where ω is the angular frequency and ϕ is the phase of harmonic oscillations, and the amplitude is set to 1. A pair of function values xt and xt+τ, where τ is the time lag, is mapped into a point xt;xt+τ in the phase plane X1OX2, where X1 and X2 are the coordinates of the delay coordinate plane. Two-dimensional time delay embedding maps the harmonic function into an ellipse; the equation of that ellipse reads
(1)X2=X1cosτ+1−X12sinτ.The geometrical shape of the ellipse can be exploited for the quantification of the quality of the embedding. Such geometric approach (for a two-dimensional phase plane) was firstly proposed by Buzug and Pfister in 1992 [9]. The radiuses of the embedded ellipse r1 and r2 can be directly expressed from (1):
(2)r1=sinωτ1−cosωτ,r2=sinωτ1+cosωτ.Thus, the area of the ellipseE reads
(3)E=π⋅r1r2=πsinωτ.The maximum area of the ellipse isπ (when the ellipse becomes a circle). Now, the function representing the quality of embedding into a two-dimensional delay coordinate space can be defined as a ratio between the area of the ellipse and the area of the circle:
(4)Qτ,ω=sinωτ.


## 2.2.m-Dimensional Delay Coordinate Space

### 2.2.1. The Embedding Quality Function

The geometrical approach for embedding a harmonic function into a two-dimensional phase plane is generalized for them-dimensional delay coordinate space in [11]. The coordinates of the reconstructed point in the m-dimensional delay coordinate space read
(5)xt,xt+τ1,xt+τ1+τ2,…,xt+τ1+⋯+τm−1,where τ1,τ2,…,τm−1 is the time delay vector; ∑k=1m−1τk is the total embedding window. It can be observed that there exist mm−1/2 planar projections of the embedded ellipse. The function representing the quality of embedding into the m-dimensional delay coordinate space is constructed as an arithmetic average of all quality functions is all possible planar projections [30]:
(6)Qτ1,τ2,…,τm−1,ω=2mm−1∑k=1m−1sinωτk+∑k=1m−1sinωτk+τk+1+⋯+sinω∑k=1m−1τk.In the case of the uniform embedding (whenτ1=τ2=⋯=τm−1=τ), (6) reads
(7)Qτ,ω=2mm−1∑k=1m−1m−ksinkωτ.The equalityQ=0 implies that the harmonic function is compressed into line segments in all possible planar projections. The maximization of Q in respect of time lags results in a maximum average area of ellipses in all possible planar projections.


### 2.2.2. A Nonharmonic Function

Let us consider a functionxt defined in the time interval t0≤t≤t1, where t0 and t1 are the limits of the observation window. Every harmonic component of xt is affected by the quality function Qτ1,τ2,…,τm−1,ω when the appropriate harmonic signal is embedded into the m-dimensional delay coordinate space. Harmonic components with frequencies where Q is small will be suppressed (in average) in all possible planar projections more than those harmonic components where Q is large. Such motivation did help to construct the optimization problem for the identification of an optimal set of time delays which do result in the richest representation of the attractor in the delay coordinate space:
(8)τ1∗,…,τm−1∗=argmaxSτ1,…,τm−1,where
(9)Sτ1,…,τm−1=π2∫0∞AωQτ1,…,τm−1,ωdω∫0∞Aωdω,τ1∗,…,τm−1∗ is the optimal set of time delays and Aω is the Fourier amplitude spectrum of the signal in the observation window [11]. The integral in the denominator is used in order to normalize the target function in respect of the signal (this integral can be computed once at the beginning of the optimization process). The multiplier π/2 is used to normalize the target function in respect of the white noise (the value of the target function for the white noise signal now becomes equal to 1 for any time delays not equal to 0) [11].


## 2.2.1. The Embedding Quality Function

The geometrical approach for embedding a harmonic function into a two-dimensional phase plane is generalized for them-dimensional delay coordinate space in [11]. The coordinates of the reconstructed point in the m-dimensional delay coordinate space read
(5)xt,xt+τ1,xt+τ1+τ2,…,xt+τ1+⋯+τm−1,where τ1,τ2,…,τm−1 is the time delay vector; ∑k=1m−1τk is the total embedding window. It can be observed that there exist mm−1/2 planar projections of the embedded ellipse. The function representing the quality of embedding into the m-dimensional delay coordinate space is constructed as an arithmetic average of all quality functions is all possible planar projections [30]:
(6)Qτ1,τ2,…,τm−1,ω=2mm−1∑k=1m−1sinωτk+∑k=1m−1sinωτk+τk+1+⋯+sinω∑k=1m−1τk.In the case of the uniform embedding (whenτ1=τ2=⋯=τm−1=τ), (6) reads
(7)Qτ,ω=2mm−1∑k=1m−1m−ksinkωτ.The equalityQ=0 implies that the harmonic function is compressed into line segments in all possible planar projections. The maximization of Q in respect of time lags results in a maximum average area of ellipses in all possible planar projections.


## 2.2.2. A Nonharmonic Function

Let us consider a functionxt defined in the time interval t0≤t≤t1, where t0 and t1 are the limits of the observation window. Every harmonic component of xt is affected by the quality function Qτ1,τ2,…,τm−1,ω when the appropriate harmonic signal is embedded into the m-dimensional delay coordinate space. Harmonic components with frequencies where Q is small will be suppressed (in average) in all possible planar projections more than those harmonic components where Q is large. Such motivation did help to construct the optimization problem for the identification of an optimal set of time delays which do result in the richest representation of the attractor in the delay coordinate space:
(8)τ1∗,…,τm−1∗=argmaxSτ1,…,τm−1,where
(9)Sτ1,…,τm−1=π2∫0∞AωQτ1,…,τm−1,ωdω∫0∞Aωdω,τ1∗,…,τm−1∗ is the optimal set of time delays and Aω is the Fourier amplitude spectrum of the signal in the observation window [11]. The integral in the denominator is used in order to normalize the target function in respect of the signal (this integral can be computed once at the beginning of the optimization process). The multiplier π/2 is used to normalize the target function in respect of the white noise (the value of the target function for the white noise signal now becomes equal to 1 for any time delays not equal to 0) [11].


## 3. The Proposed Algorithm

The technique proposed in [11] enables a fast and effective determination of the optimal set of time delays by assessing the geometrical properties of the embedded attractor. However, the optimization problem in [11] does not asses the phase spectrum of the Fourier transform. In other words, [11] is an approximate algorithm for the determination of optimal time delays when every discrete harmonic component of the Fourier amplitude spectrum is treated as a separate individual harmonic component. Nevertheless, such an approximating approach appears to be useful in practical applications—especially in time series forecasting algorithms based on fuzzy neural networks [31–33]. However, it is clear that the Fourier phase spectrum does have an impact on the geometrical representation of the reconstructed attractor in the phase space. In other words, a more accurate assessment of the geometrical properties of the reconstructed attractor is required. Let us consider a continuous function ft in the observation window 0≤t≤T and the m-dimensional delay coordinate space with the set of time lags τ=τ1,τ2,…,τm−1. Then, instead of computing any projections of the embedded attractor into a 2-dimensional phase planes, we compute the average distance of the attractor points to the origin of the embedding frame:
(10)Dτ1,τ2,…,τm−1=1Tm∫0Tf2t+f2t+τ1+f2t+τ1+τ2+⋯+f2t+∑k=1m−1τkdt.


### 3.1. Properties ofDτ1,τ2,…,τm−1: Embedding a Harmonic Function

Let us consider that the embedding is uniformτ1=τ2=⋯=τm−1=τ;T→∞ and the embedded function is a harmonic function ft=asinω0t+ϕ, where a, ω0, and ϕ are the amplitude, the circular frequency, and the phase, accordingly. Then,
(11)Dτ=amlimT→∞1T∑k=0m−1sin2ω0t+kτ+ϕdt.The change of variablesω0t+ϕ=x yields
(12)Dτ=aω0mlimT→∞1T∫ϕω0T+ϕ∑k=0m−1sin2x+ω0kτdx=aω0m12π∫02π∑k=0m−1sin2x+ω0kτdx.Now,(13)Dτ=aω0m12π∫02πsin2x∑k=0m−112+12cos2ω0kτ+12sin2x∑k=0m−1sin2ω0kτ+cos2x∑k=0m−112−12cos2ω0kτ1/2dx=aω02m12π∫02πm−cos2x∑k=0m−1cos2ω0kτ+sin2x∑k=0m−1sin2ω0kτdx.Let us denotez=ei2ω0τ. Then,
(14)Dτ=aω02m12π∫02πm−cos2x∑k=0m−1Rezk+sin2x∑k=0m−1Imzkdx=aω02m12π∫02πm−cos2xRe1−zm1−z+sin2xIm1−zm1−zdx.Note that1−zm/1−z=0 when τ=π/mω0. Thus,
(15)Dπmω0=aω02m12π∫02πmdx=a2ω0is the maximal value for the harmonic function. On the opposite, the pole z=1 is a removable singularity at τ=π/ω0 and
(16)Dπω0=aω0212π∫02π1−cos2xdx=2aω0πis the minimal value for the harmonic function.


### 3.2. The Comparison betweenDτ and Sτ

We continue with the harmonic functionft=asinω0t+ϕ and uniform embedding. Then, Aω=δω−ω0 and (9) yields Sτ=π/2Qτ,ω0, where δ is the Dirac delta function. Without the loss of generality, we assume that a=1 and ω0=1. Then, Dmin=2/π≈0.637; Dmax=1/2≈0.707. Graphs of Dτ and Sτ are illustrated in Figure 1.Figure 1
The comparison of functionsDτ (solid line) and Sτ (dashed line). The parameter m is set to 2, 3, 4, and 20 in (a), (b), (c), and (d), respectively.
(a)
(b)
(c)
(d)The first observation is that graphs ofDτ and Sτ are similar at m=2 (Figure 1(a)). The best embedding of the harmonic signal into the 2-dimensional delay coordinate space is produced at τ=π/2+πk, k∈ℤ (then the harmonic function is embedded into a circle). The harmonic function is embedded into a line interval on one of the diagonals of the frame at τ=πk, k∈ℤ. The difference is that now Sπk=0, but Dπk=2/π (note that the scale of the y-axis for Sτ is shown on the left, but the scale for Dτ is shown on the right side of Figure 1(a)).The differences betweenDτ and Sτ are more clear at m=3 (Figure 1(b)). The harmonic function cannot be mapped into circles in all 3 projections at a fixed τ; therefore, the maximum value of Sτ at m=3 is lower than that at m=2 (note that the scales of the axes are the same for different m in Figure 1). But Dmax remains the same and does not depend on the dimension of the delay coordinate space. Both Dτ and Sτ yield the best embedding at π/3l, where l∈ℤ/3k, k∈ℤ.The maximum value ofSτ at m=4 becomes smaller than the one at m=3 (Figure 1(c)). But now a major difference between Dτ and Sτ can be observed in terms of the best embedding. Sτ yields the best embedding at ±π/4+πk, but Dτ yields the best embedding at π/4l, where l∈ℤ/4k, k∈ℤ. Clearly, Sτ misses the time delay τ=π/2+πk. That can be explained by the fact that Sτ measures the average area of the mapped attractor in all possible plane projections of the m-dimensional delay coordinate space. On the contrary, Dτ directly measures the geometrical volume occupied by the attractor without considering any projections of that attractor.The maximum value ofSτ tends to 1 when m→∞ (Figure 1(d)). Note that Sτ1,…,τm−1=1 when the embedded time series is the white noise [11]. However, as mentioned previously, Dmax is an invariant value in respect of m. What is more important, Dτ yields the best embedding at π/ml, where l∈ℤ/mk, k∈ℤ, which is in stark contrast to Sτ (Figure 1(d)).


### 3.3. Embedding a Discrete Time Series

Let us consider a finite scalar discrete time seriesx1,x2,…,xn, where n is the number of observation points uniformly spread in the time axis. Let the sampling rate of the time series (the time interval between adjacent observations) be δ. Then, all time delays are multiple values of δ:τk=ikδ,ik∈ℕ,k=1,2,…,m−1. Let us denote Δ=∑k=1m−1ik. Then, (10) reads
(17)Di1δ,i2δ,…,im−1δ=1n−Δm∑k=1n−Δxk2+xk+i12+xk+i1+i22+⋯+xk+Δ2.Note thatδ does not play any role on the right-hand side of (17)—because the definite integral in (10) is replaced by the limit sum in (17). Now, the problem of the identification of the optimal set of time delays is formulated as follows:
(18)max1≤i1,…,im−1≤L1n−Δm∑k=1n−Δxk2+xk+i12+xk+i1+i22+⋯+xk+Δ2,where L is the upper limit for every individual time delay. The total embedding window (the time delay between the first and the last time delay coordinates) is Δδ, and it could be also used for the estimation of L. It is clear that n must be greater than Δ. However, n must be at least 10 times greater than m−1L in order to produce a meaningful optimization result (18).


### 3.4. Properties ofDτ1,τ2,…,τm−1: Embedding the Random Noise

It is well known that embedding of the random noise is invariant to time delays. In other words, the geometric shape of the reconstructed attractor does not depend on time delays (if only the embedded time series is an ergodic random process and the realization of that process is long enough). This geometric feature is exploited in one of the tests used to measure the quality of random number generators. For example, if a random number generator yields a sequence of numbers distributed uniformly in the interval0,1, then the density of the embedded points should be the same in any location of a two-dimensional square 0,1×0,1 for a sufficiently long sequence and any time delay greater than 1 [34]. It is clear that Di1δ,…,im−1δ must be invariant to time delays if only the embedded time series is random. We demonstrate the properties of Di1δ,…,im−1δ for the random time series distributed uniformly in the interval a,b,a≥0 and for the Gaussian time series. Of course, the invariance in respect to time delays does hold, but some features of Di1δ,…,im−1δ are interesting and are presented below.


#### 3.4.1. Embedding the Uniform Noise

Let us consider a discrete random variableX distributed uniformly in the interval a,b,a≥0. Initially, let us consider m=1 and n→∞. In fact, this situation cannot be considered an embedding because all values of the time series are mapped on a single one-dimensional axis without any time delays. However, the formal result reads
(19)limn→∞D∣m=1=limn→∞1n∑k=1nxk=a+b2.Now, let us consider the situation whenm→∞ (for any n):
(20)limm→∞Di1δ,…,im−1δ=limm→∞1n−Δm∑k=1n−Δxk2+xk+i12+xk+i1+i22+⋯+xk+Δ2=limm→∞1n−Δm∑k=1n−Δm⋅a2+ab+b23=a2+ab+b23.Thus, finally,(21)a+b2<Di1δ,…,im−1δ<a2+ab+b23,if the embedded random time series is distributed uniformly in the interval a,b,a≥0.


#### 3.4.2. Embedding the Gaussian Noise

Now, let us assume that a discrete random variableX is the Gaussian noise, with the mean value μ and standard deviation σ. Firstly, let us consider m=1 and n→∞:
(22)limn→∞D∣m=1=limn→∞1n∑k=1nxk=μ.Now, let us consider the case whenm→∞ (for any n):
(23)limm→∞Di1δ,…,im−1δ=limm→∞1n−Δm∑k=1n−Δxk2+xk+i12+xk+i1+i22+⋯+xk+Δ2=limm→∞1n−Δm∑k=1n−Δm⋅μ2+σ2=μ2+σ2.Thus, finally,(24)μ<Di1δ,…,im−1δ<μ2+σ2,if the embedded random time series is the Gaussian noise, with the mean μ and the standard deviation σ.


### 3.5. Computational Experiments: Optimal Embedding of Mackey-Glass Time Series

Equation (18) is a nonlinear integer optimization problem. The full sort method could be used for the identification of the optimal set of time delays—only if m, L, and n are not very large and computational resources are not strictly limited. Otherwise, soft computing techniques could be exploited for the determination of the near-optimal set time delays [30]. Let us consider Mackey-Glass time series [35] the discrete numerical solution to chaotic Mackey-Glass delay differential equation that reads
(25)dxtdt=a⋅xt−c1+x10t−c−bxt,where a=0.2, b=0.1, and c=17; the observation window is 0≤t≤1000 and the time step δ=1 (Figure 2). The standard FNN algorithm [36] yields (5) as the minimum embedding dimension for the investigated chaotic Mackey-Glass series.Figure 2
Mackey-Glass time series ata=0.2, b=0.1, and c=17.Now, we use (18) to determine the optimal uniform time lag in the interval 1≤τ≤100 (Figure 3). The maximum value of D is reached at τ=27:D27,27,27,27=0.2248 (Figure 3). Note that this value of D can be used as a lower bound for the integer programming problem when the near-optimal set of nonuniform time delays is sought using branch-and-bound [37] or other soft optimization algorithms.Figure 3
Plot ofDτ for the Mackey-Glass time series (uniform embedding in the 5-dimensional delay coordinate space). The maximum value of Dτ is attained at τ=27.The full sort method is used to determine the optimal time delay vector for nonuniform embedding. 1004 different combinations 1≤i1,i2,i3,i4≤100 are assessed in 8.14 hours on a two-kernel 2.20 GHz processor with 8 GB RAM. As mentioned previously, soft computing techniques can be effectively exploited for the determination of the near-optimal set of time delays. Finally, the optimal set of time delays reads 23,66,16,20; D23,66,16,20=0.2264. All possible planar projections of the embedded Mackey-Glass attractor are illustrated in Figure 4.Figure 4
All possible planar projections of the embedded Mackey-Glass attractor occupy the largest area (in average) atτ1=23, τ2=66, τ3=16, and τ4=20.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)


### 3.6. Synchronization Measure Based on a Geometric Approach to Attractor Embedding

The proposed technique for the determination of optimal time delays is based on a straightforward geometric approach. Such an approach can be effectively exploited for comparing geometric properties of attractors embedded from different time series. Moreover, we will demonstrate that such an approach is working well with short time series.The basic idea of the proposed synchronization feature is based neither on the direct comparison between geometric shapes of the reconstructed attractors nor on global characteristics of those attractors (like averaged maximum Lyapunov exponent). We will use the optimal time delay as a single numerical characteristic describing the geometry of the reconstructed attractor in a 2-dimensional delay coordinate space. The reasoning for such a measure is straightforward. It follows from the properties of the reconstructed attractor (Figure4). Such synchronization measure is much more general compared to synchronization measures used to detect identical, phase, amplitude envelope, lag, or even generalized synchronization. Two attractors are considered to be similar if their optimal time lags are the same. The higher the difference between optimal time lags, the more different these attractors are considered (from the geometrical point of view).We will demonstrate that such an approach for detecting synchronization between two different time series is highly feasible—especially if those time series are far from being similar.Let us consider two paradigmatic chaotic Rössler systems with diffusive coupling:(26)dX1dt=−Y1−Z1+εX2−X1,dY1dt=X1+a1Y1,dZ1dt=b+Z1X1−c,dX2dt=−Y2−Z2+εX1−X2,dY2dt=X2+a2Y2,dZ2dt=b+Z2X2−c,where a1=0.2, a2=0.3, b=0.2, and c=5.7; ε is the coupling parameter. The initial condition is set to X1=X2=1, Y1=Y2=1, and Z1=Z2=0.Attractors of both Rössler systems are illustrated in Figure5. The difference between the uncoupled attractors can be seen by a naked eye (Figures 5(a) and 5(b)). However, the attractors become similar as the coupling increases (Figures 5(e) and 5(f)).Figure 5
Coupled Rössler attractors (26) in the phase space: (a, b) uncoupled systems at ε=0; (c, d) coupled systems at ε=0.15; (e, f) coupled systems at ε=0.95 (a1=0.2, a2=0.3, b=0.2, and c=5.7).
(a)
(b)
(c)
(d)
(e)
(f)Analogously, attractors of coupled Rössler systems (26) are visualized in the embedding space (Figure 6). Note that a short observation window (400 scalar data points) is used to reconstruct the attractor in the planar delay-coordinate space. One time delay is used for every embedding; optimal values of time delays are depicted in Figure 6. Note that both the geometric shape of the embedded attractors and the optimal time delays become closer as the coupling parameter grows (Figure 6).Figure 6
Coupled Rössler attractors (26) in the planar delay coordinate space: (a, b) uncoupled systems at ε=0; (c, d) coupled systems at ε=0.15; (e, f) coupled systems at ε=0.95 (a1=0.2, a2=0.3, b=0.2, and c=5.7).
(a)
(b)
(c)
(d)
(e)
(f)It is also interesting to compare the maximal Lyapunov exponents calculated for the system (26) and the one in the embedding space. Averaged maximal Lyapunov exponents for the uncoupled Rössler systems are computed separately: λ1max=0.0764 for the first system and λ2max=0.1228 for the second system. However, the coupled Rössler systems (26) are considered a single system. The 6-dimensional Jacobian matrix of the coupled system yields (after averaging) a single maximal Lyapunov exponent: λmax=0.1094 at ε=0.15 and λmax=0.0524 at ε=0.95. It can be observed that the single maximal Lyapunov exponent of the coupled system at ε=0.15 approaches the average of λ1max and λ2max. But the maximal Lyapunov exponent decreases (the coupled system becomes less chaotic) at ε=0.95.The maximal averaged Lyapunov exponents computed directly from time seriesX1 and X2 in the embedded space are 0.0485 and 0.0541 for ε=0, 0.0367 and 0.0461 for ε=0.15, 0.0326 and 0.0344 for ε=0.95. Note that the maximal Lyapunov exponents computed directly from time series X1 and X2 are underestimated if compared to the maximal exponents for the coupled systems in (26). This effect is well known for synchronized differential equations [38]. Nevertheless, both methods demonstrate that the coupled systems become less chaotic when the coupling parameter is increased.Figure7(a) shows the evolution of X1 and X2 and the difference between X1 and X2 when two Rössler systems are uncoupled (the length of the observation window is 8000 points). Figure 7(b) shows X1, X2, and X1−X2 when ε=0.15, and Figure 7(c) corresponds to ε=0.95. Note that the scales of the vertical axes in the graphs in Figure 7(c) are different. Time series X1 and X2 become more similar as the coupling parameter ε increases; however, the difference X1−X2 remains chaotic even at ε=0.95 (Figure 7). The whole observation window is now split into 20 equal intervals; each interval contains 400 points (both for X1 and for X2). The optimal time lag (2-dimensional delay coordinate space) is identified for time series in each individual interval according to (18); 0≤τ≤200. Optimal time delays for X1 and X2 are represented by a solid and a dashed line correspondingly in Figures 8(a)–8(c); the correlation between two sequences of optimal time delays is ρ=0.19263, ρ=0.46533, ρ=0.99861 at ε=0, ε=0.15 and ε=0.95, respectively. It appears that the proposed algorithm based on the geometric approach to optimal attractor embedding is able to detect the synchronization between two chaotic series in an effective and efficient way. Note that the difference between X1 and X2 is not equal to zero at ε=0.95 (Figure 7). However, it appears that ρ (the correlation between two sequences of optimal time delays) is a good indicator of the geometrical synchronization between two sequences. It represents a geometric measure of synchronization between different time series. But instead of comparing geometric shapes of the embedded attractors, it compares sequences of optimal time lags in a predefined number of observation windows.Figure 7
Coupled Rössler systems (26) at different values of the coupling parameter ε. (a) X1, X2, and X1−X2 when two Rössler systems are uncoupled (ε=0). (b) X1, X2, and X1−X2 at ε=0.15. (c) X1, X2, and X1−X2 at ε=0.95.
(a)
(b)
(c)Figure 8
Correlation between sequences of optimal time lags serves as an indicator of the geometric synchronization between two time sequences. (a–c) Coupled Rössler systems without the additive noise. (d–f) Coupled Rössler systems with the additive Gaussian noise (the standard deviation equal to 2). (g–i) Coupled Rössler systems with the additive Gaussian noise (the standard deviation equal to 18).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Computational experiments are continued with the coupled Rössler systems with the additive Gaussian random noise. Note that different realizations of the Gaussian random noise are added toX1 and X2. Moreover, the noise is added after the integration of (26) and thus is not smoothed by the numerical integrator. Computational experiments are repeated for different values of the coupling coefficient (ε=0, ε=0.15, and ε=0.95) and with the different parameters of the additive Gaussian random noise. The Gaussian noise with the zero mean and the standard deviation equal to 2 is used in Figure 9, and that with the zero mean and the standard deviation equal to 18 is used in Figure 10.Figure 9
Coupled Rössler systems (26) at different values of the coupling parameter ε with the additive Gaussian noise. (a) X1, X2, and X1−X2 when two Rössler systems are uncoupled (ε=0). (b) X1, X2, and X1−X2 at ε=0.15. (c) X1, X2, and X1−X2 at ε=0.95. Different realizations of the Gaussian random noise (zero mean, standard deviation equal to 2) are added to X1 and X2 after the integration of (26).
(a)
(b)
(c)Figure 10
Coupled Rössler systems (26) at different values of the coupling parameter ε with the additive Gaussian noise. (a) X1, X2, and X1−X2 when two Rössler systems are uncoupled (ε=0). (b) X1, X2, and X1−X2 at ε=0.15. (c) X1, X2, and X1−X2 at ε=0.95. Different realizations of the Gaussian random noise (zero mean, standard deviation equal to 18) are added to X1 and X2 after the integration of (26).
(a)
(b)
(c)Again, the whole observation window is split into 20 equal intervals (each interval contains 400 points both forX1 and for X2). Optimal time lags for the coupled Rössler systems with the additive noise are represented by a solid and a dashed line correspondingly in Figures 8(d)–8(f) (for the standard deviation equal to 2) and in Figures 8(g)–8(i) (for the standard deviation equal to 18).The correlation between sequences of optimal time delays for the coupled Rössler systems with the additive noise (with standard deviation equal to 2) reads −0.39219, 0.36979, 0.79285 atε=0, ε=0.15, and ε=0.95, respectively (Figure 8). The differences between X1 and X2 are large, random, and unpredictable (Figure 9). Nevertheless, the proposed algorithm is capable of identifying the synchronization between two sequences even in such a complicated computational experiment.Analogously, the correlation between sequences of optimal time delays for the coupled Rössler systems with the additive noise (with standard deviation equal to 18) reads 0.094789, 0.23079, and 0.25094 atε=0, ε=0.15, and ε=0.95, respectively (Figure 8). Now, the proposed algorithm for the identification of the geometrical synchronization between X1 and X2 reaches the limits of applicability. Simply, the signal-noise ratio is so low that it is impossible to trace the presence of any synchronization between two time series.


### 3.7. Computational Experiments: Comparisons and Extensions

This paper presents a new technique for the quantification of geometric synchronization between different time series. This technique is based on the identification of optimal time lags which maximize the state space occupied by the embedded segments of scalar time series. A natural question is whether other existing methods for the determination of optimal time lags could be used instead.Computational experiments are repeated with the coupled Rössler systems, but the autocorrelation method for the determination of optimal time lags [5, 6] is used instead of the proposed method. All parameters of the computational setup are kept the same—including the coupling parameter, system parameters, and the length of the observation window (Figure 11).Figure 11
Correlation between sequences of optimal time lags identified by means of the autocorrelation method. (a–c) Coupled Rössler systems without the additive noise. (d–f) Coupled Rössler systems with the additive Gaussian noise (the standard deviation equal to 2). (g–i) Coupled systems with the additive Gaussian noise (the standard deviation equal to 18).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)The first observation is that the autocorrelation method yields incorrect correlation for the coupled Rössler systems atε=0.15; it appears that the correlation between uncoupled systems (Figures 11(a) and 11(b)) is higher than that between coupled systems (Figures 11(d) and 11(e)). Secondly, the autocorrelation method yields wrong results for the coupled Rössler systems with the additive noise—the correlation between the coupled systems with the additive Gaussian noise is higher than that between the coupled systems without the additive noise (Figures 11(b)–11(h)).Computational experiments are now repeated with the mutual information method for the determination of optimal time lags [4] (Figure 12). Again, correlations between the coupled systems with the additive noise (Figures 12(b)–12(e)) are higher those between the coupled systems without noise (Figures 12(c)–12(f)). Also, correlations between strongly coupled systems (Figures 12(e) and 12(f)) are lower than correlations between weakly coupled systems (Figures 12(h) and 12(i)).Figure 12
Correlation between sequences of optimal time lags identified by means of the mutual information method. (a–c) Coupled Rössler systems without the additive noise. (d–f) Coupled Rössler systems with the additive Gaussian noise (the standard deviation equal to 2). (g–i) Coupled systems with the additive Gaussian noise (the standard deviation equal to 18).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)These computational experiments demonstrate that the identification of optimal time lags based on the maximization of the state space occupied by the embedded segments of scalar time series is a perfect approach for the presented synchronization measure.Computational experiments are continued with a completely different system. The proposed synchronization measure is tested using a coupled system of periodically forced nonlinear pendulums:(27)d2Xdt2+bdXdt+sinX=f1cosωt+εY−X,d2Ydt2+bdYdt+sinY=f2cosωt+εX−Y,where X and Y are angular coordinates of the coupled systems, b is the damping coefficient, f1 and f2 are the external forcing amplitudes of both systems, and ω is the frequency. A driven damped nonlinear pendulum is a paradigmatic model in nonlinear dynamics; the model yields chaotic behavior at ω=2/3, b=1, and f=2.048 [39]. The parameters of the coupled model are set as follows: ω=2/3, b=1, f1=2.048, and f2=2.049.The proposed synchronization measureρ for (27) is presented in Figure 13. At first, the coupled system of nonlinear pendulums is considered without the additive noise (Figures 13(a)–13(c)). The synchronization measure for the uncoupled system (ε=0) is ρ=−0.21944. The coupling yields ρ=0.56745 at ε=0.04 and ρ=0.98281 (complete synchronization) at ε=0.095.Figure 13
Correlation between sequences of optimal time lags for the coupled nonlinear pendulums (27). (a–c) The coupled system without the additive noise. (d–f) The coupled system with the additive Gaussian noise (the standard deviation equal to 0.6). (g–i) The coupled system with the additive Gaussian noise (the standard deviation equal to 5.6).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Computational experiments are continued by adding the additive Gaussian noise. The proposed method is still able to detect synchronization between the systems when the standard deviation of the additive Gaussian noise is equal to 0.6 (Figures13(d)–13(f)). Finally, the method is not able to detect synchronization when the standard deviation of the additive noise is set to 5.6 (Figures 13(h)–13(i)).The ability of the proposed technique to identify the coupling between two chaotic systems (even in the presence of a moderate additive noise) could be explained by the nature of the algorithm which is used to identify optimal time lags. We directly assess the geometric properties of the embedded attractor in the delay coordinate space. The proposed technique is indeed applicable to short time series; the length of the observation window must be sufficient to reconstruct at least one full orbit in the reconstructed delay coordinate space (Figure6). This is in contrast with the autocorrelation method or the mutual information method requiring much longer observation windows [4–6].


## 3.1. Properties ofDτ1,τ2,…,τm−1: Embedding a Harmonic Function

Let us consider that the embedding is uniformτ1=τ2=⋯=τm−1=τ;T→∞ and the embedded function is a harmonic function ft=asinω0t+ϕ, where a, ω0, and ϕ are the amplitude, the circular frequency, and the phase, accordingly. Then,
(11)Dτ=amlimT→∞1T∑k=0m−1sin2ω0t+kτ+ϕdt.The change of variablesω0t+ϕ=x yields
(12)Dτ=aω0mlimT→∞1T∫ϕω0T+ϕ∑k=0m−1sin2x+ω0kτdx=aω0m12π∫02π∑k=0m−1sin2x+ω0kτdx.Now,(13)Dτ=aω0m12π∫02πsin2x∑k=0m−112+12cos2ω0kτ+12sin2x∑k=0m−1sin2ω0kτ+cos2x∑k=0m−112−12cos2ω0kτ1/2dx=aω02m12π∫02πm−cos2x∑k=0m−1cos2ω0kτ+sin2x∑k=0m−1sin2ω0kτdx.Let us denotez=ei2ω0τ. Then,
(14)Dτ=aω02m12π∫02πm−cos2x∑k=0m−1Rezk+sin2x∑k=0m−1Imzkdx=aω02m12π∫02πm−cos2xRe1−zm1−z+sin2xIm1−zm1−zdx.Note that1−zm/1−z=0 when τ=π/mω0. Thus,
(15)Dπmω0=aω02m12π∫02πmdx=a2ω0is the maximal value for the harmonic function. On the opposite, the pole z=1 is a removable singularity at τ=π/ω0 and
(16)Dπω0=aω0212π∫02π1−cos2xdx=2aω0πis the minimal value for the harmonic function.


## 3.2. The Comparison betweenDτ and Sτ

We continue with the harmonic functionft=asinω0t+ϕ and uniform embedding. Then, Aω=δω−ω0 and (9) yields Sτ=π/2Qτ,ω0, where δ is the Dirac delta function. Without the loss of generality, we assume that a=1 and ω0=1. Then, Dmin=2/π≈0.637; Dmax=1/2≈0.707. Graphs of Dτ and Sτ are illustrated in Figure 1.Figure 1
The comparison of functionsDτ (solid line) and Sτ (dashed line). The parameter m is set to 2, 3, 4, and 20 in (a), (b), (c), and (d), respectively.
(a)
(b)
(c)
(d)The first observation is that graphs ofDτ and Sτ are similar at m=2 (Figure 1(a)). The best embedding of the harmonic signal into the 2-dimensional delay coordinate space is produced at τ=π/2+πk, k∈ℤ (then the harmonic function is embedded into a circle). The harmonic function is embedded into a line interval on one of the diagonals of the frame at τ=πk, k∈ℤ. The difference is that now Sπk=0, but Dπk=2/π (note that the scale of the y-axis for Sτ is shown on the left, but the scale for Dτ is shown on the right side of Figure 1(a)).The differences betweenDτ and Sτ are more clear at m=3 (Figure 1(b)). The harmonic function cannot be mapped into circles in all 3 projections at a fixed τ; therefore, the maximum value of Sτ at m=3 is lower than that at m=2 (note that the scales of the axes are the same for different m in Figure 1). But Dmax remains the same and does not depend on the dimension of the delay coordinate space. Both Dτ and Sτ yield the best embedding at π/3l, where l∈ℤ/3k, k∈ℤ.The maximum value ofSτ at m=4 becomes smaller than the one at m=3 (Figure 1(c)). But now a major difference between Dτ and Sτ can be observed in terms of the best embedding. Sτ yields the best embedding at ±π/4+πk, but Dτ yields the best embedding at π/4l, where l∈ℤ/4k, k∈ℤ. Clearly, Sτ misses the time delay τ=π/2+πk. That can be explained by the fact that Sτ measures the average area of the mapped attractor in all possible plane projections of the m-dimensional delay coordinate space. On the contrary, Dτ directly measures the geometrical volume occupied by the attractor without considering any projections of that attractor.The maximum value ofSτ tends to 1 when m→∞ (Figure 1(d)). Note that Sτ1,…,τm−1=1 when the embedded time series is the white noise [11]. However, as mentioned previously, Dmax is an invariant value in respect of m. What is more important, Dτ yields the best embedding at π/ml, where l∈ℤ/mk, k∈ℤ, which is in stark contrast to Sτ (Figure 1(d)).


## 3.3. Embedding a Discrete Time Series

Let us consider a finite scalar discrete time seriesx1,x2,…,xn, where n is the number of observation points uniformly spread in the time axis. Let the sampling rate of the time series (the time interval between adjacent observations) be δ. Then, all time delays are multiple values of δ:τk=ikδ,ik∈ℕ,k=1,2,…,m−1. Let us denote Δ=∑k=1m−1ik. Then, (10) reads
(17)Di1δ,i2δ,…,im−1δ=1n−Δm∑k=1n−Δxk2+xk+i12+xk+i1+i22+⋯+xk+Δ2.Note thatδ does not play any role on the right-hand side of (17)—because the definite integral in (10) is replaced by the limit sum in (17). Now, the problem of the identification of the optimal set of time delays is formulated as follows:
(18)max1≤i1,…,im−1≤L1n−Δm∑k=1n−Δxk2+xk+i12+xk+i1+i22+⋯+xk+Δ2,where L is the upper limit for every individual time delay. The total embedding window (the time delay between the first and the last time delay coordinates) is Δδ, and it could be also used for the estimation of L. It is clear that n must be greater than Δ. However, n must be at least 10 times greater than m−1L in order to produce a meaningful optimization result (18).


## 3.4. Properties ofDτ1,τ2,…,τm−1: Embedding the Random Noise

It is well known that embedding of the random noise is invariant to time delays. In other words, the geometric shape of the reconstructed attractor does not depend on time delays (if only the embedded time series is an ergodic random process and the realization of that process is long enough). This geometric feature is exploited in one of the tests used to measure the quality of random number generators. For example, if a random number generator yields a sequence of numbers distributed uniformly in the interval0,1, then the density of the embedded points should be the same in any location of a two-dimensional square 0,1×0,1 for a sufficiently long sequence and any time delay greater than 1 [34]. It is clear that Di1δ,…,im−1δ must be invariant to time delays if only the embedded time series is random. We demonstrate the properties of Di1δ,…,im−1δ for the random time series distributed uniformly in the interval a,b,a≥0 and for the Gaussian time series. Of course, the invariance in respect to time delays does hold, but some features of Di1δ,…,im−1δ are interesting and are presented below.


### 3.4.1. Embedding the Uniform Noise

Let us consider a discrete random variableX distributed uniformly in the interval a,b,a≥0. Initially, let us consider m=1 and n→∞. In fact, this situation cannot be considered an embedding because all values of the time series are mapped on a single one-dimensional axis without any time delays. However, the formal result reads
(19)limn→∞D∣m=1=limn→∞1n∑k=1nxk=a+b2.Now, let us consider the situation whenm→∞ (for any n):
(20)limm→∞Di1δ,…,im−1δ=limm→∞1n−Δm∑k=1n−Δxk2+xk+i12+xk+i1+i22+⋯+xk+Δ2=limm→∞1n−Δm∑k=1n−Δm⋅a2+ab+b23=a2+ab+b23.Thus, finally,(21)a+b2<Di1δ,…,im−1δ<a2+ab+b23,if the embedded random time series is distributed uniformly in the interval a,b,a≥0.


### 3.4.2. Embedding the Gaussian Noise

Now, let us assume that a discrete random variableX is the Gaussian noise, with the mean value μ and standard deviation σ. Firstly, let us consider m=1 and n→∞:
(22)limn→∞D∣m=1=limn→∞1n∑k=1nxk=μ.Now, let us consider the case whenm→∞ (for any n):
(23)limm→∞Di1δ,…,im−1δ=limm→∞1n−Δm∑k=1n−Δxk2+xk+i12+xk+i1+i22+⋯+xk+Δ2=limm→∞1n−Δm∑k=1n−Δm⋅μ2+σ2=μ2+σ2.Thus, finally,(24)μ<Di1δ,…,im−1δ<μ2+σ2,if the embedded random time series is the Gaussian noise, with the mean μ and the standard deviation σ.


## 3.4.1. Embedding the Uniform Noise

Let us consider a discrete random variableX distributed uniformly in the interval a,b,a≥0. Initially, let us consider m=1 and n→∞. In fact, this situation cannot be considered an embedding because all values of the time series are mapped on a single one-dimensional axis without any time delays. However, the formal result reads
(19)limn→∞D∣m=1=limn→∞1n∑k=1nxk=a+b2.Now, let us consider the situation whenm→∞ (for any n):
(20)limm→∞Di1δ,…,im−1δ=limm→∞1n−Δm∑k=1n−Δxk2+xk+i12+xk+i1+i22+⋯+xk+Δ2=limm→∞1n−Δm∑k=1n−Δm⋅a2+ab+b23=a2+ab+b23.Thus, finally,(21)a+b2<Di1δ,…,im−1δ<a2+ab+b23,if the embedded random time series is distributed uniformly in the interval a,b,a≥0.


## 3.4.2. Embedding the Gaussian Noise

Now, let us assume that a discrete random variableX is the Gaussian noise, with the mean value μ and standard deviation σ. Firstly, let us consider m=1 and n→∞:
(22)limn→∞D∣m=1=limn→∞1n∑k=1nxk=μ.Now, let us consider the case whenm→∞ (for any n):
(23)limm→∞Di1δ,…,im−1δ=limm→∞1n−Δm∑k=1n−Δxk2+xk+i12+xk+i1+i22+⋯+xk+Δ2=limm→∞1n−Δm∑k=1n−Δm⋅μ2+σ2=μ2+σ2.Thus, finally,(24)μ<Di1δ,…,im−1δ<μ2+σ2,if the embedded random time series is the Gaussian noise, with the mean μ and the standard deviation σ.


## 3.5. Computational Experiments: Optimal Embedding of Mackey-Glass Time Series

Equation (18) is a nonlinear integer optimization problem. The full sort method could be used for the identification of the optimal set of time delays—only if m, L, and n are not very large and computational resources are not strictly limited. Otherwise, soft computing techniques could be exploited for the determination of the near-optimal set time delays [30]. Let us consider Mackey-Glass time series [35] the discrete numerical solution to chaotic Mackey-Glass delay differential equation that reads
(25)dxtdt=a⋅xt−c1+x10t−c−bxt,where a=0.2, b=0.1, and c=17; the observation window is 0≤t≤1000 and the time step δ=1 (Figure 2). The standard FNN algorithm [36] yields (5) as the minimum embedding dimension for the investigated chaotic Mackey-Glass series.Figure 2
Mackey-Glass time series ata=0.2, b=0.1, and c=17.Now, we use (18) to determine the optimal uniform time lag in the interval 1≤τ≤100 (Figure 3). The maximum value of D is reached at τ=27:D27,27,27,27=0.2248 (Figure 3). Note that this value of D can be used as a lower bound for the integer programming problem when the near-optimal set of nonuniform time delays is sought using branch-and-bound [37] or other soft optimization algorithms.Figure 3
Plot ofDτ for the Mackey-Glass time series (uniform embedding in the 5-dimensional delay coordinate space). The maximum value of Dτ is attained at τ=27.The full sort method is used to determine the optimal time delay vector for nonuniform embedding. 1004 different combinations 1≤i1,i2,i3,i4≤100 are assessed in 8.14 hours on a two-kernel 2.20 GHz processor with 8 GB RAM. As mentioned previously, soft computing techniques can be effectively exploited for the determination of the near-optimal set of time delays. Finally, the optimal set of time delays reads 23,66,16,20; D23,66,16,20=0.2264. All possible planar projections of the embedded Mackey-Glass attractor are illustrated in Figure 4.Figure 4
All possible planar projections of the embedded Mackey-Glass attractor occupy the largest area (in average) atτ1=23, τ2=66, τ3=16, and τ4=20.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)


## 3.6. Synchronization Measure Based on a Geometric Approach to Attractor Embedding

The proposed technique for the determination of optimal time delays is based on a straightforward geometric approach. Such an approach can be effectively exploited for comparing geometric properties of attractors embedded from different time series. Moreover, we will demonstrate that such an approach is working well with short time series.The basic idea of the proposed synchronization feature is based neither on the direct comparison between geometric shapes of the reconstructed attractors nor on global characteristics of those attractors (like averaged maximum Lyapunov exponent). We will use the optimal time delay as a single numerical characteristic describing the geometry of the reconstructed attractor in a 2-dimensional delay coordinate space. The reasoning for such a measure is straightforward. It follows from the properties of the reconstructed attractor (Figure4). Such synchronization measure is much more general compared to synchronization measures used to detect identical, phase, amplitude envelope, lag, or even generalized synchronization. Two attractors are considered to be similar if their optimal time lags are the same. The higher the difference between optimal time lags, the more different these attractors are considered (from the geometrical point of view).We will demonstrate that such an approach for detecting synchronization between two different time series is highly feasible—especially if those time series are far from being similar.Let us consider two paradigmatic chaotic Rössler systems with diffusive coupling:(26)dX1dt=−Y1−Z1+εX2−X1,dY1dt=X1+a1Y1,dZ1dt=b+Z1X1−c,dX2dt=−Y2−Z2+εX1−X2,dY2dt=X2+a2Y2,dZ2dt=b+Z2X2−c,where a1=0.2, a2=0.3, b=0.2, and c=5.7; ε is the coupling parameter. The initial condition is set to X1=X2=1, Y1=Y2=1, and Z1=Z2=0.Attractors of both Rössler systems are illustrated in Figure5. The difference between the uncoupled attractors can be seen by a naked eye (Figures 5(a) and 5(b)). However, the attractors become similar as the coupling increases (Figures 5(e) and 5(f)).Figure 5
Coupled Rössler attractors (26) in the phase space: (a, b) uncoupled systems at ε=0; (c, d) coupled systems at ε=0.15; (e, f) coupled systems at ε=0.95 (a1=0.2, a2=0.3, b=0.2, and c=5.7).
(a)
(b)
(c)
(d)
(e)
(f)Analogously, attractors of coupled Rössler systems (26) are visualized in the embedding space (Figure 6). Note that a short observation window (400 scalar data points) is used to reconstruct the attractor in the planar delay-coordinate space. One time delay is used for every embedding; optimal values of time delays are depicted in Figure 6. Note that both the geometric shape of the embedded attractors and the optimal time delays become closer as the coupling parameter grows (Figure 6).Figure 6
Coupled Rössler attractors (26) in the planar delay coordinate space: (a, b) uncoupled systems at ε=0; (c, d) coupled systems at ε=0.15; (e, f) coupled systems at ε=0.95 (a1=0.2, a2=0.3, b=0.2, and c=5.7).
(a)
(b)
(c)
(d)
(e)
(f)It is also interesting to compare the maximal Lyapunov exponents calculated for the system (26) and the one in the embedding space. Averaged maximal Lyapunov exponents for the uncoupled Rössler systems are computed separately: λ1max=0.0764 for the first system and λ2max=0.1228 for the second system. However, the coupled Rössler systems (26) are considered a single system. The 6-dimensional Jacobian matrix of the coupled system yields (after averaging) a single maximal Lyapunov exponent: λmax=0.1094 at ε=0.15 and λmax=0.0524 at ε=0.95. It can be observed that the single maximal Lyapunov exponent of the coupled system at ε=0.15 approaches the average of λ1max and λ2max. But the maximal Lyapunov exponent decreases (the coupled system becomes less chaotic) at ε=0.95.The maximal averaged Lyapunov exponents computed directly from time seriesX1 and X2 in the embedded space are 0.0485 and 0.0541 for ε=0, 0.0367 and 0.0461 for ε=0.15, 0.0326 and 0.0344 for ε=0.95. Note that the maximal Lyapunov exponents computed directly from time series X1 and X2 are underestimated if compared to the maximal exponents for the coupled systems in (26). This effect is well known for synchronized differential equations [38]. Nevertheless, both methods demonstrate that the coupled systems become less chaotic when the coupling parameter is increased.Figure7(a) shows the evolution of X1 and X2 and the difference between X1 and X2 when two Rössler systems are uncoupled (the length of the observation window is 8000 points). Figure 7(b) shows X1, X2, and X1−X2 when ε=0.15, and Figure 7(c) corresponds to ε=0.95. Note that the scales of the vertical axes in the graphs in Figure 7(c) are different. Time series X1 and X2 become more similar as the coupling parameter ε increases; however, the difference X1−X2 remains chaotic even at ε=0.95 (Figure 7). The whole observation window is now split into 20 equal intervals; each interval contains 400 points (both for X1 and for X2). The optimal time lag (2-dimensional delay coordinate space) is identified for time series in each individual interval according to (18); 0≤τ≤200. Optimal time delays for X1 and X2 are represented by a solid and a dashed line correspondingly in Figures 8(a)–8(c); the correlation between two sequences of optimal time delays is ρ=0.19263, ρ=0.46533, ρ=0.99861 at ε=0, ε=0.15 and ε=0.95, respectively. It appears that the proposed algorithm based on the geometric approach to optimal attractor embedding is able to detect the synchronization between two chaotic series in an effective and efficient way. Note that the difference between X1 and X2 is not equal to zero at ε=0.95 (Figure 7). However, it appears that ρ (the correlation between two sequences of optimal time delays) is a good indicator of the geometrical synchronization between two sequences. It represents a geometric measure of synchronization between different time series. But instead of comparing geometric shapes of the embedded attractors, it compares sequences of optimal time lags in a predefined number of observation windows.Figure 7
Coupled Rössler systems (26) at different values of the coupling parameter ε. (a) X1, X2, and X1−X2 when two Rössler systems are uncoupled (ε=0). (b) X1, X2, and X1−X2 at ε=0.15. (c) X1, X2, and X1−X2 at ε=0.95.
(a)
(b)
(c)Figure 8
Correlation between sequences of optimal time lags serves as an indicator of the geometric synchronization between two time sequences. (a–c) Coupled Rössler systems without the additive noise. (d–f) Coupled Rössler systems with the additive Gaussian noise (the standard deviation equal to 2). (g–i) Coupled Rössler systems with the additive Gaussian noise (the standard deviation equal to 18).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Computational experiments are continued with the coupled Rössler systems with the additive Gaussian random noise. Note that different realizations of the Gaussian random noise are added toX1 and X2. Moreover, the noise is added after the integration of (26) and thus is not smoothed by the numerical integrator. Computational experiments are repeated for different values of the coupling coefficient (ε=0, ε=0.15, and ε=0.95) and with the different parameters of the additive Gaussian random noise. The Gaussian noise with the zero mean and the standard deviation equal to 2 is used in Figure 9, and that with the zero mean and the standard deviation equal to 18 is used in Figure 10.Figure 9
Coupled Rössler systems (26) at different values of the coupling parameter ε with the additive Gaussian noise. (a) X1, X2, and X1−X2 when two Rössler systems are uncoupled (ε=0). (b) X1, X2, and X1−X2 at ε=0.15. (c) X1, X2, and X1−X2 at ε=0.95. Different realizations of the Gaussian random noise (zero mean, standard deviation equal to 2) are added to X1 and X2 after the integration of (26).
(a)
(b)
(c)Figure 10
Coupled Rössler systems (26) at different values of the coupling parameter ε with the additive Gaussian noise. (a) X1, X2, and X1−X2 when two Rössler systems are uncoupled (ε=0). (b) X1, X2, and X1−X2 at ε=0.15. (c) X1, X2, and X1−X2 at ε=0.95. Different realizations of the Gaussian random noise (zero mean, standard deviation equal to 18) are added to X1 and X2 after the integration of (26).
(a)
(b)
(c)Again, the whole observation window is split into 20 equal intervals (each interval contains 400 points both forX1 and for X2). Optimal time lags for the coupled Rössler systems with the additive noise are represented by a solid and a dashed line correspondingly in Figures 8(d)–8(f) (for the standard deviation equal to 2) and in Figures 8(g)–8(i) (for the standard deviation equal to 18).The correlation between sequences of optimal time delays for the coupled Rössler systems with the additive noise (with standard deviation equal to 2) reads −0.39219, 0.36979, 0.79285 atε=0, ε=0.15, and ε=0.95, respectively (Figure 8). The differences between X1 and X2 are large, random, and unpredictable (Figure 9). Nevertheless, the proposed algorithm is capable of identifying the synchronization between two sequences even in such a complicated computational experiment.Analogously, the correlation between sequences of optimal time delays for the coupled Rössler systems with the additive noise (with standard deviation equal to 18) reads 0.094789, 0.23079, and 0.25094 atε=0, ε=0.15, and ε=0.95, respectively (Figure 8). Now, the proposed algorithm for the identification of the geometrical synchronization between X1 and X2 reaches the limits of applicability. Simply, the signal-noise ratio is so low that it is impossible to trace the presence of any synchronization between two time series.


## 3.7. Computational Experiments: Comparisons and Extensions

This paper presents a new technique for the quantification of geometric synchronization between different time series. This technique is based on the identification of optimal time lags which maximize the state space occupied by the embedded segments of scalar time series. A natural question is whether other existing methods for the determination of optimal time lags could be used instead.Computational experiments are repeated with the coupled Rössler systems, but the autocorrelation method for the determination of optimal time lags [5, 6] is used instead of the proposed method. All parameters of the computational setup are kept the same—including the coupling parameter, system parameters, and the length of the observation window (Figure 11).Figure 11
Correlation between sequences of optimal time lags identified by means of the autocorrelation method. (a–c) Coupled Rössler systems without the additive noise. (d–f) Coupled Rössler systems with the additive Gaussian noise (the standard deviation equal to 2). (g–i) Coupled systems with the additive Gaussian noise (the standard deviation equal to 18).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)The first observation is that the autocorrelation method yields incorrect correlation for the coupled Rössler systems atε=0.15; it appears that the correlation between uncoupled systems (Figures 11(a) and 11(b)) is higher than that between coupled systems (Figures 11(d) and 11(e)). Secondly, the autocorrelation method yields wrong results for the coupled Rössler systems with the additive noise—the correlation between the coupled systems with the additive Gaussian noise is higher than that between the coupled systems without the additive noise (Figures 11(b)–11(h)).Computational experiments are now repeated with the mutual information method for the determination of optimal time lags [4] (Figure 12). Again, correlations between the coupled systems with the additive noise (Figures 12(b)–12(e)) are higher those between the coupled systems without noise (Figures 12(c)–12(f)). Also, correlations between strongly coupled systems (Figures 12(e) and 12(f)) are lower than correlations between weakly coupled systems (Figures 12(h) and 12(i)).Figure 12
Correlation between sequences of optimal time lags identified by means of the mutual information method. (a–c) Coupled Rössler systems without the additive noise. (d–f) Coupled Rössler systems with the additive Gaussian noise (the standard deviation equal to 2). (g–i) Coupled systems with the additive Gaussian noise (the standard deviation equal to 18).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)These computational experiments demonstrate that the identification of optimal time lags based on the maximization of the state space occupied by the embedded segments of scalar time series is a perfect approach for the presented synchronization measure.Computational experiments are continued with a completely different system. The proposed synchronization measure is tested using a coupled system of periodically forced nonlinear pendulums:(27)d2Xdt2+bdXdt+sinX=f1cosωt+εY−X,d2Ydt2+bdYdt+sinY=f2cosωt+εX−Y,where X and Y are angular coordinates of the coupled systems, b is the damping coefficient, f1 and f2 are the external forcing amplitudes of both systems, and ω is the frequency. A driven damped nonlinear pendulum is a paradigmatic model in nonlinear dynamics; the model yields chaotic behavior at ω=2/3, b=1, and f=2.048 [39]. The parameters of the coupled model are set as follows: ω=2/3, b=1, f1=2.048, and f2=2.049.The proposed synchronization measureρ for (27) is presented in Figure 13. At first, the coupled system of nonlinear pendulums is considered without the additive noise (Figures 13(a)–13(c)). The synchronization measure for the uncoupled system (ε=0) is ρ=−0.21944. The coupling yields ρ=0.56745 at ε=0.04 and ρ=0.98281 (complete synchronization) at ε=0.095.Figure 13
Correlation between sequences of optimal time lags for the coupled nonlinear pendulums (27). (a–c) The coupled system without the additive noise. (d–f) The coupled system with the additive Gaussian noise (the standard deviation equal to 0.6). (g–i) The coupled system with the additive Gaussian noise (the standard deviation equal to 5.6).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Computational experiments are continued by adding the additive Gaussian noise. The proposed method is still able to detect synchronization between the systems when the standard deviation of the additive Gaussian noise is equal to 0.6 (Figures13(d)–13(f)). Finally, the method is not able to detect synchronization when the standard deviation of the additive noise is set to 5.6 (Figures 13(h)–13(i)).The ability of the proposed technique to identify the coupling between two chaotic systems (even in the presence of a moderate additive noise) could be explained by the nature of the algorithm which is used to identify optimal time lags. We directly assess the geometric properties of the embedded attractor in the delay coordinate space. The proposed technique is indeed applicable to short time series; the length of the observation window must be sufficient to reconstruct at least one full orbit in the reconstructed delay coordinate space (Figure6). This is in contrast with the autocorrelation method or the mutual information method requiring much longer observation windows [4–6].


## 4. Concluding Remarks

This paper presents a novel approach towards geometric synchronization of chaotic time series. This approach is based on the correlation between sets of time lags used for the reconstruction of attractors from different time series in delay coordinate spaces. An efficient algorithm is proposed for the identification of optimal time delays in finite observation windows. Computational experiments with the Mackey-Glass time series, coupled Rössler systems, and coupled chaotic pendulums are used to demonstrate the efficacy of the presented techniques.It is natural to expect that every method has its strengths and weaknesses. The No-Free-Lunch theorems [40] yield a conclusion that it is always possible to find a problem which will defeat any algorithm. The proposed method does not compare the shapes of the embedded attractors; instead, it compares optimal time delays resulting in a maximal expansion of the embedded attractors in the delay coordinate space. Therefore, it should be possible to design such two artificial time series which are not synchronized at all, but the presented technique would yield a positive synchronization measure (if only optimal time lags would coincide with appropriate observation windows). It is completely natural that every algorithm does have its bounds of applicability. In the present case, one should be cautious in testing the synchronization between time series of a completely different origin.The proposed synchronization measure is based on a geometric approach to attractor embedding using finite observation windows. However, a short observation window is not the main strength of the presented method. The proposed technique uses a direct geometrical estimate of the embedded attractor in the delay coordinate space. In other words, the presented method compares the geometric features of the embedded attractors (in terms of the optimal delays). Such generalized geometric approach opens new possibilities for the design of new geometric synchronization measures between different time series.



---

*Source: 8259496-2018-08-08.xml*