Model-free prediction of spatiotemporal dynamical systems with recurrent neural networks: Role of network spectral radius

A common difficulty in applications of machine learning is the lack of any general principle for guiding the choices of key parameters of the underlying neural network. Focusing on a class of recurrent neural networks—reservoir computing systems, which have recently been exploited for model-free prediction of nonlinear dynamical systems—we uncover a surprising phenomenon: the emergence of an interval in the spectral radius of the neural network in which the prediction error is minimized. In a three-dimensional representation of the error versus the time and spectral radius, the interval corresponds to the bottom region of a “valley.” Such a valley arises for a variety of spatiotemporal dynamical systems described by nonlinear partial differential equations, regardless of the structure and the edge-weight distribution of the underlying reservoir network. We also find that, while the particular location and size of the valley depend on the details of the target system to be predicted, the interval tends to be larger for undirected than for directed networks. The valley phenomenon can be beneficial to the design of optimal reservoir computing, representing a small step forward in understanding these machine-learning systems.

Figure 1

Basic structure of reservoir computing. The left blue box represents an input subsystem that maps the $M$-dimensional input data to a vector of much higher dimension $N$, where $N\gg M$. The blue circle in the middle denotes the reservoir system, which can be, for example, a complex neural network of $N$ interconnected neurons, whose connection structure is characterized by the $N\times N$ weighted matrix ${\mathbf{W}}_{\mathrm{res}}$. The dynamical state of the $i\mathrm{th}$ neuron in the reservoir is $r_i$. The blue box on the right represents the output module, which converts the $N$-dimensional state vector of the reservoir network into an $L$-dimensional output vector, where $N\gg L$. The mapping from the input module to the reservoir is described by the $N\times M$ weighted matrix ${\mathbf{W}}_{\mathrm{IR}}$, and that from the reservoir to the output module by the $L\times N$ weighted matrix ${\mathbf{W}}_{\mathrm{RO}}$. During the training phase, the three blue boxes are activated. In this case, the whole computing device is effectively a nonlinear dynamical system with external input. In the prediction phase, the external input is cut off and the output data are directly fed back to the reservoir (green box), so the system is one without any external driving.

Figure 2

Using reservoir computing to predict Akhmediev breathers and the emergence of an optimal interval in the spectral radius. (a) For spectral radius $\rho =1.5$, a successful case of prediction of Akhmediev breathers, where the top panel shows the time evolution of the true solution, the middle panel displays the predicted solution from reservoir computing after training, and the bottom panel depicts the difference between the true and the predicted solutions. The color bar indicates the scale of the spatiotemporal wave magnitude $\left|\psi \right(x,t\left)\right|$. (b) For $0<\rho \le 2$ (ordinate), time evolution of the ensemble-averaged RMSE, denoted $\langle \text{RMSE}\rangle$, where, for each fixed value of $\rho$, 100 random realizations of the reservoir system are used to calculate the average and the color bar indicates the scale of the $\langle \text{RMSE}\rangle$ values. (c) Three-dimensional view of (b). The emergence of a valley interval in $\rho$ with a minimized prediction error can be seen unequivocally in (b) and (c). (d) Detailed time evolution of $\langle \text{RMSE}\rangle$ for four specific values of $\rho$ with standard deviation. (More quantitative details can be found in Appendixes pp1 and pp2.) Other parameters of the reservoir-computing system are $\alpha =1,N=4992,k=3,M=L=64,N_t=8010,S=10$, and $\mathrm{\Gamma }=1\times {10}^{-4}$.

Figure 3

Emergence of a valley in the prediction error versus the spectral radius of the reservoir network for predicting Kuznetsov-Ma solitons and soliton collision in the NLSE. (a)–(c) True (upper panel) and predicted (middle panel) wave patterns as well as the difference (lower panel) in wave function for Kuznetsov-Ma solitons: soliton collision for $\left(a_1=0.14,a_2=0.34\right)$ and $\left(a_1=0.42,a_2=0.18\right)$, respectively. (d)–(f) Corresponding time evolution of the ensemble-averaged prediction error $\langle \text{RMSE}\rangle$ for systematically varying $\rho$ values. For each fixed $\rho$ value, 100 random reservoir systems are used to calculate the average error and the color bar indicates the scale of the $\langle \text{RMSE}\rangle$ value. For predicting the Kuznetsov-Ma solitons, if the value of $\rho$ is chosen from the valley, then all realizations lead to near-zero errors (d). However, for the case of predicting soliton collision, only about half of the realizations yield near-zero errors even when the value of $\rho$ is chosen from the valley (outside the valley, large errors arise for nearly all realizations). Other parameter values are $N=4992,k=3,M=L=64,N_t=8010,S=10,\mathrm{\Gamma }=1\times {10}^{-4},\alpha =1.0$ for (a) and (d) and $\alpha =3.0$ for (b), (c), (e), and (f). Information about the standard deviation of $\langle \text{RMSE}\rangle$ can be found in Appendixes pp1 and pp2.

Figure 4

Emergence of a valley interval in the spectral radius of the reservoir network with minimum error for predicting spatiotemporal chaotic solutions of the KSE. (a) An example of successful prediction for $\rho =0.1$. Top panel: True spatiotemporal evolution of a typical chaotic solution of the KSE. Middle panel: Predicted spatiotemporal evolution. Lower panel: Difference between the true and the predicted solutions (true minus predicted). The color bar indicates the scale of $u(x,t)$. (b) Time evolution of ensemble-averaged RMSE (over 100 random realizations of the reservoir system) for systematically varying $\rho$ values in the interval (0.0,2.0], where the color bar indicates the scale of $\langle \text{RMSE}\rangle$. Note that the evolution time is represented as ${\mathrm{\Lambda }}_{m}t$, where ${\mathrm{\Lambda }}_m\approx 0.05$ is the maximum Lyapunov exponent of the chaotic solution and one unit of ${\mathrm{\Lambda }}_{m}t$ corresponds to one Lyapunov time. (c) Three-dimensional representation of $\langle \text{RMSE}\rangle$ in the $(\rho ,{\mathrm{\Lambda }}_{m}t)$ plane, revealing the emergence of a valley interval in $\rho$ in which $\langle \text{RMSE}\rangle$ is minimized. (d) Time evolution of $\langle \text{RMSE}\rangle$ for four representative values of $\rho$ with standard deviations (more details are given in Appendix pp1). Among the four cases, the best prediction result is achieved for $\rho =0.1$, where the chaotic solution can be predicted with near-zero error for about five Lyapunov time. Other parameters of the reservoir computing system are $\alpha =1,N=4992,k=3,M=L=64,N_t=70010,S=10$, and $\mathrm{\Gamma }=1\times {10}^{-4}$.

Figure 5

Emergence of a valley interval in the spectral radius of the reservoir network with minimum error for predicting spatiotemporal chaotic solutions of the CGLE. (a), (b) Example of successful prediction of the spatiotemporal chaotic solution of the one-dimensional CGLE, for which the maximum Lyapunov exponent is ${\mathrm{\Lambda }}_m\approx 0.22$. In (a), the true and predicted real and imaginary parts (above and below the black horizontal dashed lines, respectively) of the spatiotemporal evolution of the solution, together with their difference, are shown. In (b), the true and predicted magnitudes of the complex solution as well as their difference are displayed. The color bars in (a) indicate the scale of the real and imaginary parts of $u(x,t)$. For both (a) and (b), the value of the spectral radius of the reservoir network is $\rho =0.1$. (c) Ensemble-averaged RMSE calculated from the magnitude value of the complex solution versus $\rho$ and the Lyapunov time, where 100 random reservoir systems are used for each fixed $\rho$ value. (d) Three-dimensional view of $\langle \text{RMSE}\rangle$, where the color bar indicates its scale with the cutoff value of 3.0. The existence of a valley interval in $\rho$ that minimizes the prediction error is unequivocal. (e), (f) Time evolution of $\langle \text{RMSE}\rangle$ (with standard deviation) for five $\rho$ values. (Details of the statistical behavior of $\langle \text{RMSE}\rangle$ are presented in Appendixes pp1 and pp3.) Other parameters are $\alpha =1,N=9984,k=3,M=L=64,N_t=80010,S=10$, and $\mathrm{\Gamma }=2\times {10}^{-5}$.

Figure 6

Emergence of the optimal valley interval in the spectral radius with minimum prediction error for undirected random reservoir networks. (a) For a undirected random reservoir network, time evolution of the ensemble-averaged RMSE (over 100 random reservoirs) for systematically varying values of the network spectral radius for prediction of Akhmediev breathers in the NLSE. (b) Similar plot but for undirected small-world reservoir networks, where the value of the rewiring probability for generating the small-world topology is 0.3. In both cases, a similar valley region arises, indicating that the link topology of the reservoir network, directed or undirected, has little effect on the emergence of the valley. The location and size of the valley interval, however, do depend on the link topology, where undirected networks tend to lead to a larger interval. Parameter values are $\alpha =1,N=4992,M=L=64,N_t=8010,S=10,\mathrm{\Gamma }=1\times {10}^{-4},k=3$ for (a) and $k=4$ for (b).

Figure 7

Behavior of error during the training phase. Time-averaged error $\langle E\rangle$ versus spectral radius $\rho$ for a random reservoir network with (a) a directed topology and (b) an undirected topology with the same average degree as in (a). In each case, the network structure is fixed but the link weights are adjusted to result in systematic variations in the network spectral radius. In both cases, a region of small errors arises, indicating the existence of an optimal interval of spectral radius after training. The location and size of the region are similar to the valley interval in, e.g., Fig. 6. Parameter values are $\alpha =1,N=4992,k=3,M=L=64,N_t=70010,S=10$, and $\mathrm{\Gamma }=1\times {10}^{-4}$.

Figure 8

Standard deviation associated with the ensemble-averaged error in predicting different dynamical states of the NLSE. For each fixed value of the spectral radius $\rho$, 100 random realizations of the network are used to calculate the standard deviation of the ensemble-averaged prediction error. Shown are the 3D representation of the standard deviation versus the time and $\rho$ for four distinct dynamical states of the NLSE: (a) Akhmediev breathers, (b) Kuznetsov-Ma solitons, (c) a soliton-collision state for $a_1=0.14$ and $a_2=0.34$, and (d) another soliton-collision state, for $a_1=0.42$ and $a_2=0.18$.

Figure 9

Standard deviation associated with ensemble-averaged prediction error for (a) the KSE and (b) the 1D CGLE. Legends are the same as in Fig. 8.

Figure 10

Example of a long-term prediction of Akhmediev breathers in the NLSE. (a) True spatiotemporal evolution pattern, (b) reservoir-computing predicted pattern, and (c) difference in the instantaneous state between the true and the predicted patterns.

Figure 11

Time evolution of the RMSE for different statistical realizations in predicting Akhmediev breathers in the NLSE. Four cases are shown, each for a fixed $\rho$ value. For $\rho$ inside the valley interval, most or all realizations exhibit small errors, as in (b) and (c). For $\rho$ outside the interval, almost all realizations exhibit large errors, as in (a) and (d).

Figure 12

Time evolution of the RMSE for different statistical realizations in predicting soliton collisions in the NLSE. Parameters of the NLSE solution are $a_1=0.14,a_2=0.34$. Legends are the same as in Fig. 11.

Figure 13

Time evolution of the RMSE for different statistical realizations in predicting the spatiotemporal chaotic state of the 1D CGLE. Parameter values of the CGLE are the same as for Fig. 5. Legends are the same as in Fig. 11.