# Encoding nonlinear and unsteady aerodynamics of limit cycle oscillations
using nonlinear sparse Bayesian learning

Rimple Sandhu
Department of Civil and Environmental Engineering
Carleton University
Ottawa, ON, Canada
&Brandon Robinson
Department of Civil and Environmental Engineering
Carleton University
Ottawa, ON, Canada
&Mohammad Khalil
Quantitative Modeling & Analysis Department
Sandia National Laboratories
Livermore, CA, United States
&Chris L. Pettit
Aerospace Engineering Department
US Naval Academy
Annapolis, MD, United States
&Dominique Poirel
Department of Mechanical and Aerospace Engineering
Royal Military College of Canada
Kingston, ON, Canada
&Abhijit Sarkar
Department of Civil and Environmental Engineering
Carleton University
Ottawa, ON, Canada
Currently at National Renewable Energy Laboratory, USASandia National
Laboratories is a multimission laboratory managed and operated by National
Technology and Engineering Solutions of Sandia, LLC., a wholly owned
subsidiary of Honeywell International, Inc., for the U.S. Department of
Energy’s National Nuclear Security Administration under contract DE-
NA-0003525.

###### Abstract

This paper investigates the applicability of a recently-proposed nonlinear
sparse Bayesian learning (NSBL) algorithm to identify and estimate the complex
aerodynamics of limit cycle oscillations. NSBL provides a semi-analytical
framework for determining the data-optimal sparse model nested within a
(potentially) over-parameterized model. This is particularly relevant to
nonlinear dynamical systems where modelling approaches involve the use of
physics-based and data-driven components. In such cases, the data-driven
components, where analytical descriptions of the physical processes are not
readily available, are often prone to overfitting, meaning that the empirical
aspects of these models will often involve the calibration of an unnecessarily
large number of parameters. While it may be possible to fit the data well,
this can become an issue when using these models for predictions in regimes
that are different from those where the data was recorded. In view of this, it
is desirable to not only calibrate the model parameters, but also to identify
the optimal compromise between data-fit and model complexity. In this paper,
this is achieved for an aeroelastic system where the structural dynamics are
well-known and described by a differential equation model, coupled with a
semi-empirical aerodynamic model for laminar separation flutter resulting in
low-amplitude limit cycle oscillations. For the purpose of illustrating the
benefit of the algorithm, in this paper, we use synthetic data to demonstrate
the ability of the algorithm to correctly identify the optimal model and model
parameters, given a known data-generating model. The synthetic data are
generated from a forward simulation of a known differential equation model
with parameters selected so as to mimic the dynamics observed in wind-tunnel
experiments.

_K_ eywords Aeroelasticity, unsteady aerodynamics, inverse problems, sparse
learning, Bayesian inference, nonlinear dynamics

##  1 Introduction

In this paper we demonstrate the applicability of the recently proposed
nonlinear sparse Bayesian learning (NSBL) algorithm [1, 2] to a single degree
of freedom (SDOF) aeroelastic oscillator that is undergoing low amplitude
limit cycle oscillations. This experimental setup has been studied extensively
through experimentation [3, 4], numerical modelling of laminar separation
flutter using high-fidelity large eddy simulations (LES) [5] and unsteady
Reynolds averaged Navier-Stokes (URANS) model [6]. Furthermore, the wind
tunnel experiments have provided a reliable test-bed for developing Bayesian
techniques for system identification and model selection for nonlinear
dynamical systems [7, 8, 9, 10]. The work in [7, 8, 9] use standard methods of
evidence-based Bayesian model selection, which allows for the systematic
comparison of a set of candidate models with varying degrees of complexity.
The model evidence as a criterion for model selection ensures a balance of
favouring models with superior average data-fit, while penalizing models that
are overly complex and thus prone to overfitting [11]. In this context, model
complexity is quantified by the KL-divergence of the parameter posterior
probability density function (pdf) from the parameter prior pdf. For
parameters where there exists little prior knowledge, it is typical to assign
non-informative priors, however the width of the distribution used for the
non-informative prior will influence the optimal complexity of the model. The
issue of sensitivity to prior width is addressed in [10], whereby the problem
is reposed as a sparse learning problem. Rather than non-informative priors,
parameters with little prior information are assigned Gaussian automatic
relevance determination (ARD) priors. The precision (inverse of the variance)
of these ARD priors are determined through evidence optimization. In this re-
framing of the inference problem, the optimal model is still quantified as
such based on the model evidence. In contrast to the previous approach, rather
than proposing an entire set of nested candidate models to determine the
optimal model complexity, the current paper approaches the problem as an
automatic discovery of the optimal sparse model nested within a single
(potentially) over-parameterized model. Herein lies an additional benefit of
approaching the problem as a sparse learning task; it is only necessary to
obtain the parameter posterior for a single model, whereas standard methods
require the calibration of each model in the candidate in order to then obtain
an estimate of the model evidence. The shortcoming of this approach lies in
the fact that the optimization process involves the use of Markov Chain Monte
Carlo (MCMC) sampling at each iteration. This is addressed in the current NSBL
framework, which removes the use of MCMC from within the optimization loop,
resulting in significantly improved computational efficiency.

The NSBL framework presented here is an extension of the sparse Bayesian
learning (SBL) also known as the relevance vector machine (RVM) algorithm [12,
13]. Both methods are motivated by the desire to avoid overfitting during
Bayesian inversion. SBL/RVM and the similar Bayesian compressive sensing (BCS)
algorithm [14] provide analytical expressions for a sparse parameter posterior
distribution owing to the analytical conveniences of the semi-conjugacy that
exists between the Gaussian likelihood functions, and Gaussian ARD priors that
are conditioned on hyperpriors that are Gamma distributions. The SBL
methodology is extended to be applicable to nonlinear-in-parameter models and
for non-Gaussian prior distributions, as these both commonly arise in
engineering applications. We provide the minimum required mathematical details
to understand the objectives of the algorithm and to provide a complete
account of all terms shown in the equations used in this paper. For the full
detailed derivation and additional details, we refer the reader to [1, 2].

##  2 Methodology: Nonlinear sparse Bayesian learning

The NSBL methodology is applicable to general nonlinear mappings,
$f:\boldsymbol{\phi}\mapsto\mbox{\bf y}$ where the model operator $f$ maps the
unknown model parameter vector $\boldsymbol{\phi}\in\mathbb{R}^{N_{\phi}}$ to
the observable model output $\mbox{\bf y}\in\mathbb{R}^{N_{y}}$. In this
specific application, $f$ represents the aeroelastic model,
$\boldsymbol{\phi}$ are the deterministic system parameters and the stochastic
parameters (relating to the model error), and y are the system output. Sensor
measurements of the system output y at discrete points in time are denoted as
$\mathcal{D}$. The likelihood function
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})$ can be computed for any
$\boldsymbol{\phi}$, using the observations $\mathcal{D}$, and these
observation may be noisy, sparse, and incomplete measurements of the system
state. The purpose of the algorithm is to obtain a data-optimal sparse
representation of $\boldsymbol{\phi}$ using Bayesian inversion, while removing
redundant parameters.

NSBL operates within the following Bayesian framework; we seek the posterior
distribution of the unknown model parameters $\boldsymbol{\phi}$ conditioned
on the data $\mathcal{D}$ and hyperparameters $\boldsymbol{\alpha}$,

$\mbox{p}(\boldsymbol{\phi}|\mathcal{D},\boldsymbol{\alpha})=\frac{\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}|\boldsymbol{\alpha})}{\mbox{p}(\mathcal{D}|\boldsymbol{\alpha})}=\frac{\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}|\boldsymbol{\alpha})}{\int\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}|\boldsymbol{\alpha})d\boldsymbol{\phi}}$
(1)

for given data and hyperparameters, the denominator, which represents the
model evidence (or marginal likelihood or type-II likelihood), is just a
normalization constant. The parameter prior
$\mbox{p}(\boldsymbol{\phi}|\boldsymbol{\alpha})$ is also conditional on the
hyperparameter. Though the objective is not to perform full hierarchical
Bayesian inference, we nevertheless define a prior on
$\mbox{p}(\boldsymbol{\alpha})$ (which is notably absent in the expression
above); this hyperparameter prior (or hyperprior) becomes relevant during the
optimization of $\boldsymbol{\alpha}$.

The following sections outline the three principal tasks involved in the NSBL
framework, as depicted in Figure 1. Namely:

  1. (i)

in section 2.1 we discuss the assignment of a hybrid prior, wherein we
distinguish between a priori relevant parameters and questionable parameters,
assigning known priors and ARD priors, respectively,

  2. (ii)

in section 2.2, we detail the incorporation of data and the physics-based
model through the construction of a Gaussian mixture model (GMM) over samples
generated from the product of the likelihood function and the known prior, and

  3. (iii)

in section 2.3 we discuss the optimization of the hyperparameters. The
derivation of various semi-analytical entities that enable the NSBL
methodology is outlined in A.

Figure 1: Summary of the main steps involved in the NSBL algorithm.

###  2.1 Hybrid prior pdf

The model parameter vector $\boldsymbol{\phi}$ is first decomposed as
$\boldsymbol{\phi}=\\{\boldsymbol{\phi}_{\alpha},\boldsymbol{\phi}_{\text{-}\alpha}\\}$,
distringuishing between the set of parameters that are known to be relevant a
priori, denoted as $\boldsymbol{\phi}_{\alpha}\in\mathbb{R}^{N_{\alpha}}$, and
the parameters that the modeller has deemed to be questionable, denoted
$\boldsymbol{\phi}_{\text{-}\alpha}\in\mathbb{R}^{N_{\phi}-N_{\alpha}}$. This
classification as questionable encompasses any parameter for which little or
no prior knowledge exists, where a non-informative prior with large support
would usually be used. The vector of questionable parameters are the set of
parameters among which we will induce sparsity, as a subset of these
parameters may be redundant. The mechanism for inducing sparsity follows SBL,
where $\boldsymbol{\phi}_{\alpha}$ is assigned a Gaussian ARD prior
$\mbox{p}(\boldsymbol{\phi}_{\alpha}|\boldsymbol{\alpha})=\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
0},\mbox{\bf A}^{-1})$. This prior is a normal distribution, whose mean vector
is an $N_{\alpha}\times 1$ zero vector, with a covariance matrix of $\mbox{\bf
A}^{-1}$, where A is the precision matrix. Following SBL [15], prior
independence of the questionable parameters $\boldsymbol{\phi}_{\alpha}$ is
assumed, hence, the precision matrix is diagonal, $\mbox{\bf
A}=\text{diag}(\boldsymbol{\alpha})$. Furthermore, each parameter
$\phi_{i}\in\boldsymbol{\phi}_{\alpha}$ has a unique variable precision
$\alpha_{i}$, such that we can write
$\mbox{p}(\phi_{i}|\alpha_{i})=\mathcal{N}(\phi_{i}|0,\alpha_{i}^{-1})$. The
hyperparameter $\alpha_{i}$ dictates the prior precision of parameter
$\phi_{i}$; where low precision (or high variance) reduces to a non-
informative prior, and conversely, a high precision (or low variance) results
in an informative prior with a mean of zero. In the limit where the precision
tends to infinity, the ARD prior becomes a Dirac delta function centered at
zero, effectively pruning the parameter. Hence, the motivation behind NSBL is
that optimally selecting $\boldsymbol{\alpha}$, can allow us to discover the
model having the optimal complexity given the available data. The optimization
criteria and methodology are presented later in section 2.3.

The joint prior pdf of $\boldsymbol{\phi}$ is summarized as [2]

$\mbox{p}(\boldsymbol{\phi}|\boldsymbol{\alpha})=\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})\mbox{p}(\boldsymbol{\phi}_{\alpha}|\boldsymbol{\alpha})=\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
0},\mbox{\bf A}^{-1}).$ (2)

This hybrid prior pdf enables sparse learning of questionable parameters in
$\boldsymbol{\phi}_{\alpha}$ through the use of an ARD prior
$\mbox{p}(\boldsymbol{\phi}_{\alpha}|\boldsymbol{\alpha})$ while incorporating
prior knowledge about parameters $\boldsymbol{\phi}_{\text{-}\alpha}$ through
an informative prior $\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$. The ARD
prior is a conditional Gaussian distribution, whose precision depends on the
hyperparameter $\boldsymbol{\alpha}$. The marginal hyperprior pdf
$\mbox{p}(\alpha_{i})$ is chosen to be a Gamma distribution. Given the
assumption of prior independence, the joint hyperprior
$\mbox{p}(\boldsymbol{\alpha})$ is written as

$\mbox{p}(\boldsymbol{\alpha})=\prod_{i=1}^{N_{\alpha}}\mbox{p}(\alpha_{i})=\prod_{i=1}^{N_{\alpha}}\text{Gamma}(\alpha_{i}|r_{i},s_{i})=\prod_{i=1}^{N_{\alpha}}\frac{s_{i}^{r_{i}}}{\Gamma(r_{i})}\alpha_{i}^{r_{i}-1}e^{-s_{i}\alpha_{i}},$
(3)

where $\text{Gamma}(\alpha_{i}|r_{i},s_{i})$ denotes a univariate Gamma
distribution parameterized by shape parameter $r_{i}>0$ and rate parameter
$s_{i}>0$. The use of a Gamma distribution as the hyperprior allows us to
enforce the requirement that the precision parameters $\boldsymbol{\alpha}$ be
positive. Furthermore, for specific combinations of shape and rate parameters,
$r_{i}$ and $s_{i}$, the Gamma function can assume many forms of informative
and non-informative priors. For instance, using values of $s_{i}\approx 1$ and
$r_{i}\approx 0$ gives a flat prior over $\alpha_{i}$, or values of
$s_{i}\approx 0$ and $r_{i}\approx 0$ gives Jeffreys prior, which is a flat
over $\log\alpha_{i}$. In fact, for reasons discussed in later sections, the
NSBL algorithm operates on the natural logarithm of the hyperparameters rather
than the hyperparameters directly. For this reason, we choose to use Jeffreys
prior for the numerical results section.

Using a univariate transformation of random variables [16], the hyperprior in
Eq. (3) becomes

$\mbox{p}(\log\boldsymbol{\alpha})=\prod_{i=1}^{N_{\alpha}}\mbox{p}(\log\alpha_{i})=\prod_{i=1}^{N_{\alpha}}\frac{\mbox{p}(\alpha_{i})}{\left|\frac{d}{d\alpha_{i}}\log\alpha_{i}\right|}=\prod_{i=1}^{N_{\alpha}}\frac{s_{i}^{r_{i}}}{\Gamma(r_{i})}\alpha_{i}^{r_{i}}e^{-s_{i}\alpha_{i}}.$
(4)

###  2.2 Gaussian mixture-model approximation

After defining the joint parameter prior distribution as in Eq. (2), we
substitute the resulting expression into the conditional posterior
distribution from Eq. (1), yielding [2]

$\mbox{p}(\boldsymbol{\phi}|\mathcal{D},\boldsymbol{\alpha})=\
\frac{\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
0},\mbox{\bf A}^{-1})}{\mbox{p}(\mathcal{D}|\boldsymbol{\alpha})}.$ (5)

Given that the ARD priors are normally distributed, NSBL constructs a GMM
approximation of the remaining terms in the numerator,
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$,
which enables the derivation of semi-analytical expressions for many entities
of interest. Obtaining expressions for the model evidence and objective
function (A.1), the parameter posterior (A.2), and the gradient and Hessian of
the objective function (A.3) makes the optimization of the hyperparameters
analytically tractable. Moreover,the use of a GMM enables the preservation of
non-Gaussianity in both the likelihood function and the known prior. The
construction involves the estimation of kernel parameters $a^{(k)}$,
$\boldsymbol{\mu}^{(k)}$, and $\boldsymbol{\Sigma}^{(k)}$ ,

$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})\approx\sum_{k=1}^{K}a^{(k)}\mathcal{N}(\boldsymbol{\phi}|\boldsymbol{\mu}^{(k)},\boldsymbol{\Sigma}^{(k)}),$
(6)

where $K$, $a^{(k)}$, and
$\mathcal{N}(\boldsymbol{\phi}|\boldsymbol{\mu}^{(k)},\boldsymbol{\Sigma}^{(k)})$
are the total number of kernels, the kernel coefficient ($a^{(k)}>$0), and a
Gaussian pdf with mean vector $\boldsymbol{\mu}^{(k)}\in\mathbb{R}^{N_{\phi}}$
and covariance matrix $\boldsymbol{\Sigma}^{(k)}\in\mathbb{R}^{N_{\phi}\times
N_{\phi}}$ [2]. For a Gaussian likelihood function and a Gaussian known prior,
this reduces to the case of SBL or RVM [15], and a single kernel is
sufficient. Otherwise, except in the case of a Laplace approximation [17] of
the likelihood function times the known prior, multiple kernels will generally
be required. The construction of the GMM typically involves the use of MCMC in
order to generate samples from the arbitrary distribution, which requires
repeated function evaluations for different samples of the unknown parameter
vector. The model itself operates as a black-box, thus, the analytical form of
the model does not need to be known; the model is only needed in order to
compute the likelihood function. Once samples have been generated from the
posterior distribution, the estimation of the kernel parameters in Eq. (6) can
be performed numerically using methods such as kernel density estimation (KDE)
or expectation maximization (EM) [18]. Since the construction of the GMM
involves numerous forward solves of the model, this step is the most
computationally demanding component of the algorithm. Notably, the GMM itself
is independent of the hyperparameters, thus this process only needs to be
performed once at the onset, and does not need to be repeated during the
optimization of the hyperparameters.

###  2.3 Sparse learning optimization problem

The critical step in the NSBL algorithm is the optimization of the
hyperparameter, $\boldsymbol{\alpha}$. Within a hierarchical Bayesian
framework, we seek a point estimate for the hyperparameters, rather than
obtaining posterior estimates thereof. As in SBL, we perform type-II maximum
likelihood, seeking the values $\boldsymbol{\alpha}$, which maximize the
hypeperparameter posterior,

$\mbox{p}(\boldsymbol{\alpha}|\mathcal{D})=\frac{\mbox{p}(\mathcal{D}|\boldsymbol{\alpha})\mbox{p}(\boldsymbol{\alpha})}{\mbox{p}(\mathcal{D})}\propto\mbox{p}(\mathcal{D}|\boldsymbol{\alpha})\mbox{p}(\boldsymbol{\alpha}).$
(7)

for a fixed set of data $\mathcal{D}$, the denominator $\mbox{p}(\mathcal{D})$
is a normalization constant that is analytically intractable in general. Thus
for optimization, we consider only the numerator which is the product of the
model evidence $\mbox{p}(\mathcal{D}|\boldsymbol{\alpha})$ and the hyperprior
$\mbox{p}(\boldsymbol{\alpha})$ . The type-II maximum a posteriori (MAP-II)
estimate $\boldsymbol{\alpha}^{\text{map}}$ can therefore be posed as [18]

$\boldsymbol{\alpha}^{\text{map}}=\operatorname*{arg\,max}_{\boldsymbol{\alpha}}\\{\mbox{p}(\boldsymbol{\alpha}|\mathcal{D})\\}=\operatorname*{arg\,max}_{\boldsymbol{\alpha}}\\{\mbox{p}(\mathcal{D}|\boldsymbol{\alpha})\mbox{p}(\boldsymbol{\alpha})\\}.$
(8)

Since the natural logarithm is a strictly increasing function, we re-write the
optimization problem in terms of
$\log\mbox{p}(\boldsymbol{\alpha}|\mathcal{D})$, which facilitates the
derivation of the semi-analytical expressions of the gradient and Hessian.
Furthermore, as noted above, the optimization is performed with respect to
$\log\boldsymbol{\alpha}$ instead of $\boldsymbol{\alpha}$. This helps to
account for the potentially widespread difference in the scales of the prior
precision between relevant parameters (whose $\alpha_{i}^{\text{map}}$ are
finite) and irrelevant parameters (whose $\alpha_{i}^{\text{map}}$ tend to
infinity). This has the added benefit of automatically enforcing the
positivity constraint of $\boldsymbol{\alpha}$ throughout the optimization
procedure. Restating Eq. (8) as the optimization of the log-evidence with
respect to the log of the hyperparameters $\boldsymbol{\alpha}$, the objective
function $\mathcal{L}(\log\boldsymbol{\alpha})$ becomes

$\displaystyle\log\boldsymbol{\alpha}^{\text{map}}$
$\displaystyle=\operatorname*{arg\,max}_{\log\boldsymbol{\alpha}}\\{\mathcal{L}(\log\boldsymbol{\alpha})\\}$
$\displaystyle=\operatorname*{arg\,max}_{\log\boldsymbol{\alpha}}\\{\log\mbox{p}(\log\boldsymbol{\alpha}|\mathcal{D})\\}$
$\displaystyle=\operatorname*{arg\,max}_{\log\boldsymbol{\alpha}}\\{\log\hat{\mbox{p}}(\mathcal{D}|\log\boldsymbol{\alpha})+\log\mbox{p}(\log\boldsymbol{\alpha})\\}.$
(9)

The objective function is the sum of the estimate of the model evidence and
the hyperprior (where constant terms that are independent of $\alpha_{i}$ are
discarded). As the model evidence itself will be analytically intractable in
general, it is replaced by the estimate
$\hat{\mbox{p}}(\mathcal{D}|\log\boldsymbol{\alpha})$ in Eq. (21) that is
available in terms of the $K$ kernels of the GMM. Substituting Eq. (4) for the
hyperprior gives

$\displaystyle\mathcal{L}(\log\boldsymbol{\alpha})$
$\displaystyle=\log\hat{\mbox{p}}(\mathcal{D}|\log\boldsymbol{\alpha})+\sum_{i=1}^{N_{\alpha}}\left(r_{i}\log\alpha_{i}-s_{i}\alpha_{i}\right).$
(10)

In this form, it becomes clear that when Jeffreys prior ($s_{i}\approx 0$,
$r_{i}\approx 0$) is used, the objective function reduces to the log-evidence,
resulting in the common procedure known as the type-II maximum likelihood
estimate or Emperical Bayes method [18], which may lead to a non-convex
optimization problem [17]. The numerical examples that follow will illustrate
the possibility of non-unique optima in the objective function. The
possibility of multiple optima can be addressed using a multi-start
optimization algorithm, where the optimization routine is initiated from an
array of different coordinates of $\log\boldsymbol{\alpha}$ in an effort to
discover all local optima in order to determine the global optima. Sequential
estimates of $\\{\log\boldsymbol{\alpha}_{j}\\}$ are obtained starting from
the initial iterate $\\{\log\boldsymbol{\alpha}_{0}\\}$, according to [2]

$\log\boldsymbol{\alpha}_{j+1}=\log\boldsymbol{\alpha}_{j}+\beta_{j}\mbox{\bf
p}_{j},$ (11)

where $\beta_{j}$ is the step-length [19]. The optimization itself is
expedited by the ability to derive the gradient vector $\mbox{\bf
J}(\log\boldsymbol{\alpha}_{i})$ and Hessian matrix $\mbox{\bf
H}(\log\boldsymbol{\alpha}_{i})$ of the objective function from Eq. (10), as
outlined in Eq. (29) and Eq. (31), respectively. This permits the convenient
calculation of $\mbox{\bf p}_{j}$ as the solution to

$\mbox{\bf H}(\log\boldsymbol{\alpha}_{j})\mbox{\bf p}_{j}=-\mbox{\bf
J}(\log\boldsymbol{\alpha}_{j}).$ (12)

The specific method for determining $\mbox{\bf p}_{j}$ (e.g., modified Newton
method, trust-region Newton method) must consider that the Hessian is not
guaranteed to be a positive definite matrix [19].

The optimization of the objective function (the log hyperparameter posterior)
with respect to $\log\boldsymbol{\alpha}$, provides the MAP estimate
$\log\boldsymbol{\alpha}^{\text{map}}$. As discussed previously, a large value
of $\log\alpha^{\text{map}}_{i}$ indicates high prior precision of parameter
$\phi_{i}$ and implies it is redundant. However, it is difficult to quantify
high precision directly as this is highly parameter dependent. To remove the
scale-dependence, we leverage another concept from SBL [15], wherein a
relevance indicator is defined for each parameter $\phi_{i}$, according to
$\gamma_{i}=1-\alpha_{i}P_{ii}$, where, $P_{ii}$ is the $i$th diagonal entry
of the posterior covariance matrix and $\alpha_{i}$ is the $i$th diagonal
entry in the prior precision matrix, A. In NSBL, it is possible to define a
similar metric, however, as with the semi-analytical entities described
previously, a relevance indicator will be defined for each of the $K$ Gaussian
kernels. This can be re-written as a ratio of prior to posterior precision
[2],

$\gamma_{i}^{(k)}=1-\frac{\alpha_{i}}{(P_{ii}^{(k)})^{-1}}\in[0,1],$ (13)

hence, the relevance indicator provides a normalized metric on the scale of 0
to 1, where a value of 0 indicates irrelevance, while a value of 1 implies
relevance. The intuition behind this conclusion is that for a parameter where
the ratio of prior-to-posterior precision is close to unity, the posterior
precision is dictated by the prior, rather than the likelihood, suggesting
that the parameter does not learn from the data. The converse holds as well,
where parameters that do learn from the data will tend to have higher
posterior precision compared to the prior precision, hence the expression for
the relevance indicator will approach a value of 1. Since we perform this one
a kernel-by-kernel basis, to summarize the relevance indicators for a given
parameter across all kernels in the GMM, we propose the use of a root-mean-
square value of the relevance indicator [2]

$\gamma_{i}^{\text{rms}}=\left(\frac{1}{K}\sum_{k=1}^{K}(\gamma_{i}^{(k)})^{2}\right)^{1/2}=\left(\frac{1}{K}\sum_{k=1}^{K}\left(1-\alpha_{i}P_{ii}^{(k)}\right)^{2}\right)^{1/2}.$
(14)

##  3 Application to aeroelastic oscillator

In this section, we demonstrate the performance of NSBL for identifying
sparsity in the unknown parameters of nonlinear stochastic differential
equations. We consider a single degree-of-freedom pitching airfoil undergoing
low-amplitude limit cycle oscillations. The structural model in Eq. (15a)
consists of a typical mass-spring-damper system augmented by a cubic stiffness
and a model for dry friction. The structural dynamics are coupled with a semi-
empirical aerodynamic model in Eq. (15b), which we will refer to as an
unsteady generalized Duffing-Van der Pol model [5]. The aerodynamic moment
coefficient ($C_{M}$) modelled as a first-order ODE, which allows for
unsteadiness in the aerodynamics, driven by the parameter $B$. The model
features a polynomial expansion of the pitch angle ($\theta$) and pitch rate
($\dot{\theta}$) to attempt to capture nonlinearity in the aerodynamic moment
coefficient. Finally, it contains a model error term to capture the
discrepancy between the modelled physics and the true, but ultimately unknown,
aerodynamic loads. The resulting system of coupled ODEs are written in terms
of non-dimensional time ($\tau$) as [8]

$\displaystyle
I\ddot{\theta}+C\dot{\theta}+K\theta+C_{nl}\text{sign}(\dot{\theta})+K_{nl}\theta^{3}$
$\displaystyle=\frac{1}{2}\rho U^{2}csC_{M}$ (15a)
$\displaystyle\frac{\dot{C}_{M}}{B}+C_{M}$
$\displaystyle=a_{0}+a_{1}\theta+a_{2}\dot{\theta}+a_{3}\theta^{2}+a_{4}\theta\dot{\theta}+a_{5}\dot{\theta}^{2}+\ldots+\frac{c_{6}}{B}\ddot{\theta}+\sigma\xi(\tau)$
(15b)

where the structural parameters $I$, $C$, $C_{nl}$ $K$ and $K_{nl}$, and
aerodynamic parameters $\rho$, $U$, $c$, $s$, and $c_{6}$ are known precisely,
and coefficients $B$, $\sigma$, and aerodynamic coefficients $a_{i}$ are to be
estimated.

In the current study, synthetic LCO data is used, in order to validate the use
of NSBL in discovering the data-optimal sparse model. In this scenario, the
data-generating model and the candidate model have the same underlying
analytical form, though the candidate model may contain extra parameters that
were not present (i.e. were equal to 0) when generating the data. Hence the
algorithm should ideally recover the parameters used in the data-generating
model, and identify any parameters that are not used to generate the data as
irrelevant. The synthetic LCO from which the measurements of the pitch are
generated have similar frequency and amplitude as are observed in the wind-
tunnel experimental setup. Synthetic measurements are generated by recording
the pitch deflection (representing only a partial measurement of the state).
Measurements are recorded at a frequency of 1000Hz (which is consistent with
the temporal density of observations available experimentally). These
measurements are then corrupted by additive Gaussian noise with a noise
strength that again mimics the noise levels in the wind-tunnel experiments. As
in the physical experiments, a low-pass filter is used to remove contributions
with a frequency above 25Hz (see the bottom pannel of Figure 2). In practice,
the sensors capture contributions from the wind-tunnel motor at 30Hz. The
dominant frequency of the system is 3.25Hz, so many super-harmonics remain
below the cut-off frequency.

The data-generating model is shown in Eq. (15). The structural model from Eq.
(15a) is rewritten in a standardized form, replacing the system parameters by
coefficients $c_{1},\ldots,c_{6}$. The aerodynamic model is also rewritten,
retaining only a subset of the terms listed in Eq. (15b), and replacing the
aerodynamic parameters by coefficients $e_{1},\ldots,e_{4}$. The terms
retained are limited to: (i) linear stiffness $e_{1}\theta$, (ii) linear
damping $e_{2}\dot{\theta}$, (iii) nonlinear Duffing-type stiffness
$e_{3}\theta^{3}$, and (iv) Van der Pol-type nonlinear damping term
$e_{4}\theta^{2}\dot{\theta}$. Furthermore, unsteadiness is introduced through
the coefficient $B$, and a random forcing is introduced through
$\sigma\xi(\tau)$, where $\xi(\tau)$ is a white noise process

$\displaystyle\ddot{\theta}$
$\displaystyle=c_{1}\text{sign}(\dot{\theta})+c_{2}\theta+c_{3}C_{\text{M}}+c_{4}\dot{\theta}+c_{5}\theta^{3},$
(16a) $\displaystyle\frac{\dot{C}_{\text{M}}}{B}+C_{\text{M}}$
$\displaystyle=e_{1}\theta+e_{2}\dot{\theta}+e_{3}\theta^{3}+e_{4}\theta^{2}\dot{\theta}+\frac{c_{6}}{B}\ddot{\theta}+\sigma\xi(\tau).$
(16b)

The structural parameters $c_{1},\ldots,c_{6}$, from Eq. (16a) and the
aeroedynamic coefficients $B,e_{1},\ldots,e_{4}$ and the strength of the
random forcing $\sigma$ in Eq. (16b) are summarized in Table 1. The
observations, $\mathcal{D}$ are shown in Figure 2. Next, we exploit the NSBL
algorithm from Section 2 to demonstrate the algorithm’s ability to correctly
identify the relevance or irrelevance of parameters in one-, two-, and four-
dimensional sparse learning exercises.

Structural parameters |  Aerodynamic parameters
---|---
$c_{1}=-6.875\times 10^{-5}$ |  $B=2.000\times 10^{-1}$
$c_{2}=-2.038\times 10^{-2}$ |  $e_{1}=-1.250$
$c_{3}=-3.819\times 10^{-2}$ |  $e_{2}=-1.000$
$c_{4}=-7.275\times 10^{-3}$ |  $e_{3}=1.000\times 10^{2}$
$c_{5}=1.824\times 10^{-1}$ |  $e_{4}=-5.000\times 10^{2}$
$c_{6}=-2.507\times 10^{-1}$ |  $\sigma=2.000\times 10^{-3}$
Table 1: Structural and aerodynamic model parameters used to generate the
synthetic data. Figure 2: Synthetic pitch measurements and its PSD. The data
used for ARD computations is shown in red. (Reprinted from Sandhu et al. [10],
with permission from Elsevier).

###  3.1 Unidimensional sparse learning

We first consider a unidimensional sparse learning problem, meaning the
questionable parameter vector $\boldsymbol{\phi}_{\alpha}$ and the associated
hyperparameter vector $\boldsymbol{\alpha}$ are both scalar entities. Given
the measurements $\mathcal{D}$, the proposed model for $C_{\text{M}}$ is
chosen to be the same as the data-generating model shown in Eq. (16b). The
structural equation of motion in Eq. (16a) is assumed to be known for the sake
of inverse modelling. The cubic aerodynamic stiffness parameter $e_{3}$ (or
term $e_{3}\theta^{3}$) in Eq. (16b) is treated as a questionable parameter,
whose relevance to the LCO aerodynamics needs to be determined using the NSBL
algorithm. As per Section 2.1, $\boldsymbol{\phi}$ is decomposed into the
questionable parameter $\boldsymbol{\phi}_{\alpha}=\\{e_{3}\\}$ and the a
priori relevant parameter vector $\boldsymbol{\phi}_{\text{-}\alpha}$ =
$\\{B$, $e_{1}$, $e_{2}$, $e_{4}$, $\sigma\\}$. The prior pdf
$\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$ is known, while $e_{3}$ is
assigned an ARD prior $\mathcal{N}(e_{3}|0,\alpha^{-1})$. Table 2 summarizes
this unidimensional sparse learning setup, including the known prior pdf of
$\boldsymbol{\phi}_{\text{-}\alpha}$. Note that $\mathcal{LN}(.|r,s)$
represents a log-normal distribution with median at $r$ and coefficient of
variation of $s$; and $\mathcal{U}(.|e,f)$ represents a uniform distribution
with lower bound $e$, and upper bound $f$.

Aerodynamic model | $\displaystyle\frac{\dot{C}_{\text{M}}}{B}+C_{\text{M}}=e_{1}\theta+e_{2}\dot{\theta}+e_{3}\theta^{3}+e_{4}\theta^{2}\dot{\theta}+\frac{c_{6}}{B}\ddot{\theta}+\sigma\xi(\tau)$
---|---
$\boldsymbol{\phi}$ decomposition |  $\boldsymbol{\phi}_{\alpha}$ = $\\{e_{3}\\}$ , $\boldsymbol{\phi}_{\text{-}\alpha}$ = $\\{B$, $e_{1}$, $e_{2}$, $e_{4}$, $\sigma\\}$
ARD prior, $\mbox{p}(\boldsymbol{\phi}_{\alpha}|\alpha)$ | $\mathcal{N}(e_{3}|0,\alpha^{-1})$
Known prior, $\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$ | $\mathcal{LN}(B|0.2,0.5)\mathcal{U}(e_{1}|{-2},0)\mathcal{U}(e_{2}|{-2},0)\mathcal{U}(e_{4}|{-600},0)\mathcal{LN}(\sigma|0.002,0.5)$
Table 2: Unidimensional sparse learning setup where parameter $e_{3}$ is
treated as questionable.

The likelihood function $\mbox{p}(\mathcal{D}|\boldsymbol{\phi})$ is computed
using the extended Kalman filter (EKF) as outlined in [7, 20] to propagate the
joint state pdf through the nonlinear model. Given the known prior pdf
$\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$ from Table 2, the unnormalized
pdf
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$
is sampled using the MCMC sampler. A total of 5000 stationary (i.e. post burn-
in) MCMC samples are generated from the six-dimensional pdf
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$.
These 5000 samples possess some correlation due to the Markovian nature of
MCMC sampling. In an effort to alleviate the effect of correlation, every
tenth sample is extracted to produce 500 independent and identically
distributed (iid) samples from
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$.
These 500 iid samples are then decomposed to into 10 sets of 50 iid samples
each. These iid sets are used to investigate the finite sample properties of
NSBL algorithm. Note that this sampling of
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$
is only required to be performed once since both
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})$ and
$\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$ remain unchanged during sparse
learning. Moreover, generating 5000 stationary samples from a six-dimensional
pdf is computationally cheap using MCMC algorithms such as delayed rejection
adaptive metropolis (DRAM) [21] and transitional MCMC [22]

Next, a multivariate KDE approximation is built for
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$
using each of the 10 sets of 50 iid samples. When using a Gaussian kernel, the
KDE approximation resembles the kernel-based approximation in Eq. (6).
Henceforth, the number of kernels $K$ is the same as the number of samples,
which is 50 for the current case. The coefficient $a^{(k)}$ are all equal to
one since the Gaussian kernels in KDE are all equally weighted. The mean
vector $\boldsymbol{\mu}^{(k)}$ for each kernel is equal to the corresponding
iid sample values (vector). The covariance matrix $\boldsymbol{\Sigma}^{(k)}$
is computed automatically using Scott’s rule for estimating KDE bandwidth [23]
and is the same for all 50 kernels. Figure 3 shows the marginal parameter pdf
pertaining to each of the 10 instances of KDE approximation of
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$.
The marginal pdfs show a reasonable variation across multiple instances
considering only 50 iid samples were used.

Figure 3: Marginal parameter pdf obtained using 10 instances of 50 iid samples
from
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$.
The baseline (converged) pdf is obtained by using all of the 500 iid samples.

Figure 4 shows the objective function $\mathcal{L}(\log\alpha)$ and the log-
evidence $\log\hat{\mbox{p}}(\mathcal{D}|\log\alpha)$ for varying $\log\alpha$
value, pertaining to each instance of the KDE approximation. There is little-
to-no variability observed in the objective function and log-evidence across
the different instances of KDE. The objective function is computed using
hyperprior parameters of $\log r_{i}$ = $\log s_{i}$ = -6.0 in Eq. (10). The
similarity between the two leftmost panels implies the negligible effect of
hyperprior in $\mathcal{L}(\log\alpha)$ for the range of $\log\alpha$ values
shown. Also, the variation in $\mathcal{L}(\log\alpha)$ across different KDE
instances is much less than that for the marginal pdfs shown in Figure 3. This
behaviour of $\mathcal{L}(\log\alpha)$ indicates that the relevance property
of a parameter is more robust to sampling variability than the KDE
approximation of
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$.

Figure 4: NSBL results for the relevance determination of an a priori relevant
parameter $e_{3}$.

The optimization of the objective function ($\mathcal{L}(\log\alpha)$) and the
associated value of the relevance indicator ($\gamma^{\text{rms}}$) are shown
as a function of Newton iteration count in the two right-most panels of Figure
4, respectively. The $\log\alpha$ iterates during the multistart Newton
iteration, and the corresponding value of the relevance indicator
$\gamma_{i}^{\text{rms}}$ from Eq. (14) are shown. The $\log\alpha$ iterates
and $\gamma_{i}^{\text{rms}}$ are shown for each of the 10 KDE instances of
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$.
The Newton iteration is initiated from three $\log\alpha$ values of -20.0,
-10.0 and 0.0. The analytical solution for the gradient and Hessian of
$\mathcal{L}(\log\alpha)$ employed in the Newton iteration was first validated
using a finite-difference scheme (results not reported here for brevity). All
the multistart Newton iterations converge to a unique optimum of
$\alpha^{\text{map}}$ = -9.1. The relevance indicator $\gamma^{\text{rms}}$
converges to a value of one, implying the strong relevance of parameter
$e_{3}$ to the LCO aerodynamics. This relevance of $e_{3}$ is insensitive to
the choice of tolerance $\gamma^{\text{tol}}$ as the optimal relevance
indicator is very close to one. The posterior pdf of $\boldsymbol{\phi}$ can
be obtained as per A.2 using the optimal ARD prior
$\mathcal{N}(e_{3}|0,\text{exp}(9.1))$. Since this optimal prior has large
variance, the posterior parameter pdf obtained from NSBL is the same as that
obtained by sampling
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$
(the exact pdf shown in Figure 3). In other words, the optimal ARD prior did
not bias the posterior pdf of $e_{3}$.

Next, consider the unidimensional sparse learning setup in Table 3 where the
aerodynamic model is augmented by a higher-order nonlinear stiffness term
$e_{5}\theta^{5}$. Here, parameter $e_{5}$ is treated as questionable. Given
that this term was not used to generate the data, we should expect the $e_{5}$
parameter to be deemed irrelevant by NSBL.

Aerodynamic model | $\displaystyle\frac{\dot{C}_{\text{M}}}{B}+C_{\text{M}}=e_{1}\theta+e_{2}\dot{\theta}+e_{3}\theta^{3}+e_{4}\theta^{2}\dot{\theta}+e_{5}\theta^{5}+\frac{c_{6}}{B}\ddot{\theta}+\sigma\xi(\tau)$
---|---
$\boldsymbol{\phi}$ decomposition |  $\boldsymbol{\phi}_{\alpha}$ = $\\{e_{5}\\}$ , $\boldsymbol{\phi}_{\text{-}\alpha}$ = $\\{B$, $e_{1}$, $e_{2}$, $e_{3}$, $e_{4}$, $\sigma\\}$
ARD prior, $\mbox{p}(\boldsymbol{\phi}_{\alpha}|\boldsymbol{\alpha})$ | $\mathcal{N}(e_{5}|0,\alpha^{-1})$
Known prior, $\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$ | $\mathcal{LN}(B|0.2,0.5)\,\mathcal{U}(e_{1}|{-2},0)\,\mathcal{U}(e_{2}|{-2},0)\,\mathcal{U}(e_{3}|{-250},250)\mathcal{U}(e_{4}|{-600},0)\,\mathcal{LN}(\sigma|0.002,0.5)$
Table 3: Unidimensional sparse learning setup where parameter $e_{5}$ is
treated as questionable.

The two left-most pannels of Figure 5 show the log-evidence
$\log\hat{\mbox{p}}(\mathcal{D}|\log\alpha)$ and the objective function
$\mathcal{L}(\log\alpha)$ using the hyperprior parameters of $\log r_{i}$ =
$\log s_{i}$ = -6. Unlike the previous case, both these entities possess large
sampling variability across multiple KDE instances of
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$.
However, the relevance of $e_{5}$ is determined by the optimum of
$\mathcal{L}(\log\alpha)$, which remains same across multiple KDE instances.
This finite sample property of $\mathcal{L}(\log\alpha)$ reiterates the notion
that relevance is well-determined using limited iid samples from
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$,
even when the KDE instances possess large sampling variability. Also notice
that, unlike log-evidence, the objective function $\mathcal{L}(\log\alpha)$ is
free from flat regions for large $\log\alpha$ values. The Hessian matrix is
singular in flat regions. Therefore, the presence of flat regions in log-
evidence makes it unsuitable for the application of Newton’s method. The
absence of flat regions in $\mathcal{L}(\log\alpha)$ is due to the influence
of hyperprior $\mbox{p}(\alpha)$ for large $\log\alpha$ values. This ensures
$\mathcal{L}(\log\alpha)$ is amenable to optimization by Newton’s method
wherein the Hessian is non-singular at all times. The shape and rate
parameters of the Gamma hyperprior must have positive values, hence to assign
Jeffreys prior, they are given values that approach zero from above. If these
values are sufficiently small (here we use $\log r_{i}$ = $\log s_{i}$ =
-6.0), they should have a negligible effect of the relevance of parameters,
but from Eq. (10), it is understood that they will help regularize the sparse
learning optimization problem for large values of $\boldsymbol{\alpha}$.

Figure 5: NSBL results for the relevance determination of an a priori
irrelevant parameter $e_{5}$.

The two rightmost plots of Figure 5 show the $\log\alpha$ iterates and the RMS
relevance indicator $\gamma^{\text{rms}}$; when initiating at $\log\alpha$
values of -15.0, -5.0 and 5.0. The multistart Newton iteration converges to a
unique optimum of $\alpha^{\text{map}}$ = 0 and $\gamma^{\text{rms}}$ = 0 for
all multistart iterations and for all KDE instances of
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$.
Note that since the objective function is relatively flat close to the
optimum, the estimated optimum can vary with a varying function tolerance set
for terminating the Newton iteration. However, the key entity to monitor is
the relevance indicator. Even with a small variation in the optimum
hyperparameter value, the relevance indicator converges to the value of zero.
According to Eq. (14), a $\gamma^{\text{rms}}$ value of zero implies
irrelevance since the posterior pdf is entirely dictated by the Dirac-delta
ARD prior centered at zero. This change in the posterior pdf of $e_{5}$
following NSBL is shown in the right pannel of Figure 6.

Figure 6: NSBL results for the relevance determination of an a priori
irrelevant parameter $e_{5}$.

Figure 6 shows the posterior pdf of both nonlinear stiffness parameters
$e_{3}$ and $e_{5}$ before and after the inclusion of optimal ARD prior
$\mathcal{N}(e_{5}|0,\text{exp}(0.0))$. Notice that the posterior pdf of
$e_{3}$ following the removal of $e_{5}$ has a lower uncertainty and is in
close agreement with the true value. This behaviour can be explained by
realizing that $e_{3}$ and $e_{5}$ both quantify the nonlinear aerodynamic
stiffness. When both the parameters are present in the model, overfitting
occurs. When $e_{5}$ is removed, overfitting is remedied, and $e_{3}$ is
estimated with an increased accuracy using $\mathcal{D}$. This numerical
exercise demonstrates the benefit of sparse learning in preventing overfitting
and enabling robust probabilistic predictions outside the regime of
measurements. The two unidimensional cases studied here instills confidence in
NSBL as a practical tool to perform sparse learning among physics-based
models. Next, we study two-dimensional sparse learning cases to gain a more
in-depth understanding of NSBL.

###  3.2 Two-dimensional sparse learning

In this section, we consider two-dimensional sparse learning cases where the
number of questionable parameters is two (i.e. $N_{\alpha}$ = 2). The
hyperparameter vector is $\boldsymbol{\alpha}$ =
$\\{\alpha_{1},\alpha_{2}\\}$. NSBL is executed with the same settings as
Section 3.1. Hence, the likelihood function
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})$ is computed using EKF, followed by
the generation of 500 iid samples from
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$
using the DRAM MCMC algorithm. NSBL is executed for each of the 10 KDE
instances of
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$
containing 50 iid samples each. Also, the hyperprior parameters are chosen as
$\log r_{i}$ = $\log s_{i}$ = -6.0 for all the cases reported in this section.
We consider the three possible combinations of relevant and irrelevant
parameters: where the questionable parameters are both relevant, where one is
relevant and the other is irrelevant, and where both are irrelevant.

We first consider the Bayesian setup in Table 4 where the proposed model is
same as the data-generating model from Eq. (16b), and parameters $e_{3}$ and
$e_{4}$ are treated as questionable. Figure 7 shows the
$\log\boldsymbol{\alpha}$ iterates during Newton iteration initiated from
{-20, -20}, {-5, -20}, {-20, -5} and {-5, -5}, for each of the 10 KDE
instances. Notice that all of the multistart Newton iterations converges to a
unique optimum of $\log\boldsymbol{\alpha}^{\text{map}}$ = {-9.1, -12.4}.
Also, convergence to a unique optimum for varying KDE instances of
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$
demonstrates that the curvature of objective function
$\mathcal{L}(\log\boldsymbol{\alpha})$ remains fairly insensitive to sampling
variability. The two rightmost panels of Figure 7 show the variation in RMS
relevance indicator $\gamma^{\text{rms}}$ for parameter $e_{3}$ and $e_{4}$,
respectively. Both the relevance indicators converge to one, implying the
relevance of both the parameters to the LCO physics. NSBL converges in less
than 10 Newton iterations, demonstrating the power of the Hessian informed
optimizer.

Aerodynamic model | $\displaystyle\frac{\dot{C}_{\text{M}}}{B}+C_{\text{M}}=e_{1}\theta+e_{2}\dot{\theta}+e_{3}\theta^{3}+e_{4}\theta^{2}\dot{\theta}+\frac{c_{6}}{B}\ddot{\theta}+\sigma\xi(\tau)$
---|---
$\boldsymbol{\phi}$ decomposition |  $\boldsymbol{\phi}_{\alpha}$ = $\\{e_{3},e_{4}\\}$ , $\boldsymbol{\phi}_{\text{-}\alpha}$ = $\\{B$, $e_{1}$, $e_{2}$, $\sigma\\}$
ARD prior, $\mbox{p}(\boldsymbol{\phi}_{\alpha}|\alpha)$ |  $\mathcal{N}(e_{3}|0,\alpha_{1}^{-1})$ $\mathcal{N}(e_{4}|0,\alpha_{2}^{-1})$
Known prior, $\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$ | $\mathcal{LN}(B|0.2,0.5)\mathcal{U}(e_{1}|{-2},0)\mathcal{U}(e_{2}|{-2},0)\mathcal{LN}(\sigma|0.002,0.5)$
Table 4: Two-dimensional sparse learning setup where parameters $e_{3}$ and
$e_{4}$ are treated as questionable. Figure 7: NSBL results for the relevance
determination of parameter $e_{3}$ and $e_{4}$.

Next, consider the Bayesian setup in Table 5 where parameters $e_{3}$ and
$e_{5}$ are treated as questionable. The data-generating model in Eq. (16b)
contains the cubic aerodynamic stiffness coefficient $e_{3}$, while the fifth-
order aerodynamic stiffness coefficient $e_{5}$ is absent. Figure 8 shows the
$\log\boldsymbol{\alpha}$ iterates pertaining to starting
$\log\boldsymbol{\alpha}$ values of {-20, -20}, {-5, -20}, {-20, -5} and {-5,
-5}. The multistart Newton iteration converges to a unique optimum of
$\log\boldsymbol{\alpha}^{\text{map}}$ = {-9.1, 0.0}. Recall that these are
the two values obtained when these parameters were studied independently in
section 3.1. The two right-most panels of Figure 8 shows the RMS relevance
indicators for parameter $e_{3}$ and $e_{5}$, respectively. Parameter $e_{3}$
is determined to be relevant since $\gamma^{\text{rms}}_{1}$ converges to
exactly one. On the contrary, parameter $e_{5}$ is determined to be irrelevant
as $\gamma^{\text{rms}}_{2}$ converges to zero. These conclusion are true for
all multistart locations and for all KDE instances of
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$.

Aerodynamic model | $\displaystyle\frac{\dot{C}_{\text{M}}}{B}+C_{\text{M}}=e_{1}\theta+e_{2}\dot{\theta}+e_{3}\theta^{3}+e_{4}\theta^{2}\dot{\theta}+e_{5}\theta^{5}+\frac{c_{6}}{B}\ddot{\theta}+\sigma\xi(\tau)$
---|---
$\boldsymbol{\phi}$ decomposition |  $\boldsymbol{\phi}_{\alpha}$ = $\\{e_{3},e_{5}\\}$ , $\boldsymbol{\phi}_{\text{-}\alpha}$ = $\\{B$, $e_{1}$, $e_{2}$, $e_{4}$, $\sigma\\}$
ARD prior, $\mbox{p}(\boldsymbol{\phi}_{\alpha}|\alpha)$ |  $\mathcal{N}(e_{3}|0,\alpha_{1}^{-1})$ $\mathcal{N}(e_{5}|0,\alpha_{2}^{-1})$
Known prior, $\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$ | $\mathcal{LN}(B|0.2,0.5)\mathcal{U}(e_{1}|{-2},0)\mathcal{U}(e_{2}|{-2},0)\mathcal{U}(e_{4}|{-600},0)\mathcal{LN}(\sigma|0.002,0.5)$
Table 5: Two-dimensional sparse learning setup where parameters $e_{3}$ and
$e_{5}$ are treated as questionable. Figure 8: NSBL results for the relevance
determination of parameter $e_{3}$ and $e_{5}$.

Next, consider the Bayesian setup in Table 6, where parameters $e_{5}$ and
$e_{6}$ are treated as questionable. Both these parameters were absent from
the data-generating model in Eq. (16b). Figure 9 shows the
$\log\boldsymbol{\alpha}$ iterates for the multistart Newton iteration
beginning from {-20, -20}, {-5, -20}, {-20, -5} and {-5, -5}. The Newton
iterations converge to two different optimums of {0.0, -20.5} and {0.0, 0.0}.
Based on the $\mathcal{L}(\log\boldsymbol{\alpha})$ values at these optimums,
{0.0, -20.5} was found to the global optimum (for the range of $\alpha_{1}$
and $\alpha_{2}$ considered here). This ability of the multistart Newton
iteration to capture multiple optima was further validated using the case of
bimodal likelihood functions.

Aerodynamic model | $\displaystyle\frac{\dot{C}_{\text{M}}}{B}+C_{\text{M}}=e_{1}\theta+e_{2}\dot{\theta}+e_{3}\theta^{3}+e_{4}\theta^{2}\dot{\theta}+e_{5}\theta^{5}+e_{6}\theta^{4}\dot{\theta}+\frac{c_{6}}{B}\ddot{\theta}+\sigma\xi(\tau)$
---|---
$\boldsymbol{\phi}$ decomposition |  $\boldsymbol{\phi}_{\alpha}$ = $\\{e_{5},e_{6}\\}$ , $\boldsymbol{\phi}_{\text{-}\alpha}$ = $\\{B$, $e_{1}$, $e_{2}$, $e_{3}$, $e_{4}$, $\sigma\\}$
ARD prior, $\mbox{p}(\boldsymbol{\phi}_{\alpha}|\alpha)$ |  $\mathcal{N}(e_{5}|0,\alpha_{1}^{-1})$ $\mathcal{N}(e_{6}|0,\alpha_{2}^{-1})$
Known prior, $\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$ | $\mathcal{LN}(B|0.2,50)\mathcal{U}(e_{1}|{-2},0)\mathcal{U}(e_{2}|{-2},0)\mathcal{U}(e_{3}|{-250},{250})\mathcal{U}(e_{4}|-1\times 10^{4},0)\mathcal{LN}(\sigma|0.002,0.5)$
Table 6: Two-dimensional sparse learning setup where parameters $e_{5}$ and
$e_{6}$ are treated as questionable.

Figure 9 shows the RMS relevance indicator for parameter $e_{5}$ and $e_{6}$
during multistart Newton iterations. Parameter $e_{5}$ is rendered irrelevant
irrespective of the optimum as $\gamma^{\text{rms}}_{1}$ converges to zero for
both the optima. The parameter $e_{6}$ requires a special attention as
$\gamma^{\text{rms}}_{2}$ in Figure 9 converges to a value in the range
[0.75,0.90] for the case of the global optimum. As demonstrated previously in
the unidimensional setting, $\gamma_{i}^{\text{rms}}$ value converges to
exactly one for relevant parameters. Also, based on our experience, the
convergence of relevance indicator $\gamma_{i}^{\text{rms}}$ to one is a
necessary condition for a parameter to be deemed relevant. This condition
implies that the posterior pdf of a relevant parameter should be largely
dictated by the likelihood function and not the data-informed ARD prior. Based
on this principle, the parameter $e_{6}$ should be deemed irrelevant.
Alternatively, one could define a tolerance $\gamma^{\text{tol}}$ for
determining the relevance of such parameters. A $\gamma^{\text{tol}}$ value of
$0.5$ implies parameter $e_{6}$ as relevant, whereas a $\gamma^{\text{tol}}$
value of $0.9$ implies an irrelevant $e_{6}$. This choice of an appropriate
$\gamma^{\text{tol}}$ imparts greater flexibility to the modeller, thereby
allowing the sparse learning process to align with the modelling goals.

Figure 9: NSBL results for the relevance determination of parameter $e_{5}$
and $e_{6}$.

Figure 10 shows the marginal parameter pdf before and after the inclusion of
the optimal ARD prior for questionable parameters $e_{5}$ and $e_{6}$ at both
optima identified in Figure 9. The pre-NSBL pdf, shown in blue, is obtained
from the KDE representation of
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$
using 500 iid samples. The post-NSBL pdfs for the case where $e_{5}$ is
irrelavant, but $e_{6}$ are relevant is shown in red. The case where both
$e_{5}$ and $e_{6}$ are irrelevant are shown in green. These posterior pdfs
are obtained using the analytical solution derived in A.2. First, in red,
where only $e_{5}$ is correctly identified as irrelevant,the posterior pdf of
$e_{5}$ reduces to a Dirac-delta function centered at zero. This eliminates
all of the uncertainty in the estimate for $e_{5}$, and also can be shown to
reduce the uncertainty in the estimate for $e_{3}$, which is also a nonlinear
stiffness parameter. Second, in green, where both $e_{5}$ and $e_{6}$ are
correctly identified as irrelevant and are both reduced to Dirac-delta
function centered at zero, there is additional reduction in uncertainty
observed for both parameters $e_{4}$ and $e_{2}$, which are nonlinear and
linear damping coefficients, respectively. The removal of these higher-order
nonlinear terms, however, is shown here to have a minimal effect on the
posterior estimate for $\sigma$. This behaviour of the marginal posterior pdf
of $\sigma$ indicates that the removal of redundant parameters during NSBL
does not affect the data-fit property of the model. Had NSBL identified an
overly-simplistic model, the model error strength would be expected to
increase. Notice the decrease in posterior uncertainty of parameters $e_{3}$
and $e_{4}$ following the removal of $e_{5}$ and $e_{6}$. This decrease is
attributed to the remediation of overfitting by the sparse learning process.
Figure 10 shows the corresponding marginal posterior pdfs for the case when
$e_{5}$ is deemed irrelevant and $e_{6}$ is deemed relevant.

Figure 10: Marginal posterior pdf of model parameters before and after sparse
learning, for the case when $e_{5}$ and $e_{6}$ are deemed irrelevant.

The objective function $\mathcal{L}(\log\boldsymbol{\alpha})$ was observed to
be dictated by log-evidence for the range of $\boldsymbol{\alpha}$ values
shown in Figure 9, meaning both log-evidence and
$\mathcal{L}(\log\boldsymbol{\alpha})$ are bimodal. To our best knowledge,
this is the first reporting of a multimodal log-evidence for a physics-based
Bayesian inverse modelling. This multimodality also demonstrates the issue of
global identifiability in the model (or hyperparameter) space, while the
inverse problem may or may not be globally identifiable in the likelihood (or
model parameter) space. This multimodality also warrants careful consideration
of evidence-based Bayesian model comparison while dealing with nested models
with closely-related nonlinear terms.

NSBL was also executed for cases with 100, 500, and 1000 iid samples per KDE
instance of
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$.
These results are not reported here since they were identical to those
reported in this section but with even less variation in
$\log\boldsymbol{\alpha}$ iterates and relevance indicator
$\gamma^{\text{rms}}$ across multiple instances. This conclusion was true for
the three numerical cases studied in this section.

###  3.3 Four-dimensional sparse learning

Consider the four-dimensional sparse learning problem detailed in Table 7.
Here we only consider a single case, where among four questionable parameters,
two parameters are relevant, and two are irrelevant. NSBL is executed to
determine the relevance of questionable parameters $e_{3}$, $e_{4}$, $e_{5}$
and $e_{6}$. The NSBL setup is the same as Section 3.1. Figure 11 shows the
$\log\boldsymbol{\alpha}$ iterates and relevance indicator
$\gamma_{i}^{\text{rms}}$ for multistart Newton iterations beginning from
{-20, -20, -20, -20}, {-20, -20, -5, -5}, {-5,-5, -20, -20} and {-5, -5, -5,
-5}. Notice that all but $\log\alpha_{4}$ converges to a unique optimum. The
$\log\alpha_{4}$ converges to two different optimums depending on the starting
point. This type of behaviour has been previously reported in Figure 9. Based
on the $\mathcal{L}(\log\boldsymbol{\alpha})$ values at these optimums,
$\log\alpha^{\text{map}}_{4}$ = 0.0 was found to the global optimum. This
produces the solution $\log\boldsymbol{\alpha}^{\text{map}}$ = {-9.1, -12.4,
0.0, 0.0}, and the corresponding value of relevance indicator as {1.0, 1.0,
0.0, 0.0}. This optimal value of $\gamma_{i}^{\text{rms}}$ indicates that
parameter $e_{3}$ and $e_{4}$ are important, while parameter $e_{5}$ and
$e_{6}$ are redundant. This conclusion is in agreement with the data-
generating model in Eq. (16b) that was used to generate $\mathcal{D}$. Figure
12 shows the marginal posterior pdf before and after the inclusion of the
optimal ARD prior computed using NSBL. The pre-NSBL pdf is obtained using 500
iid samples from
$\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$,
while the post-NSBL pdf is obtained using the kernel-based analytical solution
derived in A.2. Once again, in Figure 12, we show the pre-NSBL results in
blue, the local optimum (where $e_{5}$ is relevant and $e_{6}$ is relevant) in
red, and the global optimum (where both $e_{5}$ and $e_{6}$ are irrelevant) in
green. Notice the decrease in uncertainty in relevant parameters following the
removal of parameter $e_{5}$ and $e_{6}$. Note that the level of sparsity in
physical models such as the LCO models is significantly lower than data-based
models due to the presence of physics-based parameters. Significant benefits
are realized in terms of decrease in uncertainty in the posterior pdf even
with the removal of two parameters.

Aerodynamic model | $\displaystyle\frac{\dot{C}_{\text{M}}}{B}+C_{\text{M}}=e_{1}\theta+e_{2}\dot{\theta}+e_{3}\theta^{3}+e_{4}\theta^{2}\dot{\theta}+e_{5}\theta^{5}+e_{6}\theta^{4}\dot{\theta}+\frac{c_{6}}{B}\ddot{\theta}+\sigma\xi(\tau)$
---|---
$\boldsymbol{\phi}$ decomposition |  $\boldsymbol{\phi}_{\alpha}$ = $\\{e_{3},e_{4},e_{5},e_{6}\\}$ , $\boldsymbol{\phi}_{\text{-}\alpha}$ = $\\{B$, $e_{1}$, $e_{2}$, $\sigma\\}$
ARD prior, $\mbox{p}(\boldsymbol{\phi}_{\alpha}|\boldsymbol{\alpha})$ |  $\mathcal{N}(e_{3}|0,\alpha_{1}^{-1})\mathcal{N}(e_{4}|0,\alpha_{2}^{-1})\mathcal{N}(e_{5}|0,\alpha_{3}^{-1})\mathcal{N}(e_{6}|0,\alpha_{4}^{-1})$,
Known prior, $\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})$ | $\mathcal{LN}(B|0.2,0.5)\,\mathcal{U}(e_{1}|{-2},0)\,\mathcal{U}(e_{2}|{-2},0)\mathcal{LN}(\sigma|0.002,0.5)$
Table 7: Four-dimensional sparse learning setup. Figure 11: NSBL results for
the four-dimensional sparse learning. Figure 12: Marginal posterior pdf of
model parameters before the after the inclusion of optimal ARD prior.

##  4 Conclusion

With the use of synthetic data, we have demonstrated that NSBL is a promising
approach to computationally efficient model selection for nonlinear stochastic
dynamical systems. In all cases studied herein, the algorithm has shown its
ability to recapture the data-generating model parameters and correctly
identify the relevance of the parameters. Note that for this case study, the
data is not particularly sparse temporally, nor is it particularly noisy,
though, the data only represents a partial observation of the state. For the
given data, even an overly complex model may not necessarily yield poor
predictions of the response, nonetheless, it has been shown that the parameter
estimates can be significantly improved by the removal of unnecessary
parameters.

The specific aspects of NSBL that make it particularly well suited to
engineering systems is the ability to use non-Gaussian known priors, and the
applicability to systems where there is a nonlinear mapping from the uncertain
parameters to the observed outputs. In the aeroelastic example, the underlying
physics and/or laws of statistics prevented certain parameters from assuming
non-positive values, hence a Gaussian distributions would not be proper
choices of prior. Likewise, the LCO being modelled is a nonlinear phenomenon,
which naturally leads to a nonlinear-in-parameter model in the machine
learning context. Both of these qualities are enabled by the use of a GMM
approximation of the product of the likelihood function and the known prior.
This GMM also enables the semi-analytical machinery that enables the efficient
computation of estimates of many quantities of interest including the
parameter posterior distribution, the model evidence, the objective function,
Jacobian/ gradient vector and Hessian matrix.

Finally, we have illustrated the need for robust optimization scheme as there
is the potential for multiple optima in the objective function, corresponding
to different locally optimal values of the hyperparameter, corresponding to
different classifications of the parameter as relevant versus irrelevant.

## References

  * [1] Rimple Sandhu.  Model comparison and sparse learning of nonlinear physics-based models using Bayesian inference.  PhD thesis, Carleton University, 2020.
  * [2] Rimple Sandhu, Mohammad Khalil, Chris Pettit, Dominique Poirel, and Abhijit Sarkar.  Nonlinear sparse bayesian learning for physics-based models.  Journal of Computational Physics, 426:109728, 2021.
  * [3] D. Poirel, Y. Harris, and A. Benaissa.  Self-sustained aeroelastic oscillations of a NACA0012 airfoil at low-to-moderate Reynolds numbers.  Journal of Fluids and Structures, 24(5):700 – 719, 2008.
  * [4] Daniel Rudmin, Azemi Benaissa, and Dominique Poirel.  Detection of laminar flow separation and transition on a naca-0012 airfoil using surface hot-films.  Journal of Fluids Engineering, 135(10), 2013.
  * [5] D. Poirel and W. Yuan.  Aerodynamics of laminar separation flutter at a transitional Reynolds number.  Journal of Fluids and Structures, 26(7-8):1174–1194, October 2010\.
  * [6] D Poirel, V Metivier, and G Dumas.  Computational aeroelastic simulations of self-sustained pitch oscillations of a naca0012 at transitional reynolds numbers.  Journal of Fluids and Structures, 27(8):1262–1277, 2011.
  * [7] Rimple Sandhu, Mohammad Khalil, Abhijit Sarkar, and Dominique Poirel.  Bayesian model selection for nonlinear aeroelastic systems using wind-tunnel data.  Computer Methods in Applied Mechanics and Engineering, 282:161 – 183, 2014.
  * [8] Rimple Sandhu, Dominique Poirel, Chris Pettit, Mohammad Khalil, and Abhijit Sarkar.  Bayesian inference of nonlinear unsteady aerodynamics from aeroelastic limit cycle oscillations.  Journal of Computational Physics, 316:534 – 557, 2016.
  * [9] Philippe Bisaillon, Rimple Sandhu, Chris Pettit, Mohammad Khalil, Dominique Poirel, CS Manohar, and Abhijit Sarkar.  Combined selection of the dynamic model and modeling error in nonlinear aeroelastic systems using bayesian inference.  Journal of Sound and Vibration, 522:116418, 2022.
  * [10] Rimple Sandhu, Chris Pettit, Mohammad Khalil, Dominique Poirel, and Abhijit Sarkar.  Bayesian model selection using automatic relevance determination for nonlinear dynamical systems.  Computer Methods in Applied Mechanics and Engineering, 320:237 – 260, 2017.
  * [11] Matthew Muto and James L Beck.  Bayesian updating and model class selection for hysteretic structural models using stochastic simulation.  Journal of Vibration and Control, 14(1-2):7–34, 2008.
  * [12] Michael E. Tipping.  Sparse bayesian learning and the relevance vector machine.  J. Mach. Learn. Res., 1:211–244, sep 2001.
  * [13] Anita C. Faul and Michael E. Tipping.  Analysis of sparse bayesian learning.  In Advances in Neural Information Processing Systems 14, pages 383–389. MIT Press, 2001.
  * [14] S. D. Babacan, R. Molina, and A. K. Katsaggelos.  Bayesian compressive sensing using laplace priors.  IEEE Transactions on Image Processing, 19(1):53–63, Jan 2010.
  * [15] Michael E. Tipping.  Sparse bayesian learning and the relevance vector machine.  Journal of Machine Learning Research, 1:211–244, September 2001\.
  * [16] Athanasios Papoulis.  Probability, random variables, and stochastic processes.  McGraw Hill, fourth edition, 2002.
  * [17] Christopher M. Bishop.  Pattern Recognition and Machine Learning (Information Science and Statistics).  Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2006.
  * [18] Kevin P. Murphy.  Probabilistic Machine Learning: An introduction.  MIT Press, 2022.
  * [19] Jorge Nocedal and Stephen J. Wright.  Numerical Optimization.  Springer, New York, NY, USA, second edition, 2006.
  * [20] Philippe Bisaillon, Rimple Sandhu, Mohammad Khalil, Chris Pettit, Dominique Poirel, and Abhijit Sarkar.  Bayesian parameter estimation and model selection for strongly nonlinear dynamical systems.  Nonlinear Dynamics, 82(3):1061–1080, 2015.
  * [21] Heikki Haario, Marko Laine, Antonietta Mira, and Eero Saksman.  Dram: efficient adaptive mcmc.  Statistics and computing, 16(4):339–354, 2006.
  * [22] Jianye Ching and Yi-Chu Chen.  Transitional markov chain monte carlo method for bayesian model updating, model class selection, and model averaging.  Journal of engineering mechanics, (7):816–832, 2007.
  * [23] D. W. Scott.  Multivariate Density Estimation: Theory, Practice, and Visualization.  John Wiley & Sons, New York, Chicester, 1992.

##  Appendix A Semi-analytical calculation of Bayesian entities

###  A.1 Model evidence

The model evidence was first introduced as the normalization factor in Eq.
(1). Substituting the expression for the joint prior pdf in Eq. (2), the model
evidence can be rewritten as [2]

$\mbox{p}(\mathcal{D}|\boldsymbol{\alpha})=\int\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}|\boldsymbol{\alpha})d\boldsymbol{\phi}=\int\mbox{p}(\mathcal{D}|\boldsymbol{\phi})\mbox{p}(\boldsymbol{\phi}_{\text{-}\alpha})\mbox{p}(\boldsymbol{\phi}_{\alpha}|\boldsymbol{\alpha})d\boldsymbol{\phi}.$
(17)

Now, substituting the GMM approximation of the product of the likelihood
function and the known prior from Eq. (6), and the expression for the ARD
prior, $\mbox{p}(\boldsymbol{\phi}|\boldsymbol{\alpha})$ =
$\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf 0},\mbox{\bf A}^{-1})$, Eq.
(17) becomes

$\displaystyle\hat{\mbox{p}}(\mathcal{D}|\boldsymbol{\alpha})$
$\displaystyle=\int\left\\{\sum_{k=1}^{K}a^{(k)}\mathcal{N}(\boldsymbol{\phi}|\boldsymbol{\mu}^{(k)},\boldsymbol{\Sigma}^{(k)})\right\\}\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
0},\mbox{\bf A}^{-1})d\boldsymbol{\phi}$
$\displaystyle=\sum_{k=1}^{K}a^{(k)}\int\mathcal{N}(\boldsymbol{\phi}|\boldsymbol{\mu}^{(k)},\boldsymbol{\Sigma}^{(k)})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
0},\mbox{\bf A}^{-1})d\boldsymbol{\phi}.$ (18)

The two Gaussian distributions in the integral expression will generally be
mismatched in dimensions ($N_{\boldsymbol{\phi}}$ and
$N_{\boldsymbol{\phi}_{\alpha}}$). In order to facilitate the integration, we
rewrite the $k$th kernel of the GMM in an expanded form as

$\mathcal{N}(\boldsymbol{\phi}|\boldsymbol{\mu}^{(k)},\boldsymbol{\Sigma}^{(k)})=\mathcal{N}\left(\begin{Bmatrix}\boldsymbol{\phi}_{\alpha}\\\
\boldsymbol{\phi}_{\text{-}\alpha}\end{Bmatrix}\Bigg{|}\begin{Bmatrix}\boldsymbol{\mu}_{\alpha}^{(k)}\\\
\boldsymbol{\mu}_{\text{-}\alpha}^{(k)}\end{Bmatrix},\begin{bmatrix}\boldsymbol{\Sigma}_{\alpha}^{(k)}&\mbox{\bf
C}^{(k)}\\\ (\mbox{\bf
C}^{(k)})^{T}&\boldsymbol{\Sigma}_{\text{-}\alpha}^{(k)}\end{bmatrix}\right),$
(19)

where $\boldsymbol{\mu}_{\alpha}^{(k)}$ and
$\boldsymbol{\Sigma}_{\alpha}^{(k)}$ are the mean and covariance of the $k$th
kernel pertaining to $\boldsymbol{\phi}_{\alpha}$, likewise
$\boldsymbol{\mu}_{\text{-}\alpha}^{(k)}$ and
$\boldsymbol{\Sigma}_{\text{-}\alpha}^{(k)}$ are the mean and covariance of
the $k$th kernel pertaining to $\boldsymbol{\phi}_{\text{-}\alpha}$, and
$\mbox{\bf C}^{(k)}$ is the cross-covariance of $\boldsymbol{\phi}_{\alpha}$
and $\boldsymbol{\phi}_{\text{-}\alpha}$. These entities are obtained from the
construction of the GMM (whether by KDE, EM, or other means). Now, rewriting
this joint distribution of $\boldsymbol{\phi}_{\alpha}$ and
$\boldsymbol{\phi}_{\text{-}\alpha}$ as the product of the conditional
distribution of $\boldsymbol{\phi}_{\text{-}\alpha}$ given
$\boldsymbol{\phi}_{\alpha}$ and the marginal distribution of
$\boldsymbol{\phi}_{\alpha}$, we obtain [2]

$\displaystyle\mathcal{N}(\boldsymbol{\phi}|\boldsymbol{\mu}^{(k)},\boldsymbol{\Sigma}^{(k)})$
$\displaystyle=\mathcal{N}(\boldsymbol{\phi}_{\text{-}\alpha}|\tilde{\boldsymbol{\mu}}_{\text{-}\alpha}^{(k)},\tilde{\boldsymbol{\Sigma}}_{\text{-}\alpha}^{(k)})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\boldsymbol{\mu}_{\alpha}^{(k)},\boldsymbol{\Sigma}_{\alpha}^{(k)}),$
(20a) $\displaystyle\tilde{\boldsymbol{\mu}}_{\text{-}\alpha}^{(k)}$
$\displaystyle=\boldsymbol{\mu}_{\text{-}\alpha}^{(k)}+(\mbox{\bf
C}^{(k)})^{T}(\boldsymbol{\Sigma}_{\alpha}^{(k)})^{-1}(\boldsymbol{\phi}_{\alpha}-\boldsymbol{\mu}_{\alpha}^{(k)}),$
(20b) $\displaystyle\tilde{\boldsymbol{\Sigma}}_{\text{-}\alpha}^{(k)}$
$\displaystyle=\boldsymbol{\Sigma}_{\text{-}\alpha}^{(k)}-(\mbox{\bf
C}^{(k)})^{T}(\boldsymbol{\Sigma}_{\alpha}^{(k)})^{-1}\mbox{\bf C}^{(k)},$
(20c)

where $\tilde{\boldsymbol{\mu}}_{\text{-}\alpha}^{(k)}$ and
$\tilde{\boldsymbol{\Sigma}}_{\text{-}\alpha}^{(k)}$ are the mean and
covariance of the $k$th kernel, pertaining to
$\boldsymbol{\phi}_{\text{-}\alpha}$ conditioned on a given value of
$\boldsymbol{\phi}_{\alpha}$. Consequently, the estimate of the model evidence
becomes [2]

$\hat{\mbox{p}}(\mathcal{D}|\boldsymbol{\alpha})=\sum_{k=1}^{K}a^{(k)}\int\mathcal{N}(\boldsymbol{\phi}_{\text{-}\alpha}|\tilde{\boldsymbol{\mu}}_{\text{-}\alpha}^{(k)},\tilde{\boldsymbol{\Sigma}}_{\text{-}\alpha}^{(k)})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\boldsymbol{\mu}_{\alpha}^{(k)},\boldsymbol{\Sigma}_{\alpha}^{(k)})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
0},\mbox{\bf A}^{-1})d\boldsymbol{\phi},$ (21)

where the product of the two left-most Gaussians can be rewritten as

$\displaystyle\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\boldsymbol{\mu}_{\alpha}^{(k)},\boldsymbol{\Sigma}_{\alpha}^{(k)})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
0},\mbox{\bf A}^{-1})$
$\displaystyle=\mathcal{N}(\boldsymbol{\mu}_{\alpha}^{(k)}|\mbox{\bf
0},\mbox{\bf
B}_{\alpha}^{(k)})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
m}_{\alpha}^{(k)}\mbox{\bf P}_{\alpha}^{(k)})$ (22a) $\displaystyle\mbox{\bf
m}_{\alpha}^{(k)}$
$\displaystyle=\boldsymbol{\mu}_{\alpha}^{(k)}-\boldsymbol{\Sigma}_{\alpha}^{(k)}(\mbox{\bf
B}_{\alpha}^{(k)})^{-1}\boldsymbol{\mu}_{\alpha}^{(k)}$ (22b)
$\displaystyle\mbox{\bf P}_{\alpha}^{(k)}$
$\displaystyle=\boldsymbol{\Sigma}_{\alpha}^{(k)}-\boldsymbol{\Sigma}_{\alpha}^{(k)}(\mbox{\bf
B}_{\alpha}^{(k)})^{-1}\boldsymbol{\Sigma}_{\alpha}^{(k)}$ (22c)
$\displaystyle\mbox{\bf B}_{\alpha}^{(k)}$
$\displaystyle=\boldsymbol{\Sigma}_{\alpha}^{(k)}+\mbox{\bf A}^{-1}$ (22d)

where $\mbox{\bf m}_{\alpha}^{(k)}$ and $\mbox{\bf P}_{\alpha}^{(k)}$ are the
mean and covariance of the $k$th kernel of the posterior of
$\boldsymbol{\phi}_{\alpha}$. Note that
$\mathcal{N}(\boldsymbol{\mu}_{\alpha}^{(k)}|\mbox{\bf 0},\mbox{\bf
B}_{\alpha}^{(k)})$ is now independent of the parameters. Substituting and
integrating over all parameters gives [2]

$\displaystyle\hat{\mbox{p}}(\mathcal{D}|\boldsymbol{\alpha})$
$\displaystyle=\sum_{k=1}^{K}a^{(k)}\int\mathcal{N}(\boldsymbol{\phi}_{\text{-}\alpha}|\tilde{\boldsymbol{\mu}}_{\text{-}\alpha}^{(k)},\tilde{\boldsymbol{\Sigma}}_{\text{-}\alpha}^{(k)})\mathcal{N}(\boldsymbol{\mu}_{\alpha}^{(k)}|\mbox{\bf
0},\mbox{\bf
B}_{\alpha}^{(k)})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
m}_{\alpha}^{(k)}\mbox{\bf P}_{\alpha}^{(k)})d\boldsymbol{\phi}$
$\displaystyle=\sum_{k=1}^{K}a^{(k)}\mathcal{N}(\boldsymbol{\mu}_{\alpha}^{(k)}|\mbox{\bf
0},\mbox{\bf B}_{\alpha}^{(k)}).$ (23)

###  A.2 Posterior parameter pdf

As the estimate of the model evidence is now available in Eq. (21), it is now
possible to obtain an estimate of the parameter posterior pdf
$\hat{\mbox{p}}(\boldsymbol{\phi}|\mathcal{D},\boldsymbol{\alpha})$ making the
same substitutions from above to obtain [2]

$\hat{\mbox{p}}(\boldsymbol{\phi}|\mathcal{D},\boldsymbol{\alpha})=\sum_{k=1}^{K}\cfrac{a^{(k)}\mathcal{N}(\boldsymbol{\phi}|\boldsymbol{\mu}^{(k)},\boldsymbol{\Sigma}^{(k)})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
0},\mbox{\bf A}^{-1})}{\hat{\mbox{p}}(\mathcal{D}|\boldsymbol{\alpha})}.$ (24)

Substituting
$\mathcal{N}(\boldsymbol{\phi}|\boldsymbol{\mu}^{(k)},\boldsymbol{\Sigma}^{(k)})$
from Eq. (20a) reduces Eq. (24) to [2]

$\displaystyle\hat{\mbox{p}}(\boldsymbol{\phi}|\mathcal{D},\boldsymbol{\alpha})$
$\displaystyle=\sum_{k=1}^{K}\cfrac{a^{(k)}\mathcal{N}(\boldsymbol{\phi}_{\text{-}\alpha}|\tilde{\boldsymbol{\mu}}_{\text{-}\alpha}^{(k)},\tilde{\boldsymbol{\Sigma}}_{\text{-}\alpha}^{(k)})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\boldsymbol{\mu}_{\alpha}^{(k)},\boldsymbol{\Sigma}_{\alpha}^{(k)})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
0},\mbox{\bf A}^{-1})}{\hat{\mbox{p}}(\mathcal{D}|\boldsymbol{\alpha})}$
$\displaystyle=\sum_{k=1}^{K}\underbrace{\left(\cfrac{a^{(k)}\mathcal{N}(\boldsymbol{\mu}_{\alpha}^{(k)}|\mbox{\bf
0},\mbox{\bf
B}_{\alpha}^{(k)})}{\sum_{r=1}^{K}a^{(r)}\mathcal{N}(\boldsymbol{\mu}_{\alpha}^{(r)}|\mbox{\bf
0},\mbox{\bf
B}_{\alpha}^{(r)})}\right)}_{w^{(k)}}\mathcal{N}(\boldsymbol{\phi}_{\text{-}\alpha}|\tilde{\boldsymbol{\mu}}_{\text{-}\alpha}^{(k)},\tilde{\boldsymbol{\Sigma}}_{\text{-}\alpha}^{(k)})\mathcal{N}(\boldsymbol{\phi}_{\alpha}|\mbox{\bf
m}_{\alpha}^{(k)},\mbox{\bf P}_{\alpha}^{(k)})$ (25)
$\displaystyle=\sum_{k=1}^{K}w^{(k)}\mathcal{N}\left(\begin{Bmatrix}\boldsymbol{\phi}_{\alpha}\\\
\boldsymbol{\phi}_{\text{-}\alpha}\end{Bmatrix}\Bigg{|}\begin{Bmatrix}\mbox{\bf
m}_{\alpha}^{(k)}\\\ \mbox{\bf
m}_{\text{-}\alpha}^{(k)}\end{Bmatrix},\begin{bmatrix}\mbox{\bf
P}_{\alpha}^{(k)}&{\mbox{\bf D}}^{(k)}\\\ ({\mbox{\bf D}}^{(k)})^{T}&\mbox{\bf
P}_{\text{-}\alpha}^{(k)}\end{bmatrix}\right)$ (26)

where $0\leq w^{(k)}\leq 1$, $\sum_{k}w^{(k)}=1$ is the weight coefficient of
kernel $k$, and $\mbox{\bf m}^{(k)}$ and $\mbox{\bf P}^{(k)}$ are the mean and
covariance of the $k$th kernel of the the posterior of $\boldsymbol{\phi}$,
respectively. Note that $\mbox{\bf m}_{\alpha}^{(k)}$ and $\mbox{\bf
P}_{\alpha}^{(k)}$ are known from Eq. (22b), and Eq. (22c). Furthermore,
$\mbox{\bf m}_{\text{-}\alpha}^{(k)}$ and $\mbox{\bf
P}_{\text{-}\alpha}^{(k)}$ are the mean and covariance of the $k$th kernel of
the posterior of $\boldsymbol{\phi}_{\text{-}\alpha}$ and ${\mbox{\bf
D}}^{(k)}$ is the cross-covariance of $\boldsymbol{\phi}_{\alpha}$ and
$\boldsymbol{\phi}_{\text{-}\alpha}$. These entities are given by

$\displaystyle\mbox{\bf m}_{\text{-}\alpha}^{(k)}$
$\displaystyle=\boldsymbol{\mu}_{\text{-}\alpha}^{(k)}-(\mbox{\bf
C}^{(k)})^{T}(\mbox{\bf
B}_{\alpha}^{(k)})^{-1}\boldsymbol{\mu}_{\alpha}^{(k)}$ (27a)
$\displaystyle\mbox{\bf P}_{\text{-}\alpha}^{(k)}$
$\displaystyle=\boldsymbol{\Sigma}_{\text{-}\alpha}^{(k)}-(\mbox{\bf
C}^{(k)})^{T}(\mbox{\bf B}_{\alpha}^{(k)})^{-1}\mbox{\bf C}^{(k)}$ (27b)
$\displaystyle{\mbox{\bf D}}^{(k)}$ $\displaystyle=\mbox{\bf
C}^{(k)}-\boldsymbol{\Sigma}_{\alpha}^{(k)}(\mbox{\bf
B}_{\alpha}^{(k)})^{-1}\boldsymbol{\mu}_{\alpha}^{(k)}((\mbox{\bf
C}^{(k)})^{T}(\mbox{\bf
B}_{\alpha}^{(k)})^{-1}\boldsymbol{\mu}_{\alpha}^{(k)})^{T}$ (27c)

###  A.3 Gradient vector and Hessian matrix

Let $J_{i}(\log\boldsymbol{\alpha})$ denote the $i$th element of the Jacobian
vector $\mbox{\bf J}(\log\boldsymbol{\alpha})$. Differentiating the objective
function in Eq. (10) with respect to the $\log(\alpha_{i})$ gives [2]

$\displaystyle J_{i}(\log\boldsymbol{\alpha})$
$\displaystyle=\frac{\partial\mathcal{L}(\log\boldsymbol{\alpha})}{\partial\log\alpha_{i}}=\frac{\partial}{\partial\log\alpha_{i}}\left\\{\log\hat{\mbox{p}}(\mathcal{D}|\log\boldsymbol{\alpha})+\sum_{i=1}^{N_{\alpha}}\left(r_{i}\log\alpha_{i}-s_{i}\alpha_{i}\right)\right\\}$
$\displaystyle=\frac{\partial\log\hat{\mbox{p}}(\mathcal{D}|\log\boldsymbol{\alpha})}{\partial\log\alpha_{i}}+r_{i}-s_{i}\alpha_{i}\;$
$\displaystyle=\frac{1}{\hat{\mbox{p}}(\mathcal{D}|\log\boldsymbol{\alpha})}\frac{\partial}{\partial\log\alpha_{i}}\left\\{\sum_{k=1}^{K}a^{(k)}\mathcal{N}(\boldsymbol{\mu}_{\alpha}^{(k)}|\mbox{\bf
0},\mbox{\bf B}_{\alpha}^{(k)})\right\\}+r_{i}-s_{i}\alpha_{i}$
$\displaystyle=\frac{1}{\hat{\mbox{p}}(\mathcal{D}|\log\boldsymbol{\alpha})}\sum_{k=1}^{K}a^{(k)}\frac{\partial}{\partial\log\alpha_{i}}\left\\{\exp\left(\log\mathcal{N}(\boldsymbol{\mu}_{\alpha}^{(k)}|\mbox{\bf
0},\mbox{\bf B}_{\alpha}^{(k)})\right)\right\\}+r_{i}-s_{i}\alpha_{i}$
$\displaystyle=\sum_{k=1}^{K}\underbrace{\left(\frac{a^{(k)}\mathcal{N}(\boldsymbol{\mu}_{\alpha}^{(k)}|\mbox{\bf
0},\mbox{\bf
B}_{\alpha}^{(k)})}{\hat{\mbox{p}}(\mathcal{D}|\log\boldsymbol{\alpha})}\right)}_{w^{(k)}}\underbrace{\left(\frac{\partial\log\mathcal{N}(\boldsymbol{\mu}_{\alpha}^{(k)}|\mbox{\bf
0},\mbox{\bf
B}_{\alpha}^{(k)})}{\partial\log\alpha_{i}}\right)}_{v_{i}^{(k)}}+r_{i}-s_{i}\alpha_{i}$
(28) $\displaystyle=\sum_{k=1}^{K}w^{(k)}v_{i}^{(k)}+r_{i}-s_{i}\alpha_{i},$
(29)

where $w^{(k)}$ is known from Eq. (25), and factor $v_{i}^{(k)}$ is given as

$\displaystyle v_{i}^{(k)}$
$\displaystyle=-\frac{1}{2}\left\\{-1+\alpha_{i}P_{ii}^{(k)}+\alpha_{i}(m_{i}^{(k)})^{2}\right\\}.$
(30)

Now, let $H_{ij}(\log\boldsymbol{\alpha})$ denote element $(i,j)$ of the
Hessian matrix $\mbox{\bf H}(\log\boldsymbol{\alpha})$. Differentiating the
$j$th element of the Jacobian in Eq. (29) with respect to the $i$th element
$\boldsymbol{\alpha}$ gives [2]

$\displaystyle H_{ij}(\log\boldsymbol{\alpha})$
$\displaystyle=\frac{\partial^{2}\mathcal{L}(\log\boldsymbol{\alpha})}{\partial\log\alpha_{i}\partial\log\alpha_{j}}=\frac{\partial
J_{j}(\log\boldsymbol{\alpha})}{\partial\log\alpha_{i}}=\frac{\partial}{\partial\log\alpha_{i}}\left\\{\sum_{k=1}^{K}w^{(k)}v_{j}^{(k)}+r_{j}-s_{j}\alpha_{j}\right\\}$
$\displaystyle=\sum_{k=1}^{K}\left\\{w^{(k)}\frac{\partial
v_{j}^{(k)}}{\partial\log\alpha_{i}}+v_{j}^{(k)}\frac{\partial
w^{(k)}}{\partial\log\alpha_{i}}\right\\}-\delta_{ij}s_{i}\alpha_{i},$ (31)

with

$\displaystyle\frac{\partial v_{j}^{(k)}}{\partial\log\alpha_{i}}$
$\displaystyle=\alpha_{i}\alpha_{j}\left(\frac{(P_{ij}^{(k)})^{2}}{2}+m_{i}^{(k)}m_{j}^{(k)}P_{ij}^{(k)}\right)+\delta_{ij}\left(v_{i}^{(k)}-\frac{1}{2}\right),$
(32) $\displaystyle\frac{\partial w^{(k)}}{\partial\log\alpha_{i}}$
$\displaystyle=w^{(k)}\left(v_{i}^{(k)}-\bar{v}_{i}\right),$ (33)

where $\delta_{ij}$ is the Kronecker delta function, and $\bar{v}_{i}$ is the
average of $v_{i}^{(k)}$ over all $K$ kernels.