# A benchmark for the Bayesian inversion of coefficients in partial
differential equations

David Aristoff Department of Mathematics, 1874 Campus Delivery, Colorado State
University; Fort Collins, CO 80523; USA (aristoff@math.colostate.edu).
Wolfgang Bangerth Department of Mathematics, Department of Geosciences, 1874
Campus Delivery, Colorado State University; Fort Collins, CO 80523; USA
(bangerth@colostate.edu).

###### Abstract

Bayesian methods have been widely used in the last two decades to infer
statistical properties of spatially variable coefficients in partial
differential equations from measurements of the solutions of these equations.
Yet, in many cases the number of variables used to parameterize these
coefficients is large, and obtaining meaningful statistics of their values is
difficult using simple sampling methods such as the basic Metropolis-Hastings
(MH) algorithm – in particular if the inverse problem is ill-conditioned or
ill-posed. As a consequence, many advanced sampling methods have been
described in the literature that converge faster than MH, for example by
exploiting hierarchies of statistical models or hierarchies of discretizations
of the underlying differential equation.

At the same time, it remains difficult for the reader of the literature to
quantify the advantages of these algorithms because there is no commonly used
benchmark. This paper presents a benchmark Bayesian inverse problem – namely,
the determination of a spatially-variable coefficient, discretized by 64
values, in a Poisson equation, based on point measurements of the solution –
that fills the gap between widely used simple test cases (such as
superpositions of Gaussians) and real applications that are difficult to
replicate for developers of sampling algorithms. We provide a complete
description of the test case, and provide an open source implementation that
can serve as the basis for further experiments. We have also computed $2\times
10^{11}$ samples, at a cost of some 30 CPU years, of the posterior probability
distribution from which we have generated detailed and accurate statistics
against which other sampling algorithms can be tested.

###### AMS:

65N21, 35R30, 74G75

##  1 Introduction

Inverse problems are parameter estimation problems in which one wants to
determine unknown, spatially variable material parameters in a partial
differential equation (PDE) based on measurements of the solution. In the
deterministic approach, one in essence seeks that set of parameters for which
the solution of the PDE would best match the measured values; this approach is
widely used in many applications. On the other hand, the Bayesian approach to
inverse problems recognizes that all measurements are subject to measurement
errors and that models are also inexact; as a consequence, we ought to pose
the inverse problem as one that seeks a probability distribution describing
how likely it is that parameter values lie in a given interval or set. This
generalization of the perspective on inverse problems has long roots, but
first came to the attention of the wider scientific community through a 1987
book by Tarantola [49]. It was later followed by a significantly revised and
more accessible version by the same author [50] as well as numerous other
books on the subject; we mention [35] as one example, along with [3, 2, 17]
for tutorial-style introductions to the topic. The Bayesian approach to
inverse problems has been used in a wide variety of inverse applications: Too
many to mention in detail, but including acoustics [12], flow in the Earth’s
mantle [59], laminar and turbulent flow [13], ice sheet modeling [42],
astronomy [15], chemistry [26, 39], and groundwater modeling [33].

From a computational perspective, the primary challenge in Bayesian inverse
problems is that after discretizing the spatially variable parameters one
seeks to infer, one generally ends up with trying to characterize a finite-
but high-dimensional probability distribution $\pi(\theta)$ that describes the
relative likelihood of parameters $\theta$. In particular, we are typically
interested in computing the mean and standard deviation of this probability
distribution (i.e., which set of parameters $\left<\theta\right>$ on average
fits the measured data best, and what we know about its variability given the
uncertainty in the measured data). Computing these integral quantities in high
dimensional spaces can only be done through sampling methods such as Monte
Carlo Markov Chain (MCMC) algorithms. On the other hand, sampling in high-
dimensional spaces often suffers from a number of problems: (i) Long burn-in
times until a chain finally finds the region where the probability
distribution $\pi(\theta)$ has values substantially different from zero; (ii)
long autocorrelation length scales if $\pi(\theta)$ represents elongated,
curved, or “ridged” distributions; (iii) for some inverse problems, multi-
modality of $\pi(\theta)$ is also a problem that complicates the
interpretation of the posterior probability distribution; (iv) in many ill-
posed inverse problems, parameters have large variances that result in rather
slow convergence to reliable and accurate means.

In most high-dimensional applications, the result of these issues is then that
one needs very large numbers of samples to accurately characterize the desired
probability distribution. In the context of inverse problems, this problem is
compounded by the fact that the generation of every sample requires the
solution of the forward problem – i.e., generally, the expensive numerical
solution of a PDE. As a consequence, the solution of Bayesian inverse problems
is computationally exceptionally expensive.

The community has stepped up to this challenge over the past two decades:
Numerous algorithms have been developed to make the sampling process more
efficient. Starting from the simplest sampler, the Metropolis-Hastings
algorithms with a symmetric proposal distribution [32], ideas to alleviate
some of the problems include nonsymmetric proposal distributions [44], delayed
rejection [53], non-reversible samplers [21], piecewise deterministic Markov
processes [11, 54] including Hamiltonian Monte Carlo [38], adaptive methods
[31, 6, 43], randomize-then-optimize methods [9, 8, 7], affine invariant
samplers [29, 24], and combinations thereof [16, 30].

Other approaches introduce parallelism (e.g., the differential evolution
method and variations [51, 52, 56]), or hierarchies of models (see, for
example, the survey by Peherstorfer et al. [41] and references therein, as
well as [45, 23, 14, 57, 25]). Yet other methods exploit the fact that
discretizing the underlying PDE gives rise to a natural multilevel hierarchy
(see [20] among many others) or that one can use the structure of the
discretized PDE for efficient sampling algorithms [55, 36].

Many of these methods are likely vastly faster than the simplest sampling
methods that are often used. Yet, the availability of a whole zoo of possible
methods and their various possible combinations has also made it difficult to
assess which method really should be used if one wanted to solve a particular
inverse problem, and there is no consensus in the community on this topic.
Underlying this lack of consensus is that there is no widely used benchmark
for Bayesian inverse problems: Most of the papers above demonstrate the
qualities of their particular innovation using some small but artificial test
cases such as a superposition of Gaussians, and often a more elaborate
application that is insufficiently well-described and often too complex for
others to reproduce. As a consequence, the literature contains few examples of
comparisons of algorithms using test cases that reflect the properties of
actual inverse problems.

Our contribution here seeks to address this lack of widely used benchmarks. In
particular:

  * •

We provide a complete description of a benchmark that involves characterizing
a posterior probability distribution $\pi(\theta)$ on a 64-dimensional
parameter space that results from inverting data for a discretized coefficient
in a Poisson equation.

  * •

We explain in detail why this benchmark is at once simple enough to make
reproduction by others possible, yet difficult enough to reflect the real
challenges one faces when solving Bayesian inverse problems.

  * •

We provide highly accurate statistics for $\pi(\theta)$ that allow others to
assess the correctness of their own algorithms and implementations. We also
provide a performance profile for a simple Metropolis-Hastings sampler as a
baseline against which other methods can be compared.

To make adoption of this benchmark simpler, we also provide an open-source
implementation of the benchmark that can be adapted to experimentation on
other sampling methods with relative ease.

The remainder of this paper is structured as follows: In Section 2, we provide
a complete description of the benchmark. In Section 3, we then evaluate highly
accurate statistics of the probability distribution that solves the benchmark,
based on $2\text{\times}{10}^{11}$ samples we have computed. Section 4
provides a short discussion of what we hope the benchmark will achieve, along
with our conclusions. An appendix presents details of our implementation of
the benchmark (Appendix A), discusses a simple 1d benchmark for which one can
find solutions in a much cheaper way (Appendix B), and provides some
statistical background relevant to Section 3 (in Appendix C).

##  2 The benchmark for sampling algorithms for inverse problems

###  2.1 Design criteria

In the design of the benchmark described in this paper, we were guided by the
following principle:

> A good benchmark is neither too simple nor too complicated. It also needs to
> reflect properties of real-world applications.

Specifically, we design a benchmark for inferring posterior probability
distributions using sampling algorithms that correspond to the Bayesian
inversion of coefficients in partial differential equations. In other words,
cases where the relative posterior likelihood is computed by comparing
(functionals of) the forward solution of partial differential equations with
(simulations of) measured data.

The literature has many examples of papers that consider sampling algorithms
for such problems (see the references in the introduction). However, they are
typically tested only on cases that fall into the following two categories:

  * •

Simple posterior probability density functions (PDFs) that are given by
explicitly known expressions such as Gaussians or superpositions of Gaussians.
There are advantages to such test cases: (i) the probability distributions are
cheap to evaluate, and it is consequently possible to create essentially
unlimited numbers of samples; (ii) because the PDF is exactly known, exact
values for statistics such as the mean, covariances, or maximum likelihood
(MAP) points are often computable exactly, facilitating the quantitative
assessment of convergence of sampling schemes. On the other hand, these test
cases are often so simple that any reasonable sampling algorithm converges
relatively quickly, making true comparisons between different algorithms
difficult. More importantly, however, such simple test cases do not reflect
real-world properties of inverse problems: Most inverse problems are ill-
posed, nonlinear, and high-dimensional. They are often unimodal, but with PDFs
that are typically quite insensitive along certain directions in parameter
space, reflecting the ill-posedness of the underlying problem. Because real-
world problems are so different from simple artificial test cases, it is
difficult to draw conclusions from the performance of a new sampling algorithm
when applied to a simple test case.

  * •

Complex applications, such as the determination of the spatially variable oil
reservoir permeability from the production history of an oil field, or the
determination of seismic wave speeds from travel times of earthquake waves
from their source to receivers. Such applications are of course the target for
applying advanced sampling methods, but they make for poor benchmarks because
they are very difficult to replicate by other authors. As a consequence, they
are almost exclusively used in only the single paper in which a new sampling
algorithm is first described, and it is difficult for others to compare this
new sampling algorithm against previous ones, since they have not been tested
against the same, realistic benchmark.

We position the benchmark in this paper between these extremes. Specifically,
we have set out to achieve the following goals:

  * •

Reflect real properties of inverse problems: Our benchmark should reflect
properties one would expect from real applications such as the permeability or
seismic wave speed determinations mentioned above. We do not really know what
these properties are, but intuition and knowledge of the literature suggest
that they include very elongated and nonlinear probability distributions,
quite unlike Gaussians or their superpositions. In order for our benchmark to
reflect these properties, we base it on a partial differential equation.

  * •

High dimensional: Inverse problems are originally infinite dimensional, i.e.,
we seek parameters that are functions of space and/or time. In practice, these
need to be discretized, leading to finite- but often high-dimensional
problems. It is well understood that the resulting curse of dimensionality
leads to practical problems that often make the Bayesian inverse problem
extremely expensive to solve. At the same time, we want to reflect these
difficulties in our benchmark.

  * •

Computable with acceptable effort: A benchmark needs to have a solution that
is known to an accuracy sufficiently good to compare against. This implies
that it can’t be so expensive that we can only compute a few thousand or tens
of thousands of samples of the posterior probability distribution. This rules
out most real applications for which each forward solution, even on parallel
computers may take minutes, hours, or even longer. Rather, we need a problem
that can be solved in at most a second on a single processor to allow the
generation of a substantial number of samples.

  * •

Reproducible: To be usable by anyone, a benchmark needs to be completely
specified in all of its details. It also needs to be simple enough so that
others can implement it with reasonable effort.

  * •

Available: An important component of this paper is that we make the software
that implements the benchmark available as open source, see Appendix A. In
particular, the code is written in a modular way that allows evaluating the
posterior probability density for a given set of parameter values – i.e., the
key operation of all sampling methods. The code is also written in such a way
that it is easy to use in a multilevel sampling scheme where the forward
problem is solved with a hierarchy of successively more accurate
approximations.

###  2.2 Description of the benchmark

Given the design criteria discussed in the previous subsection, let us now
present the details of the benchmark. Specifically, we seek (statistics of) a
non-normalized posterior probability distribution $\pi(\theta|\hat{z})$ on a
parameter space $\theta\in\Theta={\mathbb{R}}^{64}$ of modestly high dimension
64 – large enough to be interesting, while small enough to remain feasible for
a benchmark. Here, we think of $\hat{z}$ as a set of measurements made on a
physical system that is used to infer information about the internal
parameters $\theta$ of the system. As is common in Bayesian inverse problems,
$\pi(\theta|\hat{z})$ is defined as the product of a likelihood times a prior
probability:

$\displaystyle\pi(\theta|\hat{z})\propto
L(\hat{z}|\theta)\,\pi_{\text{pr}}(\theta).$ (1)

Here, $L(z|\theta)$ describes how likely it would be to measure values $z$ if
$\theta$ were the “true” values of the internal parameters. $\pi_{\text{pr}}$
is a (not necessarily normalized) probability distribution encoding our prior
beliefs about the parameters. A complete description of the benchmark then
requires us to describe the values of $z$ and ways to evaluate the functions
$L$ and $\pi_{\text{pr}}$. We will split the definition of $L$ into a
discussion of the forward model and a statistical model of measurements in the
following.

####  2.2.1 The forward model

The setting we want to pursue is as follows: Let us imagine a membrane
stretched over a frame that bounds a domain $\Omega$ which, for simplicity we
assume to be the unit square $\Omega=(0,1)^{2}$. The membrane is subject to an
external vertical force $f(\mathbf{x})$ which for the purpose of this
benchmark we choose constant as $f(\mathbf{x})=10$. Furthermore, the membrane
has a spatially variable resistance $a(\mathbf{x})$ to deflection (for
example, it may have a variable thickness or may be made from different
materials). In this benchmark, we assume that $a(\mathbf{x})$ is piecewise
constant on a uniform $8\times 8$ grid as shown in Fig. 1, with the 64 values
that parameterize $a(\mathbf{x})$ given by the elements of the vector
$\theta_{0},\ldots,\theta_{63}$ as also indicated in the figure. In other
words, there is a 1:1 relationship between the vector $\theta$ and the
piecewise constant coefficient function
$a(\mathbf{x})=a^{\theta}(\mathbf{x})$.

.
0816243240485619172533414957210182634425058311192735435159412202836445260513212937455361614223038465462715233139475563
$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$

$\\!\small 0$

$\\!\small 13$

$\\!\small 26$

$\\!\small 39$

$\\!\small 52$

$\\!\small 65$

$\\!\small 78$

$\\!\small 91$

$\\!\small 104$

$\\!\small 117$

$\\!\small 130$

$\\!\small 143$

$\\!\small 156$

$\\!\small 1$

$\\!\small 14$

$\\!\small 27$

$\\!\small 40$

$\\!\small 53$

$\\!\small 66$

$\\!\small 79$

$\\!\small 92$

$\\!\small 105$

$\\!\small 118$

$\\!\small 131$

$\\!\small 144$

$\\!\small 157$

$\\!\small 2$

$\\!\small 15$

$\\!\small 28$

$\\!\small 41$

$\\!\small 54$

$\\!\small 67$

$\\!\small 80$

$\\!\small 93$

$\\!\small 106$

$\\!\small 119$

$\\!\small 132$

$\\!\small 145$

$\\!\small 158$

$\\!\small 3$

$\\!\small 16$

$\\!\small 29$

$\\!\small 42$

$\\!\small 55$

$\\!\small 68$

$\\!\small 81$

$\\!\small 94$

$\\!\small 107$

$\\!\small 120$

$\\!\small 133$

$\\!\small 146$

$\\!\small 159$

$\\!\small 4$

$\\!\small 17$

$\\!\small 30$

$\\!\small 43$

$\\!\small 56$

$\\!\small 69$

$\\!\small 82$

$\\!\small 95$

$\\!\small 108$

$\\!\small 121$

$\\!\small 134$

$\\!\small 147$

$\\!\small 160$

$\\!\small 5$

$\\!\small 18$

$\\!\small 31$

$\\!\small 44$

$\\!\small 57$

$\\!\small 70$

$\\!\small 83$

$\\!\small 96$

$\\!\small 109$

$\\!\small 122$

$\\!\small 135$

$\\!\small 148$

$\\!\small 161$

$\\!\small 6$

$\\!\small 19$

$\\!\small 32$

$\\!\small 45$

$\\!\small 58$

$\\!\small 71$

$\\!\small 84$

$\\!\small 97$

$\\!\small 110$

$\\!\small 123$

$\\!\small 136$

$\\!\small 149$

$\\!\small 162$

$\\!\small 7$

$\\!\small 20$

$\\!\small 33$

$\\!\small 46$

$\\!\small 59$

$\\!\small 72$

$\\!\small 85$

$\\!\small 98$

$\\!\small 111$

$\\!\small 124$

$\\!\small 137$

$\\!\small 150$

$\\!\small 163$

$\\!\small 8$

$\\!\small 21$

$\\!\small 34$

$\\!\small 47$

$\\!\small 60$

$\\!\small 73$

$\\!\small 86$

$\\!\small 99$

$\\!\small 112$

$\\!\small 125$

$\\!\small 138$

$\\!\small 151$

$\\!\small 164$

$\\!\small 9$

$\\!\small 22$

$\\!\small 35$

$\\!\small 48$

$\\!\small 61$

$\\!\small 74$

$\\!\small 87$

$\\!\small 100$

$\\!\small 113$

$\\!\small 126$

$\\!\small 139$

$\\!\small 152$

$\\!\small 165$

$\\!\small 10$

$\\!\small 23$

$\\!\small 36$

$\\!\small 49$

$\\!\small 62$

$\\!\small 75$

$\\!\small 88$

$\\!\small 101$

$\\!\small 114$

$\\!\small 127$

$\\!\small 140$

$\\!\small 153$

$\\!\small 166$

$\\!\small 11$

$\\!\small 24$

$\\!\small 37$

$\\!\small 50$

$\\!\small 63$

$\\!\small 76$

$\\!\small 89$

$\\!\small 102$

$\\!\small 115$

$\\!\small 128$

$\\!\small 141$

$\\!\small 154$

$\\!\small 167$

$\\!\small 12$

$\\!\small 25$

$\\!\small 38$

$\\!\small 51$

$\\!\small 64$

$\\!\small 77$

$\\!\small 90$

$\\!\small 103$

$\\!\small 116$

$\\!\small 129$

$\\!\small 142$

$\\!\small 155$

$\\!\small 168$

.

Fig. 1: Left: Numbering of the 64 cells on which the parameters are defined.
Right: Numbering and locations of the $13^{2}=169$ evaluation points at which
the solution is evaluated.

Then, an appropriate model to describe the vertical deflection $u(\mathbf{x})$
of the membrane would express $u$ as the solution of the following partial
differential equation that generalizes the Poisson equation:

$\displaystyle-\nabla\cdot\left[a(\mathbf{x})\nabla u(\mathbf{x})\right]$
$\displaystyle=f(\mathbf{x})$ $\displaystyle\text{in $\Omega$},$ (2)
$\displaystyle u(\mathbf{x})$ $\displaystyle=0$ $\displaystyle\text{on
$\partial\Omega$}.$ (3)

This model is of course not exactly solvable. But its solution can be
approximated using discretization. The way we define the likelihood $L$ then
requires us to specify exactly how we discretize this model. Concretely, we
define $u_{h}(\mathbf{x})$ as the solution of a finite element discretization
of (2)–(3) using a uniform $32\times 32$ mesh and a $Q_{1}$ (bilinear)
element. Because $f(\mathbf{x})$ is given, and because there is a 1:1
relationship between $\theta$ and $a(\mathbf{x})$, this discretized model then
implies that for each $\theta$ we can find a
$u_{h}(\mathbf{x})=u_{h}^{\theta}(\mathbf{x})$ that can be thought of as being
parameterized using the $1089$ degrees of freedom of the $Q_{1}$
discretization on the $32\times 32$ mesh. (However, of these $1089$ degrees of
freedom, 128 are on the boundary and are constrained to zero.) In other words,
using the $Q_{1}$ shape functions $\varphi_{k}(\mathbf{x})$, we can express
$u_{h}(\mathbf{x})=\sum_{k=0}^{1088}U_{k}\varphi_{k}(\mathbf{x})$. It is
important to stress that the mapping $\theta\mapsto
u_{h}^{\theta}(\mathbf{x})$ (or equivalently, $\theta\mapsto U^{\theta}$) is
nonlinear.

The function $u_{h}^{\theta}(\mathbf{x})$ can be thought of as the predicted
displacement at every point $\mathbf{x}\in\Omega$ if $\theta$ represented the
spatially variable stiffness coefficient of the membrane. In practice,
however, we can only measure finitely many things, and consequently define a
measurement operator ${\cal M}:u_{h}\mapsto z\in{\mathbb{R}}^{169}$ that
evaluates $u_{h}$ on a uniform $13\times 13$ grid of points
$\mathbf{x}_{k}\in\Omega$ so that
$\mathbf{x}_{k}=\left(\frac{i}{13+1},\frac{j}{13+1}\right),1\leq i,j\leq 13$
with $k=13(i-1)+(j-1)$. The locations of these points are also indicated in
Fig. 1. This last step then defines a linear mapping. Because of the
equivalence between the function $u_{h}$ and its nodal vector $U$, the
linearity of the measurement operator implies that we can write $z=MU$ with a
matrix $M\in{\mathbb{R}}^{169\times 1089}$ that is given by
$M_{kl}=\varphi_{l}(\mathbf{x}_{k})$.

In summary, a parameter vector $\theta\in{\mathbb{R}}^{64}$ then predicts
measurements $z^{\theta}\in{\mathbb{R}}^{169}$ using the following chain of
maps:

$\displaystyle\theta\mapsto a^{\theta}(\mathbf{x})\mapsto U^{\theta}\mapsto
z^{\theta}.$ (4)

The mapping $\theta\mapsto z^{\theta}$ is commonly called the “forward model”
as it predicts measurements $z^{\theta}$ if we knew the parameter values
$\theta$. The “inverse problem” is then of course the inverse operation: to
infer the parameters $\theta$ that describe a system based on measurements $z$
of its state $u$.

All of the steps of the forward model have been precisely defined above and
are easily computable with some basic knowledge of finite element methods (or
using the code discussed in Appendix A). The expensive step is to solve for
the nodal vector $U^{\theta}$, as this requires the assembly and solution of a
linear system of size $1089$.

######  Remark 1.

The $32\times 32$ mesh to define the forward model is chosen sufficiently fine
to resolve the exact solution $u$ reasonably well. At the same time, it is
coarse enough to allow for the rapid evaluation of the solution – even a
rather simple implementation should yield a solution in less than a second,
and a highly optimized implementation such as the one discussed in Appendix
A.1 will be able to do so in less than 5 milliseconds on modern hardware. As a
consequence, this choice of mesh allows for computing a large number of
samples, and consequently accurate quantitative comparisons of sampling
algorithms.

We also mention that the $32\times 32$ mesh for $u_{h}(\mathbf{x})$ is twice
more globally refined than the $8\times 8$ mesh used to define $a(\mathbf{x})$
in terms of $\theta$. It is clear to practitioners of finite element
discretizations of partial differential equations that the mesh for $u_{h}$
must be at the very least as fine as the one for the coefficient $a^{\theta}$
to obtain any kind of accuracy. On the other hand, these choices then leave
room for a hierarchy of models in which the forward model uses $8\times 8$,
$16\times 16$, and $32\times 32$ meshes; we expect that multilevel sampling
methods will use this hierarchy to good effect.

######  Remark 2.

In our experiments, we will choose the values of $\theta$ (and consequently of
$a(\mathbf{x})$) clustered around one. With the choice $f=10$ mentioned above,
this leads to a solution $u(\mathbf{x})$ with values in the range $0\ldots
0.95$. This then also implies that we should think of the numerical magnitude
of our measurements $z^{\theta}_{k}$ as ${\cal O}(1)$.

####  2.2.2 The likelihood $L(z|\theta)$

Given the predicted measurements $z^{\theta}$ that correspond to a given set
of parameters $\theta$, the likelihood $L(z|\theta)$ can be thought of as
expressing the (non-normalized) probability of actually obtaining $z$ in a
measurement if $\theta$ were the “correct” set of parameters. This is a
statement that encodes the measurement error of our measurement device.

For the purposes of this benchmark, we assume that these measurement errors
are identical and independently distributed for all 169 measurement points.
More specifically, we define the likelihood as the following (non-normalized)
probability function:

$\displaystyle
L(z|\theta)=\exp\left(-\frac{\|z-z^{\theta}\|^{2}}{2\sigma^{2}}\right)=\prod_{k=0}^{168}\exp\left(-\frac{(z_{k}-z^{\theta}_{k})^{2}}{2\sigma^{2}}\right),$
(5)

where we set $\sigma=0.05$ and where $z^{\theta}$ is related to $\theta$ using
the chain (4).

######  Remark 3.

We can think of (5) as encoding our belief that our measurement system
produces a Gaussian-distributed measurement $z_{k}\sim
N(z^{\theta}_{k},\sigma)$. Given that $z^{\theta}_{k}={\cal O}(1)$,
$\sigma=0.05$ implies a measurement error of 5%. This is clearly much larger
than the accuracy with which one would be able to determine the deflection of
a membrane in practice. On the other hand, we have chosen $\sigma$ this large
to ensure that the Bayesian inverse problem does not lead to a probability
distribution $\pi(\theta|\hat{z})$ that is so narrowly centered around a value
$\bar{\theta}$ that the mapping $\theta\mapsto z^{\theta}$ can be linearized
around $\bar{\theta}$ – in which case the likelihood $L(\hat{z}|\theta)$ would
become Gaussian, as also discussed in Appendix B. We will demonstrate in
Section 3.4 that indeed $\pi(\theta|\hat{z})$ is not Gaussian and, moreover,
is large along a curved ridge that can not easily be approximated by a
Gaussian either.

####  2.2.3 The prior probability $\pi_{\text{pr}}(\theta)$

Our next task is to describe our prior beliefs for the values of the
parameters. Given that the 64 values of $\theta$ describe the stiffness
coefficient of a membrane, it is clear that they must be positive.
Furthermore, as with many mechanical properties that can have values over vast
ranges,111For example, the Young’s modulus that is related to the stiffness of
a membrane, can range from 0.01 GPa for rubber to 200 GPa for typical steels.
Similarly, the permeability of typical oil reservoir rocks can range from 1 to
1000 millidarcies. reasonable priors are typically posed on the “order of
magnitude” (that is, the logarithm), not the size of the coefficient itself.
We express this through the following (non-normalized) probability
distribution:

$\displaystyle\pi_{\text{pr}}(\theta)=\prod_{i=0}^{63}\exp\left(-\frac{(\ln(\theta_{i})-\ln(1))^{2}}{2\sigma_{\text{pr}}^{2}}\right),$
(6)

where we choose $\sigma_{\text{pr}}=2$. We recognize the prior density of
$\ln(\theta_{k})$ as a Gaussian with mean $\sigma_{\text{pr}}^{2}$ and
standard deviation $\sigma_{\text{pr}}$.

Because this prior distribution is posed on the logarithm of the parameters,
the prior on the parameters themselves is very heavy-tailed, with mean values
$\left<\theta_{k}\right>_{\pi_{\text{pr}}}$ for each component much larger
than the value at which $\pi_{\text{pr}}$ takes on its maximum (which is at
$\theta_{k}=1$): Indeed, the mean of each $\theta_{k}$ with respect to
$\pi_{\textup{pr}}$ is about $403.43$.

We note that this prior probability is quite weak and, in particular, does not
assume any (spatial) correlation between parameters as is often done in
inverse problems [55, 48, 35]. The correlations we will observe in our
posterior probability (see Section 3.3) are therefore a consequence of the
likelihood function only.

####  2.2.4 The “true” measurements $\hat{z}$

The last piece necessary to describe the complete benchmark is the choice of
the “true” measurements $\hat{z}$ that we want to use to infer the statistical
properties of the parameters $\theta$. For the purposes of this benchmark, we
will use the 169 values for $\hat{z}$ given in Table 1.

Table 1: The “true” measurement values $\hat{z}_{k},k=0,\ldots,168$ used in the benchmark. The values are also available in the electronic supplemental material and are shown in full double precision accuracy to allow for exact reproduction of the benchmark. $\hat{z}_{0}$ | $0.060\,765\,117\,622\,593\,69$ | $\hat{z}_{60}$ | $0.623\,530\,057\,415\,691\,7$ | $\hat{z}_{120}$ | $0.514\,009\,175\,452\,694\,3$
---|---|---|---|---|---
| $0.096\,019\,101\,208\,484\,81$ |  | $0.555\,933\,270\,404\,593\,5$ |  | $0.555\,933\,270\,404\,596\,9$
| $0.123\,885\,251\,783\,858\,4$ |  | $0.467\,030\,499\,447\,417\,8$ |  | $0.567\,744\,369\,374\,330\,4$
| $0.149\,518\,411\,737\,520\,1$ |  | $0.349\,980\,914\,381\,1$ |  | $0.547\,825\,166\,529\,545\,3$
| $0.184\,159\,612\,754\,978\,4$ |  | $0.196\,882\,637\,462\,94$ |  | $0.489\,575\,966\,490\,898\,2$
| $0.217\,452\,502\,826\,112\,2$ |  | $0.217\,452\,502\,826\,125\,3$ |  | $0.410\,964\,174\,301\,917\,1$
| $0.225\,099\,616\,089\,869\,8$ |  | $0.412\,232\,953\,784\,340\,4$ |  | $0.395\,727\,260\,284\,338$
| $0.219\,795\,476\,900\,299\,3$ |  | $0.577\,945\,241\,983\,156\,6$ |  | $0.377\,894\,932\,200\,473\,4$
| $0.207\,469\,569\,837\,092\,6$ |  | $0.685\,968\,374\,925\,437\,2$ |  | $0.359\,626\,827\,185\,712\,4$
| $0.188\,999\,647\,766\,301\,6$ |  | $0.737\,310\,833\,139\,606\,3$ |  | $0.219\,125\,026\,894\,894\,8$
$\hat{z}_{10}$ | $0.163\,272\,253\,215\,372\,6$ | $\hat{z}_{70}$ | $0.745\,881\,198\,317\,824\,6$ | $\hat{z}_{130}$ | $0.163\,272\,253\,215\,368\,3$
| $0.127\,678\,248\,003\,818\,6$ |  | $0.727\,896\,802\,240\,655\,9$ |  | $0.285\,039\,780\,666\,332\,5$
| $0.077\,118\,459\,157\,893\,12$ |  | $0.690\,479\,353\,535\,775\,1$ |  | $0.373\,006\,008\,206\,081$
| $0.096\,019\,101\,208\,485\,52$ |  | $0.636\,917\,645\,271\,028\,8$ |  | $0.432\,532\,550\,635\,420\,7$
| $0.200\,058\,953\,336\,798\,3$ |  | $0.567\,744\,369\,374\,321\,5$ |  | $0.467\,030\,499\,447\,431\,5$
| $0.338\,559\,259\,195\,176\,6$ |  | $0.478\,473\,876\,486\,586\,7$ |  | $0.478\,473\,876\,486\,602\,3$
| $0.393\,430\,002\,464\,780\,6$ |  | $0.360\,219\,063\,282\,326\,2$ |  | $0.467\,712\,268\,759\,904\,1$
| $0.404\,022\,389\,246\,154\,1$ |  | $0.203\,179\,205\,473\,732\,5$ |  | $0.434\,171\,688\,106\,105\,5$
| $0.412\,232\,953\,784\,309\,2$ |  | $0.225\,099\,616\,089\,881\,8$ |  | $0.388\,186\,479\,011\,099$
| $0.410\,048\,009\,154\,555\,4$ |  | $0.410\,048\,009\,154\,578\,7$ |  | $0.377\,894\,932\,200\,460\,2$
$\hat{z}_{20}$ | $0.394\,915\,163\,718\,996\,8$ | $\hat{z}_{80}$ | $0.555\,561\,595\,613\,713\,7$ | $\hat{z}_{140}$ | $0.363\,336\,256\,718\,736\,4$
| $0.369\,787\,326\,479\,123\,2$ |  | $0.656\,123\,536\,696\,093\,8$ |  | $0.346\,445\,726\,190\,539\,9$
| $0.334\,018\,262\,359\,24$ |  | $0.711\,655\,887\,807\,071\,5$ |  | $0.209\,636\,232\,136\,565\,5$
| $0.285\,039\,780\,666\,338\,2$ |  | $0.727\,896\,802\,240\,657$ |  | $0.127\,678\,248\,003\,814\,8$
| $0.218\,426\,003\,247\,867\,1$ |  | $0.712\,192\,867\,867\,018\,7$ |  | $0.218\,426\,003\,247\,863\,4$
| $0.127\,112\,115\,635\,095\,7$ |  | $0.671\,218\,739\,142\,872\,9$ |  | $0.282\,169\,498\,339\,525\,2$
| $0.123\,885\,251\,783\,861\,1$ |  | $0.613\,915\,777\,559\,149\,2$ |  | $0.324\,831\,514\,891\,553\,5$
| $0.338\,559\,259\,195\,181\,9$ |  | $0.547\,825\,166\,529\,538\,1$ |  | $0.349\,980\,914\,381\,109\,7$
| $0.711\,928\,516\,276\,647\,5$ |  | $0.467\,712\,268\,759\,903\,1$ |  | $0.360\,219\,063\,282\,333\,3$
| $0.817\,571\,286\,175\,642\,8$ |  | $0.358\,765\,491\,100\,084\,8$ |  | $0.358\,765\,491\,100\,079\,9$
$\hat{z}_{30}$ | $0.683\,625\,411\,657\,810\,5$ | $\hat{z}_{90}$ | $0.205\,073\,429\,167\,591\,8$ | $\hat{z}_{150}$ | $0.353\,438\,997\,477\,926\,8$
| $0.577\,945\,241\,983\,115\,7$ |  | $0.219\,795\,476\,900\,309\,4$ |  | $0.364\,264\,009\,018\,228\,3$
| $0.555\,561\,595\,613\,689\,7$ |  | $0.394\,915\,163\,719\,015\,7$ |  | $0.359\,626\,827\,185\,69$
| $0.528\,518\,156\,173\,671\,9$ |  | $0.528\,518\,156\,173\,691\,1$ |  | $0.346\,445\,726\,190\,529\,5$
| $0.491\,439\,702\,849\,224$ |  | $0.621\,319\,720\,186\,747\,1$ |  | $0.326\,072\,895\,342\,464\,3$
| $0.440\,936\,749\,485\,328\,2$ |  | $0.674\,517\,904\,909\,440\,7$ |  | $0.180\,670\,595\,355\,394$
| $0.373\,006\,008\,206\,077\,2$ |  | $0.690\,479\,353\,535\,786$ |  | $0.077\,118\,459\,157\,892\,44$
| $0.282\,169\,498\,339\,521\,4$ |  | $0.671\,218\,739\,142\,878\,7$ |  | $0.127\,112\,115\,635\,096\,3$
| $0.161\,017\,673\,385\,773\,9$ |  | $0.617\,840\,828\,935\,951\,4$ |  | $0.161\,017\,673\,385\,775\,7$
| $0.149\,518\,411\,737\,525\,7$ |  | $0.545\,360\,502\,723\,788\,3$ |  | $0.183\,460\,041\,273\,014\,4$
$\hat{z}_{40}$ | $0.393\,430\,002\,464\,792\,9$ | $\hat{z}_{100}$ | $0.489\,575\,966\,490\,909$ | $\hat{z}_{160}$ | $0.196\,882\,637\,462\,944\,3$
| $0.817\,571\,286\,175\,656\,2$ |  | $0.434\,171\,688\,106\,127\,8$ |  | $0.203\,179\,205\,473\,735\,4$
| $0.943\,915\,462\,552\,765\,3$ |  | $0.353\,438\,997\,477\,945\,6$ |  | $0.205\,073\,429\,167\,588\,5$
| $0.801\,590\,411\,509\,512\,8$ |  | $0.208\,322\,749\,696\,134\,7$ |  | $0.208\,322\,749\,696\,124\,5$
| $0.685\,968\,374\,925\,402\,4$ |  | $0.207\,469\,569\,837\,099$ |  | $0.217\,959\,990\,927\,999\,8$
| $0.656\,123\,536\,696\,059\,9$ |  | $0.369\,787\,326\,479\,136\,6$ |  | $0.219\,125\,026\,894\,882\,2$
| $0.621\,319\,720\,186\,731\,5$ |  | $0.491\,439\,702\,849\,241\,2$ |  | $0.209\,636\,232\,136\,555\,1$
| $0.575\,361\,131\,500\,004\,9$ |  | $0.575\,361\,131\,500\,020\,3$ |  | $0.180\,670\,595\,355\,388\,7$
| $0.514\,009\,175\,452\,682\,3$ |  | $0.623\,530\,057\,415\,701\,7$ | $\hat{z}_{168}$ | $0.106\,796\,555\,001\,001\,3$
| $0.432\,532\,550\,635\,416\,5$ |  | $0.636\,917\,645\,271\,049\,7$ |  |
$\hat{z}_{50}$ | $0.324\,831\,514\,891\,548\,2$ | $\hat{z}_{110}$ | $0.613\,915\,777\,559\,157\,9$ |  |
| $0.183\,460\,041\,273\,008\,6$ |  | $0.545\,360\,502\,723\,793\,5$ |  |
| $0.184\,159\,612\,754\,991\,7$ |  | $0.433\,660\,492\,961\,285\,1$ |  |
| $0.404\,022\,389\,246\,183\,2$ |  | $0.410\,964\,174\,301\,931\,2$ |  |
| $0.683\,625\,411\,657\,843\,9$ |  | $0.388\,186\,479\,011\,124\,5$ |  |
| $0.801\,590\,411\,509\,539\,6$ |  | $0.364\,264\,009\,018\,259\,2$ |  |
| $0.787\,011\,956\,114\,497\,7$ |  | $0.217\,959\,990\,928\,014\,5$ |  |
| $0.737\,310\,833\,139\,580\,8$ |  | $0.188\,999\,647\,766\,301\,1$ |  |
| $0.711\,655\,887\,807\,046\,3$ |  | $0.334\,018\,262\,359\,246\,1$ |  |
| $0.674\,517\,904\,909\,428\,3$ |  | $0.440\,936\,749\,485\,338\,1$ |  |

In some sense, it does not matter where these values come from – we could have
measured them in an actual experiment, and used these values to infer the
coefficients of the system we measured on. On the other hand, for the purposes
of a benchmark, it might be interesting to know whether these “true
measurements” $\hat{z}$ correspond to a “true set of parameters”
$\hat{\theta}$ against which we can compare statistics such as the mean
$\left<\theta\right>$ of the posterior probability $\pi(\theta|\hat{z})$.

|  |
---|---|---
|  |
Fig. 2: Top left: The $8\times 8$ grid of values $\hat{\theta}$ used in
generating the “true measurements” $\hat{z}$ via the forward model discussed
in Section 2.2.4. The color scale matches that in Fig. 5. Top right: The
solution of the Poisson equation corresponding to $\hat{\theta}$ on the
$256\times 256$ mesh and using $Q_{3}$ finite elements. “True” measurements
$\hat{z}$ are obtained from this solution by evaluation at the points shown in
the right panel of Fig. 1. Bottom left: For comparison, the solution obtained
for the same set of parameters $\hat{\theta}$, but using the $32\times 32$
mesh and $Q_{1}$ element that defines the forward model. Bottom right: The
solution of this discrete forward model applied to the posterior mean
$\left<\theta\right>_{\pi(\theta|\hat{z})}$ that we will compute later; the
values of $\left<\theta\right>_{\pi(\theta|\hat{z})}$ are tabulated in Table
2, and visualized in Fig. 5.

Indeed, this is how we have generated $\hat{z}$: We chose a set of parameters
$\hat{\theta}$ that corresponds to a membrane of uniform stiffness
$a(\mathbf{x})=1$ except for two inclusions in which $a=0.1$ and $a=10$,
respectively. This set up is shown in Fig. 2.222This set up has the accidental
downside that both the set of parameters $\hat{\theta}$ and the set of
measurement points $\mathbf{x}_{k}$ at which we evaluate the solution are
symmetric about the diagonal of the domain. Since the same is true for our
finite element meshes, the exact solution of the benchmark results in a
probability distribution that is invariant to permutations of parameters about
the diagonal as well, and this is apparent in Fig. 5, for example. A better
designed benchmark would have avoided this situation, but we only realized the
issue after expending several years of CPU time. At the same time, the
expected symmetry of values allows for a basic check of the correctness of
inversion algorithms: If the inferred mean value
$\left<\theta_{7}\right>_{\pi(\theta|\hat{z})}$ is not approximately equal to
$\left<\theta_{63}\right>_{\pi(\theta|\hat{z})}$ – see the numbering shown in
the left panel of Fig. 1 – then something is wrong.  Using $\hat{\theta}$, we
then used the series of mappings as shown in (4) to compute $\hat{z}$.
However, to avoid an inverse crime, we have used a $256\times 256$ mesh and a
bicubic ($Q_{3}$) finite element to compute
$\hat{\theta}\mapsto\hat{u}_{h}\mapsto\hat{z}={\cal M}\hat{u}_{h}$, rather
than the $32\times 32$ mesh and a bilinear ($Q_{1}$) element used to define
the mapping $\theta\mapsto u_{h}\mapsto z^{\theta}={\cal M}u_{h}$.

As a consequence of this choice of higher accuracy (and higher computational
cost), we can in general not expect that there is a set of parameters $\theta$
for which the forward model of Section 2.2.1 would predict measurements
$z^{\theta}$ that are equal to $\hat{z}$. Furthermore, the presence of the
prior probability $\pi_{\text{pr}}$ in the definition of $\pi(\theta|\hat{z})$
implies that we should not expect that either the mean
$\left<\theta\right>_{\pi(\theta|\hat{z})}$ nor the MAP point
$\theta_{\text{MAP}}=\arg\max_{\theta}\pi(\theta|\hat{z})$ are equal or even
just close to the “true” parameters $\hat{\theta}$.

##  3 Statistical assessment of $\pi(\theta|\hat{z})$

The previous section provides a concise definition of the non-normalized
posterior probability density $\pi(\theta|\hat{z})$. Given that the mapping
$\theta\mapsto z^{\theta}$ is nonlinear and involves solving a partial
differential equation, there is no hope that $\pi(\theta|\hat{z})$ can be
expressed as an explicit formula. On the other hand, all statistical
properties of $\pi(\theta|\hat{z})$ can of course be obtained by sampling, for
example using algorithms such as the Metropolis-Hastings sampler [32].

In order to provide a useful benchmark, it is necessary that at least some
properties of $\pi(\theta|\hat{z})$ are known with sufficient accuracy to
allow others to compare the convergence of their sampling algorithms. To this
end, we have used a simple Metropolis-Hastings sampler to compute
$2\text{\times}{10}^{11}$ samples that characterize $\pi(\theta|\hat{z})$, in
the form of $N=$2000$$ Markov chains of length $N_{L}=10^{8}$ each. (Details
of the sampling algorithm used to obtain these samples are given in Appendix
A.2.) Using the program discussed in Appendix A.1, the effort to produce this
many samples amounts to approximately 30 CPU years on current hardware. On the
other hand, we will show below that this many samples are really necessary in
order to provide statistics of $\pi(\theta|\hat{z})$ accurately enough to
serve as reference values – at least, if one insists on using an algorithm as
simple as the Metropolis-Hastings method. In practice, we hope that this
benchmark is useful in the development of algorithms that are substantially
better than the Metropolis-Hastings method. In addition, when assessing the
convergence properties of a sampling algorithm, it is of course not necessary
to achieve the same level of accuracy as we obtain here.

In the following let us therefore provide a variety of statistics computed
from our samples, along with an assessment of the accuracy with which we
believe that we can state these results. In the following, we will denote by
$0\leq L<N=2000$ the number of the chain, and $0\leq\ell<N_{L}=10^{8}$ the
number of a sample $\theta_{L,\ell}$ on chain $L$. If we need to indicate one
of the 64 components of a sample, we will use a subscript index $k$ for this
purpose as already used in Section 2.2.1.

###  3.1 How informative is our data set?

While we have $N=2000$ chains, each with a large number $N_{L}=10^{8}$ of
samples per chain, a careful assessment needs to include an evaluation how
informative all of these samples really are. For example, if the samples on
each chain had a correlation length of $10^{7}$ because our Metropolis-
Hastings sampler converges only very slowly, then each chain really only
contains approximately ten statistically independent samples of
$\pi(\theta|\hat{z})$. Consequently, we could not expect great accuracy in
estimates of the mean value, covariance matrices, and other quantities
obtained from each of the chains. Similarly, if the “burn-in” time of the
sampler is a substantial fraction of the chain lengths $N_{L}$, then we would
have to throw away many of the early samples.

To assess these questions, we have computed the autocovariance matrices

$\displaystyle AC_{L}(s)$
$\displaystyle=\left<[\theta_{L,\ell}-\left<\theta\right>_{L}][\theta_{L,\ell-s}-\left<\theta\right>_{L}]^{T}\right>_{L}$
$\displaystyle=\frac{1}{N_{L}-s-1}\sum_{\ell=s}^{N_{L}-1}[\theta_{L,\ell}-\left<\theta\right>_{L}][\theta_{L,\ell-s}-\left<\theta\right>_{L}]^{T}$
(7)

between samples $s$ apart on chain $L$. We expect samples with a small lag $s$
to be highly correlated (i.e., $AC_{L}(s)$ to be a matrix that is large in
some sense), whereas for large lags $s$, samples should be uncorrelated and
$AC_{L}(s)$ should consequently be small. A rule of thumb is that samples at
lags $s$ can be considered decorrelated from each other if $AC_{L}(s)\leq
10^{-2}AC_{L}(0)$ entrywise; see Appendix C.

Fig. 3: Decay of the trace of the autocovariance matrices $AC_{L}(s)$ with the
sample lag $s$. The light blue curves show the traces of $AC_{L}(s)$ for
twenty of our chains. The red curve shows the trace of the averaged
autocovariance matrices, $AC(s)=\frac{1}{N}\sum_{L=0}^{N-1}AC_{L}(s)$. The
thick yellow line corresponds to the function $4\cdot 10^{6}\cdot e^{-s/5300}$
and shows the expected, asymptotically exponential decay of the
autocovariance. The dashed green line indicates a reduction of the
autocovariance by roughly a factor of 100 compared to its starting values.

Fig. 3 shows the trace of these autocovariance matrices for several of our
chains. (We only computed the autocovariance at lags $s=0,100,200,300,\ldots$
up to $s=20,000$ because of the cost of computing $AC_{L}(s)$.) The curves
show that the autocovariances computed from different chains all largely
agree, and at least asymptotically decay roughly exponentially with $s$ as
expected. The data also suggest that the autocorrelation length of our chains
is around $N_{AC}=10^{4}$ – in other words, each of our chains should result
in approximately $N_{L}/N_{AC}=10^{4}$ meaningful and statistically
independent samples.

To verify this claim, we estimated the integrated autocovariance [47] using

$I\\!AC\approx\frac{1}{N}\sum_{L=0}^{N-1}100\sum_{s=-200}^{200}AC_{L}(|100s|).$
(8)

The integrated autocovariance is obtained by summing up the autocovariance.
(The factor of $100$ appears because we only computed $AC_{L}(s)$ at lags that
are multiples of $100$.) We show in Appendix C that the integrated
autocovariance leads to the following estimate of the effective sample size:

$\textup{Effective sample size per
chain}\approx\frac{N_{L}}{\lambda_{\max}(C^{-1}\cdot I\\!AC)}\approx 1.3\times
10^{4},$ (9)

where $\lambda_{\max}$ indicates the maximum eigenvalue and
$C=\frac{1}{N}\sum_{L=0}^{N-1}AC_{L}(0)$ is the covariance matrix; see also
(10) below. This is in good agreement with Fig. 3 and the effective sample
size derived from it.

Fig. 4: Two perspectives on the length of the burn-in phase of our sampling
scheme. Left: Values of components $k=0,9,36,54$, and $63$ of samples
$\theta_{\ell}$ for $\ell=0,\ldots,100,000$, for one, randomly chosen chain.
(For the geometric locations of the parameters that correspond to these
components, see Fig. 1.) Right: Across-chain averages
$\frac{1}{N}\sum_{L=0}^{N-1}\theta_{L,\ell}$ for the same components $k$ as
above. Both images also shows the mean values
$\left<\theta\right>_{\pi(\theta|\hat{z})}$ for these five components as
dashed lines.

This leaves the question of how long the burn-in period of our sampling scheme
is. Fig. 4 shows two perspectives on this question. The left panel of the
figure shows several components $\theta_{\ell,k}$ of the samples of one of our
chains for the first few autocorrelation lengths. The data shows that there is
at least no obvious “burn-in” period on this scale that would require us to
throw away a substantial part of the chain. At the same time, it also
illustrates that the large components of $\theta$ are poorly constrained and
vary on rather long time scales that make it difficult to assess convergence
to the mean. The right panel shows across-chain averages of the $\ell$th
samples, more clearly illustrating that the burn-in period may only last
around 20,000 samples – that is, that only around 2 of the approximately
10,000 statistically independent samples of each chain are unreliable.

Having thus convinced ourselves that it is safe to use all $10^{8}$ samples
from all chains, and that there is indeed meaningful information contained in
them, we will next turn our attention towards computing statistical
information that characterizes $\pi(\theta|\hat{z})$ and against which other
implementations of sampling methods can compare their results.

###  3.2 The mean value of $\pi(\theta|\hat{z})$

The simplest statistic one can compute from samples of a distribution is the
mean value. Table 2 shows the 64 values that characterize the mean

$\displaystyle\left<\theta_{k}\right>_{\pi(\theta|\hat{z})}=\frac{1}{N}\sum_{L=0}^{N-1}\left(\frac{1}{N_{L}}\sum_{\ell=0}^{N_{L}-1}\theta_{L,\ell,k}\right).$

A graphical representation of $\left<\theta\right>_{\pi(\theta|\hat{z})}$ is
shown Fig. 5 and can be compared to the “true” values $\hat{\theta}$ shown in
Fig. 2.333By comparing Fig. 2 and the data of Table 2 and Fig. 5, it is clear
that for some parameters, the mean
$\left<\theta_{k}\right>_{\pi(\theta|\hat{z})}$ is far away from the value
$\hat{\theta}_{k}$ used to generate the original data $\hat{z}$ – principally
for those parameters that correspond to large values $\hat{\theta}_{k}$, but
also the “white cross” below and left of center. For the first of these two
places, we can first note that the prior probability $\pi_{\text{pr}}$ defined
in (6) is quite heavy-tailed, with a mean far larger than where its maximum is
located. And second, by realizing that a membrane that is locally very stiff
is not going to deform in a substantially different way in response to a force
from one that is even stiffer in that region – in other words, in areas where
the coefficient $a^{\theta}(\mathbf{x})$ is large, the likelihood function (5)
is quite insensitive to the exact values of $\theta$, and the posterior
probability will be dominated by the prior $\pi_{\text{pr}}$ with its large
mean. For the “white cross”, one can make plausible that the likelihood is
uninformative and that, consequently, mean value and variances are again
determined by the prior. To understand why this is so, one could imagine by
analogy what would happen if one could measure the solution $u(\mathbf{x})$ of
(2)–(3) exactly and everywhere, instead of only at a discrete set of points.
In that case, we would have $u(\mathbf{x})=z(\mathbf{x})$, and we could infer
the coefficient $a(\mathbf{x})$ by solving (2)–(3) for the coefficient instead
of for $u$. This leads to the advection-reaction equation $-\nabla
z(\mathbf{x})\cdot\nabla a(\mathbf{x})-(\Delta
z(\mathbf{x}))a(\mathbf{x})=f(\mathbf{x})$, which is ill-posed and does not
provide for a stable solution $a(\mathbf{x})$ at those places where $\nabla
z(\mathbf{x})=\nabla u(\mathbf{x})\approx 0$. By comparison with Fig. 2, we
can see that at the location of the white cross, we could not identify the
coefficient at one point even if we had measurements available everywhere, and
not stably so in the vicinity of that point. We can expect that this is also
so in the discrete setting of this benchmark – and that consequently, at this
location, only the prior provides information.

Table 2: Sample means $\left<\theta\right>_{\pi(\theta|\hat{z})}$ for the 64
parameters, along with their estimated 2-sigma uncertainties.

$\left<\theta\right>_{0}=$ |  $76.32$ $\pm$ | $0.30$
---|---|---
|  $1.2104$ $\pm$ | $0.0094$
|  $0.977\,380$ $\pm$ | $0.000\,051$
|  $0.882\,007$ $\pm$ | $0.000\,039$
|  $0.971\,859$ $\pm$ | $0.000\,048$
|  $0.947\,832$ $\pm$ | $0.000\,064$
|  $1.085\,29$ $\pm$ | $0.000\,11$
|  $11.39$ $\pm$ | $0.10$
$\left<\theta\right>_{8}=$ |  $1.119$ $\pm$ | $0.011$
|  $0.093\,721\,5$ $\pm$ | $0.000\,002\,7$
|  $0.115\,799\,2$ $\pm$ | $0.000\,003\,9$
|  $0.5815$ $\pm$ | $0.0022$
|  $0.9472$ $\pm$ | $0.0079$
|  $6.258$ $\pm$ | $0.079$
|  $9.334$ $\pm$ | $0.090$
|  $1.081\,51$ $\pm$ | $0.000\,11$
$\left<\theta\right>_{16}=$ |  $0.977\,449$ $\pm$ | $0.000\,052$
|  $0.115\,796\,2$ $\pm$ | $0.000\,003\,8$
|  $0.461$ $\pm$ | $0.020$
|  $267.01$ $\pm$ | $0.55$
|  $30.87$ $\pm$ | $0.19$
|  $7.189$ $\pm$ | $0.089$
|  $12.39$ $\pm$ | $0.11$
|  $0.949\,863$ $\pm$ | $0.000\,073$

$\left<\theta\right>_{24}=$ |  $0.881\,977$ $\pm$ | $0.000\,039$
---|---|---
|  $0.5828$ $\pm$ | $0.0020$
|  $267.72$ $\pm$ | $0.62$
|  $369.35$ $\pm$ | $0.64$
|  $234.59$ $\pm$ | $0.53$
|  $13.29$ $\pm$ | $0.14$
|  $22.36$ $\pm$ | $0.16$
|  $0.988\,806$ $\pm$ | $0.000\,074$
$\left<\theta\right>_{32}=$ |  $0.971\,900$ $\pm$ | $0.000\,049$
|  $0.9509$ $\pm$ | $0.0079$
|  $30.76$ $\pm$ | $0.19$
|  $233.93$ $\pm$ | $0.52$
|  $1.169$ $\pm$ | $0.012$
|  $0.8327$ $\pm$ | $0.0057$
|  $88.52$ $\pm$ | $0.33$
|  $0.987\,809$ $\pm$ | $0.000\,079$
$\left<\theta\right>_{40}=$ |  $0.947\,816$ $\pm$ | $0.000\,065$
|  $6.260$ $\pm$ | $0.076$
|  $7.119$ $\pm$ | $0.087$
|  $13.20$ $\pm$ | $0.13$
|  $0.8327$ $\pm$ | $0.0035$
|  $176.73$ $\pm$ | $0.44$
|  $283.38$ $\pm$ | $0.58$
|  $0.914\,212$ $\pm$ | $0.000\,077$

$\left<\theta\right>_{48}=$ |  $1.085\,21$ $\pm$ | $0.000\,11$
---|---|---
|  $9.386$ $\pm$ | $0.089$
|  $12.44$ $\pm$ | $0.12$
|  $22.50$ $\pm$ | $0.17$
|  $88.57$ $\pm$ | $0.33$
|  $283.41$ $\pm$ | $0.57$
|  $218.65$ $\pm$ | $0.49$
|  $0.933\,451$ $\pm$ | $0.000\,087$
$\left<\theta\right>_{56}=$ |  $11.35$ $\pm$ | $0.11$
|  $1.081\,43$ $\pm$ | $0.000\,11$
|  $0.949\,869$ $\pm$ | $0.000\,074$
|  $0.988\,770$ $\pm$ | $0.000\,074$
|  $0.987\,866$ $\pm$ | $0.000\,083$
|  $0.914\,247$ $\pm$ | $0.000\,077$
|  $0.933\,426$ $\pm$ | $0.000\,087$
|  $1.599\,84$ $\pm$ | $0.000\,30$

To assess how accurately we know this average, we consider that we have
$N=2000$ chains of length $N_{L}=10^{8}$ each, and that each of these has its
own chain averages

$\displaystyle\left<\theta_{k}\right>_{L}=\frac{1}{N_{L}}\sum_{\ell=0}^{N_{L}-1}\theta_{L,\ell,k}.$

The ensemble average $\left<\theta_{k}\right>_{\pi(\theta|\hat{z})}$ is of
course the average of the chain averages $\left<\theta_{k}\right>_{L}$ across
chains, but the chain averages vary between themselves and we can compute the
standard deviation of these chain averages as

$\displaystyle\text{stddev}\left(\left<\theta_{k}\right>_{L}\right)=\left[\frac{1}{N}\sum_{L=0}^{N-1}\left(\left<\theta_{k}\right>_{L}-\left<\theta_{k}\right>_{\pi(\theta|\hat{z})}\right)^{2}\right]^{1/2}.$

Under standard assumptions, and assuming that the posterior is Gaussian, we
can then estimate that we know the ensemble averages
$\left<\theta_{k}\right>_{\pi(\theta|\hat{z})}$ to within an accuracy of
$\pm\frac{1}{\sqrt{N}}\text{stddev}\left(\left<\theta_{k}\right>_{L}\right)$
with 68% (1-sigma) certainty, and with an accuracy of
$\pm\frac{2}{\sqrt{N}}\text{stddev}\left(\left<\theta_{k}\right>_{L}\right)$
with 95% (2-sigma) certainty. In reality, the posterior is not Gaussian (see
Section 3.4), and the argument is not true as stated; however, computing
2-sigma values for all parameters is still a useful metric for how accurately
we know each of the parameters.

This 2-sigma accuracy is also provided in Table 2. For all but parameter
$\theta_{18}$ (for which the relative 2-sigma uncertainty is $4.5\%$), the
relative uncertainty in $\left<\theta\right>$ is between $0.003\%$ and
$1.3\%$. In other words, the table provides nearly two certain digits for all
but one parameter, and four digits for at least half of all parameters.

|
---|---
Inferred mean values $\left<\theta_{k}\right>_{\pi(\theta|\hat{z})}$. |  Variances $C_{kk}$.
Fig. 5: Left: Visualization on an $8\times 8$ grid of mean values
$\left<\theta\right>_{\pi(\theta|\hat{z})}$ obtained from our samples. This
figure should be compared with the values $\hat{\theta}$ used as input to
generate the “true measurements” $\hat{z}$ and shown in Fig. 2. Right:
Variances $C_{kk}$ of the parameters. Dark colors indicate that a parameter is
accurately known; light colors that the variance is large. Standard deviations
(the square roots of the variances) are larger than the mean values in some
cases because of the heavy tails in the distributions of parameters.

###  3.3 The covariance matrix of $\pi(\theta|\hat{z})$ and its properties

The second statistic we demonstrate is the covariance matrix,

$\displaystyle\begin{split}C_{L}&=\frac{1}{N_{L}-1}\sum_{\ell=0}^{N_{L}-1}\left(\theta_{L,\ell}-\left<\theta\right>_{\pi(\theta|\hat{z})}\right)\left(\theta_{L,\ell}-\left<\theta\right>_{\pi(\theta|\hat{z})}\right)^{T},\\\
C&=\frac{1}{N}\sum_{L=0}^{N-1}C_{L}.\end{split}$ (10)

While conceptually easy to compute, in practice it is substantially harder to
obtain accuracy in $C$ than it is to compute accurate means
$\left<\theta\right>_{\pi(\theta|\hat{z})}$: While we know the latter to two
or more digits of accuracy, see Table 2, there is substantial variation
between the matrices $C_{L}$.444For diagonal entries $C_{L,kk}$, the standard
deviation of the variation between chains is between 0.0024 and 37.7 times the
corresponding entry $C_{kk}$ of the average covariance matrix. The variation
can be even larger for the many small off-diagonal entries. On the other hand,
the average (across chains) difference $\|C_{L}-C\|_{F}$ is $0.9\|C\|_{F}$.
This would suggest that we don’t know very much about these matrices, but as
shown in the rest of the section, qualitative measures can be extracted
robustly. The remainder of this section therefore only provides qualitative
conclusions we can draw from our estimate of the covariance matrix, rather
than providing quantitive numbers.

First, the diagonal entries of $C$, $C_{kk}$, provide the variances of the
statistical distribution of $\theta_{k}$, and are shown on the right of Fig.
5; the off-diagonal entries $C_{k\ell}$ suggest how correlated parameters
$\theta_{k}$ and $\theta_{l}$ are and are depicted in Fig. 6.

Fig. 6: Left: A visualization of the covariance matrix $C$ computed from the
posterior probability distribution $\pi(\theta|\hat{z})$. The sub-structure of
the matrix in the form of $8\times 8$ tiles represents that geographically
neighboring – and consequently correlated – parameters are either a distance
of $\pm 1$ or $\pm 8$ apart. Right: Correlation matrix
$D_{ij}=\frac{C_{ij}}{\sqrt{C_{ii}}\sqrt{C_{jj}}}$.

In the context of inverse problems related to partial differential equations,
it is well understood that we expect the parameters to be highly correlated.
This can be understood intuitively given that we are thinking of a membrane
model: If we increased the stiffness value on one of the $8\times 8$ pixels
somewhat, but decreased the stiffness value on a neighboring correspondingly,
then we would expect to obtain more or less the same global deformation
pattern – maybe there are small changes at measurement points close to the
perturbation, but for measurement points far away the local perturbation will
make little difference. As a consequence, we should expect that
$L(\hat{z}|\theta)\approx L(\hat{z}|\tilde{\theta})$ where $\theta$ and
$\tilde{\theta}$ differ only in two nearby components, one component of
$\tilde{\theta}$ being slightly larger and the other being slightly smaller
than the corresponding component of $\theta$. If the changes are small, then
we will also have that $\pi(\theta|\hat{z})\approx\pi(\tilde{\theta}|\hat{z})$
– in other words, we would expect that $\pi$ is approximately constant in the
secondary diagonal directions
$(0,\ldots,0,+\varepsilon,0,\ldots,0,-\varepsilon,0,\ldots 0)$ in $\theta$
space.

On the other hand, increasing (or decreasing) the stiffness value in both of
two adjacent pixels just makes the membrane overall more (or less) stiff, and
will yield different displacements at all measurement locations. Consequently,
we expect that the posterior probability distribution $\pi(\theta|\hat{z})$
will strongly vary in the principal diagonal directions
$(0,\ldots,0,+\varepsilon,0,\ldots,0,+\varepsilon,0,\ldots 0)$ in $\theta$
space.

We can illustrate this by computing two-dimensional histograms of the samples
for parameters $\theta_{k}$ and $\theta_{l}$ corresponding to neighboring
pixels – equivalent to a two-dimensional marginal distribution. We show such
histograms in Fig. 7. These also indicate that the posterior probability
distribution $\pi(\theta|\hat{z})$ is definitely not Gaussian – see also
Remark 3.

. .

Fig. 7: Pairwise marginal probabilities for $\theta_{45}$ and $\theta_{46}$
(left) and for $\theta_{53}$ and $\theta_{54}$ (right), respectively. See Fig.
1 for the relative locations of these parameters. These marginal distributions
illustrate the anti-correlation of parameters: If one is large, the other is
most likely small, and vice versa.

A better way to illustrate correlation is to compute a singular value
decomposition of the covariance matrix $C$. Many inverse problems have only a
relatively small number of large singular values of $C$ [55, 22, 12, 59, 42,
13], suggesting that only a finite number of modes is resolvable with the data
available – in other words, the problem is ill-posed. Fig. 8 shows the
singular values of the covariance matrix $C$ for the current case. The data
suggests that from the 169 measured pieces of (noisy) data, a deterministic
inverse problem could only recover some 25-30 modes of the parameter vector
$\theta\in{\mathbb{R}}^{64}$ with reasonable accuracy.555The figure shows the
spread of each of the eigenvalues of the within-chain matrices $C_{L}$ in
blue, and the eigenvalues of the across-chain matrix $C$ in red. One would
expect the latter to be well approximated by the former, and that is true for
the largest and smallest eigenvalues, but not for the ones in the middle.
There are two reasons for this: First, each of the $C_{L}$ is nearly singular,
but because each chain is finite, the poorly explored directions are different
from one chain to the next. At the same time, it is clear that the sum of
(different) singular matrices may actually be “less singular”, with fewer
small eigenvalues, and this is reflected in the graph. A second reason is that
we computed the eigenvalues of each of the $C_{L}$ and ordered them by size
when creating the plot, but without taking into account the associated
eigenspaces. As a consequence, if one considers the, say, 32nd largest
eigenvalue of $C$, the figure compares it with the 32nd largest eigenvalues of
all of the $C_{L}$, when the correct comparison would have been with those
eigenvalues of the matrices $C_{L}$ whose eigenspace is most closely aligned;
this may be an eigenvalue elsewhere in the order, and the effect will likely
be the most pronounced for those eigenvalues whose sizes are the least well
constrained. The conclusions to be drawn from Fig. 8 are therefore not the
actual sizes of eigenvalues, but the number of “large” eigenvalues. This
observation is robust, despite the inaccuracies in our determination of $C$.

Fig. 8: Singular values of the covariance matrix $C$ as defined in (10). Red
bars indicate eigenvalues of the across-chain averaged covariance matrix. Blue
boxes correspond to the 25th and 75th percentiles of the corresponding
eigenvalues of the covariance matrices $C_{L}$ of individual chains; vertical
bars extend to the minimum and maximum across chains for the $k$th eigenvalue
of the matrices $C_{L}$; blue bars in the middle of boxes indicate the median
of these eigenvalues.

###  3.4 Higher moments of $\pi(\theta|\hat{z})$

In some sense, solving Bayesian inverse problems is not very interesting if
the posterior distribution $p(\theta|\hat{z})$ for the parameters is Gaussian,
or at least approximately so, because it can be done much more efficiently by
computing the maximum likelihood estimator through a deterministic inverse
problem, and then computing the covariance matrix via the Hessian of the
deterministic (constrained) optimization problem. For example, [55] provides
an excellent overview of the techniques that can be used in this case. Because
of these simplifications, it is of interest to know how close the posterior
density of this benchmark is to a multi-dimensional Gaussian.

Fig. 9: Histograms (marginal distributions) for all 64 components of $\theta$,
accumulated over all $2\text{\times}{10}^{11}$ samples. The histograms for
those parameters whose values were 0.1 for the purposes of generating the
vector $\hat{z}$ are shown in red, those whose values were 10 in orange, and
all others in blue. The figure highlights histograms for some of those
components $\theta_{k}$ whose marginal distributions are clearly neither
Gaussian nor log-Gaussian. Note that the histograms were generated with bins
whose size increases exponentially from left to right, so that they are
depicted with equal size in the figure given the logarithmic $\theta_{k}$
axis. Histograms with equal bin size on a linear scale would be more heavily
biased towards the left, but with very long tails to the right. See Fig. 7 for
(pair) histograms using a linear scale.

To evaluate this question, Fig. 9 shows histograms of all of the parameters,
using 1000 bins that are equally spaced in logarithmic space; i.e., for each
component $k$, we create 1000 bins between -3 and +3 and sort samples into
these bins based on $\log_{10}(\theta_{k})$. It is clear that many of the
parameters have heavy tails and can, consequently, not be approximated well by
Gaussians. On the other hand, given the prior distribution (6) we attached to
each of the parameters, it would make sense to conjecture that the logarithms
$\log(\theta_{k})$ might be Gaussian distributed.

If that were so, the double-logarithmic plot shown in the figure would consist
of histograms in the form of parabolas open to the bottom, and again a simpler
– and presumably cheaper to compute – representation of $\pi(\theta|\hat{z})$
would be possible. However, as the figure shows, this too is clearly not the
case: While some parameters seem to be well described by such a parabola, many
others have decidedly non-symmetric histograms, or shapes that are simply not
parabolic. As a consequence, we conclude that the benchmark is at least not
boring in the sense that its posterior distribution could be computed in some
comparably much cheaper way.

###  3.5 Rate of convergence to the mean

The data provided in the previous subsections allows checking whether a
separate implementation of this benchmark converges to the same probability
distribution $\pi(\theta|\hat{z})$. However, it does not help in assessing
whether it does so faster or slower than the simplistic Metropolis-Hastings
method used herein. Indeed, as we will outline in the next section, we hope
that this work spurs the development and evaluation of methods that can
achieve the same results without needing more than $10^{11}$ samples.

To this end, let us here provide metrics for how fast our method converges to
the mean $\left<\theta\right>_{\pi(\theta|\hat{z})}$ discussed in Section 3.2.
More specifically, if we denote by
$\left<\theta\right>_{L,n}=\frac{1}{n}\sum_{\ell=0}^{n-1}\theta_{L,\ell}$ the
running mean of samples zero to $n-1$ on chain $L$, then we are interested in
how fast it converges to the mean. We measure this using the following error
norm

$\displaystyle e_{L}(n)$
$\displaystyle=\left\|\text{diag}\,(\left<\theta\right>_{\pi(\theta|\hat{z})})^{-1}\left(\left<\theta\right>_{L,n}-\left<\theta\right>_{\pi(\theta|\hat{z})}\right)\right\|$
$\displaystyle=\left[\left(\left<\theta\right>_{L,n}-\left<\theta\right>_{\pi(\theta|\hat{z})}\right)^{T}\text{diag}\,(\left<\theta\right>_{\pi(\theta|\hat{z})})^{-2}\left(\left<\theta\right>_{L,n}-\left<\theta\right>_{\pi(\theta|\hat{z})}\right)\right]^{1/2}.$
(11)

The weighting by a diagonal matrix containing the inverses of the estimated
parameters $\left<\theta\right>_{\pi(\theta|\hat{z})}$ (given in Table 2 and
known to sufficient accuracy) ensures that the large parameters with their
large variances do not dominate the value of $e_{L}(n)$. In other words,
$e_{L}(n)$ corresponds to the “root mean squared relative error”.666A possibly
better choice for the weighting would be to use the inverses of the diagonal
entries of the covariance matrix – i.e., the variances of the recovered
marginal probability distributions of each parameter. However, these are only
approximately known – see the discussion in Section 3.3 – and consequently do
not lend themselves for a concise definition of a reproducible benchmark.

Fig. 10 shows the convergence of a few chains to the ensemble average. While
there is substantial variability between chains, it is clear that for each
chain $e_{L}(n)^{2}\rightarrow 0$ and, furthermore, that this convergence
follows the classic one-over-$n$ convergence of statistical sampling
algorithms. Indeed, averaging $e_{L}(n)^{2}$ over all chains,

$\displaystyle e(n)^{2}=\frac{1}{N}\sum_{L=0}^{N-1}e_{L}(n)^{2},$

the behavior of this decay of the “average” square error $e(n)^{2}$ can be
approximated by the following formula that corresponds to the orange line in
the figure:

$\displaystyle e(n)^{2}\approx\frac{$1.9\text{\times}{10}^{8}$}{n}.$ (12)

While we have arrived at the factor $1.9\text{\times}{10}^{8}$ by fitting a
curve “by eye”, it turns out – maybe remarkably – that we can also
theoretically support this behavior: using the Markov chain central limit
theorem [34] (see Appendix C for details), we can estimate the mean of
$ne(n)^{2}$ by

$\textup{tr}\left(\mathop{\operator@font
diag}\nolimits(\langle\theta\rangle_{\pi(\theta|\hat{z})})^{-1}\cdot
I\\!AC\cdot\mathop{\operator@font
diag}\nolimits(\langle\theta\rangle_{\pi(\theta|\hat{z})})^{-1}\right)\approx
1.9\times 10^{8},$

where the matrix $I\\!AC$ is defined in (8).

Fig. 10: Convergence of the square of the relative error $e_{L}(n)^{2}$
between the running mean up to sample $n$ and the true mean
$\left<\theta\right>_{\pi(\theta|\hat{z})}$, measured in the weighted norm
(11), for a subset of chains. The thick red curve corresponds to the average
of these squared errors over all chains. This average squared error is
dominated by some chains with large errors and so lies above the curves
$e_{L}(n)^{2}$ of most chains. The thick orange line corresponds to the decay
$\frac{$1.9\text{\times}{10}^{8}$}{n}$ and represents an approximately average
convergence behavior of the chains.

If one measures computational effort by how many times an algorithm evaluates
the probability distribution $\pi(\theta|\hat{z})$, then $n$ in (11) can be
interpreted as work units and (12) provides an approximate relationship
between work and error. Similar relationships can be obtained experimentally
for other sampling algorithms that we hope this benchmark will be used for,
and (12) therefore allows comparison among other algorithms as well as against
the one used here.

##  4 Conclusions and what we hope this benchmark achieves

As the data presented in the previous section illustrates, it is possible to
obtain reasonably accurate statistics about the Bayesian solution of the
benchmark introduced in Section 2, even using a rather simple method: The
standard Metropolis-Hastings sampler. At the same time, using this method, it
is not at all trivial to compute posterior statistics accurately: We had to
compute $2\text{\times}{10}^{11}$ samples, and expended 30 CPU years on this
task (plus another two CPU years on postprocessing the samples).

But all of this also makes this a good benchmark: Simple algorithms, with
known performance, can solve it to a reasonable accuracy, and more advanced
algorithms should be able to do so in a fraction of time without making the
test case trivial. For example, it is not unreasonable to hope that advanced
sampling software [19, 1, 40, 37], using multi-level and multi-fidelity
expansions [41, 20, 45, 23], and maybe in conjunction with methods that
exploit the structure of the problem to approximate covariance matrices [55],
might be able to reduce the compute time by a factor of 100 to 1000, possibly
also running computations in parallel. This would move characterizing the
performance of such algorithms for the case at hand to the range of a few
hours or days on moderately parallel computers; practical computations might
not actually need the same level of accuracy and could be solved even more
rapidly.

As a consequence of these considerations, we hope that providing a benchmark
that is neither too simple nor too hard, and for which the solution is known
to good accuracy, spurs research in the development of better sampling
algorithms for Bayesian inverse problems. Many such algorithms of course
already exist, but in many cases, their performance is not characterized on
standardized test cases that would allow a fair comparison. In particular,
their performance is often characterized using probability distributions whose
characteristics have nothing to do with those that result from inverse
problems – say, sums of Gaussians. By providing a standardized benchmark that
matches what we expect to see in actual inverse problems – along with an open
source implementation of a code that computes the posterior probability
function $\pi(\theta|\hat{z})$ (see Appendix A) – we hope that we can
contribute to more informed comparisons between newly proposed algorithms:
Specifically, that their performance can be compared with the relationship
shown in (12) and Fig. 10 to provide a concrete factor of speed-up over the
method used here.

### Acknowledgments

W. Bangerth was partially supported by the National Science Foundation under
award OAC-1835673 as part of the Cyberinfrastructure for Sustained Scientific
Innovation (CSSI) program; by award DMS-1821210; by award EAR-1925595; and by
the Computational Infrastructure in Geodynamics initiative (CIG), through the
National Science Foundation under Award No. EAR-1550901 and The University of
California – Davis.

W. Bangerth also gratefully acknowledges the discussions and early experiments
with Kainan Wang many years ago. These early attempts directly led to the
ideas that were ultimately encoded in this benchmark. He also appreciates the
collaboration with Mantautas Rimkus and Dawson Eliasen on the SampleFlow
library that was used for the statistical evaluation of samples. Finally, the
feedback given by Noemi Petra, Umberto Villa, and Danny Long are acknowledged
with gratitude.

D. Aristoff gratefully acknowledges support from the National Science
Foundation via awards DMS-1818726 and DMS-2111277.

##  Appendix A An open source code to sample $\pi(\theta|\hat{z})$

We make a code that implements this benchmark available as part of the “code
gallery” for deal.II at
https://dealii.org/developer/doxygen/deal.II/code_gallery_MCMC_Laplace.html,
using the name MCMC-Laplace, and using the Lesser GNU Public License (LGPL)
version 2.1 or later as the license. deal.II is a software library that
provides the basic tools and building blocks for writing finite element codes
that solve partial differential equations numerically. More information about
deal.II is available at [5, 4]. The deal.II code gallery is a collection of
programs based on deal.II that were contributed by users as starting points
for others’ experiments.

The code in question has essentially three parts: (i) The forward solver that,
given a set of parameters $\theta$, produces the output $z^{\theta}$ using the
map discussed in Section 2.2.1; (ii) the statistical model that implements the
likelihood $L(z|\theta)$ and the prior probability $\pi_{\text{pr}}(\theta)$,
and combines these to the posterior probability $\pi(\theta|\hat{z})$; and
(iii) a simple Metropolis-Hastings sampler that draws samples from
$\pi(\theta|\hat{z})$. The second of these pieces is quite trivial,
encompassing only a couple of functions; we will therefore only comment on the
first and the third piece below.

###  A.1 Details of the forward solver

The forward solver is a C++ class whose main function performs the following
steps:

  1. 1.

It takes a 64-dimensional vector $\theta$ of parameter values, and interprets
it as the coefficients that describe a piecewise constant field
$a(\mathbf{x})$;

  2. 2.

It assembles a linear system that corresponds to the finite element
discretization of equations (2)–(3) using a $Q_{1}$ (bilinear) element on a
uniformly refined mesh;

  3. 3.

It solves this linear system to obtain the solution vector $U^{\theta}$ that
corresponds to the function $u_{h}^{\theta}(\mathbf{x})$; and

  4. 4.

It evaluates the solution $u_{h}^{\theta}$ at the measurement points
$\mathbf{x}_{k}$ to obtain $z^{\theta}$.

It then returns $z^{\theta}$ to the caller for evaluation with the statistical
model.

Such a code could be written in deal.II with barely more than 100 lines of C++
code, and this would have been sufficient for the purpose of evaluating new
ideas of sampling methods. However, we wanted to draw as large a number of
samples as possible, and consequently decided to see how fast we can make this
code.

To this end, we focused on accelerating three of the operations listed above,
resulting in a code that can evaluate $\pi(\hat{z}|\theta)$ in 4.5
$\mathrm{m}\mathrm{s}$ on an Intel Xeon E5-2698 processor with 2.20GHz (on
which about half of the samples used in this publication were computed), 3.1
$\mathrm{m}\mathrm{s}$ on an AMD EPYC 7552 processor with 2.2
$\mathrm{G}\mathrm{H}\mathrm{z}$ (the other half), and 2.7
$\mathrm{m}\mathrm{s}$ on an Intel Core i7-8850H processor with 2.6
$\mathrm{G}\mathrm{H}\mathrm{z}$ in one of the authors’ laptops.

The first part of the code that can be optimized for the current application
uses the fact that the linear system that needs to be assembled is the sum of
contributions from each of the cells of the mesh. Specifically, the
contribution from cell $K$ is

$\displaystyle A^{K}=P_{K}^{T}A^{K}_{\text{local}}P_{K}$

where $P_{K}$ is the restriction from the global set of degrees of freedom to
only those degrees of freedom that live on cell $K$, and
$A_{\text{local}}^{K}$ is – for the $Q_{1}$ element used here – a $4\times 4$
matrix of the form

$\displaystyle(A^{K}_{\text{local}})_{ij}=\int_{K}a^{\theta}(\mathbf{x})\nabla\varphi_{i}(\mathbf{x})\cdot\nabla\varphi_{j}(\mathbf{x})\,\text{d}x.$

This suggests that the assembly, including the integration above that is
performed via quadrature, has to be repeated every time we consider a new set
of parameters $\theta$. However, since we discretize $u_{h}$ on a mesh that is
a strict refinement of the one used for the coefficient
$a^{\theta}(\mathbf{x})$, and because $a^{\theta}(\mathbf{x})$ is piecewise
constant, we can note that

$\displaystyle(A^{K}_{\text{local}})_{ij}=\theta_{k(K)}\underbrace{\int_{K}\nabla\varphi_{i}(\mathbf{x})\cdot\nabla\varphi_{j}(\mathbf{x})\,\text{d}x}_{(A_{\text{local}})_{ij}},$

where $k(K)$ is the index of the element of $\theta$ that corresponds to cell
$K$. Here, the matrix $A_{\text{local}}$ no longer depends on $\theta$ and
can, consequently, be computed once and for all at the beginning of the
program. Moreover, $A_{\text{local}}$ does not actually depend on the cell $K$
as long as all cells have the same shape, as is the case here. We therefore
have to store only one such matrix. This approach makes the assembly
substantially faster since we only have to perform the local-to-global
operations corresponding to $P_{K}$ on every cell for every new $\theta$, but
no longer any expensive integration/quadrature.

Secondly, we have experimented with solving the linear systems so assembled as
fast as possible. For the forward solver used for each sample, the size of
these linear systems is $1089\times 1089$, with at most 9 entries per row.
Following a substantial amount of experimentation, we found that a sparse
direct solver is faster than any of the other approaches we have tried, and we
use the UMFPACK [18] interfaces in deal.II for this purpose. In particular,
this approach is faster than attempting to use an algebraic multigrid method
as a solver or preconditioner for the Conjugate Gradient method. We have also
tried to use a sparse decomposition via UMFPACK as a preconditioner for the CG
method, updating the decomposition only every few samples – based on the
assumption that the samples change only relatively slowly and so a
decomposition of the matrix for one sample is a good preconditioner for the
matrix corresponding to a subsequent sample. However, this turned out to be
slower than using a new decomposition for each sample.

The linear solver described above consumes about 90% of the time necessary to
evaluate each sample. As a consequence, there is certainly further room for
improvements. After all numerical results had been generated for this
publication, we have followed up on a suggestion by Martin Kronbichler to
replace the linear solver by a conjugate gradient method preconditioned by an
incomplete LU decomposition. This accelerates the computations by about a
factor of three, from $2.7\text{\,}\mathrm{m}\mathrm{s}$ to less than
$0.9\text{\,}\mathrm{m}\mathrm{s}$ per sample on the Intel Core i7-8850H
processor mentioned above. It is this accelerated version that is available at
the website mentioned above; we have verified that the new version results in
the same results up to at least 11 digits in computing the posterior
probability using the techniques mentioned in Appendix A.4.

Finally, evaluating the solution of a finite element field $u_{h}(\mathbf{x})$
at arbitrary points $\mathbf{x}_{k}$ is an expensive operation since one has
to find which cell $K$ the point belongs to and then transform this point into
the reference coordinate system of the cell $K$. On the other hand, the point
evaluation is a linear and bounded operation, and so there must exist a vector
$m_{k}$ so that $u_{h}(\mathbf{x}_{k})=m_{k}\cdot U$ where $U$ is the vector
of coefficients that describe $u_{h}$. This vector
$m_{k}=(\varphi_{i}(\mathbf{x}_{k}))_{i=0}^{1089-1}$ can be computed once and
for all. The computation of $z=(u_{h}(\mathbf{x}_{k}))_{k=0}^{169-1}$ can then
be facilitated by building a matrix $M$ whose rows are the vectors $m_{k}$,
and then the evaluation at all measurement points reduces to the operation
$z=MU$. $M$ is a sparse matrix with at most 4 entries per row, making this a
very economical approach.

The code with all of these optimizations is not very large – it contains 197
semicolons.777Counting semicolons is a commonly used metric in C and C++
programs. It roughly coincides with the number of declarations and statements
in a program, and is a better metric for code size than the number of lines of
code, as the latter also includes comments and empty lines used to help the
readability of a code.

###  A.2 Details of the Metropolis-Hastings sampler

The steps described at the start of this appendix yield an algorithm that,
given a sample $\theta$, can evaluate $\pi(\theta|\hat{z})$. We use this
functionality to drive a Metropolis-Hastings sampler to obtain a large number
of samples characterizing the posterior probability distribution.

While the basic algorithm of the Metropolis-Hastings sampler is well known
[32, 35], its practical implementation depends crucially on a number of
details that we will describe in the following.

First, we start the sampling process with a fixed sample
$\theta_{0}=(1,\ldots,1)^{T}$, corresponding to a coefficient
$a^{\theta}(\mathbf{x})=1$.

Secondly, an important step in the Metropolis-Hastings algorithm is the
generation of a “trial” sample $\tilde{\theta}_{\ell}$ based on the current
sample $\theta_{\ell}$. To this end, we use the following strategy: We define
the components of $\tilde{\theta}_{\ell,k}$ of $\tilde{\theta}_{\ell}$ as

$\displaystyle\tilde{\theta}_{\ell,k}=e^{\ln(\theta_{\ell,k})+\xi_{k}}=\theta_{\ell,k}e^{\xi_{k}},$

where $\xi_{k}$, $k=0,\ldots,63$, are iid Gaussians with mean $0$ and standard
deviation $\sigma_{\text{prop}}$. In other words, the “proposal distribution”
for the trial samples is an isotropic Gaussian ball centered at
$\theta_{\ell}$ in log space. This has the effect that all elements of samples
always stay positive, as one would expect given that they correspond to
material stiffness coefficients. The use of a ball in log space is also
consistent with the description of our prior probability distribution in
Section 2.2.3, which is also defined in log space.

To compute the Metropolis-Hastings acceptance probability, we first need to
compute the proposal probability density. By definition

$\displaystyle{\mathbb{P}}(\tilde{\theta}_{\ell,k}\leq\tilde{t}_{k}|\theta_{\ell,k}=t_{k})$
$\displaystyle={\mathbb{P}}(t_{k}e^{\xi_{k}}\leq\tilde{t}_{k})$
$\displaystyle={\mathbb{P}}(\xi_{k}\leq\log\tilde{t}_{k}-\log t_{k})$
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma_{\textup{prop}}^{2}}}\int_{-\infty}^{\log\tilde{t}_{k}-\log
t_{k}}\exp\left(\frac{-x^{2}}{2\sigma_{\textup{prop}}^{2}}\right)\,dx.$

The probability density
$p_{\textup{prop}}(\tilde{\theta}_{\ell,k}|\theta_{\ell,k})$ of proposing
$\tilde{\theta}_{\ell,k}$ given $\theta_{\ell,k}$ is then the derivative of
this expression with respect to $\tilde{t}_{k}$, with
$\tilde{\theta}_{\ell,k}$ and ${\theta}_{\ell,k}$ in place of $\tilde{t}_{k}$
and $t_{k}$:

$\displaystyle
p_{\textup{prop}}(\tilde{\theta}_{\ell,k}|\theta_{\ell,k})=\frac{1}{\tilde{\theta}_{\ell,k}}\frac{1}{\sqrt{2\pi\sigma_{\textup{prop}}^{2}}}\exp\left(\frac{-(\log\tilde{\theta}_{\ell,k}-\log\theta_{\ell,k})^{2}}{2\sigma_{\textup{prop}}^{2}}\right).$

By definition, the components $\tilde{\theta}_{\ell,k}$ of the proposal vector
$\tilde{\theta}_{\ell}$ are independent conditional on the current state
$\theta_{\ell}$. Thus the joint probability density,
$p_{\textup{prop}}(\tilde{\theta}_{\ell}|\theta_{\ell})$, of proposing vector
$\tilde{\theta}_{\ell}$ given vector $\theta_{\ell}$ is the product of the
probabilities above:

$p_{\textup{prop}}(\tilde{\theta}_{\ell}|\theta_{\ell})=\frac{1}{\prod_{k=0}^{63}\tilde{\theta}_{\ell,k}}\times\frac{1}{(2\pi\sigma_{\textup{prop}}^{2})^{32}}\exp\left(-\frac{\sum_{k=0}^{63}(\log\tilde{\theta}_{\ell,k}-\log\theta_{\ell,k})^{2}}{2\sigma_{\textup{prop}}^{2}}\right).$

The Metropolis-Hastings acceptance probability,
$A(\tilde{\theta}_{\ell}|\theta_{\ell})$, to accept proposal
$\tilde{\theta}_{\ell}$ given the current state $\theta_{\ell}$, is then

$\displaystyle A(\tilde{\theta}_{\ell}|\theta_{\ell})$
$\displaystyle=\min\left\\{1,\frac{\pi(\tilde{\theta}_{\ell}|\hat{z})}{\pi({\theta}_{\ell}|\hat{z})}\times\frac{p_{\textup{prop}}({\theta}_{\ell}|\tilde{\theta}_{\ell})}{p_{\textup{prop}}(\tilde{\theta}_{\ell}|\theta_{\ell})}\right\\}$
$\displaystyle=\min\left\\{1,\frac{\pi(\tilde{\theta}_{\ell}|\hat{z})}{\pi({\theta}_{\ell}|\hat{z})}\times\prod_{k=0}^{63}\frac{\tilde{\theta}_{\ell,k}}{{\theta}_{\ell,k}}\right\\}.$

As usual, with probability $A(\tilde{\theta}_{\ell}|\theta_{\ell})$ the
proposal $\tilde{\theta}_{\ell}$ is accepted, in which case it becomes the
next sample $\theta_{\ell+1}=\tilde{\theta}_{\ell}$; otherwise the proposal is
rejected and we keep the current sample, $\theta_{\ell+1}={\theta}_{\ell}$.

In our experiments, we use $\sigma_{\text{prop}}=0.0725$, corresponding to
changing $\theta$ by a factor $e^{\xi}$ that with 65% probability lies within
the range
$[e^{-\sigma_{\text{prop}}},e^{+\sigma_{\text{prop}}}]=[0.93,1.075]$. This
results in an acceptance rate for the Metropolis-Hastings algorithm of just
under 24%. This is close to the recommended value of 0.234 for Metropolis-
Hastings sampling algorithms that can be derived for specific probability
distributions that are generally simpler than the one we are interested in
here [27, 46]; the theory guiding the derivation of the 0.234 value may not be
applicable here (see [10]), but absent better guidance, we stuck with it.

Finally, all steps in the Metropolis-Hastings algorithms that require the use
of a random number use the MT19937 random number generator, as implemented by
C++11 compilers in the std::mt19937 class.

###  A.3 Implementations of the benchmark in alternative languages

We consider the C++ implementation discussed above as the “reference
implementation” of the benchmark. However, we recognize that it is rather
heavy-weight in the sense that it requires the installation of the deal.II
library. While this readily facilitates otherwise non-trivial modifications
(e.g., for multilevel sampling schemes that require the solution of the
forward problem on coarser meshes), it is clear that for simple experiments,
it would be nice to have stand-alone implementations of the benchmark.

As a consequence, we have also developed Matlab and Python versions of the
benchmark. These are available from the same website from which the C++
implementation is available – see the link at the top of the appendix. These
alternative implementations provide everything one needs to build a sampler:
Namely, the functionality to provide an input vector $\theta$ and to compute
the prior $\pi_{\text{pr}}(\theta)$ and the likelihood $L(\hat{z}|\theta)$
from such an input. Using (1), these can then be used to compute the posterior
probability $\pi(\theta|\hat{z})$ that forms the basis of most sampling
algorithms. The Matlab version includes a basic Metropolis-Hastings sampler
with parallel functionality.

On a recent laptop, the Matlab version is able to compute the posterior
probability for a sample in about $4\text{\,}\mathrm{m}\mathrm{s}$, not
substantially slower than the C++ version we have used for our results. The
Python version – although an almost literal transcription of the Matlab
version – requires approximately $25\text{\,}\mathrm{m}\mathrm{s}$. One
imagines that it could be optimized further (for example, a substantial part
of the run time is spent in the insertion of the cell-local $4\times 4$
matrices into the global matrix), but we have not attempted to do so.

###  A.4 Testing of alternative implementations

In order to facilitate testing of alternative implementations such as the ones
discussed in the previous section (or testing modifications made to the C++
implementation itself), the website from which the benchmark can be obtained
(see the top of this appendix) also contains a set of known input/output
pairs. Specifically, it contains files for ten different input vectors
$\theta$, along with the corresponding outputs $z^{\theta}$, likelihood
$L(\hat{z}|\theta)$, and prior $\pi_{\text{pr}}(\theta)$ for each input
vector. The latter two can be combined via (1) to the posterior probability
$\pi(\theta|\hat{z})$ associated with the input vector $\theta$.

We have used these known input/output pairs obtained from our reference
implementation to verify that the alternative implementations of the benchmark
discussed in the previous sub-section are correct. For example, the Matlab
implementation provides the vectors $z^{\theta}$ to relative errors on the
order of $10^{-13}$; log priors and likelihoods are computed to relative
errors less than $10^{-11}$. The Python version achieves the same level of
accuracy.

##  Appendix B One dimensional version of the benchmark

Many of the features of the posterior probability on the 64-dimensional
parameter space that we have experimentally observed in Section 3 match those
that one would expect for inverse problems of the kind discussed herein. In
particular, the fact that we have a large spread between the large and small
eigenvalues of the covariance matrix, the non-Gaussianity of the probability
distribution, and the anticorrelation of parameters defined on neighboring
cells did not come as a surprise. Yet, strict proofs for these properties are
hard to come by.

At the same time, we can investigate an analogous situation in a one-
dimensional model with parameters $\theta_{k}$, for which one can derive the
posterior probability distribution analytically. In this Appendix, we work out
some details of a one-dimensional version of the benchmark. Consider the
generalized Poisson equation

$\displaystyle-\frac{d}{dx}\left(a(x)\frac{du}{dx}(x)\right)$
$\displaystyle=f(x),$ $\displaystyle 0<x<1,$ (13) $\displaystyle u(x)$
$\displaystyle=0,$ $\displaystyle x\in\\{0,1\\}.$ (14)

Again we assume $a(x)=a^{\theta}(x)$ is parametrized by
$\theta=(\theta_{0},\ldots,\theta_{N-1})$ where

$a(x)=\theta_{k},\qquad\text{if}\quad\frac{k}{N}<x<\frac{k+1}{N},$ (15)

and we take $f(x)\equiv 1$. The solution to (13) is then of the form

$u(x)=-\int_{0}^{x}\frac{y+c}{a(y)}\,dy,$ (16)

where $c$ is constant determined by requiring $u(1)=0$. Due to the piecewise
constant form of $a=a(x)$, this rewrites as

$u(x)=u_{k}(x):=-\frac{1}{\theta_{k}}\left(\frac{x^{2}}{2}+cx\right)+d_{k},\qquad\text{if}\quad\frac{k}{N}<x<\frac{k+1}{N}.$
(17)

There are $N-1$ continuity constraints and two boundary conditions, namely,

$\displaystyle
u_{k-1}\left(\frac{k}{N}\right)=u_{k}\left(\frac{k}{N}\right),\qquad
k=1,\ldots,N-1,$ (18) $\displaystyle u_{0}(0)=0,\quad u_{N-1}(1)=0.$ (19)

This translates, via (17), into $N+1$ linear equations for the $N+1$
coefficients $c$,$d_{0},\ldots,d_{N-1}$. The boundary conditions show that

$d_{0}=0,\qquad
c=-\frac{\int_{0}^{1}\frac{y}{a(y)}\,dy}{\int_{0}^{1}\frac{1}{a(y)}\,dy}=-\frac{1}{2N}\frac{\sum_{k=0}^{N-1}\theta_{k}^{-1}\left(2k+1\right)}{\sum_{k=0}^{N-1}\theta_{k}^{-1}}.$
(20)

The remaining $N-1$ equations for $d_{1},\ldots,d_{N-1}$ can be solved using
the continuity constraints, which give the equations

$d_{k-1}-d_{k}=\left(\frac{1}{2}\left(\frac{k}{N}\right)^{2}+c\frac{k}{N}\right)\left(\theta_{k-1}^{-1}-\theta_{k}^{-1}\right),\qquad
k=1,\ldots,N-1,$ (21)

which have solution

$d_{k}=-\sum_{j=1}^{k}\left(\frac{1}{2}\left(\frac{j}{N}\right)^{2}+c\frac{j}{N}\right)\left(\theta_{j-1}^{-1}-\theta_{j}^{-1}\right),\qquad
k=1,\ldots,N-1.$

Let us specifically consider the case with $N=2$ parameters to understand some
qualitative features of the benchmark. In this case,
$c=-\frac{1}{4}(3\theta_{0}+\theta_{1})/(\theta_{0}+\theta_{1})$ and

$\displaystyle u_{0}(x)$
$\displaystyle=-\frac{x^{2}}{2\theta_{0}}+\frac{3\theta_{0}+\theta_{1}}{4\theta_{0}(\theta_{0}+\theta_{1})}x,$
$\displaystyle\quad u_{1}(x)$
$\displaystyle=-\frac{x^{2}}{2\theta_{1}}+\frac{3\theta_{0}+\theta_{1}}{4\theta_{1}(\theta_{0}+\theta_{1})}x+\frac{\theta_{1}-\theta_{0}}{4\theta_{1}(\theta_{0}+\theta_{1})}.$

This solution is shown in Fig. 11.

Let us assume that we have measurements at $x_{0}=0.25$ and $x_{1}=0.75$. Then
the “exact” measurements we would get are888Unlike in Section 2, these values
are computed using the exact solution, instead of a finite element
approximation.

$\displaystyle\hat{z}_{0}$
$\displaystyle=-\frac{x_{0}^{2}}{2\theta_{0}}+\frac{3\theta_{0}+\theta_{1}}{4\theta_{0}(\theta_{0}+\theta_{1})}x_{0},$
$\displaystyle\qquad\hat{z}_{1}$
$\displaystyle=-\frac{x_{1}^{2}}{2\theta_{1}}+\frac{3\theta_{0}+\theta_{1}}{4\theta_{1}(\theta_{0}+\theta_{1})}x_{1}+\frac{\theta_{1}-\theta_{0}}{4\theta_{1}(\theta_{0}+\theta_{1})}.$

We can use these values to exactly and cheaply evaluate the posterior
probability (without sampling), and Fig. 12 shows $\pi(\theta|\hat{z})$ in
this simple case, using “true” parameters $\hat{\theta}_{0}=0.1$ and
$\hat{\theta}_{1}=1$ to define $\hat{z}$. We use the prior and likelihoods
defined before, with the prior standard deviation $\sigma_{\textup{pr}}=2$ and
two different values for the likelihood standard deviations used in (5):
$\sigma=0.01$ and $\sigma=0.1$.

Fig. 11: Exact solution to (13)–(14) when $a(x)=0.1$ for $0<x<1/2$ and
$a(x)=1$ for $1/2<x<1$, corresponding to “true” parameters
$\hat{\theta}_{0}=0.1$ and $\hat{\theta}_{1}=1$, with measurements at
$x_{0}=0.25$ and $x_{1}=0.75$.

Fig. 12: Left: Non-normalized posterior distribution $\pi(\theta|\hat{z})$
when $\sigma=0.01$ in the 1D model. Notice the distribution is approximately
Gaussian for this very small $\sigma$. Right: When using $\sigma=0.1$ in the
1D model. For this larger value of $\sigma$, the posterior distribution
extends to larger values of $\theta_{0}$ and $\theta_{1}$, and has a non-
Gaussian, “banana” shape, reminiscent of the pair histograms in Fig. 7.

The figure illustrates more concisely the strong correlation between
parameters that we have experimentally observed in Fig. 7. It also illustrates
that if $\sigma$ in the likelihood (5) is chosen small (i.e., if the
measurement error is small), then the posterior is approximately Gaussian. The
observation here therefore validates our choice of a relatively large
$\sigma$, see also Remark 3.

##  Appendix C Estimating the essential sample size

Section 3.1 assessed how much information is actually present in the many
samples we have computed, and we have also used the results shown there in
providing theoretical support for the key cost-accuracy estimate (12). Many of
the statistical errors described there can be estimated from the
autocovariance $AC_{L}(s)$ defined in (7). The basis for this is the Markov
chain central limit theorem. Informally, the Markov chain central limit
theorem says that, for large $n$ and $N_{L}$, the running means
$\langle\theta\rangle_{L,n}$ are approximately normally distributed with mean
$\langle\theta\rangle_{\pi(\theta|\hat{z})}$ and covariance

$\textup{Cov}\left(\langle\theta\rangle_{L,n}\right)\approx\frac{I\\!AC}{n},$
(22)

where $I\\!AC$ is the integrated autocovariance, obtained by summing up the
autocovariance $AC_{L}(s)$; we estimate $I\\!AC$ from the data in (8).
Equation (22) then justifies the formula for the scaling of $e(n)^{2}$ below
(12).

To establish the effective sample size formula (9), we cite the “Delta method”
in statistics [58]. Informally, the Delta method states that, for large $n$
and $N_{L}$, a continuously differentiable function $f$ of the running means
is nearly normally distributed with variance

$\textup{Var}(f(\langle\theta\rangle_{L,n}))\approx\frac{1}{n}\nabla
f(\langle\theta\rangle_{\pi(\theta|\hat{z})})^{T}\cdot I\\!AC\cdot\nabla
f(\langle\theta\rangle_{\pi(\theta|\hat{z})}),$ (23)

provided $\nabla f(\langle\theta\rangle_{\pi(\theta|\hat{z})})$ is nonzero.
The formula (23) is obtained by Taylor expanding $f$ and applying the Markov
chain central limit theorem. For a Markov chain in which successive samples
are all independent, $I\\!AC$ is the covariance matrix, which we estimate from
the data in (10) and denote by $C$. Using a standard result on generalized
Rayleigh quotients,

$\displaystyle\frac{\nabla
f(\langle\theta\rangle_{\pi(\theta|\hat{z})})^{T}\cdot I\\!AC\cdot\nabla
f(\langle\theta\rangle_{\pi(\theta|\hat{z})})}{\nabla
f(\langle\theta\rangle_{\pi(\theta|\hat{z})})^{T}\cdot C\cdot\nabla
f(\langle\theta\rangle_{\pi(\theta|\hat{z})})}\leq\max_{x\neq
0}\frac{x^{T}\cdot I\\!AC\cdot x}{x^{T}\cdot C\cdot
x}=\lambda_{\max}(C^{-1}\cdot I\\!AC),$

where $\lambda_{\max}$ denotes the largest eigenvalue. This means that the
variance (23) is at most $\lambda_{\max}(C^{-1}\cdot I\\!AC)$ times the value
it would take if all of the samples on chain $L$ were independent. In other
words, the minimum number of “effectively independent samples” is
approximately $N_{L}/\lambda_{\max}(C^{-1}\cdot I\\!AC)$ for each of our
chains of length $N_{L}$.

It is quite standard to quantify statistical error in Markov chain Monte Carlo
simulations by using the integrated autocovariance [47, 28]. In the
literature, it is also common to determine an effective sample size by
checking to see where autocovariances cross a certain threshold – such as a
small percentage of its initial value – as we have done in Figure 3 above,
where we used a threshold of $1\%$.

## References

  * [1] B. M. Adams et al., Dakota, a multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantification, and sensitivity analysis: Version 6.13 user’s manual, Tech. Report SAND2020-12495, Sandia National Laboratories, 2020.
  * [2] O. Aguilar, M. Allmaras, W. Bangerth, and L. Tenorio, Statistics of parameter estimates: A concrete example, SIAM Review, 57 (2015), pp. 131–149.
  * [3] M. Allmaras, W. Bangerth, J. M. Linhart, J. Polanco, F. Wang, K. Wang, J. Webster, and S. Zedler, Estimating parameters in physical models through bayesian inversion: A complete example, SIAM Review, 55 (2013), pp. 149–167.
  * [4] D. Arndt, W. Bangerth, B. Blais, T. C. Clevenger, M. Fehling, A. V. Grayver, T. Heister, L. Heltai, M. Kronbichler, M. Maier, P. Munch, J.-P. Pelteret, R. Rastak, I. Thomas, B. Turcksin, Z. Wang, and D. Wells, The deal.II library, version 9.2, Journal of Numerical Mathematics, 28 (2020), pp. 131–146.
  * [5] D. Arndt, W. Bangerth, D. Davydov, T. Heister, L. Heltai, M. Kronbichler, M. Maier, J.-P. Pelteret, B. Turcksin, and D. Wells, The deal.II finite element library: design, features, and insights, Computers & Mathematics with Applications, 81 (2021), pp. 407–422.
  * [6] Y. F. Atchadé and J. S. Rosenthal, On adaptive Markov Chain Monte Carlo algorithms, Bernoulli, 11 (2005), pp. 815–828.
  * [7] J. M. Bardsley, Computational Uncertainty Quantification for Inverse Problems, SIAM, 2018.
  * [8] J. M. Bardsley, A. Seppänen, A. Solonen, H. Haario, and J. Kaipio, Randomize-then-optimize for sampling and uncertainty quantification in electrical impedance tomography, SIAM/ASA J. Uncertainty Quantification, 3 (2015), pp. 1136–1158.
  * [9] J. M. Bardsley, A. Solonen, H. Haario, and M. Laine, Randomize-then-optimize: A method for sampling from posterior distributions in nonlinear inverse problems, SIAM J. Sci. Comput., 36 (2014), pp. A1895–A1910.
  * [10] M. Bédard, Optimal acceptance rates for Metropolis algorithms: Moving beyond 0.234, Stochastic Processes and their Applications, 118 (2008), pp. 2198–2222.
  * [11] J. Bierkens, P. Fearnhead, and G. O. Roberts, The zig-zag process and super-efficient sampling for Bayesian analysis of big data, The Annals of Statistics, 47 (2019), pp. 1288–1320.
  * [12] T. Bui-Thanh and O. Ghattas, Analysis of the hessian for inverse scattering problems: II. Inverse medium scattering of acoustic waves, Inverse Problems, 28 (2012), p. 055002.
  * [13] P. Chen, U. Villa, and O. Ghattas, Taylor approximation and variance reduction for PDE-constrained optimal control under uncertainty, Journal of Computational Physics, 385 (2019), pp. 163–186.
  * [14] J. A. Christen and C. Fox, Markov chain monte carlo using an approximation, Journal of Computational and Graphical Statistics, 14 (2005), pp. 795–810.
  * [15] I. Craig and J. Brown, Inverse problems in astronomy, in Bayesian Astrophysics, A. Asensio Ramos and I. Arregui, eds., Cambridge University Press, 1986, pp. 31–61.
  * [16] E. Darve, D. Rodríguez-Gómez, and A. Pohorille, Adaptive biasing force method for scalar and vector free energy calculations, The Journal of chemical physics, 128 (2008), p. 144120.
  * [17] M. Dashti and A. M. Stuart, The bayesian approach to inverse problems, 2015.
  * [18] T. A. Davis, Algorithm 832, ACM Transactions on Mathematical Software, 30 (2004), pp. 196–199.
  * [19] B. Debusschere, K. Sargsyan, C. Safta, and K. Chowdhary, Uncertainty quantification toolkit (UQTk), in Handbook of Uncertainty Quantification, R. Ghanem, D. Higdon, and H. Owhadi, eds., Springer International Publishing, Cham, 2017, pp. 1807–1827.
  * [20] T. J. Dodwell, C. Ketelsen, R. Scheichl, and A. L. Teckentrup, A hierarchical multilevel Markov Chain Monte Carlo algorithm with applications to uncertainty quantification in subsurface flow, SIAM/ASA Journal on Uncertainty Quantification, 3 (2015), pp. 1075–1108.
  * [21] A. B. Duncan, T. Lelievre, and G. A. Pavliotis, Variance reduction using nonreversible Langevin samplers, Journal of Statistical Physics, 163 (2016), pp. 457–491.
  * [22] 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.
  * [23] C. M. Fleeter, G. Geraci, D. E. Schiavazzi, A. M. Kahn, and A. L. Marsden, Multilevel and multifidelity uncertainty quantification for cardiovascular hemodynamics, Computer Methods in Applied Mechanics and Engineering, 365 (2020), p. 113030.
  * [24] Daniel Foreman-Mackey, David W Hogg, Dustin Lang, and Jonathan Goodman, emcee: the mcmc hammer, Publications of the Astronomical Society of the Pacific, 125 (2013), p. 306.
  * [25] C. Fox and G. Nicholls, Sampling conductivity images via MCMC, in Proceedings of the Leeds Annual Statistical Research Workshop (LASR), K.V. Mardia, C.A. Gill, and R.G. Aykroyd, eds., 1997, pp. 91–100.
  * [26] N. Friedman, M. Linial, I. Nachman, and D. Pe’er, Using Bayesian networks to analyze expression data, Journal of Computational Biology, 7 (2000), pp. 601–620.
  * [27] A. Gelman, W. R. Gilks, and G. O. Roberts, Weak convergence and optimal scaling of random walk Metropolis algorithms, Annals of Applied Probability, 7 (1997), pp. 110–120.
  * [28] C. J. Geyer, Practical Markov Chain Monte Carlo, Statistical science, (1992), pp. 473–483.
  * [29] Jonathan Goodman and Jonathan Weare, Ensemble samplers with affine invariance, Communications in applied mathematics and computational science, 5 (2010), pp. 65–80.
  * [30] H. Haario, M. Laine, A. Mira, and E. Saksman, DRAM: efficient adaptive MCMC, Statistics and Computing, 16 (2006), pp. 339–354.
  * [31] H. Haario, E. Saksman, and J. Tamminen, An adaptive Metropolis algorithm, Bernoulli, 7 (2001), pp. 223 – 242.
  * [32] W. K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika, 57 (1970), pp. 97–109.
  * [33] Y. Jiang and A. D. Woodbury, A full-Bayesian approach to the inverse problem for steady-state groundwater flow and heat transport, Geophysical Journal International, 167 (2006), pp. 1501–1512.
  * [34] G. L. Jones, On the Markov chain central limit theorem, Probability surveys, 1 (2004), pp. 299–320.
  * [35] J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Springer-Verlag, 2005.
  * [36] J. Martin, L. C. Wilcox, C. Burstedde, and O. Ghattas, A stochastic newton MCMC method for large-scale statistical inverse problems with application to seismic inversion, SIAM Journal on Scientific Computing, 34 (2012), pp. A1460–A1487.
  * [37] D. McDougall, N. Malaya, and R. D. Moser, The parallel C++ statistical library for Bayesian inference: QUESO, in Handbook of Uncertainty Quantification, R. Ghanem, D. Higdon, and H. Owhadi, eds., Springer International Publishing, Cham, 2017, pp. 1829–1865.
  * [38] R. Neal, MCMC using Hamiltonian dynamics, in Handbook of Markov Chain Monte Carlo, S. Brooks, A. Gelman, G. L. Jones, and X.-L. Meng, eds., CRC Press, 2011, pp. 113–162.
  * [39] R. M. Neal, Bayesian learning via stochastic dynamics, in Advances in neural information processing systems, 1993, pp. 475–482.
  * [40] M. Parno, A. Davis, P. Conrad, and Y. M. Marzouk, Mit uncertainty quantitification (muq) library.  http://muq.mit.edu/, 2021.
  * [41] B. Peherstorfer, K. Willcox, and M. Gunzburger, Survey of multifidelity methods in uncertainty propagation, inference, and optimization, SIAM Review, 60 (2018), pp. 550–591.
  * [42] N. Petra, J. Martin, G. Stadler, and O. Ghattas, A computational framework for infinite-dimensional bayesian inverse problems, part II: Stochastic Newton MCMC with application to ice sheet flow inverse problems, SIAM Journal on Scientific Computing, 36 (2014), pp. A1525–A1555.
  * [43] G. O. Roberts and J. S. Rosenthal, Examples of adaptive MCMC, Journal of Computational and Graphical Statistics, 18 (2009), pp. 349–367.
  * [44] G. O. Roberts and R. L. Tweedie, Exponential convergence of Langevin distributions and their discrete approximations, Bernoulli, 2 (1996), pp. 341–363.
  * [45] J. Seo, C. Fleeter, A. M. Kahn, A. L. Marsden, and D. E. Schiavazzi, Multi-fidelity estimators for coronary artery circulation models under clinically-informed data uncertainty, International Journal for Uncertainty Quantification, 10 (2020), pp. 449–466.
  * [46] C. Sherlock and G. O. Roberts, Optimal scaling of the random walk Metropolis on elliptically symmetric unimodal targets, Bernoulli, 15 (2009), pp. 774–798.
  * [47] A. Sokal, Monte Carlo methods in statistical mechanics: foundations and new algorithms, in Functional Integration, Springer, 1997, pp. 131–192.
  * [48] A. M. Stuart, Inverse problems: A bayesian perspective, Acta Numerica, 19 (2010), pp. 451–559.
  * [49] A. Tarantola, Inverse Problem Theory, Elsevier, Amsterdam, New York, 1987.
  * [50] , Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics, Jan. 2005.
  * [51] C. J. F. Ter Braak, A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces, Statistics and Computing, 16 (2006), pp. 239–249.
  * [52] C. J. F. ter Braak and J. A. Vrugt, Differential Evolution Markov Chain with snooker updater and fewer chains, Statistics and Computing, 18 (2008), pp. 435–446.
  * [53] L. Tierney and A. Mira, Some adaptive Monte Carlo methods for Bayesian inference, Statistics in Medicine, 18 (1999), pp. 2507–2515.
  * [54] P. Vanetti, A. Bouchard-Côté, G. Deligiannidis, and A. Doucet, Piecewise-deterministic Markov Chain Monte Carlo, arXiv preprint arXiv:1707.05296, (2017).
  * [55] U. Villa, N. Petra, and O. Ghattas, hIPPYlib: An extensible software framework for large-scale inverse problems governed by PDEs. Part I: Deterministic inversion and linearized bayesian inference, ACM Transactions on Mathematical Software, (2021, accepted).
  * [56] J. A. Vrugt, C. J. F. ter Braak, C. G. H. Diks, B. A. Robinson, J. M. Hyman, and D. Higdon, Accelerating Markov Chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling, International Journal of Nonlinear Sciences and Numerical Simulation, 10 (2009), pp. 273–290.
  * [57] D. Watzenig and C. Fox, A review of statistical modelling and inference for electrical capacitance tomography, Measurement Science and Technology, 20 (2009), p. 052002.
  * [58] K. Wolter, Introduction to variance estimation, Springer Science & Business Media, 2007.
  * [59] J. Worthen, G. Stadler, N. Petra, M.l Gurnis, and O. Ghattas, Towards adjoint-based inversion for rheological parameters in nonlinear viscous mantle flow, Physics of the Earth and Planetary Interiors, 234 (2014), pp. 23–34.