# Bayesian Spline Learning for Equation Discovery of Nonlinear Dynamics with
Quantified Uncertainty

Luning Sun
University of Notre Dame
<EMAIL_ADDRESS>
Daniel Zhengyu Huang11footnotemark: 1
California Institute of Technology
<EMAIL_ADDRESS>
Hao Sun
Renmin University of China
<EMAIL_ADDRESS>
Jian-Xun Wang
University of Notre Dame
<EMAIL_ADDRESS>
Equal contributionCorresponding author

###### 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.

##  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 understood, preventing researchers from
developing first-principled models. On the other hand, the ever-increasing
data availability opens up new avenues for scientists to identify predictive
models from enormous observation data, a process known as _system
identification_ (SI). Recent advances in deep learning have prompted the rapid
development of powerful SI models for high-dimensional problems using deep
neural networks (DNNs). Many DNN-based SI models have been proposed to learn
differential operators for complex (spatio)temporal physics from data and
shown good potential in terms of data reproduction and state prediction [1, 2,
3]. 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 pruning
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 [7, 8, 9] or coupled manner [10, 11, 12, 13,
14]. For uncertainty quantification, the dictionary-based equation discovery
algorithms have been recently extended to Bayesian settings [15, 16, 17],
based on the idea of sparse Bayesian learning pioneered by Tipping and co-
workers [18, 19, 20, 21, 22].

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 parsimonious 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 a spline basis, on which a collection
of derivatives can be reliably computed for the library construction. The
posterior distributions of spline-based model parameters are 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 to 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.111The code will
be available at https://github.com/luningsun/SplineLearningEquation

##  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 [4, 5, 24, 25, 26, 6, 27, 16]. 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 [10, 11, 12, 13, 14, 7].
For example, automatic differentiation (AD) is effective for computing
derivatives analytically from sparse, noisy data using a point-wise multi-
layer perceptron (MLP) [29, 30, 13]. 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 enforcing 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-informed neural networks (PINN) in the context of solving PDEs
[31] and discovering ODEs [14].

However, these models fail to simultaneously 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 (PiSL) [14] that performs very well for ODE
discovery cannot handle PDE problems due to the limitation of the B-spline
basis 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 is natural to use Bayesian
framework to quantify and analyze the prediction uncertainty. To enable the
posterior computation for high-dimensional trainable parameters that are
intractable by traditional Bayesian inference, people have resorted to various
approximation strategies, e.g., variational inference [32, 33, 34], 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-Gaussian (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 a predictive model for nonlinear
dynamics forecast. People have recently integrate deep learning into DA to
improve online prediction performance [40, 41, 42]. However, the predictive
model in DA is assumed to 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 us consider a dynamical system, which is governed by a parametric ODE/PDE
system in the following general form

$\frac{d\mathbf{u}}{dt}=\mathcal{F}(\mathbf{u}),$ (1)

where $\mathbf{u}$ denotes the state vector, for ODE systems,
$\mathbf{u}=[u_{1}(t),u_{2}(t),\cdots,u_{d}(t)]^{T}\in\mathbb{R}^{d}$ depends
only on time $t$, for PDE systems, $\mathbf{u(x,}$
$t)=[u_{1}(\mathbf{x},t),u_{2}(\mathbf{x},t),\cdots,u_{d}(\mathbf{x},t)]^{d}\in\mathbb{R}^{d}$
depends on both time $t$ and space $\mathbf{x}$, and
$\mathcal{F}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ represents unknown nonlinear
functions. The states are observed at discrete times $\\{t_{i}\\}_{i=1}^{n}$
and at spatial locations $\\{\mathbf{x}_{j}\\}_{j=1}^{s}$ for PDE systems. Let
$\tilde{\mathbf{u}}$ denote the noisy observation vector, the observation set
is
$\tilde{\mathbf{U}}=\\{\tilde{\mathbf{u}}(t_{1}),\tilde{\mathbf{u}}(t_{2}),\cdots,\tilde{\mathbf{u}}(t_{n})\\}^{T}\in\mathbb{R}^{n\times
d}$ and
$\tilde{\mathbf{U}}=\\{\tilde{\mathbf{u}}(\mathbf{x}_{1},t_{1}),\tilde{\mathbf{u}}(\mathbf{x}_{2},t_{1}),\cdots,\tilde{\mathbf{u}}(\mathbf{x}_{s},t_{n})\\}^{T}\in\mathbb{R}^{(n\times
s)\times d}$ for ODE/PDE systems, respectively. Our goal here is to explicitly
discover the parsimonious form of $\mathcal{F}(\cdot)$ from a library of
candidate basis functions and quantify the associated uncertainty given noisy
observation data.

Figure 1: Overview of method. Black box: overall work flow. Blue box: a demo
for spline basis with derivatives in 1D and 2D cases. Red box: a sketch for
the Bayesian ADO training process. Green box: the data assimilation enhanced
uncertainty quantification.

###  3.1 Overview

Given the dataset $\tilde{\mathbf{U}}$, the SI problem can be solved by sparse
regression techniques with a predefined library $\mathbf{\Phi}(\mathbf{u})$ of
$m$ basis functions,

$\mathbf{\Phi}(\mathbf{u})=[\phi_{1}(\mathbf{u}),\phi_{2}(\mathbf{u}),\cdots,\phi_{m}(\mathbf{u})]\in\mathbb{R}^{m},$
(2)

where $\phi_{i}:\mathbb{R}^{d}\to\mathbb{R},1\leq i\leq m$ denotes a basis
function, which, for instance, can be the polynomial or trigonometric
function. Hence, the matrix of library terms evaluated at observed states is
defined as,

ODE :
$\displaystyle\mathbf{\Phi}(\tilde{\mathbf{U}})=\bigg{[}\mathbf{\Phi}\big{(}\tilde{\mathbf{u}}(t_{1})\big{)}^{T},\mathbf{\Phi}\big{(}\tilde{\mathbf{u}}(t_{2})\big{)}^{T},\cdots,\mathbf{\Phi}\big{(}\tilde{\mathbf{u}}(t_{n})\big{)}^{T}\bigg{]}^{T}\in\mathbb{R}^{n\times
m},$ (3) PDE :
$\displaystyle\mathbf{\Phi}(\tilde{\mathbf{U}})=\bigg{[}\mathbf{\Phi}\big{(}\tilde{\mathbf{u}}(\mathbf{x}_{1},t_{1})\big{)}^{T},\mathbf{\Phi}\big{(}\tilde{\mathbf{u}}(\mathbf{x}_{2},t_{1})\big{)}^{T},\cdots,\mathbf{\Phi}\big{(}\tilde{\mathbf{u}}(\mathbf{x}_{s},t_{n})\big{)}^{T}\bigg{]}^{T}\in\mathbb{R}^{(n\times
s)\times m}.$

Given the library, the overview of the BSL framework is shown in Fig. 1.
Firstly, a spline-based model is constructed to represent state variables and
their derivatives by denoising the observation data, which is illustrated in
the blue box and discussed in Sec.3.2. Then the spline-based model and sparse
regression are trained simultaneously by alternating direction optimization
(ADO), as shown in the red box. Specifically, a 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 in the library as defined in Eq. LABEL:eq:matrixlib, and then
updates the number of relevant terms in 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.

###  3.2 Spline-based learning

Several previous works have already shown the potential of spline-based
learning and demonstrated the advantages compared with the classical DL
structures, e.g.,(MLP, 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
spline-reconstructed fields. The B-spline curves, based on the De Boor’s
algorithm, are defined in a recursive way as:

$\displaystyle\begin{aligned} N_{s,0}(t)&=\begin{cases}1,&\text{if
$\tau_{s}\leq t\leq\tau_{s+1}$}.\\\ 0,&\text{otherwise}.\end{cases}\\\
N_{s,k}(t)&=\frac{t-\tau_{s}}{\tau_{s+k}-\tau_{s}}N_{s,k-1}(t)+\frac{\tau_{s+k+1}-t}{\tau_{s+k+1}-\tau_{s+1}}N_{s+1,k-1}(t),\end{aligned}$
(4)

where $\tau_{s}$ is the location of knots, $k$ 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)=\Sigma_{s=0}^{r+k-1}N_{s,k}(t)\theta_{s}$. And the number of control
points (trainable weights) $\boldsymbol{\theta}\in\mathbb{R}^{r+k}$ is chosen
empirically. It can be proved that the derivative of $p$ order B-spline basis
is a function of $p-1$ order B-spline, written as:

$\frac{d}{dt}N_{s,k}(t)=\frac{k}{\tau_{s+k}-\tau_{s}}N_{s,k-1}(t)-\frac{k}{\tau_{s+k+1}-\tau_{s+1}}N_{s+1,k-1}(t).$
(5)

The proof of Eq. 5 is attached in Appendix A.1. The higher-order derivative
can be calculated by recursively using Eq. 5. The analytical derivatives of
spline basis are very beneficial for PDE discovery tasks since it always
involves constructing library terms containing high order derivatives. With a
proper order $p$, the first $p$ derivatives are accurate, and there is no
error introduced during the derivation, as opposed to using the numerical
methods to approximate derivatives. The spline function and its derivatives
are defined in a one-dimensional scenario. It is straightforward to extend to
n-dimension by direct using tensor-product. For example, for a two-dimensional
problem with spatial-temporal fields, the basis can be defined as

$N_{s_{1},s_{2}}^{k_{1},k_{2}}(t_{i},x_{j})=N_{s_{1},k_{1}}(t_{i})N_{s_{2},k_{2}}(x_{j}).$
(6)

Here the two-dimensional basis is denoted by a different style to write it
compactly. Similarly, the partial derivative for two-dimensional basis is
defined as:

$\frac{\partial^{(q_{1}+q_{2})}N_{s_{1},s_{2}}^{k_{1},k_{2}}}{\partial
t^{(q_{1})}\partial
x^{(q_{2})}}(t_{i},x_{j})=\frac{d^{(q_{1})}N_{s_{1},k_{1}}}{dt^{(q_{1})}}(t_{i})\frac{d^{(q_{2})}N_{s_{2},k_{2}}}{dx^{(q_{2})}}(x_{j}).$
(7)

With the definition of sparse system identification and spline reconstruction,
the whole spline learning can be stated as follows: given noisy measurement
data $\tilde{\mathbf{U}}$, find the best sets of weights $\boldsymbol{\theta}$
and $\mathbf{W}$ so that data fitting loss and the weakly physics-informed
loss can be minimized under sparsity constraints , as shown in Eq. 8:

$\displaystyle\\{\boldsymbol{\theta},\mathbf{W}\\}$
$\displaystyle=\operatorname*{arg\,min}_{\boldsymbol{\theta}^{\prime},\mathbf{W}^{\prime}}\frac{1}{\mathcal{N}_{m}}||\mathbf{N}_{m}\boldsymbol{\theta}^{\prime}-\tilde{\mathbf{U}}||_{L_{2}}+\frac{1}{\mathcal{N}_{c}}||\mathbf{\Phi}(\mathbf{N}_{c}\boldsymbol{\theta}^{\prime})\mathbf{W}^{\prime}-\dot{\mathbf{N}}_{c}\boldsymbol{\theta}^{\prime}||_{L_{2}}+\lambda||\mathbf{W}^{\prime}||_{L_{\alpha}}.$
(8)

Here, $\mathbf{N}_{m}$ and $\mathbf{N}_{c}$ denote the spline basis matrices
evaluated at measurement and collocation locations. Moreover,
$\mathcal{N}_{m}$ and $\mathcal{N}_{c}$ are numbers of measurement data and
collocation points. Furthermore, the alternating direction optimization (ADO)
shown in previous works [13, 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 a 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): (1) Observation error, where the data is
noisy, and the smoothed data is reconstructed through spline-based learning,
as shown in Eq. 9

$\tilde{\mathbf{U}}=\mathbf{N}_{m}\boldsymbol{\theta}^{\prime}+\boldsymbol{\epsilon}_{1},$
(9)

where $\boldsymbol{\epsilon}_{1}$ 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 as,

$\dot{\mathbf{u}}=\boldsymbol{\Phi}(\mathbf{u})\boldsymbol{W^{\prime}}+\boldsymbol{\epsilon}_{2},$
(10)

where $\boldsymbol{\epsilon}_{2}$ represents model form error. And Eq. (10) is
evaluated on collocation points.

In this work, these error terms are modeled as zero-mean multivariate Gaussian
random variables:
$\boldsymbol{\epsilon_{1}}\sim\mathcal{N}(\mathbf{0},\mathbf{B})$ and
$\boldsymbol{\epsilon_{2}}\sim\mathcal{N}(\mathbf{0},\mathbf{P})$. And we
further assume that these error covariance matrices $\mathbf{B},\mathbf{P}$
are diagonal matrix with diagonal terms $\\{b_{k}\\}$ and $\\{p_{k}\\}$ with
$1\leq k\leq d$ and they are learn-able parameters during training.

According to Bayes’ rule, the posterior can be written as:

$p(\boldsymbol{\theta},\mathbf{W},\mathbf{B},\mathbf{P}|\tilde{\mathbf{U}},\dot{\mathbf{U}})\propto
p(\boldsymbol{\theta},\mathbf{W},\mathbf{B},\mathbf{P})p(\tilde{\mathbf{U}},\dot{\mathbf{U}}|\mathbf{W},\boldsymbol{\theta},\mathbf{B},\mathbf{P}).$
(11)

Here $\dot{\mathbf{U}}=\dot{\mathbf{N}}_{c}\boldsymbol{\theta}$ denotes the
derivative estimation on all the collocation points. The prior is further
decomposed by

$p(\mathbf{W},\boldsymbol{\theta},\mathbf{B},\mathbf{P})\propto
p(\mathbf{W}|\alpha)p(\alpha)p(\boldsymbol{\theta}|\beta)p(\beta)p(\mathbf{B})p(\mathbf{P}).$
(12)

Currently, we specify the prior for the linear coefficient matrix as a zero
mean Gaussian distribution
$p(\mathbf{W}|\alpha)=\mathcal{N}(\mathbf{W}|0,\alpha^{-1}\mathbf{I})$ with
the hyper prior as a Gamma distribution
$p(\alpha)=\text{Gamma}(\alpha|a_{0},b_{0})$. Similarly, we also define the
prior for the spline trainable parameters as zero mean Gaussian distribution
with the hyper prior as another Gamma distribution
$p(\mathbf{\boldsymbol{\theta}}|\beta)=\mathcal{N}(\boldsymbol{\theta}|0,\beta^{-1}\mathbf{I})$
and $p(\beta)=\text{Gamma}(\beta|a_{1},b_{1})$. To account for the data
uncertainty and process uncertainty, the diagonal covariance matrices
$\mathbf{B}$ and $\mathbf{P}$ are also set as learn-able during the training.
And hence improper uniform priors are used for $\mathbf{B}$ and $\mathbf{P}$.
The likelihood consists of two parts, as:

$p(\tilde{\mathbf{U}},\dot{\mathbf{U}}|\mathbf{W},\boldsymbol{\theta},\mathbf{B},\mathbf{P})=p(\tilde{\mathbf{U}}|\boldsymbol{\theta},\mathbf{B})p(\dot{\mathbf{U}}|\mathbf{W},\mathbf{P},\boldsymbol{\theta}),$
(13)

where

$\displaystyle p(\tilde{\mathbf{U}}|\boldsymbol{\theta},\mathbf{B})$
$\displaystyle\propto\exp\\{-\frac{1}{2}(\tilde{\mathbf{U}}-\mathbf{N}_{m}\boldsymbol{\theta})^{T}\mathbf{B}^{-1}(\tilde{\mathbf{U}}-\mathbf{N}_{m}\boldsymbol{\theta})\\},$
(14) $\displaystyle
p(\dot{\mathbf{U}}|\mathbf{W},\mathbf{P},\boldsymbol{\theta})$
$\displaystyle\propto\exp\\{-\frac{1}{2}(\dot{\mathbf{N}}_{c}\boldsymbol{\theta}-\boldsymbol{\Phi}(\mathbf{N}_{c}\boldsymbol{\theta})\boldsymbol{W})^{T}\mathbf{P}^{-1}(\dot{\mathbf{N}}_{c}\boldsymbol{\theta}-\boldsymbol{\Phi}(\mathbf{N}_{c}\boldsymbol{\theta})\boldsymbol{W})\\}.$
(15)

Traditional Bayesian sampling approaches (i.e., Markov chain Monte Carlo
methods) can be usually 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 a subspace by finding the PCA
components of the SWAG trajectories [44] and then draw samples in the subspace
instead. In terms of the loss function in the probabilistic model, we chose to
maximize the log form of the posterior density Eq. 11, also by leveraging the
ADO algorithm. The sub-process I minimizes the log density function defined by
Eq. 11-13 and the sub-process II adopts the Bayesian variants of SINDy
algorithm modified from [19, 17]. The likelihood function in sub-process II
has the same form as in Eq. 15 but with a different sparsity promoting prior
$p(\mathbf{W}|\mathbf{A})={\displaystyle\prod_{j=1}^{m}\mathcal{N}(\mathbf{W}_{j}|0,{\alpha^{\prime}}_{j}^{-1})}$,
where
$\boldsymbol{A}=[\alpha^{\prime}_{1},\alpha^{\prime}_{2},...,\alpha^{\prime}_{m}]^{T}$.
In a single ADO iteration, the sub-process I provides the updated $\theta$ and
$\mathbf{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 $\mathbf{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 in 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 a 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 $\mathbf{B}$ and evolution error
covariance matrix $\mathbf{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.

##  4 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 can 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 is
defined as $\frac{dx}{dt}=y$ and $\frac{dy}{dt}=\mu(1-x^{2})y-x$ with
$\mu=0.5$. The data is corrupted with $5\%$ noise, and the library consists of
polynomials of state variables up to $3rd$ order. The parsimonious model
structure can be correctly identified, and the result is shown in the upper
red box of Fig. 2. 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
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\boldsymbol{-}}}$)
and the truth equation coefficients (black stars
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{*}}}$),
where the horizontal axis represents term indices. For example, the Van der
Pol system has four different terms, and the x-axis ranges from 1 to 4. The
two sub-figures in the lower left part show probability density distributions
(PDF) of two identified coefficients, where truth equation coefficients (red
lines
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\boldsymbol{-}}}$)
fall within the confidence interval with high probability. The right sub-
figures show the propagated ensemble results (blue lines
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\boldsymbol{-}}}$)
based on the discovered equations, measurement data (green circles
$\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}{\mathbf{\circ}}$),
and the true state trajectories of $x$ and $y$ (red dashed lines
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\boldsymbol{--}}}$).
The result clearly shows that the ensemble predictions can cover the data, and
the uncertainty range fluctuates around the truth state values.

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
is $\frac{dx}{dt}=y$ and $\frac{dy}{dt}=-x-0.5x^{2}y+0.5y$. For the Lorenz96
system, the compact form of the governing equation is
$\frac{dX_{i}}{dt}=(X_{i+1}-X_{i-2})X_{i-1}-X_{i}+F$ with periodic boundary
conditions.

The second ODE example is a chaotic system, Lorenz 96. It is a simplified
mathematical model for atmospheric convection, defined as:
$\frac{dX_{i}}{dt}=(X_{i+1}-X_{i-2})X_{i-1}-X_{i}+F$, $i=1,2,...n$ with
periodic boundary conditions $X_{-1}=X_{n-1},X_{0}=X_{n}$ and $X_{n+1}=X_{1}$.
We chose $n={4}$ in current case and $F=8$ for the forcing terms. The
measurement states variables are corrupted with $10\%$ Gaussian noise, and the
library consists of polynomials of state variables up to $3rd$ order. The
discovery result for these systems is shown in the lower blue box in Fig. 2
following the same layout as 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, 6 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 a large phase difference compared with truth
trajectories due to 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), a classical DA
technique, to assimilate noisy measurements with the identified system. The
observation error matrix $\boldsymbol{B}$ and evolution error matrix
$\boldsymbol{P}$ required for the EnKF scheme can be directly passed from the
Bayesian spline learning framework. The right subplots show the DA results for
the 3 out of 6 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
($\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}{\mathbf{\circ}}$)).
And the red solid-line
($\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\boldsymbol{-}}$)
is the true L96 states. Finally, the blue curves
(${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\boldsymbol{-}}}$)
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 is only about 1 second. More discovery results
are attached in Appendix A.5.

h

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 $x$.

Then, we use the proposed BSL framework to identify classical PDE systems and
evaluate the performance of two of them in the main text. They include
advection equation $\frac{\partial u}{\partial t}+\frac{\partial u}{\partial
x}=0$ and Burgers’ equation $\frac{\partial u}{\partial t}+u\frac{\partial
u}{\partial x}-0.5\frac{\partial^{2}u}{\partial x^{2}}=0$. The qualitative
results are shown in Fig. 3 inside the 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 red lines
(${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\boldsymbol{-}}}$)
represent the truth value, blue lines
(${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\boldsymbol{-}}}$)
are obtained by ensemble predictions, and green circles
(${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}{\mathbf{\circ}}}$)
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.

Finally we benchmark our proposed method with several state-of-the-art
discovery algorithms (PINN-SR [13], SINDy [5] and RVM [45]). The error metric
is defined as:

rmse $\displaystyle=\frac{||\mathbf{C}_{\rm Discovery}-\mathbf{C}_{\rm
True}||_{2}}{||\mathbf{C}_{\rm True}||_{2}}$ (16)
$\displaystyle{\mathbf{M_{P}}}$ $\displaystyle={\frac{||\mathbf{C}_{\rm
Discovery}\odot\mathbf{C}_{\rm True}||_{0}}{||\mathbf{C}_{\rm
Discovery}||_{0}}}$ (17) $\displaystyle{\mathbf{M_{R}}}$
$\displaystyle={\frac{||\mathbf{C}_{\rm Discovery}\odot\mathbf{C}_{\rm
True}||_{0}}{||\mathbf{C}_{\rm True}||_{0}}}$ (18)

where $\mathbf{C}_{\rm Discovery}$ are the non-zero mean prediction from the
posterior distribution and $\mathbf{C}_{\rm True}$ are the true coefficients
of the governing equations. If the method fails to converge or cannot identify
the correct parsimonious form, we will report the final result as Fail. The
errors are scaled by $\times 10^{-3}$ to have a clear comparison. Two
additional metrics, precision $M_{P}$ and recall $M_{R}$, are also defined,
where the $\odot$ represents element-wise product of vectors and the $l_{0}$
norm is the non-zero terms in a vector. It can be seen from Table 1 that our
BSL method always performs best, when the noise is significant ($>5\%$),
demonstrating its robustness to noise.

The PINN-SR can discover the PDE equation with corrupted data set, but it
fails to predict accurate time trajectories for ODE systems. This is also
reported in a relevant paper [46], that the plain MLP structure is 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\%$) |  rmse (large222Large noise for different cases: Van der Pol: $5\%$, Lorenz 96: $10\%$, Advection: $20\%$, Burgers: $10\%$, Burgers’ with source: $20\%$, Heat: $15\%$, Poisson: $5\%$) | $\mathbf{M_{P}}$ |  $\mathbf{M_{R}}$333$\mathbf{M_{P}},\mathbf{M_{R}}$ are only reported for the largest noise cases |  Training Cost444All cases are running on a Nvidia 2070 Ti GPU card
---|---|---|---|---|---|---
Van der Pol Oscillator
BSL(Ours) | $\mathbf{0.2}$ | $2.82$ | $\mathbf{18.04}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 133(+3)s$
PINN-SR |  Fail555Fail means failure in discovery of the parsimonious ODE/PDE forms. | Fail | Fail | $0.214$ | $0.75$ | $\sim 1213s$
SINDy | $1.0$ | $\mathbf{1.93}$ | Fail | $0.267$ | $1.0$ | $\sim 10s$
RVM | $1.0$ | $2.54$ | $27.46$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 10s$
Lorenz 96
BSL(Ours) | $\mathbf{0.269}$ | $1.47$ | $\mathbf{13.0}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 1654(+438)s$
PINN-SR | Fail | Fail | Fail | $0.5$ | $0.22$ | $\sim 10788s$
SINDy | $0.4$ | $0.64$ | Fail | $0.75$ | $\mathbf{1}$ | $\sim 10s$
RVM | $0.4$ | $\mathbf{0.6}$ | $49.7$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 25s$
Advection Equation
BSL(Ours) | $\mathbf{0.26}$ | $\mathbf{1}$ | $\mathbf{1.9}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 946(+233)s$
PINN-SR | $5.9$ | $4.5$ | $30.4$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 650s$
SINDy | $2.3$ | $8.2$ | $38.9$ | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 10s$
RVM | $0.77$ | $6.76$ | Fail | ${0.2}$ | ${\mathbf{1}}$ | $\sim 4s$
Burgers’ Equation
BSL(Ours) | $3.62$ | $4.13$ | $\mathbf{6.38}$ | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 117(+74)s$
PINN-SR | $10.2$ | $\mathbf{3.3}$ | $10.3$ | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 512s$
SINDy | $0.826$ | Fail | Fail | ${\mathbf{1}}$ | ${0.5}$ | $\sim 10s$
RVM | $\mathbf{0.754}$ | Fail | Fail | ${0.1429}$ | ${0.5}$ | $\sim 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 us to accurately interpolate solution surfaces and
analytically compute derivatives to form the candidate library, outperforming
other benchmark methods based on finite difference (FD) or 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 sparsity-promoting ADO iteration strategy is
proposed 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 to 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-art methods.

Admittedly, this work still relies on a pre-defined library of candidate
terms, and thus the identified 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 is adopted for multi-dimensional spline constructions
as shown in Tab. 7. These choices are not optimal and have notable limitations
for high-dimensional 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 in high-
dimensional settings [47, 48], while defining continuous spline kernels and
expanding basis in subdomains will allow a 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 a negative societal impact if it is used abusively, particularly
for predictive modeling of high-consequence natural systems, where caution
should be taken for decision making.

## Acknowledgments and Disclosure of Funding

The authors would like to acknowledge the funds from National Science
Foundation, United States of America under award numbers CMMI-1934300 and
OAC-2047127, the Air Force Office of Scientific Research (AFOSR), United
States of America under award number FA9550-22-1-0065, and startup funds from
the College of Engineering at University of Notre Dame in supporting this
study

## 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, G. 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–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.
  * [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] S. H. Rudy, S. L. Brunton, J. L. Proctor, J. N. Kutz, Data-driven discovery of partial differential equations, Science advances 3 (4) (2017) 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] S. 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 a 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 AI-based reconstruction algorithms for magnetic particle tracking, Plos one 16 (7) (2021) e0254051.
  * [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-informed learning of governing equations from scarce data, Nature Communications 12 (2021) 6136.
  * [14] F. Sun, Y. Liu, H. Sun, Physics-informed spline learning for nonlinear dynamics discovery, in: Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence (IJCAI-21), 2021, pp. 2054–2061.
  * [15] S. Zhang, G. Lin, Robust data-driven discovery of governing physical laws with error bars, Proceedings 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.07788 (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. E. 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.3838 (2013).
  * [21] A. C. Faul, M. E. Tipping, A variational approach to robust regression, in: International Conference on Artificial Neural Networks, Springer, 2001, pp. 95–102.
  * [22] A. C. Faul, M. E. 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 in 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. L. 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) 948–972.
  * [28] M. Tancik, P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ramamoorthi, 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-informed 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-informed hermite-spline cnns, arXiv preprint arXiv:2109.07143 (2021).
  * [32] C. Blundell, J. Cornebise, 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. Elhoseiny, T. Darrell, M. Rohrbach, Uncertainty-guided continual learning with bayesian neural networks, in: International Conference on 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-Conference Track Proceedings, Vol. 6, International Conference on Representation Learning, 2018.
  * [38] A. G. Wilson, P. Izmailov, Bayesian deep learning and a probabilistic perspective of generalization, 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, J. Smith, D. Cremers, M. Brenner, S. Hoyer, Variational data assimilation with a learned inverse observation operator, in: International Conference on Machine Learning, PMLR, 2021, pp. 3449–3458.
  * [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 A 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, F. Weichert, H. Müller, Splinecnn: Fast geometric deep learning with continuous B-spline kernels, in: Proceedings of the IEEE Conference on 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-time integration of parametric evolution equations with physics-informed deeponets, arXiv preprint arXiv:2106.05384 (2021).
  * [47] P. Laube, M. O. Franz, G. Umlauf, Deep learning parametrization for B-spline curve approximation, 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, ${O}(h^{2n+2-l})$ bounds on some spline interpolation errors, Bulletin of the American Mathematical Society 74 (6) (1968) 1072–1078.
  * [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 to splines, Vol. 27, springer-verlag New York, 1978.

## Checklist

Please do not 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. 1.

For all authors…

    1. (a)

Do the main claims made in the abstract and introduction accurately reflect
the paper’s contributions and scope? [Yes]

    2. (b)

Did you describe the limitations of your work? [Yes] See Section 5

    3. (c)

Did you discuss any potential negative societal impacts of your work? [Yes]

    4. (d)

Have you read the ethics review guidelines and ensured that your paper
conforms to them? [Yes]

  2. 2.

If you are including theoretical results…

    1. (a)

Did you state the full set of assumptions of all theoretical results? [Yes]
See Section 3

    2. (b)

Did you include complete proofs of all theoretical results? [Yes] See
supplemental materials

  3. 3.

If you ran experiments…

    1. (a)

Did you include the code, data, and instructions needed to reproduce the main
experimental results (either in the supplemental material or as a URL)? [Yes]
It is included as a URL

    2. (b)

Did you specify all the training details (e.g., data splits, hyperparameters,
how they were chosen)? [Yes] See supplemental materials

    3. (c)

Did you report error bars (e.g., with respect to the random seed after running
experiments multiple times)? [Yes] See Section 4 and supplemental materials

    4. (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. 4.

If you are using existing assets (e.g., code, data, models) or
curating/releasing new assets…

    1. (a)

If your work uses existing assets, did you cite the creators? [Yes]

    2. (b)

Did you mention the license of the assets? [No] It is general open source code
available in public for research purpose

    3. (c)

Did you include any new assets either in the supplemental material or as a
URL? [Yes]

    4. (d)

Did you discuss whether and how consent was obtained from people whose data
you’re using/curating? [N/A]

    5. (e)

Did you discuss whether the data you are using/curating contains personally
identifiable information or offensive content? [N/A]

  5. 5.

If you used crowdsourcing or conducted research with human subjects…

    1. (a)

Did you include the full text of instructions given to participants and
screenshots, if applicable? [N/A]

    2. (b)

Did you describe any potential participant risks, with links to Institutional
Review Board (IRB) approvals, if applicable? [N/A]

    3. (c)

Did you include the estimated hourly wage paid to participants and the total
amount spent on participant compensation? [N/A]

##  Appendix A Appendix

###  A.1 Proof: analytical form of spline derivatives

######  Theorem 1.

. The derivative of B-spline curve can be analytically evaluated as:

$\dot{y}(t)=\Sigma_{s=0}^{r+k-1}\dot{N}_{s,k}(t)\theta_{s},$ (19)

where

$\dot{N}_{s,k}(t)=\frac{k}{\tau_{s+k}-\tau_{s}}N_{s,k-1}(t)-\frac{k}{\tau_{s+k+1}-\tau_{s+1}}N_{s+1,k-1}(t).$
(20)

###### Proof.

We will prove by the induction. For the base case $k=1$, we have

$\dot{N}_{s,1}(t)=\begin{cases}\frac{1}{\tau_{s+1}-\tau_{s}}&t\in[\tau_{s},\tau_{s+1}]\\\
-\frac{1}{\tau_{s+2}-\tau_{s+1}}&t\in[\tau_{s+1},\tau_{s+2}]\end{cases}=\frac{1}{\tau_{s+1}-\tau_{s}}N_{s,0}(t)-\frac{1}{\tau_{s+2}-\tau_{s+1}}N_{s+1,0}(t).$
(21)

Let’s suppose that the formula holds for $k$ up to $n$. We will prove that
this formula also holds for $k=n+1$. By the definition in Eq. 4 and the chain
rule, we can get that:

$\displaystyle\dot{N}_{s,n+1}(t)$
$\displaystyle=\frac{t-\tau_{s}}{\tau_{s+n+1}-\tau_{s}}\dot{N}_{s,n}(t)+\frac{N_{s,n}(t)}{\tau_{s+n+1}-\tau_{s}}+\frac{\tau_{s+n+2}-t}{\tau_{s+n+2}-\tau_{s+1}}\dot{N}_{s+1,n}(t)-\frac{N_{s+1,n}(t)}{\tau_{s+n+2}-\tau_{s+1}}$
(22)
$\displaystyle=\frac{t-\tau_{s}}{\tau_{s+n+1}-\tau_{s}}\bigg{(}\frac{n}{\tau_{s+n}-\tau_{s}}N_{s,n-1}(t)-\frac{n}{\tau_{s+n+1}-\tau_{s+1}}N_{s+1,n-1}(t)\bigg{)}$
$\displaystyle+\frac{\tau_{s+n+2}-t}{\tau_{s+n+2}-\tau_{s+1}}\bigg{(}\frac{n}{\tau_{s+n+1}-\tau_{s+1}}N_{s+1,n-1}(t)-\frac{n}{\tau_{s+n+2}-\tau_{s+2}}N_{s+2,n-1}(t)\bigg{)}$
$\displaystyle+\frac{1}{\tau_{s+n+1}-\tau_{s}}N_{s,n}(t)-\frac{1}{\tau_{s+n+2}-\tau_{s+1}}N_{s+1,n}(t)$
$\displaystyle=\frac{n+1}{\tau_{s+n+1}-\tau_{s}}N_{s,n}(t)-\frac{n+1}{\tau_{s+n+2}-\tau_{s+1}}N_{s+1,n}(t).$

∎

###  A.2 Spline representation

In this section, we give error bounds for spline representation. For
simplicity, we consider 1D scenario and assume the target function
$u:[0,1]\rightarrow R$ is periodic and defined on the unit interval
$\Omega=[0,1]$. Consider a set of uniform knots
$\Gamma:0=\tau_{0}\leq\tau_{1}\leq\cdots\leq\tau_{r+k}=1$ with . The space of
$k$th degree $\\{N_{s,k}\\}_{s=0}^{r+k-1}$ splines is

$\displaystyle S^{k}(\Omega,\Gamma)=\\{p|p(t)\textrm{ is a polynomial of
degree $k$ in each }(\tau_{i},\tau_{i+1})\\}\cap\mathcal{C}^{k-1}(\Omega).$

Spline interpolation seeks $\hat{u}\in S^{k}(\Omega,\Gamma)$ that satisfies
$u(\tau_{i})=\hat{u}(\tau_{i})\quad\forall\quad 0\leq r+k$.

######  Theorem 2 (Spline interpolation error bounds [49, 50]).

Assume that $u$ is periodic, for odd degree spline interpolation $k$, we have

$\displaystyle\lVert
u-\hat{u}\rVert_{\infty}=\Delta\tau^{k+1}\Big{(}C_{k+1}\lVert
u^{(k+1)}\rVert_{\infty}+\mathcal{O}(\omega(u^{(k+1)},\Delta\tau))\Big{)},$
(23) $\displaystyle\lVert
u^{(l)}-\hat{u}^{(l)}\rVert_{\infty}=\Delta\tau^{k+1-l}\Big{(}D_{k+1-l}\lVert
u^{(k+1)}\rVert_{\infty}+\mathcal{O}(\omega(u^{(k+1)},\Delta\tau))\Big{)}\quad\forall
1\leq l\leq k-1.$

Here $C_{k+1}$ and $D_{k+1-l}$ are constant parameters, which are independent
of $\Delta\tau$ and $u$, and
$\omega(u^{(k+1)},\Delta\tau)=\sup_{|x-y|\leq\Delta\tau}|u^{(k+1)}(x)-u^{(k+1)}(y)|$.

In the present work, we focus on using spline for smoothing noisy data. It is
essentially a nonparameteric estimation of function $u$ from noisy data
$\\{t_{i},\tilde{u}_{i}\\}_{i=1}^{n}$,

$\displaystyle\tilde{u}_{i}=u(t_{i})+\eta_{i}\quad\textrm{ with
}\quad\eta_{i}\sim\mathcal{N}(0,\delta_{\eta}^{2}).$

We further assume that $\\{t_{i}\\}$ are uniformly sampled from $\Omega$. The
unregularized smoothing process estimates spline basis coefficients

$\hat{\boldsymbol{\theta}}=(\mathbf{N}_{m}^{T}\mathbf{N}_{m})^{-1}\mathbf{N}_{m}^{T}\cdot\tilde{\mathbf{U}}$

through minimizing
$\lVert\tilde{\mathbf{U}}-\mathbf{N}_{m}\cdot\boldsymbol{\theta}\rVert$. Here
$\tilde{\mathbf{U}}=\bigl{[}\tilde{u}_{1},\tilde{u}_{2},\cdots,\tilde{u}_{n}\bigr{]}^{T}$
and $\mathbf{N}_{m}$ is the spline basis matrix evaluated at these measurement
locations. The optimal $\boldsymbol{\theta}^{\rm opt}$ satisfies
$\hat{u}(t)=\mathbf{N}(t)\cdot\boldsymbol{\theta}^{\rm opt}$, where $\hat{u}$
is the spline interpolation of $u$, and
$\mathbf{N}(t)=[N_{0,k}(t),N_{1,k}(t),\cdots,N_{r+k-1,k}(t)]^{T}$ denotes
spline basis function vector. Following [51], we have spline fitting error
bounds, as following.

######  Theorem 3 (Spline fitting error bounds).

Assume that $u$ is periodic and the number of data $n$ is sufficient large,
for odd degree spline interpolation $k$, we have

$\displaystyle\lVert\mathbb{E}\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{\rm
opt}\lVert_{2}\leq
C_{1}\Delta\tau^{k+1}\qquad\lVert\mathrm{Cov}\hat{\boldsymbol{\theta}}\lVert_{2}\leq
C_{2}\frac{\delta_{\eta}^{2}}{n},$

where $C_{1}$ and $C_{2}$ are constant and independent of $n$.

###### Proof.

Let us denote

$\displaystyle\mathbf{U}=[u(t_{1}),u(t_{2}),\cdots,u(t_{n})]^{T}\quad\hat{\mathbf{U}}=[\hat{u}(t_{1}),\hat{u}(t_{2}),\cdots,\hat{u}(t_{n})]^{T}$
$\displaystyle\boldsymbol{\eta}=[\eta_{1},\eta_{2},\cdots,\eta_{n}]^{T}\quad\mathbf{e}=\mathbf{U}-\hat{\mathbf{U}}.$

We have

$\begin{split}\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{opt}&=(\mathbf{N}_{m}^{T}\mathbf{N}_{m})^{-1}\mathbf{N}_{m}^{T}\cdot\tilde{\mathbf{U}}-(\mathbf{N}_{m}^{T}\mathbf{N}_{m})^{-1}\mathbf{N}_{m}^{T}\mathbf{N}_{m}{\boldsymbol{\theta}}^{\rm
opt}\\\
&=(\mathbf{N}_{m}^{T}\mathbf{N}_{m})^{-1}\mathbf{N}_{m}^{T}\cdot(\mathbf{U}+\boldsymbol{\eta}-\hat{\mathbf{U}})\\\
&=(\mathbf{N}_{m}^{T}\mathbf{N}_{m})^{-1}\mathbf{N}_{m}^{T}\cdot(\mathbf{e}+\boldsymbol{\eta}).\end{split}$
(24)

We will first prove that

$\displaystyle\lVert\mathbf{N}_{m}^{T}\mathbf{N}_{m}\rVert_{2}=\mathcal{O}(n)\qquad\lVert(\mathbf{N}_{m}^{T}\mathbf{N}_{m})^{-1}\rVert_{2}=\mathcal{O}(\frac{1}{n}).$
(25)

For any $\boldsymbol{\theta}$,
$\frac{1}{n}\boldsymbol{\theta}^{T}\mathbf{N}_{m}^{T}\mathbf{N}_{m}\boldsymbol{\theta}$
is the Monte Carlo approximation of
$\int\bigl{(}\mathbf{N}(t)^{T}\cdot\boldsymbol{\theta}\bigr{)}^{2}dt$, and
hence

$\displaystyle\frac{1}{n}\boldsymbol{\theta}^{T}\mathbf{N}_{m}^{T}\mathbf{N}_{m}\boldsymbol{\theta}=\int(\mathbf{N}(t)^{T}\boldsymbol{\theta})^{2}dt+\mathcal{O}(\frac{1}{\sqrt{n}}).$
(26)

Bringing the following property of B-splines [52]

$\displaystyle
M_{1}\boldsymbol{\theta}^{T}\boldsymbol{\theta}\leq\int(\mathbf{N}(t)^{T}\boldsymbol{\theta})^{2}dt\leq
M_{2}\boldsymbol{\theta}^{T}\boldsymbol{\theta}\qquad\exists M_{1},M_{2}>0$

into Eq. (26) leads to Eq. (25). Then we prove that

$\displaystyle\lVert\mathbf{N}_{m}^{T}\mathbf{e}\rVert_{2}=\mathcal{O}(n\Delta\tau^{k+1}).$
(27)

Since $\sum_{s}N_{s,k}(t)=1$ and
$\lVert\mathbf{e}\rVert_{\infty}=\mathcal{O}(\Delta\tau^{k+1})$, we have

$\displaystyle\lVert\mathbf{N}_{m}^{T}\mathbf{e}\rVert_{2}\leq\lVert\mathbf{N}_{m}^{T}\mathbf{e}\rVert_{1}=\mathcal{O}(n\Delta\tau^{k+1}).$

Finally, combining Eq. (24), Eq. (25) and Eq. (27) leads to

$\displaystyle\lVert\mathbb{E}\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{\rm
opt}\lVert_{2}=\lVert(\mathbf{N}_{m}^{T}\mathbf{N}_{m})^{-1}\mathbf{N}_{m}^{T}\cdot\mathbf{e}\lVert_{2}\leq\lVert(\mathbf{N}_{m}^{T}\mathbf{N}_{m})^{-1}\lVert_{2}\lVert\mathbf{N}_{m}^{T}\cdot\mathbf{e}\lVert_{2}\leq
C_{1}\Delta\tau^{k+1}$
$\displaystyle\lVert\mathrm{Cov}\hat{\boldsymbol{\theta}}\lVert_{2}=\sigma_{\eta}^{2}\lVert(\mathbf{N}_{m}^{T}\mathbf{N}_{m})^{-1}\lVert_{2}\leq
C_{2}\frac{\delta_{\eta}^{2}}{n}$

∎

In our sparse Bayesian regression, in stead of solving the aforementioned
minimization problem, we have additional regularization terms.

###  A.3 Algorithms

In this section, we present detailed algorithms used in the present work,
which include Bayesian Alternative Direction Optimization (ADO) Learning 1,
Sequential Threshold Sparse Bayesian Learning 2, and Ensemble Kalman Filter 3.

Input : Library $\mathbf{\Phi}$, spline basis $\mathbf{N}$, time derivative of
spline basis matrix $\dot{\mathbf{N}}_{c}$, negative log form of equation Eq.
11 $\mathcal{L}$

Output : Mean estimation: $\boldsymbol{\theta}_{\text{SWA}}$,
$\mathbf{W}_{\text{SWA}}$, $\mathbf{B}_{\text{SWA}}$,
$\mathbf{P}_{\text{SWA}}$

Samples from posterior distributions: $\mathring{\boldsymbol{\theta}}$,
$\mathring{\mathbf{W}}$, $\mathring{\mathbf{B}}$, $\mathring{\mathbf{P}}$

Pretrain:

for  _$i=1:T_{\rm Pretrain}$_ do

SDG optimization with the fixed library

$\\{\boldsymbol{\theta}_{i+1},\mathbf{W}^{\prime}_{i+1},\mathbf{B}_{i+1},\mathbf{P}^{\prime}_{i+1}\\}=\operatorname*{arg\,min}{\mathcal{L}}$

Update library

$\mathbf{W}_{i+1},\mathbf{P}_{i+1}=\text{STSparseBayesian}(\mathbf{\Phi},\dot{\mathbf{U}}_{i+1}=\dot{\mathbf{N}}_{c}\boldsymbol{\theta}_{i+1},\mathbf{W}^{\prime}_{i+1},\mathbf{P}^{\prime}_{i+1})$

if
_$\mathcal{L}(\boldsymbol{\theta}_{i+1},\mathbf{W}_{i+1},\mathbf{B}_{i+1},\mathbf{P}_{i+1})
<\mathcal{L}^{\star}$_ then

$\mathcal{L}^{\star}=\mathcal{L}(\boldsymbol{\theta}_{i+1},\mathbf{W}_{i+1},\mathbf{B}_{i+1},\mathbf{P}_{i+1})$

$\boldsymbol{\theta}^{\star},\mathbf{W}^{\star},\mathbf{B}^{\star},\mathbf{P}^{\star}=\boldsymbol{\theta}_{i+1},\mathbf{W}_{i+1},\mathbf{B}_{i+1},\mathbf{P}_{i+1}$

else

break

end for

Stochastic Weight Averaging-Gaussian (SWAG) for posterior approximation:

$\boldsymbol{\theta}_{\text{SWA}},\mathbf{W}_{\text{SWA}},\mathbf{B}_{\text{SWA}},\mathbf{P}_{\text{SWA}}=\boldsymbol{\theta^{\star}},\mathbf{W}^{\star},\mathbf{B}^{\star},\mathbf{P}^{\star}$

With a constant learning rate for  _$i=1:T_{\rm SWAG}$_ do

SGD update
$\boldsymbol{\theta}_{i},\mathbf{W}_{i},\mathbf{B}_{i},\mathbf{P}_{i}$

$\boldsymbol{\theta}_{\text{SWA}},\mathbf{W}_{\text{SWA}},\mathbf{B}_{\text{SWA}},\mathbf{P}_{\text{SWA}}=\frac{i\boldsymbol{\theta}_{\text{SWA}}+\boldsymbol{\theta}_{i}}{i+1},\frac{i\boldsymbol{W}_{\text{SWA}}+\boldsymbol{W}_{i}}{i+1},\frac{i\boldsymbol{B}_{\text{SWA}}+\boldsymbol{B}_{i}}{i+1},\frac{i\boldsymbol{P}_{\text{SWA}}+\boldsymbol{P}_{i}}{i+1}$

end for

Compute low-rank square root of empirical covariance matrices
$\mathbf{\Lambda}_{\theta},\mathbf{\Lambda}_{W},\mathbf{\Lambda}_{B},\mathbf{\Lambda}_{P}$
from
$\\{\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{\text{SWA}}\\},\\{\boldsymbol{W}_{i}-\boldsymbol{W}_{\text{SWA}}\\},\\{\boldsymbol{B}_{i}-\boldsymbol{B}_{\text{SWA}}\\},\\{\boldsymbol{P}_{i}-\boldsymbol{P}_{\text{SWA}}\\}$

Sampling:

$\mathring{\boldsymbol{\theta}}=\boldsymbol{\theta}_{\text{SWA}}+\mathbf{\Lambda}_{\theta}\mathring{\boldsymbol{z}}_{\theta}\qquad\mathring{\boldsymbol{W}}=\boldsymbol{W}_{\text{SWA}}+\mathbf{\Lambda}_{W}\mathring{\boldsymbol{z}}_{W}\qquad\mathring{\boldsymbol{B}}=\boldsymbol{B}_{\text{SWA}}+\mathbf{\Lambda}_{B}\mathring{\boldsymbol{z}}_{B}\qquad\mathring{\boldsymbol{P}}=\boldsymbol{P}_{\text{SWA}}+\mathbf{\Lambda}_{P}\mathring{\boldsymbol{z}}_{P}$

where $\mathring{\boldsymbol{z}}_{*}$ are the random samples from
$\mathcal{N}(0,I)$.

Algorithm 1 Bayesian Alternative Direction Optimization (ADO) Learning

Input : Spline trainable parameter $\boldsymbol{\theta}$, library
$\mathbf{\Phi}(\boldsymbol{\theta})$, approximated derivative
$\dot{\mathbf{U}}$, library weight $\mathbf{W}$, and process error matrix
$\mathbf{P}$

Output : Best solution library $\mathbf{\Phi}^{\star}$,$\mathbf{W}^{\star}$,
$\mathbf{P}^{\star}$

Initialize : Threshold $\epsilon$, number of library terms $p_{\text{old}}$,
and Flag = $\mathbf{True}$

while  _Flag is True_ do

while  _not converged_ do

1\. Compute the relevance variable $\eta_{i}={q_{i}}^{2}-s_{i}$ as defined in
[19]

2\. Update library $\boldsymbol{\Phi}^{\star}$, weight $\mathbf{W}^{\star}$,
and process error matrix $\mathbf{P}^{\star}$, as shown in [19]

end while

for  _$j=1:p_{\text{old}}$_ do

$\mathbf{W}^{\star}(j)=0$ If $|\mathbf{W}^{\star}(j)|\leq\epsilon$

end for

Find the nonzero entries in $\mathbf{W}^{\star}$, record the index as
$\mathbf{I}$, update $\mathbf{\Phi}=\mathbf{\Phi}(:,\mathbf{I})$ and
$p_{\text{new}}$ = length of $\mathbf{I}$;

if  _$\mathbf{p_{new}=p_{old}}$_ then

Flag = $\mathbf{False}$

end while

Algorithm 2 Sequential Threshold Sparse Bayesian Learning

Input : ensemble number $J$, sampled weights $\\{\mathbf{W}^{j}\\}_{j=1}^{J}$,
discovered dynamical model
$\mathbf{M}(\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}{}\
;\ \mathbf{W})$, process noise covariance $\mathbf{P}$, observation model
$\mathbf{h}$, observation $\tilde{\mathbf{U}}$, observation noise covariance
$\mathbf{B}$

Output : Analysis ensemble trajectories $\mathbf{U}_{a}^{j}(t)$

Forecast :
$\mathbf{U}_{f}^{j}(t_{i+1})=\mathbf{M}(\mathbf{U}_{a}^{j}(t_{i});\mathbf{W}^{j})+\boldsymbol{\epsilon}^{j}_{2},\quad\boldsymbol{\epsilon}_{2}^{j}\sim\mathcal{(}\mathbf{0},\mathbf{P})$

$\overline{\mathbf{U}}_{f}(t_{i+1})=\frac{1}{J}\sum_{j=1}^{J}{\mathbf{U}_{f}^{j}(t_{i+1})}$

Analysis :
${\mathbf{U}}_{h}^{j}(t_{i+1})=\mathbf{h}({\mathbf{U}}_{f}^{j}(t_{i+1}))\qquad\overline{\mathbf{U}}_{h}(t_{i+1})=\frac{1}{J}\sum_{j=1}^{J}{\mathbf{U}_{h}^{j}(t_{i+1})}$
$\mathbf{C}^{fh}(t_{i+1})=\frac{1}{{J}-1}\Sigma_{j=1}^{J}(\mathbf{U}^{i}_{f}(t_{i+1})-\overline{\mathbf{U}}_{f}(t_{i+1}))(\mathbf{U}^{i}_{h}(t_{i+1})-\overline{\mathbf{U}}_{h}(t_{i+1}))^{T}$
$\mathbf{C}^{hh}(t_{i+1})=\frac{1}{{J}-1}\Sigma_{j=1}^{J}(\mathbf{U}^{i}_{h}(t_{i+1})-\overline{\mathbf{U}}_{h}(t_{i+1}))(\mathbf{U}^{i}_{h}(t_{i+1})-\overline{\mathbf{U}}_{h}(t_{i+1}))^{T}+\mathbf{B}$
$\mathbf{K}(t_{i+1})=\mathbf{C}^{fh}\bigl{(}\mathbf{C}^{hh}\bigr{)}^{-1}$

${\mathbf{U}}_{a}^{j}(t_{i+1})={\mathbf{U}}_{f}^{j}(t_{i+1})+\mathbf{K}(t_{i+1})\bigg{(}\tilde{\mathbf{U}}(t_{i+1})-{\mathbf{U}}_{h}^{j}(t_{i+1})-\boldsymbol{\epsilon}_{1}^{j}\bigg{)},\quad\boldsymbol{\epsilon}_{1}^{j}\sim\mathcal{(}\mathbf{0},\mathbf{B})$

Algorithm 3 Ensemble Kalman Filter

###  A.4 Training Details

Additional training hyper parameters used in Sec. 4 is shown in the Tab. 2.

Table 2: Training Details Case | Van der Pol | Lorenz 96 | Advection | Burgers’
---|---|---|---|---
ADO Iter | $5$ | $5$ | $1$ | $5$
ADO Epoch | $20K$ | $50K$ | $20K$ | $20K$
Post Epoch | $1K$ | $65K$ | $2K$ | $0.5K$
SWAG Epoch | $1.5K$ | $80K$ | $0.5K$ | $0.5K$
LR | $1\times 10^{-2}$ | $1\times 10^{-2}$ | $1\times 10^{-2}$ | $1\times 10^{-2}$
SWAG LR | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$ | $1\times 10^{-3}$

###  A.5 Additional Result: ODE

We list additional discovery and UQ results in this section. Fig. 5 shows 4
distributions of the coefficients for Van der Pol system in red box and Lorenz
96 system in blue box. Fig. 4 shows additional UQ result from the identified
L96 systems without incorporating the data assimilation process. The truth
trajectory is marked by red. The measurement is marked by green dots and the
ensemble trajectories are marked by blue. Although the system has been
identified with high accuracy, as shown in Tab. 3, the predicted ensembles of
the state variables still become chaotic after several seconds. It is
inevitable since the chaotic nature of the underlying system, which means any
small perturbation in any parameters would significantly influence the future
trajectories. Fortunately, the predicted covariance matrix of the Bayesian
framework makes it easy to incorporate the data assimilation with the
identified systems. With the identified distribution of system coefficients,
the data assimilation can be used to predict the future states with reduced
uncertainty, given noisy measurement data in the past. Fig. 6 shows additional
UQ result for all the 6 state variables for Lorenz 96 system, incorporating
EnKF algorithms.

Figure 4: Additional discovery results for ODE systems; the red box shows
results coefficients distribution for the Van der Pol system, and the blue box
shows the result for Lorenz 96 system. Figure 5: Additional UQ results for ODE
systems; the blue box shows the all the states prediction for Lorenz 96 system
without ensemble Kalman filter. Figure 6: Additional UQ results for ODE
systems; the blue box shows the all the states prediction for Lorenz 96 system
with ensemble Kalman filter.

###  A.6 Additional Result: PDE

In this section, we attached the qualitative result for PDE discovery and the
uncertainty quantification. The contour plot and the cross section result are
shown in Fig 7 (for advection equation and Burgers equation) and Fig. 8 (for
Burgers equation with source). The analytical form of the mentioned PDEs are
listed in Tab. 5. Probability distribution for the PDE coefficient are shown
in Fig. 9. Additional UQ prediction result for advection and Burgers’ equation
are shown in and Fig. 10.

Figure 7: 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 $x$. Figure 8:
Additional UQ results for PDE; the black box shows the cross section UQ
results for Burgers’ equation with source. Figure 9: Additional discovery
results for PDE systems; the red box shows results coefficients distribution
for the advection equation, and the blue box shows the result for Burgers’
equation.

Figure 10: Additional UQ results PDE; the red box shows the cross section UQ
results for Advection equation. The blue box shows the cross section UQ
results for Burgers’ equation.

###  A.7 Additional Discovery result

In this section, we list the full table that includes all the experiments made
for current work, as attached in Tab. 3.

Table 3: ODE and PDE discovery comparison Name |  rmse($0\%$) |  rmse($1\%$) |  rmse (large666Large noise for different cases: Van der Pol: $5\%$, Lorenz 96: $10\%$, Advection: $20\%$, Burgers: $10\%$, Burgers’ with source: $20\%$, Heat: $15\%$, Poisson: $5\%$) | $\mathbf{M_{P}}$ |  $\mathbf{M_{R}}$777$\mathbf{M_{P}},\mathbf{M_{R}}$ are only reported for the largest noise cases |  Training Cost888All cases are running on a Nvidia 2070 Ti GPU card
---|---|---|---|---|---|---
Van der Pol Oscillator
BSL(Ours) | $\mathbf{0.2}$ | $2.82$ | $\mathbf{18.04}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 133(+3)s$
PINN-SR |  Fail999Fail means failure in discovery of the parsimonious ODE/PDE forms. | Fail | Fail | $0.214$ | $0.75$ | $\sim 1213s$
SINDy | $1.0$ | $\mathbf{1.93}$ | Fail | $0.267$ | $1.0$ | $\sim 10s$
RVM | $1.0$ | $2.54$ | $27.46$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 10s$
Lorenz 96
BSL(Ours) | $\mathbf{0.269}$ | $1.47$ | $\mathbf{13.0}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 1654(+438)s$
PINN-SR | Fail | Fail | Fail | $0.5$ | $0.22$ | $\sim 10788s$
SINDy | $0.4$ | $0.64$ | Fail | $0.75$ | $\mathbf{1}$ | $\sim 10s$
RVM | $0.4$ | $\mathbf{0.6}$ | $49.7$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 25s$
Advection Equation
BSL(Ours) | $\mathbf{0.26}$ | $\mathbf{1}$ | $\mathbf{1.9}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 946(+233)s$
PINN-SR | $5.9$ | $4.5$ | $30.4$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 650s$
SINDy | $2.3$ | $8.2$ | $38.9$ | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 10s$
RVM | $0.77$ | $6.76$ | Fail | ${0.2}$ | ${\mathbf{1}}$ | $\sim 4s$
Burgers’ Equation
BSL(Ours) | $3.62$ | $4.13$ | $\mathbf{6.38}$ | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 117(+74)s$
PINN-SR | $10.2$ | $\mathbf{3.3}$ | $10.3$ | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 512s$
SINDy | $0.826$ | Fail | Fail | ${\mathbf{1}}$ | ${0.5}$ | $\sim 10s$
RVM | $\mathbf{0.754}$ | Fail | Fail | ${0.1429}$ | ${0.5}$ | $\sim 4s$
Name |  rmse($0\%$) |  rmse($0.1\%$) |  rmse (large) | $\mathbf{M_{P}}$ | $\mathbf{M_{R}}$ | Training Cost
Burgers’ with Source
BSL(Ours) | $11$ | $\mathbf{12.4}$ | $\mathbf{13.4}$ | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 396(+340)s$
PINN-SR | $\mathbf{10.5}$ | $15$ | $34.6$ | ${1}$ | ${1}$ | $\sim 600s$
SINDy | $26.2$ | Fail | Fail | ${\mathbf{1}}$ | ${0.33}$ | $\sim 10s$
RVM | $27.6$ | Fail | Fail | ${0.5}$ | ${0.67}$ | $\sim 10s$
Heat Equation
BSL(Ours) | $19$ | $19$ | $\mathbf{38.9}$ | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 71(+8)s$
PINN-SR | Fail | Fail | Fail | $0$ | $0$ | $\sim 285s$
SINDy | $\mathbf{1.9}$ | $\mathbf{17}$ | Fail | ${0.25}$ | ${1}$ | $\sim 10s$
RVM | $2.8$ | $6$ | Fail | ${0}$ | ${0}$ | $\sim 6s$
Poisson Equation
BSL(Ours) | $\mathbf{1.18\times 10^{-2}}$ | $\mathbf{0.133}$ | $\mathbf{16.7}$ | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 92(+9)s$
PINN-SR | Fail | Fail | Fail | ${0}$ | ${0}$ | $\sim 3737s$
SINDy | $1.15$ | 87 | 962 | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 10s$
RVM | $1.15$ | $232$ | $968$ | ${\mathbf{1}}$ | 1 | $\sim 10s$

###  A.8 Analytical forms of the discovered system

Table 4: Analytical forms of ODE Name |  $\epsilon$ (large)
---|---
Vander Pol Oscillator
True |  $\frac{dx}{dt}=y$,$\frac{dy}{dt}=-x-0.5x^{2}y+0.5y$
BSL(Ours) |  $\frac{dx}{dt}=1.0096(\pm 0.024)y$,$\frac{dy}{dt}=-0.9858(\pm 0.037)x-0.4801(\pm 0.114)x^{2}y+0.4889(\pm 0.111)y$
PINN-SR | $\frac{dx}{dt}=1.3079-0.1151x+0.9982y+0.1939x^{2}-0.4101xy+0.8559y^{2}-0.3577x^{3}$
|  $+0.6764x^{2}y-0.1207xy^{2}+0.6261y^{3}$,
| $\frac{dy}{dt}=-0.6718x+1.6035y+1.6339y^{2}+0.4803y^{3}$
SINDy | $\frac{dx}{dt}=0.1197+0.169x+0.9975y-0.0292x^{2}+0.0232xy-0.0274y^{2}$
|  $-0.0393x^{3}+0.0208x^{2}y-0.0463xy^{2}$,
|
$\frac{dy}{dt}=-1.0909x+0.149y+0.0294x^{3}-0.4x^{2}y-0.0356xy^{2}+0.0762y^{3}$
RVM |  $\frac{dx}{dt}=0.9957(\pm 0.073)y$,$\frac{dy}{dt}=-0.9943(\pm 0.0597)x-0.4777(\pm 0.0809)x^{2}y+0.4714(\pm 0.1204)y$
Lorenz 96
True |  $\frac{dX_{1}}{dt}=(X_{2}-X_{5})X_{6}-X_{1}+8$, $\frac{dX_{2}}{dt}=(X_{3}-X_{6})X_{1}-X_{2}+8$,
|  $\frac{dX_{3}}{dt}=(X_{4}-X_{1})X_{2}-X_{3}+8$,
$\frac{dX_{4}}{dt}=(X_{5}-X_{2})X_{3}-X_{4}+8$,
|  $\frac{dX_{5}}{dt}=(X_{6}-X_{3})X_{4}-X_{5}+8$,
$\frac{dX_{6}}{dt}=(X_{1}-X_{4})X_{5}-X_{6}+8$
BSL(Ours) |  $\frac{dX_{1}}{dt}=1.0033(\pm 2.37\times 10^{-4})X_{2}X_{6}-0.9926(\pm 2.03\times 10^{-4})X_{5}X_{6}$,
|  $-0.9991(\pm 9.04\times 10^{-4})X_{1}+7.9773(\pm 2.08\times 10^{-2})$,
| $\frac{dX_{2}}{dt}=0.9963(\pm 5.54\times 10^{-5})X_{1}X_{3}-0.9942(\pm
2.47\times 10^{-4})X_{1}X_{6}$
|  $-1.0054(\pm 2.31\times 10^{-4})X_{2}+7.8719(\pm 7.05\times 10^{-4})$,
| $\frac{dX_{3}}{dt}=1.0106(\pm{2.19\times
10^{-4}})X_{2}X_{4}-1.0029(\pm{2.03\times 10^{-4}})X_{1}X_{2}$
|  $-0.9979(\pm 5.3\times 10^{-4})X_{3}+7.9938(\pm{1.48\times 10^{-2}})$,
| $\frac{dX_{4}}{dt}=1.0067(\pm 1.01\times 10^{-4})X_{3}X_{5}-1.0103(\pm
9.09\times 10^{-5})X_{2}X_{3}$
|  $-0.9926(\pm 2.36\times 10^{-4})X_{4}+8.055(\pm{2.41\times 10^{-3}})$,
| $\frac{dX_{5}}{dt}=0.9953\pm(4.39\times
10^{-4})X_{4}X_{6}-0.9922\pm(3.88\times 10^{-4})X_{3}X_{4}$
|  $-1.0072\pm(4.1\times 10^{-4})X_{5}+7.8538(\pm{2.1\times 10^{-3}})$,
| $\frac{dX_{6}}{dt}=1.0095(\pm 2.71\times 10^{-4})X_{1}X_{5}-0.9967(\pm
2.83\times 10^{-4})X_{4}X_{5}$
| $-0.9772(\pm 2.46\times 10^{-4})X_{6}+7.9964(\pm 4.15\times 10^{-2})$
PINN-SR | N/A
SINDy | $\frac{dX_{1}}{dt}=8.0985-0.9423X_{1}-0.1007X_{4}-0.1659X_{6}+0.9582X_{2}X_{6}-0.944X_{5}X_{6}$
|
$\frac{dX_{2}}{dt}=7.9372-0.1268X_{1}-0.9607X_{2}+0.9636X_{1}X_{3}-0.9558X_{1}X_{6}$
|
$\frac{dX_{3}}{dt}=8.1523-0.0983X_{1}-0.1652X_{2}-0.9952X_{3}-0.9541X_{1}X_{2}+0.9651X_{2}X_{4}$
|
$\frac{dX_{4}}{dt}=7.4958-0.9754X_{4}+0.1515X_{5}-0.973X_{2}X_{3}+0.9206X_{3}X_{5}$,
|
$\frac{dX_{5}}{dt}=7.8556-0.1146X_{4}-0.9457X_{5}-0.9559X_{3}X_{4}+1.0023X_{4}X_{6}$
|
$\frac{dX_{6}}{dt}=7.6026+0.0934X_{2}-0.927X_{6}+0.9711X_{1}X_{5}-0.9754X_{4}X_{5}$
RVM |  $\frac{dX_{1}}{dt}=0.9613(\pm 1.21)X_{2}X_{6}-0.962(\pm 1.27)X_{5}X_{6}$,
|  $-0.8983(\pm 2.67)X_{1}+7.518(\pm 5.88)$,
| $\frac{dX_{2}}{dt}=0.9578(\pm 1.31)X_{1}X_{3}-0.9693(\pm 1.24)X_{1}X_{6}$
|  $-0.9217(\pm 2.8)X_{2}+7.6442(\pm 6.1575)$,
| $\frac{dX_{3}}{dt}=0.9535(\pm{1.48})X_{2}X_{4}-0.9803(\pm{1.38})X_{1}X_{2}$
|  $-0.9468(\pm 2.96)X_{3}+7.5323(\pm 6.56)$,
| $\frac{dX_{4}}{dt}=0.9365(\pm 1.24)X_{3}X_{5}-0.9748(\pm 1.15)X_{2}X_{3}$
|  $-0.9114(\pm 2.59)X_{4}+7.67(\pm 5.49)$,
| $\frac{dX_{5}}{dt}=0.9981\pm(1.26)X_{4}X_{6}-0.9667\pm(1.2)X_{3}X_{4}$
|  $-0.9216\pm(2.75)X_{5}+7.6146(\pm 5.96)$,
| $\frac{dX_{6}}{dt}=0.9677(\pm 1.42)X_{1}X_{5}-0.9757(\pm 1.38)X_{4}X_{5}$
| $-0.8986(\pm 2.86)X_{6}+7.6657(\pm 6.57)$

Table 5: Analytical forms of unsteady PDE Name |  $\epsilon$ (large)
---|---
Advection Equation
True | $u_{t}=-u_{x}$
BSL(Ours) | $u_{t}=-0.9988(\pm 0.024)u_{x}$
PINN-SR | $u_{t}=-0.997u_{x}$
SINDy | $u_{t}=-0.9961u_{x}$
RVM | $u_{t}=-0.5148(\pm 0.106)u_{x}-0.9797(\pm 0.384)uu_{x}$
| $-0.018(\pm 0.208)u^{2}u_{x}-0.075(\pm 0.148)u^{2}u_{xxx}+0.049(\pm
0.111)u^{3}u_{xxx}$
Burgers’ Equation
True | $u_{t}=-uu_{x}+0.5u_{xx}$
BSL(Ours) | $u_{t}=-0.9929(\pm 0.086)uu_{x}+0.4993(\pm 0.005)u_{xx}$
PINN-SR | $u_{t}=-1.0103uu_{x}+0.5051u_{xx}$
SINDy | $u_{t}=-0.8179uu_{x}$
RVM | $-0.0809(\pm 0.041)u_{x}-1.6684(\pm 0.2618)uu_{x}+4.1835(\pm 0.571)u^{2}u_{x}$
| $-3.9068(\pm 0.391)u^{3}u_{x}+0.1916(\pm 0.064)uu_{xx}$
| $-1.0314(\pm 0.197)u^{2}u_{xx}+1.5504(0.156)u^{3}u_{xx}$
Burgers’ Equation with Source
True | $u_{t}=-uu_{x}+0.1u_{xx}+sin(x)sin(t)$
BSL(Ours) | $-0.9882(\pm 0.246)uu_{x}+0.105(\pm 0.022)u_{xx}+0.9859(\pm 0.005)sin(x)sin(t)$
PINN-SR | $u_{t}=-0.9576uu_{x}+0.1168u_{xx}+1.0179sin(x)sin(t)$
SINDy | $u_{t}=0.8052sin(x)sin(t)$
RVM | $-0.0234(\pm 0.145)uu_{x}+0.8318(\pm 0.142)sin(x)sin(t)$
| $-0.0789(\pm 0.105)sin(x)+0.3558(\pm 0.156)sin(x)cos(t)$
Table 6: Analytical forms of steady PDE Name |  $\epsilon$ (large)
---|---
Heat Equation
True | $u_{yy}=-u_{xx}$
BSL(Ours) | $u_{yy}=-0.9611(\pm 0.059)u_{xx}$
PINN-SR | $u_{yy}=0.5544u_{x}$
SINDy | $u_{yy}=-0.069u_{xx}+13.8988uu_{x}-19.493u_{x}+0.2468uu_{xx}$
RVM | $u_{yy}=-7.5861(\pm 202.6)u_{x}$
Poisson Equation
True | $u_{yy}=-u_{xx}-sin(x)sin(y)$
BSL(Ours) | $u_{yy}=-0.9788(\pm 8.75\times 10^{-4})u_{xx}-0.9897(\pm 4.63\times 10^{-4})sin(x)sin(y)$
PINN-SR | $u_{yy}=0.13611752uu_{x}+0.29748484uu_{xxx}$
SINDy | $u_{yy}=0.2316u_{xx}-0.4221sin(x)sin(y)$
RVM | $u_{yy}=0.2284(\pm 0.007)u_{xx}-0.3954(\pm 0.01)sin(x)sin(y)$

###  A.9 Implementation detail for the spline

In the current work, we apply direct tensor product to extend spline for
solving spatial-temporal field and the relevant statistics are attached in
Tab. 7, where the numerber of control points (trainable weights) $\theta$ in
1-d is marked by red and the total number of control points for the 2-d
scenario is listed in the last column. We only store the non-zero elements for
two-dimensional basis to leverage the sparsity (local support) of the spline.
However, we must claim that it is _not_ an optimal way to extend spline for
higher spatial dimensions. In that case, a spline kernel can be defined and
the tensor-product is only processed in the subdomain, as shown in [43].

Table 7: Direct tensor-product spline for PDE basis $x$ |  basis $t(y)$ |  basis$(x,t(y))$ | sparsity | trainable params
---|---|---|---|---
Advection Equation
$50\times{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{54}}$ | $50\times{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}54}$ | $62500$ | $0.0086$ | $2916$
Burgers’ Equation
$128\times{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}13}$ | $101\times{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}19}$ | $92872$ | $0.058$ | $247$
Burgers’ Equation with Source
$201\times{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}103}$ | $101\times{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}103}$ | $251415$ | $0.0012$ | $10609$
Heat Equation
$51{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\times 11}$ | $51\times{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}11}$ | $41209$ | $0.1309$ | $121$
Poisson
$101\times{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}53}$ | $101\times{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}53}$ | $162409$ | $0.0057$ | $2809$

###  A.10 Relevant Terminologies

Table 8: Terminologies Name | Symbol | Explanation
---|---|---
Control points | $\mathbf{\theta}$ |  Trainable weight $\theta$ for spline basis.
Knots | $\mathbf{\tau}_{s}$ | Location of control points.
Measurement points |  | Sparse spatio-temporal points with labels.
Collocation points |  | Dense spatio-temporal points without labels.
| $\mathbf{N_{m}}$ | Spline basis evaluated at measurement points.
| $\mathbf{N_{c}}$ | Spline basis evaluated at collocation points.
Library candidates | $\mathbf{\Phi}$ | A collection of polynomial terms that
|  | the system identification algorithm can choose parsimonious terms from it.
|  |  e.g., $\\{x,y,x^{2}y,...\\}$ (for ODE) or $\\{u,uu_{x},u_{xx}...\\}$ for (PDE).
ADO iteration |  | Alternating direction optimization to update the trainable parameters
|  |  including control points $\mathbb{\theta}$, weight of library candidates $\mathbf{W}$ and covariance matrices.
Aleatoric Uncertainty |  | Due to intrinsic randomness by nature, which is irreducible.
Epistemic Uncertainty |  | Because of a lack of knowledge, which can be reduced by adding more information.

###  A.11 Real world application: predator-prey system

In this section, we would test our proposed BSL model on one real-world case,
predator-prey system. The real data set is obtained online and it depicts the
population of hares and lynx from 1900 to 1920 from Hudson Bay Company. The
data is presented in Tab. 11: The reference governing equation by mathematical
analysis is:

$\displaystyle\frac{dx}{dt}$ $\displaystyle=0.4807x-0.0248xy$ (28)
$\displaystyle\frac{dy}{dt}$ $\displaystyle=-0.9272y+0.0276xy$ (29)

We test the 4 methods on this data set and the result can be found in Tab 9
.We have made assumptions about constructing the libraries. We assume the
predator (lynx) only feeds on the prey (hares). Meantime, the prey (hares)
only has one predator (lynx). Therefore, the change rate of these two species
can only depend on themselves ($x$,$y$) and some higher order correlations
between them ($xy,x^{2}y,xy^{2}$). The discovered forms of the 4 methods are
listed in Tab. 10. Finally, the UQ results are shown in Fig 11. In short, only
the proposed method work on the real sparse and noisy dataset. And the UQ
prediction covers more measurement points than the reference model Eq. 28,
which helps to better explain the data set.

Table 9: ODE discovery comparison Name |  rmse (large) | $\mathbf{M_{P}}$ | $\mathbf{M_{R}}$ | Training Cost
---|---|---|---|---
Predator-prey
BSL(Ours) | ${\mathbf{30.4}}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 2278(+6)s$
PINN-SR | Fail | $0.5$ | $0.25$ | $\sim 2988s$
SINDy | Fail | $0.6$ | $0.75$ | $\sim 10s$
RVM | Fail | $0.8$ | $1$ | $\sim 10s$
Table 10: Analytical forms of ODE Name |  $\epsilon$ (large)
---|---
Predator-prey (Lotka-Volterra)
True |  $\frac{dx}{dt}=0.4807x-0.0248xy$,$\frac{dy}{dt}=-0.9272y+0.0276xy$
BSL(Ours) | $\frac{dx}{dt}=-0.5124(\pm 0.028)x-0.0266(\pm 8.72\times 10^{-4})xy$
| $\frac{dy}{dt}=-0.9258(\pm 0.065)y+0.0279(\pm 1.47\times 10^{-3})xy$
PINN-SR | $\frac{dx}{dt}=-13.9238y$
| $\frac{dy}{dt}=-0.1144y$
SINDy |  $\frac{dx}{dt}=0.5813x-0.0261xy$,
| $\frac{dy}{dt}=0.2549x-0.2702y$
RVM |  $\frac{dx}{dt}=0.5732(\pm 0.6488)x-0.2386(\pm 0.4643)y-0.0253(\pm 0.1432)xy$,
| $\frac{dy}{dt}=-0.8018(\pm 0.9459)y+0.0226(\pm 0.1481)xy$
Table 11: Lynx-Hares population Year |  Hares($\times 1000$) |  Lynx($\times 1000$)
---|---|---
1900 | 30 | 4
1901 | 47.2 | 6.1
1902 | 70.2 | 9.8
1903 | 77.4 | 35.2
1904 | 36.3 | 59.4
1905 | 20.6 | 41.7
1906 | 18.1 | 19
1907 | 21.4 | 13
1908 | 22 | 8.3
1909 | 25.4 | 9.1
1910 | 27.1 | 7.4
1911 | 40.3 | 8
1912 | 57 | 12.3
1913 | 76.6 | 19.5
1914 | 52.3 | 45.7
1915 | 19.5 | 51.1
1916 | 11.2 | 29.7
1917 | 7.6 | 15.8
1918 | 14.6 | 9.7
1919 | 16.2 | 10.1
1920 | 24.7 | 8.6
Figure 11: UQ results for predator-prey system; left sub-fig shows the result
for state $x$,the hares population; the right sub-fig shows the result for
state $y$, the lynx population. The green dot is the sparse and biased
measurement data. The red line is the reference model and the blue curves are
the ensemble predictions from our proposed model.

###  A.12 Experiment on different smoothing algorithms

This section we will test the effect of different smoothing algorithms and its
impact to the SINDy result. These methods include smoothing based polynomial
interpolation, convolutional smoother, smoothing with Tikhonov regularization
and smoothing with spline fitting. We performed comparison of SINDy with these
smoothing methods on a representative ODE system (Van der Pol system) and PDE
system (Poisson equation) studied in this work. As shown in the Table 12 and
Table 13 below, when the data noise is above $5\%$, although data is
preprocessed using smoothing and uniform resampling, none of these method
work. Basically, the SINDy still failed to discover the correct model forms
with different smoothing schemes, and the identified systems are different
from the true. In contrast, our proposed approach is very robust and superior
to handling corrupted data, thanks to the spline learning in Bayesian
settings.

Table 12: ODE and PDE discovery comparison Name | $\textbf{rmse}(\epsilon=5\%)$ | $\mathbf{M_{P}}$ | $\mathbf{M_{R}}$ | Training Cost
---|---|---|---|---
Van der Pol Oscillator
BSL(Ours) | $\mathbf{18.04}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\sim 133(+3)s$
SINDy(No smoother) | Fail | $0.33$ | $\mathbf{1}$ | $\sim 10s$
SINDy(Poly) | Fail | $0.6$ | $0.75$ | $\sim 10s$
SINDy(Conv) | Fail | $0.4$ | $\mathbf{1}$ | $\sim 10s$
SINDy(Tikhonov) | Fail | $0.21$ | $\mathbf{1}$ | $\sim 10s$
SINDy(Spline) | Fail | $0.25$ | $1$ | $\sim 10s$
Poisson Equation
BSL(Ours) | $\mathbf{16.7}$ | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 92(+9)s$
SINDy(No smoother) | Fail | $0.5$ | $\mathbf{1}$ | $\sim 10s$
SINDy(Poly) | $962$ | ${\mathbf{1}}$ | ${\mathbf{1}}$ | $\sim 10s$
SINDy(Conv) | Fail | $0.66$ | $\mathbf{1}$ | $\sim 10s$
SINDy(Tikhonov) | Fail | $0.5$ | $\mathbf{1}$ | $\sim 10s$
SINDy(Spline) | Fail | $0.5$ | $\mathbf{1}$ | $\sim 10s$
Table 13: Smoothing algorithm effects Name | Analytical form
---|---
Van der Pol Oscillator
True | $\frac{dx}{dt}=y$
| $\frac{dy}{dt}=-x-0.5x^{2}y+0.5y$
SINDy(No smoother) | $\frac{dx}{dt}=0.0277-0.217x+1.4y-0.025x^{2}$
|  $+0.0542x^{3}-0.1285x^{2}y+0.1039xy^{2}-0.0896y^{3}$,
| $\frac{dy}{dt}=-0.9835x+0.3462y-0.4675x^{2}y+0.0325y^{3}$
SINDy(Poly) | $\frac{dx}{dt}=0.1197+0.169x+0.9975y-0.0292x^{2}$
|  $+0.0232xy-0.0274y^{2}-0.0393x^{3}+0.0208x^{2}y-0.0463xy^{2}$,
|
$\frac{dy}{dt}=-1.0909x+0.149y+0.0294x^{3}-0.4x^{2}y-0.0356xy^{2}+0.0762y^{3}$
SINDy(Conv) |  $\frac{dx}{dt}=0.1518+0.9978y-0.0486x^{2}+0.0274xy-0.0327y^{2}$,
| $\frac{dy}{dt}=-1.1154x-0.2143y+0.0348x^{3}-0.4374x^{2}y+0.0613y^{3}$
SINDy(Tikhonov) |  $\frac{dx}{dt}$ $=-0.2572+0.1486x+1.1299y+0.0982x^{2}-0.0721xy$,
| $+0.0549y^{2}-0.0478x^{2}y-0.0816xy^{2}-0.0542y^{3}$
| $\frac{dy}{dt}=0.2338-1.3282x+0.1603y-0.0718x^{2}$
| $+0.058xy-0.0564y^{2}+0.1069x^{3}-0.0842x^{2}y+0.1373xy^{2}-0.0415y^{3}$
SINDy(Spline) | $\frac{dx}{dt}=0.1916+1.5304y-0.0614x^{2}+0.0215xy$
| $-0.0325y^{2}-0.155x^{2}y+0.0796xy^{2}-0.1168y^{3}$
| $\frac{dy}{dt}=-0.35-1.0127x+0.6756y+0.0695x^{2}$
| $+0.0741y^{2}-0.5493x^{2}y+0.0523xy^{2}-0.0412y^{3}$
Poisson equation
True | $u_{yy}=-u_{xx}-sin(x)sin(y)$
SINDy(No smoother) | $u_{yy}=0.5635u_{xx}+1.683uu_{x}+0.088uu_{xx}-0.17sin(x)sin(y)$
SINDy(Poly) | $u_{yy}=0.2316u_{xx}-0.4221sin(x)sin(y)$
SINDy(Conv) | $u_{yy}=0.47u_{xx}-0.2649sin(x)sin(y)+0.073sin(x)cos(y)$
SINDy(Tikhonov) | $u_{yy}=-0.17u_{xx}-0.087uu_{x}-0.458sin(x)sin(y)+0.02sin(x)cos(y)$
SINDy(Spline) | $u_{yy}=0.1562u_{xx}+0.1751uu_{x}-0.3951sin(x)sin(y)-0.0911sin(x)cos(y)$