Reservoir computing with the minimum description length principle

Reservoir computing is a branch of machine learning, in which only the responses of otherwise untrained dynamical systems are learned. Reservoir computing is attributed to Jaeger (2001) and Maass (2002) who initially introduced the concept as a simplified framework for training recurrent neural networks (Nakajima and Fischer, 2021). Since then, reservoir computing has maintained relevance at least in part because of its conceptual simplicity (Nakajima and Fischer, 2021) and low learning cost (Tanaka , 2019).

Despite simple training methods, reservoir computing has demonstrated great effectiveness in modeling temporal problems and has been used widely for the prediction of chaotic systems. In recent years, much work has been done to incorporate greater complexity in the readouts of reservoir computers while retaining the simplicity of their training. It has been shown by Bollt (2021) that reservoirs with linear activation functions may be improved by incorporating quadratic terms of reservoir component states. Quadratic terms are also used by Pathak (2018a) to build a readout for an ensemble of reservoirs and model spatio-temporal data. In Ren (2024), higher-order terms are included to construct a Hilbert space of polynomials of reservoir components.

While this recent work in echo-state network (ESN) readout augmentation has demonstrated the ability to significantly improve reservoir performance, it does so while substantially increasing readout size. We propose using the minimum description length principle (MDL) as a criterion for selecting subsets from these larger readouts. MDL is a powerful technique for comparison between models of different sizes that was first introduced by Rissanen (1978). We show that when combined with the subset selection algorithm presented by Judd and Mees (1995). MDL improves ESN performance. One possible explanation for such performance improvements is the reduction of linear dependence in the readout, which reduces the need for regularization.

We will first introduce reservoir computing with a particular focus on ESNs (Sec. II A), then we will introduce the MDL principle in Sec. II B. We look specifically at ESNs in this paper, though our implementation of MDL theory is applicable to all forms of reservoir computing including physical reservoirs. In Sec. II C, we provide methods for performance comparison between ESNs and present the dynamical systems that we use for prediction tasks. In Sec. III, we demonstrate the difference between MDL and ridge regression and compare their performance on benchmark tasks. In Sec. IV, we use the MDL principle to select components from a readout augmented with sub-graph states. Finally, in Sec. V, we provide an explanation in terms of linear dependence and consistency for the improved performance and reduced need for regularization in MDL readout subsets.

Reservoir computing is a method for time-series modeling with high-dimensional dynamical systems called reservoirs. These reservoirs are driven by an input signal, which they map onto a high-dimensional space of measurable components. After some transient period, these components become complex functions of historical input states if the echo-state property holds (Yildiz , 2012). If the echo-state property does hold, then the space of reservoir components is an embedding space (Duan , 2023), which converts a temporal problem into a spatial problem (Schrauwen , 2007). All that is needed to map the state of reservoir components back to the state of the input is a readout trained on pattern recognition (Tanaka , 2019). This readout consists of a weighted sum of reservoir components and is commonly trained using ordinary least squares regression or ridge regression (Sun , 2022). In theory, the latter helps prevent over-fitting, which is particularly important in reservoir computing since the reservoir state is high-dimensional by design.

In practice, an infinite transient period is not possible and the space of reservoir components may contain some amount of noise that is not functionally dependent on the input history. This can be demonstrated by driving two copies of a reservoir from different initial states with an identical time series and measuring the consistency capacity between them (Jüngling , 2021). Consistency capacity describes the effective dimension of the reservoir—the dimension of the subspace that is purely functionally dependent on the input history. A more detailed explanation of consistency capacity can be found in the  Appendix.

Echo-state networks (Jaeger, 2001), which pioneered the field of reservoir computing alongside the liquid state machine (Maass , 2002), remain one of the most commonly used reservoir frameworks. They can be conceptualized as consisting of a network with nodes ${ x i } i = 1 n$⁠, which are connected to a time-dependent input signal $u(t)$⁠. We denote the adjacency matrix of this network by $M$ and denote the vector of connections between the input and network nodes by $W$⁠. In order to introduce non-linearity, the network node states can be passed through an activation function $f: R→ R$⁠.

Minimum description length is derived from the idea of code-length. It is a method for comparing models of different sizes within the same model class by approximating their ability to compress data (Nakamura , 2003). The aim when building a model according to MDL theory is to minimize the sum of the code-lengths $L(⋅)$ of both the model parameters $L(θ)$ and the model errors $L(e |θ)$⁠.

However, an optimal model does not necessarily include all variables. An algorithmic approach that tries to optimize parameter subsets was given by Judd and Mees (1995). In this approach, which we will refer to as the breathing algorithm, parameter $θ i n$ is added to the optimal parameter subset of size $k$ increasing it to size $k+1$ and then the parameter with the smallest fitted coefficient $θ o u t$ is dropped. The added parameter $θ i n$ corresponds to the component with the strongest absolute correlation to the model errors made by prediction with the parameter subset of size $k$⁠. In other words, if $e k$ is the vector of errors obtained when predicting with the current set of parameters, then $θ i n$ corresponds to the component $X i n$ which maximizes $| X i n⋅ e k |$⁠. This is repeated until $θ i n$ and $θ o u t$ are the same parameter at which point $θ i n$ is included in the optimal subset and $k→k+1$⁠. After a model has been found for each subset size using this process, the description length for each is determined. The readout subset corresponding to the lowest description length is selected. This is illustrated in Fig. 1. Extensions to this algorithm have been proposed by Nakamura (2007).

Both the subset selection algorithm and the iteration of Newton’s method for optimizing coefficient precisions $δ i$ increase the computational complexity of the training process in most cases. If the subset selection algorithm is allowed to terminate early, then it is possible for the method detailed above to improve computational complexity by avoiding the inversion of large matrices. In our implementation using python code and an off-the-shelf laptop, ESNs with up to $1000$ nodes required less than a minute to train.

We use these metrics to compare ESN performance at modeling the standard Lorenz (⁠ $(σ,ρ,β)=(10,28, 8 3)$⁠), Rössler (⁠ $(a,b,c)=(0.1,0.1,14)$⁠), and Thomas (⁠ $b=0.208186$⁠) systems over a range of forecasting horizons. The equations for these systems are provided in the  Appendix. Since the characteristic timescales of each system are different, forecasting horizons are measured in Lyapunov time which is the inverse of the leading Lyapunov exponent. We estimate the leading Lyapunov exponents of the Lorenz, Rössler, and Thomas systems with the coefficients above to be roughly $0.9121±0.0044$⁠, $0.0660±0.0068$⁠, and $0.0378±0.0040$⁠, respectively. Accepted values for these systems are around $0.9057$ for Lorenz (Viswanath, 1998) and $0.0714$ for Rössler (Sprott, 2023). Our calculations are made using the method described by Sandri (1996) and the provided error margins are three times the standard deviation of bootstrap sample means. The Lyapunov times (⁠ $T L$⁠) we use are $1.096$⁠, $15.14$⁠, and $26.48$⁠. Given the difference between the timescale of the Lorenz attractor and those of Rössler and Thomas, we use $δt=0.02 T L$ for each system (This equates to $0.02$s for Lorenz, $0.3$s for Rössler, and $0.53$s for Thomas). A time delay embedding of each attractor is illustrated in Fig. 2.

In our experiments, we use the hyper-parameters outlined below:

1. Network size: We use $200$-node ESNs unless otherwise specified.

2. Activation function: We use the hyperbolic tangent $tanh(⋅)$ activation function unless specified.

3. Spectral radius: $ρ=1$ or $ρ=0.9$ when linear activation is used.

4. Input weight matrix: We draw input weights from a zero mean Gaussian distribution $W∼N( 0 →,0.3 I)$⁠.

5. Adjacency matrix: We create the adjacency $M$ matrix by first taking the adjacency matrix $A$ of a Watts–Strogatz random network with average degree $6$ and re-wire probability $10 network size$⁠. We then replace non-zero entries in $A$ with independent samples from a zero-mean normal distribution and then scale the resulting matrix by its spectral radius.

6. Time interval between inputs: $δt=0.02 T L$⁠, where $T L$ is the Lyapunov time of the driving system.

7. Testing and training data: We use $5000$ points from each time series to train the ESNs. They are then tested on a new time series from the same system with 25 000 points.

8. Leakage: $β=0$⁠. Non-zero leakage is beneficial for longer forecasting horizons. However, we avoid using leakage for simplicity.

Typically in reservoir computing, ridge regression is used for regularization to reduce model over-fitting. The penalty term introduced by ridge regression reduces the magnitude of coefficients overall and improves numerical stability in the case of perfect or imperfect multicollinearity. However, ridge regression includes all model parameters and does not encourage sparsity as strongly as alternatives like Lasso regression (Tibshirani, 1996). There is a slightly different rationale behind selecting readout subsets using MDL and the breathing algorithm described in Sec. II B. Reservoir components that have little influence on model performance are completely neglected. Multicollinearity in the matrix of output states $X$ is avoided by dropping components that are highly correlated to any others already contained within the optimal component subset. When ridge regression is used in conjunction with MDL subset selection, smaller ridge penalties $α$ tend to be required since a reduction in linear dependence in the design matrix $X$ results in a better-conditioned matrix inverse $( X T X + α I ) − 1$⁠.

MDL also encourages sparsity by limiting readout size. The size of the readout chosen using MDL is bounded above according to the amount of data used to train the model. This is because a model cannot achieve a minimum total description length if its parameters require the specification of more information than the target dataset. In reality, the optimal subset size will fall well below this limit for a good model.

If we compare this measurement to the size of the ESN, we find that the effective readout size continues to increase with the readout size when ridge regression is used on its own (Fig. 3). Perhaps counter-intuitively, readouts trained with weaker regularization (smaller $α$⁠) appear more sparse. However, this is because the matrix inverse $( X T X + α I ) − 1$ may become poorly conditioned for small $α$ causing some coefficients to blow up and account for a large part of the norm of the coefficient vector. On the other hand, the size of the MDL readout is substantially smaller than that produced by ridge regression for large ESNs. The size of the MDL readout reaches a plateau for some sufficient ESN size.

Despite the sparsity of the MDL trained readout, it performs better than ridge regression at the task of $n$-step ahead prediction for the Lorenz, Rössler, and Thomas systems. In the experiment below (Fig. 4), we compare the mean square prediction error of ESNs trained using either the entire readout (shown with dashed lines) or MDL readout subsets (shown with solid lines). MSE was found for three different forecasting horizons across a range of ridge regression penalty sizes. Results were averaged across $30$ trials and error bars indicate the $5th$ and $95th$ percentile of bootstrap sample means. We find that the benefit of MDL subset selection is most apparent on short forecasting horizons. This is likely because model over-fitting is easier on the comparatively simpler task of short-term forecasting. This would also explain why stronger regularization (larger $α$⁠) is also needed for shorter horizons.

Improvements from MDL subset selection at $n$-step ahead prediction are also apparent with $1000$ node ESNs (Fig. 5) but not on small (⁠ $50$-node) ESNs (Fig. 6). It is not surprising that MDL provides no benefit to small reservoirs since the size of the optimal subset according to MDL is the same as the size of the full readout (Fig. 3) for ESNs with less than around $100$ nodes. However, performance is better in the larger (⁠ $200$ or $1000$⁠)-node ESNs where MDL readout subset selection is beneficial.

We also test the robustness to noise of MDL trained readouts in Fig. 7 and their ability for reconstruction in Fig. 8. Robustness to noise is similar in ESNs trained with and without MDL subset selection. Improved performance in the MDL readouts with low levels of noise and short forecasting horizons just follows from results in Fig. 4. Performance in ESNs that use the entire readout is improved by low levels of Gaussian noise. This is because low levels of Gaussian noise reduce linear dependence in $X$ and result in a better-conditioned inverse $( X T X + α I ) − 1$⁠. However, performance is not improved by low levels of Gaussian noise (with variance above $10 − 6$⁠) when MDL subsets are selected since linear dependence in $X$ is reduced to begin with.

In autonomous reconstruction tasks, like in $n$-step ahead prediction tasks, ESNs trained using MDL readout subsets produce a similar performance to those trained using the whole readout when ridge penalties are large (Fig. 8). When penalties are smaller, performance appears to be qualitatively better with MDL subsets. However, MDL is arguably less relevant for autonomous reconstruction since ideal ridge parameters are larger.

Much work has been done lately to improve the performance of ESNs by incorporating higher-order products of reservoir components in their otherwise linear readouts. The performance improvements of this readout augmentation are particularly striking when linear activation functions are used in place of more conventional non-linear activation functions like $tanh⁡(⋅)$⁠. As demonstrated by Bollt (2021), ESNs with a linear activation function become Vector Auto-Regressive (VAR) models whose performances may be improved with quadratic terms. These quadratic terms are introduced by taking the product of pairs of network nodes. In general, it has been shown that by introducing non-linearity in just the readout of an otherwise linear dynamical system, a wide range of signals may be approximated (Boyd and Chua, 1985). It has also been shown that quadratic terms can help to break harmful symmetries that lead to mirror attractors in autonomous reconstruction tasks (Herteux and Räth, 2020).

If we build from the quadratic readout in Bolt’s paper and include network sub-graphs of arbitrary size by taking a direct product of sub-graph nodes, then this VAR becomes a polynomial model.

Readout weights are determined for augmented readouts as they otherwise would be except that columns representing sub-graph states are added to the design matrix $X$⁠. For example, if a triangle involving nodes $1$⁠, $3$⁠, and $14$ is added to the readout of a $200$-node ESN, then the entry $X 201 , t$ would take the value $X 1 , t⋅ X 3 , t⋅ X 14 , t$⁠.

In Fig. 9, we demonstrate that the inclusion of sub-graph states can improve performance substantially when a linear activation function is used. These improvements are perhaps unsurprising since without higher-order sub-graph states, the reservoir output is a purely linear combination of historical inputs and linear approximations of dynamical systems are only suitable over very short time horizons. However, if the whole readout is trained and small ridge penalties are used performance decreases with readout augmentation.

It is also true that performance is improved by readout augmentation for MDL readout subsets with a non-linear (hyperbolic-tangent) activation function. Performance improvements occur from augmentation but are smaller when the entire readout is trained. This is demonstrated in Fig. 10 where tests performed in Fig. 4 are repeated with and without edge, triangle and four-cycle states included in the readout.

In Figs. 11 and 12, we re-perform the noise robustness and autonomous reconstruction experiments from Sec. III. The best performance in noise robustness tests using the whole readout occurs with a larger amount of noise when sub-graph states are included. This is likely because noise has a regularizing effect and more regularization is needed for the larger readout. Augmentation appears to reduce robustness to high degrees of Gaussian noise for ESNs trained using MDL readout subsets.

Reconstruction performance with and without augmentation is harder to compare due to its qualitative nature. However, the absence of dark bands off the attractor when $α= 10 − 12$ indicates an improvement from Fig. 8. Again, no successful reconstruction was achieved using the entire readout with $α= 10 − 16$⁠.

An important finding from Secs. III and IV is that MDL readout subsets require less regularization. For example, we can see from $n$-step ahead prediction and autonomous reconstruction performance that large ridge penalties are needed when training the entire readout. However, performance on both tasks remains reasonable when smaller ridge penalties are used to train MDL readout subsets. Regularization from low levels of added Gaussian noise is beneficial for training the full readout but not for training MDL readout subsets.

Regularization in the form of a ridge regression penalty or added Gaussian noise helps stabilize the matrix inverse $( X T X + α I ) − 1$ if the columns of $X$ [which are the time series of node or sub-graph states $( x i( t 0),… x i( t T))$] are linearly dependent. We argued in Sec. III that this linear dependence is often avoided when using MDL readout subsets. We can test this theory with consistency capacity (see Subsection 3 of the  Appendix), which measures the number of linearly independent components of the driving signal contained within the reservoir (Lymburn , 2019).

The consistency spectrum, different to the consistency capacity, is a ranked set of reservoir component consistencies. The consistency of each reservoir component indicates its degree of functional dependence on the driving signal and their sum is the consistency capacity. By looking at the consistency spectrum in Fig. 13, we can see that additional independent signal components are provided by the higher-order sub-graph states.

We can also measure an effective consistency capacity for readout subsets (including or not including higher-order states) by following the calculation in Subsection 3 of the  Appendix but including only the columns of $X$ and $X ′$ that correspond to elements of the readout. Rather than being a measure of the number of linearly independent signal components contained within the entire reservoir state, this effective consistency capacity measures the number of such components contained within the sub-state of the reservoir that gets mapped to the output. If coefficients are only fitted for a subset of readout components, then it is this effective consistency that matters. We normalize the effective consistency capacity of a readout subset by the subset size to provide a measure of the proportion of the readout that corresponds to linearly independent components of the driving signal.

We demonstrate in Fig. 14 that using the breathing algorithm with MDL as a subset selection criterion improves readout consistency. This means that the proportion of readout components in the MDL readout subset that are linearly independent and functionally dependent on the driving signal is higher than in the whole readout.

One mechanism for this improvement is the reduction in readout size. Smaller readouts tend to be more consistent even with random subset selection. In theory, the optimal subset size according to MDL theory should occur while consistency is still high. Once a subset size is reached where new components are less relevant to the signal or similar to those already contained, any reduction in $L(e)$ is less likely to be worth the increase in $L(θ)$⁠.

The breathing algorithm also tends to select component subsets with disproportionately high consistency capacity when readouts are small. This is because, at each subset size, the algorithm chooses the predictor $θ in$ most strongly related to the prediction error. With large enough subsets, the breathing algorithm may eventually choose less consistent readouts than random sampling—as is the case when the readout is augmented in Fig. 14. However, this is only the case for readout sizes larger than the optimal size according to MDL and occurs because there may be signal components that are consistent but not useful. Irrelevant signal components will not be related to prediction error.

In this paper, we demonstrated the use of minimum description length in reservoir training. In doing so, we found that using MDL readout subsets is beneficial for $n$-step ahead prediction tasks and results in qualitatively similar attractor reconstructions. One explanation for the advantages of MDL readout subsets is that they contain less linear dependence which we demonstrated using consistency capacity. This reduction in linear dependence also made MDL readout subsets less sensitive to the size of the ridge penalty.

We argued that decreased linear dependence in MDL readout subsets was only partly explained by their smaller size. The other explanation was that the breathing algorithm selected subsets with higher consistency than random sampling.

The benefits of using the MDL principle in ESN training that we found did not apply to small (⁠ $50$-node) ESNs. However, these small ESNs produced inferior performance to their larger counterparts which were improved by MDL subset selection.

As such, we believe that our findings provide a strong case for the consideration of MDL in practical reservoir training tasks and that this case is strengthened by recent attention to readout augmentation.

A.M. was supported by the Australian Government RTP scholarship at the University of Western Australia. M.S. and D.M.W. are supported by the ARC Discovery Grant (No. DP200102961), funded by the Australian Government. M.S. also acknowledges the support of the Australian Research Council through the Centre for Transforming Maintenance through Data Science (Grant No. IC180100030).

The authors have no conflicts to disclose.

Antony Mizzi: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal). Michael Small: Conceptualization (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). David M. Walker: Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

The data that support the findings of this study are openly available in MDL_ESNs, at https://github.com/antony-mizzi/MDL_ESNs.git (Mizzi, 2025).

We generate time series from the Lorenz, Rössler, and Thomas systems using an Eulerian integration with an interval $dt=0.002$ for Lorenz and $dt=0.02$ for Rössler and Thomas. We then sample points at a rate of $δt=0.02$ for Lorenz and $δt=0.2$ for Rössler and Thomas. The set of sampled points is used for reservoir prediction tasks. The update equations for each system are given below:

Allow us to denote $X$ to be the state matrix with rows $x i(t)$ that contain the state of the $ith$ node at times $t∈{0,δt,…,(T−1)δt}$⁠. We require $X$ to be normalized such that $x i(t) x i T(t)= x j(t) x j T(t)∀i,j$⁠. If $B⊂{1,2,..,N}$ is a node subset then $X B$ is the sub-matrix with rows corresponding to nodes in $B$⁠. Initially, we set $B=∅$⁠.

1. Fit coefficients $θ B= ( X B X B T + α I ) − 1 X B TY$ for nodes in the subset $B$ which has size $k$⁠. Denote the errors of this fit as $e B=Y− X B Tθ$⁠.

2. Pick the row $i$ that corresponds to the maximum entry of $X T e B$⁠.

3. Fit coefficients $θ B ′$ for the matrix $X B ′$ with rows corresponding to nodes in $B ′=B∪{i}$⁠. Find the parameter $j$ with the smallest absolute entry in $θ B ′$⁠.

4. If $i=j$ (the newly included row has the smallest absolute parameter) then include $i$ in the optimal parameter subset (⁠ $B→B∪{i}$⁠, $k→k+1$⁠).

5. If $i≠j$ then repeat the process using $B→{j}∪B∩ { i } ¯$⁠.

Consistency capacity is a measure of the functional dependence of the reservoir state on the input signal. In fully consistent reservoirs, the echo-state property holds completely. Consistency capacity is measured as follows.