On the Use of Interval Extensions to Estimate the Largest Lyapunov Exponent from Chaotic Data

Amethod to estimate the (positive) largest Lyapunov exponent (LLE) from data using interval extensions is proposed.Themethod differs from the ones available in the literature in its simplicity since it is only based on three rather simple steps. Firstly, a polynomial NARMAX is used to identify amodel from the data under investigation. Secondly, interval extensions, which can be easily extracted from the identified model, are used to calculate the lower bound error. Finally, a simple linear fit to the logarithm of lower bound error is obtained and then the LLE is retrieved from it as the third step. To illustrate the proposed method, the LLE is calculated for the following well-known benchmarks: sine map, Rössler Equations, and Mackey-Glass Equations from identified models given in the literature and also from two identified NARMAXmodels: a chaotic jerk circuit and the tent map. In the latter, a Gaussian noise has been added to show the robustness of the proposed method.


Introduction
One of the most applied techniques to confirm whether a system is chaotic or not is based on the estimation of the Lyapunov exponents [1].The calculation of Lyapunov exponents is usually based on the average divergence or convergence of close trajectories along certain directions in state space.In chaotic systems, two trajectories separate exponentially with time despite having very similar initial conditions.Since the work of [2] several numerical methods to estimate Lyapunov exponents were proposed [3][4][5][6][7][8], just to mention a few.A comparison between estimation methods can be found in [9], for instance.The authors in [5] developed the first practical algorithm to calculate the Lyapunov exponents by estimating the growth of the corresponding set of vectors as the system evolves.This method allows the estimation of the complete spectrum of Lyapunov exponents.Amongst these exponents, the (positive) largest Lyapunov exponent (LLE) is the exponent considered to be the main reason for the separation rate.Therefore the estimation of such an exponent is used to build up the chaotic nature of the data under scrutiny.The authors in [1,7] independently used statistical properties of the local divergence rates of nearby trajectories of the systems under investigation to estimate the LLE.These properties come in handy when the goal is to get high exactness of numerical estimates of the LLE, as shown in the work [10].Recent and different applications using the LLE can be found in [11][12][13][14].
In [15] a straightforward method to compute the LLE using interval extensions and the original equations of motion was proposed.The exponent is estimated from the slope of the line derived from the lower bound error when considering two interval extensions of the original system [16,17].The authors of [15] demonstrated that the method is quick and simple to be used and could be considered as an alternative to other algorithms available in the literature.However, it requires the original equations of motion of the system and it can not be applied directly to time series.In order to overcome this limitation, a system identification methodology based on the well-known polynomial NAR-MAX (Nonlinear AutoRegressive Moving Average model with eXogenous inputs) [18] approach is introduced.Such an approach (system representation) has been shown to adequately represent a large number of chaotic systems [18,19].A similar method using the same system representation was proposed in [20] but, in their procedure, the Jacobian matrix must be found and its values estimated.For polynomial NARMAX it was argued that the Jacobian matrix is easily obtained but that may not be the case when other types of representation, such as Neural Networks [20], are considered.Using the interval extensions, it is shown that it is possible to put together the advantages pointed out in [20] without the need to calculate the Jacobian matrix to estimate the LLE.Moreover, as NARMAX is based on a polynomial structure, interval extensions are easily generated by simple arithmetic manipulations, such as commutative and associative.In addition to that, as indicated in [15], the proposed method requires similar or less points to achieve convergence when compared to other methods available in the literature.
The rest of this paper is laid out as follows.Preliminary concepts are introduced in Section 2 and it gives the necessary background material to understand the main idea of the work.Then, the proposed method is presented in Section 3. Section 4 presents the results obtained by the proposed approach when applying it to three different systems that have been identified in the literature: sine map, Rössler Equations, and Mackey-Glass Equations.Moreover, two identified polynomial NARMAX models have been considered: a chaotic jerk circuit [21] and the tent map.In the end, a Gaussian noise has been added to demonstrate the robustness of the proposed technique.Finally, the conclusions are given in Section 5.

Preliminary Concepts
In this section a brief review on the calculation of the LLE is given along with the recent result on the lower bound error and the interval extensions.

The LLE Calculation.
To review the basic ideas laid out in [7], consider the following equation: where U  is the set of all other delay vectors in an neighbourhood of the vector   (data from trajectories of the system under investigation) and |U  | is the number of elements in U  .The initial and final point of each window over the time series is given by  and , respectively.The LLE can be estimated by searching for a linear scaling in plot (Δ) versus Δ.In order to find the linear scaling, it is necessary to tune the parameter , to find the neighbourhood, and to define the number of elements in U  .This will not be necessary for the method proposed in Section 3; however there is the need of finding the equations of motion of the system under investigation and this will be accomplished by using a system identification methodology as reviewed in the next section.

The Polynomial NARMAX.
A NARMAX model can be written as follows [22]: where , , and  are, respectively, the output, input, and noise terms at the discrete time  ∈ N. The parameters   ,   , and   are the maximum lag considered for output, input, and noise, respectively.Noise terms are frequently included in the parameter estimation process to avoid bias in the estimates.
Here F ℓ [⋅] is assumed to be a polynomial with nonlinearity degree ℓ ∈ Z + .Nonlinear systems are frequently modelled using particular cases of the polynomial NARMAX, such as NAR, NARX, and NARMA [16,[23][24][25].In this paper, they all are treated under the name polynomial NARMAX.Further details of this mathematical representation can be found in [18].The great advantage of using the polynomial NARMAX representation is that the problem of finding the relevant terms and the estimates of their coefficients can be cast into a linear regression problem and therefore the well-known least squares method can be applied.In the eighties, Billings [18] devised an efficient method to take the advantage of the classical methods for solving the least squares problem to include the structure selection using the so-called ERR (Error Reduction Ratio) which attempts to measure the contribution of each term to the overall variance of the output signal.Such a method was called OLS (Orthogonal Least Squares method) and has been successfully used to obtain dynamically valid models from data generated by nonlinear systems behaving chaotically (see, e.g., [26,27]).Additionally, it is straightforward to generate interval extensions from the polynomial NARMAX, as seen in Definition 1.

The Lower Bound
Error.This subsection is an adaptation from [16,17] on the lower bound error applied to continuous nonlinear systems.
Let  ∈ N, a metric space  ⊂ R; the relation where  :  → , is a recursive function or a map of a state space into itself and   denotes the state at the discrete time .The sequence {  } obtained by iterating (3) starting from an initial condition  0 is called the orbit of  0 [28,29].
Definition 1 (see [16]).An interval extension of  is an interval valued function  of an interval variable , with the property where an interval is a closed set of a real number  ∈ R, where Equations ( 6) are mathematically equivalent but they present different sequence of their basic arithmetic operations.On floating point standard neither commutative nor distributive properties hold [30,31].Therefore, (6) may exhibit different results when a computer-based approach is used to estimate their values.We have exploited these features in this work to estimate the LLE.Definition 3. Let  ∈ N represent a pseudo-orbit, which is defined by an initial condition, an interval extension of , some specific hardware and software, numerical precision standard, and discretization scheme.A pseudo-orbit is an approximation of an orbit and can be represented by { x, } or x,0 , x,1 , x,2 , . .., such that where  , ∈ R is the error and  , ≥ 0.
A pseudo-orbit characterises an interval where the genuine orbit lies.Henceforth an interval related to each estimation of a pseudo-orbit is given by From ( 7) and ( 8) it is clear that Theorem 4 [16] establishes that at least one of two pseudoorbits must have an error greater than or equal to the lower bound error.

The Proposed Method
The method proposed in this work can be summarized by the following steps: (1) Fit a model to the data series.Polynomial NARMAX models were the choice in this work.
(2) Choose two different interval extensions of the system under investigation following Definition 1.
(3) With exactly the same initial conditions simulate the two interval extensions.
Table 1: LLE for sine map investigated in [34].The reference value is calculated using the original equations of the systems, as proposed in [15].
[15] Our method LLE 1.133 1.114 (4) Use the least squares method to fit a line to the slope of the logarithm curve of the absolute value of the lower bound error defined in Section 2.3.The slope of the line is the LLE.
To the interested reader, a detailed description and explanation of why this simple method works may be found in [15,32].

Results
In this section, some examples are given to illustrate the usefulness of the proposed methods.

Example 1: Sine Map.
A unidimensional sine map is defined as The LLE from the original equation ( 10) was calculated using the method proposed in [15].To this end, two interval extensions given by the following equations were considered: The curve of the logarithm of the error bound is depicted in Figure 1(a) and slope of the linear regression, that is, the LLE, is shown to be 1.133.This value is in good agreement with the values found in the literature such as the one given in [33].
From the data generated by iterating the model in (10) the following polynomial NARMAX given in [34] was identified: Considering two equivalent interval extensions of the identified model (12) we have Figure 1(b) shows the logarithm of the lower bound error for the two interval extensions given by (13).Table 1 summarizes the results.Note that the similarity of the LLE can be also used to validate the model.13), with results for (  ) and (  ) and the same initial condition; that is,  0 = 0.1.LLE calculated is the inclination of the curve given by 1.114.

Mathematical Problems in Engineering
This system has been identified in [33] and the resulting polynomial NARMAX model is given by the following equation: for instance.Several other interval extensions could be determined but two of them would suffice for the purposes of finding the LLE.
As can be seen the interval extensions used to simulate the systems are mathematically equivalent, but they differ in the sequence of arithmetical operations.Figure 2 shows the logarithm of the lower bound error of the simulation using the two specified interval extensions.Table 2 summarizes the result of LLE for Rössler Equations.The literature suggests a value of 1.242 and the estimated value using the model procedure suggested in [33] is 1.566.The estimated value seems slightly different, but if the confidence interval to calculate the LLE using data is taken into consideration, one could consider 1.242 ± 0.399 instead and therefore the value of LLE calculated in [33] is in good agreement with the value in the literature.Similar reasoning can be applied to the value estimated using the method presented in this paper, that is, 1.530.Using a polynomial NARMAX and interval extension, our method has estimated the LLE to be equal to 1.530, while the authors in [33] calculated it as 1.566.

Example 3: Mackey-Glass
Equations.The Mackey-Glass is an interesting system used in many papers as an example of a chaotic and infinite dimensional system, since it is a timedelay system [36,37].The system equation is given by where  = 0.2,  = 0.1,  = 10, and  = 30.The authors in [20] identified the following model to this system using a sequence of 234 data points sample at     Figure 3: LLE calculation for Mackey-Glass Equations.The -axis is the simulation time and the -axis is the log of the lower bound error.Using a polynomial NARMAX and interval extension, our method has estimated the LLE to be equal to 0.006, the same value found in the literature using the dynamical equations.The interval extension used in this example is obtained just by replacing the term by which means that only the order of multiplication has been changed.Figure 3 presents the logarithm of the error for the Mackey-Glass Equation.Using the method proposed here the LLE is found to be exactly the same as the value provided by [20], as reported in Table 3.It is worth emphasizing that there is no need to calculate the Jacobian matrix of the model in the proposed method.

Example 4: A Chaotic Jerk Circuit.
A chaotic jerk circuit given in [21] was also used to further illustrate the proposed approach.The circuit is composed of resistors, capacitors, and operational amplifiers as shown in Figure 4.In this work, a polynomial NAR with nonlinearity degrees ℓ = 4 and   = 4 was used to identify the variable  of the chaotic jerk circuit, as shown in Figure 5.The data was produced by means of the simulation of the electronic circuit using LTspice.The obtained model has 67 terms, which have been omitted in this manuscript for the sake of simplicity.
Table 4: LLE for the variable  of the chaotic jerk circuit investigated in [21].The reference value was computed using circuit acquired data.Using a polynomial NARMAX and interval extension, our method has estimated the LLE to be equal to 0.0787, which represents a good agreement with the value suggested by [21].
For this case, the interval extensions can be found by changing the order in which the model terms were computed.The two interval extensions were generated by the simple order commutation of elements in the polynomial, as illustrated in Example 2. The interval extensions used to simulate the systems are mathematically equivalent, but are not computationally equivalent.Figure 5 shows the logarithm of the lower bound error of the simulation using the two specified interval extensions and Table 4 summarizes the results of the LLE for the chaotic jerk circuit.For comparison purposes [21] shows that, by using experimental data, the value of LLE can be estimated as 0.0735.On the other hand, the proposed methodology estimates the LLE value as 0.0787  6: LLE calculation for the tent map.The -axis is the number of iterates and the -axis is the log of the lower bound error.Using a polynomial NARMAX and interval extension, our method has estimated the LLE to be equal to 1.0461, while the correct value is 0.9993.which is in a good agreement with the value for the original system.

Example 5:
The Tent Map.To verify the robustness of the proposed approach under the presence of noise, a tent map (see (22)) was used as a system to be identified, within an additive zero mean Gaussian noise with standard deviation  = 0.02. () = 1 − 1.999 | ( − 1) − 0.5| .
A polynomial NAR was identified with ℓ = 4 and   = 2, yielding the model  () = 0.1133 ( − 1) + 14.0488 ( − 1) 4 + 0.5214 The interval extensions can be found by changing the order in which two terms were computed.In this example, the interval extensions were generated by the simple position change of the terms 0.1133( − 1) and +0.5214.Figure 6 presents the logarithm of the error for the identified model whereas Figure 7 shows the first return map of both identified and original tent map.
Despite the presence of noise affecting the system, the LLE found by the method proposed in this paper was in good agreement with the one analytically computed from the tent map, as shown in Table 5.

Conclusion
This paper presents a system identification approach to calculate the LLE based on the interval extensions.The method proposed here can be seen as an extension of the method developed in [15] in which the LLE is calculated by taking into consideration the properties of the exponential growing of the lower bound error [16].The method is designed specifically for the calculation of the LLE from chaotic data, which was not possible to undertake with method proposed in [15].A polynomial NARMAX model was used to identify the system, but other representations could also be used in the first step.Nevertheless, the polynomial NARMAX represents a very suitable mathematical representation to produce interval extensions, which can be obtained by simple manipulation upon the terms of the identified polynomial model.Compared to existing approaches, such as the one indicated in [20], the proposed method does not require the calculation of the Jacobian matrix and, more importantly, is very simple to use.We believe that the proposed method is also a good alternative for computing Lyapunov exponents from limited data.However, the computation of Lyapunovexponent spectrum, as performed in [38], is left for future investigations.
Simulation of(11), with results for (  ) and (  ) and the same initial condition; that is,  0 = 0.1.LLE calculated is the inclination of the curve given by 1.133.

Figure 1 :
Figure 1: LLE calculation of sine map: (a) original equations and (b) identified model.The -axis is the number of iterates and the -axis is the log of the lower bound error.

Figure 2 :
Figure2: LLE calculation for Rössler's Equations.The -axis is the simulation time and the -axis is the log of the lower bound error.Using a polynomial NARMAX and interval extension, our method has estimated the LLE to be equal to 1.530, while the authors in[33] calculated it as 1.566.

Figure 7 :
Figure 7: First return map for the tent map (∘, in black) and for the identified model (△, in red).

Table 5 :
LLE for tent map.