Bayesian model inversion using stochastic spectral embedding

Consider some non-observable parameters $\bm{X}\in{\mathcal{D}}_{\bm{X}}$ and the observables $\bm{Y}\in{\mathcal{D}}_{\bm{Y}}$. Furthermore, let ${{\mathcal{Y}}=\{\bm{y}^{(1)},\dots,\bm{y}^{(N)}\}}$ be a set of $N$ measurements, i.e., noisy observations of a set of realizations of $\bm{Y}$. Statistical inference consists in drawing conclusions about $\bm{X}$ using the information from ${\mathcal{Y}}$ (Gelman et al., 2014). These measurements can be direct observations of the parameters ($\bm{Y}=\bm{X}$) or some quantities indirectly related to $\bm{X}$ through a function or model ${\mathcal{M}}:{\mathcal{D}}_{\bm{X}}\to{\mathcal{D}}_{\bm{Y}}$. One way to conduct this inference is through Bayes’ theorem of conditional probabilities, a process known as Bayesian inference.

Denoting by $\pi(\cdot)$ a probability density function (PDF) and by $\pi(\cdot|\bm{x})$ a PDF conditioned on $\bm{x}$, Bayes’ theorem can be written as

where $\pi(\bm{x})$ is known as the prior distribution of the parameters, i.e., the distribution of $\bm{X}$ before observing the data ${\mathcal{Y}}$. The conditional distribution $\pi({\mathcal{Y}}|\bm{x})$, known as likelihood, establishes a connection between the observations ${\mathcal{Y}}$ and a realization of the parameters $\bm{X}=\bm{x}$. For a given realization $\bm{x}$, it returns the probability density of observing the data ${\mathcal{Y}}$. Under the common assumption of independence between individual observations, $\{\bm{y}^{(i)}\}_{i=1}^{N}$, the likelihood function takes the form:

The likelihood function is a map ${\mathcal{D}}_{\bm{X}}\rightarrow\mathbb{R}_{+}$, and it attains its maximum for the parameter set with the highest probability of yielding ${\mathcal{Y}}$. With this, Bayes’ theorem from Eq. (1) can be rewritten as:

where $Z$ is a normalizing constant often called evidence or marginal likelihood. On the left-hand side, $\pi(\bm{x}|{\mathcal{Y}})$ is the posterior PDF, i.e., the distribution of $\bm{X}$ after observing data ${\mathcal{Y}}$. In this sense, Bayes’ theorem establishes a general expression for updating the prior distribution using a likelihood function to incorporate information from the data.

Bayesian model inversion describes the application of the Bayesian inference framework to the problem of model inversion (Beck and Katafygiotis, 1998; Kennedy and O’Hagan, 2001; Jaynes, 2003; Tarantola, 2005). The two main ingredients needed to infer model parameters within the Bayesian framework are a prior distribution $\pi(\bm{x})$ of the model parameters and a likelihood function ${\mathcal{L}}$. In practical applications, prior information about the model parameters is often readily available. Typical sources of such information are physical parameter constraints or expert knowledge. Additionally, prior inversion attempts can serve as guidelines to assign informative prior distributions. In cases where no prior information about the parameters is available, so-called non-informative or invariant prior distributions (Jeffreys, 1946; Harney, 2016) can also be assigned. The likelihood function serves instead as the link between model parameters $\bm{X}$ and observations of the model output ${\mathcal{Y}}$. To connect these two quantities, it is necessary to choose a so-called discrepancy model that gives the relative probability that the model response to a realization of $\bm{X}=\bm{x}$ describes the observations. One common assumption for this probabilistic model is that the measurements are perturbed by a Gaussian additive discrepancy term $\bm{E}\sim{\mathcal{N}}(\bm{\varepsilon}|\bm{0},\bm{\Sigma})$, with covariance matrix $\bm{\Sigma}$. For a single measurement $\bm{y}^{(i)}$ it reads:

This discrepancy between the model output ${\mathcal{M}}(\bm{X})$ and the observables $\bm{Y}$ can result from measurement error or model inadequacies. By using this additive discrepancy model, the distribution of the observables conditioned on the parameters $\bm{Y}|\bm{x}$ is written as:

where ${\mathcal{N}}(\cdot|\bm{\mu},\bm{\Sigma})$ denotes the multivariate Gaussian PDF with mean value $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$. The likelihood function ${\mathcal{L}}$ is then constructed using this probabilistic model $\pi(\bm{y}^{(i)}|\bm{x})$ and Eq. (2). For a given set of measurements ${\mathcal{Y}}$ it thus reads:

With the fully specified Bayesian model inversion problem, Eq. (3) directly gives the posterior distribution of the model parameters $\pi(\bm{x}|{\mathcal{Y}})$. In the setting of model inversion, the posterior distribution represents therefore the state of belief about the true data-generating model parameters, considering all available information: computational forward model, discrepancy model and measurement data (Beck and Katafygiotis, 1998; Jaynes, 2003).

Often, the ultimate goal of model inversion is to provide a set of inferred parameters, with associate confidence measures/intervals. This is often achieved by computing posterior statistics (e.g., moments, mode, etc.). Propagating the posterior through secondary models is also of interest. So-called quantites of interest (QoI) can be expressed by calculating the posterior expectation of suitable functions of the parameters $h(\bm{x}):\mathbb{R}^{M}\to\mathbb{R}$, with $\bm{X}|{\mathcal{Y}}\sim\pi(\bm{x}|{\mathcal{Y}})$, as in:

Depending on $h$, this formulation encompasses posterior moments ($h(\bm{x})=x_{i}$ or $h(\bm{x})=(x_{i}-{\mathbb{E}}\left[X_{i}\right])^{2}$ for the first and second moments, respectively), posterior covariance ($h(\bm{x})=x_{i}x_{j}-{\mathbb{E}}\left[X_{i}\right]{\mathbb{E}}\left[X_{j}\right]$) or expectations of secondary models ($h(\bm{x})={\mathcal{M}}^{\star}(\bm{x})$).

The idea of SLE is to use sparse PCE to find a spectral representation of the likelihood function ${\mathcal{L}}$ occurring in Bayesian model inversion problems (see Eq. (2)). We present here a brief introduction to the method and the main results of Nagel and Sudret (2016).

Likelihood functions can be seen as scalar functions of the input random vector $\bm{X}\sim\pi(\bm{x})$. In this work we assume priors of the type $\pi(\bm{x})=\prod_{i=1}^{M}\pi_{i}(x_{i})$, i.e. with independent marginals, their spectral expansion then reads:

where the explicit dependence on ${\mathcal{Y}}$ was dropped for notational simplicity.

Upon computing the basis and coefficients in Eq. (8), the solution to the inverse problem is converted to merely post-processing the coefficients $a_{\bm{\alpha}}$. The following expressions can be derived for the individual quantities:

Evidence

The evidence emerges as the coefficient of the constant polynomial $a_{\bm{0}}$

$$Z=\int_{{\mathcal{D}}_{\bm{X}}}{\mathcal{L}}(\bm{x})\pi(\bm{x})\,{\rm d}\bm{x}\approx\left<{\mathcal{L}}_{\mathrm{SLE}},1\right>_{\pi}=a_{\bm{0}}.$$

(11)

Posterior

Upon computing the evidence $Z$, the posterior can be evaluated directly through

$$\pi(\bm{x}|{\mathcal{Y}})\approx\frac{{\mathcal{L}}_{\mathrm{SLE}}(\bm{x})\pi(\bm{x})}{Z}=\frac{\pi(\bm{x})}{a_{\bm{0}}}\sum_{\bm{\alpha}\in{\mathcal{A}}}a_{\bm{\alpha}}\Psi_{\bm{\alpha}}(\bm{x}).$$

(12)

Posterior marginals

We split the random vector $\bm{X}$ into two vectors $\bm{X}_{{\bm{\mathsf{u}}}}$ with components $\{X_{i}\}_{i\in{\bm{\mathsf{u}}}}\in\mathcal{D}_{\bm{X}_{{\bm{\mathsf{u}}}}}$ and $\bm{X}_{{\bm{\mathsf{v}}}}$ with components $\{X_{i}\}_{i\in{\bm{\mathsf{v}}}}\in\mathcal{D}_{\bm{X}_{{\bm{\mathsf{v}}}}}$, where ${\bm{\mathsf{u}}}$ and ${\bm{\mathsf{v}}}$ are two non-empty disjoint index sets such that ${\bm{\mathsf{u}}}\cup{\bm{\mathsf{v}}}=\{1,\dots,M\}$. Denote further by $\pi_{{\bm{\mathsf{u}}}}(\bm{x}_{{\bm{\mathsf{u}}}})\stackrel{{\scriptstyle\text{def}}}{{=}}\prod_{i\in{\bm{\mathsf{u}}}}\pi_{i}(x_{i})$ and $\pi_{{\bm{\mathsf{v}}}}(\bm{x}_{{\bm{\mathsf{v}}}})\stackrel{{\scriptstyle\text{def}}}{{=}}\prod_{i\in{\bm{\mathsf{v}}}}\pi_{i}(x_{i})$ the prior marginal density functions of $\bm{X}_{{\bm{\mathsf{u}}}}$ and $\bm{X}_{{\bm{\mathsf{v}}}}$ respectively. The posterior marginals then read:

$$\pi_{{\bm{\mathsf{u}}}}(\bm{x}_{{\bm{\mathsf{u}}}}|{\mathcal{Y}})=\int_{{\mathcal{D}}_{\bm{X}_{{\bm{\mathsf{v}}}}}}\pi(\bm{x}|{\mathcal{Y}})\,{\rm d}\bm{x}_{{\bm{\mathsf{v}}}}\approx\frac{\pi_{{\bm{\mathsf{u}}}}(\bm{x}_{{\bm{\mathsf{u}}}})}{a_{\bm{0}}}\sum_{\bm{\alpha}\in{\mathcal{A}}_{{\bm{\mathsf{v}}}=0}}a_{\bm{\alpha}}\Psi_{\bm{\alpha}}(\bm{x}_{{\bm{\mathsf{u}}}}),$$

(13)

where ${\mathcal{A}}_{{\bm{\mathsf{v}}}=0}=\{\bm{\alpha}\in\mathcal{A}:\alpha_{i}=0\Leftrightarrow i\in{\bm{\mathsf{v}}}\}$. The series in the above equation constitutes a subexpansion that contains non-constant polynomials only in the directions $i\in{\bm{\mathsf{u}}}$.

Quantities of interest

Finally, it is also possible to analytically compute posterior expectations of functions that admit a polynomial chaos expansion in the same basis of the form $h(\bm{X})\approx\sum_{\bm{\alpha}\in{\mathcal{A}}}b_{\bm{\alpha}}\Psi_{\bm{\alpha}}(\bm{X})$. Eq. (7) then reduces to the spectral product:

$${\mathbb{E}}\left[h(\bm{X})|{\mathcal{Y}}\right]=\frac{1}{a_{\bm{0}}}\sum_{\bm{\alpha}\in{\mathcal{A}}}a_{\bm{\alpha}}b_{\bm{\alpha}}.$$

(14)

The quality of these results depends only on the approximation error introduced in Eq. (10). The latter, in turn, depends mainly on the chosen PCE truncation strategy (Blatman and Sudret, 2011; Nagel and Sudret, 2016) and the number of points used to compute the coefficients (i.e., the experimental design). It is known that likelihood functions typically have quasi-compact supports (i.e., ${\mathcal{L}}(\bm{X})\approx 0$ on a majority of ${\mathcal{D}}_{\bm{X}}$). Such functions require a very high polynomial degree to be approximated accurately, which in turn can lead to the need for prohibitively large experimental designs.

Viewing the likelihood as a function of a random variable $\bm{X}$ with independent marginals, we can directly use Eq. (15) to write down its SSLE representation

where the variable $\bm{X}$ is distributed according to the prior distribution $\pi(\bm{x})$ and, consequently, the local basis used to compute $\widehat{\mathcal{R}}_{S}^{k}(\bm{X})$ is orthonormal w.r.t. that distribution.

Due to the local spectral properties of the residual expansions, the SSLE representation of the likelihood function retains all of the post-processing properties of SLE (Section 3.1):

Evidence

The normalization constant $Z$ emerges as the sum of the constant polynomial coefficients weighted by the prior mass:

$$Z=\sum_{k\in{\mathcal{K}}}\sum_{\bm{\alpha}\in{\mathcal{A}}^{k}}a_{\bm{\alpha}}^{k}\int_{{\mathcal{D}}_{\bm{X}}^{k}}\Psi_{\bm{\alpha}}^{k}(\bm{x})\pi(\bm{x})\,{\rm d}\bm{x}=\sum_{k\in{\mathcal{K}}}\mathcal{V}^{k}a_{\bm{0}}^{k},\quad\text{where}\quad\mathcal{V}^{k}=\int_{{\mathcal{D}}_{\bm{X}}^{k}}\pi(\bm{x})\,{\rm d}\bm{x}.$$

(18)

Posterior

This allows us to write the posterior as

$$\pi(\bm{x}|{\mathcal{Y}})\approx\frac{{\mathcal{L}}_{\mathrm{SSLE}}(\bm{x})\pi(\bm{x})}{Z}=\frac{\pi(\bm{x})}{\sum_{k\in{\mathcal{K}}}\mathcal{V}^{k}a_{\bm{0}}^{k}}\sum_{k\in{\mathcal{K}}}\bm{1}_{\mathcal{D}_{\bm{X}}^{k}}(\bm{x})\widehat{\mathcal{R}}^{k}_{S}(\bm{x}).$$

(19)

Posterior marginal

Utilizing again the disjoint sets ${\bm{\mathsf{u}}}$ and ${\bm{\mathsf{v}}}$ from Eq. (13) it is also possible to analytically derive posterior marginal PDFs as

$$\pi_{{\bm{\mathsf{u}}}}(\bm{x}_{{\bm{\mathsf{u}}}}|{\mathcal{Y}})=\int_{{\mathcal{D}}_{\bm{X}_{{\bm{\mathsf{v}}}}}}\pi(\bm{x}|{\mathcal{Y}})\,{\rm d}\bm{x}_{{\bm{\mathsf{v}}}}\approx\frac{\pi_{{\bm{\mathsf{u}}}}(\bm{x}_{{\bm{\mathsf{u}}}})}{\sum_{k\in{\mathcal{K}}}\mathcal{V}^{k}a_{\bm{0}}^{k}}\sum_{k\in{\mathcal{K}}}\bm{1}_{{\mathcal{D}}_{\bm{X}_{{\bm{\mathsf{u}}}}}^{k}}(\bm{x}_{{\bm{\mathsf{u}}}})\widehat{\mathcal{R}}^{k}_{S,{\bm{\mathsf{u}}}}(\bm{x}_{{\bm{\mathsf{u}}}})\mathcal{V}^{k}_{{\bm{\mathsf{v}}}}$$

(20)

where

$$\widehat{\mathcal{R}}^{k}_{S,{\bm{\mathsf{u}}}}(\bm{x}_{{\bm{\mathsf{u}}}})=\sum_{\bm{\alpha}\in{\mathcal{A}}_{{\bm{\mathsf{v}}}=0}^{k}}a_{\bm{\alpha}}^{k}\Psi_{\bm{\alpha}}^{k}(\bm{x}_{{\bm{\mathsf{u}}}})\quad\text{and}\quad\mathcal{V}^{k}_{{\bm{\mathsf{v}}}}=\int_{{\mathcal{D}}_{\bm{X}_{{\bm{\mathsf{v}}}}}^{k}}\pi_{{\bm{\mathsf{v}}}}(\bm{x}_{{\bm{\mathsf{v}}}})\,{\rm d}\bm{x}_{{\bm{\mathsf{v}}}}.$$

(21)

$\widehat{\mathcal{R}}^{k}_{S,{\bm{\mathsf{u}}}}(\bm{x}_{{\bm{\mathsf{u}}}})$ is a subexpansion of $\widehat{\mathcal{R}}^{k}_{S}(\bm{x})$ that contains only non-constant polynomials in the directions $i\in{\bm{\mathsf{u}}}$. Note that, as we assumed that the prior distribution has independent components, the constants $\mathcal{V}^{k}$ and $\mathcal{V}^{k}_{{\bm{\mathsf{v}}}}$ are obtained as products of univariate integrals which are available analytically from the prior marginal cumulative distribution functions (CDFs).

Quantities of interest

Expected values of function $h(\bm{x})=\sum_{\bm{\alpha}\in{\mathcal{A}}^{k}}b_{\bm{\alpha}}^{k}\Psi_{\bm{\alpha}}^{k}(\bm{x})$ for $k\in{\mathcal{K}}$ under the posterior can be approximated by:

$\displaystyle\begin{split}\mathbb{E}\left[h(\bm{X})|{\mathcal{Y}}\right]&=\int_{{\mathcal{D}}_{\bm{X}}}h(\bm{x})\pi(\bm{x}|{\mathcal{Y}})\,{\rm d}\bm{x}\\ &=\frac{1}{Z}\sum_{k\in{\mathcal{K}}}\sum_{\bm{\alpha}\in{\mathcal{A}}^{k}}a_{\bm{\alpha}}^{k}\int_{{\mathcal{D}}_{\bm{X}}^{k}}h(\bm{x})\Psi_{\bm{\alpha}}^{k}(\bm{x})\pi(\bm{x})\,{\rm d}\bm{x}\\ &=\frac{1}{Z}\sum_{k\in{\mathcal{K}}}\sum_{\bm{\alpha}\in{\mathcal{A}}^{k}}a_{\bm{\alpha}}^{k}b_{\bm{\alpha}}^{k},\end{split}$

(22)

where $b_{\bm{\alpha}}^{k}$ are the coefficients of the PCE of $h$ in the $\mathrm{card}({\mathcal{K}})$ bases $\{\Psi_{\bm{\alpha}}^{k}\}_{\bm{\alpha}\in{\mathcal{A}}^{k}}$. This can also be used for computing posterior moments like mean, variance or covariance.

These expressions can be seen as a generalization of the ones for SLE detailed in Section 3.1. For a single-level global expansion (i.e., $\mathrm{card}({\mathcal{K}})=1$) and consequently $\mathcal{V}^{(0,1)}=1$, they are identical.

The original algorithm for computing an SSE was presented in Marelli et al. (2021). It recursively partitions the input domain $\mathcal{D}_{\bm{X}}$ and constructs truncated expansions of the residual. We reproduce it below for reference but replace the model $\mathcal{M}$ with the likelihood function $\mathcal{L}$. We further simplify the algorithm by choosing a partitioning strategy with $N_{S}=2$.

1. 1.

   Initialization:

   1. (a)

      $\ell=0$, $p=1$

   2. (b)

      $\mathcal{D}_{\bm{X}}^{\ell,p}=\mathcal{D}_{\bm{X}}$

   3. (c)

      ${\mathcal{R}}^{\ell}(\bm{X})=\mathcal{L}(\bm{X})$

2. 2.

   For each subdomain $\mathcal{D}_{\bm{X}}^{\ell,p}$, $p=1,\cdots,P_{\ell}$:

   1. (a)

      Calculate the truncated expansion $\widehat{\mathcal{R}}_{S}^{\ell,p}(\bm{X}^{\ell,p})$ of the residual ${\mathcal{R}}^{\ell}(\bm{X}^{\ell,p})$ in the current subdomain

   2. (b)

      Update the residual in the current subdomain ${\mathcal{R}}^{\ell+1}(\bm{X}^{\ell,p})={\mathcal{R}}^{\ell}(\bm{X}^{\ell,p})-\widehat{\mathcal{R}}_{S}^{\ell,p}(\bm{X}^{\ell,p})$

   3. (c)

      Split the current subdomain $\mathcal{D}_{\bm{X}}^{\ell,p}$ in $2$ subdomains $\mathcal{D}_{\bm{X}}^{\ell+1,\{s_{1},s_{2}\}}$ based on a partitioning strategy

   4. (d)

      If $\ell<L$, $\ell\leftarrow\ell+1$, go back to 2a, otherwise terminate the algorithm

3. 3.

   Termination

   1. (a)

      Return the full sequence of $\mathcal{D}_{\bm{X}}^{\ell,p}$ and $\widehat{\mathcal{R}}_{S}^{\ell,p}(\bm{X}^{\ell,p})$ needed to compute Eq. (15).

In practice, the residual expansions $\widehat{\mathcal{R}}_{S}^{\ell,p}(\bm{X}^{\ell,p})$ are computed using a fixed experimental design $\mathcal{X}$ and corresponding model evaluations $\mathcal{L}(\mathcal{X})$. The algorithm then only requires the specification of a partitioning strategy and a termination criterion, as detailed in Marelli et al. (2021).

Likelihood functions are typically characterized by a localized behaviour: Close to the data-generating parameters they peak while quickly decaying to $0$ in the remainder of the prior domain. This means that in a majority of the domain the likelihood evaluation is non-informative. Directly applying the original algorithm is then expected to waste many likelihood evaluations.

We therefore modify the original algorithm by adding an adaptive sampling scheme (Section 4.2.1) that includes the termination criterion and introducing an improved partitioning strategy (Section 4.2.2) that is especially suitable for finding compact support features. The rationale for these modifications is presented next.

In this first case study, the goal is the inference of a single unknown parameter with a multimodal posterior distribution. This problem is difficult to solve with standard MCMC methods due to the presence of a probability valley, i.e. low probability region, between essentially disjoint posterior peaks. It also serves as an illustrative example of how the proposed adaptive algorithm constructs an adSSLE. The problem is fabricated and uses synthetic data, but is presented in the context of a relevant engineering problem.

Consider the oscillator displayed in Figure 3 subject to harmonic (i.e., sinusoidal) excitation. Assume the prior information about its stiffness $X\stackrel{{\scriptstyle\text{def}}}{{=}}k$ is that it follows a lognormal distribution with $\mu=0.8~{}N/m$ and $\sigma=0.1~{}N/m$. Its true value shall be determined using measurements of the oscillation amplitude at the location of the mass $m$. The known properties of the oscillator system are the oscillator mass $m=1~{}kg$, the excitation frequency $\omega=1~{}rad/s$ and the viscous damping coefficient $c=0.1~{}Ns/m$. The oscillation amplitude is measured in five independent oscillation events and normalized by the forcing amplitude yielding the measured amplitude ratios ${\mathcal{Y}}=\{9.01,8.67,8.84,9.22,8.54\}$.

This problem is well known in mechanics and in the linear case (i.e., assuming small deformations and linear material behavior) can be solved analytically with the amplitude of the frequency response function. This function returns the ratio between the steady state amplitude of a linear oscillator and the amplitude of its excitation. It is given by

We assume a discrepancy model with known discrepancy standard deviation $\sigma$. In conjunction with the available measurements ${\mathcal{Y}}$, this leads to the following likelihood function:

We employ the adSSLE algorithm to approximate this likelihood function with $N_{\mathrm{ref}}=10$. A few iterations from the solution process are shown in Figure 4. The top plots show the subdomains ${\mathcal{D}}_{X}^{\ell,p}$ constructed at each refinement step, highlighting the terminal domains ${\mathcal{T}}$. The middle plots display the residual between the true likelihood and the approximation at the current iteration, as well as the adaptively chosen experimental design $\mathcal{X}$. The bottom plots display the target likelihood function and its current approximation.

The initial global approximation of the first iteration in Figure LABEL:sub@fig:cs1:SSEBehaviour:1 is a constant polynomial based on the initial experimental design. By the third iteration, the algorithm has identified the subdomain ${\mathcal{D}}_{X}^{2,2}$ as the one of interest and proceeds to refine it in subsequent steps. By the $8$th iteration both likelihood peaks have been identified. Finally, by the $10$th iteration in Figure LABEL:sub@fig:cs1:SSEBehaviour:4, both likelihood peaks are approximated well by the adSSLE approach.

The last iteration shows how the algorithm splits domains and adds new sample points. There is a clear clustering of subdomains and sample points near the likelihood peaks at $X=0.95$ and $X=1.05$.

The results from Eq. (22) show that without further computations it would be possible to directly extract the posterior moments by post-processing the SSLE coefficients. In the present bimodal case, however, the posterior moments are not very meaningful. Instead, the available posterior approximation gives a full picture of the inferred parameter $X|{\mathcal{Y}}$. It is shown together with the true posterior and the original prior distribution in Figure 5.

For this case study, non-adaptive experimental design approaches like the standard SSLE (Marelli et al., 2021) and the original SLE algorithm (Nagel and Sudret, 2016) will almost surely fail for the considered experimental design of $N_{\mathrm{ED}}=100$. In numerous trial runs these approaches did not manage to accurately reconstruct the likelihood function due to a lack of informative samples near the likelihood peaks.

This case study is a complex engineering problem that can be modified to exhibit high posterior concentration in the prior domain. It was originally presented in Nagel and Sudret (2016) and solved there using SLE. We again solve the same problem with SSLE and compare the performance of SLE (Nagel and Sudret, 2016), the original non-adaptive SSLE (Marelli et al., 2021) and the proposed adSSLE approach presented in Section 4.2. To investigate the performance of the algorithm in the case of high posterior concentrations (e.g. due to a large data-set), two instances of the problem with different discrepancy parameters are investigated.

Consider the diffusion-driven stationary heat transfer problem sketched in Figure LABEL:sub@fig:cs3:setup. It models a 2D plate with a background matrix of constant thermal conductivity $\kappa_{0}$ and $6$ inclusions with conductivities $\bm{\kappa}\stackrel{{\scriptstyle\text{def}}}{{=}}(\kappa_{1},\dots,\kappa_{6})^{\intercal}$. The diffusion driven steady state heat distribution is described by a heat equation in Euclidean coordinates $\bm{r}\stackrel{{\scriptstyle\text{def}}}{{=}}(r_{1},r_{2})^{\intercal}$ of the form

where the thermal conductivity field is denoted by $\kappa$ and the temperature field by $\tilde{T}$. The boundary conditions of the plate are given by a no-heat-flux Neumann boundary conditions on the left and right sides ($\partial\tilde{T}/r_{1}=0$), a Neumann boundary condition on the bottom (${\kappa_{0}\partial\tilde{T}/r_{2}=q_{2}}$) and a temperature $\tilde{T}_{1}$ Dirichlet boundary condition on the top.

We employ a finite element (FE) solver to solve the weak form of Eq. (31) by discretizing the domain into approximately $10^{5}$ triangular elements. A sample solution returned by the FE-solver is shown in Figure 6b.

In this example we intend to infer the thermal conductivities $\bm{\kappa}$ of the inclusions. We assume the same problem constants as in Nagel and Sudret (2016) (i.e., $q_{2}=2{,}000~{}\text{W/m}^{2}$, $\tilde{T}_{1}=200~{}\text{K}$, $\kappa_{0}=45\text{W/m/K}$). The forward model ${\mathcal{M}}$ takes as an input the conductivities of the inclusions $\bm{\kappa}$, solves the finite element problem and returns the steady state temperature $\tilde{\bm{T}}\stackrel{{\scriptstyle\text{def}}}{{=}}(\tilde{T}_{1},\dots,\tilde{T}_{20})^{\intercal}$ at the measurement points, i.e., ${\mathcal{M}}:\bm{\kappa}\mapsto\tilde{\bm{T}}$.

To solve the inverse problem, we assume a multivariate lognormal prior distribution with independent marginals on the inclusion conductivities, i.e. $\pi(\bm{\kappa})=\prod_{i=1}^{6}\mathcal{LN}(\kappa_{i}|\mu=30~{}\text{W/m/K},\sigma=6~{}\text{W/m/K})$. We further assume an additive Gaussian discrepancy model, which yields the likelihood function

with a discrepancy standard deviation of $\sigma$.

As measurements, we generate one temperature field with $\hat{\bm{\kappa}}\stackrel{{\scriptstyle\text{def}}}{{=}}(32,36,20,24,40,28)^{\intercal}~{}\text{W/m/K}$ and collect its values at $20$ points indicated by black dots in Figure LABEL:sub@fig:cs3:setup. We then perturb these temperature values with additive Gaussian noise and use them as the inversion data ${\mathcal{Y}}\stackrel{{\scriptstyle\text{def}}}{{=}}\bm{T}=(T_{1},\dots,T_{20})^{\intercal}$.

We look at two instances of this problem that differ only by the discrepancy parameter $\sigma$ from Eq. (32). The prior model response has a standard deviation of approximately $0.3~{}\text{K}$, depending on the measurement point $T_{i}$. We therefore solve the problem first with a large value $\sigma=0.25~{}\text{K}$ and second with a small value $\sigma=0.1~{}\text{K}$. As the discrepancy standard deviation determines how peaked the likelihood function is, the first problem has a likelihood function with a much wider support and is in turn significantly easier to solve than the second one. It is noted here that in practice, the peakedness of the likelihood function is either increased by a smaller discrepancy standard deviation, or the inclusion of additional experimental data.

To monitor the dependence of the algorithms on the number of likelihood evaluations, we solve both problems with a set of maximum likelihood evaluations $N_{\mathrm{ED}}=\{1{,}000;2{,}000,5{,}000;10{,}000;30{,}000\}$. The number of refinement samples is set to $N_{\mathrm{ref}}=1{,}000$.

As a benchmark, we use reference posterior samples generated by the affine-invariant ensemble sampler MCMC algorithm (Goodman and Weare, 2010) with $30{,}000$ steps and $50$ parallel chains, requiring a total of $N_{\mathrm{ED}}=1.5\cdot 10^{6}$ likelihood evaluations. Based on numerous heuristic convergence tests and due to the large number of MCMC steps, the resulting samples can be considered to accurately represent the true posterior distributions.

The results of the analyses are summarized in Figure 7, where the error measure $\eta$ is plotted against the number of likelihood evaluations for the large and small standard deviation case. For the large discrepancy standard deviation case, both SSLE approaches clearly outperform standard SLE w.r.t. the error measure $\eta$. This is most significant at mid-range experimental designs ($N_{\mathrm{ED}}={5{,}000;10{,}000}$), where SLE does not reach the required high degrees and fails to accurately approximate the likelihood function. At larger experimental designs SLE catches up to non-adaptive SSLE but is still outperformed by the proposed adSSLE approach. The real strength of the adaptive algorithm shows for the case of a small discrepancy standard deviation, where the limitations of fixed experimental designs become obvious. When the likelihood function is nonzero in a small subdomain of the prior, the global SLE and non-adaptive SSLE approach will fail in practice because of the insufficient number of samples placed in the informative regions. The adSSLE approach, however, works very well in these types of problems. It manages to identify the regions of interest and produces a likelihood approximation that accurately reproduces the posterior marginals.

Tables 1 and 2 show the convergence of the adSSLE method moment estimate (mean and variance) to the reference solution for a single run. In brackets next to the moment estimates $\xi$, the relative error $\epsilon\stackrel{{\scriptstyle\text{def}}}{{=}}|\xi_{\mathrm{MCMC}}-\xi_{\mathrm{SSLE}}|/\xi_{\mathrm{MCMC}}$ is also shown. Due to the non-strict positivity of the SSLE estimate, one variance estimate computed with Eq. (22) is negative and is therefore omitted from Table 2.