Peiyun Yin, Xiaoge Guo
College of Informatics, Huazhong Agricultural University, Wuhan, China.
DOI: 10.4236/jcc.2026.145003   PDF    HTML   XML   29 Downloads   139 Views   Citations

Abstract

Sparse identification of nonlinear dynamical systems is an important project, directly addressing the physics community’s long-standing goal of data-driven discovery. Although many effective methods have been developed to enhance the ability for data-driven discovery of governing equations, critical challenges, including derivative computation and optimization methods, demanding more robust theoretical foundations and innovative computational strategies. This paper proposes a novel sparse optimization strategy for learning nonlinear dynamics from noisy data. The strategy is based on the Frank-Wolfe algorithm (also known as the conditional gradient method), an iterative first-order optimization technique for constrained convex optimization. Unlike traditional methods that rely on projection, the algorithm minimizes a linear approximation of the objective function over the feasible set to determine the search direction, naturally yielding sparse solutions, which is consistent with prior work in sparse identification of nonlinear dynamics (SINDy). However, our approach differs fundamentally: while SINDy enforces sparsity through user-defined thresholding, the Frank-Wolfe algorithm inherently produces sparse solutions by minimizing a convex objective over a compact set. This distinction endows our method with greater scalability to large-scale dynamics and noisy data. To further mitigate instabilities in numerical differentiation, particularly under noisy conditions, we incorporate an integral form to replace derivative estimation. The proposed framework thus offers a robust, scalable, and interpretable alternative for data-driven modeling of complex systems. Numerical experiments demonstrate the effectiveness and superiority of our proposed method, which combines subsampling and data normalization, in identifying nonlinear dynamical systems.

Keywords

Nonlinear Dynamic System, Sparse Identification, Subsampling, Data Normalization, The Frank-Wolfe Algorithm

Share and Cite:

1. Introduction

With the exponential growth of data and computing power, data-driven discovery of governing equations has emerged as a transformative technique in science and engineering, enabling the discovery of dynamical system behaviors from observational data [1]-[7]. These rediscovered equations, such as the Navier-Stokes equations [8] [9] and the Lorentz equations [5], provide fundamental mathematical tools for analyzing and predicting system evolution. Recently, this technique has been extended to complex systems, including multiscale systems [10], hybrid systems [11], and complex networks [12]. Despite these advancements, critical challenges remain, including the generation of candidate functions, the computation of derivatives and the optimization method, which call for more robust theoretical foundations and innovative computational strategies.

Currently, many research efforts focused on the deep learning framework for the discovery of governing equations [3] [13] [14]. A notable example is the Physics-Informed Neural Networks (PINNs) proposed in [3]. PINNs learn the solutions to nonlinear partial differential equations (PDEs) by leveraging deep neural networks (DNNs) to capture hidden physical laws from data. Since both the solution function and the system dynamics are approximated by deep neural networks, the selection of basis functions can be avoided, and the derivatives can be evaluated at machine precision using symbolic or automatic differentiation. Authors of [15] generalized PINN to the fractional PINNs (fPINNs) where the residual in the loss function is introduced by using both automatic differentiation for the integer-order operators and numerical discretization for the fractional operators. Building on the PINN, some other versions also are developed, including physics-informed DeepONets [14], Kolmogorov-Arnold networks [16] and Fourier neural operator [17]. However, as with other neural networks, the trained PINN can act as a “black box”. It is often difficult to interpret the model’s internal workings or extract clear physical insights.

Instead of using deep neural network to represent the system dynamics, pioneering studies of data-driven discovery, proposed by [1] and [2], introduced genetic programming-based symbolic regression as a paradigm-shifting method. This evolutionary algorithm searches the space of mathematical expressions to identify nonlinear differential equations that optimally balance model complexity (quantified by term count) and predictive accuracy. By doing so, it fundamentally circumvents the rigid equation form constraints imposed by conventional regression techniques [18], which require predefined symbolic structures. Symbolic regression further addresses a long-standing challenge in physics: the simultaneous optimization of both structural and parametric representations. This capability directly aligns with the community’s goal of data-driven discovery, as demonstrated in studies on predictive modeling [19] and robust parameter estimation [20]. However, symbolic regression methods suffer from a key limitation, high computational costs and overfitting, especially for high-dimensional nonlinear dynamics. For this, some attempts have been made to improve the weakness [21]-[24]. In [24], authors combine dimensional analysis with the AI Feynman symbolic regression algorithm to show that dimensional analysis significantly improves the accuracy of the recovered equation. Therein, PINNs also are used to uncover hidden terms in differential equations.

In contrast, SINDy (Sparse Identification of Nonlinear Dynamics) is another inspiring data-driven method that leverages sparse regression to extract governing equations from observational data [5]. SINDy has demonstrated remarkable versatility in applications ranging from fluid dynamics [25] [26] to biological systems [27] [28]. The core of SINDy lies in its ability to achieve model interpretability through sparse constraints. To do this, this method constructs a large library of candidate functions (e.g., polynomials, trigonometric functions) and combines them with sequential threshold optimization techniques to automatically identify the nonlinear terms that significantly contribute to system dynamics. Unlike conventional least-squares regression [18] and the aforementioned symbolic regression that may overfit noisy observations, SINDy’s sparse regularization intrinsically suppresses irrelevant terms, yielding parsimonious models that generalize better.

Nevertheless, since the sparse regularization works by empirical pruning of the user-defined thresholds [29], this method encounters challenges in scenarios where the true underlying dynamics are unknown or highly complex. Additionally, the method’s performance can be sensitive to the choice of basis functions and the estimation of derivatives, which may introduce biases or errors in the identification process. Therefore, innovative algorithms are needed to enhance the performance of the SINDy method [30]-[33].

In this paper, we proposed a different sparse optimization strategy to learn the nonlinear dynamics from noisy data. The optimization strategy uses the classical Frank-Wolfe algorithm, also known as the conditional gradient method [34] [35]. It is an iterative first-order optimization algorithm used for constrained convex optimization. Instead of performing a projection, the algorithm minimizes a linear approximation of the objective function over the feasible set to determine the next search direction. Notably, this method produces sparse solutions which is consistent with the developed method [5] [7]. The difference is that the Frank-Wolfe algorithm is a mathematical optimization method, which aims at minimizing a convex, differentiable objective function over a compact convex set and then naturally obtains a sparse solution, while the SINDy method ensures the sparsity via the user-customized thresholding. Thus, our proposed method has the larger scalability to large-scale dynamics and noisy data. Additionally, to ease the instability of numerical differentiation especially in the noise case, the integral form, initially developed by [36], also is used to replace the estimation of derivatives. Thus, the proposed framework offers a robust, scalable, and interpretable alternative for data-driven modeling of complex systems.

This paper is organized as follows: in the second section, the problem studied in this paper is described as a sparse leaning problem; in the third section, the sub-sampling and data normalization are, and the used Frank-Wolfe algorithm to solve the sparse regression problem is introduced. In the fourth section, several numerical experiments are presented to verify the performance of the proposed method. The fifth section summarizes the text.

2. Problem Statement

Consider a nonlinear dynamic system with the following form:

$\begin{array}{r}\hfill \dot{x}=f\left(x\left(t\right)\right),x\left(0\right)=x_0,\end{array}$ (1)

where the vector $x\left(t\right)={\left[x_1\left(t\right),x_2\left(t\right),\cdots ,x_n\left(t\right)\right]}^{\text{T}}\in {ℝ}^n$ is the state of the system at time $t$ , the function $f(\cdot ):{ℝ}^n\to {ℝ}^n$ represents the dynamic constraints that mathematically defines the equations of motion governing the system, $x_0$ represents the initial state of the system, and $\dot{x}$ represents the derivative of the state variable of the system with respect to time. In general, when the function $f$ is known and satisfies Lipschitz continuity, the system (1) has a unique solution [36]. However, it is a challenge for the system modeling community to specify the $f$ in terms of the extensive physical knowledge. This work focuses on discovering the dynamic function $f$ from the measurement data containing noise.

For this, we first assume that the data for the system (1) is available, which is denoted by ${\left\{x\left(t_k\right)\right\}}_{k=1}^K$ at discrete time points $t_1,t_2,\cdots ,t_K$ and $x\left(t_0\right)=x_0$ . To discover the dynamics function $f$ , a widely adopted method is to represent the function $f\left(x\right)$ as a linear combination of a set of trial functions. More specifically, the $i$ -th ($i=1,2,\cdots ,n$ ) component $f_i\left(x\right)$ of the function $f\left(x\right)$ can be expressed as a linear combination of $N$ trial functions ${\theta }_j\left(x\right)$ , where the trial function can be the polynomial, trigonometric function, to name a few. Based on this, $f_i\left(x\right)$ can be expressed as follows:

$f_i\left(x\right)={\displaystyle {\sum }_{j=1}^N{c}_{i,j}{\theta }_j}\left(x\right),$ (2)

and Equation (1) combined with Equation (2) can be rewritten into

${\dot{x}}_i\left(t\right)={\displaystyle {\sum }_{j=1}^N{c}_{i,j}{\theta }_j}\left(x\left(t\right)\right),$ (3)

where $c_{i,j}$ is the to-be-identified coefficient. Since Equation (3) holds at any time, one can further write Equation (3) as a linear system in the matrix form:

$\begin{array}{c}{\dot{x}}_i=\Theta \left(x\right)c_i,\end{array}$ (4)

where

$\begin{array}{c}{x}_i={\left[x_i\left(t_1\right),x_i\left(t_2\right),\cdots ,x_i\left(t_K\right)\right]}^{\text{T}},\end{array}$

$\begin{array}{c}{\dot{x}}_i={\left[{\dot{x}}_i\left(t_1\right),{\dot{x}}_i\left(t_2\right),\cdots ,{\dot{x}}_i\left(t_K\right)\right]}^{\text{T}},\end{array}$

$\begin{array}{c}\Theta \left(x\right)=\left[\begin{array}{cccc}{\theta }_1\left(x\left(t_1\right)\right)& {\theta }_2\left(x\left(t_1\right)\right)& \cdots & {\theta }_N\left(x\left(t_1\right)\right)\\ {\theta }_1\left(x\left(t_2\right)\right)& {\theta }_2\left(x\left(t_2\right)\right)& \cdots & {\theta }_N\left(x\left(t_2\right)\right)\\ ⋮& ⋮& & ⋮\\ {\theta }_1\left(x\left(t_K\right)\right)& {\theta }_2\left(x\left(t_K\right)\right)& \cdots & {\theta }_N\left(x\left(t_K\right)\right)\end{array}\right],\end{array}$

and

$\begin{array}{c}{c}_i={\left[c_{i,1},c_{i,2},\cdots ,c_{i,N}\right]}^{\text{T}}.\end{array}$

Throughout this work, the candidate functions ${\theta }_j\left(x\right)$ are chosen as polynomial monomials of the state variables. For a system with state dimension $n$ , a library of total degree $d$ includes all monomials of the form $x_1^{p_1},x_2^{p_2},\cdots ,x_n^{p_n}$ with $0\le {\displaystyle {\sum }_{i=1}^n{p}_i}\le d$ . The maximum degree $d$ and the specific set of monomials may vary from one example to another, depending on the expected dominant nonlinearities of the target system. This adaptive choice improves computational efficiency and identification accuracy. For each numerical experiment in Section 4, we explicitly state the selected library in the corresponding subsection. The same library is used for all state variables within that experiment. As a consequence, one can identify the unknown coefficients $c_i$ by solving Equation (4). However, it is difficult to estimate the derivatives from the noisy data due to the sensitivity of numerical computing.

To accurately infer $f$ from data, the integral form can be adopted to avoid direct differentiation of noisy data [36]. Specifically, the integral form of Equation (4) is as follows:

$\begin{array}{c}{x}_i\left(t\right)-x_i\left(0\right)={\displaystyle {\int }_0^t\Theta }\left(x\left(\tau \right)\right)\text{d}\tau c_i.\end{array}$ (5)

For each element ${\theta }_j\left(x\left(t_k\right)\right)$ of the matrix $\Theta \left(x\left(t\right)\right)$ , we can approximate it using piecewise constant quadrature and represent it by

$\begin{array}{c}{\phi }_j\left(x\left(t_k\right)\right)={\displaystyle {\int }_0^{t_k}{\theta }_j\left(x\left(\tau \right)\right)\text{d}\tau }\approx \text{Δ}t{\displaystyle {\sum }_{l=1}^K{\theta }_j}\left(x\left(t_l\right)\right).\end{array}$ (6)

From this, we construct the feature matrix

$\begin{array}{c}\Phi \left(x\right)=\left[{\Phi }_1\left(x\right),{\Phi }_2\left(x\right),\cdots ,{\Phi }_N\left(x\right)\right]\\ =\left[\begin{array}{cccc}{\phi }_1\left(x\left(t_1\right)\right)& {\phi }_2\left(x\left(t_1\right)\right)& \cdots & {\phi }_N\left(x\left(t_1\right)\right)\\ {\phi }_1\left(x\left(t_2\right)\right)& {\phi }_2\left(x\left(t_2\right)\right)& \cdots & {\phi }_N\left(x\left(t_2\right)\right)\\ ⋮& ⋮& & ⋮\\ {\phi }_1\left(x\left(t_K\right)\right)& {\phi }_2\left(x\left(t_K\right)\right)& \cdots & {\phi }_N\left(x\left(t_K\right)\right)\end{array}\right].\end{array}$ (7)

and collect all the displacements of the $i$ -th state of the system from the initial time $t=0$ to the time $t_k\left(k=1,2,\cdots ,K\right)$ into a column vector:

$\begin{array}{c}{X}_i=x_i\left(t\right)-x_i\left(0\right)=\left[\begin{array}{c}{x}_i\left(t_1\right)-x_i\left(0\right)\\ x_i\left(t_2\right)-x_i\left(0\right)\\ ⋮\\ x_i\left(t_K\right)-x_i\left(0\right)\end{array}\right].\end{array}$ (8)

By combining Equation (5)-Equation (8), we get the following discrete linear system:

$\begin{array}{c}{X}_i=\Phi \left(x\right)c_i,\end{array}$ (9)

where $X_i\in {ℝ}^K$ is the displacement vector, $\Phi \in {ℝ}^{K\times N}$ is the feature matrix, and $c_i\in {ℝ}^N$ is the unknown sparse vector. Therefore, the identification problem of nonlinear dynamic systems can be transformed into a sparse regression problem with a solution coefficient vector $c_i$ .

3. Sparse Recognition of Nonlinear Dynamical Systems Based on Frank-Wolfe Algorithm

This section mainly introduces a data-driven method for identifying nonlinear dynamical systems. In this method, the vector $X_i$ and the feature matrix Φ in Equation (9) are subsampled [37], and further are normalized to obtain a new linear system. For the new linear system obtained after subsampling, the Frank-Wolfe algorithm [38] is used to solve the sparse coefficient vector $c_i$ .

3.1. Subsampling and Data Normalization

This subsection introduces the used subsampling and the random subsampling. For system subsampling, we first determine the sampling interval $m_1$ and the total sample $M_1$ , and select position “1” as the sampling starting point to generate the sampling row index. For random subsampling, the total sample $M_2$ is determined, and the row index of unduplicated samples is randomly generated. The method of subsampling is to combine the row index generated by system subsampling and random subsampling and apply it to the vector $X_i$ and the matrix Φ at the same time [8]. We mark the process of subsampling as $\mathcal{S}$ , then the linear system obtained after subsampling is

$\begin{array}{c}\mathcal{S}{X}_i=\mathcal{S}\Phi c_i,\end{array}$ (10)

where $\mathcal{S}{X}_i\in {ℝ}^{M_1+M_2}$ , $\mathcal{S}\Phi \in {ℝ}^{\left(M_1+M_2\right)\times N}$ and the subsampling process is shown in Figure 1.

Figure 1. Take row-identical subsamples of vector $X_i$ and matrix Φ in the linear system (9).

In the process of integrating the linear system (4), there will be a large difference of orders of magnitude, which will lead to poor identification effect. In order to solve this problem, we need to normalize the data to maintain data stability and improve the generalization ability of the model. For the vector $\mathcal{S}{X}_i$ in the linear system (10), the maximum absolute is used to scale to it. Then, we have

$\begin{array}{c}{Y}_i=\frac{\mathcal{S}{X}_i}{\text{max}\left(\left|\mathcal{S}{X}_i\right|\right)}.\end{array}$ (11)

For matrix $\mathcal{S}\Phi$ , applying the $L_2$ norm normalization to $\mathcal{S}\Phi$ column-wise:

$\begin{array}{c}{D}_j=\frac{\mathcal{S}{\Phi }_j}{{\Vert \mathcal{S}{\Phi }_j\Vert }_2},\end{array}$ (12)

then all normalized column vectors $D_j$ are represented as columns of the new matrix $D$ by:

$\begin{array}{c}D=\left[D_1,D_2,\cdots ,D_N\right].\end{array}$ (13)

Therefore, the linear system after subsampling and normalization is as follows:

$\begin{array}{c}{Y}_i=D{c}_i,\end{array}$ (14)

where $Y_i\in {ℝ}^{M_1+M_2}$ , $D\in {ℝ}^{\left(M_1+M_2\right)\times N}$ , $c_i\in {ℝ}^N$ . Using this method to preprocess data can help us to identify the governing equations of nonlinear dynamical systems more effectively in terms of reducing calculation cost, processing unbalanced data and improving data stability.

However, it is also worth noting that in practical applications, due to the fact that the data contains noise, Equation (14) does not strictly hold so that

$\begin{array}{c}{Y}_i\approx D{c}_i.\end{array}$ (15)

This approximation makes it very difficult for us to solve the inverse problem. In order to deal with this challenge, we need to find an optimal ${\widehat{c}}_i$ , which makes the error between $Y_i$ and $D{c}_i$ extremely small, and promotes the sparsity of the coefficient vector $c_i$ . For this, we use the Lasso Regression [39]-[41] to solve this problem:

${\widehat{c}}_i=\underset{c_i\in {ℝ}^N}{\arg \min }\frac{1}{2}{\Vert D{c}_i-Y_i\Vert }_2^2+\lambda {\Vert c_i\Vert }_1,$ (16)

where $\lambda$ is the regularization parameter.

3.2. Frank-Wolfe Algorithm

In this subsection, we propose to use the Frank-Wolfe algorithm [34] [35] to solve the Lasso problem corresponding to Equation (16). Frank-Wolfe algorithm is an iterative method for solving convex optimization problems, especially when the objective function is a differentiable convex function and the constraint is convex set. Therefore, it is especially effective when dealing with convex optimization problems with $L_1$ constraints.

According to [39], Equation (16) can be equitably expressed as:

$\underset{c_i\in {ℝ}^N}{\min }\frac{1}{2}{\Vert D{c}_i-Y_i\Vert }_2^2\text{s}\text{.t}.{\Vert c_i\Vert }_1\le \alpha ,$ (17)

where $\alpha$ is the regularization parameter. Problem (17) is equivalent to (16) via the KKT conditions: the Lagrange multiplier for the ${\ell }_1$ -ball constraint acts as the regularization parameter $\lambda$ . Thus, a larger $\lambda$ (stronger sparsity) corresponds to a smaller $\alpha$ . In practice, we treat $\alpha$ as the main hyperparameter and select it empirically as described in Section IV. Then, we solve the convex optimization problem (17) to find the optimal sparse coefficient vector $c_i$ by using the Frank-Wolfe algorithm.

We define the objective function as

$g\left(c_i\right)=\frac{1}{2}{\Vert D{c}_i-Y_i\Vert }_2^2.$ (18)

Since the Frank-Wolfe algorithm is also called the conditional gradient algorithm, which determines the optimization direction by calculating the gradient, we need to compute the gradient of the objective function $g\left(c_i^{\left(p\right)}\right)$ at the $p$ -th iteration:

$\nabla g\left(c_i^{\left(p\right)}\right)=D^{\text{T}}\left(D{c}_i^{\left(p\right)}-Y_i\right),$ (19)

where $p=1,2,\cdots ,p_{\max }$ ，$p_{\max }$ is the maximum number of iterations. During each iteration, for a given gradient $\nabla g\left(c_i^{\left(p\right)}\right)$ , we need to find a direction $s^{\left(p\right)}$ that minimizes the following linear subproblem:

$s^{\left(p\right)}=\underset{{\Vert s\Vert }_1\le \alpha }{\arg \min }\langle \nabla g\left(c_i^{\left(p\right)}\right),s\rangle .$ (20)

This problem can be simplified by using the dual norm. For the $L_1$ norm constraint problem, the dual norm is the $L_{\infty }$ norm [42], i.e.

${\Vert \nabla g\left(c_i^{\left(p\right)}\right)\Vert }_{\infty }=\underset{j}{\max }\left|{\left(\nabla g\left(c_i^{\left(p\right)}\right)\right)}_j\right|.$ (21)

To find the optimal direction $s^{\left(p\right)}$ , we need to minimize the inner product $\langle \nabla g\left(c_i^{\left(p\right)}\right),s\rangle$ . This is equivalent to moving in the direction of the maximum absolute component of the gradient under the constraints of the $L_1$ norm, where the maximum absolute component is given by [38]

$j^{\left(p\right)}=\underset{j}{\arg \max }\left|{\left(\nabla g\left(c_i^{\left(p\right)}\right)\right)}_j\right|.$ (22)

Since we need the inner product $\langle \nabla g\left(c_i^{\left(p\right)}\right),s\rangle$ to be minimum and ${\Vert s\Vert }_1\le \alpha$ , we need to choose the direction of s to be opposite to the direction of the gradient and take the largest possible step size $\alpha$ in the direction of the gradient

component with the largest absolute value [38], i.e.

$s^{\left(p\right)}=-\alpha sign{\left({\left(\nabla g\left(c_i^{\left(p\right)}\right)\right)}_{j^{\left(p\right)}}\right)}_{e_{j^{\left(p\right)}}},$ (23)

where $e_{j^{\left(p\right)}}$ is a unit vector that takes the value “1” in the direction of the maximum absolute component of the gradient.

Next, we perform an exact linear search between the current solution $c_i^{\left(p\right)}$ and the direction $s^{\left(p\right)}$ to find the optimal step size $\widehat{\gamma }$ to minimize the objective function $g\left(c_i\right)$ :

$\widehat{\gamma }=\underset{\gamma \in \left[0,1\right]}{\arg \min }g\left(c_i^{\left(p\right)}+\gamma d^{\left(p\right)}\right),$ (24)

where $d^{\left(p\right)}=s^{\left(p\right)}-c_i^{\left(p\right)}$ indicates the search direction, and

$g\left(c_i^{\left(p\right)}+\gamma d^{\left(p\right)}\right)=\frac{1}{2}{\Vert \gamma D\left(s^{\left(p\right)}-c_i^{\left(p\right)}\right)+D{c}_i^{\left(p\right)}-Y_i\Vert }_2^2$ (25)

To find the optimal step size $\widehat{\gamma }$ , we take the derivative of Equation (25) and make the derivative zero, then we have

${\widehat{\gamma }}^{\left(p\right)}=\frac{{\left(q^{\left(p\right)}\right)}^{\text{T}}\left(D{c}_i^{\left(p\right)}-Y_i\right)}{{\left(q^{\left(p\right)}\right)}^{\text{T}}{q}^{\left(p\right)}},$ (26)

where $q^{\left(p\right)}=D\left(s^{\left(p\right)}-c_i^{\left(p\right)}\right)$ .

Once the optimal step size ${\widehat{\gamma }}^{\left(p\right)}$ in each iteration is calculated, we can update the current solution $c_i^{\left(p\right)}$ by the optimal step size to get the next solution $c_i^{\left(p+1\right)}$ , i.e.

$c_i^{\left(p+1\right)}=c_i^{\left(p\right)}+{\widehat{\gamma }}^{\left(p\right)}\left(s^{\left(p\right)}-c_i^{\left(p\right)}\right)$ (27)

We define the dual gap $h^{\left(p\right)}$ :

$h^{\left(p\right)}=-{\left(d^{\left(p\right)}\right)}^{\text{T}}\nabla g\left(c_i^{\left(p\right)}\right).$ (28)

If the dual gap $h^{\left(p\right)}$ is less than the preset threshold, the algorithm is considered to have converged and the iteration is stopped. Additionally, the used algorithm also will stop the iterations when the maximum number of iterations is reached.

Finally, the normalized solution vector $c_i$ is restored to the original scale:

$c_i=\frac{c_i^{\left(p\right)}\max \left(\left|\mathcal{S}{X}_i\right|\right)}{{\left(\mathcal{S}{\Phi }_{norm}\right)}^{\text{T}}},$ (29)

where we define

$\mathcal{S}{\Phi }_{norm}=\left[{\Vert \mathcal{S}{\Phi }_1\Vert }_2,{\Vert \mathcal{S}{\Phi }_2\Vert }_2,\cdots ,{\Vert \mathcal{S}{\Phi }_N\Vert }_2\right].$ (30)

for nonlinear dynamical system identification based on the Frank-Wolfe algorithm is summarized in Algorithm 1.

4. Numerical Experiments and Results

In order to verify the effectiveness of our proposed identification method, several numerical experiments are presented in this section.

Algorithm 1. Frank-Wolfe algorithm for identification of nonlinear dynamical systems.

Input:	Displacement vector $X_i\in {ℝ}^K$ , feature matrix $\Phi \in {ℝ}^{K\times n}$ , maximum number of iterations $p_{\max }$ , parameter $\alpha$ and $\tau$ , threshold value $tol$ .
Output:	Vector $c_i$ such that $X_i=\Phi c_i$ .
1:	According to Equation (11) and Equation (12), the displacement Vector $X_i$ and the feature matrix Φ are normalized by using maximum absolute scaling and $L_2$ norm respectively to get $Y_i$ and $D$ , so that $Y_i=D{c}_i$ .
2:	Initialise $c_i=0$ 。
3:	for $p=1,2,\cdots ,p_{\max }$ do:
4:	Computed search direction $d^{\left(p\right)}:d^{\left(p\right)}=s^{\left(p\right)}-c_i{}^{\left(p\right)}$ , $s^{\left(p\right)}=\underset{{\Vert s\Vert }_1\le \alpha }{\arg \min }\langle \nabla g\left(c_i^{\left(p\right)}\right),s\rangle .$
5:	Determine the optimal step size using precise line search: ${\widehat{\gamma }}^{\left(p\right)}=\frac{{\left(q^{\left(p\right)}\right)}^{\text{T}}\left(D{c}_i^{\left(p\right)}-Y_i\right)}{{\left(q^{\left(p\right)}\right)}^{\text{T}}{q}^{\left(p\right)}},$ where $q^{\left(p\right)}=D\left(s^{\left(p\right)}-c_i^{\left(p\right)}\right)$ .
6:	Update solution：$c_i^{\left(p+1\right)}=c_i^{\left(p\right)}+{\widehat{\gamma }}^{\left(p\right)}\left(s^{\left(p\right)}-c_i^{\left(p\right)}\right)$ .
7:	if dual gap $h^{\left(p\right)}=-{\left(d^{\left(p\right)}\right)}^{\text{T}}\nabla g\left(c_i^{\left(p\right)}\right)$ then:
8:	break.
9:	end if
10:	end for
11:	Denormalization：$c_i=\frac{c_i^{\left(p\right)}\max \left(\left|\mathcal{S}{X}_i\right|\right)}{{\left(\mathcal{S}{\Phi }_{norm}\right)}^{\text{T}}}$ .
12:	for $j=1,2,\cdots ,N$ do
13:	if $\frac{c_{i,j}^2}{{\Vert c_i\Vert }_2^2}<\tau$ then:
14:	$c_{i,j}=0.$
15:	end if
16:	end for
17:	return Sparse vector $c_i$ .

4.1. Logistic Equation

Consider the following Logistic equation

$\dot{x}\left(t\right)=1.6x\left(t\right)-{\left(x\left(t\right)\right)}^2.$ (31)

This is a common model of population dynamics. We use the ode45 function in MATLAB to generate the simulation data $x\left(t\right)$ of the Logistic equation under the initial condition $x\left(0\right)=8\times {10}^{-8}$ , the time interval $t\in \left[0,10\right]$ , and the time-step $\Delta t=0.001$ . The generated data trajectory is shown by the solid red line in Figure 2. However, to model the practical case, we add the noise to the data via $\tilde{x}\left(t\right)=x\left(t\right)+\eta \left(t\right)$ , where $\eta \left(t\right)~N\left(0,{\sigma }^2\right)$ , here we choose $\sigma ={10}^{-5}$ .

After obtaining the noisy measurement data $\tilde{x}\left(t\right)$ containing the noise, we can combine the initial condition $x\left(0\right)=8\times {10}^{-8}$ and Equation (8) to get the displacement vector $X$ , and then define the matrix Θ as follows

$\begin{array}{c}\Theta =\left[\begin{array}{ccccc}|& |& |& \cdots & |\\ 1& x& x^2& \cdots & x^9\\ |& |& |& \cdots & |\end{array}\right].\end{array}$

The high maximum degree (9) is intentionally chosen to test the method’s ability to suppress irrelevant higher-order terms under noise. By combining Equation (6) and Equation (7), we can construct the feature matrix Φ.

Next, we use the subsampling method proposed in this paper, which combines systematic subsampling and random subsampling, to simultaneously conduct subsampling on the displacement vector $X$ and the feature matrix Φ. When we perform systematic subsampling, the sampling interval $m_1=100$ and the total sample $M_1=12$ are selected, and the position “1” is selected as the sampling starting point to generate the sampling row index. When conducting random subsampling, the total sample $M_2=30$ is used, and the row index of unduplicated samples was randomly generated. The row index generated during the subsampling process is then applied to the displacement vector $X$ and the feature matrix Φ to obtain $\mathcal{S}X$ and $\mathcal{S}\Phi$ . The subsampling process is shown in Figure 1.

Then, we apply the maximum absolute value scaling in Equation (11) and the $L_2$ norm normalization in Equation (12) to normalize the vector $\mathcal{S}X$ and the matrix $\mathcal{S}\Phi$ , respectively, and then get the linear system: $Y=Dc$ . Finally, we can use the Frank-Wolfe algorithm proposed in this paper to solve the sparse vector $c$ . In the process of applying the Frank-Wolfe algorithm, we set the maximum number of iterations $p_{\text{max}}=40000$ , the regularization parameter $\alpha =3$ , the threshold value $tol={10}^{-8}$ , and the parameter $\tau ={10}^{-4}$ . The value $\alpha =3$ was chosen empirically after testing a small set of candidate values to achieve a desirable trade‑off between sparse coefficient support and data fitting accuracy.

With the noise intensity $\sigma ={10}^{-5}$ , the identification results for the Logistic equation are shown in Table 1. In addition, we also plotted the data trajectories of the learned governing equation, as shown by the dashed black line in Figure 2.

By comparing the true coefficients with the learned coefficients in Table 1, and analyzing the data trajectory of the exact system without noise and the identified system with noise intensity $\sigma ={10}^{-5}$ in Figure 2, it can be demonstrated that the Frank-Wolfe algorithm proposed in this paper, which combines subsampling and data normalization, accurately identifies the correct activity items. This also demonstrates the effectiveness of our approach. In addition, we can see from Fig. that the trajectory of the system identified by the method proposed in this paper almost coincides with the real trajectory, which also shows that the method is robust to the measurement noise. The experimental MSE is 1.4545 × 10⁻⁷ and the RMSE is 3.8137 × 10⁻⁴.

Figure 2. Data trajectories of logistic equation under given conditions. The solid red line represents the true trajectory of the exact system. The dashed black line represents the data trajectory of the sparse identified system under the condition of noise intensity $\sigma ={10}^{-5}$ .

Table 1. Identification results of Logistic Equation (31). The “Basis” represents the used basis functions. The “True” indicates the true coefficients. The “Learning” represents the parameters identified using the proposed method.

Basis	$\text{d}x/\text{d}t$
	True	Learning
1	0	0
$x$	1.6	1.6002
$x^2$	−1	−1.0232
$x^3$	0	0
$⋮$	$⋮$	$⋮$
$x^8$	0	0
$x^9$	0	0

4.2. Lokta-Volterra System

Consider the Lokta-Volterra system

$\{\begin{array}{l}{\dot{x}}_1\left(t\right)=2/3x_1\left(t\right)-4/3x_1\left(t\right)x_2\left(t\right),\\ {\dot{x}}_2\left(t\right)=-x_2\left(t\right)+x_1\left(t\right)x_2\left(t\right).\end{array}$ (32)

The Lotka-Volterra system, also known as the predatorprey model, is a mathematical model used to describe the dynamics of biological populations. Under the initial condition $\left[x_1\left(0\right),x_2\left(0\right)\right]=\left[{10}^{-5},5\times {10}^{-4}\right]$ , the time interval $t\in \left[0,200\right]$ , and the time-step $\Delta t=0.01$ , the simulation data of Lotka-Volterra system is generated. The generated data trajectory is shown by the solid red line in Figure 3(a). Similarly, in the presence of noise, we set the noise intensity $\sigma ={10}^{-5}$ .

We calculate the displacement vector $X_i\left(i=1,2\right)$ according to Equation (8). We construct a polynomial candidate library of total degree 4, including all monomials of the form $x_1^{p_1},x_2^{p_2}$ with $p_1,p_2\ge 0$ and $0\le p_1+p_2\le 4$ . We define the matrix Θ as follows:

$\Theta =\left[\begin{array}{cccccccc}|& |& |& |& |& |& \cdots & |\\ 1& x_1& x_2& x_1^2& x_1{x}_2& x_2^2& \cdots & x_2^4\\ |& |& |& |& |& |& \cdots & |\end{array}\right]$

Although the true system (32) contains only linear and bilinear terms $\left\{x_1,x_2,x_1{x}_2\right\}$ , the higher‑degree monomials are included to test the ability of the Frank‑Wolfe algorithm to suppress irrelevant terms and recover the correct sparse support under noisy measurements. Then, according to Equation (6) and Equation (7), the feature matrix Φ is obtained.

After obtaining the displacement vector $X_i\left(i=1,2\right)$ and the feature matrix Φ, we need to subsample for $X_i\left(i=1,2\right)$ and Φ. In this case, we need to subsample the displacement vectors $X_1$ and $X_2$ and their corresponding feature matrices Φ₁ and Φ₂ respectively. The feature matrices Φ₁ and Φ₂ are the same here, and we have added subscripts for the convenience of the subsampling process.

For the displacement vector $X_1$ and the feature matrix Φ₁, we select the sampling interval $m_{11}=30$ , the total sample $M_{11}=28$ , and the position “1” as the sampling starting point to generate the sampling row index. When conducting random subsampling, we select the total sample $M_{12}=6$ and randomly generate the row index with no duplicate samples. Then the row index generated by the subsampling is applied to both the displacement vector $X_1$ and the feature matrix Φ₁ to obtain $\mathcal{S}{X}_1$ and $\mathcal{S}{\Phi }_1$ . For the displacement vector $X_2$ and the feature matrix Φ₂, we select the sampling interval $m_{21}=150$ , the total sample $M_{21}=110$ , and the position “1” as the sampling starting point to generate the sampling row index. When conducting random subsampling, we select the total sample $M_{22}=8$ and randomly generate the row index with no duplicate samples. Then the row index generated by the subsampling is applied to both the displacement vector $X_2$ and the feature matrix Φ₂ to obtain $\mathcal{S}{X}_2$ and $\mathcal{S}{\Phi }_2$ .

Next, Equation (11) is used to perform maximum absolute value scaling on the vector $\mathcal{S}{X}_i\left(i=1,2\right)$ , and Equation (12) is used to perform $L_2$ norm normalization on the matrix $Y_i=D_i{c}_i\left(i=1,2\right)$ . Thus, the linear system: $Y_i=D_i{c}_i\left(i=1,2\right)$ is obtained. Finally, we use the Frank-Wolfe algorithm proposed in this paper to solve the convex optimization problem with $L_1$ constraints in Equation (16), and then obtain the sparse vectors $c_1$ and $c_2$ .

In the process of applying the Frank-Wolfe algorithm to solve the sparse vector $c_i\left(i=1,2\right)$ , for the sparse vector $c_1$ , we set the maximum number of iterations $p_{1_{\max }}=100000$ and the regularization parameter ${\alpha }_1=500$ . For the sparse vector $c_2$ , we set the maximum number of iterations $p_{2_{\max }}=10200$ and the regularization parameter ${\alpha }_2=500$ . They share a threshold value $tol={10}^{-8}$ and a parameter $\tau ={10}^{-6}$ . Both ${\alpha }_1=500$ and ${\alpha }_2=500$ were selected empirically based on several trials, where larger values of $\alpha$ (weaker sparsity) led to overfitting and smaller values (stronger sparsity) eliminated relevant terms. The value 500 gave the best compromise for both equations under the same subsampling setting. The identification results are shown in Table 2.

Additionally, we also plotted the true trajectories of each component of the Lotka-Volterra system and the identified trajectories of each component of the learned governing equations, as shown in Figure 3(b) and Figure 3(c).

By comparing and analyzing the true and estimated coefficients in Table 2, as well as the true and the identified trajectories in Figure 3, we verify the effectiveness of the proposed method for the identification of the Lokta-Volterra system. However, it can be seen from Figure 3(a) that there are still some errors in the identification of the entire system. In addition, we also calculated the MSE and RMSE between the true parameters and the estimated parameters of the Lotka-Volterra system, as shown in Table 3.

Table 2. Identification results of Lokta-Volterra system (32).

Basis	$\text{d}x/\text{d}t$		$\text{d}y/\text{d}t$
	True	Learning	True	Learning
1	0	0	0	0
$x$	2/3	0.6666	0	0
$y$	0	0	-1	−1.0001
$x^2$	0	0	0	0
$xy$	−4/3	−1.3334	1	1.0009
$y^2$	0	0	0	0
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
$y^4$	0	0	0	0

Figure 3. (a) Phase diagram of Lotka-Volterra system. The solid red line represents the true trajectory of the exact system. The dashed black line represents the data trajectory of the sparse identified system under the condition of noise intensity $\sigma ={10}^{-5}$ . (b) The true trajectory of the state component $x_1\left(t\right)$ is represented by the solid red line. The identified trajectory is represented by dashed black line. (c) The true trajectory of the state component $x_2\left(t\right)$ is represented by the solid blue line. The identified trajectory is represented by dashed black line.

Table 3. MSE and RMSE between the true parameters and the estimated parameters of the Lotka-Volterra system.

	State component x	State component y
MSE	1.7767	8.8032 × 10−1
RMSE	1.3329	9.3826 × 10−1

4.3. Two-Dimensional Damped Oscillator

Consider the two-dimensional damped harmonic oscillator with cubic dynamics

$\{\begin{array}{l}{\dot{x}}_1\left(t\right)=-0.1x_1^3\left(t\right)+2x_2^3\left(t\right),\\ {\dot{x}}_2\left(t\right)=-2x_1^3\left(t\right)-0.1x_2^3\left(t\right).\end{array}$ (33)

The two-dimensional damped harmonic oscillator generates the simulation data under the initial condition $\left[x_1\left(0\right),x_2\left(0\right)\right]=\left[1,1.8\right]$ , the time interval $t\in \left[0,10\right]$ , and the time-step $\Delta t=0.001$ . The data trajectory is shown in Figure 4(a). In the presence of noise, we set the noise intensity $\sigma ={10}^{-5}$ .

After obtaining the noisy measurement data $\tilde{x}$ , we can obtain the displacement vector $X_i\left(i=1,2\right)$ and the matrix Θ, where we construct a polynomial candidate library of total degree 4, including all monomials in $x_1$ and $x_2$ with$0\le p_1+p_2\le 4$ . We define the matrix Θ as follows:

$\Theta =\left[\begin{array}{cccccccc}|& |& |& |& |& |& \cdots & |\\ 1& x_1& x_2& x_1^2& x_1{x}_2& x_2^2& \cdots & x_2^4\\ |& |& |& |& |& |& \cdots & |\end{array}\right]$

Although the true dynamics (33) contain only the cubic terms $x_1^3$ and $x_2^3$ , we include all lower- and higher-degree monomials to evaluate whether the Frank-Wolfe algorithm can correctly identify the sparse support and reject the irrelevant terms in the presence of noise. Then we can obtain the feature matrix Φ by Equation (6) and Equation (7).

When subsampling the displacement vector $X_i\left(i=1,2\right)$ and the feature matrix Φ, we use the same method as in the two examples above. We choose $m_{11}=100$ , $M_{11}=30$ , $M_{12}=20$ and $m_{21}=100$ , $M_{21}=30$ , $M_{22}=20$ . And then we have the vector $\mathcal{S}{X}_i\left(i=1,2\right)$ and the matrix $\mathcal{S}{\Phi }_i\left(i=1,2\right)$ . Then we use the method proposed in this paper to normalize the vector $\mathcal{S}{X}_i\left(i=1,2\right)$ and the matrix $\mathcal{S}{\Phi }_i\left(i=1,2\right)$ respectively, and we can get the linear system: $Y_i=D_i{c}_i\left(i=1,2\right)$ .

Finally, we can use the Frank-Wolfe algorithm to solve the sparse coefficient vector $c_i\left(i=1,2\right)$ . When solving the sparse coefficient vector $c_1$ , we set the maximum number of iterations $p_{1_{\max }}=3000$ and the regularization parameter ${\alpha }_1=3$ . When solving the sparse coefficient vector $c_2$ , we set the maximum number of iterations $p_{2_{\max }}=26000$ and the regularization parameter ${\alpha }_2=10$ . They share a threshold value $tol={10}^{-8}$ and a parameter $\tau =2\times {10}^{-5}$ .

These values were determined by empirical tuning: ${\alpha }_1=3$ and ${\alpha }_2=10$ each produced a sparse coefficient vector whose nonzero entries matched the true support while keeping the residual ${\Vert D{c}_i-Y_i\Vert }_2$ low. Different $\alpha$ values for the two equations reflect their different sensitivities to the ℓ₁ constraint. The identification results are shown in Table 4.

Table 4. Identification results of the two-dimensional damped oscillator (32).

Basis	$\text{d}x/\text{d}t$		$\text{d}y/\text{d}t$
	True	Learning	True	Learning
1	0	0	0	0
$x$	0	0	0	0
$y$	0	0	0	0
$x^2$	0	0	0	0
$xy$	0	0	0	0
$y^2$	0	0	0	0
$x^3$	−0.1	−0.1002	−2	−2.0042
$x^{2}y$	0	0	0	0
$x{y}^2$	0	0	0	0
$y^3$	2	2.0001	−0.1	−0.1000
$x^4$	0	0	0	0
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
$y^4$	0	0	0	0

We also plotted the true trajectories of each component of the two-dimensional damped harmonic oscillator and the identified trajectories of each component of the learned governing equation, as shown in Figure 4(b).

Figure 4. Phase diagram of the two-dimensional damped oscillator and the identified trajectories of state components x and y.

By comparing the true coefficients and learned coefficients in Table 3, and analyzing the true trajectories and identified trajectories in Figure 4(a) and Figure 4(b), it can be effectively demonstrated that the Frank-Wolfe algorithm proposed in this paper, which combines subsampling and data normalization, can correctly identify active items and approximate coefficients. The MSE and RMSE between the true parameters and the estimated parameters are shown in Table 5.

Table 5. MSE and RMSE between the true parameters and the estimated parameters of the two-dimensional damped oscillator.

	State component x	State component y
MSE	1.7299 × 10−5	1.6199 × 10−5
RMSE	4.1592 × 10−3	4.0248 × 10−3

4.4. Lorenz System

Consider the following Lorenz system

$\{\begin{array}{l}{\dot{x}}_1\left(t\right)=-10x_1\left(t\right)+10x_2\left(t\right),\\ {\dot{x}}_2\left(t\right)=28x_1\left(t\right)-x_2\left(t\right)-x_1\left(t\right)x_3\left(t\right),\\ {\dot{x}}_3\left(t\right)=-8/3x_3\left(t\right)+x_1\left(t\right)x_2\left(t\right).\end{array}$ (34)

The Lorenz system is a three-dimensional nonlinear dynamical system that is one of the classic examples of chaos theory, demonstrating complex dynamical behavior. The Lorenz system generates the simulation data under the initial condition $\left[x_1\left(0\right),x_2\left(0\right),x_3\left(0\right)\right]=\left[-152,150,300\right]$ , the time interval $t\in \left[0,10\right]$ , and the time-step $\Delta t=0.001$ . The data trajectory is shown in Figure 5. In the presence of noise, we set the noise intensity $\sigma ={10}^{-5}$ .

From the noisy measurement data $\tilde{x}$ , we can get the displacement vector $X_i\left(i=1,2,3\right)$ and the matrix Θ, where we construct a polynomial candidate library of total degree 3, including all monomials in $\left\{x_1,x_2,x_3\right\}$ with $0\le p_1+p_2+p_3\le 3$ . The matrix Θ is defined as follows:

$\Theta =\left[\begin{array}{ccccccccc}|& |& |& |& |& |& |& \cdots & |\\ 1& x_1& x_2& x_3& x_1^2& x_1{x}_2& x_1{x}_3& \cdots & x_3^3\\ |& |& |& |& |& |& |& \cdots & |\end{array}\right]$

Although the true Lorenz system (34) contains only linear and bilinear terms $\left\{x_1,x_2,x_3,x_1{x}_2,x_1{x}_3\right\}$ , the higher-degree monomials are included to rigorously test the ability of the Frank-Wolfe algorithm to suppress irrelevant terms and recover the exact sparse support in a chaotic system under noise. Then the feature matrix Φ can be obtained by Equation (6) and Equation (7).

Next, the displacement vector $X_i\left(i=1,2,3\right)$ and the feature matrix Φ are subsampled. In the subsampling process, $m_{11}=190$ , $M_{11}=20$ , $M_{12}=5$ are selected for the displacement vector $X_1$ and the feature matrix ${\Phi }_1$ . For the displacement vector $X_2$ and the feature matrix ${\Phi }_2$ , we choose $m_{21}=149$ , $M_{21}=52$ , $M_{22}=8$ . For the displacement vector $X_3$ and the feature matrix ${\Phi }_3$ , we choose $m_{31}=33$ , $M_{31}=6$ , $M_{32}=36$ .Then we normalize the vector $\mathcal{S}{X}_i\left(i=1,2,3\right)$ and the matrix $\mathcal{S}{\Phi }_i\left(i=1,2,3\right)$ according to Equation (11) and Equation (12), and obtain the linear system: $Y=Dc$ . Finally, the sparse coefficient vector $c_i\left(i=1,2,3\right)$ is solved by the Frank-Wolfe algorithm proposed in this paper. When solving the sparse coefficient vector $c_1$ , we set the maximum number of iterations $p_{1_{\max }}=70000$ and the regularization parameter ${\alpha }_1=3.8$ . When solving the sparse coefficient vector $c_1$ , we set the maximum number of iterations $p_{2_{\max }}=49600$ and the regularization parameter ${\alpha }_2=14$ . When solving the sparse coefficient vector $c_3$ , we set the maximum number of iterations $p_{3_{\max }}=68000$ and the regularization parameter ${\alpha }_3=9.59$ . They share a threshold value $tol={10}^{-8}$ and a parameter $\tau ={10}^{-6}$ . Each ${\alpha }_i$ was chosen empirically by trying a set of candidate values and selecting the one that produced the correct sparsity pattern with the smallest fitting error. The distinct values (3.8, 14, 9.59) are necessary because the three equations have different magnitudes and sensitivities to the ${\ell }_1$ constraint. The identification results are shown in Table 6. In addition, we also plotted the true trajectories of each state component of the Lorenz system, as shown in Figure 6. The MSE and RMSE between the true parameters and the estimated parameters are shown in Table 7.

Table 6. Identification results of Lorenz system (34).

Basis	$\text{d}x/\text{d}t$				$\text{d}y/\text{d}t$		$\text{d}z/\text{d}t$
	True	Learning		True		Learning	True	Learning
1	0		0		0	0	0	0
$x$	−10		−10.1408		28	28.0949	0	0
$y$	10		10.0281		−1	−1.0010	0	0
$z$	0		0		0	0	−8/3	−2.6615
$x^2$	0		0		0	0	0	0
$xy$	0		0		0	0	1	0.9965
$xz$	0		0		−1	−1.0050	0	0
$y^2$	0		0		0	0	0	0
$yz$	0		0		0	0	0	0
$z^2$	0		0		0	0	0	0
$x^3$	0		0		0	0	0	0
$⋮$	$⋮$		$⋮$		$⋮$	$⋮$	$⋮$	$⋮$
$z^3$	0		0		0	0	0	0

Figure 5. Identified trajectories of the Lorenz system.

Figure 6. Trajectory of state component x, y, z.

Table 7. MSE and RMSE between the true parameters and the estimated parameters of the Lorenz system.

	State component x	State component y	State component z
MSE	6.4427 × 101	9.1607 × 101	6.3988 × 101
RMSE	8.0267	9.5712	7.9992

By comparing the true coefficients with the learned coefficients in Table 5, and analyzing the true trajectories and identified trajectories in Figure 5 and Figure 6, we can see that the method proposed in this paper is very effective for the identification of complex Lorenz system. Additionally, it can be seen from Figure 5(a) and Figure 5(b) that the considered Lorentz system is sensitive to the system parameters. When the identified parameters contain a very small error, the true trajectory of the Lorenz system also exhibits slight deviations.

5. Conclusion

In this work, we proposed a different sparse optimization strategy to learn the parsimonious governing equation from noisy data. The proposed method combines the Frank-Wolfe algorithm with the integral strategy to enhance the robustness in the noisy condition. Additionally, the subsampling is also used to reduce the computational complexity. Several numerical experiments demonstrate that the proposed method can accurately capture the system equations from noisy data.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Bongard, J. and Lipson, H. (2007) Automated Reverse Engineering of Nonlinear Dynamical Systems. Proceedings of the National Academy of Sciences, 104, 9943-9948.[CrossRef] [PubMed]
[2]	Schmidt, M. and Lipson, H. (2009) Distilling Free-Form Natural Laws from Experimental Data. Science, 324, 81-85.[CrossRef] [PubMed]
[3]	Raissi, M., Perdikaris, P. and Karniadakis, G.E. (2019) Physics-Informed Neural Networks: A Deep Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differential Equations. Journal of Computational Physics, 378, 686-707.[CrossRef]
[4]	Cai, S., Mao, Z., Wang, Z., Yin, M. and Karniadakis, G.E. (2021) Physics-Informed Neural Networks (PINNs) for Fluid Mechanics: A Review. Acta Mechanica Sinica, 37, 1727-1738.[CrossRef]
[5]	Brunton, S.L., Proctor, J.L. and Kutz, J.N. (2016) Discovering Governing Equations from Data by Sparse Identification of Nonlinear Dynamical Systems. Proceedings of the National Academy of Sciences, 113, 3932-3937.[CrossRef] [PubMed]
[6]	Brunton, S.L., Kutz, J.N. and Proctor, J.L. (2017) Data Driven Discovery of Governing Physical Laws. SIAM News, 50.
[7]	Pan, W., Yuan, Y., Gonçalves, J. and Stan, G. (2016) A Sparse Bayesian Approach to the Identification of Nonlinear State-Space Systems. IEEE Transactions on Automatic Control, 61, 182-187.[CrossRef]
[8]	Rudy, S.H., Brunton, S.L., Proctor, J.L. and Kutz, J.N. (2017) Data-Driven Discovery of Partial Differential Equations. Science Advances, 3, e1602614.[CrossRef] [PubMed]
[9]	Schaeffer, H. (2017) Learning Partial Differential Equations via Data Discovery and Sparse Optimization. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473, Article 20160446.[CrossRef] [PubMed]
[10]	Brunton, S.L. and Kutz, J.N. (2019) Methods for Data-Driven Multiscale Model Discovery for Materials. Journal of Physics: Materials, 2, Article 044002.[CrossRef]
[11]	Yuan, Y., Tang, X., Zhou, W., Pan, W., Li, X., Zhang, H., et al. (2019) Data Driven Discovery of Cyber Physical Systems. Nature Communications, 10, Article No. 4894.[CrossRef] [PubMed]
[12]	Gao, T. and Yan, G. (2022) Autonomous Inference of Complex Network Dynamics from Incomplete and Noisy Data. Nature Computational Science, 2, 160-168.[CrossRef] [PubMed]
[13]	Raissi, M., Yazdani, A. and Karniadakis, G.E. (2020) Hidden Fluid Mechanics: Learning Velocity and Pressure Fields from Flow Visualizations. Science, 367, 1026-1030.[CrossRef] [PubMed]
[14]	Chen, Z., Liu, Y. and Sun, H. (2021) Physics-Informed Learning of Governing Equations from Scarce Data. Nature Communications, 12, Article No. 6136.[CrossRef] [PubMed]
[15]	Pang, G., Lu, L. and Karniadakis, G.E. (2019) fPINNs: Fractional Physics-Informed Neural Networks. SIAM Journal on Scientific Computing, 41, A2603-A2626.[CrossRef]
[16]	Liu, Z., Wang, Y., Vaidya, S., Ruehle, F., Halverson, J., Soljacic, M., Hou, T.Y. and Tegmark, M. (2025) KAN: Kolmogorov-Arnold Networks. The Thirteenth International Conference on Learning Representations (ICLR 2025), Singapore, 21-23 April 2025, 1-47. https://openreview.net/forum?id=Ozo7qJ5vZi
[17]	Li, Z., Kovachki, N.B., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A. and Anandkumar, A. (2021) Fourier Neural Operator for Parametric Partial Differential Equations. 9th International Conference on Learning Representations (ICLR 2021), Virtual Event, Austria, 3-7 May 2021, 1-16.
[18]	Ljung, L. (1998) System Identification. In: Procházka, A., Uhlíř, J., Rayner, P.W.J. and Kingsbury, N.G., Eds., Applied and Numerical Harmonic Analysis, Birkhäuser Boston, 163-173.[CrossRef]
[19]	Quade, M., Abel, M., Shafi, K., Niven, R.K. and Noack, B.R. (2016) Prediction of Dynamical Systems by Symbolic Regression. Physical Review E, 94, Article 012214.[CrossRef] [PubMed]
[20]	Reinbold, P.A.K., Kageorge, L.M., Schatz, M.F. and Grigoriev, R.O. (2021) Robust Learning from Noisy, Incomplete, High-Dimensional Experimental Data via Physically Constrained Symbolic Regression. Nature Communications, 12, Article 3219.[CrossRef] [PubMed]
[21]	Vaddireddy, H., Rasheed, A., Staples, A.E. and San, O. (2020) Feature Engineering and Symbolic Regression Methods for Detecting Hidden Physics from Sparse Sensor Observation Data. Physics of Fluids, 32, Article 015113.[CrossRef]
[22]	Petersen, B.K., Larma, M.L., Mundhenk, T.N., Santiago, C.P., Kim, S.K. and Kim, J.T. (2021) Deep Symbolic Regression: Recovering Mathematical Expressions from Data via Risk Seeking Policy Gradients. International Conference on Learning Representations (ICLR 2021), Virtual Event, Austria, 3-7 May 2021. https://www.iclr.cc/virtual/2021/poster/2578
[23]	Cornforth, T. and Lipson, H. (2012) Symbolic Regression of Multiple-Time-Scale Dynamical Systems. Proceedings of the 14th Annual Conference on Genetic and Evolutionary Computation, Philadelphia, 7-11 July 2012, 735-742.[CrossRef]
[24]	Podina, L., Darooneh, D., Grewal, J. and Kohandel, M. (2024) Enhancing Symbolic Regression and Universal Physics Informed Neural Networks with Dimensional Analysis. arXiv:2411.15919. https://arxiv.org/abs/2411.15919
[25]	Fukami, K., Murata, T., Zhang, K. and Fukagata, K. (2021) Sparse Identification of Nonlinear Dynamics with Low-Dimensionalized Flow Representations. Journal of Fluid Mechanics, 926, A10.[CrossRef]
[26]	Zhang, M., Zeng, F., Yousif, M.Z., Zhu, Z., Lee, J.S. and Lim, H. (2025) Convolutional Autoencoder-Augmented Sparse Identification of Noisy Nonlinear Fluid Dynamics from Unlabelled Data. Physics of Fluids, 37, Article 095177.[CrossRef]
[27]	Mangan, N.M., Brunton, S.L., Proctor, J.L. and Kutz, J.N. (2016) Inferring Biological Networks by Sparse Identification of Nonlinear Dynamics. IEEE Transactions on Molecular, Biological and Multi-Scale Communications, 2, 52-63.[CrossRef]
[28]	Wu, X., McDermott, M. and MacLean, A.L. (2025) Data-Driven Model Discovery and Model Selection for Noisy Biological Systems. PLOS Computational Biology, 21, e1012762.[CrossRef] [PubMed]
[29]	Zhang, L. and Schaeffer, H. (2019) On the Convergence of the Sindy Algorithm. Multiscale Modeling & Simulation, 17, 948-972.[CrossRef]
[30]	Fasel, U., Kutz, J.N., Brunton, B.W. and Brunton, S.L. (2022) Ensemble-SINDy: Robust Sparse Model Discovery in the Low-Data, High-Noise Limit, with Active Learning and Control. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 478, Article 20210904.[CrossRef] [PubMed]
[31]	Kaheman, K., Kutz, J.N. and Brunton, S.L. (2020) Sindy-PI: A Robust Algorithm for Parallel Implicit Sparse Identification of Nonlinear Dynamics. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 476, Article 20200279.[CrossRef] [PubMed]
[32]	Tran, G. and Ward, R. (2017) Exact Recovery of Chaotic Systems from Highly Corrupted Data. Multiscale Modeling & Simulation, 15, 1108-1129.[CrossRef]
[33]	Mangan, N.M., Kutz, J.N., Brunton, S.L. and Proctor, J.L. (2017) Model Selection for Dynamical Systems via Sparse Regression and Information Criteria. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473, Article 20170009.[CrossRef] [PubMed]
[34]	Levitin, E.S. and Polyak, B.T. (1966) Constrained Minimization Methods. USSR Computational Mathematics and Mathematical Physics, 6, 1-50.[CrossRef]
[35]	Frank, M. and Wolfe, P. (1956) An Algorithm for Quadratic Programming. Naval Research Logistics Quarterly, 3, 95-110.[CrossRef]
[36]	Schaeffer, H. and McCalla, S.G. (2017) Sparse Model Selection via Integral Terms. Physical Review E, 96, Article 023302.[CrossRef] [PubMed]
[37]	Ma, C., Huang, C., Cheng, C. and Li, X. (2024) Deterministic-Like Data-Driven Discovery of Stochastic Differential Equations via the Feynman-Kac Formalism. The European Physical Journal Special Topics, 234, 721-741.[CrossRef]
[38]	Braun, G., Carderera, A., Combettes, C.W., Hassani, H., Karbasi, A., Mokhtari, A., and Pokutta, S. (2022) Conditional Gradient Methods. arXiv:2211.14103
[39]	Hastie, T., Tibshirani, R. and Friedman, J. (2009) The Elements of Statistical Learning: Data Mining, Inference, and Prediction. 2nd Edition, Springer.
[40]	James, G., Witten, D., Hastie, T. and Tibshirani, R. (2013) An Introduction to Statistical Learning: With Applications in R (Springer Texts in Statistics, Vol. 103). Springer.
[41]	Tibshirani, R. (1996) Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58, 267-288.[CrossRef]
[42]	Balasubramanian, K. and Ghadimi, S. (2018) Zeroth Order (Non) Convex Stochastic Optimization via Conditional Gradient and Gradient Updates. Advances in Neural Information Processing Systems, 31, 3459-3468.