{'fraction_non_alphanumeric': 0.04553952520371089, 'fraction_numerical': 0.01733747222124237, 'mean_word_length': 4.321850948000751, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Power of Linear Recurrent Neural Networks

Recurrent neural networks are a powerful means to cope with time series. We show how linear, i.e., linearly activated recurrent neural networks (LRNNs) can approximate any time-dependent function f(t) given by a number of function values. The approximation can effectively be learned by simply solving a linear equation system; no backpropagation or similar methods are needed. Furthermore, the size of an LRNN can be reduced significantly in one step, after inspecting the eigenvalues of the network transition matrix, by taking only the most relevant components. Therefore, in contrast to others, we do not only learn network weights but also the network architecture. LRNNs have interesting properties: They end up in ellipse trajectories in the long run and allow the prediction of further values and compact representations of functions. We demonstrate this by several experiments, among them multiple superimposed oscillators (MSO), robotic soccer, and predicting stock prices. LRNNs outperform the previous state-of-the-art for the MSO task with a minimal number of units.


Introduction
Deep learning in general means a class of machine learning algorithms that use a cascade of multiple layers of nonlinear processing units for feature extraction and transformation (Deng and Yu, 2014). The tremendous success of deep learning in diverse fields of artificial intelligence, such as computer vision and natural language processing, seems to depend on a bunch of ingredients: artificial, possibly recurrent neural networks (RNNs), with nonlinearly activated neurons, convolutional layers, and iterative training methods like backpropagation (Goodfellow et al., 2016). But which of these components are really essential for machine learning tasks such as time-series analysis?
Research in time series analysis and hence modeling dynamics of complex systems has a long tradition and is still highly active due to its crucial role in many real-world applications (Lipton et al., 2015) like weather forecast, stock quotations, comprehension of trajectories of objects and agents, or solving number puzzles (Ragni and Klein, 2011;Glüge and Wendemuth, 2013). The analysis of time series allows, among others, data compression, i.e., compact representation of time series, e.g., by a function f (t), and prediction of further values.
Numerous research addresses these topics by RNNs, in particular variants of networks with long short-term memory (LSTM) (Hochreiter and Schmidhuber, 1997). In the following, we consider an alternative, simple, yet very powerful type of RNNs: linear recurrent neural networks (LRNNs). Thus we only use linear activation, which allows us to minimize the network size in one step, namely by inspecting the eigenvalues of the network transition matrix (cf. Section 4.3). Therefore, in contrast to others, we do not only learn network weights but also the network architecture.
The rest of this paper is structured as follows: First, we briefly review related works (Section 2). We then formally introduce LRNNs as a special and simple kind of RNNs together with their properties, including the general network dynamics and their long-term behavior (Section 3). Afterwards, learning LRNNs is explained (Section 4). It is a relatively straightforward procedure which allows network size reduction; no backpropagation or gradient-descent method is needed. We then discuss results and experiments (Section 5), before we end up with conclusions (Section 6).

Echo State Networks
Echo state networks (ESNs) (Jaeger and Haas, 2004;Jaeger, 2007) play a significant role in RNN research as they provide an architecture and supervised learning principle for RNNs. They do so by driving a random, large, fixed RNN, called reservoir in this context, with the input signal which then induces in each neuron within this reservoir network a nonlinear response signal. They combine a desired output signal by a trainable linear combination of all response signals, allowing dimensionality reduction by so-called conceptors (Jaeger, 2014(Jaeger, , 2017. Xue et al. (2007) propose a variant of ESNs that work with several independent (decoupled) smaller networks. ESN-style (random) initialization has been shown effective for training RNNs with Hessian-free optimization (Martens and Sutskever, 2011). The latter paper addresses the problem of how to effectively train recurrent neural networks on complex and difficult sequence modeling problems which may contain long-term data dependencies. This can also be done with LRNNs (cf. MSO benchmark, Section 5.1). Tiňo (2018) considers the effect of weight changes in linear symmetric ESNs on (Fisher) memory of the network. From a solely theoretical point of view, Couillet et al. (2016) investigate the asymptotic performance of linear ESNs.

Recurrent Neural Networks
Simple RNNs are proposed by Elman (1990). By allowing them to accept sequences as inputs and outputs rather than individual observations, RNNs extend the standard feedforward multilayer perceptron networks. As shown in many sequence modeling tasks, data points such as video frames, audio snippets, and sentence segments are usually highly related in time. This results in RNNs being used as the indispensable tools for modeling such temporal dependencies. Linear RNNs and some of their properties (like short-term memory) are network architecture and hence its size is reduced in one step by analyzing the network transition matrix (cf. Section 4.3).

Autoregression
An autoregressive model is a representation of a type of random process (Akaike, 1969). It specifies that the output variable or a vector thereof depends linearly on its own previous values and on a stochastic term (white noise). In consequence, the model is in the form of a stochastic differential equation as in general (physical) dynamic systems (Colonius and Kliemann, 2014). An LRNN is also linearly activated, but its output does not only depend on its own previous values and possibly white noise but on the complete state of the possibly big reservoir whose dynamics is explicitly dealt with. In addition, the size of the network might be reduced in the subsequent process. We will continue the comparison between autoregression and LRNNs later (in Sections 4. 1 and 5.4).
Autoregressive frameworks are common in current machine-learning applications like language modeling, e.g., in the generative pre-trained transformer GPT-3, by next word prediction (Brown et al., 2020). The model applied in this context, however, is very complex (175 billion parameters) and non-linear. Related to autoregression is autoencoding of sequences. For this, Pasa and Sperduti (2014) show that linear autoencoders can be used for pre-training of RNNs, while we establish completely linear RNNs here. Furthermore, Sperduti (2006) gives an exact closed-form solution for the weights of the linear autoencoder, which is related to the approximation theorem for LRNNs (cf. Property 7).

Linear Recurrent Neural Networks
RNNs often host several types of neurons, each activated in a different manner (Elman, 1990;Hochreiter and Schmidhuber, 1997). In contrast to this, we here simply understand a homogeneous interconnected group of standard neurons as a neural network which may have arbitrary loops, akin to biological neuronal networks. We adopt a discrete time model, i.e., input and output can be represented by a time series and is processed stepwise by the network.
Definition 1 (time series). A time series is a series of data points in d dimensions S(0), . . . , S(n) ∈ R d where d ≥ 1 and n ≥ 0.
Definition 2 (recurrent neural network). A recurrent neural network (RNN) is a directed graph consisting of altogether N nodes, called neurons. x(t) denotes the activation of the neuron x at (discrete) time t. We may distinguish three groups of neurons (cf. Figure 1): • input neurons (usually without incoming edges) whose activation is given by an external source, e.g., a time series, • output neurons (usually without outgoing edges) whose activation represents some output function, and • reservoir or hidden neurons (arbitrarily connected) that are used for auxiliary computations. In ESNs, only output weights are trained and the hidden layer is also called reservoir.
x 1 x 3 x 2 1 0 2 1 0 1 1 1 1 Figure 2: LRNN for f (t) = t 2 with time step τ = 1. The input/output neuron x 1 is marked by a thick border. The initial values of the neurons at time t 0 = 0 are written in the nodes. The weights are annotated at the edges.
The edges of the graph represent the network connections. They are usually annotated with weights which are compiled in the transition matrix W . An entry w ij in row i and column j denotes the weight of the edge from neuron j to neuron i. If there is no connection, then w ij = 0. The transition matrix has the form containing the following weight matrices: • input weights W in (weights from the input and possibly the output to the reservoir), • output weights W out (all weights to the output and possibly back to the input), and • reservoir weights W res (weights within the reservoir, a matrix of size N res ×N res where N res is the number of reservoir neurons).
Let us now define the network activity: The initial configuration of the neural network is given by a column vector s with N components, called start vector. It represents the network state at the start time t = t 0 . Because of the discrete time model, the activation of a (non-input) neuron x at time t + τ (for some time step τ > 0) from the activation at time t of the neurons x 1 , . . . , x n that are connected to x with the weights w 1 , . . . , w n is computed as follows: x(t + τ ) = g w 1 x 1 (t) + · · · + w n x n (t) This has to be done simultaneously for all neurons of the network. Here g is a (real-valued) activation function. Although we will not make use of it, g may be different for different parts of the network. Usually, a nonlinear, bounded, strictly increasing sigmoidal function g is used, e.g., the logistic function, the hyperbolic tangent (tanh), or the softplus function (Goodfellow et al., 2016, Sect. 3.10). In the following, we employ simply the (linear) identity function and can still approximate arbitrary time-dependent functions (cf. Property 7).
Definition 3 (linear recurrent neural network). A linear recurrent network (LRNN) is an RNN with the following properties: 1. For the start time, it holds that t 0 = 0, and τ is constant, often τ = 1.
2. The initial state S(0) of the given time series constitutes the first d components of the start vector s.
3. For all neurons we have linear activation, i.e., everywhere g is the identity.
4. The weights in W in and W res are taken randomly, independently, and identically distributed from the standard normal distribution, i.e., the Gaussian distribution with mean µ = 0 and standard deviation σ = 1, whereas the output weights W out are learned (cf. Section 4.1).
5. There is no clear distinction between input and output but only one joint group of d input/output neurons. They may be arbitrarily connected like the reservoir neurons. We thus can imagine the whole network as a big reservoir because the input/output neurons are not particularly special.
LRNNs can run in one of two modes: either receiving input or generating (i.e., predicting) output. In output generating mode, Equation 2 is applied to all neurons including the input/output neurons, whereas in input receiving mode the activation of every input/output neuron x at time t is always overwritten with the respective input value at time t given by the time series S.

Examples
Example 1. The function f (t) = t 2 can be realized by an LRNN (in output generating mode) with three neurons (cf. Figure 2). The respective transition matrix W and start vector s are: Consequently, starting at t = 0 with time step τ = 1, we have: • x 3 (0) = 1, x 3 (t + 1) = x 3 (t), and hence x 3 (t) = 1 in general.
Thus, in the neuron x 1 , the function f (t) is computed. This is the only input/output neuron.

Remark 1.
LRNNs are well suited to represent differential equations and to solve them numerically. To see this, consider the homogeneous linear differential equation where c k ∈ R are constant coefficients, x (k) (t) is the k-th derivative of the function x with respect to time t, and n > 0. For a function x and its derivativeẋ, we havė for small time steps τ . If we consider the difference of Equation 3 between the times t + τ and t, we obtain, after dividing by τ and applying Equation 4 to x (k) (t) for k ≥ 0: Therefore, we can solve differential equations approximately by LRNNs with start vector s, which should satisfy Equation 3, and the following transition matrix: Example 2. The exponential function exp(t) = e t can be defined by the differential equationẋ(t) = x(t), i.e., we have c 0 = 1 and c 1 = −1 in Equation 3. In consequence, according to Remark 1 and because exp(0) = 1, the transition matrix W and start vector s of the corresponding LRNN are: Induction over time yields immediately x(t) =ẋ(t) = (1 + τ ) t/τ ≈ e t for small τ (according to Euler) as expected.
The strong relationship between RNNs and differential equations is already known (Kruse et al., 2016, Sect. 9) as well as the extraction of eigenvalues to describe dynamical systems (Strogatz, 2015, Sect. 5). Nevertheless, as we will show in the rest of this paper, the combination of both provides an effective method for network size reduction (cf. Section 4.3) and therefore seems to be worthwhile to be considered by the machine learning community in more detail.

Network Dynamics
An LRNN runs through network states f (t) for t ≥ 0. It holds (in output generating mode) and hence simply f (t) = W t · s. Property 1. Let W = V · J · V −1 be the Jordan decomposition of the transition matrix W where J is the direct sum, i.e., a block diagonal matrix, of one or more Jordan blocks (Horn and Johnson, 2013, Sect. 3 in general with different sizes m × m and different eigenvalues λ. Then we have: If we decompose V into matrices v of size N × m and the column vector V −1 · s into a stack of column vectors w of size m, corresponding to the Jordan blocks in J, then f (t) can be expressed as a sum of vectors u = v · J m (λ) t · w where the Jordan block powers are upper triangular Toeplitz matrices, i.e., in which each descending diagonal from left to right is constant, with: (Horn and Johnson, 2013, Sect. 3.2.5) Remark 2. Although the parameter t is discrete, i.e., a nonnegative integer number, the values of f (t) = W t · s can also be computed for t ∈ R and are always real. For this, we consider the Jordan block powers from Equation 5: • The definition of the binomial coefficient t is applicable for real and even complex t and nonnegative integer k. For negative k, we have t k = 0.
• For real matrices W , there are always complex conjugate eigenvalue pairs λ and λ and corresponding complex coefficients c and c (resulting from the vectors u in Property 1). With c = |c| e iψ and λ = |λ| e iω , we get c λ t +c λ t = |c| |λ| t cos(ωt+ψ), applying Euler's formula. This obviously is defined for all t ∈ R and always yields real-valued f (t).
• Negative real eigenvalues, i.e., the case λ < 0, should be treated in a special way, namely by replacing λ t by |λ| t cos(πt). Both terms coincide for integer t, but only the latter is real-valued for all t ∈ R. The powers of positive real eigenvalues λ are always positive and real and hence need no further consideration.
A Jordan decomposition exists for every square matrix W (Horn and Johnson, 2013, Theorem 3.1.11). But if W has N distinct eigenvectors, there is a simpler decomposition, called eigendecomposition. The transition matrix W is diagonalizable in this case, i.e., similar to a diagonal matrix D, and the network dynamics can be directly described by means of the eigenvalues and eigenvectors of W : Property 2. Let W = V · D · V −1 be the eigendecomposition of the transition matrix W with column eigenvectors v 1 , . . . , v N in V and eigenvalues λ 1 , . . . , λ N , on the diagonal of the diagonal matrix D, sorted in decreasing order with respect to their absolute values. Like every column vector, we can represent the start vector s as linear combination of the eigenvectors, namely as s =

Long-Term Behavior
Let us now investigate the long-term behavior of an RNN (run in output generating mode) by understanding it as an (autonomous) dynamic system (Colonius and Kliemann, 2014;Strogatz, 2015). We will see (in Property 4) that the network dynamics may be reduced to a very small number of neurons in the long run. They determine the behavior for t → ∞. Nevertheless, for smaller t, the use of many neurons is important for computing short-term predictions.
Property 3. In none of the N dimensions f (t) = W t · s grows faster than a polynomial and only single-exponential in t.
Proof. Let f k (t) denote the value of the k-th dimension of f (t), λ be the eigenvalue of W with maximal absolute value and m be the maximal (geometric) multiplicity of the eigenvalues of the transition matrix W . Then, from Property 1, we can easily deduce as asymptotic behavior for large t.
Property 4. Consider an RNN whose transition matrix W is completely real-valued, has (according to Property 2) an eigendecomposition W = V · D · V −1 with unit spectral radius (i.e., the largest absolute value of all eigenvalues is 1), and all eigenvalues are distinct, e.g., a pure random reservoir. Then, almost all terms x k λ k t v k in Equation 6 vanish for large t because for all eigenvalues λ k with |λ k | < 1 we have lim t→∞ λ k t = 0. Although a general real matrix can have more than two complex eigenvalues which are on the unit disk, for a pure random reservoir as considered here, almost always, i.e., with probability close to 1, only the (largest) eigenvalues λ 1 and possibly λ 2 have the absolute values 1. In consequence, we have one of the following cases: 1. λ 1 = +1. In this case, the network activity contracts to one point, i.e., to a singularity: lim This means we have an oscillation in this case. The dynamic system alternates between the two points ±x 1 v 1 .
3. λ 1 and λ 2 are two (properly) complex eigenvalues with absolute value 1. Since W is a real-valued matrix, the two eigenvalues as well as the corresponding eigenvectors v 1 and v 2 are complex conjugate with respect to each other. Thus, for large t, we have an ellipse trajectory We consider now the latter case in more detail, converting all complex values into real ones (cf. Section 3.4): In this case, the matrixṼ consists of the two complex conjugate eigenvectors v = v 1 and v = v 2 . From all linear combinations of both eigenvectors (κ) = κv+κ v for κ ∈ C with |κ| = 1, we can determine the vectors with extremal lengths. Since we only have two real-valued dimensions because clearly the activity induced by the real matrix W remains in the real space, there are two such vectors. Let now κ = e iϕ = cos(ϕ) + i sin(ϕ) (Euler's formula) and v and v be the real and imaginary parts of v, respectively, i.e., v = v + iv and v = v − iv . Then, for the square of the vector length, it holds: To find out the angle ϕ with extremal vector length of (κ), we have to investigate the derivative of the latter term above with respect to ϕ and compute its zeros. This yields 2 −2 sin(2ϕ)(v 2 −v 2 ) − 2 cos(2ϕ)(2v ·v ) = 0 and thus tan(2ϕ) = −2v ·v v 2 −v 2 . Because of the periodicity of the tangent function there are two main solutions for κ that are orthogonal to each other: κ 1 = e iϕ and κ 2 = e i(ϕ+π/2) . They represent the main axes of an ellipse. All points the dynamic system runs through lie on this ellipse in the long run. The length ratio of the ellipse axes is µ = (κ 1 ) / (κ 2 ) . We normalize both vectors to unit length and put them in the matrixV = (κ 1 ) / (κ 1 ) (κ 2 ) / (κ 2 ) . We now build a matrixD, similar toD but completely real-valued which states the ellipse rotation. The rotation speed can be derived from the eigenvalue λ 1 . In each step of length τ , there is a rotation by the angle ωτ where ω is the angular frequency which can be determined from the equation λ 1 = |λ 1 | e iωτ . The two-dimensional ellipse trajectory can be stated by two (co)sinusoids: f (t) = a cos(ω t) b sin(ω t) with a, b > 0. Applying the addition theorems of trigonometry, we get: Figure 3: Dynamic system behavior in two dimensions (qualitative illustration): In the long run, we get an ellipse trajectory (dotted/red), though the original data may look like random (dots/black). Projected to one dimension, we have pure sinusoids with one single angular frequency, sampled in large steps (solid/blue).
From this, we can read off the desired ellipse rotation matrixD as indicated above and µ = a/b. Finally, we can determine the corresponding two-dimensional start vectorx by solving the equationV ·D ·x =Ṽ ·D ·x. In summary, we have for large t. Every RNN with many neurons can thus be approximated by a simple network with at most two neurons, defined by the matrixD and start vectorx. The output values can be computed for all original dimensions by multiplication with the matrixV . They lie on an ellipse in general. Nonetheless, in the beginning, i.e., for small t, the dynamics of the system is not that regular (cf. Figure 3). But although Property 4 states the asymptotic behavior of random RNNs with unit spectral radius, interestingly the network dynamics converges relatively fast to the final ellipse trajectory: The (Euclidean) distance between the actual value f (t) (according to Equation 6) and its approximation by the final ellipse trajectory (Equation 7) is almost zero already after a few hundred steps (cf. Figure 4). Of course this depends on the eigenvalue distribution of the transition matrix (Tao et al., 2010). The long-term behavior of LRNNs is related to that of ESNs. For the latter, usually the activation function is tanh and the spectral radius is smaller than 1. Then reservoirs with zero input collapse because of |tanh(z)| ≤ z for all z ∈ R but the convergence may be rather slow. This leads to the so-called echo state property (Manjunath and Jaeger, 2013): Any random initial state of a reservoir is forgotten such that, after a washout period, the current network state is a function of the driving input. In contrast to ESNs, LRNNs have linear activation and usually a spectral radius of exactly 1. But as we have just shown, there is a similar effect in the long run: The network activity reduces to at most two dimensions which are independent from the initial state of the network.

Real-Valued Transition Matrix Decomposition
For real-valued transition matrices W , it is possible to define a decomposition that, in contrast to the ordinary Jordan decomposition in Property 1, solely makes use of realvalued components, adopting the so-called real Jordan canonical form (Horn and Johnson, 2013, Sect. 3.4.1) of the square matrix W . For this completely real-valued decomposition, the Jordan matrix J is transformed as follows: 1. A Jordan block with real eigenvalue λ remains as is in J.
2. For complex conjugate eigenvalue pairs λ = λ + iλ and λ = λ − iλ with λ , λ ∈ R, the direct sum of the corresponding Jordan blocks J m (λ) and J m (λ) is replaced by a real Jordan block: This procedure yields the real Jordan matrix J. In consequence, we have to transform V also into a completely real-valued form. For a complex conjugate eigenvalue pair λ and λ, the corresponding two eigenvectors in V could be replaced by the vectors inV (cf. Section 3.3). The subsequent theorem shows a more general way: The matrix V from Property 1 is transformed into a real matrix A and, what is more, the start vector s can be replaced by an arbitrary column vector y with all non-zero entries.
Property 5. Let W = V · J · V −1 be the (real) Jordan decomposition of the transition matrix W and s the corresponding start vector. Then for all column vectors y of size N with all non-zero entries, there exists a square matrix A of size N × N such that for all integer t ≥ 0 we have: Proof. We first prove the case where the Jordan matrix J only contains ordinary Jordan blocks as in Property 1, i.e., possibly with complex eigenvalues on the diagonal. Since J is a direct sum of Jordan blocks, it suffices to consider the case where J is a single Jordan block because, as the Jordan matrix J, the matrices A and also B (see below) can be obtained as direct sums, too.
In the following, we use the column vectors y = y 1 · · · y N with all non-zero entries, Since the main diagonal of the left matrix contains no 0s because y N = 0 by precondition, there always exists a solution for b (Horn and Johnson, 2013, Sect. 0.9.3). Then A = V · B does the job: The generalization to the real Jordan decomposition is straightforward by applying the fact that for complex conjugate eigenvalue pairs λ and λ the matrix M from above in a real Jordan block is similar to the diagonal matrix D =

Learning LRNNs
Functions can be learned and approximated by LRNNs in two steps: First, as for ESNs (Jaeger and Haas, 2004), we only learn the output weights W out (cf. Section 4.1); all other connections remain unchanged. The input weights W in and reservoir weights W res are arbitrary random values that are not changed (cf. Definition 3). Nevertheless, in order to obtain better numerical stability during the computation, they are adjusted as follows: • In the presence of linear activation, the spectral radius of the reservoir weights matrix W res is set to 1. Otherwise, with increasing t, the values of W t explode if the spectral radius is greater or vanish if the spectral radius is smaller than 1 (cf. Section 3.3).
• The norms of the vectors in W in and W res should be balanced (Koryakin et al., 2012). To achieve this, we initialize the reservoir neurons such that the reservoir start vector r 0 (with N res components) has unit norm by setting: It is part of the start vector s = S(0) r 0 (cf. Section 3).
• We usually employ fully connected graphs, i.e., all, especially the reservoir neurons are connected with each other because the connectivity has nearly no influence on the best reachable performance (Koryakin et al., 2012).
Second, if possible, we reduce the network size (cf. Section 4.3); this often leads to better generalization and avoids overfitting. Thus, in contrast to many other approaches, the network architecture is changed during the learning process. In contrast to other approaches, we do not do this by incremental derivation from the original network but in only one step.

Learning the Output Weights
To learn the output weights W out , we run the input values from the time series S(0), . . . , S(n) through the network (in input receiving mode), particularly through the reservoir. This means, we build the sequence of corresponding reservoir states R(0), . . . , R(n) where the reservoir start vector r 0 in principle can be chosen arbitrarily but with all non-zero entries (cf. Property 5): We want to predict the next input value S(t + τ ), given the current input and reservoir states S(t) and R(t). To achieve this, we comprise all but the last input and reservoir states in one matrix X with: Each output value corresponds to the respective next input value S(t + τ ). Therefore, we compose another matrix Y out = S(1) · · · S(n) consisting of the predicted next values of the time series S where the first value S(0) clearly has to be omitted because it cannot be predicted. We compute Y out (t) = S(t + τ ) from X(t) by assuming a linear dependency: Its solution can easily be determined as W out = Y out /X, where / denotes right matrix division, i.e., the operation of solving a linear equation system, possibly applying the least squares method in case of an overdetermined system, as implemented in many scientific programming languages like Matlab (Higham and Higham, 2017) or Octave (Eaton et al., 2017). Prediction of further values is now possible (in output generating mode) as follows: with W as in Equation 1  (12) Property 6 (treatment of multiple sequences). It is also possible to learn from multiple sequences at once. For this, let several time series S 1 , . . . , S K in d dimensions with (not necessarily identical) lengths n 1 , . . . , n K be given. For each S k with 1 ≤ k ≤ K, we determine: • the sequence of corresponding reservoir states R k (according to Equation 8), taking always the same reservoir start vector r 0 , • the corresponding input matrix X k (according to Equation 9), and • the corresponding predicted output matrix Y out k (according to Equation 10).
We aggregate the input and output matrices to X = X 1 · · · X K and Y out = Y out 1 · · · Y out K with n 1 + · · · + n K columns each. Solving the linear matrix equation Y out = W out · X (identical with Equation 11) finally yields the output weight matrix W out .
This first phase of the learning procedure is related to a linear autoregressive model (Akaike, 1969). However, one important difference to an autoregressive model is that for LRNNs the output does not only depend on its own previous values and possibly white noise but on the complete state of the possibly big reservoir whose dynamics is explicitly dealt with in the reservoir matrix W res . The reservoir effectively allows us to do arbitrary auxiliary computation such that any (non-linear) function f (t) can be approximated by an LRNN (cf. Property 7).

An Approximation Theorem
Property 7. From a function f (t) in d ≥ 1 dimensions, let a series of function values f (t 0 ), . . . , f (t n ) be given. Then there is an LRNN with the following properties: 1. It runs exactly through all given n + 1 function values, i.e., it approximates f (t).
2. It can effectively be learned by the above-stated solution procedure (Section 4.1).
Proof. The procedure for learning output weights (cf. Section 4.1) uses the reservoir state sequence as part of the coefficient matrix X which reduces to at most two dimensions in the long run however -independent of the number of reservoir neurons (cf. Section 3.3).
Thus, the rank of the coefficient matrix X is not maximal in general and in consequence the linear equation system from Equation 11 often has no solutions (although we may have an equation system with the same number of equations and unknowns). A simple increase of the number of reservoir neurons does not help much. Therefore, we shall apply the learning procedure in a specific way, learning not only the output weights W out as in ESNs but the complete transition matrix W as follows: First, we take the series of function values f (t 0 ), . . . , f (t n ) and identify them with the time series S(0), . . . , S(n). Let then R(0) · · · R(n) be a random reservoir state sequence matrix for N res reservoir neurons, considered as additional input in this context. If all elements of this matrix are random numbers, independently and identically distributed from the standard normal distribution, its rank is almost always maximal. Let now d be the rank of the matrix S(0) · · · S(n − 1) and: We now have to solve the linear matrix equation W · X = Y with X as in Equation 9. If w 1 , . . . , w N and y 1 , . . . , y N denote the row vectors of the matrices W and Y , respectively, then this is equivalent to simultaneously solving the equations w k · X = y k , for 1 ≤ k ≤ N . We must ensure that there exists almost always at least one solution. This holds if and only if the rank of the coefficient matrix X is equal to the rank of the augmented matrix M k = X y k , for every k. We obtain the equation rank(X) = min(d + N res , n) = rank(M k ) = min(d + N res + 1, n). From this, it follows that reservoir neurons have to be employed to guarantee at least one exact solution for the w k . By construction, the neural network with transition matrix W runs exactly through all given n + 1 function values f (t).
Remark 3. The just sketched proof of Property 7 suggests a way to learn the input and reservoir weights. This topic has also been investigated by Palangi et al. (2013) for ESNs with nonlinear activation function in the reservoir. However, for LRNNs, the given input and reservoir weights W in and W res together with the learned output weights W out already provide the best approximation of the function f (t). There is no need to learn the input and reservoir weights, simply because LRNNs are completely linearly activated RNNs (including the reservoir). If one tries to learn W in and W res taking not only the output time series S but additionally the reservoir state time series R into account, then exactly the given input and reservoir weights are learned (if Equation 13 holds). Only with nonlinear activation there is a learning effect. Nonetheless, W in and W res can be learned by the procedure sketched in the proof, if they are not given in advance, starting with a random reservoir state sequence. But our experiments indicate that this procedure is less numerically stable than the one with given, i.e., predefined input and reservoir weights and unit spectral radius for the reservoir (cf. Section 4.1).
Property 7 is related to the universal approximation theorem for feedforward neural networks (Hornik, 1991). It states that a (non-recurrent) network with a linear output layer and at least one hidden layer activated by a nonlinear, sigmoidal function can approximate any continuous function on a closed and bounded subset of the R n from one finite-dimensional space to another with any desired non-zero amount of error, provided that the network is given enough hidden neurons (Goodfellow et al., 2016, Sect. 6.4.1). Since RNNs are more general than feedforward networks, the universal approximation theorem also holds for them (Maass et al., 2002). Any measurable function can be approximated with a (general) recurrent network arbitrarily well in probability (Hammer, 2000).
Because of the completely linear activation, LRNNs cannot compute a nonlinear function f (x) from the (possibly multi-dimensional) input x. Nevertheless, they can approximate any (possibly nonlinear) function over time f (t), as Property 7 shows. Another important difference between LRNNs and nonlinearly activated feedforward neural networks is that LRNNs can learn the function f (t) efficiently. No gradient-descent method like backpropagation is required; we just have to solve a linear equation system. Hence learning is as easy as learning a single-layer perceptron, which however is restricted in expressibility because only linearly separable functions can be represented.
Remark 4. Any time series S(0), . . . , S(n) can be generated by employing a backward shift matrix, i.e., a binary matrix with 1s on the subdiagonal and 0s elsewhere (Horn and Johnson, 2013, Sect. 0.9.7), as transition matrix W and s = S(0) · · · S(n) as start vector. But such a network clearly would have no ability to generalize to future data. Fortunately, this does not hold for a transition matrix W learned by the procedure in Section 4.1. Furthermore, the eigenvalue spectrum of the backward shift matrix is empty, whereas that of the learned W is not, which is important for network size reduction introduced in Section 4.3.

Network Size Reduction
To approximate a function exactly for sure, we need a large number N res of reservoir neurons in Property 7. It is certainly a good idea to lower this number. One could do this by simply taking a smaller number of reservoir neurons, but then a good approximation cannot be guaranteed. In what follows, we therefore reduce the dimensionality of the transition matrix W in a more controlled way -after learning the output weights. Our procedure of dimensionality reduction leads to smaller networks with sparse connectivity. In contrast to other approaches, we do not learn the new network architecture by incremental derivation from the original network, e.g., by removing unimportant neurons or weights, but in only one step by inspecting the eigenvalues of the transition matrix.
For ESNs, dimensionality reduction is considered, too, namely by means of so-called conceptors (Jaeger, 2014(Jaeger, , 2017Krause et al., 2021). These are special matrices which restrict the reservoir dynamics to a linear subspace that is characteristic for a specific pattern. However, as in principal component analysis (PCA) (Jolliffe, 2011), conceptors reduce only the spatial dimensionality of the point cloud of the given data. In contrast to this, for LRNNs, we reduce the transition matrix W and hence take also into account the temporal order of the data points in the time series. By applying insights from linear algebra, the actual network size can be reduced and not only the subspace of computation as with conceptors.
Property 8. By Property 5, the function f (t) = W t · s can be rewritten by means of the Jordan matrix of the transition matrix W as A · J t · y, where the start vector can be chosen as non-zero constant, e.g., y = 1 · · · 1 . Furthermore, by Property 1, f (t) can be expressed as a sum of vectors u = v · J m (λ) t · w where w is constant because it is part of the start vector y. Then it follows from Property 4 that for large t the contribution of a Jordan component vanishes if v ≈ 0 and/or |λ| 1. In consequence, we can omit all Jordan components causing only small errors, until a given threshold is exceeded. The error E of a network component corresponding to a Jordan block J m (λ) can be estimated by the root-mean-square error (RMSE) normalized to the number of all sample components between input x and predicted output y: In practice, we omit all network components with smallest errors as long as the RMSE is below a given threshold θ. Network components that are not omitted are considered relevant. Thus, from A, J, and y (according to Property 5), we successively derive reduced matrices A and J and the vector y as follows: • From A, take the rows corresponding to the input/output components and the columns corresponding to the relevant network components.
• From J, take the rows and columns corresponding to the relevant network components.
• From y, take the rows corresponding to the relevant network components.
Note that the dimensionality reduction does not only lead to a smaller number of reservoir neurons, but also to a rather simple network structure: The transition matrix J (which comprises the reservoir weights W res of the reduced network) is a sparse matrix with nonzero elements only on the main and immediately adjacent diagonals. Thus, the number of connections is in O(N ), i.e., linear in the number of reservoir neurons, not quadratic -as in general. Figure 5 summarizes the overall learning procedure for LRNNs including network size reduction. It has been implemented by the authors in Octave. Note that, although the Jordan matrix J (cf. Property 1) may contain eigenvalues with multiplicity greater than 1, Octave does not always calculate exactly identical eigenvalues in this case. Therefore, we cluster the computed eigenvalues as follows: If the distance in the complex plane between eigenvalues is below some given small threshold δ, they are put into the same cluster which eventually is identified with its centroid. Thus it is a kind of single linkage clustering (Gower and Ross, 1969). The complete implementation of the learning procedure together with some case studies (cf. Section 5) is available under the following link: https://bitbucket.org/oliverobst/decorating 4.4 Complexity and Generalization of the Procedure Property 9. In both learning steps, it is possible to employ any of the many available fast and constructive algorithms for linear regression and eigendecomposition (Demmel et al., 2007). Therefore, the time complexity is just O(N 3 ) for both output weights learning and % d-dimensional function, given sampled, as time series S = [f (0) . . . f (n)] % random initialization of reservoir and input weights W in = randn(N, d) W res = randn(N res , N res ) % learn output weights by linear regression X = W t · s t=0,...,n Y out = S(1) · · · S(n) W out = Y out /X % transition matrix and its decomposition % network size reduction y = 1 · · · 1 Y = J t · y t=0,...,n A = X/Y with rows restricted to input/output dimensions reduce(A, J, y) to relevant components such that RMSE(S,Out) < θ Figure 5: Pseudocode for learning LRNNs including network size reduction. In the last statement, actually a binary search algorithm is employed for determining the relevant network components with smallest errors. For this, the network components are sorted by their error E.
dimensionality reduction. In theory, if we assume that the basic numerical operations like + and · can be done in constant time, the asymptotic complexity is even a bit better. In practice, however, the complexity depends on the bit length of numbers in floating point arithmetic, of course, and may be worse hence. The size of the learned network is in O(N ) (cf. Section 4.3).
Note that, in contrast, feedforward networks with three threshold neurons already are NP-hard to train (Blum and Rivest, 1992). This results from the fact that the universal approximation theorem for feedforward networks differs from Property 7 because the former holds for multi-dimensional functions and not only time-dependent input. In this light, the computational complexity of O(N 3 ) for LRNNs does not look overly expensive. It dominates the overall time complexity of the whole learning procedure because it is not embedded in a time-consuming iterative learning procedure (like backpropagation) as in other state-ofthe-art methods. In this interval, all values of the parabola (solid/blue) are greater than or equal to those of the sinusoid (dashed/red).
Remark 5. We observe that most of the results presented in this paper still hold if the transition matrix W contains complex numbers. This means in particular that also complex functions can be learned (from complex-valued time series) and represented by LRNNs (Property 7). Nonetheless, the long-term behavior of networks with a random complex transition matrix W differs from the one described in Section 3.3 because then there are no longer pairs of complex conjugate eigenvalues.

Experiments
In this section, we demonstrate evaluation results for learning and predicting time series, approximating them by a function f (t) represented by an LRNN, for several tasks. We consider the following benchmarks: multiple superimposed oscillators, number puzzles, robot soccer simulation, and predicting stock prices. All experiments are performed with a program written by the authors in Octave (Eaton et al., 2017) (cf. Section 4.3). Let us start with an example that illustrates the overall method.
To investigate this, we sample both graphs for t ∈ [0, 1] with time step τ = 0.01. After that, we learn the output weights W out (cf. Section 4.1), starting with a large enough reservoir consisting of up to N res = 100 neurons (cf. Property 7). Finally, we reduce the size of the overall transition matrix W with precision threshold θ = 0.01 and cluster threshold δ = 0.03 (cf. Section 4.3). Minimal LRNNs consist of N 1 = 3 neurons for the parabola (cf. Example 1) and N 2 = 2 neurons for the sinusoid (cf. Section 3.3). The networks of minimal size are learned already with N res = 40 reservoir neurons before network size reduction in about 77% (parabola) or 99% (sinusoid) of the trials (cf. Figure 7). Learning the parabola is more difficult because the corresponding transition matrix W (cf. Example 1) has no proper eigendecomposition according to Property 2 but only a Jordan decomposition according to Property 1.
Various publications have investigated the MSO problem with different numbers of sinusoids. We concentrate here solely on the most complex case K = 8 whose graph is shown in Figure 8 because in contrast to other approaches it is still easy to learn for LRNNs. Applying the LRNN learning procedure with precision threshold θ = 0.5, we arrive at LRNNs with only N = 16 reservoir neurons and an RMSE less than 10 −5 , if we start with a large enough reservoir (cf. Figure 9). Since two neurons are required for each frequency (cf. Section 3.3), 2 K = 16 is the minimal size. Thus LRNNs outperform the previous state-of-the-art for the MSO task with a minimal number of units. Koryakin et al. (2012) report N res = 68 as the optimal reservoir size for ESNs, but in contrast to our approach, this number is not further reduced. In general, an LRNN with 2K neurons suffices to represent a signal, which might be a musical harmony (Stolzenburg, 2017), consisting of K sinusoids. It can be learned by the LRNN learning procedure with dimension reduction.

Solving Number Puzzles
Example 5. Number series tests are a popular task in intelligence tests. The function represented by a number series can often be learned also by artificial neural networks, in particular RNNs. Glüge and Wendemuth (2013) list 20 number puzzles (cf. Ragni and Klein, 2011). Among them are the series: 28,33,31,36,34,39,37,42] 3,6,12,24,48,96,192,384] [6,9,18,21,42,45,90,93] f (t) = 2f (t − 2) + 4.5 + 1.5(−1) t−1 8,12,16,20,24,28,32,36] We apply the LRNN learning procedure to all 20 number puzzles taking small reservoirs because the number series are short. As a side effect, this leads to learning more general functions, which seems to be fully adequate because number puzzles are usually presented to humans. The first 7 of 8 elements of each series are given as input. In each trial, we repeatedly generate LRNNs, until the RMSE is smaller than θ = 0.1. Then the last (8-th) element of the respective series is predicted (according to Equation 12) and rounded to the nearest integer because all considered number series are integer. Table 1 lists the percentages of correct predictions of the last element for different settings. Here the series with definitions recurring to f (t − 2) but not f (t − 1), e.g., S 8 and S 15 , turned out to be the most difficult. If we now just add the previous values of the time series, i.e., f (t − 2), as clue to the input, then the correctness of the procedure increases significantly: For 19 of 20 number puzzles, the most frequently predicted last element (simple majority) is the correct one. It is predicted in 76.5% on average over all trials and number puzzles. Let us remark that the whole evaluation with altogether 20 · 50 · 1000 = 1 000 000 trials including possibly repeated network generation and network size reduction ran in only a few minutes on standard hardware.

Replaying Soccer Games
RoboCup (Kitano et al., 1997) is an international scientific robot competition in which teams of multiple robots compete against each other. Its different leagues provide many sources of robotics data that can be used for further analysis and application of machine learning. A soccer simulation game lasts 10 mins and is divided into 6000 time steps where the length of each cycle is 100 ms. Logfiles contain information about the game, in particular about the current positions of all players and the ball including velocity and orientation for each cycle.   1.5% 0.5% 0.6% 0.5% 57.2% Table 1: Percentages of correct predictions of the last element for 20 number puzzles (Ragni and Klein, 2011;Glüge and Wendemuth, 2013) in 1000 trials for different settings: (a) with fixed reservoir size N res = 3, 4, 5; (b) with network size reduction starting with N res = 7 reservoir neurons; (c) same procedure but in addition always the previous series value is used as clue. Michael et al. (2019) describe a research dataset using some of the released binaries of the RoboCup 2D soccer simulation league (Chen et al., 2003;Gabel et al., 2017) from 2016 and 2017 . In our experiments we evaluated ten games of the top-five teams (available from https://bitbucket.org/oliverobst/robocupsimdata), considering only the (x, y)-coordinates of the ball and the altogether 22 players for all time points during the so-called "play-on" mode. For LRNN learning, we use only every 10-th time step of each game with d = 2 + 2 · 22 = 46 input dimensions and start with a reservoir consisting of N res = 500 neurons. We repeat the learning procedure until the RMSE is smaller than 1; on average, already two attempts suffice for this. This means, if we replay the game by the learned LRNN (in output generating mode) then on average the predicted positions deviate less than 1 m from the real ones (Euclidean distance) -over the whole length of the game (cf. Figure 10). Network size reduction leads to significantly less neurons compared to the original number N = 46 + 500 = 546 -on average 29.2% if we concentrate on the relevant components for the ball trajectory (cf.   (2) after dimensionality reduction where θ = 1 m. The network size can be reduced significantly -29.2% on average (last column).
(cf. Property 9). Thus the LRNN model is definitely smaller than the original time series representation of a game. The complete learning procedure runs in less than a minute on standard hardware. Property 6 shows how we can learn from multiple time series at once. This is also helpful here because by this procedure we can investigate the overall behavior of a specific robot soccer agent. As example for this, we consider the trajectories of the goalkeeper of the RoboCup simulation team FRA-UNIted during the seeding and the qualifying round of RoboCup Japan Open 2020 (see http://bit.ly/japanopen2020ssim). For learning one LRNN from this, we employ a reservoir with N res = 1000 neurons, adopt again a maximum threshold for the RMSE of θ = 1 m, and only use every 20-th step of each of the 7 games. The overall trajectory of the FRA-UNIted goalkeeper can be learned easily then (cf. Figure 11). From this, one may conclude that the goalkeeper takes up three basic positions in front of the goal, does not approach the centre line more than about 30 m and hardly leaves the centre line.

Predicting Stock Prices
Stock price prediction is a topic that receives a considerable amount of attention in finance. Complexity of markets resulting in multiple and sudden changes in trends of stock prices make their prediction a challenge. Consequently, a number of different approaches and methods have been developed, with the autoregressive integrated moving average (ARIMA) model (cf. Section 2.3) being a popular choice. Here we analyze 30 different stocks using the closing stock prices 2016-2019 from https://de.finance.yahoo.com/ and compare their ARIMA prediction with LRNNs (cf. Litz, 2020).
A standard ARIMA(p, d, q) model consists of autoregression AR(p) and a moving average MA(q). The parameter p describes the history length (lag order) used to predict the stock price at time t. We have f (t) = c 1 f (t − 1) + · · · + c p f (t − p) + e t where c 1 , . . . , c p are (real-valued) autocorrelation coefficients and e t , the residuals, are Gaussian white noise.
In the moving average process MA(q), the value q specifies the number of past residuals considered for prediction. An underlying trend of the stock is modeled by using a drift, i.e., a constant that extends the term. This procedure is particularly well-suited for stationary time series. Most financial time series, however, exhibit non-stationary behavior and thus require a transformation to make them stationary. This is achieved by investigating the derivatives of the series, with the order of this process given by the parameter d.
In our case, we predict each stock using an ARIMA model, implemented with the auto.arima package from Hyndman and Khandakar (2008), with seasonal patterns modeled using Fourier terms (Hyndnman and Athanasopoulos, 2013, p. 321). To fit the model, we split stock prices into training and test sets, with the first 80% of each series for training and the final 20% used for evaluation. Figure 12 shows predictions of the ARIMA model for Airbus (AIR.DE). We use the first 610 points from this series to build an ARIMA(2,1,2) model with drift and the subsequent 152 data points to evaluate our prediction using the RMSE as metric. Figure 13 shows the prediction using an LRNN for the same task. In this example, the RMSE for LRNNs is E test = 8.07 Euro and better than the RMSE of the ARIMA model with E test = 9.05 Euro.
For a representative comparison, we calculate the RMSE of predictions on every stock in the set. The average RMSE using LRNNs is E test = 18.40 Euro, lower than the average   (Litz, 2020). With an average stock price of 286.71 Euro of all stocks in the set, the average deviation is only 6.1%. Apart from the good prediction results, LRNNs have the advantage that they allow the prediction of multiple stocks at the same time. An LRNN can read in 30 stocks and predict each of them concurrently. For a concurrent forecast for 60 steps, LRNNs using 600 neurons achieved an average RMSE of E test = 30.02 Euro with 240 training steps. Compared to ARIMA, LRNNs have also an advantage when it comes to the number of hyper-parameters that have to be tuned. The LRNN model is robust when it comes to choosing the number of neurons, whereas the ARIMA model requires the adjustment of many parameters (e.g., for seasonal patterns). With an increase in number of hyper-parameters the compute time for ARIMA also increases. In our example, LRNNs are computed about 15 times faster than the ARIMA models with the selected number of Fourier terms.
The review by Shah et al. (2019) shows that predicting stock prices remains a challenging problem, especially for a longer timeframe, which we have investigated here. Contemporary research often uses complex models, ranging from LSTM RNNs (Nelson et al., 2017;Roondiwala et al., 2017) to attention-based models that use further information about the events that drive the stock prices, e.g., news texts from social media (Liu et al., 2019). These models and also LRNNs yield accurate results, but mainly only in the short run.

Conclusions
In this paper, we have introduced LRNNs. The major innovation in this work is a closedform approach to network size reduction (cf. Section 4.3) that learns both architecture and parameters of linearly activated RNNs. No backpropagation, gradient-descent, or other iterative procedure with several epochs is required, and the approach leads to significantly smaller, sparsely connected, and in many cases even minimal size networks.
We have shown that despite its simplicity of using only linear activation in the recurrent layer, LRNNs are a powerful approach to model time-dependent functions. The training procedure only uses standard matrix operations and is thus quite fast. In contrast to ESNs, also, no washout period is required. Any function can be approximated directly from its first step, with an arbitrary starting vector (cf. Property 5).
Although any time-dependent function can be approximated with arbitrary precision (cf. Property 7), not any function can be implemented by RNNs, in particular functions increasing faster than single-exponential (cf. Property 3) like 2 2 t (double-exponential) or t! (factorial function). Nevertheless, experiments with reasonably large example and network sizes can be performed successfully within seconds on standard hardware. However, if thousands of reservoir neurons are employed, the procedure may become numerically unstable, at least our Octave implementation. The likelihood of almost identical eigenvectors and eigenvalues with absolute values greater than 1 in the learned transition matrix W is increased then. Nonetheless, the underlying major problem here seems to be that existing scientific programming libraries do not calculate the eigenvalues of large matrices accurately enough.
A particularly interesting application of our approach reducing the network size is in hardware implementations of neural networks, e.g., for neuromorphic or reservoir computing (Mead, 1990;Indiveri et al., 2011;Liao and Li, 2017). Neuromorphic computing refers to new hardware that mimics the functioning of the human brain, and neuromorphic hardware results from the exploration of unconventional physical substrates and nonlinear phenomena. Future work shall include improving predictive and memory capacity of LRNNsanalyzed for small networks by Marzen (2017) and to some extent also by Couillet et al. (2016) -taking inspiration from convolutional networks (Goodfellow et al., 2016). Last but not least, other machine learning tasks besides prediction shall be addressed in more detail, including classification and reinforcement learning (Sutton and Barto, 2018;Pong et al., 2017).