# Detection of Chaos Using Reservoir Computing Approach

Ali Rida, Slavisa Jovanovic, S. Petit-Watelot, Hassan Rabah

► To cite this version:

Ali Rida, Slavisa Jovanovic, S. Petit-Watelot, Hassan Rabah. Detection of Chaos Using Reservoir Computing Approach. IEEE Access, 2022, 10, pp.52686-52699. 10.1109/ACCESS.2022.3174960. hal-03981818

HAL Id: hal-03981818

https://hal.science/hal-03981818

Submitted on 10 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Received April 11, 2022, accepted May 2, 2022, date of publication May 13, 2022, date of current version May 20, 2022.

Digital Object Identifier 10.1109/ACCESS.2022.3174960

# Detection of Chaos Using Reservoir Computing Approach

ALI RIDA ISMAIL<sup>10</sup>, SLAVISA JOVANOVIC<sup>10</sup>, (Member, IEEE), SÉBASTIEN PETIT-WATELOT, AND HASSAN RABAH<sup>10</sup>, (Senior Member, IEEE)

Jean Lamour Institute (JL), University of Lorraine, UMR 7198, 54000 Nancy, France

Corresponding author: Slavisa Jovanovic (slavisa.jovanovic@univ-lorraine.fr)

This work was supported by the French National Research Agency through CHIPMuNCS under Contract ANR-17-CE24-0008.

**ABSTRACT** 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.

**INDEX TERMS** Chaos detection, time series, reservoir computing, conceptors.

## 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 time-series) 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 non-linear 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 Fourier 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 reservoir-computing 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 non-integer (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 conceptor-driven 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.

## II. BACKGROUND

### 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
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 Rosenberg'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*<sub>j-ave</sub>(0) is taken as the average of Euclidean distances of all points to the reference one (*j* ∈ [1, *M*]):

$$d_{j-ave}(0) = E([d_{j1}(0), d_{j2}(0), \dots, d_{jn}(0)]) \quad (1)$$

where *E* denotes the mean, and *d*<sub>*jn*</sub> is the distance separating the *n*<sup>th</sup> nearest neighbor *X*<sub>*jn*</sub> from the reference point *X*<sub>*j*</sub>:

$$d_{jn}(0) = ||X_{jn} - X_{j}|| \quad (2)$$

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 $\lambda_{max}$ as follows:

$$d(t) = d_{ave}(0)e^{\lambda_{max}t} \quad (3)$$

where $d_{ave}(0)$ is the initial separation of the neighboring points. Using Equation 3, one can write for the *j*<sup>th</sup> pair:

$$d_{j-ave}(i) \approx d_{ave}(0)e^{\lambda_{max}(i\Delta t)} \quad (4)$$

Applying a logarithmic function on both sides, we get:

$$\ln[d_{j-ave}(i)] \approx \ln[d_{ave}(0)] + \lambda_{max}(i\Delta t) \quad (5)$$

The above equation describes *M* parallel lines whose slopes are relatively proportional to $\lambda_{max}$. Finally, $\lambda_{max}$ can be accurately estimated by applying the least-squares fit on the averaged line defined as follows:

$$y(i) = \frac{1}{\Delta t} E[d_{j-ave}(i)] \quad (6)$$

where $E[d_{j-ave}(i)]$ is the average over all pairs of *j*.

## B. 0-1 TEST

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 $\phi(j)$ for $j \in [1, N]$, the translation variables $p_c(n)$ and $q_c(n)$ are computed as follows:

$$p_c(n) = \sum_{j=1}^{n} \phi(j) \cos(jc) \quad \& \quad q_c(n) = \sum_{j=1}^{n} \phi(j) \sin(jc) \quad (7)$$

for $n \in [1, N]$ and $c \in (0, \pi)$. 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 multidimensional 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:

$$D_c(n) = \lim_{N \to \infty} \frac{1}{N} \sum_{r=1}^{N} \{[p_c(j+n) - p_c(j)]^2 \\ + [q_c(j+n) - q_c(j)]^2\} \\ - V_{osc}(c, n) \quad (8)$$

for $n \le n_{cut}$ where $n_{cut} \le 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:

$$V_{osc}(c, n) = \lim_{N \to \infty} \frac{1}{N} \sum_{j=1}^{N} \phi(j)^2 \frac{1 - \cos nc}{1 - \cos c} \quad (9)$$

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:

$$K_c = \lim_{n \to \infty} \frac{\log \tilde{D}_c(n)}{\log(n)} \quad (10)$$

where $\tilde{D}_c(n) = D_c(n) - \min\{D_c(n)\}_{n=1}^{n_{cut}}$, 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, \pi)$. 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 \notin [0, 1]$.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. General structure of conceptor-driven network.</figcaption></figure>

### 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^{v}(n + 1) = \tanh(W^{*}x^{v}(n) + W^{in}v(n) + b) \quad (11)$$

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:

$$\sum_{n=1}^{L} (W^{out} x^{v}(n) - v(n))^2 \quad (12)$$

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:

$$x^j(n + 1) = \tanh(W^*x^j(n) + W^{in}p^j(n) + b) \quad (13)$$

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:

$$\sum_{j=1}^{K} \sum_{n=n_0}^{L} \|W^*x^j(n) + W^{in}p^j(n) - Wx^j(n)\|^2 \quad (14)$$

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)$:

$$x(n + 1) = C^j \tanh(Wx(n) + b) \quad (15)$$

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:

$$L(C^j) = E \|C^j x^j(n) - x^j(n)\|^2 + (\gamma^j)^{-2} \|C^j\|_{fro}^2 \quad (16)$$

where E is the mathematical mean, $\|\cdot\|_{fro}^2$ is the squared Frobenius matrix norm, and $\gamma^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 $\gamma^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 \approx I$, the first term converges to its minimal value. By contrast, the second terms vanishes in the case $C^j \approx O$. The minimization of the quadratic loss above leads to the following optimal solution of $C^j$:

$$C^j = R^j(R^j + (\gamma^j)^{-2}I)^{-1} \quad (17)$$

$R^j = E[x^j(n)x^j(n)^T]$ 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:

$$y(n) = W^{out}x(n) \approx p(n) \quad (18)$$

Figure 2 summarizes the described main steps of the ConDN process in flowchart form.
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Flowchart of time series prediction process by
conceptor-driven network.</figcaption></figure>

### 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 \in [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 \in s$, $i \in [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| \ge 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 \in s$ of different size $|q^1| > |q^2| > \dots > |q^{N_p}|$. A ConDN network must be fed with the signals of the same size. Thus, the size equalizing of the signals $q^i$, $i \in s$ from the first phase is mandatory. The size of the smallest signal $q^i$, $i \in s$ corresponding to the highest sampling factor $s_{i|s|} = N_p$ must be used for this operation. As the result of the equalizing operation, we have $N_S$ signals $p^i$, $i \in s$ of the same length, which can be presented in a matrix form $\mathbf{P} = [p^1 p^2 \dots p^{N_p}]$ as shown in Algorithm 1, line 11. In the next phase, the prepared input signals matrix $\mathbf{P}$ is fed to an already installed ConDN network

<caption>**Algorithm 1: ConDN-Based Chaos Detection (p, s)**</caption>

1. **ConDN-based Chaos Detection (p, s)**
   inputs: Time series to analyze **p**;
   Ordered sampling vector **s** = [1, 2, ..., $N_p$];
   output: Dynamics variation **Δ** or fitting parameter **a**
   // Time series sampling
2. **for** $i \in s$ **do**
3.   **[** $q^i \leftarrow$ Sampling (**p**, **i**)
4.   $\Rightarrow |q^1| > |q^2| > \dots > |q^{N_p}|$
5.   // Sampled time series length equalizing
6.   **for** $i \in s$ **do**
7.   **[** $p^i \leftarrow q^i$[1:|$q^{N_p}|$]
8.   $\Rightarrow |p^1| = |p^2| = \dots = |p^{N_p}|$
9.   // ConDN input matrix **P**
10.  **P** = [$p^1 p^2 \dots p^{N_p}$]
11.  // C - overall conceptor matrix
12.  **C** $\leftarrow$ ConDN (**P**), where **C** = [$C^1 C^2 \dots C^{N_p}$]
13.  // conceptor matrix $C^i$ corresponding to $p^i$, $i \in s$
14.  // S - overall SVD matrix
15.  **S** $\leftarrow$ SVD (**C**), where **S** = [$S^1 S^2 \dots S^{N_p}$]
16.  // SVD matrix $S^i$ corresponding to $C^i$, $i \in s$
17.  **T** $\leftarrow$ Trace (**S**), where **T** = [$t^1 t^2 \dots t^{N_p}$]
18.  $t^i = \text{tr}(S^i) = \sum_j^N S_{jj}^i$, $i \in s$
19.  **if** (one-shot) **then**
20.    **[** $\begin{cases}
21.      \text{ // } s = [1, 2] \\
22.      \Delta_1 = t^{s_2} - t^{s_1}; \\
23.      \Delta = [\Delta_1] \\
24.      \text{ return } \Delta;
\end{cases}$]
25.  **else**
26.    **[** $\begin{cases}
27.      \text{ // } s = [1, 2, \dots, N_p] \\
28.      \text{for } i \in [1, |s| - 1] \text{ do} \\
29.        \text{[} \Delta_i = t^{s_{i+1}} - t^{s_1}; \\
30.        \Delta = [\Delta_1, \Delta_2, \dots, \Delta_{|s|-1}] \\
31.        (a, \tau) \leftarrow \text{Fitting}(\Delta, s) \\
        \text{return } a;
\end{cases}$]
32.  **]$$

(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 $\mathbf{C} = [C^1 C^2 \dots C^{N_p}]$ are generated, where $C^i$, $i \in s$ is an $N \times N$ conceptor matrix corresponding to the input signal $p^i$, $i \in s$. 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 $\mathbf{C}^j$ undergoes a singular value decomposition (SVD) as shown in Algorithm 1, line 16:

$$C^i = U^i S^i U^{iT}, \quad i \in s \qquad (19)$$

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 \in s$ contain the precious information about the dynamics of the input signals $p^i$, $i \in s$ and are used to determine their either regular or chaotic behaviour. For each matrix $S^i$, $i \in s$, we compute then the value of
trace $t^i = tr(S^i)$, $i \in 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 \in 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 \in s$ are expected to increase with higher sampling factors $s_i$, $i \in 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 $\Delta$ can be computed. The elements of $\Delta$ are computed from the trace vector $T = [t^1 t^2 \dots 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 $\Delta$ is $N_S - 1$ where $\Delta_j = t^{s_{j+1}} - t^1$, $j \in [1, N_S - 1]$.

The results obtained in the dynamics variation vector $\Delta$ 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 $\Delta(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 $2p^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 $\Delta$ can be obtained. $\Delta$ has only one element $\Delta_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 $\Delta_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 $\Delta_1$.

In the multiple-shots sampling method, the sampling vector $s$ has more than two values $s = [1, 2, \dots, N_p]$ and is used for a much deeper analysis of the evolution of the obtained dynamics variation vector $\Delta$ with sampling $s$. For $N_S = |s|$ different sampling factors in $s$, the proposed approach will generate a dynamics variation vector $\Delta$ of size $N_S - 1$, as shown in line 29 of Algorithm 1. The evolution of $\Delta$ 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 $\tau$ in the latter case due to the rapid increase in the value of the trace in correspondence with the increase of the sampling factor.

$$ \Delta = a(1 - e^{-(s-1)/\tau}) \qquad (20) $$

## 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:

* Lorenz System:

$$ \begin{aligned} \dot{x} &= \alpha(y - x) \\ \dot{y} &= x(\beta - z) - y \\ \dot{z} &= xy - \gamma z \end{aligned} \qquad (21) $$

* Chua System:

$$ \begin{aligned} \dot{x} &= \alpha(y - x + \epsilon x - W(w)x) \\ \dot{y} &= x - y + z \\ \dot{z} &= -\beta y - \gamma z \\ \dot{w} &= x \end{aligned} \qquad (22) $$

where $W(w) = f(\delta, \zeta, w(t))$

* Rossler System:

$$ \begin{aligned} \dot{x} &= -y - z \\ \dot{y} &= x + \alpha y \\ \dot{z} &= \beta + z(x - \gamma) \end{aligned} \qquad (23) $$

* Chen System:

$$ \begin{aligned} \dot{x} &= \alpha(y - x) \\ \dot{y} &= \epsilon x - xz + \gamma y \\ \dot{z} &= xy - \beta z \end{aligned} \qquad (24) $$

* Lu System:

$$ \begin{aligned} \dot{x} &= \alpha(y - x) \\ \dot{y} &= -xz + \gamma y \\ \dot{z} &= xy - \beta z \end{aligned} \qquad (25) $$

* FO Rossler System:

$$ \begin{aligned} \frac{d^q x}{dt^q} &= -y - z \\ \frac{d^q y}{dt^q} &= x + \alpha y \\ \frac{d^q z}{dt^q} &= \beta + z(x - \gamma) \end{aligned} \qquad (26) $$

where $q \in [0, 1]$ is the derivation order.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. 3D presentation of the attractor of the conventional chaotic systems: (a) Lorenz, (b) Chua, (c) Rossler, (d) Chen, (e) Lu, and (f) the fractional-order systems (FO) Rossler, and (g) FO Chen.</figcaption></figure>

• FO Chen System:

$$ \begin{aligned} \frac{d^q x}{dt^q} &= \alpha(y - x) \\ \frac{d^q y}{dt^q} &= \epsilon x - xz + \gamma y \\ \frac{d^q z}{dt^q} &= xy - \beta z \end{aligned} \qquad (27) $$

where $q \in [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 is 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

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Multilayer composition of the studied pseudo-spin valve, with layers' thickness in nm in parentheses.</figcaption></figure>

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) = \pm 1$.
<figure><img src="image_6.png" /><figcaption>FIGURE 5. 2D plot of the trajectory of the NCVO's vortex gyration in the three different modes: (a) non-modulated, (b) modulated, and (c) chaotic.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Artificially reconstructed $x(t)$ signal without jitter effect for the (a) non-chaotic, (b) half-chaotic, and (c) chaotic cases.</figcaption></figure>

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_z(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><caption>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.</caption><thead><tr><th>System</th><th>Parameters</th><th>Chaotic Mode</th><th>Regular Mode</th></tr></thead><tbody><tr><td>Lorenz</td><td>{α, β, γ}</td><td>{10, 28, 8/3}</td><td>{10, 28, 3}</td></tr><tr><td>Chua</td><td>{α, β, γ, ε, δ, ζ}</td><td>{0.3, 0.8, 10, 13, 0.1, 1.5}</td><td>{0.3, 0.8, 30, 30, 0.25, 1}</td></tr><tr><td>Rossler</td><td>{α, β, γ}</td><td>{0.2, 0.4, 5.7}</td><td>{0.2, 0.4, 4.1}</td></tr><tr><td>Chen</td><td>{α, β, γ, ε}</td><td>{35, 3, 28, -7}</td><td>{20, 3, 28, -7}</td></tr><tr><td>Lu</td><td>{α, β, γ}</td><td>{36, 3, 20}</td><td>{10, 3, 10}</td></tr></tbody></table>

<table><caption>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.</caption><thead><tr><th rowspan="2">System</th><th rowspan="2">Parameters</th><th colspan="2">Order</th></tr><tr><th>Chaotic</th><th>Non-Chaotic</th></tr></thead><tbody><tr><td>FO Rossler</td><td>{α, β, γ} = {0.4, 0.4, 10}</td><td>o = 0.95</td><td>o = 0.87</td></tr><tr><td>FO Chen</td><td>{α, β, γ, ε} = {35, 3, 28, -7}</td><td>o = 0.96</td><td>o = 0.853</td></tr></tbody></table>

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 \in 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 $\Delta = [\Delta_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 $\Delta_1$ (close to 0), whereas in the chaotic case $\Delta_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 non-chaotic) evaluated in the last section, a half-chaotic one was
<figure><img src="image_8.png" /><figcaption>**FIGURE 7.** One-shot sampling: variation of trace values $t^i$, $i \in s = [1, 2]$ versus data length, for the chaotic (blue) and non-chaotic (orange) behaviours of the systems: (a) FO Rossler, (b) FO Chen, (c) Lorenz, (d) Chua, (e) Rossler, (f) Chen, and (g) Lu.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 8.** One-shot sampling: dynamics vector variation $\Delta$ for the chaotic (blue) and non-chaotic (orange) behaviours of the seven studied conventional systems.</figcaption></figure>

<table><caption>**TABLE 3.** Multiple-shots sampling: values of the parameters used in the ordinary differential equations of the three studied integer-order systems (Eqns. 23 to 25) in the chaotic, half-chaotic, and non-chaotic modes.</caption><thead><tr><th rowspan="2">System</th><th rowspan="2">Common Parameters</th><th colspan="3">Variable Parameters</th></tr><tr><th>Chaotic</th><th>Half-Chaotic</th><th>Non-Chaotic</th></tr></thead><tbody><tr><td>Rossler</td><td>{a, b} = {0.2, 0.4}</td><td>c = 2</td><td>c = 5.1</td><td>c = 5.7</td></tr><tr><td>Chen</td><td>{b, c, d} = {3, 28, -7}</td><td>a = 20</td><td>a = 27.7</td><td>a = 35</td></tr><tr><td>Lu</td><td>{b, c} = {3, 20}</td><td>a = 23</td><td>a = 23.5</td><td>a = 35</td></tr></tbody></table>

<table><caption>**TABLE 4.** Multiple-shots sampling: values of the parameters and derivation order used in the ordinary differential equations of the two studied fractional-order systems (Eqns. 26 and 27) in the chaotic, half-chaotic, and non-chaotic modes.</caption><thead><tr><th rowspan="2">System</th><th rowspan="2">Parameters</th><th colspan="3">Order</th></tr><tr><th>Chaotic</th><th>Half-Chaotic</th><th>Non-Chaotic</th></tr></thead><tbody><tr><td>FO Rossler</td><td>{α, β, γ} = {0.4, 0.4, 10}</td><td>o = 0.95</td><td>o = 0.9</td><td>o = 0.87</td></tr><tr><td>FO Chen</td><td>{α, β, γ, ε} = {35, 3, 28, -7}</td><td>o = 0.96</td><td>o = 0.86</td><td>o = 0.853</td></tr></tbody></table>

<table><caption>**TABLE 5.** Multiple-shots sampling: number of neurons of the ConDNs receiving the first component ($x(t)$) of the systems: FO Rossler, FO Chen, Rossler, Chen, and Lu; as well as the non-filtered and filtered signal $x(t)$ of NCVO, each of three cases: chaotic, half-chaotic, and non-chaotic.</caption><thead><tr><th rowspan="2">System</th><th colspan="3">Reservoir Size (N)</th></tr><tr><th>Chaotic</th><th>Half-Chaotic</th><th>Non-Chaotic</th></tr></thead><tbody><tr><td>FO Rossler</td><td>120</td><td>100</td><td>30</td></tr><tr><td>FO Chen</td><td>100</td><td>80</td><td>20</td></tr><tr><td>Rossler</td><td>120</td><td>100</td><td>50</td></tr><tr><td>Chen</td><td>120</td><td>50</td><td>30</td></tr><tr><td>Lu</td><td>120</td><td>120</td><td>50</td></tr><tr><td>Non-Filtered NCVO</td><td>400</td><td>200</td><td>80</td></tr><tr><td>Filtered NCVO</td><td>100</td><td>80</td><td>80</td></tr></tbody></table>

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 $\Delta(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 $\Delta$ 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
<figure><img src="image_10.png" /><figcaption>**FIGURE 9.** Multiple-shots sampling: dynamics vector variation $\Delta$ as a function of sampling factor $s$, with its fitting function, for (a) FO Rossler and (b) FO Chen.</figcaption></figure>

variation $\Delta(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 $\Delta(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 $\tau$ 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} \approx 1$) whereas at other times it is recognized as a non-chaotic signal ($K_{median} \approx 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

<figure><img src="image_11.png" /><figcaption>**FIGURE 10.** Multiple-shots sampling: dynamics vector variation $\Delta$ as a function of sampling factor $s$, with its fitting function, for (a) Rossler, (b) Chen, and (c) Lu.</figcaption></figure>

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 $\Delta(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 $\Delta$ in the chaotic case typically decreases with the addition of noise. In contrast, the value of $\Delta$ 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 $\Delta$ 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}$
<figure><img src="image_12.png" /><figcaption>FIGURE 11. Multiple-shots sampling: dynamics vector variation $\Delta$ as a function of sampling factor $s$ with its fitting function, for (a) simulated and (b) artificially generated (no jitter effect) signal $x(t)$ (horizontal position of vortex) of NCVO.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12. Multiple-shots sampling: slope values $a$ (see Eq. 20) for the Rossler, Chen, and Lu systems, as well as the simulated and artificially generated (or filtered) signal $x(t)$ of the NCVO, in three cases: chaotic (blue), half-chaotic (orange), and non-chaotic (green).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 13. 0-1 Test: value of $K_{median}$ for the FO Rossler, FO Chen, Rossler, Chen, and Lu systems, as well as the simulated and artificially generated (or filtered) signal $x(t)$ of the NCVO, in three cases: chaotic (blue), half-chaotic (green), and non-chaotic (orange).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 14. One-shot sampling: dynamics vector variation $\Delta$ versus SNR of added noise (dB), for the chaotic (blue) and non-chaotic (orange) behaviours of the systems: (a) FO Rossler, (b) FO Chen, (c) Lorenz, (d) Chua, (e) Rossler, (f) Chen, and (g) Lu.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 15. 0-1 Test: Variation of $K_{median}$ versus SNR of added noise (dB), for the chaotic (blue) and non-chaotic (orange) behaviours of the systems: (a) FO Rossler, (b) FO Chen, (c) Lorenz, (d) Chua, (e) Rossler, (f) Chen, and (g) Lu.</figcaption></figure>

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 = 30\,\text{dB}$ and $14\,\text{dB}$. 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 corre-
sponding 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 con-
clusion drawn in Section IV-B, where a chaotic system is

<figure><img src="image_17.png" /><figcaption>**FIGURE 16.** Multiple-shots sampling: (a) dynamics vector variation Δ as a function of the sampling factor and its fitting function, as well as (b) its corresponding slope for the five experimental packages of the NCVO.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**FIGURE 17.** One-shot sampling: dynamic vector variation (Δ) for the samples of the five experimental packages of NCVO.</figcaption></figure>

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*<sub>median</sub>
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 pre-
sented in Figure 19. As shown, the packages are distributed
irregularly along the axis of the applied magnetic field μ<sub>0</sub>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
<figure><img src="image_19.png" /><figcaption>FIGURE 18. 0-1 Test: value of $K_{median}$ of the five experimental packages of the NCVO device.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 19. Largest Lyapunov exponent: value of $\lambda$(GHz) of the samples of the five experimental packages of NCVO [32].</figcaption></figure>

package D is non-chaotic whereas those in the second part are chaotic. The value of the LLE ($\lambda$) is large for the samples of package B and the second part of package D, due to their chaotic nature. However the value of $\lambda$ is always positive even for the non-chaotic samples due to the jitter effect. Moreover, the value of $\lambda$ 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 ConDN-based method surpasses the 0-1 Test, the LLE, and the single-shot 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) multiple-shots 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.

## ACKNOWLEDGMENT

The authors would like to thank Thibaut Devolder and Damien Rontani for sharing experimental data and LLG calculations: Thibaut Devolder for initially measuring the NCVO time traces used in [32] and sharing them for our analyses, and Damien Rontani, for computing the LLEs, initially used in [32], and presented in Figure 19 of this work.

## REFERENCES

[1] M. Cattani, I. L. Caldas, S. L. D. Souza, and K. C. Iarosz, "Deterministic chaos theory: Basic concepts," *Revista Brasileira de Ensino de Física*, vol. 39, no. 1, Oct. 2016.

[2] U. Parlitz, "Estimating Lyapunov exponents from time series," in *Chaos Detection and Predictability*, 2016, pp. 1-34.

[3] M. Balcerzak, D. Pikunov, and A. Dabrowski, "The fastest, simplified method of Lyapunov exponents spectrum estimation for continuous-time dynamical systems," *Nonlinear Dyn.*, vol. 94, no. 4, pp. 3053-3065, Sep. 2018.

[4] J. Awrejcewicz, A. Krysko, N. Erofeev, V. Dobriyan, M. Barulina, and V. Krysko, "Quantifying chaos by various computational methods. Part 1: Simple systems," *Entropy*, vol. 20, no. 3, p. 175, Mar. 2018.

[5] J. Pathak, Z. Lu, B. R. Hunt, M. Girvan, and E. Ott, "Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data," *Chaos, Interdiscipl. J. Nonlinear Sci.*, vol. 27, no. 12, Dec. 2017, Art. no. 121102.

[6] E. G. Nepomuceno, S. A. M. Martins, M. J. Lacerda, and E. M. A. M. Mendes, "On the use of interval extensions to estimate the largest Lyapunov exponent from chaotic data," *Math. Problems Eng.*, vol. 2018, pp. 1-8, Jan. 2018.

[7] G. A. Gottwald and I. Melbourne, "The 0-1 test for chaos: A review," in *Chaos Detection and Predictability*. Berlin, Germany: Springer, 2016, pp. 221-247.

[8] J. R. Tempelman and F. A. Khasawneh, "A look into chaos detection through topological data analysis," *Phys. D, Nonlinear Phenomena*, vol. 406, May 2020, Art. no. 132446.

[9] D. Chen and W. Han, "Prediction of multivariate chaotic time series via radial basis function neural network," *Complexity*, vol. 18, no. 4, pp. 55-66, Mar. 2013.

[10] A. S. Weigend, *Time Series Prediction: Forecasting the Future and Understanding the Past*. Evanston, IL, USA: Routledge, 2018.

[11] A. D. Pano-Azucena, E. Tlelo-Cuautle, B. Ovilla-Martinez, L. G. de la Fraga, and R. Li, "Pipeline FPGA-based implementations of ANNs for the prediction of up to 600-steps-ahead of chaotic time series," *J. Circuits, Syst. Comput.*, vol. 30, no. 9, Jul. 2021, Art. no. 2150164.

[12] E. Tlelo-Cuautle, A. M. González-Zapata, J. D. Díaz-Muñoz, L. G. de la Fraga, and I. Cruz-Vega, "Optimization of fractional-order chaotic cellular neural networks by metaheuristics," *Eur. Phys. J. Special Topics*, pp. 1-7, Jan. 2022.

[13] H. Hassani, N. Alharbi, and M. Ghodsi, "Distinguishing chaos from noise: A new approach," *Int. J. Energy Statist.*, vol. 2, no. 2, pp. 137-150, Jun. 2014.

[14] A. R. Ismail, S. Jovanovic, S. Petit-Watelot, and H. Rabah, "Application of reservoir computing for the modeling of nano-contact vortex oscillator," *Electronics*, vol. 8, no. 11, p. 1315, Nov. 2019. [Online]. Available: https://www.mdpi.com/2079-9292/8/11/1315

[15] J.-P. Eckmann, S. O. Kamphorst, D. Ruelle, and S. Ciliberto, "Liapunov exponents from time series," *Phys. Rev. A, Gen. Phys.*, vol. 34, no. 6, pp. 4971-4979, Dec. 1986.

[16] S. Mehdizadeh, "A robust method to estimate the largest Lyapunov exponent of noisy signals: A revision to the Rosenstein's algorithm," *J. Biomech.*, vol. 85, pp. 84-91, Mar. 2019.

[17] H. Jaeger, "Conceptors: An easy introduction," 2014, arXiv:1406.2671.
[18] H. Jaeger, "Using conceptors to manage neural long-term memories for temporal patterns," *J. Mach. Learn. Res.*, vol. 18, no. 13, pp. 1-43, 2017. [Online]. Available: http://jmlr.org/papers/v18/15-449.html

[19] T. Furlanello, J. Zhao, A. M. Saxe, L. Itti, and B. S. Tjan, "Active long term memory networks," 2016, arXiv:1606.02355.

[20] S. A. Deka, D. M. Stipanović, B. Murmann, and C. J. Tomlin, "Long-short term memory neural network stability and stabilization using linear matrix inequalities," in *Proc. IEEE Int. Symp. Circuits Syst. (ISCAS)*, May 2019, pp. 1-4.

[21] E. N. Lorenz, "Deterministic nonperiodic flow," *J. Atmos. Sci.*, vol. 20, no. 2, pp. 130-141, 1963, doi: 10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2.

[22] C. Stegemann, H. A. Albuquerque, R. M. Rubinger, and P. C. Rech, "Lyapunov exponent diagrams of a 4-dimensional Chua system," *Chaos, Interdiscipl. J. Nonlinear Sci.*, vol. 21, no. 3, 2011, Art. no. 033105.

[23] K. M. Ibrahim, R. K. Jamal, and F. H. Ali, "Chaotic behaviour of the Rössler model and its analysis by using bifurcations of limit cycles and chaotic attractors," *J. Phys., Conf. Ser.*, vol. 1003, May 2018, Art. no. 012099.

[24] R.-D. Ene, C. Pop, and C. Petrişor, "Systematic review of geometrical approaches and analytical integration for Chen's system," *Mathematics*, vol. 8, no. 9, p. 1530, Sep. 2020.

[25] A. D. Pano-Azucena, E. Tlelo-Cuautle, G. Rodriguez-Gomez, and L. G. de la Fraga, "FPGA-based implementation of chaotic oscillators by applying the numerical method based on trigonometric polynomials," *AIP Adv.*, vol. 8, no. 7, Jul. 2018, Art. no. 075217.

[26] C. Li and G. Chen, "Chaos and hyperchaos in the fractional-order Rössler equations," *Phys. A, Stat. Mech. Appl.*, vol. 341, pp. 55-61, Oct. 2004.

[27] J. G. Lu and G. Chen, "A note on the fractional-order Chen system," *Chaos, Solitons Fractals*, vol. 27, no. 3, pp. 685-688, Feb. 2006.

[28] S. Petit-Watelot, J.-V. Kim, A. Ruotolo, R. M. Otxoa, K. Bouzehouane, J. Grollier, A. Vansteenkiste, B. V. de Wiele, V. Cros, and T. Devolder, "Commensurability and chaos in magnetic vortex oscillations," *Nature Phys.*, vol. 8, no. 9, pp. 682-687, 2012.

[29] A. Vansteenkiste, J. Leliaert, M. Dwornik, M. Helsen, F. García-Sánchez, and B. van Waeyenberge, "The design and verification of MuMax3," *AIP Adv.*, vol. 4, Oct. 2014, Art. no. 107133.

[30] S. Petit-Watelot, R. M. Otxoa, and M. Manfrini, "Electrical properties of magnetic nanocontact devices computed using finite-element simulations," *Appl. Phys. Lett.*, vol. 100, no. 8, Feb. 2012, Art. no. 083507.

[31] J. Létang, S. Petit-Watelot, M.-W. Yoo, T. Devolder, K. Bouzehouane, V. Cros, and J.-V. Kim, "Modulation and phase-locking in nanocontact vortex oscillators," *Phys. Rev. B, Condens. Matter*, vol. 100, Oct. 2019, Art. no. 144414, doi: 10.1103/PhysRevB.100.144414.

[32] T. Devolder, D. Rontani, S. Petit-Watelot, K. Bouzehouane, S. Andrieu, J. Létang, M.-W. Yoo, J.-P. Adam, C. Chappert, S. Girod, V. Cros, M. Sciamanna, and J.-V. Kim, "Chaos in magnetic nanocontact vortex oscillators," *Phys. Rev. Lett.*, vol. 123, no. 14, Oct. 2019, Art. no. 147701.

[33] A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, "Determining Lyapunov exponents from a time series," *Phys. D, Nonlinear Phenomena*, vol. 16, no. 3, pp. 285-317, 1985.

<figure><img src="image_21.png" /></figure>

**ALI RIDA ISMAIL** received the B.S. degree in electronic engineering from Lebanese University, in 2016, the M.S. degree in control systems and information technology from Grenoble Alpes University, in 2018, and the Ph.D. degree in electronics and circuits design from the University of Lorraine, in 2022. He is currently working as a Postdoctoral Researcher with the Jean Lamour Institute, Nancy, where his research interests include machine learning, signal processing, data

simulation, and features extraction.

<figure><img src="image_22.png" /></figure>

**SLAVISA JOVANOVIC** (Member, IEEE) received the B.S. degree in electrical engineering from the University of Belgrade, Serbia, in 2004, and the M.S. and Ph.D. degrees in electrical engineering from the University of Lorraine, Nancy, France, in 2006 and 2009, respectively. From 2009 to 2012, he was with the Diagnosis and Interventional Adaptive Imaging Laboratory (IADI), Nancy, as a Research Engineer working on MRI-compatible sensing embedded systems. Then, he joined the

Faculty of Sciences and Technologies and Jean Lamour Institute (UMR 7198), University of Lorraine, where he is currently an Associate Professor. He is the author and coauthor of more than 50 papers in conference proceedings and international peer-reviewed journals, and he holds one patent. His research interests include energy harvesting circuits, neuromorphic architectures, reconfigurable network-on-chips, and algorithm-architecture matching for real-time signal processing.

<figure><img src="image_23.png" /></figure>

**SÉBASTIEN PETIT-WATELOT** received the Ph.D. degree, in 2007. His Ph.D. works on the study of spin-transfer torque realized with the SPINTEC Laboratory. He is currently a Researcher in condensed matter physics. After his Ph.D. degree, he worked as a Postdoctoral Fellow at CEA Saclay (SPEC), Nanomagnetism Group. During this period, he studied magnetotransport properties in mesoscopic devices. In 2010, he joined the Institut d'Electronique Fondamentale (now C2N)

to work on spin torque oscillators. Since 2012, he has been an Associate Professor at Lorraine University (UL) and carries out his research activities with the Nanomagnetism and Spintronic Group, Institut Jean Lamour (UL/CNRS). His field of expertise addresses spintronic, nano-magnetism, and magnetization dynamics. Taking advantage of the fundamental knowledge, he has acquired on nonlinear dynamics, he devotes a part of his research time to the study of chaotic-based systems.

<figure><img src="image_24.png" /></figure>

**HASSAN RABAH** (Senior Member, IEEE) received the M.S. degree in electronics and control engineering and the Ph.D. degree in electronics from Henri Poincaré University, Nancy, France, in 1987 and 1993, respectively. He became an Associate Professor in electronics microelectronics and reconfigurable computing with the University of Lorraine, Nancy, in 1993, and a Full Professor, in 2011. In 1997, he joined the Architecture Group of LIEN, where he super-

vised research on very-large-scale integration implementation of a parallel architecture for image and video processing. He also supervised research on the field-programmable gate array (FPGA) implementation of adaptive architectures for smart sensors in collaboration with industrial partners. He participated in several national projects for quality of service measurement and video transcoding techniques. He joined Institut Jean Lamour, University of Lorraine, where he held the position of the Head of the Measurement and Electronics Architectures Group, from 2013 to 2021. His current research interests include the design and implementation of FPGA-based embedded systems with an emphasis on power optimization, video compression-decompression and transcoding, compressive sensing, sensor networks, and machine learning. He has been a program committee member and organized special sessions for a number of conferences. He is currently the Vice-Chair of the IEEE Instrumentation and Measurement France Chapter.

<figure><img src="image_25.png" /></figure>