remarkRemark hypothesisHypothesis claimClaim Certified Dimension Reduction for
Bayesian Updating with the Cross-Entropy MethodMax Ehre, Rafael Flock, Martin
Fußeder, Iason Papaioannou, Daniel Straub

# Certified Dimension Reduction for Bayesian Updating with the Cross-Entropy
Method††thanks: Submitted to the editors March 12, 2022. We acknowledge
support by the German Research Foundation (DFG) through Grants STR 1140/11-1
and PA 2901/1-1.

Max Ehre Technical University of Munich, School of Engineering and Design,
Engineering Risk Analysis Group<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Rafael Flock Technical University of
Denmark, Department of Applied Mathematics and Computer Science, DTU Compute
<EMAIL_ADDRESS>Martin Fußeder Technical University of Munich, School of
Engineering and Design, Chair of Structural Analysis ().
<EMAIL_ADDRESS>Iason Papaioannou22footnotemark: 2 Daniel
Straub22footnotemark: 2

###### Abstract

In inverse problems, the parameters of a model are estimated based on
observations of the model response. The Bayesian approach is powerful for
solving such problems; one formulates a prior distribution for the parameter
state that is updated with the observations to compute the posterior parameter
distribution. Solving for the posterior distribution can be challenging when,
e.g., prior and posterior significantly differ from one another and/or the
parameter space is high-dimensional. We use a sequence of importance sampling
measures that arise by tempering the likelihood to approach inverse problems
exhibiting a significant distance between prior and posterior. Each importance
sampling measure is identified by cross-entropy minimization as proposed in
the context of Bayesian inverse problems in Engel et al. (2021). To
efficiently address problems with high-dimensional parameter spaces we set up
the minimization procedure in a low-dimensional subspace of the original
parameter space. The principal idea is to analyse the spectrum of the second-
moment matrix of the gradient of the log-likelihood function to identify a
suitable subspace. Following Zahm et al. (2021), an upper bound on the
Kullback-Leibler-divergence between full-dimensional and subspace posterior is
provided, which can be utilized to determine the effective dimension of the
inverse problem corresponding to a prescribed approximation error bound. We
suggest heuristic criteria for optimally selecting the number of model and
model gradient evaluations in each iteration of the importance sampling
sequence. We investigate the performance of this approach using examples from
engineering mechanics set in various parameter space dimensions.

###### keywords:

Bayesian inverse problems, high dimensions, cross-entropy method, importance
sampling, certified dimension reduction

62F15, 62L12, 62P30, 60G60, 65C05

##  1 Introduction

We consider inverse problems in the context of a computational model $f$ with
$\bm{y}=f(\bm{\theta})$. That is, we want to characterise a cause (parameters
of the computational model, $\bm{\theta}$) based on observations of the
corresponding effects or consequences of said cause (output of the
computational model $\bm{y}$). An example is a structural system represented
with a finite element model that is parameterized by loads, geometric and
material properties $\bm{\theta}\in\mathcal{X}\subseteq\mathbb{R}^{d}$ and
that produces outputs such as stresses and deflections
$\bm{y}\in\mathcal{Y}\subseteq\mathbb{R}^{m}$. In the majority of
applications, we cannot expect the inverse problem to be well defined, i.e.,
there need not be a solution, the solution may not be unique or it might be
very sensitive to the observations [54]. To further complicate the matter, in
practice, observations are often incomplete and/or contaminated with noise.
Here, we focus on the Bayesian approach to inverse problems, which offers a
consistent framework for incorporating both noisy and incomplete observations
as well as addressing ill-posedness by regularizing the problem using prior
information [29, 54]. [29, 54] discuss the Bayesian inverse problem (BIP) in
infinite-dimensional settings, while in practice, we usually retreat to the
finite-dimensional case by means of discretizing infinite-dimensional random
objects such as random fields and processes. Hence, we focus on finite-
dimensional BIPs in this work.

We represent $\bm{\theta}$ and $\bm{y}$ as real-valued random vectors
$\bm{\Theta}:\mathcal{X},\mathcal{B}(\mathcal{X})\to\mathbb{R}$ and
$\bm{Y}:\mathcal{Y},\mathcal{B}(\mathcal{Y})\to\mathbb{R}$, where
$\mathcal{B}(\cdot)$ is the Borel $\sigma$-algebra, and we assume the
probability measures $\mathbb{P}_{\bm{\Theta}}$, $\mathbb{P}_{\bm{Y}}$ to be
absolutely continuous with respect to the respective Lebesgue measures on
$\mathcal{B}(\mathcal{X})$ and $\mathcal{B}(\mathcal{Y})$. We then may use the
associated probability density functions (PDF) $p(\bm{\theta})$ and
$p(\bm{y})$ to characterize $\bm{\Theta}$ and $\bm{Y}$.

We start by placing a prior distribution on $\bm{\theta}$ by defining the
prior PDF $p_{0}(\bm{\theta}):\mathbb{R}^{d}\rightarrow\mathbb{R}_{\geq 0}$.
As the name suggests, $p_{0}(\bm{\theta})$ formalizes any information one may
have on $\bm{\theta}$ _prior_ to considering any observations. This
information may come as the outcome of an expert elicitation, selection rules
[41] and/or guiding principles to construct noninformative priors such as
Jeffreys’s priors [26] or priors satisfying the maximum entropy principle
[25]. While many of these principles rest on the idea to minimize the
influence the prior exerts on the posterior distribution and thus aim at
‘letting the data speak’, it is flat/weak priors in particular that can lead
to overly confident inference results [18]. A single layer of priors may not
do justice to complex models with a large number of unobserved variables, in
which case _hierarchical models_ with several layers of prior distributions
can be utilized [17, Section 2.8].

Next, one or several observations of $\bm{y}$ that we refer to as
$\tilde{\bm{y}}$, are represented by the _likelihood_
$L(\bm{\theta})\vcentcolon=p(\tilde{\bm{y}}|\bm{\theta}):\mathbb{R}^{d}\rightarrow\mathbb{R}_{\geq
0}$, which states how likely these observations are to occur under any given
set of parameters $\bm{\theta}$. In Bayesian inverse problems, the likelihood
will be a function of $f$, thereby facilitating the backpropagation of
information on outputs of $f$, $\bm{y}$, to its parameters $\bm{\theta}$. With
this, the posterior PDF of $\bm{\theta}$ conditional on observations
$\tilde{\bm{y}}$ follows from Bayes’ theorem as

(1)
$p_{\bm{y}}(\bm{\theta})\vcentcolon=p(\bm{\theta}|\tilde{\bm{y}})=\frac{p(\tilde{\bm{y}}|\bm{\theta})p_{0}(\bm{\theta})}{p(\tilde{\bm{y}})}=\frac{L(\bm{\theta})p_{0}(\bm{\theta})}{Z},$

where

(2)
$Z=\int_{\mathcal{X}}L(\bm{\theta})p_{0}(\bm{\theta})\mathrm{d}\bm{\theta}$

is the marginal likelihood of the data also known as the _model evidence_. We
assume the evidence is finite $Z<\infty$ and the likelihood is Borel-
measurable. In the Bayesian approach, solving the inverse problem amounts to
computing the posterior distribution of $\bm{\theta}$ and generating samples
from $p_{\bm{y}}(\bm{\theta})$. In many instances, the posterior distribution
cannot be computed exactly. Instead, sampling approaches such as importance
sampling (IS) [19] or Markov Chain Monte Carlo (MCMC) [22] are used to sample
from the posterior and construct estimates of posterior expectations. If prior
and posterior distributions differ from one another significantly,
constructing efficient biasing or proposal densities for IS or MCMC,
respectively, becomes difficult. Such problems can be addressed by repeatedly
applying sampling methods on an artificial sequence of distributions that
gradually approach the posterior starting from the prior, namely _sequential
Monte Carlo_ methods (SMC) [13, 37, 11]. In _adaptive_ SMC [30, 24, 32], the
distribution sequence is determined during runtime based on intermediate
samples.

In SMC approaches, the distributions appearing in the sequence are
characterized by samples that are obtained through resample-move steps;
samples from each previous distribution are moved via MCMC sampling to obtain
samples from the next distribution. However, MCMC produces dependent samples.
Alternatively, in _cross-entropy importance sampling_ (CE-IS) [50], a sequence
of parametrized distributions is defined such that each new distribution in
the sequence is identified through solving an optimization (cross-entropy
minimization) problem. Estimation of the target distribution is then performed
with IS using the final fitted parametric density. Hence, CE-IS avoids MCMC
sampling and dependent samples. CE-IS was introduced in the context of rare
event estimation in [48] and was recently applied to solve the BIP in [15].

Both acceptance rate and mixing time — and hence, computational efficiency —
of many MCMC algorithms deteriorate as the problem dimension $d$ increases
[46, 36]; notable exceptions include the _preconditioned Crank-Nicholson_
(pCN) sampler [3, 8] and Hamiltonian MCMC [38]. Therefore, different
approaches have been proposed to reduce the dimension of the inverse problem
by identifying low-dimensional subspaces on which the solution to the original
problem may be identified in good approximation. While their existence cannot
be guaranteed independent of the inverse problem, low-dimensional subspaces
frequently occur in BIPs as a result of $f$ being a smoothing operator applied
to the input vector $\bm{\theta}$, e.g., in the form of solutions to a set of
partial differential equations (PDEs). In [35, 57] the problem dimension is
reduced by representing the prior with a truncated Karhunen-Loéve-expansion.
In the context of linear BIPs, [16, 51] construct low-rank approximations of
the prior-preconditioned Hessian of the log-likelihood thereby exploiting
structure in both prior and likelihood. The _likelihood-informed subspace_
method of [9] extends this approach to nonlinear BIPs based on a low-rank
approximation of the posterior-preconditioned Hessian of the log-likelihood.
[7] propose a similar approach in which they identify an _active subspace_ of
the BIP, i.e., a low-rank approximation of the prior-preconditioned negative
log-likelihood gradient. Building on the idea of likelihood-informed
subspaces, [60] proposes _certified dimension reduction_ for nonlinear BIPs
and derives an upper bound on the _Kullback-Leibler-divergence_ between
reduced and full space posterior in function of the subspace dimension.

While the CE-IS approach to BIPs of [15] circumvents MCMC altogether, its
performance deteriorates with increasing parameter dimension. This is both due
to an increasing degeneracy of the IS weights that are used in the context of
CE-IS [49] as well as the rapidly growing number of parameters in the employed
distribution models. For example, in Gaussian models with full covariance
structure, the number of parameters is $p=d(d+3)/2$, implying that the number
of $f$-evaluations required to obtain an accurate fit scale quadratically with
$d$. Solving BIPs with CE-IS is therefore only suitable for low parameter
dimension. [58] uses CE-IS for estimating rare event probabilities of models
with large parameter dimension by applying certified dimension reduction.

In this work, we devise a scheme to efficiently solve nonlinear BIPs using CE-
IS and certified dimension reduction. Our method extends the approaches of
[15, 58] to address high-dimensional BIPs. Moreover, we introduce heuristic
rules for adaptively selecting the number of model and model gradient
evaluations during the simulation. Section 2 recapitulates CE-IS for BIPs
following [15]. Section 3 details the certified dimensionality reduction
approach for CE-IS and Subsection 3.4 contains a discussion on methodology,
algorithmic details and a summary of the final procedure. In Section 4, we
investigate the efficacy of our method on two structural engineering examples
both featuring large parameter dimensions. Concluding remarks are given in
Section 5.

##  2 Cross-entropy-based importance sampling for Bayesian updating

###  2.1 Importance sampling and the cross-entropy method

In this chapter, we briefly describe the _CE-based IS method for Bayesian
updating_ (CEBU) proposed in [15]. Importance sampling is a variance reduction
method for estimating expectations of a function $G(\bm{\theta})$,
$\mathbb{E}_{p}[G(\bm{\Theta})]$ [50, Chapter 5]. Note that we use lowercase
letters for deterministic variables. We use uppercase letters for matrices and
random variables/vectors with the exception of random samples, which are
denoted with lowercase letters yet treated as random variables. Throughout
this work we assume all random vectors to be real-valued, i.e.,
$\bm{X}:(\mathcal{X},\mathcal{B}(\mathcal{X}))\to\mathbb{R}$, where
$(\mathcal{X},\mathcal{B}(\mathcal{X}))$ is a measurable space consisting of
the outcome space $\mathcal{X}$ and its associated Borel $\sigma$-algebra
$\mathcal{B}(\cdot)$. Further, we assume probability measures
$\mathbb{P}_{\bm{X}}$ to be absolutely continuous with respect to the
respective Lebesgue measures on $\mathcal{B}(\mathcal{X})$ so that we may use
the associated PDFs $p(\bm{x})$ to describe $\bm{X}$. Let $q(\bm{\theta})$ be
a PDF on $\mathcal{X}$ such that $q(\bm{\theta})>0$ whenever
$p(\bm{\theta})>0$ and suppose we only know
$\psi(\bm{\theta})=cp(\bm{\theta})$ pointwise with unknown normalizing
constant $c$. Then we can write

(3)
$\mu\vcentcolon=\mathbb{E}_{p}[G(\bm{\Theta})]=\frac{1}{c}\mathbb{E}_{q}\left[\frac{G(\bm{\Theta})\psi(\bm{\Theta})}{q(\bm{\Theta})}\right]=\frac{\mathbb{E}_{q}\left[G(\bm{\Theta})w(\bm{\Theta})\right]}{\mathbb{E}_{q}\left[w(\bm{\Theta})\right]},~{}~{}~{}\cite[cite]{[\@@bibref{}{Owen2013}{}{},
\mathrm{Chapter~{}9}]}$

where $q$ is termed the importance, auxiliary, instrumental or biasing density
and $w(\bm{\theta})=\psi(\bm{\theta})/q(\bm{\theta})$ is referred to as the
likelihood ratio or IS weight. LABEL:e:IS leads to the _self-normalized IS
estimate_

(4)
$\widehat{\mu}_{\mathrm{IS},q}=\frac{1}{n\widehat{c}}\sum_{k=1}^{n}G(\bm{\theta}_{k})w(\bm{\theta}_{k}),~{}~{}~{}\bm{\theta}_{k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}q(\bm{\theta}),$

where an estimate of the normalizing constant is given as
$\widehat{c}=n^{-1}\sum_{k=1}^{n}w(\bm{\theta}_{k})$. For many problems $q$
can be chosen such that Eq. 4 has lower variance
$\widehat{\mathbb{V}}[\widehat{\mu}_{\mathrm{IS},q}]$ than the crude Monte
Carlo estimate [40, Chapter 9].

In the context of BIPs $p$ is a posterior distribution $p_{\bm{y}}$ and the
normalizing constant $c$ in Eq. 3 is the model evidence $Z$. $p_{\bm{y}}$ is
the optimal IS density to estimate the model evidence as
$\mathbb{V}[\widehat{Z}]=0$ if $q=p_{\bm{y}}$. Since sampling from the
posterior is usually difficult, we continue with a discussion of how to get a
parametric $q$ close to $p_{\bm{y}}$.

[48] proposed finding a parametric IS density $q(\bm{\theta},\bm{v})$ with
parameters $\bm{v}\in\mathcal{V}$ by minimizing the _Kullback-Leibler
divergence_ (KLD) between $q(\bm{\theta},\bm{v})$ and an optimal IS density in
the context of rare event probability estimation. [15] builds on this
principle to estimate a parametric distribution that is close to the posterior
$p_{\bm{y}}$ as follows. The KLD between the posterior and the parametric
density $D_{\mathrm{KL}}(p_{\bm{y}}(\bm{\theta})||q(\bm{\theta},\bm{v}))$ is
defined as [50]

(5)
$\begin{split}D_{\mathrm{KL}}(p_{\bm{y}}(\bm{\theta})||q(\bm{\theta},\bm{v}))&=\mathbb{E}_{p_{\bm{y}}}\left[\ln\left(\frac{p_{\bm{y}}(\bm{\Theta})}{q(\bm{\Theta},\bm{v})}\right)\right]\\\
&=\frac{1}{Z}\mathbb{E}_{p_{0}}[L(\bm{\Theta})\ln(p_{\bm{y}}(\bm{\Theta}))]\underbrace{-\frac{1}{Z}\mathbb{E}_{p_{0}}[L(\bm{\Theta})\ln(q(\bm{\Theta},\bm{v}))]}_{\mathrm{cross~{}entropy}~{}H(p_{\bm{y}},q(\cdot,\bm{v}))}.\end{split}$

The first summand on the right-hand side of Eq. 5 is not a function of
$\bm{v}$ so that minimizing
$D_{\mathrm{KL}}(p_{\bm{y}}(\bm{\theta})||q(\bm{\theta},\bm{v}))$ is
equivalent to maximizing the negative cross entropy:

(6)
$\bm{v}=\operatorname*{arg\,max}\limits_{\bm{v}\in\mathcal{V}}\mathbb{E}_{p_{0}}[L(\bm{\Theta})\ln(q(\bm{\Theta},\bm{v}))],$

which conveniently does not depend on the unknown $Z$. An approximate solution
of this optimization problem based on samples from $p$ reads

(7)
$\widehat{\bm{v}}=\operatorname*{arg\,max}\limits_{\bm{v}\in\mathcal{V}}\frac{1}{n}\sum\limits_{k=1}^{n}L(\bm{\theta}_{k})\ln(q(\bm{\theta}_{k},\bm{v})),~{}~{}~{}\bm{\theta}_{k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}p_{0}(\bm{\theta}).$

The optimization problem in Eq. 7 is usually convex, continuous and the
objective function is differentiable with respect to $\bm{v}$ such that
identifying $\widehat{\bm{v}}$ is straight-forward. Closed-form solutions of
Eq. 7 exist in various situations, e.g., if $q(\bm{\theta},\bm{v})$ is any
member of the exponential family [50, Chapter 8]. [31, 20] use a Gaussian
mixture model in order to capture several disconnected failure regions. [59]
and [43] use _von Mises-Fisher_ and _von Mises-Fisher-Nakagami_ (vMFN) mixture
models, respectively, to overcome the poor performance of Gaussian models in
high-dimensional rare event probability estimation problems. [15] test the
performance of both Gaussian mixture and vMFN mixture models in the context of
the CE method for BIPs and show that although the latter has superior
performance in certain high-dimensional settings, the former possesses higher
flexibility and is thus able to accurately describe complicated posteriors. In
all these works, different variants of expectation maximization are used to
solve for $\widehat{\bm{v}}$.

How well $q(\bm{\theta},\widehat{\bm{v}})$ approximates
$p_{\bm{y}}(\bm{\theta})$ hinges on how well samples from $p$ can inform the
objective function about $L(\bm{\theta})$. In other words, if prior and
likelihood are not close to one another, we cannot expect the solution of Eq.
7 to yield a satisfying approximation to $p_{\bm{y}}(\bm{\theta})$ independent
of the parametric model choice. This problem can be overcome by tempering the
likelihood as described in the following section.

###  2.2 Tempering the likelihood

In order to bridge the distance between prior and likelihood one may break
down the single CE problem into several smaller ones. To this end, we define a
sequence of PDFs $\\{q_{t}(\bm{\theta})\\}_{j=1}^{m}$ with

(8)
$q_{t}(\bm{\theta})\vcentcolon=\frac{L^{\beta_{t}}(\bm{\theta})p_{0}(\bm{\theta})}{Z_{t}},$

where
$Z_{t}=\int_{\mathcal{X}}L^{\beta_{t}}(\bm{\theta})p_{0}(\bm{\theta})\mathrm{d}\bm{\theta}$
and ensuring
$0=\vcentcolon\beta_{0}<\beta_{1}<\dots<\beta_{m-1}<\beta_{m}\vcentcolon=1$
such that $q_{0}(\bm{\theta})\vcentcolon=p_{0}(\bm{\theta})$ and
$q_{m}(\bm{\theta})\vcentcolon=p_{\bm{y}}(\bm{\theta})$. The idea is to start
with samples from $p_{0}(\bm{\theta})$ and select $\beta_{1}$ small enough to
facilitate an accurate estimate $\widehat{\bm{v}}_{1}$. Next, upon selecting
$\beta_{2}\in(\beta_{1},1]$, samples from
$q(\bm{\theta},\widehat{\bm{v}}_{1})$ can be used to estimate
$\widehat{\bm{v}}_{2}$. This procedure is repeated until $\beta_{m}=1$ after
$m$ steps and the CE problem is solved for the target posterior density. The
$t$-th CE minimization problem reads

(9)
$\widehat{\bm{v}}_{t}=\operatorname*{arg\,max}\limits_{\bm{v}\in\mathcal{V}}\frac{1}{n}\sum\limits_{k=1}^{n}\ln(q(\bm{\theta}_{k},\bm{v}))w_{t}(\bm{\theta}_{k}),~{}~{}~{}\bm{\theta}_{k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}q(\bm{\theta},\widehat{\bm{v}}_{t-1}),$

with
$w_{t}(\bm{\theta})=L^{\beta_{t}}(\bm{\theta})p_{0}(\bm{\theta})/q(\bm{\theta},\widehat{\bm{v}}_{t-1})$.

In Eq. 9, the likelihood ratio or weight $w_{t}(\bm{\theta})$ accounts for the
fact that the $t$-th PDF parameter estimate $\widehat{\bm{v}}_{t}$ is based on
samples from the $(t-1)^{\mathrm{th}}$ PDF
$q(\bm{\theta},\widehat{\bm{v}}_{t-1})$. The variance of
$\widehat{\bm{v}}_{t}$ depends on the variance of the weights
$w_{t}(\bm{\theta})$. In particular, if the numerator PDF of $w$ has fatter
tails than its denominator PDF, the weight variance blows up and the parameter
estimate $\widehat{\bm{v}}_{t}$ deteriorates. The normalized _effective sample
size_ (nESS) is a common performance metric of IS that is directly related to
the variance of the weights [40, Chapter 9]:

(10)
$n_{\mathrm{eff}}=\frac{1}{1+\delta_{w}^{2}}~{}~{}~{}\mathrm{with}~{}~{}~{}\delta_{w}=\frac{\sqrt{\mathbb{V}[w(\bm{\Theta})]}}{\mathbb{E}[w(\bm{\Theta})]}$

the coefficient of variation of the weights. Therefore in [15], $\beta_{t}$ is
computed adaptively in each step such as to achieve a target nESS
$n_{\mathrm{eff}}^{*}$ by utilizing a sample-based estimate of the coefficient
of variation of the weights $\widehat{\delta}_{w}$:

(11) $\beta_{t}=\operatorname*{arg\,min}\limits_{\beta_{j-1}<\beta\leq
1}\left(n_{\mathrm{eff}}^{*}-\frac{1}{1+\widehat{\delta}_{w}(\beta)^{2}}\right)^{2}=\operatorname*{arg\,min}\limits_{\beta_{j-1}<\beta\leq
1}\left(n_{\mathrm{eff}}^{*}-\frac{(\sum_{k=1}^{n}w(\bm{\theta}_{k};\beta))^{2}}{\sum_{k=1}^{n}w^{2}(\bm{\theta}_{k};\beta)}\right)^{2}.$

The weights on the right-hand side of Eq. 11 can be evaluated approximately by
assuming $q(\bm{\theta},\widehat{\bm{v}}_{t-1})=q_{t-1}(\bm{\theta})$ in each
step, such that $w(\bm{\theta})\propto L(\bm{\theta})^{\beta_{t}-\beta_{t-1}}$
(the factor $Z_{t-1}$ required here for equality cancels out in Eq. 11 and is
immaterial to its solution).

###  2.3 Method

In [15], CEBU is implemented in the $d$-dimensional standard-normal space
$(\mathcal{U},\mathcal{B}(\mathcal{U}),\mathbb{P}_{\bm{U}})$ with
$\mathcal{U}=\mathbb{R}^{d}$, so that the standard-normal random vector
$\bm{U}\sim\varphi_{d}(\bm{u})$, where $\varphi_{d}$ denotes the
$d$-dimensional standard-normal PDF. Under a suitable isoprobabilistic
transformation $T:\bm{\Theta}\rightarrow\bm{U}$, e.g., using the inverse CDF
transform, the Rosenblatt transform [47] or copula models [34, 39, 55],
arbitrary priors $p_{0}(\bm{\theta})$ are transformed to $\varphi_{d}(\bm{u})$
while the likelihood in standard-normal space is given as
$\tilde{L}\vcentcolon=(L\circ T^{-1})(\bm{u})$. Similarly, we define
$\tilde{p}_{\bm{y}}\vcentcolon=\tilde{L}(u)\varphi(\bm{u})/Z$ where the
evidence $Z$ is invariant under the transformation $T$ [15, Appendix A].

The CEBU loop terminates once $\beta_{t}=1$. Then, a set of $n$ samples is
drawn from the final parametric density corresponding to $\beta_{t}=1$. These
samples are subsequently reweighted to generate samples from the true
posterior distribution $\tilde{p}_{\bm{y}}$. To this end, a final set of
weights
$\bm{\mathrm{w}}_{\mathrm{final}}(\bm{u})=\tilde{L}(\bm{u})\varphi_{d}(\bm{u})/q(\bm{u},\widehat{\bm{v}}_{t})$
is computed as the likelihood ratio of the unnormalized posterior in standard-
normal space $\tilde{L}(\bm{u})\varphi_{d}(\bm{u})$ and the parametric density
corresponding to $\beta_{t}=1$, $q(\bm{u},\widehat{\bm{v}}_{t})$. The evidence
can be written as

(12)
$Z=\mathbb{E}_{\varphi_{d}}[\tilde{L}(\bm{U})]=\mathbb{E}_{q(\bm{u},\widehat{\bm{v}}_{t})}[\tilde{L}(\bm{U})\varphi_{d}(\bm{U})/q(\bm{U},\widehat{\bm{v}}_{t})]=\mathbb{E}_{q(\bm{u},\widehat{\bm{v}}_{t})}[\bm{\mathrm{w}}_{\mathrm{final}}],$

which suggests estimating $Z$ as

(13)
$\widehat{Z}=\frac{1}{n}\sum_{k=1}^{n}w_{\mathrm{final}}(\bm{u}_{k}),~{}~{}~{}\bm{u}_{k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}q(\bm{u},\widehat{\bm{v}}_{t}).$

A desired number of $N$ weighted posterior samples
$\\{\bm{u}_{k}\\}_{k=1}^{N}$ may then be obtained by resampling the last set
of $n$ samples corresponding to $\beta_{t}=1$ with replacement and weighted
with the normalized final weights
$\\{\bm{\mathrm{w}}_{\mathrm{final}}(\bm{u}_{k})/(n\widehat{Z})\\}_{k=1}^{n}$.
In [15], the authors use a stratified version of this resampling step based on
[14]. In a final step, these samples are transformed back to
$\bm{\Theta}$-space through applying the inverse transform $T^{-1}$. The
entire procedure is summarized in Alg. 1.

Algorithm 1 CE-BU

Input Likelihood $L$, transform $T$, target nESS $n_{\mathrm{eff}}^{*}$, #
post. samples $N$, # samples/level $n$
Output posterior samples $\mathbf{\Theta}_{\mathrm{post}}$, estimated evidence
$\widehat{Z}$

1: Set $t\leftarrow 0$, $\beta_{0}\leftarrow 0$, $\widehat{\bm{v}}_{0}$ (so
that $q_{0}(\bm{u})=\varphi_{d}(\bm{u})$)

2: while $\beta_{t}<1$ do

3: $t\leftarrow t+1$

4: Sample $\bm{\mathrm{U}}\in\mathbb{R}^{n\times
d}\leftarrow\\{\bm{u}_{k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}q(\bm{u},\widehat{\bm{v}}_{t-1})\\}_{k=1}^{n}$

5: Evaluate $\bm{\mathrm{\ell}}\in\mathbb{R}^{n\times 1}\leftarrow(L\circ
T^{-1})(\bm{\mathrm{U}})$

6: Compute $\beta_{t}$ with Eq. 11

7: Evaluate $\bm{\mathrm{w}}\in\mathbb{R}^{n\times
1}\leftarrow\bm{\mathrm{\ell}}^{\beta_{t}-\beta_{t-1}}$

8: Compute $\widehat{\bm{v}}_{t}$ with Eq. 9

9: end while

10: Sample $\bm{\mathrm{U}}\in\mathbb{R}^{n\times
d}\leftarrow\\{\bm{u}_{k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}q(\bm{u},\widehat{\bm{v}}_{t})\\}_{k=1}^{n}$

11: Evaluate $\bm{\mathrm{w}}_{\mathrm{final}}\in\mathbb{R}^{n\times
1}\leftarrow\left\\{\frac{(L\circ
T^{-1})(\bm{u}_{k})\varphi_{d}(\bm{u}_{k})}{q(\bm{u}_{k},\widehat{\bm{v}}_{t})}\right\\}_{k=1}^{n}$

12: Estimate evidence
$\widehat{Z}\leftarrow\frac{1}{n}\sum_{k=1}^{n}\mathrm{w}_{\mathrm{final},k}$

13: Normalize weights
$\bar{\bm{\mathrm{w}}}_{\mathrm{final}}\leftarrow\bm{\mathrm{w}}_{\mathrm{final}}/(n\widehat{Z})$

14: $\mathbf{U}_{\mathrm{post}}\leftarrow$ Resample (with replacement) $N$
times from $\mathbf{U}$ with weighting
$\bar{\bm{\mathrm{w}}}_{\mathrm{final}}$

15: $\mathbf{\Theta}_{\mathrm{post}}\leftarrow
T^{-1}(\mathbf{U}_{\mathrm{post}})$

16: return $\mathbf{\Theta}_{\mathrm{post}}$, $\widehat{Z}$

##  3 CEBUred: Certified dimension reduction for CEBU

###  3.1 Linear subspaces of $\mathcal{U}$

[15] test both Gaussian and vMFN mixture models for $q$. A $K$-component
Gaussian mixture requires fitting $Kd(d+3)/2+K-1$ parameters whereas a
$K$-component vMFN mixture features only $K(d+3)+K-1$ parameters. In spite of
the more advantageous linear scaling in $d$ offered by vMFN mixtures, the
required number of samples per CE-level can quickly exceed the computational
budget if $d$ is large.

In [15], it is shown that the CE method with the Gaussian mixture model is
able to obtain accurate representations of posterior densities in various
problem settings. However, the Gaussian mixture model performs poorly in high-
dimensional IS [20]. This is due to the fact that IS weights with respect to
Gaussian densities tend to degenerate in high dimensions. Further, the number
of parameters of the GM model increases quadratically with the input
dimension. The latter implies that the required number of samples per CE-level
to obtain accurate parameter estimates becomes prohibitively large in high
dimensions. To alleviate these problems, we draw on the ideas presented in
[60] to determine a low-dimensional linear subspace of $\mathcal{X}$ in which
an effective IS density can be constructed. The resulting approach can be
viewed as an extension of the CE method with failure-informed dimension for
rare event estimation, proposed in [58].

In each step of CEBU, for the tempered posterior distribution
$\tilde{p}_{\bm{y},\beta}(\bm{u})={Z}^{-1}_{\beta}\tilde{L}^{\beta}(\bm{u})\varphi_{d}(\bm{u})$
with
${Z}_{\beta}=\int_{\mathbb{R}^{d}}\tilde{L}^{\beta}(\bm{u})\varphi_{d}(\bm{u})\mathrm{d}\bm{u}$
we seek an approximation of the form

(14)
$\tilde{p}_{\bm{y},\beta}^{(r)}\propto(g\circ\mathbf{P}_{r})(\bm{u})\varphi_{d}(\bm{u}),$

where $g:\mathbb{R}^{d}\rightarrow\mathbb{R}_{>0}$ is a Borel-measurable
function referred to as _profile function_ in the following.
$\mathbf{P}_{r}\in\mathbb{R}^{d\times d}$ is a rank-$r$ projection matrix,
i.e., $\mathbf{P}_{r}\circ\mathbf{P}_{r}=\mathbf{P}_{r}$. Any
$\bm{u}\in\mathbb{R}^{d}$ can be decomposed as
$\bm{u}=\mathbf{P}_{r}\bm{u}+\mathbf{P}_{\perp}\bm{u}=\bm{u}_{r}+\bm{u}_{\perp}$
with the complementary projection
$\mathbf{P}_{\perp}\vcentcolon=\mathbf{I}_{d}-\mathbf{P}_{r}$ satisfying
$\mathrm{Im}(\mathbf{P}_{\perp})=\mathrm{Ker}(\mathbf{P}_{r})$. We call
$\mathcal{U}_{r}\vcentcolon=\mathrm{Im}(\mathbf{P}_{r})$ the _likelihood-
informed subspace_ (LIS) and
$\mathcal{U}_{\perp}\vcentcolon=\mathrm{Im}(\mathbf{P}_{\perp})$ the
_complementary subspace_ (CS). The LIS and CS are at this point still subsets
of the ambient space $\mathbb{R}^{d}$ so that no effective dimension reduction
is achieved by their introduction. However, they correspond to lower-
dimensional spaces we refer to as local LIS $\bar{\mathcal{U}}_{r}$ and local
CS $\bar{\mathcal{U}}_{\perp}$, where in standard-normal space
$\bar{\mathcal{U}}_{r}=\mathbb{R}^{r}$ and
$\bar{\mathcal{U}}_{\perp}=\mathbb{R}^{d-r}$. We discuss the mapping to these
local subspaces in more detail in Subsection 3.4. The profile function $g$ is
only a function of $\bm{u}_{r}\in\mathcal{U}_{r}$ and is constant in
$\bm{u}_{\perp}\in\mathcal{U}_{\perp}$. Following [58], we first define an
optimal $g$ for the tempered posterior distributions of CEBU given a
projection $\mathbf{P}_{r}$ in Subsection 3.2. Next, we identify the
projection that minimizes the KLD between full and low-rank posterior in
Subsection 3.3 and lay out the certified dimensionality reduction for CEBU in
Subsection 3.4.

###  3.2 Optimal profile function $g$

[60] show that for a given projection matrix $\mathbf{P}_{r}$, the optimal
profile function $g_{\beta}^{\star}(\bm{u})$ that minimizes the KLD
$D_{\mathrm{KL}}(\tilde{p}_{\bm{y},\beta}||\tilde{p}_{\bm{y},\beta}^{(r)})$,
is the following conditional expectation

(15)
$\mathbb{E}_{p}[\tilde{L}^{\beta}(\bm{U})|\mathbf{P}_{r}\bm{u}]:\bm{u}\rightarrow\int_{\mathbb{R}^{d-r}}\tilde{L}^{\beta}(\mathbf{P}_{r}\bm{u}+\bm{\Phi}_{\perp}\bar{\bm{u}}_{\perp})p_{\perp}(\bar{\bm{u}}_{\perp}|\mathbf{P}_{r}\bm{u})\mathrm{d}\bar{\bm{u}}_{\perp},$

where $\bm{\Phi}_{\perp}\in\mathbb{R}^{d\times d-r}$ such that
$\mathrm{span}(\bm{\Phi}_{\perp})=\mathrm{Ker}(\mathbf{P}_{r})$ and
$\bar{\bm{u}}_{\perp}\in\mathbb{R}^{d-r}$. The conditional PDF
$p_{\perp}(\bar{\bm{u}}_{\perp}|\mathbf{P}_{r}\bm{u})$ reads

(16)
$p_{\perp}(\bar{\bm{u}}_{\perp}|\mathbf{P}_{r}\bm{u})=\frac{\varphi_{d}(\mathbf{P}_{r}\bm{u}+\bm{\Phi}_{\perp}\bar{\bm{u}}_{\perp})}{\int_{\mathbb{R}^{d-r}}\varphi_{d}(\mathbf{P}_{r}\bm{u}+\bm{\Phi}_{\perp}\bar{\bm{u}}_{\perp}^{\prime})\mathrm{d}\bm{u}_{\perp}^{\prime}},$

which, by convention, equals zero whenever the denominator of Eq. 16 equals
zero. Following from the optimality of Eq. 15, the optimal reduced posterior
reconstruction in standard-normal space reads

(17)
$\tilde{p}_{\bm{y},\beta}^{(r,\star)}\propto\mathbb{E}_{p}[\tilde{L}^{\beta}(\bm{U})|\mathbf{P}_{r}\bm{u}]\varphi_{d}(\bm{u}).$

[60] also remarks that the conditional expectation Eq. 15 is not only optimal
with respect to the KL divergence but also minimizes the mean-square
reconstruction error of the likelihood function with respect to the prior
measure
$\mathbb{E}_{p}[(\tilde{L}^{\beta}(\bm{U})-(g\circ\mathbf{P}_{r})(\bm{U}))^{2}]$.

###  3.3 Optimal projection $\mathbf{P}_{r}$

Under assumptions on the prior distribution that hold in the standard-normal
setting [60, Example 2.6], the _subspace logarithmic Sobolev inequality_ in
[60, Theorem 2.9] states that $\int_{\mathbb{R}^{d}}\lVert\nabla
h(\bm{u})\lVert^{2}\varphi_{d}(\bm{u})\mathrm{d}\bm{u}\leq\infty$ for any
continuously differentiable function
$h:\mathbb{R}^{d}\rightarrow\mathbb{R^{d}}$ and for any projection
$\mathbf{P}_{r}\in\mathbb{R}^{d\times d}$,

(18)
$\int_{\mathbb{R}^{d}}h^{2}(\bm{u})\ln\left(\frac{h^{2}(\bm{u})}{\mathbb{E}_{\varphi_{d}}[h(\bm{U})|\mathbf{P}_{r}\bm{u}]}\right)\varphi_{d}(\bm{u})\mathrm{d}\bm{u}\leq
2\int_{\mathbb{R}^{d}}\lVert(\mathbf{I}-\mathbf{P}_{r}^{\mathrm{T}})\nabla
h(\bm{u})\lVert^{2}\varphi_{d}(\bm{u})\mathrm{d}\bm{u}.$

By choosing $h^{2}(\bm{u})=Z_{\beta}^{-1}\tilde{L}^{\beta}(\bm{u})$ we obtain
the KLD
$D_{\mathrm{KL}}(\tilde{p}_{\bm{y},\beta}||\tilde{p}_{\bm{y},\beta}^{(r,\star)})$
on the left-hand side of Eq. 18. With $\nabla
h(\bm{u})=\frac{1}{2}(Z_{\beta}^{-1}\tilde{L}^{\beta}(\bm{u}))^{\frac{1}{2}}\nabla\ln\tilde{L}^{\beta}(\bm{u})$,
an upper bound on
$D_{\mathrm{KL}}(\tilde{p}_{\bm{y},\beta}||\tilde{p}_{\bm{y},\beta}^{(r,\star)})$
emerges on the right-hand side of Eq. 18 as

$\displaystyle
D_{\mathrm{KL}}(\tilde{p}_{\bm{y},\beta}||\tilde{p}_{\bm{y},\beta}^{(r,\star)})$
$\displaystyle\leq\frac{1}{2}\int_{\mathbb{R}^{d}}\lVert(\mathbf{I}-\mathbf{P}_{r}^{\mathrm{T}})\nabla\ln\tilde{L}^{\beta}(\bm{u})\lVert^{2}\tilde{p}_{\bm{y},\beta}(\bm{u})\mathrm{d}\bm{u}$
$\displaystyle=\frac{1}{2}\int_{\mathbb{R}^{d}}\mathrm{tr}\left[(\mathbf{I}-\mathbf{P}_{r}^{\mathrm{T}})\beta^{2}\nabla\ln\tilde{L}(\bm{u})(\nabla\ln\tilde{L}(\bm{u}))^{\mathrm{T}}(\mathbf{I}-\mathbf{P}_{r})\right]\tilde{p}_{\bm{y},\beta}(\bm{u})\mathrm{d}\bm{u}$
$\displaystyle=\frac{1}{2}\mathrm{tr}\left[(\mathbf{I}-\mathbf{P}_{r}^{\mathrm{T}})\mathbf{H}(\mathbf{I}-\mathbf{P}_{r})\right]=\vcentcolon\frac{1}{2}\mathcal{R}(\mathbf{P}_{r},\mathbf{H}),$

where $\mathcal{R}(\mathbf{P}_{r},\mathbf{H})$ is the mean-squared error
incurred by approximating $\nabla\ln\tilde{L}(\bm{U})$ with
$\mathbf{P}_{r}^{\mathrm{T}}\nabla\ln\tilde{L}(\bm{U})$ when
$\bm{U}\sim\tilde{p}_{\bm{y},\beta}(\bm{u})$ and we define

(19)
$\mathbf{H}\vcentcolon=\beta^{2}\mathbb{E}_{\tilde{p}_{\bm{y},\beta}}\left[\nabla\ln\tilde{L}(\bm{U})(\nabla\ln\tilde{L}(\bm{U}))^{\mathrm{T}}\right].$

Our goal is to find the rank-r projection that minimizes
$\mathcal{R}(\mathbf{P}_{r},\mathbf{H})$. [60, Proposition 2.11] states that a
minimizer of $\mathcal{R}(\mathbf{P}_{r},\mathbf{H})$ over all viable
projections of rank $r$ is given by the $r$ eigenvectors of $\mathbf{H}$
corresponding to its $r$ leading eigenvalues. Let the solutions of the
eigenproblem $\mathbf{H}\bm{\phi}_{i}=\bm{\phi}_{i}\lambda_{i}$,
$\\{\bm{\phi}_{i},\lambda_{i}\\}_{i=1}^{d}$, be ordered so that
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{d}$, then collecting
$\bm{\Phi}_{r}\vcentcolon=[\bm{\phi}_{1},\bm{\phi}_{2},\dots,\bm{\phi}_{r}]\in\mathbb{R}^{d\times
r}$, the optimal projector is given as
$\mathbf{P}_{r}=\bm{\Phi}_{r}\bm{\Phi}_{r}^{\mathrm{T}}$. With this definition
of the projection and since the standard-normal prior satisfies inequality Eq.
18, the accuracy of the reduced posterior can be controlled using a tolerance
$\epsilon$ as

(20)
$D_{\mathrm{KL}}(\tilde{p}_{\bm{y},\beta}||\tilde{p}_{\bm{y},\beta}^{(r,\star)})\leq\frac{\beta^{2}}{2}\sum_{i=r+1}^{d}\lambda_{i}\leq\epsilon.$

Upon selecting $\epsilon$ and computing the eigenpairs of $\mathbf{H}$, we
choose $r$ as small as possible so that Eq. 20 holds. Efficient dimension
reduction is therefore contingent on a sharp decay of the
$\mathbf{H}$-spectrum, which is a property of the computational model $f$ and
the observation model $p(\bm{u},\bm{y})$ (i.e., prior and likelihood).

###  3.4 Method

$\bm{\Phi}_{r}$ is the eigenspace of the symmetric matrix $\mathbf{H}$ and
thus is an orthogonal basis of $\mathcal{U}_{r}$. $\bm{\Phi}_{r}$ maps the
ambient LIS coordinate $\bm{u}_{r}\in\mathcal{U}_{r}$ to its local counterpart
$\bar{\bm{u}}_{r}\in\bar{\mathcal{U}}_{r}=\mathbb{R}^{r}$ as
$\bar{\bm{u}}_{r}=\bm{\Phi}_{r}^{\mathrm{T}}\bm{u}_{r}$. In the same way, we
define
$\bm{\Phi}_{\perp}\vcentcolon=[\bm{\phi}_{r+1},\bm{\phi}_{r+2},\dots,\bm{\phi}_{d}]\in\mathbb{R}^{d\times
d-r}$, so that the ambient CS coordinate $\bm{u}_{\perp}$ is mapped to its
local counterpart
$\bar{\bm{u}}_{\perp}\in\bar{\mathcal{U}}_{\perp}=\mathbb{R}^{d-r}$ as
$\bar{\bm{u}}_{\perp}=\bm{\Phi}_{\perp}^{\mathrm{T}}\bm{u}_{\perp}$. Thus, we
can write any $\bm{u}\in\mathbb{R}^{d}$ as
$\bm{u}=\bm{\Phi}_{r}\bar{\bm{u}}_{r}+\bm{\Phi}_{\perp}\bar{\bm{u}}_{\perp}$
and collect the ambient coordinate of $\mathbb{R}^{d}$ with respect to the
basis defined by $\bm{\Phi}=[\bm{\Phi}_{r},\bm{\Phi}_{\perp}]$ as

(21) $\bar{\bm{u}}=\left[\begin{matrix}\bar{\bm{u}}_{r}\\\
\bar{\bm{u}}_{\perp}\end{matrix}\right]=\bm{\Phi}^{\mathrm{T}}\cdot\bm{u}=\left[\begin{matrix}\bm{\Phi}_{r}^{\mathrm{T}}\\\
\bm{\Phi}_{\perp}^{\mathrm{T}}\end{matrix}\right]\cdot\bm{u}.$

The LIS and CS in local and ambient coordinates are illustrated in Fig. 1.

$\mathbf{P}_{\perp}\bm{u}={\color[rgb]{0.0314,0.3961,0.0392}\bm{u}_{\perp}}$$\mathbf{P}_{r}\bm{u}={\color[rgb]{0.1216,0.3137,0.7020}\bm{u}_{r}}$$\bm{u}$${\color[rgb]{0.1216,0.3137,0.7020}\mathcal{U}_{r}}\subseteq\mathbb{R}^{d}$${\color[rgb]{0.0314,0.3961,0.0392}\mathcal{U}_{\perp}}\subseteq\mathbb{R}^{d}$$\mathcal{U}=\mathcal{U}_{r}\oplus\mathcal{U}_{\perp}$${\color[rgb]{0.1216,0.3137,0.7020}\bar{\bm{u}}_{r}}$${\color[rgb]{0.0314,0.3961,0.0392}\bar{\bm{u}}_{\perp}}$${\color[rgb]{0.1216,0.3137,0.7020}\bar{\mathcal{U}}_{r}}\subseteq\mathbb{R}^{r}$${\color[rgb]{0.0314,0.3961,0.0392}\bar{\mathcal{U}}_{\perp}}\subseteq\mathbb{R}^{(d-r)}$$\bm{\Phi}_{r}^{\mathrm{T}}$$\bm{\Phi}_{r}$$\bm{\Phi}_{\perp}^{\mathrm{T}}$$\bm{\Phi}_{\perp}$
Figure 1: Ambient space $\mathbb{R}^{d}$ along with the LIS $\mathcal{U}_{r}$
and CS $\mathcal{U}_{\perp}$ as defined by $\mathbf{P}_{r}$ and
$\mathbf{P}_{\perp}$ (left) and their local counterparts
$\bar{\mathcal{U}}_{r}$ and $\bar{\mathcal{U}}_{\perp}$ (right).

Due to orthogonality of $\bm{u}_{r}$ and $\bm{u}_{\perp}$ and the rotatinal
symmetry of the standard-normal PDF, we may factorize the prior as
$\varphi_{d}(\bm{u})=\varphi_{r}(\bar{\bm{u}}_{r})\varphi_{d-r}(\bar{\bm{u}}_{\perp})$.
With this local coordinate prior, the reduced tempered posterior Eq. 17 reads

(22)
$p_{\bm{y},\beta}^{(r,\star)}(\bar{\bm{u}})\propto\underbrace{\mathbb{E}_{\varphi_{d}}[\tilde{L}^{\beta_{t}}(\bm{U})|\bm{\Phi}_{r}\bar{\bm{u}}_{r}]\varphi_{r}(\bar{\bm{u}}_{r})}_{\begin{subarray}{c}\mathrm{reduced~{}tempered}\\\
\mathrm{posterior}\end{subarray}}\underbrace{\varphi_{d-r}(\bar{\bm{u}}_{\perp}).}_{\begin{subarray}{c}\mathrm{complementary}\\\
\mathrm{prior}\end{subarray}}$

By analogy with CEBU, in the $t$-th step of CEBUred, we approximate the
reduced tempered posterior in Eq. 22 with a parametric model
$q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r,t}):\bar{\mathcal{U}}_{r}\rightarrow\mathbb{R}_{>0}$.
The parametric, tempered posterior is

(23)
$q(\bar{\bm{u}},\bm{v}_{r,t})=q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r,t})\varphi_{d-r}(\bar{\bm{u}}_{\perp}).$

Following [58], we select a Gaussian model for
$q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r,t})$, although more complicated PDFs such
as mixture models may be used as well here. The parameter set
$\bm{v}_{r,t}=\\{\bm{\mu}_{r,t}\in\mathbb{R}^{r},\bm{\Sigma}_{r,t}\in\mathbb{R}^{r\times
r}\\}$ contains the mean vector $\bm{\mu}_{r,t}$ and covariance matrix
$\bm{\Sigma}_{r,t}$ of the Gaussian model.

In the $t$-th step of CEBUred, the new temperature $\beta_{t}$ is computed
according to Eq. 11. If $t>1$, the likelihood in ambient space is evaluated by
plugging samples from the previous’ step’s reduced biasing density
$\bar{\bm{u}}_{r,k}\sim q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r,t})$ and the
complementary prior
$\bar{\bm{u}}_{\perp,k}\sim\varphi_{d-r}(\bar{\bm{u}}_{\perp})$ in
$\bm{u}_{k}=\bm{\Phi}_{r,t-1}\bar{\bm{u}}_{r,k}+\bm{\Phi}_{\perp,t-1}\bar{\bm{u}}_{\perp,k}$.
Thereafter, the gradient covariance matrix $\mathbf{H}$ of the likelihood
function with respect to the tempered posterior is estimated to determine the
current LIS and CS projections In each step but the first ($t>1$), a self-
normalized IS estimate of $\mathbf{H}_{t}$ based on samples from the previous
biasing density $q(\bar{\bm{u}},\widehat{\bm{v}}_{t-1})$ is computed as

(24)
$\widehat{\mathbf{H}}_{t}=\frac{\beta_{t}^{2}\sum_{k=1}^{n_{\mathbf{H}}}\nabla\ln
L(\bm{u}_{k})(\nabla\ln
L(\bm{u}_{k}))^{\mathrm{T}}w_{t}(\bar{\bm{u}}_{k})}{\sum_{k=1}^{n_{\mathbf{H}}}w_{t}(\bar{\bm{u}}_{k})}~{}~{}~{}~{}~{}\begin{cases}\bar{\bm{u}}_{r,k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r,t-1})\\\
\bar{\bm{u}}_{\perp,k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\varphi_{d-r}(\bar{\bm{u}}_{\perp})\end{cases}.$

If $t=1$, the weights equal 1 and samples are drawn from the $d$-dimensional
prior in ambient space, $\varphi_{d}(\bm{u})$ by setting
$\bm{\Phi}_{r}=\mathbf{I}_{d}$ and $\bm{\Phi}_{\perp}=\bm{0}_{d}$. For any
$t>1$ the weights are computed as

(25)
$w_{t}(\bar{\bm{u}})=\frac{\tilde{L}^{\beta_{t}}(\bm{\Phi}_{r}\bar{\bm{u}}_{r}+\bm{\Phi}_{\perp}\bar{\bm{u}}_{\perp})\varphi_{r}(\bar{\bm{u}}_{r})\varphi_{d-r}(\bar{\bm{u}}_{\perp})}{q(\bar{\bm{u}},\bm{v}_{r,t-1})}=\frac{\tilde{L}^{\beta_{t}}(\bm{\Phi}_{r}\bar{\bm{u}}_{r}+\bm{\Phi}_{\perp}\bar{\bm{u}}_{\perp})\varphi_{r}(\bar{\bm{u}}_{r})}{q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r,t-1})}.$

Upon computing the spectrum of $\widehat{\mathbf{H}}_{t}$, the LIS-dimension
$r$ is selected according to Eq. 20. Once the projections $\bm{\Phi}_{r}$ and
$\bm{\Phi}_{\perp}$ are defined, the parameters of
$q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r,t})$ are computed by minimizing the KLD
$D_{\mathrm{KL}}(\tilde{p}_{\bm{y},\beta}^{(r,\star)}(\bar{\bm{u}})||q(\bar{\bm{u}},\bm{v}_{t}))$.
As in Eq. 5, this is equivalent to maximizing the negative cross-entropy
between the two distributions, i.e.,

(26)
$\begin{split}\bm{v}_{r,t}&=\operatorname*{arg\,min}\limits_{\bm{v}\in\mathcal{V}}D_{\mathrm{KL}}\left(\tilde{p}_{\bm{y},\beta}^{(r,\star)}\lVert
q(\cdot,\bm{v})\right)=\operatorname*{arg\,max}\limits_{\bm{v}\in\mathcal{V}}-H\left(\tilde{p}_{\bm{y},\beta}^{(r,\star)},q(\cdot,\bm{v})\right)\\\
&=\operatorname*{arg\,max}\limits_{\bm{v}_{r}\in\mathcal{V}_{r}}\int_{\bar{\mathcal{U}}_{r}}\int_{\bar{\mathcal{U}}_{\perp}}\mathbb{E}_{\varphi_{d}}[\tilde{L}^{\beta_{t}}(\bm{U})|\bm{\Phi}_{r}\bar{\bm{u}}_{r}]\ln\left(q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r})\right)\varphi_{r}(\bar{\bm{u}}_{r})\varphi_{d-r}(\bar{\bm{u}}_{\perp})\mathrm{d}\bar{\bm{u}}_{\perp}\mathrm{d}\bar{\bm{u}}_{r}\\\
&=\operatorname*{arg\,max}\limits_{\bm{v}_{r}\in\mathcal{V}_{r}}\int_{\bar{\mathcal{U}}_{r}}\mathbb{E}_{\varphi_{d}}[\tilde{L}^{\beta_{t}}(\bm{U})|\bm{\Phi}_{r}\bar{\bm{u}}_{r}]\ln\left(q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r})\right)\varphi_{r}(\bar{\bm{u}}_{r})\mathrm{d}\bar{\bm{u}}_{r}\\\
&=\operatorname*{arg\,max}\limits_{\bm{v}_{r}\in\mathcal{V}_{r}}\mathbb{E}_{\varphi_{d}}\left[\tilde{L}^{\beta_{t}}(\bm{\Phi}_{r}\bar{\bm{U}}_{r}+\bm{\Phi}_{\perp}\bar{\bm{U}}_{\perp})\ln\left(q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r})\right)\right].\end{split}$

Throughout Eq. 26 the normalization constant $\tilde{Z}_{t}$ has been dropped
as it is irrelevant for solving the optimization problem. The final equality
in Eq. 26 is a consequence of the factorized prior in standard-normal space,
i.e.,

(27)
$\begin{split}\mathbb{E}_{\varphi_{d}}[\tilde{L}^{\beta_{t}}(\bm{U})|\bm{\Phi}_{r}\bar{\bm{u}}_{r}]&=\int_{\bar{\mathcal{U}}_{\perp}}\tilde{L}^{\beta_{t}}(\bm{\Phi}_{r}\bar{\bm{u}}_{r}+\bm{\Phi}_{\perp}\bar{\bm{u}}_{\perp})\frac{\varphi_{r}(\bar{\bm{u}}_{r})\varphi_{d-r}(\bar{\bm{u}}_{\perp})}{\int_{\mathbb{R}^{d-r}}\varphi_{r}(\bar{\bm{u}}_{r})\varphi_{d-r}(\bar{\bm{u}}_{\perp}^{\prime})\mathrm{d}\bar{\bm{u}}_{\perp}^{\prime}}\mathrm{d}\bar{\bm{u}}_{\perp}\\\
&=\int_{\bar{\mathcal{U}}_{\perp}}\tilde{L}^{\beta_{t}}(\bm{\Phi}_{r}\bar{\bm{u}}_{r}+\bm{\Phi}_{\perp}\bar{\bm{u}}_{\perp})\varphi_{d-r}(\bar{\bm{u}}_{\perp})\mathrm{d}\bar{\bm{u}}_{\perp}.\end{split}$

An IS estimate of $\bm{v}_{r,t}$ based on samples from
$q(\bar{\bm{u}},\widehat{\bm{v}}_{r,t-1})$ reads

(28)
$\widehat{\bm{v}}_{r,t}=\operatorname*{arg\,max}\limits_{\bm{v}_{r}\in\mathcal{V}_{r}}\frac{1}{n}\sum\limits_{k=1}^{n}\ln\left(q^{(r)}(\bar{\bm{u}}_{r,k},\bm{v}_{r})\right)w_{t,\text{adj}}(\bar{\bm{u}}_{r,k},\bar{\bm{u}}_{\perp,k}),~{}~{}~{}\begin{cases}\bar{\bm{u}}_{r,k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r,t-1})\\\
\bar{\bm{u}}_{\perp,k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\varphi_{d-r}(\bar{\bm{u}}_{\perp})\end{cases}$

and requires the computation of the adjusted weights $w_{t,\text{adj}}$

(29)
$w_{t,\text{adj}}(\bar{\bm{u}}_{r},\bar{\bm{u}}_{\perp})=\frac{\tilde{L}^{\beta_{t}}(\bm{\Phi}_{r,t}\bar{\bm{u}}_{r}+\bm{\Phi}_{\perp,t}\bar{\bm{u}}_{\perp})\varphi_{r}(\bar{\bm{u}}_{r})\varphi_{d-r}(\bar{\bm{u}}_{\perp})}{q^{(d)}(\bar{\bm{u}},\widehat{\bm{v}}_{t,\text{adj}})}.$

Therein,
$\widehat{\bm{v}}_{t,\text{adj}}=\\{\bm{\mu}_{t,\text{adj}}\in\mathbb{R}^{d},\bm{\Sigma}_{t,\text{adj}}\in\mathbb{R}^{d\times
d}\\}$ represents the parameters of the $d$-dimensional Gaussian density
$q(\bar{\bm{u}},\widehat{\bm{v}}_{t-1})$ expressed with respect to the updated
orthogonal basis $\bm{\Phi}_{t}$. Computing adjusted weights with the
transformed parameters is necessary to address non-matching bases in the
numerator and denominator of Eq. 29. That is, if the new basis $\bm{\Phi}_{t}$
differs from the basis $\bm{\Phi}_{t-1}$, the complementary prior will no
longer be standard-normal with respect to $\bm{\Phi}_{t}$. The transformation
from $\bm{\Phi}_{t-1}$ to $\bm{\Phi}_{t}$ is linear whereby
$q(\bar{\bm{u}},\widehat{\bm{v}}_{t,\text{adj}})$ is Gaussian again and its
parameters with respect to $\bm{\Phi}_{t}$ can be expressed as

(30)
$\bm{\mu}_{t,\text{adj}}=\bm{\Phi}_{t}^{\mathrm{T}}\underbrace{\bm{\Phi}_{r,t-1}\bm{\mu}_{r,t-1}}_{\bm{\mu}_{t-1}},~{}~{}~{}\bm{\Sigma}_{t,\text{adj}}=\bm{\Phi}_{t}^{\mathrm{T}}\underbrace{[\bm{\Phi}_{r,t-1}\bm{\Sigma}_{r,t-1}\bm{\Phi}_{r,t-1}^{\mathrm{T}}+\bm{\Phi}_{\perp,t-1}\bm{\Phi}_{\perp,t-1}^{\mathrm{T}}]}_{\bm{\Sigma}_{t-1}}\bm{\Phi}_{t}.$

$\bm{\mu}_{t-1}$ and $\bm{\Sigma}_{t-1}$ are the mean and covariance vector in
ambient space that are subsequently transformed to the reduced spaces given
the novel basis $\bm{\Phi}_{t}$. This transformation between local and global
and subsequent subspaces in steps $t-1$ and $t$ is illustrated in Fig. 2.

${\color[rgb]{0.1216,0.3137,0.7020}\bm{P}_{r,t}\bm{u}}$$\bm{\Phi}_{\perp,t-1}^{\mathrm{T}}$$\bm{\Phi}_{\perp,t}^{\mathrm{T}}$$\bm{\Phi}_{r,t-1}^{\mathrm{T}}$$\bm{\Phi}_{r,t}^{\mathrm{T}}$$\bm{u}$${\color[rgb]{0.1216,0.3137,0.7020}\bm{P}_{r,t-1}\bm{u}}$${\color[rgb]{0.0314,0.3961,0.0392}\bm{P}_{\perp,t-1}\bm{u}}$${\color[rgb]{0.0314,0.3961,0.0392}\bm{P}_{\perp,t}\bm{u}}$$\bar{\mathcal{U}}_{r}$$q^{(r)}(\bar{u}_{r},\bm{v}_{r,t-1})$$\bar{\mathcal{U}}_{\perp}$$\varphi_{d-r}(\bar{u}_{\perp})$$\mathcal{U}_{\perp,t-1}$$\mathcal{U}_{\perp,t}$$\mathcal{U}_{r,t-1}$$\mathcal{U}_{r,t}$
Figure 2: Left: Mapping between two subsequent orthogonal bases
$\bm{\Phi}_{t-1}$ and $\bm{\Phi}_{t}$ in ambient space. Right: Mapping from
the two ambient bases to the local LIS (upper right) and CS (lower right).

###  3.5 Choosing $n_{\mathbf{H}}$ and $n$ adaptively

$n_{\mathbf{H}}$ is the number of log-likelihood gradient evaluations used to
compute $\widehat{H}_{t}$ in Eq. 24. $n$ on the other hand is the number of
direct evaluations of the tempered likelihood used to estimate the parameters
of the $t$-th biasing density in Eq. 28. In the absence of $f$-solvers that
are specifically geared towards efficient gradient evaluation such as adjoint
solvers [1], computing $\nabla\ln L(\bm{u})$ is considerably more expensive
than evaluating $L(\bm{u})$.

[6] suggests a heuristic for determining $n_{\mathbf{H}}$ when estimating the
second-moment matrix of the gradient $\nabla f$ of a Lipschitz-continuous
function $f$, i.e., $\nabla f\leq a$, in order to discover an _active
subspace_ of $f$. They use work of [23] on the spectrum of sums of
($n_{\mathbf{H}}$) random matrices to establish bounds on the relative
accuracy of the estimated spectrum of $\mathbf{C}=\mathbb{E}_{p}[\nabla
f(\nabla f)^{\mathrm{T}}]$. In the context of CEBUred, we have
$f=\beta\log\tilde{L}$ and $p=\tilde{p}_{\bm{y},\beta}$.

[6, Corollary 3.5] states that for $\varepsilon\in(0,1]$,
$\mathbb{P}[|\widehat{\lambda}_{r}-\lambda_{r}|/\lambda_{r}\leq\varepsilon]\leq
2d^{-b}$ if the spectrum of $\widehat{\mathbf{H}}$,
$\\{\widehat{\lambda}_{i}\\}_{i=1}^{d}$, is computed with $n_{\mathbf{H}}\geq
4(b+1)a\lambda_{1}\ln(d)/(\lambda_{r}\varepsilon)^{2}$ log-likelihood gradient
samples. Drawing on a matrix Bernstein inequality in [56], [6, Corollary 3.8]
states that for $\varepsilon\in(0,1]$,
$\mathbb{P}[\mathbf{H}-\widehat{\mathbf{H}}\lVert_{2}/\lVert\mathbf{H}\lVert_{2}\leq\varepsilon]\leq
2m^{1-3c/8}$ (the $2$-norm of a matrix here is its spectral norm, which also
corresponds to its largest singular value) when estimating
$\widehat{\mathbf{H}}$ with at least $n_{\mathbf{H}}\geq
ca^{2}\ln(d)/(\lambda_{1}\varepsilon^{2})$ samples. Finally, choosing
$\varepsilon$ such that
$\varepsilon\leq(\lambda_{r}-\lambda_{r+1})/(5\lambda_{1})$ and using this
last lower bound on $n_{\mathbf{H}}$, the distance between the image of the
local estimated and true LIS projections is bounded with high probability as
well: the distance as measured with the spectral norm
$d(\mathrm{Im}(\bm{\Phi}_{r}),\mathrm{Im}(\widehat{\bm{\Phi}}_{r}))=\lVert\bm{\Phi}_{r}\bm{\Phi}_{r}^{\mathrm{T}}-\widehat{\bm{\Phi}}_{r}\widehat{\bm{\Phi}}_{r}^{\mathrm{T}}\lVert_{2}=\lVert\bm{\Phi}_{r}^{\mathrm{T}}\widehat{\bm{\Phi}}_{\perp}\lVert_{2}$
is bounded as
$\mathbb{P}[\lVert\bm{\Phi}_{r}^{\mathrm{T}}\widehat{\bm{\Phi}}_{\perp}\lVert_{2}\leq
4\lambda_{1}\varepsilon/(\lambda_{r}-\lambda_{r+1})]\leq 2m^{1-3c/8}$
according to [6, Corollary 3.10]. [6] translates this bound into a heuristic
on account of $c$, $a$ and the true spectrum $\\{\lambda_{i}\\}_{i=1}^{d}$
being unknown in many use cases involving numerical/simulation models $f$. The
heuristic emerges by summarizing all unknown constants in a fudge factor
$\alpha_{\mathbf{H}}$ resulting in

(31) $n_{\mathbf{H}}=\alpha_{\mathbf{H}}r\ln(d),$

where [6] recommends $\alpha_{H}\in[2,10]$ and the target rank $r$ corresponds
to the smallest eigenvalue $\lambda_{r}$ that shall be estimated with the
desired relative accuracy $\varepsilon$. Remarkably, the effort scales
logarithmically with the ambient space dimension $d$ suggesting that we can
hope to estimate $\mathbf{H}$ with relatively few log-likelihood gradient
samples even in very high-dimensional settings. As the target rank $r$ is not
known a priori, we detail an iterative procedure to jointly determine $r$ and
$n_{\mathbf{H}}$ in Alg. 2.

The required number of samples in each level of the CEBU procedure in turn
depends on the adaptively selected LIS-rank $r$ through the number of
parameters that have to be fitted in the Gaussian reduced biasing density
$n_{\mathrm{par}}$. In particular, an $r$-variate Gaussian requires fitting
$n_{\mathrm{par}}=r(r+3)/2$ parameters. To select the number of samples
required to accurately estimate $\widehat{\bm{v}}_{r,t}$ in Eq. 28, we use the
following heuristic:

(32)
$n(r)={\underbrace{\vphantom{\frac{1}{2}}\alpha_{\mathrm{par}}}_{\begin{subarray}{c}\mathrm{~{}fudge~{}}\\\
\mathrm{factor}\end{subarray}}}~{}{\underbrace{\frac{1}{2}r(r+3)}_{\begin{subarray}{c}\mathrm{number~{}of}\\\
\mathrm{parameters}\end{subarray}}}~{}{\underbrace{\vphantom{\frac{1}{2}}(1+\delta_{w}^{2})}_{\begin{subarray}{c}\mathrm{inverse}\\\
\mathrm{nESS}\end{subarray}}},$

where [6] recommend to chose $\alpha_{\mathrm{par}}\in[2,10]$. In case an
adjoint solver is used for $f$, the estimation of $\widehat{H}_{t}$ as in Eq.
24 will return $n_{\mathbf{H}}$ likelihood evaluations as a byproduct that can
be utilizied in estimating $\widehat{\bm{v}}_{r,t}$ so that only an effective
$n(r)-n_{\mathbf{H}}$ new samples need to be drawn and evaluated at each
level.

Algorithm 2 adapt_H

Input Likelihood and log-gradient $\tilde{L}$ and $\nabla\ln\tilde{L}$,
reduced biasing density $q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r})$, local LIS & CS
projections $\bm{\Phi}_{r}$ & $\bm{\Phi}_{\perp}$, fudge factor
$\alpha_{\mathbf{H}}$, error tolerance $\epsilon$, temperature $\beta_{0}$
Output Subspace samples
$\bar{\bm{\mathrm{U}}}_{r},\bar{\bm{\mathrm{U}}}_{\perp}$, Likelihood samples
$\bm{\ell}$, local LIS & CS projections $\bm{\Phi}_{r}$
& $\bm{\Phi}_{\perp}$, temperature $\beta$, LIS-rank $r$, # of
$\mathbf{H}$-samples $n_{\mathbf{H}}$

1: Set $r\leftarrow 1$, $d\leftarrow\mathrm{rowdim}(\bm{\Phi}_{r})$,
$n_{\mathbf{H}}\leftarrow\alpha_{\mathbf{H}}\ln(d)$, $\Delta n\leftarrow
n_{\mathbf{H}}$,
$\bar{\bm{\mathrm{U}}}_{r},\bar{\bm{\mathrm{U}}}_{\perp},\bm{\ell},d\bm{\ell}\leftarrow[~{}]$

2: while $\Delta n>0$ do

3: Sample $\bar{\bm{\mathrm{U}}}_{r,\mathrm{add}}\in\mathbb{R}^{\Delta n\times
r}\leftarrow\\{\bar{\bm{u}}_{r,k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r})\\}_{k=1}^{\Delta n}$

4: Sample $\bar{\bm{\mathrm{U}}}_{\perp,\mathrm{add}}\in\mathbb{R}^{\Delta
n\times(d-r)}\leftarrow\\{\bar{\bm{u}}_{\perp,k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\varphi_{d-r}(\bar{\bm{u}}_{\perp})\\}_{k=1}^{\Delta n}$

5: Append vertically
$\bar{\bm{\mathrm{U}}}_{r}\leftarrow[\bar{\bm{\mathrm{U}}}_{r},\bar{\bm{\mathrm{U}}}_{r,\mathrm{add}}],~{}\bar{\bm{\mathrm{U}}}_{\perp}\leftarrow[\bar{\bm{\mathrm{U}}}_{\perp},\bar{\bm{\mathrm{U}}}_{\perp,\mathrm{add}}]$

6: Compute
$\bm{\ell}_{\mathrm{add}}\leftarrow\tilde{L}(\bar{\bm{U}}_{r,\mathrm{add}}^{\mathrm{T}}\bm{\Phi}_{r}+\bar{\bm{U}}_{\perp,\mathrm{add}}^{\mathrm{T}}\bm{\Phi}_{\perp})$
and
$d\bm{\ell}_{\mathrm{add}}\leftarrow\nabla\tilde{L}(\bar{\bm{U}}_{r,\mathrm{add}}^{\mathrm{T}}\bm{\Phi}_{r}+\bar{\bm{U}}_{\perp,\mathrm{add}}^{\mathrm{T}}\bm{\Phi}_{\perp})$

7: Append vertically
$\bm{\ell}\leftarrow[\bm{\ell},\bm{\ell}_{\mathrm{add}}],d\bm{\ell}\leftarrow[d\bm{\ell},d\bm{\ell}_{\mathrm{add}}]$

8: Evaluate $\beta$ and $\mathbf{w}$ in function of $\beta_{0}$,
$\bar{\bm{U}}_{r}$, $\bar{\bm{U}}_{\perp}$, $\bm{\ell}$,
$q^{(r)}(\bar{\bm{U}}_{r},\bm{v}_{r})$ with Eq. 11 & Eq. 25

9: Evaluate $\widehat{\mathbf{H}}$ in function of $\beta$, $\mathbf{w}$ and
$d\bm{\ell}$ with Eq. 24

10: Evaluate
$\\{\bm{\phi}_{i},\lambda_{i}\\}_{i=1}^{d}\leftarrow\operatorname*{solve}\\{\bm{\phi}\in\mathbb{R}^{d},\lambda\in\mathbb{R}:\widehat{\mathbf{H}}\bm{\phi}=\bm{\phi}\lambda\\}$

11: Select $r$ in function of $\epsilon$ and
$\\{\bm{\phi}_{i},\lambda_{i}\\}_{i=1}^{d}$ with Eq. 20

12: Set $\Delta n=\alpha_{H}r\ln(d)-n_{\mathbf{H}}$

13: end while

14: Define $\bm{\Phi}_{r}\leftarrow[\bm{\phi}_{1},\dots,\bm{\phi}_{r}]$,
$\bm{\Phi}_{\perp}\leftarrow[\bm{\phi}_{r+1},\dots,\bm{\phi}_{d}]$

15: return $\bar{\bm{\mathrm{U}}}_{r}$, $\bar{\bm{\mathrm{U}}}_{\perp}$,
$\bm{\ell}$, $\bm{\Phi}_{r}$, $\bm{\Phi}_{\perp}$, $\beta$, $r$,
$n_{\mathbf{H}}$

The CEBUred algorithm is summarized in Alg. 3.

Algorithm 3 CE-BU-red

Input Likelihood and log-gradients $\tilde{L}$ and $\nabla\ln\tilde{L}$,
transform $T$, parameters
$n_{\mathrm{eff}}^{*},\alpha_{\mathbf{H}},\alpha_{\mathrm{par}},\epsilon$, #
post. samples $N$,
Output posterior samples $\mathbf{\Theta}_{\mathrm{post}}$, estimated evidence
$\widehat{Z}$

1: Set $t\leftarrow 0$, $\beta_{0}\leftarrow 0$,
$\bm{\Phi}_{r,0}\leftarrow\mathbf{I}_{d\times d}$,
$\bm{\Phi}_{\perp,0}\leftarrow 0$,
$\widehat{\bm{v}}_{r,0}=\\{\bm{0}_{d},\mathbf{I}_{d\times d}\\}$,

2: while $\beta_{t}<1$ do

3: $t\leftarrow t+1$

4:
$\bar{\bm{\mathrm{U}}}_{r,0},\bar{\bm{\mathrm{U}}}_{\perp,0},\bm{\ell}_{0},\bm{\Phi}_{r,t},\bm{\Phi}_{\perp,t},\beta_{t},r_{t},n_{\mathbf{H}}$$\leftarrow\mathrm{adapt\\_H}(\tilde{L},\nabla\ln\tilde{L},q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r,t-1}),\bm{\Phi}_{r,t-1},\bm{\Phi}_{\perp,t-1},\beta_{t-1},\alpha_{\mathbf{H}})$

5: Compute the required number of samples $n$ in function of $r_{t}$ with Eq.
32

6: Sample
$\bar{\bm{\mathrm{U}}}_{r,\mathrm{add}}\in\mathbb{R}^{(n-n_{\mathbf{H}})\times
r_{t}}\leftarrow\\{\bar{\bm{u}}_{r,k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}q^{(r)}(\bar{\bm{u}}_{r},\bm{v}_{r,t-1})\\}_{k=1}^{n-n_{\mathbf{H}}}$

7: Sample
$\bar{\bm{\mathrm{U}}}_{\perp,\mathrm{add}}\in\mathbb{R}^{(n-n_{\mathbf{H}})\times(d-r_{t})}\leftarrow\\{\bar{\bm{u}}_{\perp,k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\varphi_{d-r_{t}}(\bar{\bm{u}}_{\perp})\\}_{k=1}^{n-n_{\mathbf{H}}}$

8: Compute $\bm{\ell}_{\mathrm{add}}\in\mathbb{R}^{(n-n_{\mathbf{H}})\times
1}\leftarrow\tilde{L}(\bar{\bm{U}}_{r,\mathrm{add}}^{\mathrm{T}}\bm{\Phi}_{r,t}+\bar{\bm{U}}_{\perp,\mathrm{add}}^{\mathrm{T}}\bm{\Phi}_{\perp,t})$

9: Join vertically
$\bar{\bm{\mathrm{U}}}_{r}\leftarrow[\bar{\bm{\mathrm{U}}}_{r,0},\bar{\bm{\mathrm{U}}}_{r,\mathrm{add}}],~{}\bar{\bm{\mathrm{U}}}_{\perp}\leftarrow[\bar{\bm{\mathrm{U}}}_{\perp,0},\bar{\bm{\mathrm{U}}}_{\perp,\mathrm{add}}],~{}\bm{\ell}\leftarrow[\bm{\ell}_{0},\bm{\ell}_{\mathrm{add}}]$

10: Compute
$\widehat{\bm{v}}_{t,\text{adj}}(\widehat{\bm{v}}_{r,t-1},\bm{\Phi}_{r,t-1},\bm{\Phi}_{r,t},\bm{\Phi}_{\perp,t-1},\bm{\Phi}_{\perp,t})$
according to Eq. 30

11: Compute the adjusted weights
$\mathbf{w}_{t,\text{adj}}\in\mathbb{R}^{n\times 1}\leftarrow
w_{t,\text{adj}}(\bar{\bm{\mathrm{U}}}_{r},\bar{\bm{\mathrm{U}}}_{\perp},\bm{\ell},\widehat{\bm{v}}_{t,\text{adj}})$
with Eq. 29

12: Compute $\widehat{\bm{v}}_{r,t}$ with Eq. 28

13: end while

14: Evaluate $\bm{\mathrm{w}}_{\mathrm{final}}\in\mathbb{R}^{n\times
1}\leftarrow\left\\{\frac{\tilde{L}(\bar{\mathbf{U}}_{r}^{\mathrm{T}}\bm{\Phi}_{r,t}+\bar{\mathbf{U}}_{\perp,t}^{\mathrm{T}}\bm{\Phi}_{\perp})\varphi_{r,t}(\bar{\mathbf{U}}_{r}^{\mathrm{T}}\bm{\Phi}_{r,t})}{q^{(r)}(\bar{\mathbf{U}}_{r}^{\mathrm{T}}\bm{\Phi}_{r,t},\widehat{\bm{v}}_{r,t})}\right\\}_{k=1}^{n}$

15: Estimate evidence
$\widehat{Z}\leftarrow\frac{1}{n}\sum_{k=1}^{n}\mathrm{w}_{\mathrm{final},k}$

16: Normalize weights
$\bar{\bm{\mathrm{w}}}_{\mathrm{final}}\leftarrow\bm{\mathrm{w}}_{\mathrm{final}}/(n\widehat{Z})$

17: $\mathbf{U}_{\mathrm{post}}\leftarrow$ Resample (with replacement) $N$
times from
$\bar{\mathbf{U}}_{r}^{\mathrm{T}}\bm{\Phi}_{r,t}+\bar{\bm{U}}_{\perp,t}^{\mathrm{T}}\bm{\Phi}_{\perp}$
with weighting $\bar{\bm{\mathrm{w}}}_{\mathrm{final}}$

18: $\mathbf{\Theta}_{\mathrm{post}}\leftarrow
T^{-1}(\mathbf{U}_{\mathrm{post}})$

19: return $\mathbf{\Theta}_{\mathrm{post}}$, $\widehat{Z}$

##  4 Experimental results

We perform two numerical examples to demonstrate the capability and test for
potential limitations of CEBUred. In the first example, we compare the
computational cost and accuracy of CEBUred and CEBU in dependency of the
ambient space dimension and verify the results with an analytical solution. In
the second example, we examine the performance of CEBUred for different error
thresholds as defined by Eq. 20. Both methods are implemented with a Gaussian
model as parametric IS density. In both examples, we infer a material
parameter random field based on model output observations. We measure the
quality of posterior random field approximations $\widehat{Y}$ against a
reference solution $Y$ (either analytical or numerical) in terms of the
following spatially averaged relative mean and variance errors:

(33)
$\varepsilon_{\mu_{Y}}=\frac{\lVert\mu_{Y}(\bm{x})-\hat{\mu}_{Y}(\bm{x})\rVert_{2}}{\lVert\mu_{Y}(\bm{x})\rVert_{2}}\quad\text{and}\quad\varepsilon_{\sigma^{2}_{Y}}=\frac{\lVert\sigma^{2}_{Y}(\bm{x})-\hat{\sigma}^{2}_{Y}(\bm{x})\rVert_{2}}{\lVert\sigma^{2}_{Y}(\bm{x})\rVert_{2}}.$

###  4.1 1D Cantilever beam

####  4.1.1 Problem description

We consider an Euler-Bernoulli beam with one clamped and one free end. It has
length $L=5~{}\mathrm{m}$ and a point load of $P=20~{}\mathrm{kN}$ acting on
the free end (Fig. 4). Its bending moment $M(x)$ can be obtained from the
Euler-Bernoulli equation and reads [5]

(34) $M(x)=-E(x)I(x)\frac{d^{2}w(x)}{dx^{2}},$

where $E(x)$ is the beam’s Young’s modulus and $I(x)$ is its moment of
inertia. Both can be summarized as the beam’s axial flexibility
$F(x)=1/E(x)I(x)$. The bending moment of the cantilever beam is computed as
$M(x)=-P(L-x)$. Hence, the vertical deformation is given by

(35) $w(x,F(x))=P\int_{0}^{x}\int_{0}^{s}(L-x)F(x)\mathrm{d}t\mathrm{d}s.$

The axial flexibility $F(x)$ is considered uncertain and spatially variable
along the beam axis. We assign a homogeneous Gaussian prior random field with
mean $\mu_{F}=10^{-4}\mathrm{kN}^{-1}\mathrm{m}^{-2}$, standard deviation
$\sigma_{F}=3.5\cdot 10^{-5}\mathrm{kN}^{-1}\mathrm{m}^{-2}$ and exponential
autocorrelation kernel

(36) $\rho(x,x^{\prime};l)=\exp{\left(-\frac{\lVert
x-x^{\prime}\rVert_{1}}{l}\right)},$

where $l$ is the correlation length. The correlation length of the random
field of the axial flexibility is $l_{F}=2\mathrm{m}$.

Figure 3: Cantilever beam with point load.

Figure 4: True deformation and measurements.

The forward model is given by a finite element (FE) model employing $100$
Euler-Bernoulli beam elements with cubic shape functions. The goal is to
obtain samples of the posterior distribution of the axial flexibility given
$n_{\mathrm{obs}}=50$ equally spaced measurements
$\\{x_{\mathrm{meas},i}\\}_{i=1}^{n_{\mathrm{obs}}}$ of the vertical
deformation along the beam axis (see Fig. 4). Adjoint methods [1] are a
computationally efficient tool for obtaining the model gradients required to
compute $\nabla\tilde{L}$ as long as the number of model outputs of which
derivatives are computed ($n_{\mathrm{obs}}$) is smaller than the number of
model inputs ($d$) with respect to which derivatives are computed. In the
context of this example, the adjoint method is thus used for the setting
$d=100$ only and the direct method is used in all other settings.

We assume the measurements to be corrupted by the additive, centered Gaussian
noise vector $\bm{\eta}\sim\mathcal{N}(\bm{0},\bm{\Sigma}_{\eta\eta})$. The
noise covariance matrix is defined as
$[\bm{\Sigma}_{\eta\eta}]_{ij}=\sigma_{\eta}^{2}\rho(x_{\mathrm{meas},i},x_{\mathrm{meas},j})$
with noise standard deviation $\sigma_{\eta}=0.001\mathrm{m}$, exponential
correlation kernel $\rho(\cdot,\cdot;l_{\eta})$ and correlation length
$l_{\eta}=1\mathrm{m}$. The random vector describing the vertical deformations
in data space $\mathbb{R}^{n_{\mathrm{obs}}}$, i.e., at the $n_{\mathrm{obs}}$
measurement locations $\\{x_{\mathrm{meas},i}\\}_{i=1}^{n_{\mathrm{obs}}}$ is
defined as

(37) $\tilde{\bm{w}}=\bm{w}+\bm{\eta}.$

Given a set of realizations of $\tilde{\bm{w}}$, i.e., observational data
$\tilde{\bm{y}}$, the likelihood function reads

(38)
$L(\bm{F})=\frac{1}{\sqrt{(2\pi)^{n_{\mathrm{obs}}}\det(\bm{\Sigma}_{\eta\eta})}}\exp{\left([\tilde{\bm{y}}-\mathcal{G}(\bm{F})]\bm{\Sigma}_{\eta\eta}^{-1}[\tilde{\bm{y}}-\mathcal{G}(\bm{F})]^{\mathrm{T}}\right)},$

where $\mathcal{G}(\cdot)$ represents the FE-model and returns the vertical
deformations of the beam at the $n_{\mathrm{obs}}$ measurement locations.

The measurements for this example are obtained by generating a single random
realization of the prior random field of the axial flexibility, solving Eq. 35
numerically at $1001$ equally spaced discretization points and then adding
randomly generated noise according to Eq. 37 to the solutions at the locations
of the measurements $\\{x_{\mathrm{meas},i}\\}_{i=1}^{n_{\mathrm{obs}}}$. By
using the analytical expression instead of the FE-model for the generation of
the measurements, we avoid the so-called ’inverse crime’ [29].

####  4.1.2 Analytical posterior

The following derivations closely follow [57] where the example is
investigated as well. Since $F(x)$ is Gaussian and $w(x,F(x))$ is a linear
function of $F(x)$ Eq. 34, the prior distribution of $w(x)$ is also Gaussian.
Its mean and covariance read

(39a)
$\mu_{w}(x)=P\int_{0}^{x}\int_{0}^{s}(L-t)F(t)\mathrm{d}t\mathrm{d}s=\frac{P\mu_{F}}{6}x^{2}(3L-x)\quad\text{and}$
(39b)
$\Sigma_{ww}(x,x^{\prime})=P\int_{0}^{x^{\prime}}\int_{0}^{x}\int_{0}^{s^{\prime}}\int_{0}^{s}(L-t)(L-t^{\prime})\Sigma_{FF}(t,t^{\prime})\mathrm{d}t\mathrm{d}t^{\prime}\mathrm{d}s\mathrm{d}s^{\prime}.$

The explicit expression of Eq. 39b is obtained using a computer algebra system
and omitted here due to its tedious form.

An analytical solution of the posterior of the axial flexibility can be
derived, since both the prior and the likelihood are Gaussian [45]. To this
end, the Gaussian random vector $\bm{F}^{\prime}=[\bm{F};\tilde{\bm{w}}]$ is
considered, which contains the discretized random flexibility field,
$\bm{F}\in\mathbf{R}^{n}$ and the $n_{\mathrm{obs}}$ deformation measurements
$\tilde{\bm{w}}\in\mathbf{R}^{n_{\mathrm{obs}}}$. The mean vector and
covariance matrix of $\bm{F}^{\prime}$ may be partitioned as

(40) $\bm{\mu}_{F^{\prime}}=\begin{bmatrix}\bm{\mu}_{F}\\\
\bm{\mu}_{\tilde{w}}\end{bmatrix}\quad\text{and}\quad\bm{\Sigma}_{F^{\prime}F^{\prime}}=\begin{bmatrix}\bm{\Sigma}_{FF}&\bm{\Sigma}_{F\tilde{w}}\\\
\bm{\Sigma}_{F\tilde{w}}^{\mathrm{T}}&\bm{\Sigma}_{\tilde{w}\tilde{w}}\end{bmatrix}.$

As $\bm{F}^{\prime}$ is jointly Gaussian, the posterior
$\bm{F}|\tilde{\bm{y}}$ is Gaussian as well and has PDF

(41)
$p(\bm{f}|\tilde{\bm{y}})=\frac{1}{\sqrt{(2\pi)^{n}\det(\bm{\Sigma}_{FF|\tilde{y}})}}\exp{\left(-\frac{1}{2}[\bm{f}-\bm{\mu}_{F|\tilde{y}}]^{\mathrm{T}}\bm{\Sigma}_{FF|\tilde{y}}^{-1}[\bm{f}-\bm{\mu}_{F|\tilde{y}}]\right)}.$

The posterior mean and covariance matrix are equal to the following
conditional mean $\bm{\mu}_{F|\tilde{y}}$ and covariance matrix
$\bm{\Sigma}_{FF|\tilde{y}}$:

(42)
$\bm{\mu}_{F|\tilde{y}}=\bm{\mu}_{F}+\bm{\Sigma}_{F\tilde{w}}\bm{\Sigma}_{\tilde{w}\tilde{w}}^{-1}(\tilde{\bm{y}}-\bm{\mu}_{\tilde{w}})\quad\text{and}\quad\bm{\Sigma}_{FF|\tilde{y}}=\bm{\Sigma}_{FF}-\bm{\Sigma}_{F\tilde{w}}\bm{\Sigma}_{\tilde{w}\tilde{w}}^{-1}\bm{\Sigma}_{F\tilde{w}}^{\mathrm{T}}.$

All quantities in Eq. 42 are computed within the partition in Eq. 40 except
from $\bm{\mu}_{\tilde{w}}$, which is obtained by
$\mathbb{E}[\tilde{\bm{w}}]=\mathbb{E}[\bm{w}+\bm{\eta}]=\bm{\mu}_{w}$.

####  4.1.3 Parameters of numerical study

The flexibility random field is discretized in space using a midpoint method
[12] with $d$ collocation points. $d$ is therefore the ambient space dimension
of the Bayesian inverse problem, where scenarios $d=\\{10,25,50,100\\}$ are
investigated. We use CEBU and CEBUred to obtain samples from the
$d$-dimensional posterior distribution of the axial flexibility given a set of
$n_{\mathrm{obs}}=50$ measurements. We use $\delta_{v,target}=1.5$ and
$n=1000$ samples per level for all CEBU runs. For CEBUred we chose
$\delta_{v,target}=1.5$, $\alpha_{\mathrm{par}}=4$, $\alpha_{\mathbf{H}}=6$
and $\varepsilon=1.0$. Results are averaged over $54$ repeated runs of both
CEBU and CEBUred.

####  4.1.4 Discussion of results

Fig. 5 shows the posterior flexibility fields obtained with both CEBU and
CEBUred at varying ambient dimension. At $d=5$, the results obtained with both
CEBU and CEBUred coincide with the analytical reference posterior as indicated
by the almost congruent scatter points in the top left panel of Fig. 5.
However, discretizing the flexibility field with only $5$ subparts does not
allow for an accurate representation of the posterior field at the clamping.
There, the axial flexibility exerts the strongest influence on the beam
deformation thus requiring a finer discretization.

Figure 5: Axial flexibility posterior field: mean and 95% credible intervals.
CEBU and CEBUred solutions are plotted at the random field collocation points
(midpoints).

As $d$ increases, the results obtained by CEBU deteriorate, as indicated by
both an increasing deviation of the CEBU solution from its analytical
counterpart in both mean and 95% posterior credible bounds. At $d>11$, CEBU
has too few samples available to accurately fit all biasing density parameters
in ambient space. CEBUred, on the other hand, agrees closely with the
analytical solution if the chosen discretization is fine enough.

Figure 6: Left: Average number of selected LIS dimension with increasing
CEBUred step index. Right: Average number of evaluated samples per step with
increasing CEBUred step index. The solid and the dashed lines represent the
total number of model and model gradient evaluations, respectively.

Figure 7: Relative posterior mean error for different combinations of
$\alpha_{\mathbf{H}}$ and $\alpha_{\mathrm{par}}$ at $d=100$.

Fig. 7 shows the number of selected ranks $r$ (corresponds to the dimension of
the effectively used subspace (LIS) in CEBUred) and number of samples plotted
over the CEBUred step index. The number of LIS dimensions reduces to $r=1$
within the first step for all tested $d$. The number of runs per number of
steps for different $d$ are shown in Table 1. At $d=5$, one of the $54$ runs
terminated after $6$ steps, whereas all other simulations terminated after a
maximum of $5$ steps. Table 1 suggests that for the given FE-discretization,
finer random field discretizations tend to stabilize the simulation in the
sense that most runs require the same number of steps.

Table 1: Number of runs broken down according to required number of CEBUred steps at varying $d$. # of steps | $d=5$ | $d=25$ | $d=50$ | $d=100$
---|---|---|---|---
2 | $1$ | $0$ | $0$ | $0$
3 | $15$ | $5$ | $1$ | $1$
4 | $29$ | $43$ | $42$ | $46$
5 | $8$ | $6$ | $11$ | $7$
6 | $1$ | $0$ | $0$ | $0$
Table 2: Likelihood and gradient evaluations for the beam example per run of CEBU and CEBUred (CEBUred: averaged over 54 runs). Problem | CEBU | CEBUred | CEBUred
---|---|---|---
dimension | (Likelihood calls) | (Likelihood calls) | (Gradient calls)
$d=5$ | 3796 | 170.0 | 54.3
$d=10$ | 3870 | 181.0 | 115.8
$d=25$ | 3833 | 184.6 | 143.5
$d=100$ | 3833 | 182.7 | 168.6

For all investigated ambient dimensions $d$ and in each CEBUred step, the beam
problem possesses very low-dimensional (likelihood-informed) subspaces, in
which the inverse problem can be solved efficiently ($r=1-2$). In this LIS,
significantly less samples are required to accurately characterize biasing
densities compared to CEBU, which operates in $d$-dimensional ambient space.
Fig. 7 (solid lines) illustrates the correspondence of the number of required
samples in CEBUred with the reduced space dimension (rank $r$). As $d$
increases, more gradients are evaluated (Fig. 7, right, dashed lines). This is
due to the factor $\ln d$ in Eq. 31.

Table 3: Relative posterior mean and variance error at varying $d$ averaged over 54 repeated runs of CEBUred. error | $d=5$ | $d=25$ | $d=50$ | $d=100$
---|---|---|---|---
$\varepsilon_{\mu_{F}}$ | $0.0588$ | $0.0219$ | $0.0167$ | $0.0142$
$\varepsilon_{\sigma^{2}_{F}}$ | $0.2168$ | $0.1039$ | $0.1090$ | $0.0907$

In Table 2 we list the average number of required likelihood and likelihood
gradient evaluations for both CEBU and CEBUred. The number of required
likelihood evaluations within CEBUred remains approximately constant across
all investigated dimensions and is more than an order of magnitude lower
compared to number of evaluations required by CEBUred. The numbe of likelihood
gradient evaluations grows with $d$ but remains below the number of likelihood
evaluations. Depending on the method of evaluating these gradients, a gradient
call may however be considerably more expensive than a likelihood call.

$\mu_{F}(\bm{x})$ and $\sigma^{2}(\bm{x})$ are the analytical posterior mean
and variance, respectively, evaluated at the discretization points and
$\hat{\mu}_{F}(\bm{x})$ and $\hat{\sigma}^{2}_{F}(\bm{x})$ are their sample-
based counterparts obtained with CEBUred. Table 3 shows the relative error of
the mean and the variance for the different dimensions $d$. What is not
immediately obvious from the plots in Fig. 5 is that a finer discretization
indeed leads to smaller relative errors. However, the decrease slows down from
$d=25$ and is rather small between $d=50$ and $d=100$. Fig. 7 shows the
relative posterior mean error $\varepsilon_{\mu_{F}}$ for different
combinations of the fudge factors $\alpha_{H}$ and $\alpha_{\mathrm{par}}$ at
$d=100$. The relative posterior mean error decreases signifcantly between
$\alpha_{\mathbf{H}}=2$ and $\alpha_{\mathbf{H}}=6$ and remains constant as
$\alpha_{\mathbf{H}}$ is increased from $6$ to $10$. This is likely due to the
fixed error threshold of $\epsilon=1.0$, which prescribes the approximation
quality of the optimal projector. Once the relevant part of the
$\mathbf{H}$-spectrum (corresponding to the choice of $\epsilon$) is estimated
accurately, increasing $\alpha_{\mathbf{H}}$ bears no further effect. In this
case, larger values of $\alpha_{\mathrm{par}}$ only lead to a larger number of
model evaluations, which in turn can lead to a better fit of the biasing
density. At $\alpha_{\mathbf{H}}=2$, the error decreases with increasing
$\alpha_{\mathrm{par}}$, whereas it remains approximately constant when
increasing $\alpha_{\mathrm{par}}$ at $\alpha_{\mathbf{H}}\geq 6$. In the
latter case, the large number of gradient evaluations (each of which also
yields a model evaluation) are already sufficient to accurately estimate the
parameters of the reduced biasing density such that increasing
$\alpha_{\mathrm{par}}$ will not further reduce the error.

###  4.2 2D plate in plane stress

The example was first presented in [33] in the context of uncertainty
quantification. We consider the adapted version from [57]. Through this
example, we investigate how the accuracy of the resulting posterior improves
by using different error thresholds $\epsilon$ as defined by Eq. 20.

####  4.2.1 Problem description

We consider a 2D square steel plate in plane-stress with side length
$32~{}\mathrm{cm}$, thickness $t=1~{}\mathrm{cm}$ and a hole with radius
$r=2~{}\mathrm{cm}$ located at its center (Fig. 8). The plate is clamped at
the left-hand side and loaded with a constant line load
$q=6~{}\mathrm{kN}~{}\mathrm{cm}^{-2}$ acting on its right-hand side. The
plate has density $\rho=7850~{}\mathrm{kg}~{}\mathrm{m}^{-3}$, which is
required to account for body forces (oriented in negative $x_{2}$-direction),
and the Poisson ratio is $\nu=0.29$.

Figure 8: Left: FE-model of the plate. The red squares mark the positions of
the strain gauges. Center and right: True fields of the strains in $x_{1}$\-
and $x_{2}$-direction.

Assuming plane stress, the displacement field
$\bm{u}(x_{1},x_{2})=[u_{x_{1}}(x_{1},x_{2}),u_{x_{2}}(x_{1},x_{2})]^{\mathrm{T}}$
can be computed implicitly based on elasticity theory through a set of
elliptic PDEs (Cauchy-Navier equations) [28]:

(43)
$G(x_{1},x_{2})\nabla^{2}\bm{u}(x_{1},x_{2})+\frac{E(x_{1},x_{2})}{2(1-\nu)}\nabla(\nabla\cdot\bm{u}(x_{1},x_{2}))+\bm{B}=0.$

$G(x_{1},x_{2}):=E(x_{1},x_{2})/(2(1+\nu))$ is the shear modulus,
$E(x_{1},x_{2})$ is Young’s modulus, and
$\bm{B}=[b(x_{1}),b(x_{2})]^{\mathrm{T}}$ is the vector of body forces acting
on the plate. In order to solve Eq. 43, an FE model with $282$ eight-noded
quadrilateral finite elements is used (Fig. 8).

In this example, the plate’s Young’s modulus is considered uncertain and
spatially variable. We assign a homogeneous random field prior with log-normal
marginal distributions with mean $\mu_{E}=2\cdot
10^{4}~{}\mathrm{kN}~{}\mathrm{cm}^{-2}$ and standard deviation
$\sigma_{E}=3\cdot 10^{3}~{}\mathrm{kN}~{}\mathrm{cm}^{-2}$. The mean and
standard deviation of the underlying Gaussian field $\ln E(x_{1},x_{2})$
follow as $\mu_{\ln E}=9.89$ and $\sigma_{\ln E}=0.15$, respectively, and its
correlation structure is modelled with the exponential kernel of Eq. 36 and
correlation length $l_{\ln E}=10~{}\mathrm{cm}$.

We discretize $\ln E$ by means of a Karhunen- Loève-expansion (KL-expansion).
To this end, we solve the following homogeneous Fredholm integral equation of
the second kind [21]

(44) $\sigma_{\ln E}^{2}\int_{D}\rho(\bm{x},\bm{x}^{\prime};l_{\ln
E})\phi_{k}(\bm{x}^{\prime})\mathrm{d}\bm{x}^{\prime}=\lambda_{k}\phi_{k}(\bm{x})$

for the covariance kernel’s set of eigenpairs $\\{\lambda_{k},\phi_{k}\\}$.
Consequently, we can express the log-normal Young’s modulus prior as the
$\exp$ of a KL-expansion [21] like

(45)
$E(x_{1},x_{2};\bm{\theta})=\exp{\left[\mu_{E}+\sum_{k=1}^{\infty}\sqrt{\lambda_{k}}\phi_{k}(x_{1},x_{2})\theta_{k}\right]},$

where the coefficients $\theta_{k}$ are pairwise independent standard-normal
Gaussian random variables.

We estimate the set of eigenpairs $\\{\lambda_{k},\phi_{k}\\}$ for the KL-
expansion by solving Eq. 44 using the Nyström method on a grid of $160\times
160$ Gauss-Legendre quadrature points. The eigenfunctions are interpolated at
the numerical integration points of the elements of the FE-model [44].

Truncating the KL-expansion Eq. 45 after $M$ terms results in an $M$-order KL-
approximation of $E$, which we denote as $\hat{E}(\bm{x};\bm{\theta})$. This
approximation recovers the random field mean exactly, however is associated
with an under-representation of its variance $\ln\sigma_{E}^{2}$. This under-
representation is often measured with the global relative variance error of
the $M$-order KL-approximation:

(46)
$\bar{\epsilon}_{\ln\sigma^{2}}=\frac{1}{|D|}\int_{D}\left|\frac{\mathbb{V}[E(\bm{x};\bm{\theta})]-\mathbb{V}[\hat{E}(\bm{x};\bm{\theta})]}{\mathbb{V}[E(\bm{x};\bm{\theta})]}\right|\mathrm{d}\bm{x}=1-\frac{1}{|D|\cdot\ln\sigma_{E}^{2}}\sum_{k=1}^{M}\lambda_{k}.$

Therein, $D$ is the spatial domain of the random field $E$. The inference task
for the example consists in learning the Young’s modulus’ posterior
distribution based on strain measurements at $n_{\mathrm{obs}}=10$ positions
on the plate. At each position, two gauges measure the strain in $x_{1}$\- and
$x_{2}$-direction (Fig. 8, left, red squares), respectively. Hence, a set of
$20$ measurements are available to solve the inference task. We generate the
measurements artificially by using a FE-model on a finer mesh of 779 elements,
in order to once again avoid the ’inverse crime’ [29]. The true strains are
depicted in the center and right plot in Fig. 8.

The strain measurements are generated by solving the forward problem based on
a single Young’s modulus prior random field realization. This realization is
generated using a midpoint method discretized at the numerical integration
(Gauss) points of the plate FE model rather than a KL-approximation.
Consequently, noise is added to the computed strains at the measurement
locations. We model the noise as a centered Gaussian random vector
$\bm{\eta}\sim\mathcal{N}(\bm{0},\bm{\Sigma}_{\eta\eta})$. The noise standard
deviation is set to $\sigma_{\eta}=10^{-4}$ and the autocorrelation of both
$x_{1}$\- and $x_{2}$-strain measurements is modelled with the exponential
kernel Eq. 36 using a correlation length of $l_{\eta}=10~{}\mathrm{cm}$. The
cross-correlation function between $x_{1}$\- and $x_{2}$-strain measurements
is taken as the autocorrelation function multiplied by a cross-correlation
coefficient of $0.25$.

####  4.2.2 Numerical reference posterior

We use adaptive Bayesian Updating with Subset Simulation (aBUS-SuS) to verify
the solution obtained with CEBUred. aBUS-SuS has been tested on a variety of
engineering applications, e.g., in [52, 53, 4, 27]. aBUS-SuS recasts the
Bayesian inverse problem as a structural reliability problem [57]. Subset
simulation (SuS) [2] is a robust and efficient method for solving such
structural reliability problems and within aBUS-SuS, SuS is employed to solve
general Bayesian inverse problems. SuS itself requires carrying out an MCMC
sampling step for which we use a pCN sampler [8] with adaptive scaling [42].

####  4.2.3 Parameters of numerical study

We discretize the Young’s modulus random field by a KL-approximation Eq. 45
with 879 terms producing a Bayesian inverse problem with ambient dimension
$d=879$. The chosen number of terms accounts for at least $97\%$ of the
spatial variance of the random field meaning the average variance error Eq. 46
is $\leq 3\%$.

The inference task is solved by using CEBUred with different error tolerances
$\epsilon$. We choose $\epsilon=\\{1.0,10^{-1},10^{-2},10^{-3}\\}$. The
remaining parameters are set as $\delta_{v,target}=1.5$,
$\alpha_{\mathrm{par}}=3$, and $\alpha_{\mathbf{H}}=4$. The gradients of the
likelihood function required at each step of CEBUred are evaluated with the
adjoint method [1] (derived for this particular problem in [58, Appendix A]).
In the final step, we draw $N=10000$ samples from the approximate posterior
distribution. The parameters for the numerical reference posterior generated
with aBUS-SuS are $n=20000$ for both samples per subset level and final
samples of the approximated posterior and intermediate conditional probability
$p=0.1$. Except for the reference posterior, which is computed once only, we
repeat each analysis with CEBUred $40$ times and average all results over the
individual runs.

####  4.2.4 Discussion of results

Fig. 9 compares the posterior fields obtained with CEBUred using an error
threshold of $\epsilon=10^{-2}$ and the reference posterior obtained with
aBUS-SuS along six sections across the plate. Along each section, the CEBUred-
based posterior means are in good agreement with the reference posterior mean.
For all other tested error thresholds ($\epsilon\in\\{1,10^{-1},10^{-2}\\}$),
similar results are obtained for the average posterior means and variances
taken over 40 repeated CEBUred runs (see Table 5). The coefficients of
variation of the posterior mean and variance estimates are rather large for
$\epsilon\geq 10^{-1}$, but decrease significantly between $\epsilon=10^{-1}$
and $\epsilon=10^{-2}$ based on the results given in Table 5.

The numbers of required likelihood and likelihood gradient evaluations for
aBUS-SuS and CEBUred at the four tested error threshold are listed in Table 4.
At $\epsilon=10^{-2}$, CEBUred reduces the number of required likelihood calls
by roughly two orders of magnitude compared to the reference aBUS-SuS run,
which comes at the cost of 642 additional gradient calls.

Figure 9: Posterior Young’s modulus random field obtained with aBUS-SuS and
with CEBUred using $\epsilon=10^{-2}$: Means (solid lines) and 95% credible
intervals (CI, dashed lines) at three vertical (top row) and three horizontal
(bottom row) cross-sections.

Fig. 11 shows the mean ranks and the corresponding number of model and
gradient evaluations. To ensure an accurate construction of the LIS, the
eigenvectors corresponding to the $r$ largest eigenvalues of $\mathbf{H}$ must
be reasonably well estimated. The number of samples required for the
estimation is determined with the heuristic formula given by Eq. 31. According
to this formula, the number of samples for the estimation of $\mathbf{H}$
linearly depends on $r$. Therefore, the lines in the left plot and the dashed
lines in the right plot in Fig. 11 are linearly dependent. For $\epsilon\geq
10^{-1}$, we obtain LIS dimensions of $r\leq 5$ in all steps.

Table 4: Likelihood and gradient evaluations for the plate example per run of aBUS-SuS and CEBUred (CEBUred: averaged over 40 runs). | aBUS-SUS | CEBUred | CEBUred | CEBUred | CEBUred
---|---|---|---|---|---
|  |  ($\epsilon=1$) |  ($\epsilon=10^{-1}$) |  ($\epsilon=10^{-2}$) |  ($\epsilon=10^{-3}$)
Llikelihood calls | 120000 | 76 | 276 | 1313 | 10621
Gradient calls | - | 105.3 | 254.5 | 642.0 | 1975.1
Table 5: Relative posterior mean and variance error at varying $\epsilon$ averaged over 40 runs with coefficient of variation in round brackets error | $\epsilon=1$ | $\epsilon=10^{-1}$ | $\epsilon=10^{-2}$ | $\epsilon=10^{-3}$
---|---|---|---|---
$\varepsilon_{\mu_{E}}$ (CoV) | $0.017(1.66)$ | $0.019(0.98)$ | $0.003(0.32)$ | $0.002(0.17)$
$\varepsilon_{\sigma^{2}_{E}}$ (CoV) | $0.174(0.52)$ | $0.201(0.75)$ | $0.039(0.33)$ | $0.027(0.10)$

Fig. 11 (right) depicts the number of gradient and model evaluations per
CEBUred step for varying $\epsilon$ exposing that for $\epsilon=1$, the number
of gradient evaluations performed to estimate the LIS exceeds the overall
number model evaluations required to perform the parameter update of the
biasing density. This is due to the fact that the number of required model
evaluations quadratically depends on the LIS dimension $r$ through Eq. 32,
where in turn $r$ decreases with increasing $\epsilon$.

Figure 10: Left: Average number of selected LIS dimension with increasing
CEBUred step index (95% CI indicated as shaded area). Right: Average number of
evaluated samples per step with increasing CEBUred step index. The solid and
the dashed lines represent the total number of model and model gradient
evaluations, respectively.

Figure 11: Mean of the normalized effective sample size (nESS) (95% CI
indicated as shaded area) plotted over the tested error thresholds along with
the target nESS (solid red line).

Fig. 11 displays the nESS of the final posterior sample along with the target
nESS, $n_{\mathrm{eff}}^{*}=0.3077$, that is related to the target coefficient
of variation of the IS weights $\delta_{w}=1.5$ through Eq. 10. The
simulations with $\epsilon=1.0$ and $\epsilon=10^{-1}$ exhibit small nESS well
below the target. In these cases, the dimensionality of the LIS is too low to
accurately represent the posterior of the random field and, consequently, the
IS weights have large variance. This leads to larger posterior mean and
variance errors and larger associated coefficients of variation of these error
measures for $\epsilon=1.0$ and $\epsilon=10^{-1}$ as documented in Table 5.

Compared to the beam application, where $n_{\mathrm{eff}}^{*}$ was achieved
for all investigated ambient dimensions with $\epsilon=1.0$, the plate
obviously requires a stricter error threshold. According to Table 5 and Fig.
11, $\epsilon=10^{-2}$ is a good choice for the present example. However, if
the threshold is chosen very small, e.g., $\epsilon=10^{-3}$, no significant
improvement is observed. In this case, the marginally increased accuracy will
not justify the additional computational expenses incurred by reducing the
threshold.

Table 6: Number of runs broken down according to required number of CEBUred steps at varying $\epsilon$. # of steps | $\epsilon=1$ | $\epsilon=10^{-1}$ | $\epsilon=10^{-2}$ | $\epsilon=10^{-3}$
---|---|---|---|---
2 | $3$ | $0$ | $0$ | $0$
3 | $30$ | $39$ | $40$ | $40$
4 | $7$ | $1$ | $0$ | $0$
Figure 12: Top: $\mathbf{H}$-spectra plotted versus CEBUred step number for
each error threshold $\epsilon$. Bottom: Upper bound of the KLD between the
tempered, full posterior and the tempered, optimally reduced posterior as
defined by Eq. 20. For each $\epsilon$ and CEBUred step, the largest scatter
point corresponds to $r=1$, the second-largest corresponds to $r=2$ and so on
to $r=d-1$. Solid lines indicate the four tested values of the error threshold
$\epsilon$.

Fig. 12 (top) depicts the spectrum of $\mathbf{H}$ at each step of CEBUred and
each predefined error threshold. As shown in Table 6, all runs with
$\epsilon\leq 10^{-2}$ required three steps, whereas only few runs with
$\epsilon\geq 10^{-1}$ required four steps. For this reason, no spectra appear
in the rightmost panels of Fig. 12 at step 4 for $\epsilon\leq 10^{-2}$. All
displayed spectra share two dominant eigenvalues of comparable magnitude
followed by a sharp decay. The main difference amongst spectra associated with
different $\epsilon$ is the number of samples used to estimate $\mathbf{H}$,
$n_{\mathbf{H}}$. While $n_{\mathbf{H}}$ has no influence on the dominant
eigenvalues, it bears some effect on a gap in the center of the spectrum, that
becomes narrower and eventually closes as $n_{\mathbf{H}}$ increases. This
effect, however, is negligible for the LIS construction as it takes place at
eigenvalue magnitudes well below that of the smallest eigenvalue whose
corresponding eigenvector is included in the LIS for any choice of $\epsilon$.

In Fig. 12 (bottom), we show the criterion Eq. 20 for each possible choice of
$r=1\dots d-1$ along with the tested error thresholds
$\epsilon=\\{1.0,10^{-1},10^{-2},10^{-3}\\}$. At any given step $t$, we can
directly compare the scatter points belonging to different error thresholds
since they have equal $\beta_{t}$ on average. The number of scatter points
that lie above the horizontal solid lines indicating $\epsilon$ determine the
LIS dimension $r$. The upper scatter points are approximately equal for any
choice of $\epsilon$ and any $t$ as the corresponding summations are dominated
by their leading term.

Using the heuristic given in Eq. 31 leads to a good estimate of the desired
first eigenvalues and vectors. This is important because the accuracy of the
LIS depends on these eigenvectors. However, for the computation of the upper
bound of the KLD between the full posterior and the optimal reduced posterior
Eq. 20 and the subsequent determination of the rank, all eigenvalues are
needed. For problems with rapidly decaying eigenvalue spectra of their
$\mathbf{H}$ matrix, as we see in this example, this is not a concern in
practice, since their smallest eigenvalues have little effect on the
computation of Eq. 20.

##  5 Concluding remarks

We present CEBUred (Cross-Entropy-based IS method for Bayesian Updating in
reduced space), an algorithm for approximating posterior distributions that
are the solutions of nonlinear Bayesian inverse problems. Such problems often
arise in the context of finding inverse solutions to computationally expensive
numerical models and solvers. Thus, computational efficiency is of the
essence, which translates to minimizing the number of required samples
(evaluations of the numerical model) to approximate the sought posterior
distribution at a prescribed accuracy. We address high-dimensional problem
settings that arise, e.g., if the inference target is represented by random
fields or processes, by identifying low-dimensional linear subspaces [60] in
which we perform cross-entropy-based importance sampling [15]. These subspaces
are obtained as truncated eigenspaces of the second-moment matrix of the
gradient of the log-likelihood $\mathbf{H}$.

We investigate CEBUred using two benchmark problems from engineering
mechanics. In the first example, the material parameter random field of a
cantilever beam subject to a point load is inferred from noise-distorted
deflection measurements. We examine the performance of CEBU versus CEBUred
versus a known analytical posterior reference solution at varying dimension of
the material parameter random field discretization. We find that the
dimensionality reduction is vital to ensure the posterior approximation
accuracy is independent of the problem dimension by comparing CEBU and
CEBUred. The second problem consists of inferring the material parameter
random field of a clamped steel plate under load from a strain measurement at
10 locations on the plate. The inference problem is set in a $879$-dimensional
space as the parameter random field is discretized with a $879$-term KLE. We
find that the quality of the posterior approximation produced by CEBUred is
closely connected to the choice of the error threshold $\epsilon$ that
controls the number of dimensions retained in reduced space. As $\epsilon$
decreases, CEBUred is able to recover the reference posterior solution
accurately both in mean and credible intervals. We further compare the
performance of CEBUred to that of aBUS-SUS (adaptive Bayesian updating with
subset simulation, [4]), which is a well-established method for nonlinear BIPs
in high-dimensions.

The results of our numerical investigations show that CEBUred is a powerful
method for solving high-dimensional nonlinear BIPs if the underlying
computational model admits a low-dimensional representation, i.e., if the
spectrum of $\mathbf{H}$ exhibits fast decay. From a computational
perspective, CEBUred is particularly useful if the model allows for the cheap
evaluation of gradients, e.g., if an adjoint solver is used and the the number
of BIP inputs exceeds the number of available observations. In these cases,
CEBUred achieves the same accuracy as aBUS-SuS at considerably lower
computational expense.

Conversely, if an adjoint solver is not available or not efficient in the
sense described above, the gain in computational efficiency provided by
dimensionality reduction may be overcompensated by expensive gradient
evaluations. In such case, the number of required gradient evaluations could
be significantly reduced by using ’data-free likelihood-informed dimension
reduction‘ as recently proposed in [10]. There, $\mathbf{H}$ is constructed in
expectation over the data such that no knowledge about the posterior density
is needed. Consequently, the upper bound between the exact and approximated
posterior is controlled in expectation over the data. Following this method,
the algorithm of CEBUred could be modified such that $\mathbf{H}$ is
constructed at the beginning and therefore, only at this stage model gradient
information would be required. Alternatively, one may turn to gradient-free
supervised dimension reduction methods in order to identify sutaible subspaces
to solve the Bayesian inverse problem.

## Acknowledgments

We acknowledge support by the German Research Foundation (DFG) through Grants
STR 1140/11-1 and PA 2901/1-1.

## References

  * [1] J. S. Arora and E. J. Haug, Methods of design sensitivity analysis in structural optimization, AIAA Journal, 17 (1979), pp. 970–974, https://doi.org/10.2514/3.61260.
  * [2] S.-K. Au and J. L. Beck, Estimation of small failure probabilities in high dimensions by subset simulation, Probabilistic engineering mechanics, 16 (2001), pp. 263–277.
  * [3] A. Beskos, G. Roberts, A. Stuart, and J. Voss, MCMC methods for diffusion bridges, Stochastics and Dynamics, 08 (2008), pp. 319–350, https://doi.org/10.1142/S0219493708002378.
  * [4] W. Betz, I. Papaioannou, J. L. Beck, and D. Straub, Bayesian inference with subset simulation: Strategies and improvements, Computer Methods in Applied Mechanics and Engineering, 331 (2018), pp. 72–93, https://doi.org/https://doi.org/10.1016/j.cma.2017.11.021.
  * [5] A. F. Bower, Applied mechanics of solids, CRC press, 2009.
  * [6] P. G. Constantine, Active Subspaces, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015, https://doi.org/10.1137/1.9781611973860, https://arxiv.org/abs/https://epubs.siam.org/doi/pdf/10.1137/1.9781611973860.
  * [7] P. G. Constantine, C. Kent, and T. Bui-Thanh, Accelerating markov chain monte carlo with active subspaces, SIAM Journal on Scientific Computing, 38 (2016), pp. A2779–A2805, https://doi.org/10.1137/15M1042127.
  * [8] S. L. Cotter, G. O. Roberts, A. M. Stuart, and D. White, MCMC Methods for Functions: Modifying Old Algorithms to Make Them Faster, Statistical Science, 28 (2013), pp. 424 – 446, https://doi.org/10.1214/13-STS421.
  * [9] T. Cui, J. Martin, Y. M. Marzouk, A. Solonen, and A. Spantini, Likelihood-informed dimension reduction for nonlinear inverse problems, Inverse Problems, 30 (2014), p. 114015, https://doi.org/10.1088/0266-5611/30/11/114015, https://doi.org/10.1088/0266-5611/30/11/114015.
  * [10] T. Cui and O. Zahm, Data-free likelihood-informed dimension reduction of bayesian inverse problems, Inverse Problems, 37 (2021), p. 045009, https://doi.org/10.1088/1361-6420/abeafb.
  * [11] P. Del Moral, A. Doucet, and A. Jasra, Sequential monte carlo samplers, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68 (2006), pp. 411–436, https://doi.org/https://doi.org/10.1111/j.1467-9868.2006.00553.x.
  * [12] A. Der Kiureghian and J.-B. Ke, The stochastic finite element method in structural reliability, Probabilistic Engineering Mechanics, 3 (1988), pp. 83–91, https://doi.org/https://doi.org/10.1016/0266-8920(88)90019-7, https://www.sciencedirect.com/science/article/pii/0266892088900197.
  * [13] A. Doucet, A. Smith, N. de Freitas, and N. Gordon, Sequential Monte Carlo Methods in Practice, Information Science and Statistics, Springer New York, 2001.
  * [14] V. Elvira, L. Martino, D. Luengo, and M. F. Bugallo, Novel weighting and resampling schemes in Population Monte Carlo, in XXVIème Colloque GRETSI, Juan-les-Pins, France, Sept. 2017, https://hal.archives-ouvertes.fr/hal-01684862.
  * [15] M. Engel, O. Kanjilal, I. Papaioannou, and D. Straub, Bayesian updating and marginal likelihood estimation by cross entropy based importance sampling.  Manuscript, 2021.
  * [16] H. P. Flath, L. C. Wilcox, V. Akçelik, J. Hill, B. van Bloemen Waanders, and O. Ghattas, Fast algorithms for bayesian uncertainty quantification in large-scale linear inverse problems based on low-rank partial hessian approximations, SIAM Journal on Scientific Computing, 33 (2011), pp. 407–432, https://doi.org/10.1137/090780717.
  * [17] A. Gelman, J. Carlin, H. Stern, D. Dunson, A. Vehtari, and D. Rubin, Bayesian Data Analysis, Third Edition, Chapman & Hall/CRC Texts in Statistical Science, Taylor & Francis, 2013.
  * [18] A. Gelman and Y. Yao, Holes in bayesian statistics, Journal of Physics G: Nuclear and Particle Physics, 48 (2020), p. 014002, https://doi.org/10.1088/1361-6471/abc3a5, https://doi.org/10.1088/1361-6471/abc3a5.
  * [19] J. Geweke, Bayesian inference in econometric models using monte carlo integration, Econometrica, 57 (1989), pp. 1317–39.
  * [20] S. Geyer, I. Papaioannou, and D. Straub, Cross entropy-based importance sampling using gaussian densities revisited, Structural Safety, 76 (2019), pp. 15–27.
  * [21] R. G. Ghanem and P. D. Spanos, Spectral stochastic finite-element formulation for reliability analysis, Journal of Engineering Mechanics, 117 (1991), pp. 2351–2372.
  * [22] W. Gilks, S. Richardson, and D. Spiegelhalter, Markov Chain Monte Carlo in Practice, Chapman & Hall/CRC Interdisciplinary Statistics, Taylor & Francis, 1 ed., 1995.
  * [23] A. Gittens and J. A. Tropp, Tail bounds for all eigenvalues of a sum of random matrices, 2011, https://arxiv.org/abs/1104.4513.
  * [24] A. Jasra, D. A. Stephens, A. Doucet, and T. Tsagarakis, Inference for Lévy-driven stochastic volatility models via adaptive sequential monte carlo, Scandinavian Journal of Statistics, 38 (2011), pp. 1–22.
  * [25] E. Jaynes, On the rationale of maximum-entropy methods, Proceedings of the IEEE, 70 (1982), pp. 939–952, https://doi.org/10.1109/PROC.1982.12425.
  * [26] H. Jeffreys, Theory of Probability, Oxford, Oxford, England, third ed., 1961.
  * [27] S.-H. Jiang, I. Papaioannou, and D. Straub, Bayesian updating of slope reliability in spatially variable soils with in-situ measurements, Engineering Geology, 239 (2018), pp. 310 – 320.
  * [28] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Dover Books on Mathematics Series, Dover Publications, Incorporated, 2012, https://books.google.de/books?id=PYXjyoqy5qMC.
  * [29] J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Springer, 2005, https://doi.org/10.1007/b138659.
  * [30] P. Koutsourelakis, A multi-resolution, non-parametric, bayesian framework for identification of spatially-varying model parameters, Journal of Computational Physics, 228 (2009), pp. 6184–6211, https://doi.org/https://doi.org/10.1016/j.jcp.2009.05.016, https://www.sciencedirect.com/science/article/pii/S0021999109002708.
  * [31] N. Kurtz and J. Song, Cross-entropy-based adaptive importance sampling using gaussian mixture, Structural Safety, 42 (2013), pp. 35 – 44.
  * [32] J. Latz, I. Papaioannou, and E. Ullmann, Multilevel sequential2 monte carlo for bayesian inverse problems, Journal of Computational Physics, 368 (2018), pp. 154–178, https://doi.org/https://doi.org/10.1016/j.jcp.2018.04.014, https://www.sciencedirect.com/science/article/pii/S0021999118302286.
  * [33] P. Liu and K. Liu, Selection of random field mesh in finite element reliability analysis, Journal of Engineering Mechanics, 119 (1993), pp. 667–680, https://doi.org/10.1061/(ASCE)0733-9399(1993)119:4(667).
  * [34] P.-L. Liu and A. Der Kiureghian, Multivariate distribution models with prescribed marginals and covariances, Probabilistic Engineering Mechanics, 1 (1986), pp. 105–112, https://doi.org/https://doi.org/10.1016/0266-8920(86)90033-0, https://www.sciencedirect.com/science/article/pii/0266892086900330.
  * [35] Y. M. Marzouk and H. N. Najm, Dimensionality reduction and polynomial chaos acceleration of bayesian inference in inverse problems, Journal of Computational Physics, 228 (2009), pp. 1862–1902, https://doi.org/https://doi.org/10.1016/j.jcp.2008.11.024.
  * [36] J. C. Mattingly, N. S. Pillai, and A. M. Stuart, Diffusion limits of the random walk Metropolis algorithm in high dimensions, The Annals of Applied Probability, 22 (2012), pp. 881 – 930, https://doi.org/10.1214/10-AAP754.
  * [37] R. M. Neal, Annealed importance sampling, Statistics and Computing, 11 (2001), p. 125–139, https://doi.org/10.1023/A:1008923215028, https://doi.org/10.1023/A:1008923215028.
  * [38] R. M. Neal, MCMC using Hamiltonian dynamics, in Handbook of Markov Chain Monte Carlo, Chapman; Hall/CRC, 2010, pp. 113–162.
  * [39] R. B. Nelsen, An Introduction to Copulas, Springer Publishing Company, Incorporated, 2010.
  * [40] A. B. Owen, Monte Carlo theory, methods and examples, 2013.
  * [41] A. O’Hagan, Expert knowledge elicitation: Subjective but scientific, The American Statistician, 73 (2019), pp. 69–81, https://doi.org/10.1080/00031305.2018.1518265, https://doi.org/10.1080/00031305.2018.1518265, https://arxiv.org/abs/https://doi.org/10.1080/00031305.2018.1518265.
  * [42] I. Papaioannou, W. Betz, K. Zwirglmaier, and D. Straub, Mcmc algorithms for subset simulation, Probabilistic Engineering Mechanics, 41 (2015), pp. 89–103.
  * [43] I. Papaioannou, S. Geyer, and D. Straub, Improved cross entropy-based importance sampling with a flexible mixture model, Reliability Engineering & System Safety, 191 (2019), p. 106564.
  * [44] W. H. Press, W. T. Vetterling, S. A. Teukolsky, and B. P. Flannery, Numerical recipes, vol. 818, Cambridge university press Cambridge, 1986.
  * [45] C. Robert, The Bayesian choice: from decision-theoretic foundations to computational implementation, Springer Science & Business Media, 2007.
  * [46] G. O. Roberts and J. S. Rosenthal, Optimal scaling for various Metropolis-Hastings algorithms, Statistical Science, 16 (2001), pp. 351 – 367, https://doi.org/10.1214/ss/1015346320.
  * [47] M. Rosenblatt, Remarks on a multivariate transformation, The annals of mathematical statistics, 23 (1952), pp. 470–472.
  * [48] R. Y. Rubinstein, Optimization of computer simulation models with rare events, European Journal of Operational Research, 99 (1997), pp. 89–112.
  * [49] R. Y. Rubinstein and P. W. Glynn, How to deal with the curse of dimensionality of likelihood ratios in monte carlo simulation, Stochastic Models, 25 (2009), pp. 547–568, https://doi.org/10.1080/15326340903291248.
  * [50] R. Y. Rubinstein and D. P. Kroese, Simulation and the Monte Carlo method, Springer, 3 ed., 2017.
  * [51] A. Spantini, A. Solonen, T. Cui, J. Martin, L. Tenorio, and Y. Marzouk, Optimal low-rank approximations of bayesian linear inverse problems, SIAM Journal on Scientific Computing, 37 (2015), pp. A2451–A2487, https://doi.org/10.1137/140977308.
  * [52] D. Straub and I. Papaioannou, Bayesian updating with structural reliability methods, Journal of Engineering Mechanics, 141 (2015), https://doi.org/10.1061/(ASCE)EM.1943-7889.0000839, http://ascelibrary.org/doi/abs/10.1061/%28ASCE%29EM.1943-7889.0000839.
  * [53] D. Straub, I. Papaioannou, and W. Betz, Bayesian analysis of rare events, Journal of Computational Physics, 314 (2016), pp. 538–556.
  * [54] A. M. Stuart, Inverse problems: A bayesian perspective, Acta Numerica, 19 (2010), p. 451–559, https://doi.org/10.1017/S0962492910000061.
  * [55] E. Torre, S. Marelli, P. Embrechts, and B. Sudret, A general framework for data-driven uncertainty quantification under complex input dependencies using vine copulas, Probabilistic Engineering Mechanics, 55 (2019), pp. 1–16, https://doi.org/https://doi.org/10.1016/j.probengmech.2018.08.001.
  * [56] J. A. Tropp, User-friendly tail bounds for sums of random matrices, Foundations of Computational Mathematics, 12 (2012), pp. 389–434, https://doi.org/10.1007/s10208-011-9099-z.
  * [57] F. Uribe, I. Papaioannou, W. Betz, and D. Straub, Bayesian inference of random fields represented with the Karhunen-Loève expansion, Computer Methods in Applied Mechanics and Engineering, 358 (2020), p. 112632, https://doi.org/https://doi.org/10.1016/j.cma.2019.112632.
  * [58] F. Uribe, I. Papaioannou, Y. M. Marzouk, and D. Straub, Cross-entropy-based importance sampling with failure-informed dimension reduction for rare event simulation, SIAM/ASA Journal on Uncertainty Quantification, 9 (2021), pp. 818–847, https://doi.org/10.1137/20M1344585.
  * [59] Z. Wang and J. Song, Cross-entropy-based adaptive importance sampling using von mises-fisher mixture for high dimensional reliability analysis, Structural Safety, 59 (2016), pp. 42–52.
  * [60] O. Zahm, T. Cui, K. Law, A. Spantini, and Y. Marzouk, Certified dimension reduction in nonlinear bayesian inverse problems, 2021, https://arxiv.org/abs/1807.03712.