Classification via Tensor Decompositions
of Echo State Networks
This work introduces a tensor-based method to perform supervised classification on spatiotemporal data processed in an echo state network. Typically when performing supervised classification tasks on data processed in an echo state network, the entire collection of hidden layer node states from the training dataset is shaped into a matrix, allowing one to use standard linear algebra techniques to train the output layer. However, the collection of hidden layer states is multidimensional in nature, and representing it as a matrix may lead to undesirable numerical conditions or loss of spatial and temporal correlations in the data.
This work proposes a tensor-based supervised classification method on echo state network data that preserves and exploits the multidimensional nature of the hidden layer states. The method, which is based on orthogonal Tucker decompositions of tensors, is compared with the standard linear output weight approach in several numerical experiments on both synthetic and natural data. The results show that the tensor-based approach tends to outperform the standard approach in terms of classification accuracy.
Echo State Networks (ESNs), first introduced in  and  under the name Liquid State Machines, have been shown to be effective at performing classification and temporal prediction tasks on spatiotemporal data, including such diverse tasks as speech recognition [20, 26, 31, 11, 24], chaotic time-series prediction , and forming the objective function in reinfocement learning methods [3, 28]. ESNs are a special type of recurrent neural network (RNN). Typically the hidden layer in an RNN is trained using a computational expensive back propagation method . In ESNs however, the weights of the hidden layer are randomly assigned with sparse, random connections among the nodes. The hidden layer of an ESN is often called a reservoir, with the collection of node values called the reservoir states. The advantages ESNs offer over traditional RNNs include a much faster training time, and a configuration that does not require retraining to use in new applications on different datasets.
In ESNs, only the output layer is trained for a particular task, and the training generally produces linear output weights using a regularized least squares approach [6, 9, 18, 19]. Recently, the output layer of an ESN was replaced with a classification scheme based on the principal components of the reservoir states [23, 24]. This approach showed promise in improving classification accuracy and in being more robust to noisy perturbations in the input data than the traditional linear trained output weights approach. However, the collection of reservoir states generated by an ESN is multidimensional in nature. Both trained linear output weights and the principal components approaches require one to superficially flatten the reservoir data into matrices, potentially eliminating or weakening spatial and temporal correlations present in the data in the process. Additionally, flattening the reservoir data may result in a very overdetermined system, especially for the trained linear output weights approach, which may yield overtrained or overly sensitive results. For example, see .
To this end, this work proposes a tensor-based supervised classification method for use with reservoir states of ESNs. Tensors, or multidimensional arrays, are natural structures for representing collections of reservoir data. Rather than using the raw reservoir data, the tensors will be approximated using a Tucker decomposition , with each mode having a factor matrix of smaller rank along with an associated core tensor that extracts the higher order features of the data. This decomposition enables one to reduce the complexity of the training tensor data by keeping only the most significant contributions while preserving the multidimensional correlations among the features.
The following notation appears in this work. Tensors will be written in a script capital letters, e.g. , matrices will appear as capital latin letters, e.g. , and vectors and scalars will be lower case latin or greek letters, e.g. or . Elements of an array will be given in ‘matlab’ notation. Thus the element of matrix , that is the element in the row and column, will be denoted by . The vector determined by extracting the column from is denoted by . Similar notation holds for lower or higher-order arrays. The usual matrix-matrix multiplication will be represented by writing the matrices adjacent to one another. This is in contrast to modal tensor-matrix multiplication, which will be defined in Section II. Superscripts will denote indices in a set, not power-wise multiplication, except for when the base is a set as in . The Hermitian transpose of a matrix is given by . Finally, will denote a vector, of length clear from context, with a ‘1’ in the k-th position and zeros elsewhere.
The remainder of this paper is organized as follows. Section II discusses background information for ESNs and relevant tensor decompositions. Section III describes the proposed classification method using tensor decompositions on the reservoir states. The results of several numerical experiments are presented in Section IV. Finally, Section V contains conclusions and a dicussion of future work.
In this section, background information is given on tensor decompositions and echo state networks.
Ii-a Tensor Decompositions
A tensor is a higher-order analogue of a vector; A vector is a first-order tensor, a matrix is a second-order tensor, and so on. Tensors are natural structures to represent and investigate multidimensional data. For example, video may be considered a third-order tensor, with the first two modes describing the and pixel coordinates of a single frame, and the third mode representing time variations. Tensors have been used to represent and explore relationships among data in diverse research areas and applications, including video processing , multiarray signal processing , independent component analysis , and others.
One challenge to employing tensor methods is the volume of the data, which suffers from the so-called ‘curse of dimensionality’ . That is, the amount of data increases exponentially with each additional mode, and naive tensor methods quickly become intractable. To help alleviate this challenge various approaches have been proposed, including several types of low-rank and sparse representations . In this work, we will employ the orthogonal Tucker-2 decomposition.
Before discussing the decomposition, some definitions need to be introduced. Let be an -order tensor, with mode having dimension , and let be a matrix. The mode product of by is the -order tensor whose entries are given by
Modal products of a tensor by matrices are commutative, provided the modes are distinct. That is, given the tensor and the matrices and , then
A tensor may be represented as a matrix through the process of unfolding [12, 13], similar to the way a matrix may be represented as a vector by vectorization. The matrix unfolding of in the third mode is the matrix with elements
One may think of as the concatenated matrix
Unfoldings of higher-order tensors or in different modes follow an analogous procedure.
An inner product may be defined on tensors as follows. Let
be defined by
This inner product induces the tensor Frobenius norm,
Now equipped with the requisite definitions, we are ready to build the orthogonal Tucker-2 tensor decomposition.
where and are the factor matrices, typically with , and is the core tensor. If the matrices and each have orthogonal columns, then (1) is called an orthogonal tucker decomposition, and is expressed as
A variant of the orthogonal Tucker decomposition, called the Tucker-2 decomposition, is of the form
with core tensor . That is, the original tensor of order greater than two is written using only two factor matrices. An illustration of the decomposition (3) is shown in Figure 1. This type of decomposition was used in  to perform feature extraction and classification of a collection of 2-dimensional signals, where the third mode of the tensors and correspond to the individual samples in the training dataset. Note that in the decomposition (3), the basis matrices and are universal across the entire tensor . In classification tasks, this method does not generate different basis matrices for distinct classes in the data set. The decomposition (3) and classification methods using it will be discussed further in Section III for use with Echo State Network data.
Ii-B Echo State Networks
Spatiotemporal data are processed in an echo state network as follows. Let be an input of spatial dimension and temporal length . The values of the hidden layer nodes, also referred to as the reservoir states of the input , denoted by , are determined by the recursion
In (4), is a nonlinear activation function, are the fixed input weights, are the fixed, randomly assigned reservoir weights, is a bias value, and is the leaking rate. This work does not use output feedback in the ESN recursion as is sometimes used in literature , since it is incompatible with using the proposed tensor approach.
For supervised classification tasks, one will have a collection of training samples. Suppose the training data are partitioned into classes, and the class contains examples. Say
is the entire collection of training inputs and
is the collection of training inputs from the class. Suppose each of the inputs are processed in the ESN (4) with the same weights and parameters. Denote the reservoir states from the input as . All of the reservoir states from may be concatenated along the third mode into the tensor , where
A representation of the tensor is shown in Figure 2. Similarly, all of the reservoir states may be concatenated along the third mode into the tensor .
Traditionally, a collection of linear output weights is trained on the unfolded collection of training reservoir tensors. Let be the unfolding of along the third mode. That is, may be written as the contatenation
Let be a matrix whose columns indicate class membership. Say the columns of satisfy .
In the standard approach  to performing classification using ESNs, first ones finds a collection of output weights so that
Commonly a solution to (5) is found using a regularization method, such as
If a new input belongs to the class at time , and describes the data well, then . This observation drives the classification scheme. Say the input from the test set with reservoir states is predicted to belong to class at time if the the maximal element of the vector is in the entry. Similarly, is predicted to belong to the class overall if the vector
is maximized in the entry. Note that a single output weight matrix (6) is found and applied to all test samples at all time steps.
Training output weights as in (6) is fast, however this approach does have weakensses. The system (5) is typically overdetermined. Although the regularization (6) tries to overcome this, the matrix may have several orders of magnitude more columns than rows in practice and may poorly represent the data, even with regularization. This may be controlled by the model by selecting only a subset of times at which to sample the reservoir. One common method is to use only a single point , however the accuracy results may suffer greatly . Even though the reservoir states hold some ‘memory’ of previous states, using only a subset of the data generally results in information loss. Moreover, the unfolding procedure loses some temporal correlations of reservoir node behavior. Finally, using a single linear output weight may simply be insufficient to separate the classes in the dataset well, as shown in .
Iii Tensor Decompositions of Reservoir States
To alleviate the deficiencies encountered using method (7), this paper proposes using a classification method based on the decomposition of the tensor of training reservoirs. To this end, let approximations of the orthogonal Tucker decompositions of the tensors and be
where the factors have orthogonal columns with and , and and are core tensors. The core tensors may be intrepreted as the entry describing the strength of the feature in the reservoir states of the input captured by the interaction of the bases and .
To approximate the Tucker-2 decomposition of the form (8), we use the Higher-Order Orthogonal Iteration (HOOI) algorithm, first introduced in  and explored in . For completeness, we include the pseudocode as Algorithm 1. A similar procedure can be used to obtain the decompositions (9).
In , the factor matrices are initialized as the dominant left singular subspace of the unfolded tensors and . In this work, we randomly initialize them for two reasons. First, just finding the dominant subspaces of the unfolded matrices is a very computationally intensive task if is large. In practice, this step may be intractable, even if steps 6 and 11 in Algorithm 1 are computable. Second, it was noted in experiments that randomly initializing the factors results in only a few additional iterations. Algorithm 1 uses a stopping criterion based on the convergence of the singular values found in Steps 6 and 11, this is to avoid problems in difference in signs when using a criterion based on the factor matrices.
One may perform classification using the decompositions (8) and (9) obtained from collections of reservoir states of a training set. To do so, suppose is a new input signal with reservoir states . Although is a matrix rather than a three dimensional tensor, it can still be expressed in terms of the factor matrices from the Tucker decompositions. Say
Indeed, each frontal slice of is of the form (12), with a collection of reservoir states from the training set in place of . It is reasonable to assume that inputs from the same class have similar reservoir responses, and therefore also produce similar frontal slices in the core tensors. Therefore, one may predict that an input belongs to the class if the slices from describe well. That is, say is in the class if
and . Similarly, one could predict belongs to the class if
Iv Experimental Results
In this section, the results of numerical experiments are given, comparing the classification accuracy using ESNs with the linear output weight approach (6) with the proposed tensor-based classification methods (14), and (15). Three datasets are used. The first dataset uses inputs that randomly switch between sine wave and square wave segments. The second dataset is a subset of the USPS collection of handwritten digits. The final dataset is a collection of cepstrum coefficients from audio recordings of speakers saying the Japanese vowel ‘ae’.
All experiments are performed in MATLAB 2017a on a PC with 16GB RAM. Several parameter combinations are considered for each dataset, with experiments repeated a number of times for each combination. The randomizations in the weights and and in generating training and testing datasets are reselected for each experiment, however they are held constant for all classification methods within a single experiment.
Iv-a Sine vs. Square Wave
In this collection of experiments, the input signals are formed by randomly placing sine and square wave segments of length and period 100, and paired with an indicator matrix where
A sample training input is shown in Figure 3. A dataset of this type has been studied in previous ESN work, including  for use with ESN matrix principal component output layer classification methods, and [21, 33] for study using photonic reservoirs.
The training data is formed by generating 20 distinct input patterns, each containing 100 segments of randomly placed sine and square waves. The test set is generated similarly, but contains 50 distinct patterns. The ESN parameters used in the simulations are and . For each , the experiments are repeated 50 times with new randomizations in and the training and testing sets. For tensor-based classification, we used the method (15) with ranks . The choice is justified because the the dataset is rather simple; only one type of square wave and one type of sine wave are used in generating the samples.
For tensor-based classification, the data generated by the ESNs on the training set are partitioned into two tensors and , where contains the reservoir states corresponding to the ‘sine’ inputs, and contains the reservoir states corresponding to the ‘square’ inputs. Each tensor is decomposed as in (9):
For each test set element, the segments of length 100 are classified according to (15). That is, let be a new input pattern from the test set with reservoir states . Let be the collection of states corresponding to the input segment of . Say that the segment is classified as a ‘sine’ wave if
and as a ‘square’ wave segment if the inequality sign is flipped.
For trained linear output weights-based classification, we generate a single matrix via Equation (6), and perform classification on the test set both pointwise and block-wise on each input segment as in Equation (7).
The mean and standard deviation of the percent classification accuracy results using these methods over 50 simulations for several parameter choices are displayed in Table I. The tensor-based classification method, in columns labeled ‘Tensor’, achieved 100% accuracy in every simulation on both the training and testing datasets. The pointwise and block-wise output weight classification methods, in columns labeled ‘Weights (pt)’ and ‘Weights (bk)’ respectively, achieved good classification accuracy for some parameter choices, but poor results for others. The block-wise method is sensitive to the number of nodes in the ESN and the bias choice, and in particular achieved near-perfect test set classification accuracy when was chosen well. On the other hand, the pointwise output weights method achieved near perfect results only when was chosen well and was sufficiently large.
|N||Weights (pt)||Weights (bk)||Tensor||Weights (pt)||Weights (bk)||Tensor|
|10||52.01 (0.75)||52.00 (10.88)||100.00 (0.00)||51.88 (0.67)||49.28 (6.77)||100.00 (0.00)|
|20||53.22 (1.05)||48.90 (13.03)||100.00 (0.00)||53.11 (0.59)||47.68 (7.88)||100.00 (0.00)|
|50||56.41 (1.56)||60.70 (17.41)||100.00 (0.00)||56.51 (0.95)||61.96 (14.64)||100.00 (0.00)|
|10||88.40 (1.80)||100.00 (0.00)||100.00 (0.00)||87.99 (1.91)||100.00 (0.00)||100.00 (0.00)|
|20||91.35 (1.72)||100.00 (0.00)||100.00 (0.00)||91.09 (2.05)||100.00 (0.00)||100.00 (0.00)|
|50||99.51 (0.10)||100.00 (0.00)||100.00 (0.00)||99.49 (0.20)||99.48 (3.68)||100.00 (0.00)|
|10||51.77 (0.57)||48.40 (11.49)||100.00 (0.00)||51.86 (0.58)||50.44 (7.67)||100.00 (0.00)|
|20||53.09 (0.79)||48.40 (10.62)||100.00 (0.00)||53.27 (0.86)||50.52 (7.45)||100.00 (0.00)|
|50||56.25 (1.66)||61.90 (18.98)||100.00 (0.00)||56.39 (1.11)||63.88 (15.10)||100.00 (0.00)|
|10||88.00 (2.13)||100.00 (0.00)||100.00 (0.00)||87.25 (3.62)||100.00 (0.00)||100.00 (0.00)|
|20||92.24 (1.50)||100.00 (0.00)||100.00 (0.00)||92.55 (1.53)||100.00 (0.00)||100.00 (0.00)|
|50||99.51 (0.11)||100.00 (0.00)||100.00 (0.00)||99.46 (0.14)||100.00 (0.00)||100.00 (0.00)|
Although (15) is an instance-based classifier, 100% classification accuracy on the test set is not guaranteed. Individual sine or square segments are indeed identical whether from the training or testing set. However, the resulting reservoir states from these segments are all distinct due to the memory of the reservoir. That is, the reservoir is not in a resting state when accepting segments in a sequence, and the initial state will continue propagating through the reservoir for some time. The ouput weight includes this contamination when it is trained, however the tensor decomposition method can capture the contribution from the initial state in only a small number of factors, while focusing primarily on the reservoir behavior stemming from the input itself. Overall, the tensor-based approach outperformed the output weights method in that it achieved a higher classification accuracy on the test set for all parameter choices.
Iv-B USPS Handwritten Digits
In this collection of experiments, classification is performed on grayscale images of handwritten digits. The dataset is partitioned into 10 classes, representing digits ‘0’ through ‘9’. Some samples of these images are shown in Figure 4. The images are treated as spatiotemporal signals, with coordinates corresponding to the spatial dimension and coordinates corresponding to the temporal dimension.
For the training dataset, 100 images from each class are randomly selected, and paired with an indicator matrix where if the input belongs to the class. The test set is formed similarly, but with a distinct 100 images from each class so the test set and training set have no overlap in samples.
The ESN parameters used are , , , and the modal ranks for the tensor-based approach are and . For each triplet , forty simulations are performed with different randomizations in the training and test set selections, as well as the ESN input and reservoir weight matrices for each simulations. However, for each simulation the same selections are used with each of the classification methods.
For tensor-based classification, , the tensor of reservoir states of the training inputs, is decomposed via Algorithm 1 into the form
is minimized for some training input in the class.
For trained linear output weights-based classification, a single output weight matrix is generated as in Equation (6) for each simulation. Classification is performed on the entire collection of reservoir states for each test input, as in Equation (7).
The results of these simulations are presented in Figure 1. Results using the HOOI Algorithm 1 are displayed in blue, and results using linear output weights are shown in red. The maximal mean accuracy over all pairs are shown for each for the tensor-based approach. Note that the tensor-based approach consistently yields higher classification accuracy than the trained linear output weight approach. Although the results are not competitive with state-of-the-art on this particular dataset, they do show that the tensor-based classification method yields higher results than standard ESN techniques.
Iv-C Japanese Vowels
In this collection of experiments, speaker identification is performed using samples of audio recordings. The dataset contains 640 samples of nine male speakers saying the Japanese vowel ‘ae’. The data is split into 270 training samples of 30 utterances by each speaker, and 370 test samples of 24-88 utterances by each speaker. Each sample is an array of cepstrum coefficients, where is the temporal length of the sample, using 12 cepstrum coefficients and two bias terms. The dataset, first appearing in  and obtained via , is popular in machine learning and ESN literature [1, 8, 11, 24, 27]. A test accuracy of 100% was reported in  using an ensemble classifier of 1000 ESNs of four leaky-integrator nodes.
In the examples below, we use a single ESN with or nodes, a nonlinear activation function and bias . The classifiers are found using (6) and Algorithm 1 with (8) from the collection of reservoir states corresponding to the training inputs. The test inputs are modified by adding Gaussian noise with . Classification is then performed on the resulting reservoir states using (7) for trained linear output weights or (14) for the tensor method. For each pair , 20 simulations were performed with new randomizations of and for each simulation.
The test classification accuracy results are displayed in Figure 6. In the figure, the blue lines correspond to the tensor-based method, and the red lines correspond to the linear output weight method. The individual lines within each method correspond to different levels of noise added to the test inputs. The -axis is the number of reservoir nodes . The points on the lines give the mean accuracy over the 20 simulations, while the vertical lines represent one standard deviation from the mean.
In the figure, classification accuracy initially decreases as increases, but eventually improve for large enough . This is consistent with results published in . Both methods degrade as the level of added noise increases, however the tensor-based method consistently yielded better accuracy results. Not only does the tensor-based approach have higher mean classification accuracy for all parameter choices, but the standard deviation is smaller indicating that the results are less sensitive to the randomizations in and .
This work introduced a tensor-based method to perform supervised classification on spatiotemporal data processed in an ESN. The numerical experiments demonstrate that the proposed method may outperform the traditional trained linear output weights method in terms of classification accuracy. Future directions include investigating other types of tensor decompositions, including low rank polyadic and sparse decompositions, as well as using other types of non instance-based classifiers on the resulting decompositions.
This work was cleared for public release by Wright Patterson Air Force Base Public Affairs on 15 Aug 2017. Case Number: 88ABW-2017-3910.
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.
-  D. Barber. Dynamic Bayesian networks with deterministic latent tables. In Proc. NIPS 2003, 2003.
-  R. Bellman. Dynamic Programming. Princeton University Press, Princeton, NJ, USA, 1 edition, 1957.
-  K. Bush and C. Anderson. Modeling reward functions for incomplete state representations via echo state networks. In 2005 IEEE International Joint Conference on Neural Networks (IJCNN), 2005.
-  A. Cichocki, D. P. Mandic, A. H. Phan, C. F. Caiafa, G. Zhou, Q. Zhao, and L. De Lathauwer. Tensor decomopsitions for signal processing applications: From two-way to multiway component analysis. IEEE Signal Processing Magazine, 32(2):145–163, 2015.
-  L. De Lathauwer, B. De Moor, and J. Vandewalle. On the best rank- and rank- approximation of higher-order tensors. SIAM J. Matrix Anal. Appl., 21(4):1324–1342, 2000.
-  A. Goudarzi, P. Banda, M. Lakin, C. Teuscher, and D. Stefanovic. A comparitive study of reservoir computing for temporal signal processing. Technical report, University of New Mexico, 2014.
-  N. Goyal, S. Vempala, and Y. Xiao. Fourier PCA and robust tensor decomposition. ACM, 2014.
-  P. Guerts. Pattern extraction for time series classification. In Proc. PKDD 2001, pages 115–127, 2001.
-  H. Jaeger. The ‘echo state’ approach to analysing and training recurrent neural networks - with an erratum note. Technical Report GMD Report Number 148, Fraunhofer Institute for Autonomous Intelligent Systems, 2011.
-  H. Jaeger. Long short-term memory in echo state networks: Details of a simulation study. Technical Report 27, Jacobs University Bremen, 2012.
-  H. Jaeger, M. Lukosevicius, D. Popovici, and U. Siewert. Optimization and applications of echo state networks with leaky-integrator neurons. Neural Networks, 20(3):335–352, 2007.
-  T. G. Kolda and B. W. Bader. Tensor decompositions and applications. 51(3):455–500, 2009.
-  P. Kroonenberg. Three-mode principal component analysis: Theory and applicaitons. DSWO Press, 1983.
-  M. Kudo, J. Toyama, and B. Li. Multidimensional curve classification using passing-through regions. Pattern Recognition Letters, 20:11–13.
-  Q. Li, A. Prater, L. Shen, and G. Tang. Overcomplete tensor decomposition via convex optimization. In Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), IEEE International Workshop on, page 53, 2015.
-  M. Lichman. Japanese vowels data set. http://archive.ics.uci.edu/ml/datasets/Japanese+Vowels, Accessed: July 2016.
-  L.-H. Lim and P. Comon. Multiarray signal processing: Tensor decomposition meets compressed sensing. Comptes Rendus Mecanique, 338(6):311–320, 2010.
-  M. Lukoševičius. A practical guide to applying echo state networks. In G. Montavon and et al., editors, NN: Tricks of the Trade, pages 650–686. Springer-Verlag Berlin Heidelberg, 2 edition, 2012.
-  M. Lukoševičius and H. Jaeger. Reservoir computing approaches to recurrent neural network training. Comp. Sci. Rev., 3:127–149, 2009.
-  W. Maass, H. Natschlager, and H. Markram. Real-time computing without stable states: a new framework for neural computation based on perturbations. Neural Comput., 14:2531–2560, 2002.
-  Y. Paquot, F. Duport, A. Smerieri, J. Dambre, B. Schrauwen, M. Haelterman, and S. Massar. Optoelectronic reservoir computing. Sci. Rep., 2, 2012.
-  A. H. Phan and A. Cichocki. Tensor decompositions for feature extraction and classification of high dimensional datasets. Nonlinear theory and its applications, IEICE, 1(1):37–68, 2010.
-  A. Prater. Comparison of echo state network output layer classification methods on noisy data. In 2017 International Joint Conference on Neural Networks (IJCNN), 2017.
-  A. Prater. Spatiotemporal signal classification via principal components of reservoir states. Neural Networks, 91:66–75, 2017.
-  P. Shah, N. Rao, and G. Tang. Sparse and low-rank tensor decomposition. In Advances in Neural Information Processing Systems (NIPS), Proceedings of, 2015.
-  M. D. Skowronski, H. Natschlager, and H. Markram. Minimum mean squared error time series classification using an echo state network prediction model. In Proceedings of ICSAS, 2006.
-  M. Strickert. Self-organizing neural networks for sequence processing. PhD thesis, Univ. of Osnabrück, Department of Computer Science, 2004.
-  I. Szita, V. Gyenes, and A. Lorincz. Reinforcement learning with echo state networks. In ICANN, Part I LNCS, pages 830–839, 2006.
-  L. R. Tucker. Some mathematical notes on three-mode factor analysis. Psychometrika, 31(3):279–311, 1966.
-  K. Vandoorne and et al. Toward optical signal processing using photonic reservoir computing. Opt. Express, 16:11182–92, 2008.
-  D. Verstraeten, B. Schrauwen, D. Stroobandt, and J. Van Capenhout. Isolated word recognition with the liquid state machine: A case study. Information Processing Letters: Special issue on applications of spiking neural networks, 96:521–528, 2005.
-  P. Werbos. Backpropagation through time: What it does and how to do it. Proc. IEEE, 78:1550–1560, 1990.
-  H. Zhang, X. Feng, B. Li, and et al. Integrated photonic reservoir compuing based on hierarchical time-multiplexing structure. Opt. Express, 22, 2014.