Takens-inspired neuromorphic processor: a downsizing tool for random recurrent neural networks via feature extraction

TL;DR

This paper introduces a Takens-inspired method to reduce the size of random recurrent neural networks for time series prediction by leveraging feature extraction and attractor reconstruction, resulting in a more efficient hybrid processor.

Contribution

It presents a novel technique combining Takens-based attractor reconstruction with machine learning to minimize neural network size for time series prediction.

Findings

Reduced neural network size by a factor of 15.

Demonstrated stabilization of an arrhythmic neural model.

Derived scaling laws for prediction error.

Abstract

We describe a new technique which minimizes the amount of neurons in the hidden layer of a random recurrent neural network (rRNN) for time series prediction. Merging Takens-based attractor reconstruction methods with machine learning, we identify a mechanism for feature extraction that can be leveraged to lower the network size. We obtain criteria specific to the particular prediction task and derive the scaling law of the prediction error. The consequences of our theory are demonstrated by designing a Takens-inspired hybrid processor, which extends a rRNN with a priori designed delay external memory. Our hybrid architecture is therefore designed including both, real and virtual nodes. Via this symbiosis, we show performance of the hybrid processor by stabilizing an arrhythmic neural model. Thanks to our obtained design rules, we can reduce the stabilizing neural network's size by a factor of 15 with respect to a standard system.

Equations

$${\mathbf{x}}_{n+1}=f_{NL}(\mu W\cdot {\mathbf{x}}_n+\alpha W^{IN}\cdot y_{n+1}^{IN}+W^{off}\cdot b),$$

$$y_{n+1}^{out}=W^{out}\cdot {\mathbf{x}}_{n+1}.$$

$$W_{op}^{out}=\underset{W^{out}}{\min }\Vert W^{out}\cdot {\mathbf{x}}_{n+1}-y_{n+1}^T\Vert ,$$

$$y_{n+1}=y{}_n+\delta \left(\frac{\vartheta y_{{\tau }_m}}{1+(y_{{\tau }_m}{)}^{\nu }}-\psi y_n\right),$$

$$(1-{ϵ}_1)\Vert {\mathbf{y}}_{n+1}^{IN}-{\mathbf{y}}_n^{IN}\Vert \le \Vert \phi ({\mathbf{y}}_{n+1}^{IN})-\phi ({\mathbf{y}}_n^{IN})\Vert \le (1+{ϵ}_2)\Vert {\mathbf{y}}_{n+1}^{IN}-{\mathbf{y}}_n^{IN}\Vert .$$

$${\text{x}}_{n+1}=f_{NL}(\mu W\cdot {\text{x}}_n+\alpha W^{IN}\cdot y_{n+1}^{IN}+W^{off}\cdot b),$$

$$y_{n+1}^{out}=W^{out}\cdot ({\text{x}}_{n+1},{\text{x}}_{n+1+{\tau }_T})$$

$$ϵ\frac{dv(t)}{dt}=v(t)[v(t)-g][1-v(t)]-w+I+\xi (t),$$

$$\frac{dw(t)}{dt}=v(t)-Dw(t)-H,$$

$${\mathbf{y}}_n^{IN}=(y_n^{IN},y_{n+{\tau }_0}^{IN},\dots ,y_{n+(M-1){\tau }_0}^{IN}),$$

$${\mathbf{y}}_{n+1}^{IN}=(y_{n+1}^{IN},y_{(n+1)+{\tau }_0}^{IN},\dots ,y_{(n+1)+(M-1){\tau }_0}^{IN}).$$

$$\phi ({\mathbf{y}}_n^{IN})=[{\phi }_{l_1}({\mathbf{y}}_n^{IN}),{\phi }_{l_2}({\mathbf{y}}_n^{IN}),\dots ,{\phi }_{l_h}({\mathbf{y}}_n^{IN})],$$

$$\phi ({\mathbf{y}}_{n+1}^{IN})=[{\phi }_{l_1}({\mathbf{y}}_{n+1}^{IN}),{\phi }_{l_2}({\mathbf{y}}_{n+1}^{IN}),\dots ,{\phi }_{l_h}({\mathbf{y}}_{n+1}^{IN})];$$

$$\frac{\Vert \phi ({\mathbf{y}}_{n+1}^{IN})-\phi ({\mathbf{y}}_n^{IN})\Vert }{\Vert {\mathbf{y}}_{n+1}^{IN}-{\mathbf{y}}_n^{IN}\Vert }\in [(1-{ϵ}_1),(1+{ϵ}_2)].$$

$$\frac{\Vert \phi ({\mathbf{y}}_{n+1}^{IN})-\phi ({\mathbf{y}}_n^{IN}){\Vert }_{min}}{\Vert {\mathbf{y}}_{n+1}^{IN}-{\mathbf{y}}_n^{IN}\Vert }={ϵ}_{min};$$

$$\frac{\Vert \phi ({\mathbf{y}}_{n+1}^{IN})-\phi ({\mathbf{y}}_n^{IN}){\Vert }_{max}}{\Vert {\mathbf{y}}_{n+1}^{IN}-{\mathbf{y}}_n^{IN}\Vert }={ϵ}_{max},$$

$$\Vert \phi ({\mathbf{y}}_{n+1}^{IN})-\phi ({\mathbf{y}}_n^{IN}){\Vert }_{min}=$$

$${\left(\sum _{l_g=l_{min}}^{l_{max}}[{\phi }_{l_g}({\mathbf{y}}_{n+1}^{IN})-{\phi }_{l_g}({\mathbf{y}}_n^{IN}){]}_{min}^2\right)}^{1/2},$$

$$\Vert \phi ({\mathbf{y}}_{n+1}^{IN})-\phi ({\mathbf{y}}_n^{IN}){\Vert }_{max}=$$

$${\left(\sum _{l_g=l_{min}}^{l_{max}}[{\phi }_{l_g}({\mathbf{y}}_{n+1}^{IN})-{\phi }_{l_g}({\mathbf{y}}_n^{IN}){]}_{max}^2\right)}^{1/2},$$

$$\Vert {\mathbf{y}}_{n+1}^{IN}-{\mathbf{y}}_n^{IN}\Vert =\sqrt{(y_{n+1}^{IN}-y_n^{IN}{)}^2+\cdots +(y_{(n+1)+(M-1){\tau }_0}^{IN}-y_{n+(M-1){\tau }_0}^{IN}{)}^2}.$$

Full text

PACS numbers

May be entered using the \pacs{#1} command.

††preprint: APS/123-QED

I Introduction

Artificial neural networks (ANNs) are systems prominently used in computational science as well as investigations of biological neural systems. In biology, of particular interest are recurrent neural networks (RNNs), whose structure can be compared among others with nervous system’s networks of advanced biological species Maass (2016). In computation, RNNs have been used to solve highly complex tasks which pose problems to other classical computational approaches Jaeger and Haas (2004); Graves et al. (2009, 2013); Sak et al. (2014); Li and Wu (2015). Their recurrent architecture allows the generation of internal dynamics, and consequently RNNs can be studied utilizing principles of dynamical systems theory. Therefore, the network’s nonlinear dynamical properties are of major importance to its information processing capacity. In fact, optimal computational performances are often achieved in a stable equilibrium network’s state, yet near criticality Langton (1990); Natschlager et al. (2005).

Among the recurrent networks reported in current literature, rRNNs are popular models for investigating fundamental principles of information processing. In these models the synaptic neural links are randomly weighted, typically following a uniform Marquez et al. (2018) or Gaussian distribution Bruckstein et al. (2009); Ganguli and Sompolinsky (2012); Babadi and Sompolinsky (2004); Schottdorf et al. (2015). Recently, there is an increasing interest in some particular types of random recurrent networks with a simplified design, where just the output layers are trained using a supervised learning rule. Such rRNNs are typically referred to as Reservoir Computers, which include echo state networks Jaeger and Haas (2004) and also liquid state machines Maass et al. (2002). Reservoir Computing provides state-of-the-art performance for challenging problems like high-quality long term prediction of chaotic signals Jaeger and Haas (2004); Antonik et al. (2017).

Prediction corresponds to estimating the future developments of a system based on knowledge about its past. Chaotic time series prediction is of importance to a large variety of fields, including the forecasting of weather Lorenz (1963), the evolution of some human pathologies Mackey and Glass (1977), population density growth Ott (1993), or dynamical control as found in the regulation of chaotic physiological functions FitzHugh (1955, 1961); Nagumo et al. (1962). In order to build a predictor for chaotic systems, most common techniques can be divided into the following groups Weigend and Gershenfeld (1993): (i) linear and nonlinear regression models such as Autoregressive-Moving-Average, Multi Adaptive Regression Spline Zarandi et al. (2013), and Support Vector Machine Celikyilmaz and Turksen (2007); (ii) state space based techniques for prediction of continuous-time chaotic systems, which utilize attractor reconstruction and interactions between internal degrees of freedom to infer the future Kantz and Schreiber (1997); Weigend and Gershenfeld (1993); Farmer and Sidorowich (1987); Alparslan et al. (1998); Kugiumtzis et al. (1998). Attractor reconstruction method is based on the embedding of the original state space in a delay-coordinate space Ott et al. (1990). And (iii) the connectionist approach, including recurrent and feedforward Rojas (1996); Gurney (1997), deep LeCun et al. (2015), and convolutional ANNs LeCun et al. (1998). This approach usually comprehend the design of ANNs using large amounts of neurons to process information Goodfellow et al. (2016).

The high-dimensionality of the ANNs’ hidden layer is commonly translated in a computationally expensive problem when considering the optimization of such networks to solve a task. Elseways, in the reservoir computing approach such training efforts are reduced due to the training is done on the output layer only. However, the more neurons in the hidden layer which are connected to the output layer, the higher the computational cost of the training step. Introducing a novel methodology, we develop a state space-based concept which guide the downsizing of the rRNNs’ hidden layer. To achieve this objective, we describe rRNNs and state space-based models within the same framework. From where we show state space patterns revealed by spontaneous reconstructions inside the high dimensional space of our random recurrent network. Furthermore, we introduce a novel methodology based on the Takens embedding theorem to identify the embedding dimensions of the input system’s spontaneous reconstruction, and their relevance to the system’s prediction performance.

We immediately exploit our insight and devise a new hybrid Takens-inspired ANN concept, in which a network is extended by an a-priori designed delay external memory. The delay term is used to virtually extend the size of the network by introducing virtual nodes which exist in the delay path. We use this new design to first validate our interpretation, and then devise an advanced hybrid rRNN to stabilize a non-periodic neuronal model which requires $15$ times less neurons than a bench-mark rRNN. As this system is driven by a stochastic signal, we show how our approach can leverage properties of the underlying deterministic system even for the case of a stochastic drive.

II Random Recurrent Networks for prediction

A rRNN is illustrated in Fig. 1, indicating the temporal flow of information received by each neuron or node. Nodes are represented by $\bigoplus$. The rRNN consists of a reservoir of $m$ nodes in state $\mathbf{x}_{n}$ at integer time $n$. Nodes are connected through random, uniformly distributed internal weights defined as coefficients in the matrix $W$ of dimensionality $m\times m$. The resulting randomly connected network is injected with 1D input data $y^{IN}_{n+1}$ according to input weights defined as random coefficients in the vector $W^{IN}$ of dimensionality $m\times 1$. The time-discrete equation that governs the network is Jaeger and Haas (2004)

[TABLE]

where $\mu$ is the bifurcation parameter, $\alpha$ the input gain, $f_{NL}(\cdot)$ is a nonlinear sigmoid-like activation function, and $b$ the constant phase offset injected through offset weights, defined as random coefficients in the vector $W^{off}$ of dimensionality $m\times 1$. In practice, we construct a network with $m=1000$, using the Matlab routine random. Connection weights $W_{i,j}$ are distributed around zero. In order to set our recurrent network in its steady state, $W$ is normalized by its largest eigenvalue. The network’s connectivity is set to one, hence it is fully connected.

An output layer creates the solution $y^{out}$ to the prediction task. In this step the network approximates the underlying deterministic law that rules the evolution of the input system. The output layer provides the computational result according to

[TABLE]

The output weights vector $W^{out}$ is calculated according to a supervised learning rule, using a reference teacher/target signal $y^{T}_{n+1}$ Jaeger (2001). We calculate the optimal output weights vector $W^{out}_{op}$ by

[TABLE]

via its pseudo-inverse (using singular value decomposition) with the Matlab routine pinv. Equation (3) therefore minimizes the error between output $W^{out}\cdot\mathbf{x}_{n+1}$ and teacher $y^{T}_{n+1}$. As training error measure we use the normalized mean squared error (NMSE) between output $y^{out}_{n+1}$ and target signal $y^{T}_{n+1}$, normalized by the standard deviation of teacher signal $y^{T}_{n+1}$.

When we use a rRNN for time series prediction of a chaotic oscillator such as the Mackey-Glass (MG) time-delayed system Mackey and Glass (1977). We can achieve good long-term prediction performances using a network of 1000 neurons. Here, long-term predictions are defined as predictions far beyond one step in the future. The task is to predict future steps of the chaotic MG system in its discrete-time version:

[TABLE]

where $y_{\tau_{m}}=y_{(n-\tau_{m}/\delta)}$, $\tau_{m}=17$ as the time delay, and $\delta=1/10$ is the stepsize indicating that the time series is subsampled by 10. The other parameters are set to $\vartheta=0.2,\nu=10,\psi=0.1$. For any prediction task in this paper, we consider 20-network model via different initializations of $\{W,W^{IN},W^{off}\}$. The prediction horizon is estimated to be 300 time step, defined by the inverse of the largest Lyapunov exponent of the MG system ($\lambda_{max}^{MG}=0.0036$). For such prediction horizon, and for predicting the value of 20 different MG sequences, we obtain the average of all NMSEs, resulting in $0.091\pm 0.013$. This performance was obtained for $b=0.2$, $\alpha=0.8$ and $\mu=1.1$, which was found to offer the best prediction performance in a range of $\mu\in[0.1,1.3]$. Moreover, the network was trained with $3000$ values of the MG system, with a teacher signal given by $y^{T}_{n+1}=y^{IN}_{n+1}$. We subtracted the average of the MG time series before injection into the rRNN, which is a common practice Jaeger (2001). Then we discarded the first $1000$ points of its response to avoid initial transients. Right after training, where we determined $W_{out}$, we connected $y^{out}_{n+1}$ to $y^{IN}_{n+1}$, and we left the network running freely 300 time steps, indicated by the prediction horizon.

Given that good long-term prediction performances are obtained, we wonder how input information is represented in the network’s space. Such qualitative investigations are for instance common practices in image classification tasks, where the extracted features are often identified and illustrated together with the hidden layers that have generated them Taigman et al. (2014). In the following, we introduce a technique that allow us to identify feature representations of the input information in the rRNN’s high dimensional space, which are linked to good prediction performances.

III A method for feature extraction in random recurrent networks

As shown previously, our rRNN is able to predict the future values of 1D input chaotic data $y^{IN}$. Such kind of data comes from a continuous-time chaotic system. As it is known in chaos theory, a minimum of three dimensions are required in a continuous-time nonlinear system to generate chaotic solutions. Typically, continuous-time chaotic solutions come from models consisting of a system of at least three nonlinear ordinary differential equations (ODEs). However, there are other ways to obtain such chaotic dynamics. One of those ways includes the introduction of an explicit temporal variable to an ODE, such as a time delay with respect to the main temporal variable. We define these models as delay differential equations (DDEs), where a DDE is in fact equivalent to an infinite-dimensional system of differential equations Demidenko et al. (2009). The number of solutions to a DDE is in theory infinite due to the infinite amount of initial conditions in the continuous rank required to solve the equation. Each initial condition initializes an ODE from the infinite-dimensional system of equations. Thus, the introduction of a time delay in an ODE, resulting in a DDE, provides sufficient dimensionality to allow for the existence of chaotic solutions. For instance, this is how the MG system can develop chaotic solutions. In such case, one just have access to a time series from a single accessible variable $y_{n+1}$ in Eq. (4), while all others remain hidden. Nevertheless, hidden variables are participating in the development of the global dynamics as well as the accessible variables.

In order to approximate a full dimensional representation of these oscillators, we could embed the 1D sequence $y_{n+1}$ into a high dimensional space. Consequently, reconstructing its state-space. Among the most practiced methods to embed 1D information, we highlight state space reconstruction techniques such as delay reconstruction (Whitney and Takens embedding theorems Whitney (1936); Takens (1981)) or through the Hilbert transform Pikovsky et al. (2001).

Delay reconstruction is a widely used method to complete missing state space information. According to the Takens embedding theorem, the time delayed version of a time series suffices to reveal the structure of the state space trajectory. Let us represent the data in a $M$-dimensional space by the vectors $\mathbf{y}_{n}=[y_{n},\;y_{(n+\tau_{0})},\;\ldots,\;y_{(n+(M-1)\tau_{0})}]^{\dagger}$, where $[\cdot]^{\dagger}$ is a transpose matrix and $y_{n}$ is the original time series. The essential parameters for state-space reconstruction are $M$ and time delay $\tau_{0}$. Embedding delay $\tau_{0}$ is often estimated by applying autocorrelation analysis or time delayed mutual information to the signal $y_{n}$. The temporal position of the first zero Kantz and Schreiber (1997) of either method maximizes the possibility to extract additional information which contains independent observations from the original signal, and hence obtain trajectories or dynamic motion along potentially orthogonal state-space dimensions. If successful, information of such maximized linear independence can enable the inference of the missing degrees of freedom Eckmann and Ruelle (1985); Packard et al. (1980); Takens (1981). On the other hand, we estimate the minimum amount of required embedding dimensions $M$ by using the method of false nearest neighbors Kantz and Schreiber (1997). Crucially, the following concepts are not restricted to a particular method of determining $\tau_{0}$ or $M$.

The autocorrelation function (ACF) is often employed to identify the temporal position $\tau_{0}$ used to reconstruct missing coordinates that define any continuous-time dynamical system. In order to show how Takens-based attractor reconstruction performs, we use a dataset of $1\times 10^{4}$ values provided by Eq. (4). The outcome of the autocorrelation analysis reveals that the ACF has its first minimum at lag $\tau_{0}=-12$. In Fig. 2 we show the absolute value of the ACF together with three examples of 2D delay-reconstructions for three different values of the embedding lag $\tau_{0}$. As it can be seen, $2$D-reconstructions based on lags $\tau_{0}=-7$ and $-17$ also unfold the geometrical object within a state space. However, aiming at a maximally orthogonal embedding, we base our analysis on the attractors reconstructed by exactly the first minimum of the ACF, $\tau_{0}=-12$. According to the false nearest neighbor analysis, the minimum dimensions $M$ required to reconstruct the MG attractor is $4$. The Takens scheme therefore provides a set of coordinates $\mathbf{y}_{n}=\{y_{n},y_{n-12},y_{n-24},y_{n-36}\}$ which reconstructs the state space object.

As it can be seen, this is a classic method to extract feature representations $\mathbf{y}_{n}$ from the original and unknown feature $y_{n}$. Next, we want to identify comparable attractor reconstruction inside our rRNN’s state space.

III.1 A Takens-inspired feature extraction technique

The network’s state space is defined by the set of orthogonal vectors in $W$. For that, we analyze the injected signal representation $y^{IN}_{n+1}$ via network’s node responses. As input information $y^{IN}_{n+1}$ is being randomly injected into the network’s high dimensional space, we search for possible spatial representations of the 1D input, where network nodes serve as embedding dimensions.

With that goal in mind, we proceed with an analysis comparable to the ACF used in delay embedding, based on the maximum absolute value of the cross-correlation function $|CC(x^{i},y^{IN})|_{max}$ between all node responses $\{x^{i}\}$ and the input data $y^{IN}$. As our aim is to identify nodes providing observations approximately orthogonal to $y^{IN}$, we record for every network node such cross-correlation maxima and the temporal position, or lag, of this maxima. Figure 3(a) shows the cross-correlation analysis (CCAs) for $\mu=1.1$, where node correlation lags and maxima correspond to the abscissa and ordinate, respectively. In this case, a distribution of node lags, where an extension of the range to $\{l_{min},l_{max}\}=\{-49,45\}$ is shown. Thus, the rRNN’s nodes reveal a strong cross-correlation at time lags covering all Takens embedding delays $\{0,\tau_{0},2\tau_{0},$ $3\tau_{0}\}$. Nodes lagged around $\{-36,-24,-12,0\}$ should well approximate the Takens embedded attractor, and a state space representation of the input sequence is presented within the rRNNs’ space.

In Fig. 3(b), we illustrate the numerous possible extracted features from the originally injected attractor embedded in the rRNN’s space for $\mu=1.1$. The first column of the figure shows simple $2$D-projections of the delay-reconstructed attractor by using the set of lags $\{-24,-18,-12,-6\}$. The attractors reconstructed by network nodes lagged at $\{-24,-18,-12,-6\}$ are shown in the three next columns. Attractors reconstructed by network nodes with maximum cross-correlation values at lags $\{-18,-6\}$ do not belong to the set of original Takens delay-coordinates. We included these additional delay-dimensions to better illustrate the breadth of features provided by the rRNN.

As shown above, the CCA also finds node lags and correlations that do not agree with the Takens framework. In the following, we introduce a methodology to exclude such additional delay-dimensions from our predictor. As an starting point, we suppress nodes with specific CCA-lag positions during the training step of the output layer, such that they will not be available to the readout matrix but still take part in the rRNN’s state evolution. To that end, we select nodes for which their CCA-lag positions are within windows of width $\delta\tau_{0}^{net}$, centered at integer multiples of $\tau_{0}^{net}$, where $\tau_{0}^{net}$ represents the time lag used for delay reconstructions in the Takens’ scheme. The windows width $\delta\tau_{0}^{net}$ defines a time-lag uncertainty associated to the identified rRNNs’ delay coordinates. All nodes with CCA-lag positions not inside the set of ($n\tau_{0}^{net}\pm\delta\tau_{0}^{net})$, $n\in\mathbb{Z}$ will not be available to the readout layer.

To illustrate our method and its effect, we show the non-filtered CCA of the rRNN when driven by the MG signal and with bifurcation parameter $\mu=1.1$ in Fig. 4(a). Using a constant CCA-windows width of $\delta\tau_{0}^{net}=3$, as the minimum uncertainty found associated with good performances, we now scan the position of these CCA-windows by changing $\tau_{0}^{net}$. We restrict the number of windows to $n\in\{-M,-M+1,\dots,M\}$. For $\tau_{0}^{net}\in\{-17,-12,-7\}$, in panels (b-d) of Fig. 4, we show examples of filtered CCAs where just rRNN nodes available for the readout layer are present.

Based on such movable CCA-filters, we can estimate the relevance of different CCA lags on the rRNN’s capacity to predict a particular temporal sequence. We define 20-network model via different initializations of $\{W,W^{IN},W^{off}\}$, and for each we obtain the NMSE for predicting the value of 20 different MG sequences at 300 time steps into the future. In Fig. 4(e), we show the resulting NMSE, averaging over all system combinations and for $-25\leq\tau_{0}^{net}\leq-1$. Here, the average NMSE is given by the solid line, and the standard deviation (stdev) interval by the dashed black lines. In panel (f) we show the percentage of rRNN’s for which prediction diverged from the target, i.e. NMSE$>1$. Panel (g) shows how many nodes are available to the network’s output layer. The constant red curves present in all panels show averaged performances obtained for the non-filtered rRNN, namely the classical RC, again with the stdev interval given by the dashed red lines.

Restricting the system’s output based on the CCA-filter windows has a strong and systematic impact onto the rRNN’s prediction performance. Performance is optimized for very characteristic filter positions, i.e. for $\tau_{0}^{net}\in[-16,-5]$. For $\tau_{0}^{net}\in[-12,-11]$ the embedding available to the network’s output closely corresponds with the lags of the original Takens attractor embedding. The performance achieved by setting $\tau_{0}^{net}\simeq\tau_{0}$ in our approach reduces the NMSE more than 1-fold, even though the system has in average significantly less network nodes available to the output layer ($\sim 500$), see Fig. 4(g). For $\tau_{0}^{net}=-7$, the performance is higher by one order of magnitude, accessed by using approximately the same initial set of nodes available to the output ($\sim 1000$), see Fig. 4(g). We are therefore able to identify node families that shows attractor embedding features in the network’s space for the first time, based on their CCA-lag.

In summary, it can be seen that for filters based on $\tau_{0}^{net}\in[-16,-5]$ the prediction performances are either coincident or better than the non-filtered CCA case. This result is in agreement with Fig. 2, where we show that lags in such interval seem to still unfold the object in the state space.

Finally, some further aspects are seen in our data. For $\tau_{0}^{net}>-4$ the CCA-filter windows overlap and nodes with a lag inside such positions of overlap are assigned to multiple windows. This artificially increases the number of nodes available to the system’s output beyond $m=1000$. The resulting NMSE strongly increases beyond the one for the original rRNN. We attribute this characteristic to over-fitting during learning, where there are just repetitions of the same few delay-coordinates.

III.2 Characteristics of the extracted features

Our previous analysis identifies and harness meaningful features related to good prediction performances. As this was realized by randomly connected networks, attractor reconstruction was achieved by randomly mapping the originally 1D-input onto the high-dimensional rRNN’s space. Such random mapping is treated by the framework of random projections theory (RPT) Johnson and Lindenstrauss (1984); Indyk and Motwani (1998); Frankl and Maehara (1988); Dasgupta and Gupta (2002); Sivakumar (2002). An active area of study within this field treats the question if original input data is randomly mapped onto the dimensions of the projection space, the structural damage to the original object is minimized.

To determine the degree of potential structural distortions to the original input attractor after random mapping onto the network’s high-dimensional space according to $\mathbf{y}^{IN}_{n}\rightarrow\varphi(\mathbf{y}^{IN}_{n})$, we measure distances between consecutive states (interstate distances) of the original $||\mathbf{y}_{n+1}^{IN}-\mathbf{y}_{n}^{IN}||$ and the projected $||\varphi(\mathbf{y}_{n+1}^{IN})-\varphi(\mathbf{y}_{n}^{IN})||$ objects. Following these steps, we take inspiration from RPT extended to nonlinear mapping Tenenbaum et al. (2000). Under such mapping, the interstate distances of the original attractor $\|\mathbf{y}^{IN}_{n+1}-\mathbf{y}^{IN}_{n}\|$, and of the projected attractors $\|\varphi(\mathbf{y}^{IN}_{n+1})-\varphi(\mathbf{y}^{IN}_{n})\|$ in the rRNN’s space, are bound to the range $[(1-\epsilon_{1}),(1+\epsilon_{2})]$ according to:

[TABLE]

where $\{\varphi(\mathbf{y}^{IN}_{n}),\varphi(\mathbf{y}^{IN}_{n+1})\}$ are states built by rRNN node responses, and assigned to an embedding dimension using the CCA. Details in the estimation of $\{\epsilon_{1},\epsilon_{2}\}$ are added in Appendix A.

As $\tau_{0}^{net}=-12$ was found to be associated to a good prediction performance using approximately half of the initial set of nodes, we show in Fig. 5(a) the estimation of $\{\epsilon_{1}(\mu),\epsilon_{2}(\mu)\}$ for the rRNN for which nodes’ CCA-lag positions are within windows of width $\delta\tau_{0}^{net}=3$, centered at integer multiples of $\tau_{0}^{net}=-12$. Here, we present the average of the statistical distribution that includes $20$-network models and predicting 20 different MG time series at 300 time steps into the future. In panel (b), we schematically illustrate the relevant geometrical properties of the attractors mapped onto the rRNN’s space. Such study considered $\tau_{0}=-12$ with which we obtain the set of coordinates $\mathbf{y}_{n}=\{y_{n},y_{n-12},y_{n-24},y_{n-36}\}$ that we use to unfold the state space object; where $\mathbf{y}^{IN}_{n}=\mathbf{y}_{n}$.

The consequence of an increasing $\mu$ in Fig. 5(a) can be explained with the graphic representation of the limits shown by Fig. 5(b). Here, we illustrate the three general cases (I, II, III) connected to their corresponding ranges in $\mu$ in Fig. 5(a). The evolution of the system’s trajectory is illustrated along a curve sampled at the positions of the big black dots. Gray dots are network’s neighbor states to the state $y^{IN}_{n+1}$. The first case (I) corresponds to neighbors which form a dense cloud of samples since $\epsilon_{1}\sim 1$, but are arranged closely around the original sample since $\epsilon_{1}\lesssim 1$ and $\epsilon_{2}<1$. The neighbor samples insufficiently enhance diversity in feature extraction and the network cannot predict the system’s future evolution. This is confirmed by Fig. 5(c,d), which show bad prediction performance and unity divergence: the RNN cannot predict the system.

Case (II) includes the values of $\mu$ where $\epsilon_{1}\simeq 1$ and $\epsilon_{2}\geq 1$. Within this parameter range, the system’s prediction performance strongly increases until reaching the lowest prediction error. Our analysis reveals the following mechanism behind this improvement. According to $\epsilon_{2}$, the maximum interstate distance possible inside the rRNN’s space is twice the interstate distance of the original trajectory. As a consequence, the rRNN samples neighbors to state $\mathbf{y}^{IN}_{n+1}$. Hence, as the state neighborhood is broader, there now is a sufficient random scanning of the attractor’s vicinity, such that the network can use the different features to solve prediction. The network can therefore use the projected objects to predict, which is confirmed by excellent performance according to Fig. 5(c,d).

The last case (III) appears for $\mu>1.3$, where all approximated distances of the embedded attractor are much larger than the original distance. The rRNN’s autonomous dynamics therefore enlarge the sampling distance such that no dense nearest neighbors are anymore available for prediction. The result are the distortions caused by the autonomous network’s dynamics typically found to be chaotic for $\mu\gtrsim 1.4$, already identified in our previous work Marquez et al. (2018). Consequently, the network’s folding property distorts the projected features. Therefore, $\epsilon_{1}$ becomes undefined, meaning that information about the structure of the embedded trajectory is lost. As a consequence, prediction performance strongly reduces, see Fig. 5(c,d).

The process described in this section allowed us to identify and use relevant features that benefit good long-term prediction performances with a downsized rRNN. Additionally, we found that the neighborhood generated by the inter-state distances of such features has an meaningful impact on the network’s ability to predict at all. At this point, we show how to simplify even more our process and package all the above described steps that allow us to generate features related to good prediction performances.

IV Hybrid Takens-inspired random recurrent networks

In this section, we directly exploit our newly gained insight and introduce a modified version of the classical random neural network for time series prediction. Here, we design a system which aims to only take into account such Takens-like dimensions that we found to be relevant for prediction. As it was previously described, actions provided by the nodes of the rRNN can be interpreted in the light of delay-embedding. We consequently modify the classical rRNN by including a Takens-inspired external memory:

[TABLE]

where $\textbf{x}_{n+1+\tau_{T}}$ is a delayed term added to the output layer, see Fig. 6(a). All elements of the reservoir layer have been copied and then time shifted by a delay term $\tau_{T}$ that could be the Takens embedding delay $\tau_{0}$. This process allows us to add virtual nodes to our network, which are distributed in the delay lines. This Takens rRNN (TrRNN) combines nonvolatile external memory (virtual nodes) with a neural network (real nodes) and therefore shares functional features of the recently introduced hybrid computing concept Graves et al. (2016). Yet, our concept makes any additional costly optimization unnecessary.

We start our analysis by identifying the embedding delay related to the best prediction performance. Here, we fix $\mu=0.1$, and modify then the delay term $\tau_{T}\in[-20,-1]$. For $\mu=0.1$, the delay-coordinates found in the network’s space only span approximately two Takens-like embedding dimensions of the MG system with delays $\{2\tau_{0},0\}$, when $\tau_{0}=-12$, see Fig. 6(b). Furthermore, most node responses are distributed along the columns centered in lags $\{2\tau_{0},0\}$. Consequently, as shown in Fig. 5(c,d), the prediction performance is almost the lowest possible due to insufficient dimensionality to get attractor-like features. Additionally, we set the number of nodes to $350$, which is $70\%$ of the nodes that the best CCA-windows filtered rRNN had to disposal, see Sec. III.A; and it uses $35\%$ of the nodes used by the non-filtered classical rRNN.

According to Fig. 7(a,b), the best average prediction performance, for 20-network model over 20 MG different time series at 300 time steps into the future, is found for $\tau_{T}=-12$, belonging to an interval $\tau_{T}\in[-13,-10]$ with the lowest NMSE values and divergent rates. The delay $\tau_{T}$ therefore agrees with the one identified for Takens attractor embedding $\tau_{0}$ and is almost identical to one of the lags $\tau_{0}^{net}$ found optimal in the CCA-window filtering. Furthermore, here the system only has as many CCA-windows at its disposal as dimensions required to embed the MG attractor. This lifts the disambiguate present in the CCA-window filtering analysis, and consequently the optimum performance is found only for an TrRNN embedding exactly along lags according to Takens embedding. In comparison to the classical rRNN, our TrRNN achieves the same performance, simultaneously reducing the amount of nodes in the reservoir layer from $1000$ to $350$ in the output when compared to the CCA filtered system. Compared to the pristine rRNN, we obtain one order of magnitude better performance with a network three times smaller.

Figure 7(c) shows the estimation of $\{\epsilon_{1},\epsilon_{2}\}$ with the variation of $\tau_{T}$. As it can be seen, $\epsilon_{2}\approx 1$ is associated to good prediction performances, found for $\tau_{T}\in[-13,-10]$. This result agrees with the results provided by the classical random network in Sec. III.B. The CCA for $\tau_{T}=-12$ is shown by Fig. 7(d), where we can find the set of nodes with all delay-coordinates required to fully reconstruct the MG attractor. In the cases where prediction was not possible, the CCA identifies the non-adequacy of the rRNN delay embedding as the reason, see Fig. 7(e) for $\tau_{T}=-3$.

IV.1 Application: control an arrhythmic neuronal model

We directly utilize our TrRNN as a part of an efficient feedback control mechanism in an arrhythmic excitable system. We task the TrRNN to aid stabilizing a system which models the firing behavior of a noise-driven neuron. It consists in the FitzHugh-Nagumo (FHN) neuronal model Longtin (1993); Christini and Collins (1995),

[TABLE]

where $v(t)$ and $w(t)$ are voltage and recovery variables. $I=0.3$ is a tonic activation signal, $\xi$ is Gaussian white noise with zero mean and standard deviation $\sim 0.02$, $\epsilon=0.005$, $g=0.5$, $D=1.0$, and $H=0.15$. These equations have been solved by the Euler-Maruyama algorithm for stochastic differential equation’s integration. In its resting state, the neuron’s membrane potential is slightly negative. Once the membrane voltage $v(t)$ is sufficiently depolarized through an external stimuli, the neuron spikes due to the rise of the action potential Enoka (2015); Tirozzi et al. (2007). The time between consecutive spikes are defined as interspike intervals (ISIs). These random ISIs are shown by dots in black in Fig. 8(a), clearly showing that spiking is non-regular.

We aim to control this random spiking behavior of the FHN neuronal model by proportional perturbation feedback (PPF) method Garfinkel et al. (1992) and by using either a TrRNN or a rRNN. The PPF method consists in the application of perturbations to locate the system’s unstable fixed point onto a stable trajectory Garfinkel et al. (1992); Schiff et al. (1994). This method is used to fit instabilities in the FHN neuronal model through the design of a control equation. In our case, the goal of using the PPF method is to build a control subsystem, which applies an external stimuli to trigger spiking and reduce the degree of chaos.

In our approach, the past information provided by voltage $v(t)$ in the FHN model is used to determine two things: (i) the parameters to design the control equation, and (ii) training parameters for rRNN ($\mu=1.1$) and TrRNN ($\mu=0.1$, and $\tau_{T}=-166$ obtained via the ACF minimum of $v(t)$ as in the Takens’ scheme) to predict future values of $v(t)$. The predicted $v(t)$ is used to calculate the full control signal with which we stabilize the neuron’s spiking activity. This methodology allows us to replace the quantity under control $v(t)$ by a reconstructed signal, which in control theory is related with the replacement of sensors in an exemplary control system. To train the network, we inject $1\times 10^{5}$ values and we let the network run freely for other $4\times 10^{6}$ steps, allowing us to stabilize $5619$ ISI points. We then evaluate the quality of the stabilization for networks ranging from $11$ to $340$ nodes. In Fig. 8(a), the set of blue dots along a constant line shows how the TrRNN can stabilize the ISI activity. As it can be seen, the network can control the random ISI starting from $n=7465$. The excellent stabilization was achieved with TrRNN with only $12$ nodes.

Figure 8(b) shows the full comparison between rRNN and TrRNN. The mean value of ISI is calculated for the different sizes of rRNN and TrRNN and then normalized by the mean value of the random ISI. The TrRNN starts inferring the inner dynamics of the FHN system for an extremely small network containing just $12$ nodes, from which point on it is always capable to correctly stabilize the ISI. In contrast, the classical rRNN does not predict at all until its architecture has at least $80$ nodes, but performance remains poor in comparison with TrRNN. For $200$ nodes the rRNN starts predicting the dynamic of the FHN system more or less correctly, allowing the control signal to fully stabilize the ISI. Yet, for more than $200$ nodes the good performance still can fluctuate, even significantly dropping again. This indicates that in general the stabilization via a classical rRNN is not robust. Furthermore, with the TrRNN one can reduce the number of nodes to $15$ times less than the classical rRNN. This stark difference in performance highlights (i) the difficulty of the task, and (iI) the excellent efficiency that the addition of a simple, linear delay term adjusted to the Takens embedding delay brings to the system. Our TrRNN therefore is not only a interesting novel ANN-concept for the prediction of complex systems. It also once more confirms our theory: embedding inside ANNs is essential if the systems are to be used for prediction.

V Comparison with classical state space prediction

Up to now, we have been describing a method for downsizing the hidden layer of rRNNs which can be summarized by the same steps that should be followed if we attempt to solve continuous-time chaotic signal prediction in the state space framework. This framework can be divided according to three fundamental aspects Weigend and Gershenfeld (1993); Kantz and Schreiber (1997); Farmer and Sidorowich (1987): (i) insufficient information to represent the complete state space trajectory of the chaotic system. This problem originates from the fact that in many cases one does not have access to all state space dimensions. In this case a reconstruction of the dynamics along the chaotic system’s missing degrees of freedom is required. The knowledge of all dimensions allows us to design predictors based on the full state space trajectories. (ii) The second problem is related to the sampling resolution. All information that is acquired, be it from simulations or from experiments, comes with a particular resolution. To minimize the divergence between a prediction and the correct value, the sampling resolution has to be maximized. This is of particular importance for prediction of chaotic systems as these by definition show exponential divergence. (iii) For deterministic chaotic systems, future states of a given trajectory can in principle be approximated from the exact knowledge of the present state. Therefore, the final step towards prediction is approximating the underlying deterministic law ruling the dynamical system’s evolution.

Therefore, step (i) is fulfilled by the random mapping which takes place in the high dimensional space of the network. Here, RPT supports the fact that the original input data are randomly mapped onto the dimensions of the projection space, and then the structural damage to the original object is minimized. Step (ii) is fulfilled by the analysis made in Section III.B, where the sampling resolution is maximized to cover the region between the states $\mathbf{y}^{IN}_{n}$ and $\mathbf{y}^{IN}_{n+1}$. Finally, (iii) relates to the training itself of the rRNN, where $W^{out}$ has to be determined via regression.

VI Conclusion

We have introduced a new method of rRNNs analysis which demonstrates how prediction is potentially achieved in high dimensional nonlinear dynamical systems. Random recurrent networks and prediction of a specific signal can consequently be described via a common methodology. Quantifying measures such as the memory related cross-correlation analysis and the feature extraction are quantitatively interpretable. We therefore significantly extend the toolkit previously available for random neural network analysis. Tools developed in the article might be comparable to the utilization of the t-SNE van der Maaten and Hinton (2012) technique for analyzing ANNs during a classification task.

Our scheme has numerous practical implications. The most direct is motivating the development and analysis of new learning strategies. Furthermore, we already designed a novel hybrid-computer which includes both, virtual and real nodes, that efficiently predicts via a priori defined external memory access rules. This approach allows us to improve the design of our neural network in order to reduce the number of nodes and connections required to solve prediction. Finally, our work partially removes the black-box property of random recurrent networks for prediction, possibly giving translational insight into how such tasks can be solved in comparable systems.

Acknowledgements

Funding for B.A.M. and B.J.S. was provided by the 2019 Queen’s postdoctoral fellowship fund, the Natural Sciences and Engineering Research Council of Canada (NSERC), and the Queen’s Research Initiation Grant (RIG).

Appendix A Estimation of $\{\epsilon_{1},\epsilon_{2}\}$

Each state in Takens space is described by $M$ delay-coordinates,

[TABLE]

The second step is to define the corresponding two arbitrary consecutive states $\{\varphi(\mathbf{y}^{IN}_{n}),$ $\varphi(\mathbf{y}^{IN}_{n+1})\}\in\mathbb{R}^{h}$, where $h$ depends on $\mu$. The value of $h$ is determined from the CCA, where we approximately assign the mapped objects dimensionality to the number of elements found in the interval $[l_{min},l_{max}]$ for each $\mu$, see Fig. 2(b,c). In order to construct those projected states, we use all the delay-coordinates provided by the network, i.e. the full range $[l_{min},l_{max}]$ for each value of $\mu$, as follows

[TABLE]

where $\{\varphi_{l_{1}}(\mathbf{y}^{IN}_{n}),\varphi_{l_{2}}(\mathbf{y}^{IN}_{n}),\ldots\}$ are node responses lagged at $[l_{min},l_{max}]$. The size of the interval $[l_{min},l_{max}]$ depends on the value of $\mu$, as it was shown by Fig. 2(b,c), where we find a broader distribution of delay-coordinates for higher values of $\mu$.

The interstate distances $\|\mathbf{y}^{IN}_{n+1}-\mathbf{y}^{IN}_{n}\|$ and $\|\varphi(\mathbf{y}^{IN}_{n+1})-\varphi(\mathbf{y}^{IN}_{n})\|$, have to be bounded in the interval $[(1-\epsilon_{1}),(1+\epsilon_{2})]$ according to

[TABLE]

Under these conditions, we can claim that the transformation by the rRNN agrees with a nonlinear random projections. Estimating limits $\{\epsilon_{1},\epsilon_{2}\}$ requires to find the inferior $\epsilon_{min}$, and superior $\epsilon_{max}$ interstate distance limits:

[TABLE]

[TABLE]

where $\epsilon_{1}$ and $\epsilon_{2}$ are calculated by isolating these constants from $\epsilon_{min}=(1-\epsilon_{1})$, and $\epsilon_{max}=(1+\epsilon_{2})$. These limits contain information about the minimum and maximum distortions that we can find in order to get the best neighbors in the rRNN. $\|\varphi(\mathbf{y}^{IN}_{n+1})-\varphi(\mathbf{y}^{IN}_{n})\|_{min}$ and $\|\varphi(\mathbf{y}^{IN}_{n+1})-\varphi(\mathbf{y}^{IN}_{n})\|_{max}$ are calculated by using Euclidean distance under minimum and maximum norms,

[TABLE]

[TABLE]

where $\varphi_{l_{g}}(\mathbf{y}^{IN}_{n})$ are node responses lagged at $l_{g}\in[l_{min},l_{max}]$, $\forall g=1,2,\ldots,h$.

Here, we therefore identify the smallest and largest distances $[\varphi_{l_{g}}(\mathbf{y}^{IN}_{n+1})-\varphi_{l_{g}}(\mathbf{y}^{IN}_{n})]_{min,max}$ along each delay coordinate. Then, it is true that these smallest and largest distances bound the Euclidean distance of these. Finally, we determine $\|\mathbf{y}^{IN}_{n+1}-\mathbf{y}^{IN}_{n}\|$ via

[TABLE]