Echo State Networks with Self-Normalizing Activations on the Hyper-Sphere
Pietro Verzelli Cesare Alippi Lorenzo Livi


[
    [
    Received ; accepted 
========================





empty



§ INTRODUCTION


Although the use of RNN in machine learning is boosting, also as effective building blocks for deep learning architectures, a comprehensive understanding of their working principles is still missing <cit.>.
Of particular relevance are ESN, introduced by Jaeger <cit.> and independently by Maass et al.<cit.> under the name of LSM, which emerge from RNN due to their training simplicity.
The basic idea behind ESN is to create a randomly connected recurrent network, called reservoir, and feed it with a signal so that the network will encode the underlying dynamics in its internal states. The desired – task dependent – output is then generated by a readout layer (usually linear) trained to match the states with the desired outputs. Despite the simplified training protocol, ESN are universal function approximators <cit.> and have shown to be effective in many relevant tasks <cit.>.

These networks are known to be sensitive to the setting of hyper-parameters like the SR, the input scaling and the sparseness degree <cit.>, which critically affect their behavior and, hence, the performance at task.
Fine tuning of hyper-parameters requires cross-validation or ad-hoc criteria for selecting the best-performing configuration. Experimental evidence and some results from the theory show that ESN performance is usually maximized in correspondence of a very narrow region in hyper-parameter space called EoC <cit.>.
However, we comment that beyond such a region ESN behave chaotically, resulting in useless and unreliable computations. 
At the same time, it is everything but trivial configuring the hyper-parameters to lie on the EoC still granting a non-chaotic behavior.
A very important property for ESN is the ESP, which basically asserts that their behavior should depend on the signal driving the network only, regardless of its initial conditions <cit.>.
Despite being at the foundation of theoretical results <cit.>, the ESP in its original formulation raises some issues, mainly because it does not account for multi-stability and is not tightly linked with properties of the specific input signal driving the network <cit.>.

In this context, the analysis of the memory capacity (as measured by the ability of the network to reconstruct or remember past inputs) of input-driven systems plays a fundamental role in the study of ESN <cit.>.
In particular, it is known that ESN are characterized by a memory–nonlinearity trade-off <cit.>, in the sense that introducing nonlinear dynamics in the network degrades memory capacity.
Moreover, it has been recently shown that optimizing memory capacity does not necessarily lead to networks with higher prediction performance <cit.>.

In this paper, we propose an ESN model that eliminates critical dependence on hyper-parameters, resulting in models that cannot enter a chaotic regime. In addition to this major outcome, such networks denote nonlinear behavior in phase space characterized by a large memory of past inputs: the proposed model generates dynamics that are rich-enough to approximate nonlinear systems typically used as benchmarks.
Our contribution is based on a nonlinear activation function that normalizes neuron activations on a hyper-sphere.
We show that the spectral radius of the reservoir, which is the most important hyper-parameter for controlling the ESN behavior, plays a marginal role in influencing the stability of the proposed model, although it has an impact on the capability of the network to memorize past inputs. Our theoretical analysis demonstrates that this property derives from the impossibility for the system to display a chaotic behavior: in fact, the maximum Lyapunov exponent is always null. An interpretation of this very important outcome is that the network always operates on the EoC, regardless of the setting chosen for its hyper-parameters.
We tested the memory of our model on a series of benchmark time-series, showing that its performance for memorization tasks is comparable with that exposed by linear networks. We also explored the memory–nonlinearity trade–off following <cit.>, showing that our model outperforms networks with linear and hyperbolical activations when both memory and nonlinearity are required by the task at hand. 



§ RESULTS




 §.§ Echo state networks


An ESN without output feedback connections is defined as:

 

    
    a⃗_k = W x⃗_k-1 + W_ins⃗_k
       x⃗_k = ϕ(a⃗_k) 
       y⃗_k = W_outx⃗_k



where x⃗_k ∈ℝ^N is the system state at time-step k, N is the number of hidden neurons composing the network reservoir, and W ∈ℝ^N × N is the connection matrix of the reservoir. 
The signal value at time k, s⃗_k ∈ℝ^N_in, is processed by the input-to-reservoir matrix W_in∈ℝ^N × N_in. 
The activation function ϕ usually takes the form of the hyperbolic tangent function, for which the network is a universal function approximator <cit.>.
Also linear networks (i.e., when ϕ is the identity) are commonly studied, both for the proven impact on applications and the very interesting results that can be derived in closed-form<cit.>.
Other activation functions have been proposed in the neural networks literature, including those that normalize the activation values on a closed hyper-surface<cit.> and those based on non-parametric estimation via composition of kernel functions<cit.>.

The output y⃗∈ℝ^N_out is generated by the matrix W_out, whose weights are learned: generally by ridge regression or lasso<cit.> but also with online training mechanisms<cit.>. In ESN, the training phase does not affect the dynamics of the system, which are de facto task-independent. 
It follows that once a suitable reservoir rich in dynamics is generated, the same reservoir may serve to perform different tasks.
A schematic representation of and ESN in depicted in Fig. <ref>





 §.§ Self-normalizing activation function

Here, we propose a new model for ESN characterized by the use of a particular self-normalizing activation function that provides important features to the resulting network.
Notably, the proposed activation function allows the network to exhibit nonlinear behaviors and, at the same time, provides memory properties similar to those observed for linear networks. 
This superior memory capacity is linked to the fact that the network never displays a chaotic behavior: we will show that the maximum Lyapunov exponent is always zero, implying that the network operates on the EoC.
The proposed activation function guarantees that the SR of the reservoir matrix (whose value is used as a control parameter) can vary in a wide range without affecting the network stability.

The proposed self-normalizing activation function is
 

    
    a⃗_k = W x⃗_k-1 + W_ins⃗_k
       x⃗_k = r a⃗_k/‖a⃗_k‖



and leaves the readout (<ref>) unaltered.
The normalization in Eq. <ref> projects the network pre-activation a_k onto an (N-1)-dimensional hyper-sphere 𝕊^N-1_r := {p⃗∈ℝ^N,  ‖p⃗‖ = r } of radius r.  
Fig. <ref> illustrates the normalization operator applied to the state.
Note that the operation (<ref>) is not element-wise like most of activation functions as its effect is global, meaning that a neuron's activation value depends on all other values.





















 §.§ Universal function approximation property

The fact that recurrent neural networks are universal function approximators has been proven in previous works <cit.> and some results on the universality of reservoir-based computation are given in<cit.>. Recently, Grigoryeva and Ortega <cit.> proved that ESN share the same property.
Here, we show that the universal function approximation property also holds for the proposed ESN model (<ref>). 
To this end, we define a squashing function as a map f: ℝ→ [-1,1] that is non decreasing and saturating, i.e., lim_x →±∞ f(x) = ± 1.
We note that ϕ_i(x⃗) := x_i/ ‖x⃗‖ can be intended as a squashing function for each ith component. 
In the same work<cit.>, the authors show that an ESN in the form (<ref>) is a universal function approximator provided it satisfies the contractivity condition: ‖ϕ(x_1) - ϕ(x_2) ‖≤‖ x_1 - x_2 ‖ for each x_1 and x_2, when ϕ is a squashing function.
Jaeger and Haas<cit.> showed that this condition is sufficient to grant the ESP, implying that ESN are universal approximators.

We prove in the Methods section that (<ref>) grants the contractivity condition, provided that:


    σ_min(W) ≥ 1 + ‖ W_in‖‖s⃗_max‖/r


where σ_min(W) is the smallest singular value of matrix W and ‖s⃗_max‖ denotes the largest norm associated to an input.
Eq. <ref> results in a sufficient yet not necessary condition that may be understood as requiring that the input will never be strong enough to contrast the expansive dynamics of the system, leading the network state inside the hyper-sphere of radius r.
In fact, unless the signal is explicitly designed for violating such a condition, it will very likely not bring the system inside the hyper-sphere as long as the norm of W is large enough compared to the signal variance.






 §.§ Network state dynamics: the autonomous case


We now discuss the network state dynamics in the autonomous case, i.e., in the absence of input.
This allows us to prove why the network cannot be chaotic.

From now on, we assume r=1 as this does not affect the dynamics, provided that condition (<ref>) is satisfied. 
From (<ref>), the system state dynamics in the autonomous case reads:


    ϕ(x⃗) := Wx⃗/‖ W x⃗‖


At time-step n the system state is given by

    x⃗_n =  ϕ(x⃗_n-1) = Wx⃗_n-1/‖ Wx⃗_n-1‖ = 
        W/‖ W Wx⃗_n-2/‖ W x⃗_n-2‖‖Wx⃗_n-2/‖ W x⃗_n-2‖ =
        W^2x⃗_n-2/‖ W^2 x⃗_n-2‖


By iterating this procedure, one obtains:


    x⃗_n := ϕ^n(x⃗_0)
        =W^nx⃗_0/‖ W^nx⃗_0‖

where x⃗_0 is the initial state. 
This implies that, for the autonomous case, a system evolving for n time-steps coincides with updating the state by performing n matrix multiplications and projecting the outcome only at the end.

It is worth to comment that this holds also if matrix W changes over time.
In fact, let W_n:=W(n) be W at time time n. Then, the evolution of the dynamical system reads:


    x⃗_n := ϕ^n(x⃗_0)
        =W_n W_n-1… W_2 W_1 x⃗_0/‖ W_n W_n-1… W_2 W_1 x⃗_0‖


Furthermore, note that a system described by matrix W and a system characterized by W' = aW coincide. In turn, this implies that the SR of the matrix does not alter the dynamics in the autonomous case. 




 §.§ Edge of criticality


When tuning the hyper-parameters of ESN, one usually tries to bring the system close to the EoC, since it is in that region that their performance is optimal<cit.>.
This can be explained by the fact that, when operating in that regime, the system introduces rich dynamics without denoting chaotic behavior.

Here, we show that the proposed recurrent model (<ref>) cannot enter a chaotic regime.
Notably, we prove that, when the number of neurons in the network is large, the maximum (local) Lyapunov exponent cannot be positive, hence neglecting the possibility to introduce a chaotic behavior.
To this end, we determine the Jacobian matrix of (<ref>) and then show that, since its spectral radius tends to 1, the maximum LLE must be null.
The Jacobian matrix of (<ref>) reads:


    J_ij =
        ∂/∂ x_iϕ_j(x⃗)
           =
        ∑_l^N
        ∂ϕ_j/∂ a_l∂ a_l/∂ x_i
        = 
        W_ij/‖a⃗‖( 1 - a_i a_j/‖a⃗‖^2)


where the time index k is omitted to ease the notation.
We observe that, asymptotically for large networks (N →∞), we have that a_i/ ‖a⃗‖→ 0 for each i, meaning that the Jacobian matrix reduces to J(x⃗) = W / ‖ Wx⃗‖. As we are considering the case with r=1, we know that ‖x⃗‖ = 1. 

This allows us to approximate the norm of W with its SR ρ = ρ(W),


    J 
        ≈W/‖ Wx⃗‖≈W/ρ.


Under this approximation (<ref>), the largest eigenvalue of J must be 1 as the SR ρ is the largest absolute value among the eigenvalues of W.
We thus characterize the global behavior of (<ref>) by considering the maximum LLE<cit.>, which is defined as:


    Λ := lim_n →∞1/nlog( ∏_k^n ρ_k )


where ρ_k is the spectral radius of the Jacobian at time-step k.
Eq. <ref> implies that Λ = 0 as n→∞, hence proving our claim.

In order to demonstrate that Λ = 0 holds also for networks with a finite number N of neurons in the recurrent layer, we numerically compute the maximum LLE by considering the Jacobian in (<ref>). The results are displayed in Fig. <ref>.
Fig. <ref> panel (a) shows the average value of the maximum LLE with the related standard deviation obtained for different values of SR.
Results show that the LLE is not significantly different from zero. In Fig. <ref> panel (b), instead, we show the Lyapunov spectrum of a network with N = 100 neurons in the recurrent layer, obtained for different SR values.
Again, our results show that the maximum LLE of (<ref>) is zero for finite-size networks, regardless of the values chosen for the SR.







 §.§ Network state dynamics: input-driven case


Let us define the “effective input” contributing to the neuron activation as u⃗ := W_ins⃗.
Accordingly, Eq. <ref> takes on the following form:


    a⃗_k = Wx⃗_k-1+ u⃗_k,

where u⃗_k operates as a time-dependent bias vector.
Let us define the normalization factor as:


    N_k 
    = ‖ W x⃗_k-1 + u⃗_k‖,


Then, as shown in the Methods section, the state at time-step n can be written as:

    x⃗_n = M^(n, 0)x⃗_0 + ∑^n_k=1 M^(n, k)u⃗_k

where

    M^(n, k) := W^n-k/∏_l = k^n  N_l.


By looking at (<ref>), it is possible to note that each u⃗_k is multiplied by a time-dependent matrix, i.e., the network's final state x⃗_n is obtained as a linear combination of all previous inputs. 





 §.§ Memory of past inputs


Here, we analyze the ability of the proposed model (<ref>) to retain memory of past inputs in the state.
In order to perform such analysis, we define a measure of network's memory that quantifies the impact of past inputs on current state x⃗_n.
Our proposal shares similarities with a memory measure called Fisher memory, first proposed by Ganguli et al. <cit.> and then further developed in <cit.>.
However, our measure can be easily applied also to non-linear networks like the one Eq. <ref> proposed here, justifying the analysis discussed in the following.

Considering one time-step in the past, we have:


    x⃗_n = W/N_nx⃗_n-1 + 1/N_nu⃗_n


All past history of the signal is processed in x⃗_n-1.
Note that (<ref>) keeps the same form for every n.
Going backward in time one more step, we obtain:

    x⃗_n = W/N_n(W/N_n-1x⃗_n-2 + 1/N_n-1u⃗_n-1)+ 1/N_nu⃗_n


As ‖x⃗_n‖ = ‖x⃗_n-1‖= ‖x⃗_n-2‖ = 1, we see that there is a sort of recursive structure in this procedure, where each incoming input is summed to the current state and then their sum is normalized. This is a key feature of the proposed architecture. In fact, it guarantees that each input will have an impact on the state, since the norm of the activations will not be too big or too small, because of the normalization.
We now express this idea in a more formal way.

The norm of (<ref>) can be written as

    ‖ M^(n, k)‖ =
        1/∏_l=k^n  N_l‖ W^n-k‖n-k ≫ 1≈1/∏_l=k^n  N_lρ^n-k

where the approximation holds thanks to the Gelfand's formula[The Gelfand's formula <cit.> states that for any matrix norm ρ(A)=lim_k →∞A^k ^1/k.].


If the input is null, we expect each N_l to be of the order of ρ as ‖ x_l ‖ = 1.
So we write:


    N_l = ‖ W x⃗_l-1 + u⃗_l‖
    =ρ + δ_l


where δ_l denotes the impact of the l-th input on the l-th normalization factor. Its presence is due to the fact that without any input N_l would be approximately equal to ρ, while the  input modifies the state and so the normalization value will be modified accordingly.
The value δ_l is introduced to account for this fact: δ_ := N_l - ρ.

If we assume that each input produces a similar effect on the state (i.e. δ_l = δ for every l), we finally obtain:


    ‖ M^(n, k)‖≈ρ^n-k/(ρ+δ)^n-k


We note that such assumption is reasonable for inputs that are symmetrically distributed around a mean value with relatively small variance, e.g. for Gaussian or uniformly distributed inputs (as demonstrated by our results).
However, our argument might not hold for all types of inputs and a more detailed analysis is left as future research.

We define the memory of the network ℳ_α(k|n) as the influence of input at time k on the network state at time n. 
More formally, having defined α := ρ / δ (since (<ref>) only depends on this ratio), we use (<ref>) to define the memory as:


    ℳ_α(k|n) :=
        (α/α+1)^n-k
        =(1 - 1/α + 1)^n-k


Eq. <ref> goes to 0 (i.e., no impact of the input on the states) as α→ 0 and tends to 1 for α→∞, implying that for an infinitely large SR the network perfectly remembers its past inputs, regardless of their occurrence in the past.
Note that (<ref>) does not depend on k and n individually, but only on their difference (elapsed time): the larger the difference, the lower the value of ℳ_α(k|n), meaning that far away inputs have less impact on the current state.




 §.§ Memory loss

By using (<ref>), we define the memory loss between state x⃗_n and x⃗_m of the signal at time-step k (with m>n>k and m = n+ τ) as follows:


    Δℳ(k|m,n) :=   ℳ_α(k|m)  - ℳ_α(k|n) =
    (α/α+1)^m-k - (α/α+1)^n-k
    
    =   (α+1) ^k/α^ kα^ n/(α+1) ^n(  α^τ/(α+1) ^τ - 1) 
    
    =   (1 - 1/α+1)^n (1 + 1/α)^k ( (1 -1/α+1)^τ-1 ) ≤ 0


For very large values of α, we have Δℳ(k|m,n) → 0, meaning that in our model (<ref>) larger spectral radii eliminate memory losses of past inputs.
Now, we want to assess if inputs in the far past have higher/lower impact on the state more than recent inputs.
To this end, we selected n>k>l and defined k = n - a and l = n - b, leading to b > a > 0 and b = a + δ.
Define:


    Δℳ(k, l|n)  =   ℳ(n - a|n) - ℳ(n-a - δ |n) = α^a/(α+1)^a - α^a +δ/(α+1)^a + δ = 
     =
           α^a/(α+1)^a(1 - α^δ/(α+1)^δ) =
        (1 - 1/α + 1)^a (1 - (1 - 1/α + 1)^δ) ≥ 0


We have that lim_δ→∞Δℳ(n - a, n - a - δ |n) = (1 - 1/α + 1)^a, showing how the impact of an input that is infinitely far in the past is not null compared with one that is only a<∞ steps behind the current state. This implies that the proposed network is able to preserve in the network state (partial) information from all inputs observed so far.



  §.§.§ Linear networks
 In order to assess the memory of linear models, we perform the same analysis for linear RNN (i.e., x⃗_n = a⃗_n) by using the definitions given in the previous section.
In this case, an expression for the memory can be obtained straightforwardly and reads:

    ℳ_L(k|n) := ρ^n-k

It is worth noting that there is no dependency on δ and, therefore, on the input.
Just like before, we have the memory loss of the signal at time-step k between state x⃗_n and x⃗_m, as:


    Δℳ_L(k|m,n)    =ℳ_L(k|m) - ℳ_L(k|n)
         
        = ℳ_L(k|n+τ) - ℳ_L(k|n)
         = 
         (ρ^τ -1) ρ^n-k≤ 0


where we set m>n>k and m = n+ τ.
As before, we also discuss the case of two different inputs observed before time-step n.
By selecting time instances n>k>l, k = n - a and l = n - b, we have b > a > 0 allowing us to write b = a + δ.

As such:


    Δℳ_L(k,l|n)    =ℳ_L(k|n) - ℳ_L(l|n)
         
         = 
         (1 - ρ^δ ) ρ^a≥ 0.

We see that in both (<ref>) and (<ref>) the memory loss tends to zero as the spectral radius tends to one, which is the critical point for linear systems.
So, according to our analysis, linear networks should maximize their ability to remember past inputs when their SR in close to one and, moreover, their memory should be the same disregarding the particular signal they are dealing with. 
We will see in the next section that both these claims are confirmed by our simulations.



  §.§.§ Hyperbolic tangent
 
We now extend the analysis to standard ESN, i.e., those using a tanh activation function. Define ϕ:=tanh (applied element-wise). Then, for generic scalars a and b, when |b| ≪ 1 the following approximation holds:


    ϕ(a⃗+b⃗) = ϕ(a⃗) + ϕ(b⃗)/1 + ϕ(a⃗) ϕ(b⃗)≈ϕ(a⃗) + ϕ(b⃗)


When a⃗ is the state and b⃗ is the input, our approximation can be understood as a small-input approximation <cit.>.
Then, the state-update reads:

    x⃗_n    = ϕ(W x⃗_n-1 + u⃗_n) ≈ϕ(W x⃗_n-1) + ϕ(u⃗_n) 
       = ϕ(Wϕ(W x⃗_n-2 + u_n-1)) + ϕ(u⃗_n) 
        ≈ϕ(Wϕ(W x⃗_n-2)) +ϕ(W ϕ(u⃗_n-1))+ ϕ(u⃗_n)


Thus, applying the same argument used in the previous cases, it is possible to write:


    ‖x⃗_n ‖   =  S_ρ^(0)(ϕ(‖u⃗_n) ‖)+S_ρ^(1)(ϕ(‖u⃗_n-1‖)) + … + S_ρ^(n)(ϕ(‖u⃗_0‖))
        
        =
        ∑^n_k=0 S_ρ^(n-k)(ϕ(‖u⃗_k‖))


where we defined:


    S_ρ^(n)(ξ) := ϕ(ρ·…ϕ (ρ·ϕ(_n  timesξ)))



We see that, differently from the previous cases, the final state x⃗_n is a sum of nonlinear functions of the signal u⃗_k.
Each signal element u⃗_k is encoded in the network state by first multiplying it by ρ and then passing the outcome through the nonlinearity ϕ(·).
This implies that, for networks equipped with hyperbolic tangent function, larger spectral radii favour the forgetting of inputs (as we are in the non-linear region of ϕ(·)).
On the other hand, for small spectral radii the network operates in the linear regime of the hyperbolic tangent and the network behaves like in the linear case.

A plot depicting the decay of the memory of input at time-step k for all three cases described above is shown in Fig. <ref>.
The trends demonstrate that, in all cases, the decay is consistent with our theoretical predictions.







§ PERFORMANCE ON MEMORY TASKS

We now perform experiments to assess the ability of the proposed model to reproduce inputs seen in the past, meaning that the desired output at time-step k is y⃗_k = u⃗_k-τ, where τ ranges from 0 to 100. We compare our model with a linear network and a network with hyperbolic tangent (tanh) activation functions.
We use fixed, but reasonable hyper-parameters for all networks, since in this case we are only interested in analyzing the network behavior on different tasks. In particular, we selected hyper-parameters that worked well in all cases taken into account; in preliminary experiments, we noted that different values did not result in substantial (qualitative) changes of the results.
The number of neuron N is fixed to 1000 for all models.
For linear and nonlinear networks, the input scaling (a constant scaling factor of the input signal) is fixed to 1 and the SR equals ρ = 0.95. 
For the proposed model (<ref>), the input scaling is chosen to be 0.01, while the SR is ρ = 15. 
For the sake of simplicity, in what follows we refer to ESN resulting from the use of (<ref>) as “spherical reservoir”.

To evaluate the performance of the network, we use the accuracy metric defined as γ = max{1 - NRMSE, 0}, where the NRMSE is defined as:


    NRMSE := √(⟨‖y⃗_k - ŷ⃗_k‖^2 ⟩/⟨‖y⃗_k - ⟨y⃗_k ⟩‖^2 ⟩).


Here, ŷ⃗_k denotes the computed output and ⟨·⟩ represents the time average over all time-steps taken into account.
In the following, we report the accuracy γ while varying τ in the past for various benchmark tasks.
Each result in the plot is computed using the average over 20 runs with different network initializations. 
The networks are trained on a data set of length L_train=5000 and the associated performance is evaluated using a test set of length L_test=2000.
The shaded area represents the standard deviation. All the considered signals were normalized to have unit variance.






  §.§.§ White noise

In this task, the network is fed with white noise uniformly distributed in the [-1,1] interval. Results are shown in Fig.<ref>, panel (a). We note that networks using the spherical reservoir have a performance comparable with linear networks, while tanh networks do not correctly reconstruct the signal when k exceeds 20.
It is worth highlighting the similarity of the plots shown here with our theoretical predictions about the memory (Eqs. (<ref>), (<ref>), and (<ref>)) shown in Fig. <ref>.



  §.§.§ Multiple superimposed oscillators

The network is fed with Multiple Superimposed Oscillators MSO with 3 incommensurable frequencies:

    u(k) = sin(0.2 k) + sin (0.311 k) + sin (0.42 k)

Results are shown in Fig.<ref>, panel (b). This task is difficult because of the multiple time scales characterizing the inputs<cit.>.
We note the performance of the linear and the spherical reservoirs are again similar and both outperform the network using the hyperbolic tangent.
The accuracy peak when k ≈ 60 is due to the fact that the autocorrelation of the signal reaches its maximum value at that time-step.



  §.§.§ Lorenz series

    The Lorenz system is a popular mathematical model, originally developed to describe atmospheric convection.
    
    
    ẋ= σ (y - x), 
       ẏ= x (ρ - z) - y, 
       ż= x y - β z.

    
    It is well-known that this system is chaotic when σ = 10, β = 8/3 and ρ = 28.
    We simulated the system with these parameters and then fed the network with the x coordinate only.
    Results are shown in Fig.<ref>, panel (c). Also in this case, while the accuracy for spherical and linear networks does not seem to be affected by k, the performance of networks using the tanh activation dramatically decreases when k is large. This stresses the fact that non-linear networks are significantly penalized when they are requested to memorize past inputs.





  §.§.§ Mackley-Glass system

    The Mackley-Glass system is given by the following delayed differential equation:
    
    ẋ(t) = -β x(t) + α x(t - λ)/ 1 + x(t - λ)

    In our experiments, we simulated the system using standard parameters, that is, λ = 17, α = 0.2 and β = 0.1. 
    Results are shown in Fig.<ref>, panel (d). Note that in this case the performance of the network with spherical reservoir is comparable with the one obtained using the hyperbolic tangent and both of them are outperformed by the linear networks.
    


  §.§.§ Santa Fe laser dataset

    The Santa Fe laser dataset is a chaotic time series obtained from laser experiments.
    The results are shown in Fig.<ref>, panel (e). Also in this case the hyperbolic tangent networks do not manage to remember the signal, while the other systems show the usual behavior.





§ MEMORY/NON-LINEARITY TRADE-OFF

Here, we evaluate the capability of the network to deal with tasks characterized by tunable memory and non-linearity features <cit.>. 
The network is fed with a signal where the u_ks are univariate random variables drawn from a uniform distribution in [-1,1]. The network is then trained to reconstruct the desired output of y_k = sin(ν· u_k-τ). We see that τ accounts for the memory required to correctly solve the task, while ν controls the amount non-linearity involved in the computation.
For each configuration of τ and ν chosen in suitable ranges, we run a grid search on the range of admissible values of SR and input scaling. Notably, we considered 20 equally-spaced values of the SR and for the input scaling.
For networks using hyperbolic tangent, the SR varies in [0.2,3]; [0.2,1.5] for linear networks, and [0.2,10] for networks with spherical reservoir. The input scaling always ranges in [0.01,2].

[This choice of exploring large values of SR for the hyperspherical case is motivated by the fact that condition (<ref>) must be satisfied in order for the network to work properly: choosing a large SR will also lead to larger σ_min(W), as discussed in the `Universal function approximation property' section.]

We then choose the hyper-parameters configuration that minimizes (<ref>) on a training set of length L_train=500 and then assess the error on a test set with L_test=200.

In Figures <ref>, <ref>, and <ref>, we show the NRMSE for the task described above for different values of ν and τ, and the ranges of the input scaling factor and the spectral radius which performed best for the hyperbolic tangent, linear and spherical activation function, respectively. 
The results of our simulations agree with those reported in <cit.> and, most importantly, confirm our theoretical prediction: the proposed model possess memory of past inputs that is comparable with the one of linear networks but, at the same time, it is also able to perform nonlinear computations.
This explains why the proposed model denotes the best performance when the task requires both memory and nonlinearity, i.e., when both τ and ν are large.
Predictions obtained for a specific instance of this task requiring both features are given in Fig <ref>, showing how the proposed model outperforms the competitors.

In order to explore the memory-nonlinearity relationship we followed the experimental design proposed in <cit.>: our goal is to study the memory property of the proposed model and not to develop a model specifically designed to maximize the compromise between memory and nonlinearity. The reader interested in models that aim at explicitly tackling the memory–nonlinearity problem may refer to recently-developed hybrid systems <cit.> which try to deal with the memory–nonlinearity trade–off by combining, in different ways, linear and non-linear units so that the network can exploit a combination of them according to the problem at hand. These approaches introduce new hyper-parameters which, basically, allow to control the memory–nonlinearity trade–off. Future works will investigate such a trade–off. 














§ DISCUSSION


In this work, we studied the properties of ESN equipped with an activation function that projects at each time-step the network state on a hyper-sphere.
The proposed activation function is global, contrarily to most activation functions used in the literature, meaning that the activation value of a neuron depends on the activations of all other neurons in the recurrent layer.
Here, we first proved that the proposed model is a universal approximator just like regular ESN and gave sufficient conditions to support this claim.
Our theoretical analysis shows that the behavior of the resulting network is never chaotic, regardless of the setting of the main hyper-parameters affecting its behavior in phase space.
This claim was supported both analytically and experimentally by showing that the maximun Lyapunov exponent is always zero, implying that the proposed model always operates on the EoC.
This leads to networks which are (globally) stable for any hyper-parameter configurations, regardless of the particular input signal driving the network, hence solving stability problems that might affect ESN.

The proposed activation function allows the model to display a memory capacity comparable with the one of linear ESN, but its intrinsic nonlinearity makes it capable to deal with tasks for which rich dynamics are required.
To support this claim, we developed a novel theoretical framework to analyze the memory of the network. By taking inspiration from the Fisher memory proposed by Ganguli et al. <cit.>, we focused on quantifying how much past inputs are encoded in future network states. The developed theoretical framework allowed us to account for the nonlinearity of the proposed model and predicted a memory capacity comparable with the one of linear ESN. This theoretical prediction was validated experimentally by studying the behavior of the memory loss using white noise as inputs as well as well-known benchmarks used in the literature.
Finally, we experimentally verified that the proposed model offers an optimal balance of nonlinearity and memory, showing that it outperforms both linear and ESN with tanh activation functions in those tasks where both features are required at the same time.




§ METHODS




 §.§ Contractivity condition










The universal approximation property, as exhaustively discussed in <cit.>, can be proved for an ESN of the form (<ref>) provided that it has the ESP. 
To prove the ESP, it is sufficient to show that the network has the contractivity condition, i.e., that for each x⃗ and y⃗ in the domain of the activation function ϕ, it must be true that:

    ‖ϕ(x⃗) - ϕ(y) ‖≤‖x⃗ - y⃗‖

Let us introduce the the notation x⃗̂⃗ :=x⃗/ ‖x⃗‖ and ŷ⃗ :=y⃗/ ‖y⃗‖. 
Taking the square of the norm, one gets:

    ‖ϕ(x⃗) - ϕ(y⃗)  ‖^2 = 
      r^2 ( 2 - 2 (x⃗̂⃗·ŷ⃗)) 
        ≤ r^2 ( ( ‖x⃗‖/‖y⃗‖ + ‖y⃗‖/‖x⃗‖) - 2 x⃗̂⃗·ŷ⃗)


Where a⃗·b⃗ is the scalar product between two vectors and the inequality a/b + b/a≥ 2 follows from the fact that (a-b)^2 = a^2 + b^2 - 2ab > 0 for all a,b > 0.
Now, we assume ‖x⃗‖, ‖y⃗‖ > r and show that:


    r^2 ·(‖x⃗‖/‖y⃗‖ + ‖y⃗‖/‖x⃗‖ - 2 x⃗̂⃗·ŷ⃗)
        ≤‖x⃗‖‖y⃗‖(‖x⃗‖/‖y⃗‖ + ‖y⃗‖/‖x⃗‖ - 2 x⃗̂⃗·ŷ⃗) = ‖x⃗‖ ^2 + ‖y⃗‖^2 - 2(‖x⃗‖) (‖y⃗‖) x⃗̂⃗·ŷ⃗_ 2x⃗·y⃗ =
        ‖x⃗ - y⃗‖_2^2


proving the contractivity condition (<ref>).

We see that the only condition needed is ‖x⃗‖, ‖y⃗‖ > r, which means that the linear part of the update (<ref>) must map states outside the hyper-sphere of radius r.
Finally, by applying properties of norms, we observe that:

    ‖W·x⃗_k+ W_in·s⃗_k‖≥‖W·x⃗_k ‖ - ‖ W_in·s⃗_k‖≥σ_min(W)r- ‖ W_in‖‖s⃗_max‖

and asking this to be larger than r results in the condition (<ref>).




 §.§ Input-driven dynamics


Here we show how to derive (<ref>).
Consider a dynamical system evolving according to (<ref>): if we explicitly expand the first steps from the initial state x⃗_0 we obtain:

    x_1   = Wx⃗_⃗0⃗+u⃗_1/‖ W x⃗_⃗0⃗+u⃗_1‖ = Wx⃗_⃗0⃗+u⃗_1/N_1=
    W/N_1x⃗_0 + 1/N_1u⃗_1
    
    x_2   = Wx⃗_⃗1⃗+u⃗_1/N_2 = W^2/N_2 N_1x⃗_0 + W /N_2 N_1u⃗_1 + 1/N_2u⃗_2
    
    x_3   = Wx⃗_⃗2⃗+u⃗_3/N_3 = W^3/N_3 N_2 N_1x⃗_0 + W^2/N_3 N_2 N_1u⃗_1 + W/N_3 N_2u⃗_2 + 1/N_3u⃗_3


So that the general case can be written as:

    x_n   = Wx⃗_n-1+u⃗_n/N_n = W^n/N_n N_n-1… N_1x⃗_0 + W^n-1/N_n N_n-1… N_1u⃗_1 + W^n-2/N_n N_n-1… N_2u⃗_2
    + … + 1/N_nu⃗_n


and in a more compact form as follows:


    x_n = M^(n, 0)x⃗_0 + ∑^n_k=1 M^(n, k)u⃗_k


where


    M^(n, k) := W^n-k/∏_l = k^n  N_l


and N_0 = 1.






§ ACKNOWLEDGEMENTS


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.



§ AUTHOR CONTRIBUTIONS STATEMENT

L.L. outlined the research ideas. P.V. conceived the methods and performed the experiments. C.A. contributed to the technical discussion. All authors developed the theoretical framework, took part to the paper writing and approved the final manuscript.



§ ADDITIONAL INFORMATION

Competing interests The authors declare no competing  interests.





















    bESNbinary ESN
	EoCEdge of Criticality
	ESNEcho State Network
	ESPEcho State Property
	FIMFisher Information Matrix
	LLELocal Lyapunov Exponent
	LSMLiquid State Machine
	MFTMean Field Theory
	MLMachine Learning
	MSOMultiple Superimposed Oscillator
	NRMSENormalized Root Mean Squared Error
	PDFProbability Density Function
	RBNRandom Boolean Network
	RNNRecurrent Neural Network
	RPRecurrency Plot
	SRSpectral Radius




