{'fraction_non_alphanumeric': 0.040116201711979424, 'fraction_numerical': 0.01215239277952721, 'mean_word_length': 4.4667232391805065, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Explaining the Surprising Success of Reservoir Computing Forecaster of Chaos? The Universal Machine Learning Dynamical System with Contrasts to VAR and DMD

Machine learning has become a widely popular and successful paradigm, including in data-driven science and engineering. A major application problem is the forecasting of future states from a complex dynamical, given a cache of data representing time ordered observation of states from the system. Artificial neural networks (ANN) have evolved as a clear leader amongst many machine learning approaches, and recurrent neural networks (RNN) are especially well suited due to an aspect of memory associated with the concept, even if the major step of training the RNN to data which typically involves backpropagation and optimization becomes computationally especially intensive. In this setting, the echo state networks (ESN) or reservoir computer (RC) have emerged for their simplicity since they are a special case of an RNN, which are not fully trained to the data. Instead only the readout weights are trained but read-in weights and internal weights are simply selected randomly. However, the read-out weights can be trained by a simple least squares step and so it is simple and efficient to train. There remains an unanswered question as to why and how an RC works at all, with randomly selected weights. To this end this work analyzes a further simplified RC. Specifically we are then able to connect the RC to the well developed time-series literature on vector autoregressive averages (VAR) that includes theorems on representability through the WOLD theorem. Further we can associate this paradigm with the now widely popular dynamic mode decomposition (DMD), and thus these three are in a sense different faces of the same thing. Our simplification is not presented for sake of tuning or improving an RC, but rather for sake of analysis the surprise being not that it doesn't work better but that such random methods work at all.


Introduction
Artificial neural networks (ANN) have emerged as a core and powerful technology in machine learning [22,46,56,57,59] that is well suited for the supervised learning in data-driven science and engineering, specifically including for forecasting problems in complex dynamical systems [25,42,36,47,13,44,34]. However, the most straight forward feedforward ANN with back propagation for training concepts can be extremely expensive to optimize to the data, even considering important recent innovations such as stochastic gradient descent or hardware break throughs such as GPU-based processing. Recurrent neural network concepts, (RNN) are especially suitable for temporal data from a dynamical system [25,41,5,6,20,68], as they naturally embed temporal information, and especially the long short term memory (LSTM) approach demonstrate excellent fidelity, [82,35,81,20,18,86], but these are especially expensive to fully train, [63].
The reservoir computing (RC) [38,52,80] and the closely related echo state network (ESN) [37,51] and liquid state machine (LSM) [54,31] have emerged as a special variant of an RNN, where only the output layer is trained rather than the entire network of weights. As such, this requires only a simple and efficient least squares estimation, rather than the more expensive full nonlinear optimization associated with a fully training an RNN. Nonetheless, and perhaps a most surprising outcome is that despite this gross simplification, the forecasting capability can still be competitive even for chaotic or spatiotemporally complex problems [82,64,88,50,16,17,27]. Specifically, an RC thrives when a full state observation is available, while fuller and more expensive variants of RNN, especially the LSTM would considered higher perming, especially when only a reduced variable set is available [81,82]. Still, the RC are popular, surely because of their simplicity to train, and perhaps in part because of their undeniable even if surprising fidelity.
The purpose of this work is to bridge the theoretical gap, to offer at least a partial explanation as to how an RC can be such a successful and general universal dynamical system for forecasting such a wide array of systems despite randomly "trained" read-in and inner weights. The purpose of this paper is not specifically to build a new method, or improve the current method, but to explain what is perhaps a surprising that the RC method works at all. In so doing, we challenge the concept with a simplified linear activation function version for sake that this allows our simplified analysis throughout, and even if this version has reduced fidelity, we show it does still have theoretic reasons it still works. There have been few significant explanations as to why it works so well, but notably [29,15,21]. Usually instead we find in the literature a collection of descriptions as to how to choose random networks as the inner layers, regarding sparsity [75,30], or regarding design of the spectral radius for linear stability [17,39,30] and the echo property, [26].
In the spirit of still incomplete theoretical basis as to the underlying success of RC, we allow a simplified version of RC with linear activation functions for which we are able to more fully identify the inner workings of how the RC can be a universal forecasting machine, even for time-series from complex and chaotic dynamical systems. We show that by this simplification the RC still works, albeit with reduced performance quality, but nonetheless the purpose here being theoretical explanation of how such a simple system of only training the readout is possible. We offer this variant as a theoretical construction. By this interpretation, we will also be able to connect the RC concept to other theoretically more matured theories. Specifically, the theory of autoregression (AR) from time-series analysis and moving averages (MA), and together called ARMA [62,19,12,66,78,74] , are founded on the Wold theorem [84] that we show are directly related to the RC concept. The vector formulation of these [67,53], called vector autogregression (VAR) and vector moving averages (VMA) are also connected by a corresponding Wold theorem. Further, we describe a relationship to the recently highly popular dynamic mode decomposition (DMD) [72,83,43,47,8], which is an empirical formulation of Koopman spectral theory, [8,3,10]. So while we do not offer this simplified RC for performance over other approaches, we hope that this work will serve to shed light on how the simplified RC approach is capable of providing useful time-series forecasts, and likewise as a suggestion as to how the general RC is successful.
This paper is arranged as follows. In Sec. 2 we describe the nature of the data as derived from a stochastic process. In Sec. 3, we review the standard RC concept, and we demonstrate it already with time-series data from a Mackey-Glass differential delay equation. In Sec. 4, is the heart of this paper, where we first present that a linear activiation function allows the RC to be stated as a linear recursion, and therefore fitting just the readout weights can proceed in a manner such that despite random readin and inner weights, there is a well posed problem. Then, in this form we are able to directly relate the linear RC solution to the classical VAR(k) solution. As such, we are then able to enlist statistical time-series forecasting theory for forecasting stochastic processes, so that the Wold theorem that guarantees a VMA can then be translated to a VAR. Furthermore the associated companion form of a VAR(1) usefully states the full vectorized problem. In Sec. 6, we note that the companion form of the VAR(1) is reminiscent of prior work for another famous concept in data-driven dynamical systems which is the time-delay formulation of a DMD-Koopman analysis. Finally in the examples Sec. 7, we present two classical examples, the Mackey-Glass differential delay equation and the Lorenz63 ordinary differential equation, with examples comparing aspects of a full nonlinear RC and the linear variant of an RC.  56), has become a standard example for time-series forecasting, for benchmarking data-driven methods since it is dynamical rich and high-dimensional and therefore challenging. (Top) Time-series, index. (Bottom) Three-dimensional projection in delay coordinates, (x(t), x(t − τ ), x(t − 2 * τ )), τ = 20. A sample of N = 10, 000 data points is chosen as the training data set.

The Data as Sampled From a Stochastic Process
For data-driven forecasting problems, we require data from a process, including from a deterministic or otherwise from a stochastic dynamical system, [11]. A process, stated, is in terms of a collection of random variables, X t on a common probability space, (Ω, B, P ), where Ω is the sample space, B the σ-alebra, and P a corresponding probability measure. T is a "time" index set and commonly it is chosen as either R, or Z or subsets. For sake of discussing finite samples of data, we emphasize maps, which may well be from discretely sampling a flow. A data set from such a process samples x ti of X ti , stated as a time sorted sample, {x i } N i=1 , t 1 < t 2 < ... < t N , using indexing notation, x i := x ti . Uniform timing is also a simplifying assumption, Data derived from a flow, say, may be collected by stroboscopic map, Suppressing the stroboscopic time t, this is a discrete time map F, and likewise other Poincare' maps may be useful for flight between surface of section, and random dynamical systems may also be relevant [73,11]. An underlying principle here is that the data should be "long enough", and likewise a general failing of any data-driven machine learning method for forecasting a stochastic process will tend to do much better in terms of interpolation than extrapolation. Generalizing, to allow for out of sample forecasts will tend to fare much better when the point to be forecasts is close to other observed inputs. Said another way, the quality of results can be brittle, depending as much upon curating a representative data set as the details of the method used to avoid that struggle between fitting between observations and overfitting and too far out of sample.
As a matter of presenting examples, we will highlight two classic problems that remain popular in benchmarking for machine learning in recent literature. These will be, • The Mackey-Glass differential delay equations, Eq. (56), and • The Lorenz63 system, Eq. (57), both of which will be presented in fuller detail in Sec. 7. In Fig. 1 we show early in this presentation for sake of context, a time-series data set of the Mackey-Glass system, from Eq. (56), to stand in as a typical data set. This problem is a useful benchmark, and it is often used as such [55,32,2,60,9,24,85], perhaps because it is a well known chaotic process, but also for sake of dimensional complexities that we recall in Sec. 7.1.

Review of The Traditional RC With Nonlinear Sigmoidal Activation Function
In this section we review the standard and fully nonlinear RC method, by which we mean, including the use of a nonlinear activation function q(s). In this context, q(s) is usually taken to be a sigmoidal function such as the hyperbolic tangent function. However, in the next section we will challenge these steps including simplifying to the identity function, q(s) = s. Assuming the training data, {x i } N i=1 ⊂ R dx , the reservoir computing RNN is stated, The hidden variable r i ∈ R dr is generally taken to be of a much higher dimension d r > d x , by a linear lifting transformation, and W in is a randomly selected matrix d r × d x of weights. See Fig. 2. A is also a linear transformation, as randomly chosen square matrix d r × d r of weights, that should be designed with certain properties such as spectral radius for convergence [17,39], or sparsity, [50,64,82] or otherwise consideration of the "echo-state" property, [15]. Likewise, the readout is by a linear transformation, using a d x × d r matrix of weights W out . However, W out , and only W out , is trained to the data, allowing for forecasts y i given data x i ∈ R dx , which is the major simplify aspect of RC since it can be done by a simple and cheap least squares computation. Finally q : R → R is an "activation" function, using the phrasing from machine learning in the neural network community to mimic the concept of a biological network that fires when a voltage has reached a threshold. Popular choices include q(s) = tanh(s), meaning a componentwise application of the scalar hyperbolic tangent function when s is multivariate. Other activations are popular in general neural network theory, including other sigmoidal functions, and also the ReLu function in certain contexts but not so commonly in RC, [27]. 0 ≤ α ≤ 1 serves to slow down the RC, to moderate stability of the fitting, but we will restrict to α = 0 in this paper as outside the purpose of challenging the concept of explaining how the RC may work in a special case of identity q in which case, nonzero α can be considered as absorbed into the random A; (1 − α)r + αAr = ((1 − α)I + αA)r and since A is chosen randomly, then ((1 − α)I + αA) may be an alternative random selection. Finally, b serves as an offset for activation, that is useful in some contexts, but it is also not relevant for our needs for the same reason we choose α = 0, and we choose b = 0.
What is remarkable about RC is that the usual hard work of optimally developing a full RNN is almost entirely skipped. Instead of learning W in and A optimally fitted to the data, these seemingly very important matrices are simply picked randomly. This is an enormous savings over what would usually be inherently hard to handle since the parameters are composed within the nonlinear activation q and require at least a gradient descent optimization of back-propagation in a high dimensional and likely multi-peaked optimization space. Almost any matrix distribution may plausibly due, but several different recipes are suggested. We say "recipe" rather than algorithm since these are descriptions of successful observations in practice, rather than a product of mathematical theory that is still not complete. Here, we choose the entries of A uniformly, A i,j ∼ U (−β, β), with β to scale the spectral radius, but other choices are common, notably for sparsity. The read in matrix is also chosen uniformly randomly, W in i,j ∼ U (0, γ), with γ > 0 chosen to scale the inner variables r. Reservoir trained across the data set, and 500 samples are shown for clarity where we see the error is sufficiently small that the one-time-step forecasts and the true data are almost the same so that the plot is indistinguishable (both shown, but curves overlay). Regularity is chosen to be, λ = 1.e − 6. (Bottom) Some randomly selected 7 of the (usually hidden) d r = 500 activation functions illustrate the general appearance. Contrast to forecasting into the future as shown in Fig. 4, and linear method in Fig. 5.
The crucial aspect of the simplification that makes reservoir computing so easy and computationally efficient, is that training to the output becomes just a linear process. The cheap and simple least squares solution is easily handled directly by matrix computations. Let, Notation here is standard that · F denotes the Frobenius-norm of the matrix, which is the least squares equivalent of the least squares matrix parameter estimation problem. The data are the forecasts to X to be optimized in least squares by W out , processed through the RC, While k = 1 is allowable, here for theoretical development in subsequent sections, we allow for larger k ≥ 1, describing memory. In practice a ridge regression (Tikhonov regularization with least squares regularity, [27,64,28,7]) is used to mitigate overfitting, the solution of which may be written formally, Notation includes · T is the matrix transpose, I is the identity matrix, and the choice of regularity parameter is λ ≥ 0. We will write a regularized pseudo-inverse with the notation, In Appendix 11 we review the matrix theory as to how to form regularized pseudo-inverses such as R † λ by a regularized singular value decomposition (SVD) in terms of regularized singular values such as σ i /(σ 2 i + λ) obtained from the singular values σ i from the SVD of R.
In Fig. 3, we show an example of an RC machine obtained from data obtained from the Mackey-Glass differential delay equations, Eq. (56). We see fitting for N = 10, 000 data points x(t), d x = 1, regularizing parameter λ = 1.0 × 10 −8 , and fitting for constant time offset. Fit and true data are shown to be so close that in fact the blue fit curve hides the red true data curve. Also shown are several (7) of the d r = 500 hidden variables r(t). The fit matrix A is randomly chosen with entries from a uniform distribution, and then scaled so that the spectral radius ρ(A) = 1. The random random matrix W in is also chosen uniformly, scaled so that x values lead to r in [−0.6, 0.6]. In Fig. 4, the trained RC are used to forecast into the future. We see small errors grow in scale, as illustrated by the bottom error curve. Results from an RC forecasting for the Lorenz63 system are presented in Sec. 7.2, and notably the forecasting quality degrades more quickly in part due to known large Lyapunov instability of that system.
What is amazing is that despite that RC may seemingly be a gross oversimplification of the RNN concept, it still seems to work quite well. Also from experience, it is generally stable in that it is somewhat insensitive to the parameters and hyperparameters of the fitting process, even if the level of quality does depend on these. Furthermore, once it starts to make larger errors the kind of dynamics it produces are still plausible alternative wave forms of the process. Nonetheless there are some parameter choices to make, notably, d r > d x must be "large enough," but how big is not well understood. Furthermore, the nature of the underlying distribution of matrices W in and A is not fully understood. We hope to contribute some general perspective as to why an RC may work at all, even though our goal here will not be to handle aspects of improving performance.

RC With A Fully Linear Activation, q(s) = s, Yields a VAR(k)
Now we attempt to challenge a central typical assumption of the RC method. Instead of choosing the activation function to be a sigmoid function, instead, we use the identity function, q(x) = x. With this assumption, we can show that the resulting linear RC machine is equivalent to a vector autoregressive process (VAR) [67,33], which is extremely popular and successful in the timeseries forecasting field, particularly in econometrics [1]. With this simplification, we find that not only can the linear RC still make useful forecasts, but we are able to connect the RC concept to this well established theory associated with VAR time-series analysis, notably the existence of representation WOLD theorem, [84,62]. However, while this gives some explanation as to why a standard nonlinear RC may work despite the seemingly oversimplification of a full RNN, we show that that the linear RC does still performs and furthermore, now with theoretical underpinnings, even if the full nonlinear RC may still perform better. So it is for the theoretical connections that we make this simplification, rather than a suggestion that it may be a new or simpler method.
Before proceeding with a discussion of q(s) = s, notice that r is related to the scale of the read-in matrix, W in . Proceed by initializing the process, by Eq. (5), Consider that since W in is randomly chosen, and we choose uniformly W in ∼ U (0, γ), then the parameter γ > 0 moderates the subsequent scale of terms u i and then r i . See for example Fig. 3, where the native data x from the Mackey-Glass system is translated to scaled internal variables. Recall the power series of the nonlinear activation function, Clearly for s << 1, then q(s) ∼ s even if chosen as a sigmoid, and the choice of read-in scale could be designed to put us in this regime as long as A is designed to keep us in this regime. However, in the following we proceed to study the consequences of stating the activation exactly as the identity, q(s) = s.
With this assumption, the first several iterations follow from Eq. (4) and Eq. (11) as a forward propagation, for which we explicitly observe the following recursion.
using notation, A 0 = I, the identity matrix. Since the readout of this process is by Eq. (4), y i = W out r i , then we may rewrite the final equation, Eq. (18), by left multiplying by W out .
with notation, Each of these coefficients a j are d x × d x matrices. This follows simply by Eq. (20), collecting products between d x × d r to d r × d r and then d r × d x matrices and notation By Eq. (19), a linear RC yields a classical VAR(k), (a vector autoregression model of k-delays) that in a general form is [67], In this writing, c allows for a general offset term, a d x × 1 vector that here we do not pursue. The ξ k+1 is underlying "noise" of the stochastic process which is part of the stability theory we review in the next section, must be assumed to come from a covariance stationary process. This relationship between an RC and a VAR(k) allows us to relate to the corresponding theoretical discussions of relevant alternative forms and stability and convergence from the stochastic process time-series literature, that we will also expand upon in the next section.
Considering the complete data set of vector time-series, Restating this as a single linear equation, Again, remembering that x i are d x × 1 vectors and that a i are d Notice that we have stated the target values, which are chosen from the data in practice, y k+1 = x k+1 . Formally, minimizing in least squares, with regularization, with Y being the target output of the right hand side of Eq. (22) by best fitted matrix a * . The solution of this regularized least squares problem may be written in its matrix form, where the symbol † refers to the Penrose pseudo-inverse, with notation described in detail in Eqs. (62)-(63), when formulating the "ridge" Tikhonov regularized pseudo-inverse X † λ .

Decomposing the VAR(k) Solution Explicitly Relates to RC
Now, we will further decompose the derived VAR(k) coefficients found in Eq. (25), to emphasize the training of just the output matrix W out of an associated RC, in terms of randomly pre-choosing A and W in . Referring to Eqs.
with the matrix defined, This is a combination of exponents of the random d r × d r matrix A, and the random d r × d x matrix W in , and so it is itself a kd r × d x random matrix. Interestingly, considering just one column at a time of the W in l , l = 1, 2.., d r , A k−1 can be understood as a collection of columns from a Krylov space and this entire process can be discussed as an Arnoldi-iteration, which is something we will explore further in Section 6.
Consider that the least squares objective Eq. (29) can be expanded to split, to emphasize that since if we pre-choose A and W in , then only the readout matrix v is a free parameter, Optimizing for v yields, Comparing this equation with Eq. (25), defining a, we see (XX † ) formally appears in both expressions. Only the associative property of matrix multiplication is needed to emphasize the role of A. More importantly, this expression Eq. (30) for W out is written so as to emphasize that the reservoir computing process is designed with A and X. Combined through the iteration, as (AX) is the data that results from Eq. (17), This is written naturally, R = (AX).
by the simple way Eq. (30) uses a matrix identity of pseudo-inverses, [28], Associativity emphasizes that since A is deterministically defined, once A and W in are chosen, and separately from the data X, then the fitting of only the parameters of W out are sufficient. If we want the VAR(k) parameters, we could either ignore the prior knowledge of choice of A and W in , and compute a directly from Eq. (29), or from Eq. (30), defining, We summarize that these manipulations concluding with Eq. (34) serve directly as the connection between the RC fitted readout and the coefficient matrices of a VAR(k). The roles of pre-choosing A and W in relate directly to W out coefficients, or indirectly to the fitted data. Concluding this section with the an example, we simplify the nonlinear RC of the Mackey-Glass data from Figs. 1, 3, to a purely linear RC fit shown in Fig. 5 which clearly is not as well performing but it does still make some forecast into the future. Further discussion of this example and also likewise a Lorenz63 example in Sec. 7.

VAR(k) Theory Suggests Convergence with k
Since the VAR(k) model of vector autoregression appears naturally in our discussion from the simplified activation function q(x) = x, as summarized by Eqs. (19), and (21), we now recall some of the classical underlying theory from the statistical time-series analysis literature [84,67] that describes sufficient conditions under which we expect existence of a VAR(k) representation.
The Wold theorem plays a central role in time-series analysis as it describes existence of a vector moving average (VMA) model representation, which then under further assumptions for invertibility, is equivalent to a VAR. Assumptions require a stationary process as a sum of two components: 1) a stochastic component consisting of "linear" combinations of lags from a white noise process, and 2) a deterministic component that is uncorrelated with the stochastic component. First we recall definitions. A d-dimensional stochastic process ξ t of zero mean, E(ξ t ) = 0, is derived from a white noise stochastic process, written with zero mean ξ t = [ξ 1,t , ξ 2,t , ..., ξ d,t ] ∼ W N (0, Ω) if E(ξ t ) = 0 and E(ξ t1 ξ T t2 ) = 0, for t 1 = t 2 , but E(ξ t ξ T t ) = Ω is symmetric positive semi-definite. A stochastic process is covariance stationary if all terms of the sequence have the same mean, and any two terms depend only on their relative positions. That is, E(ξ t ) = E(ξ t ), for all t , and for all t ≥ 0, there exists γ t ∈ R such that, Cov(ξ t , ξ t−t ) = γ t , for all t > t , meaning depending on t − t rather than the t or t . With these definitions, we can state the central theorem of this section that we recall: Theorem 1 (Wold Decomposititon Theorem, [84,67]) A zero mean covariance stationary vector process {x t } admits a representation, where C(L) = ∞ i=0 C i L i is a polynomial delay operator polynomial, the C i are the moving average matrices, and L i (ξ t ) = ξ t−i . The term C(L)ξ is the stochastic part of the decomposition. The µ t term is the deterministic (perfectly predictable) part as a linear combination of the past values of X t . Furthermore, • µ t is a d-dimensional linearly deterministic process.
• For each t, µ t is called the innovation for X t , in that µ t = X t − P (X t |X t−1 , X t−1 , ...). The so-called shock is fundamental.
Clarifying notation of the delay operator polynomial, with an example, let, so if for example, x t ∈ R 2 , For interpretation and definition, consider: • If µ t = 0, then this is called a "regular" process, and therefore there is a purely vector moving average (VMA) representation. If C i = 0 for i > p for some finite p > 0 then it is called a VMA(p) or otherwise it is a VMA(∞) representation.
• If X t is regular then the representation is unique. Now to our point to relate a Wold VMA representation to our discussion following the linear RC where we saw a VAR(k) results Eqs.
Existence of this inverse implies that the Wold implied VMA process has a VAR representation, In practice, when an infinite order vector moving average process, VMA(∞) corresponds to an infinite order vector autoregressive process, VAR(∞), then recursion of expanding Eq. (39) and matching term by term yields, Though, a VAR representation may be found from a VMA through several methods, including a method of moments leading to the Walker-Yule equations, [65], or a least squares method in the case of finite presentations. Often, for parsimonious efficiency of presentation, a mixed form of a p-step AR and a q-step MA model might make a suitable approximation, for what is called a ARMA(p,q) model. While not allowing ourselves to be drawn entirely into the detailed theory of econometrics and statistical time-series analysis, pursuing stronger necessary conditions, we wish to point out some already apparent relevant points from the stated special sufficient conditions. Remark 1 Summary statements. If a vector stochastic process satisfies the hypothesis of a Wold theorem, then it: • Can be written either as a VMA or a VAR, when Eq. (40) of C(L) is invertible, Eq.(41).
• Furthermore, in practice a least squares estimate of a VAR(k) may be used for finite k, which relates to an RC by the least squares fit, Eqs. (30), (34).
Finally, we separate from the above technical points, the following fundamental remark to distinguish existence versus uniqueness of a representation, Remark 2 While a stochastic process may have a VMA representation and if through invertibility, a corresponding VAR, which is a linear descriptions of the process, it may not taken to be "the" unique physical underlying description since nonlinear descriptions certainly may exist.

Remark 3
The processes that we may be interested in, such as those derived from Eq. (1), may describe the evolution of a (chaotic) dynamical system and these may allow a representation, Eq. (3), [11,76,45]. However, in many of these natural examples, the "color" or even the nature of the noise may well not be conforming to the white noise assumption of the Wold theorem 1. Certainly contrasting samples from an invariant measure from a chaotic dynamical system to a white noise process is a well studied [69,40], but still undecided topic. While existence of the VMA and corresponding VAR representation by referring to the Wold theorem does depend on that hypothesis, nonetheless, successful constructive fitting of a VAR(k) by regression, even if implicitly through an RC, seems to proceed successfully in practice in a wide array of examples.
With this last remark, we admit that while the details of the rigor guaranteeing existence may in practice break down, due to inability to check all hypothesis, as often such gaps occur between mathematics, applied mathematics, and practice as related to real world data, we feel that the concept is still highly instructive as underlying explanation, despite strong sufficient assumptions used to extend a rigorous theory.

Stability of the VAR(k) and relationship to a VAR(1)
To discuss stability, we recall [67] the fact that a VAR(k), Eqs. (19), (21), can be stated as a VAR(1) in terms of "stacked" (delayed) variables, called the companion system. This idea is familiar in dynamical systems as we see it is related to stating time-delay variables and the Taken's embedding theorem [77,61,71,58,9,87]. Define, where since a i are each d x × d x matrices,then A is kd x × kd x , and X k is kd x × 1. For discussion in the next section, it will be convenient to consider for contrast to Eq. (50), a matrix of all the data, are kd x × (N − k − 1). Notice that the data in X, also from Eq. (22).
It follows that analysis of stability of a VAR(1) sufficiently describes the stability of a VAR(k). If there is even a small offset c, whether by a bias or imperfection of fit, then follows the recursion, This relates the VAR(1) back to a VMA(l) form. Clearly, even a small constant disturbance C is successively influenced by the (delay) matrix A. In the limit l → ∞, recall the geometric series of matrices, converges if the spectral radius is strictly contained in the complex unit disc, Equivalently, a general VAR(k), Eq. (42), is stable if and only if a characteristic polynomial, has all its roots outside the unit disc. Under this stability assumption we conclude that, which relates a VAR(1) form to a Wold form through a VMA(∞). Since by Eq. (20), each matrix a j = W out A j−1 W in , then the magnitude of entries in the matrix A and the read in matrix W in each moderate the magnitudes of entries of a j . So considerations by the Gershgorin disc theorem, [28] relates these magnitudes to the magnitudes of z. Generally sparsity of A, magnitude of the spectrum of A and magnitudes of W in can be reduced for stability, and to moderate the "memory" associated with converges with k, and these magnitudes were already discussed for sake of a regime where the usual sigmoidal q would be close to the identity.
6 Is There a Connection to DMD-Koopman?
To briefly answer the question titling this section, the answer is yes, there is a connection between VAR and DMD, and so to RC. The more nuanced answer is that the connection is not complete. Throughout the discussion so far, a specialized version of an RC using an identity activation function, yields a linear process that is shown to relate to a VAR that is also a linear process. In this section we ask if it also relates to the dynamic mode decomposition, DMD [72,70,43,83], a concept that is also premised on a linear process model as a finite estimation of the infinite dimensional linear action of the Koopman operator on a function space of observables [3]. In Koopman theory, instead of describing the evolution and geometry of orbits in the phase space, the transfer operator methods generally describe evolution of functions whose domain is the phase space [45,11]. Recently this approach has excited a huge trend in applied dynamical systems, with many excellent research papers [83,43,47,8], review papers [3,14], and books [43] toward theory, numerical implementation and scientific application practice. Our focus here will remain narrow, the goal being to simply identify a connection to the RC and its related VAR, as discussed above. A primary purpose of DMD methods are for modal analysis of the system to describe coherent and typical behaviors, but it also can be used for forecasting, and for this sake the analogy is drawn here. For direct comparison, first allow some minor manipulations to relate the VAR(k), Eq. (42) and Eq. (22), to a typical DMD form. A time-delay version of a linear evolution is a special case of an exact DMD written as follows, with notation used as above, where X, and X are the kd x × (N − k − 1) data matrices in Eq. (44), and K is a kd x × kd x DMD matrix approximating the action of the infinite-dimensional Koopman operator. Abusing notation slightly, the least squares problem, has the solution, which is called the "exact DMD" solution. While there are many variants of DMD, this one called exact DMD is popular for its simplicity of implementation while still useful for interpreting the system in terms of modal behaviors. Contrasting K derived by exact DMD, Eq. (50), versus A for the VAR(1) form described in Eqs. (43)-(44) reveals clear similarities since each states a linear relationship between the same data, X = KX, versus X = AX, but these are ill-posed equations and the A need not be the same as K. Closer inspection reveals that Eq. (53) allows freedom for best least squares fit considering the entire matrix K, and so differences relative to Eqs. (22), (25). Whereas, only the first k rows of A are free parameters in the regression; the subsequent rows of A are sparsely patterned with either zero's or the identity matrix, Eq. (43).
A similar, but not identical, structural difference appears when contrasting the SVD based exact DMD to the original DMD method of Schmidt [72] and also Rowley and Mezic, [70] which is an Arnoldi-like version of DMD in terms of iterations in a Krylov space [4,79]. Reviewing that Arnoldi-version of DMD, using the notation of [70], observations x k ∈ R d are assumed (fitted) to be from a linear process, but also by considering the iterations are to be fitted in the Krylov space, assuming that x m ∈ Kry m (x 0 ) = span{x 0 , Ax 0 , ..., . . .; c m−1 ] is the vector of coefficients. Then a key and clever observation was to rewrite this in terms of a companion matrix, So that results, From there, exploiting the theme of Arnoldi methods, the eigenvalues of C are related as a subset of the eigenvalues of A and with a direct linear relationship between eigenvectors of C and A, Ritz vectors, and the unknown coefficients c of C can be computed by a least squares procedure. Keeping in mind that power iterations as one does in Krylov spaces emphasize just the dominant direction, the Arnoldi methods take care to orthogonalize at each step in the algorithm for stabilization an otherwise unstable search for large sparse matrices, and these make deliberate use of QR decompositions. Our interest here is only to point out analogies between A from reservoir computing and VAR(1) and K, rather than to continue toward discussion of modal analysis as one does in DMD analysis. Summarizing, the analogy we see that the companion matrix C in Eq. (54) reminds us of the companion matrix A in Eq. (43). However the most significant difference is that while c i are scalars, that a i are d x × d x matrices.

Examples
The example Figs. 1, 3-5 already threaded in the above presentation of methods, were in terms of the Mackey-Glass differential delay equation system, which we now recall. Then in the subsequent, we will show similar figures highlighting the concepts in a different system, the famous Lorenz63 ODEs.

Example 1: Mackey-Glass Differential Delay Equation
The Mackey-Glass differential delay equation, [55], has become a now classic standard example in time-series analysis [23,48] of a high (infinite) dimensional dynamical system with a low-dimensional attractor, which we have used as a benchmark in our own previous work, [9] and including recently in presentations for machine learning [9]. The problem is physiologically relevant for describing dynamic diseases. A differential delay equations can be described as infinite dimensional dynamical systems, a concept that is more easily understandable in terms of the notion that an initial condition state advances not just a finite dimensional vector, but rather an entire interval [t 0 , t 0 + t d ] of initial values of x(t) are required. However, the MG equations have a nice property for the its practical use as a benchmark problem, which is that there is essentially an attractor whose fractal dimension varies with respect to the parameters chosen, allowing for a complexity tunable test problem. We have chosen parameters t d = 17, a = 0.2, b = 0.1, c = 10.0 for which if pursing time delay embedding gives an embedding dimension of d = 4. We use integration steps of ∆t = t d /100 throughout. We show time-series in Fig. 1, a standard nonlinear RC forecast of the system in Fig. 3, and the linear RC/VAR forecast of the system in Figs. 4-5.

Example 2: Lorenz63
The Lorenz63 system [49] is the three coupled ordinary differential equations: While these Lorenz equations may have been originally posed as time varying Fourier coefficients for describing a partial differential equation system describing convection rolls of heated fluid in an atmospheric system, they have become been a popular paradigm in the study of chaotic systems, for foundation principles of chaos historically and ongoing, as a simple and familiar benchmark problem and also in the pedagogy of dynamical systems. The chaotic attractor in the phase space {x(t), y(t), z(t)} illustrates a familiar butterfly, but we show a segment of the x(t) time series that will be used as our data set, in Fig. 6 and Fig. 7 . Also   Fig. 6. Quality is not as good, but it still works reasonably, at least for many randomly chosen A.
shown are nonlinear RC, q(x) = tanh(x) activation forecasts in Fig. 6 using the usual nonlinear reservoir computing, with excellent success. In Fig. 7, we show forecasting results using a linear q(x) = x activation RC with still good results and associated with the theory in this work connecting to VAR(k).

Conclusion
The success of machine learning and artificial neural networks has lead to a clear and overwhelming widely adopted wave across so many areas where data is relevant and patterns are of interest. Dynamical systems is no exception and forecasting a dynamical system is a specific application that is broadly relevant and of interest to us here. The RNN framework is particularly relevant for a dynamical systems since the reserve memory aspect of the concept allows for a good framework and some aspects of delay embedding. However while the RNN tends to have many many parameters to fit, with the danger of overfitting always present, and in any case the large cost of optimization in the training phase, there is a surprising short cut. The echo-state/reservoir computing concept presumes to choose the weights for input layer and internal layer entirely randomly. Then only the output layer is trained. Furthermore that output layer training will be a linear least squares solve, rather than the usual nonlinear optimization necessary for the many parameters of the full nonlinear RNN. That this would allow a huge computational savings is clear. What is not clear and perhaps even a surprise, is how can this gross simplification still lead to useful results. While there have been a number of studies describing experimentally how to choose better random processes to define the random parameters, e.g. such as to emphasize sparsity, or to control the spectral radius, and other properties, in this work we have taken a different approach, which is not specifically to improve the performance of the concept but instead to give a partial explanation as to how the concept can work at all. After all, at first glance it may seem that it would be impossible that such a simplification could work. In this work, we have simplified the RC concept, allowing for the activation function to be an identity function instead of the more typical sigmoidal function. In this case, the process is entirely linear and so easier to study. As it turns out the RC still performs well as the nonlinear activation may be important for better performance, is not necessary for some good performance, and certainly leads to easier analysis. Herein we have been able to prove that the linear RC is in factly directly related to the more matured topic of VAR, vector autoregressive time series forecasting, and with all the related theory, such as the WOLD theorem through which we recall the representation theorem also now applies to the RC. Also we are able to make a direct connection to the increasingly popular DMD theory. So now with a linear RC analyzed from the start, including the randomly selected parameters, we offer a few explicit examples that also demonstrate that the concept works. Therefore, while our simplification to a linear RC is not proved to also explain the full nonlinear RC, we feel it does give some insight as to how a random partial selection can still be compensated, and also, we argue that a scaling step further connects the linear to a linearized sample to the sigmoidal activation. We hope in the future that we might relate this work more explicitly to explain the nonlinear RC success beyond the series explanation here, and also to explain more closely the role of memory beyond spectral radius, perhaps to embedding theory, and also perhaps to sufficient dimensional observation for an approximate isometric embedding in terms of the Johnson-LIndenstrauss theorem.

Data Availability Statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments
The author received funding from the Army Research Office (N68164-EG) and also DARPA.
We review how to numerically and stably compute the pseudo-inverse by the singular value decomposition, with regularized singular values (SVD). Reviewing the matrix theory of regularized pseudo-inverses for general matrices, if, Xb = z, X n×p , b p×1 , z n×1 , then if the SVD is X = U ΣV , with orthogonal matrices, n × n, U satisfies, U U T = U T U = I and p × p, V satisfies V V T = V T V = I, and Σ is n × p "diagonal" matrix of singular values, Then, . . .
The least squares estimator of Xb = z is, and we write the ridge regression Tikhonov regularized solution, The regularized pseudo-inverse X † λ is better stated in terms of the regularized singular values, by, and then, Throughout, since we will always refer to regularized pseudo-inverses, we will not emphasize this by abusing notation allowing that b * denotes b * λ even if only a very small λ > 0 is chosen, unless otherwise stated, λ = 1.0 × 10 −8 . This mitigates the tendency of overfitting or likewise stated in terms of zero or almost zero singular values that would otherwise appear in the denominators of Σ † . The theory is similar for n = p and n < p, as well as the scenario where z is not just a vector but a matrix, and likewise as in Eq. (9) where we refer to the transpose scenario.