Data-Based Reconstruction of Chaotic Systems by Stochastic Iterative Greedy Algorithm

It is challenging to reconstruct a nonlinear dynamical system when suﬃcient observations are not available. Recent study shows this problem can be solved by paradigm of compressive sensing. In this paper, we study the reconstruction of chaotic systems based on the stochastic gradient matching pursuit (StoGradMP) method. Comparing with the previous method based on convex optimization, the study results show that the StoGradMP method performs much better when the numerical sampling period is small. So the present study enables potential application of the reconstruction method using limited observations in some special situations where limited observations can be acquired in limited time.


Introduction
Dynamical systems play a key role in many scientific disciplines including physical, biological, and social systems. e dynamical models can help us analyze the way that the present state of the system depends on the past state and predict the possible state in the future. To accomplish these goals, accurate dynamical models describing these systems are required, so the reconstruction of dynamical systems from observation of time series has become a central issue in these scientific disciplines [1]. It is challenging to reconstruct a nonlinear dynamical system without prior knowledge. e delay-coordinate embedding method [2] is the traditional model-free reconstruction method for nonlinear dynamical systems.
e mathematical foundation of the method is Takens's theorem [3] which was laid out three decades ago. e reconstruction system by the delay-coordinate embedding method can be used to estimate some import quantities of the original system such as the dimension of attractor [4], the Lyapunov exponents [5], and unstable periodic orbits that constitute the skeleton of attractor [6]. With development of deep learning, some model-free methods based on neural networks are also developed [7,8]. One disadvantage of model-free methods is that they cannot give accurate dynamical model of the underlying system. e mathematical tools of the dynamical system cannot be used to analyze the behaviors of the system.
In general, it is difficult to reconstruct an accurate dynamical model of the underlying system without any prior knowledge. But the nonlinear functions describing the dynamic of the system can be approximated by power series in many famous chaotic systems, such as the Lorenz system [9], Rössler system [10], hyperchaotic Chen system [11], and Hénon map [12]. In this case, the reconstruction of an accurate dynamical model can be converted to parameter identification and can be achieved by some parameter identification methods such as the synchronization-based method [13], particle swarm optimization algorithm [14], genetic algorithm [15], and differential evolution algorithm [16]. In general, the number of unknown parameters to be estimated can be large, so large amount of observations may be required. However, it is difficult or expensive to acquire observations in some special situations. Can the system be constructed when sufficient observations are not available? In fact, though the number of unknown parameters is large, many of them are zero (the sparsity condition). So, the compressive sensing method [17,18] for sparse signal recovery can provide a solution to the problem. In [19], the authors firstly studied the reconstruction of chaotic systems by a compressive sensing method. e results show that reconstruction of the chaotic system is possible with limited observations. en reconstruction of dynamical networks [20][21][22] is also studied based on the compressive sensing method. However, we find that reconstruction success rate of the chaotic system is low when the numerical sampling period is small. e possible reason for this is that the algorithm gets trap in local optima when restricted isometry property (RIP) breaks down. In fact, the above reconstruction method is solved by the convex optimization methods [23]. In addition, there exist greedy methods for sparse signal recovery which allow faster computation compared to the convex optimization methods. In practical applications, the widely used ones are the orthogonal matching pursuit (OMP) [24], compressive sampling matching pursuit (CoSaMP) [25], and the iterative hard thresholding (IHT) [26]. Based on these methods, some stochastic greedy methods have also been developed, such as stochastic iterative hard thresholding (StoIHT), stochastic gradient matching pursuit (StoGradMP) [27], and stochastic compressive sampling matching pursuit (StoCoSaMP) [28]. In general, the randomness typically has a dramatic and positive impact on escaping from local optima. So, we will try to use the stochastic greedy method to reconstruct the chaotic dynamic systems, when the numerical sampling period is small.

Problem and Method
Let us consider chaotic dynamical systems described by where x � (x 1 , x 2 , . . . , x n ) ∈ R n is the state vector and f i (x)(i � 1, 2, . . . , n) are nonlinear functions which describe the dynamic of the system. Here, we assume that the accurate model of the system is not available, but f i (x)(i � 1, 2, . . . , n) can be approximated by power series: . . .
For many famous chaotic systems, such as the Lorenz system [9], Rössler system [10], and Hénon map [12], the nonlinear functions f i (x)(i � 1, 2, . . . , n) describing the dynamic of the system are just a part of power series. By the above assumption (2), the reconstruction of the chaotic system (1) can be converted to the identification of parameters a i l 1 ,l 2 ,...,l n . Moreover, assuming that observations of x and _ x at some time t 1 , t 2 , . . . , t k are available, and written as en, we obtain the following equation: where which is usually called design matrix or measurement matrix. Obviously, if we have k � (m + 1) n observations of x and _ x which make G a nonsingular matrix, then a i (i � 1, 2, . . . , n) can be solved from (3) easily, and the chaotic system (1) is reconstructed exactly at the same time. In general, the number of coefficients to be estimated is large, and large amount of observations are required, when n (dimension of system) is large. In some cases, however, we cannot get enough observations of x and _ x, i.e., k ≪ (m + 1) n , and the reconstruction of chaotic systems when k ≪ (m + 1) n is an interesting and significant problem. According to the conventional wisdom, the system (1) cannot be reconstructed when k ≪ (m + 1) n because the equation is underfitting, which has a nonunique solution. But most elements of a i are zero for reconstruction of the chaotic system (i.e., the a i is sparse), so the equation (3) can be solved by the compressive sensing method. In [19], the authors firstly studied the reconstruction of chaotic systems by the compressive method.
e results show that the reconstruction of the chaotic system is possible with limited observed data, when the nonlinear function described by the dynamic is unknown. e constructed systems are accurate, which can be used to predict the crisis of chaotic systems. However, we find that success rate of reconstruction is low, when the numerical sampling period is small. e reason for this may be that the algorithm gets trap in local optima that may not be global when RIP condition breaks down. Actually, the reconstruction method in [19] is to solve the following basis pursuit problem: Above basis pursuit can be recast into a linear constraint minimization problem, and corresponding MATLAB coding (l1magic) can be found at site [29]. In fact, convex optimization problem (4) is relaxation of the following nonconvex problem: where ‖·‖ 0 is l 0 norm which counts the number of nonzero elements, and K i is sparsity level. Recently, some greedy matching pursuit algorithms have been proposed to solve this problem, such as OMP [24] and CoSaMP [25]. e OMP method usually offers a much faster runtime than the convex approach but lacks strong recovery guarantees. CoSaMP offers both the advantage of fast runtime and strong recovery guarantees as a convex method [27]. e CoSaMP methods have been used remarkably in the signal processing community due to their simplicity and computational efficiency. However, we find the success rate of reconstruction is very low, when we tried to reconstruct chaotic systems by the CoSaMP method. e reason for this may due to get trap in local optima when the initial condition is not properly chosen. Basing on CoSaMP methods, some stochastic greedy methods have been developed, such as StoGradMP [27] and StoCoSaMP [28]. ough a comprehensive analysis of the phenomenon of escaping local minima in these stochastic greedy methods is difficult, balancing randomness and greediness carefully has a dramatic and positive impact on escaping from local optima. So, we will try to use the StoGradMP method to reconstruct the chaotic dynamic systems, when the numerical sampling period is small. e algorithm basing on StoGradMP is described in Algorithm 1, which consists of following steps at each iteration: (1) Randomize: Randomly select a subset S of M � 1, 2, 3, . . . , k { } with specified probability, where k is sampling number.
(2) Proxy: Form a signal proxy r j . r j � l∈S (y il − <g l , a j i > )g l is the gradient of r j � l∈S (y il − <g l , a j i > ) 2 /2, where y il is the l th element of y l , and g l is the l th row of measurement matrix of G.
(3) Identify: Locate the largest components of the proxy. e operation approx 2K i (r j )gets the indices of 2K i largest entries (in magnitude) of r j to compose of Γ. (4) Merge: Form the subspace Γ � Γ ∪ Λ, where Γis the set of newly identified components in the identify step, and Λ is the set of components that appear in the current approximation.
is is a least squares problem for present problem. According to theoretical results in [27], the above method is linear convergence in expectation when the measurement matrix Gsatisfies RIP condition. ough it is an open problem to prove whether a general matrix G satisfies the RIP condition, the numerical results show the present method is effective for reconstruction of chaotic systems. Indeed, the present algorithm is a stochastic version of CoSaMP. Comparing with original CoSaMP, the main difference is gradient computation in step proxy. e present method do not need to compute the gradient of l∈M (y il − <g l , a j i > ) 2 . Instead, it only needs to compute the gradient of l∈S (y il − <g l , a j i > ) 2 , where S is a random chosen subset of M.
is property not only makes the algorithm efficient in large scale problem, but also introduces some random perturbation which usually has positive impact on escaping from local optima. In the previous paper [28], the authors also proposed a stochastic version of the CoSaMP method. Comparing the original CoSaMP method, a greedy or a stochastic step at every iteration is randomly executed. During a stochastic step, choose 2K atoms randomly to merge into the support. In a greedy step, signal proxy is formed as the CoSaMP method. So this method is only applied for a problem where the least square loss is used.
e StoGradMP method can be applied for loss function which does not exhibit quadratic structure. is type of loss function will present in our further study, so we choose the present StoGradMP-based method.

Numerical Experiments
In this section, some numerical simulations are performed to show that the StoGradMP method is more effective than the convex method used in [19], when the numerical sampling period is small. To show the importance of random perturbation, the reconstruction results by the CoSaMP method are also presented. ere are two input parameters for StoGradMP: the sparsity level K and halting criterion. ough there are several strategies for choosing K for some particular application, these methods are not suitable for present problem. For reconstruction of chaotic systems, the number of nonzero coefficients is very small, so we run the algorithm with iteratively adding K in some range and select the best approximation. For stopping criterion, we use fixed number of iterations. Some other simple stopping criterions can be found in [25], which may also be useful in practice. In a randomize step, we simply select all the subsets of M with the same probability. In an estimate step, we use the conjuncture gradient method to solve the least square problem. To solve convex problem (4), we use the method in [23], and corresponding MATLAB coding (l1magic) can be found at site [29].
In simulation, we do some processing on original observed data. e observed data y i (i � 1, 2, . . . , n) and measure matrix G are centralized and normalized, i.e., let Mathematical Problems in Engineering where y ij is jthe lement of y i , g l is the l th column of measurement matrix G, and g lj is the jth element of g l . Because all elements of the first column of G are one, the centralization step can avoid the estimation of the constant term in function f i (x)(i � 1, 2, . . . , n). Each element of G is constructed by x l 1 1 x l 2 2 , . . . , x l n n ; thus, the difference of amplitude of g l may be large. e normalization step can remove impact of the amplitude of g l on reconstruction. Moreover, the normalization may make the matrix G ′ satisfy RIP condition, if G is constructed properly [19]. en, instead of solving equation (3), we solve where G ′ � [g 2 ′ , g 3 ′ , , g ′ (m+1) n ] and a i ′ � [a i2 ′ , a i3 ′ , . . . , a i(m+1) n ′ ] T . After the above equation is solved, a i can be acquired by where a il is the lth element of a i . a i1 the constant term of f i (x) can be computed by (3), when a il (l � 2, 3, . . . , (m + 1) n ) are known.
To measure the success of reconstruction, the following scales are defined: where n is the dimension of system, and m is an input parameter which should be large enough such that nonlinear functions f i (x) (i � 1, 2, . . . , n,) can be presented by power series (2). Obviously, if e i < ε (small constant) for all i, the unknown system is reconstructed successfully. To show the performance of each construction method, the success rate is computed. e success rate is the ratio of number of successful reconstructions to total number of independent realizations. In this section, the famous Lorenz system [9], Rössler system [10], hyperchaotic Chen system [11], and chaotic memristive circuits [30] are performed as illustrative examples. For these four systems, the highest power of f i (x) is 3, so we set m � 3. In numerical simulation, the fourthorder Runge-Kuttamethod with a time step 0.001 is adopted to create simulation data.

Example 1: Lorenz
System. e famous Lorenz system can be described as follows: (11) e Lorenz system is derived from a two-dimensional fluid layer uniformly warmed from below and cooled from above. x, y, z describe the rate of change of convection, the horizontal temperature variation, and the vertical temperature variation with respect to time. σ, r, b are positive system parameters, and the system exhibits chaotic behavior for some value of these parameters. e Lorenz equations also arise in simplified models for lasers, dynamos, thermosyphons, and so on. In this example, let σ � 10, r � 24.1, and b � 8/3. Under this group of parameters, nonlinear system (11) has chaotic attractor. For a different numerical sampling periodτ, the success rates of reconstruction by the convex method (solved by l1magic) and StoGradMP are shown in Figures 1 and 2 separately. In these two figures, we set ε � 10 − 5 , and 1000 independent realizations are computed to obtain the success rate. e Input: sparse level K i and stopping criterion Initialize: a 0 i ← 0, j ← 0, Λ ← While the stop criterion not satisfy do Randomize: randomly select a set S ⊆ M � 1, 2, 3, . . . , K { } Proxy: r j ← l∈S (y il − < g l , a j i >)g l Identify: Mathematical Problems in Engineering numerical results show that the StoGradMP method performs much better when the numerical sampling period is small. In Figure 3, the reconstruction results by the CoSaMP method are also presented. Comparing with Figure 2, one can see that random perturbation has dramatic and positive impact on improving success rate of reconstruction.

Example 2: Rössler
System. e famous Rössler system can be described as follows: Here x, y, z are the states of the system, and a, b, c are positive system parameters. Rössler designed the Rössler attractor in 1976, but the original theoretical equations were later found to be useful in modeling in chemical reactions. e original Rössler paper says the Rössler attractor was intended to behave similarly to the Lorenz attractor, but it is simpler and has only one manifold. In this example, let a � 0.1, b � 0.1, and c � 14.0. Under this group of parameters, the nonlinear system (12) has chaotic attractor. All the other parameters are the same as in example 1. e corresponding numerical results are shown in Figures 4-6. ese results also show that the StoGrandMP method is much better when the numerical sampling period is small, and random perturbation has dramatic and positive impact on improving success rate of reconstruction.

Example 3: Hyperchaotic Chen
System. e hyperchaotic Chen system can be described as follows Herex, y, z, u are the states of system, and a, b, c, d, rare positive system parameters. e positive Lyapunov exponent has an important impact on the behaviors of chaotic systems. In some applicable problems, one may face systems that have more than one positive Lyapunov exponent. is may cause the system dynamics expand in several different directions at the same time. In this example, let a � 35, b � 3, c � 12, d � 7, and r � 0.5. Under this group of parameters, the nonlinear system (13) has hyperchaotic attractor (two positive Lyapunov exponents). All the other parameters are the same as in example 1. e corresponding numerical results are shown in Figures 7-9. ese results also show that the StoGrandMP method is much better when the numerical        sampling period is small, and random perturbation has dramatic and positive impact on improving success rate of reconstruction.

Example 4: Chaotic Memristic Circuits.
In this example, we consider a simplest memristic circuit [30] which consist of an inductor, capacitor, and nonlinear memristor. e system can be described as follows: where x is the voltage across capacitor, y is the current through inductor, and z is the internal state of memristor. e parameters are C � I s C n T s , L � L n T s /I s , β � β 5kpot /R, and α � 1/T s C f α 10kpot . We choose I s � 10000, T s � 10 5 C n � 1nF, L n � 300mH, R � 1KΩ, C f � 10nF,β 5kpot � 1.5KΩ, and α 10kpot � 5KΩ. Using these component values from the circuit, we get the following values for parameters in the system (14): L � 3.0, C � 1.0, β � 1.5, and α � 0.6. Under this group of parameters, nonlinear system (14) has chaotic attractor. All the other parameters are the same as in example 1. e corresponding numerical results are shown in Figures 10-12. ese results also show that the StoGrandMP method is much better when the numerical sampling period is small, and random perturbation has dramatic and positive impact on improving success rate of reconstruction.

Conclusions
Reconstruction of the nonlinear dynamical system from observed time series is interesting and significant in many scientific disciplines. In general, it is difficult to reconstruct the exact dynamical model when sufficient observed data are not validated. Recently, some studies show that this problem can be solved by the compressive sensing method which is widely used in sparse signal recovery. But we find that the reconstruction success rate of previous reconstruction methods is low, when the numerical sampling period is small. e reason for this may due to get trap in local optima that may not be global when RIP breaks down. In this paper, we try to solve the problem by stochastic gradient matching pursuit developed recently. In general, randomness has positive impact on escaping from local optima. Numerical simulation results show the stochastic gradient matching pursuit method is much more effective than the previous reconstruction method when the numerical sampling period is small. So the present study enables potential application of the reconstruction method based on limited observations in some special situations where limited data just can be acquired in limited time.

Data Availability
No data were used to support this study.