An Echo State Network (ESN) is a type of recurrent neural network, where the internal layer has randomly assigned connections and weights among the nodes. ESNs, first introduced simultaneously in  and in  under the name Liquid State Machines (LSM), have shown promise in performing classification and pattern prediction on spatio-temporal datasets. ESNs and LSMs, collectively studied under the name Reservoir Computing , have been used in diverse applications such as speech recognition [4, 5, 6], noise modeling , event detection , robot control [9, 10], as well as chaotic time-series prediction .
Many variations of the behavior and topology of reservoirs have been studied, including both fully connected  and sparse  random connections, or cyclical time-delay reservoirs [13, 14], and the reservoir nodes themselves may be analog, leaky  or spiking 
. The hallmark of all reservoir designs is that the internal layer is fixed and is not trained for specific datasets or applications. Non-reservoir recurrent-neural networks are typically trained by using a computationally expensive backpropagation method. By fixing the internal layer of a reservoir, all training occurs only at the output layer. Commonly the output layer consists of a linear readout of the internal states, where the readout weights are trained using a regularized least squares approach. Although this method is computationally fast, it may not produce the most accurate results, particularly in the presence of noisy data [18, 19]
. Specialized ESNs have been proposed to handle noisy data, in particular embedded in support vector machines or in a Bayesian framework , but done only in the context of time-series prediction tasks.
In contrast, this study compares the performance of comparatively basic reservoir output methods on noisy data for classification tasks. This side-by-side study of the basic approaches will guide the development of better specialized methods for specific applications. The three basic methods included within are:
Using regularized least squares to train the linear readout weights.
Using the Dantzig selector to train sparse, linear readout weights.
Finding a low-rank approximation of the reservoir states for each class to determine membership of test data.
Method A, using regularized least squares, was introduced in the seminal works, and still appears the most commonly in literature. Its appeal is that it typically produces good results and is computationally fast. However, it may be sensitive to reservoir parameters and noise. Three variants of Method A are considered: (1) computing a new output weight matrix for each timestep, (2) computing a single output weight matrix based on the final timestep, and (3) computing a single global output weight matrix to fit all timesteps.
Method B finds a sparse collection of readout weights. This approach may be beneficial when the reservoir has significantly more nodes than necessary to fully distinguish the different classes of the input data. Sparse readout weights were used with ESNs using computationally intense pruning methods  and efficient elastic nets for time-series forecasting , but does not appear to be a widely adopted approach. In this work, the sparse readout weights are obtained using the Dantzig selector , a regularized optimization problem similar to LASSO . It will be slightly slower to train than Method A, but will be equally fast in testing.
Finally, method C does not use a linear transformation of the internal layer data, but forms collections of low-rank principal components of the reservoir responses to the input data of each class. The reservoir responses of a new test signal are compared to these collections. The signal is assigned to the class whose principal components best describe the variation in the reservoir response. Using principal components to perform classification has been widely adopted in other contexts, but has only recently appeared for use with reservoirs. This method should perform very well with noisy data, but may be slower to produce results.
All methods are compared theoretically and experimentally, with experimental results performed on synthetic and real-world data with artificial noise applied. The results will show that in the noise-free case, the investigated methods typically yield similar classification accuracy. For ESNs with a small number of hidden-layer nodes, the regularized least squares approaches in Method A tend to give the most accurate results, although accuracy degrades quickly as the noise level increases. According to the experiments, better accuracy can be achieved on noisy data by increasing the number of nodes and using the low-rank approximation technique as in Method C. Only simple classification is performed; we do not employ ensemble classifiers.
This work is organized as follows: Section II describes the behavior of an ESN in detail, and discusses implementation of the three output methods explored in this work. Section III reports on the performance of the three methods on two different datasets. The article concludes with a summary of the experimental findings in Section IV.
The following notation is used. The set of nonnegative integers is denoted by . Let denote a vector, of length clear from the context, of all zeros except for a "1" in the entry. The norm of a matrix is denoted by , for
. The normal distribution with meanis given by
, and the uniform distribution on the intervalis given by . The transpose of a matrix is denoted by .
Ii ESN Models and Output Classifiers
The behavior of ESNs is described in detail here. Let be an input pattern, and let denote its value at time . Suppose is discretized and parameterized so . Suppose is fed into an random ESN with nodes, possibly leaky, with no output feedback. Then the internal-layer reservoir state at time is given by the vector , whose value can be determined by
Here, the ESN model appearing in  is followed. A representation of this design of an ESN is shown in Figure 1. Throughout this article we set and , since nonzero feedback is inconsistent with the third output classification method in this study. appearing in (2) may be chosen as either or a concatenation such as . For notational convenience, is used throughout this section, but the results hold similarly for with appropriate changes to matrix dimensions.
The parameters and are the feedback strength and input gain, respectively. The leaking rate is given by . The input weights are a fixed matrix, where is the dimension of the input patterns at each time step. The reservoir weights are a fixed matrix. The output weights are learned matrices, where is the number of classes of the input data.
Ii-a Output weights trained by regularized least squares
Traditionally in reservoir computing literature, input signals are classified by using a linear combination of the reservoir states, where the readout weights are trained using a least squares method [1, 12, 26, 27]. Three variations of the trained linear output weights are explored: temporal pointwise weights, single weights trained against the final timestep, and global weights trained against all timesteps. To determine the readout weights using any of these variations, first let
be a collection of distinct training patterns, each belonging to one ofclasses. Compute the reservoir states of each at time according to Equation (1). Collect these states into matrices by concatenating columns. Define a matrix indicating class membership of the training set. Say if is in the class and either (or possibly ) otherwise.
A1. Pointwise output weights
The input patterns may display large variation with respect to the time variable , even within the same class. In this case, a new output weight matrix should be computed for each . That is, the output weight matrices are determined so that . This can be done using least squares, but the result may be very sensitive to small purturbations . To obtain a more robust result, the regularized least squares method (3) is employed: For each , let be obtained via
The value appearing in Equation (3) is the regularization parameter, and helps to minimize overfitting of the training data. A larger will force dampening of the coefficients in . A smaller forces the approximation to be closer. An appropriate must be chosen to balance these two goals, and may be user-specified or determined using a cross-validation technique as in [20, 21]. Equation (3) can be solved explicitly by
After determining the collection of output weights, a newly encountered pattern is classified depending on how well the weights project reservoir states onto columns of . That is, let be the test pattern, and be its reservoir states determined by Equation (1). If belongs to the class, then theoretically for most . More precisely, let
Then is predicted to be in the class if .
This approach should work well on datasets with patterns that display similar intra-class temporal behavior, but has a somewhat high computational cost.
A2. Final output weights
In this approach, a single weight matrix is found based on the final temporal state of the reservoir. Suppose is the collection of final reservoir states for all patterns in the training set. The weight matrix is found so that . Here, the matrix is found using the following regularized optimization problem:
A newly encountered pattern is classified based on its final reservoir state . Let . Then is predicted to be in the class if .
This approach has a lower computational cost than the pointwise approach in (A1), and should work well if the final reservoir states are well separated for different classes. However, since classification is based on just a single point in time, accuracy may degrade if noise propagates through the reservoir and sufficiently contaminates the final reservoir states.
A3. Global output weights
In this approach, a single weight matrix is found based upon all of the reservoir responses. The weight matrix is found so that for all . Here, the matrix is found using regularized least squares:
A newly encountered pattern is classified similarly to the previous two approaches. Let be the reservoir state for the input at time , and let
Then is predicted to be in the class if .
In essence, this approach separates the average reservoir state over time of the classes. It should perform well on datasets with well separated class means and low intra-class variation.
Ii-B Sparse trained output weights
This method produces sparse linear output weights. The output weights determined by methods A1, A2 and A3 will in general be dense, with every reservoir node making a nonzero contribution to the output layer. However, it is reasonable to assume that only a subset of the reservoir nodes can be used for performing classification, especially if the number of nodes in the reservoir is very large. To this end, a sparse matrix recovery method may be used, such as LASSO  or the Dantzig selector , to determine sparse output weights .
We employ the Dantzig selector, which works as follows. Let be the collection of reservoir nodes of the training set, and let be the matrix indicating class membership. For each , compute
The value is a fixed regularization parameter and is a diagonal scaling matrix with entries equal to the norm of the rows of . Similar to Equation (3), a larger forces to be sparser with small nonzero coefficients, while a smaller makes the approximation better. The parameter must be chosen to balance accuracy and sparsity. Equations (3) and (7) are similar in that they both have terms to encourage close approximations, however, the regularization term of Equation (7) controls the overall sparsity of , rather than the regularization term of (3) that tends to make the entries of more uniform in magnitude. Equation (7) can be solved using the algorithms appearing in  or .
The test phase for this method follows that in Subsection II-A1.
Ii-C Principal components of reservoir states
This method performs classification based on the principal components of the reservoir states. This approach is motivated by the fact that often inputs from the same class have similar underlying reservoir responses. Since ESNs operate at the “edge of chaos" , perturbations of inputs within the same class may lead to noisy discrepancies in the corresponding reservoir states, however, a principal component based method will identify any strong unifying underlying reservoir response pattern if it exists.
The method employed below is a variant of that found in . In the training phase, for each time , create matrices each one equaling the collection of reservoir states from inputs in the class of the training set, determined using Equation (1). For some fixed positive integer , let be the collection of the first principal components of .
In the testing phase, let be a newly encountered input signal and its reservoir responses. If belongs to the class, then the columns of should capture the behavior of well at each . Therefore the class membership of may be predicted as follows. Let be a -vector with
Say is in the class if .
Since noise is a high-frequency component overlaying the true signal, by considering the low-rank approximations, this method should adapt to noisy data well.
Iii Experimental Results
In this section, the output classification methods are compared experimentally. Two datasets are used: A synthetic dataset with signals randomly alternating between sine waves and square waves, and a real-world dataset consisting of samples of speakers uttering the Japanese vowel ‘ae’.
The sine/square wave dataset has been used in reservoir computing experiments, notably in  and  for study with photonic reservoirs. The Japanese vowel dataset (JV) was used in  with ensemble classifiers built from several small leaky ESNs.
The output layer methods are explored using these datasets in a noise-free case, as well as when the data are modified by applying gaussian noise to the input for various levels . The reservoirs are initialized as in Equation (1) for various . For each pair , a number of simulations are performed, with any randomizations in the matrices and , the signals , or in the application of noise redetermined for each simulation.
For both datasets, two regularization parameters are used with the output layer method A1 of pointwise regularized least squares, and , which are labeled as ‘(A1) ’ and ‘(A1) ’ respectively in the tables and figures of results. For the output layer methods A2, A3, and B, the regularization parameter is used. These results are labeled as ‘(A2) Endpoints’, ‘(A3) Global’, and ‘(B) Sparse’ respectively. Finally, the results using method C using low-rank approximations are labeled as ‘(C) PCA’.
A small number of regularization parameters are used, providing a snapshot comparison of their performance for various noise levels. The results show that a larger regularization parameter tends to yield higher classification accuracy for data with a higher level of noise contamination.
The sine/square wave experiments are performed in MATLAB 2015b on a PC at 2.40 GHz and 16 GB RAM. The JV experiments are performed in MATLAB 2016a on the Condor Cluster hosted at the RIT DOD HPC Affiliated Resource Center.
Iii-a Sine vs. Square Wave
In this collection of experiments, the training dataset is a signal with randomly placed sine and square wave segments, along with an indicator function , where
A sample of a typical training set is shown in Figure 2. The individual sine and square wave segments have the same period. The reservoir responses are found as in Equation (1), with both and randomly choosing entries from , but scaled so that its spectral radius is less than . The parameters are chosen as , , and . The reservoir size is selected from .
The test set is generated by randomly combining the sine and square wave segments, then adding noise , for to the inputs of the test set. No artificial noise is applied to the training set.
|(A1)||(A1)||(A2) Endpoints||(A3) Global||(B) Sparse||(C) PCA|
Table I displays the mean and standard deviation of the test accuracy for the 50 simulations for each pair . At the noise-free level, methods A1, A2, A3, and C perform similarly, with the endpoint method A3 giving the best results. This is expected, due to the large discrepancy in the final value of the inputs for each class ( for square waves vs for sines). As the noise increases, the low-rank method C outperforms the others as expected. The performance of the sparse method B improves as the reservoir size increases, for then it has more dimensions from which to choose the few most significant nodes contributing to classification. The performance of the other methods do not depend on
, probably due to the simple behavior of the input patterns. Finally, note that the single global output weight method B performs poorly. This is likely explained by the fact that both the sine and square waves have an average value of 0.
Iii-B Japanese Vowels
In this collection of experiments, speaker identification is performed using the Japanese vowel dataset. The dataset originally appeared in  and was obtained via . JV is a temporal set with samples recorded from nine male speakers saying the Japanese vowel ‘ae’. Each sample forms an real-valued matrix, where
is an integer between 7 and 29 (inclusive) depending on the duration of the utterance, and 12 is the number of cepstrum coefficient features observed for each speaker. The dataset includes 640 samples; of which a fixed selection of 270 are used for training (with 30 utterances by each of the nine speakers), and the remaining are used for testing (24-88 utterances by each of the nine speakers). The test dataset was modified by adding gaussian white noise, with .
Results were obtained using a simulated ESN, with input nodes accounting for the 12 cepstrum coefficients plus two bias terms, output nodes representing the 9 classes (one for each speaker), and N reservoir nodes with . The ESN parameters were fixed at , and , chosen to match those appearing in .
Table II displays the mean and standard deviation of the test classification accuracy over 100 simulations for each pair . For low noise levels, the methods perform similarly, especially for larger , with the endpoint method A2 outperforming the rest. As expected, the endpoint method A2 appears to be the least tolerant to noise. The single global weight method A3 is somewhat tolerant to noise corruption, but the low-rank method C is shown to be the most robust to noise. The table suggests that to achieve the best classification accuracy on low noise data, one should use the endpoint method A2 with a moderate reservoir size, and to achieve the best accuracy on higher noise data, one should use the low-rank approach C with a large reservoir size.
|(A1)||(A1)||(A2) Endpoints||(A3) Global||(B) Sparse||(C) PCA|
Figures 3 and 4 display more detailed classification accuracy results. Figure 3 fixes the reservoir size at , and displays the mean and standard deviation of the test accuracy over all 100 simulations at each noise level . Figure 4 fixes the noise level at , and displays the mean test accuracy at each reservoir size. The right-hand plot of Figure 4 displays the ‘boxed’ area of the left-hand plot, and includes standard deviation information as vertical bars. The figures show that the low-rank approximation method C performs well, tending to have the highest classification accuracy and lowest variation of results. Notice that most of the methods have similar classification accuracy for small reservoir sizes as in the right-hand plot of Figure 4, but the low-rank method C and single global weights method A3 have much lower variation than the others.
In this study, several approaches for performing classification with ESNs were compared. The numerical experiments suggest that the classification accuracy depends strongly on the characteristics of the data, the classification method used, and the design of the ESN. It was observed that in both synthetic and real-world datasets with low noise corruption, employing an output method using linear readout weights trained on the final temporal state of a reservoir, as in Equation (5), with moderate number of nodes produces the most accurate results. On the other hand, if the data are noisy then the findings suggest that one should design a reservoir with a larger number of nodes and use the classification method of finding a low-rank approximation of the reservoir states, as in Equation (8
). In general, we did not observe evidence to support using the classification method using sparse linear readout weights in practice. Future directions include investigating on additional datasets and with other reservoir types and output methods, including time delay reservoirs, multilayer perceptrons, and a global variation of the low-rank approximation method, rather than the pointwise approach studied here.
This research was partially supported by Air Force Office of Scientific Research [LRIR:15RICOR122].
Approved for public release by USAF 88 ABW on 18 NOV 2016. Case Number: 88ABW–2016–5889. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the view of the United States Air Force.
-  H. Jaeger, The ‘echo state’ approach to analysing and training recurrent neural networks - with an erratum note, Tech. Rep. GMD Report Number 148, Fraunhofer Institute for Autonomous Intelligent Systems (2011).
-  W. Maass, H. Natschläger, H. Markram, Real-time computing without stable states: a new framework for neural computation based on perturbations, Neural Comput. 14 (2002) 2531–2560. doi:10.1162/089976602760407955.
-  D. Verstraeten, B. Schrauwen, M. D’Haene, D. Stroobandt, An experimental unification of reservoir computing methods, Neural Netw. 20 (2007) 391–403. doi:10.1016/j.neunet.2007.04.003.
-  W. Maass, H. Natschläger, H. Markram, A model for real-time computation in generic neural microcircuits, in: Proceedings of NIPS, 2003, pp. 229–236.
-  M. D. Skowronski, K. Meier, J. Schemmel, Minimum mean squared error time series classification using an echo state network prediction model, in: Proceedings of ISCAS, 2006.
-  D. Verstraeten, B. Schrauwen, D. Stroobandt, J. Van Campenhout, Isolated word recognition with the liquid state machine: A case study, Inform. Process. Lett. 96 (2005) 521–528.
-  H. Jaeger, H. Haas, Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communications, Science 304 (2004) 78–80, doi: 10.1126/science.1091277.
J. Hertzberg, H. Jaeger, F. Schönherr, Learning to ground fact symbols in behavior-based robots, in: ECAI’02 Proceedings of the 15th European Conference on Artificial Intelligence, 2002, pp. 708–712.
-  H. Burgsteiner, On learning with recurrent spiking neural networks and their applications to robot control with real-world devices, Ph.D. thesis, Graz University of Technology (2005).
-  P. Joshi, W. Maass, Movement generation and control with generic neural microcircuits, in: Biologically Inspired Approaches to Advanced Information Technology, Vol. 3141, 2004, pp. 258–273.
H. Jaeger, Long short-term memory in echo state networks: Details of a simulation study, Tech. Rep. 27, Jacobs University Bremen (2012).
-  M. Lukoševičius, A practical guide to applying echo state networks, in: G. Montavon, et al. (Eds.), NN: Tricks of the Trade, 2nd Edition, Springer-Verlag Berlin Heidelberg, 2012, pp. 650–686.
-  L. Grigoryeva, J. Henriques, L. Larger, J. Ortega, Optimal nonlinear information processing capacity in delay-based reservoir computers, Sci. Rep. 5 (12858). doi:10.1038/srep12858.
-  Y. Paquot, F. Duport, A. Smerieri, J. Dambre, B. Schrauwen, M. Haelterman, S. Massar, Optoelectronic reservoir computing, Sci. Rep. 2. doi:10.1038/srep00287.
H. Jaeger, M. Lukosevicius, D. Popovici, U. Siewert, Optimization and applications of echo state networks with leaky-integrator neurons, Neural Netw. 20 (3) (2007) 335–352.
-  A. Almassian, Information representation and computation of spike trains in reservoir computing systems with spiking neurons and analog neurons, Ph.D. thesis, Portland State University (2016).
-  P. Werbos, Backpropagation through time: What it does and how to do it, Proc. IEEE 78 (1990) 1550–1560. doi:10.1109/5.58337.
-  A. Prater, Reservoir computing for spatiotemporal signal classification without trained output weights, arXiv: 1604.03073 (2016).
-  L. N. Trefethen, D. Bau, Numerical Linear Algebra, SIAM, 1997.
-  Z. Shi, M. Han, Support vector echo-state machine for chaotic time-series prediction, IEEE Trans. Neural Netw. 18, doi: 10.1109/TNN.2006.885113.
-  D. Li, M. Han, J. Wang, Chaotic time series prediction based on a novel robust echo state network, IEEE Trans. Neural Netw. Learn. Syst. 23, doi: 10.1109/TNNLS.2012.2188414.
-  X. Dutoit, B. Schrauwen, J. Van Campenhout, D. Stroobandt, H. Van Brussel, M. Nuttin, Pruning and regularization in reservoir computing, Neurocomputing 72 (2009) 1534–1546, doi: 10.1016/j.neucom.2008.12.020.
-  F. M. Bianchi, S. Scardapane, A. Uncini, A. Rizzi, A. Sadeghian, Prediction of telephone calls load using echo state networks with exogenous variables, Neural Netw. 71 (2015) 204–214, doi: 10.1016/j.neunet.2015.08.010.
E. Candes, T. Tao, The Dantzig selector: Statistical estimation whenis much larger than , Ann. Statist. 35 (2007) 2392–2404.
-  R. Tibshirani, Regression shrinkage and selection via the lasso, J. Royal Statist. Soc. B 58 (1996) 267–288.
-  A. Goudarzi, P. Banda, M. Lakin, C. Teuscher, D. Stefanovic, A comparitive study of reservoir computing for temporal signal processing, Tech. rep., University of New Mexico (2014).
-  M. Lukoševičius, H. Jaeger, Reservoir computing approaches to recurrent neural network training, Comp. Sci. Rev. 3 (2009) 127–149.
-  A. Prater, L. Shen, B. Suter, Finding Dantzig selectors with a proximity operator based fixed-point algorithm, Comput. Stat.Data Anal. 90 (2015) 36–46, doi: 10.1016/j.csda.2015.04.005.
-  Z. Lu, T. K. Pong, Y. Zhang, An alternating direction method for finding Dantzig selectors, Comput. Stat. Data Anal. 56 (2012) 4037–4046, doi: 10.1016/j.csda.2012.04.019.
-  R. Legenstein, W. Maass, Edge of chaos and prediction of computational performance for neural circuit models, Neural Netw. 20 (2007) 323–334. doi:10.1016/j.neunet.2007.04.017.
-  H. Zhang, X. feng, B. Li, Y. Wang, K. Cui, L. F., W. Dou, Y. Huang, Integrated photonic reservoir computing based on hierarchical time-multiplexing structure, Opt. Express 22 (2014) 31356–70, doi: 10.1364/OE.22.031356.
-  M. Kudo, J. Toyama, B. Li, Multidimensional curve classification using passing-through regions, Pattern Recogn. Lett. 20 (1999) 1103–1111.
-  M. Lichman, Japanese vowels data set, http://archive.ics.uci.edu/ml/datasets/Japanese+Vowels, accessed: July 2016.