Minimising Stochastic Complexity with Ridge Regression

1. Introduction

We approach the problem of selecting a regularisation strength for ridge regression with the intent of producing the simplest generalised linear model (GLM). This complexity-based approach is particularly relevant since the popularity of GLMs stems, in part, from their simplicity as a model class [1,2,14,22,43]. They are often used to identify and characterise relationships from data in fields like medicine [6], dietary science [7], economics [8,9] and sociology [10,11] . Such applications lie in accordance with Occam’s razor principle, which tells us to seek simple explanations, or in this case, simple model classes.

In these endeavours, it is common to use some form of regularisation in the fitting process, like an L1 ridge penalty. In the case of ridge regularisation, the penalty term stabilises the matrix inverse under multicollinearity of input variables [12,13,14] and improves generalisation to new data by reducing model variance at the expense of increased model bias [15,16,17]. The application of regularisation also further reduces the complexity of linear models. This makes it possible to choose the regularisation strength that minimises model complexity.

In this paper, we measure complexity in terms of stochastic complexity [18], which is a refined variant of the minimum description length principle [19] (MDL). The stochastic complexity of a model is directly related to its flexibility. When used for model selection, it suggests that we should hesitate to choose a model class that fits our dataset well if it also provides good fits to unrelated datasets using the same predictors. While regularisation constrains flexibility and produces simpler models, this is not directly accounted for by stochastic complexity. As such, the purpose of this paper, as well as deriving a ridge strength criterion, is to demonstrate a method for determining the effect of regularisation on the stochastic complexity of a model. Our method makes use of a novel interpretation of stochastic complexity’s underlying coding scheme in order to directly account for the constraint imposed by regularisation on model flexibility.

Our work is organised as follows: In Section 2, we provide an overview of stochastic complexity, including an explanation of the typical underlying coding scheme. In Section 3 we discuss our method for modifying this scheme to account for regularisation. In Section 4 we derive an analytical form for the stochastic complexity of a linear model fitted by ridge regression using our modified scheme. We also present our ridge penalty selection criterion. In Section 5 we compare our criterion with several other information-theoretic criteria. In particular, we note that our work is closely related to that of Dwivedi et al [20] who recently used luckiness normalised maximum likelihood [21] (LNML) for ridge strength selection. Their approach produces a different selection criterion, which stems from a different interpretation of regularisation’s effect on complexity. In Section 6 we use our criteria to train linear readouts with the example of reservoir computing. In Section 7, we conclude with a discussion and summary.

2. Background into Stochastic Complexity and the Minimum Description Length Principle

The premise behind MDL is that models compress information. This is done by making use of a parametric model class to derive a coding scheme for possible data sets. The shorter the encoding of the dataset we wish to model under the scheme we derive, the better the compression achieved by the model. For example, in classical two-part MDL, the scheme is to first encode our model parameters and then to encode the errors that the model makes using these parameters. The best models, according to two-part MDL, minimise the combined lengths of their error and parameter codes.

The coding scheme assumed by stochastic complexity is more abstract. No one set of model parameters is used. Instead, a probability $P^{*}(Y;\mathcal{M},X)$ is assigned to each potential dataset Y from the set of possibilities $Y\in \mathcal{Y}$ using a weighted ensemble of parameterisations from the same model class $\mathcal{M}:(X,\theta )\to Y$. If we denote the weighting associated with each parametrisation $\theta$ in the ensemble as $\pi \left(\theta \right)$, then the probability assigned by the resulting ensemble would be equivalent to the marginal probability under $\pi$ as a prior

$$P^{*}(Y;\mathcal{M},X)=\sum _{\theta }P\left(Y\right|\theta ;\mathcal{M},X)\pi \left(\theta \right).$$

Importantly, the ensemble weightings are chosen deterministically based on the model class. This is what enables the resulting code to be deciphered and makes stochastic complexity a valid description length. Under the stochastic complexity scheme, ensemble weights $\pi \left(\theta \right)$ are chosen before the dataset Y is known in the way that minimises regret for the worst possible scenario of Y. This regret $R(Y,\pi )$ is the difference between the length of code assigned to Y by the weighted ensemble (which we denote by $L_{\mathrm{sc}}\left(Y\right)$) and the shortest length that could be used to encode Y with any other ensemble. Under the lower bound for efficiency given by Shannon’s source coding theorem, these code lengths are negative log probabilities $L_{\mathrm{sc}}(Y;\mathcal{M})=-log\left(P^{*}\left(Y\right)\right)$. Since the best ensemble for encoding any particular dataset Y is just its maximum likelihood estimate ${\widehat{\theta }}_Y$, the worst-case regret which we seek to minimise by choice of weightings is the following

$$R_0\left(Y\right):=\underset{Y\in \mathcal{Y}}{max}\{-log\left(P^{*}(Y;\mathcal{M},X)\right)+log\left(P\left(Y\right|{\widehat{\theta }}_Y;\mathcal{M},X)\right)\}.$$

The weightings that minimise this worst-case regret are those that induce proportionality to the maximum likelihood estimate ($P^{*}(Y;\mathcal{M})={\displaystyle \frac{1}{c}}P\left(Y\right|{\widehat{\theta }}_Y;\mathcal{M})$). This means that we can work directly with $P^{*}$ without ever actually determining the optimal weightings. Under proportionality, regret is constant across all datasets $R(Y,\pi )=R(Y^{\prime },\pi )=-log\left(c\right)$. The constant of proportionality c must be the sum of all maximum likelihood estimates in order for it to induce a probability distribution that integrates to one. The probability of Y under $P^{*}$ is therefore given as follows

$$P^{*}(Y;\mathcal{M})={\displaystyle \frac{P\left(Y\right|{\theta }_Y^{*};\mathcal{M},X)}{{\sum }_{Y^{\prime }\in \mathcal{Y}}P\left(Y^{\prime }\right|{\theta }_{Y^{\prime }}^{*};\mathcal{M},X)}}.$$

When we work with continuous data, and therefore an infinite set of possible datasets $\mathcal{Y}$, the ensemble weights $\pi \left(\theta \right)$ become difficult to determine. They have been shown [22,23], to approach the Jeffreys prior [24] for large datasets under certain conditions [25]. However, we are not required to determine the ensemble weights since proportionality still holds, and we can again work directly with $P^{*}$. In the case of continuous data, the stochastic complexity makes use of the integral version for $P^{*}$

$$L_{\mathrm{SC}}(Y;\mathcal{M})=-log{\displaystyle \frac{P\left(Y\right|{\theta }_Y^{*};\mathcal{M},X)}{{\int }_{Y^{\prime }\in \mathcal{Y}}P\left(Y^{\prime }\right|{\theta }_{Y^{\prime }}^{*};\mathcal{M},X)d{Y}^{\prime }}}.$$

(1)

Notice that stochastic complexity can be broken into two terms, the first being the model’s negative log likelihood. The other term measures the complexity of the parametric model class $\mathcal{M}$

$$\mathrm{Penalty}=log\left({\int }_{Y^{\prime }\in \mathcal{Y}}P\left(Y^{\prime }\right|{\theta }_{Y^{\prime }}^{*};\mathcal{M},X)d{Y}^{\prime }\right),$$

and functions as the penalty term in model selection. This penalty directly measures flexibility and it lends stochastic complexity an important intuition beyond its meaning as a code length. It tells us not to trust models that fit the data well if they can also provide good fits to unrelated data using the same predictors.

3. The effect of Regularisation on Complexity

The role of model regularisation on stochastic complexity can potentially be interpreted in different ways. Our interpretation is to directly consider the effect of regularisation on the flexibility of the model class. This follows on from the intuitive meaning of stochastic complexity described in the previous section: that we should have more confidence in our regularised model if the regularisation applied prevents the model from easily fitting unrelated datasets using the same predictors.

Rather than building a weighted ensemble of maximum likelihood estimates, we build a weighted ensemble of regularised fits. The weightings of this ensemble will depend on the regularisation strength $\alpha$. We select the weightings such that the probability assigned by the ensemble $P_{\alpha }^{*}(Y;\mathcal{M})$ is proportional to the probability assigned by the regularised fit $P\left(Y\right|{\theta }_{Y,\alpha }^{*};\mathcal{M})$, which again minimises worst-case regret. This regret, which we now denote by $R_{\alpha }\left(Y\right)$, because of its dependence on the regularisation strength $\alpha$, is still the length of code that could have been saved had the parameter fit for the dataset been known ahead of time. The only difference is that this parameter fit is now the set of regularised parameters and not the maximum likelihood estimate

$$R_{\alpha }:=\underset{Y\in \mathcal{Y}}{max}\{-log\left(P_{\alpha }^{*}(Y;\mathcal{M},X)\right)+log\left(P\left(Y\right|{\widehat{\theta }}_{Y,\alpha };\mathcal{M},X)\right)\}.$$

Assuming that the encoder and decoder also know the regularisation strength ahead of time, the resulting ensemble weights can be calculated by the decoder, allowing them to decipher any subsequent encoding. Under this interpretation, which therefore describes a valid description length, the stochastic complexity of a dataset Y becomes the following

$$L_{\mathrm{SC}}(Y;\mathcal{M},\alpha )=-log{\displaystyle \frac{P\left(Y\right|{\theta }_{Y,\alpha }^{*};\mathcal{M},X)}{{\int }_{Y^{\prime }\in \mathcal{Y}}P\left(Y^{\prime }\right|{\theta }_{Y^{\prime },\alpha }^{*};\mathcal{M},X)d{Y}^{\prime }}}.$$

(2)

One benefit of this approach is its generalisability. It can be applied directly to other methods for fitting and regularising parameters. Another benefit is the fairness of the comparison it draws between regularisation techniques. One technique can be compared with another based purely on how much each restricts the flexibility of the model class. The caveat, which is not unique to our approach, is that to derive an analytical value for stochastic complexity, we still need an analytical form for the model likelihood when regularisation is applied.

4. Deriving a Ridge Penalty to Minimise Stochastic Complexity

In the case of regression with a known state matrix X, the parametric model class $\mathcal{M}$ is the multivariate normal distribution $\mathcal{N}(X\theta ,{ϵ}^2\mathbb{I})$. This is the model class for which the regression solution $\theta ={\widehat{\theta }}_Y$ provides the maximum likelihood estimate. From a Bayesian perspective, the ridge penalty comes from maximising the posterior with the added assumption of a normal prior over the parameters. However, it is not strictly necessary to assume a prior over the parameters in order to calculate stochastic complexity. We only need to know the probability of each dataset $Y\in \mathcal{Y}$ given the regularised parameter fit ${\theta }_{Y,\alpha }^{*}$. According to our model class, that probability is the following

$$P\left(Y\right|{\theta }_{Y,\alpha }^{*};\mathcal{M},X)={\left(2\pi ϵ\right)}^{-n/2}\exp \left(-{\displaystyle \frac{\left|\right|Y-\mathcal{M}(X,{\theta }_{Y,\alpha }^{*}){\left|\right|}^2}{2{ϵ}^2}}\right).$$

Fortunately, an analytical form for this probability exists. First, we can express the ridge regression solution ${\theta }_{Y\alpha }^{*}$ in terms of the penalty $\alpha$, the state matrix X and the dataset Y. It is given by ${\theta }_{Y\alpha }^{*}={(X^{T}X+\alpha \mathbb{I})}^{-1}{X}^{T}Y$. This means that an analytical form for the mean square error, which we denote by ${\sigma }_{Y,\alpha }^2$, also exists

$$\begin{array}{cc}\hfill {\sigma }_{Y,\alpha }^2={\displaystyle \frac{1}{n}}{Y}^T{Q}_{\alpha }Y,\mathrm{where}Q& ={(P_{\alpha }-\mathbb{I})}^T(P_{\alpha }-\mathbb{I})\hfill \\ \hfill \mathrm{and}{P}_{\alpha }& =X{(X^{T}X+\alpha \mathbb{I})}^{-1}{X}^T.\hfill \end{array}$$

In the expression above, n refers to the length of Y.

To find the stochastic complexity penalty term, we integrate the resulting probability ${\left(2\pi ϵ\right)}^{-n/2}\exp (-{\displaystyle \frac{1}{2{ϵ}^2}}{Y}^T{Q}_{\alpha }Y)$ over the space of all possible datasets. For continuous datasets of length n, this space is ${\mathbb{R}}^n$ and the penalty is the following

$$log{\int }_{{\mathbb{R}}^n}{\left(2\pi {ϵ}^2\right)}^{-n/2}exp(-{\displaystyle \frac{Y^T{Q}_{\alpha }Y}{2{ϵ}^2}})dY.$$

We can write the value of this integral in terms of the determinant of $Q_{\alpha }$

$$-{\displaystyle \frac{1}{2}}log(det\left(Q_{\alpha }\right)),$$

or in terms of the eigenvalues of the design matrix $X^{T}X$, which we denote by ${\lambda }_i$

$$-{\displaystyle \frac{1}{2}}log(det\left(Q_{\alpha }\right))={\displaystyle \sum _{i=1}^m}log(1+{\displaystyle \frac{{\lambda }_i}{\alpha }}).$$

This penalty describes how well we can fit various possible datasets using regression with a ridge constant of $\alpha$ and a gram matrix $X^{T}X$ that has the given spectrum $\{{\lambda }_1,{\lambda }_2,...{\lambda }_m\}$. Lastly, we need the other component of the stochastic complexity, which is the log likelihood of Y under its own parameter fit $-logP\left(Y\right|{\theta }_{Y,\alpha }^{*};\mathcal{M})$. Together, the stochastic complexity becomes the following function of the ridge penalty:

$$L_{\mathrm{SC}}(Y;\alpha ,X)={\displaystyle \frac{n}{2{ϵ}^2}}{\sigma }_{Y,\alpha }^2+{\displaystyle \sum _{i=1}^m}log\left(1+{\displaystyle \frac{{\lambda }_i}{\alpha }}\right)+\mathrm{const}.$$

(3)

The constant term, which we ignore because it is irrelevant for our purpose, is $nlog\left(\sqrt{2\pi ϵ}\right)$.

Calculating the stochastic complexity requires choosing a variance ${ϵ}^2$ for the multivariate Gaussian used in our coding scheme. This is also the case for various other ridge penalty criteria [26,27]. A typical method for doing so in the under-parameterised case ($m<<n$) is to estimate $ϵ$ in terms of the adjusted mean square error of the linear regression solution ${ϵ}^2=n{\sigma }_{Y,0}^2/(n-m)$. A more robust method, which is not limited to the under-parametrised case, is provided by Liu et al [28]. They suggest approximating ${ϵ}^2$ by the following

$${ϵ}^2={\displaystyle \frac{Y^T(\mathbb{I}-P_{\alpha })Y}{n-\mathrm{Tr}\left(P_{\alpha }\right)}}.$$

Our formulation for optimising the ridge penalty passes two sanity checks that we can construct from our understanding of ridge regression. Firstly, the optimal penalty scales correctly with the terms of the state matrix. If X is multiplied by a scalar $X\to c*X$, the optimal ridge penalty should scale with its square $\alpha \to c^2\alpha$, since this leaves the projection matrix $X{(X^{T}X-\alpha \mathbb{I})}^{-1}{X}^T$ unchanged. This is the case. Multiplying the state matrix X by a scalar adds a trivial constant to the penalty term and maps ${\sigma }_{Y,\alpha }^2$ to ${\sigma }_{Y,c^2\alpha }^2$. Secondly, since the penalty term is a measure of model flexibility, it should only depend on the unique, non-zero columns in X. This is also true. Duplicate columns, or columns of zeros, add eigenvalues of zero to the spectrum and these do not contribute to the penalty term.

5. Alternative Approaches

In Section 3 we accounted for regularisation by devising a coding scheme from an ensemble of regularised parameter fits. An alternative approach is to build the ensemble from maximum posterior likelihood estimates after assuming a prior distribution over parameter space. This is similar to the way that luckiness normalised maximum likelihood (LNML) assumes a so-called `luckiness function’ $p_{\mathrm{luck}}\left(\theta \right)$ to derive a model score

$$\mathrm{LNML}\left(Y\right|\mathcal{M}):={\displaystyle \frac{{max}_{\theta }\left[P\left(Y\right|\theta ;\mathcal{M})p_{\mathrm{luck}}\left(\theta \right)\right]}{{\int }_{Y^{\prime }\in \mathcal{Y}}{max}_{\theta }\left[P\left(Y^{\prime }\right|\theta ;\mathcal{M})p_{\mathrm{luck}}\left(\theta \right)\right]d{Y}^{\prime }}}.$$

(4)

Dwivedi et al [20] chose $p_{\mathrm{luck}}$ to be a Gaussian density for ridge regression

$$p_{\mathrm{luck}}\left(\theta \right)\propto exp(-{\displaystyle \frac{\alpha }{2{ϵ}^2}}{\theta }^T\theta ).$$

This choice is meaningful because a Gaussian prior over the parameters produces a posterior likelihood which is maximal at the ridge regression solution. The numerator and denominator therefore both measure posterior likelihoods. The corresponding description length assigned by Dwivedi et al, which we denote by $L_{\mathrm{Dwivedi}}$ and provide below, is similar to our own, but includes an additional log parameter likelihood term

$$L_{\mathrm{Dwivedi}}(\alpha ;X,Y)={\displaystyle \frac{n}{2{ϵ}^2}}{\sigma }_{Y,\alpha }^2+{\displaystyle \frac{\alpha {\left({\theta }_{Y,\alpha }^{*}\right)}^T{\theta }_{Y,\alpha }^{*}}{2{ϵ}^2}}+{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^m}log(1+{\displaystyle \frac{{\lambda }_i}{\alpha }}).$$

The approach taken by Dwivedi et al can be generalised to any regularisation techniques that offer a Bayesian interpretation. The luckiness score would simply be replaced with whatever prior distribution is implicitly assumed over the parameters.

The MDL-inspired approach taken by Silhavy et al [29] makes use of a two-part coding scheme. They take the description length found for linear regression by Giurcaneanu et al [30] and consider the additional cost of encoding the ridge parameter (or Lasso parameter). It should be noted, that Giurcaneanu et al found the description length for linear regression under a different framing of the model class. They encode the error variance ${\sigma }_{Y,\alpha =0}^2$ and then use the variance to encode the errors $Y\sim \mathcal{N}(X\theta ,{\sigma }_{Y,\alpha =0}^2\mathbb{I})$. As such, Silhavy et al’s method for choosing an MDL ridge parameter is quite different from our own. They seek to minimise the following criterion

$$L_{\mathrm{Silhavy}}(\alpha ;X,Y):={\displaystyle \frac{n-k}{2}}log\left({\sigma }_{Y,\alpha }^2\right)+{\displaystyle \frac{1}{2}}log\left(\alpha \right)+{\displaystyle \frac{k}{2}}log\left({\theta }_{\alpha }^T{X}^{T}X{\theta }_{\alpha }\right)+\mathrm{const}.$$

(5)

It should also be noted that the major objective of Silhavy et al’s work was to apply MDL to multi-criteria decision analysis. This means seeking a tradeoff between model complexity and a combination of various performance measures. In this paper, we only deal with a single objective, which is model log likelihood (effectively mean square model error).

Other information criteria, like the Akaike [31] and Bayesian [32] Information Criteria (AIC and BIC respectively), can also be used to choose a ridge strength if the number of model parameters is replaced with the number of effective degrees of freedom which has dependence on $\alpha$ [33,34]. The commonly used measure of effective degrees of freedom is the trace of the projection matrix $P_{\alpha }=X{(X^{T}X+\alpha \mathbb{I})}^{-1}{X}^T$ [35,36]. This results in the modified AIC and BIC criteria provided below

$${\mathrm{AIC}}_{\mathrm{eff}}(\alpha ;X,Y)={\displaystyle \frac{n}{{ϵ}^2}}{\sigma }_{Y,\alpha }^2+2{\displaystyle \sum _{i=1}^m}(1+{\displaystyle \frac{\alpha }{{\lambda }_i}}),\mathrm{and}$$

$${\mathrm{BIC}}_{\mathrm{eff}}(\alpha ;X,Y)={\displaystyle \frac{n}{{ϵ}^2}}{\sigma }_{Y,\alpha }^2+log\left(n\right){\displaystyle \sum _{i=1}^m}(1+{\displaystyle \frac{\alpha }{{\lambda }_i}}).$$

6. Fitting a Reservoir Computer

As a proof of concept, we use our ridge penalty criterion to select a penalty strength for fitting reservoir computer readouts. Reservoir computers are high-dimensional dynamical systems that become conditionally stable when driven by an input signal. Their dynamic response is mapped to the future state of the input with a linear readout. This mapping is possible because the reservoir state-space becomes a time-delay embedding-space of the driving time series [37] according to the echo-state property [38,39], which reservoirs must satisfy. Reservoir readouts are often fitted with ridge regression [40,41,42] and their complexity is a subject of interest [43,44,45,46]. In particular, we use echo-state networks [47] as our reservoirs and perform time series prediction on Lorenz-96 [48]. Our exact method is detailed in the Appendix along with a description of echo-state networks.

The average generalisation error for one-step prediction of Lorenz-96 is compared with network size in Figure 1. The ridge penalty chosen according to stochastic complexity (in black) produces smaller generalisation error on average than any fixed penalty magnitude. This is true across a range of network sizes and suggests that our method leads to sensible choices for the ridge parameter. Our criterion also achieves comparable results to some of the alternative approaches mentioned in Section 5, and performs comparatively well on large networks. This is demonstrated in Figure 2. It should be noted that the relative strength of the model complexity term is dependent on the value of $ϵ$ for these methods, which we choose to be the smaller of the two values described in Section 4

$${ϵ}^2=min({\displaystyle \frac{n{\sigma }_{Y,0}^2}{n-k}},{\displaystyle \frac{Y^T(\mathbb{I}-P_{\alpha })Y}{n-\mathrm{Tr}\left(P_{\alpha }\right)}}).$$

The first of these values is always positive and defined, since the length of the training time series is 1200 points and the largest network has 1096 nodes.

Qualitative comparisons of reconstruction ability are also provided in Figure 3 and Figure 4. This time, echo-state networks are trained for one-step-ahead prediction of the Lorenz system. Individual trajectories in each subplot are generated by continuously driving echo-state networks with their own past predictions. These predictions depend on output coefficients, which in turn depend on the ridge strength selected by each criterion. All criteria tend to produce reconstructions that are qualitatively similar to the Lorenz attractor when networks and prediction step sizes are small. This indicates that the ridge strengths they select tend to be appropriate for capturing the underlying dynamics. For example, the reconstruction composites in Figure 3 are made using 200 node networks trained at a sampling rate of $0.02$s. On the other hand, ridge strengths selected by stochastic complexity appear to produce reasonable reconstructions more often when networks and step sizes are larger. This is the case in Figure 4 where predictions are made using 800 node networks and a step size of $0.1$s.

7. Summary and Comments

We have proposed a selection criterion for the choice of penalty strength in ridge regression. In the process, we have demonstrated an application of stochastic complexity to a common problem. This required that we devise a model encoding scheme in accordance with the typical stochastic complexity scheme, but which is sensitive to the application of regularisation. Our scheme suggests that we can consider the reduction in model complexity induced by regularisation as a constraint on model flexibility. We believe that this result is consistent with the intuitive message espoused by stochastic complexity - that we should be hesitant in choosing a model class that fits our data well if it provides good fits to unrelated datasets using the same predictors. Regularisation should ease our hesitancy by preventing our model from easily fitting these unrelated datasets.

The ridge penalty selection criterion that we derived produced promising results when used to fit linear readouts for echo-state networks. Moreover, the interpretation of stochastic complexity from which it was derived is generalisable to other regularisation techniques and different model classes.

Appendix A.1

Echo-state networks (ESNs) are a specific type of reservoir. They consist of a network and a set of input weights and behave like the recurrent part of a recurrent neural network. At each time step, the network is driven according to the value of the time series $u\left(t\right)$. The network state $x\left(t\right)$ updates according to the governing equation

$$x(t+1)=f(Ax\left(t\right)+W^{in}u\left(t\right)+b).$$

Unlike recurrent neural networks, the adjacency matrix A, input connections $W^{in}$ and node biases b all remain untrained. The activation function f, which in our case is the hyperbolic tangent $f(·)=tanh(·)$, introduces non-linearity. Predictions at each time step are given by the linear combination of node values ${\theta }^{T}x\left(t\right)$ with parameters $\theta$ determined from ridge regression ${\theta }_{Y\alpha }^{*}={(X^{T}X+\alpha \mathbb{I})}^{-1}{X}^{T}Y$. We train echo-state networks for one-step-ahead prediction. In this case, X is the matrix with a $t^{\mathrm{th}}$ row of $x\left(t\right)$ and Y is the vector with a $t^{th}$ entry of $u(t+1)$.

In the experiments performed in Section 6, every sample involves a new ESN architecture and a newly generated time series. We generate the network adjacency A as a weighted Watts-Strogatz network with a connectivity of 6 weightings drawn from the standard normal. The adjacency matrix is then normalised to have a spectral radius of $\rho$ which is sampled uniformly between $3/4$ and $3/2$. Input connection strengths are sampled from a Gaussian with zero mean and a standard deviation of $0.3$. We do not use node biases ($b=\overrightarrow{0}$).

For Figure 1 and Figure 2 the driving time series consists of 1200 points generated from a Runge-Kutta (RK4) integration of Lorenz-96 with a sampling rate of $0.1$s. Lorenz96 is a high-dimensional dynamical system with a Lyapunov time that we calculate to be about $0.68$s using the method in [49]. It has a governing set of equations detailed below in eq Appendix A.1, and we use the hyperparameters $(R,n)=(8,15)$.

$${\dot{x}}_i=(x_{i+1}-x_{i-2})x_{i-1}-x_i+R,$$

Note that in this equation, $x_0$ is the same as $x_n$.

For the time delay embedding plots (Figure 3 and Figure 4), we instead train the reservoir at one-step prediction using 1200 points from the Lorenz system. The governing equations of this system are provided below. We use standard hyperparameters $(\sigma ,\rho ,\beta )=(10,28,{\displaystyle \frac{8}{3}})$ for which we calculate the Lyapunov time to be around $1.1$ seconds.

$$\begin{array}{cc}\hfill \dot{x}& =\sigma (y-x)\hfill \\ \hfill \dot{y}& =x(\rho -z-y)\hfill \\ \hfill \dot{z}& =xy-\beta z\hfill \end{array}$$

The trajectories shown in Figure 3 and Figure 4 are generated autonomously. This means that rather than driving ESNs by an external signal, we instead drive them at each time step by their previous prediction. Effectively, their predictions are an Eulerian integration of the following system of equations

$$\begin{array}{cc}\hfill \widehat{u}\left(t\right)& ={\theta }^{T}x\left(t\right),\hfill \\ \hfill x(t+1)& =f(Ax\left(t\right)+W^{in}\widehat{u}\left(t\right)+b).\hfill \end{array}$$

This autonomous prediction is continued for a fixed duration of 200 time steps. For each predicted trajectory $\widehat{u}$, We graph $\widehat{u}\left(t\right)$ against $\widehat{u}(t+5)$ in Figure 3 and $\widehat{u}\left(t\right)$ against $\widehat{u}(t+1)$ in Figure 4. Both plots are comprised of 50 different trajectories, each of which is produced by a different ESN architecture trained on a different time series. The same 50 architectures and time series are used in each subplot.