# Bayesian Dynamical System Identification with Unified Sparsity Priors and
Model Uncertainty

Prem Ratan Mohan Ram
Institut für Dynamik und Schwingungen
Technische Universität Braunschweig
38106 Braunschweig, Germany
<EMAIL_ADDRESS>
&Ulrich Römer
Institut für Dynamik und Schwingungen
Technische Universität Braunschweig
38106 Braunschweig, Germany
<EMAIL_ADDRESS>
&Richard Semaan
Institute of fluid mechanics
Technische Universität Braunschweig
38106 Braunschweig, Germany
<EMAIL_ADDRESS>

###### Abstract

This work is concerned with uncertainty quantification in dynamical system
identification. Dynamical systems are ubiquitous in design and control
applications and recent efforts focus on their data-driven construction. Our
starting point is the sparse-identification of nonlinear dynamics (SINDy)
framework, which reformulates system identification as a regression problem,
where unknown functions are approximated from a sparse subset of an underlying
library. In this manuscript, we formulate this system identification method in
a Bayesian framework to handle parameter and structural model uncertainties.
We present a general approach to enforce sparsity, which builds on the
recently introduced class of neuronized priors. We perform comparisons between
different variants such as Lasso, horseshoe, and spike and slab priors, which
are all obtained by modifying a single activation function. We also outline
how state observation noise can be incorporated with a probabilistic state-
space model. The resulting Bayesian regression framework is robust and simple
to implement. We apply the method to two generic numerical applications, the
pendulum and the Lorenz system, and one aerodynamic application employing
experimental measurements.

##  1 Introduction

Developing reliable dynamical models is crucial to all scientific disciplines
ranging from epidemiology, neuroscience, finance, to turbulence. Dynamic
modeling for the long-term features is a key enabler for physical
understanding, state estimation from limited sensors signals, prediction,
control, and optimization. Dynamical system modeling has seen tremendous
progress in the last decades, driven by algorithmic advances, accessibility to
large data, and hardware speedups. One breakthrough in system identification
was reported by Bongard and Lipson [1] using symbolic regression. The method
performs a heuristic search of the best equation that describes the dynamics.
Symbolic regression is however expensive and not easily scalable to large
systems. This limitation may be bypassed by black-box techniques. These
include Volterra series [2], autoregressive models [3] (e.g., ARX, ARMA, and
NARMAX), eigensystem realization algorithm (ERA) [4], and neural network (NN)
models [5]. These approaches, however, have limited interpretability and
provide little physical insights. Some (e.g. NN) require large volumes of data
and long training time, luxuries that are not always at hand.

Recently, sparse identification of nonlinear dynamical systems (SINDy) was
introduced [6]. The approach relies on sparse regression for system
identification from time-series data and shows remarkable performance. The
regression allows the determination of the system coefficients, which multiply
a set of candidate basis functions. Sparsity is achieved through a
regularization step that reduces complexity and the risk of overfitting.
However, both the regularization coefficient and the regularization norm are
usually selected in a heuristic or ad hoc manner, without a general rigorous
approach. Moreover, the deterministic approach of SINDy does not consider the
uncertainty in the inferred model parameters nor in the library of basis
functions.

This shortcoming can be addressed by embedding the SINDy algorithm into a
Bayesian framework, which is a subject of current great interest. In [7], the
authors discuss the analogy between the SINDy formulation and a Bayesian MAP
estimate. A Gaussian prior is used in a variational Bayesian approach to
quantify uncertainties in the model. The approach directly employs derivative
data, whereas the state is assumed to be perfectly known. A similar setting
has been considered in [8], where a Bayesian thresholding algorithm was
employed to obtain a SINDy variant with quantified uncertainty. The use of
standard Gaussian priors in both cases simplifies the algorithmic treatment,
however, it has limited capabilities in enforcing sparsity. Additionally, such
a prior is mainly suited to handle parametric uncertainty in the dynamical
model. Here, we emphasize the concept of model uncertainty, understood as the
uncertainty in the structure of the identified model. When considering model
uncertainty, a possible approach is to introduce a latent binary vector
$\gamma_{i}\in\\{0,1\\}$, which indicates exclusion/inclusion of the
associated basis function. Bayesian inference then allows to estimate the
exclusion/inclusion probabilities $P(\gamma_{i}=0,1)$ together with the
marginal distributions of the parameter. The subsequent prediction can then be
based on model averaging, where multiple models are considered according to
their probabilities. Another possibility uses a single model, for instance,
the highest posterior probability model or the median model [9]. The
difficulty in using a model indicator vector is the exponential growth with
the size of the library, which results in a complicated algorithmic treatment
for complex models. This aspect of model uncertainty has been discussed in the
context of the Bayesian Lasso in [10].

In a recent contribution [11], spike and slab priors, which are among the most
efficient sparsity priors in Bayesian regression, have been introduced in the
dynamical system learning context. The authors derive sparse models relying on
median model selection and demonstrate an improved sparsity in the inferred
model compared to the relevance vector machine. Here, again noise has been
considered in the model equation, while the system states were assumed to be
perfectly known. A nonlinear Bayesian learning approach, which is however not
related to the SINDy framework, was put forth in [12] employing automatic
relevance determination priors. Finally, a general formulation of Bayesian
identification of dynamical system was recently presented in [13]. There, the
authors formulate a hidden Markov model of a time-discrete system, including
process, observation and parameter uncertainty. The SINDy framework and
dynamic mode decomposition are recovered as MAP estimates under specific
choices for measurement and process noise, as well as the observation
operator. Based on nonlinear Kalman approximations of the marginal Likelihood
(marginalized over the inferred hidden states), the parameter posterior
distribution is obtained with a dedicated MCMC algorithm.

In this study, we propose a flexible and automatable Bayesian framework for
sparse dynamical system identification. The Bayesian formulation provides a
rigorous approach for the choices of the residual and regularization terms in
the sparse system regression framework. The method is enabled by the recently-
proposed neuronized priors [14] that provide a unified formulation. In
particular, by choosing a suitable activation function, neuronized priors can
represent highly efficient shrinkage priors such as discrete spike and slab,
Lasso, and Horseshoe priors. Hence, we can explicitly target model selection
and model uncertainty in addition to quantifying the uncertainty in model
parameters. We outline how this flexible prior framework can be combined with
the general setting of [13], which accounts for both process and observation
noise. Considering process noise provides a natural means to account for
numerical differentiation errors, which need to be taken into account if no
state derivative data are available. Including observation noise, in turn,
seems to be a natural, yet currently often omitted, step in the SINDy
framework. Our approach is compared against the least-squares regression with
thresholding and applied on two generic problems and one experimental dataset
using different shrinkage priors and different error minimization functions.
The proposed method is accurate and flexible.

The manuscript is structured as follows. The original SINDy method is recalled
in Section 2. The proposed stochastic formulation with neuronized priors is
introduced in Section 3 and applied on three different applications in Section
4. The study is summarized in Section 5.

##  2 Sparse system identification

Before we introduce the Bayesian formulation, it is helpful to first review
the deterministic SINDy algorithm. Following [6], we consider dynamical
systems of the form

$\dot{\mathbf{x}}=\frac{d\mathbf{x}}{dt}=\mathbf{f}(\mathbf{x}(t)),\quad\mathbf{x}(0)=\mathbf{x}_{0},$

where $\mathbf{x}\in\mathbb{R}^{n}$ refers to the state of the system with
initial condition $\mathbf{x}_{0}$. Our aim is to learn the function
$\mathbf{f}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ from time-series data
assuming a set of basis functions. The SINDy paradigm formulates this learning
task as a parameter estimation problem. Without loss of generality, let us
introduce the problem setting with a simple one-dimensional state $x$.
Possible library function candidates for $\mathbf{f}$ could be a linear and
quadratic basis, hence,

$\dot{x}(t)\approx\xi_{1}x(t)+\xi_{2}x^{2}(t),\quad x(0)=x_{0}.$ (1)

The parameters $\bm{\mathbf{\xi}}=(\xi_{1},\xi_{2})^{\top}$ are estimated as
the solution of a regression problem. In vector notation, (1) can be expressed
as

$\dot{\mathbf{x}}\approx\underbrace{[\mathbf{x}\
\mathbf{x}^{2}]}_{=\Theta(\mathbf{x})}\bm{\mathbf{\xi}}\,,$ (2)

where the basis functions $\mathbf{x}=(x(t_{1}),\ldots,x(t_{m}))^{\top}$ and
$\mathbf{x}^{2}=(x^{2}(t_{1}),\ldots,x^{2}(t_{m}))^{\top}$ are vectors
evaluated over discrete time $t_{i},i=1,\ldots,m$. The regression problem
reads

$\displaystyle\bm{\mathbf{\xi}}^{*}$
$\displaystyle=\mathrm{argmin}_{\bm{\mathbf{\xi}}}\sum_{i=1}^{m}(\dot{x}(t_{i})-\xi_{1}x(t_{i})-\xi_{2}x^{2}(t_{i}))^{2}$
(3)
$\displaystyle=\mathrm{argmin}_{\bm{\mathbf{\xi}}}\|\dot{\mathbf{x}}-\Theta(\mathbf{x})\bm{\mathbf{\xi}}\|^{2}_{2},$
(4)

where $\|\cdot\|_{2}$ refers to the Eucledian norm. Crucial for the success of
this approach is the choice of a library of basis functions. In particular, we
assume that the basis library is sufficiently rich to capture all relevant
dynamics. Note that since the derivative is rarely available, a common
approach employs a numerical approximation, i.e., a difference quotient or
total variation denoising. After identification, the dynamics of each
component can be expressed as

$\dot{x}\approx\Theta(x(t))\bm{\mathbf{\xi}},$ (5)

where $\Theta(x(t))$ contains the library of functions, as given for instance
in (1). To keep the notation simple, we maintain the case of a one-dimensional
state variable for the time being and provide details on the general case only
later. We note that also in the general case with higher state dimensions,
each elements of the vector $\mathbf{f}$ are learned independently, resulting
in an individual parameter vector $\bm{\mathbf{\xi}}$ for each dynamical
equation.

So far, we have assumed that both the state and the time derivative of the
state are known, which is uncommon in real-life applications. With a suitable
approximation, such as a forward Euler different quotient, we can rewrite (5)
as

$\frac{x(t_{l})-x(t_{l-1})}{\Delta
t}\approx\Theta(x(t_{l-1}))\bm{\mathbf{\xi}},\quad l=2,\ldots,m,$ (6)

where we have assumed a uniform time step size $\Delta t$, for simplicity. The
regression problem is now expressed as,

$\bm{\mathbf{\xi}}^{*}=\mathrm{argmin}_{\bm{\mathbf{\xi}}}\sum_{l=2}^{m}\left(\frac{x(t_{l})-x(t_{l-1})}{\Delta
t}-\Theta(x(t_{l-1}))\bm{\mathbf{\xi}}\right)^{2}.$ (7)

This approach assumes that the state is perfectly known, i.e., observation
error is neglected, which is a potential source of bias. Consider for instance
a data model of the form

$y(t)=x(t)+e(t),$ (8)

which expresses that the true state $x$ cannot be observed, instead we only
have access to the surrogate (sometimes referred to as the observable) $y$. In
a general parameter estimation setting, it has been reported [15] that
employing $y$ instead of $x$ results in a biased estimate. A possible remedy
is using a pseudo-least-squares estimate, where first estimators of the
solution and the derivative are introduced as $\hat{x},d\hat{x}/dt$, which are
then used to solve the regression problem. The properties of the pseudo-least
squares estimator, as well as the connection to error in measurement models
have been investigated in [15].

We proceed by introducing a matrix vector version of (6), by setting
$z_{l-1}=(x(t_{l})-x(t_{l-1}))/\Delta t$,
$\mathbf{z}=(z_{1},\ldots,z_{m-1})^{\top}$ and the $\mathbb{R}^{m-1\times p}$
matrix $\mathbf{D}$, where the row $j$ of $\mathbf{D}$ is given by
$\Theta(\mathbf{x}(t_{j}))$. Then, the regression model can be expressed as

$\mathbf{z}\approx\mathbf{D}\bm{\mathbf{\xi}},$ (9)

which is a standard linear model. One main assumption of the SINDy approach is
the sparsity of the parameter vector $\bm{\mathbf{\xi}}$. Sparsity is often
present in mathematical models of physical systems, which are frequently
amenable to reduced order and low-rank modeling. Specifically, sparsity refers
to the fact that many coefficients of $\bm{\mathbf{\xi}}$ are zero. Hence, a
dedicated (linear) regression approach needs to incorporate some
regularization technique. A regularized least-squares approach to solve (9)
minimizes the objective function

$\bm{\mathbf{\xi}}^{*}=\mathrm{argmin}_{\bm{\mathbf{\xi}}}\|\mathbf{z}-\mathbf{D}\bm{\mathbf{\xi}}\|_{2}^{2}+\lambda^{2}\|\bm{\mathbf{\xi}}\|^{2}.$
(10)

Different choices are available concerning the norm of the regularization term
$\|\bm{\mathbf{\xi}}\|$; popular choices include the $l^{2}$-norm (ridge
regression), the $l^{1}$-norm (LASSO regression) and the $l^{0}$-norm
(counting the number of non-zero elements). The choice of a Lasso least-square
approach has already been suggested in the original SINDy paper [6]. However,
eventually, a standard least-square algorithm together with iterative
thresholding was recommended. The SINDy procedure can be compactly summarized
as in Algorithm 1 [16].

Result: Approximation of $\bm{\mathbf{\xi}}^{*}$

$\bm{\mathbf{\xi}}^{(0)}=\mathbf{D}^{\dagger}\mathbf{z}$ (pseudoinverse
$\mathbf{D}^{\dagger}$);

$S=\\{1,\ldots,p\\}$ ;

for  _$k=1,\ldots,k_{\text{max}}$_ do

if  _$|\xi_{j}^{(k-1)}|\geq\lambda,j=1,\ldots,p,$_ then

Set $S=S\setminus\\{j\\}$;

end if

$\bm{\mathbf{\xi}}^{(k)}=\mathrm{argmin}_{\bm{\mathbf{\xi}}_{S}}\|\mathbf{D}\bm{\mathbf{\xi}}_{S}-\mathbf{z}\|\quad(\bm{\mathbf{\xi}}_{S,i}=0,\forall
i\notin S)$ ;

end for

Algorithm 1 SINDy

Hence, in the first step, an ordinary least-squares problem is solved. It can
be shown [16] that this iterative thresholding converges to a local minimum of
(10). Hence, the minimization can serve as a common starting point for a
Bayesian SINDy formulation.

Before moving to the Bayesian setting, we briefly discuss how to handle the
general case, with multiple state variables. To this end, we introduce the
matrices

$\displaystyle\mathbf{X}=\left(\begin{array}[]{c}\mathbf{x}^{\top}(t_{1})\\\
\mathbf{x}^{\top}(t_{2})\\\ \vdots\\\
\mathbf{x}^{\top}(t_{m})\end{array}\right),\qquad\dot{\mathbf{X}}=\left(\begin{array}[]{c}\dot{\mathbf{x}}^{\top}(t_{1})\\\
\dot{\mathbf{x}}^{\top}(t_{2})\\\ \vdots\\\
\dot{\mathbf{x}}^{\top}(t_{m})\end{array}\right),\qquad\mathbf{X},\dot{\mathbf{X}}\in\mathbb{R}^{m\times
n}.$

The library for approximating the unknown function $\mathbf{f}$ is then
constructed as

$\Theta(\mathbf{X})=\left[\mathbf{1}\quad\mathbf{X}\quad\mathbf{X}^{2}\cdots\right],\quad\Theta(\mathbf{X})\in\mathbb{R}^{m\times
p},$ (11)

where, for instance, row number $i$ of $\mathbf{X}^{2}$ is given by the vector
representation of $\mathbf{x}(t_{i})\otimes\mathbf{x}(t_{i})$. With this
notation at hand, we seek for an approximation

$\dot{\mathbf{X}}\approx\Theta(\mathbf{X})\bm{\mathbf{\xi}},$ (12)

with coefficients $\bm{\mathbf{\xi}}\in\mathbb{R}^{p\times n}$, to be
determined.

##  3 Bayesian SINDy with neuronized priors

To incorporate measurement noise and uncertainty, the linear model (9) is
modified as

$\mathbf{z}=\mathbf{D}\bm{\mathbf{\xi}}+\bm{\mathbf{\eta}},$ (13)

where $\mathbf{D}$ represents the design matrix. We assume, for the time
being, perfect knowledge of the state and hence of $\mathbf{D}$ and consider
uncertainty only in the data derivative. A Bayesian approach first formulates
a prior density $p(\bm{\mathbf{\xi}})$ for the unknown model parameter vector.
The likelihood expresses the probability of obtaining the data $\mathbf{z}$ as

$p(\mathbf{z}|\bm{\mathbf{\xi}})\sim\mathcal{N}(\mathbf{D}\bm{\mathbf{\xi}},\bm{\mathbf{\Sigma}}_{\bm{\mathbf{\epsilon}}}),$
(14)

where $\bm{\mathbf{\Sigma}}_{\bm{\mathbf{\epsilon}}}$ denotes the covariance
matrix of the Gaussian noise $\bm{\mathbf{\eta}}$. Bayes’ theorem then
expresses the posterior density as

$p(\bm{\mathbf{\xi}}|\mathbf{z})\propto
p(\mathbf{z}|\bm{\mathbf{\xi}})p(\bm{\mathbf{\xi}}).$ (15)

From (15), a MAP estimate is obtained by solving

$\bm{\mathbf{\xi}}^{*}=\mathrm{argmin}_{\bm{\mathbf{\xi}}}\left(-\log\left(p(\bm{\mathbf{\xi}}|\mathbf{z})\right)\right)=\mathrm{argmin}_{\bm{\mathbf{\xi}}}\left(\|\mathbf{z}-\mathbf{D}\bm{\mathbf{\xi}}\|_{\bm{\mathbf{\Sigma}}_{\bm{\mathbf{\epsilon}}}^{-1}}^{2}-\log\left(p(\bm{\mathbf{\xi}})\right)\right),$
(16)

with the $\bm{\mathbf{\Sigma}}_{\bm{\mathbf{\epsilon}}}^{-1}$-weighted
$l^{2}$-norm $\|\cdot\|_{\bm{\mathbf{\Sigma}}_{\bm{\mathbf{\epsilon}}}^{-1}}$.
A comparable formulation to (10) is obtained by choosing a Laplace prior
$p(\bm{\mathbf{\xi}})\propto\exp(-\lambda\|\bm{\mathbf{\xi}}\|_{1})$, whereas
a normal prior yields a ridge-regression. Many more sparsity-inducing priors
can be found in the literature. Examples include the horseshoe prior and the
spike and slab priors. The latter results in a MAP estimates with tight
connections to $l^{0}$-regularization, which completes the analogy between
Bayesian and standard sparse regression [17].

A more general approach employs the probabilistic state space model
formulation

$\displaystyle x_{i+1}$ $\displaystyle=x_{i}+\Delta
t\Theta(x_{i})\bm{\mathbf{\xi}}+\eta_{i},$
$\displaystyle\eta_{i}\sim\mathcal{N}(0,\sigma_{\eta}^{2}),$ (17)
$\displaystyle y_{j}$ $\displaystyle=x_{j}+\varepsilon_{j},$
$\displaystyle\varepsilon_{j}\sim\mathcal{N}(0,\sigma_{\varepsilon}^{2}),$
(18)

where $i=1,\ldots,m-1$ and $j=1,\ldots,k$. The measurements can be either
acquired at every discrete time step ($k=m$), or at subsets. The assumption of
independent, zero mean Gaussian observation noise is quite common and often
justifiable. The process noise $\eta_{i}$ can represent the discretization
error associated with the numerical derivative in our case. In this case, the
Gaussian assumption is more debatable since the numerical error is typically
structured. We maintain the independent normality assumption nevertheless and
refer to [18] for treatment in this direction. Recent studies on Bayesian
generalizations of SINDy consider $\eta_{i}$, but ignore $\varepsilon_{j}$. In
fact, the standard deterministic SINDy formulation is obtained as a MAP
estimate, when setting $\sigma_{\eta}=1$ and $\sigma_{\varepsilon}=0$, see
Theorem 4.2 of [13].

Following [19], we proceed by recasting (17) and (18) as

$\displaystyle\bm{\mathbf{\xi}}$ $\displaystyle\sim p(\bm{\mathbf{\xi}}),$
$\displaystyle\text{(prior)},$ (19) $\displaystyle x_{i}$ $\displaystyle\sim
p(x_{i}|x_{i-1},\bm{\mathbf{\xi}}),$ $\displaystyle\text{(dynamics)},$ (20)
$\displaystyle y_{j}$ $\displaystyle\sim p(y_{j}|x_{j}),$
$\displaystyle\text{(observation)},$ (21)

where we assume perfect knowledge of the initial conditions, for simplicity.
We now collect the measurements in a vector
$\mathbf{y}=(y_{1},\ldots,y_{k})^{\top}$. Accounting for observation noise now
requires to jointly estimate $(\mathbf{x},\bm{\mathbf{\xi}})$. Here, we are
mainly interested in learning the dynamics, encoded in $\bm{\mathbf{\xi}}$,
which is guided by the marginal posterior density

$p(\bm{\mathbf{\xi}}|\mathbf{y})=\int
p(\mathbf{x},\bm{\mathbf{\xi}}|\mathbf{y})\ \mathrm{d}\mathbf{x}.$

Computing the marginal posterior, which is mainly based on a formula for the
marginal likelihood $p(\mathbf{y}|\bm{\mathbf{\xi}})$, is commonly carried out
recursively, including a prediction and an update step, exploiting the
Markovian structure of the dynamical model. The main motivations being
computational simplicity and the sequential character of new measurements.
This has been covered in [13] and is not presented in any detail here.

Our goal is to obtain a quantification of uncertainty in the inferred model
parameters and models, where the latter is tightly related to the inclusion
probabilities of the individual library functions. Hence, we have to go beyond
MAP estimates and consider the posterior distribution. There exists a large
body of literature addressing the challenges in sparse Bayesian regression.
Here, we recall two main classes, continuous global-local and spike and slab
priors. A continuous global-local prior is given as

$p(\xi_{j}|\lambda_{j},\tau)\sim\mathcal{N}(0,\tau\lambda_{j}^{2}),$ (22)

where $\tau,\lambda_{j}$ are called the global and local shrinkage parameter
or variance components, respectively. This type of model is discussed in
detail, for instance in [20]. For simplicity, we consider a zero mean in the
discussion of different priors. The global shrinkage parameter $\tau$ pushes
the coefficients towards zero, whereas some parameters are allowed to assume
large values if needed, by choosing a heavy-tailed distribution for
$\lambda_{j}$. The model (22) is hierarchical and different choices for the
distributions of $\tau,\lambda$ give rise to different priors. Popular choices
are the Horseshoe prior [21] or the Bayesian Lasso [22].

On the other hand, spike and slab priors employ a mixture distribution

$p(\xi_{j}|c,\delta,\gamma_{j})\sim\gamma_{j}\mathcal{N}(0,\tau
c^{2})+(1-\gamma_{j})\mathcal{N}(0,\tau\delta^{2}),$ (23)

where the Bernoulli variable $\gamma_{j}$ expresses the inclusion probability
of the coefficient $\xi_{j}$. By setting $\delta\ll c$ the second distribution
on the right-hand-side of (23) is concentrated around zero (spike), whereas
$\mathcal{N}(0,\tau c^{2})$ is less concentrated (slab). Spike and slab priors
often provide excellent sparsity in regression [23]. Whereas (23) describes
the continuous version, a discrete counterpart is obtained by replacing
$\mathcal{N}(0,\tau\delta^{2})$ with a delta function at the origin as

$p(\xi_{j}|c,\delta,\gamma_{j})\sim\gamma_{j}\mathcal{N}(0,\tau
c^{2})+(1-\gamma_{j})\delta_{0}.$ (24)

In this study, we combine multiple state-of-the-art techniques for dynamical
system identification with uncertainty quantification. To this end, we employ
a recently proposed formulation of priors, the so-called neuronized priors
[14], which provide a unified treatment. Moreover, we utilize an efficient
MCMC algorithm [14] to infer posterior distributions of the model parameters.
A neuronized prior is formulated as,

$\displaystyle p(\xi_{j}|\alpha_{0})\sim T(\alpha_{j}-\alpha_{0})w_{j},$ (25)

where $\alpha_{0}$ is a fixed hyperparameter and $T$ is an activation function
such as the ReLu function. The notion of an activation function is borrowed
from artificial neural networks, whence the name neuronized prior. The model
is hierarchical, i.e.,

$\displaystyle p(\alpha_{j})$ $\displaystyle\sim N(0,1),$ (26) $\displaystyle
p(w_{j}|\tau_{w})$ $\displaystyle\sim N(0,\tau_{w}^{2}).$ (27)

The main interest of the formulation is its generality, i.e., by choosing
different activation functions, a spike and slab, Horseshoe, and Lasso prior
are recovered, at least asymptotically. Additionally, the formulation allows
obtaining efficient algorithms, avoiding the necessity to use latent binary
indicator variables, as pointed out in [24]. Different versions of neuronized
priors, together with their descriptions are presented in Figure 1.

Lasso$\xi$

Horseshoe$\xi$

ReLu$\xi$

Figure 1: Examples of neuronized priors. All priors enforce sparsity by
allocating probability mass near (or at) zero. Left: Lasso prior for
$\tau_{w}=0.1$, which is obtained by setting $T(\xi)=\xi$, $\alpha_{0}=0$.
Middle: Approximate horseshoe prior with
$T(\xi)=\exp(0.37\mathrm{sign}(\xi)\xi^{2}+0.89\xi+0.08)$, $\alpha_{0}=0$ and
$\tau_{w}=0.05$. The distribution is only asymptotically equivalent to the
horseshoe prior. Right: Discrete spike and slab prior with
$T(\xi)=\max(0,\xi)$, which corresponds to the ReLu activation function. In
this case $\alpha_{0}=0.5$ plays the role of the sparsity parameter and the
slab variance is chosen as $\tau_{w}=0.1$. The vertical axes-range differs
among the three figures to allow for a better visualization.

Having formulated the prior, two different strategies can be found in the
literature for model selection. On the one hand, variable selection methods
aim to identify a single model without necessarily quantifying uncertainty. We
mention for instance the expectation-maximization algorithm presented in [25].
This path is typically chosen, if the sampling from the full posterior
distribution is too expensive. Since the problems we consider here, in
particular the size of the data and the library, are not too large, we employ
an MCMC approach instead. The algorithm to approximate the marginal posterior
distribution
$p(\bm{\mathbf{\alpha}},\bm{\mathbf{w}}|\mathbf{y},\tau_{w},\alpha_{0})$ is
summarized in Algorithm 2, which is adopted from [14]. This algorithm
successively updates $\mathbf{w}$ and $\mathbf{\alpha}$, where $\mathbf{w}$
can be sampled from a normal distribution with covariance matrix

$\tilde{\bm{\mathbf{\Sigma}}}(\bm{\mathbf{\alpha}})=\sigma_{\eta}^{2}(D_{\bm{\mathbf{\alpha}}}\mathbf{D}^{\top}\mathbf{D}D_{\bm{\mathbf{\alpha}}}+\sigma_{\eta}^{2}\tau_{w}^{-2}\mathbf{I}_{p\times
p})^{-1},\
D_{\bm{\mathbf{\alpha}}}=\mathrm{diag}(\alpha_{1},\ldots,\alpha_{p})$

and mean value
$\tilde{\bm{\mathbf{\mu}}}(\mathbf{y},\bm{\mathbf{\alpha}})=\tilde{\bm{\mathbf{\Sigma}}}(\bm{\mathbf{\alpha}})D_{\bm{\mathbf{\alpha}}}\mathbf{D}^{\top}\mathbf{y}$.
The sampling distribution for $\mathbf{\alpha}$ depends on the chosen
activation function and is therefore, not explicitly stated.

Based on the MCMC approximation, we obtain the posterior of the state as

$p(x(t)|\mathbf{y},\tau_{w},\alpha_{0})=\int
p(x(t)|\mathbf{w},\bm{\mathbf{\alpha}},\tau_{w},\alpha_{0})p(\mathbf{w},\bm{\mathbf{\alpha}}|\mathbf{y},\tau_{w},\alpha_{0})\
\mathrm{d}\mathbf{\bm{\mathbf{\alpha}}}\mathrm{d}\mathbf{\mathbf{w}}.$ (28)

We will comment on the choice of the remaining hyperparameters
$\tau_{w},\alpha_{0}$ in the numerical results section. As previously
mentioned, estimating the marginal Likelihood requires a Kalman filter
together with marginalization over the state. In the simple case where
$\sigma_{\varepsilon}=0$, this is not required and we obtain the simplified
expression

$p(\mathbf{y}|\bm{\mathbf{\alpha}},\bm{\mathbf{w}})\propto\mathrm{exp}\left(-\frac{1}{2}\frac{\|\mathbf{D}\bm{\mathbf{\xi}}(\bm{\mathbf{\alpha}},\bm{\mathbf{w}})-\mathbf{z}\|_{2}^{2}}{\sigma_{\eta}^{2}}\right),$

which will be used in the following examples.

Input: Hyper parameter $\tau_{w},\alpha_{0}$ (omitted in formulas below)

Input: Estimate of data marginal Likelihood
$\hat{p}(\mathbf{y}|\bm{\mathbf{\alpha}},\bm{\mathbf{w}})\approx
p(\mathbf{y}|\bm{\mathbf{\alpha}},\bm{\mathbf{w}})$

Result: Sample $(\mathbf{w}^{(i)},\bm{\mathbf{\alpha}}^{(i)})$ of posterior
density $p(\bm{\mathbf{\alpha}},\bm{\mathbf{w}}|\mathbf{y})$

for  _$i=1,\ldots,N$_ do

Sample $\mathbf{w}^{(i)}$ from
$\mathbf{w}|\mathbf{y},\bm{\mathbf{\alpha}}\sim\mathcal{N}(\tilde{\bm{\mu}}(\mathbf{y},\bm{\mathbf{\alpha}}),\tilde{\bm{\Sigma}}(\bm{\mathbf{\alpha}}))$;

Sample $\bm{\mathbf{\alpha}}^{(i)}$ based on
$\hat{p}(\bm{\mathbf{\alpha}}|\mathbf{y},\bm{\mathbf{w}}^{(i)})$ with a random
walk Metropolis Hastings algorithm;

end for

Algorithm 2 MCMC with neuronised prior

##  4 Results

The Bayesian SINDy algorithm with neuronised MCMC is applied on three systems;
two generic examples, the linear pendulum and the Lorenz system, and one
aerodynamic application, for which experimental measurement data are used.

###  4.1 Linear Pendulum

$x_{1}$$t$

$x_{1}$$t$

Figure 2: Results for the identification of the linear pendulum with a Lasso
prior. Reference trajectory, posterior mean and credible intervals covering
$90\%$ of the posterior samples for noise level $\sigma_{\eta}=10^{-1}$
(left), and noise level $\sigma_{\eta}=10^{-2}$ (right).

$\dot{x}_{1}$$t$

$\dot{x}_{1}$$t$

Figure 3: Results for the identification of the linear pendulum with a Lasso
prior. Time derivative data $\mathbf{z}$ and posterior of
$\Theta(\mathbf{X})\bm{\mathbf{\xi}}$. The $90\%$ credible intervals decrease
together with the noise level $\sigma_{\eta}=10^{-1}$ (left) and
$\sigma_{\eta}=10^{-2}$ (right).

The simple pendulum example, adapted from [13], enables easy testing of the
method. The dynamical system is given by

$\frac{d}{dt}\left(\begin{array}[]{c}x_{1}\\\
x_{2}\end{array}\right)=\left(\begin{array}[]{cc}0&1\\\
-g/L&0\end{array}\right)\left(\begin{array}[]{c}x_{1}\\\
x_{2}\end{array}\right),$ (29)

where $g=9.81\,$ms-1 and $L=1\,$m. The system is simulated in time with
MATLAB’s built-in Runge Kutta ODE45 time integrator over a time interval
$[0,4]\,$s. Using this example, we investigate different activation functions
(hence, different priors) and the influence of the data noise. Data are
generated at $m=50$ uniformly spaced time intervals. In particular, the data
vector is created as

$\mathbf{z}_{j}=\dot{\mathbf{x}}(t_{j})+\sigma_{\eta}\bm{\mathbf{\eta}}_{j},\quad
j=1,\ldots,m$

where $t_{j}$ denotes the discrete time, $\sigma_{\eta}$ denotes the data
noise amplitude, and
$\bm{\mathbf{\eta}}_{j}=(\eta_{j,1},\eta_{j,2})^{\top},\eta_{j,1},\eta_{j,2}\sim\mathcal{N}(0,1)$
are independently drawn. The library in this case is chosen as

$\Theta(\mathbf{X})=\left[\mathbf{X}\ \mathbf{X}^{2}\right].$

The two equations defining the system are identified individually with the
neuronized MCMC algorithm. The hyperparameters are chosen, as described in
Section 3. For each MCMC run a sample size of $M=10^{5}$ elements is chosen,
which yields a Geweke score

$Z=\frac{\mu_{A}-\mu_{B}}{\sqrt{\frac{\sigma_{A}}{n_{A}}+\frac{\sigma_{B}}{n_{B}}}}$
(30)

smaller than $2$. In equation (30), $\mu_{A/B},\sigma_{A/B}$ refer to the
moments of Markov sub-chains of size $n_{A}=0.1M,n_{B}=0.5M$ extracted from
the beginning and the end, respectively. We always neglect the first 20$\%$ of
the Markov Chain to remove the burn-in phase.

In Figure 2, we present the first component of the reference solution along
with the identified system with the Lasso prior-based posterior distribution.
The latter is obtained by constructing a family of dynamical models according
to (5), and solving them with the same Runge Kutta method. The posterior mean
represents the arithmetic mean over all posterior trajectories and the $90\%$
credible interval bounds cover $90\%$ of all trajectories. We observe the
larger uncertainty band, reflected by the posterior credible interval on the
left ($\sigma_{\eta}=10^{-1}$), compared to the right
($\sigma_{\eta}=10^{-2}$). Since for the Lasso prior $\alpha_{0}=0$, only a
single hyperparameter needs to be selected, which is chosen according to the
signal-to-noise ratio as [14],

$\tau_{w}=\|\bm{\mathbf{\xi}}\|_{2}/(\sqrt{p\mathbb{E}[T^{2}(\alpha)]}),$ (31)

where we draw $\alpha\sim\mathcal{N}(0,1)$ and approximate the expected value
with the Monte Carlo method. This yields $\tau_{w}=0.0039$ and
$\tau_{w}=0.0148$ for the first and second equation, respectively. The
posterior distributions for the Horshoe and discrete spike and slab prior are
qualitatively similar. They are reported in Appendix A. In Figure 3, we report
the data $\mathbf{z}$ and the posterior credible intervals for the linear
model $\Theta(\mathbf{X})\bm{\mathbf{\xi}}$, for each equation separately. We
again observe a decreasing uncertainty, when the noise level is reduced from
$\sigma_{\eta}=10^{-1}$ on the left, to $\sigma_{\eta}=10^{-2}$ on the right.
The results for the other available prior choices are largely similar and are
not reported here.

Our investigations are not complete without examining the accuracy of the
different prior methods with regard to estimating the regression coefficients.
In Figure 4, the true parameter values are shown together with the posterior
distributions of the estimated coefficients. The posterior distribution is
visualized with box plots. The Lasso and ReLu prior deliver visually the best
estimates, whereas the scatter in the Horseshoe prior is larger. However,
comapred to the Lasso, the ReLu prior improves the sparsity in the solution.
Whenever there is a real interest in selecting one “winning” model, a common
choice is selection of the median model, i.e., a basis function is included if
$P(\gamma_{i}=1)>0.5$. In this regard, both Lasso and Horseshoe priors
typically yield the full set of of variables. In other words, most
coefficients are small but non-zero with probability greater than $0.5$. Here,
an additional threshold would need to be applied, in order to recover a sparse
model. Only the ReLu prior delivers a sparse model directly based on the
median model selection criterion, where only the coefficient associated to
$x_{1}x_{2}$ is erroneously included if the noise is too large. In view of
these findings, we consider ReLu activation functions for the remaining, more
challenging, examples to follow.

Lasso$\xi_{i}$$i$

Horseshoe$\xi_{i}$$i$

ReLu$\xi_{i}$$i$

Figure 4: True parameter values for the first equation of the pendulum and box
plots for posterior sample of the estimated coefficients. The library basis
functions, indexed with $i$, are $x_{1}$ ($i=1$), $x_{2}(i=2)$,
$x_{1}x_{2}(i=3)$, $x_{1}^{2}(i=4)$, $x_{2}^{2}(i=5)$. The results are shown
for a noise level $\sigma_{\eta}=0.1$ and $m=50$ data points are used.

###  4.2 Lorenz System

The second application employing the neuronised MCMC algorithm is the Lorenz
system. Here too, the reference system dynamics are known. The Lorenz system
is particularly chosen because of its chaotic nature. It is driven by a set of
three coupled nonlinear ordinary differential equations (ODEs) given by

$\begin{split}\dot{x}_{1}&=c_{1}(x_{2}-x_{1}),\\\
\dot{x}_{2}&=x_{1}(c_{2}-x_{3})-x_{2},\\\
\dot{x}_{3}&=x_{1}x_{2}-c_{3}x_{3}.\end{split}$ (32)

The selected parameters are $c_{1}=10$, $c_{2}=28$ and $c_{3}=2.667$ with
initial conditions $(x_{1}^{0},x_{2}^{0},x_{3}^{0})=(-8,8,27)$. The simulation
is performed with a time step $\Delta t=0.1$ for a total time of $T=100$. The
initial state $(x_{1}^{0},x_{2}^{0},x_{3}^{0})^{\top}$ is propagated in time
using MATLAB’s ODE45 function, and the resulting time series
$[\mathbf{x_{1}}\mathbf{x_{2}}\mathbf{x_{3}}]$ is used to compute the
derivatives $\mathbf{\dot{x}_{1}},\mathbf{\dot{x}_{2}}$ and
$\mathbf{\dot{x}_{3}}$. To simulate numerical differentiation noise, normally
distributed random noise $\bm{\eta}\sim\mathcal{N}(\bm{0}_{m\times
1},\mathbf{I}_{m\times m})$ is separately added to
$\mathbf{\dot{x}_{1}},\mathbf{\dot{x}_{2}}$ and $\mathbf{\dot{x}_{3}}$.

Figure 5: The reference trajectory, posterior mean and credible intervals
covering 90% of the posterior samples of the autocorrelation distributions
$\mathbf{x}_{1},\mathbf{x}_{2}$ and $\mathbf{x}_{3}$.

Figure 6: Each row in the above figure corresponds to one particular equation
of the Lorenz system. Each row contains the histograms of the coefficients of
each of the basis functions comprising $\mathbf{D}$. The values of the
coefficients used to generate the data vectors are indicated with a red star.

With the data vectors computed, the neuronised MCMC algorithm is applied once
on each of the data vectors. For all the three vectors, a discrete spike and
slab/ReLu prior with $\alpha_{0}=0.5$ is used. As per equation (31),
$\tau_{w}$ the parameter that determines the variance of the prior’s slab
part, is set to 12.462, 24.791 and 2.518, for the three components
$\mathbf{\dot{x}},\mathbf{\dot{y}}$, and $\mathbf{\dot{z}}$, respectively,
whereas $\alpha_{0}$ is set to 0.5 to obtain a high degree of sparsity. While
computing $\tau_{w}$’s, the least squares estimate
$\bm{\xi}^{i}=(\mathbf{D}^{\top}\mathbf{D})^{-1}\mathbf{D}^{\top}\mathbf{z}^{i}$
is used, with $i=1\dots 3$ for the 3 Lorenz equations. Note that the
$\tau_{w}$ values are set before starting the MCMC for each of the equations.
In this example and the following one, the error variance $\sigma_{\eta}^{2}$
was also updated within the MCMC, by imposing an inverse gamma prior Inv-
Gam($a_{0},b_{0}$) on $\sigma_{\eta}^{2}$, with both $a_{0}$ and $b_{0}$
assigned a value of 1. On the imposition of such a prior, the posterior of
$\sigma_{\eta}^{2}$ too has an inverse gamma form (see [14]). In order to
facilitate efficient sampling of $\bm{\alpha}$ for the case of the ReLu prior,
$\bm{\alpha}$ is directly sampled from a mixture of two truncated Gaussians
[14], instead of the Random Walk Metropolis Hastings step in Algorithm 2. This
sampling strategy is performed for both the current and the next application.
The following basis library:
$\mathbf{D}(\mathbf{x_{1}},\mathbf{x_{2}},\mathbf{x_{3}})=\left[\mathbf{x_{1}},\mathbf{x_{2}},\mathbf{x_{3}},\mathbf{x_{1}x_{2}},\mathbf{x_{1}x_{3}},\mathbf{x_{2}x_{3}}\right]$
is employed. The MCMC is run for a total of $M=10000$ iterations and the first
3000 sample points are discarded as burn-in sample points. This yields a
Geweke score smaller than 2.

The posterior distributions of the coefficients $\bm{\xi}$ are presented in
Figure 6. Each row corresponds to one equation of the Lorenz system. The
reference coefficients are marked with red dots. As the figure shows, the
histograms of coefficients corresponding to insignificant bases are associated
with spikes at zero, whereas those of the non-zero coefficients are correctly
centered around their corresponding reference values. With the obtained
posteriors, using a median model selection criterion, the correct basis
functions for each of the equations in equation (32) are recovered. In Table
1, we compare the single point estimates obtained using the neuronised prior
algorithm and the SINDy algorithm against the reference values. For the
neuronised prior estimates, the posterior mean of the coefficients included in
the median model is used. For SINDy, the estimates obtained for $\lambda=0.1$
are shown. It is to be noted that the SINDy estimates for other values of
$\lambda$, that is $\lambda=0.3,0.5,0.7$ and 0.9, were also computed and they
were found to be identical to the estimate for $\lambda=0.1$, and hence they
are not reported here. It can be seen that both the SINDy and the neuronised
prior estimates match very closely with the reference values.

$\dot{x}_{1}=10(x_{2}-x_{1})$ | $\dot{x}_{2}=x_{1}(28-x_{3})-x_{2}$ | $\dot{x}_{3}=x_{1}x_{2}-8/3x_{3}$
---|---|---
Bases | Ref. | Neu. Pr. | SINDy | Ref. | Neu. Pr. | SINDy | Ref. | Neu. Pr. | SINDy
$\lambda=0.1$ | $\lambda=0.1$ | $\lambda=0.1$
$x_{1}$ | -10 | -10.0076 | -10.0076 | 28 | 28.0295 | 28.0296 | 0 | 0 | 0
$x_{2}$ | 10 | 10.0101 | 10.0102 | -1 | -1.0097 | -1.0097 | 0 | 0 | 0
$x_{3}$ | 0 | 0 | 0 | 0 | 0 | 0 | -2.6667 | -2.6679 | -2.6679
$x_{1}x_{2}$ | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0.9999 | 0.9999
$x_{1}x_{3}$ | 0 | 0 | 0 | -1 | -1.0007 | -1.0007 | 0 | 0 | 0
$x_{2}x_{3}$ | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0

Table 1: Single point estimates comparison between SINDy and neuronized prior
system dynamics identification. For the latter the posterior means of the
median model are shown.

With the posteriors of $\bm{\mathbf{\xi}}^{i}$ determined, the states
$\mathbf{x}_{1},\mathbf{x}_{2}$ and $\mathbf{x}_{3}$ are propagated through
equation (32). Since a direct comparison of time series is meaningless for
chaotic systems, a quantitative assessment of the model requires an
alternative metric. In this study, we compare the autocorrelation function,
computed as
$R_{\mathbf{x}\mathbf{x}}(t_{1},t_{2})=\mathbb{E}[\mathbf{x}_{t_{1}}\mathbf{x}_{t_{2}}]$,
of the reference data and that of the model. Figure 5 presents window-averaged
auto-correlation distributions for the states $\mathbf{x}_{1}$,
$\mathbf{x}_{2}$ and $\mathbf{x}_{3}$. Forty time-windows of length 2.5 units
are used for the averaging. It can be observed that the reference curves match
the posterior mean distributions very well.

###  4.3 Aerodynamic Application

The final application employs measured data from an aerodynamic experiment.
The experiment consists of an airfoil equipped with a pitching flap. Here, we
are interested in identifying a model for the unsteady lift coefficient
$C_{l}$ as a function of the pitch angle $\delta$ and enriched features
thereof. Both $C_{l}$ and $\delta$ are experimentally acquired in real-time
with surface pressure sensors and position encoders, respectively. The
particular test case is purposely selected to exhibit challenging dynamics
through dynamic stall. Details on the experimental setup as well as the
aerodynamic problem are found in [26] and [27].

The dynamics is learned from the following basis library

$\mathbf{D}(\mathbf{C}_{l},\bm{\delta})=[\mathbf{C}_{l},\dot{\bm{\delta}},\bm{\delta},\bm{\delta}\dot{\bm{\delta}},\dot{\bm{\delta}}^{2},\bm{\delta}^{2},\dot{\bm{\delta}}^{3},\bm{\delta}^{3},\dot{\bm{\delta}}_{\frac{1}{4}},\bm{\delta}_{\frac{1}{4}},\dot{\bm{\delta}}_{\frac{1}{2}},\bm{\delta}_{\frac{1}{2}},\dot{\bm{\delta}}_{\frac{1}{8}},\bm{\delta}_{\frac{1}{8}}],$

consisting of 14 basis functions, where $\bm{\delta}\dot{\bm{\delta}}$ refers
to an element-wise multiplication of both vectors. Here, the subscripts
$\frac{1}{8}$, $\frac{1}{4}$ and $\frac{1}{2}$ refer to corresponding time-
delayed basis functions by fractions of a time period. The derivatives
$\dot{\mathbf{C}}_{l}$ and $\dot{\bm{\delta}}$ are obtained by first smoothing
$\mathbf{C}_{l}$ and $\bm{\delta}$ and then using MATLAB’s gradient function.
A noise vector $\bm{\eta}\sim\mathcal{N}(\bm{0}_{m\times 1},0.01^{2}\
\mathbf{I}_{m\times m})$ is added to $\dot{\mathbf{C}}_{l}$ to obtain
$\mathbf{z}$ for the subsequent steps. The model coefficients $\bm{\xi}$ are
estimated with the neuronised MCMC algorithm employing a discrete spike and
slab/ReLu prior with $\alpha_{0}=1.5$ and $\tau_{w}=227.468$. The hyper-
parameter $\alpha_{0}$ is assigned this value to obtain a high degree of
sparsity, and the parameter $\tau_{w}$ is set as per equation (31). As
described in the Lorenz example, $\tau_{w}$ is set before starting the MCMC
and the least squares estimate of $\bm{\xi}$ is used in its computation. Since
we are mainly interested in the unsteady lift dynamics, both $\mathbf{D}$ and
the data vector are first de-meaned. To simplify the sparse optimisation
procedure, we additionally normalize $\mathbf{D}$ columnwise with the
respective standard deviations of each basis vector. Similar to the Lorenz
example, the MCMC is run for $M=10000$ iterations, which results in a Geweke
score of less than 2. The first 3000 iterations are discarded as burn-in
sample points. The magnitudes of $a_{0}$ and $b_{0}$ pertaining to the inverse
gamma prior for $\sigma^{2}$ are set to 1, whereas the initial value of
$\sigma$ is set to 0.1. The value of $\sigma$ gets updated as part of the MCMC
procedure.

The resulting posterior distributions of the coefficients $\bm{\xi}$ are shown
in Figure 7. As the figure shows, if we consider the median model, only five
basis functions contribute to the dynamics,
$\dot{\bm{\delta}},\bm{\delta}\dot{\bm{\delta}},\bm{\delta}^{2},\dot{\bm{\delta}}^{3}$,
and $\bm{\delta}^{3}$. This yields a model with the following structure,

$\dot{\mathbf{C}_{l}}\approx\xi_{1}\dot{\bm{\delta}}+\xi_{2}\bm{\delta}\dot{\bm{\delta}}+\xi_{3}\bm{\delta}^{2}+\xi_{4}\dot{\bm{\delta}}^{3}+\xi_{5}\bm{\delta}^{3}$
(33)

where the coefficients can assume values according to their respective
posteriors. The identified model offers an intriguing structure. While the
pitch angle rate $\dot{\bm{\delta}}$ is an expected basis in unsteady
aerodynamics modeling, the remaining ones are not. Their role is to model
nonlinear interactions deviating from harmonic input-output response. This
nonlinearity is caused by the so-called dynamic stall, where the rapid
pitching causes a strong vortex to be shed from the leading edge that travels
downstream over the wing. Hence, these basis functions are compensating for
the assumed restrictive model structure. Airfoils with dynamic stall require
more complex dynamical models, such as a Duhamel integral for the indical
response approach [28], or an additional state variable for the ONERA modeling
approach [29].

Figure 7: Posterior distributions of the model coefficients $\bm{\xi}$ for the
aerodynamic example. Starting with a library of 14 basis functions, the
algorithm identifies five bases
$\dot{\bm{\delta}},\bm{\delta}\dot{\bm{\delta}},\bm{\delta}^{2},\dot{\bm{\delta}}^{3}$
and $\bm{\delta}^{3}$ that contribute to the lift dynamics.

Figure 8: (Left) the lift coefficient $C_{l}$ and (right) its time derivative
$\dot{C_{l}}$ over one pitching period $T$. The three distributions are the
reference trajectory, posterior mean, and credible intervals covering 90% of
the posterior samples.

As done in the previous examples, the entire posterior of $\bm{\xi}$ is used
to reconstruct the dynamics $\dot{\mathbf{C}_{l}}\approx\mathbf{D}\bm{\xi}$,
which is subsequently used to reconstruct $\mathbf{C}_{l}$. The reference lift
coefficient distribution along with the distribution of its time derivative
over a time period of one pitching period $T$ are shown in Figure 8. Despite
the restrictive model structure, both distributions compare favorably to the
reference data.

##  5 Summary and conclusions

In this paper, we have presented an uncertainty-aware method for sparse
dynamical system identification and applied it to both generic model examples
and an aerodynamic application with real data. Here, the main contribution to
Bayesian dynamical system identification was the adaption of a unified prior
family, the neuronized prior. In particular, by choosing different activation
functions we were able to recover and implement well-known sparsity priors,
such as spike and slab, Horseshoe and Lasso priors in particular, in a common
framework. We then employed dedicated MCMC methods to infer dynamical reduced
order models together with parameter and structural model uncertainties. The
framework also allows to obtain a single model with a model selection
strategy, which is desired from an application point of view to gain physical
insights from the data. Several directions are to be further explored. We have
limited the considerations to moderate basis dimensions and additional efforts
are needed to address large-scale problems. Moreover, including model
parameters into the framework would allow to answer relevant tasks in control
applications. Finally, more general noise structure should be considered to
reflect the structure of complex data.

## References

  * [1] Josh Bongard and Hod Lipson.  Automated reverse engineering of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 104(24):9943–9948, June 2007.
  * [2] F. C. Fu and J. B. Farison.  On the Volterra series functional evaluation of the response of non-linear discrete-time systems.  International Journal of Control, 18(3):553–558, July 1973.
  * [3] Chris Chatfield.  Time-series forecasting.  CRC press, 2000.
  * [4] Jer-Nan Juang.  Applied system identification.  Prentice-Hall, Inc., 1994.
  * [5] Tong Wang, Huijun Gao, and Jianbin Qiu.  A combined adaptive neural network and nonlinear model predictive control for multirate networked industrial process control.  IEEE Transactions on Neural Networks and Learning Systems, 27(2):416–425, 2016.
  * [6] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the national academy of sciences, 113(15):3932–3937, 2016.
  * [7] Robert K Niven, Ali Mohammad-Djafari, Laurent Cordier, Markus Abel, and Markus Quade.  Bayesian identification of dynamical systems.  Multidisciplinary Digital Publishing Institute Proceedings, 33(1):33, 2020.
  * [8] Sheng Zhang and Guang 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):20180305, 2018.
  * [9] Maria Maddalena Barbieri, James O Berger, et al.  Optimal predictive model selection.  The annals of statistics, 32(3):870–897, 2004.
  * [10] Chris Hans.  Model uncertainty and variable selection in Bayesian lasso regression.  Statistics and Computing, 20(2):221–229, 2010.
  * [11] R Fuentes, R Nayek, P Gardner, N Dervilis, T Rogers, K Worden, and EJ Cross.  Equation discovery for nonlinear dynamical systems: A Bayesian viewpoint.  Mechanical Systems and Signal Processing, 154:107528, 2021.
  * [12] Rimple Sandhu, Mohammad Khalil, Chris Pettit, Dominique Poirel, and Abhijit Sarkar.  Nonlinear sparse Bayesian learning for physics-based models.  Journal of Computational Physics, page 109728, 2020.
  * [13] Nicholas Galioto and Alex Gorodetsky.  Bayesian system id: Optimal management of parameter, model, and measurement uncertainty.  Nonlinear Dynamics, 102:241–267, 2020.
  * [14] Minsuk Shin and Jun S Liu.  Neuronized priors for Bayesian sparse linear regression.  Journal of the American Statistical Association, pages 1–43, 2021\.
  * [15] Hua Liang and Hulin Wu.  Parameter estimation for differential equation models using a framework of measurement error in regression models.  Journal of the American Statistical Association, 103(484):1570–1583, 2008.
  * [16] Linan Zhang and Hayden Schaeffer.  On the convergence of the sindy algorithm.  Multiscale Modeling & Simulation, 17(3):948–972, 2019.
  * [17] Nicholas G Polson and Lei Sun.  Bayesian l 0-regularized least squares.  Applied Stochastic Models in Business and Industry, 35(3):717–731, 2019.
  * [18] Patrick R Conrad, Mark Girolami, Simo Särkkä, Andrew Stuart, and Konstantinos Zygalakis.  Statistical analysis of differential equations: introducing probability measures on numerical solutions.  Statistics and Computing, 27(4):1065–1082, 2017.
  * [19] Simo Särkkä.  Bayesian filtering and smoothing.  Number 3. Cambridge University Press, 2013.
  * [20] Nicholas G Polson and James G Scott.  Shrink globally, act locally: Sparse Bayesian regularization and prediction.  Bayesian statistics, 9(501-538):105, 2010.
  * [21] Carlos M Carvalho, Nicholas G Polson, and James G Scott.  The Horseshoe estimator for sparse signals.  Biometrika, 97(2):465–480, 2010.
  * [22] Trevor Park and George Casella.  The Bayesian lasso.  Journal of the American Statistical Association, 103(482):681–686, 2008.
  * [23] Juho Piironen, Aki Vehtari, et al.  Sparsity information and regularization in the horseshoe and other shrinkage priors.  Electronic Journal of Statistics, 11(2):5018–5051, 2017.
  * [24] Lizhen Nie and Veronika Ročková.  Bayesian bootstrap spike-and-slab lasso.  arXiv preprint arXiv:2011.14279, 2020.
  * [25] Veronika Ročková and Edward I George.  Emvs: The em approach to Bayesian variable selection.  Journal of the American Statistical Association, 109(506):828–846, 2014.
  * [26] Johannes Pohl and Richard Semaan.  Quantification and modeling of dynamic lift on a DLR-F15 research airfoil with active trailing-edge flap.  In AIAA Paper 2020-0556. American Institute of Aeronautics and Astronautics, 2020.
  * [27] Johannes Pohl, Richard Semaan, and Anya R. Jones.  Dynamic lift measurements on an airfoil with periodic flap motion at high Reynolds number.  In AIAA Paper 2019-1396. American Institute of Aeronautics and Astronautics, 2019.
  * [28] Wagner Herbert.  Über die Entstehung des dynamischen Auftriebes von Tragflügeln.  Zeitschrift fur Angewandte Mathematic und Mechanik, 5(1):17–35, 1925.
  * [29] K. W. McAlister, O. Lambert, and D. Petot.  Application of the ONERA Model of Dynamic Stall.  Technical Report NASA-A-9824, NASA, Moffett field, CA, Ames Research center, 1984.

##  Appendix A Additional Results for the Pendulum

In Figure 9 we report the posterior results for the Horseshoe prior with
activation function $T(x)=\exp(0.5\mathrm{sign}(x)x^{2}+0.733x)$. Again
$\alpha_{0}=0$, however, a different value for $\tau_{w}$ needs to be chosen
to yield good results. In [14], it is recommended to choose $\tau_{w}$ such
that $1-\mathbb{E}[1/(1+T(\alpha)^{2}\tau_{w}^{2})]$,
$\alpha\sim\mathcal{N}(0,1)$ equals the expected portion of zero values, which
is $\approx 0.2$ in our case. In this regard, a value $\tau_{w}$ seems
appropriate and indeed yields satisfactory results for the case
$\sigma_{\eta}=0.1$. However, for $\sigma_{\eta}=0.2$ we observed that larger
values are needed. The results reported in Figure 9 (right) have been obtained
with $\tau_{w}=10,20$ for the first and second equation respectively.

$x_{1}$$t$

$x_{1}$$t$

Figure 9: Results for the identification of the linear pendulum with a
Horseshoe prior. Reference trajectory, posterior mean and credible intervals
covering $90\%$ of the posterior samples. These samples are obtained by
computing the time trajectory for posterior samples of $\bm{\mathbf{\xi}}$.
Left: data size $n=10$ and noise level $\sigma_{\eta}=10^{-1}$. Right: data
size $n=50$ and noise level $\sigma_{\eta}=10^{-2}$.

Finally, in Figure 10 the results for the ReLu activation function are
reported, which corresponds to a discrete spike and slab prior. Here, the
hyperparameter $\tau_{w}$ is again chosen according to the signal-to-noise
ratio, as described in Section 4.1, where also the concrete numerical values
are given. In case of the ReLu activation function, a second hyperparameter
$\alpha_{0}$ has to be chosen. In [14] this parameter was related to the
expected fraction of non-zero coefficients $\eta$ as
$\alpha_{0}=\Phi^{-1}(\eta)$, where $\Phi$ refers to the CDF of the standard
normal distribution. Note, that the parameter $\eta$ is commonly employed as
the Bernoulli-hyperparameter in the spike and slab formulation with a latent,
model-inclusion, variable.

$x_{1}$$t$

$x_{1}$$t$

Figure 10: Results for the identification of the linear pendulum with a
discrete spike and slab prior. Reference trajectory, posterior mean and
credible intervals covering $90\%$ of the posterior samples. These samples are
obtained by computing the time trajectory for posterior samples of
$\bm{\mathbf{\xi}}$. Left: data size $n=10$ and noise level
$\sigma_{\eta}=10^{-1}$. Right: data size $n=50$ and noise level
$\sigma_{\eta}=10^{-2}$.