11institutetext: Cyprus University of Technology, Limassol Cyprus
11email<EMAIL_ADDRESS>
22institutetext: National and Kapodistrian University of Athens, Athens,
Greece33institutetext: Aalborg University, Denmark

# Variational Conditional Dependence Hidden Markov Models for Skeleton-Based
Action Recognition††thanks: This work has received funding from the European
Union’s Horizon 2020 research and innovation program under grant agreement No
872139, project aiD.

Konstantinos P. Panousis 11 Soritios Chatzis 11 Sergios Theodoridis 2233

###### Abstract

Hidden Markov Models (HMMs) comprise a powerful generative approach for
modeling sequential data and time-series in general. However, the commonly
employed assumption of the dependence of the current time frame to a single or
multiple immediately preceding frames is unrealistic; more complicated
dynamics potentially exist in real world scenarios. This paper revisits
conventional sequential modeling approaches, aiming to address the problem of
capturing time-varying temporal dependency patterns. To this end, we propose a
different formulation of HMMs, whereby the dependence on past frames is
dynamically inferred from the data. Specifically, we introduce a hierarchical
extension by postulating an additional latent variable layer; therein, the
(time-varying) temporal dependence patterns are treated as latent variables
over which inference is performed. We leverage solid arguments from the
Variational Bayes framework and derive a tractable inference algorithm based
on the forward-backward algorithm. As we experimentally show, our approach can
model highly complex sequential data and can effectively handle data with
missing values.

###### Keywords:

Hidden Markov Models Approximate Inference Temporal Dependence.

##  1 Introduction

Modeling sequential data, typically encountered in many real-world
applications such as bioinformatics and computer vision, remains a fundamental
task in the field of machine learning. Hidden Markov Models have been a
popular approach for modeling such data, but more recently, have largely been
replaced by their “deep” variants. Despite their success, these approaches
exhibit significant drawbacks such as over-parameterization even in simple
tasks as image classification[20]. Nevertheless, both models aim to learn
hidden representations for sequential data.

Commonly, first-order HMMs are usually considered; this allows for simplicity
and low computational complexity. Then, model temporal dynamics are restrained
to a simple step back, undermining model effectiveness. This compromise is
especially detrimental in many applications where longer temporal dynamics may
be present. Researchers have considered alternatives to alleviate this
restriction, for example by introducing extended temporal dependencies in the
form of second or higher order Markovian dynamics and had successful
applications in many domains [2, 17, 13]. However, the increased order
introduces three significant drawbacks: (i) additional complexity to the model
that may render it unusable in real-world scenarios; (ii) unnecessary burden
to the researchers, namely the need to determine the best postulated order for
each specific application and dataset; and (iii) unnecessary complex models
with strong overfitting tendencies, which may negate the benefits and
flexibility of HMMs.

Lastly, a further drawback of conventional HMM approaches is the static and
homogeneity assumptions [3, 10, 8], where the possible dynamic temporal
dependencies of the data are ignored. The same effect applies to hidden semi-
Markov models (HSMMs) [26]; even though the relaxed temporal assumptions allow
for more flexible modeling of the dynamics, potential non-homogeneous temporal
dynamics in the data are still ignored. Drawing inspiration from there
results, researchers have developed HMMs with variable order Markov chains,
e.g., [5]. The resulting models have been shown to be effective in a diverse
range of applications, nevertheless exhibiting significant drawbacks such as
the inability to model continuous observations [12].

Recently, [29] introduced a variant of a simple HSMM model, dubbed
Hierarchical Dynamic Model (HDM). Specifically, the proposed model constitutes
a hierarchical extension by leveraging on the Bayesian framework to increase
the capacity of HSMMs to account for temporal and spatial variations in the
Human Action Recognition task. The produced generative model is shown to be
more robust to the natural variations of the data; it exhibits increased
generalization capabilities, while at the same time requiring less data to
train.

This paper draws from these results and attempts to offer a principled way of
modeling sequential data with complex dynamics. To this end, we develop a
novel variant of HMMs that is able to capture complex temporal dependencies.
The considered approach constitutes a hierarchical model: by postulating an
additional latent first-order Markov Chain, called the dependence-generator
layer, the model can alter the effective temporal dynamics of the conventional
observation-emitting Markov Chain. In this way, the proposed approach can
effectively infer which past state more strongly affects the current time
frame. The proposed inferential construction is enhanced by a fully Bayesian
treatment under the Variational Inference framework. We posit that the
proposed hierarchical variant, combined with the flexibility and effectiveness
of Variational Bayes, can greatly increase the capacity of the resulting
architecture to model sequential data that exhibit a complex combination of
spatial and temporal variations.

We formulate efficient training and inference algorithms for our approach by:
(i) deriving a variant of the classical forward-backward algorithm used in
HMMs, (ii) relying on Variational Inference and conjugate priors for closed-
form solutions. We dub our approach Variational Bayesian Conditional
Dependence Hidden Markov Model (VB-CD-HMM) and evaluate its performance in the
Human Action Recognition domain. Human Action Recognition poses one of the
most challenging tasks in the computer vision community. Researchers have
devoted significant effort to address this particular area [25] and remarkable
progress has been made with the popularization of action recognition through
3D data [1]. In this context, various approaches have been employed in order
to model the activity recognition task, such as HMM-based approaches [29],
LSTM-based [28] and Gaussian Processes [6, 7]. However, the significant
spatial and temporal variations arising in the execution of an action, along
with variations due to the capturing process, e.g., camera occlusion, leave
room for needed improvements in existing approaches. We provided strong
empirical evidence for our approach though a series of experiments dealing
with Human Action Recognition benchmarks. Our evidence shows competitive
predictive accuracy, persistent even in the presence of missing values.

##  2 Proposed Model

In this work, we consider an extension of a conventional first-order HMM,
whereby a hierarchical structure is employed consisting of two layers: (i) a
dependence-generator layer; this comprises a simple first-order Markov Chain
that determines the steps-back taken into account in the second layer; and
(ii) a chain of observation-emitting hidden states, where the temporal
dependencies are determined from the output of the dependence-generator layer.
The temporal dependencies are inferred in a data-driven fashion. An outline of
the proposed Conditional Dependence HMM (CD-HMM) model is provided in Fig.
1(a).

$z_{t-1}$$z_{t}$$z_{t+1}$$x_{t-1}$$x_{t}$$x_{t+1}$$\dots$$\dots$$x_{\scaleto{t-K+1}{3.5pt}}$$x_{t-K}$$y_{t-1}$$y_{t}$$y_{t+1}$$y_{\scaleto{t-K+1}{3.5pt}}$$y_{t-K}$$\dots$

(a) A graphical illustration of the Conditional Dependence Hidden Markov
Model.

$z_{t-1}$$z_{t}$$z_{t+1}$$\dots$$z_{\scaleto{t-K+1}{3.5pt}}$$z_{t-K}$$\dots$$z_{1}$$\hat{\pi}$$\hat{A}$$\alpha_{0}$$\alpha$$x_{t-1}$$x_{t}$$x_{t+1}$$\dots$$\dots$$x_{\scaleto{t-K+1}{3.5pt}}$$x_{t-K}$$\dots$$x_{1}$$\pi$$A$$\eta_{0}$$\eta$$y_{t-1}$$y_{t}$$y_{t+1}$$\dots$$y_{\scaleto{t-K+1}{3.5pt}}$$y_{t-K}$$\dots$$y_{1}$$\scaleto{\phi_{im}}{7.5pt}$$\lambda,m,\eta,S$

(b) The resulting hierarchical VB-CD-HMM.

###  2.1 Model Definition

Let $Y=\\{\boldsymbol{y}_{t}\\}_{t=1}^{T}\in\mathbb{R}^{D}$ be an observed
data sequence with $T$ time frames and $D$ features. Following the definition
of conventional HMMs, let us assume an $N$-state emitting state sequence
denoted as $X=\\{x_{t}\\}_{t=1}^{T},x_{t}\in[1,\dots,N]$ where $x_{t}$
indicates the state from which the $t^{th}$ observation is emitted. Each
emission density is modeled by an $M$-component finite mixture model; the
mixture component indicators are denoted as
$L=\\{l_{t}\\}_{t=1}^{T},l_{t}\in[1,\dots,M]$ with $l_{t}$ indicating which
component generated the $t^{th}$ observation.

The hierarchical construction postulates an additional layer, the dependence-
generator layer, such that the latent data associated with each sequence are
augmented with the supplementary sequence of temporal dependence indicators,
$Z=\\{z_{t}\\}_{t=1}^{T},z_{t}\in[1,\dots,K]$. These denote the current
temporal dependencies between the observation-emitting states at time $t$ and
times $t-1,\dots,t-K$ as shown in Fig. 1(a). Thus, different pairwise state
dependencies $(x_{t},x_{t-z_{t}})$ are valid at each time point $t$, dictated
by the inferred latent variables $z_{t}$.

The model parameters comprise $\theta=\\{\phi,\psi\\}$, where $\phi$ contains
the parameters of the emission distributions of the model and $\psi$ the
effective parameters of the latent processes.

For the emission distributions, finite mixture models of multivariate Gaussian
distributions are considered,
$p(\boldsymbol{y}_{t}|x_{t}=i)=\sum_{m=1}^{M}c_{im}\mathcal{N}(\boldsymbol{y}_{t}|\boldsymbol{\mu}_{im},R_{im}),\forall
i$ where $\mathcal{N}(\cdot|\boldsymbol{\mu},R)$ is a multivariate Normal
distribution with mean $\boldsymbol{\mu}\in\mathbb{R}^{D}$ and precision
matrix $R\in\mathbb{R}^{D\times D}$ with distinct parameters for each state
and mixture; $\boldsymbol{c}_{i}\in\mathbb{R}^{M}$ are the mixture component
coefficients for state $x_{t}=i$. Hence, the effective parameters for the
model distributions are
$\phi=\\{\boldsymbol{\mu}_{im},R_{im},c_{im}\\}_{i,m=1}^{N,M}$.

We proceed with the definition of the parameters and their relationships
pertaining to the postulated latent processes. Specifically, the first,
dependence-generator, layer is a first-order Markov Chain and is described by
the following initial and transition probabilities:

$\displaystyle\hat{\pi}_{k}$ $\displaystyle\triangleq p(z_{1}=k),\quad\forall
k$ (1) $\displaystyle\hat{A}_{kk^{\prime}}$ $\displaystyle\triangleq
p(z_{t}=k^{\prime}|z_{t-1}=k),\quad\forall t>1,k,k^{\prime}$ (2)

For the second, observation-emitting, layer process, the initial probability
of the emitting state reads:

$\displaystyle\pi_{i}\triangleq p(x_{1}=i),\quad\forall i$ (3)

To model the dependency between the first and second layers, we define the set
of conditional dependence transition probability matrices
$\\{A^{k}\\}_{k=1}^{K}$ on the second-layer process, where
$A^{k}\triangleq[A_{ij}^{k}]_{i,j=1}^{N}$,

$\displaystyle\begin{split}A_{ij}^{k}&\triangleq
p(x_{t}=j|x_{t-k}=i,z_{t}=k)\end{split}$ (4)

Thus, based on the inferred dependence of the dependence-generator layer, we
consider different temporal dependencies in the second layer at each time
point; these are regulated through the $z_{t}$ variables, which are imposed a
first-order Markov chain themselves. Intuitively, the conditional dependence
mechanism determines the strength of the dependence of the current state to
individual past states. The strongest inferred dependence is retained from the
model, and is encoded in the $z_{t}$ values; thus, the model dynamically
alters the transition dynamics for the observation-emitting state at each time
point.

Having fully defined the parameters for the latent processes, the
corresponding set comprises
$\psi=\\{\hat{\boldsymbol{\pi}},\hat{A},\boldsymbol{\pi},\\{A^{k}\\}_{k=1}^{K}\\}$.
The definition of the model is now concluded.

The joint distribution reads:

$\displaystyle\begin{split}&p(Y,X,Z)=\hat{\pi}_{z_{1}}\pi_{x_{1}}\prod_{t=1}^{T-1}\hat{A}_{z_{t},z_{t+1}}\prod_{t>1}A_{x_{t-z_{t}},x_{t}}^{z_{t}}\prod_{t=1}^{T}p(\boldsymbol{y}_{t}|x_{t})\end{split}$

###  2.2 Model Training

We train the CD-HMM model presented in the previous section considering a
variational Bayesian treatment. The Variational Bayes (VB) approach was chosen
over MCMC inference methods, e.g., Gibbs Sampling, considering the negligible
performance differences and the significantly lower computational complexity.
Sampling-based approaches encompass large computational complexity; this
restricts the applicability of the model to either more complex applications
and datasets or different settings (e.g., streaming).

The VB treatment of the model comprises the introduction of appropriate prior
distributions over all model parameters and maximization of the resulting
Evidence Lower Bound (ELBO) expression. To obtain closed-form solutions for
the updates of the parameters, we employ conjugate priors; a choice
accompanied by the benefits of lower computational complexity and
interpretability [23].
Priors and Evidence Lower Bound. Dirichlet distributions are imposed as the
priors of the initial state probabilities for both layer processes
$p(\boldsymbol{\pi})=\mathcal{D}(\boldsymbol{\pi}|\boldsymbol{\eta}_{0})$ and
$p(\hat{\boldsymbol{\pi}})=\mathcal{D}(\hat{\boldsymbol{\pi}}|\boldsymbol{\alpha}_{0})$.
Analogous Dirichlet priors are imposed on the rows of the state transition
probabilities and the mixture coefficients of the emission distributions,
i.e.,

$\displaystyle p(\boldsymbol{c}_{i,:})$
$\displaystyle=\mathcal{D}(\boldsymbol{c_{i,:}}|\boldsymbol{w}_{i}),\
\boldsymbol{w}_{i}\in\mathbb{R}^{M},\ i=1,\dots,N$ (5) $\displaystyle
p(\hat{\boldsymbol{A}}_{k,:})$
$\displaystyle=\mathcal{D}(\hat{\boldsymbol{A}}_{k,:}|\boldsymbol{\alpha}_{k}),\boldsymbol{\alpha}_{k}\in\mathbb{R}^{K},k=1,\dots,K$
(6) $\displaystyle p(\boldsymbol{A}_{i,:}^{k})$
$\displaystyle=\mathcal{D}(\boldsymbol{A}_{i,:}^{k}|\boldsymbol{\eta}_{i}^{k}),\boldsymbol{\eta}_{i}^{k}\in\mathbb{R}^{N},i=1,\dots,N,\forall
k$ (7)

For the prior over the means and precision matrices of the mixture components,
we assign a Normal-Wishart distribution with hyperparameters
$\lambda_{ij},\boldsymbol{m}_{ij},\eta_{ij},S_{ij}$,

$\displaystyle
p(\boldsymbol{\mu}_{ij},R_{ij})=\mathcal{NW}(\boldsymbol{\mu}_{ij},R_{ij}|\lambda_{ij},\boldsymbol{m}_{ij},\eta_{ij},S_{ij}),\
\forall i,j$ (8)

This concludes the formulation of the prior specification for the VB treatment
of the CD-HMM model, henceforth referred to as VB-CD-HMM. A graphical
illustration of the fully proposed model can be found in Fig. 1(b).

Let $\theta=\\{\phi,\psi\\}$ be the sets of latent process and emission
distribution parameters where a conjugate exponential prior has been imposed.
Under the VB treatment, we introduce an arbitrary distribution $q(\theta)$ and
derive the Evidence Lower Bound (ELBO) for the model. Maximizing the ELBO is
equivalent to the minimization of the KL divergence between the true and the
postulated posterior. We introduce the mean-field (posterior-independence)
assumption on the joint variational posterior, such that $q(\theta)$
factorizes over all latent variables and model parameters. Then, the ELBO can
be written as

$\displaystyle\begin{split}\text{ELBO}&=\mathbb{E}_{q(\theta),X,Y,Z}\Bigg{[}\log\left(\hat{\pi}_{z_{1}}\pi_{x_{1}}\prod_{t=1}^{T-1}\hat{A}_{z_{t},z_{t+1}}\prod_{t>1}^{T}A_{x_{t},x_{t-z_{t}}}^{z_{t}}\prod_{t=1}^{T}p(\boldsymbol{y}_{t}|x_{t})\right)\Bigg{]}\\\
&\quad+\mathbb{E}_{q(\theta)}\left[\log
p(C,\pi,A,\hat{\pi},\hat{A},\boldsymbol{\mu},\boldsymbol{R})\right]-\mathbb{E}_{q(\theta)}\left[\log
q(\theta)\right]\\\ \end{split}$ (9)

Variational Posteriors. For each variable, the optimal member of the
exponential family111By imposing conjugate priors, we obtain posteriors of the
same functional form as their corresponding priors. can be found by maximizing
(9). This maximization is performed in an EM-like fashion; in the E-step the
parameters of the posterior distributions of the latent processes are updated,
while in the M-step the rest of the variational distributions; the resulting
iterative procedure is guaranteed to monotonically increase the ELBO.
M-step. Starting from Eq. (9), we collect all relevant terms pertaining to the
transition matrix of the first layer process $\hat{A}$; we then maximize the
resulting expression, yielding the following variational posterior:

$\displaystyle
q(\hat{A})=\prod_{i=1}^{K}\mathrm{Dir}(\hat{A}_{i1},\dots,\hat{A}_{iK}|\omega_{i1}^{\hat{A}},\dots,\omega_{iK}^{\hat{A}})$
(10)

where $\omega_{ij}^{\hat{A}}=\alpha_{ij}+\sum_{t}\gamma_{ijt}^{\hat{A}}$,
$\gamma_{ijt}^{\hat{A}}\triangleq q(z_{t}=i,z_{t-1}=j)$.

We follow the same procedure for the initial state probabilities, $\hat{\pi}$,
the corresponding posteriors for the second layer process parameters $\pi$ and
$A=\\{A^{k}\\}_{k=1}^{K}$, and the mixture component weights $C$. Due to
conjugate priors, the emission distributions yield Normal-Wishart
distributions with variational parameters
$\\{\tilde{\lambda}_{ij},\tilde{\boldsymbol{m}}_{ij},\tilde{\eta}_{ij},\tilde{\boldsymbol{S}}_{ij}\\}_{i,j=1}^{N,M}$;
these are similarly obtained via standard computations, i.e., maximizing the
corresponding ELBO expression.
E-step. Turning to the E-step of the iterative procedure, the joint
variational posterior optimizer for the latent processes $Z$ and $X$, and the
mixture component indicators $L$ is analogous to the expression of the
conditional probability defined in the context of conventional HMMs, e.g., see
[9]. Thus, the variational posterior entails the calculation of the
corresponding responsibilities that comprise said posterior $q(z_{1}=i)$,
$q(x_{1}=1)$, $q(z_{t}=i,z_{t-1}=j)$ and $q(x_{t}=j,x_{t-k}=i)$; these can
easily be computed by means of of the well-known forward-backward algorithm.
Thus, in the following, we derive a variant of the forward-backward algorithm
for our model.
Forward-Backward Algorithm. For calculating the forward-backward algorithm, we
first define the forward messages as:

$\displaystyle\alpha_{t}(\\{x_{\tau}\\}_{\tau=t-K+1}^{t},z_{t})\triangleq
p(\\{\boldsymbol{y}_{\tau}\\}_{\tau=1}^{t},\\{x_{\tau}\\}_{\tau=t-K+1}^{t},z_{t})$
(11)

The defined messages can be computed recursively, using the following
initialization:

$\displaystyle\alpha_{1}(x_{1},z_{1})=\begin{cases}\pi_{i}^{*}p^{*}(\boldsymbol{y}_{1}|x_{1}=i),&z_{1}=1\\\
0,&z_{1}>1\end{cases}$ (12)

and the induction step reads:

$\displaystyle\begin{split}\alpha_{t}(\\{x_{\tau}\\}_{\tau=t-K+1}^{t},z_{t})&=p^{*}(\boldsymbol{y}_{t}|x_{t})\sum_{x_{t-K}}\sum_{z_{t-1}}\hat{A}_{z_{t-1},z_{t}}^{*}A_{x_{t-k},x_{t}}^{*z_{t}}\alpha_{t-1}(\\{x_{\tau}\\}_{\tau=t-K}^{t-1},z_{t-1})\end{split}$
(13)

The backward messages are analogously defined as:

$\displaystyle\beta_{t}(\\{x_{\tau}\\}_{\tau=t-K+1}^{t},z_{t})\triangleq
p(\\{\boldsymbol{y}_{\tau}\\}_{\tau=t+1}^{T}|\\{x_{\tau}\\}_{\tau=t-K+1}^{t},z_{t})$
(14)

The corresponding initialization reads:

$\displaystyle\beta_{T}(\\{x_{\tau}\\}_{\tau=T-K+1}^{T},z_{T}=k)=1,\quad\forall
k$ (15)

while the recursion becomes:

$\displaystyle\begin{split}\beta_{t}(\\{x_{\tau}\\}_{\tau=t-K+1}^{T},z_{t})&=\sum_{x_{t+1}}\sum_{z_{t+1}}\hat{A}_{z_{t},z_{t+1}}^{*}p^{*}(\boldsymbol{y}_{t+1}|x_{t+1})\times\\\
&A_{x_{t-k+1},x_{t+1}}^{*z_{t+1}}\beta_{t+1}(\\{x_{\tau}\\}_{\tau=t-K+2}^{T},z_{t+1})\end{split}$
(16)

Responsibilities. By utilizing the forward-backward messages, the necessary
responsibilities for updating the parameters of the model can now be computed.

For the first layer process, the marginal initial state responsibilities
yield:

$\displaystyle\begin{split}\gamma_{t}^{z}(k)&\triangleq
p(z_{t}=k|\\{\boldsymbol{y}_{\tau}\\}_{\tau=1}^{T})\propto\sum_{X}\alpha_{t}(\\{x_{\tau}\\}_{\tau=t-K+1}^{t},z_{t})\beta_{t}(\\{x_{\tau}\\}_{\tau=t-K+1}^{t},z_{t})\end{split}$

Analogously, for the observation-emitting states, we have:

$\displaystyle\begin{split}\gamma_{t}^{x}(i)&\triangleq
p(x_{t}=i|\\{\boldsymbol{y}_{\tau}\\}_{\tau=1}^{T})\propto\sum_{X^{\prime},z_{t}}\alpha_{t}(\\{x_{\tau}\\}_{\tau=t-K+1}^{t},z_{t})\beta_{t}(\\{x_{\tau}\\}_{\tau=t-K+1}^{t},z_{t})\end{split}$

where
$X^{\prime}=\\{x_{\tau}\\}_{\tau=t-K+1}^{t}\setminus\\{x_{t}\\}=\\{x_{\tau}\\}_{\tau=t-K+1}^{t-1}$.

The state transitions responsibilities for the temporal dependence indicators:

$\displaystyle\begin{split}\gamma_{t}^{z}(k,k^{\prime})&\triangleq
p(z_{t}=k^{\prime},z_{t-1}=k|\\{\boldsymbol{y}_{\tau}\\}_{\tau=1}^{T})\\\
&\propto\sum_{X^{\prime}}\alpha_{t-1}(\\{x_{\tau}\\}_{\tau=t-K}^{t-1},z_{t-1})\hat{A}^{*}_{z_{t-1},z_{t}}\
A_{x_{t-k},x_{t}}^{*z_{t}}p^{*}(\boldsymbol{y}_{t}|x_{t})\beta_{t}(\\{x_{\tau}\\}_{\tau=t-K}^{t},z_{t})\end{split}$

Finally, the dependent state transition responsibilities:

$\displaystyle\begin{split}\gamma_{t}^{x}(i,&j,k)\triangleq
p(x_{t}=j,x_{t-k}=i|\\{\boldsymbol{y}_{\tau}\\}_{\tau=1}^{T})\propto\sum_{X^{\prime\prime},z_{t}}\alpha_{t}(\\{x_{\tau}\\}_{\tau=t-K+1}^{t},z_{t})\beta_{t}(\\{x_{\tau}\\}_{\tau=t-K+1}^{t},z_{t})\end{split}$

where
$X^{\prime\prime}=\\{x_{\tau}\\}_{\tau=t-K+1}^{t}\setminus\\{x_{t},x_{t-k}\\}$.

###  2.3 Inference

Let us consider a VB-CD-HMM model, fully trained with data $Y$. We seek to
calculate the predictive density of the test sequence $Y^{\mathrm{test}}$:

$\displaystyle p(Y^{\text{test}}|Y)=\int d\hat{\theta}\
p(\hat{\theta}|Y)P(Y^{\text{test}}|\hat{\theta})$ (17)

Using the introduced variational posterior in place of the unknown true
posterior $p(\hat{\theta}|Y)$ and similar to [9], we can compute the necessary
quantities by utilizing the forward algorithm variant defined in Section 2.2.
Specifically, the density of a given test set
$Y^{\text{test}}=\\{\boldsymbol{y}_{t}^{\text{test}}\\}_{t=1}^{T}$ yields:

$\displaystyle
p(Y^{\text{test}}|\hat{\theta})=\sum_{z_{T}}\sum_{X}\alpha_{T}(x_{T-K+1},\dots,x_{T},z_{T})$
(18)

##  3 Experimental Evaluation

###  3.1 Experimental Details

We use four popular benchmark datasets for the human action recognition task,
namely MSR Action 3D (MSRA) [15], UTD-MHAD (UTD) [11], Gaming 3D (G3D) [4] and
UPenn Action (Penn) [27]. Only skeletal data are used for training the models.
The preprocessing and augmentation of the data are the same as in [29]. We
extract the motion of the joints for every pair of joints between two
consecutive frames. PCA is then employed for dimensionality reduction.

In addition, we further consider evaluation with a more complex dataset,
namely NTU [21]. This comprises $\approx 57000$ samples of $60$ different
actions. We adopt the two data splits suggested in the original publication:
(i) Cross-View and (ii) Cross-Subject.

###  3.2 Parameter Initialization & Hyperparameter Selection

The resulting ELBO expression is non-convex with respect to the variational
posterior, and potentially many local maxima exist; the obtained solution is
thus dependent on parameter initialization. To avoid the rather time-consuming
procedure of multiple runs from different initial values, we employ a common
initialization scheme: K-Means is used to initialize the posterior parameters
of the second layer process and the means and variances of the distributions.
The parameters of the first postulated layer are initialized to random values,
while for the hyperparameters we use ad-hoc uninformative values, similar to
[9]. Specifically, the hyperparameters of the Dirichlet priors are set to very
small values ($10^{-3}$). For the parameters of the emission distributions, we
set $\lambda_{ij}=\tilde{\lambda}_{ij}=0.25$ and
$\eta_{ij}=\tilde{\eta}_{ij}=D+2$, where $D$ is the dimensionality of the
data.

###  3.3 States, Mixtures and Temporal Dependencies

We trained different models with varying number of states and mixtures to
investigate their effect on model effectiveness. We observed that the
recognition accuracy is not considerably sensitive to the chosen
configuration; thus, we run multiple configurations and retain the one with
the higher ELBO. For all experiments, we set $K=2$, as we observed that the
inferred temporal dependence posteriors rarely assigned high probability (if
any) to $p(z_{t}>2|z_{t-1})$; this behavior is reasonable considering the
physical meaning of the used data features. Nevertheless, we did not observe
any drop in accuracy when setting $K>2$. For $K=1$, the model reduces to a
simple HMM, and exhibits a significant drop in recognition accuracy. This
outcome strongly corroborates the benefits of our dynamic temporal dependency
inference approach.

###  3.4 Experimental Results

We follow the train and test splits as suggested by the authors of the
respective datasets. We train one model per class and assign the test data to
the class with the maximum predictive posterior (Eq. (18)). We compare our
approach to related approaches, along with a comparison with HDM [29]. Out of
the HDM model variants therein, we focus on the BV variant for comparability,
whereby variational inference is utilized for computing the predictive
likelihood.

Table 1: Recognition Accuracy ($\%$) for individual dataset experiments.

Model | MSRA | UTD | G3D | Penn | Avg.
---|---|---|---|---|---
HMM | 67.8 | 82.8 | 68.1 | 82.3 | 75.3
HMM2 | 80.2 | 83.1 | 82.6 | 84.4 | 82.6
HSMM | 66.3 | 82.3 | 77.5 | 78.9 | 76.35
LSTM | 74.7 | 77.0 | 82.2 | 90.3 | 81.1
HDM-PI | 70.3 | 84.4 | 79.4 | 89.8 | 81.0
HDM-PL | 80.6 | 90.2 | 87.7 | 91.6 | 87.5
HDM-BV | 82.1 | 91.4 | 87.7 | 90.8 | 88.0
VB-CD-HMM | $\mathbf{82.5}$ | $\mathbf{92.7}$ | $\mathbf{90.6}$ | $\mathbf{92.0}$ | $\mathbf{89.45}$

Individual Datasets. We first consider the MRSA dataset, where the recognition
rates of our approach can be found in the first column of Table 1. Compared to
baseline models such as HMMs, HSMMs and LSTMs, our model performs much better
in classification accuracy with an average improvement of $12.9\%$. Even
though the recently proposed HDM [29] consistently improves over the
considered alternatives, it falls short compared to the proposed VB-CD-HMM
model. This behavior is consistent across all the considered datasets,
resulting in an average classification accuracy of $89.45\%$, outperforming
HDM by $1.45\%$. Additionally, we observe that our proposed approach
consistently and significantly outperforms a second-order HMM (HMM2) evaluated
under the same experimental and modeling setup. Even though HMM2 improves over
simple first-order HMMs, by considering more complex dynamics, it still falls
short to the flexibility of modeling dependencies over long horizons, that VB-
CD-HMM offers. Note also that the prowness of HMM2 comes at the cost of a
considerable increase in the resulting computational complexity; this is in
contrast to our model, which still relies on (a variant of) the forward-
backward algorithm.

We additionally compare the recognition accuracy of our approach to other
state-of-the-art methods. The corresponding performances are illustrated in
Table 2. As we observe, our approach yields clearly improved accuracy compared
to the competition for UTD, Penn and NTU datasets. For the G3D dataset, our
method slightly outperforms LRBM [19], but R3DG [24] performs better; this is
due to the sophisticated feature engineering and combination of several
approaches in contrast to our simple joints locations and motion features.

Table 2: Recognition accuracy for all the considered datasets compared to
state-of-the-art.

Dataset
---
Method $\parallel$ Acc
MSRA | UTD | G3D | Penn | NTU (x-view) | NTU (x-subject)
ST-LSTM[16] $\parallel$ 94.8 | Fusion[11] $\parallel$ 79.1 | LRBM [19] $\parallel$ 90.5 | Actemes[27] $\parallel$ 86.5 | ST-TSL [22] $\parallel$ 92.4 | ST-TSL[22] $\parallel$ 84.8
B-GC-LSTM[28] $\parallel$ 94.5 | SOS-CNN[14] $\parallel$ 87.0 | R3DG[24] $\parallel$ 91.1 | AOG[18] $\parallel$ 84.8 | B-GC-LSTM [28] $\parallel$ 89.0 | B-GC-LSTM [28] $\parallel$ 81.8
VB-CDHMM $\parallel$ 92.5 | VB-CDHMM $\parallel$ 92.7 | VB-CDHMM $\parallel$ 90.6 | VB-CDHMM $\parallel$ 92.0 | VB-CDHMM $\parallel$ 93.0 | VB-CDHMM $\parallel$ 85.0

The main focus of this Section is to compare and prove the efficacy of our
approach within the family of more general HMM-based methods in the tackled
application domain. Thus, we did not explore complicated model-based and
specific-task tailored feature engineering; this allows for both transparency
of comparison to recent related published work, and interpretability of our
results. Nevertheless, as we shall see in the next section, our method clearly
outperforms R3DG, when we randomly omit observations; this behavior vouches
for the improved robustness of our approach. Likewise, ST-LSTM [16]
incorporates complicated feature extraction mechanisms, thus explaining the
performance gap compared to our approach. For the NTU dataset [21], to the
best of our knowledge, there are no results for recent HMM-based methods. As
we observe, compared to the state-of-the-art methods ST-LSTM [16] and B-GC-
LSTM [28], VB-CD-HMM yields considerable improvements in both splits, despite
its simplicity and lower computational complexity.

###  3.5 Missing Values

Generative models come with the additional benefit of robustness to missing
values. This advantage is of great importance in the human action recognition
field, especially when using skeletal data, where the dataset may be corrupted
with missing observations, e.g., due to occlusion. As an HMM variant, the
proposed model is such a generative model; we thus assess the efficacy of our
approach when missing values are present. To this end, we train the considered
architecture in the UTD, MSRA and G3D datasets, with randomly missing values
from both the train and test data. Analogous experiments have been performed
in [29], and we adopt the same procedure, including the proportion of missing
values with respect to the original data.

Table 3: Recognition Accuracy ($\%$) with missing values.

Dataset | UTD | MSRA | G3D
---|---|---|---
Missing Portion | $10\%$ | $30\%$ | $50\%$ | $10\%$ | $30\%$ | $50\%$ | $10\%$ | $30\%$ | $50\%$
Model | R3DG [24] | $81.5$ | $74.0$ | $72.0$ | $78.0$ | $72.0$ | $70.0$ | $87.0$ | $86.0$ | $83.0$
DLSTM [30] | $70.5$ | $66.0$ | $63.0$ | $68.0$ | $63.0$ | $61.0$ | $81.0$ | $76.0$ | $73.0$
HDM [29] | $91.0$ | $90.5$ | $90.0$ | $80.5$ | $78.0$ | $76.0$ | $90.0$ | $89.0$ | $88.0$
VB-CD-HMM | $\mathbf{92.55}$ | $\mathbf{91.6}$ | $\mathbf{90.2}$ | $\mathbf{81.7}$ | $\mathbf{80.1}$ | $\mathbf{79.1}$ | $\mathbf{90.2}$ | $\mathbf{89.3}$ | $\mathbf{88.5}$

Thus, in Table 3, we report the recognition rates when we randomly omit
$10\%,30\%$ and $50\%$ of the observations. As we observe, our method clearly
outperforms the R3DG [24] and DLSTM [30] methods by a clear margin. The same
behavior is observed compared to the more relative HDM approach [29]. Note
that even though the reported HDM accuracies are obtained through Gibbs
sampling (BG variant), our approach exhibits the smallest decrease in
recognition accuracy relative to the increase of missing values.

Further Insights. Finally, we examine the patterns of the first layer process,
that is the dependence generator layer of the VB-CD-HMM model. We gain further
insights to the behavior of the model and assert that the temporal
dependencies do not collapse to simple first order dynamics (which would
effectively reduce the model to a conventional HMM). To this end, we focus on
a trained model on the UTD dataset and choose a subset of ten actions to
investigate their posterior parameters concerning the generation of
dependencies. As is clearly shown in Fig. 2, the distribution of temporal
dependencies is essentially different for each action, providing strong
empirical evidence that the introduced mechanism can capture complex varying
temporal dynamics in the data.

Figure 2: The inferred temporal dependence probabilities $A_{kk^{\prime}}$ for
the first $10$ actions in the UTD dataset. Black denotes very high
probability, while white very low. Actions with similar patterns, e.g., 1
(swipe left), 2 (swipe right), 5 (throw), i.e., comprising sequential
movements to “one” direction, are dominantly dependent on the previous frame;
in contrast, more complicated movements, e.g., 3 (wave), 4 (clap), mainly
depend on two steps back, possibly dominated by the reversal in direction.
Even more complicated actions, as 7 (basketball shoot) and 8 (draw x) exhibit
more complex temporal dependencies, that are captured via the resulting
inferred patterns.

##  4 Conclusions

Our experimental results have provided strong empirical evidence for the
efficacy of our approach, yielding competitive recognition accuracy in the
Human Action Recognition task. Compared to the commonly considered LSTM-based
implementations, our proposed architecture:

  * •

offers the flexibility to model complex sequential data accompanied by the
interpretability of its discrete bottlenecks; this is in stark contrast to
LSTM-based methods that assume continuous states that essentially entail no
probabilistic interpretation.

  * •

enjoys closed-form updates that allow for significantly lower computational
and training complexity, devoid of complications, e.g., vanishing gradients,
and

  * •

can effectively handle data with missing values, exhibiting the smallest
decrease in recognition performance relative to other approaches.

We emphasize that the proposed mechanism for inferring the dependencies
between the states in our proposed model allows for analyzing the inferred
temporal dynamics. Thus, we can gain significant insights and examine how the
dynamics adapt to the complexity and variations of different train/test
sequences. Therefore, this capacity can be of great benefit to the current
effort of the community towards interpretable machine learning models.

## References

  * [1] Aggarwal, J., Xia, L.: Human activity recognition from 3d data: A review. Pattern Recognition Letters 48 (10 2014)
  * [2] Aycard, O., Mari, J., Washington, R.: Learning to automatically detect features for mobile robots using second-order hidden markov models. IJARS (2005)
  * [3] Begleiter, R., El-Yaniv, R., Yona, G.: On prediction using variable order markov models. J. Artif. Int. Res. 22 (Dec 2004)
  * [4] Bloom, V., Makris, D., Argyriou, V.: G3d: A gaming action dataset and real time action recognition evaluation framework. In: CVPRW (June 2012)
  * [5] Bühlmann, P., Wyner, A.J.: Variable length markov chains. Ann. Statist. (1999)
  * [6] Chatzis, S.P., Demiris, Y.: Nonparametric mixtures of gaussian processes with power-law behavior. IEEE Transactions on Neural Networks and Learning Systems 23(12), 1862–1871 (2012)
  * [7] Chatzis, S.P., Kosmopoulos, D.: A latent manifold markovian dynamics gaussian process. IEEE Transactions on Neural Networks and Learning Systems (2015)
  * [8] Chatzis, S., Kosmopoulos, D., Papadourakis, G.: A nonstationary hidden markov model with approximately infinitely-long time-dependencies. IJAIT 25 (2016)
  * [9] Chatzis, S.P.: A variational bayesian methodology for hidden markov models utilizing student’s-t mixtures. Pattern Recognition 44 (2011)
  * [10] Chatzis, S.P.: Margin-maximizing classification of sequential data with infinitely-long temporal dependencies. Expert Systems with Applications 40 (2013)
  * [11] Chen, C., Jafari, R., Kehtarnavaz, N.: Utd-mhad: A multimodal dataset for human action recognition utilizing a depth camera and a wearable inertial sensor. In: Procs. ICIP (Sep 2015)
  * [12] Dimitrakakis, C.: Bayesian variable order markov models. In: Procs. ICAIS 13
  * [13] Engelbrecht, H., du Preez, J.: Efficient backward decoding of high-order hidden markov models. Pattern Recognition 43 (2010)
  * [14] Hou, Y., Li, Z., Wang, P., Li, W.: Skeleton optical spectra-based action recognition using convolutional neural networks. IEEE TCSVT 28(3), 807–811 (2018)
  * [15] Li, W., Zhang, Z., Liu, Z.: Action recognition based on a bag of 3d points. In: IEEE CVPRW 2010 (07 2010)
  * [16] Liu, J., Shahroudy, A., Xu, D., Kot, A., Wang, G.: Skeleton-based action recognition using spatio-temporal lstm network with trust gates. IEEE Transactions on Pattern Analysis and Machine Intelligence 40, 3007–3021 (12 2018)
  * [17] Nel, E., du Preez, J.A., Herbst, B.M.: Estimating the pen trajectories of static signatures using hidden markov models. IEEE Transactions on Pattern Analysis and Machine Intelligence 27 (Nov 2005)
  * [18] Nie, B.X., Xiong, C., Zhu, S.: Joint action recognition and pose estimation from video. In: Procs. CVPR. pp. 1293–1301 (June 2015)
  * [19] Nie, S., Wang, Z., Ji, Q.: A generative restricted boltzmann machine based method for high-dimensional motion data modeling. Comput. Vis. Image Underst. (2015)
  * [20] Panousis, K., Chatzis, S., Theodoridis, S.: Nonparametric Bayesian deep networks with local competition. Procs ICML, PMLR (June 2019)
  * [21] Shahroudy, A., Liu, J., Ng, T., Wang, G.: Ntu rgb+d: A large scale dataset for 3d human activity analysis. In: Procs. CVPR. pp. 1010–1019 (2016)
  * [22] Si, C., Jing, Y., Wang, W., Wang, L., Tan, T.: Skeleton-Based Action Recognition with Spatial Reasoning and Temporal Stack Learning: 15th EC, Munich, Germany, September, 2018, Proceedings, Part I, pp. 106–121 (09 2018)
  * [23] Theodoridis, S.: Machine Learning: A Bayesian and Optimization Perspective. Academic Press, Inc., USA, 2nd edn. (2020)
  * [24] Vemulapalli, R., Arrate, F., Chellappa, R.: Human action recognition by representing 3d skeletons as points in a lie group. In: Procs. CVPR (June 2014)
  * [25] Weinland, D., Ronfard, R., Boyer, E.: A survey of vision-based methods for action representation, segmentation and recognition. Comput. Vis. Image Underst. (2011)
  * [26] Yu, S.Z.: Hidden semi-markov models. Artif. Intell. 174 (2010)
  * [27] Zhang, W., Zhu, M., Derpanis, K.G.: From actemes to action: A strongly-supervised representation for detailed action understanding. In: ICCV (Dec 2013)
  * [28] Zhao, R., Wang, K., Su, H., Ji, Q.: Bayesian graph convolution lstm for skeleton based action recognition. In: Procs. ICCV. pp. 6881–6891 (2019)
  * [29] Zhao, R., Xu, W., Su, H., Ji, Q.: Bayesian hierarchical dynamic model for human action recognition. In: Procs. CVPR (June 2019)
  * [30] Zhu, W., Lan, C., Xing, J., Zeng, W., Li, Y., Shen, L., Xie, X.: Co-occurrence feature learning for skeleton based action recognition using regularized deep lstm networks. In: Procs. AAAI (2016)


#  Adaptive Gaussian Process Regression for Bayesian inverse problems††thanks:
This work has been supported by Bundesministerium für Bildung und Forschung –
BMBF, project number 05M20ZAA (siMLopt) and by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) – 436400679.

Paolo Villani Zuse Institut Berlin<EMAIL_ADDRESS>Jörg Unger
Bundesanstalt für Materialforschung und -prüfung<EMAIL_ADDRESS>Martin
Weiser22footnotemark: 2

###### Abstract

We introduce a novel adaptive Gaussian Process Regression (GPR) methodology
for efficient construction of surrogate models for Bayesian inverse problems
with expensive forward model evaluations. An adaptive design strategy focuses
on optimizing both the positioning and simulation accuracy of training data in
order to reduce the computational cost of simulating training data without
compromising the fidelity of the posterior distributions of parameters. The
method interleaves a goal-oriented active learning algorithm selecting
evaluation points and tolerances based on the expected impact on the Kullback-
Leibler divergence of surrogated and true posterior with a Markov Chain Monte
Carlo sampling of the posterior. The performance benefit of the adaptive
approach is demonstrated for two simple test problems.

Keywords: Gaussian process regression, Bayesian inverse problems, surrogate
models, parameter identification, active learning

MSC 2010: 60G15, 62F15, 62F35, 65N21

##  1 Introduction

The inverse problem of inferring the posterior probability of parameters
$p\in\mathbb{R}^{d}$ in a forward model $y(p)$ from measurements
$y^{m}\in\mathbb{R}^{m}$ is often addressed by sampling with Markov Chain
Monte Carlo (MCMC) methods [5]. The large number of forward evaluations
required for a faithful representation of the posterior density renders this
inapplicable in case of computationally expensive forward models such as large
finite element (FE) simulations. The forward model is thus often replaced by a
fast surrogate model when sampling the posterior. Here, we focus on the
efficient construction of Gaussian Process Regression (GPR) surrogates.

Surrogate models are learned from values $y(p_{i})$ at specific evaluation
points $p_{i}$ as training data. The accuracy of the resulting surrogate
depends on the number and position of the sample points. Constructing an
accurate surrogate model can become computationally expensive when a large
number of evaluations is required. Consequently, strategies for selecting
near-optimal evaluation points have been proposed for various settings [11]. A
priori point sets [4, 10] are effectively supplemented by adaptive designs [3,
6, 8, 16] selecting the most beneficial evaluation points $p_{i}$.

When using FE simulations for computing training data, the evaluations of
$y(p_{i})$ are affected by discretization and truncation errors. The trade-off
between accuracy and cost has been investigated using different low and high
fidelity models [9], and by an adaptive choice of evaluation tolerances [12,
13, 14] in different settings. Here, we extend [13] from an offline training
for maximum posterior point estimates to an interleaved posterior sampling and
surrogate training driven by a goal-oriented approach.

##  2 Gaussian Process regression

Gaussian process regression is a regression technique which allows to
approximate any function, naturally fits the Bayesian framework, and provides
an uncertainty estimate of its prediction.

We consider a forward model $y:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}$, which
we assume to be a realisation of a Gaussian process $\mathcal{G}$ with mean
$\mu_{0}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}$ and covariance kernel
$k:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{m\times m}$ to be
defined later.

For training points $(p_{i},y_{i})_{i=1,\dots,s}$ with $y_{i}\approx y(p_{i})$
of accuracy $\tau_{i}\geq 0$, we are interested in a prediction of
$y_{s+1}\approx y(p_{s+1})$ for any $p_{s+1}$. The GPR posterior covariance
block matrix is $\Gamma=(K^{-1}+T^{-2})^{-1}\in\mathbb{R}^{m(s+1)\times
m(s+1)}$ with prior covariance blocks $K_{ij}=k(p_{i},p_{j})$ and formally
likelihood covariance $T=\mathrm{diag}(\tau_{1}I,\dots,\tau_{s}I,\infty I)$.
The GPR posterior mean is $\bar{Y}=\Gamma(K^{-1}M_{0}+T^{-2}Y)$ with
$Y=(y_{1},\dots,y_{s},0)$. Then, the GPR prediction is the marginal normal
distribution $y_{s+1}\sim\mathcal{N}(\bar{Y}_{s+1},\Gamma_{s+1,s+1})$. As
$p_{s+1}\in\Omega$ is arbitrary, this defines mean
$\bar{y}:\Omega\to\mathbb{R}^{m}$ and covariance
$\Gamma:\Omega\to\mathbb{R}^{m\times m}$ on the whole parameter space. We
refer to [11, 13] for a more detailed exposition.

##  3 Bayesian surrogate-based parameter identification

We consider the forward model
$y:\Omega\subset\mathbb{R}^{d}\to\mathbb{R}^{m}$, which cannot be evaluated
directly, but can be approximated through a numerical procedure $y_{\tau}$
with arbitrary precision in exchange of computational work: We assume that for
any $\tau>0$, we obtain an evaluation $y_{\tau}(p)\sim\mathcal{N}(y(p),\tau
I)$, with cost $W_{\tau}$.

We assume measurements $y^{m}$ to be random variables generated by a linear
additive Gaussian noise model

$y^{m}=y(p)+\eta$ (1)

with $\eta\sim\mathcal{N}(0,\Sigma_{l})$. For simplicity, we consider a
diagonal covariance structure
$\Sigma=\text{diag}(\sigma_{1},\dots\sigma_{n})$, corresponding to independent
noise components. The conditional distribution of the measurements is then
$y^{m}\mid p\sim\mathcal{N}(y(p),\Sigma_{l})$,

$\pi(y^{m}\mid
p)=(2\pi)^{-m/2}\det(\Sigma_{l})^{-1/2}\exp\Big{(}-\frac{1}{2}\lVert
y^{m}-y(p)\rVert_{\Sigma_{l}^{-1}}^{2}\Big{)}$

is the likelihood of the problem. Evaluating the likelihood requires
evaluating the forward model $y$, which we assume to be computationally
expensive.

To reduce costs, we assume that $y$ is a realisation of a GP, and introduce a
GP surrogate model $\mathcal{G}$ of predictive mean $\bar{y}$ and variance
$\Gamma$. For simplicity, we consider a surrogate with independent output
components, i.e. diagonal covariance $\Gamma(p)$. The training points for this
GP are given by numerical evaluations $y_{\tau_{i}}(p_{i})$ of the forward
model. These points and the corresponding evaluation tolerances $\tau_{i}$
form the training design $\mathcal{D}$. We postpone the question of how to
build training designs to the next section.

To evaluate the likelihood, we could substitute the forward model with the
mean estimate $\bar{y}$, obtaining

$\pi_{\text{plug-in}}(y^{m}\mid
p,\bar{y})=(2\pi)^{-m/2}\det(\Sigma_{l})^{-1/2}\exp\Big{(}-\frac{1}{2}\lVert
y^{m}-\bar{y}(p)\rVert_{\Sigma_{l}^{-1}}^{2}\Big{)}.$ (2)

This, from a decision-theoretic point of view, corresponds to the minimisation
of the $L^{1}$ loss [7], but ignores the uncertainty estimate given by the
predictive variance: since $y$ is assumed to be a realisation of
$\mathcal{G}$, the measurement noise model (1) becomes
$y^{m}=\mathcal{G}(p)+\eta$. Marginalizing over GP realizations results in a
different conditional distribution of the measurements $y^{m}\mid
p,\mathcal{G}\sim\mathcal{N}(\bar{y}(p),\Sigma_{l}+\Gamma(p))$ and in a
marginal likelihood:

$\pi_{\mathcal{D}}(y^{m}\mid
p,\mathcal{D})=(2\pi)^{-m/2}\det\left(\Sigma_{l}+\Gamma(p)\right)^{-\frac{1}{2}}\exp\left(-\frac{1}{2}\lVert
y_{m}-\bar{y}(p)\rVert_{\left(\Sigma_{l}+\Gamma(p)\right)^{-1}}^{2}\right),$
(3)

see, e.g., [2]. Note that the conditional distribution is still Gaussian due
to the normality of both the noise and the GP. Moreover, the likelihood
$\pi_{\mathcal{D}}$ is closely related to the $L^{2}$ loss [7, 15]. Including
the GP variance into the likelihood can be important for avoiding
overconfident yet wrong posterior approximations by surrogated forward models,
see Fig. 1 for an illustration.

By adopting a Bayesian point of view, we express prior belief on the parameter
by assigning a prior distribution $\pi(p)$. Then, by Bayes’ theorem, we obtain
a true posterior distribution

$\pi(p\mid y^{m})=\frac{\pi(p)\ \pi(y^{m}\mid p)}{\pi(y^{m})},$ (4)

corresponding to the true likelihood $\pi(y^{m}\mid p)$ and an approximate
posterior

$\pi(p\mid y^{m},\mathcal{D})=\frac{\pi(p)\ \pi(y^{m}\mid
p,\mathcal{D})}{\pi(y^{m}\mid\mathcal{D})},$ (5)

corresponding to the likelihood approximation $\pi(y^{m}\mid p,\mathcal{D})$
as given in (2) and (3), respectively.

Figure 1: Impact of the likelihoods (2) and (3) on the posterior for an
illustrative inverse problem problem with forward model $y(p)=p^{2}sin(p)$ and
uniform prior on parameter space $[0,1]$. The horizontal lines show the actual
measurement and the $2\sigma$ range of measurement noise. The marginal
likelihood (3) is wider due to including the GP variance, and avoids
overconfident posteriors.

In both cases, the normalising constant $\pi(y^{m})$ or
$\pi(y^{m}\mid\mathcal{D})$, respectively, will not be computationally
available, as it requires integration over the parameter space $\Omega$:
fortunately, it is not needed for posterior sampling by Markov-Chain Monte
Carlo (MCMC) methods.

##  4 Posterior-oriented surrogate model

As in [15], we do not aim at building a surrogate which is globally accurate
on the whole parameter space $\Omega$, but at finding a design $\mathcal{D}$
such that the approximate posterior is accurate, i.e. $\pi(p\mid
y^{m})\approx\pi(p\mid y^{m},\mathcal{D})$. Repeatedly selecting training
points randomly sampled from $\pi(p\mid y^{m},\mathcal{D})$, updating
$\mathcal{G}$ and then iterating is sufficient for convergence of $\pi(p\mid
y^{m},\mathcal{D})$ to $\pi(p\mid y^{m})$ in the Hellinger metric [2]. Here,
we also aim at finding a design $\mathcal{D}$ which incurs a small
computational cost of evaluating training data $y_{\tau_{i}}$.

We measure the deviation of the surrogated and the true posterior densities by
the Kullback-Leibler (KL) divergence

$\displaystyle D_{\rm KL}\left(\pi(\cdot\mid y_{m})\mid\pi(\cdot\mid
y_{m},\mathcal{D})\right)$ $\displaystyle=\mathbb{E}_{\pi(p\mid
y^{m})}\left[\log\frac{\pi(p\mid y_{m})}{\pi(p\mid y_{m},\mathcal{D})}\right]$
$\displaystyle=\int_{\Omega}\pi(p\mid y_{m})\log\frac{\pi(p\mid
y_{m})}{\pi(p\mid y_{m},\mathcal{D})}\ dp.$ (6)

Since computing the KL divergence requires evaluating the full model, we
derive a numerical approximation which relies on the surrogate only. Using the
marginal likelihood (3) and the posteriors (4) and (5), their logarithmic
ratio can be written as

$\log\frac{\pi(p\mid y^{m})}{\pi(p\mid
y^{m},\mathcal{D})}=\log\frac{\pi(y^{m}\mid p)}{\pi_{\mathcal{D}}(y^{m}\mid
p,\mathcal{D})}-\log\frac{\pi(y^{m})}{\pi(y^{m}\mid\mathcal{D})}.$

The first term, the logarithmic ratio of true and surrogated likelihood,
equals

$\displaystyle\log$ $\displaystyle\frac{\pi(y^{m}\mid
p)}{\pi_{\mathcal{D}}(y_{m}\mid p,\mathcal{D})}$
$\displaystyle=\frac{1}{2}\Bigg{(}\log\frac{\text{det}\left(\Sigma_{l}+\Gamma(p)\right)}{\text{det}\left(\Sigma_{l}\right)}-\lVert
y(p)-y^{m}\rVert_{\Sigma_{l}^{-1}}^{2}+\lVert\bar{y}(p)-y^{m}\rVert_{\left(\Sigma_{l}+\Gamma(p)\right)^{-1}}^{2}\Bigg{)}.$

As $\Sigma_{l}^{-1}-\left(\Sigma_{l}+\Gamma(p)\right)^{-1}\succeq 0$, we can
upper bound the difference between norms by

$\displaystyle-\lVert y(p)$
$\displaystyle-y^{m}\rVert_{\Sigma_{l}^{-1}}^{2}+\lVert\bar{y}(p)-y^{m}\rVert_{(\Sigma_{l}+\Gamma(p))^{-1}}^{2}$
$\displaystyle\leq-\lVert
y(p)-y^{m}\rVert_{\Sigma_{l}^{-1}}^{2}+\lVert\bar{y}(p)-y^{m}\rVert_{\Sigma_{l}^{-1}}^{2}$
$\displaystyle=-\lVert
y(p)-\bar{y}(p)\rVert_{\Sigma_{l}^{-1}}^{2}-2\left(\bar{y}(p)-y^{m}\right)^{T}\Sigma_{l}^{-1}\left(y(p)-\bar{y}(p)\right).$

By assuming that $y$ is a realisation of $\mathcal{G}$,
$\mathbb{E}\left[\left(y^{(i)}(p)-\bar{y}^{(i)}(p)\right)^{2}\right]=\Gamma^{(i,i)}(p)$
and therefore $\lVert
y(p)-\bar{y}(p)\rVert_{\Sigma_{l}^{-1}}^{2}\approx\mathop{\mathrm{missing}}{tr}\left(\Sigma_{l}^{-1}\Gamma(p)\right)$
hold. Defining
$v=\Sigma_{l}^{-1}\sqrt{\mathop{\mathrm{missing}}{diag}\left(\Gamma(p)\right)}\in\mathbb{R}^{m}$,
we obtain

$\displaystyle-\lVert
y(p)-y^{m}\rVert_{\Sigma_{l}^{-1}}^{2}+\lVert\bar{y}(p)-y^{m}\rVert_{(\Sigma_{l}+\Gamma(p))^{-1}}^{2}\lesssim-\mathop{\mathrm{missing}}{tr}\left(\Sigma_{l}^{-1}\Gamma(p)\right)+2\left|\bar{y}(p)-y^{m}\right|^{T}v.$

We therefore define the local error quantity

$\displaystyle e_{\mathcal{D}}(p)$
$\displaystyle:=\frac{1}{2}\left(\log\det(I+\Sigma_{l}^{-1}\Gamma(p))-\mathop{\mathrm{missing}}{tr}\left(\Sigma_{l}^{-1}\Gamma(p)(I-\Sigma_{l}^{-1})\right)+2\left|\bar{y}(p)-y^{m}\right|^{T}v\right)$
(7) $\displaystyle\gtrsim\log\frac{\pi(y^{m}\mid
p)}{\pi_{\mathcal{D}}(y_{m}\mid p,\mathcal{D})}$

as an approximate upper bound on the log ratio of true and surrogated
likelihood.

By optimistically assuming that the normalisation factors are similar
independent of the design $\mathcal{D}$, and thus
$\log\frac{\pi(y^{m})}{\pi(y^{m}\mid\mathcal{D})}\approx 0$, we substitute (7)
into (6) and obtain the global error quantity

$\displaystyle E(\mathcal{D})=\int_{\Omega}e_{\mathcal{D}}(p)\pi(p\mid y^{m})\
dp.$ (8)

To create an optimal surrogate model, we aim at a training design
$\mathcal{D}$ which minimises $E(\mathcal{D})$ under a computational work
constraint. By denoting the computational work needed to realize $\mathcal{D}$
by $W(\mathcal{D})$ , for a given budget $W$ we aim at solving the
optimisation problem

$\min_{\mathcal{D}}E(\mathcal{D})\ \text{ subject to }\ W(\mathcal{D})\leq W.$
(9)

##  5 Sequential design of experiments

It is far from trivial to predict a priori how design choices impact the error
quantity $E$, especially when a large budget $W$ is available or the initial
surrogate is unreliable. Fortunately, an exact solution of (9) is not needed –
an approximate solution will do, even if it yields a slightly less efficient
design. We follow [13, 14] and adopt a greedy sequential approach, where the
budget $W=\sum_{j=1}^{J}\Delta W_{j}$ is partitioned and sequentially spent.

We start from an initial design $\mathcal{D}_{0}$ and then, for $j=1,\dots,J$,
aim at solving

$\min_{\mathcal{D}_{j}\leq\mathcal{D}_{j-1}}E(\mathcal{D}_{j})\quad\text{s.t.}\quad
W(\mathcal{D}_{j}\mid\mathcal{D}_{j-1})\leq\Delta W_{j}.$ (10)

We write $\mathcal{D}\leq\mathcal{D}_{j-1}$ for any design $\mathcal{D}$ which
refines $\mathcal{D}_{j-1}$ in the sense that it includes all evaluation
points $p_{i}$ contained in $\mathcal{D}_{j-1}$ with lesser or equal
tolerances $\tau_{i}$. We write
$W(\mathcal{D}\mid\mathcal{D}_{j-1})=W(\mathcal{D})-W(\mathcal{D}_{j-1})$ for
the work needed to obtain $\mathcal{D}$ from $\mathcal{D}_{j-1}$.

Even this sequential formulation is highly non-linear and non-convex. An
accurate solution would require a considerable amount of computational work,
possibly exceeding the savings in computational budget possible with a better
design. Consequently, we adopt the heuristic approach of separating the
selection of new candidate evaluation points from the optimisation of the
evaluation tolerances. In the latter, we also decide about the actual
inclusion of the new points in the training set.

_Candidate points._ We choose points where spending computational budget is
likely to reduce the error most. In order to do so, we look at the sensitivity
of the global error $E$ with respect to a reduction of training error at a
candiate position $p^{\prime}$ [14]. This is given by

$\displaystyle\frac{dE({\mathcal{D}})}{dW(p^{\prime})}$
$\displaystyle=\int_{\Omega}\frac{de_{{\mathcal{D}}}(p)}{dW(p^{\prime})}\pi_{\mathcal{D}}(p\mid
y_{m})\,dp$
$\displaystyle=\int_{\Omega}\frac{de_{\mathcal{D}}(p)}{d\Gamma(p)}\frac{d\Gamma(p)}{d\tau(p^{\prime})}\bigg{|}_{\tau=\tau^{\prime}}\frac{d\tau(p^{\prime})}{dW(p^{\prime})}\bigg{|}_{\tau=\tau^{\prime}}\pi_{\mathcal{D}}(p\mid
y_{m})\ dp,$ (11)

where the linearization tolerance $\tau^{\prime}$ is the current GP standard
deviation at point $p^{\prime}$. We adopt (11) as a utility function and
select local minimizers of $\frac{dE(\mathcal{D}_{j-1})}{dW}$ as next
candidate points.

The optimisation problem is solved approximately via a multistart pattern
search. Quadrature is performed by Monte Carlo integration on samples
$\mathcal{S}_{j}$ to be defined in Sec. 6 below. This results in the numerical
utility function

$\frac{dE({\mathcal{D}_{j-1}})}{dW(p^{\prime})}\approx\frac{1}{|\mathcal{S}_{j}|}\sum_{p\in\mathcal{S}_{j}}\frac{de_{\mathcal{D}_{j-1}}(p)}{dW(p^{\prime})}.$

If more than $c_{j}$ local maxima are found, the best $c_{j}$ ones are
selected as candidates; if less are found, all of them are included. A larger
number of candidates allows more points to be considered, but results in a
harder accuracy optimisation problem.

_Evaluation tolerances._ Let
$\mathcal{D}_{j}=\big{\\{}(p_{i}^{j},\tau_{i}^{j})\mid
i=1,\dots,s_{j}\big{\\}}$ be the set of training points at step $j$. By the
selection of candidate points, $s_{j}\geq s_{j-1}$ and $p_{i}^{j}=p_{i}^{j-1}$
for $i=1,\dots,s_{j-1}$ hold.

Optimal tolerances $\tau_{i}^{j}$ are given by the solution of (10) as a
function of the tolerances. In order to be able to solve the problem, we
ignore the shifts in the mean $\bar{y}$ as they cannot be predicted before
evaluating the model. Consequently, we only consider the impact of evaluation
tolerances on the predictive variance and, for evaluation tolerances
$\tau^{j}=(\tau_{1}^{j},\dots,\tau_{s_{j}}^{j})$, write $E(\tau^{j})$. As
already spent computational budget cannot be recovered by forgetting
previously acquired information, we impose the constraint
$\tau_{i}^{j}\leq\tau_{i}^{j-1}$ for $i=1,\dots,s_{j-1}$.

This results in the problem

$\min_{\tau^{j}\in\mathcal{T}_{j}}E(\tau^{j})\quad\text{subject to}\quad
W_{\tau^{j}|\mathcal{D}_{j-1}}\leq\Delta W_{j},$ (12)

where the set of admissible tolerances is

$\mathcal{T}_{j}=\\{(\tau_{1},\dots,\tau_{s_{j}})\in(\mathbb{R}^{+}\cup\\{+\infty\\})^{s_{j}}\mid\tau_{i}\leq\tau_{i}^{j-1}\text{
for }i\leq s_{j-1}\\}.$

If after optimization $\tau_{i}^{j}=+\infty$ holds for some $i>s_{j-1}$,
$p_{i}^{j}$ is excluded from the training set.

Before we can numerically solve the problem, we need to notice that
computational costs are not available before the evaluation is performed, such
that we need to resort to a priori work models. Following [13, 18], we make
use of established a priori asymptotic estimates for finite elements of degree
$r$ in space dimension $l$ and an optimal solver such as multigrid, and define

$W(\tau)=\tau^{-l/r}.$ (13)

This estimate is asymptotic for $\tau\rightarrow 0$. Consequently, despite
being inaccurate for low-accuracy evaluations, it is usually accurate for the
expensive high-accuracy ones.

Problem (12) is solved by multistart gradient descent with projection and
backtracking linesearch. The integral in $E$ is approximated again by Monte
Carlo integration on the samples $\mathcal{S}_{j}$, resulting in a numerical
objective

$E(\tau^{j})\approx\frac{1}{|\mathcal{S}_{j}|}\sum_{p\in\mathcal{S}_{j}}e_{\tau^{j}}(p).$

To implement gradient descent with projection, we adopt the coordinate change

$\tau^{j}=\left(\tau_{1},\dots,\tau_{s_{j}}\right)\mapsto\left(\tau_{1}^{-l/r},\dots,\tau_{s_{j}}^{-l/r}\right)=W^{j},$

such that the constraint in (12) becomes linear, transforming the set of
admissible tolerances $\mathcal{T}^{j}$ into a simplex and enabling efficient
projection.

##  6 Solution of the inverse problem

The previous sections established the inverse problem (4) and the sequential
approach (10) to surrogate model training. Similar to [17], we combine them to
an interleaved strategy given as pseudocode in Alg. 1.

Both the global error quantity (8) and the utility function (11) require
integration with respect to the posterior $\pi(p\mid y^{m})$. We perform the
integration through an MCMC sampling of the posterior, which is is at the same
time the ultimate goal of the inversion.

We start with an empty sample chain $\mathcal{S}_{0}=\emptyset$. At iteration
$j$, we draw a number $n_{j}$ of samples form $\pi(p\mid
y^{m},\mathcal{D}_{j-1})$, append them to $\mathcal{S}_{j-1}$, and remove the
oldest $h_{j}<n_{j}$ elements of the chain, as they have been drawn from a
less accurate posterior approximation. This results in the sample chain
$\mathcal{S}_{j}$, which is used to evaluate the integrals involved in the
training problem (10) at step $j$.

As the sample size $|\mathcal{S}_{j}|$ may be is too large for an efficient
evaluation of the integrals in (8) and (11), we use a sufficiently large
randomly extracted subset of $\mathcal{S}_{j}$ instead of the whole chain for
Monte Carlo integration.

When the computational budget is exhausted, the training of the surrogate
model terminates. A last round of samples is added to the chain, obtaining the
final set of samples from the posterior.

Algorithm 1 Surrogate-based Bayesian inversion

0: $\mathcal{D}_{0}$ initial design, $W$ budget

1: $\mathcal{S}_{0}\leftarrow\emptyset$

2: $W_{\mathcal{D}}\leftarrow 0$

3: $j\leftarrow 1$

4: while $W_{\mathcal{D}}\leq W$ do

5: decide: $n_{j}$ samples to draw, $h_{j}$ samples to remove

6: remove $h_{j}$ samples from $\mathcal{S}_{j-1}$

7: draw $n_{j}$ samples $\mathcal{S}$ from $\pi(p\mid y^{m},\mathcal{D})$

8: $\mathcal{S}_{j}\leftarrow\mathcal{S}_{j-1}\cup\mathcal{S}$

9: decide: $\Delta W_{j}$ iteration budget, $c_{j}$ number of candidates

10: obtain $c_{j}$ candidates

11: optimize accuracies $\tau^{j}$, update $\mathcal{D}$

12: evaluate forward model for decreased tolerances

13: $W_{\mathcal{D}}\leftarrow W_{\mathcal{D}}+\Delta W_{j}$

14: $j\leftarrow j+1$

15: end while

16: draw $n_{j}$ samples $\mathcal{S}$ from $\pi(p\mid y^{m},\mathcal{D})$

17: $\mathcal{S}_{j}\leftarrow\mathcal{S}_{j-1}\cup\mathcal{S}$

##  7 Numerical experiments

We present two illustrative experiments based on a Python implementation of
Alg. 1, where GPR is implemented with PyTorch. We adopt a separable kernel
with diagonal output structure and a Gaussian kernel as base [1]. The
hyperparameters are tuned by marginal likelihood maximisation using PyTorch’s
Adam optimiser, with the kernel’s correlation length scale constrained to
$[0,0.15]$.

As a benchmark, the results are compared with a non-adaptive space filling
approach, Latin Hypercube Sampling, and the position-adaptive-only training
strategy given by candidate point selection according to (11), i.e. all
candidates are accepted and evaluated with a fixed accuracy. For comparing the
approaches, the approximation errors (6) are computed numerically with MCMC
sampling utilising the true forward model. The implementation used for these
examples is available at
Zenodo111https://zenodo.org/doi/10.5281/zenodo.11066159.

###  7.1 1D analytical experiment

The first experiment is performed on a one-dimensional parameter space, with
$m=2$ measurements. We consider an analytical forward model
$y:\mathopen{]}0,1\mathclose{[}\to\mathbb{R}^{2}$ given by

$\displaystyle y(p)=\left[\frac{1}{2}p+\frac{1}{2}p^{2}\text{
exp}\left(\frac{1}{3}\text{sin}(12p-i)\right)\right]_{i=0,1}.$

This mimics the evaluation of a FE model on a 2D domain with quadratic
elements, i.e. $l/r=1$. The discretization error is simulated via a zero mean
Gaussian noise and the measurement likelihood is
$\Sigma_{l}=10^{-4}\mathrm{diag}(\frac{16}{9},\frac{4}{9})$.

A budget of 500 is considered: at each iteration two candidate points are
considered and a budget of 20 is assigned to each point. With the work model
(13), this results in a default tolerance of $0.05$ per point in the non-
adaptive strategies and a total of 12 iterations.

The number $n_{j}$ of new samples added into $\mathcal{S}_{j}$ is gradually
increased from 200 samples at the first iteration to 2000 in the last,
according to $n_{j}=200+1800\left(\frac{j}{12}\right)^{2}$. Similarly, the
number of discarded samples ranges from 200 to 1000, with $h_{1}=0$ as in the
first iteration the chain is empty, and
$h_{j}=200+800\left(\frac{j}{12}\right)^{2}$ for $j>1$.

Figure 2: Kullback-Leibler divergence of surrogated posterior and true
posterior for different training designs over the computational work spent in
the 1D example.

The obtained accuracies in terms of the Kullback-Leibler divergence between
true posterior $\pi(p\mid y^{m})$ and surrogated posterior $\pi(p\mid
y^{m},\mathcal{D})$ are shown in Fig. 2. Optimizing evaluation tolerances
provides a significant performance improvement over both other strategies.

###  7.2 2D analytical experiment

The second experiment considers a parameter space of two dimensions and $m=3$
measurements. The forward model
$y:\mathopen{]}-0.5,0.5\mathclose{[}^{2}\to\mathbb{R}^{3}$ is again
analytical, given by

$\displaystyle y(p)=\bigg{[}$ $\displaystyle\text{sin}(10k)(p_{1}-p_{2})\text{
exp}\left(\frac{1}{3}\text{sin}(8p_{2})\right)$
$\displaystyle+\text{cos}(10k)(p_{1}+p_{2})\text{
exp}\left(\frac{1}{3}\text{sin}(8p_{1})\right)\bigg{]}_{k\in\\{0,2,3\\}}.$

Figure 3: Reduction of surrogate standard deviation of the $y_{1}$, i.e.
$k=0$, component (left) and change of posterior distribution (right) between
iterations 7 and 9. The computational work for each point is represented by
its size. New points are added and some of the old points are refined. The
true parameter used for creating the artificial measurements $y^{m}$ is
denoted by a green star.

The underlying model is assumed to be a quadratic FE scheme on a 3D domain,
i.e. $l/r=1.5$. The discretization error is again simulated via zero mean
Gaussian noise and the measurement likelihood is
$\Sigma_{l}=10^{-4}\mathrm{diag}(1,1,4)$.

A working budget of 3600 is considered: at each iteration, 3 candidate points
are considered and a fixed budget of 100 corresponding to a fixed tolerance
$\tau=0.046$ is assigned to each point in the non-adaptive strategies for a
total of 12 iterations.

Figure 4: Kullback-Leibler divergence of surrogated posterior and true
posterior for different training designs over the computational work spent in
the 2D example.

The number $n_{j}$ of new samples added into $\mathcal{S}_{j}$ is gradually
increased from 200 samples at the first iteration to 4000 in the last,
according to $n_{j}=200+\lfloor 26.4j^{2}\rfloor$. Similarly, the number of
discarded samples ranges from 200 to 2000, with $h_{1}=0$ as in the first
iteration the chain is empty, and then $h_{j}=200+\lfloor 12.5j^{2}\rfloor$
for $j>1$. The error reduction by adding new points and decreasing tolerances
is illustrated in Fig. 3 for a single iteration. The performance in terms of
the Kullback-Leibler divergence between true and surrogated posteriors over
computational work is shown in Fig. 4. Again, a substantial performance
improvement is achieved by optimizing evaluation tolerances in addition to the
evaluation positions.

## Conclusions

When learning GPR surrogate models with numerically simulated training data as
a replacement for the true forward model in posterior sampling, significant
reductions of computational effort can be achieved with adaptive approaches.
With numerical forward models that allow exploiting accuracy-work trade-offs,
such as finite element simulations, the goal-oriented adaptive selection of
simulation tolerances appears to be particularly effective.

## References

  * [1] M.A. Álvarez, L. Rosasco, and N.D. Lawrence.  Kernels for vector-valued functions: A review.  Foundations and Trends in Machine Learning, 4(3):195–266, 2012\.
  * [2] T. Bai, A.L. Teckentrup, and K.C. Zygalakis.  Gaussian processes for Bayesian inverse problems associated with linear partial differential equations.  Technical report, arXiv:2307.08343, 2023.
  * [3] K. Crombecq, E. Laermans, and T. Dhaene.  Efficient space-filling and non-collapsing sequential design strategies for simulation-based modeling.  European Journal of Operational Research, 214:683–696, 2011.
  * [4] A. Giunta, S. Wojtkiewicz, and M. Eldred.  Overview of modern design of experiments methods for computational simulations (invited).  In 41st Aerospace Sciences Meeting and Exhibit, AIAA 2003-649, pages 1–17, 2003.
  * [5] P.J. Green, K. Łatuszyński, M. Pereyra, and C.P. Robert.  Bayesian computation: a summary of the current state, and samples backwards and forwards.  Stat. Comput., 25:835–862, 2015.
  * [6] V. Joseph and Y. Hung.  Orthogonal-maximin latin hypercube designs.  Statistica Sinica, 18:171–186, 2008.
  * [7] M. Järvenpää, M. U. Gutmann, A. Vehtari, and P. Marttine.  Parallel Gaussian process surrogate Bayesian inference with noisy likelihood evaluations.  Bayesian Analysis, 16, pp. 147–178., 2021.
  * [8] R. Lehmensiek, P. Meyer, and M. Müller.  Adaptive sampling applied to multivariate, multiple output rational interpolation models with application to microwave circuits.  International Journal of RF and Microwave Computer-Aided Engineering, 12(4):332–340, 2002.
  * [9] J. Nitzler, J. Biehler, N. Fehn, P.-S. Koutsourelakis, and A. Wall.  A generalized probabilistic learning approach for multi-fidelity uncertainty quantification in complex physical simulations.  Comp. Meth. Appl. Mech. Eng., 400:115600, 2022.
  * [10] N. Queipo, R. Haftka, W. Shyy, T. Goel, R. Vaidyanathan, and P. Tucker.  Surrogate-based analysis and optimization.  Progress in Aerospace Sciences, 41(1):1–28, 2005.
  * [11] C. Rasmussen and C.K.I. Williams.  Gaussian Processes for Machine Learning.  MIT Press, 2006.
  * [12] G. Sagnol, H.-C. Hege, and M. Weiser.  Using sparse kernels to design computer experiments with tunable precision.  In Proceedings of COMPSTAT 2016, pages 397–408, 2016.
  * [13] P. Semler and M. Weiser.  Adaptive Gaussian process regression for efficient building of surrogate models in inverse problems.  Inverse Problems, 39:125003, 2023.
  * [14] P. Semler and M. Weiser.  Adaptive gradient enhanced gaussian process surrogates for inverse problems.  In Proceedings of the MATH+ Thematic Einstein Semester 2023, 2024 (submitted).
  * [15] M. Sinsbeck and W. Nowak.  Sequential Design of Computer Experiments for the Solution of Bayesian Inverse Problems.  SIAM/ASA Journal on Uncertainty Quantification, 5:1, 640-664., 2017\.
  * [16] M. Sugiyama.  Active learning in approximately linear regression based on conditional expectation of generalization error.  Journal of Machine Learning Research, 7:141––166, 2006.
  * [17] Z. Wang and M. Broccardo.  A novel active learning-based Gaussian process metamodelling strategy for estimating the full probability distribution in forward UQ analysis.  Struct. Safety, 84:101937, 2020.
  * [18] M. Weiser and S. Ghosh.  Theoretically optimal inexact spectral deferred correction methods.  Commu. Appl. Math. Comp. Sci., 13(1):53–86, 2018.