A Robust SINDy Approach by Combining Neural Networks and an Integral
  Form
The discovery of governing equations from data has been an active field of
research for decades. One widely used methodology for this purpose is sparse
regression for nonlinear dynamics, known as SINDy. Despite several attempts,
noisy and scarce data still pose a severe challenge to the success of the SINDy
approach. In this work, we discuss a robust method to discover nonlinear
governing equations from noisy and scarce data. To do this, we make use of
neural networks to learn an implicit representation based on measurement data
so that not only it produces the output in the vicinity of the measurements but
also the time-evolution of output can be described by a dynamical system.
Additionally, we learn such a dynamic system in the spirit of the SINDy
framework. Leveraging the implicit representation using neural networks, we
obtain the derivative information -- required for SINDy -- using an automatic
differentiation tool. To enhance the robustness of our methodology, we further
incorporate an integral condition on the output of the implicit networks.
Furthermore, we extend our methodology to handle data collected from multiple
initial conditions. We demonstrate the efficiency of the proposed methodology
to discover governing equations under noisy and scarce data regimes by means of
several examples and compare its performance with existing methods.
Ali Forootani, Pawan Goyal, Peter Benner
# A Robust SINDy Approach by Combining Neural Networks and an Integral Form

###### Abstract

The discovery of governing equations from data has been an active field of research for decades. One widely used methodology for this purpose is sparse regression for nonlinear dynamics, known as SINDy. Despite several attempts, noisy and scarce data still pose a severe challenge to the success of the SINDy approach. In this work, we discuss a robust method to discover nonlinear governing equations from noisy and scarce data. To do this, we make use of neural networks to learn an implicit representation based on measurement data so that not only it produces the output in the vicinity of the measurements but also the time-evolution of output can be described by a dynamical system. Additionally, we learn such a dynamic system in the spirit of the SINDy framework. Leveraging the implicit representation using neural networks, we obtain the derivative information--required for SINDy--using an automatic differentiation tool. To enhance the robustness of our methodology, we further incorporate an integral condition on the output of the implicit networks. Furthermore, we extend our methodology to handle data collected from multiple initial conditions. We demonstrate the efficiency of the proposed methodology to discover governing equations under noisy and scarce data regimes by means of several examples and compare its performance with existing methods.

Keywords:Spare regression, discovering governing equations, neural networks, nonlinear system identification, Runge-Kutta scheme. +
Footnote †: journal: Journal of Computational and Graphical Statistics

## 1 Introduction

System identification is a crucial aspect of understanding and modeling the dynamics of various physical, chemical, and biological systems. Over the years, various powerful and efficient system identification techniques have been developed, and these methods have been applied in a wide range of applications, see, e.g., [1, 2, 3]. Traditionally, system identification techniques rely on prior model hypotheses. With a linear model hypothesis, several methodologies have been proposed; see, e.g., [1, 2].
However, for nonlinear system identification, defining a prior is challenging, and it is often done with the help of practitioners. Despite several earlier works [4, 5, 6], nonlinear system identification is still an active and exciting research field. Towards automatic nonlinear system identification, generic algorithms and symbolic regression have shown their effectiveness and promises in discovering governing nonlinear equations using measurement [7, 8]. However, their computational expenses remain undesirable.

Instead of building suitable functions in the spirit of symbolic regression, there has been a focus on sparsity-promoting approaches for nonlinear system identification [9, 10, 11]. They rely on the assumption that nonlinear dynamics can be defined by a few nonlinear basis functions from a dictionary with a large collection of nonlinear basis functions. Such a technique enables the discovery of interpretable, parsimonious, and generalizable models that balance precision and performance. It is nowadays widely referred to as SINDy[11]. SINDy has been employed for a handful number of challenging model discovery problems such as fluid dynamics [12], plasma dynamics [13], turbulence closures [14], mesoscale ocean closures [15], nonlinear optics [16], computational chemistry [17], and numerical integration [18]. Moreover, the results of the SINDy have been extended widely to many applications, such as nonlinear model predictive control [19], rational functions [20, 21], enforcing known conservation laws and symmetries [12], promoting stability [22], generalizations for stochastic dynamics [23], from Bayesian perspective [24].

Often, SINDy approaches require a reliable estimate of the derivative information, making them very challenging for noisy and scare data regimes. Blending numerical methods [21, 25, 26, 27, 28] and weak formulations of differential equations [28] avoid these requirements, but their performance still deteriorates for low signal-to-noise measurements. In addition, the method in [28] relies on the choice of basis functions that allow to write differential equations in a weak formulation. The work in [29] utilizes the concepts of an ensemble to improve the predictions, but it still requires reliable estimates of derivatives to some extent. To discover governing equations from noisy data, the authors in [30] proposed a scheme that aims to decompose the noisy signals into clean signals and the noise using a Runge-Kutta-based integration method. However, the scheme explicitly estimates the noise, making it harder to scale, and requires all the dependent variables to be available at the same time grid.

Recently, applications of deep neural networks (DNN) have received attention in sparse regression model discovery methods. For instance, in [31], a deep learning-based discovery algorithm has been employed to identify underlying (partial) differential equations. However, therein, only a single trajectory is considered to recover governing equations, but in many complex processes, we might require data for different parameters and initial conditions to explore rich dynamics, thus, the reliable discovery of governing equations. Furthermore, the work [31] discovers governing equations based on estimating derivative information using automatic differential tools. However, we know that differential equations can also be written in the integral form, whereby the numerical approaches can employed as well, see, e.g., [21, 32].

In this paper, we discuss an approach, namely, iNeural-SINDy, for the data-driven discovery of nonlinear dynamical systems using noisy data from the tens of SINDy. For this, we make use of DNN to learn an implicit representation based on the given data set so that the network outputs denoised data, which is later utilized for the sparse regression to discover governing equations. To solve the sparse regression problem, we make use of not only automatic differential tools but also integral forms for differential equations. As a result, we observe a robust discovery of governing equations. We note that such a concept has recently been used in the context of neural ODEs in [33] to learn black-box dynamics using noisy and scarce data. We further discuss how to incorporate the data coming from multiple initial conditions.

The rest of this paper is organized as follows. Section 2 briefly recalls the SINDy approach [11]. In Section 3, we propose a novel methodology for sparse regression to learn underlying governing equations by making use of DNN, automatic differential tools, and numerical methods. Furthermore, in Section 4, we discuss its extension for multiple initial conditions and different parameters. In Section 5, we demonstrate the proposed framework by means of various synthetic noisy measurements and present a comparison with the current state-of-the-art approaches. Finally, Section 6 concludes the paper with a brief summary and future research avenues.
## 2 A Brief Recap of SINDy

The SINDy algorithm is a nonlinear system identification approach which is based on the hypothesis that governing equations of a nonlinear system can be given by selecting a few suitable basis functions, see, e.g., [11]. Precisely, it aims at identifying a few basis functions from a dictionary, containing a large number of candidate basis functions. In this regard, sparsity-promoting approaches can be employed to discover parsimonious nonlinear dynamical systems to have a good trade-off for model complexity and accuracy [34, 35]. Consider the problem of discovering nonlinear systems of the form:

\[\dot{\mathbf{x}}(t)=\mathbf{f}(\mathbf{x}(t)), \tag{1}\]

where \(\mathbf{x}(t)=\left[\mathbf{x}_{1}(t),\mathbf{x}_{2}(t),\dots,\mathbf{x}_{n}( t)\right]^{\top}\in\mathbb{R}^{n}\) denotes the state at time \(t\), and \(\mathbf{f}(\mathbf{x}):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}\) is a nonlinear function of the state \(\mathbf{x}(t)\).

Towards discovering the function \(\mathbf{f}\) in (1), which defines the vector field or dynamics of the underlying system, we start by collecting time-series data of state \(\mathbf{x}(t)\). Let us further assume to have time derivative information of the state. If it is not readily available, we can approximate it using numerical methods, e.g., a finite difference scheme. Thus, consider that the data \(\{\mathbf{x}(t_{0}),\dots,\mathbf{x}(t_{\mathcal{N}})\}\) and its derivative \(\{\dot{\mathbf{x}}(t_{0}),\dots,\dot{\mathbf{x}}(t_{\mathcal{N}})\}\) are given. In the next step, we assemble the data in matrices as follows:

\[\mathbf{X}=\begin{bmatrix}\mathbf{x}(t_{1})^{\top}\\ \mathbf{x}(t_{2})^{\top}\\ \vdots\\ \mathbf{x}(t_{\mathcal{N}})^{\top}\end{bmatrix}=\begin{bmatrix}\mathbf{x}_{1} (t_{1})&\mathbf{x}_{2}(t_{1})&\cdots&\mathbf{x}_{n}(t_{1})\\ \mathbf{x}_{1}(t_{2})&\mathbf{x}_{2}(t_{2})&\cdots&\mathbf{x}_{n}(t_{2})\\ \vdots&\vdots&\ddots&\vdots\\ \mathbf{x}_{1}(t_{\mathcal{N}})&\mathbf{x}_{2}(t_{\mathcal{N}})&\cdots& \mathbf{x}_{n}(t_{\mathcal{N}})\end{bmatrix}, \tag{2}\]

where each row represents a snapshot of the state. Similarly, we can write the time derivative as follows:

\[\dot{\mathbf{X}}=\begin{bmatrix}\dot{\mathbf{x}}(t_{1})^{\top}\\ \dot{\mathbf{x}}(t_{2})^{\top}\\ \vdots\\ \dot{\mathbf{x}}(t_{\mathcal{N}})^{\top}\end{bmatrix}=\begin{bmatrix}\dot{ \mathbf{x}}_{1}(t_{1})&\dot{\mathbf{x}}_{2}(t_{1})&\cdots&\dot{\mathbf{x}}_{n }(t_{1})\\ \dot{\mathbf{x}}_{1}(t_{2})&\dot{\mathbf{x}}_{2}(t_{2})&\cdots&\dot{\mathbf{x}} _{n}(t_{2})\\ \vdots&\vdots&\ddots&\vdots\\ \dot{\mathbf{x}}_{1}(t_{\mathcal{N}})&\dot{\mathbf{x}}_{2}(t_{\mathcal{N}})& \cdots&\dot{\mathbf{x}}_{n}(t_{\mathcal{N}})\end{bmatrix}. \tag{3}\]

The next key building block in the SINDy algorithm is the construction of a dictionary \(\Theta(\mathbf{y})\), containing candidate basis functions (e.g., constant, polynomial or trigonometric functions). For instance, our dictionary matrix can be given as follows:

\[\Theta(\mathbf{X})=\begin{bmatrix}\begin{array}{cccccc}\vline&\vline& \vline&\vline&\vline&\vline&\vline&\vline&\vline&\vline\\ \vline&\vline&\vline&\vline&\vline&\vline&\vline&\vline&\vline& \vline\\ \vline&\vline&\vline&\vline&\vline&\vline&\vline&\vline&\vline&\vline& \vline&\vline\\ \end{array}, \tag{4}\]

assume \(\Theta(\mathbf{X})\in\mathbb{R}^{m\times D}\), and in the above formulations polynomial terms are denoted by \(\mathbf{X}^{\mathbb{P}_{2}}\) or \(\mathbf{X}^{\mathbb{P}_{3}}\); to be more descriptive \(\mathbf{X}^{\mathbb{P}_{2}}\) denotes the quadratic nonlinearities of the state \(\mathbf{X}\) as follows:

\[\mathbf{X}^{\mathbb{P}_{2}}=\begin{bmatrix}\mathbf{x}_{1}^{2}(t_{1})&\mathbf{ x}_{1}(t_{1})\mathbf{x}_{2}(t_{1})&\cdots&\mathbf{x}_{2}^{2}(t_{1})&\mathbf{x}_{2}(t_{1} )\mathbf{x}_{3}(t_{1})&\cdots&\mathbf{x}_{2}^{2}(t_{1})\\ \mathbf{x}_{1}^{2}(t_{2})&\mathbf{x}_{1}(t_{2})\mathbf{x}_{2}(t_{2})&\cdots& \mathbf{x}_{2}^{2}(t_{2})&\mathbf{x}_{2}(t_{2})\mathbf{x}_{3}(t_{2})&\cdots& \mathbf{x}_{n}^{2}(t_{2})\\ \vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\ \mathbf{x}_{1}^{2}(t_{\mathcal{N}})&\mathbf{x}_{1}(t_{\mathcal{N}})\mathbf{x}_ {2}(t_{\mathcal{N}})&\cdots&\mathbf{x}_{2}^{2}(t_{\mathcal{N}})&\mathbf{x}_{2}( t_{\mathcal{N}})\mathbf{x}_{3}(t_{\mathcal{N}})&\cdots&\mathbf{x}_{n}^{2}(t_{ \mathcal{N}})\end{bmatrix}.\]

In this setting, each column of the dictionary \(\Theta(\mathbf{x})\) denotes a candidate function in defining the function \(\mathbf{f}(\mathbf{x})\) in (1). We are interested in identifying a few candidate functions from the dictionary \(\Theta\) so that a weighted sum of these selected functions can describe the function \(\mathbf{f}\). For this, we can set up a sparse regression formulation to achieve this goal. Precisely, we seek to identify a sparse vector \(\Xi=[\xi_{1},\ \xi_{2},\dots,\ \xi_{n}]\), where \(\xi_{i}^{\top}\in\mathbb{R}^{m}\) with \(m\) denoting the number of columns in \(\Theta\), that determines which features from the dictionary are active and their corresponding coefficients.

The SINDy algorithm formulates the sparse regression problem as an optimization problem as follows. Given a set of observed data \(\mathbf{X}\) and the corresponding time derivatives \(\dot{\mathbf{X}}\), the goal is to find the sparsest matrix \(\Xi\) that fulfills the following:

\[\dot{\mathbf{X}}=\Theta(\mathbf{X})\Xi.\]
However, finding such a matrix is an NP hard problem. Therefore, there is a need to come up a sparsity promoting regularization, and in this category, LASSO is a widely known approach [36, 37]. Despite its success, it is unable to yield matrix which is the sparsest, or the approaches, e.g., discussed in [38, 39], require a prior information about how many non-zeros elements are expected in the matrix \(\Xi\), which is not known. On the other hand, the authors in [11] discuss a sequential thresholding approach, where simple least squares problems are solved iteratively, and at each step, coefficients below a given tolerance are pruned. Analysis of such an algorithm is discussed in [40]. We summarize the SINDy approach in Algorithm 1. Moreover, we mention that other regularization schemes or heuristics are discussed in [21, 22, 25, 41, 42], but in this work, we focus only on the sequential thresholding approach, similar to in Algorithm 1 due to its simplicity.

```
0: Dictionary \(\Theta\), time-series data \(\mathbf{X}\), time derivative information \(\mathbf{X}\), threshold value tol, and maximum iterations max-iter.
0: Estimated coefficients \(\Xi\) that define governing equations for nonlinear systems.
1:\(\Xi=(\Theta^{\top}\Theta)\backslash\Theta^{\top}\mathbf{\hat{X}}\)\(\triangleright\) For initial guess, solving a least-squares problem
2:\(k=1\)
3:while\(k<\texttt{max-iter}\)do
4:\(\texttt{small\_inds}=(\texttt{abs}(\Xi)<\texttt{tol})\)\(\triangleright\) identifying small coefficients
5:\(\Xi(\texttt{small\_inds})=0\)\(\triangleright\) excluding small coefficients
6: Solve \(\Xi=(\Theta^{\top}\Theta)\backslash\Theta^{\top}\mathbf{\hat{X}}\) subject to \(\Xi(\texttt{small\_inds})=0\)
7:\(k=k+1\)
```

**Algorithm 1** SINDy algorithm [29]

## 3 iNeural-SINDy: Neural Networks and Integrating Schemes Assisted SINDy Approach

A challenge in the classical SINDy approach discussed in the previous section is the availability of an accurate estimate of the derivative information. If the derivative information is inaccurate, the resulting sparse model may not accurately capture the underlying system dynamics.

In this section, we present an approach that combines SINDy framework with a numerical integration scheme and neural networks in a particular way so that a robust discovery of governing equations can be made amid poor signal-to-noise ratio and irregularities in data. The methodology is inspired by the work [33]. The main components of the methodology are as follows. For given noisy data, we aim to learn an implicit representation using a neural network based on the noisy data so that the network yields denoised data but still in the vicinity of the noisy collected data, and governing equations describing the dynamics of the denoised data can be obtained by employing SINDy. For SINDy, we utilize automatic differential tools to obtain the derivative information via the network and also make use of an integral form of the differential equations. In the following, we make these discussions more precise.

Consider noisy data \(\mathbf{y}(t)\in\mathbb{R}^{n}\) at the time instances \(\{t_{0},\ldots,t_{\mathcal{N}}\}\), i.e., \(\{\mathbf{y}(t_{0}),\ldots,\mathbf{y}(t_{\mathcal{N}})\}\). Moreover, \(\mathbf{y}(t)=\mathbf{x}(t)+\epsilon(t)\), where \(\mathbf{x}(t)\) and \(\epsilon(t)\) denote clean data and noise, respectively. Under this setting, we aim to discover the structure of vector field \(\mathbf{f}\) by identifying the most active terms in the dictionary \(\Theta\) so that it satisfies as follows:

\[\dot{\mathbf{x}}(t)=\mathbf{f}(\mathbf{x}). \tag{5}\]

Note that we do not know \(\epsilon\)'s. In order to learn \(\mathbf{f}\) from \(\mathbf{y}\), we blend three ingredients together, which are discussed in the following.

1. **Sparse regression assumption:** In our setting, we utilize the principle of SINDy, which we discussed in the previous section. This means that the system dynamics (or vector field defining dynamics) can be represented by a few suitable terms from a dictionary of candidate functions. This allows us to obtain a parsimonious representation of dynamical systems and reduces the model complexity, leading to better generalization and interpretability of the models.
2. **Automatic differential to estimate derivative information:** As mentioned earlier, SINDy algorithm requires accurate derivative information for the system, which can be challenging
to obtain from experiments or to estimate using numerical methods. To cope with this issue, we make use of neural networks with its automatic differentiation (AD) feature, which is a technique used to estimate derivative information. The use of a DNN in combination with the SINDy algorithm was earlier discussed in [31], where it has been shown that the discovery of nonlinear system dynamics without explicit need of accurate derivative information [31]. We make use of a DNN to parameterize a nonlinear mapping from time \(t\) to the dependent variable \(\mathbf{y}(t)\). To that end, let us denote a DNN by \(\mathcal{G}_{\theta}\), where \(\theta\) contains DNN parameters. The input to \(\mathcal{G}_{\theta}\) is time \(t\), and its output is \(\mathbf{y}(t)\), i.e., \(\mathbf{y}(t)=\mathcal{G}_{\theta}(t)\). However, in the case of noisy measurement \(\mathbf{y}(t)\) at time \(\{t_{0},\ldots,t_{\mathcal{N}}\}\), we expect \(\mathcal{G}_{\theta}\) to predict outputs in the proximity of \(\mathbf{y}\), i.e., \[\mathbf{y}(t)\approx\mathbf{x}(t)=\mathcal{G}_{\theta}(t),\quad t=\{t_{0}, \ldots,t_{\mathcal{N}}\}.\] With the sparse regression hypothesis, we aim to learn a dynamical model for \(\mathbf{x}\), as it can be seen as a denoised version of \(\mathbf{y}\). For this, we construct a dictionary of possible candidate functions using \(\mathbf{x}\), which we denote by \(\Theta\big{(}\mathbf{X}(t)\big{)}\). Next, we require the derivative information of \(\mathbf{x}\) with respect to time \(t\). Since we have an implicit representation of \(\mathbf{x}(t)\) using a DNN, we can employ AD to obtain the required information. Having dictionary and derivative information, we set up a sparse regression problem as follows: \[\dot{\mathbf{X}}(t)=\Theta\big{(}\mathbf{X}(t)\big{)}\Xi,\] (6) where \(\Xi\) is the sparsest possible matrix, which selects the most active terms from the dictionary to define dynamics. Finding the sparsest solution is computationally infeasible; we, thus, utilize the sequential thresholding approach as discussed in Algorithm 1 with minor modifications. Instead of solving least-squares problems in Steps 1 and 5 in Algorithm 1, we have a loss function as follows: \[\mathcal{L}:=\min_{\theta,\Xi}\sum_{i=0}^{\mathcal{N}}\lambda_{1}\|\mathbf{y}( t_{i})-\mathbf{x}(t_{i})\|+\lambda_{2}\|\dot{\mathbf{x}}(t)-\Theta\big{(} \mathbf{x}(t)\big{)}\Xi\|,\] (7) where \(\mathbf{x}(t_{i}):=\mathcal{G}_{\theta}(t_{i})\), and \(\lambda_{1}\) and \(\lambda_{2}\) are hyperparameters.
3. **Numerical integration scheme:** A dynamical system is a particle or an ensemble of particles whose state varies over time and thus obeys differential equations involving time derivatives [43]. To predict the evolution of the dynamical system, it is necessary to have an analytical solution of such equations or their integration over time through computer simulations. Therefore, we aim to incorporate the information contained in the form of integration of dynamical systems while discovering governing equations via sparse regression, which is expected to make the process of discovering equations robust to the noise and scarcity of data. When differential equations are written in an integral form, then we do not require derivative information as well; however, the resulting optimization problem involves an integral form. In this regard, one can employ the principle of Neural-ODEs [44] to solve efficiently such optimization problems. One can also approximate the integral form using suitable integrating schemes [45], and recently, fourth-order Runge-Kutta (RK4) scheme [21] and linear multi-step methods [46] are combined with SINDy. In this work, we make use of the RK4 scheme to approximate an integral. Following [21], our goal is to predict the state of a dynamic system \(\mathbf{x}(t_{k+1})\) at time \(t=t_{k+1}\) from the state \(\mathbf{x}(t_{k})\) at time \(t=t_{k}\), where \(k\in\{0,1,\ldots,\mathcal{N}-1\}\). By employing the RK4 scheme, \(\mathbf{x}(t_{k+1})\) can be computed as a weighted sum of four components that are the product of the time-step and gradient field information \(\mathbf{f}(\cdot)\) at the specific locations. These components are computed as follows: \[\mathbf{x}(t_{k+1})\approx\mathbf{x}(t_{k})+\frac{1}{6}h_{k}(\mathbf{a}_{1}+2 \cdot\mathbf{a}_{2}+2\cdot\mathbf{a}_{3}+\mathbf{a}_{4}),\quad h_{k}=t_{k+1}-t _{k},\] (8) where, \[\mathbf{a}_{1}=\mathbf{f}(\mathbf{x}(t_{k})),\quad\mathbf{a}_{2}=\mathbf{f} \Big{(}\mathbf{x}(t_{k})+h_{k}\frac{\mathbf{a}_{1}}{2}\Big{)},\quad a_{3}= \mathbf{f}\Big{(}\mathbf{x}(t_{k})+h_{k}\frac{\mathbf{a}_{2}}{2}\Big{)},\quad \mathbf{a}_{4}=\mathbf{f}\Big{(}\mathbf{x}(t_{k})+h_{k}\mathbf{a}_{3}\Big{)}.\]
For the sake of simplicity with a slight abuse of a notation, the right-hand side of (8) is denoted by \(\mathcal{F}_{\texttt{Rk4}}\big{(}f,\mathbf{x}(t_{k}),h_{k}\big{)}\), i.e.,

\[\mathbf{x}(t_{k+1})=\mathbf{x}(t_{k}+h_{k})\approx\mathcal{F}_{\texttt{Rk4}} \big{(}\mathbf{f},\mathbf{x}(t_{k}),h_{k}\big{)}. \tag{9}\]

Like the SINDy algorithm, we collect samples from the dynamical system at time \(t=\{t_{0},\ldots,t_{\mathcal{N}}\}\) and define the time step as \(h_{k}:=t_{k+1}-t_{k}\).

With sparse regression assumption, we can write \(\mathbf{f}(\mathbf{x})=\Theta(\mathbf{x})\Xi\), where \(\Theta(\mathbf{x})\) is a dictionary and \(\Xi\) is a sparse matrix. Then, we can set up a sparse regression as follows. We seek to identify the sparsest matrix \(\Xi\) so that the following is minimized:

\[\sum_{k}\left\|\mathbf{x}(t_{k+1})-\mathcal{F}_{\texttt{Rk4}}\big{(}\Theta( \mathbf{x})\Xi,\mathbf{x}(t_{k}),h_{k}\big{)}\right\|.\]

When the RK4 scheme is merged with the previously discussed DNN framework, we apply a one-time ahead prediction based on RK4-SINDy to the output of our DNN, i.e.,

\[\mathbf{x}_{\texttt{Rk4}}(t_{k+1})\approx\mathcal{F}_{\texttt{Rk4}}\Big{(} \mathbf{f},\mathbf{x}(t_{k}),h_{k}\Big{)}.\]

Having all these ingredients, we combine them to define a loss function to train our DNN structure, as well as to discover governing equations describing underlying dynamics. To that end, we have the following loss function:

\[\mathcal{L}=\mu_{1}\mathcal{L}_{\texttt{MSE}}+\mu_{2}\mathcal{L}_{\texttt{ deri}}+\mu_{3}\mathcal{L}_{\texttt{Rk4}},\quad\mu_{1},\mu_{2},\mu_{3}\in[0,1], \tag{10}\]

where \(\mathcal{L}_{\texttt{MSE}}\) is the mean square error (MSE) of the output of the DNN \(\mathcal{G}_{\theta}\) (denoted by \(\mathbf{\hat{x}}\)) with respect to the collected data \(\mathbf{y}\), and \(\{\mu_{1},\mu_{2},\mu_{3}\}\) are positive constants, determining the weight of different losses in the total loss function. It is given as

\[\mathcal{L}_{\texttt{MSE}}=\frac{1}{\mathcal{N}}\sum_{k=1}^{\mathcal{N}} \left\|\mathbf{y}(t_{k})-\mathbf{x}(t_{k})\right\|_{2}^{2}. \tag{11}\]

It forces the DNN to produce output in the vicinity of the measurements, and \(\mu_{1}\) is its weight. \(\mathcal{L}_{\texttt{deri}}\) is inspired by the sparse regression and aims to compute the sparse coefficient matrix \(\Xi\). It is computed as follows:

\[\mathcal{L}_{\texttt{deri}}=\frac{1}{\mathcal{N}}\sum_{k=1}^{\mathcal{N}} \left\|\hat{\mathbf{x}}(t_{k})-\Theta\big{(}\mathbf{x}(t_{k})\big{)}\Xi \right\|_{2}^{2}, \tag{12}\]

The term \(\mathcal{L}_{\texttt{Rk4}}\) encodes the capabilities of the vector field to predict the state at the next time step. This is the MSE of the output of the RK4 scheme and the output of DNN, given as follows:

\[\mathcal{L}_{\texttt{Rk4}}=\frac{1}{\mathcal{N}-1}\sum_{k=1}^{\mathcal{N}} \left\|\frac{1}{h_{k}}\left(\mathbf{x}(t_{k+1})-\mathcal{F}_{\texttt{Rk4}} \left(\Theta\big{(}\mathbf{x}(t_{k})\big{)}\Xi,\mathbf{x}(t_{k}),h_{k}\big{)} \right)\right\|_{2}^{2}. \tag{13}\]

It is worth highlighting that the coefficient matrix \(\Xi\) will be updated alongside the weights and biases of the DNN, and the dictionary terms are calculated by (4). Furthermore, after a certain number of epoch training, we employ sequential thresholding on \(\Xi\) to remove small coefficients as sketched in Algorithm 1, and update the remaining parameters thereafter. We summarize the procedure in Algorithm 2. Additional steps in Algorithm 2 are as follows. We train our network for initial iterations (denoted by init-iter) without employing sequential thresholding; this helps the DNN to learn the underlying dynamics of the dataset. Afterward, we employ sequential thresholding every \(q\) iterations. In the rest of the paper, the proposed methodology is referred to as iNeural-SINDy.
```
0: Data set \(\{\mathbf{y}(t_{0}),\mathbf{y}(t_{1}),\ldots,\mathbf{y}(t_{\mathcal{N}})\}\), tol for sequential thresholding, a dictionary containing candidate functions \(\Theta\), a neural network \(\mathcal{G}_{\theta}\) (parameterized by \(\theta\)), initial iteration (init-iter), maximum iterations max-iter, and parameters \(\{\mu_{1},\mu_{2},\mu_{3}\}\).
0: Estimated coefficients \(\Xi\), defining governing equations.
1: Initialize the DNN module parameters, and the coefficients \(\Xi\)
2:\(k=1\)
3:while\(k<\texttt{max-iter}\)do
4: Feed time \(t_{i}\) as an input to the DNN (\(\mathcal{G}_{\theta}\)) and predict output \(\mathbf{x}\).
5: Compute the derivative information \(\hat{\mathbf{x}}\) using automatic differentiation.
6: Compute the cost function (10).
7: Update the parameters of DNN (\(\theta\)) and the coefficient \(\Xi\) using gradient descent.
8:for\(k\%q=0\) & \(k>\texttt{init-iter}\)do\(\triangleright\) Employing sequential thresholding after \(q\) iterations
9:\(\texttt{small\_inds}=(\texttt{abs}(\Xi)\)\(<\texttt{tol})\)\(\triangleright\) identifying small coefficients
10:\(\Xi(\texttt{small\_inds})=0\)\(\triangleright\) excluding small coefficients
11: Update the parameters of DNN (\(\theta\)) and the coefficient \(\Xi\) using gradient descent, while ensuring \(\Xi(\texttt{small\_inds})\) remains zero.
12:\(k=k+1\).
```

**Algorithm 2** iNeuralSINDy: SINDy combined with neural network and integral scheme for nonlinear system identification.

## 4 Extension to Multi-trajectories Data

Thus far, we have presented the discovery of governing equations using a single trajectory time series data set using a single initial condition. However, for complex dynamical processes, a single trajectory is not sufficient to describe underlying dynamics completely. Therefore, it is necessary to collect data using multiple trajectories; hence, we need to adopt our proposed methodology to account for multiple trajectories.

To achieve this goal, we augment the input time \(t\) with an initial condition so that a DNN can capture the nonlinear behavior of the system with respect to different initial conditions. To that end, let us consider \(\mathcal{M}\) different trajectories with initial conditions \(y_{0}^{[j]}\), where \(j\in\{1,\ldots,\mathcal{M}\}\). To reflect the multi-trajectories in our framework, we modify the architecture of the DNN, which now takes \(t_{k}\) and \(y_{0}^{[j]}\) as inputs, and intend to predict \(y_{k}^{[j]}\)--that is, the state at time \(t_{k}\) with respect to the initial condition \(y_{0}^{[j]}\). Then, we also adapt our loss function (10) as follows:

\[\mathcal{L}=\mu_{1}\sum_{j=1}^{\mathcal{M}}\mathcal{L}_{\texttt{BSE}}^{[j]}+ \mu_{2}\sum_{j=1}^{\mathcal{M}}\mathcal{L}_{\texttt{deri}}^{[j]}+\mu_{3}\sum _{j=1}^{\mathcal{M}}\mathcal{L}_{\texttt{RRA}}^{[j]},\quad\mu_{1},\mu_{2},\mu _{3}\in[0,1], \tag{14}\]

where

\[\mathcal{L}_{\texttt{BSE}}^{[j]} =\frac{1}{\mathcal{M}\cdot\mathcal{N}}\sum_{j=1}^{\mathcal{M}} \sum_{k=1}^{\mathcal{N}}\left\|\mathbf{y}^{[j]}(t_{k})-\mathbf{x}^{[j]}(t_{k} )\right\|_{2}^{2},\] \[\mathcal{L}_{\texttt{deri}}^{[j]} =\frac{1}{\mathcal{M}\cdot\mathcal{N}}\sum_{j=1}^{\mathcal{M}} \sum_{k=1}^{\mathcal{N}}\left\|\Theta\big{(}\mathbf{x}^{[j]}(t_{k})\big{)} \hat{\Xi}-\hat{\mathbf{x}}(t_{k})\right\|_{2}^{2},\] \[\mathcal{L}_{\texttt{RRA}}^{[j]} =\frac{1}{h}\frac{1}{\mathcal{M}\cdot\mathcal{N}}\sum_{j=1}^{ \mathcal{M}}\sum_{k=1}^{\mathcal{N}}\left\|\mathbf{x}^{[j]}(t_{k+1})-\mathbf{x }_{\texttt{RRA}}^{[j]}(t_{k})\right\|_{2}^{2}\ \ \text{with}\ \ h=t_{k+1}-t_{k}.\]

We depict a schematic diagram of our proposed approach in Figure 1 for such a case.

## 5 Numerical Experiments

In this section, we demonstrate the proposed methodology, the so-called iNeural-SINDy, by means of several numerical examples and present a comparison with existing methodologies. For the comparison, we primarily consider two approaches, namely DeePyMoD[31], and RK4-SINDy[21]. DeePyMoD utilizes only automatic-differential tools to estimate derivative information by constructing an implicit
representation of the noisy data, while RK4-SINDy embeds a numerical integration scheme to avoid computation of derivative information. The proposed methodology iNeural-SINDy can be viewed as a combination of DeePyMoD and RK4-SINDy. For the chaotic Lorenz example, we also present a comparison with Weak-SINDy [28].

To quantify the performance of the considered methodologies, we define the following coefficient error measure for each state variable \(\mathbf{x}_{i}\):

\[\mathcal{E}(\mathbf{x}_{i})=\left\|\Xi_{\mathbf{x}_{i}}^{\texttt{truth}}-\Xi_{ \mathbf{x}_{i}}^{\texttt{set}}\right\|_{1}, \tag{15}\]

where \(\Xi_{\mathbf{x}_{i}}^{\texttt{truth}}\) and \(\Xi_{\mathbf{x}_{i}}^{\texttt{set}}\) are, respectively, the true and estimated coefficients, corresponding to the state variable \(\mathbf{x}_{i}\), and \(\|\cdot\|_{1}\) denotes the \(l_{1}\)-norm. A motivation to quantity each state variable separately is that their dynamics can be of different scales; thus, their coefficients might also be in a different order. Therefore, to better understand the quality of the discovered models, we analyze them separately. Furthermore, to observe the performance of the methodologies under the noisy data, which is often the case in real-world scenarios, we artificially generate noisy data by corrupting the clean data. For this, we use a white Gaussian noise \(\mathcal{N}(\mu,\sigma^{2})\) with a zero mean \(\mu=0\) and variance \(\sigma^{2}\), where \(\sigma\) denotes the standard deviation. The noise level in the data is controlled by \(\sigma\), i.e., larger \(\sigma\) implies more noise present in the data. Additionally, since iNeural-SINDy and DeePyMoD both involve neural networks, we also compare their performance sensitivity in two scenarios as follows:

* Scene_A: In the first scenario, we consider having a single initial condition and a fixed number of neurons in the hidden layers but vary the amount of training data and noise levels.
* Scene_B: In the second one, we consider having a single initial condition and a fixed number of training data but vary the number of neurons in the hidden layers and noise levels.

In addition, in the following, we further clarify common implementation and reproducibility details that are considered for all the examples.

**Data generation.** We have generated the data synthetically by using solve_ivp function from scipy.integrate package to solve a given set of differential equations and produce the data set. When an identification approach terminates based on the considered methodologies (e.g., iNeural-SINDy, DeePyMoD, RK4-SINDy, or Weak-SINDy), we multiply the dictionary \(\Theta\) by estimated coefficient matrix

Figure 1: A schematic diagram of the approach iNeural-SINDy. (a) noisy measurement data, (b) feeding the initial condition \(\left(\mathbf{y}_{1,0},\ \mathbf{y}_{2,0}\right)\) and the time \(t\) to the DNN, (c) using the output of the DNN, construct a polynomial dictionary, (d) estimating the parameters of the DNN and sparse vector \(\Xi\) by considering a loss function.
\(\Xi^{\texttt{set}}\) to obtain the discovered governing equations. We then make use of the solve_ivp function from scipy.integrate to obtain time-evolution dynamics.

Moreover, we perform a data-processing step before feeding to a neural network by mapping the minimum and maximum values to \(-1\) and \(1\), respectively. The hyper-parameters \(\mu\)'s in (14) are set to \(\mu_{1}=1\), \(\mu_{2}=0.1\) and \(\mu_{3}=0.1\) for iNeural-SINDy. Note that we can drive RK4-SINDy and DeePyMoD approaches by setting \(\mu_{3}=0\) and \(\mu_{2}=0\), respectively, in (14).

Architecture.We use multi-layer perception networks with periodic activation functions, namely, SIREN[47], to learn an implicit representation based on measurement data. The numbers of hidden layers and neurons will be discussed for each example separately.

Hardware.For training neural networks and parameter estimations for discovering governing equations, we have used Nvida(r)RTX A4000 GPU with 16 GB RAM, and for CPU computations (e.g., for generating data), we have used a 12th Gen Intel(r) Core(tm)i5-12600K processor with 32 GB RAM.

Training set-up.We use the Adam optimizer [48] to update the coefficient matrix \(\Xi\) that is trained alongside the DNN parameters. The threshold value (tol), learning rate of the optimizer, maximum iterations(max-iter), initial iterations (init-iter), the iteration \(q\) for employing sequential thresholding for Algorithm 2 will be mentioned for each example separately.

However, we note that after each thresholding step in Algorithm 2, we reset the learning rate \(5\times 10^{-6}\) for DNN parameters and \(1\times 10^{-2}\) for the coefficient matrix \(\Xi^{\texttt{set}}\) except for the Lorenz example, which is explicitly mentioned in the Lorenz example.

### Two-dimensional damped oscillators

In our first example, we consider the discovery of a two-dimensional linear oscillatory damped system using data. The dynamics of the oscillator can be given by

\[\begin{split}\dot{\mathbf{x}}_{1}(t)&=-0.1\mathbf{x }_{1}(t)+2.0\mathbf{x}_{2}(t),\\ \dot{\mathbf{x}}_{2}(t)&=-2.0\mathbf{x}_{1}(t)-0.1 \mathbf{x}_{2}(t).\end{split} \tag{16}\]

Simulation setup:To generate the training data set, we consider three initial conditions in the range \([-2,2]\) for \(\mathbf{x}_{1}\) and \(\mathbf{x}_{2}\), and for each initial condition, we take \(400\) equidistant points in the time interval \(t\in[0,\ 10]\). Our DNN architecture has three hidden layers, each having \(32\) neurons. We set the number of epochs \(\texttt{max-iter}=15,000\) and threshold value \(\texttt{tol}=0.05\). The initial iteration init-iter is set to \(5,000\) with the learning rate of \(10^{-4}\) for the DNN parameters and \(10^{-3}\) for the coefficient matrix \(\Xi^{\texttt{set}}\), and after \(q=2,000\) iterations, we employ the sequential thresholding. Moreover, we construct a dictionary containing polynomials of degrees up to two.

Results:Figure 2 demonstrates the performance of different algorithms in the presence of noise. We consider additive white Gaussian noise with different standard variances \(\sigma=\{0,\ 0.02,\ 0.04,\ 0.08\}\). It shows that as we increase the noise level, the RK4-SINDy fails to estimate the coefficients. However, iNeural-SINDy and DeePyMoD are robust in discovering the underlying equations accurately, even for high noise levels, and both exhibit similar performances. In Table 1 (in the appendix), we also report learned governing equations from data with various noise levels, which again illustrate that both iNeural-SINDy and DeePyMoD have similar performance, and RK4-SINDy fails to recover governing equations from highly noisy data. Furthermore, in Figure 3, the convergence of the non-zero coefficients for the different methods is shown as the training progresses. It can be seen that iNeural-SINDy has a faster convergence rate compared to DeePyMoD and RK4-SINDy. Next, we discuss the performance of iNeural-SINDy and DeePyMoD for Scene_A and Scene_B.

* Scene_A: We consider a DNN architecture with three hidden layers, each having \(32\) neurons. For comparison, we consider noise levels with standard variance \(\sigma=\{0,\ 0.02,\ 0.04,\ 0.06\}\), and take the number of samples \(\{30,\ 40,\ 50,\ 100,\ 200,\ 300,\ 400\}\) in the time interval \([0,10]\) for a single initial condition \((\mathbf{x}_{1}(0),\mathbf{x}_{2}(0))=(5,2)\). The rest of the settings are the same as mentioned earlier in the simulation setup. By varying the noise levels and the number of samples, we report
the quality of the learned governing equations in Figure 4. Note that the error criterion defined in (15) is used. Each cell shows the error corresponding to sample sizes and noise levels. By comparing the simulation results, we notice that DeePyMoD performs better for low data regime, but as the number of data is increased, both iNeural-SINDy and DeePyMoD perform similarly.
* Scene_B: In this case, we consider a DNN architecture with three hidden layers but vary the number of neurons at each layer from 2 to 64. Again, we consider various noise levels. We take 400 samples in the time interval \([0,10]\) for a single arbitrary initial condition \((\mathbf{x}_{1}(0),\mathbf{x}_{2}(0))=(5,2)\). The rest of the settings are the same as mentioned earlier in the simulation setup. By varying the noise levels and number of neurons, we report a comparison between iNeural-SINDy and DeePyMoD in Figure 5, where each cell shows the error, corresponding to a specific number of neurons and noise level. These comparisons again show that both methodologies perform comparably and learn correct coefficients with similar performance for a large number of neurons as the DNN has more capacities to capture the dynamics present in the data. More interesting, we would like to highlight that both methods do not over-fit as the capacity of the DNN is increased.

### Cubic damped oscillator

The cubic oscillatory system is given by the following equation:

Figure 3: Linear oscillator: Estimated coefficients during the training loop for iNeural-SINDy, DeePyMoD and RK4-SINDy.

Figure 2: Linear oscillator: A comparison of the learned equations using different methods under various noise levels present in measurement with the ground truth.
\[\dot{\mathbf{x}}_{1}(t) =-0.1\mathbf{x}_{1}^{3}(t)+2.0\mathbf{x}_{2}^{3}(t), \tag{17}\] \[\dot{\mathbf{x}}_{2}(t) =-2.0\mathbf{x}_{1}^{3}(t)-0.1\mathbf{x}_{2}^{3}(t).\]

The system consists of two coupled, non-linear differential equations describing the time evolution of two variables, \(\mathbf{x}_{1}\) and \(\mathbf{x}_{2}\). Given noisy data, we aim to recover the governing equations and perform a similar analysis as done for the previous example.

Simulation setup:To generate the training data set, we consider two initial conditions \((\mathbf{x}_{1}(0),\mathbf{x}_{2}(0))=\{(2,2),(-2,-2)\}\) and collect 800 points in the time interval \(t\in[0,10]\). Our DNN architecture has three hidden layers, each having 32 neurons. We set the number of epoch max-iter\(=30,000\), threshold value \(\texttt{tol}=0.05\). The initial training iteration (init-iter) is set to \(15,000\) with the learning rate \(10^{-4}\) for the DNN parameters and \(10^{-3}\) for the coefficient matrix \(\Xi^{\texttt{next}}\). After the initial training, for

Figure 4: Linear oscillator: A comparison of iNeural-SINDy and DePyMoD under Scene_A.

Figure 5: Linear oscillator: A comparison of iNeural-SINDy and DePyMoD under Scene_B
every \(q=5,000\) iterations afterward, we employ the sequential thresholding and update the DNN parameters and \(\Xi^{\texttt{next}}\). The dynamical system is estimated in the space of polynomials up to order three.

Results:To see the performance of these different methodologies under the presence of noise, we consider a Gaussian noise with the standard variance \(\sigma=\{0,\ 0.02,\ 0.04,\ 0.06\}\). We report the obtained results in Figure 6 and in Table 2 (see Appendix), and we notice that RK4-SINDy performs poorly for high noise levels, but iNeural-SINDy and DeePyMoD have competitive performance. Further, in Figure 7, we plot the convergence of the non-zero coefficients as the training progresses for the noise-free case. Here, we again observe a faster convergence for iNeural-SINDy as compared to the other two approaches. Next, we investigate performances of iNeural-SINDy and DeePyMoD for Scene_A and Scene_B, which are discussed in the following.

* Scene_A: We fix a DNN architecture with three hidden layers, each having 32 neurons. We consider a set of noise level with \(\sigma=\{0,\ 0.02,\ 0.04,\ 0.06\}\) and a set of sample size \(\{30,40,\ 50,\ 100,\ 200,\ 300,\ 400\}\). The data are collected using a random initial condition in the interval \([1,4]\) for \(\{\mathbf{x}_{1},\mathbf{x}_{2}\}\). The rest of the settings are the same as mentioned earlier in the simulation setup. The results are shown in Figure 8, where we notice that for a smaller data set, iNeural-SINDy performs slightly better as compared to DeePyMoD, whereas for larger data set, it is otherwise.
* Scene_B: For this case, we fix the sample size to 400 but consider a DNN architecture with three hidden layers with the number of neurons ranging from 2 to 64. Furthermore, we consider a set of noise levels with \(\sigma=\{0,\ 0.02,\ 0.04,\ 0.06\}\). The data are generated as in Scene_A and the training setting is also to be as above. The results are depicted in Figure 9, where we observe that iNeural-SINDy performs better as compared to DeePyMoD for fewer neurons, and as we increase the number of neurons, both methods perform similarly.
\[\dot{\mathbf{x}}_{1}(t) =1.0\mathbf{x}_{1}(t)-1.0\mathbf{x}_{2}(t)-\frac{1}{3}\mathbf{x}_{1} ^{3}(t)+0.1, \tag{18}\] \[\dot{\mathbf{x}}_{2}(t) =0.1\mathbf{x}_{1}(t)-0.1\mathbf{x}_{2}(t).\]

Simulation setup:For this simulation example, we consider two initial conditions \((\mathbf{x}_{1}(0),\mathbf{x}_{2}(0))=\{(2,1.5),(1.5,2)\}\) and take \(400\) data points in the time interval \(t\in[0,\ 200]\). The DNN architecture has three hidden layers with \(32\) neurons. We set the number of epoch max-iter\(=50,000\) and threshold value tol\(=0.05\). The number of iterations for the initial training is set to \(15,000\) with the learning rate \(10^{-4}\) for the DNN parameters and \(10^{-3}\) for the coefficient matrix \(\Xi^{\texttt{next}}\). After the initial training, we employ the sequential thresholding after each \(q=5,000\) iterations. We aim to learn the underlying governing equations in the space of polynomials with degrees up to order three.

Results:Converse to the results that we earned in the previous two examples, for the FHN, iNeural-SINDy has a slower convergence rate compared to DeePyMoD and RK4-SINDy, see Figure 10.

For this example, we again make a similar observation (see Figure 11, and Table 3 in Appendix), where we notice that iNeural-SINDy and DeePyMoD exhibit similar performances for lower noise levels, but for the higher noise values, (see the results for \(\sigma=0.08\) Table 3), iNeural-SINDy tends to outperform DeePyMoD. Moreover, RK4-SINDy clearly fails for high noise levels. However, converse to the results reported in the previous two examples, for this example, we notice a slower convergence

Figure 8: Cubic oscillator: A comparison of iNeural-SINDy and DeePyMoD under Scene_A.

Figure 7: Cubic oscillator: Estimated coefficients during the training loop for iNeural-SINDy, DeePyMoD and RK4-SINDy.
of iNeural-SINDy compared to DeePyMoD and RK4-SINDy; see, Figure 10. But we highlight that iNeural-SINDy can identify governing equations for highly noisy data, as stated earlier. Next, we compare the performances of iNeural-SINDy and DeePyMoD under Scene_A and Scene_B.

* Scene_A: In this case, we consider a fixed DNN architecture with three hidden layers, each consisting of 32 neurons. The different noise levels \(\{0.0,\ 0.02,\ 0.04,0.06\}\) are considered, while the sample size is considered in the range from 150 to 450 with an increment of 50. Here, a single initial condition is used for data collection; that is, \((\mathbf{x}_{1}(0),\mathbf{x}_{2}(0))=(3,2)\). The rest of the settings are the same as mentioned earlier for this example. The results are shown in Figure 12, where we notice that DeePyMoD outperforms iNeural-SINDy and has a better performance.
* Scene_B: Here, we conduct a study where we keep the number of samples fixed at 400, obtained

Figure 10: Fitz-Hugh Nagumo: Estimated coefficients during the training loop for iNeural-SINDy, DeePyMoD and RK4-SINDy.

Figure 9: Cubic oscillator: A comparison of iNeural-SINDy and DeePyMoD under Scene_B
using the initial condition \((\mathbf{x}_{1}(0),\mathbf{x}_{2}(0))=(3,2)\). The DNN architecture is designed to have three hidden layers. We aim to explore how iNeural-SINDy and DeePyMoD perform under different combinations of neurons for each layer and noise level. The training settings for each case remain the same, as mentioned earlier. The outcomes are presented in the heat-map depicted in Figure 13, where we notice that both DeePyMoD and iNeural-SINDy almost have the same performance in all the settings.

### Chaotic Lorenz system

The chaotic Lorenz system is a set of three differential equations as follows [49]:

\[\dot{\mathbf{x}}_{1}(t) =\gamma\big{(}\mathbf{x}_{2}(t)-\mathbf{x}_{1}(t)\big{)}, \tag{19a}\] \[\dot{\mathbf{x}}_{2}(t) =\mathbf{x}_{1}(t)\big{(}\rho-\mathbf{x}_{3}(t)\big{)}-\mathbf{x} _{2}(t),\] (19b) \[\dot{\mathbf{x}}_{3}(t) =\mathbf{x}_{1}(t)\mathbf{x}_{2}(t)-\beta\mathbf{x}_{3}(t), \tag{19c}\]

Figure 11: Fitz-Hugh Nagumo: Comparison of the estimation with different techniques and noise level

Figure 12: Fitz-Hugh Nagumo: A comparison of iNeural-SINDy and DeePyMoD under Scene_A.
where the parameters \(\gamma\), \(\rho\), and \(\beta\) are positive constants with associated standard values \(\gamma=10,\ \rho=28,\ \beta=\frac{8}{3}\). The Lorenz system is a classic example of a chaotic system, which means that small differences in the initial conditions can lead to vastly different outcomes over time. It is a widely used benchmark example for discovering governing equations [11].

Simulation setup:We collect our data in the time interval \(t\in[0,10]\) with a sample size of \(200\) for three different initial conditions \((\mathbf{x}_{1}(0),\mathbf{x}_{2}(0),\mathbf{x}_{3}(0))=\{(-8,7,27),(-6,6,25),(-9,8,22)\}\). The DNN architecture has three hidden layers, each having \(64\) neurons. We set the number of iterations max-iter\(=35,000\), and threshold value \(\mathtt{tol}=0.2\). We set the number of iterations for the initial training to init-iter\(=10,000\), and the learning rate \(7\cdot 10^{-4}\) for the DNN parameters and \(10^{-2}\) for the coefficient matrix \(\Xi^{\mathtt{est}}\). After finishing the initial iterations, we employ sequential thresholding after every \(q=3,000\) iterations. Moreover, after each sequential thresholding, we reset the learning rate for DNN parameters to \(5\cdot 10^{-6}\) and for the coefficient matrix \(\Xi^{\mathtt{est}}\) to \(10^{-2}\). The governing equations are estimated by constructing a dictionary with polynomials up to degree two.

Since the magnitude of \(\{\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3}\}\) for the Lorenz example can be large, we consider scaling the \(\mathbf{x}\)'s using a scaling factor \(\alpha\). Note that such scaling does not affect the interaction between different \(\mathbf{x}\)'s; thus, the sparsity pattern remains the same as well. However, it is observed that improving the condition number of the dictionary matrix enhances the estimate of the coefficients and helps us to determine the right governing equations.

Results:We conduct experiments using a scaling factor \(\alpha=0.1\). Further, we aim to learn governing equations from the noisy data with noise levels of \(\sigma=\{0,\ 0.04,\ 0.1,\ 0.2,\ 0.4\}\). We report the obtained results in Table 4, where we notice that iNeural-SINDy and DeePyMoD yield similar performance except for the case of higher noise level (e.g., see the results for \(\sigma=0.4\)), where iNeural-SINDy recovered the equations better. However, RK4-SINDy performs poorly for the higher noise levels. Next, we conduct a performance analysis of iNeural-SINDy and DeePyMoD for Scene_A and Scene_B. We note that in both scenarios, the training data are generated using a single initial condition \((\mathbf{x}_{1}(0),\mathbf{x}_{2}(0),\mathbf{x}_{3}(0))=(-8,\ 7,\ 27)\).

* Scene_A: We compare iNeural-SINDy and DeePyMoD under Scene_A. We also investigate the effect of the scaling factor \(\alpha\) and consider two values of it, i.e., \(\alpha=\{0.1,\ 1\}\). We fix the DNN architecture to have three hidden layers, each having \(64\) neurons. We consider different sample sizes and noise levels to compare the performance of iNeural-SINDy and DeePyMoD. For \(\alpha=0.1\),

Figure 13: Fitz-Hugh Nagumo: A comparison of iNeural-SINDy and DeePyMoD under Scene_B
we show the results in Figure 14, where we notice that iNeural-SINDy outperforms DeePyMoD in most cases. A similar observation is made for \(\alpha=1\), which is reported in Figure 15. For these experiments, it is hard to conclude the effect of the scaling factors, as we notice that in some cases, the scaling improves the performance, and in some cases, it is not the case.
* Scene_B: In this case, we fix the sample size to \(400\). We also fix the number of hidden layers for the DNN architecture to three but vary the number of neurons in each layer. We also conduct experiments to see the effect of the scaling factor in this case as well. The results for \(\alpha=1\) and \(\alpha=0.1\) are shown in Figure 16 and Figure 17, respectively. These heat maps indicate the outperformance of iNeural-SINDy in most cases. We also observe that for a larger number of neurons, the scaling factor \(\alpha=0.1\) slightly performs better compared to scaling factor \(\alpha=1\) in both iNeural-SINDy as well as DeePyMoD.

A comparison of iNeural-SINDy with WEAK-SINDy:Beside our previous comprehensive study, we next compare iNeural-SINDy with Weak-SINDy, which also does not require any estimate of

Figure 16: Lorenz example: A comparison of iNeural-SINDy and DeePyMoD under Scene_B with the scaling factor \(\alpha=0.1\).

Figure 14: Lorenz example: A comparison of iNeural-SINDy and DeePyMoD under Scene_A with the scaling factor \(\alpha=0.1\)

Figure 15: Lorenz example: A comparison of iNeural-SINDy and DeePyMoD under Scene_A with the scaling factor \(\alpha=1\)
derivatives using noisy data; for more details on Weak-SINDy, we refer to [28].

For this study, we again consider the same initial condition as used for Scene_A and Scene_B. We take 2000 data points in the time interval \([0,10]\), which are corrupted using different noise levels \(\sigma=\{0,0.02,0.08,0.1\}\). To discover the governing equations, we consider a dictionary of polynomials up to degree two. For training iNeural-SINDy we use the same setting as discussed in Section 5.4. For Weak-SINDy, we consider the code provided by the authors1. We report the results in Table 5, which indicates that iNeural-SINDy outperforms Weak-SINDy in the presence of high noise.

Footnote 1: [https://github.com/MathBioCU/WSINDy_ODE/tree/master](https://github.com/MathBioCU/WSINDy_ODE/tree/master)

## 6 Conclusions

In this work, we proposed a methodology, namely iNeural-SINDy, to discover governing equations using noisy and scarce data. It consists of three main components--these are: (a) learning an implicit representation based on given noisy data using a deep neural network, (b) setting up a sparse regression problem inspired by SINDy[11], discovering governing equations, and (c) utilize an integral form of differential equations. We have combined all these components innovatively to learn governing equations from noisy data. Particularly, we highlight that we leverage the implicit representation using neural networks to estimate the derivative using automatic differential to avoid any numerical derivative estimation using noisy data. We have shown how iNeural-SINDy can be employed when data are collected using multiple trajectories. Furthermore, we have presented an extensive comparison of the proposed methodology with RK4-SINDy[21] and DeePyMoD[31], where we noticed that iNeural-SINDy clearly out-performed RK4-SINDy, and in many cases, iNeural-SINDy also yielded better or comparable results as compared to DeePyMoD, expect for the FHN example. We also compared iNeural-SINDy with Weak-SINDy using the Lorenz example, where we noticed a better performance of iNeural-SINDy. In the future, we would like to extend the proposed framework to the identification of parametric and control-driven dynamical systems. We also like to combine the idea of the ensemble discussed in [29] to further improve the quality of learned governing equations.
