<figure><img src="image_1.png" /></figure>

ifip

International Federation for Information Processing

This document is the original author manuscript of a paper submitted to an IFIP
conference proceedings or other IFIP publication by Springer Nature. As such, there
may be some differences in the official published version of the paper. Such
differences, if any, are usually due to reformatting during preparation for publication or
minor corrections made by the author(s) during final proofreading of the publication
manuscript.
# Principal Component Modes of Reservoir Dynamics in Reservoir Computing

José María Enguita<fnref n="1" />(✉), Ignacio Díaz<fnref n="1" />, Diego García<fnref n="1" />, Abel Alberto Cuadrado<fnref n="1" />, and José Ramón Rodríguez<fnref n="2" />

**Abstract.** In Reservoir Computing, signals or sequences are fed into a set of interconnected non-linear units (neurons) with capabilities for storing information (reservoir). The reservoir generates an expanded representation of the input, which is subsequently mapped onto the desired output using a trained output layer (readout). However, despite their success in various experimental tasks, the dynamics of the reservoir are not yet well understood. In this paper we introduce a new technique, based on the well known Singular Value Decomposition (SVD), to obtain the main dynamic modes of the reservoir when excited with an input signal. We conduct experiments using Echo State Networks (ESN) to demonstrate the technique's potential and its ability to decompose input signals into Principal Component Modes as expanded by the reservoir. We expect that this approach will open new possibilities in its application to the field of visual analytics in process state visualisation, determination of attribute vectors, and detection of novelties. Furthermore, this technique could serve as a foundation for a better understanding of the reservoir's dynamic state that could help in other areas of research, such as domain shift or continual learning.

**Keywords:** Reservoir Computing · Echo State Networks · Reservoir Dynamics · Singular Value Decomposition · Visual Analytics

## 1 Introduction

### 1.1 Reservoir Computing

Reservoir computing (RC) is a computational framework based on *Recurrent Neural Network* (RNN) models which extends the concept of earlier neural network architectures such as Liquid-State Machines and Echo State Networks [15]. RC is highly effective for processing sequential or temporal data, such as the behaviour of dynamic systems or time series. The input signals are fed into a set

<footnote n="1">¹ Department of Electrical Engineering, University of Oviedo, 33204 Gijón, Spain
{jmenguita,idiaz,garciaperdiego,cuadradoabel}@uniovi.es</footnote>

<footnote n="2">² SUPPRESS Research Group, University of León, 24007 León, Spain
jrodro@unileon.es
http://isa.uniovi.es/GSDPI/</footnote>
of interconnected neurons, referred to as the ‘reservoir’, which generates a high-dimensional expansion [13]. The reservoir contains feedback (recurrent) loops that provide memory and exhibits a complex and highly non-linear dynamic behaviour [9]. A simple readout system is trained to map the internal state of the reservoir to the desired output.

RC systems offer a key advantage over other recurrent neural networks in that they simplify the learning process. Specifically, as the reservoir remains fixed, only the readout layer is trained. This method requires a simple least-squares estimation rather than the more expensive and computationally intensive non-linear optimisation required to train an RNN, reducing the computational cost and facilitating rapid learning. Remarkably, the forecasting capability of RC systems can still be competitive even when dealing with chaotic or spatio-temporally complex problems [1].

Despite all the very convincing experimental results on a variety of tasks, a thorough understanding of the dynamics of the reservoir is still lacking [16]. This contribution aims to extend this understanding by presenting a simple yet effective approach of characterising the reservoir dynamics by means of the *Singular-Value Decomposition* (SVD), which has already been effectively employed in various related applications such as identifying underlying patterns and detecting outliers in time series [10], configuring the weight matrix of a reservoir with specific spectral characteristics [11], and obtaining an optimal low-rank representation of a time series expanded by Dynamic Mode Decomposition [3,4]. In addition, Gallicchio et al. [6] used a similar approach to investigate redundancy in the reservoir.

Our proposal uses the SVD over the evolution in time of the activation values of the neurons in the reservoir when excited by an input signal. We show that this technique yields an interpretable description of the reservoir dynamics and provides a decomposition in terms of principal component modes in a lower-dimensional space.

# 2 Methods

## 2.1 Recurrent Neural Networks

*Recurrent Neural Networks* (RNNs) are rooted on the works by Rumelhart, Hinton and Williams [14] back in 1986. RNNs are a type of artificial neural network in which connections introduce a non-linear feedback by connecting the node outputs with their inputs, thus enabling the network to display temporal dynamic behaviour.

The evolution of the states of the hidden layer neurons $\mathbf{x}(k)$ and the output of the network $\mathbf{y}(k)$ for a given input signal $\mathbf{u}(k)$, follows:

$$ 
\begin{aligned}
\mathbf{x}(k) &= f(\mathbf{x}(k-1), \mathbf{u}(k)) \\
\mathbf{y}(k) &= h(\mathbf{x}(k)).
\end{aligned} 
\quad (1) $$

where $f$ and $h$ are general functions, usually non-linear.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Representation of an ESN.</figcaption></figure>

## 2.2 Echo State Networks

Although the methods presented in this paper are general to reservoir computing systems, we implement the results on the well known *Echo State Networks* (ESN) architecture. Nevertheless, it should be easily adapted to other models.

Proposed by Jaeger and Hass back in the beginning of the 21st century [7–9], ESNs use an sparsely connected RNN with random weights as reservoir.

An ESN, in its basic form, allows the modelling of non-linear systems by supervised learning. ESNs, see Fig. 1, are a particular case of equation (1), based on the concept of *non-linear expansion* [13], in which we consider a high-dimensional state vector $\mathbf{x}(k) \in \mathbb{R}^n$, a non-linear model in the state equation and a linear model of the output obtained from the state<fnref n="1" />:

$$ 
\begin{align} 
\mathbf{x}(k) &= \sigma(\mathbf{W}_{\text{res}} \mathbf{x}(k-1) + \mathbf{W}_{\text{in}} \mathbf{u}(k)) \tag{2} \\ 
\mathbf{y}(k) &= \mathbf{W}_{\text{out}} \mathbf{x}(k). 
\end{align} 
$$

The model represented by equation (2) includes the reservoir matrix $\mathbf{W}_{\text{res}} \in \mathbb{R}^{n \times n}$, the input matrix $\mathbf{W}_{\text{in}} \in \mathbb{R}^{n \times p}$, and the output matrix $\mathbf{W}_{\text{out}} \in \mathbb{R}^{q \times n}$ as parameters. A non-linear function $\sigma$, typically sigmoidal, is used in the state equation, and sparsity is ensured in $\mathbf{W}_{\text{res}}$ using a small fraction of the possible feedback connections among nodes. The matrices $\mathbf{W}_{\text{res}}$ and $\mathbf{W}_{\text{in}}$ are generated with random values, with some modifications applied to bring the system to the edge of stability. Only the matrix $\mathbf{W}_{\text{out}}$ is learned to reproduce specific time patterns.

In order for the ESN network to work, the reservoir must satisfy the *echo state property*. This property links the asymptotic characteristics of the excited reservoir dynamics to the driving signal, so that the initial state of the reservoir

<footnote n="1">¹ The original version is given here. Numerous variants exist in the literature incorporating, for example, a bias term added to $\mathbf{u}(k)$, direct effect of $\mathbf{u}(k)$ on the output, feedback term (inclusion of $\mathbf{y}(k)$) in the equation of state, application of a low-pass filter to the states, use of non-linear regression models to obtain $\mathbf{y}(k)$ from $\mathbf{x}(k)$, etc. A detailed description of many of these variants can be found in [12].</footnote>
is asymptotically 'washed out'. This property is satisfied, for additive-sigmoid neuron reservoirs, if the spectral radius $\rho$ of the reservoir matrix (maximum of the absolute value of its eigenvalues) is less than 1 [2], although some authors reported good results with spectral radius slightly over this value [12].

In the above equations, both $\mathbf{u}(k)$ and $\mathbf{y}(k)$ are vectors because multiple inputs and outputs can be considered. However, to keep the notation simpler and without loss of generality, we will restrict to only single-input and single-output systems for the rest of the paper, so that $\mathbf{y}$ and $\mathbf{u}$ are time series, $u(k)$, $y(k)$ are elements of each series, and $\mathbf{W}_{out} \in \mathbb{R}^{1 \times n}$.

## 2.3 Analysis of the Reservoir Dynamics

By applying the model in equation (2) recursively, an input or excitation sequence $\mathbf{u}$ generates a sequence of state vectors containing the activation values of the $n$ neurons of the reservoir, $\mathbf{x}(k)$ which, for the case of $m$ samples, can be grouped in a matrix:

$$ \mathbf{X} \in \mathbb{R}^{n \times m} = \begin{pmatrix} | & | & | \\ \mathbf{x}(1) & \mathbf{x}(2) \cdots \mathbf{x}(m) \\ | & | & | \end{pmatrix}. \quad (3) $$

The rows of the matrix $\mathbf{X}$ contain the temporal evolution of the reservoir with a great richness and variety of dynamic behaviours as a response to an input signal. Each state vector $\mathbf{x}(k)$ can be considered an expanded set of descriptors of the dynamics of the input signal at instant $k$.

In this work, we propose to analyse the dynamics of the reservoir through the SVD of the matrix $\mathbf{X}$:

$$ \mathbf{X} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T = \sigma_1 \mathbf{u}_1 \mathbf{v}_1^T + \dots + \sigma_n \mathbf{u}_n \mathbf{v}_n^T, \quad (4) $$

where $\mathbf{u}_i$ is the $i$-th column of matrix $\mathbf{U}$, $\mathbf{v}_i^T$ is the $i$-th row of matrix $\mathbf{V}^T$.

Singular values $\sigma_1, \sigma_2, \dots, \sigma_n$ represent the weights of the principal modes of $\mathbf{X}$, and may be used as descriptors of the dynamic evolution of the reservoir, when excited by the input signal.

The columns of $\mathbf{U}$ are the eigenvectors of $\mathbf{X}\mathbf{X}^T$, the columns of $\mathbf{V}$ are the eigenvectors of $\mathbf{X}^T\mathbf{X}$. As they are ordered by decreasing values of $\sigma_i$, the columns of $\mathbf{U}$ and $\mathbf{V}$ are hierarchically ordered according to the correlation they capture from the columns and rows of $\mathbf{X}$. That's why the former captures the spatial correlation (across the $n$ components of the state vector) and the latter the temporal correlation (across the $m$ time steps). Both are orthonormal bases for the space to which the columns and rows of $\mathbf{X}$ belong, respectively.

If we select only the $r$-first singular values, an approximation of $\mathbf{X}$ can be obtained as $\tilde{\mathbf{X}} = \mathbf{U}_r \mathbf{\Sigma}_r \mathbf{V}_r^T$. In fact, this is the best rank-$r$ approximation of $\mathbf{X}$ in the $L_2$-sense, and it provides a reduced-order latent space in which to study the dynamic evolution of the reservoir when the system is excited with an input signal $\mathbf{u}$.
This space can also be projected onto a lower dimension space, such as 2D, for visualisation purposes. This enables enhanced analysis and interpretation of the system's behaviour over time.

## 2.4 Principal Component Modes

It should be noted that, thus far, the ESN can be analysed for its dynamic behaviour without the need for training. However, a trained ESN can offer interesting possibilities. According to the SVD, in line with (4), the matrix $X$ can be decomposed in $X_i = u_i \sigma_i v_i^T$, so that $X = \sum_i X_i$. The target time series $y$ is approximated by the ESN using $W_{out}$ as:

$$y \approx W_{out}X = W_{out} \sum_{i} u_i \sigma_i v_i^T = \sum_{i} \alpha_i v_i^T = \sum_{i} m_i, \quad (5)$$

where $\alpha_i$ are scalar values, since $W_{out} \in \mathbb{R}^{1,n}$, $u_i \in \mathbb{R}^{n,1}$, and $\sigma_i$ is a scalar. In sum, taking $\alpha_i = W_{out}u_i\sigma_i$ we have that the target time series can be approximated as a sum of modes $m_i$.

The decomposition can also be applied locally on sliding submatrices $X(k) \in \mathbb{R}^{n \times m'}$, with the columns between $k$ and $k+m' - 1$ of $X$, to track the evolution of the reservoir when excited by the input signal. This allows to keep track of any variation in the modes $m_i$, for instance, in process visual monitoring.

There are two primary configurations of operation depending on how the readout is trained:

1. **Autoencoder:** The ESN is trained using the input signal. Then, by utilising the trained readout, the input signal can be decomposed into a hierarchy of dynamic modes as observed by the reservoir:

$$W_{out}^{AE} = \arg \min_{W} \|WX - u\| \quad \text{and} \quad m_i = W_{out}^{AE} X_i. \quad (6)$$

2. **Transcoder:** The ESN is trained using the output signal. In this scenario, it is the output signal which is decomposed into its dynamic modes:

$$W_{out}^{TC} = \arg \min_{W} \|WX - y\| \quad \text{and} \quad m_i = W_{out}^{TC} X_i. \quad (7)$$

# 3 Results and Discussion

## 3.1 Second Order System

In this experiment we simulate a generic second order system:

$$\frac{d^2 y(t)}{dt^2} + 2\zeta\omega_n \frac{dy(t)}{dt} + \omega_n^2 y(t) = K\omega_n^2 u(t), \quad (8)$$

with gain $K = 1$, damping ratio $\zeta = 0.5$ (underdamped) and natural frequency $\omega_n = 1$, for 1200 s seconds. Figure 2 shows an example of hierarchical decomposition in transcoder configuration of the main principal components. The system
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Plot of the first 4 modes ($\mathbf{m}_1, \dots, \mathbf{m}_4$) in which the response of the system to a square signal is decomposed by the ESN. X-axis is time in seconds.</figcaption></figure>

is excited with a square signal and the output is used to train an ESN with 300 neurons in the hidden layer<fnref n="2" />. Figure 3 shows the reconstruction of the output using the first modes in the decomposition.

To analyse the possibilities for local decomposition and tracking, the damping ratio is varied, taking the values 0.2, 0.3, 0.4 and 0.5 during 300 s each. The response to a square signal is shown in Fig. 4.

An ESN is fed with this response signal using a sliding window of size $m' = 800$ samples and a stride value of 20 samples (as we are analysing the reservoir only, the readout doesn't need to be trained). As the ESN has $n = 300$ neurons in its hidden layer, the signal is decomposed into 300 modes, corresponding to the nonzero singular values $\sigma_1, \dots, \sigma_{300}$, at each step. Projecting these 300-dimensional vectors of singular values into a 2D plot with a simple principal component analysis (PCA) algorithm, the trajectories in the latent space (from the high state of the square signal on the right, to its low state on the left) corresponding to each damping-ratio value can be clearly observed (see Fig. 5, left).

In this example, the input signal has a period of 500 samples. As the window size approaches that value, the trajectories become shorter and better separated until they collapse into 4 distinct positions in the projected map (see Fig. 5, right). Some transition dots still appear in the plot, corresponding to sections where the sliding window partially overlaps with the output of two different system states. These trajectories show the transitions between states, but are, in fact, an artefact of the sliding window mechanism.

<footnote n="2">² Other parameters are: $\rho = 0.95$, sparsity = 0.01, leaking rate = 0.025, input scaling = 10 and warm up = 20. A bias is added to the inputs.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Reconstruction of the response to a square signal by the ESN using the first 2, 8, 15 and 50 modes (x-axis is time in seconds). The original signal is marked with dots for better reference.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Response to a square signal. The damping ratio is changed every 300 s, so that $\zeta = 0.2, 0.3, 0.4, 0.5$.</figcaption></figure>

Thanks to the fact that each window corresponds to a single period of the signal, each stride of the window effectively performs a circular shift of the columns of the matrix $\mathbf{X}$. Any row permutation can be represented as $\mathbf{X}_{perm} = \mathbf{X}P$, where $\mathbf{P}$ is a permutation matrix, therefore orthogonal. If $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T$ is a valid SVD decomposition of $\mathbf{X}$, then $(\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T)\mathbf{P} = \mathbf{U}\mathbf{\Sigma}(\mathbf{V}^T\mathbf{P})$ is a valid SVD decomposition of $\mathbf{X}_{perm}$. As a result, the singular values remain unchanged, and the rows of matrix $\mathbf{V}^T$ are permuted by $\mathbf{P}$. All values in both $\mathbf{V}^T$ and $\mathbf{U}$ that correspond to non zero singular values remain unaltered, except for a potential change in sign, as both matrices are real.

As a final note in this example, the magnitude of the main singular value, $\sigma_1$, can differentiate among the states of the system, as is shown in Fig. 6.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Latent space trajectories projected to 2D for the different damping ratio values (both axis represent coordinates in the latent space). Left, with a sliding window of size 800 samples (high state of the input signal corresponds to the right and low state to the left). Right, the same trajectories with a window size equal to the period of the signal (500 samples). Four distinct positions, related to the four modes of the system, can be clearly seen.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Plot of the singular value $\sigma_1$ when the sliding window size equals the period of the signal. The different dynamics corresponding to different damping ratios can be clearly observed.</figcaption></figure>

## 3.2 Dynamic Modes of Vibration Signals in an Induction Motor Under Test

For this experiment we used an available dataset [5] containing vibration and current data of a 4kW induction motor with 6306-2Z/C3 bearings that rotates at 1500 rpm (25 Hz) with a supply frequency 50 Hz. This machine has been subjected to nine different tests (Table 1), for which five operating variables have been measured at a sampling frequency of 5000 Hz: three vibrations (measured as accelerations $a_c$, $a_x$ and $a_y$) and two phase currents ($i_R$ and $i_S$). We have selected the horizontal vibration ($a_x$) for this experiment.

The vibration data were fed into the same ESN as before, but with an input scaling factor of 1. Data were analysed with a sliding window with a length of 1000 samples, and a stride of 100 samples.

The projection of the singular values over a 2D space quickly showed the operation modes and the trajectories indicated how the system varied from one test to the next. This can be seen in Fig. 7. Tests 8 and 9 consisted in gradual variations from one electrical fault state to another, which can be seen in their
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Latent space trajectories projected to 2D. The different locations of the operational states in tests 1 to 7 were labelled for reference.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Top: input signal from the horizontal accelerometer. Bottom: evolution of the singular value $\sigma_2$ throughout the tests.</figcaption></figure>

respective trajectories. The former starts at normal operating mode, moves passing through all the electrical fault modes, and gets back to normal operation. The latter, however, starts from normal operation and just moves through all the electrical fault modes. A visual analytics application could help the user to keep track of the operation of the motor just by observing the trajectories. Moreover, the singular value $\sigma_2$ is able to differentiate between all the states (see Fig. 8), meaning it could be used for automatic fault diagnosis.

By training the ESN in autoencoder configuration, it is possible to decompose the vibration signal into Principal Component Modes. The result is shown in Fig. 9. A portion of the signal of 200 ms (which equals the used window size) is
<table><caption>**Table 1.** Description of performed tests.</caption><thead><tr><td>Number</td><td>Performed test</td><td>Duration (s)</td></tr></thead><tbody><tr><td>1</td><td>Mechanical fault (eccentric mass on pulley)</td><td>4</td></tr><tr><td>2</td><td>Combined electrical and mechanical fault</td><td>4</td></tr><tr><td>3</td><td>Normal operation</td><td>4</td></tr><tr><td>4</td><td>Electrical fault (15 Ω resistor in phase R)</td><td>4</td></tr><tr><td>5</td><td>Electrical fault (20 Ω resistor in phase R)</td><td>4</td></tr><tr><td>6</td><td>Electrical fault (25 Ω resistor in phase R)</td><td>4</td></tr><tr><td>7</td><td>Electrical fault (5 Ω resistor in phase R)</td><td>4</td></tr><tr><td>8</td><td>Gradual electrical fault (resistance increases and decreases)</td><td>4</td></tr><tr><td>9</td><td>Gradual electrical fault (gradual increase of resistance)</td><td>8</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** First Principal Component Modes of the vibration for different tests. A section of 200 ms of the signal (the used window size) at each test segment is displayed. Modes have been centred for clarity.</figcaption></figure>

shown along with the first three modes. Interestingly, mode 3 carries significant information about the vibration caused by each type of fault. It can be observed that mode 3 captures the most representative harmonics of each vibration state. In the case of a mechanical fault, a frequency 25 Hz originating from the eccentric mass is observed, whereas an electrical fault produces a 100 Hz harmonic as a result of electrical asymmetry. It is noteworthy that the combination of both electrical and mechanical faults yields similar proportions of the aforementioned harmonics.
# 4 Conclusion

The SVD decomposition of the temporal evolution of the reservoir provides valuable information about its dynamic evolution, creating a latent space of dimension $r \le n$ that can be used to observe trajectories, search for attribute vectors, detect novelties, or project onto a visualisation space to apply visual analytics techniques. Unlike Fourier analysis, which uses an imposed basis of sine wave signals, the SVD performs a decomposition into *principal* components, those that contain the most energy, potentially resulting in more relevant and interpretable information.

This idea can be applied to quasi-periodic signals or variant systems, as long as the echo state property is fulfilled. This does not result in any additional restrictions, as such property must be fulfilled for the ESN to function properly.

Training the output layer of the network is not a requirement, however, if this layer is trained, it becomes possible to reconstruct the principal modes of the temporal series dynamics from the reservoir's perspective. Two different modes have been defined: the *autoencoder* configuration, where the network is trained with the input signal itself, and the *transcoder* configuration, where the network is trained to learn a dynamical system between the input excitation and a target output signal. In both cases, the obtained information can be used in multiple ways, such as searching for the principal modes for monitoring, predictive maintenance, or novelty detection in a process.

Although its main purpose is analytical, we believe that this technique can also serve as a basis for developing new methods that allow analysing the dynamic characteristics of the reservoir and determining its suitability for specific problems. Moreover, the location of the operating point in the latent space or the information obtained from certain principal modes can be used to address other research issues in this type of networks, such as domain shift or continual learning, both of which are topics of interest in the field of machine learning in general and in reservoir computing in particular.

**Acknowledgements.** This work was supported by the Ministerio de Ciencia e Inno-
vación / Agencia Estatal de Investigación (MCIN/AEI/ 10.13039/501100011033) grant
[PID2020-115401GB-I00].

# References

1. Bollt, E.: On explaining the surprising success of reservoir computing forecaster of chaos? The universal machine learning dynamical system with contrast to VAR and DMD. Chaos: Interdiscip. J. Nonlinear Sci. **31**, 013108 (2021). https://doi.org/10.1063/5.0024890. http://aip.scitation.org/doi/10.1063/5.0024890

2. Buehner, M., Young, P.: A tighter bound for the echo state property. IEEE Trans. Neural Netw. **17**, 820–824 (2006). https://doi.org/10.1109/TNN.2006.872357

3. Dylewsky, D., Barajas-Solano, D., Ma, T., Tartakovsky, A.M., Kutz, J.N.: Stochastically forced ensemble dynamic mode decomposition for forecasting and analysis of near-periodic systems. IEEE Access **10**, 33440–33448 (2022). https://doi.org/10.1109/ACCESS.2022.3161438
4. Dylewsky, D., Kaiser, E., Brunton, S.L., Kutz, J.N.: Principal component trajectories for modeling spectrally continuous dynamics as forced linear systems. Phys. Rev. E **105**, 015312 (2022). https://doi.org/10.1103/PHYSREVE.105.015312/FIGURES/10/MEDIUM. https://journals.aps.org/pre/abstract/10.1103/PhysRevE.105.015312

5. Díaz Blanco, I., Cuadrado Vega, A.A., Muñiz, A.G., García Pérez, D.: Dataicann: datos de vibración y corriente de un motor de inducción. https://digibuo.uniovi.es/dspace/handle/10651/53461 (2019)

6. Gallicchio, C., Micheli, A.: A Markovian characterization of redundancy in echo state networks by PCA. In: Proceedings of the 18th European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning (ESANN-2010), pp. 321–326 (2010)

7. Jaeger, H.: The “echo state” approach to analysing and training recurrent neural networks. GMD Report 148, GMD - German National Research Institute for Computer Science (2001). http://www.faculty.jacobs-university.de/hjaeger/pubs/EchoStatesTechRep.pdf

8. Jaeger, H.: Echo state network. Scholarpedia **2**(9), 2330 (2007)

9. Jaeger, H., Haas, H.: Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication. Science **304**, 78–80 2004). https://doi.org/10.1126/science.1091277. https://www.science.org/doi/10.1126/science.1091277

10. Khoshrou, A., Pauwels, E.J.: Data-driven pattern identification and outlier detection in time series. Adv. Intell. Syst. Comput. **858**, 471–484 (2018). https://doi.org/10.1007/978-3-030-01174-1-35. http://arxiv.org/abs/1807.03386

11. Li, F., Wang, X., Li, Y.: Effects of singular value spectrum on the performance of echo state network. Neurocomputing **358**, 414–423 (2019). https://doi.org/10.1016/j.neucom.2019.05.068. https://linkinghub.elsevier.com/retrieve/pii/S092523121930774X

12. Lukoševičius, M.: A practical guide to applying echo state networks. In: Montavon, G., Orr, G.B., Müller, K.-R. (eds.) Neural Networks: Tricks of the Trade. LNCS, vol. 7700, pp. 659–686. Springer, Heidelberg (2012). https://doi.org/10.1007/978-3-642-35289-8_36

13. Lukoševičius, M., Jaeger, H.: Reservoir computing approaches to recurrent neural network training. Comput. Sci. Rev. **3**(3), 127–149 (2009)

14. Rumelhart, D.E., Hinton, G.E., Williams, R.J.: Learning representations by back-propagating errors. Nature **323**, 533–536 (1986). https://doi.org/10.1038/323533a0. https://www.nature.com/articles/323533a0

15. Tanaka, G., et al.: Recent advances in physical reservoir computing: a review. Neural Netw. **115**, 100–123 (2019). https://doi.org/10.1016/J.NEUNET.2019.03.005

16. Verstraeten, D., Schrauwen, B.: On the quantification of dynamics in reservoir computing. In: Alippi, C., Polycarpou, M., Panayiotou, C., Ellinas, G. (eds.) ICANN 2009. LNCS, vol. 5768, pp. 985–994. Springer, Heidelberg (2009). https://doi.org/10.1007/978-3-642-04274-4_101