Bayesian optimization of hyper-parameters in reservoir computing

Echo state networks were introduced independently by Herbert Jaeger [22] and Wolfgang Maass (under the name Liquid State Machines [23]). The idea is to use a recurrent neural network which is fixed, i.e., its connection weights and any biases are not trained. One can drive this network using an input, which will change the state of the network (i.e., the value at each of the nodes). The output is then constructed using a linear combination of all, or a subset of, the node values. The weights of this linear combination are most commonly obtained using linear regression. The values of the nodes are updated according to the rule:

$$\bm{x}\left(n+1\right)=\sigma\left(\bm{Wx}(n)+\bm{W}^{\textit{in}}\bm{u}(n+1)+\bm{b}\right),$$

(1)

where $\bm{x}(n)$ is the $N$-dimensional state vector of the network at time $n$, $\bm{u}(n)$ is the $K$-dimensional input vector at time $n$, $\bm{W}$ is the $N\times N$ reservoir weight matrix, $\bm{W}^{\textit{in}}$ is the $N\times K$ input weight matrix, $\bm{b}$ is a constant bias term and $\sigma(.)$ is a sigmoid function (in our case we use the $\tanh$ function).

The reservoir weight matrix is rescaled as $\bm{W^{\prime}}=\rho\bm{W}/\left|\lambda_{\textit{max}}\right|$, where $\left|\lambda_{\textit{max}}\right|$ is the spectral radius of the network (i.e., the largest eigenvalue of $\bm{W}$) and $\rho$ is a scaling parameter (effectively the spectral radius of the rescaled network). This rescaling operation is a popular choice in the reservoir computing literature. It, however, does not guarantee the echo state property [24]. See, e.g., [25, 26] for rigorous discussions on sufficient conditions for the echo state property.

The output is given by:

$$\bm{y}(n)=\bm{U}^{\textit{out}}\bm{x}(n),$$

(2)

where the output weights matrix $\bm{U}^{\textit{out}}$ is obtained using, e.g., linear regression. We performed linear regression using ridge regression (a.k.a. Tikhonov regularization), which is a modified version of the linear regression equations [7]:

$$\bm{U}^{\textit{out}}=\left(\bm{X}^{T}\bm{X}+\lambda^{2}\bm{I}\right)^{-1}\bm{X}^{T}\bm{y},$$

(3)

where $\bm{I}$ is the identity matrix, $\lambda$ is the regularization parameter, and $\bm{X}$ is the matrix of all reservoir states. Performance is measured with the normalized mean squared error (NMSE):

$$\text{NMSE}=\frac{\langle\lVert\hat{y}(t)-y(t)\rVert^{2}\rangle}{\langle\lVert y(t)-\langle y(t)\rangle\rVert^{2}\rangle},$$

(4)

where $\hat{y}(t)$ is the predicted value and $y(t)$ is the target value. The denominator is simply the variance of the test set under consideration. The root mean squared version (NRMSE) is:

$$\text{NRMSE}=\sqrt{\frac{\langle\lVert\hat{y}(t)-y(t)\rVert^{2}\rangle}{\langle\lVert y(t)-\langle y(t)\rangle\rVert^{2}\rangle}}.$$

(5)

The concept of a nonlinear delay node as a reservoir was introduced in [27]. It can be implemented in hardware with optical and electronic components [28, 29, 30, 12, 31, 8, 9, 10, 11, 12]. We use two different types of NDNs.

The Santa Fe data set consists of the output of a chaotic laser system [40]. We use the A.cont data set, which can be obtained from [41]. The goal is to predict the next step of the time series.

The Nonlinear Auto-Regressive Moving Average (NARMA) of order 10 was originally introduced for use as a timeseries prediction benchmark in [42]. It is generated by

$$y(t+1)=0.3y(t)+0.05y(t)\sum_{i=0}^{9}y(t-i)+1.5s(t-9)s(t)+0.1,$$

(10)

with $s(t)$ a random term with uniform distribution in $[0,0.5]$. The goal is to do a one-step ahead prediction of $y(t)$.

This task was first introduced in [43] and has since been used as a benchmark several times [44, 29, 45]. The data set is created as follows: An i.i.d. sequence $d(t)$ is generated by randomly choosing values from $\{-3,-1,1,3\}$. This signal is then passed through a linear channel described by:

	$\displaystyle q(n)=$	$\displaystyle 0.08d\left(n+2\right)-0.12d\left(n+1\right)+d\left(n\right)+0.18d\left(n-1\right)-0.1d\left(n-2\right)+0.091d\left(n-3\right)$
		$\displaystyle-0.05d\left(n-4\right)+0.04d\left(n-5\right)+0.03d\left(n-6\right)+0.01d\left(n-7\right),$

which in turn is passed through a nonlinear channel:

$$u(n)=q(n)+0.036q(n)^{2}-0.011q(n)^{3}+\nu(n),$$

where $\nu(n)$ is an i.d.d. Gaussian noise with zero mean and a standard deviation adjusted to yield signal-to-noise ratios (SNR) ranging from $12$ to $32$ dB. The objective is to reconstruct $d(n-2)$ given $u(n)$, for several values of the SNR (so $u(n)$ is the input and the target is $y(n)=d(n-2)$). Following [44, 46], we shifted $u(n)$ by +30. The quality measure of the algorithm is given by the Symbol Error Rate (SER), which is the fraction of incorrect symbols obtained.

If $p$ is uneven, the MG dynamics Eq. (6) is unstable for negative values. We therefore rescale $u(k)$ so that it is always positive. We furthermore add an offset parameter $\delta\geq 1$ to the mask $M\in\{-1+\delta,1+\delta\}^{N}$.

For comparison purposes, the number of nodes $N$ is fixed. The (time) distance between the nodes is denoted by $\theta$. A larger $\theta$ gives a longer relaxation time between the nodes. By definition one has that $\tau=N\theta$.

We employ Euler integration to solve the MG equation. The number of integration steps between adjacent nodes determines the precision of the obtained solution. However, the method used to solve the differential equation does not necessarily influence the quality of the reservoir. Put differently, one does not need an exact solution of the differential equation to have a high-quality reservoir. The essential property of the reservoir is a “good” nonlinear projection to a high-dimensional space, which is not necessarily the exact solution of the differential equation. We perform Euler integration with either several integration steps between adjacent nodes (number of steps $2|\lfloor\theta/0.1-1\rfloor|+2$), or with one integration step between adjacent nodes. We refer to these methods as, respectively, multi-step integration (MSI) and one-step integration (OSI). Both choices give good results, with the latter method obviously being faster.

For all NDN implementations, ridge regression is performed with sci-kit learn [47].

To implement the echo state networks we made use of the software library Oger [48], which in turn builds on the Modular toolkit for Data Processing (MDP) [49].

For the implementation of the GP we used the Spearmint library, available at https://github.com/HIPS/Spearmint. A fork of this framework, as well as a step-by-step tutorial detailing both installation and a practical example (NDN MG on NARMA 10), is available at https://bitbucket.org/uhasseltmachinelearning/spearmint/. This code replaces the MongoDB backend with an SQLite one. This removes the need for a background server process, which makes it easier to deploy in a cluster environment or as part of a Jupyter Notebook.

For the MG NDN, we compare our results with those of [50]. Details of the implementation can be found in A. As in [50], we take $N=200$. Before the other parameters are optimized, ten different masks are generated, and the best performing one is selected. There are six parameters to optimize: $p$, $\gamma$, $\eta$, $\theta$, $\delta$, and the regularization parameter $\lambda$. Spearmint runs are performed for both OSI and MSI. One needs to specify the boundaries of the hyper-parameter search space. We therefore start by running Spearmint for different boundary settings, to find reasonable areas for the hyper-parameter search. Correct order of magnitudes for the boundaries can also be estimated from previously published results and physical arguments. For different boundary values, the optimum for $p$ always converged to a value close to 1. Because a smaller number of variables implies a smaller search space, especially in high dimensions, we set $p=1$ in our final run.

The final results are shown in Table 1. The NMSE for the MSI is an order of magnitude lower than [50], while the OSI result is four times lower. A plot of the NMSE on the validation set as a function of the number of Spearmint iterations is shown in Figure 4. If we start with good hyper-parameter boundaries, the first runs show a strong decrease in error, after which further improvements come in steps. This behaviour qualitatively resembles the exponential convergence speed found in analytical work, cf. Section 3. It should be noted that different runs with small changes in the boundaries of the search space lead to parameter sets that are quite different, but have the same performance.

Some error surfaces are shown in B. One sees that the error surface is smooth around the optimum parameter values for both OSI and MSI. The error surfaces of the OSI and MSI differ significantly. The error surface is non-stationary. For example, the variation in the $\gamma$ direction is low near the optimum, Figure 10, while it is high away from the optimum, Figure 12. A difference in spatial variation between different parameters can be seen in Figure 11, where $\gamma$ exhibits a slow variation and $p$ a faster variation. This difference in variation is taken care of by Spearmint, as discussed in Section 3. Note that $\delta$ has little to no influence on the performance. Given the smoothness of the error surface, we expect good convergence behaviour from the GP, as is witnessed in the results.

We now compare our optimal parameters to those from reference [50]. One should keep in mind that the integration schemes are slightly different from the implementation in [50], even though Euler integration is also used there. Our programs are written in Python, while Matlab is used in [50]. For reference, our MSI method with Appeltant’s parameters gives NMSE = 0.028, which is similar to his result NMSE = 0.02 (note that we don’t know all the parameters of [50], some are educated guesses). All parameters differ at most one order of magnitude. It appears that higher values of $\gamma$ and $\eta$ are useful, although this seems to be specific to the Santa Fe task. A more important conclusion is that the value of $\theta=0.2$ is suboptimal. Larger $\theta$ values allow for better performance in most considered benchmark tasks. This is important, because the choice $\theta=0.2$ has become popular in the literature [28, 12, 51, 52, 53, 30, 54, 55, 56, 57]. A larger $\theta$ causes a longer calculation time for the MSI. The OSI calculation time is independent of $\theta$, and has the same order of magnitude improvement in NMSE.

In addition to the NDN, we studied an ESN to compare the performance of the Bayesian optimization algorithm to the results from [58], who use binary PSO to optimize the ESN architecture. They use $6400$ samples for training, $1700$ for validation and another $1700$ for testing. The optimized parameters were the input scaling (which is a scaling parameter of $\bm{W^{\textit{in}}}$ in Eq. (1)), spectral radius $\alpha$ and the ridge regression parameter $\lambda$. The number of nodes in the network was set to $N=200$. The validation error was averaged over 10 independent realizations of the ESN. Using Bayesian optimization on these parameters we obtained a NMSE of $0.00694\pm 0.00087$, which should be compared to $0.0284$ in [58] for their standard ESN, and $0.0143$ after applying their optimization algorithm. This is a significant improvement compared to their results. During testing an irregularity in the data set was found, on which we elaborate in A. In the partitioning we study, the irregularity is part of the validation set. It should be noted that the irregularity may have driven the optimization to a hyper-parameter set that minimizes the error caused by the irregularity, which isn’t necessarily the optimal set for prediction on the test set, as no such deviations occur there. In order to check the performance of the ESN on the Santa Fe laser task without the complications of the irregularity, we also performed some tests using the Appeltant division of the data set, as was used to test the performance of the delay node algorithm. The resulting NMSE was $0.00693\pm 0.00097$. Even though the error is approximately the same as for the other data set division, the result is better since only 1000 points instead of 6400 points are used for training.

We used the same methods for the NARMA 10 time series. We train and validate on two fixed time series. The test error shows quite some variance, so we average the test error over 15 randomly generated time series of 2000 points each. We found an NRMSE $=0.08$ for the OSI scheme, which compares favorably to NRMSE $=0.12$ of Appeltant [50]. The histogram of test errors is shown in Figure 5. Most values lie around the average, but there are some outliers with higher NRMSE. We don’t know if such an average was taken in [50]. The MSI scheme gives NRMSE $=0.146$, which is worse than the result from [50]. Note that $\theta\approx 0.2$ for the optimal MSI scheme (compared to $\theta=0.908$ for OSI). The MSI scheme for NARMA 10 was used in [50] to arrive at the choice of $\theta=0.2$.

Because NARMA 10 is deterministic given the ten previous input steps, the linear regression weights of the output are extremely large (a factor 100 larger than Santa Fe, which itself has large weights because it is deterministic). Large regression weights occur when overfitting, but in this case they give good results. The regularization parameter should be set to zero $\lambda=0$, i.e., no regularization is needed. This “overfitting” issue also leads to a strong sensitivity on noise in the implementation, making hardware implementations inherently difficult [50]. It is possible that rounding errors in the MSI scheme are the reason for the worse performance, which makes the calculation of reservoir states less stable.

From the plots of the error surfaces in C one sees that there is a very low sensitivity for $\delta$. The error surface is less smooth than the Santa Fe error surface, but still rather smooth. While the MSI scheme is indeed noisier compared to the OSI scheme, the amount of noise seems small enough for Spearmint to handle (this amount of noise was not problematic for the demo runs, see Section 7). The fact that the MSI error surface is noisier corroborates our idea that the worse performance of the MSI scheme is because of rounding error in the integration scheme. Away from the optimal parameter values, the error surface can be very noisy, cf. Figure 21. We finally note that the error surfaces are non-stationary here as well.

We now assess the performance on the nonlinear channel equalization task (Section 4.3). In order to compare with the results found in [29] we use the delay node system with the sine nonlinearity, as described in Section 2.2.2. Following the setup in [29] we generate $10$ different realizations of the data set, for each of the $6$ SNRs (ranging from $12$ to $32$ dB). The delay node was trained using an initial washout of $200$ samples, followed by $3000$ points of training, and was validated on $10^{6}$ validation points. This was done for each of the $10$ realizations of the data set, and the validation error was the average of the error on the individual realizations. To produce the output, ridge regression was used, after which the resulting values were binned to the values $\{-3,-1,1,3\}$. The reservoir consists of $N=50$ nodes.

The parameters that were optimized using Bayesian optimization are $\alpha$, $\beta$, $\gamma$ and the regularization parameter $\lambda$. These were optimized separately for each of the SNRs. The results are shown in Table 3. The results of the SER are plotted in Figure 6. For low noise levels, Bayesian optimization finds parameters which result in a performance that is 2 orders of magnitude better than the literature. For high noise the performance is on par with [29]. Presumably, the error surface at high levels of noise becomes less well behaved, and the assumption of smoothness breaks down, reducing the algorithm to random sampling. This is apparent at SNR 16, where the $\alpha$ parameter is larger by 2 orders of magnitude when compared to the other values of the SNR.

We also used an ESN on this data set, and compare our results with [44]. Mimicking their setup we used an ESN with 47 nodes, used $200$ washout points, $5000$ training points and validated on $10^{6}$ points. The validation error was an average over the error on $10$ independent realizations of the network. Once the hyper-parameters were determined on the validation set, a test was performed on a test set of $10^{6}$ points, averaged over $100$ separate realizations of the network. This was done for $10$ instantiations of the data set, and the final result is the average of the performance on each of these data sets. The result is shown in Figure 6. We observe a performance which is an order of magnitude better than [44] for high values of the SNR. Similar to the results with the NDN architecture with sine nonlinearity, our method converges to the results obtained in other works for the lower values of the SNR .