{'fraction_non_alphanumeric': 0.0403535195657734, 'fraction_numerical': 0.012873341100912216, 'mean_word_length': 4.393930104230534, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interpretable Design of Reservoir Computing Networks using Realization Theory

The reservoir computing networks (RCNs) have been successfully employed as a tool in learning and complex decision-making tasks. Despite their efficiency and low training cost, practical applications of RCNs rely heavily on empirical design. In this paper, we develop an algorithm to design RCNs using the realization theory of linear dynamical systems. In particular, we introduce the notion of $\alpha$-stable realization, and provide an efficient approach to prune the size of a linear RCN without deteriorating the training accuracy. Furthermore, we derive a necessary and sufficient condition on the irreducibility of number of hidden nodes in linear RCNs based on the concepts of controllability and observability matrices. Leveraging the linear RCN design, we provide a tractable procedure to realize RCNs with nonlinear activation functions. Finally, we present numerical experiments on forecasting time-delay systems and chaotic systems to validate the proposed RCN design methods and demonstrate their efficacy.


I. INTRODUCTION
The reservoir computing network (RCN) is a bio-mimetic computational tool that is increasingly used in a variety of applications to solve complex decision making problems [1]- [3]. Essentially, the RCN is a class of recurrent neural networks (RNNs), which is composed of one hidden layer, typically with a large number of sparsely interconnected neurons, and a linear output layer. In contrast to the classical RNN, a distinct feature of the RCN is that all of its connections in the hidden-layer are randomly pre-determined and fixed. Hence, the training process of the RCN involves only learning the weights of its linear output-layer in a supervised learning framework.
The existing supervised learning approach to training an RCN was proposed in [1]. Subsequently, the RCN was successfully employed for forecasting time-series with applications in finance [4], [5], wireless communication [6], speech recognition [7], and robot navigation [8]. Notwithstanding its efficient training procedure, the major limitation of the RCN lies in the fact that its practical application relies heavily on empirical design of the hyper-parameters of the network, including its size [9].
Recently, there has been a renewed interest in developing tractable methods for designing neural networks that are suitably deployed in diverse scenarios [10], [11]. In this context, deriving rigorous and systematic techniques to design neural networks, especially establishing principled strategies for selecting their hyper-parameters that yield a desired performance, is compelling but challenging. In this paper, we propose a tractable approach to design RCNs that warrant effective functioning for given datasets. In particular, we focus on the application of RCNs to learn dynamic models of dynamical systems from their time-series measurement data, and develop rigorous design principles to prune the number of hidden-layer nodes in RCNs without deteriorating the training performance. Leveraging the notions of controllability and observability matrices, we derive a necessary and sufficient condition on the irreducibility of RCNs with linear activation functions. This in turn results in an interpretable RCN pruning algorithm, where the RCNs' controllability and observability matrices inform on its size and irreducibility. Furthermore, we illustrate that the developed irreducible linear realization of the RCN not only sufficiently represents the underlying dynamics inherited in the time-series data, but also contributes to a tractable design of general RCNs with nonlinear activation functions.
The paper is organized as follows. In Section II, we provide a brief review of related works that motivate the need of our developments. In Section III, through tailoring existing results on the learnability of RNNs, we motivate our realizationtheoretic RCN design principles by illustrating how an RCN achieves τ -step ahead forecast of time-series associated with a dynamical system. In Section IV, we introduce realizationtheoretic aspects from systems theory to facilitate a comprehensive RCN design, and then establish a necessary and sufficient condition on the irreducibility of a linear RCN that achieves desired training accuracy. This result in turn forms the basis to an educated design of RCNs with nonlinear activation functions. In Section V, we present several results using an RCN to forecast time-series and learn chaotic systems to demonstrate the applicability of the proposed RCN design framework.

II. RELATED WORKS AND MOTIVATION
In this section, we briefly review the existing works on RCNs and point out the specific problems that we address in this paper.
The computational framework of an RCN and its training procedure were first proposed in [1], where two conditions were hypothesized as requirements for successful applications of the RCNs -the echo state property (ESP) and a general compactness assumption on the training signal. In addition, a mathematical definition of the ESP was also provided in [12], [13]. Intuitively, the ESP implies that the state of the RCN is uniquely determined by its input history rather than the initial condition of the network.
Thereafter, several results explaining the principles of the RCN, especially in evaluating some of its features, such as the ESP [14], [15], the memory capacity [16]- [18], and the stability [19], have been reported. In addition to investigating the fundamental properties of the RCN, multiple attempts addressing its design were also reported. The performance of the RCN in relation to the complexity of its network topology was analyzed in [20]. Using the interpretation of contracting maps, a discussion on the architecture of the RCN was presented in [21]. Supported by heuristic analyses, it was shown in [22] that the small-world network topology improved the performance (e.g., forecasting accuracy or memory capacity) of the RCN. In spite of the prescribed results, existing applications using RCNs rely on randomly generated connection matrices. Using the random matrix theory, an explanation on why such a design of the RCN, in general, leads to acceptable performance was presented [23].
Recently, there has been a rising tide of interest in analyzing RCNs, especially RCNs with linear activation functions, using control-theoretic approaches [24], [25]. For instance, the connectivity patterns of RCNs were studied in [26] leveraging the concept of controllability matrix. Furthermore, it was proved in [27] that the memory capacity of a linear RCN can be characterized by the rank of its controllability matrix. Nevertheless, due to the inherited randomness and usage, the design process of the RCNs are still based on empirical strategies, and designing an RCN with minimum size to ensure a desired performance is compelling but remains elusive.
In this work, we focus on establishing schematic design principles to realize RCNs with linear and nonlinear activation functions, which achieve guaranteed training accuracy. The main contributions of this paper include: (1) establishing the notion of α-stable realization for designing the weight/connection matrix in an RCN with desired training accuracy for a given dataset; (2) devising an algorithm based on realization theory to prune the size of linear RCNs with quantifiable training accuracy; and (3)

III. ROLE OF TAKENS EMBEDDING IN RCN FRAMEWORKS
In this section, we provide details of the RCN dynamics and its training procedure. We tailor existing results on the learnability of RNNs, in particular, the Takens embedding theorem, to render a comprehensive analysis of how an RCN learns the underlying dynamics of a time-series. Based on the analyses, we discuss the applications of RCNs for the time series forecasting problem, which motivates the design of RCNs using linear realization theory in Section IV. We begin with a brief introduction to the Takens theorem and discuss its role in understanding RCNs.

A. Takens theorem and its implication on time-series forecasting
We consider a time-dependent variable s(t), evolving on an m-dimensional manifold M ⊂ R p , following the dynamicṡ s(t) = F s (s(t)), where F s : M → R p is a smooth vector field. Let v : M → R be an observation function, and in practice, we measure a discrete sequence of observations, say v[s(t i )], where t i for i = 0, 1, 2, . . . denoting the sampling instants. We can then define the propagation map φ : M → M describing the flow of s(t) at time t i by s(t i+1 ) = φ(s(t i )). Theorem 1 (Takens theorem [28]). Let M be a com- generic' means open and dense in C 1 topology.

B. The RCN dynamics and its training procedure
In this part, we introduce the dynamics of the RCN and a sufficient condition for the RCN to possess the ESP. With the guarantee of the ESP and using Takens theorem, we show that an RCN can 'learn' to forecast the data generated by a dynamical system on a compact manifold.
Consider the RCN described by where x[k] ∈ R N denotes the state of the RCN and N is the number of nodes; u[k] ∈ R p is the input to the RCN; y[k] ∈ R q is the output of the RCN; and A ∈ R N ×N and B ∈ R N ×p are pre-determined matrices denoting the connections in the RCN; C ∈ R q×N is the weight matrix in output-layer to be trained, with c i , i = 1, . . . , N as its i th column; α ∈ (0, 1) is called the 'leakage rate', and σ is an activation function that is applied to a vector component-wise. A well-known result (see [1], [30]) for the RCN in (1) to possess ESP is provided as the following lemma. for k > T . The canonical way to do this is to first select the size of the RCN (i.e., N ), and randomly generate the matrices A and B of appropriate dimensions. Then, the sequence u[0; T − 1] is fed into the RCN dynamics (1) as an input to generate a sequence of the RCN states x[1; T ]. A fixed positive integer w is selected as the 'washout' length, and only the sequence after the w th step, i.e., x[w; T ], is collected. Finally, the coefficient in the output layer, C, is trained to minimize the error between the RCN output y[w; T ] and the reference sequenceỹ[w; T ], that is, Following this supervised learning procedure, when the input to the RCN is u[k] (for k > T ), the RCN outputs a value that approximatesỹ[k + 1]. In this sense, the training process enables the RCN to learn the underlying dynamics of the reference sequenceỹ[k].
In the following section, we explain in detail on how the RCN encodes the underlying dynamics of the reference sequence during the training process through the lens of Takens embedding theorem.

C. Learning dynamics using an RCN
To begin with, we consider the task of 1-step ahead forecast of a given time-series and explain how the training process enables the RCN to perform this task. For ease of exposition, we illustrate the idea with one-dimensional time-series, i.e., u[k], y[k] ∈ R, and the framework is directly applicable to the multi-dimensional cases since the Takens theorem holds regardless of the dimension of the time-series. 1) 1-step ahead forecast: Suppose a time-series {u[k]} ⊂ R is generated by a dynamical system on a compact manifold of dimension m. The sequence u[0; T ] is used as input to the RCN, and the 1-step shifted sequence u[1; T + 1] is provided as the training reference. Denote the solution for the state equation in (1) as , then by the uniqueness of the solution of a dynamical system, we have As a result of Lemma 1, when A 2 < 1 L , the RCN acquires the ESP. This implies that there exists a k 0 ∈ N such that, after the RCN evolves for k 0 steps, the effect of the initial condition on the solution trajectory becomes negligible. Therefore, for a fixed 'washout' length w > k 0 , we can define for j = 0, . . . , T − w. We observe from (3) that when the RCN is endowed with the ESP, the training procedure is equivalent to learning the map between u[w + j] and the window of historical data u[j; w + j − 1]. The existence of such a map is guaranteed by Takens theorem. In particular, with w in (3) greater than 2m + 1, there exists a smooth map ψ w : . Therefore, training an RCN can also be interpreted as approximating the map ψ w using nonlinear functions ξ (i,w) for i = 1, . . . , N .
2) Multi-step ahead forecasting: Based on the idea of 1step ahead forecast, we can extend the RCN to accomplish τ -step ahead forecast. Specifically, if the training reference is A similar argument to 1-step ahead forecasting holds for τ -step forecasting since the output of the RCN can be configured to approximate ψ τ Hence, in this way, the RCN training can be viewed as approximating the last component of ψ τ w by using N i=1 c i ξ (i,w) . Remark 1. The above explanation on multi-step ahead forecasting based on time-delay embedding provides a much better understanding on why the RCN successfully forecasts a chaotic system as reported in [31], [32]. In particular, for a chaotic system with an attractor, such as the Lorenz is uniquely determined by a window of historical data. Therefore, the training process of the RCN, in this case, can be interpreted as using N i=1 c i ξ (i,w) to approximate the nonlinear mapping G (similar to (3)). This illustrates the ability of the RCN to forecast a wide variety of time-series by simply changing the training referenceỹ[k].

4) Choices of the activation function:
When the activation function σ is Lipschitz continuous with the Lipschitz constant L, the condition A 2 < 1 L is sufficient to ensure the ESP for (1) as a result of Lemma 1. For instance, the commonly used hyperbolic tangent function, tanh, has a Lipschitz constant L = 1 so that we need A 2 < 1 to ensure the ESP for (1). Since the above explanation on how RCNs learn dynamical systems does not post any restriction on the activation function σ, we indeed have much freedom on the choice of σ. In fact, the activation function can be as simple as a linear function.
The trade-off between the RCN performance and the choice of activation function for the RCN is well-documented. Although it is believed, in general, that nonlinear activation functions perform better than linear ones, it is reported in the literature (see [33]- [36]) that this conclusion is indeed based on a case-by-case study. One major advantage of using linear activation functions, as we shall see in Section IV, is that a linear activation function enables a thorough analysis of the RCN using system-theoretic tools, and allows for an explicit design of the RCN, which can eventually be used as a baseline for designing an RCN with a nonlinear activation function, e.g., tanh or sigmoid, as widely used in the literature.

IV. REALIZATION THEORY FOR THE RCN DESIGN
In this section, we present a realization-theoretic framework for systematic design of RCNs. We illustrate the main idea and conduct the analysis for RCNs with linear activation functions. Specifically, we first show that for given inputoutput sequences u[0; T − 1] ⊂ R p andỹ[1; T ] ⊂ R q , there exists an RCN that approximates the input-output relation in the data. Then, we provide a systematic scheme to prune the size of RCNs while maintaining the same training error. At the end of this section, we illustrate how these results can be carried over to the design of RCNs with general nonlinear activation.

A. Realization of RCNs with linear activation function
We begin by introducing the notion of realization from systems theory [37] and defining it in the context of the RCN, which will form the basis of the proposed RCN design framework.

Definition 1 (Realization of linear systems). Given two se-
Furthermore, if ǫ = 0, then such a realization is called an N-dimensional perfect realization.
For simplicity, we will refer to an N -dimensional realization using (Ã,B,C) N , and this will denote the linear dynamical system in (4). To facilitate the delineation between two realizations, we introduce Markov parameters, equivalent and irreducible realizations as follows.
Definition 2 (Markov parameter). The k th Markov parameter of a realization (Ã,B,C) N is a matrix of real numbers γ k ∈ R q×p defined by γ k =CÃ kB .  [37]. Since in this work, we consider the reduction of ǫ-error realizations of RCNs, Definition 4 is named as 'irreducible realization' to avoid ambiguity.

Definition 3 (Equivalent realizations). Two realizations
Next, we will describe the RCN training procedure through the use of a realization (Ã,B,C) N , and then establish the realization framework tailored for the RCN that explicitly accounts for the ESP, an important and necessary property for the functioning of the RCN.
Consider the RCN as given in (1) with a linear activation function, e.g., σ is the identity function, given by where I is the identity matrix of appropriate dimension. Let A := (1 − α)I + αA ∈ R N ×N andB = αB ∈ R N ×p , then the RCN dynamics can be expressed as Now. we introduce the connection between the RCN training procedure and the notion of realization theory. Specifically, the first step of RCN training is to fix the dimension N , the leakage rate α, and the randomly generated matricesÃ andB. Then, training the output layer of the RCN constructed using A andB is equivalent to findingC * such that Let us consider an N 1 -dimensional ǫ-error realization, denoted (Ã 1 ,B 1 ,C 1 ) N1 , and suppose (Ã 2 ,B 2 ,C 2 ) N2 is the irreducible equivalent realization of (Ã 1 , which implies that if we train the RCN constructed using the N 2 -dimensional matricesÃ 2 ∈ R N2×N2 andB 2 ∈ R N2×p , the training error will be bounded above by ǫ. Therefore, finding the irreducible equivalent realization to a given RCN enables quantifying a smaller size of the RCN that provides a desired training accuracy.
To adopt this realization-theoretic idea for the RCN design, an additional constraint has to be imposed on the matrix A in order for the RCN to be equipped with the ESP. To achieve this, we propose the notion of α-stable realization.
For instance, we know, by Lemma 1, that an RCN in (1) using a linear activation function σ with the Lipschitz constant L = 1 possesses the ESP when A 2 < 1. Therefore, in this case, for any α ∈ ( 1 2 , 1), having Ã 2 < 2α − 1 in (5) is sufficient to guarantee the ESP. This is because As a consequence, finding a (2α − 1)-stable realization (Ã,B,C) N will ensure that the corresponding RCN possesses the ESP. Therefore, in the remainder of this section, we will consider a fixed leakage rate α ∈ ( 1 2 , 1), so that 2α − 1 ∈ (0, 1). For simplicity, we use 'stable realization' in place of '(2α − 1)stable realization'.

B. The irreducible stable realization
Theoretically, if there exists one stable realization that achieves a desired training error for a given input-output sequence, one can construct many different realizations with the same training error. A consequent question of paramount practical importance to ask is how to prune the size of a realization as much as possible while maintaining the training error tolerance of the given stable realization. The answer to this question is pertinent to the concept of fundamental properties of a control system.

respectively.
Controllability and observability properties then lead to the characterization of a irreducible realization (see [38]). Note that Lemma 2 poses no constraints on the matrix norm of the realization, and thus an irreducible realization may be unstable. As a result, modifications have to be made in order to construct a stable irreducible realization resulting in an RCN with ESP. In the following, we develop a systematic scheme to construct an equivalent stable realization of RCN with reduced size. . Note that by construction V is an orthonormal matrix so that V −1 = V ⊺ . Therefore, we have Note that each column ofÃV 1 lies in R(V 1 ). Since columns of V 2 are in the orthogonal complement of R(V 1 ) by construction, it holds that V ⊺ 2Ã V 1 = 0 . Therefore, V −1Ã V can be re-written as Moreover, since every column ofB lies in R(W N ), we have V ⊺ 2B = 0 by the construction of V 2 so that Now we denoteCV = [C 1 , C 2 ], where C 1 consists of the firstN columns ofCV , and C 2 is formed by the remaining N −N columns. Then, it can be observed that holds for all k = 0, 1, . . .. Furthermore, due to the structure provided by (7) and (8), it can be verified that which implies that (A 11 , B 1 , C 1 )N is anN -dimensional realization that is equivalent to (Ã,B,C) N . Now it remains to show that (A 11 , B 1 , C 1 )N is α-stable, given that (Ã,B,C) N is α-stable. By the definition of matrix 2-norm, we have Let y = V 1 x, then y is a vector in R N satisfying Hence, (9) can be bounded by y ⊺Ã y ≤ Ã 2 < α, which concludes the proof.
From a dual perspective, we also have the following lemma regarding the observability matrix.  Proof. We prove the theorem by proving the contraposition, i.e., (Ã,B,C) N is not irreducible if and only if rank W N G N < N .
On the one hand, if (Ã,B,C) N is not irreducible, then by Lemma 2, either rank W N < N or rank G N < N . Therefore, On the other hand, if rank(W N G N ) < N , then, by Sylvester's rank inequality, it holds that which implies that rank W N + rank G N < 2N and that either rank (W N ) < N or rank (G N ) < N . Hence, from Lemma 3 (or Lemma 4), it holds that (Ã,B,C) N is not irreducible.
As a consequence of Lemma 3, Lemma 4, and Theorem 2, the procedure for finding the irreducible stable realization that is equivalent to a given stable realization (Ã,B,C) N is described in Algorithm 1, where Orth(A) returns an orthonormal basis of R(A), and dimÃ returns the dimension of the matrixÃ.
Remark 3. It is worthwhile to mention that as proved in [27], a linear RCN attains maximal memory capacity when its weight matrices (A, B) constitutes a full-rank controllability matrix. As a consequence, Algorithm 1 not only returns an irreducible linear realization of RCN, but also provides an RCN that reaches maximum memory capacity.

Remark 4. Each iteration in Algorithm 1 consists of computing the eigen-decomposition of the controllability or the observability matrix, which has a time-complexity of O(N 3 ).
In the worst case, Algorithm 1 may take N iterations to terminate, which results in a total time-complexity of O(N 4 ). Nevertheless, in all our numerical experiments, we observe that the number of iterations for Algorithm 1 to terminate is of order much smaller than N , i.e., around 10 to 20 iterations for the cases when N = 500, 1000 or even N = 2000, so that we empirically expect that the average time-complexity of Algorithm 1 is O(N 3 ). On the other hand, model selection procedures for designing an initial learning model for the given data typically involves evaluating the performance of models with various hyper-parameters and choosing the model that yields the best result [20], [39]. One of the main features of the proposed approach is that the results of Theorem 2 can be used to evaluate the irreducibility of the RCN and Algorithm 1 can be used to prune the RCN size, irrespective of how the initial RCN model is selected.

C. Further implications of irreducible stable realizations
The development of irreducible realization in the previous section has a lot to offer for designing linear and nonlinear RCNs, as well as understanding the underlying dynamics in the training dataset. We start by explaining how the size of the irreducible realization is related to Takens embedding, which provides a criterion to characterize the complexity of the underlying dynamics determined by u[0; T −1] andỹ[1; T ].
From the theory of linear dynamical system (see [40]), it is a known fact that for any N -dimensional realization (Ã,B,C) N , there exists an invertible matrix P ∈ R N ×N such that (P −1Ã P, P −1B ,CP ) N is in the observable canonical form given by where θ j , j = 0, 1, . . . , N − 1 are arbitrary constants and γ j , j = 0, 1, · · · , N − 1 are the first N Markov parameters of the realization (Ã,B,C) N . It is not hard to verify that for any invertible matrix P ∈ R N ×N , we have (P −1Ã P, P −1B ,CP ) N to be equivalent to (Ã,B,C) N . Therefore, without loss of generality, we assume that (Ã,B,C) N is in the observable canonical form as in (10). In this case, the dynamics of each component of x[k] is given by . .
. Therefore, we havẽ When the RCN possesses the ESP, the effect of x N [k + 1] oñ y[k+N ] is negligible, so that from (11),ỹ[k+N ] is determined by u[k], . . . , u[k + N − 1], which is the input data history of size N . From Takens theorem, the time evolution on a compact manifold of dimension m can be represented by any timedelay embedding longer than 2m + 1. Therefore, if a linear realization describes the underlying dynamics determined by u[0; T − 1] andỹ[1; T ] perfectly, the dimension of such realization satisfies N > 2m + 1 by Takens theorem. On the other hand, if a linear realization of size N attains a training error of ǫ, then it implies that there exists a dynamical system evolving on a manifold of dimension 1 2 (N − 1) that approximately represents the underlying dynamics up to ǫerror. This analysis provides a bound on the size of an RCN representing the dynamics of the underlying dynamical system generating the given input-output data sequences.
In addition, it is worth mentioning that since we are using linear dynamics to approximate the map ψ 2m+1 in Takens embedding, the bound on m mentioned above can be improved through the use of a nonlinear realization. It is intuitive to argue that using a nonlinear realization (e.g., with tanh or sigmoid as the activation function) to design an RCN may result in better approximation compared to a linear realization of the same dimension. However, the explicit solution of RCN dynamics with a nonlinear activation function is in general unavailable. In this case, the linear realization of the RCN offers a guideline towards designing a nonlinear RCN realization. A reasonable design approach for the RCN with a nonlinear activation function is to first design a linear realization, say (Ã,B,C) N , for the given input-output data sequences; then construct an RCN with nonlinear activation function using the matrices A, and B as in (1) through and train the readout layer again.

V. NUMERICAL EXPERIMENTS
In this section, we present several numerical examples to illustrate the developed tractable realization-theoretic approach with training error guarantees, for which the irreducible size of RCNs with respect to desired training errors can be explicitly quantified. Based on linear stable realizations, we further design nonlinear RCNs with the canonical tanh or sigmoid activation functions that achieve desired training error performance.

A. The irreducible realization of time-delay systems
In this example, we design an RCN using the proposed approach for forecasting a time-delay system to elucidate the intimate connection between the irreducible linear realization and Takens embedding.
We first introduce how the training data was generated and how the RCN was trained. In this experiment, we varied τ from 1 to 50. For each τ , we randomly generated a N 0 -dimensional 0.001-error realization with N 0 > τ and computed its irreducible realization using Algorithm 1. Then, using the irreducible realization, we forecast the sequence τ -step ahead for 2000 time steps. When training RCNs, we fixed training length as T = 1000, and leakage rate as α = 0.9. Figure 1a demonstrates the size of irreducible linear realization N min with respect to τ with N 0 selected as N 0 = 4τ . Each point in the figure is the average plus/minus the standard deviation of 10 independent experiments under the same setup. Figure 1a shows that using Algorithm 1, we can trim a large RCN (dashed green line) into much smaller size (solid blue line) with the same performance. This result is consistent with our analysis in Section IV-C that the minimum size of an RCN realization can be used as a criterion to quantify the length of time-delay embedding (dashed red line) associated with the training dataset. Figure 1b provides the averaged mean squared errors (MSE) of the irreducible linear realization with respect to τ .

B. Time-evolution forecast for chaotic systems
In this example, we show the use of an RCN to learn the temporal evolution of a chaotic system, and analyze how the network configuration affects the performance of the designed RCN. Specifically, we consider the Rössler system, given by where x, y, z are state variables and a, b, c are constant parameters selected as a = 0.5, b = 2.0, and c = 4.0. Initial conditions were selected as x(0) = 0, y(0) = 0, and z(0) = 1.
The Rössler system was simulated for t ∈ [0, 1], with 8000 sampling points collected in this time window. The first 5000 points were used as the training input, and the sample points between 30 and 5030 were used as the reference sequence to train the RCN for 30-step ahead forecasting. Then, we recorded the output of the RCN for another 3000 steps as a forecasting sequence. Figure 2a provides a demonstration of forecasting the x-component of Rössler system for 30-steps ahead.
We first randomly generate two RCNs with all hyperparameters to be the same, but one with tanh activation function, the other one with linear activation function. The hyper-parameters of the RCN were fixed as follows: number of nodes N = 500, leakage rate α = 0.8, length of training data t r = 5000, length of forecasting data t s = 3000, and washout length w = 500. If not specifically mentioned, the matrix A was generated under a uniform distribution on [0, 1] of sparsity 63.2%, and then normalized to have a matrix 2norm equals to 0.9. The output-layer was trained via ridge regression [41] with λ = 10 −8 .
Figures 2b presents the results of forecasting the time-series provided in Figure 2a with a mean-squared-error (MSE) of 1.49 × 10 −3 . Under the same setup as in Figure 2b, we varied the leakage rate α and the activation function of the RCN, and present the corresponding results of MSE versus timedelay (τ ) in Figures 2c and 2d, respectively. From Figure  2c, we observed that α = 0.05 resulted in the best MSE across different cases (of τ ); and from Figure 2d, we observed that the RCN with linear activation function achieved similar performance as the RCNs with a nonlinear activation function when the time-delay τ was small. Therefore, we used the same RCN as in Figure 2b, but changed the leakage rate to α = 0.05 and the activation function into linear function to forecast the same time-series. The corresponding results are presented in Figure 2e with an MSE of 2.97 × 10 −2 .
We further applied Algorithm 1 on the above RCN with linear activation function, yielding a irreducible linear realization of size 31 with a MSE of 2.54 × 10 −2 . Figure 2f presents the result of using the irreducible RCN to forecast the same time-series generated by Rössler system for 30-steps ahead. Based on our empirical analysis in Section IV-C, the results in Figure 2f implies that the underlying dynamics of the given input-output sequences can be well-approximated by a linear dynamics with a time-delay of no longer than 31 steps, which is evident by the experiment setups.

C. Analysis of linear and nonlinear activation functions
In this part, we further investigate the difference in performance between linear and nonlinear RCNs. The results in this section support our idea of using the linear realization of the RCN to help design a nonlinear RCN with a guaranteed training error, as mentioned in Section IV-C. Table I presents the results of forecasting time-delay systems using linear and nonlinear RCNs. We varied the size of the RCN from 50 to 500. For each N , we generated a time-delay system with N delay steps (as in Section V-A) and simulated 2000 independent experiments of forecasting (N − 1)steps ahead. The training length was fixed as t r = 5000 and the forecast length was fixed as t s = 2000. To make a fair comparison, in each experiment, we generated three RCNs of size N , one with linear activation function, one with tanh activation function, and another with sigmoid activation function, using the same randomly generated matricesÃ and B. Other hyper-parameters of the three RCNs were set to be the same as in Section V-A. The average and standard deviation of training and forecast MSE are reported in Table  I. As we observe from the last column in Table I, the nonlinear RCNs outperform the linear RCN in terms of the training MSE in most cases. Table II reports the results of forecasting the Rössler system using RCNs with different activation function. Using the same dataset as the one in V-B, we compared the performance of linear and nonlinear RCN at different scales. Similar to the previous table, we varied the size of the RCN from 50 to 500 and conducted 2000 independent experiments to forecast the Rössler system for 10-steps ahead. The training length was fixed as t r = 5000 and the forecast length was fixed as t s = 2000. In each experiment, we generated three RCNs of size N , one with linear activation function, one with tanh activation function, and another with sigmoid activation function, using the same randomly generated matricesÃ andB. Other hyperparameters of the three RCNs were set to be the same as in Section V-B. The average and standard deviation of training and forecast MSE are reported in Table II. As we observe from the last column in Table II, in this experiment, the training MSE of nonlinear RCNs was always smaller than that of the linear RCNs.

VI. CONCLUSIONS
In this paper, we provided a detailed analysis and a holistic description of the training procedure and the operation of the RCN. With the help of Takens embedding theorem, we derived the delay embedding map, which an RCN potentially learns during the training process from the given input-output data. This provided insights into the role that the linear activation function and other hyper-parameters play in the design and working of RCNs in applications such as forecasting a timeseries. Furthermore, based on the notions of linear realization theory, we provided a systematic approach to trim RCNs with guaranteed training accuracy. In this context, we introduced the idea of α-stable realizations for designing stable RCNs that achieve the desired training objective with reduced size, and established a tractable design algorithm to synthesize RCNs with nonlinear activation. The numerical experiments on forecasting time-delay systems and the Rössler system were used to substantiate our proposed design approach for interpretable RCNs. We observed from the experiments that the nonlinear RCNs with both reduced size and guaranteed training accuracy can be attained based on the minimum realization for linear RCNs. These results suggested that the proposed approach offers an informed and interpretable design methodology to devise nonlinear RCNs for a given dataset to decode the underlying dynamics in the data. N Training MSE (mean ± std) Forecasting MSE (mean ± std) P(ǫ linear > max(ǫ tanh , ǫ sigmoid )) linear (×10 −5 )