Reservoir computing approaches for representation and classification of
multivariate time series
Filippo Maria Bianchia,∗, Simone Scardapaneb , Sigurd Løksea , Robert Jenssena
a Machine

arXiv:1803.07870v1 [cs.NE] 21 Mar 2018

b Department

Learning Group, UiT the Arctic University of Norway, Hansine Hansens veg 18, 9019 Tromsø, Norway.
of Information Engineering, Electronics and Telecommunications (DIET), Sapienza University of Rome, Via
Eudossiana 18, 00184 Rome, Italy

Abstract
Classification of multivariate time series (MTS) has been tackled with a large variety of methodologies and
applied to a wide range of scenarios. Among the existing approaches, reservoir computing (RC) techniques,
which implement a fixed and high-dimensional recurrent network to process sequential data, are computationally efficient tools to generate a vectorial, fixed-size representation of the MTS, which can be further
processed by standard classifiers. Building upon previous works, in this paper we describe and compare
several advanced RC-based approaches to generate unsupervised MTS representations, with a specific focus
on their capability of yielding an accurate classification. Our main contribution is a new method to encode
the MTS within the parameters of a linear model, trained to predict a low-dimensional embedding of the
reservoir dynamics. We also study the combination of this representation technique when enhanced with
a more complex bidirectional reservoir and non-linear readouts, such as deep neural networks with both
fixed and flexible activation functions. We compare with state-of-the-art recurrent networks, standard RC
approaches and time series kernels on multiple classification tasks, showing that the proposed algorithms
can achieve superior classification accuracy, while being vastly more efficient to train.
Keywords: Reservoir computing, model space, time series classification, recurrent neural networks

1. Introduction
In this paper we consider the problem of classifying multivariate time series (MTS), i.e., assigning each
MTS to one of a fixed number of classes. This is a fundamental task in many applicative scenarios, including
(but not limited to) health monitoring [1, 2, 3, 4], civil engineering [5], action recognition [6], and speech
analysis [7, 8]. The problem has been tackled by a wealth of different approaches, spanning from the definition
of tailored distance measures over MTS [9], to the identification and modeling of short-time patterns in the
form of dictionaries or shapelets [10]. In this paper, we focus on recurrent neural networks (RNNs) [7, 11]
for MTS classification. While there are several variations of RNNs, they all share the same fundamental
principle: the MTS is first processed sequentially by a dynamic (possibly adaptable) model, and then the
sequence of its internal states generated over time is exploited to perform classification [12].
Reservoir computing (RC) is a family of RNN models whose recurrent part is kept fixed and is either
generated randomly, or by means of custom topologies for facilitating the information flow [13, 14, 15, 16].
Despite this strong architectural simplification, the recurrent part of the model (the reservoir) provides a
rich pool of dynamic features which are suitable for solving a large variety of tasks. Indeed, RC models can
achieve excellent performance in many fields, including time series forecasting [17, 18, 19, 20, 21], process
modelling [15], and speech analysis [8]. Since the reservoir is fixed, one needs to train only the readout,
∗ Corresponding

author.
Email addresses: filippo.m.bianchi@uit.no (Filippo Maria Bianchi), simone.scardapane@uniroma1.it (Simone
Scardapane), sigurd.lokse@uit.no (Sigurd Løkse), robert.jenssen@uit.no (Robert Jenssen)

1

which provides an instant mapping between the internal representation of the reservoir and the task-specific
output [14, 16]. In machine learning, RC techniques were originally introduced under the name echo state
networks (ESNs) [13]; in this paper, we use the two terms interchangeably.
Several authors have shown that RC models can be a particularly efficient solution for MTS classification [22, 23, 24, 25]. The MTS is fed to the reservoir that generates a sequence of states over time and,
in the simplest case, the last reservoir state becomes a representation of the input, which is subsequently
processed by a classification algorithm for vectorial data [22]. An alternative approach originally proposed
in [26] and later applied to time series classification and fault detection [27, 28, 29], advocates to map the
inputs in a “model-based” feature space where they are represented by statistical properties that potentially
yield better class separability. In this formulation, for each MTS a model is trained to predict its next input
from the current reservoir state, and the parameters of this model become the MTS representation. Since
prediction is performed with ridge regression [26], classification with model space representations is still order
of magnitudes faster than fully adaptable recurrent models trained with gradient descent.
Contributions of the paper
We introduce and evaluate several advanced RC architectures for MTS classification, with a strong focus
on providing training efficiency without hindering the classification performance.
As main contribution, we propose an unsupervised representation of the input MTS that extends the
model space criterion. We call it “reservoir model space” and it consists in the parameters of the one-stepahead predictor that estimates the future reservoir state, as opposed to the future MTS input proposed in
the original formulation [26]. The prediction model in this case must account for all the dynamics in the
reservoir to solve its task and we argue that provides a more accurate characterization of the input MTS,
thanks to the generalization capability of the reservoir. To ensure tractability and computational efficiency
in our approach, we apply dimensionality reduction to the reservoir states, to train the prediction model on
a low-dimensional embedding of the original dynamics. This is accomplished by applying to the reservoir
states sequence a modified version of principal component analysis (PCA) that keeps separated the modes
of variation among time step and data samples.
The proposed representation is at the core of a RC-based classifier, whose reservoir and readout can be
enhanced as follows. First, we extend our preliminary work [30] by considering a bidirectional reservoir to
process in parallel the input sequence both in forward and backward time directions. When implementing
the reservoir model space with a bidirectional reservoir, the prediction model must learn also a memorization
task, as it yields at each time step both the previous and past reservoir state. Second, we consider replacing
the simple linear readout with a support vector machine (SVM) classifier or, still inspired by [30], with a
more sophisticated deep (feedforward) model, i.e., a multilayer perceptron (MLP). While we showed in [30]
that a MLP can improve the classification performance, we only experimented a standard neural network
architecture. In this paper, following recent trends in the deep learning literature we also investigate endowing
the deep readout with flexible nonlinear activation functions, namely Maxout [31] and kernel activation
functions [32].
We evaluate how bidirectional reservoirs and deep readouts affect the performance (both in training time
and classification accuracy) of RC-based classifiers, which use as representations the last reservoir state,
the output model space or the proposed reservoir model space. Experiments are performed on a large set
of experimental benchmarks and a real-world dataset of medical MTS. The RC-based classifiers are also
compared against fully trained RNNs and SVM classifiers configured with precomputed kernels for MTS.
We show that the features extracted using the reservoir model space criterion are fast to compute and highly
informative, as they yield superior classification accuracy.
Structure of the paper
Sec. 2 introduces the problem of classification of MTS with RNN models and the modern architectures
that improve memory capabilities. Sec. 3 is focused on the RC models, introduces the output model
space representation and several advanced techniques for building the classification framework, including
bidirectional reservoirs and deep readouts. Sec. 4 describes the proposed reservoir model space, the adopted

2

dimensionality reduction procedure, and how we extend the bidirectional reservoir formulation in this context.
Sec. 5 performs an extensive experimental evaluation of the techniques, and Sec. 6 provides conclusions.
Notation
We denote variables as lowercase letters (x); constants as uppercase letters (X); vectors as boldface
lowercase letters (x); matrices as boldface uppercase letters (X); tensors as calligraphic letters (X ). All
vectors are assumed to be columns. The operator k·kp is the standard `p norm in Euclidean spaces. The
notation x(t) indicates time step t and x[n] sample n in the dataset.
2. Classification and representation learning with recurrent neural networks
We consider classification of generic F -dimensional MTS observed for T time instants, whose observation
at time t is denoted as x(t) ∈ RF . We represent a MTS in a compact form as a T × F matrix X =
T
[x(1), . . . , x(T )] . The problem of assigning a class label y represented with one-hot encoding to the sequence
X can be framed as a generic density estimation problem:
p (y|x(T ), x(T − 1), . . . , x(1)) .

(1)

A commonly adopted approach in machine learning is to build the sequence model as the combination of
an encoding function and a decoding function. The encoder is used to generate a representation of the
input, while the decoder is a discriminative (or predictive) model that computes the posterior probability of
the output given the representation provided by the encoder. Among all possible choices for the encoding
method, RNNs [33] are a type of artificial neural network particularly suitable to model sequential data. An
RNN is governed by the following state-update equation
h(t) = f (x(t), h(t − 1); θrnn ) ,

(2)

where h(t) is the internal state of the RNN at time t that depends on its previous value h(t − 1) and the
current input x(t), f (·) is a nonlinear activation function (usually a sigmoid or hyperbolic tangent), and θrnn
are the adaptable weights of the RNN. The simplest (vanilla) formulation reads:
h(t) = f (Win x(t) + Wr h(t − 1)) ,

(3)

with θrnn = {Win , Wr }. The matrices Win and Wr are the weights of the input and recurrent connections,
respectively, and are learned with gradient descent in order to generate optimal representations for the task
at hand.
T
From the sequence of the RNN states generated over time, described by the matrix H = [h(1), . . . , h(T )] ,
it is possible to define an encoding (representation) r(H) = rX of the input sequence X. Rather than
accounting for the whole sequence of RNN states, however, it is common to represent the MTS only with a
subset of H. A common choice is to take rX = h(T ), i.e. discard all states except the last one. Thanks to
the ability of capturing temporal dependencies, the RNN can embed into its last state all the information
required to reconstruct the original input [34].
Such a state becomes a fixed-size vectorial representation of the MTS and can be processed by standard
machine learning algorithms. Specifically, the decoder maps the input representation rX into the output
space, which contains all class labels y in a classification task:
y = g(rX ; θdec ).

(4)

where θdec are trainable parameters. In practice, g(·) is either implemented by another neural network or
by a simpler linear model. Its parameters can be jointly learned with the ones contained inside the RNN

3

N

component within a common training procedure. To this end, given a set of MTS {X[n]}n=1 and associated
N
labels {y[n]}n=1 , we can train the model by minimizing an empirical cost:1
∗
∗
θrnn
, θdec
= arg min
θrnn ,θdec


N 


1 X
,
l y[n], g r f (X[n])
N n=1

(5)

where l(·, ·) is a generic loss function (e.g., cross-entropy over the labels). If the encoder has adaptable
weights, their gradient of (5) with respect to θrnn and θdec can be computed by back-propagation through
time [12].
To regularize the parameters values during the training phase, a common procedures is to add an `2
norm penalty to both the weights in the encoding and decoding functions, r(·) and g(·). The penalty
term is controlled by a scalar parameter λ. Furthermore, in all our experiments we also consider dropout
regularization, that randomly drops connection weights with probability pdrop during training [35]. In the
encoding function, we apply a dropout specific for recurrent architectures [36].
2.1. Gated RNN architectures
Despite the theoretical capability of basic RNNs to model any dynamical system, in practice their effectiveness is hampered by the difficulty of training their parameters [37, 38]. To ensure stability, the derivative
of the recurrent function in an RNN must not exceed unity. However, as an undesired effect, the gradient
of the loss shrinks when back-propagated in time through the network. Using RC models (described in
the next section) is one way of avoiding this problem. Another common solution is the long short-term
memory (LSTM) network [39]; differently from the vanilla RNNs in (3), LSTM exploits gating mechanisms
to maintain its internal memory unaltered for long time intervals. However, LSTM flexibility comes at the
cost of a higher computational and architectural complexity. A popular variant is the gated recurrent unit
(GRU) [40], that provides a better memory conservation by using less parameters than LSTM, but its state
update requires an additional operation, hence a higher computational cost. Both LSTM and GRU still
require to back-propagate through time the gradient of their loss.
3. Reservoir computing classifiers
In order to avoid the costly operation of back-propagating through time, the standard ESN architecture
takes a radical different direction: it still implements the encoding function in (3), but with the fundamental
difference that the weights parameters Win and Wr are randomly generated and left untrained (or, possibly,
implemented according to a prefixed topology [41]). To compensate this lack of adaptability, when processing
the input sequence the recurrent layer (reservoir) generates a rich pool of heterogeneous dynamics, from which
the ones useful to solve many different tasks can be drawn. The generalization capabilities of the reservoir
mainly depend on three ingredients: (i) a large number of processing units in the recurrent layer, (ii) sparsity
of the recurrent connections, and (iii) a spectral radius of the connection weights matrix Wr , set to bring
the system to the edge of stability [42]. The behaviour of the reservoir can therefore be controlled by simply
modifying the following structural hyperparameters instead of training the internal weight matrices: the
spectral radius ρ; the percentage of non-zero connections β; and the number of hidden units R. Another
important hyperparameter is the scaling ω of the values in Win , which controls the amount of nonlinearity
in the processing units and, jointly with ρ, can shift the internal dynamics from a chaotic to a contractive
regime. Finally, a Gaussian noise with standard deviation ξ can be added to the argument of the state
update function (3) for regularization purposes [13].
In ESNs, the decoder (commonly referred as readout) is usually a linear model:
y = g(rX ) = Vo rX + vo

(6)

1 Note that MTS may have different lengths. For readability, in the paper we refer with T to the length of a single MTS,
implicitly assuming that T is a function of the MTS itself.

4

2

The weight matrix Vo and biases vo can be learned by means of ridge regression by setting l(a, b) = ka − bk2
in (5), which provides a closed-form solution [16].
The combination of a (untrained) reservoir and a linear readout defines the basic ESN model [13]. In the
remainder of this section, we present multiple extensions (e.g., more sophisticated reservoirs and readouts),
which have been proposed in the literature over the last years to improve the basic architecture. Most of
these extensions can be considered independently one from the other: a complete overview of all the possible
combinations considered in this paper is given later on in Fig. 2 at the end of Sec. 4.
3.1. Output model space
As stated in the introduction, rather than representing the MTS by simply taking the last reservoir state
h(T ), [26] proposed to process each MTS with a common reservoir and then fit an additional linear model
(one for each time series), whose parameters become the representation of the MTS. In particular, a ridge
regression model is trained to implement an output function that performs one step-ahead prediction of each
input MTS:
x(t + 1) = Uo h(t) + uo
(7)
The parameters in Uo and uo completely characterize the linear model and θo = [vec(Uo ); uo ] ∈ RF (R+1)
becomes the representation rX of the MTS, which is, in turn, processed by the classifier in (4).
It is possible to see that rX characterize the following discriminative model
p (x(t + 1)|h(t); rX ) .

(8)

At the same time, the reservoir encoding that operates according to (3) can be expressed as
p (h(t)|h(t − 1), x(t)) .

(9)

The state h(t − 1) in (9) depends, in turn, from x(t − 1) and h(t − 2), and so on. By expressing those
dependencies explicitly with the chain rule of probability and then plugging (9) in (8) one obtains
T
Y

p (x(t + 1)|x(t), x(t − 1), . . . , x(1); rX ) = p (x(T ), x(T − 1), . . . , x(1); rX ) ,

(10)

t=1

Eq. (10) is as a generative model [43], which provides a powerful characterization of the inputs and
also induces a metric relationship between samples [44]. As a result, classification with the model space
criterion can be cathegorized as a hybrid discriminative/generative method, and the representation rX (learnt
unsupervisedly) can be exploited in a variety of different tasks. Indeed, the model view characterization of
a MTS has proven effective in anomaly detection [27], classification [28, 29], and to build a kernel similarity
matrix [26].
3.2. Bidirectional reservoir
Bidirectional architectures have been successfully applied in RNNs to extract from the input sequences
temporal features that account also for dependencies very far in time [7]. Within the RC framework, a
bidirectional reservoir has been used in the context of time series prediction to collect future information,
which is provided at training stage to improve the accuracy of the model [15]. The future information is not
available during test and the model uses a standard reservoir to extract temporal features. On the other
hand, in a classification setting the whole time series is given at once, in both training and test stages. In
[30] the last state of a bidirectional reservoir is used as a fixed-size representation for classifying the input
MTS.
Bidirectionality is implemented by feeding into the same reservoir an input sequence both in straight
and reverse order, so that the following states are generated


~h(t) = f Win x(t) + Wr ~h(t − 1) ,


(11)
~ = f Win x(t)
~ − 1) ,
~ + Wr h(t
h(t)
5

~ = x(T − t). The final states sequence of the bidirectional
where x(t)
reservoir
is obtained by concatenating
i
h
b
~
~
at each time interval the backward and forward states h (t) = h(t); h(t) ∈ R2R .
The representation of the input as the last state generated by the reservoir in this case would become
rX = hb (T ). The main advantage of using a bidirectional reservoir in conjunction with the last state
representation is that, compared to unidirectional architectures, temporal dependencies spanning longer time
intervals can be captured. Since the reservoir trades its internal stability with a vanishing memory of the
past inputs [42], at time T the state ~h(T ) maintains scarce information about the first inputs processed. On
~ ) is more influenced by the first inputs and, therefore, hb (T ) summarizes well recent and
the other hand, h(T
past information. This representation improves classification results especially when important information
are contained also in the first part of the input.
3.3. Dimensionality reduction of the states
Dimensionality reduction provides a mapping from a high dimensional input space to a new lower dimensional one, usually relying on unsupervised criteria. Since the reservoir is characterized by a large
number of neurons, dimensionality reduction applied on top of h(t) yields a more compact representation,
which can provide a regularization to the model that enhances its generalization capability and simplifies
its training [45]. The role of dimensionality reduction becomes particularly important in the presence of a
bidirectional reservoir, as the state dimension is doubled.
In the context of RC, commonly used algorithms for reducing the dimensionality of the reservoir are PCA
and kernel PCA, which project data on the first D eigenvectors of a covariance matrix [46]. For example,
in [47] we dealt with a prediction task and we applied dimensionality reduction to a single sequence of
reservoir states generated by the input MTS. On the other hand, in [30] we tackled a classification task with
multiple MTS, each one associated to a different sequence of states. Since we used as representation the last
reservoir states, those were stacked into a matrix to which dimensionality reduction was applied.
3.4. Nonlinear readouts
In a standard ESN, the output layer is a linear readout that is quickly trained by solving a convex
optimization problem. However, a simple linear model might not possess sufficient representational power
for modeling the high-level embeddings derived from the reservoir states. For this reason, several authors
proposed to replace the standard linear decoding function g(·) in (6) with a nonlinear model, such as support
vector machines (SVMs) [18, 20] or MLPs [48, 49, 50].
In particular, MLP is an universal function approximator that can learn complex representations of
the input by stacking multiple layers of neurons configured with non-linear activations, e.g., rectified linear
units (ReLUs). Deep MLPs are known for their capability of disentangling factors of variations from highdimensional spaces [51], and therefore can be more powerful and expressive in their instantaneous mappings
from the representation to the output space than linear readouts.
Readouts implemented as MLPs accomplished only modest results in the earliest works on RC [14].
However, nowadays MLPs can be trained much more efficiently by means of sophisticated initialization
procedures [52], regularization techniques [35] and more expressive activation functions [31, 32]. In [30] we
showed that combining ESNs with MLPs trained with modern techniques can provide a substantial gain
in performance as compared to a linear formulation. Also in this paper, we consider nonlinear readouts
performed by means of a MLP, whose number of layers determines a “feedforward” depth to the RNN [11].
4. Reservoir model space
In this section we introduce the main contribution of this paper, i.e., the reservoir model space for
representing MTS. To make the model tractable, we extend the approach from Sec. 3.3 for reducing the
dimensionality of reservoir features, by computing PCA on data represented as matrices rather than vectors.
In addition, we discuss the effect of a bidirectional reservoir (see Sec. 3.2) on the proposed representation.

6

4.1. Proposed MTS representation
The great generalization capability of the reservoir is grounded on the large amount of different dynamical
features it generates from the input time series. Indeed, the readout selects and combines the ones which are
useful to accomplish a specific task. In the case of prediction, different combinations of these dynamics are
accounted depending on the forecast horizon of interest. Therefore, by fixing a specific prediction step (e.g.,
1 step-ahead) we implicitly ignore all those dynamics that are not particularly useful to solve the task. We
argue that this potentially introduces a bias in the model space induced by output prediction, since some
dynamical features that are not important for the prediction task can still be useful to characterize the MTS.
Therefore, we propose a variation of the output model space (see Sec. 3.1), where each MTS is represented
as the parameters of a linear model that predicts the next reservoir state, rather than the output sample.
The prediction model in this case must account for all the dynamics in the reservoir to solve its task and we
argue that provides a more accurate characterization of the input MTS. More formally, the following linear
model is trained to implement the prediction of the next reservoir state
h(t + 1) = Uh h(t) + uh ,

(12)

and, thus, rX = θh = [vec(Uh ); uh ] ∈ RR(R+1) is our proposed representation of the input MTS.
Contrarily to the output model space representation presented in Sec. 3.1, the prediction model here
describes a generative model of the reservoir sequence, rather than of the input sequence. The learned feature
space can capture both the data and their generative process. As for the output model space, a classifier
that processes the reservoir model representation combines the explanatory capability of generative models
with the classification power of the discriminative methods.
Due to the potentially large size of the reservoir, the parameters of the prediction model in (12) would
be too much to make this idea tractable in a naı̈ve fashion. Besides the undesired effects of producing a
representation of the input in a large dimensional space without enforcing sparsity constraints, evaluating
the ridge regression solution for each MTS would require too many computational resources, which is in
disagreement with our purpose of designing an efficient classification framework. As a solution, in the next
section we propose to reduce the number of features in the reservoir, by keeping separated the modes of
variation among time in the different data samples.
4.2. Dimensionality reduction for reservoir states tensor
To ensure efficiency, we fit the prediction model on an embedding of the reservoir states sequence given
by a dimensionality reduction procedure, which extends the methods previously discussed. In particular, in
both cases described in Sec. 3.3, dimensionality reduction was performed by applying regular PCA or kPCA
to a dataset represented as a matrix (tensor of order two). However, our proposed MTS representation
no longer coincides with the last reservoir state, but derives from the whole sequence of states generated
over time. Therefore, we conveniently describe our dataset as a 3-mode tensor H ∈ RN ×T ×R and require a
procedure to map R → D s.t. D  R, while maintaining the other dimensions unaltered. We note that for
MTS of varying length, zero-padding (or a more elaborate interpolation procedure) is required to build the
tensor.
Dimensionality reduction on high-order tensors can be achieved through Tucker decomposition [53], which
decomposes a tensor into a core tensor (the lower-dimensional representation) multiplied by a matrix along
each mode. When only one dimension of H is modified, Tucker decomposition becomes equivalent to applying
a two-dimensional PCA on a particular matricization of H [54]. In particular, to reduce the third dimension
(R) one computes the mode-3 matricization of H by arranging the mode-3 fibers (high-order analogue of
matrix rows/columns) to be the rows of a resulting matrix H(3) ∈ RN T ×R . Then, standard PCA projects the
rows of H(3) on the eigenvectors associated to the D largest eigenvalues of the covariance matrix C ∈ RR×R ,
defined as
C=

NT
X

T
1
hi − h̄ hi − h̄ .
N T − 1 i=1

7

(13)

3/7/2018

rmESN.xml
^
H

H

N

h(t)

...
R

dimensionality
reduction

h(t + 1)

model space
embedding

...

D
θh [N ]

U h [1], uh [1]

θh [1]

T

Figure 1: Schematic depiction of the procedure to generate the reservoir model space representation. For
each input MTS X[n] a sequence of states H[n] is generated by a fixed reservoir. Those are the frontal
slices (dimension N ) of H, but notice that in the figure lateral slices (dimension T ) are shown. By means
of dimensionality reduction, the reservoir features from R become D. A distinct linear model is trained to
predict the columns of each frontal slice Ĥ of Ĥ. The parameters θh [n] of the model trained on Ĥ[n] become
the representation of X[n]
PN T
In (13), hi is the i-th row of H(3) and h̄ = N1 i hi . As a result of the concatenation of the first
two dimensions in H, C evaluates the variation of the components in the reservoir states across all samples
and time steps at the same time. Consequently, both the original structure of the dataset and the temporal
orderings are lost, as reservoir states relative to different samples and generated in different time steps are
mixed together. This may lead to a potential loss in model accuracy, as the existence of modes of variation
in time courses within individual samples is ignored [55]. To address this issue, we consider as individual
samples the matrices Hn ∈ RT ×H , obtained by slicing H across its first dimension. The sample covariance
matrix in this case reads
N
T

1 X
S=
Hn − H̄
Hn − H̄ .
(14)
N − 1 n=1
The first D leading eigenvectors of S are stacked in a matrix E ∈ RR×D and the desired tensor of reduced
dimensionality is obtained as Ĥ = H ×3 E, where, ×3 denotes the 3-mode product.
Like C, S ∈ RR×R describes the variations of the variables in the reservoir. However, the whole time
sequence of reservoir states generated by a MTS is considered as an observation. Accordingly, states pertaining to different MTS are grouped together and their temporal ordering is preserved. We notice that an
analogous approach was adopted to reduce the dimensionality of images [56], which are characterized by
multiple spatial dimensions, while the dimensions in our samples represent space and time, respectively.
After dimensionality reduction, the model in (7) becomes
ĥ(t + 1) = Uh ĥ(t) + uh ,

(15)

where ĥ(·) are the columns of a frontal slice Ĥ of Ĥ, Uh ∈ RD×D , and bh ∈ RD . The representation will now
coincide with the following vector of parameters rX = θh = [vec(Uh ); uh ] ∈ RD(D+1) , whose dimensionality
is controlled by the described dimensionality reduction procedure. A schematic description of the proposed
unsupervised procedure to derive the reservoir model representation is shown in Fig. 1
4.3. Reservoir model space with bidirectional reservoirs
In the reservoir model space formulation a bidirectional reservoir modifies the representation rX , as the
linear model in (12) becomes
h
i
h
i
~ + 1) = Ub ~h(t); h(t)
~
h(t + 1); h(t
+ ubh ,
(16)
h
8

3/7/2018

main_schema.xml

dimensionality
reduction

...

MLP

last state

unidirectional
reservoir

...

h(T )

...

output model
lin reg
h(t)

X

¯
H

H

bidirectional
reservoir

rX

...

...

x(t + 1)

t

SVM
reservoir model

...
t

y

θo

...
h(t)

h(t + 1)

θh

Figure 2: Proposed variants of the RC-based classifiers. Several distinct models are obtained by selecting
different options for each processing block.
where Ubh ∈ R2R×2R and ubh ∈ R2R are the new set of parameters. In this case, the linear model is
trained to solve two different tasks at the same time, which are predicting the next state h(t + 1) and
~ + 1) (or equivalently their low-dimensional embeddings). We argue that
reproducing the previous one h(t
such a model provides a more accurate representation of the input, by modeling at the same time the
temporal dependencies in both time directions. We notice that the bidirectional reservoir produces a similar
effect also in the output model space (see Sec. 3.1), where the linear model learns to reproduce the previous
input and to predict the next one.
When using the bidirectional reservoir, the state size is doubled and, accordingly, also the size of the
representation increases. In particular, when the last state is used as representation its size is doubled, rX ∈
R2R . In the output and reservoir space representations instead, we have rX ∈ RF (2R+1) and rX ∈ R2R(2R+1) ,
respectively. Also in this case, by applying dimensionality reduction as described in the previous section, we
limit the computational complexity and reduce the risk of overfitting
4.4. Summary of the RC classifier variants
The presented variants of the RC-based classifiers are schematized in Fig. 2. The input MTS X is
processed by a reservoir, which is either unidirectional or bidirectional, and it generates the sequence of
states H over time. An optional dimensionality reduction step can be applied to reduce the number of
reservoir features, and a new sequence H̄ is obtained. Three different approaches can be chosen to generate
a representation rX of the input from the sequence of reservoir features: the last state of the sequence h(T ),
the output state model θo , or the proposed reservoir state model θh . The obtained representation rX is finally
processed by a decoder (readout), which is implemented by a classifier for real-valued vectors. As previously
discussed, we focus on three different classifiers, which are the standard linear readout implemented by ridge
regression, a SVM classifier and a deep neural network (MLP). The classifier outputs the class label y to be
assigned to the MTS X.
The implementation of the RC classifier, with all the presented variants, is publicly available online.2
2 https://github.com/FilippoMB/Reservoir-model-space-classifier

9

5. Experiments
In this section we test the different RC-based architectures for MTS classification presented in the paper,
including our proposed classifiers with the reservoir model space embedding. We also compare against the
performance obtained by RNNs classifiers trained with gradient descent (LSTM and GRU) and with SVM
classifiers using pre-trained kernels evaluated by means of similarities for MTS. Depending if the input MTS
in the RC-based model is represented by the last reservoir state (rX = h(T )), or by the output space model
(Sec. 3.1), or by the reservoir space model (Sec. 4), we refer to the models as lESN, omESN and rmESN,
respectively. Whenever we use a bidirectional reservoir, a deep readout or a SVM readout we add the prefix
“bi-”, “dr-”, and “svm-”, respectively (e.g., bi-lESN or dr-bi-rmESN ).
First, we introduce the MTS datasets for classification that we consider and we define our experimental
setup. In Sec. 5.1, we compare the performance obtained on several benchmark datasets, by the different
representations yielded by the RC architectures and by the fully trained RNNs. In Sec. 5.2, we analyze
how performance change in the RC-based models, when the representations are first generated using a
bidirectional reservoir and then are processed with MLP readouts. Finally, in Sec. 5.3 we process a realworld dataset and we further compare the classification performance with two time series kernels, using a
SVM classifier.
Benchmark datasets. To provide an extensive evaluation of the performance of each classification model,
we consider ten benchmark classification datasets of univariate and multivariate time series taken from the
UCR and UCI repositories,3 whose details are reported in Tab. 1.
Table 1: Time series benchmark datasets details. Column 2 to 5 report the number of variables (#V ),
samples in training and test set, and number of classes (#C), respectively. Tmin is the length of the shortest
MTS in the dataset and Tmax the longest MTS.
Dataset

#V Train

Test

#C

T min T max Source

Swedish Leaf
Chlorine Concentration
DistPhal
ECG
Libras
Ch.Traj.
Wafer
Jp.Vow.
Arab. Dig.
Auslan

1
1
1
2
2
3
6
12
13
22

625
3840
139
100
180
2558
896
370
2200
1425

15
3
3
2
15
20
2
9
10
95

128
166
80
39
45
109
104
7
4
45

500
467
400
100
180
300
298
270
6600
1140

128
166
80
152
45
205
198
29
93
136

UCR
–
UCR
UCR
UCI
UCI
UCR
UCI
UCI
UCI

Blood samples dataset. As real-world case study, we analyze MTS of blood measurements obtained from
electronic health records of patients undergoing a gastrointestinal surgery at the University Hospital of North
Norway in 2004–2012.4 Each patient is represented by a MTS of 10 blood sample measurements collected for
20 days after surgery. We consider the problem of classifying patients with and without surgical site infections
from their blood samples, collected 20 days after surgery. The dataset consists of 883 MTS, of which 232
pertain to infected patients. The original MTS contain missing data, corresponding to measurements not
collected for a given patient at certain time intervals, which are replaced by zero-imputation in a preprocessing
step.
3 www.cs.ucr.edu/

~eamonn/time_series_data, archive.ics.uci.edu/ml/datasets.html
dataset has been published in the AMIA Data Competition 2016: https://groups.google.com/forum/#!topic/
ml-news/MQtVkxizrrU
4 The

10

0.34

88
0.32

87
86

0.3

85
84
25

50

75

100

125

Training time

Classification accuracy

89

0.28
150

D

Figure 3: Classification accuracy and training time for different dimensions D of the space with reduced
dimensionality.
Experimental setup. For each dataset, we train the models 10 times using independent random parameters
initializations and each model is configured with the same hyperparameters in all experiments. Indeed, the
reservoir is usually very sensitive to hyperparameter setting and, therefore, a proper configuration in RC
models is usually much more critical than in networks trained with gradient descent, such as LSTM and
GRU. Nevertheless, we show that by using “default” hyperparameters values (set according to common rules
of thumb), rather than tuning them with an independent validation procedure for each task, the RC-based
approaches, especially the proposed rmESN, is robust and achieves superior performance. In a later section
we also consider fine-tuning for all architectures.
To provide a more significant comparison, lESN, omESN and rmESN always share the same randomly
generated reservoir, which is configured with the following hyperparameters: internal units R = 800; spectral
radius ρ = 0.99; non-zero connections percentage β = 0.25; input scaling ω = 0.15; noise level ξ = 0.01. When
classification is performed with a ridge regression readout, we set the regularization parameter λ = 1.0. The
ridge regression prediction models, used to generate the model-space representation in omESN and rmESN,
are instead configured with λ = 5.0.
In each RC-based method, we apply the dimensionality reduction procedure proposed in Sec. 4, as
it provides important computational advantages (in terms of both memory and CPU time), as well as a
regularization that improves the overall generalization capability and robustness of the models. To determine
the optimal value D to project H into H̄, we evaluate how the training time and the average classification
accuracy (computed with a k-fold cross-validation procedure) of the RC classifiers varies on the benchmark
dataset in Tab. 1. We report the average results in Fig. 3. We notice that while the training time increases
approximately linearly with D, it is possible to identify an “elbow” in the classification accuracy for D = 75.
We, therefore, select such value as it provides a compromise between high accuracy and low training time.
Finally, LSTM and GRU are configured with H = 30 hidden units; the decoding function is implemented
as a neural network with a dense layer of 20 hidden units followed by a softmax layer; the dropout probability
is pdrop = 0.1; the `2 regularization parameter is λ = 0.0001; gradient descent is implemented with the Adam
algorithm [57] and we train the models for 5000 epochs.
5.1. Comparison of baseline architectures
In this experiment we compare the classification accuracy obtained on the representations yielded by
the RC models, lESN, omESN and rmESN, and by the fully trainable RNNs, implementing either GRU or
LSTM cells. Results are evaluated on the benchmark datasets in Tab. 1. In this experiment, the decoder
is implemented by linear regression in the RC models and by a non-linear function in LSTM and GRU, as
described at the beginning of this section. Since all the other parameters in LSTM and GRU are learned with
a non linear optimization technique, the non-linearities in the decoding function do not result in additional
computational costs. Results are reported in Fig. 4. The first panel reports the mean classification accuracy
and standard deviation obtained by 10 independent runs on all benchmark datasets. rmESN achieves the

11











5DQNLQJ

7UDLQLQJWLPHPLQV

&ODVVLILFDWLRQDFFXUDF\











*58

/670

UP(61

RP(61

(61

*58

/670

UP(61

RP(61

(61

*58

/670

UP(61

RP(61

(61



Figure 4: Comparison of the average results obtained on all benchmark datasets.
best results, both in terms of high mean accuracy and low standard deviation, and is followed by omESN.
The other methods yield higher standard deviation and lESN obtains the worst average accuracy overall.
The second panel shows the average training time of each model. As expected, LSTM and GRU are much
slower (more than 10 minutes versus few seconds for the RC models), due to the time consuming gradient
descent optimization procedures. On the other hand, the RC architectures exploit linear models that can
be trained efficiently due to their closed form solution. Interestingly, despite the high number of parameters
in rmESN, mainly due to the large size of the Uh matrix, its training time is still low and is comparable to
the other, less complex, RC architectures.
Finally, the third panel shows a ranking score computed as follows: a model is assigned with 5 points
when it achieves the best results in a benchmark problem, 4 points if it reaches the second best result, and
so on. The sum of these points obtained across each task is the final ranking reported. As we can see, the
gradient-based RNNs achieve better results than lESN and omESN according to this metric than to the
average classification accuracy; the reason is that on few datasets LSTM and GRU perform poorly, with a
consequent drop on the average accuracy. Nonetheless, the proposed rmESN yields superior representations,
obtaining the highest ranking score even when processed by a simple linear readout. The result outlined by
the ranking score is particularly interesting because, while LSTM and GRU exploit supervised information to
learn the representations rX , in rmESN the encoding procedure is completely unsupervised. Furthermore,
the discriminative classifier implemented by LSTM and GRU is non-linear, and thus it provides a more
powerful modeling capability.
5.2. Experiments with bidirectional reservoir and deep-readout
In this experiment we investigate the effects of a bidirectional reservoir and a deep-readout, implemented
by a MLP, in terms of classification accuracy and training time of the RC-based classifiers. To further increase
the flexibility of the readout, apart from using standard rectified linear units in the MLP, we also employ more
sophisticated transfer functions, namely Maxout [31] and kernel activation functions (KAFs) [32]. Thanks
to their adaptable parameters, trained jointly with the other MLP weights, these functions can learn more
complicated relationships and improve the expressive capability of the model at every layer. We refer the
reader to the original publications for details on their formulation. The deep readout is implemented with
3 layers of 20 neurons each and is trained for 5000 epochs, using a dropout probability pdrop = 0.1 and `2
regularization parameter λ = 0.001.
Once again, we evaluate the models on the all the benchmark datasets and in Fig. 5 we report results
in terms of classification accuracy and training time. We can see that both the bidirectional reservoir
and deep readout improve, to different extents, the classification accuracy of each RC classifier. The largest
improvement occurs in lESN when implemented with a bidirectional reservoir. Indeed, the representation rX
provided by this model is the last state, which depends mostly on the last observed values of the input MTS.
12

(61
EL(61
GU(615H/8
GU(610D[
GU(61.$)
RP(61
ELRP(61
GURP(615H/8
GURP(610D[
GURP(61.$)
UP(61
ELUP(61
GUUP(615H/8
GUUP(610D[
GUUP(61.$)

7UDLQLQJWLPHPLQV

13

GUUP(61.$)

GUUP(610D[

GUUP(615H/8

ELUP(61

UP(61

GURP(61.$)

GURP(610D[

GURP(615H/8

ELRP(61

RP(61

GU(61.$)

GU(610D[

GU(615H/8

EL(61

(61

&ODVVLILFDWLRQDFFXUDF\






















Figure 5: Classification accuracy, ranking and training time when using rmESN with a bidirectional reservoir
and deep readouts, configured with ReLUs, KAFs, and Maxout activations.

In all those cases where the most relevant information is contained at the beginning of the input sequence
or when the MTS are too long and the short-term memory of the reservoir prevents to capture long-term
dependencies, the bidirectional architecture greatly improves the lESN representation. The bidirectional
reservoir slightly improves the performance also in omESN and rmESN. We recall that in these cases, rather
than learning only a model for predicting the next output/state, when using a bidirectional reservoir we learn
a model that also performs a memorization task. Therefore, its parameters provide a further characterization
of the input. However, the improvement of the classification performance in these model is limited, probably
because the representations obtained with a unidirectional reservoir characterize the input sufficiently well.
A deep-readout also enhance the accuracy in each model; improvements are larger in lESN and more
limited in omESN and rmESN. Once again, this underlines that the weaker lESN representation benefits
by adding more complex mechanisms in the classification framework. Even more than the bidirectional
reservoir, a deep-readout trades better modeling capabilities with more computational resources, especially
when implemented with adaptive activation functions. Remarkably, when using Maxout functions rather
than standard ReLUs, the training time is slightly higer, but we obtain significant improvements in the
average classification accuracy. In particular, dr-omESN (Maxout) obtains almost the same performance of
the basic version of rmESN, while dr-rmESN (Maxout) reaches a mean classification accuracy beyond 90%.
Another interesting result obtained by both Maxout and KAF is a reduction in the standard deviation of
the accuracy, meaning that they implement a more robust classifier.
As a final consideration, we observe that the improvements provided by the bidirectional reservoirs and
the deep-readouts to rmESN are smaller than for lESN and omESN. This underlines how the rmESN
representations are already very powerful in their basic formulation.
5.3. Classification of blood samples MTS
As last experiment, we evaluate the performance of the RC classifiers on the blood sample MTS. In
this test we consider only omESN and rmESN, as they provide an optimal compromise between training
efficiency and classification accuracy, and are the main focus in this work. For this test, we include in the
comparison also two state-of-the-art kernels for MTS. The first is the learned pattern similarity (LPS) [58],
which identifies segments-occurrence within the MTS by means of regression trees. Those are used to generate
a bag-of-words type compressed representation, on which the similarity scores are computed. The second
method is the time series cluster kernel (TCK) [9], which is based on an ensemble learning approach wherein
the clustering results of many Gaussian mixture models, fit several times on random subsets of the original
dataset, are joined to form the final kernel.
For LPS and TCK, SVM is configured with the pre-computed kernels returned by the two procedures,
while for omESN and rmESN we build a RBF kernel with bandwidth γ. Note that this corresponds to
implementing the decoding function (4) with a non-linear SVM, as discussed in Sec. 3.4. We optimize on
a validation set the SVM hyperparameters, which are the smoothness of the decision hyperplane, c, and
bandwidth, γ (only omESN and rmESN ). The hyperparameter space is explored with a grid search, by
varying c in [0.1, 5.0] with resolution 0.1 and γ in [0.01, 1.0] with resolution 0.01. LPS is configured using 200
regression trees and maximum segments length 10. TCK is configured with 40 different random initializations
and 30 maximum mixtures for each partition. RC classifiers use the same hyperparameters as in the previous
experiments.
To compute the performance of the models, those are evaluated 15 times with independent random
initializations and randomly splitting the original dataset into training, validation, and test set, containing
70%, 10% and 20% of the original samples respectively, mixed in a random order. Each time, we normalize
the data by subtracting the mean and dividing by the standard deviation of each variable in the training
set. Mean and standard deviations are computed excluding the imputed values.
The results are depicted in Fig. 6. For completeness, we report also the classification results obtained
on this task by omESN and rmESN, with g(·) implemented as a linear readout. Also in this case, rmESN
outperforms omESN either when it is configured with a linear or a SVM readout. As for the deep-readout,
we notice that the more powerful decoding function improves the classification accuracy in rmESN only
slightly, while the increment in omESN is much larger. Nevertheless, svm-rmESN manages to outperform

14

7UDLQLQJWLPHPLQV

&ODVVLILFDWLRQDFFXUDF\









 

VYPUP(61
VYPRP(61
VYP/36
VYP7&.
UP(61
RP(61

RP(61

UP(61

VYP7&.

VYP/36

VYPRP(61

VYPUP(61



Figure 6: Classification accuracy obtained with SVM using different precomputed kernels. We also report
the results obtained by rmESN and omESN on the same problem.
the SVM classifiers configured with LPS and TCK kernels. We notice standard deviations in all methods are
quite high, since the train-validation-test splits are generated randomly at every iteration and, therefore, the
classification task changes each time. TCK yields results with the lowest standard deviation and is followed
by the two versions of rmESN. Using a SVM readout increases the training time of the RC models, especially
rmESN, but is still much lower than computing the TCK and LPS kernels.
6. Conclusions and future work
In this work we investigated several alternatives to build a classifier based on reservoir computing, focusing
on unsupervised procedures to learn fixed-size representations of the input time series. As main contribution,
we proposed a classifier based on the ESN reservoir model space representation, which can be categorized
as a hybrid generative-discriminative approach. Specifically, the parameters of a model that predict the
next reservoir states characterize the generative process of high-level dynamical features of the inputs. Such
parameters are, in turn, processed by a discriminative decoder that classifies the original time series.
In a hybrid generative-discriminative approach where data are assumed to be generated by a parametric
distribution, the subsequent discriminative model cannot be specified independently from the generative
model type, otherwise classification accuracy could be biased [59]. However, in our case the reservoir is
flexible and generic, as it can express a large variety of dynamical features of the input. Therefore, the
reservoir captures relationships among the data, without posing constraints on the particular model underlying the data distribution. This provides two advantages: (i) different discriminative models can be used in
conjunction with the same reservoir model space representation and (ii) the same reservoir can be used to
model data generated by different distributions.
To make the reservoir model space tractable we applied an unsupervised dimensionality reduction procedure, suitable for datasets represented as high-order tensors. Such dimensionality reduction can greatly
reduce computational time and memory usage, without hampering the final classification accuracy and can
also be applied to generate other ESN-based representations.
We evaluate the proposed model on several benchmarks problems for classification of multivariate time
series, showing that it can outperform other ESN models and fully trained RNNs with a lower amount of
fine-tuning and a much smaller training time. We also considered a real-world case study of time series
pertaining to blood samples data and we compared our method with kernels for multivariate time series.
Results show that our approach yields a slightly superior accuracy with a significantly lower training time.
Interestingly, combining the proposed representation with more sophisticated architectures for the reservoir
15

and the readout provides only a minimal contribution to the accuracy, pointing to the strong informative
content of this representation in terms of discriminative power.
In future works, we are interested in extending this model to handle additional architectures for the
reservoir (e.g., fixed topologies [41] or leaky neurons [60]), and to extend the reservoir model to manage
missing data in the time series leveraging its generative properties.
Acknowledgments
This work was partially funded by the Norwegian Research Council FRIPRO grant no. 239844 on
developing the Next Generation Learning Machines.
The authors would like to thank Arthur Revhaug, Rolv-Ole Lindsetmo and Knut Magne Augestad, all
clinicians currently or formerly affiliated with the gastrointestinal surgery department at the University
Hospital of North Norway, for preparing the blood samples dataset. The authors would also like to acknowledge Karl Øyvind Mikalsen, UiT, and Cristina Soguero-Ruiz, University of Rey Juan Carlos, for discussions
regarding this data set.
References
References
[1] A. Kampouraki, G. Manis, C. Nikou, Heartbeat time series classification with support vector machines,
IEEE Transactions on Information Technology in Biomedicine 13 (4) (2009) 512–518.
[2] P. Buteneers, D. Verstraeten, B. Van Nieuwenhuyse, D. Stroobandt, R. Raedt, K. Vonck, P. Boon,
B. Schrauwen, Real-time detection of epileptic seizures in animal models using reservoir computing,
Epilepsy Research 103 (2-3) (2013) 124–134.
[3] G. Clifford, C. Liu, B. Moody, L. Lehman, I. Silva, Q. Li, A. E. Johnson, R. G. Mark, AF classification from a short single lead ECG recording: The Physionet computing in cardiology challenge 2017,
Computing in Cardiology 44 (2017) 1–4.
[4] K. Ø. Mikalsen, F. M. Bianchi, C. Soguero-Ruiz, S. O. Skrøvseth, R.-O. Lindsetmo, A. Revhaug,
R. Jenssen, Learning similarities between irregularly sampled short multivariate time series from EHRs,
in: Proc. 3rd International Workshop on Pattern Recognition for Healthcare Analytics at ICPR 2016,
2016.
[5] E. Carden, J. Brownjohn, Arma modelled time-series classification for structural health monitoring of
civil infrastructure, Mechanical Systems and Signal Processing 22 (2) (2008) 295–314.
[6] D. Hunt, D. Parry, Using echo state networks to classify unscripted, real-world punctual activity, in:
Engineering Applications of Neural Networks, Springer, 2015, pp. 369–378.
[7] A. Graves, J. Schmidhuber, Framewise phoneme classification with bidirectional LSTM and other neural
network architectures, Neural Networks 18 (5-6) (2005) 602–610.
[8] E. Trentin, S. Scherer, F. Schwenker, Emotion recognition from speech signals via a probabilistic echostate network, Pattern Recognition Letters 66 (2015) 4–12.
[9] K. Ø. Mikalsen, F. M. Bianchi, C. Soguero-Ruiz, R. Jenssen, Time series cluster kernel for learning
similarities between multivariate time series with missing data, Pattern Recognition 76 (2018) 569 –
581.
[10] A. Bagnall, J. Lines, A. Bostrom, J. Large, E. Keogh, The great time series classification bake off:
a review and experimental evaluation of recent algorithmic advances, Data Mining and Knowledge
Discovery 31 (3) (2017) 606–660.
16

[11] R. Pascanu, C. Gulcehre, K. Cho, Y. Bengio, How to construct deep recurrent neural networks, in:
Proc. 2017 International Conference on Learning Representations (ICLR), 2014.
[12] A. Graves, A.-r. Mohamed, G. Hinton, Speech recognition with deep recurrent neural networks, in:
Proc. 2013 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP),,
IEEE, 2013, pp. 6645–6649.
[13] H. Jaeger, The “echo state” approach to analysing and training recurrent neural networks-with an
erratum note, GMD Technical Report 148 (34).
[14] M. Lukoševičius, H. Jaeger, Reservoir computing approaches to recurrent neural network training,
Computer Science Review 3 (3) (2009) 127–149.
[15] A. Rodan, A. Sheta, H. Faris, Bidirectional reservoir networks trained using SVM+ privileged information for manufacturing process modeling, Soft Computing 21 (22) (2017) 6811–6824.
[16] S. Scardapane, D. Wang, Randomness in neural networks: an overview, Wiley Interdisciplinary Reviews:
Data Mining and Knowledge Discovery 7 (2).
[17] Z. Shi, M. Han, Support vector echo-state machine for chaotic time-series prediction, IEEE Transactions
on Neural Networks 18 (2) (2007) 359–372.
[18] D. Li, M. Han, J. Wang, Chaotic time series prediction based on a novel robust echo state network,
IEEE Transactions on Neural Networks and Learning Systems 23 (5) (2012) 787–799.
[19] A. Deihimi, H. Showkati, Application of echo state networks in short-term electric load forecasting,
Energy 39 (1) (2012) 327–340.
[20] F. Bianchi, S. Scardapane, A. Uncini, A. Rizzi, A. Sadeghian, Prediction of telephone calls load using
Echo State Network with exogenous variables, Neural Networks 71 (2015) 204–213.
[21] F. M. Bianchi, E. De Santis, A. Rizzi, A. Sadeghian, Short-term electric load forecasting using echo
state networks and PCA decomposition, IEEE Access 3 (2015) 1931–1943.
[22] M. Skowronski, J. Harris, Minimum mean squared error time series classification using an echo state
network prediction model, in: Proc. 2006 IEEE International Symposium on Circuits and Systems
(ISCAS), IEEE, 2006.
[23] Q. Ma, L. Shen, W. Chen, J. Wang, J. Wei, Z. Yu, Functional echo state network for time series
classification, Information Sciences 373 (2016) 1–20.
[24] F. Palumbo, C. Gallicchio, R. Pucci, A. Micheli, Human activity recognition using multisensor data
fusion based on reservoir computing, Journal of Ambient Intelligence and Smart Environments 8 (2)
(2016) 87–107.
[25] P. Tanisaro, G. Heidemann, Time series classification using time warping invariant echo state networks,
in: Proc. 2016 15th IEEE International Conference on Machine Learning and Applications (ICMLA),
IEEE, 2016, pp. 831–836.
[26] H. Chen, F. Tang, P. Tino, X. Yao, Model-based kernel for efficient time series analysis, in: Proceedings
of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining, ACM,
2013, pp. 392–400.
[27] H. Chen, P. Tino, A. Rodan, X. Yao, Learning in the model space for cognitive fault diagnosis, IEEE
Transactions on Neural Networks and Learning Systems 25 (1) (2014) 124–136.

17

[28] W. Aswolinskiy, R. Reinhart, J. Steil, Time series classification in reservoir-and model-space: a comparison, in: IAPR Workshop on Artificial Neural Networks in Pattern Recognition, Springer, 2016, pp.
197–208.
[29] W. Aswolinskiy, R. Reinhart, J. Steil, Time series classification in reservoir-and model-space, Neural
Processing Letters (2017) 1–21.
[30] F. M. Bianchi, S. Scardapane, S. Løkse, R. Jenssen, Bidirectional deep-readout echo state networks, in:
European Symposium on Artificial Neural Networks, 2018.
[31] I. J. Goodfellow, D. Warde-Farley, M. Mirza, A. C. Courville, Y. Bengio, Maxout networks., Proc. 30th
International Conference on Machine Learning (ICML).
[32] S. Scardapane, S. Van Vaerenbergh, S. Totaro, A. Uncini, Kafnets: kernel-based non-parametric activation functions for neural networks, arXiv preprint arXiv:1707.04035.
[33] F. M. Bianchi, E. Maiorino, M. C. Kampffmeyer, A. Rizzi, R. Jenssen, Recurrent Neural Networks for
Short-Term Load Forecasting: An Overview and Comparative Analysis, Springer, 2017.
[34] I. Sutskever, O. Vinyals, Q. V. Le, Sequence to sequence learning with neural networks, in: Advances
in neural information processing systems, 2014, pp. 3104–3112.
[35] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, R. Salakhutdinov, Dropout: A simple way to
prevent neural networks from overfitting, The Journal of Machine Learning Research 15 (1) (2014)
1929–1958.
[36] W. Zaremba, I. Sutskever, O. Vinyals, Recurrent neural network regularization, arXiv preprint
arXiv:1409.2329.
[37] Y. Bengio, P. Simard, P. Frasconi, Learning long-term dependencies with gradient descent is difficult,
IEEE Transactions on Neural Networks 5 (2) (1994) 157–166.
[38] R. Pascanu, T. Mikolov, Y. Bengio, On the difficulty of training recurrent neural networks, in: International Conference on Machine Learning, 2013, pp. 1310–1318.
[39] S. Hochreiter, J. Schmidhuber, Long short-term memory, Neural computation 9 (8) (1997) 1735–1780.
[40] K. Cho, B. Van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, Y. Bengio, Learning
phrase representations using rnn encoder-decoder for statistical machine translation, arXiv preprint
arXiv:1406.1078.
[41] A. Rodan, P. Tino, Minimum complexity echo state network, IEEE Transactions on Neural Networks
22 (1) (2011) 131–144.
[42] F. M. Bianchi, L. Livi, C. Alippi, Investigating echo-state networks dynamics by means of recurrence
analysis, IEEE Transactions on Neural Networks and Learning Systems 29 (2) (2018) 427–439.
[43] A. Y. Ng, M. I. Jordan, On discriminative vs. generative classifiers: A comparison of logistic regression
and naive bayes, in: Advances in neural information processing systems, 2002, pp. 841–848.
[44] T. Jaakkola, D. Haussler, Exploiting generative models in discriminative classifiers, in: Advances in
neural information processing systems, 1999, pp. 487–493.
[45] M. Belkin, P. Niyogi, Laplacian eigenmaps for dimensionality reduction and data representation, Neural
computation 15 (6) (2003) 1373–1396.
[46] S. Mika, B. Schölkopf, A. J. Smola, K.-R. Müller, M. Scholz, G. Rätsch, Kernel pca and de-noising in
feature spaces, in: Advances in neural information processing systems, 1999, pp. 536–542.
18

[47] S. Løkse, F. M. Bianchi, R. Jenssen, Training echo state networks with regularization through dimensionality reduction, Cognitive Computation 9 (3) (2017) 364–378.
[48] W. Maass, T. Natschläger, H. Markram, Real-time computing without stable states: A new framework
for neural computation based on perturbations, Neural computation 14 (11) (2002) 2531–2560.
[49] K. Bush, C. Anderson, Modeling reward functions for incomplete state representations via echo state
networks, in: Proc. International Joint Conference on Neural Networks (IJCNN), Vol. 5, IEEE, 2005,
pp. 2995–3000.
[50] Š. Babinec, J. Pospı́chal, Merging echo state and feedforward neural networks for time series forecasting,
in: International Conference on Artificial Neural Networks, Springer, 2006, pp. 367–375.
[51] I. Goodfellow, H. Lee, Q. V. Le, A. Saxe, A. Y. Ng, Measuring invariances in deep networks, in:
Advances in neural information processing systems, 2009, pp. 646–654.
[52] X. Glorot, Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, in:
Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, 2010,
pp. 249–256.
[53] L. R. Tucker, Some mathematical notes on three-mode factor analysis, Psychometrika 31 (3) (1966)
279–311.
[54] T. G. Kolda, B. W. Bader, Tensor decompositions and applications, SIAM review 51 (3) (2009) 455–500.
[55] C. F. Beckmann, S. M. Smith, Tensorial extensions of independent component analysis for multisubject
fmri analysis, Neuroimage 25 (1) (2005) 294–311.
[56] D. Zhang, Z.-H. Zhou, (2d) 2pca: Two-directional two-dimensional pca for efficient face representation
and recognition, Neurocomputing 69 (1-3) (2005) 224–231.
[57] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980.
[58] M. G. Baydogan, G. Runger, Time series representation and similarity based on local autopatterns,
Data Mining and Knowledge Discovery 30 (2) (2016) 476–509.
[59] K. H. Brodersen, T. M. Schofield, A. P. Leff, C. S. Ong, E. I. Lomakina, J. M. Buhmann, K. E. Stephan,
Generative embedding for model-based classification of fmri data, PLoS computational biology 7 (6)
(2011) e1002079.
[60] U. Siewert, W. Wustlich, Echo-state networks with band-pass neurons: Towards generic time-scaleindependent reservoir structures, Internal status report, PLANET intelligent systems GmbH.

19

