Integer Echo State Networks: Efficient Reservoir Computing for Digital Hardware

TL;DR

This paper introduces integer Echo State Networks (intESN), a hardware-efficient reservoir computing model using integer vectors and cyclic shifts, achieving high energy efficiency with minimal performance loss.

Contribution

The paper presents a novel intESN architecture that replaces traditional matrix operations with cyclic shifts, enabling efficient digital hardware implementation.

Findings

Significant reduction in memory footprint and computational complexity.

Confirmed energy efficiency improvements on FPGA hardware.

Maintained competitive performance on typical reservoir computing tasks.

Abstract

We propose an approximation of Echo State Networks (ESN) that can be efficiently implemented on digital hardware based on the mathematics of hyperdimensional computing. The reservoir of the proposed integer Echo State Network (intESN) is a vector containing only n-bits integers (where n<8 is normally sufficient for a satisfactory performance). The recurrent matrix multiplication is replaced with an efficient cyclic shift operation. The proposed intESN approach is verified with typical tasks in reservoir computing: memorizing of a sequence of inputs; classifying time-series; learning dynamic processes. Such architecture results in dramatic improvements in memory footprint and computational efficiency, with minimal performance loss. The experiments on a field-programmable gate array confirm that the proposed intESN approach is much more energy efficient than the conventional ESN.

Figures

Tables

Table 1. TABLE I: Details of datasets for time-series classification.

Univariate datasets from UCR
Name	#V	Train	Test	#C
Swedish Leaf	1	500	625	15
Distal Phalanx	1	139	400	3
ECG	1	100	100	2
Wafer	1	1000	6164	2
Multivariate datasets from UCI
Character Trajectories	3	300	2558	20
Spoken Arabic Digit	13	6600	2200	10
Japanese Vowels	12	270	370	9

Table 2. TABLE II: Area utilization of the hardware architectures.

$N$	LUT	FF	BRAM	DSP48
ESN
$100$	$8317$	$8724$	$31.5$	$47$
$200$	$8368$	$8770$	$87.5$	$47$
$300$	$15950$	$8920$	$135.5$	$47$
intESN
$100$	$4577$	$5488$	$7.5$	$5$
$200$	$4602$	$5492$	$11.5$	$6$
$300$	$4671$	$5497$	$12.5$	$6$

Table 3. TABLE III: Number of clock cycles (time) involved in the sequence recall task.

$N$	software ESN	hardware ESN	hardware intESN
$100$	$3.0\text{×}{10}^8$ (3.04 $s$)	$5.8\text{×}{10}^8$ (5.77 $s$)	$1.4\text{×}{10}^8$ (1.42 $s$)
$200$	$1.0\text{×}{10}^9$ (9.99 $s$)	$1.9\text{×}{10}^9$ (18.74 $s$)	$2.8\text{×}{10}^8$ (2.82 $s$)
$300$	$2.1\text{×}{10}^9$ (21.16 $s$)	$3.9\text{×}{10}^9$ (38.90 $s$)	$4.5\text{×}{10}^8$ (4.46 $s$)

Table 4. TABLE IV: Speed-up comparison with the use of pipelining.

	ESN		intESN
Pipelining	✓	✗	✓	✗
Clock cycles	$1,297,072$	$51,948$	$148,241$	$2,856$
Speed-up	$24.97$		$51.91$

Table 5. TABLE V: XPE overall power consumption (Watts)

$N$	software ESN	hardware ESN	hardware intESN
$100$	$1.53$	$1.73$	$1.59$
$200$	$1.53$	$1.78$	$1.60$
$300$	$1.53$	$1.95$	$1.60$

Equations

$$\text{x}(n)=\tanh (\rho \text{W}\text{x}(n-1)+\beta {\text{W}}^{\text{in}}\text{u}(n)+\beta {\text{W}}^{\text{back}}\text{y}(n-1)),$$

$$\hat{\text{y}}(n)=g({\text{W}}^{\text{out}}[\text{x}(n);\text{u}(n)]),$$

$$\text{dist}={\text{x}}^{\top }\text{y}$$

$$f_{\kappa }(x)=\{\begin{array}{ll}-\kappa & x\le -\kappa \\ x& -\kappa <x<\kappa \\ \kappa & x\ge \kappa \end{array}$$

$$\text{x}(n)=f_{\kappa }(\text{Sh}(\text{x}(n-1),1)+{\text{u}}^{\text{HD}}(n)+{\text{y}}^{\text{HD}}(n-1)).$$

$$\text{u}(n{)}_q=\lfloor 200\text{u}(n)\rceil /200,$$

$$y(n)=0.5\sin (n/4).$$

$$y(n{)}_q=\lfloor 100y(n)\rceil /100.$$

Full text

I Introduction

Recent work in Reservoir Computing [1, 2] (RC) illustrates how a Recurrent Neural Network with fixed connectivity can memorize and generate complex spatio-temporal sequences. RC has been shown to be a powerful tool for modeling and predicting dynamic systems, both living [3] and technical [4, 5]. Recently, it has been shown in [6] that RC is able to predict large chaotic systems.

Recent work on feed-forward networks shows that the binarization of filters in Convolutional Neural Networks can lead to enormous gains in memory and computational efficiency [7]. Reducing the memory allocated to each neuron or synapse from a 32-bit float to a few bits or binary saves computation with minimal loss in performance (see, e.g., [8, 9]). The increase in efficiency broadens the range of applications for these networks.

This article addresses two important research directions in RC: training reservoir networks and implementing networks efficiently. We discovered several direct functional similarities between the operations in RC and those of hyperdimensional computing (HDC) [10]. HDC [11] or more generally Vector Symbolic Architectures [12] are frameworks for neural symbolic representation, computation, and analogical reasoning. The distinction from the traditional computing is that all entities (objects, phonemes, symbols) are represented by random vectors of very high dimensionality – several thousand dimensions. Complex data structures and analogical reasoning are implemented by simple arithmetical operations (binding, addition/bundling, and permutation) and a well-defined similarity metric [11]. Specifically, RC and HDC are connected by the following core principles:

• •

Random projections of input values onto a reservoir (which in essence is a high-dimensional vector) matches random HDC representations stored in a superposition;

• •

The update of the reservoir by a random recurrent connection matrix is similar to HDC binding/permutation operation;

• •

The nonlinearity of the reservoir can be approximated with the thresholded addition of integers in HDC.

We exploited these findings to design integer Echo State Networks (intESN), which perform like Echo State Networks (ESN) but with smaller memory footprint and computational cost.

In the proposed architecture, the reservoir of the network contains only constrained integers for each neuron, reducing the memory of each neuron from a 32-bit float to only a few bits. The recurrent matrix multiply update is replaced by a permutation (or even a cyclic shift), which results in the dramatic boosting of the computational efficiency. We validate the architecture on several tasks common in the RC literature. All examples demonstrate satisfactory approximation of performance of the conventional ESN; while the implementation on field-programmable gate array confirms the amenability of intESN for digital hardware.

The article is structured as follows. Background and related work are presented in Section II. The main contribution – integer Echo State Networks are described in Section III. The performance evaluation follows in Section IV. Section V presents the experiments on digital hardware. Sections VI and VII present a discussion and conclusions.

II Background and related work

There are many practical tasks that require the history of inputs to be solved. In the area of artificial neural networks (ANN), such tasks require working memory. This could be implemented by recurrent connections between neurons of an RNN. Training RNNs is much harder than that of feed-forward ANNs (FFNNs) due to vanishing gradient problem [13].

The challenge of training RNNs was addressed from two approaches. One approach eliminates the vanishing gradient problem through neurons with special memory gates [14]. Another approach is to reformulate the training process by learning only connections to the last readout layer while keeping the other connections fixed. This approach originally appeared in two similar architectures: Liquid State Machines [15] and ESNs [16], now referred to as reservoir computing[1].

It is interesting to note, that similar ideas were conceived in the area of FFNNs, which can be seen as an RNN without memory, and are known under the name of Extreme Learning Machines (ELMs) [17]. ELMs are used to solve various machine learning problems including classification, clustering, and regression [18].

Important applications of RC are the modeling and predicting of complex dynamic systems. Generating and predicting chaotic systems was an important use-case from the beginning [19], for example, ESNs were used for chaotic time-series from low-order aberrations caused by turbulence [20]. A thorough study on emulating chaotic systems was recently presented in [21]. It was also shown that ESNs can be used for forecasting of Electroencephalography signals and for solving classification problems in the context of Brain-Computer Interfaces [22]. There are different classification strategies and readout methods when performing classification of time-series with ESN. In [23] three classification strategies and three readout methods were explored under the conditions that testing data is purposefully polluted with noise. Interestingly, different readout methods are preferable in different noise conditions. Recent work in [24] also studied classification of multivariate time-series with ESN using standard benchmarking datasets. The work covered several advanced approaches, which extend the conventional ESN architecture, for generating a representation of a time-series in a reservoir.

Another recent research area is binary RC with cellular automata (CARC) which started as a interdisciplinary research within three areas: cellular automata, RC, and HDC. CARC was initially explored in [25] for projecting binarized features into high-dimensional space. Further in [26], it was applied for modality classification of medical images. The usage of CARC for symbolic reasoning is explored in [27]. The memory characteristics of a reservoir formed by CARC are presented in [28]. Work [29] proposed the usage of coupled cellular automata in CARC. Examples of recent RC developments also include advanced architectures such as Laplacian ESN [30], learning of reservoir’s size and topology [31], new tools for investigating reservoir dynamics [32] and determining its edge of criticality [33].

The design of ESNs has been an important research area (see, e.g., [34, 35, 36]). One important aspect of the design is, of course, the choice of network’s parameters for a given task. Another important aspect considered in this study is the computational complexity. One of the ways of reducing computational costs would be to use quantized reservoir states. It was explored in Fractal Prediction Machines [37] and Neural Prediction Machines [38, 39] RC models. Another way of reducing computational costs involves modifications of network’s connectivity structure. In this respect, the approach which is ideologically closest to our intESN, was presented in [40]. The authors demonstrated that a simple cycle reservoir (referred to as the ring-based ESN) can be used to achieve a performance similar to the conventional ESN. Similar conclusions about the ring-based ESN were obtained in [36] when studying different design strategies for reservoir connection matrices in four typical RC tasks. While the ring-based ESN explored reservoir update solution, which is similar to one of our optimizations, the technical side is very different from our approach as intESN strives at using only integers as neurons activation values.

While in this article the main focus is on reservoir states comprised of integers only, it is worth mentioning related works considering the general problem of reducing the computational complexity of the conventional ESN. For example, several optimizations were used in [41] in order to deploy an ESN on a resource-constrained device for anomaly detection. These optimizations include sparse matrix algebra via compressed row storage for weights of connections between between the input layer neurons and the reservoir; single floating point precision; and an activation function, which resembles $\tanh()$ function but has lower complexity. Similarly, in works [42, 43] ESNs were used in the context of user movements prediction on resource-constrained devices, therefore, the authors studied how the parameters of the network will affect computational costs and task performance. In particular, the varied parameters were sparsity of reservoir connection matrix, number of bits per weight, number of neurons in reservoir.

II-A Echo State Networks

This subsection summarizes the functionality of the conventional ESN, it follows the description in [44] for a special case of leaky integration when $\alpha=1$111For the detailed tutorial on ESNs diligent readers are referred to [44].. Fig. 1 depicts the architectural design of the conventional ESN, which includes three layers of neurons. The input layer with $K$ neurons represents the current value of input signal denoted as $\textbf{u}(n)$. The output layer ($L$ neurons) produces the output of the network (denoted as $\textbf{y}(n)$) during the operating phase. The reservoir is the hidden layer of the network with $N$ neurons, with the state of the reservoir at time $n$ denoted as $\textbf{x}(n)$.

In general, the connectivity of ESN is described by four matrices. $\textbf{W}^{\text{in}}$ describes connections between the input layer neurons and the reservoir, and $\textbf{W}^{\text{back}}$ does the same for the output layer. Both matrices project the current input and output to the reservoir. The memory in ESN is due to the recurrent connections between neurons in the reservoir, which are described in the reservoir matrix W. Finally, the matrix of readout connections $\textbf{W}^{\text{out}}$ transforms the current activity levels in the input layer and reservoir ($\textbf{u}(n)$ and $\textbf{x}(n)$, respectively) into the network’s output $\textbf{y}(n)$.

Note that three matrices ($\textbf{W}^{\text{in}}$, $\textbf{W}^{\text{back}}$, and W) are randomly generated at the network initialization and stay fixed during the network’s lifetime. Thus, the training process is focused on learning the readout matrix $\textbf{W}^{\text{out}}$. There are no strict restrictions for the generation of projection matrices $\textbf{W}^{\text{in}}$ and $\textbf{W}^{\text{back}}$. They are usually randomly drawn from either normal or uniform distributions and scaled as shown below. The reservoir connection matrix, however, is restricted to posses the echo state property. This property is achieved when the spectral radius of the matrix W is less or equal than one. For example, W can be generated from a normal distribution and then normalized by its maximal eigenvalue. Unless otherwise stated, in this article an orthogonal matrix was used as the reservoir connection matrix; such a matrix was formed by applying QR decomposition to a random matrix generated from the standard normal distribution. Also, W can be scaled by a feedback strength parameter, see (1).

The update of the network’s reservoir at time $n$ is described by the following equation:

[TABLE]

where $\beta$ and $\rho$ denote projection gain and the feedback strength, respectively. Note that it is assumed that the spectral radius of the reservoir connection matrix W is one. Note also that at each time step neurons in the reservoir apply $\tanh()$ as the activation function. The nonlinearity prevents the network from exploding by restricting the range of possible values from -1 to 1. The activity in the output layer is calculated as:

[TABLE]

where the semicolon denotes concatenation of two vectors and $g()$ the activation function of the output neurons, for example, linear or Winner-take-all.

II-A1 Training process

This article only considers training with supervised-learning when the network is provided with the ground truth desired output at each update step. The reservoir states $\textbf{x}(n)$ are collected together with the ground truth $\textbf{y}(n)$ for each training step. The weights of the output layer connections are acquired by solving the regression problem which minimizes the mean square error between predictions (2) and the ground truth. While this article does not focus on the readout training task, it should be noted that there are many alternatives reported in the literature including the usage of regression with regularization, online update rules, etc. [44].

II-B Fundamentals of hyperdimensional computing

In a localist representation, which is used in all modern digital computers, a group of bits is needed in its entirety to interpret a representation. In HDC, all entities (objects, phonemes, symbols, items) are represented by vectors of very high dimensionality – thousands of bits. The information is spread out in a distributed representation, which contrary to the localist representations, any subset of the bits can be interpreted. Computing with distributed representations utilizes statistical properties of vector spaces with very high dimensionality, which allow for approximate, noise-tolerant, highly parallel computations. Item memory (also referred to as clean-up memory) is needed to recover composite representations assigned to complex concepts. There are several flavors of HDC with distributed representations, differentiated by the random distribution of vector elements, which can be real numbers [45, 46, 47, 48], complex numbers [49], binary numbers [11, 50], or bipolar [46, 51].

We rely on the mathematics of HDC with bipolar distributed representations to develop intESN. Kanerva [11] proposed the use of distributed representations comprising $N=10,000$ binary elements (referred to as HD vectors). The values of each element of an HD vector are independent equally probable, hence they are also called dense distributed representations. Similarity between two binary HD vectors is characterized by Hamming distance, which (for two vectors) measures the number of elements in which they differ. In very high dimensions Hamming distances (normalized by the dimensionality $N$) between any arbitrary chosen HD vector and all other vectors in the HD space are concentrated around 0.5. Interested readers are referred to [11] and [52] for comprehensive analysis of probabilistic properties of the high-dimensional representational space.

The binary HD vectors can be equivalently mapped to the case of bipolar representations, i.e., where each vector’s element is encoded as “-1” or “+1”. This definition is sometimes more convenient for purely computational reasons. The distance metric for the bipolar case is a dot product:

[TABLE]

Basic symbols in HDC are referred to as atomic HD vectors. They are generated randomly and independently, and due to high dimensionality will be nearly orthogonal with very high probability, i.e., similarity (dot product) between such HD vectors is approximately 0. An ordered sequence of symbols can be encoded into a composite HD vector using the atomic HD vectors, the permutation (e.g., cyclic shift as a special case of permutation) and bundling operations. This vector encodes the entire sequence history in the composite HD vector and resembles a neural reservoir.

Normally in HDC, the recovery of component atomic HD vectors from a composite HD vector is performed by finding the most similar vectors stored in the item memory222It is not common to do such decoding in RC. Normally, in the scope of RC a readout matrix is learned. In this article, we follow this standard RC approach to extracting information back from a reservoir.. However, as more vectors are bundled together there is more interference noise and the likelihood of recovering the correct atomic HD vector declines.

Our recent work [10] reveals the impact of interference noise, and shows that different flavors of HDC have universal memory capacity. Thus, the different flavors of HDC can be interchanged without affecting performance. From these insights, we are able to design much more efficient networks for reservoir computing for digital hardware.

III Integer Echo State Networks

This section presents the main contribution of the article – an architecture for integer Echo State Network. The architecture is illustrated in Fig. 2. The proposed intESN is structurally identical to the the conventional ESN (see Fig. 1) with three layers of neurons: input ($\textbf{u}(n)$, $K$ neurons), output ($\textbf{y}(n)$, $L$ neurons), and reservoir ($\textbf{x}(n)$, $N$ neurons). It is important to note from the beginning that training the readout matrix $\textbf{W}^{\text{out}}$ for intESN is the same as for the conventional ESN (Section II-A1).

However, other components of intESN differs from the conventional ESN. First, activations of input and output layers are projected into the reservoir in the form of bipolar HD vectors [47] of size $N$ (denoted as $\textbf{u}^{\text{HD}}(n)$ and $\textbf{y}^{\text{HD}}(n)$). For problems where input and output data are described by finite alphabets and each symbol can be treated independently, the mapping to $N$-dimensional space is achieved by simply assigning a random bipolar HD vector to each symbol in the alphabet and storing them in the item memory [11, 53]. In the case with continuous data (e.g., real numbers), we quantized the continuous values into a finite alphabet. The quantization scheme (denoted as $Q$) and the granularity of the quantization are problem dependent. Additionally, when there is a need to preserve similarity between quantization levels, distance preserving mapping schemes are applied (see, e.g., [54, 55]), which can preserve, for example, linear or nonlinear similarity between levels. An example of a discretization and quantization of a continuous signal as well as its HD vectors in the item memory is illustrated in Fig. 2. Continuous values can be also represented in HD vectors by varying their density. For a recent overview of several mapping approaches readers are referred to [56]. Also, an example of applying such mapping is presented in Section IV-A2. Another feature of intESN is the way the recurrence in the reservoir is implemented. Rather than a matrix multiply, recurrence is implemented via the permutation of the reservoir vector. Note that permutation of a vector can be described in matrix form, which can play the role of W in intESN. Note that the spectral radius of this matrix equals one. However, an efficient implementation of permutation can be achieved for a special case – cyclic shift (denoted as $\text{Sh}()$). It is important to note that we have shown in [10] that the recurrent weight matrix W creates key-value pairs of the input data. Note that W is chosen randomly and kept fixed, and this always leads to the same properties. Moreover, there is no advantage of the fully connected random recurrent weight matrix over the simple cyclic shift operation for storing the input history. Thus, the use of the cyclic shift in place of a random recurrent weight matrix does not limit intESN’s ability to produce linearly separable representations. Fig. 2 shows the recurrent connections of neurons in a reservoir with recurrence by cyclic shift of one position. In this case, vector-matrix multiplication $\textbf{W}\textbf{x}(n)$ is equivalent to $\text{Sh}(\textbf{x}(n),1)$.

Finally, to keep the integer values of neurons, intESN uses different nonlinear activation function for the reservoir – clipping (4). Note that the simplest bundling operation is an elementwise addition. However, when using the elementwise addition, the activity of a reservoir (i.e., a composite HD vector) is no longer bipolar. From the implementation point of view, it is practical to keep the values of the elements of the HD vector in the limited range using a threshold value (denoted as $\kappa$).

[TABLE]

The clipping threshold $\kappa$ is regulating nonlinear behavior of the reservoir and limiting the range of activation values. Note that in intESN the reservoir is updated only with integer bipolar vectors, and after clipping the values of neurons are still integers in the range between $-\kappa$ and $\kappa$. Thus, each neuron can be represented using only $\log_{2}(2\kappa+1)$ bits of memory. For example, when $\kappa=7$, there are fifteen unique values of a neuron, which can be stored with just four bits. We have also shown recently that the usage of the clipping might be beneficial when implementing resource-efficient alternatives of Self-Organizing Maps [57].

Summarizing the aforementioned differences, the update of intESN is described as:

[TABLE]

IV Performance evaluation

In this section, the proposed intESN approach is verified and compared to the conventional ESN and the ring-based ESN [40] on a set of typical RC tasks. In particular, three aspects are evaluated: short-term memory, classification of time-series, and modeling of dynamic processes. Short-term memories are compared using the trajectory association task [49], introduced in the area of holographic reduced representations [45]. Additionally, an approach for storing and decoding analog values using intESN is demonstrated on image patches. Classification of time-series is studied using the standard datasets from UCI and UCR. Modeling of dynamic processes is tested on two typical cases. First, the task of learning a simple sinusoidal function is considered. Next, networks are trained to reproduce a complex dynamical system produced by a Mackey-Glass series. Unless otherwise stated, ridge regression (the regularization coefficient is denoted as $\lambda$) with the Moore-Penrose pseudo-inverse was used to learn the readout matrix $\textbf{W}^{\text{out}}$. Values of input neurons $\textbf{u}(n)$ were not used for training the readout in any of the experiments below.

IV-A Short-term memory

IV-A1 Sequence recall task

The sequence recall task includes two stages: memorization and recall. At the memorization stage, a network continuously stores a sequence of tokens (e.g., letters, phonemes, etc). The number of unique tokens is denoted as $D$ ($D=27$ in the experiments), and one token is presented as input each timestep. At the recall stage, the network uses the content of its reservoir to retrieve the token stored $d$ steps ago, where $d$ denotes delay. In the experiments, the range of delay varied between 0 and 15.

For the conventional and ring-based ESNs, the dictionary of tokens was represented by a one-hot encoding, i.e. the number of input layer neurons was set to the size of the dictionary $K=D=27$. The same encoding scheme was adopted for the output layer, $L=27$. The input vector was projected to the reservoir by the projection matrix $\textbf{W}^{in}$ where each entry was independently generated from the uniform distribution in the range $[-1,1]$, the projection gain was set to $\beta=0.1$. The reservoir connection matrix W for the conventional ESN was first generated from the standard normal distribution and then orthogonalized. The reservoir connection matrix W for the ring-based ESN was generated as a permutation matrix. The feedback strength of both reservoir connection matrices was set to $\rho=0.94$.

For intESN, the item memory was populated with $D$ random high-dimensional bipolar vectors. The threshold for the clipping function was set to $\kappa=3$. The output layer was the same as in ESN with $L=27$ and one-hot encoding of tokens. It is worth noting that $\rho$ and $\kappa$ were chosen in such a way that the accuracy curves would resemble each other as close as possible. The diligent readers are kindly referred to the Supplementary materials (Fig. S.1) where for the case $N=200$ the curves for the range of $\rho$ and $\kappa$ values are presented.

For each value of of the delay $d$ a readout matrix $\textbf{W}^{\text{out}}$ was trained, producing 16 matrices in total. The training sequence presented 2000 random tokens to the network, and only the last 1500 steps were used to compute the readout matrices. The regularization parameter for ridge regression was set to $\lambda=0$. The training sequence of tokens delayed by the particular $d$ was used as the ground truth for the the activations of the output layer. During the operating phase, both the inclusion of a new token into the reservoir and the recall of the delayed token from the reservoir were simultaneous. Experiments were performed for three different sizes of the reservoir: $N=100$, $N=200$, and $N=300$.

The memory capacity of the network is characterized by the accuracy of the correct decoding of tokens for different values of the delay. Fig. 3 depicts the accuracy for all networks conventional ESN (solid lines), ring-based ESN (dash-dotted line) and intESN (dashed lines). The capacities of all the networks grow with the increased number of neurons in the reservoir. Since the capacities of the conventional ESN and the ring-based ESN are almost identical, which is in line with [36], for the rest of this subsection we assume both of these networks when using the term ESN.

The capacities of ESN and intESN are comparable for small $d$, i.e., for the most recent tokens. For the increased delays the curves featured slightly different behaviors. With increase of the value of $d$ the performance of intESN started to decline faster compared to ESN. Eventually, all curves converge to the value of the random guess which equals $1/D$. Moreover, the information capacity of a network is characterized by the amount of the decoded from the reservoir information. This amount is determined using the amount of information per token ($\log_{2}D$), the probability of correctly decoding a token at each delay value, and the concept of mutual information. We calculated the amount of information for all networks in Fig. 3 in the considered delay range. For 100 neurons intESN preserved 19.3% less information, for 200 and 300 neurons 21.7% less.

These results highlight a very important trade-off: the performance versus a complexity of implementation. While the performance of intESN is somewhat poorer in this task, one has to bear in mind its memory footprint. With the clipping threshold $\kappa=3$ only 3-bit are needed to represent the state of a neuron compared to 32-bit per neuron (the size of type float) in ESN. In other words, intESN allowed lowering the memory footprint of the reservoir by an order of magnitude by sacrificing only a fraction of the performance with respect to the information capacity. Thus, we conjecture that some reduction in the performance for ten folds memory footprint reduction is an acceptable price in applications on resource-constrained computing devices. On the other hand, we can check the performance of the networks with equal memory footprints. For this we increased the number of neurons in intESN so that the total memory consumed by the reservoir with the same clipping threshold $\kappa=3$ would match that of the conventional or ring-based ESN. This network is denoted as “intESN-large”. Since $\kappa=3$ requires only 3-bit, in order to get the memory footprint corresponding to ESN, intESN could use more than ten times more neurons. Thus, the memory footprint of intESN with $1000$ neurons corresponds to ESN with $100$ neurons; while intESN with $2000$ and $3000$ neurons correspond to ESN with $200$ and $300$ neurons respectively. The results for this case are presented in Fig. 4 (the training sequence was prolonged to 9000 random tokens). With such settings intESN-large has clearly higher information capacity. In particular, for ESN memory footprint with 100 neurons the decoded amount of information has increased 2.2 times while for 200 and 300 neurons it increased 1.6 and 1.3 times respectively. It is important to note, however, that while the memory consumed by the reservoir of intESN-large was comparable to the corresponding ESN, the readout matrix for intESN-large was larger and more computationally demanding than the ESN readout matrix since the size of a readout matrix is proportional to the number of neurons in the reservoir.

IV-A2 Storage of analog values in intESN

This subsection presents the feasibility of storing analog values in intESN using image patches as a showcase. It is important to emphasize that this subsection does not go into detailed comparisons with other methods as the main purpose here is the principal demonstration of the possibilities of storing continuous data in reservoirs consisting of integers in a limited range. In other words, with this showcase, we are aiming at demonstrating the feasibility of using integer approximation of neuron states in intESN to work with analog representations. A value of a pixel (in an RGB channel) can be treated as an analog value in the range between 0 and 1. For each pixel it is possible to generate a unique bipolar HD vector. The typical approach to encode an analog value is to multiply all elements of the HD vector by that value. The whole sequence is then represented using the bundling operation on all scaled HD vectors. The result of bundling can be used as an input to a reservoir. However, the resultant composite HD vector will not be in the integer range anymore. We address this problem by using sparsity. Instead of scaling elements of an HD vector, we propose to randomly set the fraction of elements of the HD vector to zeros, i.e., the HD vector will become ternary. The proportion of zero elements is determined by the pixel’s analog value. Pixels with values close to zero will have very sparse HD vectors while pixels with values close to one will have dense HD vectors, but all entries will always be [-1, 0, or +1]. The result of bundling of such HD vectors (i.e., HD vector for an image) will still have integer values. Such representational scheme allows keeping integer values in the reservoir but it still can effectively store analog values.

The examples of results are presented in Fig. 5. Top row depicts original images stored in the reservoir. The other rows depict images reconstructed from the reservoir. The following parameters of intESN were used (top to bottom): $N=64000$, $\kappa=11$; $N=32000$, $\kappa=8$; $N=16000$, $\kappa=6$; $N=8000$, $\kappa=4$. The values of $\kappa$ were optimized for a particular $N$. Columns correspond to the delay values (i.e., how many steps ago an image was stored in the reservoir) as in the previous experiment. As one would anticipate, the quality of the reconstructed images is improving for larger reservoir sizes. At the same time, the quality of the reconstructed images is deteriorating for larger delay values, i.e., the worst quality of the reconstructed image could be observed in the bottom right corner while the best reconstruction is located in the top left corner. Nevertheless, the main observation for this experiment is that it is possible to project analog values into the reservoir with integer values using the mapping via varying sparsity and then retrieve the values from the reservoir. Moreover, we have shown recently [58] that the mapping via varying sparsity could even be helpful when solving classification problems with a feed-forward variant of the ESN.

IV-B Classification of time-series

In this section, ESN (conventional and ring-based) and intESN networks are compared in terms of classification accuracy obtained on standard time-series datasets. Following [24] we used several (four) univariate datasets from UCR333UCR. Time Series Classification Archive [online], 2018. – Available online: https://www.cs.ucr.edu/%7Eeamonn/time_series_data_2018/. and several (three) multivariate datasets from UCI444UCI. Machine Learning Repository [online], 2019. – Available online: http://archive.ics.uci.edu/ml/datasets.html.. Details of datasets are presented in Table I. For each dataset, the table includes the name, number of variables (#V), number of classes (#C), and the number of examples in training and testing datasets.

Configurations of the networks were kept fixed for all datasets. In fact, the configuration of the conventional and ring-based ESNs were set in accordance to [24]: reservoir size was set to $N=800$, projection gain was set to $\beta=0.25$, the feedback strength was set to $\rho=0.99$. The regularization parameter for the ridge regression was set to $\lambda=1.0$. The intESN was also trained with the same $\lambda$. The clipping threshold for the intESN was set to $\kappa=7$. Also for intESN the quantized values of time-series were mapped to bipolar vectors using scatter codes [56, 59]. The input signal u(n) was quantized as:

[TABLE]

where $\lfloor*\rceil$ denotes rounding to the the closest integer. Two sizes of intESN’s reservoir were used. The first size corresponded to the size of the conventional and ring-based ESNs, i.e., $N=800$. The second size (“intESN-large”) corresponded to the same memory footprint555Except for the Japanese Vowels dataset where such reservoir size seemed to significantly overfit the training data. In that case, the number of neurons was increased twice. required for ESN reservoir assuming that one ESN neuron requires 32-bit while one intESN neuron requires 4-bit (when $\kappa=7$). Thus, “intESN-large” had $N=6400$ neurons.

The output layers of networks were representing one-hot encodings of classes in a dataset, i.e., for the particular dataset $L=$#C of that dataset. The readout layers of all networks were trained using time-series from a training dataset in the so-called endpoints mode [23] when only final temporal reservoir states for each time-series are used for training a single readout matrix.

The experimental accuracies obtained from the networks for the considered datasets are presented in Fig. 6 and  7. Fig. 6 presents the results for univariate datasets while Fig. 7 presents the results for multivariate datasets. The figures depict mean and standard deviation values across ten independent random initializations of the networks. Similar to subsection IV-A1, the accuracy of the conventional ESN and the ring-based ESN are almost identical, thus, for the rest of this subsection we assume both of these networks when using the term ESN.

The obtained results strongly depend on the characteristics of the data. However, it was generally observed that intESN with the memory footprint equivalent to ESN demonstrated higher classification accuracy. On the other hand, the classification accuracy of intESN with the same number of neurons as in ESN was similar to ESN’s performance for all considered datasets but two (“Swedish Leaf” and “Character Trajectories”) for which the accuracy degradation was sensible. We, therefore, conjecture that in a general case, one cannot guarantee the same classification accuracy as for ESN. The empirical evidence, however, shows that it is not infeasible. Since placing the reported results into the general context of time-series classifications is outside the scope of this article we do not further elaborate on fine-tuning of hyperparameters of intESN for the best classification performance. However, the interested readers are kindly referred to the Supplementary materials (Fig. S.2 and S.3) where several different values of $N$, $\kappa$, and $\rho$ were examined for each dataset.

IV-C Modeling of dynamic processes

IV-C1 Learning Sinusoidal Function

The task of learning a sinusoidal function [60] is an example of a learning simple dynamic system with the constant cyclic behavior. The ground truth signal was generated as follows:

[TABLE]

In this task, the input layer was not used, i.e. $K=0$ but the network projected the activations of the output layer back to the reservoir using $\textbf{W}^{\text{back}}$. The output layer had only one neuron ($L=1$). The reservoir size was fixed to $N=1000$ neurons. The length of the training sequence was 3000 (first 1000 steps were discarded from the calculation). For ESN, the feedback strength for the reservoir connection matrix was set to $\rho=0.8$, for both networks $\lambda$ was set to 0. A continuous value of the ground truth signal was fed-in to ESN during the training.

For intESN, in order to map the input signal to a bipolar vector the quantization was used. The signal was quantized as:

[TABLE]

The item memory for the projection of the output layer was populated with bipolar vectors preserving linear (in terms of dot product) similarity between quantization levels [55]. The threshold for the clipping function was set to $\kappa=3$.

In the operating phase, the network acted as the generator of the signal feeding its previous prediction (at time $n-1$) back to the reservoir. Fig. 8 demonstrates the behavior of intESN during the first 100 prediction steps. The ground truth is depicted by dashed line while the prediction of intESN is illustrated by the shaded area between 10% and 90% percentiles (100 simulations were performed). The figure does not show the performance of the conventional ESN as it just followed the ground truth without visible deviations. intESN clearly follows the values of the ground truth but the deviation from the ground truth is increasing with the number of prediction steps. It is unavoidable for the increasing prediction horizon but, in this scenario, it is additionally accelerated due to the presence of the quantization error at each prediction step. It is worth noting, however, that the quality of predictions for this task could be improved by increasing the value of $\kappa=3$, i.e., at the cost of extract memory allocated for each neuron. The diligent readers are kindly referred to the Supplementary materials (Fig. S.4) where several different values of $\kappa$ were examined. The next subsection will clearly demonstrate effects caused by the quantization process. The error is accumulated because every time when feeding the prediction back to the reservoir of intESN it should be quantized in order to fetch a vector from the item memory.

IV-C2 Mackey-Glass series prediction

A Mackey-Glass series is generated by the nonlinear time delay differential equation. It is commonly used to assess the predictive power of an RC approach. In this scenario, we followed the preprocessing of data and the parameters of ESN described in [4]. The parameters of intESN (including quantization scheme) were the same as in the subsection above. The interested readers are kindly referred to the Supplementary materials (Fig. S.5) where several different values of $N$ and $\kappa$ were examined. The length of the training sequence was 3000 (first 1000 steps were discarded from the calculation). Fig. 9 depicts results for the first 300 prediction steps. The results were calculated from 100 simulation runs. The figure includes four panels. Each panel depicts the ground truth, the mean value of predictions as well as areas marking percentiles between 10% and 90%. The lower right corner corresponds to intESN while three other panels show performance of ESN in three different cases related to the quantization of the data.

In these scenarios ESN was trained to learn the model from the quantized data in order to see to which extent it affects the network. The upper left corner corresponds to ESN without data quantization. In this case, the predictions precisely follow the ground truth. The upper right corner corresponds to ESN trained on the quantized data but with no quantization during the operational phase. In such settings, the network closely follows the ground truth for the first 150 steps but then it often explodes. The lower left corner corresponds to ESN where the data was quantized during both training and prediction. In this scenario, the network was able produce to produce good prediction just for the first few dozens of steps and then entered the chaotic mode where even the mean value does not reflect the ground truth. These cases demonstrate how the quantization error could affect the predictions especially when it is added at each prediction step. Note that intESN operated in the same mode as the third ESN. Despite this fact, its performance rather resembles that of the second ESN where the speed deviation of the ground truth is faster. At the same time, the deviation of intESN grows smoothly without a sudden explosion in contrast to ESN.

V Digital hardware experiments

In order to demonstrate the amenability of intESN for digital hardware, we used a field-programmable gate array (FPGA) and implemented three different architectures: software ESN as well as hardware accelerated ESN and intESN. All architectures were implemented on a ZedBoard FPGA666ZedBoard. Hardware User’s Guide [online], 2014. – Available online: http://zedboard.org/sites/default/files/documentations/ZedBoard˙HW˙UG˙v2˙2.pdf., which contains a dual core ARM Cortex A9 CPU interfacing with a programmable logic fabric. The Xilinx Vivado Design Suite777Xilinx. Vivado Design Suite User Guide [online], 2018. – Available online: https://www.xilinx.com/support/documentation/sw˙manuals/xilinx2018˙1/ug910-vivado-getting-started.pdf. and Vivado SDK888Xilinx. Generating Basic Software Platforms [online], 2018. – Available online: https://www.xilinx.com/support/documentation/sw˙manuals/xilinx2018˙3/ug1138-generating-basic-software-platforms.pdf. were used to design the hardware architectures and program the FPGA board, respectively. The efficiency evaluation (e.g., energy consumption) of the architectures was based on the recall stage of the sequence recall task as described in Section IV-A1.

V-A Architectures design

The software ESN was implemented using only CPU on the FPGA board. For the hardware acceleration experiments, we programed the CPU to generate the inputs and feed them into a hardware architecture (either for ESN or intESN), and then retrieve the outputs once the computation is over. The hardware ESN and intESN architectures were designed using Vivado High Level Synthesis (HLS) using the C programming language, and later on synthesized into hardware Intellectual Property components to be integrated in a larger hardware system comprised of: the Zynq Processing System, AXI Direct Memory Access, AXI Interconnect, ESN/intESN architecture, and various other peripherals used for clocking and resetting mechanisms. It is important to note that for resource limitations on the ZedBoard, no pipelining or hardware optimization directives have been used on the HLS designs in order to cope with the resource usage for growing reservoir sizes and remain within the board’s capacity. However, we still provide a comparison of a speed-up expected from the pipelining for the two hardware architectures.

V-B Evaluation methodology

For each architecture, the reservoir size was varied between $100$, $200$, and $300$ neurons. In the following, we report the number of clock cycles required to accomplish the sequence recall task, the power consumption, and the area utilization (i.e., resources) required for the hardware designs. The number of clock cycles was measured using a hardware timer incorporated into the designs. The area utilization is reported from the hardware synthesis reports provided by Vivado. All three architectures use the same clock frequency of $100$ MHz. The ZedBoard contains a $10$ m$\Omega$ shunt-resistor in series with the input supply of the whole board, which could be used to obtain the overall power consumption by measuring the voltage across it. However, at the desired scale of changes voltage fluctuation and measurement sensitivity made it almost impossible to properly perform a precise comparison. Therefore, instead, we used the Xilinx Power Estimator (XPE) tool999Xilinx. Power Estimator User Guide [online], 2018. – Available online: https://www.xilinx.com/support/documentation/sw˙manuals/xilinx2018˙3/ug440-xilinx-power-estimator.pdf. provided by the Vivado Design Suite to estimate the power consumption of each architecture, which is a standard option [61, 62, 63].

V-C Efficiency evaluation results

Table II presents the area utilization of the hardware architectures for different sizes of reservoir. It is clear that the resource utilization of the hardware ESN is always larger than that of hardware intESN. This is an empirical manifestation of the facts that intESN: a) requires a lower memory footprint ($\kappa=3$) and b) its machinery uses simpler operations, e.g., the clipping instead of $\tanh()$ and the cyclic shift instead of vector-matrix multiplication. It is important to note, that the drastic increase of LUT utilization when the reservoir size of ESN was set to $300$ is due to resource constraints on the FPGA board since the total number of usable BRAM units is 140. Thus, the increased number of LUTs was used as an alternative way of increasing the memory capacity of the board.

Table II presents the number of clock cycles (time in seconds) necessary to perform the operating phase of the sequence recall task for each architecture. The number of clock cycles was measured with the help of the hardware timer. The time was calculated using the known frequency of the clock ($100$ MHz) As expected, for each architecture the operation time increases with the increased reservoir size. However, for any reservoir size, intESN is several times faster than both implementations of ESN (at least $2.1$ against the software ESN when $N=100$). The gain is increasing with the increased reservoir size so that when $N=300$ the hardware intESN is $8.7$ times faster compared to the hardware ESN.

An interesting remark in Table III is that the hardware ESN seems slower than the software ESN. It is counterintuitive since the hardware architecture is expected to accelerate the computations. In the consider case, this fact is explained by the absence of pipelining (cf. Section V-A) what prevents further hardware optimizations. The pipelining would certainly increase the speed at the price of the increased area utilization, which would make the hardware designs being inconceivable on the target board. In order to evaluate the effect of the pipelining on the speed of the hardware architectures, we have used the Vivado HLS synthesis report to estimate the number of clock cycles per iteration when $N=300$. The results are presented in Table IV. It shows that the pipelining significantly decreases the number of clock cycles, which makes the hardware ESN much faster than the software one. At the same time, the speed-up obtained for the hardware intESN is still higher than that of the hardware ESN, which makes intESN even more efficient than ESN.

Table V presents the power consumption of each architecture. The first observation is that the consumption of the software ESN remains constant for different $N$. This is because the hardware itself remains fixed, and the software computations that are performed on the Cortex CPU are the same for each configuration. Moreover, the software ESN has the lowest power consumption because it is only comprised of the CPU and the hardware timer, whereas the other two include other hardware components in addition to those. However, for the hardware implementation of ESN one could see a noticeable increase for a larger reservoir sizes. There is also an increase for the hardware intESN but it is slower compared to the hardware ESN. Moreover, for the fixed reservoir size, intESN’s power consumption remains lower than that of ESN.

Fig. 10 depicts the overall energy consumption of each architecture for different sizes of reservoir. The energy consumption was calculated as a product of the overall power consumption reported in Table V and the operating time reported in Table III. Fig. 10 clearly demonstrates that the energy efficiency of the hardware intESN is much better than that of both the hardware and software ESNs. Scaling the reservoir size leads to an increase of the energy consumption of all architectures. However, the slope of the curve for intESN is lower than for ESN’s architectures. Therefore, the energy saved by the use of intESN drastically increases with increasing size of reservoir (e.g., for $N=300$ it needs $10.6$ times lower energy than for the hardware ESN), which is a strong argument in favor of intESN.

VI Discussion

VI-1 Efficiency of intESN

In principle, evaluation of the computational efficiency of the proposed intESN could be done even in a simplistic manner using only a computer. For example, we have performed the initial assessment by simplified execution time measurements with Matlab. We used our Matlab implementations of both networks for the trajectory association task (Section IV-A1) with $N=300$ neurons to compare the times of projecting data and executing a reservoir. ESN was implemented using 32-bit float type (type single in Matlab) while intESN was implemented using 8-bits integer type (type int8 in Matlab). On average, the time required by intESN was $3.9$ times less than that of ESN. Thus, even the training time needed for convergence with intESN is shorter than that of ESN because intESN is much faster when it comes to projecting data and executing a reservoir while the time for estimating the readout matrix would be the same for both networks.

Section V presented the proper evaluation of the computational efficiency of the proposed intESN approach and the conventional ESN using the FPGA board. The hardware implementation of intESN was compared to two reference implementations: software ESN and hardware ESN. The detailed benchmarking tests have supported our claims about the efficiency of intESN. It is also worth noting that the efficiency gains of the intESN would be less relative to the ring-based ESN as essentially both networks use the same efficient mechanism for reservoir’s connectivity.

VI-2 Hyperparameters

In comparison to ESN, intESN has fewer hyperparameters. The common hyperparameter is the reservoir size $N$. Two specific intESN hyperparameters are clipping threshold $\kappa$ and the mapping to the reservoir; while for ESN we have to choose $\rho$, $\beta$, and $\alpha$. When $N$ is fixed then $\kappa$ in intESN has an effect on the network’s memory similar to $\rho$ and $\beta$ in ESN (see the next subsection for details). However, the difference is that $\kappa$ takes only positive integer values. This has its pros and cons. It could be much easier to optimize a single hyperparameter in the integer range. On the other hand, having real-valued hyperparameters one can get a configuration providing much finer tuning of network’s memory for a given task.

With respect to the mapping (projection) to the reservoir, it is probably the most non-trivial part on intESN. Especially, when data to be projected are real numbers. In this article, we mention three different strategy for mapping real numbers to bipolar or ternary vectors: puncturing of non-zero elements (Section IV-A2), mapping preserving linear similarity between vectors [55], and nonlinear approximate mapping using scatter codes [56, 59]. We suggest that it is a useful heuristic to try each of the approaches and choose the one performing best.

VI-3 On equivalence of ESN and intESN in terms of forgetting time constants

Section IV-A1 presented the experimental comparison of storage capabilities of intESN and ESN. An analytical approach to the treatment of memory capacity of reservoir was presented in [10] (Section 2.4.4). The work introduced an analytical tool called the forgetting time constant (denoted as $\tau$). The forgetting time constant is a scalar characterizing the memory decay in a network. In the case of intESN, it can be calculated analytically using the clipping threshold $\kappa$. For ESN, currently only special cases can be analyzed. For example, when reservoir neurons are linear the feedback strength $\rho$ will determine $\tau$. The other example is when reservoir neurons use $\tanh$, the feedback strength is fixed to one, and the reservoir update rule is modified to $f(\textbf{x})=\gamma\tanh(\textbf{x}/\gamma)$, where $\gamma$ denotes a gain parameter. This parameter affects $\tau$, which could be calculated numerically. Fig. 11 presents the implicit comparison of this special case of ESN and intESN via their forgetting time constants which are determined by $\gamma$ and $\kappa$ respectively. Both parameters similarly (not far from linear in logarithmic coordinates) affect $\tau$. It allows arguing that the networks are close to being functionally identical in terms of storage capabilities. The development of an approach for estimating the forgetting time constant for ESN in a general case is a part of our agenda for the future work.

VI-4 Training the readout in a generator mode

In the experiments generating time-series, we used the teacher forcing approach for training a readout matrix. This, however, does not have to be the compulsory choice for intESN. We do not foresee any complications for applying other approaches allowing to modify the network’s behavior for producing complex target functions. In particular, FORCE method [2], which uses error-based modification of readout weights during the training process, can be used as it has a mode in which only weights of a readout matrix are changed while leaving the rest of the network fixed.

VII Conclusions

In this article we proposed an architecture for integer approximation of the reservoir computing, which is based on the mathematics of hyperdimensional computing. The neurons in the reservoir are described by integers in the limited range and the update operations include only addition, permutation (cyclic shift), and clipping. Therefore, the integer Echo State Network has substantially smaller memory footprint and higher computational efficiency compared to the conventional Echo State Network with the same number of neurons in the reservoir. The actual number of bits for representing a neuron depends on the clipping threshold $\kappa$, but can be significantly lower than 32-bit floats in Echo State Network. For example, in our experiments the results were obtained with $\kappa=3$ and $\kappa=7$, which effectively makes it sufficient to represent a neuron with only three or four bits respectively. We demonstrated that the performance of the integer Echo State Network is comparable to the conventional Echo State Network in terms of memory capacity, potential capabilities for classification of time-series and modeling dynamic systems. The better performance was observed when the memory footprint of reservoir of the integer Echo State Network was set to that of the conventional Echo State Network. The experiment on the digital hardware have validated the amenability of the integer Echo State Network for significantly improving the energy efficiency of computations. Further improvements can be made by optimization of the parameters and better quantization schemes for handling continuous values. Naturally, due to the peculiarity of input data projection into the integer Echo State Network, the performance of the network in tasks for modeling dynamic systems is to a certain degree lower than that of the conventional Echo State Network. This, however, does not undermine the importance of integer Echo State Networks, which are extremely attractive for memory and power savings, and in the general area of approximate computing, where errors and approximations are becoming acceptable as long as the outcomes have a well-defined statistical behavior.