Interpreting recurrent neural networks behaviour via excitable network attractors
Andrea Ceni Peter Ashwin Lorenzo Livi


Spectrum Matching in Licensed Spectrum Sharing
    M. Majid Butt, Senior Member, IEEE, Irene Macaluso, Eduard A. Jorswieck, Senior Member, IEEE, 
Julie Bradford, Nicola Marchetti, Senior Member, IEEE, and Linda Doyle


        
Irene Macaluso, Nicola Marchetti and Linda Doyle are with CONNECT center for future networks, Trinity College Dublin, Ireland. Email:{macalusoi, nicola.marchetti, linda.doyle}@tcd.ie.
M. Majid Butt is with school of engineering, University of Glasgow, UK. Email: majid.butt@glasgow.ac.uk.
Eduard A. Jorswieck is with Department of Electrical Engineering and Information Technology, TU Dresden, Germany. Email: eduard.jorswieck@tu-dresden.de.
Julie Bradford is with Real Wireless, UK. Email: julie.bradford@realwireless.biz
The project ADEL acknowledges the financial support of the Seventh Framework Programme for Research of the
European Commission under grant number: 619647. We also acknowledge support from the Science Foundation Ireland
under grants No. 13/RC/2077 and No. 10/CE/i853.












    December 30, 2023
======================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




Introduction: Machine learning provides fundamental tools both for scientific research and for the development of technologies with significant impact on society.
It provides methods that facilitate the discovery of regularities in data and that give predictions without explicit knowledge of the rules governing a system.
However, a price is paid for exploiting such flexibility: machine learning methods are typically black-boxes where it is difficult to fully understand what the machine is doing or how it is operating. This poses constraints on the applicability and explainability of such methods.
Methods: Our research aims to open the black-box of recurrent neural networks, an important family of neural networks used for processing sequential data. We propose a novel methodology that provides a mechanistic interpretation of behaviour when solving a computational task.
Our methodology uses mathematical constructs called excitable network attractors, which are invariant sets in phase space composed of stable attractors and excitable connections between them.
Results and Discussion: As the behaviour of recurrent neural networks depends both on training and on inputs to the system, we introduce an algorithm to extract network attractors directly from the trajectory of a neural network while solving tasks.
Simulations conducted on a controlled benchmark task confirm the relevance of these attractors for interpreting the behaviour of recurrent neural networks, at least for tasks that involve learning a finite number of stable states and transitions between them.







§ INTRODUCTION


Artificial recurrent neural networks (RNNs) are widely used to solve tasks involving temporal data, e.g., speech <cit.> and handwriting recognition <cit.>, audio classification <cit.> or time series forecasting <cit.>.
RNNs are characterised by the presence of feedback connections in a hidden layer, which allows generating a state-space representation that equips the network with short-term memory capability.
RNNs are universal approximators of dynamical systems <cit.> meaning that, given enough neurons in the hidden layer, it is possible to fine tune the weights to achieve any desired level of accuracy.
Nevertheless, training via back-propagation through time is difficult due to the vanishing/exploding gradient problem <cit.>. This has led to the development of new and faster techniques for training RNNs, including a different paradigm known as reservoir computing <cit.>.
Echo state networks (ESNs) <cit.> constitute an important example of reservoir computing, where a recurrent layer (called a reservoir) is composed of a large number of neurons with randomly initialised connections that are not fine-tuned via gradient-based optimisation mechanisms.
The main idea behind ESNs is to exploit the rich dynamics generated by the reservoir with an output layer, the read-out that is optimised to solve a specific task.




 §.§ Problem statement and research hypothesis


The high-dimensional and non-linear nature of RNNs complicates interpretability of their internal dynamics, which are characterised by complex, input-dependent spatio-temporal patterns of activity <cit.>.
This poses constraints on understanding the behaviour of RNNs: they are usually viewed as black-boxes from which it is hard to extract useful knowledge about their inner workings. As highlighted by recent research efforts <cit.>, similar interpretability issues affect many other machine learning methods. Furthermore, an increasing societal need to develop accountability and explainability of decision making by AI <cit.> is driving the development of methodologies for explaining the behaviour of such methods.

Our aim in this paper is to develop effective models that capture the essential dynamical behaviour of RNNs on computational tasks as input-driven responses of a dynamical system, while neglecting microscopic details of the RNN dynamics in phase space (i.e., the space of all possible neuron activations).
To this end, we hypothesise that RNNs can undertake computations by exploiting (i) transient dynamical regimes and (ii) excitable connections to switch between different stable attractors, depending on input and current state.
The RNN behaviour depends on the task at hand: for example, long transients are mostly exploited in time series forecasting problems, while switching between attractors is mostly exploited in classification problems or tasks requiring to learn a finite set of memory states. These mechanisms are not mutually exclusive and can be exploited synergistically to realise more complex computations.




 §.§ Contribution and paper organisation


In order to test our hypothesis, we develop a theoretical framework based on network attractors of dynamical systems.
An attractor is a subset of the state space of a dynamical system (i.e., the phase space), where the state will asymptotically converge for “usual” choices of initial conditions.
Network attractors are special kinds of attractor which can be thought of as directed graphs, where nodes correspond to local attractors and directed edges corresponds to particular trajectories that connect (or almost connect) those local attractors. 
Since for this paper, the local attractors of interest are all stable fixed points, we will use the term local attractor synonymously with stable fixed point. A stable fixed point is the simplest possible attractors: this is a point in phase space where all variables of the dynamical system assume constant values, and all close enough initial conditions converge to this point.
However, fixed points need not be stable: they can be partially stable (saddle points), or totally unstable (repellers).
According to the nature of the trajectories connecting local attractors it is possible to consider both heteroclinic network attractors composed of heteroclinic connections between saddles (i.e. partially stable fixed points), and excitable network attractors (ENAs) composed by connections that are excitable, i.e. that require a small initial perturbation to move between stable attractors <cit.>.

Heteroclinic network attractors have already been suggested as models to describe transient computation <cit.>, here, conversely, we focus attention on ENAs which have the advantage of being more robust to perturbations.
We focus particularly on ENAs between stable fixed points <cit.>, although ENAs can be theoretically conceived between any kind of local attractor: limit cycles, limit tori or strange attractors.
More precisely, we show how to construct ENAs describing RNN behaviour for tasks that involve switching between a finite set of attractors.
Interestingly, from a directed graph (representing the underlying network attractor of a dynamical system) it is possible to reverse engineer <cit.> a set of ordinary differential equations ruling, in our case, the RNN behaviour.
This, in turn, opens the way to a more analytic description of how RNNs solve computational tasks.

In this paper, we focus on the flip-flop task with two bits <cit.>, since it provides a controlled workbench to test our hypothesis. The input to discrete-time RNN is assumed to be a discrete temporal sequence with values from {+1,0,-1}.
In general, we consider discrete time k varying input, denoted as u[k], and output, denoted as z[k], of a system related by

    ( x[k] ,  z[k] ) = R ( u[k], x[k-1]),

where R represents an RNN with evolving internal state, namely x[k], and connection weights that are optimised during the training phase.
We consider response to inputs where, independently for any k, u[k] is (0,0) with probability 1-p or otherwise chosen to be one of the four inputs (± 1,0) or (0,± 1) with equal probability. This generates a sequence of inputs that remain at (0,0) for an exponentially distributed period of time, but occasionally takes one of the four values (0,± 1), (± 1,0).
The target output to be learned by the RNN is a two-dimensional vector z[k] that can assume four possible configurations: (1, 1), (1,-1), (-1,1), and (-1,-1).
The system (<ref>) is considered to successfully accomplish the task if the network is able to reliably and accurately reproduce all 20 possible actions described in Table <ref>.



We use RNNs that are ESNs subjected to supervised training with a perturbation matrix obtained by injecting the output into the ESN dynamics.
This choice does not limit the validity of our hypothesis: we claim that our results are general and can be used to describe the behaviour of more general discrete-time, input-driven RNNs. The contributions of this paper can be summarised as follows:


  * We provide a theoretical framework to describe the behaviour of RNNs by means of ENAs. The proposed theoretical framework is general and covers several types of computational tasks. We present a specific instance of such a framework applied to describe RNN behaviour on tasks requiring to learn a finite set of memory states;


  * We use bifurcation analysis to manually design (low-dimensional) ESNs that give rise to ENAs able to solve the flip-flop task with any number of bits. This allows us to justify our choice of modelling framework and suggests its validity in the context of RNNs;


  * As ESNs are driven by inputs and subject to training via perturbation matrices, bifurcation analysis alone is not sufficient to explain changes to their behaviour. In fact, inputs and perturbation matrices can affect ESN behaviour in a non-trivial way.
To this end, we introduce an algorithm to extract the ENA describing the ESN behaviour for the task at hand.
The algorithm analyses an ESN trajectory and constructs a directed graph encoding the underlying ENA on which the dynamics takes place. The vertices (nodes) of such a graph are associated with the fixed points of the dynamics and directed edges describe excitable connections between them. We apply this algorithm to an ESN that is trained to solve the flip-flop task and show that the resulting ENA is able to explain how ESNs perform computations in a detailed and mechanistic way;


  * We define a notion of excitability threshold for this high-dimensional, non-linear dynamical system driven by inputs. We propose a method for computing this excitability threshold that accounts for inputs and can directly be applied to trajectories generated by a neural network while solving a task;


  * Our simulations suggest three important findings.
First, as already noted by recent research <cit.>, the dynamics of high-dimensional RNNs takes place in a much lower-dimensional phase space region that is determined by the structure introduced with training and inputs.
Here, we observe that the dynamics is indeed low-dimensional, but highlight the fact that additional dimensions are used occasionally to switch between stable states according to control inputs.
This suggests that, for example, a simplistic use of methods based on explained variance to reduce dimensions needs to be avoided.
Second, we show how ENA models describing RNN behaviour can be exploited to provide a mechanistic interpretation of errors occurring during the computation undertaken by RNNs.
Finally, we show how excitability thresholds of the extracted ENAs allow us to assess the robustness of RNNs to noise-induced perturbations.




The remainder of this paper is organised as follows.
Section <ref> introduces the essential background material needed in this paper. Section <ref> shows how to design low-dimensional ESN models that give rise to ENAs able to solve the flip-flop task.
In Section <ref> we propose an algorithm for automatically extracting ENAs directly from an ESN trajectory generated while solving a task.
In Section <ref>, we present and discuss results of the simulations.
Finally, Section <ref> draws conclusions and points to future research directions.
We include three appendices: Appendix <ref> reviews notions of linear stability used throughout the paper.
Appendix <ref> discuss bifurcations of fixed points for low-dimensional ESN maps.
Appendix <ref> provides details of the procedure used to determine fixed points.




§ BACKGROUND





 §.§ Echo State Networks


We consider a specific system of the form (<ref>), corresponding to a discrete-time ESN state-update and related output:

    𝐱[k] =   ϕ(𝐖_r𝐱[k-1] + 𝐖_in𝐮[k]), 
    𝐳[k] =   𝐖_o𝐱[k].

𝐱[k]∈ℝ^N_r is the state, 𝐮[k]∈ℝ^N_i and 𝐳[k]∈ℝ^N_o denote input and output, respectively. The activation function ϕ(·) is applied component-wise; without loss of generality, we consider ϕ = tanh : ℝ⟶ (-1,1).
It is worth mentioning <cit.> that (<ref>) is often written 𝐱[k] = 𝐖_rϕ(𝐱[k-1]) + 𝐖_in𝐮[k]. 
However, <cit.> proved that the two formulations are equivalent up to a change of coordinates; they produce the same discrete-time dynamics.
In this paper, we build on (<ref>) and consider a network of discrete-time leaky-integrator neurons <cit.> of the form:

    𝐱[k] = (1 - α)𝐱[k-1] + αϕ(𝐖_r𝐱[k-1] + 𝐖_in𝐮[k] + ϵ).

Here, α∈(0, 1] is called a leak rate and explicitly sets the time-scale of the ESN <cit.>.
The ϵ term represents additive white Gaussian noise with spherical covariance matrix and unit standard deviation.

The reservoir 𝐖_r∈ℝ^N_r× N_r and input 𝐖_in∈ℝ^N_r× N_i matrices are usually random with i.i.d. entries drawn from uniform or Gaussian distributions <cit.>.
However, in the literature it is possible to find reservoirs with different connection patterns, including deterministic topologies <cit.> and those exploiting the norm-preserving property of orthogonal matrices <cit.>. In our case the read-out matrix 𝐖_o ∈ℝ^N_o× N_r is optimised for the task at hand.
Relevant hyperparameters directly affecting ESN performance include the number of neurons and sparseness of their connections, the spectral radius of the reservoir matrix, and leak rate α <cit.>.
The so-called echo-state property (ESP) <cit.> guarantees the existence and uniqueness of a global attracting trajectory for any input sequence in a compact set. The ESP, although useful in some tasks like forecasting tasks, is in practice difficult to verify and it is usually formulated only for ESNs with state-update of the form shown in (<ref>) and (<ref>).
Therefore, it is not suitable to ESN models and tasks we discuss here; see the following subsection.




 §.§ Training ESNs with low-rank perturbation matrices



Training of RNNs is typically implemented by means of stochastic gradient descent or variations of thereof <cit.>.
Learning long-term dependencies in RNNs with gradient descent is known to be problematic, as a consequence of the so-called vanishing/exploding gradient problem <cit.>. To this end, different approaches have been proposed that can be summarised in two categories: (i) methods using gating mechanisms (such as Long Short-Term Memory <cit.> and Gated Recurrent Unit <cit.> networks) and (ii) those based on unitary matrices and constant-slope activations <cit.>.
On the other hand, training of the recurrent layer in ESNs is typically realised by perturbing a randomly-initialised reservoir with a low-rank, deterministic matrix.
This is conventionally accomplished by feeding back the ESN output to the recurrent layer <cit.> or, as <cit.> recently proposed, by designing the reservoir directly as 𝐖_r=𝐗 + 𝐃, where 𝐗 is a random matrix and 𝐃 is a deterministic, low-rank matrix encoding the task of interest.



In this paper, we use a supervised learning algorithm which exploits the feedback of the ESN output as a mechanism for training the recurrent layer.
The state-update equation (<ref>) takes the following form:


    𝐱[k] = (1 - α)𝐱[k-1] + αϕ(𝐖_r𝐱[k-1] + 𝐖_in𝐮[k] + 𝐖_fb𝐲[k-1]  + ϵ),

where 𝐖_fb∈ℝ^N_r× N_o is a matrix with i.i.d. random coefficients usually drawn from a uniform or Gaussian distribution.
Depending on whether training is performed batch or online, 𝐲[k-1] in (<ref>) takes the form of either the target signal or the output produced by the ESN (<ref>), respectively.
In batch mode, it is possible to distinguish two main phases; see Fig. <ref> for an illustration. 
First, the reservoir is fed with an auxiliary input, i.e. the target signal 𝐲.
We construct a matrix 𝐗∈ℝ^N× N_r containing the states 𝐱[k] of the ESN generated in response to target and input signals. Finally, the weights 𝐖_o of the read-out are determined by solving a regularised least-squares problem, 𝐖_o = (( 𝐗^⊤𝐗 + λ^2 𝐈)^-1𝐗^⊤𝐲)^⊤, where 𝐈 is an N_r× N_r identity matrix and λ≥0 is a regularisation parameter.
Successively, the target signal is replaced by the generated output 𝐳; this “closed-loop phase” corresponds to the test phase of the trained ESN. An analysis of the stability of the transition from open- to closed-loop can be found in <cit.>.



We call the trained reservoir the following matrix

    𝐌 := 𝐖_r + 𝐖_fb𝐖_o.


Once the read-out matrix 𝐖_o is optimised, by imposing 𝐲[k]=𝐳[k] and expanding in (<ref>) with (<ref>), we obtain:

    𝐱[k] =    (1 - α)𝐱[k-1] + αϕ(𝐖_r𝐱[k-1] + 𝐖_in𝐮[k] + 𝐖_fb𝐖_o𝐱[k-1] + ϵ)= 
      
     =    (1 - α)𝐱[k-1] + αϕ(𝐌𝐱[k-1]  + 𝐖_in𝐮[k] + ϵ).


The trained reservoir (<ref>) is obtained by adding a matrix 𝐖_fb𝐖_o∈ℝ^N_r× N_r, which is low-rank N_o≪ N_r relative to the randomly initialised reservoir matrix 𝐖_r. Therefore the reservoir is in some sense trained using output feedback connections.




Inputs 𝐮[k] in (<ref>) play the role of control inputs and are typically constant or impulsive signals.
The ESN read-out matrix 𝐖_o is conventionally determined by solving a regularised least-squares problem.
Nonetheless, we note that also online training schemes have been developed, e.g., the FORCE learning algorithm originally introduced in <cit.> and further extended by <cit.>.
During the test phase, regardless of the adopted training mechanism, the state-update of ESNs is described by (<ref>).
In this paper, we consider batch training via ridge regression and analyse the trajectory generated by (<ref>) during the test phase.

It is worth nothing that our theoretical framework does not rely on a particular training method or a particular RNN architecture. We note that complicated (trained) RNN models may be described using a noisy nonautonomous dynamical system, which in our system is represented by (<ref>). We focus on ESNs as they are the simplest RNN models to test our hypothesis. For this reason, the terms ESN and RNN are used interchangeably.





 §.§ Network attractors



Many common dynamical systems encountered in nature are dissipative <cit.>.
In such systems, the absence of any conservation law means that typically the system evolves towards an attracting set of dimension strictly less than the original phase space dimension; such a set (or a particular subset of it) is commonly called an attractor: formal definitions are discussed for example in <cit.>. The basin of attraction of an attractor is the set of all initial conditions from which the system evolves toward the attractor.
Attractors convey crucial information about the behaviour of the dynamical systems which have generated them.  

Here, we consider the following noise-free discrete-time dynamical system with inputs:

    𝐱[k]=𝐆(𝐱[k-1], 𝐮[k]),

where 𝐆: ℝ^N_r×ℝ^N_i⟶ℝ^N_r is related to (<ref>) as follows:

    𝐆(𝐱, 𝐮)= 
    [       (1-α)x_1 + αϕ(𝐌_(1)·𝐱 + (𝐖_in)_(1)·𝐮);       (1-α)x_2 + αϕ(𝐌_(2)·𝐱 + (𝐖_in)_(2)·𝐮);                                           ⋮; (1-α)x_N_r + αϕ(𝐌_(N_r)·𝐱 + (𝐖_in)_(N_r)·𝐮) ],

where is 𝐌 the trained reservoir matrix (<ref>), the subscript (i) denotes the ith rows of a matrix and · the usual dot product.
As mentioned before, the input signal for the flip-flop task is null most of the time, at which point it is governed by the autonomous dynamics of

    𝐱[k] = 𝐅(𝐱[k-1]) = 𝐆(𝐱[k-1], 0).

Fixed points 𝐩∈ℝ^N_r are solutions of 𝐅(𝐩) = 𝐩.
Related to the notion of attractor is the notion of limit set <cit.>, thus we introduce the following 


The ω-limit set of a point 𝐱_0 under the iterated map (<ref>) is defined by

    ω_𝐅(𝐱_0) := ⋂_n ∈ℕ{𝐅^h(𝐱_0)   |    h > n }.



The ω-limit set of a point 𝐱_0 is the set of limit points of the forward trajectory {𝐅^h(𝐱_0) }_h ∈ℕ. For a given fixed point 𝐩, its basin of attraction is formed by all points 𝐱∈ℝ^N_r such that ω_𝐅(𝐱) = {𝐩}. If such a fixed point is stable, then there exists a neighbourhood of 𝐩 whose ω-limit set corresponds to this fixed point.




More complex attractors consisting of networks of invariant sets in phase space have been proposed in the literature <cit.>.
Such models found renewed interest in neuroscience <cit.> and other fields of research, as they provide a fundamental tool to describe dynamic processes occurring on transients that explore excitable connections.
More relevant to our paper, we focus on networks of stable fixed points that are connected by excitable connections.
Following <cit.>, we say that there exists an excitable connection for amplitude δ>0 from stable fixed point 𝐩_i to 𝐩_j whenever

    B_δ(𝐩_i) ∩ W^s(𝐩_j) ≠∅,

where B_δ(𝐩_i) stands for the closed ball centred on 𝐩_i with radius δ > 0 and W^s(𝐩_j) = {𝐱∈ℝ^N_r   |   ω_𝐅(𝐱) = {𝐩_j }} denotes the basin of attraction[The basin of attraction corresponds to the stable manifold whenever 𝐩_j is hyperbolic. This has interior if 𝐩_j is an attractor.] of the fixed point 𝐩_j. 


We define excitability threshold <cit.>  (or just threshold) of the excitable connection from 𝐩_i to 𝐩_j, and denote it as δ_th(𝐩_i,𝐩_j), the following nonnegative real number:

    δ_th(𝐩_i,𝐩_j) := inf{δ>0    :    B_δ(𝐩_i)∩ W^s(𝐩_j)  ≠∅}.



The quantity (<ref>) can be informally interpreted as the fact that the fixed point 𝐩_i is δ_th(𝐩_i,𝐩_j) away from the basin of 𝐩_j; see Fig. <ref> for a visual explanation.







A set X_exc⊂ℝ^N_r is called an excitable network attractor (ENA) for amplitude δ>0 if there exists a collection of fixed points {𝐩_i }_i=1^M, such that
    
    X_exc({𝐩_i }_i=1^M, δ) := ⋃_i,j=1 
     i≠ j^M{𝐅^h( B_δ(𝐩_i) ) }_h≥ 0 ∩  W^s(𝐩_j),

where 𝐅( B_δ(𝐩_i) ) := {𝐅( 𝐱 )    |   𝐱∈ B_δ(𝐩_i) } is based on (<ref>) (see <cit.>).

The autonomous dynamics on such a set converges to one of the stable fixed points; hence, external inputs are necessary to get interesting dynamics.
If the external input is large enough, then the state will escape from the current basin of attraction and switch to a different one, until another sufficiently large input will lead to another change of basin.

The excitability threshold (<ref>) is defined as the (Euclidean) distance between a given stable fixed point, 𝐩_i, and the basin of attraction of another fixed point, say 𝐩_j.
Such a quantity measures the minimum distance in phase space necessary to escape from the basin of 𝐩_i and converge towards 𝐩_j.
Nevertheless, if the dynamical system is high dimensional, then there will be a large number of possible escaping directions that could be exploited by inputs.
Therefore, excitability thresholds (<ref>) alone may not be representative of nonautonomous systems driven by inputs.
For this purpose, in order to take into account the action of inputs on the dynamics, we introduce the notion of input-driven excitability threshold of an excitable connection.
Considering K as a compact subspace of ℝ^N_i, we define 𝒢(𝐱_0 ; K) := ⋃_𝐮∈ K𝐆(𝐱_0 , 𝐮 ), where 𝐆 is the function in (<ref>) defining the nonautonomous dynamical system which describes the trained neural network.
The set 𝒢(𝐱_0 ; K ) contains all states reachable from 𝐱_0 under the action of input values 𝐮∈ K.
Importantly, the presence of noise has the effect to perturb away the internal state from the exact location of the stable point in the deterministic counterpart of the dynamics. Therefore, instead of 𝒢(𝐩_i; K) we rather observe the following set.


We call local input response set of the stable point 𝐩_i the subset of phase space defined by 
    
    𝒢(B_r(𝐩_i); K ) := ⋃_𝐱∈ B_r(𝐩_i)𝒢(𝐱; K) ,

where the radius r can be shrunk according to the amplitude of the noise and K is a subset of admissible input values for the task at hand.


Continuity of 𝐆 implies that, if r → 0^+, then the monotonically decreasing sequence of sets {𝒢(B_r(𝐩_i); K ) }_r≥ 0 converges to 𝒢(𝐩_i; K ) regardless of K ⊂ℝ^N_i.
Furthermore, we note that 𝒢(B_r(𝐩_i); K ), as a subspace of ℝ^N_r, is compact if K⊂ℝ^N_i is compact.

If K represents all possible inputs of a particular task[For example, in the flip-flop task we have K={ (0,0),   (1,0),   (-1,0),  (0,1),  (0,-1) }], then we can drop it and denote (<ref>) simply as 𝒢(B_r(𝐩_i) ). 
Therefore, the subregion of the phase space represented by the local input response set 𝒢(B_r(𝐩_i)), encodes the action exercised by the input when the internal state of the RNN is nearby the stable point 𝐩_i.


We define input-driven excitability threshold of the excitable connection from 𝐩_i to 𝐩_j, and denote it as δ_inp(𝐩_i,𝐩_j), the following nonnegative real number:

    δ_inp(𝐩_i,𝐩_j) := inf{δ>0    :   𝒢(𝐩_i)∩ B_δ(𝐩_i)∩ W^s(𝐩_j)  ≠∅}.



The excitability threshold in (<ref>) has a similar meaning to the one defined in (<ref>), although it considers only the subspace exploited by the action of inputs nearby 𝐩_i, where the excitable connection starts from, see Fig. <ref>.
Finally, from the fact that 𝒢(𝐩_i)∩ B_δ(𝐩_i)∩ W^s(𝐩_j) ⊆  B_δ(𝐩_i)∩ W^s(𝐩_j), it follows that δ_th(𝐩_i,𝐩_j)  ≤δ_inp(𝐩_i,𝐩_j) holds for all pairs of fixed points. 







§ DESIGNING LOW-DIMENSIONAL ESNS TO SOLVE FLIP-FLOP TASKS



In this section, starting from ESN models we show how to manually design ENAs that realise computations needed for the flip-flop task. This step allows us to justify the choice of the modelling framework adopted here and its validity in the context of RNNs.
For this purpose, we rely on bifurcation theory; see Appendix <ref> for technical notions.
A bifurcation <cit.> is a qualitative change in the solution set (phase space topology) on changing a parameter of a dynamical system. On varying the parameters of the model, if such qualitative changes appear, then we say the system has undergone a bifurcation.
For this reason, the notion of bifurcation plays an important role in training RNNs and in describing their behaviour. 
Training the recurrent layer of RNNs corresponds to shaping their phase space topology, possibly inducing bifurcations that lead to a qualitative change in behaviour.
We argue that adding a low-rank perturbation matrix to the reservoir (<ref>) can induce bifurcations that lead to the creation of attracting regions in ESN phase space useful to solve the task at hand. Clearly, this can also happen with more sophisticated training algorithms.
Let us assume the origin as the only global attractor for the autonomous system 𝐱[k] = (1 - α)𝐱[k-1] + αϕ(𝐖_r𝐱[k-1] ) associated with (<ref>).
Then, the bifurcation induced by (<ref>) leads to a transition from such a trivial dynamic to one where the origin becomes unstable[The origin could in principle remain stable and other attractors appear elsewhere: an example can be found in <cit.>.] and repels trajectories towards other attracting regions in phase space, where the nonautonomous dynamics actually takes place.

In Sec. <ref>, we provide a low-dimensional example of an ESN that is able to solve the flip-flop task; the reservoir of such ESN is formed by two neurons. In Sec. <ref>, instead we design an ESN model with a reservoir formed by 2k neurons which is able to solve the flip-flop task with k bits. For both examples we show there is an underlying ENA that explains their dynamics.




 §.§ A minimal-dimension example to solve the two-bit flip-flop task



In order to solve the k-bit flip-flop task, one needs to learn 2^k memory states (stable fixed points) and the related switching patterns dictated by control inputs. Here, we show how to design an ESN with two neurons in the recurrent layer, giving rise to an ENA able to solve the flip-flop task with k=2 bits.
According to the analysis of Appendix <ref>, we know that it is possible to obtain, with k neurons, up to 3^k fixed points, 2^k of which are stable, 3^k - 2^k -1 are saddles, and 1 (the origin) is a repeller.
Therefore, we set the trained reservoir weights (<ref>) according to condition (<ref>).
In particular,

    𝐌(b) := ω_r[[ 1 b; b 1;   ]]

with ω_r=3 scaling the reservoir weights, and b acting as tuning parameter.
As shown in Fig. <ref>, for every 0 ≤ b < 0.47, the autonomous system 𝐱[k] = tanh(𝐌(b) 𝐱[k-1] ) has four stable attractors located near the vertices of the invariant square [-1,1]^2.



The input is injected into the ESN via the following weight matrix, 𝐖_in=ω_in𝐈_2, i.e. a scaled 2× 2 identity matrix, setting the scaling factor (called a hyperparameter in machine learning) to ω_in=6.
Let us set the remaining parameters in (<ref>) as ϵ=0 and α=1, and let the ESN output (<ref>) be the identity mapping.
Let us assume that, at time k, the system is in state 𝐩=(p_1,p_2), which is close to one of the four attractors, i.e., 𝐩_1 ≈ (1,1),  𝐩_2 ≈ (-1,1),  𝐩_3 ≈ (-1,-1),  𝐩_4 ≈ (1,-1). 
The possible, non-null inputs at time-step k are 𝐮[k]∈{(1,0), (0,1), (-1,0), (0,-1)}. 
The action of the input pulse is encoded in a vector Δ∈ [-1, 1]^2, representing the difference between the state before and after the occurrence of such a pulse, whose components are

    Δ_j ≈tanh( 3 p_j +  6 u_j  ) - tanh(  3  p_j ) ,    j=1,2,

considering ω_r = 3, ω_in = 6 and, for the sake of simplicity, b=0, which corresponds to the case of two independent neurons.
Hence, the switching mechanism between attractors is ruled by the signs of both p_j and u_j; the former encoding the position and the latter the direction where to move.
Therefore, the state changes if and only if p_j and u_j assume different signs for some j∈{1,2}.
In fact, if they have the same sign, then (<ref>) is null because tanh( (p_j) [ 3 + 6 ] ) ≈tanh( (p_j)  3 ) due to saturating activation functions. While, if they have different sign, then (<ref>) becomes close to -1 or 1, according to (p_j), because tanh( (p_j)  [ 3 - 6 ]  ) = -tanh( (p_j) 3 ). 
Clearly, it is not necessary that |ω_r-ω_in|= ω_r holds, as in our set up. 
Although we assumed b=0 for clarity of explanation, the conclusion is the same for each b ∈ [0,0.47). However, when the parameter approaches the bifurcation value b ≈ 0.47, the escaping dynamics from certain fixed points become significantly slower.

Given a set of stable fixed points, a δ>0 gives rise to a specific ENA; see definition of ENA in (<ref>).
For instance, a large δ will potentially activate all excitable connections between the attractors.
However, in order to properly solve the flip-flop task, some of the connections need not to be active, namely, the diagonal connections between 𝐩_1 and 𝐩_3, and between 𝐩_2 and 𝐩_4.
Due to symmetry, there are only three different excitability thresholds that an excitable connection can have, that is, δ_th(𝐩_2,𝐩_1) , δ_th(𝐩_1,𝐩_2), δ_th(𝐩_1,𝐩_3). In the particular configuration shown in Fig. <ref>, it holds that δ_th(𝐩_2,𝐩_1) < δ_th(𝐩_1,𝐩_2) < δ_th(𝐩_1,𝐩_3).
Therefore, the underlying ENA supporting the nonautonomous ESN dynamics is defined as X_exc ({𝐩_i }_i=1^4, δ), for every δ_th(𝐩_1,𝐩_2) < δ < δ_th(𝐩_1,𝐩_3).
We note that δ_th(𝐩_1,𝐩_3) ≈√(2), regardless of b ∈ [0,   0.47).
We can reduce the size of the basin of attraction of 𝐩_2 and 𝐩_4 by increasing the value of b, which in turn reduces the excitability threshold of the corresponding connections, until at b≈ 0.47 a bifurcation occurs making them disappear[A bifurcation where both saddles collide with the corresponding stable points 𝐩_2 and 𝐩_4, simultaneously annihilating all six fixed points.].
However, the basins of attraction of the other two stable points, 𝐩_1 and 𝐩_3, become bigger (δ_th(𝐩_1,𝐩_2) ≈ 2 - δ_th(𝐩_2,𝐩_1)), so increasing their excitability thresholds. 
Therefore, when the ESN state is close to 𝐩_1 (or 𝐩_3), the input needs to be very precise to throw back the state to the narrow basins of 𝐩_2 and 𝐩_4. Moreover, it must have a large amplitude compared to what is needed to escape from such narrow basins.
Nevertheless, setting ω_in large enough and depending on b (until a value of ω_in=6, which is enough for every 0 < b < 0.47), the system is able to reproduce the flip-flop dynamics without errors.

Unfortunately, it does not seem to be possible to design a two-dimensional ESN for the two-bit flip-flop where the excitability of each attractor can be tuned by changing the location of the nearby saddles. Nevertheless, as we will show in the next section, this becomes possible by including two additional dimensions in the model.




 §.§ A 2k-dimensional model for k-bit flip-flop tasks



In this section, we propose a 2k-dimensional model able to solve k-bit flip-flop tasks.
For the sake of clarity, but without loss of generality, we show the k=2 bit case.
The key insight is to model the switching dynamics between two fixed points in a two-dimensional space and then build a model formed by two independent systems with two-dimensional switching dynamics.

Unlike the previous example, here we want to tune the excitability of all connections. This can be obtained by imposing the condition in (<ref>) and tuning the reservoir weights to change the position of the saddles.
Therefore, we define

    𝐁 :=
    [[ 1.1   4;  -s   4;     ]],

where s is a real parameter. Hence, for 0≤ s<2.15, the autonomous dynamical system defined by 𝐱[k] = tanh(𝐁𝐱[k-1] ) has one repeller (the origin), two attractors, namely 𝐩_1≈(1,1),  𝐩_2≈ (-1,-1), and two saddles.
By increasing s within the [0, 2.15) interval, we can make the saddles closer to the respective stable points, see Fig. <ref>, hence decreasing the excitability thresholds of these attractors, until a saddle-node bifurcation occurs approximately at s=2.15, annihilating attractors with saddles.




Let us define the input matrix as follows:

    𝐖_in :=ω_in[[ 0 0; 1 0;   ]],

where ω_in is a positive real parameter.
For instance, setting s=2 and ω_in=1, the two-dimensional dynamical system defined by 𝐱[k] = tanh(𝐁𝐱[k-1] + 𝐖_in𝐮[k] ) is able to accomplish the switching mechanism between 𝐩_1 and 𝐩_2 according to inputs (1,0), (-1,0), and ignore other inputs, i.e., (0,1), (0,-1).
Of course, replacing zeros in 𝐖_in with relatively small values (compared to ω_in) does not change the results.

The complete system able to solve the two-bit flip-flop dynamics can be obtained by defining the reservoir as the following four-dimensional block diagonal matrix,

    𝐌 := [[ 𝐁 0; 0 𝐁 ]],

where 0 denotes a 2× 2 matrix with all zeros.
Starting from a given initial condition, for every 0≤ s<2.15, the state of the four-dimensional autonomous dynamical system 𝐱[k] = tanh(𝐌𝐱[k-1] ) converges to one of these four attractors (𝐩_1,𝐩_1), (𝐩_1,𝐩_2), (𝐩_2,𝐩_1), (𝐩_2,𝐩_2). 
In order to produce a two-dimensional output (<ref>) suitable for the task at hand, we define 𝐖_o := [ 1 0 0 0; 0 0 1 0 ], which basically corresponds to a projection onto the first and third component, i.e., (x_1, x_2, x_3, x_4) ⟼ (x_1, x_3).
Finally, we define the input matrix as

    𝐖_in := ω_in[[ 0 0; 1 0; 0 0; 0 1 ]].

We considered the case of a 4-dimensional system composed by two independent, two-dimensional systems.
Nevertheless, the dynamics remains qualitatively the same even if the coupling between these two-dimensional systems is weak, i.e., if the zero matrices in (<ref>) are replaced by matrices with relatively small (absolute) values.
With those definitions in mind and, as before, ϵ=0 and α=1, the ESN ruled by (<ref>) and (<ref>) correctly implements the flip-flop task with two bits.

As the dynamics of systems (x_1, x_2) and (x_3, x_4) are independent from each other, the set of fixed points of the overall dynamics turns out to be the Cartesian product of the set of fixed points of (x_1, x_2) and the fixed points of (x_3, x_4) system. This gives rise to a large number of fixed points: there are 4 stable points, 8 saddles with 1 unstable directions, 8 saddles with 2 unstable directions, 4 saddles with 3 unstable directions and 1 repeller.
However, the set of fixed points where the nonautonomous flip-flop dynamics take place is formed by 4 stable fixed points and 8 saddles with only one unstable direction; see Sec. <ref> for a detailed example.
Every stable fixed point is close to a pair of saddles and, due to symmetry, they are all at the same distance δ_th(𝐩_1,𝐩_2).
This quantity defines the excitability thresholds of the connections needed in the flip-flop task, and therefore implicitly defines the ENA ruling the behaviour of the four-dimensional ESN.







§ EXTRACTING ENAS FROM THE ESN TRAJECTORY



In this section, we describe the proposed algorithm to extract an ENA from an ESN trajectory. The proposed algorithm, which is schematically illustrated in Fig. <ref>, takes the trained reservoir matrix 𝐌 (<ref>) and a trajectory of the nonautonomous ESN (<ref>) as input and produces a weighted directed graph, representing the ENA describing the ESN behaviour, as output.
The algorithm is composed of two main steps.
In Sec. <ref>, we describe the procedure to find fixed points of the ESN dynamics (corresponding to the vertices of the graph).
Successively, in Sec. <ref>, we show how to determine the excitable connections and related thresholds between stable fixed points (corresponding to the weighted directed edges).







 §.§ Finding fixed points of the dynamics



The optimisation algorithm we have used to find fixed points is based on <cit.>; see <cit.> for an open-source Tensorflow toolbox for finding fixed points in arbitrary RNN architectures and <cit.> for an alternative method to identify fixed points.
The key idea is to define a scalar function whose minima correspond to fixed points of the ESN dynamics.



We call velocity field of the autonomous system (<ref>) the vector field 𝐐 : ℝ^N_r⟶ℝ^N_r, defined as
    
    𝐐(𝐱):= 𝐅(𝐱) - 𝐱,

with 𝐅 being the map in (<ref>).

Therefore, the velocity field takes the following form:

    𝐐(𝐱)= α[ tanh(𝐌·𝐱) - 𝐱],

where 𝐐(𝐱[k]) is the vector that needs to be added to the current state 𝐱[k] to obtain the next one.
In fact,

    𝐱[k+1]=𝐅(𝐱[k]) ⟺𝐱[k+1] - 𝐱[k]=𝐅(𝐱[k]) - 𝐱[k]  ⟺𝐱[k+1] = 𝐱[k] + 𝐐(𝐱[k]).




We define kinetic energy of the autonomous system (<ref>) to be the following scalar function,
    
    q(𝐱) :=12‖𝐐(𝐱)‖^2.


Note that fixed points 𝐱^* ∈ℝ^N_r satisfy 𝐐(𝐱^*)=0 if and only if q(𝐱^*)=0. Fixed points are hence identified as the global minima of (<ref>).
We use the quasi-Newton algorithm BFGS <cit.> to minimise (<ref>).
In order to speed-up the optimisation by several orders of magnitude, we explicitly provided the gradient of (<ref>) to BFGS, which reads:

    ∇ q(𝐱_0) =  𝐉_𝐐(𝐱_0)^T𝐐(𝐱_0)     =   α( 𝐃(𝐱_0) 𝐌 - 𝐈_N_r)^T𝐐(𝐱_0),

where 𝐈_N_r is an N_r× N_r identity matrix, 𝐉_𝐐(𝐱_0) denotes the Jacobian matrix of the velocity field (<ref>) and 𝐃(𝐱_0) is a diagonal matrix defined in (<ref>) of Appendix <ref>, both evaluated on 𝐱_0.

The initial conditions for minimising (<ref>) are determined by randomly sampling states from a trajectory of the nonautonomous ESN (<ref>). The convergence of the BFGS algorithm depends on a tolerance.
In fact, the algorithm may converge to similar solutions that, depending on chosen tolerance, are numerically different.
As these solutions represent fixed points of the dynamics, we aggregate them in a meaningful way and return a non-redundant set of fixed points. For this purpose, as post-processing step, we run a clustering algorithm on the set of all solutions and return only cluster representatives; details provided in Appendix <ref>.




 §.§ Determining excitable connections between attractors



Once stable fixed points have been identified, we determine the excitable connections between them.
As discussed before (and in more detail in Appendix <ref>) , the ESN training (<ref>) induces bifurcations that generate new fixed points; and possibly also other attracting regions in phase space that, however, are not explicitly modelled in this work.
The ESN is driven by control inputs that allow us to correctly perform switches between stable states.
As a consequence, we first need to understand how such inputs affect the dynamics from a geometric point of view.
This is done in Sec. <ref> by introducing the notion of local switching subspace (LSS), which is strictly related to the notion of local input response set, (<ref>), introduced in Sec. <ref>.
Then, in Sec. <ref> we describe how we determine all excitable connections that are relevant for the task under consideration and compute their excitability thresholds.
The method we propose is based on a grid of points lying in the LSS, accounting for the input action on the dynamics. By simulating the autonomous dynamics with initial conditions taken from such a grid, we are able to approximate input-driven excitability thresholds (<ref>) and also to quantify how likely it is that the RNN uses such connections while solving the task.




  §.§.§ Local switching subspaces





Let us consider an ESN trajectory (<ref>), 𝐱[0], 𝐱[1], 𝐱[2], …, 𝐱[k], …, obtained with inputs 𝐮[1], 𝐮[2], …, 𝐮[k+1], ….
Moreover, let us denote with 𝒦 :={k_i}_i ∈ℕ the set of indices for which 𝐮[k_i] ≠0.
We define pulse difference vector (PDV) a vector containing the difference between pre- and post-input states, namely  
    
    𝐱[k_i] - 𝐱[k_i-1],    k_i ∈𝒦.


The PDVs (<ref>) convey relevant information about the action of inputs on ESN state and we exploit such information to determine the excitable connections and related thresholds.

We remind the reader that the number of non-zero inputs is controlled by the parameter p of the exponential distribution (see (<ref>) and related discussion in the text), which in turn determines the (average) number of PDVs available for the following analysis.


In order to compute only those connections that are actually used by the ESN while solving the task, we need to focus on the action of inputs in the neighbourhood of each attractor.
To this end, we consider an Euclidean ball B_r( 𝐩 ) of radius r≥0 centred on an attractor 𝐩, and call local PDVs of 𝐩 the finite set of PDVs that originate inside B_r(𝐩).
Therefore, referring to (<ref>), the local PDVs of 𝐩 is a collection of vectors originating in B_r(𝐩) and ending in 𝒢( B_r( 𝐩 )).

 
Let 𝕃(𝐩) be the vector space obtained by means of principal components analysis of the local PDVs (<ref>) of the attractor 𝐩, retaining only l≪ N_r principal components.
We define the local switching subspace (LSS) of 𝐩, and we denote it as 𝐩 + 𝕃(𝐩), the affine space composed by attaching the vector space 𝕃(𝐩) over 𝐩; see Fig. <ref>, top panel.








  §.§.§ Estimation of excitability thresholds



The idea is to sample the LSS of a stable point 𝐩_i and describe it as a grid of points G(𝐩_i). Then, we consider those points on the grid as initial conditions for the autonomous system (<ref>), which is then iterated for a sufficiently large number of steps to ensure convergence to nearby attractors. Finally, tracing the evolution of these initial conditions allows us to estimate excitability thresholds and other relevant quantities.

The proposed algorithm for finding excitability thresholds is graphically illustrated in Fig. <ref>.
The algorithm is based on a hypercube H(𝐩_i) centred on attractor 𝐩_i that is contained within the LSS 𝐩_i + 𝕃(𝐩_i).
The length of the hypercube is such that H(𝐩_i) contains the projection of 𝒢( B_r(𝐩_i) ) on 𝐩_i + 𝕃(𝐩_i).
In order to estimate excitability thresholds, we consider a mesh with a pre-defined density of points on the hypercube H(𝐩_i), thus obtaining a grid of points, G(𝐩_i) = {𝐠_j^i}. To simplify the notation, we use a single index j to enumerate points of the grid.
Through the ω-limit set (<ref>) of the grid G(𝐩_i), that is:

    ω_𝐅(G(𝐩_i)) = ⋃_jω_𝐅(𝐠_j^i),

we can compute the following nonnegative integer c(𝐩_i):= | ω_𝐅(G(𝐩_i)) | - 1, which counts the number of excitable connections originating from 𝐩_i which are allowed to be activated through input.
Indeed, ω_𝐅(G(𝐩_i)) is composed of a set of stable fixed points, {𝐩_i_0, 𝐩_i_1, 𝐩_i_2, …, 𝐩_i_c(𝐩_i)}, determining the endpoints of the c(𝐩_i) different excitable connections.

Note that if the grid is sufficiently dense then the attractor itself is always included in such a set of fixed points. However we do not count this as an excitable connection. In what follows, we assume that the attractor is indexed by i_0, i.e., 𝐩_i_0 = 𝐩_i.


With these definitions in mind, we are ready to compute thresholds of excitable connections.
Let us denote with

    σ_i: {1,…, |G(𝐩_i)| }⟶{i_0, i_1,…, i_c(𝐩_i)},

an indexing function such that 𝐩_σ_i(j) = ω_𝐅( 𝐠_j^i ). Through the preimage σ_i^-1(·), we obtain a partition of points of the grid into the following subsets:

    G(𝐩_i; 𝐩_i_t ) :=  {𝐠_j^i ∈ G(𝐩_i)     |     j ∈σ_i^-1(i_t)  },  t = 0, 1, …, c(𝐩_i).

The points of the grid belonging to the subset defined by (<ref>) are all destined to converge to 𝐩_i_t. Therefore, we estimate the input-driven excitability threshold (<ref>) of the connection from 𝐩_i to 𝐩_i_t as follows:

    δ_inp(𝐩_i, 𝐩_i_t) = min{‖𝐠_j^i - 𝐩_i ‖    :    𝐠_j^i  ∈ G(𝐩_i; 𝐩_i_t )  },    t=1, …,  c(𝐩_i).


The excitability thresholds (<ref>) represent geometric properties of the attractors and related basins in phase space learned through training.
In order to determine the effective excitability (accounting for inputs) of each connection outgoing from 𝐩_i, we exploit the local topology of the LSS of 𝐩_i by means of the grid partition (<ref>) induced by σ_i(·).
To this end, we define the ratio of initial conditions taken from the grid that converged to attractor 𝐩_i_t as follows:

    ν_i,i_t := |G(𝐩_i; 𝐩_i_t )||G(𝐩_i )| - |G(𝐩_i; 𝐩_i)|   ∈ [0,1].


In the limit of infinite number of points in the grid, the quantity (<ref>) converges to the ratio of volumes between the portion of the hypercube belonging to the basin of 𝐩_i_t and the portion of the hypercube that does not belong to the basin of 𝐩_i.
Therefore, dense grids give ratios (<ref>) providing accurate information about how the LSS is distributed between basins of attraction of stable fixed points.
Finally, merging both (<ref>) and (<ref>) into a single expression, we define effective excitability of the connection from 𝐩_i to 𝐩_i_t as follows:

    β_i,i_t :=  ν_i,i_tδ_inp(𝐩_i, 𝐩_i_t).

Note that this quantity takes into account both the distance between the attractor 𝐩_i and the basin's boundary, and the volume of the basin itself.
A low value for β_i,i_t indicates that it is difficult to activate the connection from 𝐩_i to 𝐩_i_t during the task. This may be due (i) to the small volume occupied by the basin of the attractor 𝐩_i_t in the LSS of 𝐩_i or (ii) to a very high excitability threshold associated with such a connection.
On the other hand, β_i,i_t≫ 1 necessarily implies that such a connection has a low excitability threshold, since ν_i,i_t∈ [0,1].
As a consequence, the distance between the basin of attraction of 𝐩_i_t and 𝐩_i is small, thus the connection can be easily activated during the task.




  
Remarks on computational complexity.
There are three parameters controlling the complexity and, accordingly, the accuracy of the search in the grid: the dimension ζ_1 of the hypercube, the length ζ_2 of the hypercube's edge, and the number of points ζ_3 on each edge determining the density of the grid.
ζ_3 and ζ_2 have a linear and polynomial impact on the computational complexity of the algorithm, respectively. However, ζ_1 is more critical as it increases exponentially the number of points in the grid and, accordingly, the overall complexity.
In the simulations we always set ζ_1=l, that is, the dimension of the hypercube is equal to the dimension of the LSS which is low in our case. More generally, the dimension of LSSs depends on the complexity of the inputs and their impact on the dynamics, and this will need to be assessed on a case-by-case basis.





§ SIMULATIONS



In this section we discuss the results of our simulations and relate this to our theoretical framework.
In Sec. <ref>, in order to evaluate the correctness of the algorithm proposed in Sec. <ref>, we apply it to the manually-designed, low-dimensional ESN maps discussed in Sec. <ref> and Sec. <ref>.
Sec. <ref> applies our method to high-dimensional trained ESNs.
We show that, even though the ESN is high-dimensional, the dynamics generated by the trained reservoir (<ref>) is effectively low-dimensional. We also show the usefulness of ENAs for giving a mechanistic interpretation of prediction errors occurring during task execution. Finally, in Sec. <ref>, we show how ENA models can be used to assess the robustness to noise of trained ESNs. For all simulations, we use p=0.1 as parameter of the exponential distribution governing the occurrence of input pulses and set the tolerance of the kinetic energy (<ref>) for detecting minima to 10^-6.




 §.§ Evaluation on manually designed low-dimensional ESNs



For the grid search algorithm, we used ζ_1=2, ζ_2=4, and ζ_3=223.



  
Minimal-dimension model
The LSS determined as described in Sec. <ref> corresponds to the whole 2D phase space.
As a consequence, excitability thresholds (<ref>) match the corresponding input-dependent counterparts (<ref>).
For each attractor, the Euclidean distances from the two closest saddles are consistent with the estimation of excitability thresholds provided by our method.
As discussed in Sec. <ref> and graphically represented in Fig. <ref> bottom-centre panel, we find three excitability thresholds in the computed ENA: 0.49, 1.39, and 1.42. In Fig. <ref>, bottom-right panel, we show that undesired connections (i.e., connections corresponding to state transitions that are not defined in the flip-flop task) between fixed points (0.97,-0.97) and (-0.97, 0.97) have very low effective excitability values (<ref>), indicating that it is unlikely to use such connections during the task execution.
The low effective excitability is due to the small volume ratio (<ref>) of the basins associated with the two attractors. On the other hand, larger thresholds characterise the undesired connections from (-1,-1) to (1,1), reflecting a lack of symmetry between the basin volumes of (1,1), (-1,-1) and (0.97, -0.97), (-0.97, 0.97).






  
Four-dimensional model
The first two principal components obtained via principal component analysis are related to all identified fixed points and produce a cumulative variance ratio larger than 0.96; Fig. <ref> shows a trajectory of the map described in Sec. <ref> and related fixed points projected on the plane spanned by the two principal components.
For every attractor, the computed LSS is a two-dimensional plane; this is expected, since it is actually the plane depicted in Fig. <ref>. Furthermore, we note that none of these planes is aligned with the one where the attractors lie, stressing the importance of defining reference frameworks local to each attractor that take the action of inputs into account.
Fig. <ref>, middle-right panel, shows how the input moves the states out of the plane, using additional dimensions for the switches.
The computed ENA, weighted with excitability thresholds (<ref>) and effective excitability (<ref>), is shown in Fig. <ref>, bottom-left and bottom-right panels, respectively. Symmetries of the dynamical system are clearly present in the resulting directed graphs.
All desired connections are characterised by an excitability threshold of ≃ 0.83, while undesired ones have higher thresholds equal to ≃  1.19. 
We note that the presence of undesired connections does not imply that the ESN actually uses such connections during the task execution.
To this end, the effective excitability thresholds (<ref>), shown in the bottom-right panel of Fig. <ref>, provide us with a more realistic picture of the behaviour under the action of inputs.
Note that the effective excitability thresholds are very low for the undesired connections, implying that the LSS used by inputs is mostly occupied by basins corresponding to attractors adjacent to the end-point of desired connections.






 §.§ Application of the proposed method to high-dimensional trained ESNs



We now consider implementation using ESNs (<ref>) with a reservoir composed of 500 neurons. 
Experimenting with ESNs with different reservoir sizes confirm that one can get ESNs can be successfully trained independent of the precise number of neurons, as long as there are enough. We choose 500 neurons for hardware and computing time considerations.
For the grid search algorithm, we use ζ_1=3, ζ_2=18,and ζ_3=12.
The state-update (<ref>) is configured without leakage, α = 1 and standard deviation of noise ϵ is set to 10^-4 during training.
The entries of matrices 𝐖_in, 𝐖_fb, and 𝐖_r are i.i.d. drawn from a uniform distribution in [-1,1]; the sparseness of 𝐖_r is 95%. Moreover, the reservoir matrix was rescaled to obtain a spectral radius equal to 0.9. Finally, the training set length is always 50000 time-steps.




  
Low-dimensional dynamics
Fig. <ref>, top panel, shows the output produced by an ESN achieving high prediction performance.
The extracted ENA, shown in the bottom-right panel, reveals that undesired connections have excitability thresholds (<ref>) significantly higher than those of desired connections. This means that, in the LSS of every stable point, basins corresponding to attractors adjacent to the end-point of undesired connections stand relatively far away from the stable point compared to basins of attractors adjacent to the end-point of desired connections.
The fixed points of the dynamics lie in a two-dimensional subspace of the 500-dimensional phase space (the cumulative variance ratio is close 1).
It is observed that trajectory spends most of the time close to such a plane. Hence, based on a principal component analysis of the trajectory, we can claim that the dynamics is two-dimensional. 
Nevertheless, during the task execution, the trajectory is occasionally driven away from such a 2D plane by the inputs in order to achieve the switches between attractors, and this feature is crucial to understand how the trained neural network behaves while solving the task. 
The cumulative variance ratio for the identification of the LSS of every attractor, as described in Sec. <ref>, revealed that the switching between stable fixed points takes place in a three-dimensional subspace of the phase space, highlighting that the overall dynamics of ESNs is effectively low-dimensional after training.
It is worth stressing that such LSSs are usually not aligned with the standard coordinate system of the original phase space and the subspaces where attractors lie, suggesting that inputs operate in phase space regions that are disjoint with respect to the low-dimensional linear subspace of the attractors; this is consistent with results reported in <cit.>.






  
Computation accuracy and spurious attractors
Various measures of accuracy exist for quantifying the performance on prediction tasks. For instance, the mean-squared-error (MSE) is typically adopted in tasks involving continuous targets. The MSE is a real-valued scalar that informs us about how close the computed output is to the target one.
However, it is not possible to infer additional insights by looking only at the MSE; for instance, it is not possible to provide a mechanistic description of errors in the computation, i.e., why and how they occur.
Here, we show how the effective ENA extracted from the ESN trajectory can be used to describe how the computation takes place in phase space and, in particular, how to diagnose the nature of prediction errors.

The top panel in Fig. <ref> shows the output produced by an ESN achieving low MSE of the order of 10^-3.
The small errors observable around time-step 13200 can be explained by looking at the ENA model depicted in the bottom panels of Fig. <ref>, which is represented with excitability thresholds (<ref>), left panel, and with effective excitability values (<ref>), right panel.
The directed graphs (whose topology is clearly identical) reveal the presence of two extra stable fixed points in the ESN phase space, which are generated during training.
Nodes of these graphs are coloured according to output values produced by the ESN.
The activation of the related excitable connections brings the ESN to operate in the proximity of such superfluous states, producing inaccurate output values and hence explaining the origin of such errors. Notably, the directed edge – in the graph on the right-hand side – connecting the cyan with the yellow node has a relatively high value of effective excitability.
This means that, in the LSS of the related attractor (cyan), whose output is (-1.03, -1.07), it is relatively easy to transition to the basin of the other attractor (yellow).
This, in turn, produces an output value equal to (-0.84, -0.40), which is significantly different from the target output, i.e., (-1, -1). The prediction error, visible to the naked eye in the top panel around time-step 13200, is indeed due to the activation of that excitable connection.

Both of these spurious attractors act as a sort of surrogates for the correct ones (i.e., those producing lower prediction errors). In fact, once the internal state switches to a spurious attractor, the ESN still behaves consistently.
More precisely, spurious attractors are associated with higher effective excitability values on connections ending up in the correct attractors – see the bottom-right panel of Fig. <ref>. The existence of these spurious attractors and the unravelling of their roles in the ESN computation could not easily be inferred by looking only at the MSE or plotting the outputs produced by the ESN. This demonstrates the importance of ENA models for describing the ESN behaviour.







 §.§ Noise tolerance and effective excitability of ENAs



Here, we aim to provide a further example of the importance of ENAs for characterising the computation of ESNs.
In particular, we ask the following question: given a set of ESNs trained on the same task and achieving the same performance during training, which one will be more robust to noise during test? Typical performance measures, such as MSE, cannot be used to answer such a question and provide useful insights.
We show that an ENA model of the computation, weighted with effective excitability values (<ref>), allows us to assess robustness to noise of ESNs.

As a perturbation, we consider the usual white Gaussian noise term ϵ in the state-update (<ref>), corresponding to perturbations directly applied to all neuron pre-activation values. We stress that such perturbations are applied only during the test phase; training is implemented as described in the previous section.
By increasing the noise standard deviation, the ESN trajectory gets increasingly perturbed neglecting the possibility to reach the proximity of attractors, hence affecting the accuracy of the resulting output values.

We note that: (i) ESNs yielding ENAs with higher effective excitability values are less robust to noise perturbations than ESNs giving rise to ENAs with low effective excitability values (see <cit.>), and (ii) ESNs producing ENAs with balanced edge weights on desired outgoing connections are more robust to noise perturbations than ENAs with unbalanced outgoing weights.
Regarding (i), high excitability values imply the existence of connections with low excitability thresholds.
That is, the basin of the attractor corresponding to the end-point of such a connection is very close to the attractor associated with the origin of the connection, resulting in unnecessary switches that induce errors.
Concerning (ii), for a given attractor, unbalanced outgoing connection weights could be a symptom of unbalanced distribution of space between basins of attraction in the LSS or the existence of some basins significantly closer to the attractor than other basins.
In the latter case, a reasoning similar to (i) applies. On the other hand, if there is indeed an asymmetric distribution of volumes between basins then basins of attractors corresponding to connections possessing large volume ratios will fill most of the LSS, leaving little space for basins of those attractors related to connections with small volume ratios.
Therefore, if the ESN state is close to an attractor with unbalanced outgoing connection weights, but with similar excitability thresholds, then under noise perturbations it is more likely to switch towards an attractor reachable through a connection with high effective excitability even when such a connection should not be activated.
For these reasons, an ENA, where each node is characterised by balanced weights on outgoing connections and very low effective excitability values on undesired connections, provides a prototypical example of ESNs whose behaviour is robust to noise.

To quantify the robustness of this behaviour to noise, we selected two ESNs that solve the two-bit flip-flop task with very high and comparable accuracy – MSE during training is ≃10^-4. We test them on the same input series composed of 10^5 time-steps and recorded their MSEs by considering different noise instances with increasing standard deviation. Directed graphs representing the extracted ENAs are shown in Fig. <ref>.
Results for increasing noise standard deviation are divided in two columns: on the left column, we report results obtained by the ESN that is least tolerant to noise.
The directed graph on the left-hand side of Fig. <ref> shows the presence of two undesired connections. The edge weights of desired connections are not balanced, especially those outgoing from attractors with output values (1,1) and (-1,-1). Conversely, the directed graph on the right-hand side does not contain undesired connections and possess balanced outgoing weights.
The absence of undesired connections indicates that, in the LSS of every attractor, the basins of the attractors corresponding to undesired connections do not exist or, alternatively, they are very small and hence not detectable with the grid density used in our simulations; therefore, they are not relevant for describing the ESN behaviour according with the numerical precision we considered in our simulations.
Finally, we highlight that a performance breakdown for the ESN on the left-hand side panel is observed starting from noise standard deviation of 8×10^-2. On the other hand, the ESN on the right-hand side is significantly more robust to noise, denoting a performance breakdown for noise standard deviation of 1.4×10^-1.








§ CONCLUSIONS



In this paper, we present a novel methodology for modeling and interpreting the behaviour of RNNs driven by inputs.
In order to obtain a mechanistic model describing how RNNs solve tasks, we exploit the theoretical framework offered by excitable network attractors <cit.>, which are defined as networks of stable fixed points connected by excitable connections.
We introduce a procedure to extract excitable network attractors directly from a trajectory generated by a trained RNN.
Such a procedure is composed of two main steps: first, fixed points are computed by solving a non-linear optimisation problem <cit.> and successively excitable connections, with related thresholds, are determined by simulating the dynamics of the autonomous system.

We validate our theoretical developments by considering ESNs trained on the flip-flop task, a simple yet relevant benchmark that consists of learning a prescribed number of stable states and related switching patterns guided by control inputs.
We cannot see any particular theoretical limitations in the application of our framework to more complex RNN architectures, although for more complicated dynamical tasks a detailed computation of bifurcation behaviour may become unfeasible.
Simulation results provide several interesting insights on how RNNs solve tasks and highlight the usefulness of excitable network attractors in describing how RNN undertake computations.
We train echo state networks by means of ridge regression. Our results (not shown here) suggest that the regularisation parameter has a direct impact on the number of attracting regions in phase space generated through training: using too low values produces under-regularised models with a large number of attractors.
An interesting future perspective consists of studying the impact of different training mechanisms (e.g., via FORCE learning or similar online approaches) on the resulting network attractor.

We believe that the proposed modelling framework based on excitable network attractors will be suitable to describe the RNN behaviour for many if not all tasks requiring the learning of a number of attractors (which need not be stable fixed points) and related switching patterns. To this end, in future we will also examine classification tasks. A related next step includes the possibility to handle inputs that are not instantaneous pulses, opening the way to more interesting case studies of practical relevance.
As mentioned in the introductory sections network attractors can in principle be constructed between any type of invariant sets, including limit cycles and strange attractors. Our focus on fixed points was mainly dictated by the fact that computing fixed points from a trajectory is significantly easier than computing limit cycles, for instance. 

Clearly, fixed points are not powerful enough to accurately model all possible behaviours in phase space.
Therefore, an interesting future perspective consists in extending our modeling framework to handle networks composed of heterogeneous attractors, such as a network in phase space connecting fixed points and limit cycles. In turn, this will allow modeling more complex RNN behaviour.
Finally, future directions include embedding the directed graph representing the excitable network attractor extracted from the trajectory in a new phase space, thus producing a set of ordinary differential equations <cit.> describing the RNN behaviour for the task under consideration.








§ LINEAR STABILITY ANALYSIS



By considering the case where the neural network is autonomous, we follow (<ref>) and write the related map as 𝐱[k]=𝐅(𝐱[k-1])=𝐆(𝐱[k-1], 0).
The ith component of the ESN state, x_i, i=1,…, N_r, evolves according to the scalar map F_i(𝐱)=(1-α)x_i + αtanh(𝐌_(i)·𝐱).
By noting that 

    ∂ F_i∂ x_j(𝐱_0)= (1-α)δ_ij + α( 1- tanh^2(𝐌_(i)·𝐱_0) )m_ij,

where δ_ij is the Kronecker delta, it is possible to write the Jacobian matrix evaluated onto 𝐱_0 as 𝐉_𝐅(𝐱_0)= (1-α)𝐈_N_r + α𝐃(𝐱_0) 𝐌, where


    𝐃(𝐱_0)=
    [[   1 - tanh^2 ( 𝐌_(1)·𝐱_0 )                          0                          …                          0;                          0   1 - tanh^2 ( 𝐌_(2)·𝐱_0 )                          …                          0;                          ⋮                          ⋮                          ⋱                          ⋮;                          0                          0                          … 1 - tanh^2 ( 𝐌_(N_r)·𝐱_0 );                            ]]

is an N_r × N_r diagonal matrix representing the squashing action of tanh(·) along the saturating components of 𝐱_0.
By linearizing the network state-update (<ref>) around a given fixed point 𝐱^*, we obtain the linear system δ𝐱[k+1] = 𝐉_𝐅(𝐱^*) δ𝐱[k], where δ𝐱[k] = 𝐱[k] - 𝐱^*.
It is known <cit.> that if a fixed point 𝐱^* is hyperbolic (i.e., 𝐉_𝐅(𝐱^*) has no eigenvalues on the unit circle in the complex plane), then the linear approximation provides a bona-fide characterisation of the non-linear behaviour around that fixed point.
Therefore, the (linear) stability of 𝐱^* is completely determined by the spectral radius of 𝐉_𝐅(𝐱^*). If all eigenvalues of 𝐉_𝐅(𝐱^*) are inside the unit circle, then 𝐱^* is a stable fixed point; on the other hand, if even just one eigenvalues has norm larger than one, then the linearized map is expanding along the corresponding eigenvectors and 𝐱^* is called a saddle[If every eigenvalue has norm greater than one, then the fixed point is a repeller.].
We conclude observing that holds x_i^* =  tanh ( 𝐌_(i)·𝐱^* ) ⟺α x_i^*   = αtanh ( 𝐌_(i)·𝐱^* )  ⟺      x_i^* + α x_i^* =  x_i^* + αtanh ( 𝐌_(i)·𝐱^* )	 ⟺    x_i^*= (1-α)x_i^* + αtanh ( 𝐌_(i)·𝐱^* ), hence the number of fixed points and their positions in phase space do not change on varying α∈ (0,1].
However, their linear stability properties are directly affected by α.





§ BIFURCATION OF FIXED POINTS IN ESNS



In what follows, we perform a bifurcation analysis of one-dimensional ESN map and provide some sufficient conditions to design two-dimensional ESNs with a desired number of fixed points. Particularly, we focus on fold bifurcations[We refer to <cit.> for the terminology; the fold bifurcation is also known as saddle-node, limit point or turning point bifurcation.] of low-dimensional ESN maps, which constitute the only codimension-1 bifurcation that generate new fixed points; as discussed by <cit.>, it is the usual mechanism for creating new attractive fixed point in RNNs.
Moreover, as shown by <cit.>, some fold bifurcations are responsible for reducing the dimensions where the actual dynamics takes place, dividing the RNN parameter space in different regions of effective dimensionality.



 §.§ ESN maps in one dimension

Here, we consider the following one-dimensional map,

    x[k] = tanh(m x[k-1] + w),

where (m, w) ∈ℝ^2 are the bifurcation parameters; we simplify the state-update in (<ref>) and set leak rate α=1, ϵ=0, and remove explicit reference to inputs.
Nevertheless, studying the map (<ref>) is still useful to obtain insights about high-dimensional input-driven ESN dynamics.
In fact, in the high-dimensional case (<ref>) the activation function of the jth neuron is determined by 

    x_j[k] = tanh(m_jjx_j[k-1]     +    ∑_s≠ j m_jsx_s[k-1] + (𝐖_in)_(j)·𝐮[k] +  ε_j   ).

Hence, the parameter w in (<ref>) can be interpreted as the weighted sum of all incoming neurons, plus input and noise terms.
Fixed points of the map F(x) = tanh(m x + w) are the solutions x^* of the equation Q(x^*) = 0, where 

    Q(x) := F(x) - x    =   tanh( m x + w) -x.

It is not possible to find a closed-form expression for the fixed points.
However, it is possible to state that: (i) for every (m, w) ∈ℝ^2, there exists at least one fixed point; (ii) there can be one, two, or three fixed points.
The proof of these statements follow straightforwardly from the fact that lim_x→±∞Q(x)=∓∞ and because, if m<1 then Q is monotonic. Otherwise, there exist two critical points, x_l < x_r, namely

    x_l,r = 1m[  ±tanh^-1(√(m-1m)) - w  ],

such that the function Q(x) folds.
Moreover, since Q(x^*) = 0 ⟺  x^* =  tanh( m x^* + w) ∈ (-1,1) for every ( m, w) ∈ℝ^2, all fixed points lie in the open interval (-1,1).
Critical points in (<ref>) are solutions to Q'(x)=0, or equivalently to F'(x) = 1.

Let us assume m>0. We observe a fold bifurcation whenever a critical point assumes the zero value. The condition Q(x_l,r) = 0 gives rise to the following parametrization of the fold bifurcation curve,

    w_±(m) := ±[  m √(m-1m)   -   tanh^-1( √(m-1m))   ],      m ∈ [1, +∞),

which possesses two symmetric branches ending on a cusp; see Fig. <ref> for an illustration.
Crossing that curve in parameter space towards the region containing the semi-axis of m>1, a new fixed point is formed (x_l or x_r, depending on the branch) and splits in a pair of fixed points, one stable and one unstable[The particular case of crossing that curve through the cusp gives rise to a pitchfork bifurcation of the origin.].
This curve delimits the boundary between a dynamic regime where two stable points coexist with an unstable one, and a regime where only one fixed point exists.
Due to symmetry, in the m<0 case, there are two bifurcation branches identifying a flip or period doubling bifurcation, which is analytically described by the curve w_±(-m) with m≤ -1. Crossing that curve towards the region containing the semi-axis of m<-1, a stable fixed point loses stability and gives rise to a 2-periodic attracting trajectory surrounding it.
We note that the flip bifurcation is detrimental for the flip-flop task considered here, as it gives period rather than fixed point attractors.







 §.§ ESN maps in two dimensions

Here, we consider a two-neuron reservoir. We denote the activation functions of these neurons as x and y, so that the ESN state evolution defines a trajectory (x[k],y[k]) ∈ [-1, 1]^2, k=1, 2, ..., ruled by:

    [ x[k]=tanh(a x[k-1] + b y[k-1] )    ,; y[k]=tanh(c x[k-1] + d y[k-1] )    .;                                      ]


It is know that a fixed point can undergo only fold or flip bifurcations in one-dimensional, discrete-time dynamical systems <cit.>. Nevertheless, considering two or more dimensions, also Neimark-Sacker bifurcations could occur, where a stable point loses stability and an invariant curve surrounding it is created. These are the only possible codimension-1 bifurcations of a fixed point for a discrete-time dynamical system.
Among them, only the fold bifurcation can generate new fixed points.
The origin is always a fixed point of (<ref>). Considering (<ref>), the Jacobian matrix evaluated on the origin reads

    𝐖_r=
        [ a b; c d ].

Therefore, training ESNs via (<ref>) implies a qualitative change of the dynamics around the origin if some eigenvalue crosses the unit circle, i.e., fixed point (0,0) bifurcates.
However, adding a low-rank matrix to the reservoir can induce global bifurcations far away from the origin as well; hence, we cannot rely on the local bifurcation of the origin to deduce the global attractor structure after training.
As a counterexample, suppose λ_1, λ_2 >1. The diagonal matrix diag(λ_1, λ_2) and the upper triangular matrix [ λ_1   γ;   0 λ_2 ] share the same spectrum {λ_1, λ_2}.
Nevertheless, if the coupling is strong enough, that is, |γ|>w_+(λ_1) y^*, where w_+(λ_1) is the function (<ref>) evaluated on λ_1 and y^* is the positive stable solution of y[k]=tanh(λ_2 y[k-1]), then the fold bifurcation curve is crossed, making the dynamics for x variable trivial. As a consequence of this, the upper triangular matrix induces dynamics with just two attractors (plus two saddles and the origin is a repeller), while the diagonal matrix induces four attractors (plus four saddles and the repeller).

In order to count the number of fixed points of the map (<ref>), we can draw its nullclines.
Defining function N_α,β(η) := 1β[-αη +  tanh^-1(η) ], the nullclines are given by

    [ y=N_a,b(x),; x=N_d,c(y),;             ]

and they represent the locus of points where x-dynamic / y-dynamic is stationary.
As a consequence, the solutions of the algebraic nonlinear system (<ref>) coincide with the set of fixed points.
In the last part of this subsection, we show some sufficient conditions to control the number of fixed points assuming a,d>1, i.e., when both nullclines are folding. We refer to Fig. <ref> for a graphical illustration.

    
  * Case b c ≥ 0.

    
    |N^'_a,b(0)| < |N^'_d,c(0)|^-1 ,   i.e.  (1-a)(1-d) < bc   ⟹  there are exactly 3 fixed points.

    On the other hand, when (1-a)(1-d) ≥ bc, there could exist 5, 7 or 9 fixed points.
    Critical points of the function N_α,β(η) are η_± = ±√(α - 1/α).
    Therefore, if (1-a)(1-d) ≥ bc holds then
    
    |N_a,b( √(a - 1/a))|< √(d - 1d)   or   |N_d,c( √(d - 1/d))|<√(a - 1a)  ⟹  there are exactly 5 fixed points.

    Fig. <ref> depicts a nullcline configuration where there are 5 intersections[Note that it holds N_α,β(±√(α - 1/α))= -1/β w_±(α).].
    From that configuration, we can make the humps of y=N_a,b(x) more pronounced by increasing the a/b ratio, until a fold bifurcation occurs, giving rise to a new couple of pair of fixed points.
    The case of 7 fixed points is obtained after the fold bifurcation.
    
    
  * Case b c < 0. 

    
    |N_a,b( √(a - 1/a))|< √(d - 1d)   and   |N_d,c( √(d - 1/d))|<√(a - 1a)  ⟹  there is exactly 1 fixed point.

    From the above configuration, increasing the d/c ratio makes stretch the humps of N_d,c until a fold bifurcation occurs producing exactly three fixed points. Beyond the fold bifurcation, a new pair of fixed points is generated.
    
    |N_a,b( √(a - 1/a))|< √(d - 1/d)   and   |N_d,c( √(d - 1/d))| > 1   ⟹  there are exactly 5 fixed points.[Due to symmetry, this condition holds also when inverting the role of N_d,c, N_a,b and consequently (a,b), (d,c).]



In both cases, a simple sufficient condition to ensure the existence of 9 fixed points, see Fig. <ref>, is 

    |N_a,b( √(a - 1/a))|>1    and   |N_d,c( √(d - 1/d))|>1    ⟹  there are exactly 9 fixed points.

The maximum number of nullcline intersections is 9, setting the maximum number of fixed points that can be generated in a two-dimensional ESN map.







§ AGGREGATION OF FIXED POINTS VIA CLUSTERING



The optimization procedure described in Sec. <ref> provides a number of solutions equal to the number of initial conditions taken into account (assuming all initial conditions converge); however, such solutions might be similar up to a prescribed numerical precision.
In order to reduce the set of all solutions to a few effective fixed points, we run the k-means clustering algorithm and retain only the elements belonging the final clusters minimizing the distance w.r.t. the cluster centroids.
Instead of aggregating all solution at once, we first group such solutions according to their linear stability properties as indicated by the spectrum of the Jacobian matrix of the autonomous map (<ref>).
In particular, we form the group of linear stable fixed points, the group with one unstable direction and so on.
Finally, for each group, we run k-means with parameter k identified according with the minimum of the Davies-Bouldin index <cit.>.




§ COMPLIANCE WITH ETHICAL STANDARDS




  
Conflict of Interest
Andrea Ceni, Peter Ashwin, and Lorenzo Livi declare that they have no conflict of interest.



  
Funding
We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan Xp GPU used for this research.
LL gratefully acknowledges partial support of the Canada Research Chairs program.
PA gratefully acknowledges partial support of the EPSRC Centre for Predictive Modelling in Healthcare via grant EP/N014391/1.



  
Human and Animal Rights
This article does not contain any studies with human participants or animals performed by any of the authors.


abbrvnat



