Dynamical system analysis of a data-driven model constructed by reservoir computing
Miki U Kobayashi, Kengo Nakai, Yoshitaka Saiki, Natsuki Tsutsumi
aa r X i v : . [ m a t h . D S ] F e b Dynamical system analysis of a data-driven model constructed by reservoir computing
Miki U. Kobayashi
Faculty of Economics, Rissho University, Tokyo 141-8602, Japan
Kengo Nakai
Faculty of Marine Technology, Tokyo University of Marine Science and Technology, Tokyo 135-8533, Japan
Yoshitaka Saiki
Graduate School of Business Administration, Hitotsubashi University, Tokyo 186-8601, Japan
Natsuki Tsutsumi
Faculty of Commerce and Management, Hitotsubashi University, Tokyo 186-8601, Japan (Dated: March 1, 2021)This study evaluates data-driven models from a dynamical system perspective, such as unstablefixed points, periodic orbits, Lyapunov exponents, manifold structures, and statistical values. Thesedynamical characteristics can be reconstructed much more precisely by a data-driven model thanby computing directly from training data. With this idea, we predict the laminar lasting timedistribution of a particular macroscopic variable of chaotic fluid flow, which cannot be calculatedfrom a direct numerical simulation of the Navier–Stokes equation because of its high computationalcost.
Introduction.
Reservoir computing, a brain-inspiredmachine-learning technique that employs a data-drivendynamical system, is effective in predicting time se-ries and frequency spectra in chaotic behaviors, includ-ing fluid flow and global atmospheric dynamics [1–10].Pathak et al. [3] examined the Lorenz system and theKuramoto–Sivashinsky system and reported that thedata-driven model obtained from reservoir computingcould generate an arbitrarily long time series that mimicsthe dynamics of the original systems.The extent to which a data-driven model using reser-voir computing can capture the dynamical properties oforiginal systems should be determined. Lu et al. [11]reported that a data-driven model has an attractor sim-ilar to that of the original system under an appropri-ate choice of parameters. Nakai and Saiki [12] confirmedthat a single data-driven model could infer the time se-ries of chaotic fluid flow from various initial conditions.Zhu et al. [13] identified some unstable periodic orbits ofa data-driven model through delayed feedback control.They suggested that a data-driven model could recon-struct the attractor of the original dynamical system.This Letter clarifies that a data-driven model usingreservoir computing has richer information than that ob-tained from a training data, especially from dynamicalsystem point of view, suggesting that dynamical prop-erties of the original unknown dynamical system canbe estimated by reservoir computing from a relativelyshort time series. Besides the invariant sets, such asfixed points and periodic orbits, the dynamical proper-ties, such as Lyapunov exponents and manifold struc-tures between stable and unstable manifolds, can be re-constructed by the data-driven model through reservoircomputing, even if the system does not have structural stability.We mainly deal with the Lorenz system [14]: dxdt = 10( y − x ) , dydt = rx − y − xz, dzdt = xy − z, (1)and will be denoted as the actual Lorenz system in thisLetter. A data-driven model is constructed from a shorttime training data created from (1), the method of whichis explained later. Two different parameter values of r are considered. One of the parameters ( r = 28) has hy-perbolic dynamics, whereas the other ( r = 60) generatesdynamics with tangencies between stable and unstablemanifolds [15]. The latter property is one of the two pri-mary sources for the breaking structural stability [16],which often appears in the real-world physical phenom-ena. As an application of the obtained knowledge, thisstudy examines high-dimensional chaotic fluid flow to de-termine if the laminar lasting time distribution can bepredicted using the data-driven model constructed fromshort training time-series data. Reservoir computing.
A reservoir is a recurrentneural network whose internal parameters are not ad-justed to fit the data in the training process [17, 18].The reservoir can be trained by feeding it an input timeseries and fitting a linear function of the reservoir statevariables to the desired output time series. We do notuse a physical knowledge in constructing a model. Thedata-driven model using reservoir computing we study isthe following: ( u ( t ) = W out r ( t ) , r ( t + ∆ t ) = (1 − α ) r ( t ) + α tanh( Ar ( t ) + W in u ( t )) , (2) z xr=28, training 0 20 40 60 80 100 -30 -20 -10 0 10 20 30 z xr=60, training FIG. 1.
Poincar´e section-like plots ( r = 28 (top) and 60(bottom)). The sets of points ( x, z ): along a trajectory ofthe data-driven model using reservoir computing and alonga long trajectory of the actual Lorenz system and along ashort trajectory used for the training data are plotted when | x − y | < ǫ p , where ǫ p = 0 .
05. The time lengths of the threetrajectories are T = 10 , 10 , and 5000, respectively. where u ( t ) ∈ R M is a vector-valued variable, the com-ponent of which is denoted as an output variable; r ( t ) ∈ R N ( N ≫ M ) is a reservoir state variable; A ∈ R N × N , W in ∈ R N × M , and W out ∈ R M × N are matrices; α (0 < α ≤
1) is a coefficient; ∆ t is a time step. We de-fine tanh( q ) = (tanh( q ) , tanh( q ) , . . . , tanh( q N )) T , fora vector q = ( q , q , . . . , q N ) T , where T represents thetranspose of a vector. We determine W out using train-ing time series data { u ( l ∆ t ) } L − l =0 . More details can befound elsewhere [6].In this Letter we evaluate a data-driven model (2) con-structed using short training time series data from a dy-namical system perspective. The main focus is on theproperties in the space of output variables (correspond-ing to x , y , and z for the case of the Lorenz system),which compare them with those of the actual system. Poincar´e section-like plots.
The Poincar´e sectionof the data-driven model of the Lorenz system has beenstudied [3]. We compare the shape and size of the attrac-tor of a data-driven model (2) with those of the attractorof the actual Lorenz system (1), and also with those ofthe set of points along the training time series data. Fig-ure 1 presents their Poincar´e section-like plots for r = 28and 60. For each of the two parameter cases, a set oftrajectory points generated from the data-driven modelseem to coincide with the chaotic attractor of the actualLorenz system. Furthermore, the data-driven model hasan attractor which is significantly larger than the set oftraining data used to construct the model. Density distribution.
The density distribution of x variable along a trajectory of the data-driven model ispresented in Fig. 2. We compare the distribution withthat obtained from the trajectory of the actual Lorenz P D F x Φ res Φ actual Φ training -0.004-0.002 0 0.002 0.004-20 -10 0 10 20 d i ff e r en c e x Φ res − Φ actual Φ training − Φ actual P D F x Φ res Φ actual Φ training -0.003-0.002-0.001 0 0.001 0.002 0.003 -30 -20 -10 0 10 20 30 d i ff e r en c e x Φ res − Φ actual Φ training −Φ actual FIG. 2.
Density distributions of a variable ( r = 28 (top)and r = 60 (bottom)). The density distribution of the x variable calculated from a length T = 10 trajectory ofthe data-driven model by reservoir computing (Φ res ) is plot-ted together with that computed from the length T = 10 long trajectory of the actual Lorenz system (Φ actual ) and thelength T = 5000 short trajectory (Φ training ) used as train-ing data for constructing the data-driven model. The differ-ence in the distributions are shown on the right. The av-erage x and the standard deviation σ of the density distri-bution are as follows; ( x, σ ) = ( − . , . . , . r = 28; ( x, σ ) = ( − . , . − . , . r = 60. R | Φ res − Φ actual | dx/ R | Φ training − Φ actual | dx ≈ / r = 28,and ≈ / r = 60. system (1) and that calculated directly from the trainingdata. The distribution of the actual Lorenz system canbe captured by employing the data-driven model. Re-markably, the distribution with a singular structure [19]in r = 60 can be recovered using the data-driven model. Fixed points and their stabilities.
Fixed points,which are fundamental structures of dynamical systems,are examined. We identify fixed points in the space of theoutput variables directly, even though they were identi-fied through the fixed points in the space of the reservoirstate variables using the directional fibers method [20].We also study the stability of each fixed point in thespace of output variables. For the data-driven model weconsider a point x ∗ = ( x ∗ , y ∗ , z ∗ ) as a fixed point, whenthe following condition is satisfied: δ = max n ∈ [0 ,n ] k x ∗ − ψ x ∗ ( n ∆ t ) k l < ǫ for some ǫ sufficiently small and forsome n sufficiently large, where ψ x ∗ ( n ∆ t ) is the pointiterated n times from x ∗ by the data-driven model with L res R res O res L actual R actual O actual x ∗ − .
47 8 .
50 0 . − .
49 8 .
49 0 . y ∗ − .
47 8 .
50 0 . − .
49 8 .
49 0 . z ∗ .
04 27 .
01 0 .
54 27 .
00 27 .
00 0 . .
09 + 10 . i .
10 + 10 . i .
67 0 .
09 + 10 . i .
09 + 10 . i . . − . i . − . i − .
66 0 . − . i . − . i − . − . − . − . − . − . − . TABLE I.
Coordinates and eigenvalues of the Jacobianmatrix at each of the three fixed points. L res , R res ,and O res are fixed points of the data-driven model, whereas L actual , R actual , and O actual are fixed points of the actualLorenz system with r = 28. The coordinates ( x ∗ , y ∗ , z ∗ ) andthe eigenvalues (Λ , Λ , Λ ) of the Jacobian matrix at eachfixed point of the data-driven model are close to those of thecorresponding fixed point of the actual Lorenz system.FIG. 3. Fixed points.
The three fixed points ( x ∗ , y ∗ , z ∗ ) ofthe data-driven model (left) and those of the actual Lorenzsystem (right) are plotted together with the trajectory pointswith the time length T = 10 . The three fixed points of thedata-driven model are close to those of the actual Lorenz sys-tem, despite the fixed points being outside the training data,which is part of the actual trajectory. See the coordinates andthe eigenvalues of the Jacobian matrix at each fixed point inTable I. the time step ∆ t . For the computation of a trajectoryfrom a given point x ∗ of the data-driven model, reservoirstate vector r (0) is determined to correspond to x ∗ bythe pre-iterates. The echo state property [17] in whichour choice of parameters in the data-driven model (2)is satisfied guarantees that for each x ∗ , the correspond-ing reservoir state vector is determined uniquely. Table Ilists the obtained coordinates of the three fixed points, L res , R res and O res , together with those of the actualLorenz system. We fix ( ǫ , n ) = (0 . , L res and R res , and ( ǫ , n ) = (1 ,
30) for O res . Figure 3 showsthe fixed points together with the trajectory points. Ta-ble I also lists the eigenvalues of the Jacobian matrix ateach fixed point. The values are obtained from the esti-mated formula of the Jacobian matrix described later forcalculating the Lyapunov exponents and vectors. Periodic trajectory.
Periodic orbits are also thefundamental structures of dynamical systems. We con-firm that the data-driven model of discrete time has aperiodic orbit-like trajectory that travels near the corre-sponding periodic orbit of the actual Lorenz system (1)of continuous time. We call { ψ x (0) ( n ∆ t ) } n ∈ [0 ,n p ] a peri- -20-10 0 10 20 -30-15 0 15 30 0 10 20 30 40 50 reservoiractualx yz -20-15-10-5 0 5 10 15 20 0 1 2 3 4 5 6 x t reservoiractual FIG. 4.
A periodic orbit-like trajectory.
A periodicorbit-like trajectory obtained from the data-driven model isplotted together with the corresponding periodic orbit (period T p =5 . r = 28(left), and their time developments of the x variable (right). r λ (1)res λ (2)res λ (3)res D KY res λ (1)actual λ (2)actual λ (3)actual D KY actual
28 0 .
905 0 . − .
574 2 .
06 0 .
905 0 . − .
571 2 . .
397 0 . − .
065 2 .
09 1 .
403 0 . − .
070 2 . TABLE II.
Lyapunov exponents and Lyapunov dimen-sions.
Lyapunov exponents of the data-driven model usingreservoir computing ( λ (1)res , λ (2)res , λ (3)res ) and those of the actualLorenz system ( λ (1)actual , λ (2)actual , λ (3)actual ) are listed. The valuesare computed using the four-stage and fourth-order Runge–Kutta method with time step 4∆ t from the points along anorbit trajectory and the estimated Jacobian matrices. TheLyapunov dimensions D KY res for the data-driven model and D KY actual for the actual Lorenz system are estimated from theKaplan–Yorke formula [21]. odic orbit-like trajectory, if the following value is suffi-ciently small: δ p = max n ∈ [0 ,n p ] k x ( n ∆ t ) − ψ x (0) ( n ∆ t ) k l , where x (0) is a point on the corresponding periodic or-bit with period T p of the actual Lorenz system (1), n p is the smallest integer satisfying n p ∆ t ≥ T p . Among theperiodic orbit-like trajectories of the data-driven modelcorresponding to the 50 periodic orbits with low periods, δ p < . δ p < . δ p among the 50periodic orbits with low periods. Lyapunov exponents and Lyapunov vectors.
The Lyapunov exponents are used to evaluate the degreeof instability and estimate the Lyapunov dimension of adynamical system. In some studies, the Lyapunov ex-ponents of a data-driven model by reservoir computingwere calculated in the space of N -dimensional reservoirstate variables [3, 4, 22, 23]. To the best of the authors’knowledge, they have not been computed only in a spaceof output variables.In this Letter, an attempt is made to compute Lya-punov exponents in the space of output variables. Theresults are compared with those of the actual Lorenz sys-tem (1) for two sets of parameters. In this computa- P D F anglereservoiractual 0 0.01 0.02 0.03 0.04 0 10 20 30 40 50 60 70 80 90 P D F anglereservoiractual 0 0.01 0.02 0.03 0.04 0 10 20 30 40 50 60 70 80 90 P D F anglereservoiractual 0 0.01 0.02 0.03 0.04 0 10 20 30 40 50 60 70 80 90 P D F anglereservoiractual FIG. 5.
Distribution of the angle between stable andunstable manifolds ( r = 28 (left) and r = 60 (right)). Thedensity distribution of the manifold angles (degree) at pointsalong a trajectory is shown for a data-driven model usingreservoir computing together with that of the actual Lorenzsystem. tion the discrete time trajectory points of a data-drivenmodel are considered samples of the continuous time tra-jectory. Table II shows the agreement of the Lyapunovexponents and the Lyapunov dimensions. We also com-pute (co-variant) Lyapunov vectors, which measure thedegree of hyperbolicity by calculating the angle betweenthe stable and unstable manifolds at some trajectorypoint [24]. The Lyapunov exponents and the Lyapunovvectors for the data-driven model are obtained from theJacobian matrices estimated numerically from the trajec-tory points in the space of output variables [25]. Manifold structure and Tangency.
Using thecomputed Lyapunov vectors we investigate the manifoldstructures of the data-driven model, particularly the de-gree of hyperbolicity and the tangencies between the sta-ble and the unstable manifolds. We consider the Lorenzsystem of r = 28 without tangencies and of r = 60 withtangencies for the comparison [15]. Figure 5 shows theprobability density function of an angle between a tan-gent vector of a stable manifold and that of an unstablemanifold along an orbit trajectory for each of the actualsystem and the data-driven model. For each case of theparameters, r = 28 and r = 60, the angle distributionsare quite similar in shape, indicating that the data-drivenmodel can reconstruct the manifold structures. More-over, Fig. 5 (right) suggests that the data-driven modelcan represent a non-hyperbolic structure with tangenciesbetween stable and unstable manifolds. Laminar lasting time distribution of chaoticfluid flow.
Applying the knowledge obtained abovefor the Lorenz system, we study laminar lasting timedistribution of a macroscopic quantity of chaotic fluidflow in three dimensions under periodic boundary con-ditions [6, 12]. The distribution shown in Fig. 6 is gen-erated from the very long trajectory of the data-drivenmodel constructed by reservoir computing with a rela-
500 1000 1500 2000 2500 3000 3500laminar lasting timereservoirtraining 10
500 1000 1500 2000 2500 3000 3500laminar lasting timereservoiractual (Navier-Stokes)
FIG. 6.
Laminar lasting time distribution of a fluidflow.
The laminar lasting time normalized distribution of acertain energy function E ( t ) of a fluid variable (correspondingto ˜ E (3 , t ) in [6]) estimated from the trajectory of the data-driven model ( T = 2 × ) is shown together with that fromthe short time training time series data ( T = 0 . × T )(left) and with that by a long time actual time series data( T = 0 . × T ) (right). The training data and actual data arecalculated from a direct numerical simulation of the Navier–Stokes system. E ( t ) is the normalized variable (average 0,standard deviation 1) and we consider the state is laminarwhen | E ( t ) | < . tively low computational cost. The detailed macroscopicdynamical structures can be determined using the data-driven model constructed from time series data withoutreferring to microscopic behaviors. We hardly obtainthese structures from a direct numerical simulations ofthe Navier–Stokes equation because of the high compu-tational cost. See the discrepancy in the distributions inFig. 6 (right). Acknowledgements.
YS was supported by the JSPSKAKENHI Grant No.17K05360 and No.19KK0067. KNwas supported by the Project of President DiscretionaryBudget of TUMST. Part of the computation was sup-ported by JHPCN (jh200020), HPCI (hp200104), andthe Collaborative Research Program for Young · WomenScientists of ACCMS and IIMC, Kyoto University. [1] D. Verstraeten, B. Schrauwen, M. D’Haene, and D. A.Stroobandt, Neural Networks , 391 (2007).[2] Z. Lu, J. Pathak, B. Hunt, M. Girvan, R. Brockett, andE. Ott, Chaos , 041102 (2017).[3] J. Pathak, Z. Lu, B. Hunt, M. Girvan, and E. Ott, Chaos , 121102 (2017).[4] J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, Phys.Rev. Lett. , 024102 (2018).[5] P. Antonik, M. Gulina, J. Pauwels, and S. Massar,Phys. Rev. E , 012215 (2018).[6] K. Nakai and Y. Saiki, Phys. Rev. E , 023111 (2018).[7] T. Arcomano, I. Szunyogh, J. Pathak, A. Wikner,B. R. Hunt, and E. Ott, Geophys. Res. Lett. ,e2020GL087776 (2020). [8] S. Pandey and J. Schumacher, Phys. Rev. Fluids ,113506 (2020).[9] Y. Huang, L. Yang, and Z. Fu,Earth System Dynamics , 835 (2020).[10] L.-W. Kong, H.-W. Fan, C. Grebogi, and Y.-C. Lai, Phys.Rev. Res. , 013090 (2021).[11] Z. Lu, B. R. Hunt, and E. Ott, Chaos , 061104 (2018).[12] K. Nakai and Y. Saiki, Discrete Contin. Dyn. Syst. S ,1079 (2021).[13] Q. Zhu, H. Ma, and W. Lin, Chaos , 093125 (2019).[14] E. Lorenz, J. At. Sci. , 130 (1963).[15] Y. Saiki and M. U. Kobayashi,JSIAM Lett. , 107 (2010).[16] C. Bonatti, L. D´ıaz, and M. Viana, Dynamics BeyondUniform Hyperbolicity (Springer-Verlag, Berlin, 2005).[17] H. Jaeger, GMD Report , 13 (2001).[18] H. Jaeger and H. Haas, Science , 78 (2004). [19] S. M. Zoldi, Phys. Rev. Lett. , 3375 (1998).[20] S. Krishnagopal, G. Katz, M. Girvan, and J. Reggia, in (IEEE, 2019) pp. 1–8.[21] J. Kaplan and J. Yorke, in Functional differential equa-tions and approximation of fixed points (Proc. SummerSchool and Conf., Univ. Bonn, Bonn, 1978) , LectureNotes in Math., Vol. 730 (Springer, Berlin, 1979) pp.204–227.[22] C. Gallicchio, A. Micheli, and L. Silvestri, Neurocomput-ing , 34 (2018).[23] V. Pyragas and K. Pyragas, Phys. Lett. A , 126591(2020).[24] F. Ginelli, P. Poggi, A. Turchi, H. Chat´e, R. Livi, andA. Politi, Phys. Rev. Lett.99