Echo State Networks with Self-Normalizing Activations on the Hyper-Sphere

TL;DR

This paper introduces a self-normalizing Echo State Network model that eliminates the need for hyper-parameter tuning, prevents chaotic behavior, and maintains high memory capacity, improving robustness and performance.

Contribution

The paper proposes a novel ESN model with self-normalizing activations on the hyper-sphere, ensuring non-chaotic dynamics and high memory retention without hyper-parameter sensitivity.

Findings

Model provably avoids chaotic regimes

Networks exhibit high memory of past inputs

Experimental results confirm rich nonlinear dynamics

Abstract

Among the various architectures of Recurrent Neural Networks, Echo State Networks (ESNs) emerged due to their simplified and inexpensive training procedure. These networks are known to be sensitive to the setting of hyper-parameters, which critically affect their behaviour. Results show that their performance is usually maximized in a narrow region of hyper-parameter space called edge of chaos. Finding such a region requires searching in hyper-parameter space in a sensible way: hyper-parameter configurations marginally outside such a region might yield networks exhibiting fully developed chaos, hence producing unreliable computations. The performance gain due to optimizing hyper-parameters can be studied by considering the memory--nonlinearity trade-off, i.e., the fact that increasing the nonlinear behavior of the network degrades its ability to remember past inputs, and vice-versa. In this paper, we propose a model of ESNs that eliminates critical dependence on hyper-parameters, resulting in networks that provably cannot enter a chaotic regime and, at the same time, denotes nonlinear behaviour in phase space characterised by a large memory of past inputs, comparable to the one of linear networks. Our contribution is supported by experiments corroborating our theoretical findings, showing that the proposed model displays dynamics that are rich-enough to approximate many common nonlinear systems used for benchmarking.

Equations

$${\mathbf{a}}_k=W{\mathbf{x}}_{k-1}+W_{\text{in}}{\mathbf{s}}_k$$

$${\mathbf{x}}_k=\varphi ({\mathbf{a}}_k)$$

$${\mathbf{y}}_k=W_{\text{out}}{\mathbf{x}}_k$$

$${\mathbf{a}}_k=W{\mathbf{x}}_{k-1}+W_{\text{in}}{\mathbf{s}}_k$$

$${\mathbf{x}}_k=r\frac{{\mathbf{a}}_k}{\Vert {\mathbf{a}}_k\Vert }$$

$${\sigma }_{\text{min}}(W)\ge 1+\frac{\Vert W_{\text{in}}\Vert \Vert {\mathbf{s}}_{\text{max}}\Vert }{r}$$

$$\varphi (\mathbf{x}):=\frac{W\mathbf{x}}{\Vert W\mathbf{x}\Vert }$$

$${\mathbf{x}}_n=\varphi ({\mathbf{x}}_{n-1})=\frac{W{\mathbf{x}}_{n-1}}{\Vert W{\mathbf{x}}_{n-1}\Vert }=\frac{W}{\Vert WW\frac{{\mathbf{x}}_{n-2}}{\Vert W{\mathbf{x}}_{n-2}\Vert }\Vert }\frac{W{\mathbf{x}}_{n-2}}{\Vert W{\mathbf{x}}_{n-2}\Vert }=\frac{W^2{\mathbf{x}}_{n-2}}{\Vert W^2{\mathbf{x}}_{n-2}\Vert }$$

$${\mathbf{x}}_n:={\varphi }^n({\mathbf{x}}_0)=\frac{W^n{\mathbf{x}}_0}{\Vert W^n{\mathbf{x}}_0\Vert }$$

$${\mathbf{x}}_n:={\varphi }^n({\mathbf{x}}_0)=\frac{W_n{W}_{n-1}\dots W_2{W}_1{\mathbf{x}}_0}{\Vert W_n{W}_{n-1}\dots W_2{W}_1{\mathbf{x}}_0\Vert }$$

$$J_{ij}=\frac{\partial }{\partial x_i}{\varphi }_j(\mathbf{x})$$

$$J\approx \frac{W}{\Vert W\mathbf{x}\Vert }\approx \frac{W}{\rho }.$$

$$\Lambda :=\underset{n\to \infty }{\lim }\frac{1}{n}\log \left(\prod _k^n{\rho }_k\right)$$

$${\mathbf{a}}_k=W{\mathbf{x}}_{k-1}+{\mathbf{u}}_k,$$

$$N_k=\Vert W{\mathbf{x}}_{k-1}+{\mathbf{u}}_k\Vert ,$$

$${\mathbf{x}}_n=M^{(n,0)}{\mathbf{x}}_0+\sum _{k=1}^n{M}^{(n,k)}{\mathbf{u}}_k$$

$$M^{(n,k)}:=\frac{W^{n-k}}{{\prod }_{l=k}^n{N}_l}.$$

$${\mathbf{x}}_n=\frac{W}{N_n}{\mathbf{x}}_{n-1}+\frac{1}{N_n}{\mathbf{u}}_n$$

$${\mathbf{x}}_n=\frac{W}{N_n}\left(\frac{W}{N_{n-1}}{\mathbf{x}}_{n-2}+\frac{1}{N_{n-1}}{\mathbf{u}}_{n-1}\right)+\frac{1}{N_n}{\mathbf{u}}_n$$

$$\Vert M^{(n,k)}\Vert =\frac{1}{{\prod }_{l=k}^n{N}_l}\Vert W^{n-k}\Vert \stackrel{n-k\gg 1}{\approx }\frac{1}{{\prod }_{l=k}^n{N}_l}{\rho }^{n-k}$$

$$N_l=\Vert W{\mathbf{x}}_{l-1}+{\mathbf{u}}_l\Vert =\rho +{\delta }_l$$

$$\Vert M^{(n,k)}\Vert \approx \frac{{\rho }^{n-k}}{(\rho +\delta {)}^{n-k}}$$

$${\mathcal{M}}_{\alpha }(k|n):={\left(\frac{\alpha }{\alpha +1}\right)}^{n-k}={\left(1-\frac{1}{\alpha +1}\right)}^{n-k}$$

$$\begin{array}{rl}\Delta \mathcal{M}(k|m,n):=& {\mathcal{M}}_{\alpha }(k|m)-{\mathcal{M}}_{\alpha }(k|n)={\left(\frac{\alpha }{\alpha +1}\right)}^{m-k}-{\left(\frac{\alpha }{\alpha +1}\right)}^{n-k}\\ =& \frac{(\alpha +1{)}^k}{{\alpha }^k}\frac{{\alpha }^n}{(\alpha +1{)}^n}\left(\frac{{\alpha }^{\tau }}{(\alpha +1{)}^{\tau }}-1\right)\\ =& {\left(1-\frac{1}{\alpha +1}\right)}^n{\left(1+\frac{1}{\alpha }\right)}^k\left({\left(1-\frac{1}{\alpha +1}\right)}^{\tau }-1\right)\le 0\end{array}$$

$$\begin{array}{rl}\Delta \mathcal{M}(k,l|n)=& \mathcal{M}(n-a|n)-\mathcal{M}(n-a-\delta |n)=\frac{{\alpha }^a}{(\alpha +1{)}^a}-\frac{{\alpha }^{a+\delta }}{(\alpha +1{)}^{a+\delta }}=\\ =& \frac{{\alpha }^a}{(\alpha +1{)}^a}\left(1-\frac{{\alpha }^{\delta }}{(\alpha +1{)}^{\delta }}\right)={\left(1-\frac{1}{\alpha +1}\right)}^a\left(1-{\left(1-\frac{1}{\alpha +1}\right)}^{\delta }\right)\ge 0\end{array}$$

$${\mathcal{M}}_L(k|n):={\rho }^{n-k}$$

$$\Delta {\mathcal{M}}_L(k|m,n)$$

$$\Delta {\mathcal{M}}_L(k,l|n)$$

$$\varphi (\mathbf{a}+\mathbf{b})=\frac{\varphi (\mathbf{a})+\varphi (\mathbf{b})}{1+\varphi (\mathbf{a})\varphi (\mathbf{b})}\approx \varphi (\mathbf{a})+\varphi (\mathbf{b})$$

$$\begin{array}{rl}{\mathbf{x}}_n& =\varphi (W{\mathbf{x}}_{n-1}+{\mathbf{u}}_n)\approx \varphi (W{\mathbf{x}}_{n-1})+\varphi ({\mathbf{u}}_n)\\ & =\varphi (W\varphi (W{\mathbf{x}}_{n-2}+u_{n-1}))+\varphi ({\mathbf{u}}_n)\approx \varphi (W\varphi (W{\mathbf{x}}_{n-2}))+\varphi (W\varphi ({\mathbf{u}}_{n-1}))+\varphi ({\mathbf{u}}_n)\end{array}$$

$$\Vert {\mathbf{x}}_n\Vert$$

$$S_{\rho }^{(n)}(\xi ):=\underset{n\text{ times}}{\underbrace{\varphi (\rho \cdot \dots \varphi (\rho \cdot \varphi (}}\xi )))$$

Full text

Introduction

Although the use of Recurrent Neural Networks in machine learning is boosting, also as effective building blocks for deep learning architectures, a comprehensive understanding of their working principles is still missing [1, 2]. Of particular relevance are Echo State Networks, introduced by Jaeger [3] and independently by Maass et al.[4] under the name of Liquid State Machine (LSM), which emerge from RNNs due to their training simplicity. The basic idea behind ESNs 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, ESNs are universal function approximators [5] and have shown to be effective in many relevant tasks [6, 7, 8, 9, 10, 11, 12].

These networks are known to be sensitive to the setting of hyper-parameters like the Spectral Radius (SR), the input scaling and the sparseness degree [3], 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 ESNs performance is usually maximized in correspondence of a very narrow region in hyper-parameter space called Edge of Criticality (EoC) [13, 14, 15, 16, 17, 18, 19, 20]. However, we comment that beyond such a region ESNs 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 ESNs is the Echo State Property (ESP), which basically asserts that their behavior should depend on the signal driving the network only, regardless of its initial conditions [21]. Despite being at the foundation of theoretical results [5], 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 [21, 22, 23].

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 ESNs [24, 25, 26, 27]. In particular, it is known that ESNs are characterized by a memory–nonlinearity trade-off [28, 29, 30], 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 [31].

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 [30], 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:

[TABLE]

where $\boldsymbol{x}_{k}\in\mathbb{R}^{N}$ is the system state at time-step $k$, $N$ is the number of hidden neurons composing the network reservoir, and $W\in\mathbb{R}^{N\times N}$ is the connection matrix of the reservoir. The signal value at time $k$, $\boldsymbol{s}_{k}\in\mathbb{R}^{N_{\text{in}}}$, is processed by the input-to-reservoir matrix $W_{\text{in}}\in\mathbb{R}^{N\times N_{\text{in}}}$. The activation function $\phi$ usually takes the form of the hyperbolic tangent function, for which the network is a universal function approximator [5]. Also linear networks (i.e., when $\phi$ 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[24, 25, 31, 27, 32]. Other activation functions have been proposed in the neural networks literature, including those that normalize the activation values on a closed hyper-surface[33] and those based on non-parametric estimation via composition of kernel functions[34].

The output $\boldsymbol{y}\in\mathbb{R}^{N_{\text{out}}}$ is generated by the matrix $W_{\text{out}}$, whose weights are learned: generally by ridge regression or lasso[3, 35] but also with online training mechanisms[36]. In ESNs, 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. 1

Self-normalizing activation function

Here, we propose a new model for ESNs 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

[TABLE]

and leaves the readout (1c) unaltered. The normalization in Eq. 2b projects the network pre-activation $a_{k}$ onto an $(N-1)$-dimensional hyper-sphere $\mathbb{S}^{N-1}_{r}:=\{\boldsymbol{p}\in\mathbb{R}^{N},\lVert\boldsymbol{p}\rVert=r\}$ of radius $r$. Fig. 2 illustrates the normalization operator applied to the state. Note that the operation (2b) 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 [37, 38] and some results on the universality of reservoir-based computation are given in[4, 39]. Recently, Grigoryeva and Ortega [5] proved that ESNs share the same property. Here, we show that the universal function approximation property also holds for the proposed ESNs model (2). To this end, we define a squashing function as a map $f:\mathbb{R}\to[-1,1]$ that is non decreasing and saturating, i.e., $\lim_{x\to\pm\infty}f(x)=\pm 1$. We note that $\phi_{i}(\boldsymbol{x}):=x_{i}/\lVert\boldsymbol{x}\rVert$ can be intended as a squashing function for each $i$th component. In the same work[5], the authors show that an ESN in the form (1) is a universal function approximator provided it satisfies the contractivity condition: $\lVert\phi(x_{1})-\phi(x_{2})\rVert\leq\lVert x_{1}-x_{2}\rVert$ for each $x_{1}$ and $x_{2}$, when $\phi$ is a squashing function. Jaeger and Haas[3] showed that this condition is sufficient to grant the ESP, implying that ESNs are universal approximators.

We prove in the Methods section that (2) grants the contractivity condition, provided that:

[TABLE]

where $\sigma_{\text{min}}(W)$ is the smallest singular value of matrix $W$ and $\lVert\boldsymbol{s}_{\text{max}}\rVert$ denotes the largest norm associated to an input. Eq. 3 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 (3) is satisfied. From (2), the system state dynamics in the autonomous case reads:

[TABLE]

At time-step $n$ the system state is given by

[TABLE]

By iterating this procedure, one obtains:

[TABLE]

where $\boldsymbol{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:

[TABLE]

Furthermore, note that a system described by matrix $W$ and a system characterized by $W^{\prime}=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 ESNs, one usually tries to bring the system close to the EoC, since it is in that region that their performance is optimal[14]. 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 (2) 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 (2b) and then show that, since its spectral radius tends to $1$, the maximum Local Lyapunov Exponent (LLE) must be null. The Jacobian matrix of (2b) reads:

[TABLE]

where the time index $k$ is omitted to ease the notation. We observe that, asymptotically for large networks ($N\to\infty$), we have that $a_{i}/\lVert\boldsymbol{a}\rVert\to 0$ for each $i$, meaning that the Jacobian matrix reduces to $J(\boldsymbol{x})=W/\lVert W\boldsymbol{x}\rVert$. As we are considering the case with $r=1$, we know that $\lVert\boldsymbol{x}\rVert=1$.

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

[TABLE]

Under this approximation (9), the largest eigenvalue of $J$ must be $1$ as the SR $\rho$ is the largest absolute value among the eigenvalues of $W$. We thus characterize the global behavior of (2b) by considering the maximum LLE[14], which is defined as:

[TABLE]

where $\rho_{k}$ is the spectral radius of the Jacobian at time-step $k$. Eq. 10 implies that $\Lambda=0$ as $n\rightarrow\infty$, hence proving our claim.

In order to demonstrate that $\Lambda=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 (8). The results are displayed in Fig. 3. Fig. 3 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. 3 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 (2b) 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 $\boldsymbol{u}:=W_{\text{in}}\boldsymbol{s}$. Accordingly, Eq. 2a takes on the following form:

[TABLE]

where $\boldsymbol{u}_{k}$ operates as a time-dependent bias vector. Let us define the normalization factor as:

[TABLE]

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

[TABLE]

where

[TABLE]

By looking at (13), it is possible to note that each $\boldsymbol{u}_{k}$ is multiplied by a time-dependent matrix, i.e., the network’s final state $\boldsymbol{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 (2) 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 $\boldsymbol{x}_{n}$. Our proposal shares similarities with a memory measure called Fisher memory, first proposed by Ganguli et al. [27] and then further developed in [24, 32]. However, our measure can be easily applied also to non-linear networks like the one Eq. 2 proposed here, justifying the analysis discussed in the following.

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

[TABLE]

All past history of the signal is processed in $\boldsymbol{x}_{n-1}$. Note that (15) keeps the same form for every $n$. Going backward in time one more step, we obtain:

[TABLE]

As $\lVert\boldsymbol{x}_{n}\rVert=\lVert\boldsymbol{x}_{n-1}\rVert=\lVert\boldsymbol{x}_{n-2}\rVert=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 (14) can be written as

[TABLE]

where the approximation holds thanks to the Gelfand’s formula111The Gelfand’s formula [40] states that for any matrix norm $\rho(A)=\lim_{k\to\infty}\left\|A^{k}\right\|^{\frac{1}{k}}$.. If the input is null, we expect each $N_{l}$ to be of the order of $\rho$ as $\lVert x_{l}\rVert=1$. So we write:

[TABLE]

where $\delta_{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 $\rho$, while the input modifies the state and so the normalization value will be modified accordingly. The value $\delta_{l}$ is introduced to account for this fact: $\delta_{:}=N_{l}-\rho$.

If we assume that each input produces a similar effect on the state (i.e. $\delta_{l}=\delta$ for every $l$), we finally obtain:

[TABLE]

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 $\mathcal{M}_{\alpha}(k|n)$ as the influence of input at time $k$ on the network state at time $n$. More formally, having defined $\alpha:=\rho/\delta$ (since (19) only depends on this ratio), we use (19) to define the memory as:

[TABLE]

Eq. 20 goes to [math] (i.e., no impact of the input on the states) as $\alpha\to 0$ and tends to $1$ for $\alpha\rightarrow\infty$, implying that for an infinitely large SR the network perfectly remembers its past inputs, regardless of their occurrence in the past. Note that (20) does not depend on $k$ and $n$ individually, but only on their difference (elapsed time): the larger the difference, the lower the value of $\mathcal{M}_{\alpha}(k|n)$, meaning that far away inputs have less impact on the current state.

Memory loss

By using (20), we define the memory loss between state $\boldsymbol{x}_{n}$ and $\boldsymbol{x}_{m}$ of the signal at time-step $k$ (with $m>n>k$ and $m=n+\tau$) as follows:

[TABLE]

For very large values of $\alpha$, we have $\Delta\mathcal{M}(k|m,n)\to 0$, meaning that in our model (2b) 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+\delta$. Define:

[TABLE]

We have that $\lim_{\delta\to\infty}\Delta\mathcal{M}(n-a,n-a-\delta|n)=\left(1-\frac{1}{\alpha+1}\right)^{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<\infty$ 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 RNNs (i.e., $\boldsymbol{x}_{n}=\boldsymbol{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:

[TABLE]

It is worth noting that there is no dependency on $\delta$ and, therefore, on the input. Just like before, we have the memory loss of the signal at time-step $k$ between state $\boldsymbol{x}_{n}$ and $\boldsymbol{x}_{m}$, as:

[TABLE]

where we set $m>n>k$ and $m=n+\tau$. 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+\delta$. As such:

[TABLE]

We see that in both (24) and (25) 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 ESNs, i.e., those using a $\tanh$ activation function. Define $\phi:=\tanh$ (applied element-wise). Then, for generic scalars $a$ and $b$, when $|b|\ll 1$ the following approximation holds:

[TABLE]

When $\boldsymbol{a}$ is the state and $\boldsymbol{b}$ is the input, our approximation can be understood as a small-input approximation [41]. Then, the state-update reads:

[TABLE]

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

[TABLE]

where we defined:

[TABLE]

We see that, differently from the previous cases, the final state $\boldsymbol{x}_{n}$ is a sum of nonlinear functions of the signal $\boldsymbol{u}_{k}$. Each signal element $\boldsymbol{u}_{k}$ is encoded in the network state by first multiplying it by $\rho$ and then passing the outcome through the nonlinearity $\phi(\cdot)$. 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 $\phi(\cdot)$). 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. 4. 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 $\boldsymbol{y}_{k}=\boldsymbol{u}_{k-\tau}$, where $\tau$ ranges from [math] 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 $\rho=0.95$. For the proposed model (2), the input scaling is chosen to be $0.01$, while the SR is $\rho=15$. For the sake of simplicity, in what follows we refer to ESNs resulting from the use of (2) as “spherical reservoir”.

To evaluate the performance of the network, we use the accuracy metric defined as $\gamma=\max\{1-\text{NRMSE},0\}$, where the Normalized Root Mean Squared Error (NRMSE) is defined as:

[TABLE]

Here, $\boldsymbol{\hat{y}}_{k}$ denotes the computed output and $\langle\cdot\rangle$ represents the time average over all time-steps taken into account. In the following, we report the accuracy $\gamma$ while varying $\tau$ 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_{\text{train}}=5000$ and the associated performance is evaluated using a test set of length $L_{\text{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.5, 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. (20), (23), and (28)) shown in Fig. 4.

Multiple superimposed oscillators

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

[TABLE]

Results are shown in Fig.5, panel (b). This task is difficult because of the multiple time scales characterizing the inputs[3]. 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\approx 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.

[TABLE]

It is well-known that this system is chaotic when $\sigma=10$, $\beta=\frac{8}{3}$ and $\rho=28$. We simulated the system with these parameters and then fed the network with the $x$ coordinate only. Results are shown in Fig.5, 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:

[TABLE]

In our experiments, we simulated the system using standard parameters, that is, $\lambda=17$, $\alpha=0.2$ and $\beta=0.1$. Results are shown in Fig.5, 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.5, 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 [30]. The network is fed with a signal where the $u_{k}$s 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(\nu\cdot u_{k-\tau})$. We see that $\tau$ accounts for the memory required to correctly solve the task, while $\nu$ controls the amount non-linearity involved in the computation. For each configuration of $\tau$ and $\nu$ 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]$. 222This choice of exploring large values of SR for the hyperspherical case is motivated by the fact that condition (3) must be satisfied in order for the network to work properly: choosing a large SR will also lead to larger $\sigma_{\text{min}}(W)$, as discussed in the ‘Universal function approximation property’ section. We then choose the hyper-parameters configuration that minimizes (30) on a training set of length $L_{\text{train}}=500$ and then assess the error on a test set with $L_{\text{test}}=200$.

In Figures 7, 8, and 9, we show the NRMSE for the task described above for different values of $\nu$ and $\tau$, 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 [30] 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 $\tau$ and $\nu$ are large. Predictions obtained for a specific instance of this task requiring both features are given in Fig 6, showing how the proposed model outperforms the competitors. In order to explore the memory-nonlinearity relationship we followed the experimental design proposed in [30]: 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 [30, 42] 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 ESNs 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 ESNs 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 EoCs. 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 ESNs.

The proposed activation function allows the model to display a memory capacity comparable with the one of linear ESNs, 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. [27], 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 ESNs. 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 ESNs 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 [5], can be proved for an ESN of the form (1b) 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 $\boldsymbol{x}$ and $\boldsymbol{y}$ in the domain of the activation function $\phi$, it must be true that:

[TABLE]

Let us introduce the the notation $\boldsymbol{\hat{x}}:=\boldsymbol{x}/\lVert\boldsymbol{x}\rVert$ and $\boldsymbol{\hat{y}}:=\boldsymbol{y}/\lVert\boldsymbol{y}\rVert$. Taking the square of the norm, one gets:

[TABLE]

Where $\boldsymbol{a}\cdot\boldsymbol{b}$ is the scalar product between two vectors and the inequality $\frac{a}{b}+\frac{b}{a}\geq 2$ follows from the fact that $(a-b)^{2}=a^{2}+b^{2}-2ab>0$ for all $a,b>0$. Now, we assume $\lVert\boldsymbol{x}\rVert,\lVert\boldsymbol{y}\rVert>r$ and show that:

[TABLE]

proving the contractivity condition (34).

We see that the only condition needed is $\lVert\boldsymbol{x}\rVert,\lVert\boldsymbol{y}\rVert>r$, which means that the linear part of the update (2a) must map states outside the hyper-sphere of radius $r$. Finally, by applying properties of norms, we observe that:

[TABLE]

and asking this to be larger than $r$ results in the condition (3).

Input-driven dynamics

Here we show how to derive (13). Consider a dynamical system evolving according to (2): if we explicitly expand the first steps from the initial state $\boldsymbol{x}_{0}$ we obtain:

[TABLE]

So that the general case can be written as:

[TABLE]

and in a more compact form as follows:

[TABLE]

where

[TABLE]

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.

bESN binary ESN EoC Edge of Criticality ESN Echo State Network ESP Echo State Property FIM Fisher Information Matrix LLE Local Lyapunov Exponent LSM Liquid State Machine MFT Mean Field Theory ML Machine Learning MSO Multiple Superimposed Oscillator NRMSE Normalized Root Mean Squared Error PDF Probability Density Function RBN Random Boolean Network RNN Recurrent Neural Network RP Recurrency Plot SR Spectral Radius