Detection of Generalized Synchronization using Echo State Networks

Detection of Gen eralized S ynchronization using E cho State Networks D Ibáñez-Soria 1 , J Garcia- Ojalvo 2 , A Soria-Frisch 1 , G Ruffini 1 ,3 1 Starlab Barcelo na S.L., Neuroscience Resear ch Business Unit 2 Department of Experi mental and Health Science s , Universi tat Pompeu Fabra, Barcelona, Spain 3 Neuroelectrics Corpor ation E-mail: david.ibanez@s tarlab.es Abstract. Generalized synchronization between coupled d ynamical systems is a phenomenon of relevance in ap plication s that range from secure communications to p hysiological modelling . Her e we test the capabilities of reservoir computing and , in particular , echo state networks for the detection o f generalized s ynchronization. A nonlinear d ynamical system consisting of t wo coupled Rössler chaotic attractors is used to generate temporal series consisting of time-loc ked generalized synchronized sequences interleaved by unsynchro nized ones. Correctly tuned , echo state net works are able to efficiently d iscriminate bet ween unsynchronized and synchronized seq uences . Compared to other state - of -the-art techniques of synchronization detection, the online capabilities of the proposed ESN based meth odo logy make it a pro mising c hoice for real -time applicatio ns aiming to monitor dynamical synchronization change s in co ntinuous signals. Keywords: Reservoir Computing, Generalized Synchr onization, Echo State Network s, Dynamical Systems, Chaotic S ystems , R össler Attractor . 1. Introduc tion In its everyday use , the concept of synchronization is comm only taken to i mply identical behavio r between two i nteracting systems. However, when the affected systems are chaotic, more sophisticated synchronization for ms can exist, such as ph ase synchronization ( where only the phase but not the amplitude of two chaotic osc illations agree with each other) [1], lag synchronization (where one of the chaotic systems follow s th e other with a certain delay ) [2], or g eneralized synchronization (where the states of the two systems are functiona lly – but not identically – relat ed in a nontrivial, in general nonlinear, m anner) [3,4 ]. Generalized synchroniz ation, in p articular, h as prov en to be a relevan t feature in the ana lysis of neurologic al disorders such as Alzheim er’s disease [5] and brain tum ors [ 6 ]. Over the years, several methods have bee n proposed for the detection of generalized synchronization, including the replica method [ 7], t he synchronizat ion likelihood approach [8], and the mutual false nearest neighbor method [ 4], among others. T hese methods are usually computationally costly, cannot be applied in a cont inuous manner, and in some cases, suffer from different biases [9]. Here we propose a machine-learni ng-based approach that enables the online detection of generalized synchronization in a n effectiv e manner. The m ethod relies on a recurren t neural network to provide the necessary fading memory that allows processing dy namical signals. Unlike feedforward neural networks, in which a static input-output mapping is applied, recurrent neural network s (RNNs) hav e cyclic connection s that pr ovide mem ory, implementing a system with dynamical capabilities [ 10]. T raining RNNs has traditionally be en computationally more expensive than training feedforward networks. Echo State Networks [ 11] (ESN) and L iquid State Machines [ 12 ] (LSM), which were developed inde pendently and simultaneously, and the more recent decorrelation learning rule for RNNs [ 13 ], are complementary approaches for designing , t raining and analyz ing RNN s within a m ethodological fram ework know n as Reservoir Com puting (RC). RC is based on the principle that if the ne twork possesses c ertain alge braic pro perties the superv ised training of al l connections is not necessary. Only supervised training of readout weight is sufficient to obtain optimal classification per formance i n many tasks. Here we exam ine the capabilities o f ESNs in the task of detec ting generalized sy nchronization chang es in synthetic tem poral series based on the coup ling of two chaotic systems. We wi ll use in particular the well-know n Rössler attractors. The paper is structured as follows. In section 2 we explain t he methodology fol lowed to construct generalized synchr onized sequences using coupled chaotic R össler oscillators. In section 3 we provide an overview of the echo state network architectu re and their key training parameters . T he procedure used to study ESN generalized synchronization detection is presented i n section 4 and the obtained r esults i n section 5. We c onclude wi th a discussion in section 6. 2. Generation of in silico tim e series The main object ive of this work is to explore the capabilities of ESNs f or discriminating between generalized synchronized tim e-series and unsynchroniz ed sequences. In this se ction, we descr ibe ho w synchronized ch aotic attra ctors are cons tructed. To t hat end we follow the un idirectional ( master- slave) cou pling of two Rössler oscillators proposed by Rulkov et al [4]. A Rössler oscillator is a dynamical system defined by three n on-linear o rdinary differential equat ions that exhibit chaotic dynamics. T he Rössler oscillator d escribed by t he state v ariables  presented in (1) has been adopte d as the driving system, while  in (2) constitu tes the coupled res ponse system:   󰇗   󰇛     󰇜 (1)   󰇗           󰇗      󰇛    󰇜   󰇗    󰇛      󰇜  󰇛     󰇜 (2)   󰇗          󰇗      󰇛   󰇜 The m anifold      ,       and       contains the t rajectories of sync hronized oscillations. The conditional Lyapunov exponents (CLE) of the response system char acterize its a symptotic local stability [14]. It has bee n proven that if all CLEs are negative, the response system is stable and shows synchronization [ 15]. Reservoir computing has pr oved i ts capabilities for the assessment of positive and negative Ly apunov exponents in high dimensio nal spatio- temporal chaotic system s [16]. For a = 0.2 and b = 5.7 Rulkov et al [4] proved that trajectories of sy nchronized motions were stabl e for a coupling factor of g = 0. 2, becoming unstable for g = 0.15. Figure1A shows a plot of   󰇛󰇜    󰇛󰇜 for     , and Figure 1B for     . In A we obser ve a straight line that denotes identical synchronization between driving and response system , while i n B unsynchronized oscillations appear. By correctly tuning the coupling factor g it is thus possible to force synchronizati on between the coupled systems. To achieve generalized synchronization Rulkov proposes the nonline ar transform ation of the response system 󰇛 󰇜 into  󰇛  󰇜 presented in equat ion 3.      (3)                    Figure1C p lots   󰇛󰇜    󰇛󰇜 for    , where a straig ht line denoting ide ntical sy nchronization is no longer observed. Figure1D shows   󰇛󰇜    󰇛󰇜 for     . In C, even though a complex relationship between variables is observed, synchronization is not lost since only a non- linear transformation was ap plied. I n this case   󰇛󰇜 and   󰇛󰇜 present generalized synchron ization [4]. Figure1 Dynamics of t wo coupled Rössler oscillators: ( A)   󰇛 󰇜    󰇛󰇜 for    , ( B ):   󰇛 󰇜    󰇛󰇜 for      , ( C)   󰇛 󰇜    󰇛 󰇜 for    and ( D)   󰇛󰇜     󰇛 󰇜 for     . Using the behavior reported in Fig. 1, a continuous temporal si gnal consisting of a series of synchronized sequences int erleaved by unsynchronized ones has been cons tructed. For t his purpose, a the time-v arying square wav e coupling function g( t) is used:  󰇛  󰇜    󰇡    󰇡 󰇡     󰇢󰇢  󰇢    (4) where T defines the up or down state duration in time units. T he model has been solved using a 2 nd -3 rd order Rung e-Kutta m ethod (i mplem ented with the Matlab function ode23) with fixed integration step of 0.4 time uni ts. The synchronization d etection perfo rmance is evaluated u sing time series of length 10  time units and an up/down state duration T of 10  time units. The resultant signal thus concatenates, depending on the initial time, at least 5 generalized synchronized seq uences followed by 5 unsynchronized on es. A C B D 3. ESNs for generalized syn chr onization det ection We now discuss the ESN architecture s explored, and the tests to assess generaliz ed synchronization . Artificial Neural Networks (ANN) such as the multilayer- perceptron (MLP) present a feed-forward structure where the i nform ation flows from the input nodes, through the hidden nodes and to the output nodes [17] . Based on t he Representation T heorem, ANNs ( and in particula r MLPs) are able t o approximate an y given ar bitrary functi on. T his static input-output architecture makes them suitable for the analysis of stationary problems, but it is in general not adequate to deal with dynamical time- dependent problems. Recurrent neural networks (RNNs) i ncorporate cyclic connections that allow the system to pr ovide memory capability to the network, and therefore to encode time-dependent information. This addition transform s the network into a dynamical system . The network keeps in its internal states non-linear transformation s of the input history (fading memory) allowing it to process information w ith temporal context [18 ]. RNNs training ha ve t raditionally been com putationally expe nsive t o train becau se of their cy clic non- linear nature [22]. Reservoir Computing (RC), a methodolog ical framework to understand, train, and apply Recurrent Neural Networks (RNNs), was proposed independently and simultaneously with the development of Ec ho State Networks (ESNs) [11] and Liq uid State Machines (LSMs) [ 19]. The fundamental principle of E SN is that if the network presents a certain algebraic property known as the Echo State Property (ESP), only supervised training of readout connections is nee ded [20]. Tra ined output units combine the internal RNN dynam ics i nto desired outputs. The untrained RNN is called the dynamical reservoir (DR) and is formed by input, hidden and backpropagation connections. T he echo state prope rty ensures that the r eservoir state does not depend in the long ter m on t he initial conditions, which are thus forgotten as time passes [21 ]. The pr ecise values of input and i nterior weights are irrelevant to the echo state property. These weight values can be randomly gene rated accord ing to some par ameters [ 22], among which it is worth pointing out the so-called spectral r adius. The spectral radi us, calculated as the largest ab solute eigenvalue of the internal connections matrix, determines t he timescale of the reservoir. I n practice, if the spe ctral r adius is smaller than one, the echo state property holds for m ost application s [ 22] . Some unlikely exceptions to this rule have been proposed, how ever [23]. Additionally, in som e situations the echo state prope rty may also hold for values of the spectr al radius larg er than one [21]. The spectral radius and som e other key parameters of the ESN , mainly the size of the reservoir and t he input scaling, rule t he dynam ical beha vior of the netw ork [22 ]. A small spe ctral radius induces a faster response, while a larger spectral radius is m ore suitab le for tasks requiri ng longer fading mem ory. In general, the size of ec ho state networks can be large r compared to other recurrent neural network approaches [24]. In general and given sufficient data, the larger the network size, the better it can learn complex dynam ics. I n case of dat a shortage, large reservoirs can lead to overfitting and m ake the network present poor predictive performance. I n ESN s, inp ut scaling i s implemented multiplying every input sample by t he same scaling factor. T he input scaling determines the degree of nonlineari ty in the reservoir: while linear tasks require small i nput scaling factors, tasks with complex dynamics demand larger inpu t scaling v alues. In our case, we configured a single ESN network to have two input node s and a single output node. The signals   󰇛󰇜 and   󰇛󰇜 from the Rössler osci llator feed the input nodes. The op timal param eters of the ESN have been de termined throug h exhaustive search in a grid. Conc retely we have used following value grid s: number of in ternal un its - (5, 10, 25, 50, 75, 100, 200, 300, 400, 500), spectral radius - (1, 0.7, 0.5, 0.2, 0. 1, 0.01, 0.001), and inpu t scaling - (0.001, 0.01, 0.1, 0.5, 1, 5, 10, 25, 50, 75, 100). The number of internal units has to be large enough f or the system to l earn the com plex dynamics associate d with generalized synchroniza tion , but not too large so the system generalizes well enough. The target of the ESN is to discriminate between synchronized and unsynchron ized sequences. For this binary classi fication p roblem, during t raining the ESN output of synchronized sequences is teacher-forced to 1 and to -1 for unsy nchronized sequences. 4. Generalized Synch ronization Det ection The performance of each input scaling, network size and spectral radius tuple has bee n individual ly assessed. A Rössler signal formed by a series of 5 generalized synchronized sequences followed by unsynchronized sequences, as described in Section 2, is used for training. A different signal of the same characteristics and the same number of sa mples but starting at different initial time is used for testing the tuple performance. Hence we are using for performance evaluation a hold -out validation scheme with 50% of traini ng samples and 50% of different test ones. According to the chaot ic nat ure of the Rössler dynamics, by starting at different instants, althoug h the two attractors present a comm on pattern in spa ce state, they wil l develop trajectories that exponenti ally separates in an unsynchronized manner [25 ]. T he ESN t est output has been smoothed using a 10000-sample moving average window before perform ance evaluation. The Receiver Operating Characteristic (ROC) Curv e has been largely used as an effective m ethod to evaluate the perform ance of binary classi fication systems. In such two-class prediction probl ems outcomes are labeled as positive and negativ e class. ROC curves show the trade-off between sensitivity and speci ficity as function of a varying dec ision threshold. I n our case, the output test samples correspondi ng to synchronized sequence s are labeled as the positive class and unsynchronized ones as the negative class. The false positive rate (specifi ci ty) and true positive rate (sensitivity) of the ESN output is then calculated f or the all possible values of the decision threshold as applied on the ESN output after sm oothing. Figure 3C exemplary shows the ROC curves calculat ed for differen t smoothing window lengths. The area under the ROC curv e, or simply Area Under the Curve (AUC), measures the probability of the system of ri ght ranking the positive and negativ e clas s sam ples. An area of 1 m eans that a ll samples we re correctly classified while an area of 0.5 represents random classification. With the objective of r educing the random effect introduced by the reservoir initialization, each tuple p erformance evaluat ion pro cess is repeated 5 independent times an d i ts av erage AU C is used to ass ess the tuple’s discrim ination perform ance. The largest AUC value (0.85) was found for a spectral radius equal to 0.01, an input scaling of 25, and 500 i nternal units. The r eserv oir size appears to be a key training par ameter. Figure 2A displays the AUC score of the best spectral-radius/inpu t-scaling tuple as function of the number of i nternal units . The AUC score increase s asymptotically with the reservoi r siz e. According to these results, it is necessary to have at least 100 internal units to achieve an AUC performance larger t han 0.8. The performance of more than 500 internal units was not evaluated due to com putational restrictions. The input scaling constitute s also a training parameter that sens itively affects perf ormance. Figure 2B represents the AUC score for the best spectral radius as a function of t he input scaling. The AUC performance si gnificant ly improv es in the range [5, 50]. This large input scaling factor suggests a large non- linearity of the regression problem under evaluation. On the other hand, in t his general ized synchronization detection scenario the discrimination performance robustly behaves with respect to the sp ectral radius as shown i n Fi gure 2C, wh ere the AUC score obtained f or the best input scaling is plotted as a function of the spectral radius. Figures 2D, 2E and 2F corroborate the aforementioned observations, with the colour map r epresen ting the AUC for 100, 300 and 500 internal units respectively. Figure 2 Generalized synchronizati on detection performance parameter anal y sis. (A) AUC score of the best spectral-radius/input-scal ing t uple as function of t he reser voir size . (B) AUC score as function of the input scaling ( C) AUC score as a function of the s pectral radius. AUC ( colormap value) variation with respect to spectral radius (X-axis) and input scaling (Y -axis) for diff erent reser voir dimensi ons: (D) 100 units, (E) 300 units, and (F) 50 0 units. Figure 3A depicts an exam ple of t he output of t he outperforming ESN parameter tuple before averaging, which achieves a 0.54 AUC. The o utput has been s caled in the [-1, 1] for the sake of visualization. The dotted black line represents the aimed ESN Output , where 1 represents synchronized sequenc es and -1 unsynchronized ones. Despite of the small AUC obtained, a substantial difference between synchr onized and unsynchronized intervals can be observed. De synchronized samples present a larger high-freque ncy amplitude response. In Figure 1B , where the coupled attractors are considered to be unsynchronized, we can observe that many samples lie around the straight line correspond ing to synchronization. T he ES N discrimination o f t hese sequences appears to be more difficu lt, as within them not all samples seem to be totally unsy nchronized. T his behaviou r motivates the use of technique s that sm ooth the ESN output, in order to im prove the detection performance. To smooth the ESN output we have used a simple moving average approach implem ented as the unweig hted mean of the previous W samples, where W states for the window length in number o f sam ples. Figure 3B shows the ESN output after W =10000 samples averag ing. In this c ase w e can observe a better discrim ination between generalized synchronized and unsynchronized sequences achieving an AUC of 0.85. Figure 3C shows the receivers operator curve computed for W = 1, 100, 500, 1000, 5000 and 10000 samples. As expec ted, the AUC increase s wit h the averaging window length. Figure 3 (A) Unfiltered ESN output for spectral radi us 0.01, input scaling 25, and 5 00 internal units (B) ESN output after 1000 samples moving a verage for spectral radius 0.01, input scalin g 25 and 500 internal units . (C) Receivers operat ors curve calculated for averaging wi ndow s of 1, 50 0, 1000, 5000 and 10000 sample s. 5. Discussion and Conclusio ns We have presented a reliable, light training generalized synchronization dete ction m ethodology based on echo stat e networks. An optimal param eterization of ESNs was able to disc riminate between time- locked generalized synchronized sequences from unsynchronized ones delivering an area under the curve abov e 0.85 . Unlike other GS detection m ethods t hat canno t be applied in a continuous fashion, artificial neural networks update its output with every input sample. ESN t herefore prove to be an ideal choice to dev elop applications capable of monitoring generalized synchroniz ation changes in real-tim e. An appropriate tuning of ESN parameters has proved necessary for achieving a good discrim ination performance. The r eservoir size turns out t o be a fundam ental training parameter along with the input- scaling. Acco rding to our re sults, a minim um of 100 inte rnal units is required t o achiev e good performance and thus learning the generalized synchroniz ation complex dynamics between Rössler oscillators . Input scaling d etermines t he degree of non- linearity in t he reserv oir. An input scaling between 5 and 50 improv ed ESN discrimination capabilities. This large input scaling factor suggests the expected high non-linea rity nature of generalized sy nchronization between the two coupled chaotic systems. We expe ct from the theoretical study presented herein to be tter understand th e properties o f ESN and the r ole of its different parameters. This can increase the number of applications of th is ANN approach for the analysis of time se ries. In this context we are aiming to apply ESN for the analysis of electroencephalog raphy data, which is used for brain monitoring, as we will show in futu re communications. 6 Acknowledgmen ts We want t o thank the E uropean Union’s Horizon 2020 research and innovation programm e who funded STI PED proj ect under grant agreement No 731827. We also want to thank the support from ICREA Academia program and the Spanish Ministry of Econom y and Competitiveness and FEDER (project FIS2015-66503- C3 -1-P and Maria de Maeztu Programme for Units of Excellence in R&D, MDM-2014- 0370). We also k indly thank developers o f Oger toolbox. 7 References [1] Rosenblum, M. G., Pikovsk y, A. S., & Kurths, J . (1996). Phase s y nchr onization of chaotic oscil lators. Ph y sical review letters, 76(11), 1804. [2] Rosenblum, M. G., P ikovsky, A. S., & Kurths, J . (1 997). From phase t o lag synchroni zation in coupled chaotic oscillators. P hysical Re view Letters, 78(22), 41 93. [3] Afraimovich, V. S., Verichev, N. N., & Rabinovich, M. I. (1986). Stochastic synchronizati on of oscillations in dissipati ve systems. Radiofizika, 29, 1050 -1060. [4] Rulkov, N. F., Sushc hik, M. M., Tsimring, L. S., & Abarbanel, H. D. (1995). Generalized s yn chronizati on of chaos in directionall y coupled chaotic s y stem s. Ph y sical Re view E, 51(2), 980. [5] Stam, Cornelis J., et al. Journal of Clinical Neuroph y siolog y 19.6 (2002): 562- 574 [6] Bartolomei, Fabrice, et al. Clinical neur ophysiology 1 17.9 (2006): 2039- 2049]. [7] Abarbanel, Henr y DI, Nikolai F. Rulkov, and Mi khail M. Sushchik. Ph ysical Review E 53.5 (1996): 45 28 [8] Stam, C. J., and B. W. Van Di jk. P hysica D: Nonlinear Phe nomena 163.3 (2002): 23 6-251 [9] Pereda, E., etal. Ph y sica D: N onlinear Phenomena 14 8.1 (2001): 147 - 158].4 [10] Buonoma no, D. V., & Maass, W. ( 2009). State- dependent computati ons: s patiotemporal processin g i n cortical networks. Nature Re views Neuroscie nce, 10(2). [11] J aeger, H. (2001). T he “ech o state” appr oach to analysing and training recurren t ne ural n etworks -with an erratum note. Bonn, Germany: German National Research Center for In formati on Technolog y GMD Technical Report, 148(3 4), 13 [12] Natsch läger, T., Maass, W., & Mar kram, H. ( 2002). The" liq uid computer ": A novel strategy f or real -time computing on time series. Special issue on Foundations of Information Processing of TELEMATIK, 8(LNMC-ARTICLE -2002-005), 39- 43 [13] Steil, J. J. (2004, July ). Backpropagati on-dec orrelation: online recurrent learning with O (N) complexity. In Neural Netw orks, 2004. P roceedings. 2004 IEEE Interna tional J oint Conference on (Vol. 2, pp. 8 43 - 848). IEEE. [14] L. M. Pecora, T.L. Carr oll, Phys. Rev. Lett. 64 (1990) 821. * * [15] P arlitz, U., & Hornstein, A. (2005, August). Detecting ge neralized s y nchr onization from time series. In Fifth EUROMECH N onlinear D ynamics Conference ENOC - 2005 (pp. 1174-1181). [16] J aideep Patha k, Zhixin Lu, Brian R. Hunt, Michelle Gir van, Edward Ott . (2017). Usi ng Machine Learning to Replicate Chaotic A ttractors and Calculate Lyapunov Exponent s from Data. arXi v:1710.07313 . [17] Me hrotra, Kis han, Chil ukuri Moh an, and Sanja y Ranka. Elemen ts o f Artific ial Neural Netw orks. Boston: MIT Press, 1997. [18] Lukoše vičius, M., & Jaeg er, H. (2009). Reservoir co mputing approaches to recurrent neural network training. Computer Science Review, 3(3), 1 27- 149. [19] Maa ss, W., Natschl 󰇘 a ger, T., Markram, H.: Real-time c omputing with out stable states: a new framework for neural computation ba sed on perturbations. Neural Com putation 14( 11), 2531 – 2 560 (2002). [20] Lukoše vičius, M., Jae ger, H. , & Schrauwen, B. (2012). Reservoir computing trends . KI -Kü nstliche Intelligenz, 26(4), 365- 371. [21] Yildiz, I. B., Jae ger, H., & Kiebel, S. J. (2012). Re -visiting the ec ho state pr operty. Ne ural networ ks, 35, 1-9 [22] J aeger, H. ( 2002). Tutorial on traini ng recurrent neural net works, c overing BPP T, RTRL, EK F and the" echo state network" appr oach (V ol. 5). GMD- Forschungszentrum Informat ionstechnik. [23] Lukoše vičius, M. (2012). A prac tical guide to a pplying ec ho state netw orks. In Neural networks: Tric ks of the trade (pp. 659-6 86). Springer Berli n Heidelberg [24] Fabia n Triefenbach, Azarakhsh Jalalvand, Benjamin Schr auwen, and Jean -P ierre Martens. Phoneme recognition with large hierarc hical reservoirs. In A dvanc es in Ne ural Information Pr ocessing Sy stems 23 (NIPS 2010), pages 230 7 – 2315. MI T Press, Cambridge, MA, 201 [25] Wu, C. W. (2002). S ynchronization in coupled chaotic circuits and systems (Vol. 41). Si ngapore: World Scientific.