SINDy-PI: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics

Accurately modelling the nonlinear dynamics of a system from measurement data is a challenging yet vital topic. The sparse identification of nonlinear dynamics (SINDy) algorithm is one approach to discover dynamical systems models from data. Although extensions have been developed to identify implicit dynamics, or dynamics described by rational functions, these extensions are extremely sensitive to noise. In this work, we develop SINDy-PI (parallel, implicit), a robust variant of the SINDy algorithm to identify implicit dynamics and rational nonlinearities. The SINDy-PI framework includes multiple optimization algorithms and a principled approach to model selection. We demonstrate the ability of this algorithm to learn implicit ordinary and partial differential equations and conservation laws from limited and noisy data. In particular, we show that the proposed approach is several orders of magnitude more noise robust than previous approaches, and may be used to identify a class of ODE and PDE dynamics that were previously unattainable with SINDy, including for the double pendulum dynamics and simplified model for the Belousov–Zhabotinsky (BZ) reaction.

Parsimonious modelling has a rich history, with many scientific advances being argued on the basis of Occam's razor, that the simplest model is likely the correct one. SINDy exemplifies this principle, identifying a potentially nonlinear model with the fewest terms required to describe how the measurement data changes in time. The basic idea behind SINDy may be illustrated on a one-dimensional systemẋ = f (x); the general formulation for multidimensional dynamics will be described in the following sections. An interpretable form of the nonlinear dynamics may be learned by writing the rate of change of the state of the system x as a sparse linear combination of a few terms in a library of candidate functions, Θ(x) = [ θ 1 (x) θ 2 (x) ... θ p (x) ]: x(t) = f (x(t)) ≈ Θ(x(t))ξ , (1.1) where each θ j (x) is prescribed candidate term (e.g. x, x 2 , sin(x), . . .). The derivative of the state and the library of candidate functions may both be computed from measured trajectory data. It then remains to solve for a sparse vector ξ with non-zero entries ξ j indicating which functions θ j (x) are active in characterizing the dynamics. The resulting models strike a balance between accuracy and efficiency, and they are highly interpretable by construction. In a short time, the SINDy algorithm has been extended to include inputs and control [48], to identify partial differential equations [29,30], to incorporate physically relevant constraints [34], to include tensor bases [45], and to incorporate integral terms for denoising [49,50]. These extensions and its simple formulation in terms of a generalized linear model in (1.1) have resulted in SINDy being adopted in the fields of fluid mechanics [34,37], nonlinear optics [31], plasma physics [32], chemical reactions [33,36,39], numerical methods [41] and structural modelling [42]. The generalized linear model in (1.1) does not readily lend itself to representing implicit dynamics and rational functions, which are not naturally expressible as sum of a few basis functions. Instead, the implicit-SINDy algorithm [47] reformulates the SINDy problem in an implicit form: Θ(x,ẋ)ξ = 0. (1.2) This formulation is flexible enough to handle a much broader class of dynamics with rational function nonlinearities, such asẋ = N(x)/D(x) which may be rewritten asẋD(x) + N(x) = 0. However, the sparsest vector ξ that satisfies (1.2) is the trivial solution ξ = 0. Thus, the implicit-SINDy algorithm leverages a recent non-convex optimization procedure [51,52] to find the sparsest vector ξ in the null space of Θ(x,ẋ), which differs from other approaches [53,54] that identify the rational dynamics. For even small amounts of noise, the dimension of the null space will become prohibitively large, making this approach extremely sensitive to noise and compromising the model discovery process.
This work develops an optimization and model selection framework that recasts implicit-SINDy as a convex problem, making it as noise robust as the original non-implicit SINDy algorithm and enabling the identification of implicit ODEs and PDEs that were previously inaccessible. The key to making the implicit-SINDy algorithm robust is the realization that if we know even a single term in the dynamics, corresponding to a non-zero entry ξ j , then we can rewrite (1.2) in a non-implicit form θ j (x,ẋ) = Θ (x,ẋ)ξ , (1.3) where Θ and ξ have the jth element removed. Because none of these terms are known a priori, we sweep through the library, term by term, testing (1.3) for a sparse model that fits the data. This procedure is highly parallelizable and provides critical information for model selection. Our approach is related to the recent work of Zhang et al. [46], which also makes the implicit problem more robust by testing candidate functions individually. However, there are a number of key differences in the present approach. Our work explicitly considers rational nonlinearities to discover exceedingly complex implicit PDEs, such as a simplified model of the Belousov-Zhabotinsky (BZ) reaction. Our framework also provides several new greedy algorithms, including parallel and constrained formulations. We further extend this method to include the effect of control inputs, making it applicable to robotic systems [55], and we use this procedure to discover Hamiltonians. Finally, our approach provides guidance on model selection, a comprehensive comparison with previous methods, and a careful analysis of noise robustness.

Background
We briefly introduce the full multidimensional SINDy and implicit-SINDy algorithms, which will provide a foundation for our robust implicit identification algorithm in §3.

(a) Sparse identification of nonlinear dynamics
The goal of SINDy [28] is to discover a dynamical system d dt from time-series data of the state x(t) = [x 1 (t), . . . , x n (t)] T ∈ R n . We assume that the dynamics, encoded by the function f , admit a sparse representation in a library of candidate functions: Thus, each row equation in (2.1) may be written as where ξ k is a sparse vector, indicating which terms are active in the dynamics. We determine the non-zero entries of ξ k through sparse regression based on trajectory data. The time-series data is arranged into a matrix X = [ x(t 1 ) x(t 2 ) ··· x(t m ) ] T , and the associated time derivative matrixẊ = [ẋ(t1)ẋ(t2) ···ẋ(t m ) ] T is computed using an appropriate numerical differentiation scheme [1,28,56]. It is then possible to evaluate the library Θ on trajectory data in X so that each column of Θ(X) is a function θ j evaluated on the m snapshots in X.
It is now possible to write the dynamical system in terms of a generalized linear model, evaluated on trajectory data:Ẋ = Θ(X)Ξ . (2.4) There are several approaches to identify the sparse matrix of coefficients Ξ , including sequentially thresholded least squares (STLSQ) [ x data collection sparse regression over many candidate nonlinear terms selection and reconstruction sparse dense regression (SR3) [59,60], stepwise sparse regression (SSR) [36] and Bayesian approaches [46,61]. It is possible to augment the library to include partial derivatives for the identification of partial differential equations (PDEs) [29,30]. Similarly, it is possible to include external forcing terms in the library Θ, enabling the identification of forced and actively controlled systems [48]. To alleviate the effect of noise, it is possible to reframe the SINDy problem in terms of an integral formulation [49,50]. There are a number of factors that affect the robustness of SINDy, some of which are discussed in appendix I.

(b) Implicit sparse identification of nonlinear dynamics
The implicit-SINDy algorithm [47] extends SINDy to identify implicit differential equations f (x,ẋ) = 0, (2.5) and in particular, systems that include rational functions in the dynamics, such as chemical reactions and metabolic networks that have a separation of timescales. The implicit-SINDy generalizes the library Θ(X) in (2.4) to include functions of x andẋ: However, this approach requires solving for a matrix Ξ whose columns ξ k are sparse vectors in the null space of Θ(X,Ẋ). This approach is non-convex, relying on the alternating directions method (ADM) [47,52], and null space computations are highly ill-conditioned for noisy data [1,47,62], thus inspiring the current work and mathematical innovations.

SINDy-PI: robust parallel identification of implicit dynamics
We have developed the SINDy-PI (parallel, implicit) framework for the robust identification of implicit dynamics, bypassing the null space approach discussed in §2b (figure 1). The idea is that if even a single term θ j (x,ẋ) ∈ Θ(x,ẋ) in the dynamics (2.5) is known, it is possible to rewrite (2.6) as θ j (X,Ẋ) = Θ(X,Ẋ|θ j (X,Ẋ)ξ j , (3.1) where Θ(X,Ẋ|θ j (X,Ẋ)) is the library Θ(X,Ẋ) with the θ j column removed. equation ( [28][29][30]36,46,49,50,[59][60][61]. In particular, we solve for a sparse coefficient vector ξ j that minimizes the following loss function: where β is the sparsity promoting parameter. There are numerous relaxations of the nonconvex optimization problem in (3.2), for example the sequentially thresholded least-squares algorithm [28]. Because there is no null space calculation, the resulting algorithm is considerably more robust to noise than the implicit-SINDy algorithm [47], i.e. we longer have to deal with an ill-conditioned null space problem. In general, the entire point of SINDy is that the dynamics are not known ahead of time, and so it is necessary to test each candidate function θ j until one of the models in (3.1) admits a sparse and accurate solution. When an incorrect candidate term is used, then the algorithm results in a dense (non-sparse) model ξ j and an inaccurate model fit, and when a correct term is included, the algorithm identifies a sparse model ξ j and an accurate model fit. In this way, it is clear when the algorithm has identified the correct model. Moreover, there is a wealth of redundant information, since each term in the correct model may be used as the candidate function on the left-hand side, and the resulting models may be cross-referenced. This approach is highly parallelizable, and each candidate term may be tested simultaneously in parallel. The non-parallel formulation in (3.1) was recently introduced by Zhang et al. [46] in the context of Bayesian regression, where they also make the implicit problem more robust by testing candidate functions individually; however, they do not consider dynamics with rational function nonlinearities or control inputs. In this work, we extend the robust implicit formulation to identify several challenging implicit ODE and PDE systems with rational function nonlinearities, which are ubiquitous in engineering and natural systems, and systems with external forcing and control inputs. We also introduce the parallel formulation and model selection frameworks. Further, we will introduce a constrained optimization framework to simultaneously test all candidate functions.

(a) Model selection
For each candidate function in (3.1), we obtain one candidate model. When the candidate function θ j is not in the true dynamics, then the resulting coefficient vector ξ j will not be sparse and there will be large prediction error. In contrast, when a correct candidate function is selected, then we obtain a sparse coefficient vector ξ j and small prediction error. For an implicit dynamical system, there may be several different implicit equations that must be identified, resulting in several candidate functions that admit sparse models. The sequentially thresholded least-squares (STLSQ) algorithm that we use here, and whose convergence properties are considered by Zhang and Schaeffer [57], iteratively computes a least-squares solution to minimize θ j (X,Ẋ) − Θ(X,Ẋ|θ j (X,Ẋ))ξ j 2 and then zeros out small entries in ξ j that are below a set threshold λ. This threshold λ is a hyperparameter that must be tuned to select the model that most accurately balances accuracy and efficiency. Thus, we must employ model selection techniques to identify the implicit models that best supports the data, while remaining as simple as possible.
There are several valid approaches to model selection. To select a parsimonious yet accurate model we can also employ the Akaike information criterion (AIC) [63,64] and Bayesian information criterion (BIC) [65], as in [66]. It is also possible to sweep through the parameter λ and candidate functions θ j , and then choose the Pareto optimal model from a family of models on the Pareto front balancing accuracy and efficiency; this is the approach in the original SINDy work [28] and in earlier work leveraging genetic programming to discover dynamics [25,26]. In this work, we take a different approach, selecting models based on performance on a test dataset X t that has been withheld for model validation to automate the model selection process. For each threshold λ, the resulting model is validated on the test set X t , and the model with the lowest test error is selected. One error function is the model fit: In practice, for rational dynamics, we select based upon the predicted derivativeẊ t : For implicit dynamics where each state derivative may be written as a rational functioṅ then we restrict the candidate functions to θ j (x,ẋ) =ẋ k θ j (x) for some θ j (x) ∈ Θ(x) to identify a separate sparse model for eachẋ k . Several candidate functions may provide accurate and sparse models. These different models may further be cross-references to check that the same terms are being selected in each model, providing additional information for model selection and validation.

(b) Constrained optimization formulation
In (3.1) each candidate function was tested individually in a parallel optimization. However, each of these individual equations may be combined into a single constrained system of equations We constrain Ξ to have zero entries on the diagonal, as shown in figure 2, which is the same as removing the candidate function from the library in the separate optimization problems in (3.1). Without this constraint, the trivial solution Ξ = I p×p will provide the sparsest Ξ and the most accurate model. This may be written as a formal constrained optimization problem: This optimization is non-convex, although there are many relaxations that result in accurate and efficient proxy solutions [28,58,59]. In this work, we will use sequentially thresholded least squares, so that any entry Ξ ij < λ will be set to zero; the sparsity parameter λ is a hyperparameter, and each column equation may require a different parameter λ j . The constrained formulation in (3.7) can be solved efficiently in modern optimization packages, and we use CVX [67,68]. After solving (3.7) we have numerous candidate models, one for each column ξ k of Ξ , given by The sparse models that result in an accurate fit are candidate implicit models, and they may be assessed using the model selection approaches outlined above. These various models may be cross-referenced for consistency, as the same models will have the same sparsity pattern. This information can then be used to refine the library Θ, for example to only include the non-zero entries in the sparse columns of Ξ .

(c) Noise robustness
We now compare the noise sensitivity of SINDy-PI and implicit-SINDy on the one-dimensional Michaelis-Menten model for enzyme kinetics [47,69,70], given bẏ where x denotes the concentration of the substrate, j x denotes the influx of the substrate, V max denotes the maximum reaction time and K m represents the concentration of half-maximal reaction. We use the same parameters as in [47], with j x = 0.6, V max = 1.5, and K m = 0.3. Figure 3 shows the result of the noise robustness of SINDy-PI and implicit-SINDy. In this example, SINDy-PI is able to handle over 10 5 more measurement noise than implicit-SINDy, while still accurately recovering the correct model. Details are provided in appendix A, and key factors that limit robustness are discussed in appendix I.

(d) Data usage
The data required to correctly identify a model is a critical aspect when comparing SINDy-PI and implicit-SINDy. Many experimental datasets are limited in volume, and thus our goal is to  identify a model with as little data as possible. In this section, we compare the SINDy-PI and implicit-SINDy methods on the challenging yeast glycolysis model [47,71] given bẏ is the most challenging equation to discover in this system, and figure 4 compares the success rate of SINDy-PI and implicit-SINDy in identifying this equation. SINDy-PI uses about 12 times less data than the implicit-SINDy when identifying (3.10f). Details are provided in appendices B and D.

(e) Comparison for implicit PDE identification
We now investigate the ability of SINDy-PI to discover a PDE with rational terms, given by a modified KdV equation where γ u is a loss term and 2g 0 /(1 + u) is a gain term. We fix γ = 0.1 and vary the value of g 0 from 0 to 1. As g 0 increases, the implicit term gradually dominates the dynamics. Figure 5 shows the results of SINDy-PI and PDE-FIND [29] for different values of g 0 . For large g 0 , SINDy-PI is able to accurately identify the rational function term, while this is not possible for PDE-FIND, since this term is not in the library. Details of the identification process are given in appendix C.

Advanced examples
We will now demonstrate the SINDy-PI framework on several challenging examples, including the double pendulum, an actuated single pendulum on a cart, the Belousov-Zhabotinsky PDE, and the identification of conserved quantities. All examples are characterized by rational nonlinearities, and we were unable to identify them using SINDy or implicit-SINDy, even in the absence of noise.

modified KdV
loss gain prediction error In our first example, we use SINDy-PI to discover the equations of motion of a mounted double pendulum, shown in figure 6. The double pendulum is a classic example of chaotic dynamics [72], and was an original challenging example used to demonstrate the capabilities of genetic programme for model discovery [26]. Correctly modelling the nonlinear dynamics is vital for accurate control [72]. We simulate the double pendulum dynamics, derived from the Euler-Lagrange equations, and use SINDy-PI to re-discover the dynamics from noisy measurements of the trajectory data. The governing equations and SINDy-PI models are provided in appendix F. Because these dynamics have rational nonlinearities, the original SINDy algorithm is unable to identify the dynamics, making this a challenging test case. The state vector is given by x = [φ 1 , φ 2 ,φ 1 ,φ 2 ] T , and the parameters of the simulation are given in appendix D. The training data is generated from an initial condition x train = [π + 1.2, π − 0.6, 0, 0] T , simulated for 10 s using a time step of dt = 0.001 s. The validation data is generated from an initial condition x val = [π − 1, π − 0.4, 0.3, 0.4] T , simulated for 3 s with time step dt = 0.001 s.
To test the robustness of SINDy-PI, we add Gaussian noise to both the training and validation data. We test the resulting models using a new testing initial condition x test = [π + 0.3, π − 0.5, 0, 0] T . We construct our library Θ to include over 40 trigonometric and polynomial terms. The most challenging part of this example is building a library with the necessary terms, without it growing too large. The library cannot be too extensive, or else the matrix Θ becomes ill conditioned, making it sensitive to noise. To reduce the library size, we use one piece of expert  knowledge: the trigonometric terms should only consist of φ 1 and φ 2 , the rotational angles of the pendula.
The candidate functions are chosen as a combination of state derivatives and trigonometric functions. Figure 6 shows that SINDy-PI can identify the equations of motion for low noise. For larger noise, SINDy-PI misidentifies the dynamics, although it still has short term prediction ability.

(b) Single pendulum on a cart
We now apply SINDy-PI to identify a fractional ODE problem with control input, given by the single pendulum on a cart in figure 7. SINDy has already been extended to include control inputs [48], although the original formulation doesn't accommodate rational functions.
The dynamics are derived from the Euler-Lagrange equations. All system parameters except for gravity are chosen to be 1, as summarized in appendix D; the governing equations and SINDy-PI models are shown in appendix E. The cart position is denoted by s. The state vector is given by The library is constructed using a combination of trigonometric and polynomial terms. Around 50 different basis functions are used for the library, and around 10 terms are tested as candidate functions. We add Gaussian noise to all system states. We then test the SINDy-PI model on a testing initial condition x test = [π , 0, 0, 0] T with control input F test = −0.5 + 0.2 sin (t) + 0.3 sin (2t) for time t = 0 to t = 2. Figure 7 shows the resulting SINDy-PI models. The structure of the model is correctly identified up to a noise magnitude of 0.01. Beyond this noise level, the SINDy-PI identified model only has short term prediction ability.

(c) Simplified model of the Belousov-Zhabotinsky reaction
We now apply SINDy-PI to a challenging PDE with rational nonlinearities, a simplified model of the Belousov-Zhabotinsky (BZ) reaction. The simplified BZ reaction model is given by [73] and where x, z, s and u are dimensionless variables and = (∂ 2 /∂x 2 s ) + (∂ 2 /∂y 2 s ) denotes the Laplacian operator.
The strong coupling dynamics and implicit behaviour in (4.2a) make the data-driven discovery of the simplified BZ reaction challenging when using implicit-SINDy and PDE-FIND. However, SINDy-PI correctly identifies the simplified dynamics of the BZ-Reaction, as shown in figure 8. To generate the simplified BZ reaction data, we use a spectral method [74,75] with time horizon T = 1 and time step of dt = 0.001. We use n = 128 discretization points with spatial domain ranging from −10 to 10. The initial condition is chosen to be a mixture of Gaussian functions. Eighty per cent of the data is used for training, and the remaining 20% is used for model selection. The right-hand side library is normalized during the sparse regression process. A range of sparsity parameters λ are tested from 0.1 to 1, with increments of 0.1 The other system parameters in (4.2) are given in appendix D and the SINDy-PI model is given in appendix G.   possible to extract the physical laws directly. These equations contain important information about the system and may be more concise, useful and straightforward than the underlying ODE or PDE. For example, given a Lagrangian, we can derive the equations of motion. The most difficult aspect of using SINDy-PI to identify a physical law is how to build the library. Conservation laws may contain higher-order derivatives, such asẍ. To include all possible terms, the library may become exceedingly large. The library size will also increase if the system has many states. Large libraries make the sparse regression sensitive to noise. Thus, extracting the physical law from data using SINDy-PI is still challenging due to the lack of constraints when constructing the library function. We only show one example in our paper to demonstrate that it is possible to achieve this using SINDy-PI, but further work is required to reduce the library size so that the sparse regression is robust.
As an example, we consider the double pendulum shown in figure 9, with the system parameters given in appendix D. In this case, we also account for the friction in the pendulum joint, with friction constants of k 1 = 7.2484 × 10 −4 and k 2 = 1.6522 × 10 −4 for the pendulum arms, respectively. In this case, we extract the Lagrangian of the double pendulum [72] using SINDy-PI. To extract this Lagrangian, we simulate the system with initial condition x train = [π − 0.6, π − 0.4, 0, 0] T from t = 0 to t = 15 with time step dt = 0.001. The resulting model is shown in figure 9.

Conclusion and future work
In this paper, we develop SINDy-PI (parallel,implicit), a robust variant of the SINDy algorithm to identify implicit dynamics and rational nonlinearities. SINDy-PI overcomes the sensitivity of the previous implicit-SINDy approach, which is based on a null-space calculation, making it highly sensitive to noise. Instead, we introduce both parallel and constrained optimizations to test candidate terms in the dynamics, making the new SINDy-PI algorithm as robust as the original SINDy algorithm. We also extend the algorithm to incorporate external forcing and actuation, making it more applicable to real-world systems. We demonstrate this approach on several challenging systems with implicit and rational dynamics, including ODEs, actuated systems and PDEs. In particular, we discover the implicit dynamics for a simplified model for the BZ chemical reaction PDE, the double pendulum mechanical system and the yeast glycolysis model, which have all been challenging test cases for advanced identification techniques. Throughout these examples, we demonstrate considerable noise robustness and reductions to the data required, over the previous implicit-SINDy algorithm.
Despite the advances outlined here, there are still many important avenues of future work. One limitation of this approach, and of SINDy in general, is in the design of the library of candidate functions. The goal is a descriptive library, but the library size grows rapidly, which in turn makes the sparse regression ill-conditioned; other issues effecting robustness are discussed in appendix I. Recently, tensor approaches have been introduced to alleviate this issue, making libraries both descriptive and tractable [45], and this is a promising approach that may be incorporated in SINDy-PI as well. More generally, automatic library generation, guided by expert knowledge, is an important topic. Other research directions will involve parametrizing elements of the library, so that the algorithm simultaneously identifies the model structure and the parameters of the sparsely selected terms. Recent unified optimization frameworks, such as SR3 [59,60], may make this possible. Model selection is another key area that will required focused attention. Balancing accuracy on test data, sparsity of the model, and the potential for overfitting are all serious concerns. The sparse regression and optimization may also be improved for better noise robustness. Finally, modifying SINDy-PI to incorporate prior physical knowledge and to only model the discrepancy with an existing model [76] will be the focus of ongoing work.
Data accessibility. Our code and data could be seen in the following link: https://github.com/dynamicslab/ SINDy-PI.
Authors' contributions. K.K. performed the research, designed the code and generated data and figures. All authors analysed the data, refined the figures and wrote the paper. All authors gave final approval for publication and agree to be held accountable for the work performed therein. To compare the performance of SINDy-PI and implicit-SINDy for noisy data, we must define an evaluation criteria. We compare the performance of the best model generated by each method that has the lowest prediction error on the test data, selected according to equation (3.4). To compare the models generated by the two methods with the ground truth model, we use the concept of model discrepancy [76][77][78] and set prediction accuracy, structural accuracy and parameter accuracy as our performance criteria. A good prediction error does not guarantee the model has good structural accuracy and parameter accuracy, and vice versa, motivating multiple performance criteria.

(b) Numerical experiments
We use the Michaelis-Menten kinetics, given by equation (3.9), to compare the performance of SINDy-PI and implicit-SINDy. We performed our numerical experiments as follows: Step 1: Randomly generate 2400 different initial conditions of different magnitudes ranging from 0 to 12.5. Simulate those initial conditions using a fourth-order Runge-Kutta method with time step dt = 0.1 and time horizon T = 5. The testing data is generated using 600 random initial conditions using the same method as the training data.
Step 2: Add Gaussian noise to the training and testing data. 23 different Gaussian noise levels with magnitudes ranging from 10 −7 to 5 × 10 −1 are used. For each noise level, 30 different random noise realizations are generated, resulting in 30 different noisy datasets for each noise level.
Step 3: Compute the derivative of the noisy data. We investigate several approaches, including finite-difference and total-variation regularized difference (TVRegDiff) [56] derivatives. In all cases, SINDy-PI is several orders of magnitude more robust to noise than implicit-SINDy, and only the result of using TVRegDiff is shown in this paper. TVRegDiff generates more accurate derivatives, but also requires hyperparameter tuning and causes aliasing, so we trim the ends of the time series generated by each initial condition (first and last 30%). It is possible to add Gaussian noise to the clean derivative data to investigate robustness, although this is less relevant for real-world scenarios, where only noisy state data is available.
Step 4: Train SINDy-PI and implicit-SINDy models on noisy training data. For each noise level, we sweep through 68 different sparsity parameters λ for SINDy-PI, from 0.01 to 5. The λ is varied by a factor of 2 [52] to calculate the null space in the implicit-SINDy method.
The library for implicit-SINDy and SINDy-PI is Step 5: Due to the various parameter values, we use model selection to choose a model. We use the test data with the same noise magnitude to perform the model selection process. The ratio of training data and testing data is 8 : 2. missing from the model), and parameter error as our model performance evaluation criteria. We average the performance over 30 different noise realizations for each noise level. We then plot the distribution of structure error in figure 3.
Many parameters affect the performance of these methods: the length of training data, prediction steps to calculate prediction error, the initial conditions for training data, choice of the library, and the derivative computation. We have attempted to carefully optimize each method, although an exhaustive parameter sweep is beyond the scope of the present work. However, in all cases SINDy-PI outperforms implicit-SINDy.

Appendix B. Data usage of SINDy-PI and implicit-SINDy
Section 3d investigates the data usage of SINDy-PI and implicit-SINDy on the yeast glycolysis model in equation (3.10). The parameters of this problems are given in table 1. The data usage comparison is performed by the following steps: Step 1: Generate training data by simulating equation (3.10) with parameters in table 1 and a time step of dt = 0.1, with time horizon T = 5. We simulate the system using 900 random initial conditions with magnitude ranging from 0 to 3.
Step 2: Shuffle the training data and select j per cent of the entire training dataset at random to train the SINDy-PI and implicit-SINDy models. These training data are sampled from all trajectories, and they are not necessarily consecutive in time. No noise is added since we only care about the effect of the data length in this case. The sparsity parameter λ is fixed for both algorithms (different values); this value is selected for a single percentage j where both methods fail to identify the correct model, and we sweep through λ.
Step 3: Run the numerical experiment for 20 times for each data length and calculate the percentage of times the two algorithms yield the correct structure of the equation (3.10f).
The final comparison is shown in figure 4. Data usage requirements for other state equations are given in table 2; figure 4 shows results for the hardest equation to identify. The other equations require less data. Normalizing the SINDy-PI library improves data learning rates as well.

Appendix C. SINDy-PI and PDE-FIND on rational PDE problem
In §3e, we compared the performance of SINDy-PI and PDE-FIND on a modified KdV equation. The simulation data is obtained using a spectral method [74] with a time step of dt = 0.01 and time horizon T = 20, spatial domain L = −25 to 25, and n = 128 spatial discretization points. The library of PDE-FIND is chosen to be while the left-hand side library is chosen to be For both SINDy-PI and PDE-FIND, we used 100 different values for the sparsity parameter λ ranging from 0.1 to 10 with step size 0.1. We use 80% of the simulation data for training and 20% for testing and model selection. We calculate the normalized prediction error for all models on state u t and the model with minimum prediction error is selected as the final model.

Appendix D. Parameter values for simulations
The parameters for the double pendulum simulation in §4a are given in

Appendix E. SINDy-PI models for the single pendulum on a cart
The Lagrangian for the single pendulum on a cart with an input force on the cart is: where m is the mass at the end of the pendulum arm, M is the mass of the cart, L is the length of the pendulum arm, s is the position of the cart and φ is the pendulum angle. We do not consider damping in this case. Using the numeric values m = M = L = 1 and g = −9.81 this simplifies to The Euler-Lagrange equation of the system are d dt where F is the force applied to the pendulum cart. It is possible to isolateφ ands:   It is possible to write this as a system of four coupled first-order equations With the numerical values shown in table 5, this becomes Parameters identified by SINDy-PI under different noise magnitudes are presented in tables 6 and 7.
If we use the values in table 3  With no noise, SINDy-PI discovers the correct equations. When we add random noise with magnitude of 0.005, SINDy-PI discovers the following    of the polynomial terms in the model, etc. The relative impact of all of these factors varies with the system we are studying.
In general, training data that explores more of the phase space will result in a more robust model discovery process. Generally speaking, data that results in a better-conditioned Θ matrix will provide more robust results. Several studies have explored strategies to improve the robustness, with Wu and Xiu suggesting the use of a large ensemble of initial conditions [43]. Exciting transients is also important [28].
Another key factor that affects the condition number of Θ is the number of library elements, which is determined by the order for a polynomial library. Including higher-order terms increases the condition number, making it more difficult to accurately disambiguate which nonlinear term is responsible for the observed behaviour. The library size scales exponentially with the maximum polynomial order.
Noise robustness is also affected by the model parameters. Smaller parameters in Ξ are more likely to be removed during thresholding, making the procedure less robust to noise. To illustrate this, we change the parameter K m in (3.9) and test the maximum noise SINDy-PI can handle. We set K m equal to 0.01, 0.1, 1 and 10 separately. We generate the training data with 120 random initial conditions ranging from 0 to 10. Each initial condition is simulated until T = 5 with dt = 0.01, and the same magnitude of Gaussian noise is added. TVRegDiff is used to calculate the derivative, and the first and last 30% of data is discarded due to aliasing effect. The testing data is generated using the same process and the ratio of training and testing data is 1 : 1. The model and maximum magnitude of noise allowed for each values of K m is summarized in table 8. These results suggest that as K m increases, the maximum noise SINDy-PI can handle also increases.