Multiplex visibility graphs to investigate recurrent neural networks dynamics

2016·Arxiv

1Machine Learning Group, Department of Physics and Technology, University of Tromsø, 9019 Tromsø, Norway 2Department of Computer Science, College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter EX4 4QF, United Kingdom 3Department of Electronics, Information, and Bioengineering, Politecnico di Milano, 20133 Milan, Italy 4Faculty of Informatics, Universit`a della Svizzera Italiana, 6900 Lugano, Switzerland *ﬁlippo.m.bianchi@uit.no ABSTRACT

1Machine Learning Group, Department of Physics and Technology, University of Tromsø, 9019 Tromsø, Norway 2Department of Computer Science, College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter EX4 4QF, United Kingdom 3Department of Electronics, Information, and Bioengineering, Politecnico di Milano, 20133 Milan, Italy 4Faculty of Informatics, Universit`a della Svizzera Italiana, 6900 Lugano, Switzerland *ﬁlippo.m.bianchi@uit.no ABSTRACT

A recurrent neural network (RNN) is a universal approximator of dynamical systems, whose performance often depends on sensitive hyperparameters. Tuning them properly may be difficult and, typically, based on a trial-and-error approach. In this work, we adopt a graph-based framework to interpret and characterize internal dynamics of a class of RNNs called echo state networks (ESNs). We design principled unsupervised methods to derive hyperparameters configurations yielding maximal ESN performance, expressed in terms of prediction error and memory capacity. In particular, we propose to model time series generated by each neuron activations with a horizontal visibility graph, whose topological properties have been shown to be related to the underlying system dynamics. Successively, horizontal visibility graphs associated with all neurons become layers of a larger structure called a multiplex. We show that topological properties of such a multiplex reflect important features of ESN dynamics that can be used to guide the tuning of its hyperparamers. Results obtained on several benchmarks and a real-world dataset of telephone call data records show the effectiveness of the proposed methods.

Introduction

A current research trend aims at investigating complex time-variant systems through graph theory, by considering suitable features associated with vertices and edges1. Of particular interest are those systems that also perform a computation when driven by an external input signal. An example is that of artificial RNNs2–4, which are computational dynamical systems whose link with physics and neurosciences dates back to the ’80 with some pioneering works from Jordan5 and Amit et al.6. Nowadays, RNNs are gaining renewed interest in neuroscience due to their biological plausibility7–10 and in computer science and engineering for their modeling ability11,12. RNNs are capable to generate complex dynamics and perform inference based on current inputs and internal state, the latter maintaining a vanishing memory of past inputs13,14.

Let us consider trajectories describing the evolution of a dynamical system in state space, e.g., the space containing all possible system states. As an example, in Fig. 1 we show the trajectories of a dynamical system operating in ordered (left) and chaotic (right) regimes, whose state is defined by the values of variables . Discriminating between order and chaos is of fundamental importance for investigating the properties of a dynamical system. It emerges that, such properties are manifested in the memory of the dynamical system and in the divergence rate of its state trajectories15. Fading memory is a desirable property of a dynamical system, and is characterized by ordered (contractive) dynamics. This is also referred to as the echo state property in the reservoir computing community16 and ensures that the current state/output of the system depends only on a finite number of past states/inputs17. At the same time, a high divergence rate between state trajectories (a property of chaotic dynamics) is also a desirable feature of RNNs. Results show that RNNs operating in a chaotic regime are able to produce meaningful patterns of activity18 and a balance must be struck in order to meet both properties. As a consequence, a (computational) dynamical system has to operate on the transition between order and chaos, in a region of the controlling parameter space called “edge of criticality”. On the edge of criticality, RNN internal dynamics becomes richer, meaning that neuron activations become heterogeneous15. Such a diversity both improves memory capacity of RNNs as well as their capability to reproduce complex target dynamics19,20. The notion of edge of criticality permeates several complex systems21–23, including random Boolean networks24 and populations of spiking neurons25. It is furthermore well-known that several complex systems spontaneously adapt to operate towards the edge of criticality, according to a mechanism known as self-organizing criticality in the statistical physics literature26.

In this work, we study ESNs27,28, a class of RNNs characterized by a large, untrained recurrent layer producing outputs by means of a linear (trained) readout layer. As such, ESNs are trained by optimizing readout weights only through a fast linear

Figure 1. Trajectories in the state space of a system described by the evolution of variables over time. The behavior of the system is controlled by a parameter (not specified), which we assume to be able to produce ordered (left plot) and chaotic (right plot) dynamics, depending on its value. In the ordered regime, the trajectories converge to the same fixed point when starting from 2 different initial conditions (red dots). In this configuration, the system is characterized by a fading memory of previous states. In the chaotic regime, on the other hand, trajectories remain separated if the system starts from different initial conditions.

regression, rather than learning the weights of the recurrent layer by means of backpropagation through time, as in more general RNNs11. Please, note that more sophisticated readout layers, requiring longer training procedures, could be used as well29. ESNs trade the precision of gradient descent with the “brute force” redundancy of the random recurrent layer. This, inevitably, makes ESNs more sensitive to selection of the hyperparameters that control their internal dynamics. Therefore, the difficulty involved in finding the optimal hyperparameters is more crucial for ESNs than in other RNN models. The non-linear recurrent nature of ESNs hampers interpretability of internal state dynamics30. As such, ESNs (as any other RNN) are usually treated as black box models with hyperparameters tuned in a supervised manner through cross-validation. However, cross-validation requires evaluation of the computational performance on a validation set for each hyperparameter configuration. This might be a serious limitation in real-life applications, whenever the data are scarce and/or supervised information unavailable. Moreover, the model has to be re-trained for each hyperparameter configuration; this is an issue whenever the learning procedure is computationally demanding29.

Several unsupervised methods for ESNs tuning have been proposed in the literature, e.g., see Refs.31,32. Although unsupervised tuning seldom achieves the same level of accuracy of the supervised counterpart, it offers some insights to the interpretation of the system behavior33. Typically, unsupervised learning is achieved by observing the dynamics of the states. These methods are usually based on a statistical framework, which requires reservoir outputs to be independent and identically distributed (i.i.d.). However, internal states are not identically distributed when ESNs are driven by non-stationary inputs34. Independence, instead, is always violated both by dependent inputs and by temporal dependencies introduced by recurrent connections. Therefore, even if the input signal is i.i.d., neurons activations may be identically distributed but never independent.

With this work, we propose two unsupervised graph-based methods for tuning ESN hyperparameters, with the aim to maximize (i) prediction accuracy and (ii) memory capacity. In our framework, temporal dependencies in time series of neurons activations (states) are converted into connections of a graph representation. This step allows to relax the i.i.d. assumption of statistical methods. Then, graphs associated with all neurons activations become layers of a structure called a multiplex. A multiplex graph35–39 is a special type of multilayer network, whose vertices are replicated through each layer and connected across layers only with their replicas. The topology (i.e., the edges) of the graph in each layer can be different. In particular, here we use a multiplex visibility graph40, whose layers consist of a particular class of graphs called horizontal visibility graphs (HVGs)41. HVGs are planar graphs built from real-valued time series, whose elements have one-to-one correspondence with the vertices. Important properties, linking the structure of HVGs with features (e.g., onset of chaotic behavior) of the dynamic system underlying the analyzed time series, have been recently studied42–46. This last aspect provides an important justification for using HVGs to model neuron activations for the purpose of hyperparameter tuning.

The novelty introduced in our paper refers to the characterization of ESN dynamics through structural characteristics of the associated multiplex network. To the best of our knowledge, this is the first time ESNs are studied within the framework of multiplex networks. However, it is worth citing a loosely related paper by Zhang et al.47, where the authors set the basis for a graph-theoretical analysis of RNNs. Here, we represent the instantaneous state of an ESN through a set of vertex properties of the multiplex (e.g., degree, clustering coefficient). Edges in HVGs might cover a relevant time interval (contained within the longest period of the signal), while they are still local in terms of topology. Accordingly, the ESN state can be characterized also by information that is non-local in time. To find hyperparameters yielding highest prediction accuracy, we

search hyperparameter configurations producing neuron activations as diverse as possible. This occurs when neurons dynamics is maximally heterogeneous (critical dynamics), a characteristic that, as we will show in the paper, is well-captured by the average entropy of vertex properties of the multiplex. Successively, to quantify the amount of memory in an ESN, we check the existence of neuron activations that are “similar” to different delayed versions of the input. In fact, memory in ESNs depends on the ability to reproduce past input sequences from information kept within some neuron activations. We describe dynamics of delayed inputs and of neuron activations through unsupervised graph-based measures. Then, by evaluating the agreement between such measurements, we quantify the memory capacity. We provide experimental evidence that our methods achieve performance comparable with supervised techniques for identifying hyperparamer configurations with high prediction accuracy and large memory capacity.

Methods

Echo state networks

A standard ESN architecture consists of a large, recurrent layer of non-linear neurons sparsely interconnected by edges with randomly generated weights, called a reservoir, and a linear, feedforward readout layer that is usually trained with regularized least-squares optimization27. De facto, the recurrent layer acts as a non-linear kernel48 that maps inputs to a high-dimensional space. The time-invariant difference equations describing the ESN state-update and output are, respectively, defined as

The reservoir contains Nr neurons whose transfer/activation function fis usually a hyperbolic tangent and gis usually the identity function. At time instant t, the ESN is driven by the input signal xand produces the output y, being Ni and No input and output dimensions, respectively. The vector hdescribes the ESN (instantaneous) state. The weight matrices Wrr (reservoir connections), Wri (input-to-reservoir connections), and Wro (output-to-reservoir connections) contain real values randomly sampled from a uniform or Gaussian distribution. Weights in Woi and Wor, instead, are trained (typically by means of regularized least-squares algorithm) according to the task at hand.

Once given, the designer can control Wrr, Wro, and Wri only through scaling coefficients. For instance Wro is scaled by a multiplicative constant , which has the effect to tune the impact of the ESN output on the next state. In this study, the output feedback is removed by setting 0. Hence, the resulting state-update difference equation (1) becomes

while the output equation (2) remains unchanged.

The properties of the ESN described in (3) heavily depend on two hyperparameters. First, input weights in Wri are scaled by a multiplicative constant , which controls the amount of non-linearity introduced by neurons. Large values for saturate the non-linear activation functions. The second hyperparameter is the spectral radius (eigenvalue with largest absolute value) of Wrr, which is related to the echo-state property, discussed earlier. For a detailed discussion on the relationships between stability, performance and , we suggest to the interested reader Ref.30 and references therein. Here, it suffices to say that a widely adopted rule-of-thumb49 suggests to set to a value slightly smaller than 1 (e.g., ). However, to reach higher performance in some practical tasks, it could be necessary to pick a small value for or to breach the aforementioned “safety” bound and push its value beyond unity. Note that in this latter case asymptotic stability of the ESN might still hold, even if some assumptions are locally violated.

In this work, we focus only on the tuning of the hyperparamers and , without affecting generality of the proposed methods that can be adopted to tune also other hyperparamers. The optimal values of for the task at hand are typically identified with a cross-validation procedure, with the potential associated shortcomings mentioned above. For that reason, different unsupervised approaches have been proposed for their tuning30,31,33. For example, the effect of on the ESN computational capability can be investigated through the maximal local Lyapunov exponent (MLLE), which measures the divergence rate in state space of trajectories with similar initial conditions. In autonomous (not input-driven) systems, chaos occurs when the maximal Lyapuanov exponent becomes positive, while in input-driven systems, like ESN, one typically relies on local first-order approximations of this quantity (see30 for details). Accordingly, the onset of criticality in ESNs can be detected by checking when MLLE crosses 0. Another quantity, which was shown to be more accurate in detecting criticality in dynamic systems and well-correlated with ESN performance, is the minimal singular value (in average, over time) of the reservoir Jacobian, denoted as is unimodal and in correspondence of its maximum the dynamical system is far from singularity, has many degrees of freedom, has a good excitability, and it separates well the input signals in state space33. By assuming a null input in Eq. 3, the Jacobian matrix of the reservoir at time t is given by Jt where the diagoperator returns a diagonal matrix. In this paper, will serve as a baseline for comparison to our proposed graph-based unsupervised methods for improved ESN hyperparameter tuning, as discussed in the next section.

The HVG41 associated with a finite univariate time series x , is constructed by assigning a vertex vt to each datum x[t]. The adjacency matrix A characterizes the graph: two vertices vi and vj, i are connected by an edge (A[i, j] = 1) iff the corresponding data fulfill the criterion x. In a multivariate scenario, the data stream is composed of Nr different time series , of equal length tmax. A multivariate time series can be mapped into a multiplex visibility graph . Specifically, the lth multiplex layer is defined by the HVG Gl constructed from xl. In the multiplex, a vertex is replicated on all layers and such replicas are linked by inter-layer connections, while intra-layer connections might change in each layer. From now on, we denote vl[t] to be the vertex of Gl in layer l associated with time interval t.

In this paper we introduce a weighted HVG (wHVG), with edge values defined as A. Since self-loops are forbidden in HVGs, edge weights are always well-defined (i.e., finite). The use of weights permits to capture additional information, as it accounts for distance in time and amplitude differences of two data points connected by the visibility rule. This weighting scheme is motivated by our need to characterize and exploit the instantaneous state by means of a suitable measure of heterogeneity (discussed in the next section). To distinguish the original adjacency matrix from the one of the wHVG, we refer the former as binary adjacency matrix and the latter as weighted adjacency matrix. Algorithm 1 delivers the pseudo-code for constructing a HVG (and a wHVG) from time series x. The worst case complexity of this algorithm is , which occurs when values in x are monotonically decreasing. Instead, the best case complexity is and arises in correspondence of monotonically increasing values in x.

Vertex properties

Let us consider an ESN with a reservoir of size Nr, driven by an input of length tmax. According to (3), the ESN generates Nr time series , of neuron activations (state). The multivariate time series is represented here by a multiplex M . Indexes of the vertices on each layer l of M are associated one-to-one with the time index of the original time series. Hence, the ESN state is represented by the vertices and by vertex properties we introduce four indexes: vertex degree, clustering coefficient50, betweenness and closeness centrality51.

Table 1. Vertex characteristics used to generate different instances of the vertex property vector

The capability of an ESN to reproduce the dynamics of a target system, hence to predict its trajectory in state space, is maximized on the edge of criticality, where the internal dynamical patterns of an ESN become sufficiently rich. In the literature,

such a “richness” is usually expressed in terms of diversity of connection weights53, entropy or rank of the matrix of neuron activations15,32.

In the same spirit, in order to find hyperparameters yielding maximal prediction accuracy, here we look for those hyperparameter configurations giving rise to neuron activations that are as heterogeneous as possible. Fig. 2 provides a visual example of the concept herein discussed. As an illustration, we consider an ESN driven by a sinusoid. We depict the neuron activations and the corresponding HVGs. We select a time step (marked by a black square in the picture) and we show the correspondence between the element in each time series of the neuron activations and the vertices in the associated HVGs. When the ESN

Figure 2. Heterogeneity of the graph-based representation of the ESN instantaneous state. From each time series, we visualize only a portion of its HVG and consider a specific time step (black square). In the HVGs, we highlight in red the edges connecting each vertex that corresponds to the selected time step. For the clustering coefficient, instead, one has to consider also edges highlighted in green. The non-local nature of ESN graph-based state is evident from such representations, since vertices can be connected to other vertices that are associated to time instants far away along the sequence.

operates with a contractive dynamic (2a), neuron activations weakly depend on previous ESN states and they are all very similar to the input signal. This results in a lack of diversity among activations. Accordingly, the corresponding HVGs in the multiplex contain vertices with similar properties, e.g., similar degree and clustering coefficient. When the ESN approaches the edge of criticality, neuron activities highly depend on previous internal states, which encode information of past inputs and internal structure of neuron connections. As we can see from (2b), on the edge of criticality, the activations have different frequencies and their phases are shifted. Such a heterogeneity in the ESN instantaneous state is captured by vertex properties of the multiplex. Finally, if pushed beyond the edge of criticality, the ESN transits into a chaotic regime as shown in (2c). In this case, neuron activations become noise-like and disordered oscillations generate HVGs with vertex properties very different from previous configurations. However, the diversity of the patterns in different time series disappears, hence heterogeneity is again lost. This lack of variety is also highlighted by recurrence of similar motifs in the corresponding HVGs.

To determine the heterogeneity of neuron activations, we consider the entropy of the related vertex property distribution. In particular, heterogeneity is computed as follows:

1. At each time t, evaluate the vertex properties using one of the properties listed in Tab. 1.

2. Estimate distribution pusing a histogram with b bins.

3. Compute instantaneous entropy: Ht is the Shannon entropy.

4. Define heterogeneity as average entropy over time instants: ¯H

Figure 3. High-level illustration of the proposed methodology. The activity of an ESN that is driven by some input signal is analyzed through a multivariate time series of neuron activations. We map such time series with a multiplex formed by HVGs on the layers. In this way, we get rid of problems related to stationarity and dependence of neuron activations. Note that with such a representation, the state of an ESN is defined through graph-theoretical tools, which allows to include non-local (in time) information in the ESN instantaneous state. For instance, at time instant t4, instead of the original activations, we consider a

vertex property (e.g., vertex degrees ) of the HVG to represent the ESN state. The instantaneous state is eventually

characterized by a measure of heterogeneity of such values; here, entropy is taken into account.

The procedure described above, is visually represented in Fig. 3. By referring to the figure, we observe that, at given time t, the ESN state is represented by those vertices in the multiplex tagged with the same label across the layers. For example, the ones in red describe the instantaneous state of the ESN at time step t4 and a vector of vertex properties is computed for this set of vertices. Then, we estimate the distribution pand compute the entropy H. In order to compute ¯H, the procedure is repeated for each time step. Finally, ¯His used to characterize the current hyperparameter configuration. The average entropy depends on the specific vertex property chosen for the analysis (we taken into account four properties, e.g., vertex degree, clustering coefficient, betweenness and closeness centrality – see Tab. 1). These vertex properties lead to four different entropy values ¯HDG, ¯HCL, ¯HBC, and ¯HCC.

To summarize, given an input signal, we select hyperparameter configurations that maximize such entropy values. This criterion is inspired by the aforementioned observation linking performance of a computational dynamical system (i.e., prediction accuracy and memory) with heterogeneity of its dynamics (critical dynamics). Accordingly, in this case this information is exploited in order to derive, in an unsupervised way, the configuration yielding highest prediction accuracy.

Other multiplex complexity measures

Recently40, two measures have been proposed in order to characterize the dynamics of a system observed through a multivariate time series and represented as a multiplex composed of HVGs. Here, we consider these measures in order to evaluate whether they are useful for identifying the hyperparameter configurations yielding maximum accuracy performances or not.

The Average Edge Overlap (AEO) computes the expected number of layers of the multiplex on which an edge is present. For binary HVGs, it is defined as

where Al[i, j] equals 1 if vertices i and j are connected in layer l. The second measure is the Interlayer Mutual Information (IMI), which quantifies the correlations between the degree distributions of two different layers li and lj. It is defined, for binary HVGs, as

where pis the joint probability to find a vertex with degree kli in layer li and degree klj in layer l j, respectively. In the experiments, we use the average IMI between all pairs of multiplex layers.

Memory measures

The memory of an ESN is quantified by the capability to retain information about past inputs in its transient dynamics54. A supervised measure called memory capacity (MC) is usually adopted to quantify ESN memory27. When computing MC, the reservoir topology and weights are kept fixed and several readout layers are trained in order to reproduce delayed versions of the input at different time lags . The MC is computed as the squared correlation coefficient between desired (delayed input) and computed outputs,

In order to properly evaluate the capability of ESNs to introduce memory through recurrent connections, x[t] is defined as a stationary uncorrelated noisy signal.

In the following, we propose an unsupervised graph-based method to identify hyperparameter configurations for which an ESN achieves large memory capacity. Given the input time series x , we determine if there exists a subset of neuron activations, which is correlated with a past input sequence xaBeing Gl the HVG representing the lth layer of the multiplex, is the sequence of its vertex degrees ordered according to the time index (not to be confused with , the degrees of vertices relative to time t across the different layers). With Gx, instead, we refer to the HVG constructed over the input xis the vector of its vertex degrees.

First, we define a measure of maximum agreement between and each sequence is a similarity measure between sequences: in this paper we consider the Pearson correlation , the Spearman correlation , and the mutual information

A second measure of agreement is defined on the adjacency matrix Ql , where Qland Al[i, j] = 1; Ql[i, j] = 0 otherwise. The agreement is then computed as the largest number of edges among all possible intersection graphs, counts the number of ones in the adjacency matrix.Finally, one might consider the similarity between original input x and neuron activations hl, In this case, is directly evaluated on the time series rather than on the sequences of vertex degrees. This last measure of agreement is taken into account in order to quantitatively show (in the experiments) the benefits of using HVGs for representing neuron activations.

By referring to the illustrative example in Fig. 4, we generate the HVGs Gx, G1, relative to delayed input x15and activations of the Nr reservoir neurons, respectively. Successively, we evaluate their similarities by means of an agreement measure (in this example). The similarity measure is chosen among the three previously proposed measures: or

Results

In the following, we perform two experiments in order to evaluate the two proposed unsupervised methods for, respectively, finding hyperparameter configurations giving rise to ESNs with high prediction accuracy and large memory capacity. In the first experiment, we show that, on different real and synthetic tasks, (supervised) prediction accuracy is maximized for the same hyperparameter configurations that yields the largest heterogeneity for the vertex properties of the multiplex. In the second experiment, we show the reliability of the graph-based memory measures in identifying hyperparameters where the (supervised) MC is maximized.

Figure 4. We extract a subsequence from the history of the input signal, for example x15, and compute the related HVG Gx. We repeat the procedure for each neuron activation, ending up with G1. For each HVG, we compute the sequence of vertex degrees and evaluate the similarities . The largest similarity value determines the value of . Here, the activation of neuron n1 is clearly the most related one with the input sequence x15and therefore it

determines the value of

Test for prediction accuracy

In this experiment, we consider several prediction tasks and, for each of them, we set the forecast step to be the smallest time-lag that guarantees the measurements in a time window of size to be uncorrelated (e.g., the first zero in the autocorrelation function of the input signal). Prediction error is evaluated by Normalized Root Mean Squared Error, NRMSE is the prediction provided by the ESN and y[t] is the desired/teacher output. The Prediction accuracy is defined as In the following, we describe the datasets used in this experimental campaign. Sinusoidal input: we feed an ESN with a sinusoid yand we predict future input values with a forecast step

Mackey-Glass time series: the Mackey-Glass (MG) system is commonly used as a benchmark in chaotic time series

prediction. The input signal is generated from the MG time-delay differential equation

We adopt the standard parameters , initial condition x(0) = 1.2, and integration step equal to 0.1. The forecast step here is

Multiple superimposed oscillator: prediction of superimposed sinusoidal waves with incommensurable frequencies is a hard forecasting exercise, due to the extension of the wavelength55. The ESN is fed with the multiple superimposed oscillator (MSO)

For this task, the ESN is trained to predict future input values, with forecast horizon NARMA: the chosen Non-Linear Auto-Regressive Moving Average (NARMA) task56 consists in modeling the output of the r-order system:

x[t] is a uniform random noise signal in [0,1] and is the input of the ESN, which is trained to reproduce y[t +1]. The NARMA task is known to require a memory of at least r past time-steps, since the output is determined by inputs and outputs from the last r time-steps. For this task we set r

Polynomial task: the ESN is fed with uniform noise in and is trained to reproduce the following output

where ciis uniformly distributed in . The difficulty of prediction can be controlled by varying the polynomial degree p and the time delay d. For this task we set p

Telephone call load time series: as a last test, we consider a real-world dataset relative to the load of phone calls registered over a mobile network. The data comes from the Orange telephone dataset, published in the Data for Development (D4D) challenge58. D4D is a collection of call data records, containing anonymized events of Orange’s mobile phone users in the Ivory Coast, in a period spanning from December 1, 2011 to April 28, 2012. The dataset consists of 6 time series consisting of: number and volume of incoming calls, number and volume of outgoing calls, day and time (1 hour resolution) when the telephone activity was registered. More detailed information is available at the related website59. All 6 time series are fed into the ESN as inputs; the goal is to predict 6 hours ahead the volume of incoming calls – the profile of this latter time series is depicted in Fig. 5.

Figure 5. D4D dataset – load profile of incoming calls volume for the first 300 time intervals.

In each test, we evaluate the correlation between the average entropy of vertex properties in the multiplex and the prediction accuracy , as we vary the hyperparameters . Multiplexes are generated using both the binary and weighted version of the HVG adjacency matrix. To appreciate the effectiveness of our methodology, we also estimate the correlations of with , the minimum singular value of the Jacobian of the reservoir (see Methods). Additionally, we consider the correlation of with the two layer-based measures IMI and AEO (see Methods). Correlations are evaluated as follows. For each configuration , we have the prediction accuracy , the entropy ¯HDGk, and so on. The values assumed by these quantities by varying and generate a two-dimensional manifold. The point-wise linear (Pearson) correlation between the manifolds relative to and the other considered measures is the result we are interested in. The configurations of the hyperparameters that are examined are generated by varying in [0.5,1.3] (20 different values) and in [0.2,0.9] (10 different values). A total of 200 configurations are evaluated. Due to the stochastic nature of the ESN initialization, for each configuration we compute and all the other measures 15 different times. Successively, we compute the correlations among their average values. We use a reservoir with Nr = 100 neurons and sparsity of the internal connectivity equal to 25%. The readout is trained by standard ridge regression with regularization parameter set to 0.05. The distributions of the vertex properties are estimated by histograms with b = 50 bins. We comment that the reservoir size Nr can be increased/decreased without affecting the applicability of the proposed methodology (only the number of bins used to estimate the vertex properties distribution might be modified – see Methods). However, the number of neurons has an impact on the overall performance of the network and on the time and space complexities of the proposed method.

Fig. 6 depicts the values assumed by and four graph-based measures in the case of the MSO prediction task. As can be seen, high correlation emerges between the (average) entropy of the vertex properties and the prediction accuracy. Since our approach is fully unsupervised, the proposed graph-based measures can approximate well the accuracy , regardless of the task learned by the readout (prediction, function approximation, reconstruction of past inputs, etc). In Tab. 2, we report the average correlation values and their statistical significance (expressed by p-values) on all tasks. As we can see, the highest (and statistically significant) correlation is achieved by using one of the four average entropy measures of vertex properties. In particular, the measure based on vertex cluster coefficient distribution, ¯HCL, achieves the best results in 5 of the 6 tasks. For what concerns the D4D time series, we observe that achieves high correlation with , but still lower than the one achieved by ¯HCL. This demonstrates the effectiveness of the proposed methodology, also in the case of a real-world application. In SIN and MSO tasks, the graph-based quantifiers estimated on the weighted HVG achieve a higher degree of accuracy, with respect to the binary counterpart. In these cases, additional qualitative information relative to temporal and amplitude differences in the connected data allows to better represent the dynamics of the system. Finally, it is worth noting that IMI takes high, yet negative correlation values on both MG and POLY tasks. In such cases, results are close to the ones achieved with our approach.

The performed experiment consists in generating 100 different random reservoirs, each one characterized by an increasing value of spectral radius in the [0.1,2] interval. As varies, we evaluate the MC by training four readouts in order to reproduce different time-lagged versions of input signal x10. Then, on the output of each reservoir, we evaluate the similarities

Figure 6. In each panel, the two-dimensional gray manifold represents the values of the prediction accuracy in the MSO task, for different configurations of . The colored manifolds are the average entropy of the distributions of (a) the vertex degree ¯HDG, (b) the clustering coefficient ¯HCL, (c) the betweenness centrality ¯HCC. The manifold in (d) represents which is computed from the ESN Jacobian. From the figure, we observe strong correlations between supervised prediction accuracy and the unsupervised measures (a–c) derived from the multiplex representation of the ESN activity.

Table 2. Correlations and related p-values (in brackets) of the proposed graph-based measures with accuracy in different prediction tasks, as change. We report the correlation of the manifolds generated by the values of with the manifolds relative to ¯HDG, ¯HCL, ¯HBC, and ¯HCC, which are the average entropy values of the distributions of vertex degree, clustering coefficient, betweenness and closeness centrality. Each measure is computed on both the binary (b) and weighted (w) versions of the HVG adjacency matrix (adj). We also report the correlations of with the manifolds relative to the minimum singular value of the reservoir Jacobian over time () and the two multiplex-based measures AEO and IMI, presented in40. In each task, the highest correlations with are highlighted in bold.

, which are high if there exists at least one series of activations that is similar to the considered past input sequence. This is evaluated in such as way that the measure taken into account. Even if some neurons retain dynamics of previous input sequences, the reservoir introduces shifts in the phase and the amplitude of the input signal. To filter out these differences, in this test we consider only HVGs defined by binary adjacency matrices, which do not account for differences in the amplitude of the connected values. To evaluate the effectiveness of the proposed unsupervised memory measures, we compute the correlation between the supervised MC and varies within the chosen interval. Note that we only monitor the effect of on the dynamics, since it is the hyperparameter that mostly affects the memory capacity32. We kept the input scaling fixed, 7, while the remaining hyperparamers are configured as in the previous experiment. As before, we repeated each experiment 15 times with different and independent random initializations. In Tab. 3, we show the mean correlation values, along with the standard deviations, between the MC and the proposed unsupervised measures of memory capacity.

From the table, we observe that the best agreement with the MC is achieved by the measures derived from the HVGs. In

Table 3. Mean correlations and standard deviations of MC with the unsupervised memory quantifiers We consider 4 different sequences of past inputs x10. Values for each measure are computed as the spectral radius of the ESNs reservoirs varies from 0.1 to 2. To compute , we consider three different similarities: the Pearson correlation , the Spearman correlation , and the mutual information . Best results for each input sequence are reported in bold.

particular, configured with the Spearman rank is always highly correlated with MC and, in three of the four delayed input sequences taken into account, is the best performing one. In each setup, works better if configured with than . Instead, results obtained with are significantly worse in all cases. achieves the best results only for the first time lag taken into account, while the agreement with MC is lower in the remaining cases. Interestingly, several measures show a high degree of correlation with the MC as the size of the delay increases. , the unsupervised measure computed directly on the input time series and neuron activations, shows positive correlations with the MC, but the agreement is always lower with respect to the graph-based measures. For , the setting with works better than . Finally, also in this case by using we obtain the worst performance. In Fig. 7, we show an example of the values of MC, (configured with ), (configured with is varied within the [0.2,2] interval.

Figure 7. Values of supervised MC and of two selected unsupervised memory measures, when the input sequence x20taken into account. The smaller the residuals (red areas), the better. We can observe that the agreement of the MC with (using ) is higher than

Discussion

Experimental results show satisfactory correlations for average entropy of vertex properties with respect to prediction accuracy. Moreover, the two unsupervised graph-based memory measures that we proposed (and ) correlate well with the supervised measure of memory capacity.

We first discuss the results of the prediction accuracy test, where we analyzed topological properties of vertices in the multiplex, representing the ESN instantaneous state. On all tests taken into account, we observed a remarkable correlation between and the average entropy of the clustering coefficient distribution ¯HCL, hence suggesting that the clustering coefficient is able to describe well the heterogeneity of the activations. To explain this result, it is necessary to elaborate on the properties of the clustering coefficient CL. In the HVG literature, CLits behavior has been analyzed for time series characterized by different Hurst exponent60. Additionally, an upper bound (CL) is provided for HVGs derived from random time series41. In the following, we present an in-depth interpretation of the results by accounting for geometrical properties of the clustering coefficient.

In a HVG, CL(v) measures the inter-visibility among neighbors of v. For convex functions, it is possible to connect any two points with a straight line. This feature is also (partially) captured by the HVG. If v is contained in a convex part of the related

time series, there is a high degree of intervisibility among the neighbor vertices to which v is connected, hence CL(v) is high. Additionally, moving along the same convex part of the time series results only in minor changes of the clustering coefficient in the associated HVG vertices. Instead, if v is a local maximum of a concave part, then it is connected to points belonging to two different basins, which do not have reciprocal visibility. In this case, CL(v) is low and its value rapidly changes as one moves away from the maximum. This results in great losses of visual information. Therefore, large values of CLindicate the presence of dominating convexities, while low values characterize concavities61. Accordingly, CLcan be used to measure the length of a convex (concave) part of the time series and how fast the convexity is changing, which is a measure of the fluctuations in the time series62. In a regime characterized by contractive dynamics, convexity changes at the same (slow) rate in different neuron activations and this results in a low entropy value of the clustering coefficient distribution among vertices in different layers. On the edge of criticality, instead, convex and concave parts in the time series of activations are characterized by heterogeneous lengths and they change at different rates. This corresponds to a high degree of clustering coefficient diversity of the same HVG vertex, replicated at different layers in the multiplex. Finally, in the chaotic regime, all time series fluctuate very rapidly and their convexity changes every few time steps. In this case, in each time series of activations the lengths of convex and concave parts are always very short and hence the desired heterogeneity is again lost.

For what concerns experiments on memory, the best overall results in terms of agreement with the supervised MC are achieved by the graph-based measure . As previously discussed, such a measure evaluates the maximum similarity between the sequence of vertex degrees on the input HVG Gx and the HVG Gl of neuron activations. This measure is closely related with the degree distribution P(k), whose importance is known in the HVG literature41. For example, it has been shown that for time series generated from an i.i.d. process, P(k) follows Pand the mean degree is . As the correlations in the time series increase, the i.i.d. assumption is lost and P(k) decays faster. Furthermore, vertex degrees are key parameters to describe dynamic processes on the graph, such as synchronization of coupled oscillators, percolation, epidemic spreading, and linear stability of equilibrium in networked coupled systems63. Their role has been studied also in the HVG framework64. HVGs have been studied in the context of time series related to processes with power-law correlations60. In our case, the time series of neuron activations have short-term correlations and increments in the correlation coefficients can have opposite signs at consecutive time lags. For these cases, we are not aware of any previous study in terms of HVGs.

In networks which are inherently degree disassortative, the range of degree values increases with network size, with a consequent decrease of the assortativity value65. In such networks, the Spearman rank correlation provides a more suitable choice with respect to calculating degree-degree Pearson correlations. It is important to notice that the rank is computed through a non-linear rescaling, which is data dependent. The information on the actual values of the data is discarded as only its inherent ordering (rank) is preserved. We argue that HVGs convey the same type of information captured by the Spearman correlation. Hence, the latter should be preferred to Pearson correlation to characterize the characteristics of the vertices and related topological properties in HVGs. This fact justifies the higher agreement with memory capacity achieved by means of when configured with , which accounts for Spearman correlations between sequences of vertex degrees in the HVGs related to the input signal and the neuron activations.

Modeling ESN dynamics through a multiplex network allowed us to connect two seemingly different research fields, thus fostering multidisciplinary research in the context of recurrent neural networks. By converting a temporal problem into a topological one, we handled temporal dependencies introduced by ESNs (as well as by other types of RNNs), hence overcoming technical limitations of statistical approaches that require independence of samples. We performed and discussed several experiments that provided empirical evidence that our methodology achieves performance higher than other unsupervised methods and comparable to cross-validation techniques. These results suggest to allocate efforts to further improve the effectiveness of unsupervised learning methods in the context of ESNs and RNNs. Finally, we would like to stress that, while this paper is primarily focused on network structures in machine learning, our results might suggest new ideas for theoretical understanding of recurrent structures in biological models of neuronal networks7,66,67. In particular, we believe that it is possible to identify emergent structural patterns in the developed graph-based representations of network dynamics. This would allow to further explore and analyze the route to chaos in input-driven neural models by exploiting the language of graph theory, which is an already established framework within the neuroscience field68.

Additional information

Author contributions

LL outlined the research ideas. FMB and LL conceived methods. FMB performed experiments and generated figures. CA and

RJ contributed to the technical discussion. All authors took part to the paper writing and approved the final manuscript.

References

1. Barzel, B. & Barab´asi, A-L. Universality in network dynamics. Nature Physics 9, 673–681 (2013).

2. Hammer, B., Micheli, A., Sperduti, A. & Strickert, M. Recursive self-organizing network models. Neural Networks 17, 1061–1085 (2004).

3. Maass, W., Joshi, P. & Sontag, E. D. Computational aspects of feedback in neural circuits. PLoS Computational Biology 3, e165 (2007).

4. Reinhart, R. F. & Steil, J. J. Regularization and stability in reservoir networks with output feedback. Neurocomputing 90, 96–105 (2012).

5. Jordan, M. I. Serial order: A parallel distributed processing approach. In Donahoe, J. W. & Dorsel, V. P. (eds.) NeuralNetwork Models of Cognition: Biobehavioral Foundations, vol. 121 of Advances in Psychology, 471 – 495 (North-Holland, 1997).

6. Amit, D. J., Gutfreund, H. & Sompolinsky, H. Spin-glass models of neural networks. Physical Review A 32, 1007–1018 (1985).

7. Enel, P., Procyk, E., Quilodran, R. & Dominey, P. F. Reservoir computing properties of neural dynamics in prefrontal cortex. PLoS Computational Biology 12, e1004967 (2016).

8. Barak, O., Sussillo, D., Romo, R., Tsodyks, M. & Abbott, L. F. From fixed points to chaos: three models of delayed discrimination. Progress in Neurobiology 103, 214–222 (2013).

9. Rajan, K., Abbott, L. F. & Sompolinsky, H. Stimulus-dependent suppression of chaos in recurrent neural networks. Physical Review E 82, 011903 (2010).

10. Fusi, S., Miller, E. K. & Rigotti, M. Why neurons mix: high dimensionality for higher cognition. Current Opinion in Neurobiology 37, 66–74 (2016).

11. Schmidhuber, J. Deep learning in neural networks: An overview. Neural Networks 61, 85 – 117 (2015).

12. Barra, A., Bernacchia, A., Santucci, E. & Contucci, P. On the equivalence of Hopfield networks and Boltzmann machines. Neural Networks 34, 1–9 (2012).

13. Charles, A., Yin, D. & Rozell, C. Distributed sequence memory of multidimensional inputs in recurrent networks. arXiv preprint arXiv:1605.08346 (2016).

14. Tiˇno, P. & Rodan, A. Short term memory in input-driven linear dynamical systems. Neurocomputing 112, 58–63 (2013).

15. Legenstein, R. & Maass, W. What makes a dynamical system computationally powerful? In Haykin, S., Principe, J., Sejnowski & McWhirter (eds.) New Directions in Statistical Signal Processing: From Systems to Brain, 127–154 (MIT Press, Cambridge, 2007).

16. Yildiz, I. B., Jaeger, H. & Kiebel, S. J. Re-visiting the echo state property. Neural Networks 35, 1–9 (2012).

17. Dambre, J., Verstraeten, D., Schrauwen, B. & Massar, S. Information processing capacity of dynamical systems. Scientific Reports 2 (2012).

18. Lajoie, G., Lin, K. K., Thivierge, J.-P. & Shea-Brown, E. Encoding in balanced networks: Revisiting spike patterns and chaos in stimulus-driven systems. PLoS Computational Biology 12, 1–30 (2016).

19. Livi, L., Bianchi, F. M. & Alippi, C. Determination of the edge of criticality in echo state networks through Fisher information maximization. IEEE Transactions on Neural Networks and Learning Systems 1–12 (2017). doi:

10.1109/TNNLS.2016.2644268.

20. Mayer, N. M. Input-anticipating critical reservoirs show power law forgetting of unexpected input events. Neural Computation 27, 1102–1119 (2015).

21. Moretti, P. & Mu˜noz, M. A. Griffiths phases and the stretching of criticality in brain networks. Nature Communications 4 (2013).

22. Mora, T. & Bialek, W. Are biological systems poised at criticality? Journal of Statistical Physics 144, 268–302 (2011).

23. Scheffer, M. et al. Anticipating critical transitions. Science 338, 344–348 (2012).

24. Wang, X., Lizier, J. & Prokopenko, M. Fisher information at the edge of chaos in random Boolean networks. Artificial Life 17, 315–329 (2011).

25. Tkaˇcik, G. et al. Thermodynamics and signatures of criticality in a network of neurons. Proceedings of the National Academy of Sciences 112, 11508–11513 (2015).

26. Markovi´c, D. & Gros, C. Power laws and self-organized criticality in theory and nature. Physics Reports 536, 41–74 (2014).

27. Lukoˇseviˇcius, M. & Jaeger, H. Reservoir computing approaches to recurrent neural network training. Computer Science Review 3, 127–149 (2009).

28. Grigoryeva, L., Henriques, J., Larger, L. & Ortega, J-P. Optimal nonlinear information processing capacity in delay-based reservoir computers. Scientific Reports 5 (2015).

29. Bianchi, F. M., Scardapane, S., Uncini, A., Rizzi, A. & Sadeghian, A. Prediction of telephone calls load using echo state network with exogenous variables. Neural Networks 71, 204–213 (2015).

30. Bianchi, F. M., Livi, L. & Alippi, C. Investigating echo state networks dynamics by means of recurrence analysis. IEEE Transactions on Neural Networks and Learning Systems 1–13 (2016). doi: 10.1109/TNNLS.2016.2630802.

31. Boedecker, J., Obst, O., Lizier, J. T., Mayer, N. M. & Asada, M. Information processing in echo state networks at the edge of chaos. Theory in Biosciences 131, 205–213 (2012).

32. Ozturk, M. C., Xu, D. & Pr´ıncipe, J. C. Analysis and design of echo state networks. Neural Computation 19, 111–138 (2007).

33. Verstraeten, D. & Schrauwen, B. On the quantification of dynamics in reservoir computing. In Artificial Neural Networks– ICANN 2009, 985–994 (Springer Berlin Heidelberg, 2009).

34. Leisch, F., Trapletti, A. & Hornik, K. Stationarity and stability of autoregressive neural network processes. In Kearns, M. J., Solla, S. A. & Cohn, D. A. (eds.) Advances in Neural Information Processing Systems 11, 267–273 (MIT Press, 1999).

35. Lee, K.-M., Min, B. & Goh, K.-I. Towards real-world complexity: an introduction to multiplex networks. The European Physical Journal B 88, 1–20 (2015).

36. Boccaletti, S. et al. The structure and dynamics of multilayer networks. Physics Reports 544, 1–122 (2014).

37. Menichetti, G., Remondini, D., Panzarasa, P., Mondrag´on, R. J. & Bianconi, G. Weighted multiplex networks. PloS ONE 9, e97857 (2014).

38. Kivel¨a, M. et al. Multilayer networks. Journal of Complex Networks 2, 203–271 (2014).

39. De Domenico, M. et al. Mathematical formulation of multilayer networks. Physical Review X 3, 041022 (2013).

40. Lacasa, L., Nicosia, V. & Latora, V. Network structure of multivariate time series. Scientific Reports 5 (2015).

41. Luque, B., Lacasa, L., Ballesteros, F. & Luque, J. Horizontal visibility graphs: Exact results for random time series. Physical Review E 80, 046103 (2009).

42. Luque, B., Lacasa, L., Ballesteros, F. J. & Robledo, A. Feigenbaum graphs: A complex network perspective of chaos. PLoS One 6, e22411 (2011).

43. Luque, B., Cordero-Gracia, M., G´omez, M. & Robledo, A. Quasiperiodic graphs at the onset of chaos. Physical Review E 88, 062918 (2013).

44. Lacasa, L. & Toral, R. Description of stochastic and chaotic series using visibility graphs. Physical Review E 82, 036120 (2010).

45. Ravetti, M. G., Carpi, L. C., Gonc¸alves, B. A., Frery, A. C. & Rosso, O. A. Distinguishing noise from chaos: Objective versus subjective criteria using horizontal visibility graph. PloS ONE 9, e108004 (2014).

46. Luque, B., Lacasa, L., Ballesteros, F. J. & Robledo, A. Analytical properties of horizontal visibility graphs in the Feigenbaum scenario. Chaos: An Interdisciplinary Journal of Nonlinear Science 22, 013109 (2012).

47. Zhang, S. et al. Architectural complexity measures of recurrent neural networks. arXiv preprint arXiv:1602.08210 (2016).

48. Hermans, M. & Schrauwen, B. Recurrent kernel machines: Computing with infinite echo state networks. Neural Computation 24, 104–133 (2012).

49. Lukoˇseviˇcius, M. A Practical Guide to Applying Echo State Networks, 659–686 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012).

50. Saram¨aki, J., Kivel¨a, M., Onnela, J.-P., Kaski, K. & Kertesz, J. Generalizations of the clustering coefficient to weighted complex networks. Physical Review E 75, 027105 (2007).

51. Boccaletti, S., Latora, V., Moreno, Y., Chavez, M. & Hwang, D. Complex networks: Structure and dynamics. Physics Reports 424, 175–308 (2006).

52. Lopez-Fernandez, L., Robles, G. & Gonzalez-Barahona, J. M. Applying social network analysis to the information in CVS repositories. In Proceedings of the International Workshop on Mining Software Repositories, 101–105 (Edinburgh, UK, 2004).

53. Bertschinger, N. & Natschl¨ager, T. Real-time computation at the edge of chaos in recurrent neural networks. Neural Computation 16, 1413–1436 (2004).

54. Pascanu, R. & Jaeger, H. A neurodynamical model for working memory. Neural Networks 24, 199–207 (2011).

55. Jaeger, H. & Haas, H. Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication. science 304, 78–80 (2004).

56. Jaeger, H. Adaptive nonlinear system identification with echo state networks. In Advances in Neural Information Processing Systems, 593–600 (MIT Press, 2002).

57. Butcher, J., Verstraeten, D., Schrauwen, B., Day, C. & Haycock, P. Reservoir computing and extreme learning machines for non-linear time-series data analysis. Neural Networks 38, 76 – 89 (2013).

58. Blondel, V. D. et al. Data for Development: the D4D Challenge on Mobile Phone Data. ArXiv preprint arXiv:1210.0137 (2012).

59. Orange d4d challenge. http://http://www.d4d.orange.com/en/Accueil. Accessed: 2016-09-22.

60. Xie, W.-J. & Zhou, W.-X. Horizontal visibility graphs transformed from fractional Brownian motions: Topological properties versus the Hurst index. Physica A: Statistical Mechanics and its Applications 390, 3592–3601 (2011).

61. Costa, L. d. F., Rodrigues, F. A., Travieso, G. & Villas Boas, P. R. Characterization of complex networks: A survey of measurements. Advances in Physics 56, 167–242 (2007).

62. Turner, A., Doxa, M., O’Sullivan, D. & Penn, A. From isovists to visibility graphs: A methodology for the analysis of architectural space. Environment and Planning B: Planning and Design 28, 103–121 (2001).

63. Restrepo, J. G., Ott, E. & Hunt, B. R. Approximating the largest eigenvalue of network adjacency matrices. Physical Review E 76, 056119 (2007).

64. Fioriti, V., Tofani, A. & Di Pietro, A. Discriminating chaotic time series with visibility graph eigenvalues. Complex Systems 21 (2012).

65. Newman, M. E. J. Assortative mixing in networks. Physical Review Letters 89, 208701 (2002).

66. Marblestone, A. H., Wayne, G. & Kording, K. P. Toward an integration of deep learning and neuroscience. Frontiers in Computational Neuroscience 10 (2016).

67. Sussillo, D., Churchland, M. M., Kaufman, M. T. & Shenoy, K. V. A neural network that finds a naturalistic solution for the production of muscle activity. Nature Neuroscience 18, 1025–1033 (2015).

68. Sporns, O. Networks of the Brain (MIT press, Cambridge, MA, USA, 2011).

designed for accessibility and to further open science