Constraining Chaos: Enforcing dynamical invariants in the training of
  recurrent neural networks
Drawing on ergodic theory, we introduce a novel training method for machine
learning based forecasting methods for chaotic dynamical systems. The training
enforces dynamical invariants--such as the Lyapunov exponent spectrum and
fractal dimension--in the systems of interest, enabling longer and more stable
forecasts when operating with limited data. The technique is demonstrated in
detail using the recurrent neural network architecture of reservoir computing.
Results are given for the Lorenz 1996 chaotic dynamical system and a spectral
quasi-geostrophic model, both typical test cases for numerical weather
prediction.
Jason A. Platt, Stephen G. Penny, Timothy A. Smith, Tse-Chun Chen, Henry D. I. Abarbanel
# Constraining Chaos: Enforcing dynamical invariants in the training of recurrent neural networks

###### Abstract

Drawing on ergodic theory, we introduce a novel training method for machine learning based forecasting methods for chaotic dynamical systems. The training enforces dynamical invariants--such as the Lyapunov exponent spectrum and fractal dimension--in the systems of interest, enabling longer and more stable forecasts when operating with limited data. The technique is demonstrated in detail using the recurrent neural network architecture of reservoir computing. Results are given for the Lorenz 1996 chaotic dynamical system and a spectral quasi-geostrophic model, both typical test cases for numerical weather prediction.

## 1 Introduction

Predicting the future trajectory of a dynamical system--a time series whose evolution is governed by a set of differential equations--is crucial in fields such as weather prediction, economics, chemistry, physics and many others [1, 2]. A prediction can be generated by deriving the governing equations of motion (EOM) for the system and integrating forward in time, perhaps with data being used to determine the value of particular constants or the initial conditions. Machine learning (ML), on the other hand, allows the construction of a forecast purely from observational data in lieu of a physical model. When the EOM are expensive to evaluate numerically, ML can be used to construct a surrogate model; such models can be integrated into data assimilation [3] algorithms--such as the Kalman filter [4, 5]--a typical use case when data are noisy and the model imperfect, such as in numerical weather prediction [6].

The inclusion of physical knowledge--EOM, conservation laws and dynamical invariants--into ML algorithms has been a topic of ongoing interest [7, 8, 9, 10, 11, 12, 13, 14, 15]. Enforcing these laws effectively reduces the searchable parameter space for a workable model, decreasing the training time and increasing the accuracy of the resulting models. An ML model trained without knowledge of invariants may fail to generalize and can produce solutions that violate fundamental constraints on the physical system [16]. Many of the examples cited above involve conservation of quantities based on the symmetry of the equations of motion, such as conservation of energy and momentum [13], or the inclusion of previously derived differential equations [11] as components of the ML training. "Physics informed" neural networks [14, 15, 17, 18, 19] add the known or partially known differential equations as a soft constraint in the loss function of the neural network, but conservation laws are not necessarily enforced and the equations need to be known.

Many physical dynamical systems of interest are dissipative--_e.g._, any dynamical system containing friction--meaning they exchange energy and mass with the surrounding environment [20]. High dimensional dissipative systems are very likely to exhibit chaos [21]--making them extremely sensitive to initial conditions. Enforcing conservation of quantities such as momentum, mass, or energy [11, 13, 22, 23] for dissipative systems in isolation may not be
sufficient for generalization due to the exchange of energy/momentum at the boundaries. Problems concerning chaotic dynamics, such as weather forecasting, exhibit fractal phase space trajectories that make it difficult to write down analytic constraints [11].

With the goal of enforcing dynamical invariants, we suggest an alternative cost function for dissipative systems based on ergodicity, rather than symmetry. This has broad implications for time series prediction of dynamical systems. After formulating the invariants, we give a recurrent neural network (RNN) [24] example applied to the Lorenz 1996 system [25] and quasi-geostrophic dynamics [26] where we add soft constraints into the training of the network in order to ensure that these quantities are conserved.

## 2 Deriving Dynamical Invariants

Ergodicity is a property of the evolution of a dynamical system. A system exhibiting ergodicity, called ergodic, is one in which the trajectories of that system will eventually visit the entire available phase space [27], with time spent in each part proportional to its volume. In general, the available phase space is a subset of the entire phase space volume. For instance a Hamiltonian system will only visit the hypersurface with constant energy [20]. Ergodicity implies that time averages over the system trajectories can be replaced with spatial averages

\[\lim_{t_{f}\rightarrow\infty}\frac{1}{t_{f}}\int_{0}^{t_{f}}g(F^{t}(\mathbf{u }_{0}))dt=\int_{\mathbf{u}\in B}g(\mathbf{u})\rho_{B}(\mathbf{u})du\qquad \forall\mathbf{u}_{0}\in B \tag{1}\]

for an arbitrary function \(g\), where \(F^{t}\) is the application of the flow of the dynamical system of \(t\) iterations, and \(\rho_{B}(\mathbf{u})\) defines an invariant density over the finite set \(B\)[27]. The invariant density gives an intuitive measure of how often a trajectory visits each part of \(B\). The invariant density defines the invariant measure [27]

\[\mu(B\subset R)=\int_{B}\rho_{B}(u)du. \tag{2}\]

\(B\) will often consist of exotic geometries such as quasi-periodic orbits and strange attractors [20]. The strange attractor is a hypersurface that contains the available subspace for a chaotic dynamical system--see Fig.(1). Deterministic chaotic systems are of importance to a vast array of applications such as in numerical weather prediction [3], chemical mixing [31, 32], optics [33], robotics [34] and many other fields.

Figure 1: left) 2D slice through the strange attractor of the spectral QG model reproduced from Fig.(3) in [26] using the implementation in [28]. An attractor is a hypersurface that draws in nearby trajectories of the system, such that the system will eventually be constrained to stay on the manifold. The motion of the dynamical system on a strange attractor is chaotic with extreme sensitivity to initial conditions. Strange attractors can be analyzed through the invariant measure that describes how often the system visits each part of the attractor [29, 30]. Related quantities such as the fractal dimension and the Lyapunov exponents [30] are globally invariant through any smooth change of coordinates and thus are natural invariant quantities for chaotic systems. right) The Lyapunov exponents of the spectral QG model. There are two positive exponents making the system chaotic.
Despite being deterministic, the precise long term prediction of chaotic systems is impossible due to the exponential growth of errors, as quantified by the system's Lyapunov spectrum. The Lyapunov spectrum, composed of a system's Lyapunov exponents (LEs), characterizes a dynamical system [29, 35] by giving a quantitative measure of how a volume of phase space stretches or shrinks over time.

For the prediction of chaotic systems, we suggest that although short term predictions will inevitably diverge, long term prediction of any system must preserve the invariants of motion characterized by the invariant measure \(\mu_{B}\) Eq.(2). Furthermore, enforcing such invariants could help to generalize the training of neural networks designed to emulate dissipative chaotic systems, in much the same way that conservation of energy and momentum has for conservative systems. While any function \(g(\mathbf{u})\) integrated with the invariant density is a constant--as seen in Eq.(1)--by the multiplicative ergodic theorem [36] the LEs and the fractal dimension are both invariant under smooth coordinate transformations and have algorithms that make them feasible to compute from observed data [29].

In the next sections we provide a concrete example using the fractal dimension and LEs as invariants that must be enforced when training a neural network. We use an RNN based on the reservoir computer (RC) architecture [37, 38, 39]. We impose a loss function that takes into account the preservation of the LEs and fractal dimension and detail the benefits of doing so. We stress that the concept is not limited to RC models and can in fact apply to any neural network architecture.

## 3 Recurrent Neural Networks and Reservoir Computing

An RNN is a network composed of nonlinear elements that are connected in such a way as to enable self excitation [40]. Therefore, given a state of the network \(\mathbf{r}(t-1)\), the next state

\[\mathbf{r}(t)=F_{r}(\mathbf{r}(t-1),\mathbf{u}(t-1),\theta) \tag{3}\]

is a function of the input \(\mathbf{u}\), the RNN equations \(F_{r}\) and the internal weights \(\theta\). The label over the input data \(t\in\mathbb{Z}\)--conveniently called time--gives the natural order and allows the analysis of the RNN as a dynamical map. \(\mathbf{r}(t)\) can then be decoded by a function \(W_{\mathrm{out}}(\mathbf{r}(t))=\dot{\mathbf{u}}\). \(W_{\mathrm{out}}\) is trained so that \(\dot{\mathbf{u}}\) is as close to the target output as possible [24]. In time series prediction tasks \(\dot{\mathbf{u}}\sim\mathbf{u}(t)\) so that the driven system Eq.(3) can become autonomous (with no external input)

\[\mathbf{r}(t)=F_{r}(\mathbf{r}(t-1),W_{\mathrm{out}}(\mathbf{r}(t-1)),\theta) \tag{4}\]

and predict the future of the dynamical system.

Reservoir computing (RC) [37, 38, 41, 42, 43, 44] is a simplified form of RNN for which only the large scale parameters of the network are varied with the detailed weights selected from probability distributions. For an RC with \(\tanh\) units at the nodes the RNN equations become [45]

\[\mathbf{r}(t)=\alpha\tanh(\mathbf{Ar}(t-1)+\mathbf{W}_{\mathrm{in}}\mathbf{u} (t-1)+\sigma_{b})+(1-\alpha)\mathbf{r}(t-1). \tag{5}\]

The elements of the \(N\times N\) adjacency matrix \(\mathbf{A}\) are fixed _i.e._, not trained--in contrast to other RNN architectures--with only its overall properties chosen such as the size \(N\), density \(\rho_{A}\) and spectral radius \(\rho_{SR}\). \(\mathbf{W}_{\mathrm{in}}\in\mathbb{R}^{N\times D}\) maps the input into the high dimensional reservoir space \(\mathbf{u}\in\mathbb{R}^{D}\rightarrow\mathbf{W}_{\mathrm{in}}\mathbf{u}\in \mathbb{R}^{N}\); the elements of \(\mathbf{W}_{\mathrm{in}}\) are chosen between \([-\sigma,\sigma]\). \(\alpha\) is the leak rate and is related to the time constant of the RC [39]. \(\sigma_{b}\) is an input bias governing the fixed point of the RC and the strength of the \(\tanh\) nonlinearity. See [39] for detailed explanations of the architecture and parameter choices.

Training the RC includes training the function \(W_{\mathrm{out}}\)--often taken to be a matrix \(\mathbf{W}_{\mathrm{out}}\) and trained through linear regression--as well as finding the correct parameters \(N,\ \rho_{A},\ \rho_{SR},\ \sigma,\ \sigma_{b}\). In [6, 45, 46] it is shown how to train the RC network through a two step training procedure that takes into account both the one step prediction accuracy as well as the long term forecast skill.

RC has been shown to be extremely successful in time series prediction tasks. Its simple form allows the easy computation of the Jacobian and other quantities that can help in a dynamical systems analysis of the RNN. In Platt et al.[48] the authors showed that the RC, when well trained, can reproduce invariants of the motion such as the LEs and fractal dimension and that the reproduction of these quantities maximized the prediction time and ensured the stability of the predictions. The training procedure in those previous works does not enforce these invariants explicitly. The hope in those previous works is that the RC is both capable of reproducing these quantities and that the loss function of short and long term forecasts guides the RC towards these values by proxy. Here we reformulate the training to take into account the invariant quantities.
calculated as the time when the root mean square error \(\mathrm{RMSE(t)}=\sqrt{\frac{1}{D}\sum_{i=1}^{D}\left[\frac{u_{i}^{\mathrm{f}}(t)- u_{i}(t)}{\sigma_{i}}\right]^{2}}\) exceeds a certain value \(\epsilon\), in this case \(\epsilon=0.3\) approximately in line with [6, 45, 47]. \(D\) is the system dimension, \(\sigma\) is the long term standard deviation of the time series and \(u^{f}\) is the RC forecast. The largest LE \(1/\lambda_{1}\sim 12\) days gives the natural time scale for error growth of the system and thus can be used as a measurement for the predictive skill. The RC returns fantastic prediction times for this low resolution model. The size of the RC is N=1500 and 100,000 training steps were provided with a \(\Delta t=80min\) giving about 30 years of data, we consider this the "data rich" case.

## 4 Enforcing Invariants

The training of an RC is determined by the training data and the selection of the parameters governing the global properties of the RC: \(N,\ \rho_{A},\ \rho_{SR},\ \sigma,\ \sigma_{b}\) and a regularization coefficient \(\beta\). Once these quantities are chosen and the weights instantiated, then \(\mathbf{W}_{\mathrm{out}}\) is given

\[\mathbf{W}_{\mathrm{out}}=\mathbf{u}\mathbf{r}^{T}(\mathbf{r}\mathbf{r}^{T}+ \beta\mathbb{I})^{-1}\]

where \(\mathbf{u}\) is the \(D\times T\) matrix of input data, \(T\) is the number of time steps and \(\mathbf{r}\) is the \(N\times T\) matrix of reservoir states [45]. For an RC we can add into the selection of these parameters knowledge of the global invariants of the system \(\mathbf{u}\). Therefore we construct a loss function

\[\mathrm{Loss}=\epsilon_{1}\|\mathrm{C_{u}}-\mathrm{C_{RC}}\|^{2}+\epsilon_{2} \sum_{k=1}^{M}\sum_{t=t_{i}}^{t_{f}}\left\|\mathbf{u}_{k}^{\mathrm{f}}(t)- \mathbf{u}_{k}(t)\right\|^{2}\exp\biggl{\{}-\frac{t-t_{i}}{t_{f}-t_{i}}\biggr{\}} ;t\in\mathbb{Z} \tag{6}\]

that can be minimized over the RC parameters and with \(\epsilon_{x}\) hyperparameters. The selection of the parameters leads to the matrix \(\mathbf{W}_{\mathrm{out}}\) based on training data \(\mathbf{u}_{\mathrm{train}}\). Platt et al. [45] generated a number of long term forecasts \(\mathbf{u}^{f}(t)\) and compared them to the data \(\mathbf{u}\); with enough data this procedure often leads to a model that reproduces the correct dynamical invariants. Without the explicit enforcement of these invariants, however, the model can fail to capture the dynamics--particularly for high dimensional systems and in cases where the number of trajectories \(M\) is constrained. Here we add the dynamical invariants (\(C_{x}\)) as a constraint in order to directly train for generalizability, similar to [49]. This scheme is illustrated in Fig.(3). The global optimization routine used to minimize the cost function was the covariance matrix adaption evolution strategy (CMA-ES) [50].

We show the Lyapunov exponents and the fractal dimension as examples of dynamical invariants in order to demonstrate the technique. With the equations of motion, such as Eq.(4) for the RC, it is quite simple to calculate these quantities using well known and efficient algorithms [30]. When training directly from data--without knowledge of the underlying system--we may not know the equations of motion so these quantities must be estimated. The largest LE can often be approximated from time series data [29, 51, 52] and the fractal dimension can be calculated using various techniques [29, 53]. A calculation of the full LE spectrum is more difficult. Use of other dynamical invariants derived from the invariant measure Eq.(2) are also possible, for instance the energy density spectrum of a fluid dynamical system as a function of wavenumber.

Figure 2: Example prediction and the probability distribution function of the valid prediction time (VPT) for 200 initial conditions over the 20 dimensional QG system described in section 5.2. The valid prediction time [6, 47] is \(\epsilon\), in this case \(\epsilon=0.3\) approximately in line with [6, 45, 47]. \(D\) is the system dimension, \(\sigma\) is the long term standard deviation of the time series and \(u^{f}\) is the RC forecast. The largest LE \(1/\lambda_{1}\sim 12\) days gives the natural time scale for error growth of the system and thus can be used as a measurement for the predictive skill. The RC returns fantastic prediction times for this low resolution model. The size of the RC is N=1500 and 100,000 training steps were provided with a \(\Delta t=80min\) giving about 30 years of data, we consider this the “data rich” case.
## 5 Results

### Lorenz 1996

Our first test case for the RC is the Lorenz 1996 system (L96), a standard testbed for data assimilation applications in numerical weather prediction. L96 describes the evolution of a scalar quantity over a number of sites scattered uniformly over a periodic spatial lattice of constant latitude with quantities approximating advection and diffusion

\[\frac{\mathrm{d}u_{k}}{\mathrm{d}t}=-u_{k-1}(u_{k-2}-u_{k+1})-u_{k}+F. \tag{7}\]

In this case we take the number of sites to be \(D=10\) and forcing \(F=8\) with the purpose of making the system hyperchaotic, with three positive Lyapunov exponents Fig.(4).

Figure 4: Lyapunov exponents of the 10D Lorenz 1996 system with F=8. There are three positive LEs, a single zero exponent, and the rest negative.

Figure 3: Parameter optimization of a reservoir computer showing the introduction of dynamical invariants into the routine. The observed data is split into training, validation and testing sets with the invariants calculated from the data [29]. These quantities can then be incorporated into the loss function to improve the overall training of the reservoir computer. A general discussion of the training strategy is found in [45].
The results for \(C_{RC}=LEs\) Eq.(6) are shown in Fig.(5). When no global information is given to the RC then it can fail to generalize when presented with unseen input. Simply providing the largest LE to the RC during training enables the neural networks to

1. generalize to unseen data so that there are good predictions over the entire range of possible initial conditions
2. reconstruct the attractor as in [54] and [48] with the prediction giving the correct ergodic properties of the data even after the prediction necessarily diverges from the ground truth.

In this case providing the largest exponent was enough to improve the predictions, with no further gains coming from providing the smaller exponents. This could perhaps be due to the parameter space being constrained enough for those exponents to be matched by the RNN even though they are not directly given.

The second invariant given is the fractal dimension of the data calculated through the Kaplan-Yorke formulation

\[\text{Dimension}=\alpha+\frac{\sum_{i=1}^{\alpha}\lambda_{i}}{|\lambda_{\alpha +1}|}, \tag{8}\]

with \(D\) the system dimension, \(\lambda\) the ordered Lyapunov exponents and \(\alpha\) the smallest index where the sum of the LE doesn't cross zero [29]. There are alternate definitions and methods for calculating the fractal dimension [53].

Providing the fractal dimension as an invariant has a similar effect to providing the LEs by raising the mean valid prediction time from \(\sim 2\to 4\). The fractal dimension may not, however, be unique to a particular set of data while the LE spectrum has a much greater chance of constraining the shape of the resulting strange attractor. Therefore we see that there is not as much improvement in the forecast when providing the fractal dimension compared to the LE spectrum.

### Synoptic Scale Atmospheric Model

For more complex higher-dimensional dynamical systems, the Lyapunov spectrum or Kaplan-Yorke dimension are quite difficult, if not impossible, to calculate. However, our previous results showed that capturing the leading Lyapunov exponent (LLE) enhanced prediction skill greatly, and even with complex models this quantity can be estimated more readily either from data [51, 52] or from a model [55, 56]. We therefore explore the value of representing the LLE

Figure 5: The RC is initialized 10 times, providing k LEs and M=7 long term forecasts Eq(6); we report on the average of the distribution of predictions for the 10 dimensional Lorenz 1996 system [25] Eq(7). There are 10 total LEs—3 positive, shown in Fig.(5). When 0 LEs are provided the RC has no global information and the prediction time is poor. Giving the largest LE is sufficient for generalizable predictions and this quantity is quite easily calculated from numerical data [29]. The last column shows the result for providing the fractal dimension as another example of an invariant quantity. The size of the RC is \(N=400\) and the number of test initial conditions was 1000. The forecast time is given as an average over the 1000 trajectories and then scaled by the largest Lyapunov exponent to give the prediction in terms of the number of Lyapunov timescales.
in the more complex case of quasi-geostrophic (QG) dynamics [57], which assume that large-scale atmospheric disturbances are governed by the conservation of potential temperature and absolute potential vorticity, while the horizontal velocity is quasi-geostrophic and the pressure quasi-hydrostatic. Numerical models based on the QG approximation were a precursor to larger scale primitive equation models used for global numerical weather prediction [3] and is frequently used in data assimilation studies targeting the atmosphere and ocean [58, 59]. Here, we consider the two-layer baroclinic model of Charney and Strauss (1980) [60] used to study the planetary-scale motions of a thermally driven atmosphere in the presence of topography. We further incorporate the adaption of Reinhold and Pierrehumbert (1982) [26] to include an additional wave in the zonal direction making it highly baroclinically unstable. We use the implementation of [28], which provides a truncated 2-layer QG atmospheric model on a mid-latitude \(\beta-\)plane frictionally coupled to a mountain and a valley with a dimension of 20 in the spectral space of the model.

For the atmospheric streamfunctions \(\psi_{a}^{1}/\psi_{a}^{3}\) at heights 250/750 hPa and the vertical velocity \(\omega=\frac{\mathrm{d}p}{\mathrm{d}t}\), the equations of motion are derived to be

\[\frac{\partial}{\partial t}\big{(}\overline{\nabla^{2}\psi_{a}^{1}}\big{)}+ \overbrace{J(\psi_{a}^{1},\nabla^{2}\psi_{a}^{1})}^{\text{horizontal advection}}+ \overbrace{\beta\frac{\partial\psi_{a}^{1}}{\partial x}}^{\text{Princision force}}=-\overbrace{k_{d}^{\prime}\nabla^{2}(\psi_{a}^{1}-\psi_{a}^{3})}^{ \text{friction}}+\overbrace{\frac{f_{0}}{\Delta p}\omega}^{\text{vertical stretching}} \tag{9}\]

\[\frac{\partial}{\partial t}\big{(}\nabla^{2}\psi_{a}^{3}\big{)}+J(\psi_{a}^{3},\nabla^{2}\psi_{a}^{3})+J(\psi_{a}^{3},f_{0}h/H_{a})+\beta\frac{\partial\psi_ {a}^{3}}{\partial x}=k_{d}^{\prime}\nabla^{2}(\psi_{a}^{1}-\psi_{a}^{3})-k_{d }\nabla^{2}\psi_{a}^{3}+\frac{f_{0}}{\Delta p}\omega \tag{10}\]

with \(\nabla=\frac{\partial}{\partial x}\hat{x}+\frac{\partial}{\partial y}\hat{y}\), \(k_{d}^{\prime}\) the friction between the layers, \(k_{d}\) the friction between the atmosphere and the ground, \(h/H_{a}\) the ratio of ground height to the characteristic depth of the atmospheric layer, \(\Delta p=500\) hPa the pressure differential between the layers and \(J(g_{1},g_{2})=\frac{\partial g_{1}}{\partial x}\frac{\partial g_{2}}{ \partial y}-\frac{\partial g_{3}}{\partial y}\frac{\partial g_{2}}{\partial x}\) the Jacobian. More details are given in [26, 28].

After integrating the model forward in time we ask if an RC model is capable of predicting the dynamics given a significant amount of training data. An example forecast and distribution of predictions is shown in Fig.(2) where the RC successfully predicts the synoptic scale atmospheric dynamics for a number of months. Such significant predictive power on a low resolution QG model is an interesting result in and of itself, showcasing the ability of RNNs to resolve more realistic atmospheric dynamics.

When we reduce the amount of data and set \(M=1\)--the limited data case in Eq. (6) with a single forecast as part of the training loss--the reservoir loses its predictive power. Adding in the information contained in the LLE, which is \(\sim 0.01\), enables the model to recover a large amount of predictive capability. The calculated LLE of the RC with no provided exponents is 0.14, compared to \(\sim 0.01\) when it is provided. The mismatch between the LEs of the two systems is a clear indication that synchronization between the two is not achieved [48]. The impact of the added information provided via the LLE is clear in Figure 6, where the average VPT has extended from only a few days to multiple weeks.

## 6 Discussion and Conclusion

Chaotic dynamical systems are difficult to predict due to their sensitivity to initial conditions [61]. Better understanding and accounting for dynamical uncertainties has, however, allowed fields like numerical weather prediction to provide useful and continually improving forecasts [62]. Previous works (e.g. [13, 49]) proposed that including conserved quantities such as energy/momentum may help to improve the application of neural networks to physical systems. However, the introduction of the proposed conserved quantities is not generally applicable to dissipative chaotic dynamical systems. Thus, we instead considered dynamical invariants based on the invariant measure.

We provided a concrete example using quantities derived from the invariant measure, such as the Lyapunov exponents and the fractal dimension, to train a particular RNN architecture called reservoir computing. Previous RC training algorithms used long-term forecasts initialized from many different initial conditions in order to improve generalizability [6, 46], essentially imposing these invariant measures by proxy. Here, we imposed the invariant measures as constraints directly in the training algorithm, allowing the RC to generalize with fewer data. Fortunately, we have found that much of the value of this additional constraint is achieved through the use of the leading Lyapunov exponent. While the entire Lyapunov spectrum can be quite difficult to calculate, particularly for large systems, the leading Lyapunov exponent can be estimated by using numerical techniques such as the breeding method [63, 64] or other methods described in [21, 55, 56]. This provides an opportunity for extension of this technique to higher-dimensional systems.

Recent works from [65, 66, 67] have shown promise in producing data-driven surrogate weather models that are competitive by some metrics with conventional operational forecast models. A key property that has not yet been demonstrated
with such surrogate models is their ability to reproduce dynamical quantities such as the LEs, which indicate an average measure of the response to small errors in the initial conditions. For weather models in particular, the enforcement of LEs is crucial for the correct operation of data assimilation algorithms [3]. Platt et al. [48] demonstrated the importance of reconstructing the LE spectrum for producing a skillful deterministic forecast model. Similarly, Penny et. al. [6] indicated the ability of the RC to reproduce accurate finite-time LEs as a requirement for RC-based ensemble forecasts to produce good estimates of the forecast error covariance, which is the primary tool used in conventional data assimilation methods to project observational data to unobserved components of the system. This information can then be used to make the RC robust to sparse and noisy measurements. While this is the more realistic scenario used in online weather prediction systems, it is a fact that is rarely taken into account in neural network applications. The introduction of explicit constraints in the training cost function both improves prediction and trains the RNN to reconstruct the correctly shaped attractor [45, 54].

## 7 Acknowledgements

J.A. Platt, S.G. Penny, and H.D.I. Abarbanel acknowledge support from the Office of Naval Research (ONR) grants N00014-19-1-2522 and N00014-20-1-2580. S.G. Penny and T.A. Smith acknowledge support from NOAA grant NA20OAR4600277. T.-C. Chen acknowledges support from the NOAA Cooperative Agreement with the Cooperative Institute for Research in Environmental Sciences at the University of Colorado Boulder, NA17OAR4320101.

## 8 Source Code

The basic RC implementation used in this study is available

[https://github.com/japlatt/BasicReservoirComputing](https://github.com/japlatt/BasicReservoirComputing)
