An improved sparse identification of nonlinear dynamics with Akaike information criterion and group sparsity.

A crucial challenge encountered in diverse areas of engineering applications involves speculating the governing equations based upon partial observations. On this basis, a variant of the sparse identification of nonlinear dynamics (SINDy) algorithm is developed. First, the Akaike information criterion (AIC) is integrated to enforce model selection by hierarchically ranking the most informative model from several manageable candidate models. This integration avoids restricting the number of candidate models, which is a disadvantage of the traditional methods for model selection. The subsequent procedure expands the structure of dynamics from ordinary differential equations (ODEs) to partial differential equations (PDEs), while group sparsity is employed to identify the nonconstant coefficients of partial differential equations. Of practical consideration within an integrated frame is data processing, which tends to treat noise separate from signals and tends to parametrize the noise probability distribution. In particular, the coefficients of a species of canonical ODEs and PDEs, such as the Van der Pol, Rössler, Burgers' and Kuramoto-Sivashinsky equations, can be identified correctly with the introduction of noise. Furthermore, except for normal noise, the proposed approach is able to capture the distribution of uniform noise. In accordance with the results of the experiments, the computational speed is markedly advanced and possesses robustness.

Introduction

Across a range of engineering fields, if the form of the governing equations is clearly known, it is possible to analyze, forecast and control dynamics. However, because it is impossible to obtain all observed data due to the complexity of the nonlinear dynamics, solving the governing equations of nonlinear dynamics has become a challenging research task. Although models are initiated from first principles [1], as the latest motivation toward machine learning, the highlighted developments posit on data-driven model discovery [2], with a much broader class of methods, including Koopman mode decomposition [3-5], dynamic mode decomposition (DMD) [6, 9], neural networks [7, 8, 12, 31], equation-free modeling [10], and nonlinear Laplacian spectral analysis [11]. Advancements regarding parsimonious models, which strike a balance between accuracy and complexity, are particularly notable.

In this line of evolution, innovative achievements utilizing symbolic regression [12] were able to directly realize nonlinear dynamics from data. Recently, sparsity-promoting technologies via regularization have been employed to robustly find the sparse representation in the space of the potential functions for nonlinear dynamics, incorporating the least absolute shrinkage and selection operator (LASSO) [13], sequentially thresholded least squares (STLSQ) [14, 24], stagewise regression [15], basis pursuit denoising (BPDN) [16], and matching pursuits [17].

Superior model discovery demands a severe technology for validation. Model selection [18-23] built on Occam’s Razor filters parsimonious, explanatory and generalizable models from an unknown feature library via Pareto analysis. Furthermore, a well-selected model can be evaluated by information loss between observed data and model-generated data via information criteria, such as the Akaike information criterion (AIC) [19] or Bayesian information criterion (BIC) [20].

The sparse identification of nonlinear dynamics (SINDy) [24] algorithm is adopted in the context of data-driven model discovery and sparse regression. Its developments are discussed in Sect. 1.1, and the contributions of this paper are outlined in Sect. 1.2.

Overview of SINDy

As the SINDy algorithm directly obtains the governing equations of the nonlinear dynamics through partially known data, its extensions have existed across various engineering and science fields due to its particular architecture. It frames the problem for solving governing equations of nonlinear dynamics as the sparse regression framework via a nonlinear, predefined function library. Table 1 outlines a review of the SINDy algorithm.

Table 1

Overview of SINDy

Classification	Literature	Method	Modification/Field	Objective
Candidate function library	Brunton. et al. [25]	SINDYc (SINDy with control)	Include external actuation and control signals	Generalize the SINDy algorithm to include inputs and control
	Mangan. et al. [26]	Implicit-SINDy	Rational function nonlinearities with cross terms	Infer networked nonlinear dynamical systems
	Quade et al. [29]	Abrupt-SINDy	Either additional, deleted, or modified model terms	Learn parsimonious models of a system in response to abrupt changes
	Kaiser et al. [47]	SINDy-MPC (model prediction control)	Include the effects of actuation	Enhance data-driven control of physical systems
	Champion et al. [31]	Autoencoder-SINDy	Include intrinsic coordinates	Discover a coordinate transformation into a reduced space
	Bramburger et al. [30]	UPO-SINDy (unstable periodic orbits)	Contain parameters for Poincaré mappings	Produce parameter dependent Poincare mappings to stabilize UPOs
	Kaptanoglu et al. [32]	Constrained-SINDy	Incorporate known physical laws	Build a connection to the large-scale Galerkin in fluid mechanics
	Wei [33]	Integral-SINDy	Contain the integral form of state equations	Identify model structure and parameters of nonlinear ODEs
	Hoffmann et al. [48]	Reactive-SINDy	Incorporate concentration time series	Estimate a parsimonious reaction network
Sparse regression	Champion et al. [35]	SR3(sparse relaxed regularized regression)	Additional regularizer for auxiliary variable	Solve a relaxation of the regularized problem
	Chu et al. [49]	Lagrangian-SINDy	Transform the naive formulation as the Legendre form	Learn interpretable and sparse formulations of overall energy and the Lagrangian
	Dai et al. [50]	Extended-SINDy	Stepwise Sparse Regressor (SSR)	Detect the dynamic behavior from data
	Fasel et al. [36]	Ensemble-SINDy	Incorporate the bootstrap aggregating (bagging)	Robustify the SINDy algorithm
	Hirsh et al. [37]	UQ-SINDy (uncertainty quantification)	Add the spike-and-slab prior and the regularized horseshoe prior	Promote robustness against noise measurements for out-of-sample forecast
	Ram et al. [38]	Bayesian-SINDy	Combine with the neuronized priors	Handle parameter and structural model uncertainties
Modified structure	Brunton et al. [51]	HAVOK (Hankel alternative view of Koopman analysis)	Integrate with delay embedding and Koopman theory	Decomposition of chaotic dynamics into linear models using forcing actuation
	Loiseau et al. [52]	SINDy-KNN (K-nearest neighbors)	Incorporate the KNN to map into the low-dimensional space	Yield interpretable nonlinear models from measurement data
	Kaheman et al. [53]	Modified-SINDy	Integrate with the automatic differential for neural networks	Denoise the data and parametrize the noise probability distribution
	Dam et al. [54]	Bifurcation-SINDy	Combine with the bifurcation analysis	Learn a underlying model from simulations of a convection system
Mathematical theory	Zhang et al. [14]	SINDy	Convergence performance	Make a proof that SINDy approximates to local minimizers of an unconstrained l0-penalty least squares issue
	de Silva et al. [45]	PySINDy (Python for the SINDy)	Code package for the SINDy algorithm	A Python package for SINDy
	Messenger et al.[55]	Weak-SINDy	Formulate the problem as the weak form	Data-driven reconstruction of model coefficients from measurements
	Wu et al. [44]	SINDy	Error processing from an optimize perspective	Improve the accuracy of data
	Cortiella et al. [56]	SINDy	Consider regularized weight	Determine governing equations of nonlinear dynamics from noisy observations
Application	Jiang et al. [40]	SINDy-LM (Levenberg–Marquardt)	Disease (especially for COVID-19)	Analyze the infectious disease in both Chinese mainland and other countries
	Bhadriraju et al. [57]	SINDy	Chemistry	Achieve the adaptive model identification
	Fukami et al. [42]	CNN-LSTM-SINDy	Fluid flows	Investigate influence on the parameter choice
	Loiseau et al. [43]	SINDy	Thermal convection dynamics	Predict a bifurcation of the high-dimensional system
	Jadhav et al. [58]	SINDy	Video process	Extract related intelligence from actual videos of strongly stochastic systems

As seen in Table 1, the extensions of the SINDy algorithm have been classified by modified architecture and distinct applications. Control terms or external forcing are introduced to augment the potential nonlinear function library [25], even for rational nonlinearities [26], integrating the alternating direction method [27] to solve the implicit function in the null space. In addition, since Lorenz-like systems [28] are exceedingly sensitive to initial conditions, the dismissed terms or discontinuous points [29] may be considered.

More recently, Poincaré mappings [30], intrinsic coordinates [31], constrained physical laws [32], and integral forms of the SINDy model [33] have been incorporated into the candidate function library to obtain well-selected models.

The process of the sparse regression problem can be recast as a convex problem. Therefore, convex relaxation regularization [34, 35] is embedded in the native SINDy frame for parameter estimation except for its inadequate robustness analysis. Bootstrap aggregating is utilized to robustify the SINDy algorithm [36]. Furthermore, spike-and-slab prior, regularized horseshoe prior [37] or Bayesian inference [38] is considered to promote robustness for limited data. Additionally, clustering is integrated into the SINDy model to identify the turning point for infectious disease [39, 40]. The underlying idea behind the SINDy may be deeply exploited to exceedingly extended areas, such as biology [26], optics [41] and physical fluids [42, 43]. Additionally, Zhang et al.[14] provided the proof of the convergence of SINDy, and optimization techniques have been utilized to reduce the derivative error in the original SINDy algorithm [44]. Afterward, the Python package for SINDy was explored [45, 46].

Contributions of this work

The proposed technique in this work focuses on the improved framework of the SINDy algorithm to identify the coefficients of a class of typical nonlinear dynamics and the noise distribution. The principle dedications and innovations of this paper are illustrated as follows.

A variant method of the SINDy is developed with the Akaike information criterion to directly identify governing equations from data.

Model selection is integrated with the original SINDy framework to identify the number of terms of nonlinearities in dynamics.

The group sparsity is embedded in the sparse regression to learn the coefficients of PDEs due to the increased complexity of nonlinear dynamics.

Considering the situation in which noise exists [53], input data may be polluted by noise or other perturbation elements. To avoid generating error for derivative estimation, it is required to divide the observations into noisy estimations and estimated measurement data to simultaneously denoise the measurements and to determine the probability distribution of the noise.

Structure of the article

The remaining structure is arranged as follows: Sect. 2 provides the overall framework of the proposed method. Section 3 elaborates the mathematical theory for the proposed method, including the SINDy, the sparse model selection with the Akaike information criterion (SINDy-AIC), and the group sparsity identification. Relevant results for experiments are displayed in Sect. 4. The discussion and further research are described in Sect. 5. Finally, the meaning of symbols appeared in this work is synthesized in “Appendix A”.

Framework of the proposed approach

In this work, the SINDy algorithm is extended to integrate with the Akaike information criterion (AIC) and the group sparsity to identify a class of typical ODE and PDE models. Noise with normal and uniform distributions is introduced to test the robustness of the proposed approach. The framework and the flowchart of the suggested method in this paper are displayed in Figs. 1 and 2. The primary process is described as follows:

Fig. 1

Schematic of the proposed technique process

Fig. 2

Overall calculation flow of the adopted method

Step 1: Sparse identification of nonlinear dynamics (SINDy). A set of the time series data extracted from several ODE systems, prescribed as , is imported into the SINDy model and a predefined library is constructed on the basis of a priori physical information. Next, the sparsity threshold is utilized to iteratively regularize the nonzero entries in the matrix to obtain a sparse model .

Step 2: Sparse model selection. The SINDy algorithm provides a combinatorially large amount of candidate models in the training set . Then the Akaike information criterion (AIC) is incorporated into the SINDy to select the optimal model by hierarchically ranking the AIC scores with the inclusion of different categories of support values in the validation set .

Step 3: Sequentially thresholded group ridge regression. In the primitive SINDy architecture, the sequentially thresholded least squares is employed to learn the sparse terms incorporated in each candidate function library . Due to the complex structure of the parametric PDE dynamics, the group sparsity is used to learn the dependently parametric coefficients with function libraries and sparse vectors bundled into a group , where the ridge regression is leveraged to threshold the coefficients .

Step 4: Automatic noise identification. In general, the robustness of this approach should be considered. First, the noise estimation is achieved by presmoothing the noisy observation and obtaining the estimated value of the clean measurements and its derivative . Finally, the technique makes use of optimization and the SINDy algorithm to synchronously recognize the distribution of the noise and to denoise the observations.

Methods

SINDy

Here, we consider the form of equations for dynamical systemsThe function f(x(t)) denotes the governing equations of nonlinear dynamical systems, and the vector indicates the shape of a system at time t.

Actually, there exist only a few important terms in governing equations for physical systems of interest, so that the right-hand side of the equations is assumed to be sparse in the space of potential functions. To identify the form of Equations f from measurements, observations x(t) extracted from dynamical systems in time history and sampled at several times are arranged in the following matrix:

However, as the observed data will be employed to learn a model that can capture the entire trajectory of the motion of systems, the amount of observations should not be extremely small. In this way, the data can be augmented by coordination transformation or matrix transformation [59]. Here, the derivative calculated by numerically approximate or total variation regularization is applied as the augmentation of the data dimension, which can be represented byRealistically, X and are often contaminated with noise. Relying on the noise, it may be indispensable to filter X and . Otherwise, the vector of sparse coefficients does not hold exactly. Total variation regularization [60] is utilized to denoise the derivative to counteract the differentiation error.

As control variables constituting dynamics can be established by individual customs, it is promising to make an assumption that the control variables can be deemed known terms. Accordingly, a candidate function, which includes all the possible nonlinearities composed of control variables in the form of permutations, including constant, polynomial or trigonometric terms, can be constructed. The specific form can be defined by the form of the governing equations of nonlinear dynamics, where denotes a higher order polynomial term. For instance, represents the quadratic nonlinearities with respect to the state x, the exact form of which is shown as

At this point, Eq. (1) can be transformed into the form of Eq. (6). In this case, it is assumed that the left-hand side of the equation and the first term on the right-hand side of the equation represented as a self-set library are known, and the quantity to be obtained is the matrix of coefficients that characterize the control terms

Simply put, the problem for obtaining specific governing equations is converted to a requirement for the vector of sparse coefficients. Naturally, symbolic regression is exploited to simplify the issue. Therefore, Eq. (6) is substituted with Eq. (7).

Here, which is a numerical matrix, is a vector of the symbolic function of elements of x with respect to . Naturally, the final form of the equations to be solved is indicated by where the superscript T denotes the matrix transpose.

Each column of Eq. (6) demands a different optimization to obtain the sparse vector of coefficients for the lth row equation. To mitigate computational errors caused by inconsistent data dimensions or excessively small entries of X, columns of may be normalized. In view of presupposing that most coefficient matrices are sparse in an appropriate basis, a sparse solution to an overdetermined system with noise should be sought. Here, the sequentially thresholded least squares [24] is used, the particular form of which is given bywhere is the sparsity knob and the R(·) denotes the regularized function.

Once the coefficient matrix has been acquired, the identified model through SINDy is obtained, which can also be considered a reconstruction or prediction of the ground truth system. The detailed flow is shown in Figs. 1 and 2, and Algorithm 1 outlines its pseudo-code in Table 2.

Table 2

The pseudo-code of Algorithm 1

Sparse model selection (SINDy-AIC)

For higher accuracy, the models obtained by regression analysis may be further refined in conjunction with statistical learning methods. There are many model selection methods, including the Akaike information criterion (AIC) [19, 20], Bayesian information criterion (BIC) [21], deviance information criterion (DIC) [23] and cross-validation (CV) [22].

In this paper, the AIC is chosen as the statistical score for model comparison from many combinatorially probable models. Originally pioneered by Akaike.H. [19, 20], AIC incorporated the principles of information entropy and Kullback–Leibler (K-L) distance [18] and built upon the notion of maximum likelihood estimation to appraise parameters. Given a candidate Model i, its corresponding AIC value iswhere denotes the conditional probability function for state variable x under the condition of the estimation , and k indicates the number of free parameters. It should be noted that 2 k represents a penalty enforced on the lower bound of the AIC values. In a wide variety of cases, the sampled data are finite datasets and a correction to the AIC value is required bywhere l refers to the dimension of the observed data. In linear regression, the residual sum of squares (RSS) is generally applied to least-squares fitting for the objective function, which is embodied in the form of , where denotes the independent variable for observation, signifies the observed dependent variable, and the candidate model is defined by . Thus, Eq. (10) can be reformulated as and its value will be used in Eq. (11).

Indeed, the AIC scores prefer a rescaled criterion named by AIC with the value of the minimum AIC since there are multiple AIC values for each potential models. In addition, the relative AICc can be linked to the statistical indicator p-value, so that the relative AIC values can be straightforwardly defined as a strength-of-evidence criterion for a contrast with informative models [61]. The strong support value corresponds to the models with , and the weak support is relative to the models with , while the models with are interpreted as no support.

The observations are divided into a training set and a validation set to verify the performance of this approach. The length and the sampling frequency of the time-series traces in each validation set are the same as the homologous training data, except for the initial conditions. For each example in Sect. 4.2, to validate the noise sensitivity, Gaussian noise with mean zero and standard deviation is added to both the training sets and the validation sets. The sparse threshold appearing in the validation sets is optimized by cross-validation.

Identification of parametric PDEs with group sparsity

In this section, the partial differential equation functional identification of nonlinear dynamics (PDE-FIND) [62] is introduced firstly. The mathematical form of the PDEs with constant coefficients is formulated aswhere the unknown nonlinear function N(·) is assumed as the summation of monomial basis functions comprised by variable v and its differential values, and the subscripts indicate the partial equation.

For simplicity, identifying a constant coefficient PDE can be idealized as an overdetermined linear regression problem , its detailed formulation is illustrated byNote that measurement matrix Vt demonstrates data sampled on m time instances and n spatial grids. is a matrix composed with potential models and represents the sparse vector of coefficients.

In contrast to PDE-FIND, the coefficients of either temporally or spatially parameter-dependent PDEs are required to be learned in this work. The time-varying parametric dependent expressions are shown as

The form of Eq. (14) has similarity to Eq. (12) except for the parametric term . If the coefficient estimates are obtained in variable space, is selected to substitute .

Similar to Eq. (4), the guessing candidate function library is constructed as

It encapsulates cross terms constituted by derivatives and parameters, where the form of the series of m equations is shown by Eq. (16).

Unlike the sparse regression problems in SINDy, the constraint that all the share the same sparsity threshold is considered in this article. Therefore, the notion of group sparsity [63] is introduced to group sets of features in together. Then Eq. (13) can be transformed into the block diagonal form

Similar to Eq. (13), the class of equations represented by Eq. (16) can also be regarded as a solitary linear system. And the group sparsity will be involved in the original SINDy algorithm, then Eq. (9) will change towhere denotes the groups, for m time grids and q potential features in library, is defined as . And indicates the block diagonal matrix .

Indeed, the approach evenly divides the measurements into different time steps or spatial locations with respect to the varied coefficients. Additionally, the validation set is constructed with 20% sampled data. And 20 cross-validation experiments are used to get the optimal solution to the sparse threshold. Next the well-selected models will be evaluated with the AIC-enlightened loss function:where n represents the amount of rows in , k denotes the amount of the non-zero elements in the learned PDEs, indicates the threshold benchmark that exerts on each model to prevent over-fitting, and represents the normalized vector of that all terms in it have unit length.

Identify noise

In practice, a wide range of systems of interest can be affected by many adverse factors such as noise, which are usually ineradicable. Indeed, the challenging problems in actual application are how to identify the regime of the noise distribution and how to learn the type of noise apart from the Gaussian noise with mean zero and standard variance .

Crucially, noise can be considered as additional data added to the clean data [64]. Based on this notion, the observations u(t) = [u1(t), u2(t), …, u(t)] 1× can be classified into noiseless data x(t) = [x1(t), x2(t), …, x(t)] 1× and noisy estimation n(t) = [n1(t), n2(t), …, n(t)] 1×.

Therefore, the candidate function library can also consist of two parts: the noiseless measurements X = [x(t1); x(t2); …; x(t)] and the noisy measurements N = [n(t1); n(t2); …; n(t)], which is illustrated as

Throughout the training, the advent of noise has the implication on the correctness of the model identification. Therefore, it is essential to isolate the noise from signals as two segments, the noiseless observations and the noisy observations, whose patterns are constituted with past states and future states, in which step q denotes the degree advanced to current snapshots. Accordingly, the model is illustrated in integration form:where indicates the vector field of the system, and F(·) denotes the nonlinear dynamics.

The two parts are then distinctly solved for sparse coefficients, and the result is compromised with the summation of distinct products. Aiming to diminish the absolute error between the estimated derivative and the actual derivative , named, , which is shown as

The additional constraint would be needed to regularize the Eq. (23) to couple the optimization parameters and . Therefore, the formulation is revised bywhere is used to constrain the numerical error, h is the decay factor, which is set to 0.9, and the details can be described as . Particularly, upper triangle symbols represent the estimation value and the subscripts indicate the time step.

As the validation set is separated for noiseless and noisy estimation, the loss function is naturally composed of two elements: numerical error and simulation error . Additionally, is provided by the simulation of the vector field across the overall trajectory.

Since certain errors are bound to exist due to the extra noise, the same performance criterion advocated by Rudy et al. [65] is applied to judge the performance of the model. Specifically, these errors are the noise estimation error , the vector field error and the simulated trajectory error . Notably, the relative noise estimation error iswhich is defined by the inconsistency between the true noise and the identified noise with the mean error. The vector field error iswhich denotes the deviation from the identified vector field to the true vector field with relative squared error. The simulated trajectory error iswhich represents the discrepancy between the true trajectory and the forward simulation trajectory , especially, denotes the Frobenius norm. The exact process is shown in Figs. 1 and 2.

Results

In this section, we demonstrate the effectiveness of the above methods by several hybrid dynamical systems: the Van der Pol system, the Rössler system, the temporally dependent Burgers’ equation and the spatially dependent Kuramoto–Sivashinsky equation, respectively. Finally, a class of physical systems are provided to investigate the robustness of the proposed approach here.

Experiment I: Sparse identification of nonlinear dynamics

Van der Pol system

The Van der Pol system [66] is defined by

The initial condition is given by and the training set is acquired at interval [− 25, 25] as the snapshot t = 0.001.

Figure 3a shows a comparison of the actual trajectory with the identified trajectory for the Van der Pol system in the form of a phase diagram, where it shows that the learned trajectory is approximated to the true trajectory with the extremely small error. Figure 3b decomposes the two-dimensional phase into two dimensions: x and y. Similarly, the result of the comparison between the actual time-series data and the estimated time-series instances on two dimensions shows that the error between them is extremely small, even almost zero.

Fig. 3

Example of nonlinear dynamics to test the performance of the SINDy algorithm

Rössler system

The Rössler system [67] is also used to verify the performance of SINDy, which is governed by

This system is simulated with initial condition . The training set is sampled with the interval t = 0.001 from t = 0 to t = 25.

Figure 3c represents the actual Rössler system, while Fig. 3d is the identified Rössler system via SINDy. Comparing the two panels, it is known that SINDy accurately reconstructs the Rössler system. Simultaneously, the results of the coefficients indicate that the model error between them is within the range from 0.01% to 0.05%.

Experiment II: Sparse model selection

The Rössler system includes three state variable (m = 3) with a third-order polynomial library (p = 3). The total number of models is in the library, each of them is denoted with possible monomials. One hundred initial conditions are randomly given to generate distinct validation sets with identical dimensions as the training sets, so that 100 cross-validation experiments are enforced on the validation set to adjust the values of the hyper-parameter and to determine the optimal model to be estimated.

As shown in Fig. 4c–e, there are so sizeable candidate models selected from the hierarchically ranked relative AICc that the rescaled relative AICc is used to narrow the range of candidate models. It is found that at the position relevant to the value of 6, the value of the relative AICc falls dramatically. Simultaneously, as shown in Fig. 4e, the 6-term model is equipped with ‘strong support’. Therefore, the model involving 6 nonlinearities is the optimal model.

Fig. 4

Selected sparse models for the Rössler system with SINDy-AIC

Van der Pol oscillator

The Van der Pol system is taken for the test case, whose formulation is governed by Eq. (30), with two state variables (m = 2) and a sixth-order polynomial library (p = 6). On this occasion, there exist potential monomials, and possible models. One hundred cross-validation trials are performed for each system on a validation set of the same dimension as the training set for ranking AICc.

Only the individual trajectory for each state variable (x and y) of the Van der Pol is employed as input, as shown in Fig. 5a. SINDy is able to subselect a series of models from the candidate library, including models with 1, 2 and 4 terms, as shown in Fig. 5b. Notably, there is an abrupt drop in the value of relative AICc at the number of nonlinear terms = 4, and the support values of the selected models fall into the strongly supported range, while other models fall into the unsupported domain according to Fig. 5b. Therefore, the optimal model selected is the one containing four nonlinear terms.

Fig. 5

Selected models for the Van der Pol system with SINDy-AIC

Finally, the value of the prediction error of SINDy-AIC is less than that of SINDy within the range of noise level from to 0.5, as shown in Fig. 5c. Nevertheless, the larger the noise value is, the greater the prediction of SINDy-AIC, with a value close to 1. The robustness of SINDy-AIC can perhaps be the focus of future research.

Experiment III: Identify the coefficients of the partial differential equations

Burgers’ equation

The first example is focused on Burgers’ equation, which has a nonlinear advection term in the form of a sinusoidally oscillatory coefficient , which is defined by

Data are generated on the interval [− 8, 8] via the spectral method and the periodic boundary requirements are constrained by k = 256 snapshots and l = 256 time instances to evenly segment the time scale and to construct the time-varying dependent parametric PDE. Our library of candidate functions is fabricated by employing the cubic order derivatives of v, which can be multiplied by the fourth-order derivatives of v. It is necessary to segregate 20% of the data from each snapshot to be utilized as a validation set. The proper level of sparsity requires the optimal through cross-validation to build the sparse coefficients.

The consequent time series for the identified coefficient and the true coefficient in the case of either noiseless data or noisy data are illustrated in Fig. 6c, d, respectively. The results show that the parametric coefficients are identified elaborately with subtle error.

Fig. 6

Identified coefficients of parametric Burgers’ equation with or without the introduction of noise

Kuramoto–Sivashinsky equation

The Kuramoto–Sivashinsky equation is a nonlinear PDE that simulates it with spatially varying coefficientswhere , , and L = 20. The spatially dependent Kuramoto–Sivashinsky equation is solved mathematically employing a spectral method on a periodic domain [− 20, 20] with k = 512 snapshots and l = 512 time instances. The library consists of the product between the derivatives of v and its cubic order terms up to fourth order.

Figure 7b–d show the differences between the true coefficients and the learned coefficients in the noise-free and noisy measurements, respectively. Traditionally, the goodness-of-fit of the model is judged by the mean square error (MSE) in Eq. (19). And the related mean square errors in two cases are and . Obviously, this method can correctly determine the coefficients for Kuramoto–Sivashinsky in the former, while it does not accurately learn the coefficients with the introduce of noise.

Fig. 7

Identified coefficients of the spatial Kuramoto–Sivashinsky equation and its respective performance indices

Additionally, the AIC can be employed to evaluate the discrepancy for true coefficients and identified coefficients, as shown in Fig. 7e. When the noise magnitude is less than 1, the differences between the identified model and the true model is close to 0.

In particular, the range of the threshold tolerances, from to , is initially employed to select the potential models, where demonstrates the minimum tolerance that has consequences for the sparsity of the , and indicates the minimum tolerance that ensures all coefficients in to be zero. And is illustrated as , where . The notion of and are illustrated in Sect. 3.3, is the threshold, and I indicates the identity matrix.

Figure 7f provides the model evaluation with the loss function from Eq. (19). From 50 values of , ranging from to evenly spaced in steps of , the candidate models will be selected in both intervals of from to and from to . Furthermore, the models corresponding to within an interval are overfit to the measurements, while the models related to are too sparse to predict the coefficients. Then the AIC will be enforced to yield the optimal model for the PDEs.

Comparison results

Regarding the validity of the augmented SINDy with group sparsity, the structure error and the success rate are introduced. The amount of entries in the model that are abnormally increased or removed are denoted by the structure error. As seen in Fig. 8a, the indigo area implies no discrepancy between the discovered model and the actual model, and the violin patterns quantify the allocation of inaccurate entries with cross-validation among 20 models. It is obvious that the augmented SINDy with group sparsity can handle 30,000 more noise data than SINDy-AIC.

Fig. 8

Comparison results of the augmented SINDy and other methods

Figure 8b shows the comparison of the success rate of SINDy, SINDy-AIC and the augmented SINDy with group sparsity in discovering the coefficients of Burgers’ equation, where the success rate is denoted by the average of the results of 15 runs. The augmented SINDy employs approximately 5 times less data than SINDy-AIC, and SINDy-AIC utilizes approximately 2 times less data than SINDy across a distinct percentage of measurements.

Experiment IV: Noise identification

Van der Pol System

Notably, the Van der Pol system is taken as the test instance, and its equation is formulated as Eq. (30).

The system is simulated with two state variables (x and y), and the prediction step q is set to 1. It is essential to compute the Jacobian and to solve Eq. (26). The optimization problem in Eq. (26) utilizes TensorFlow2.7 and the Adam optimizer with a learning rate of 0.001 for calculations. To promote the sparsity of the identified model, the sequentially thresholded is applied to Eq. (26) for times.

The initial conjecture of the next iteration was set as the preceding iteration’s optimization outcome in each iteration. Additional parameters for each example, such as the noise level, the initial conditions, the number of iterations and the sparse threshold , are shown in Table 5.

Table 5

A variety of parameters for different test instances

Models	Initial condition	q	λ	Nloop	Library order	Max Adam optimizer
Duffing	[− 2, 2]	1	0.05	5	3	5000
Van der Pol	[− 2, 1]	1	0.15	8	3	5000
Lotka-Volterra	[1, 2]	1	0.18	10	3	5000
Lorenz	[5, 5, 25]	1	0.20	8	2	15,000
Rössler	[3, 5, 0]	1	0.08	8	2	15,000
SprottB	[1, 1, 1]	1	0.20	8	2	15,000

Figure 9 shows the results after applying the proposed method to the dataset with the introduction of both 15% uniform and 15% normal measurement noise. It illustrates that the noise estimation is approximate to the true measurement noise with the form of the single simulation trajectory for the time series, histogram and probability density respectively, as shown in Fig. 9a, b, c, e, respectively. The error metrics for each noise distribution shown in Table 3.

Fig. 9

Identification for two kinds of noise distribution and system reconstruction for van der pol system

Table 3

Error metrics for the Van der Pol system

Model	Type of noise distribution	Noise percentage (%)	EN	Ef	EP
Van der Pol	Normal	15	0.00053	0.00024	6.14E−05
	Uniform	15	1.04E−05	8.49E−06	5.40E−06

The Rössler system is taken as the final example, whose form is given by Eq. (31).

Data are generated by three input variables (, and ), and levels of 30% artificial noise for uniform or normal distribution are added.

Similar to the first case in Sect. 4.4.1, a learning rate of 0.001 of the Adam optimizer and TensorFlow2.7 are executed to learn the models. Auxiliary parameters, including the magnitude of noise, prediction step q, initial conditions, the amount of iterations and sparsity threshold , are listed in Table 5.

A total of 2500 time steps from t = 0 to t = 25 are generated via high-fidelity simulation and then these data are corrupted with 30% virtual noise in a normal or uniform distribution. The error metrics for the property of the algorithm are outlined in Table 4, and it is shown that the values are approximate to zero.

Table 4

Error metrics for the Rössler system

Model	Type of noise distribution	Noise percentage (%)	EN	Ef	EP
Rössler	Normal	30	0.02272	0.00046	0.01604
	Uniform	30	0.00025	1.63E−06	3.52E−05

Figure 10 displays the bias in the approximation of the measurement noise for a single trajectory with the time series, the simulation evolving into the attractor and the probability density corrupted by 30% artificial noise. In the left column, the system is contaminated by 30% normal measurement noise. Another column shows the system with 30% uniform noise against the true system.

Fig. 10

The results of identification for both noise distribution with different forms and simulation reconstructions for Rossler system

The learned results, either the simulation trajectory or the probability density of various noise distributions, approach the actual values with minimal error. Additionally, it is evident that the differences for uniform observation noise are smaller than the measurements with 30% normal noise, in such an occasion where the error for 6 iterations declines gradually to approximately 0.

Finally, a range of canonical systems are provided to verify the robustness to significant measurement noise for the algorithm, and their relevant expressions are formulated as Eqs. (30), (31), (34), (35), (36) and (37).

The results are boxed in Fig. 11, which shows that the systems are reconstructed successfully with the levels of increasing noise, ranging from 5 to 30%. More broadly, the parameters and the performance indices for these systems are detailed in Tables 5, 6 and 7.

Fig. 11

Corrupted systems and identified systems of a set of levels of artificial noise

Table 6

Error metrics for 5–15% noise of different test cases

Models	5%			10%			15%
	EN	Ef	EP	EN	Ef	EP	EN	Ef	EP
Duffing	8.48E−06	4.43E−05	0.0001	3.59E−05	0.00018	0.00042	8.33E−05	0.0004	0.0009
Van der Pol	5.41E−05	1.77E−05	5.07E−06	0.00022	8.76E−05	2.32E−05	0.00053	0.00024	6.14E−05
Lotka-Volterra	0.00073	3.01E−05	4.48E−05	0.00303	0.00016	0.00011	0.03211	0.00018	0.50404
Rössler	0.0005	6.28E−06	0.00045	0.0021	3.10E−05	0.00181	0.00488	7.84E−05	0.00401
Lorenz	0.00484	0.00011	0.48734	0.0193	0.00042	0.46719	0.05444	0.00117	0.44858
SprottB	6.53E−05	1.07E−05	1.53695	0.00026	4.44E−05	1.02298	0.0006	0.00011	1.07388

Table 7

Error metrics for 20–30% noise of different test cases

Models	20%			25%			30%
	EN	Ef	EP	EN	Ef	EP	EN	Ef	EP
Duffing	0.00015	0.00067	0.00149	0.00024	0.00101	0.00216	0.00035	0.00138	0.00287
Van der Pol	0.00099	0.00046	0.00012	0.00163	0.00078	0.00019	0.00245	0.00119	0.00029
Lotka-Volterra	0.00765	0.00037	0.00589	0.0183	0.00108	0.00074	0.02883	0.00164	0.00168
Rössler	0.00884	0.00017	0.00722	0.0147	0.00029	0.00029	0.02272	0.00046	0.01604
Lorenz	0.1267	0.00227	0.44827	1.59086	0.01819	0.40318	1.62864	0.01896	0.50842
SprottB	0.00109	0.0002	1.39116	0.00173	0.00032	1.41288	0.00255	0.00048	1.11184

As the prediction step q, the number of iterations for denoising the data and sparse threshold affect the results of this algorithm, ergodic trials are utilized to undermine the derivative error and to identify the optimum for the abovementioned parameters, with the exception of the overwhelming computational disadvantage.

Discussion and conclusion

The incorporation of mathematical methods is recommended here as follows: (i) identification of nonlinear dynamics via SINDy, (ii) optimal model selection through the AIC, (iii) integration between sparse regression and group sparsity for coefficients of PDEs and (iv) automated identification of noise via data splitting.

From the point of view of experimental results, the proposed method accurately identifies the coefficients for either the ODEs or the PDEs. More broadly, the following limitations may be attracted.

Structural analysis. The clean measurements, the proper candidate function library and the appropriate regularized sparsity threshold are regarded as three preliminary elements for the structure of the SINDy model. Additionally, the pure measurements play the most important role in the framework; in case noise is introduced, the items in the identified models are bound to be incorrect regardless of the choice of the candidate function library and the sparsity threshold. Furthermore, the sparsity threshold determines the coefficient values of the identified models, and the type of the identified coefficients is defined by potential items in the candidate function library. Due to the actual consideration, the noise is impossible to neglect. Crucially, the indispensable process in sparse regression is parameter tuning.

Model analysis. As the structure of the majority of the models in engineering or biological fields is unknown or complex, the class of the items in the candidate function library is inconclusive. Therefore, the solution to the identified model is obscure. More broadly, the threshold selection will undoubtedly be influenced by the polynomial order in the potential function library. As the magnitude of the polynomial order increases, a regularization threshold may be employed to scale down the dimensionality of the function library.

Robustness analysis. In fact, there are many kinds of probability distributions of noise, but uniform and normal distributions are showcased in this paper. Moreover, the identified model will not be reconstructed with the introduction of the noise due to the asymmetric structure of the probability distribution, such as the gamma or Rayleigh distribution. The method may be extended to integrate with machine learning models to detect the regime of various probability distributions of noise, particularly those with a nonzero mean.

In conclusion, the proposed method has consequences for data-driven model identification and signal-to-noise separation. The flexibility of the architecture of this approach yields a structure for discovering the underlying dynamics of intricate systems with interpretability. Moreover, it can perhaps be engaged with brain science to exploit epilepsy or Alzheimer’s disease.

Table: The summation of the variable

Symbol	Meaning	Symbol	Meaning
State variable		h	The decay factor in Sect. 3.4
X˙,X	The training set in Eq. (2) and Eq. (3)	ω_l	The constrained condition in Eq. (24)
X^,X^˙	The estimated measurements in Eq. (23)	μ(t)	The parametric dependencies in Eq. (14)
Y	The validation set in Sect. 2	N_a	The number of possible monomials in Sect. 4.2
U	The actually observed data in Eq. (20)	N_c	The total number of potential models in Sect. 4.2
N^	The estimated noise in Eq. (26)	p	The pth order polynomial library in Sect. 4.2
X^Q_i	The ith order of polynomial terms in Eq. (5)	m	The dimension of state variable in Sect. 4.2
Θ	The candidate function collection in Eq. (4)	λ	The threshold in Sect. 4.4
V t	The measurements for partial differential equations (PDEs) in Eq. (13)	Performance criterion
Θ~	The block diagonal matrix Θ in Eq. (17)	e_d	The absolute error of the estimated derivative and system’s vector field in Eq. (23)
Ξ	The vector of the sparse coefficients in Eq. (6)	e_s	The simulation error in Eq. (25)
ξ_i	The ith element in the vector of the sparse coefficients in Eq. (7)	E_N	The noise identification error in Eq. (27)
θ_l	The lth potential features in the candidate function library in Sect. 2	E_f	The vector field error in Eq. (28)
x_i	The state variables of distinct dimension in Eq. (30)	E_P	The prediction error in Eq. (29)
v_xx…x⏟_n	The nth order derivative in Eq. (12)	Function
t_i	The time points in Sect. 3.1	f(·)	The function of nonlinear dynamics in Eq. (1)
Parameter		F(·)	The flow map of the dynamic systems in Eq. (21)
Γ	A collection of groups in Eq. (18)	RSS	The residual sum of squares in Sect. 3.2
μ^	The estimated parameter in Eq. (10)	L(·)	The loss function in Eq. (26)
k	The degree free parameters in Eq. (10)	Δ_i	The relative AICc in Sect. 4.2
l	The dimension of measurements in Eq. (11)	AIC_min	The minimum AIC value in Sect. 3.2
ε	The noise magnitude in Sect. 3.2	R(·)	The regularization function in Eq. (9)
q	The backward or forward time steps in Eq. (22)	N(·)	The nonlinear function for the dynamics in Eq. (12)