Detection of Chaos Using Reservoir Computing Approach

Detection of chaos in time series is of utmost importance in many scientific fields. Indeed, the presence of chaos and its significance, especially in multidimensional systems, plays an essential role in the control and analysis of such systems, and in their practical use in a variety of applications. In this paper, we demonstrate a new methodology for the detection of chaos in time series using a reservoir computing (RC) paradigm called conceptor-driven network (ConDN). Case studies on the known chaotic attractors (i.e. Lorenz, Rossler, Chua) of integer (conventional) and non-integer (fractional-order) orders, as well as a physically simulated and designed spintronic device (NCVO) are used in this study to validate the proposed chaos detection approach. The proposed chaos detection approach is tested on clean and noisy time series of the mentioned attractors. It outperforms the 0–1 chaos detection test and the largest Lyapunov exponent (LLE) estimation approach especially in the high noise-level conditions. In addition, the proposed approach is capable of differentiating the time series generated by the systems whose dynamics is at the edge of chaos. The simplicity of use of the proposed chaos detection approach can be counted, as well, as one of its main advantages over traditional chaos detection methods.


I. INTRODUCTION
Chaos is a natural or artificial non-linear behaviour found in systems. Such a behaviour is seemingly unpredictable without possibility to converge to a stationary or a periodic function of time, although it is actually governed by deterministic laws. The main reason behind such unpredictability is the high sensitivity a chaotic system shows to initial conditions. Indeed, an infinitesimal deviation of the system's starting point may lead to a massive change of its dynamics. Accordingly, a chaotic system can be characterized as [1]: bounded (dynamics stay inside an orbit rather than escaping off to infinity), deterministic (dynamics can be analytically described where same initial conditions give same behaviour over time), and sensitive (small perturbations to the system are exponentially amplified).
The detection of chaos in time series in the case of high-dimensional and noisy systems is not a trivial task. The common approach of detecting chaos in such time series is the computation of the largest Lyapunov exponent (LLE) The associate editor coordinating the review of this manuscript and approving it for publication was Michael Lyu. whose sign informs about the chaotic (or periodic) behaviour. Although instability mechanisms of chaotic systems are rigorously described with the Lyapunov exponent (and vectors), its computation remains very challenging even today despite the numerous algorithms and methods proposed in the literature. One way to compute the LLE is to use the so-called direct methods [2]. These methods, by tracking the evolution of the distance between neighboring orbits, afford a graphical feedback whether local exponential divergence due to different reasons (i.e. measurement noise, length of timeseries) is properly identified or not. On the other hand, the full spectrum of Lyapunov exponents can be extracted by the dot products of perturbations vectors and their derivatives. The n Lyapunov exponents are estimated by the integration of n − 1 perturbations and by using basic computing operations such as summation, subtraction, multiplication, and division. This yields to faster and more efficient results compared to the methods involving the calculation of perturbations lengths logarithms and Gram-Schmidt orthonormalization, especially in low-dimensional problems [3]. Machine learning has also been investigated in LLE computation. An approach is demonstrated in [4], where a single layer feed-forward network of given parameters is trained on the time-series to be inspected. Then, the parameters of the initial network are varied and the network is retrained. The output signals produced by the two networks are analyzed and compared for the extraction of the entropy characteristics of the input time-series. A different methodology, based on the use of reservoir computing, is shown in [5]. A reservoir of a recurrent neural network is first built and trained on a given input time-series. In the testing phase, the reservoir becomes capable of producing an image of the input time series whose behaviour can be defined using the parameters of the neural network. Thus, the Lyapunov exponent estimation can be calculated using classical ways. The main advantages of the latter method is that it is applicable even in model-free cases, where data could be of unknown source or too noisy to be exploited; besides it is sufficient to be fed by short time-series where large ones are generated by the network on which the computation of Lyapunov exponents is done. Analytical modelling is an additional option for the computation of Lyapunov exponents [6]. This can be done by identifying a model for the available data using NARMAX (Non-linear Autoregressive Models) representation, and then extract the interval extensions that will be used to calculate the lower bound error. Eventually, the logarithm of the lower bound error is fitted by a simple linear function from which the LLE is extracted. Nevertheless, the poor reliability of these different methods when applied to noisy time series often limits their use to clean time series from theoretical systems.
Another approach for chaos detection in time series is the 0-1 test proposed by Gottwald and Melbourne in [7]. The advantage of this method over the largest Lyapunov exponent is the possibility to detect chaos even when the knowledge about state-space variables of the observed system is incomplete or unknown, which is the case of many real systems.
The 0-1 test shows its limits in noisy time series, especially when the gauge-measurements needed for a reliable detection and corresponding to the regular system behaviour, cannot be obtained (often in real world systems). Recently, Tempelman and Khasawneh presented an enhanced version of the Gottwald and Melbourne's 0-1 test where by using topological data analysis its reliability for high-noise levels time series is greatly improved [8]. Another class of chaos detection approaches is the use of non-linear prediction. This method is originally established for short-term generation of non linear time-series. However, by applying a long-term nonlinear prediction technique on a given signal and calculating the accuracy (error between real and predicted signal), the differentiation between chaotic and regular dynamics can be attained. This notion of chaos detection can be achieved using a variety of approaches found in the literature such as radial basis function [9], neural networks [10], [11]. In [11], the authors propose three ANN-based designs (MLP, NARX, and ESN) with their FPGA implementation for the chaotic time series prediction. The authors presented optimized hardware pipeline architectures of the considered models allowing high speed computation and optimal chaotic time series prediction.
The frequency analysis as well, and particularly Fourrier transform (FT), could serve for the task of chaos detection in time-series as presented recently by Tlelo-Cuautle et al. in [12]. The authors used this FT-based chaos detection to accelerate the optimization process of fractional-order chaotic networks driven by metaheuristics such as differential evolution and particle swarm optimization. They observed that chaotic time series are usually related to high amplitude FT presentations whereas regular time series are of low amplitude presentations. Indeed, chaos detection in time series is still a widely open research problem, because chaos can be induced by noise as well as some standard random processes may have deterministic origin [13].
In this paper, we propose an alternative reservoircomputing based approach for chaos detection. By introducing a time series under investigation as input without any additional parameters, the insights about the chaotic or regular behaviour of the time series are produced. Case studies on the known chaotic attractors (Lorenz, Rossler, Chua, Chen and Lu), obtained using integer (conventional) and noninteger (fractional) orders of derivations, as well as a real world spintronic device, are used in this study to validate its effectiveness. The simplicity of use of our proposed method combined with high robustness to noise are its main advantages. The preliminary results of this work has been presented at ISCAS 2021 in [14]. In this paper, a more in-depth analysis of the presented approach as well as an extended experimental section including real emergent chaotic devices has been provided.
The rest of the paper is organized as follows: Section II gives background details about the computation of the Largest Lyapunov Exponent (LLE), 0-1 test, and conceptordriven networks. The proposed chaos detection approach is presented in Section III. In Sections IV and V we demonstrate and discuss respectively the effectiveness and robustness of our proposed method on several theoretical and real worlds systems including fractional order systems. Finally, Section VI gives concluding remarks and presents some future perspectives.

A. LARGEST LYAPUNOV EXPONENT
One way to detect chaos in a time series (or in a system at its origin) is to calculate its Largest Lyapunov Exponent (LLE). The LLE measures the exponential rate of divergence of adjoining trajectories on the system's attractor. Generally, the calculation of the LLE passes through two main stages: (1) constructing an appropriate state space that captures the system's dynamics, and (2) measuring the strength of divergence of the points on the state space. For that, two random neighboring points from the state space are taken as initial points, and their separating distance is continuously measured. The common challenge that all LLE's estimation algorithms face is the influence of noise, especially when dealing with experimental data where the possibility of choosing false VOLUME 10, 2022 neighbors increases. One way to tackle this problem is to consider multiple neighboring points and average their distances from the reference ones. In addition, other factors making the computation of the LLE very challenging and which may lead to its false estimation are inherent noise, system uncertainties and perturbations as well as the inappropriate choice of parameters [15].
The Rosenstein's approach is one of the algorithms for the estimation of the LLE from small data sets [16]. It starts with the placement of M groups of neighboring points whose initial separation d j−ave (0) is taken as the average of Euclidean distances of all points to the reference one (j ∈ [1, M ]): where E denotes the mean, and d jn is the distance separating the n th nearest neighbor X jn from the reference point XĴ : where ||.|| is the Euclidean norm. Note that all nearest points have to be located from the reference point at a distance at least equal to the mean period of the time-series under inspection (i.e a nearest point should be located on a trajectory different from that of the reference point). The evolution of the average distance between neighboring points is commonly defined using the Largest Lyapunov Exponent λ max as follows: where d ave (0) is the initial separation of the neighboring points. Using Equation 3, one can write for the j th pair: Applying a logarithmic function on both sides, we get: The above equation describes M parallel lines whose slopes are relatively proportional to λ max . Finally, λ max can be accurately estimated by applying the least-squares fit on the averaged line defined as follows: where E[d j−ave (i)] is the average over all pairs of j.
The 0-1 test is a chaos detection method firstly introduced by Gottwald and Melbourne for the analysis of deterministic dynamical systems [7]. Unlike the LLE methods, the 0-1 test uses time series to distinguish between regular and chaotic dynamics with no need to reconstruct the phase space. The final result of this test is either 1 or 0 indicating the presence of chaos or its absence respectively. A recipe of the 0-1 test can be summarized in computing the following quantities: Translation Variables p c (n) and q c (n), Mean Square Displacement D c (n), and Asymptotic Growth Rate K c (n).
Given a discrete time series φ(j) for j ∈ [1, N ], the translation variables p c (n) and q c (n) are computed as follows: for n ∈ [1, N ] and c ∈ (0, π). The number of c values taken in the test, is denoted by N c . In practice, N c = 100 seems to be sufficient [7]. The translation variables are bounded if the underlying dynamics is regular (periodic or quasiperiodic) whereas they behave asymptotically like Brownian motion for mutlidimensional chaotic systems.
The mean square displacement D c (n), a second quantity needed for the 0-1 test, describes the diffuseness of a given behaviour. For regular dynamic systems, D c (n) is a bounded function of time, whereas it scales linearly with time for chaotic systems. Thus, the growth rate of D c (n) is the main criterion used for chaos detection. D c (n) is defined as follows: for n ≤ n cut where n cut ≤ N . In practice, n cut = N /10 is a common choice [7]. The term V osc (c, n) is added for less non-linearity during convergence: Finally, the asymptotic growth rate K c , which describes the strength of the correlation of D c (n) with the linear growth, is estimated. This is done by the linear regression of the log − log plot of D c (n). Accordingly, one can write: whereD , is set to eliminate the negative values. Graphically, the slope of the straight line fitting the graph of log D c (n) versus log n is an approximation of K c . This process is repeated N c times, where each time c takes a new random value in the open set (0, π). Note that this set could be reduced for the sake of test validity. Consequently, the median of all K c values is computed to obtain K : K = 0 signifies regular dynamics, whereas K = 1 indicates the presence of chaos. It should be noted that the dataset N should be of a sufficient data length so that an asymptotic behaviour of D c (n) takes place. Although its few requirements in application, the 0-1 Test is however not always reliable. The used data should not be oversampled nor downsampled, otherwise false results can be obtained. Furthermore, the parameter c has to be tuned in a convenient range to avoid resonance that could yield to K c ∈ [0, 1].

C. CONCEPTOR-DRIVEN NETWORK
Conceptor-Driven Network (ConDN) is a recent Reservoir Computing approach, firstly introduced by Jaeger in 2014 [17], [18]. It is based on the principle of long short-term memories (LTMs) [19], [20] initially established for the prediction of continuous time-series. Its general structure, as shown in Figure 1, is an enhanced version of the echo-state network (ESN) [11], also established for prediction task. A ConDN consists of four layers: a K -dimensional input layer, an N × N reservoir, a conceptor layer of dimensions N × N , and a single-unit output layer. The conceptor layer, found between the reservoir and the output layer as shown in Fig. 1, presents the main difference in such networks with respect to the conventional ESNs.
K time-varying input signals, each of length L, are introduced to the network via the input layer. An input time series, also called a pattern, is indexed by j ∈ [1, K ] and denoted by p j (n). The reservoir placed between the two layers, is a Recurrent Neural Network (RNN) at which, the connections between nodes are spontaneous and in all directions. The conceptor layer involves a group of matrices called conceptors C j (n), where each represents, if trained, a given input pattern p j (n).
In the installation phase, an N -neurons RNN is established by randomly creating a dynamic reservoir described by a tuple (W in , W * , b) where W in is the input N × K weights matrix; W * is the internal N × N connection matrix; and b is a bias vector of length N . Afterwards, the readout weights matrix W out is computed by introducing into the network, as an input, a random time-varying signal. The main reason behind this choice is to find an optimal output matrix that could perform well for any input patterns. For that, an L-dimensional white noise signal v(n) enters the network to run the system according to the following state-update equation: x The resultant state vector x v is used then to find the optimal estimation of W out by minimizing the following quadratic loss by linear regression: The second phase represents the loading phase. In this phase, the ConDN starts to recognize the input patterns p j (n). K time varying signals are separately inserted into the network to compute their corresponding intermediate state vectors x j (n) using the following state-update equation: The collected states from Equation 13 hold the dynamical characteristics of the original input patterns. Thus, it becomes possible to expel the input drivers p j (n) and build instead one matrix that involves all the states x j (n). For that, the input weights matrix W * is replaced by the input internalization weights matrix W by minimizing the quadratic loss below after the washout period n 0 ends: The intermediate states generated from Equation 13 are used as well for the computation of the conceptor matrices C j . These matrices guide the system in the testing phase rather than the input signals p j (n): In mathematical terms, C j works as an identity matrix for the term f (Wx j (n)+b), whereas it is a null matrix for the terms holding other states. Accordingly, the optimal values of the conceptor matrices C j can be computed using the following quadratic loss: where E is the mathematical mean, . 2 fro is the squared Frobenius matrix norm, and γ j is the aperture parameter that specifies the internal structure of C j . A conceptor matrix C j is close to the identity matrix I if γ j is large. Conversely, when the aperture is small, the C j is close to the null matrix O. This mutation can be evidently observed in Equation 16. When C j ≈ I , the first term converges to its minimal value. By contrast, the second terms vanishes in the case C j ≈ O. The minimization of the quadratic loss above leads to the following optimal solution of C j : being the correlation matrix whose state vector x j is already collected for the computation of the internal matrix W (Eqs. 13 and 14).
Finally, in the testing phase, any pattern p j (n) can be re-estimated by the ConDN by running first the state-update equation (Eq. 15) followed by the linear output equation below:

III. PROPOSED CHAOS DETECTION APPROACH
The proposed ConDN-based chaos-detection method is summarized in Algorithm 1. The inputs of this method are two vectors: p -representing a time series whose dynamics should be analyzed (regular or chaotic); and s -representing an ordered sampling vector whose elements s i , i ∈ [1, |s|] provide sampling factors to use for the sampling of the initial time series p. Note that, |s| denotes the length of s. The first phase of this approach is to build the database of input signals needed to feed the ConDN network. This database is built by using the time series p and sampling factors s i (s i ∈ s, i ∈ [1, |s|]). The first element of this ordered sampling vector is always 1 (s 1 = 1) meaning that the initial time series without sampling is always a part of the database of input signals. The size of the sampling vector s should be at least 2, i.e |s| ≥ 2. In the case |s| = 2, the second element of the sampling factor is s 2 = 2, meaning that the initial time series is sampled by a factor of 2. For |s| > 2, all other values of sampling factors are allowed and should be in increasing order in the vector s. Let us assume the size |s| = N S > 2 for the following explanation. The result of this first phase are N S signals q i , i ∈ s of different size A ConDN network must be fed with the signals of the same size. Thus, the size equalizing of the signals q i , i ∈ s from the first phase is mandatory. The size of the smallest signal q i , i ∈ s corresponding to the highest sampling factor s |s| = N p must be used for this operation. As the result of the equalizing operation, we have N S signals p i , i ∈ s of the same length, which can be presented in a matrix form P = [p 1 p 2 . . . p N p ] as shown in Algorithm 1, line 11. In the next phase, the prepared input signals matrix P is fed to an already installed ConDN network  (see Section II-C for more details). The training phase of the ConDN networks is carried out on these input signals.
As a result of this training phase, the conceptor matrices The size of the conceptor matrix is equal to the size of the ConDN's reservoir (see Section II-C). Afterwards, each generated conceptor matrix C j undergoes a singular value decomposition (SVD) as shown in Algorithm 1, line 16: where S i is the diagonal matrix of SVD at which the eigenvalues of C i are arranged on its diagonal in descending order. The matrices S i , i ∈ s contain the precious information about the dynamics of the input signals p i , i ∈ s and are used to determine their either regular or chaotic behaviour. For each matrix S i , i ∈ s, we compute then the value of trace t i = tr(S i ), i ∈ s corresponding to the summation of all eigenvalues on its diagonal, as illustrated in line 19 of Algorithm 1. Large values of t i , i ∈ s are expected to arise in the case of chaotic patterns, whereas small values indicate the absence or the low level of chaos. Moreover, the values of t i , i ∈ s are expected to increase with higher sampling factors s i , i ∈ s in the case of chaotic patterns, and to remain stable and low otherwise. In the last phase, by using the computed trace values from the previous phase, the dynamics variation vector can be computed. The elements of are computed from the trace vector T = [t 1 t 2 . . . t N P ] of size N S by subtracting the value of t s 1 = t 1 from all other elements starting from its second element t s 2 . Thus, the size of the vector is The results obtained in the dynamics variation vector can be analyzed in two ways: a quick estimation of the nature (regular or chaotic) of the input time series called one-shot sampling method; or a much deeper analysis based on the fitting results of the relationship (s) called multiple-shots sampling method. In the first method, the sampling vector s has only two values s = [1,2], meaning that the initial time series p is sampled only once with a factor of 2. The two signals or patterns (the initial time series p 1 and sampled by 2 p 2 ) once equalized in length are introduced to the ConDN, and allow to obtain, after the ConDN training, two conceptor matrices C 1 and C 2 corresponding to each input pattern respectively. By computing the SVD matrices S 1 and S 2 of these conceptors as well as their trace values t 1 and t 2 , the dynamics variation vector can be obtained. has only one element 1 corresponding to the difference t 2 − t 1 , as shown in line 23 of Algorithm 1. If p (or p 1 and p 2 ) is chaotic, two distinct values of t 1 and t 2 are expected, and thus a large value of 1 . On the other side, sampling should not largely affect the value of trace in the non-chaotic case, yielding to a small value of 1 .
In the multiple-shots sampling method, the sampling vector s has more than two values s = [1, 2, · · · , N p ] and is used for a much deeper analysis of the evolution of the obtained dynamics variation vector with sampling s. For N S = |s| different sampling factors in s, the proposed approach will generate a dynamics variation vector of size N S − 1, as shown in line 29 of Algorithm 1. The evolution of with sampling s should follow a capacitor charging function trend (see Equation 20) in both, regular and chaotic behaviour case, with much higher values of fitting parameters a and τ in the latter case due to the rapid increase in the value of the trace in correspondence with the increase of the sampling factor.

IV. CHAOS DETECTION RESULTS
In this section, the proposed chaos detection approach has been tested and validated on known chaotic systems described with differential equations of integer and fractional orders that can be tuned to generate regular dynamics also, as well as on an emergent non-linear spintronic device called Nano-Contact Vortex Oscillator (NCVO), which is still not fully described with analytical expressions. For some of the case studies, the proposed approach has been compared with the 0-1 Test approach detailed in Section II-B.

A. CASE STUDIES 1) KNOWN CHAOTIC ATTRACTORS
Five known chaotic systems of integer orders were considered in this study: Lorenz, Chua, Rossler, Chen, and Lu [21]- [24]. We also include the two fractional-order (FO) systems: FO Rossler and FO Chen [25]- [27]. A 3D representation of their attractors in the chaotic state is shown in Figure 3. They describe continuous-time dynamics of these systems and are defined with non-linear ordinary differential equations as follows: where W (w) = f (δ, ζ, w(t)) • Rossler System:ẋ where q ∈ [0, 1] is the derivation order.  • FO Chen System: where q ∈ [0, 1] is the derivation order. The parameters and orders shown in the ordinary differential equations above affect the shape of an attractor and it may even vanish at some values. For that, we consider for each integer-order system two sets of parameters: the first one reinforces its chaotic dynamics whereas the second one its regular (or non-chaotic) dynamics. Some parameters are kept unchanged among the two sets. For the fractional-order systems, we fix the parameters and vary the order to realize chaotic and regular dynamics. The values of the parameters as well as the orders used in our simulation will be shown later.

2) NANO-CONTACT VORTEX OSCILLATOR
Nano-Contact Vortex Oscillator (NCVO) is a highly non-linear spintronic device of spin valve structure [14], [28] (See Fig. 4). The most significant layers at which the physical phenomenon occurs are: (1) the Permalloy layer at which the magnetization alters continuously and it is called the free layer; and (2) the Cobalt layer that shows uniform magnetization so it is called the fixed layer. The two ferromagnetic layers are separated by a copper band. A metallic portion on top, called the Nano Contact (NC), allows the passage of the current into the device.
When the NCVO is subjected to a direct bias current I and a magnetic field B, a magnetic vortex is formed in the upper part of the body. This vortex revolves around the center of the device creating an alternative resistance. With the applied DC current, an oscillating voltage V GMR at the vortex gyration frequency appears across the device. Therefore, the NCVO exhibits different behaviours that can be classified in two major modes: • Modulated non-chaotic mode: the vortex is exposed to a periodic relaxation in its gyration. After fulfilling a given number of complete cycles, the vortex directs toward the center and changes its polarity and the sense of gyration.
• Chaotic mode: the vortex revolves around the nano-contact with no definite rules. Unlike the former mode, the gyration relaxation occurs in a non-periodic way although the vortex shows complete cycles.
In addition, the NCVO device exhibits also a mode called non-modulated non-chaotic mode at which the vortex rotates continuously around the nano-contact only in one direction. The different modes can be studied by visualizing the time varying signals of both the NCVO's space averaged output magnetization of the free layer M (t) and the vortex position with respect to the NanoContact s(t). These output signals are composite which makes their evaluation quite difficult. In the non-modulated mode, these time series are periodic signals whose fundamental frequency f 0 specifies the gyration speed of the vortex. In the modulation mode, the time series generated by the NCVO are quasi-periodic. This mode is also called the locking mode, where f 0 is locked to the modulation frequency f mod . The later frequency specifies the number of cycles the vortex performs before changing its polarity. In this case, the fraction f mod /f 0 is simply a rational number. In the chaotic mode, the modulation frequency is no more locked to the fundamental frequency; this leads to an irrational value of the fraction f mod /f 0 . Two types of data are generated from an NCVO: micro-magnetic simulated data and real experimental data. In simulation, M (t) and s(t) are generated using the finite difference micromagnetic solver Mumax3 [29]. This solver implements a numerical space and a time integration of the Landau-Lifshitz-Gilbert equation including the different spin transfer torque terms. More details concerning the procedure description and the values of parameters used in simulation can be found in [28], [30], [31]. The extracted signals M x (t), M y (t), and M z (t) represent the projection of M (t) on the three dimensions in space respectively. Similarly, s(t) has the components: x(t), y(t), and core polarity z(t) = ±1.  The gyration of the magnetic vortex around the NanoContact in the upper layer of the NCVO reflects its operating mode whether it is non-modulated, modulated, or chaotic (see Fig. 5). However, this circular movement representing the polarity reversal of the vortex takes places at slightly distinct positions in its trajectory. As a result, the time varying signals corresponding to the NCVO's output magnetization (M x (t), M y (t), and M y (t)) as well as the vortex position in space (x(t), y(t), and z(t)) are contaminated with a jitter effect, an additional slightly chaotic behaviour different to the chaos obtained due to the polarity reversal frequency.
The NCVO's output signals are composed of different patterns. Each pattern represents one cycle of the trajectory of the vortex gyration around the Nano-Contact. Consequently, the shape of a pattern depends on the direction of gyration of the vortex around the Nano-Contact. Due to the jitter effect, the patterns are slightly variable in length. Thus, it seems difficult to determine whether the observed complex signal arises from a chaotic or stochastic process. One way to consider only the chaotic behaviour of the NCVO device and get rid of the jitter effect is to build an artificial signal having the same patterns (approximated with sine and cosine functions) as the ones found in the NCVO's output signals which evolve in the rhythm of the vortex polarity and direction of gyration, but have the same length. This artificial signal is also used as one of the case studies in this evaluation part (see Figure 6).

B. ONE-SHOT SAMPLING
We first evaluated the one-shot sampling method described in Algorithm 1 on the seven known attractors detailed in Section IV-A1. The set of parameters of the integer-order TABLE 1. One-shot sampling: values of the parameters used in the ordinary differential equations of the five studied integer-order systems (Equations 21 to 25) in the chaotic and regular modes.

TABLE 2.
One-shot sampling: values of the parameters and derivation order used in the ordinary differential equations of the two studied fractional-order systems (Equations 26 and 27) in the chaotic and regular modes.
systems, as well as the set of parameters and orders of the fractional-order systems, used to generate signals with chaotic and regular nature, are presented in Tables 1 and 2 respectively. The first component x(t) of each of the seven systems, for a given set of parameters and orders corresponding to chaotic or regular nature, is injected beside its sampled counterpart into a separate ConDN. The reservoir size N of all networks is arbitrary set to 1000. Figure 7 presents the variation of the trace values t i , i ∈ s = [1, 2] of the chaotic and regular input patterns as a function of their data length, where i = 1, 2 denote the raw p 1 and sampled p 2 patterns respectively. The trace values of the two chaotic (initial and sampled) patterns are plotted in blue, whereas the ones of the non-chaotic (initial and sampled) patterns are in red. Generally speaking, t i of a chaotic signal is large compared to that of a non-chaotic one. This could be referred to the high instability that characterizes chaos and is at the origin of large eigenvalues captured by the ConDN's conceptor matrix. Additionally, sampling of these input patterns mainly impacts chaotic time-series. This can be noticed for all cases in Figures 7 where the two blue plots corresponding to the two (initial and sampled by 2) chaotic signals are separated by a big gap, whereas the two non-chaotic signals in red are almost stacked. To numerically express this difference for the two behaviours of each system, the dynamics vector variation = [ 1 ] is calculated. These results are shown in Figure 8. The non-chaotic signals of any of the studied conventional systems show a small value of 1 (close to 0), whereas in the chaotic case 1 is large (between 45 and 95).

C. MULTIPLE-SHOTS SAMPLING
We also evaluated the multiple-shots sampling method described in Algorithm 1. For this method, three of the five integer-order known attractors (Rossler, Chen, and Lu) with the two fractional-order systems (FO Rossler and FO Chen), all detailed in Section IV-A1, as well as the (simulated and artificial) signals describing the NCVO's vortex position in space (see Section IV-A2 for more details), were used.
For each system, to the two modes (chaotic and nonchaotic) evaluated in the last section, a half-chaotic one was VOLUME 10, 2022 added. The half-chaotic mode represents a transient state between the chaotic and non-chaotic modes. It should be mentioned that the two remaining integer-order conventional attractors, Lorenz and Chua, were discarded in this part due to the difficulties to obtain a half-chaotic behaviour for these systems (abrupt transitions from regular to chaotic cases). For the integer-order systems, the three different modes are obtained by varying the parameters of the differential equations 23 to 25. The values of these parameters used in the simulation are shown in Table 3. For the fractional-order systems, the three different modes are obtained by fixing the parameters and varying the fractional order of the differential equations 26 and 27. The values of such parameters and orders used in the simulation are shown in Table 4. On the other hand, for the NCVO device, a micromagnetic simulation is carried out at I = 17mA with different values of the applied magnetic field. The horizontal magnetic field is fixed to B x = 1mT , whereas the vertical magnetic field takes the values B z = 10mT (for non-chaotic case), B z = 1mT (for half-chaotic case), and B z = 20mT (for chaotic case). Similarly, three behavioural modes of x(t) are obtained, where x(t) denotes the horizontal projection of the time-varying position of the NCVO's vortex in space. Note that, the non-chaotic signal of the NCVO device used in this evaluation is in the modulated mode. The number of neurons used for each system is shown in Table 5. It should be noticed   that the difference in number of neurons used for different case studies is mainly related to the optimization technique we applied to accelerate simulations. Indeed, for each system and mode, this technique computes (based on the experimental trials) the number of non-zero eigenvalues of the ConDN's conceptor matrices and based on this reduces the number of neurons to use. Consequently, significant acceleration of the overall simulation time without influencing the final result is obtained. Figures 9, 10 and 11 present the variation of the dynamics vector (s) with its fitting function (see Equation 20) for the five known chaotic systems and the NCVO respectively. From Figures 9 and 10, it can be noticed that the signals of FO Rossler, FO Chen, Rossler, Chen, and Lu systems show a quick increase in for the chaotic case, almost constant for the non-chaotic case, and an intermediate growth for the half-chaotic case. The same can be noticed for the NCVO's results from Figure 11. Indeed, the three modes of the NCVO can be easily distinguished by the dynamics vector variation (s). The non-chaotic and half-chaotic modes are close compared to the chaotic one, which is mainly due to the few pattern distortions that differentiates the half-chaotic signal from the regular one. The effect of filtrating the jitter effect (see Section IV-A2) by imposing the same length to all patterns can be noticed by comparing Figures 11(a) and 11(b). Each plot, of the three modes, shows now a slower growth as function of the sampling factor (smaller slope). Moreover, the half-chaotic mode shows the same evolution of (s) as that of the non-chaotic mode after expelling the jitter effect as these two modes are distinguished by only few pattern distortions.
Using the linear regression, the fitting parameters a and τ presented in Equation 20 of all studied cases are computed. These values of a are shown in Figure 12. The values of a can be categorized in three ranges: less than 20 for a non-chaotic signal; greater than 50 for a chaotic signal; and in between for half-chaotic signals. It is clear that the presented chaos detection method shows clear distinction between these three modes. Moreover, the same case studies were introduced to the 0-1 Test for comparison. These results are shown in Figure 13. As it can be noticed, the 0-1 Test does not provide the same efficiency in chaos detection where the three behavioural modes can not be defined using the value of K median . Indeed, a half-chaotic signal is sometimes recognized as a chaotic (K median ≈ 1) whereas at other times it is recognized as a non-chaotic signal (K median ≈ 0).

V. ROBUSTNESS TO NOISE
In this section, the robustness to noise of the two proposed chaos detection methods, compared to that of the 0-1 Test and largest Lyapunov exponent approaches (detailed in Sections II-B and II-A respectively), has been evaluated. First, the robustness of one-shot sampling method was examined using additive noise. Each studied time-series, before its examination was contaminated by an additive white noise of a Signal-to-Noise Ratio (SNR) between 184dB and 11dB. The dynamics vector variation (s) for sampling factor 2 versus the SNR of the added noise for the seven conventional systems is shown in Figure 14. As noticed, the value of in the chaotic case typically decreases with the addition of noise. In contrast, the value of in the non-chaotic case varies but stay negative for any SNR value. This reveals the fact that our proposed one-shot sampling method mainly recognizes deterministic dynamics. The FO Rossler, FO Chen, Lorenz, Rossler, and Lu systems (Figures 14 (a), (b), (c), (e) and (g) respectively), can tolerate a white noise of SNR up to 11dB, where at any noise level, is always positive in the chaotic case and negative in the non-chaotic case. The Chua system can handle a SNR < 14dB (Figure 14(d)), and the Chen system can handle a SNR < 11dB (Figure 14(f)). The same seven systems were evaluated by the 0-1 Test approach. The results are shown in Figure 15. As noticed, the value of K median VOLUME 10, 2022   for any of the inspected chaotic systems remains constant at 1 for any added white noise level. In contrast, the K median of a non-chaotic system increases with the addition of noise. It seems evident that the 0-1 Test, unlike the one-shot sampling method, mixes up between deterministic and stochastic behaviours. Furthermore, the 0-1 Test does not show a better robustness against white noise, where the tolerance of the  inspected systems ranges between SNR = 30dB and 14dB. Note that, we assume the 0-1 Test is still robust if K median > 0.5 in the chaotic case and K median < 0.5 in the regular case.
The robustness of the multiple-shots sampling method was studied using highly noisy experimental traces of NCVO. For the sake of comparison, the same data was tested using one-shot sampling method as well as the 0-1 Test and largest Lyapunov exponent approaches. Wherefore, a real NCVO was experimentally run using a DC bias current I = 17mA and subjected to different magnetic field quantities. Consequently, five packages containing recorded horizontal time varying position x(t) were analyzed: A, B, C, D, and N. In the packages A and C, a synchronized gyration of the vortex around the nano-contact was established during recording resulting in that x(t) behaves as a quasi-periodic signal (regular nature). In the package B, a chaotic movement of the vortex is detected, which normally gives rise for a chaotic signal x(t). In the package N, the vortex is either hidden or static. In such case, x(t) is described as pure white noise signal. In the package D, the vortex starts gyration chaotically, and then it stabilizes with a synchronized rotation.
A signal x(t) from any package does not show clean shape due to the high contamination of the recorded signals with the experimental noise. From each package (except for the package D) two samples of x(t) of the same nature (regular or chaotic) were considered. The two samples from the package D were recorded at different behavioural modes. For all tests, a ConDN with a 100-neurons reservoir was considered.
Starting with the multiple-shots sampling, the different plots of (s), as well as their fitting plots are all shown in Figure 16(a). As it can be noticed, the three plots corresponding to the three chaotic samples (the two samples of the package B and the first one from the package D) show a very quick growth compared to the plots of the other non-chaotic samples. The values of the slope a of the fitting function for all plots are shown in Figure 16(b). As noticed, the dynamics vector variation of the three chaotic samples can be fitted to a charging capacitor function having a slope a ≥ 60. On the contrary, the non-chaotic samples are fitted with a slope a ≤ 20. These results are in agreement with the results obtained in Section IV-C, where a chaotic signal was characterized by a slope a ≥ 50, a regular signal by a slope a ≤ 20, and any signal with a slope 20 < a < 50 was related to a transitional mode between chaos and regular nature.
We also evaluated the one-sampling method using the same NCVO experimental data. Figure 17 shows the values of the dynamic variation vector of the two samples of each package after applying one-shot sampling method. As shown, of the three chaotic samples (the two samples of package B and the first sample of package D) is between 17 and 20. The remaining non-chaotic samples have lower percentage variation that ranges between 2 and 14. Although the value of is larger in the chaotic case than in the non-chaotic case, but the values still close to each other. Thus, the percentage variation derived by one-shot sampling does not show clear differentiation between chaotic and non-chaotic experimental samples. Moreover, these findings are incompatible with conclusion drawn in Section IV-B, where a chaotic system is

characterized by
> 45 whereas a non-chaotic system is characterized by ≈ 0.
In addition, the same samples were tested using the 0-1 Test and the largest Lyapunov exponent (LLE) approach. For the 0-1 Test, the results are shown in Figure 18. As it can be noticed, the 0-1 Test is not capable to differentiate between chaotic and regular experimental data, where K median ranges between 0.5 and 1 regardless of the tested sample. For the LLE, the same experimental samples were evaluated by T. Devolder et al. in [32] using two LLE methods: Wolf [33] and Ronsentein (detailed in Section II-A). The results are presented in Figure 19. As shown, the packages are distributed irregularly along the axis of the applied magnetic field µ 0 H (x-axis). Each of the packages C, D, and N are partitioned into two parts. As a reminder, the samples of the first part of VOLUME 10, 2022  package D is non-chaotic whereas those in the second part are chaotic. The value of the LLE (λ) is large for the samples of package B and the second part of package D, due to their chaotic nature. However the value of λ is always positive even for the non-chaotic samples due to the jitter effect. Moreover, the value of λ is high for the samples of package N compared to other non-chaotic samples even that in package N the samples are simply white noise signals (the vortex is either hidden or static). Consequently, the LLE method is not reliable in distinguishing between the chaotic and regular experimental data. To summarize, the multiple-shots sampling ConDNbased method surpasses the 0-1 Test, the LLE, and the singleshot ConDN-based methods in terms of robustness. The later methods do not afford a decisive differentiation between the chaotic and non-chaotic modes of a real NCVO device.

VI. CONCLUSION
In this paper, we investigated a novel methodology for chaos detection in time series using a reservoir computing (RC) paradigm called concept-driven network (ConDN). This chaos detection approach, that involves the sampling of the input series before being injected to the network, can be applied in two ways: (1) one-shot sampling, and (2) multipleshots sampling. The proposed approach has been validated on known chaos attractors driven by ordinary differential equations (i. e. Lorenz, Rossler, Chen) of both integer and fractional order, as well as on a highly non-linear spintronic oscillator called NCVO whose data are collected by micromagnetic simulation. All considered case studies involve chaotic and regular behaviors, whereas some of them involve additionally a transient behavior between them. The robustness of the proposed ConDN-based chaos detection method is also evaluated, using data contaminated by additive noise as well as highly noisy experimental data (for the NCVO).
The multiple-shots sampling method has shown to be outperforming the 0-1 Test, the largest Lyapunov exponent method, and the one-shot sampling methods in terms of efficiency and robustness to noise.