Bayesian Spline Learning for Equation Discovery of Nonlinear Dynamics with Quantified Uncertainty
Authors: Luning Sun
Daniel Zhengyu Huang
Hao Sun
Jian-Xun Wang
Venue: NeurIPS 2022 Accept (2022)
Keywords: Equation Discovery
Deep Learning
Spline Learning
Bayesian Method

Abstract: Nonlinear dynamics are ubiquitous in science and engineering applications, but the physics of most complex systems is far from being fully understood. Discovering interpretable governing equations from measurement data can help us understand and predict the behavior of complex dynamic systems. Although extensive work has recently been done in this field, robustly distilling explicit model forms from very sparse data with considerable noise remains intractable. Moreover, quantifying and propagating the uncertainty of the identified system from noisy data is challenging, and relevant literature is still limited. To bridge this gap, we develop a novel Bayesian spline learning framework to identify parsimonious governing equations of nonlinear (spatio)temporal dynamics from sparse, noisy data with quantified uncertainty. The proposed method utilizes spline basis to handle the data scarcity and measurement noise, upon which a group of derivatives can be accurately computed to form a library of candidate model terms. The equation residuals are used to inform the spline learning in a Bayesian manner, where approximate Bayesian uncertainty calibration techniques are employed to approximate posterior distributions of the trainable parameters. To promote the sparsity, an iterative sequential-threshold Bayesian learning approach is developed, using the alternative direction optimization strategy to systematically approximate L0 sparsity constraints. The proposed algorithm is evaluated on multiple nonlinear dynamical systems governed by canonical ordinary and partial differential equations, and the merit/superiority of the proposed method is demonstrated by comparison with state-of-the-art methods.

Bayesian Spline Learning for Equation Discovery of Nonlinear Dynamics with Quantified Uncertainty
Luning Sun' University of Notre Dame 1sun78nd. edu
Daniel Zhengyu Huang" California Institute of Technology
dzhuangOcaltech edu
Hao Sun
Jian-Xun Wang
Renmin University of China haosungruc. edu. cn
University of Notre Dame jwang33end. edu
Abstract
Nonlinear dynamics are ubiquitous in science and engineering applications, but the physics of most complex systems is far from being fully understood. Discovering interpretable governing equations from measurement data can help us understand and predict the behavior of complex dynamic systems. Although extensive work has recently been done in this field, robustly distilling explicit model forms from very sparse data with considerable noise remains intractable. Moreover, quanti- fying and propagating the uncertainty of the identified system from noisy data is challenging, and relevant literature is still limited. To bridge this gap. we develop novel Bayesian spline learning framework to identify parsimonious governing equations of nonlinear (tpppio)lemmoral dynamies from sparse, noisy data with quantified uncertainty. The proposed method utilizes spline basis to handle the data scarcity and measurement noise, upon which group of derivatives can be accurately computed to form a library of candidate model terms. The equation residuals are used to inform the spline learning a Bayesian manner, where approximate Bayesian uncertainty calibration techniques are employed to approxi- mate posterior distributions of the trainable parameters. To promote the sparsity, iterative sequential -threshol-thh Bayesian learning approach is developed, using the alternative direction optimization strategy to systematically approximate LO sparsity constraints. The proposed algorithm evaluated multiple nonlinear dy- namical systems governed by canonical ordinary and partial differential equations, and the merit/superiority of the proposed method demonstrated by comparison with state-of-the-a methods.
1 Introduction
In the realm of science and engineering, dynamical systems are ubiquitous. However, in actual circumstances, the governing equations behind complicated dynamics may not be completely derstood. preventing researchers from developing first- trist-rincipped models. On the other hand, the ever-increasing data availability opens up new avenues for scientists to identify predictive models from enormous observation data, process known as system identification (SI). Recent advances in deep learning have prompted the rapid development of powerful SI models for sigh-dimessional problems using deep neural networks (DNNs). Many DNN-based SI models have been proposed to
"qqual contribution Corresponning author
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

learn differential operators for complex (spatio physics from data and shown good potential in terms of data reproduction and state prediction [TH However, deep learning models usually lack interpretability and are difficult to comprehend. Furthermore, it is questionable in terms of generalizability when compared to first-principle models, as such black-box DNN models provide less insight into the underlying processes.
Instead of identifying a black-box model, we focus on extracting analytical equation forms from data, which has higher interpretability and has the potential to advance our knowledge of unknown physics. The sparse identification of non-linear dynamic (SINDy) algorithm 4 is an excellent development along this path. The central idea is to use sparse linear regression to uncover parsimonious governing equations from a dictionary of basis functions constructed by data, where the sparsity is promoted by pruming out redundant terms based on certain specified thresholds [5]. SINDy. although showing great promise, faces several challenges: (1) it heavily relies on high-quality data to extract derivative information, which is typically based on finite difference (FD) methods, making it impossible to handle incomplete, scarce, or noisy data; (2) it is formulated in a deterministic fashion and cannot account for the uncertainties from multiple sources, which is critical for real-world applications where data is frequently corrupted and can be very sparse, In the past few years, the SINDy framework has been further improved in various aspects to address these challenges, e.g., enhancing the library [6] or using deep learning for denoising and derivative computation by fitting the noisy data in a decoupled 77-9 or coupled manner [10 For uncertainty quantification, the dictionary- -bssdd equation discovery algorithms have been recently extended to Bayesian settings [7], based on
the idea of sparse Bayesian learning pioncered by Tipping and co-workers [8822
Despite recent progress and extensive work in this field, reliably distilling explicit equation forms from very sparse data with significant noise remains an unsolved challenge. There are still significant gaps in handling data scarcity and noise, quantifying and reducing multi-source uncertainties, and promoting sparsity, which most existing equation discovery techniques struggle to simultaneously address. To this end, we propose a novel Bayesian Spline Learning (BSL) approach to identify parsi- monious ordinal/partial differential equations (ODEs/PDEs) from sparse and noisy measurements; meanwhile, the associated uncertainties are quantified. To deal with data scarcity and measurement noise, the proposed BSL uses spline basis, on which collection of derivatives can be reliably computed for the library construction. The posterior distributions of spline-bascd model parameters approximated by a stochastic gradient descent (SGD) trajectory-based training scheme, where the first moment of SGD iterations is computed by stochastic weight averaging approach (23]. The proposed BSL approach is effective in two aspects. On the one hand, the spline-based representation can help to interpolate locally the spatiotemporal field and perform differentiation analytically. As a result, it considerably enhances learning efficiency in cases with sparse and noisy data. The posterior distributions of spline, library, and equation coefficients are estimated simultaneously, without adding too much overhead the training process. Besides the measurement uncertainty, the model-form uncertainty is also obtained in the proposed method, which can be used for downstream Bayesian data assimilation, where online data is assimilated to improve the predictability of the identified system for chaotic scenarios
2 Related Work
The data in real-world circumstances is frequently sparse in spatial/temporal domains and may contain considerable noise, posing significant challenges to SINDy or its variants [46)] 16, (7] Deep learning (DL) has been leveraged as a superb interpolator for concurrently generating metadata and smoothing high-frequency noise (28). effectively improving the performance in identifying equations from imperfect measurements [0.14. For example, automatic differentiation (AD) is effective for computing derivatives analytically from sparse, noisy data using a point-wise multi- layer perceptron (MLP) [3] 29 30). However, it is difficult to impose locality constraints in the point-wise formulation. Wandel et al. [31] demonstrated the merit of employing a spline basis for analytically calculating derivatives while enforeing the locality of spatiotemporal fields. Owing to the superiority of spline-based differentiation compared to numerical discretization, physics-informed spline networks remarkably outperform point-wise physics-informat neural networks (PINN) in the
context of solving PDEs 31 and discovering ODEs [14
'The code will be available at http::///ghhhcoooo///ugn///////n

---

However, these models fail to simultancously deal with ODE and PDE systems when the data is sparse and substantially corrupted. For the PDE datasets, for example, the PINN-based sparse regression (PINN-SR) 13 works admirably, but it fails to converge on ODE datasets. On the other hand, algorithms such as physics-informed spline learning (PiSI.) [4] that performs very well for ODE discovery cannot handle PDE problems due to the limitation of the B-spline hasis adopted. More importantly, all these SOTA sparse learning algorithms are formulated in deterministic settings and uncertainties introduced from data/library imperfection cannot be quantified.
For inverse problems like equation discovery, it natural to use Bayesian framework quantify and analyze the prediction uncertainty. To enable the posterior computation for high-dimension trainable parameters that are intractable by traditional Bayesian inference, people have resorted various approximation strategies, e.g., variational inference (4, Monte Carlo dropout [35], Bayes by backprop (36), Laplace approximation [37], and deep ensemble approaches [38. Although these techniques have had great success, training may still be challenging, and costs will rise dramatically as the problems become more sophisticated. Alternatively, we chose to employ an SGD trajectory- -based algorithm, Stochastic Weight Averaging- Gaasssinn (SWAG) [23] where the information in the SGD trajectory is exploited to approximate the posterior. As demonstrated empirically, SWAC is well scalable to high-dimensional problems and can accurately estimate uncertainty across many different Bayesian learning tasks 23 39].
Data assimilation (DA) has been widely used in numerical weather prediction (NWP) by fusing online sensing data into predictive model for nonlinear dynamics forecast. People have recently integrate deep learning into DA to improve online prediction performance [20. However, the predictive model in DA is assumed be known a priori, and uncertainties from the model and observations, which are required for DA, are usually hard to obtain. In our framework, the predictive model can be unknown a priori and will be identified to assimilate additional data for online forecasting. Moreover, instead of arbitrarily guessing the observation and model-form errors, this information can be learned in the proposed BSL framework, as the data and model-form uncertainties are quantified. The UQ capability of BSL naturally integrates the equation discovery with Bayesian DA techniques.
The main contributions of this work are three-fold: (1) we extended spline learning for sparse equation discovery of spatiotemporal physics governed by ODEs or PDEs; (2) we developed sparsity- promoting Bayesian learning for UQ tasks; and (3) the proposed BSL framework can seamlessly integrate with Bayesian data assimilation techniques to improve online dynamics forecasting.
3 Methodology
Let consider dynamical system, which is governed by parametric ODE/DDE system in the following general form
du dt = F(u).
(I)
where u denotes the state vector, for ODE systems, u = [us (t), uz(t), (t())) E Rd depends only on time I, for PDE systems, u(x, () = [ui (x, l), uz (x, (),-. ... ua ua(x, f)j E Rd depends on both time and space x, and F Rd + R2d represents unknown nonlinear functions. The states are observed at discrete times h= and at spatial locations {x, } for PDE systems. Let û denote the noisy observation vector, the observation set is U (ū(4), u(t2). ... ū({,)) R and Ü= [ū(x), f1), ú(x;,tt): ...,..... ü(x,, fn) E (mxs)xdd for ODE/PDE systems, respectively. Our goal here is to explicitly discover the parsimonious form of FF-) from library of candidate basis functions and quantify the associated uncertainty given noisy observation data.
3.1 Overview
Given the dataset U, the SI problem can be solved by sparse regression techniques with predefined library & (u) of m basis functions,
(u) (6 (u), d2) 2m (u)] E R.
(2)
where di : Rdd_ + R,ISis m denotes basis function, which, for instance, can be the polynomial or trigonometric function. Hence, the matrix of library terms evaluated at observed states is defined
3

---

1-D basis
2-D basis
Spline Learning
Bayesian ADO Training
UQ
L priors (Ú N B (Ú N0
ADO Itr = 1 Stage1 + Stage2
W*.B.P
+NN0 WW NN.0 OW
du, du
dt
dt
18,4; my W/N2
O U
EnKF
Stage
Stage
W
nodd
Prediction
Approximate Posteriors
m S
Epochs
x priors+ -NN 1N.0 (NN
Figure 1: Overview of method. Black box: overall work flow. Blue box: demo for spline basis with derivatives ID and 2D cases. Red box: a sketch for the Bayesian ADO training process. Green box: the data assimilation enhanced uncertainty quantification.
as,
(((()))),(((())))) ((Ú) (ū(tt))) ((ū(tt))* .M ((((tn)))
ODE:
R=xm
(3)
PDE:
((U)
o ti)x" " ((((x,,)))" , {((xx,tt))
e R (nxa)xm
Given the library, the overview of the BSL framework is shown in Fig. Firstly, a spline- based model constructed represent state variables and their derivatives by denoising the observation data, which is illustrated in the blue box and discussed Sec 3.2 Then the spline-based model and sparse regression trained simultaneously by alternating direction optimization (ADO), shown in the red box. Specifically, single ADO iteration contains two sub-processes, where sub-process I trains the spline-based model using log posterior loss and passes the trainable parameters and the noise estimations to sub-process II, which uses a Bayesian SINDy-like method to prune out redundant terms the library defined in Eq. 3] and then updates the number of relevant terms the training loss for sub-process I. After several ADO iterations, the parsimonious form of the governing equations and the posterior distribution for the coefficients will be estimated. More details about the ADO iterations in the Bayesian framework are shown in Sec. 3.3 and Appendix A.3. Finally, with the estimated posterior, the predictive uncertainty can be quantified by evaluating the identified system with an ensemble of parameters. To further improve the prediction capability, especially for chaotic systems, we propose to leverage data assimilation techniques, which is shown in the green box and discussed in Sec 3.4 and Appendix A.5.
spline-based
3.2 Spline-based learning
Several previous works have already shown the potential spline-based learning and demonstrated the advantages compared with the classical DI. structures, (MIP, CNN) 14 31 43 Therefore, we use this structure to smooth the solution fields based on the noisy measurement and then identify the underlying PDE/ODE and calculate the corresponding derivatives from the rpoosserccontttuctted

---

fields. The B-spline curves, based on the De Boor's algorithm, are defined in recursive way as:
No,((1) [1, Tf Sts s++1-
00 otherwise.
(4)
(t) = -- 1t)- Tatk+1 Ts+1
Toth+1 -t
No+1,A-1(t),
Tatk
where Tu is the location of knots, & is the degree of polynomial. When k: 3, it is the well-used Cubic-B Spline curve. With the defined basis, the spline interpolation can be write as y(t) = *x=0 2 ((,, And the number of control points (trainable weights) O E B is chosen empirically. can be proved that the derivative of p order B-spline basis is function of p 1 order B-spline, written
k
k
N (a)((1) di
Na,E-1(1) Ts+h+1 Ts+1 N541IE-1())-
(5)
Tsth
The proof of Eq. 5 is attached in Appendix A.I. The higher-order derivative be calculated by recursively using Eq. [ The analytical derivatives spline basis are very beneficial for PDE discovery tasks since it always involves constructing library terms containing high order derivatives. With a proper order the first derivatives are accurate, and there is no error introduced during the derivation, opposed to using the numerical methods to approximate derivatives. The spline function and its derivatives are defined in one-oonnenn nne-dimensiona scenario. is straightforward to extend n-dimension by direct using tensor-product. For example, for two-dimensional problem with spatial-temporal fields, the basis can be defined as
Nk (fz, 22) = Ne, (.)))),,a, $1.82
(6)
Here the two-dimensional basis is denoted by different style write it compactly, Similarly, the partial derivative for two-dimensional basis is defined as:
a 022) Nhu. 91.92 (t,,
d(er) No (t,) dtlas)
((42) Na2. (zz))
dr (12)
With the definition of sparse system identification and spline reconstruction, the whole spline learning can be stated follows: given noisy measurement data U, find the best sets of weights O and W so that data fitting loss and the weakly physies-informed loss can be minimized under sparsity constraints, as shown in Eq. 8
(0. W] arg min 0'.W" Nn
II - + No ((NO'' N'0'l A/|W'
(8)
Here, Nm, and N₂ denote the spline basis matrices evaluated at measurement and collocation locations. Moreover, Nmm and N, are numbers of measurement data and collocation points. Furthermore, the alternating direction optimization (ADO) shown in previous works [3] 14 can be adopted to minimize the loss function efficiently, and the details are attached in Appendix A.3.
3.3 Sparse system identification in Bayesian formulation
Bayesian methods provide natural probabilistic representation of uncertainty, which is crucial for model predictions. System identification from noisy and sparse measurements generally contains two types of errors (similar to the hidden Markov model): and the smoothed data is reconstructed through spline-based learning, shown in Eq. 9
(1) Observation error, where the data noisy,
N 0' +11,
where e represents the observation error. (2) Evolution error or model form error, since the discovered system cannot be exact due to library imperfection and needs to be reformulated
i = {u))
+ E2
(10)
where 62 represents model form error. And Eq. 10) is evaluated collocation points.
In this work, these error terms are modeled as zero-mean multivariate Gaussian random variables:
E (O, B) and E2 ~ (O. And we further assume that these error covariance matrices B, P
5

---

are diagonal matrix with diagonal terms {ba} and {pk} with 1 I k S d and they are learn-abl parameters during training.
According to Bayes' rule, the posterior can be written as:
p(0, W, W, P,U Ù) & p(O, W, B, W.B.P))(U,UW.OOB
(1I)
Here Ü N,, denotes the derivative estimation on all the collocation points. The prior further decomposed by
P W. O, B) x ((Waapt()pp)))))))))))
(12)
Currently, we specify the prior for the linear coefficient matrix p(W[a)= =
mean Gaussian distribution
with the hyper prior as a Gamma distribution p(a:) = Gamma (a: (0o. bo). Similarly, also define the prior for the spline trainable parameters zero mean Gaussian distribution with the hyper prior as another Gamma distribution p(0] = ( 0, 3) and p(3) (3)am b1). To account for the data uncertainty and process uncertainty, the di- agonal covariance matrices B and P are also set as learn-able during the training. And hence improper uniform priors are used for B and P. The likelihood consists of two parts, as:
P(ú U/W U,W.O,B.P) = (/UO,BB)xUUWWPOO.)
(13)
where
p( B) x expl - -0 Nm - (U- N Nm
(14)
P(Ú & expl- (N.0
- ((N.O)W) -NNN (N.O-9(N.0)W))
(15)
Traditional Bayesian sampling approaches (i.e., Markov chain Monte Carlo methods) can be usu- ally intractable and expensive, especially when the parameter dimensionality is high. Therefore, researchers tend to use alternative approximation approaches. Current work uses a stochastic gradient descent (SGD) trajectory-based approach, Stochastic Weight Averaging Gaussian (SWAG) [23) algorithm to approximately sample from the posterior distribution. This method approximates the posterior by collecting the parameters near the loss plateau after a sufficient number of the training steps. To further reduce the inference cost, we construct subspace by finding the PCA components of the SWAG trajectories 44 and then draw samples in the subspace instead. terms of the loss function in the probabilistic model, we chose to maximize the log form of the posterior density Eq [T] also by leveraging the ADO algorithm. The sub-prou sub-process I minimizes the log density function defined by Eq. 13 13 and the sub-process II adopts the Bayesian variants of SINDy algorithm modified from The likelihood function in sub-process II has the same form as in Eq. 15 but with a different
117 19.
sparsity promoting prior p( (WIA) = II (W")),,a'")) where A = [a ag. o'm In single
j=1
ADO iteration, the sub-process provides the updated 0 and P to sub-process II for constructing the library terms. While the sub-process II shrinks the library terms and passes updated library/weight structure W back to sub-process I. The whole training requires multiple ADO iterations before it reaches the final balance, where no more terms will be pruned out sub-process II. The detailed ADO algorithm is listed in Appendix A.3. After obtaining the approximated posterior distribution, predicting uncertainty can be estimated by marginalizing out the model parameters.
3.4 Data assimilation for enhanced predictability
For chaotic system, one notorious problem is that a slight perturb in any model parameters can significantly influence the prediction. For example, in numerical weather prediction (NWP) tasks where the researcher always needs to predict the behavior of chaotic weather systems, various data assimilation (DA) techniques have been developed to assimilate the available data and the known equations. Kalman filter and its variants are successful mathematical tools for data assimilation. Our Bayesian formulation, providing the observation error covariance matrix B and evolution error covariance matrix P, can be naturally incorporated with the Kalman filter frameworks to improve the predictability of chaotic systems. Current works choose the ensemble Kalman filter for the task, and more background can be found in Appendix A.3.

---

True COEFF
Predicted
-1.0
Equation Terms
Ttue Learned
True Learned
Coy (1)
C(x
x2)
True COFFF
Predicted Mean
Equation Terms
Predicted
Predicted
cooop
RO
7.50
7.95
8.00
8.05
-1. 1..02 11.000 - 998 -0.996
C(F())
C(X
Figure 2: The discovery results for ODE systems; the red box shows results for the Van der Pol (VdP) system, and the blue box shows the result for Lorenz 96 system. The layout inside each box follows the rules below. Upper left sub-fig: discovered mean for the relevant library terms; Lower left sub-figs: selected posterior distribution for the identified distributions: Right sub-figs: ensemble prediction plots for UQ. For VdP system, the governing equation de
dx = y and
dy - 0.5x² 0.5y. For the Lorenz96 system, the compact form of the governing equation is de
dx: (X,+1 X1 X, + F with periodic boundary conditions.
Experiment and Result
In this section, we first show the equation discovery and uncertainty quantification results for nonlinear ODE systems. We also show that incorporating DA techniques improve the predicting ability for chaotic systems. Finally, we present the PDE discovery results for several canonical PDE systems.
The first pedagogical ODE example is the Van der Pol oscillator, which defined as de di dt = =(1 - x²)y with = 0.5. The data is corrupted with 5% noise, and the library be
and dy
consists of polynomials of state variables up 3rd order. The parsimonious model structure correctly identified, and the result is shown in the upper red box of Fig. [ There are five sub-figures, which depict the library discovery, coefficient distribution, and the forward propagated UQ results. Specifically, the single upper left figure shows the mean of discovered coefficients (blue lines -) and the truth equation coefficients (black stars *), where the horizontal axis represents term indices. For example, the Van der Pol system has four different terms, and the x-axis ranges from to 4. The two sub-figures the lower left part show probability density distributions (PDF) of two identified coefficients, where truth equation coefficients (red lines -) fall within the confidence interval with high probability. The right sub-figures show the propagated ensemble results (blue lines -) based on the discovered equations, measurement data (green circles v), and the true state trajectories of æ and y (red dashed lines -)) The result clearly shows that the ensemble predictions can cover the data, and the uncertainty range fluctuates around the truth state values.
7

---

72
Figure 3: The discovery results for PDE systems; the red box shows results for the advection system, and the blue box shows the result for the Burgers' system. The layout inside each box follows the rules below. Leftmost sub-fig: true contour plot; Middle two sub-figs: the spatial results at different time t; Rightmost sub-figs: the temporal result at a fixed point æ.
The second ODE example is a chaotic system, Lorenz 96. It is a simplified mathematical model for atmospheric convection, defined as: dx. d! (X1+1- X X; F i= 1,2, 2....... ... with periodic boundary conditions X-1 X Xn-1, Xo = X, and X X1. We chose current case and F = 8 the forcing terms. The measurement states variables are corrupted with 10% Gaussian noise, and the library consists of polynomials of state variables to 3rd order. The discovery result for these systems is shown in the lower blue box in Fig. 2 following the same layout the Van der Pol systems. The upper left figure indicates that the sparsity structure can be identified, and the horizontal axis marks 24 parsimonious terms (4 terms for each state, states in total) out of 84 library terms. The detailed PDF plots show that posterior distributions still cover truth values with high probability. Although the discovery result is quite accurate, the forward simulations of identified Lorenz 96 system would induce large phase difference compared with truth trajectories due the chaotic nature of the underlying system (in Appendix A.5). To alleviate the chaotic behavior and improve the predictability, we seamlessly coupled the ensemble Kalman filter (EnKF), classical DA technique, to assimilate noisy measurements with the identified system. The observation error matrix B and evolution error matrix P required for the En scheme can be directly passed from the Bayesian spline learning framework. The right subplots show the DA results for the 3 out of state variables to save space. The horizontal axis marks the evolving time, and the vertical axis represents the state variables. In each subplot, the green shaded region indicates the time interval with available noisy measurement data (marked by green circles (o)). And the red solid-line (-) is the true L96 states. Finally, the blue curves (-) are the ensemble predictions. It can be seen that the ensemble predictions inside the region with measurement data almost overlap with the true trajectory. Furthermore, the ensemble uncertainty grows more significantly in the extrapolation region, but the ensemble ranges still fluctuate around the true trajectory. Note that the prediction for the chaotic system with the DA process is much better than directly forward simulating the identified system, where the predictable interval only about second. More discovery results are attached in Appendix A.5.
4 in
Then, we use the proposed BSL framework identify classical PDE systems and evaluate the performance of two of them in the main text. They include advection equation Ou du 0 and
au + Bg Burgers' equation ät red and blue boxes, respectively. Inside each box, the leftmost contour plot shows the true state value in the spatiotemporal field, and the right three sub-figures show the UQ results for different cross-sections from the left contour plot. The lines (-) represent the truth value, blue lines (-) are obtained by ensemble predictions, and green circles (o) represent sparse and noisy measurements. The prediction from the identified PDE system is accurate in the spatiotemporal field, and the ensemble fluctuates around the truth state variables.
0² = O. The qualitative results are shown in Fig. 3 inside the du 0.5 Bar

---

Finally we benchmark our proposed method with several state-of satte-of-thharre the-art discovery algorithms (PINN- SR [13, SINDy [5] and RVM 445 The error metric is defined as:
Discovery
C
rmse
(16)
C la
Diseoveery
o C
Mp
(17)
||C Discovery llu
Discooeery O True llo
MR
(18)
|[Truello
where DDicoveryy are the non-zero mean prediction from the posterior distribution and CTrne are the true coefficients of the governing equations. the method fails to converge or cannot identify correct parsimonious form, we will report the final result as Fail. The errors are scaled by x x10-3 to have clear comparison. Two additional metrics, precision Mp and recall MR. are also defined, where the @ represents element- wise product of vectors and the norm is the non-zero terms in vector. can be seen from Table that our BSL method always performs best, when the noise is significant (> 5%), demonstrating its robustness noise.
The PINN-SR can discover the PDE equation with corrupted data set, but fails to predict accurate time trajectories for ODE systems. This is also reported in relevant paper [46], that the plain MLP structure not satisfactory for predicting time series. The SINDy method can behave much better when large high-quality data exist. However, the SINDy method can easily fail when noise is significant. Our BSL requires fewer parameters and is easier to converge due to the enforcement of locality constraints by the spline basis model. Furthermore, the SINDy and RVM methods require much more data (100%) to train but are still vulnerable to data noise. These benchmark cases show the potential of our BSL for equation discovery tasks. The full table with more benchmark test cases can be found in Tab. 3 in Appendix. We also apply the proposed method on a real-world dataset, as shown in Tab. 9 and Tab. 10 and discuss the effect of smoothing method in Tab. 12 and Tab. 13.
Table 1: ODE and PDE discovery comparison
Name
rmse(0%)
rmse(1%% 1))
rmse (large
Mp
M
Training Cos
Van der Pol Oscillator
BSLrOurs)
0.2
2.82
18.04
133(+-))
PINN-SR
Fail7
Fail
Fail Fail
0.214
0.75
1213s
SINDY
1.0
1.93
0.267
1.0
10s
RVM
2.54
27.46
Lorenz 96
BSL(Ours)
0.269
1.47
13.0
165-((+438)s
PINN-SR
Fail
Fail
Fail
0.5
0.22
10788s
SINDY
0.4
0.64
Fail
0.75
1
10s
RVM
49.7
Advection Equation
BSL(Ours)
0.26
1
1.9
946(+233))
PINN-SR
5.9
4.5
30.4
1
650s
SINDy
2.3
8.2
38.9
1
10s
RVM
0.77
6.76
Fail
0.2
4s
Burgers' Equation
BSL(Ours)
3.62
4.13
6.38
117(+74)s
PINN-SR
10.2
3.3
10.3
1
1
512s
SINDY
0.826
Fail
Fail
0.5
10s
RVM
0.754
Fail
Fail
0.1429
0.5
4s

---

5 Discussion and Limitation
In this work, we developed a novel Bayesian spline learning (BSL) framework for equation discovery from sparse and noisy datasets with quantified uncertainty. The proposed framework significantly improves SI performance and contributes to the existing literature in the following aspects: Firstly, the use of spline basis enables to accurately interpolate solution surfaces and analytically compute derivatives to form the candidate library, outperforming other benchmark methods based on finite difference (FD) auto differentiation (AD), which either suffer from noisy/sparse data or overfitting issues. Secondly, the proposed Bayesian learning formulation notably enhances the robustness for large data noise and sparsity, and meanwhile, quantifies the predictive uncertainty with minimum computational overhead. Moreover, a Bayesian paaritt-promooting ADO iteration strategy is pro- posed to promote sparsity and recover the parsimonious governing equation as well as the posterior distribution of its coefficients. Last but not least, the Bayesian DA is also integrated into the BSL framework improve the online predictability of the chaotic systems, which can potentially benefit real-world tasks such as numerical weather forecasting. The proposed framework is evaluated on discovering multiple canonical ODE and PDE systems, and great superiority has been demonstrated in comparison with state-of-the ff-the--rt art methods.
Admittedly, this work still relies system is largely limited to the functional space determined by the user-specified library. In general, for library-based methods, how to design an inclusive but not unnecessarily large library a priori is important yet very challenging, which may require prior knowledge of the system to be identified and thus limit their applications for systems involving complex governing physics. Moreover, in this work, only a uniform displacement of knots is used for spline representation, and the tensor product of 1-D splines adopted for multi- diltiiddimesionall spline constructions as shown in Tab. These choices are not optimal and have notable limitations for high- -igeedimeesional problems. To tackle this issue, we propose to apply more advanced spline learning techniques to reduce the computational cost and improve the scalability. For instance, using deep learning to optimize knot size has shown great promise for lowering computational costs high-dimensional settings 447 48 while defining continuous spline kernels and expanding basis in subdomains will allow significant reduction of trainable parameters [43]. The improvement of spline learning will be explored in our future work. Lastly, similarly to all data-driven models, the proposed method could have negative societal impact is used abusively, particularly for predictive modeling of high-consequence natural systems, where caution should be taken for decision making.
pre-defined library of candidate terms, and thus the identified
Acknowledgments and Disclosure of F Funding
The authors would like acknowledge the funds from National Science Foundation, United States of America under award numbers CMMI--194330 and OAC-2047127, the Air Force Office of Scientific Research (AFOSR). United States of America under award number FA9550-22-1-008 funds from the College of Engineering at University of Notre Dame in supporting this study
and startup
References
[1] T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, P. Battaglia, Learning mesh-based simulation with graph networks, in: International Conference on Learning Representations, 2020.
[2] L. Lu, P. Jin, Pang. Z. Zhang, G. E. Karniadakis, Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nature Machine Intelligence 3 (3) (2021) 218-21122 229.
[3] X. Han, H. Gao, T. Pfaff, J.-X. Wang, L. Liu, Predicting physics in mesh-reduced space with temporal attention, in: International Conference on Learning Representations, 2022.
'Large noise for different cases: Van der Pol: 5%, Lorenz 96: 10%, Advection: 20%, Burgers: 10%, Burgers'
with source: 20%. Heat: 15%, Poisson: 5%
$ MM, MR are only reported for the largest noise cases
6 All cases are running on a Nvidia 2070 Ti GPU card
Fail means failure in discovery of the parsimonious ODE/PDE forms.
10

---

[4] S. L. Brunton, J. L. Proctor, J. N. Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the national academy of sciences 113 (15) (2016) 3932-3937.
[5] H. Rudy, S. Brunton, L.. Proctor, J. N. Kutz, Data-driven discovery of partial differential equations, Science advances (4) (2017) e e1602614.
[6] H. K. Chu, M. Hayashibe, Discovering interpretable dynamics by sparsity promotion on energy and the lagrangian, IEEE Robotics and Automation Letters 5 (2) (2020) 2154-2160.
[7] H. Rudy, J. N. Kutz, S. L.. Brunton, Deep learning of dynamics and signal noise decomposition with time-stepping constraints, Journal of Computational Physics 396 (2019) 483-506.
[8] J. Wang, X. Xie, J. Shi, W. He, Q. Chen, L. Chen, W. Gu, T. Zhou, Denoising autoencoder, a deep learning algorithm, aids the identification of novel molecular signature of lung adenocarcinoma, Genomics, proteomics & bioinformatics 18 (4) (2020) 468-480.
[9] H. Wu, P. Du, R. Kokate, J.-X. Wang, A semi-analytical solution and Al-based reconstruction algorithms for magnetic particle tracking, Plos one 16 (7) (2021) c0254051.
[10] Z. Long, Y. Lu, B. Dong, PDE-Net 2.0: Learning PDEs from data with a numeric-symbolic hybrid deep network, Journal of Computational Physics 399 (2019) 108925.
[11] S. Kim, P. Y. Lu, S. Mukherjee, M. Gilbert, L. Jing, V. Čeperić. M. Soljačić, Integration of neural network-based symbolic regression in deep learning for scientific discovery, IEEE Transactions on Neural Networks and Learning Systems 32 (9) (2020) 4166-4177.
[12] M. Corbetta, Application of sparse identification of nonlinear dynamics for physics-informed learning, in: 2020 IEEE Aerospace Conference, IEEE, 2020, pp. 1-8.
[13] Z. Chen, Y. Liu, H. Sun, Physics-informen learning of governing equations from scarce data, Nature Communications 12 (2021) 6136.
[14] Sun, Liu, H. Sun, mhysics-informed spline learning for nonlinear dynamics discovery, in: Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence (UCAI-21), 2021, pp. 2054-2061.
[15] S. Zhang, G. Lin, Robust data-driven discovery of governing physical laws with error bars, Pro- ceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 474 (2217) (2018) 20180305.
[16] S. Zhang, G. Lin, Robust data-driven discovery of governing physical laws using a new subsampling- based sparse Bayesian method to tackle four challenges (large noise, outliers, data integration, and extrapolation), arXiv preprint arXiv:1907.0770 (2019).
[17] S. M. Hirsh, D. A. Barajas-Solano, J. N. Kutz, Sparsifying priors for Bayesian uncertainty quantification in model discovery, arXiv preprint arXiv:2107.02107 (2021).
[18] M. F. Tipping, Sparse Bayesian learning and the relevance vector machine, Journal of machine learning research 1 (Jun) (2001) 211-244.
[19] M. E. Tipping. A. C. Faul, et al., Fast marginal likelihood maximisation for sparse Bayesian models., in: AISTATS, 2003.
[20] C. M. Bishop, M. Tipping, Variational relevance vector machines, arXiv preprint arXiv:1301.3883 (2013).
[21] A. C. Faul, M. Tipping, variational approach to robust regression, in: International Conference on Artificial Neural Networks, Springer, 2001, pp. 95-102.
[22] A. C. Faul, M. F. Tipping, Analysis of sparse Bayesian learning, in: Advances in neural information processing systems, 2002, pp. 383-389.

---

[23] W. J. Maddox, P. Izmailov, T. Garipov, D. P. Vetrov, A. G. Wilson, A simple baseline for Bayesian uncertainty in deep learning, Advances Neural Information Processing Systems 32 (2019).
[24] M. Quade, M. Abel, J. Nathan Kutz, S. L. Brunton, Sparse identification of nonlinear dynamics for rapid model recovery, Chaos: An Interdisciplinary Journal of Nonlinear Science 28 (6) (2018) 063116.
[25] K. Champion, B. Lusch, J. N. Kutz, S. L. Brunton, Data-driven discovery of coordinates and governing equations, Proceedings of the National Academy of Sciences 116 (45) (2019) 22445-22451.
[26] K. P. Champion, S. Brunton, J. N. Kutz, Discovery of nonlinear multiscale systems: Sampling strategies and embeddings, SIAM Journal on Applied Dynamical Systems 18 (1) (2019) 312- 333.
[27] L. Zhang, H. Schaeffer, On the convergence of the sindy algorithm, Multiscale Modeling & Simulation 17 (3) (2019) 942-
[28] M. Tancik, P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ra- mamoorthi, J. Barron, R. Ng. Fourier features let networks learn high frequency functions in low dimensional domains, Advances in Neural Information Processing Systems 33 (2020) 7537-7547.
[29] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics- inffommedd neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378 (2019) 686-707.
[30] G.-J. Both, S. Choudhury, P. Sens, R. Kusters, Deepmod: Deep learning for model discovery in noisy data, Journal of Computational Physics 428 (2021) 109985.
[31] N. Wandel, M. Weinmann, M. Neidlin, R. Klein, Spline-pinn: Approaching pdes without data using fast, physics phfoies-inforedd hermite-spline cnns, arXiv preprint arXiv:2109.07149 (2021).
[32] C. Blundell, J. Comebiiee, K. Kavukcuoglu, D. Wierstra, Weight uncertainty in neural network, in: International conference on machine learning, PMLR, 2015, pp. 1613-1622.
[33] Q Liu, D. Wang, Stein variational gradient descent: A general purpose Bayesian inference algorithm, Advances in neural information processing systems 29 (2016).
[34] D. M. Blei, A. Kucukelbir, J. D. McAuliffe, Variational inference: A review for statisticians, Journal of the American statistical Association 112 (518) (2017) 859-877.
[35] Y. Gal, Z. Ghahramani, Dropout as a Bayesian approximation: Representing model uncertainty in deep learning, in: international conference on machine learning, PMLR, 2016, pp. 1050-1059.
[36] S. Ebrahimi, M. Elhosciny, T. Darrell, M. Rohrbach, Uncertainty-guided continual learning with bayesian neural networks, in: International Conference Learning Representations, 2020.
[37] H. Ritter, A. Botev, D. Barber, A scalable laplace approximation for neural networks, in: 6th International Conference on Learning Representations, ICLR 2018 8--ooffeeeeeee Track Proceedings, Vol. 6, International Conference on Representation Learning, 2018
[38] A. G. Wilson, P. Izmailov, Bayesian deep learning and a probabilistic perspective of generaliza- tion, Advances in neural information processing systems 33 (2020) 4697-4708.
[39] M. Abdar, F. Pourpanah. S. Hussain, D. Rezazadegan, L. Liu, M. Ghavamzadeh, P. Fieguth, X. Cao, A. Khosravi, U. R. Acharya, et al., A review of uncertainty quantification in deep learning:
Techniques, applications and challenges, Information Fusion 76 (2021) 243-297.
[40] T. Frerix, D. Kochkov, Smith, D. Cremers, M. Brenner, S. Hoyer, Variational data assimilation with learned inverse observation operator, in: International Conference on Machine Learning, PMLR, 2021, pp. 3449-3458.
12

---

[41] P. Grönquist, C. Yao, T. Ben-Nun, N. Dryden, P. Dueben, S. Li, T. Hoefler, Deep learning for post-processing ensemble weather forecasts, Philosophical Transactions of the Royal Society 379 (2194) (2021) 20200092.
[42] J. McCabe, J. Brown, Learning to assimilate in chaotic dynamical systems, Advances in Neural Information Processing Systems 34 (2021) 12237-12250
[43] M. Fey, J. E. Lenssen, Weichert, H. Müller, Splinecnn: Fast geometric deep learning with continuous B-spline kernels, in: Proceedings of the IEEE Conference Computer Vision and Pattern Recognition, 2018, pp. 869-877.
[44] P. Izmailov, W. J. Maddox, P. Kirichenko, T. Garipov, D. Vetrov, A. G. Wilson, Subspace inference for Bayesian deep learning, in: Uncertainty in Artificial Intelligence, PMLR, 2020, pp. 1169-1179.
[45] M. Tipping, The relevance vector machine, Advances in neural information processing systems 12 (1999).
[46] S. Wang, P. Perdikaris, Long-ting time integration of parametric evolution equations with physics- informed deeponets, arXiv preprint arXiv:2106.05384 (2021).
[47] Laube, M. O. Franz, G. Umlauf, Deep learning parametrization for B-spline curve approxi- mation, in: 2018 International Conference on 3D Vision (3DV), IEEE, 2018, pp. 691-699.
[48] C. Durkan, A. Bekasov, I. Murray, G. Papamakarios, Neural spline flows, Advances in neural information processing systems 32 (2019).
[49] C. de Boor, On the convergence of odd-degree spline interpolation, Journal of approximation theory 1 (4) (1968) 452-463.
[50] B. Swartz, 2 1) bounds Mathematical Society 74 (6) (1968) 1072-1078.
spline interpolation errors, Bulletin of the American
[51] J. Z. Huang, C. O. Wu, L. Zhou, Polynomial spline estimation and inference for varying coefficient models with longitudinal data, Statistica Sinica (2004) 763-788.
[52] C. De Boor, C. De Boor, A practical guide splines, Vol. 27, springer-verlag New York, 1978.
13

---

Checklist
Please do modify the questions and only use the provided macros for your answers. Note that the Checklist section does not count towards the page limit. In your paper, please delete this instructions block and only keep the Checklist section heading above along with the questions/answers below.
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] See Section 5
(c) Did you discuss any potential negative societal impacts of your work? [N/A]
(d) Have you read the ethies review guidelines and ensured that your paper conforms them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Section
(b) Did you include complete proofs of all theoretical results? [Yes] See supplemental materials
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or URL)? [Yes] It is included as a URL
(b) Did you specify all the training details (c.g., data splits, hyperparameters, how they were chosen)? [Yes] See supplemental materials
(c) Did you report error bars (e.g., with respect to the random seed after running experi-
ments multiple times)? [Yes] See Section and supplemental materials
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See supplemental materials
4. If you are using existing assets (c.g., code, data, models) or curating/releasing new assets...
(a) your work uses existing assets, did you eite the creators? [Yes]
(b) Did you mention the license of the assets? [No] It is general open source code available in public for research purpose
(c) Did you include any assets either in the supplemental material URL? [Yes]
(d) Did you discuss whether and how consent was obtained from people whose data you't 're using/curating? [N/A]
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent participant compensation? [N/A]
14


Reviews:
[{"review_id": "v6k-dinHgvi", "reviewer": "Reviewer_R79g", "summary": "The paper proposes a Bayesian spline learning framework to discover Nonlinear dynamics (ODE or PDE) using sparse data, which is interpretable and is able to quantify the uncertainty. Technical proposals are made to implement this strategy.", "questions": "I have several questions regarding clarifying the contribution of the proposed method.\n\n1. In 3.3, the authors assume that Aleatoric Uncertainty is $\\epsilon_1 \\sim N(0,B)$ and Epistemic Uncertainty is $\\epsilon_2 \\sim N(0,P)$. Since two noises are independent, the sum of two noises is $N(0, P+B)$, which can be understood as Aleatoric Uncertainty only. Can you elaborate how the proposed methodology differentiate between Aleatoric Uncertainty $N(0,B)$ and Epistemic Uncertainty $N(0,P)$ and how it affects the result?\n \n2. It looks like that the authors try to prove the spline fitting error bound with respect to the mean, which is not zero. Have you considered using this upper bound as a means of the Epistemic Uncertainty? Or have you considered other distributions(e.g. uniform distribution) for $\\epsilon_2$ (Epistemic Uncertainty)?\n\n3. Can you please provide additional details on how long it takes to train the model? The proposed method in the paper uses the posterior distribution to estimate uncertainty. However, the key challenge in the Bayesian Inference is estimating posterior distribution. I’m curious to see how much computation you need in order to make the result better.\n", "limitations": "1. The authors recognize that the proposed method relies on a predefined library of candidate terms\n\n2. There are some papers considering finding the best knot size $(\\tau_{s+1} - \\tau_{s})$ as a way to improve the expressibility of spline fitting function, for instance, \"Neural Spline Flow\", Conor Durkan, Artur Bekasov, Iain Murray and George Papamakarios. I think this method could potentially reduce the total computation required.\n", "rating": 5, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The strengths of the paper are listed below:\n1. The experiment results seem promising.\n2. The proposed methodology is interpretable and is able to quantify the uncertainty.\n\nThe weakness of the paper are listed below:\n1. The suggested method relies on the predefined library.\n2. The authors did not include the training time to compare the result. \n\nApplying spline based learning with alternating direction optimization to train Nonlinear Dynamics are original, significant. However, I have further questions about the assumptions and results.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "b6VcFT_Rzh", "reviewer": "Reviewer_nVYU", "summary": "This paper extended an existing spline-based model for fitting dynamic systems (ODEs) reported in [14]. There are two main changes: 1) Use tensor-product of splines to model both time and spatial dimensions. 2) Use Bayesian regression to fit the coefficients. With the first change, authors attempted to model PDEs with one spatial dimension (avection equation and Burgers' equation). \n\nThe authors also claimed that the Bayesian estimates help with the data assimilation tasks, although the approach was not described in detail.\n\n", "questions": "1. Equation (3), for PDE, how do you model terms like $u_x$ and $u_{xx}$? Please give an example like Burger's equation. Is the representation in $R^{n\\times s\\times m}$?\n2. Equation (10), what is $U$? Same for equation (13).\n3. Table 1, what are the values?\n", "limitations": "The authors should clarify the limitations of splines in modeling higher-dimension data.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n\nThe paper contains some reasonable ideas (tricks) that have not seemed to be explored before.\n\nWeaknessses:\n\n1. The paper contains quite a few typos in equations and counter-intuitive use of symbols (e.g., using $\\hat{U}$ to represent data instead of an estimate). These problems can make the method difficult to understand. The author should have been more careful proof-reading the draft before submission.\n2. The tensor-product of 1-d splines has severe limitations, particularly when the spatial dimension is higher. The authors should have made this limitation more clear.\n3. The experiment results do not contain sufficient comparison with baseline methods (Only Table 1 contains the comparison of, I assume, the running time).  ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "VkisoyWY6lG", "reviewer": "Reviewer_7qjr", "summary": "In this paper, the authors:\n- propose a system identification algorithm that leverages spline bases to fit sparsely fit noisy data, and provides parameter and prediction uncertainty (which can be used in an active learning / Bayesian optimization setting)\n- show that the algorithm can be used to infer ODEs and PDEs, and with high level of noise\n- show that the proposed algorithm performs well against two other standard / state of the art methods, but has tradeoffs at different noise levels", "questions": "- in Figure 2, the discovered mean of the library terms are plotted against their actual values: what is the actual value in these cases for the VdP and Lorenz96 systems? If I understand correctly, the output of the algorithm are the weights of the splines at different orders, which do not relate to the actual dynamical equations for the systems at all, so how are the comparisons made? Or do the library terms CX2 and CX3 correspond to, e.g., first and second order terms such as x, and xy? It would be nice if this point could be clarified in the text\n\n- in this vein, and it's very possible that I'm confused on this point, but the motivation was to find interpretable models that match the true system. For example, with SINDY, the library can be constructed with polynomial terms up to arbitrary orders, which can then be exactly matched with the original system. In this work, the basis functions are splines, which do not allow for such interpretations\n\n- the schematic in Figure 1 demonstrates a nice effort in explaining the algorithm, but it's a bit crowded, and it's especially unclear how the loss functions in the red box relate to each other, and what the correspond to in the main? In particular, it's a bit unclear to me how the data fitting loss (L2) is combined with the logL loss?\n\n- related to my confusion on the learning procedure, are the spline parameters separately learned (as indicated by the blue box being separated)? or is the whole thing trained end to end. I apologize for the confusion here and do not really have concrete suggestions for how to improve clarity.\n\n- while I understand that the algorithm comparisons are conducted on the benchmark ODE/PDEs, it would be nice to see how the proposed method works with real noisy data that is routinely modeled with one of the systems (e.g., advection equation), or another i.e. predator-prey data with LotkaVolterra, or COVID numbers with SIR model\n\n- figure fonts are extremely small, basically unreadable in many of them\n\n- some terminology should be explained for a reader unfamiliar with SI, e.g., \"control points\"", "limitations": "a brief discussion on library choice is included, but a longer discussion is warranted regarding the splines, as well as how it performs in high noise regimes (especially in comparison with the other algorithms). Furthermore, potential negative societal impact should be discussed as such an algorithm can have broad consequences in accurately modeling natural systems.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "originality: I don't know the literature well enough to know if others have attempted similar approaches, but in comparison to standard SINDY, this work builds upon the idea of constructing the system from sparse basis functions but additionally provides a probabilistic output\n\nquality: I did not check the math/code, but I believe the results are of high quality at a superficial level, and support the overall claim (i.e., uses splines to discovery system equations, is sparse, and provides uncertainty estimates). Unclear if it would work beyond the standard benchmark synthetic systems.\n\nclarity: the description of the algorithm is a bit unclear (specific issues outlined below), but the overall goal, the components of the algorithm, and the results were clearly communicated\n\nsignificance: system identification that tries to retrieve the underlying mechanistic model is especially important in the context of scientific applications with noisy observations; has potential to be broadly impactful in many areas of ML and science", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}, {"review_id": "SO3XPh-rl8_", "reviewer": "Reviewer_3NV7", "summary": "The authors propose a Bayesian Spline Learning method for equation discovery, which uses spline basis to handle the data scarcity and noise. The Bayesian uncertainty calibration techniques are used to approximate posterior distributions of the trainable parameters. The proposed method is built upon the ADO of PINN-SR, but outperforms PINN-SR and SINDy in all the metrics in the experiments of four equations.", "questions": "In the current presentation, it can still be hard to interpret why uncertainty calibration can benefit equation discovery against measurement noise and sparsity. We only know that the spline learning can perform differentiation accurately to enhance learning efficiency and the posterior distribution of spline, library and equation coefficients are estimated simultaneously. These explain the training time superiority but lacks in interpretation why the obtained uncertainty can benefit equation discovery/coefficient estimation, which is more at the core of the question. It makes it harder to understand why the proposed method has a smaller error than PINN-SR and SINDy, especially calibration method generally succeeds in efficiency of uncertainty estimation but hardly provides significant accuracy boost based on relevant studies of previous works on other topics such as traffic prediction, etc.", "limitations": "The authors mention the limitation of pre-defined library of candidate terms. ", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The motivation is clear and the used method makes sense in the current presentation.\n- The experimental results support the claim that their method is more robust with 10% data with large noise (5% to 20%) on Van der Pol Oscillator, Lorenz 96, Advection Equation and Burgers' Equation. \n- The codes and theoretical properties of the algorithm are provided in the supplementary materials.\n\nWeaknesses:\n- It is not clear what are the analytical forms of the used four equations and what are the parsimonious equations discovered. There should be some descriptions of the discovered equations at least in the supplementary materials if not in the main manuscript.\n- The presentation about why the propose method can benefit equation discovery can be further extended and described in detail.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept