Lyapunov analysis of data-driven models of high dimensional dynamics using reservoir computing: Lorenz-96 system and fluid flow

We computed the Lyapunov spectrum and finite-time Lyapunov exponents of a data-driven model constructed using reservoir computing. This analysis was performed for two dynamics that exhibit a highly dimensionally unstable structure. We focused on the reconstruction of heterochaotic dynamics, which are characterized by the coexistence of different numbers of unstable dimensions. This was achieved by computing fluctuations in the number of positive finite-time Lyapunov exponents.


Introduction
Interest in modeling chaotic dynamics using machine learning techniques is increasing.Lu et al [1] reported that a data-driven model using reservoir computing has an attractor similar to that of the reference system.Nakai and Saiki [2,3] confirmed that a single data-driven model obtained from fixed training data can infer the time series of chaotic fluid flows from various initial conditions.They suggested that such a data-driven model can reconstruct the attractor of a chaotic fluid flow.The effectiveness of machine learning methods has been demonstrated in the fields of atmospheric and climate science.Considerable progress has been made in forecasting phenomena such as the El Niño Southern Oscillation [4,5], Asian summer monsoons [6], and hurricanes [7] at incomparably low computational costs.
Reservoir computing, a brain-inspired machine learning technique using a data-driven dynamical system, is effective in predicting time series and frequency spectra of chaotic behaviors, including fluid flow and global atmospheric dynamics [2,[8][9][10][11][12][13][14][15][16][17][18].Pathak et al [10] examined the Lorenz-63 system and Kuramoto-Sivashinsky system and reported that the data-driven model obtained from reservoir computing can generate an arbitrarily long time series that mimics the dynamics of the original systems.
In our previous study [19], we demonstrated that a data-driven model using reservoir computing obtains richer information than that obtained from the training data, particularly from a dynamical system perspective.This indicates that the dynamical properties of the original unknown system can be estimated by reservoir computing, even with a relatively short training time series.The data-driven model reconstructs invariant sets such as fixed points and periodic orbits.Furthermore, it can reconstruct Lyapunov exponents and manifold structures between stable and unstable manifolds.
The data-driven model using reservoir computing has a neural network whose dimension is significantly higher than that of the targeted dynamics.In this current paper, we analyze the dynamics of a data-driven model using Lyapunov exponents.Pathak et al [10] found that the Lyapunov exponents of the reservoir model, constructed from time series data of the Kuramoto-Shivashinsky equation mimics those of the actual systems.
Chaotic dynamics can be quite heterogeneous, meaning that some regions have more unstable directions than others.When trajectories wander between these regions, the dynamics become complicated.A chaotic invariant set is considered heterogeneous when arbitrarily close to each point in the set, there are different periodic orbits with different numbers of unstable dimensions.This type of dynamics is referred to as hetero-chaos [20,21], and is commonly observed in high-dimensional chaotic systems.Because of the movement of a typical trajectory across periodic orbits with different numbers of unstable dimensions, the number of positive finite-time Lyapunov exponents fluctuates (for each time T > 0) [20,[22][23][24].This phenomenon is known as unstable dimension variability.For heterochaotic dynamics [20], we investigate the reconstruction of unstable dimension variability in the data-driven model by focusing on the fluctuation of the number of positive finite-time Lyapunov exponents.
Our previous study [19] demonstrated that the Lyapunov exponents of the data-driven model in the coordinates of the original system (e.g.(x, y, z) for the Lorenz-63 system) and their values agree with those of the reference system.However, this study focuses on the Lyapunov analysis of a data-driven model using the coordinates of the reservoir state vectors rather than the reference system coordinates.
Furthermore, we choose the Lorenz-96 system [25] and a fluid flow as examples of high-dimensional chaotic dynamics.After providing a brief explanation of reservoir computing in section 2, we show the results for the Lorenz-96 dynamics in section 3 and for fluid flow in section 4. Section 5 concludes our study.

Reservoir computing
Reservoir computing refers to a specific type of recurrent neural network in which the internal parameters, W in and A, are not adjusted to fit the data during the training phase [26,27].A reservoir computing model can be trained by providing an input time series, and fitting a linear function of the reservoir state vector to the output time series.Notably, no prior physical knowledge was used during the model construction process.
The data-driven model is as follows: where u(t) ∈ R M is a vector-valued variable, the component of which is denoted as an output variable; r(t) ∈ R N (N ≫ M) is a reservoir state vector.A ∈ R N×N , W in ∈ R N×M , and W * out ∈ R M×N are matrices; α (0 < α ⩽ 1) is the coefficient; ∆t is the time step of training time series.We define tanh(q) = (tanh(q 1 ), tanh(q 2 ), . . ., tanh(q N )) T , for vector q = (q 1 , q 2 , . . ., q N ) T , where T represents the transpose of the vector.
We now explain how to obtain W * out in the second equation of (1).The time development of the reservoir state vector r(l∆t) is determined by the first equation of (1), with training time series data {u(l∆t)}(−L 0 ⩽ l ⩽ L), where L 0 is the transient time and L is the time length.For given matrices A and W in , we determine W out ∈ R M×N so that the following quadratic form attains the minimum: where ∥q∥ 2 = q T q and β (β > 0) is a regularization parameter.The minimizer is set as where I ∈ R N×N is the identity matrix, δR (respectively, δU) is the matrix whose l-th column is r(l∆t) (respectively, u(l∆t)).See [28] P.140 and [29] chapter 1 for details.Note that A is randomly chosen to have the spectral radius ρ [30].See [2,10] for the details of the reservoir computing.

Lorenz-96 system
Lorenz [25,31] proposed a dissipative high-dimensional ordinary differential equation system as a model of an oscillating scalar atmospheric quantity described by where the system has cyclical symmetry, indicating that x −1 = x K−1 , x 0 = x K , and x 1 = x K+1 , and f is the forcing parameter.Figure 1.Short time inference of a variable of the Lorenz-96 system.The time series of the u1 variable inferred from the data-driven model (red) is compared with that from the actual Lorenz-96 system (blue).The two panels show the time series of the single data-driven model from different initial conditions, with the corresponding actual time series.In each case, the model trajectory exhibits behavior similar to the actual trajectory for a certain period.
The model reconstructs statistical properties and short-time trajectories.See figure 1 for the short-time trajectories of the model from two different initial conditions, both of which agree well with the actual trajectories.
Before the Lyapunov analysis, we examined the autocorrelation function to investigate the validity of a statistical property computed from the model's long trajectory.The data-driven model well reconstructed the autocorrelation function.See figure 2. This is used to identify the characteristic time length of the dynamics [34].
Lyapunov exponents were used to evaluate the degree of instability and estimate the Lyapunov dimension of a dynamical system.In some studies, the Lyapunov exponents of a data-driven model using reservoir computing were calculated in the space of N-dimensional reservoir state vector [10,11,35,36].Pathak et al [10] computed Lyapunov exponents for reservoir state vector and found that they almost coincide with those of the original system for the case of a partial differential equation, whereas only positive and neutral exponents coincide with those for the Lorenz-63 system.Figure 3 shows the Lyapunov spectrum of the data-driven model.The top several Lyapunov exponents are considered to play essential roles in the reconstruction of the dynamics.In particular, the first, the second, and the third Lyapunov exponents show good agreements, which is shown in table 2. Figure 4 shows density distributions of the first and second finite-time Lyapunov exponents with 5 finite-time for the actual Lorenz-96 system and the reservoir model.Note that the zero Lyapunov exponent corresponding to the orbit direction is not considered.Figure 5 shows time evolutions of the number of positive finite-time Lyapunov exponents for the actual Lorenz-96 system and the data-driven model.The time length of the orbit for computing the finite-time Lyapunov exponents    λ (1)  λ (2)  λ (3)   Reservoir

Macroscopic fluid flow
In this section, we constructed a data-driven model of a chaotic fluid flow and performed a Lyapunov analysis of the model.The model focused on the energy in the wavenumber space, which serves as a macroscopic variable for a fluid flow.To create the training time series of the variable, we computed a direct numerical simulation of the incompressible three-dimensional Navier-Stokes equation under periodic boundary conditions: where T = R/2πZ, ν > 0 is the viscosity parameter, π (x, t) is the pressure, and v(x, t) is the velocity.In this study, we choose ν = 0.058.We employed the Fourier spectral method [37] with 9 modes in each direction to approximate the Navier-Stokes equation.The time integration was performed using the fourth-order   Runge-Kutta method at each time step.The forcing f was input into the low-frequency variables at each time step to roughly preserve the total energy.See [18] for details.The energy function E 0 (k, t) for wavenumber k ∈ N is defined as follows:  We constructed a reservoir computing model (1) using u(t) = (E(3, t), E(3, t − ∆τ ), . . ., E(3, t − 22∆τ )) as the input variable, where ∆τ is the delay time (see table 1 for the set of hyperparameters used).The model reconstructs the statistical properties and short-time trajectories.See figure 6 for the short-time trajectories of the model from various initial forecast times, all of which agree well with the actual trajectories.The data-driven model well reconstructs the autocorrelation function.See figure 7. Figure 8 shows the Lyapunov   spectrum of the reservoir model of a fluid flow.The number of positive Lyapunov exponents is 5. Figure 9 shows the fluctuation of unstable dimensions of the reservoir model over a period.It suggests that the model has heterochaos, indicating the coexistence of different unstable dimensions.The finite-time Lyapunov dimension D KY estimated using the Kaplan-Yorke formula [38] also fluctuates between 5 ⩽ D KY ⩽ 13 for finite-time T = 20, and between 4 ⩽ D KY ⩽ 12 for finite-time T = 40.

Concluding remarks
We computed the Lyapunov spectrum of the data-driven models for high-dimensional dynamics constructed using reservoir computing.For heterochaotic dynamics, we found that the data-driven model reconstructs fluctuations in the number of positive finite-time Lyapunov exponents.

Figure 2 .
Figure 2. Autocorrelation functions for the data-driven model and the actual system of the Lorenz-96 system.The autocorrelation function of the data-driven model for the u1 variable agrees well with that of the actual system.Both functions show exponential decay.

Figure 3 .
Figure 3. Lyapunov spectrum of the Lorenz-96 system.The left panel shows the full Lyapunov spectrum and the right panel shows the spectrum for the first 20 Lyapunov exponents for the data-driven model of the Lorenz-96 system.

Figure 4 .
Figure 4. Density distributions of the top two finite-time Lyapunov exponents, except for the orbit direction, for the Lorenz-96 system.The left panel shows the first finite-time Lyapunov exponents and the right panel shows the second finite-time Lyapunov exponents, each transversal to the orbit direction.The finite-time Lyapunov exponents are averaged over a time length of 5.The distribution of the actual Lorenz-96 system with f = 6 is shown in blue, and that of the data-driven model is shown in red.

Figure 5 .
Figure 5.Time evolution of the number of positive finite-time Lyapunov exponents for the Lorenz-96 system.These are computed for the actual Lorenz-96 system (left) and the data-driven model (right) using two different finite times (5 (purple) and 10 (green)).Each of the two mainly takes the values of 1 or 2. See table 3 for the time ratios of the 1D, 2D, and 3D unstable dimensions for the finite-time Lyapunov exponents.
D k := {K ∈ R 3 |k − 0.5 ⩽ |K| < k + 0.5}.To obtain the essential dynamics of the energy function, we used a bandpass filter to eliminate high-and low-frequency fluctuations.The resulting filtered energy variable is referred to as an energy function E(k, t) throughout the study.

Figure 6 .
Figure 6.Short time inference of the energy function of a fluid flow.Time series of u1 variable inferred from the data-driven model (red) is compared with that from the actual fluid flow (blue).Each of the eight panels shows the time series of a single data-driven model from a different initial condition, with the corresponding actual time series.In each case, the model trajectory exhibits behavior similar to the actual trajectory for a certain period.

Figure 7 .
Figure 7. Autocorrelation functions for the data-driven model and the actual system for a chaotic fluid flow.The autocorrelation function of the data-driven model for the E(3) variable agrees well with that of the actual system.Both functions show exponential decay.

Figure 8 .
Figure 8. Lyapunov spectrum of the fluid flow.The Lyapunov spectrum for the data-driven model is shown in the left panel.The number of positive Lyapunov exponents is 5, and the Lyapunov dimension is 9.The right panel is the enlargement of the left panel.

Figure 9 .
Figure 9.Time evolution of fluctuations of unstable dimensions.The time evolution of fluctuations in unstable dimensions for the data-driven model of a chaotic fluid flow is shown (finite-time 20 (purple) and 40 (green)).

Table 1 .
The list of parameters and their values used in the reservoir computing.The values in each column show the parameters used to model a data-driven model for the Lorenz-96 (L-96) and a fluid model (fluid).

Table 2 .
Lyapunov exponents of the Lorenz-96 system.Lyapunov exponents of the data-driven model using reservoir computing (λ res , λ considering the lowest period of the periodic orbits.The fluctuations in the number of positive Lyapunov exponents indicate the existence of heterochaos in the data-driven model for Lorenz-96 dynamics.

Table 3 .
The time ratios of 1D, 2D, and 3D unstable dimensions for finite-time Lyapunov exponents using time series data with a length of 10 4 for Lorenz-96 dynamics.