##  Chapter 16 Posterior exploration for computationally intensive forward
models

Mikkel B. Lykkegaard, Colin Fox, Dave Higdon, C. Shane Reese,
and J. David Moulton

To appear in the Handbook of Markov Chain Monte Carlo (2nd edition),
probably as Chapter 16.

###  1 Introduction

In a common inverse problem, we wish to infer about an unknown spatial field
$x=(x_{1},\ldots,x_{m})^{T}$ given indirect observations
$y=(y_{1},\ldots,y_{n})^{T}$. The observations, or data, are linked to the
unknown field $x$ through some physical system

$y=\zeta(x)+\epsilon$

where $\zeta(x)$ denotes the actual physical system and $\epsilon$ is an
$n$-vector of observation errors. Examples of such problems include medical
imaging (Kaipio and Somersalo,, 2004), geologic and hydrologic inversion
(Stenerud et al.,, 2008), and cosmology (Jimenez et al.,, 2004). When a
forward model, or simulator $\eta(x)$, is available to model the physical
process, one can model the data using the simulator

$y=\eta(x)+e,$

where $e$ includes observation error as well as error due to the fact that the
simulator $\eta(x)$ may be systematically different from reality $\zeta(x)$
for input condition $x$. Our goal is to use the observed data $y$ to make
inference about the spatial input parameters $x$ – predict $x$ and
characterize the uncertainty in the prediction for $x$.

The likelihood $L(y|x)$ is then specified to account for both mismatch and
sampling error. We will assume zero-mean Gaussian errors so that

$L(y|x)\propto\exp\\{-\mbox{\small$\frac{1}{2}$}(y-\eta(x))^{T}\Sigma_{e}^{-1}(y-\eta(x))\\},$
(1.1)

with $\Sigma_{e}$ known. It is worth noting that the data often come from only
a single experiment. So while it is possible to quantify numerical errors,
such as those due to discretization (e.g.,Kaipio and Somersalo, (2004);
Nissinen et al., (2008)), there is no opportunity to obtain data from
additional experiments for which some controllable inputs have been varied.
Because of this limitation, there is little hope of determining the sources of
error in $e$ due to model inadequacy. Therefore, the likelihood specification
will often need to be done with some care, incorporating the modeler’s
judgment about the appropriate size and nature of the mismatch term.

In many inverse problems we wish to reconstruct $x$, an unknown process over a
regular 2-d lattice. We consider systems for which the model input parameters
$x$ denote a spatial field or image. The spatial prior is specified for $x$,
$\pi(x)$ which typically takes into account modeling, and possibly
computational considerations.

The resulting posterior is then given by

$\pi(x|y)\propto L(y|\eta(x))\times\pi(x).$

This posterior can, in principle, be explored via Markov chain Monte Carlo
(MCMC). However the combined effects of the high dimensionality of $x$ and the
computational demands of the simulator make implementation difficult, and
often impossible, in practice. By itself, the high dimensionality of $x$ isn’t
necessarily a problem. MCMC has been carried out with relative ease in large
image applications (Weir,, 1997; Rue,, 2001). However, in these examples, the
forward model was either trivial, or non-existent. Unfortunately, even a
mildly demanding forward simulation model can greatly affect the feasibility
of doing MCMC to solve the inverse problem.

In this chapter we apply a standard single-site updating scheme that dates
back to Metropolis et al., (1953) to sample from this posterior. While this
approach has proven effective in a variety of applications, it has the
drawback of requiring millions of calls to the simulation model. In Section 3
we consider an MCMC scheme that uses highly multivariate updates to sample
from $\pi(x|y)$: the multivariate random walk Metropolis algorithm (Gelman et
al.,, 1996). Such multivariate updating schemes are alluring for
computationally demanding inverse problems since they have the potential to
update many (or all) components of $x$ at once, while requiring only a single
evaluation of the simulator.

Next, in Section 4, we consider augmenting the basic posterior formulation
with additional formulations based on faster, approximate simulators. We
create two faster, approximate forward models: $\eta_{\text{m}}(x)$ by
leveraging an incomplete multigrid solve; and $\eta_{\text{c}}(x)$ by
coarsening the initial “fine-scale” formulation described in the next section.
These approximate simulators can be used in a delayed acceptance scheme (Liu,,
2001; Christen and Fox,, 2005), as well as in an augmented formulation (Higdon
et al.,, 2002). Both of these recipes can be utilized with any of the above
MCMC schemes, often leading to substantial improvements in efficiency. For a
given $x$, resulting output produced by the approximate multigrid model
$\eta_{\text{m}}(x)$ is quite close to that obtained by exact solve of the
fine scale model $\eta(x)$. However, the coarser model $\eta_{\text{c}}(x)$
shows substantial, systematic differences from the output of the fine-scale
model $\eta(x)$. So we also consider the stochastic approaches for correcting
systematic errors between the coarse- and fine-scale model output, introduced
in Cui et al., (2011). We find that the corrected coarse model enables cheap
and semi-automatic generation of multivariate proposals from the single-site
proposal, with an appreciable acceptance rate in the fine-scale chain, giving
improved efficiency.

The updating schemes are illustrated with an electrical impedance tomography
(EIT) application described in the next section, where the values of $x$
denote electrical conductivity of a 2-d object. The chapter concludes with a
discussion and some general recommendations.

###  2 An inverse problem in electrical impedance tomography

Bayesian methods for EIT applications have been described in Fox and Nicholls,
(1997), Kaipio et al., (2000) and Andersen et al., (2003). A notional inverse
problem is depicted in Figure 16.1; this setup was given previously in Moulton
et al., (2008). Here a 2-d object composed of regions with differing
electrical conductivity is interrogated by $16$ electrodes. From each
electrode, in turn, a current $I$ is injected into the object and taken out at
a rate of $I/(16-1)$ at the remaining 15 electrodes. The voltage is then
measured at each of the 16 electrodes. These 16 experimental configurations
result in $n=16\times 16$ voltage observations which are denoted by the
$n$-vector $y$. The measurement error is simulated by adding iid mean 0
Gaussian noise to each of the voltage measurements. The standard deviation
$\sigma$ of this noise is chosen so that the signal to noise ratio is about
1000:3, which is typical of actual EIT measurements. The resulting simulated
data is shown in the right frame of Figure 16.1 – one plot for each of the 16
circuit configurations. In each of those plots, the injector electrode is
denoted by the black plotting symbol.

Figure 16.1:  A synthetic EIT application. A 2-d object is surrounded by
electrodes at 16 evenly spaced locations around its edge. The conductivity of
the object is 3 in the white regions, and 4 in the black regions (the units
are arbitrary since the data are invariant to scaling of the conductivity).
First, a current of $I$ is injected at electrode 1, and extracted evenly at
the other 15 electrodes. The voltage is measured at each electrode, with
respect to the mean voltage on electrodes. This data is shown in the plot
labeled 1 on the right. Similar experiments are carried out with each
electrode taking a turn as the injector. The resulting voltages are shown in
the remaining 15 plots. In each plot, the voltage corresponding to the
injector electrode is given by a black plotting symbol.

We take $s$ to denote spatial locations within the object
$\Omega=[0,1]\times[0,1]$, and take $x(s)$ to denote the electrical
conductivity at site $s$. We also take $v(s)$ to be the potential at location
$s$, and $j(s)$ to be the current at boundary location $s$. A mathematical
model for the measurements is then the Neumann boundary-value problem

$\begin{array}[c]{rlll}{-\nabla}\cdot x\left(s\right)\nabla
v\left(s\right)&=&0&\qquad s\in\Omega\\\ \displaystyle
x\left(s\right)\frac{\partial v\left(s\right)}{\partial
n\left(s\right)}&=&j\left(s\right)&\qquad s\in\partial\Omega,\end{array}$

where $\partial\Omega$ denotes the boundary of the object $\Omega$ and $n(s)$
is the unit normal vector at the boundary location $s\in\partial\Omega$. The
conservation of current requires that the sum of the currents at each of the
16 electrodes be 0. The arbitrary additive offset in the solution of this
Neumann problem is determined by requiring that the sum of voltages at
boundary electrodes is zero.

In order to numerically solve this problem for a given set of currents at the
electrodes and a given conductivity field, $x(s)$, the conductivity field is
discretized into an $m=24\times 24$ lattice. We use a standard Cholesky-based
sparse solver for this system. Later in Section 4, we also create a faster,
approximate solver based on an $8\times 8$ discretization of the conductivity
field.

Now, for any specified conductivity configuration $x$ and current
configuration, the solver produces 16 voltages. For all 16 current
configurations, the forward solver produces an $n=256$-vector of resulting
voltages $\eta(x)$. Hence the sampling model for the data $y$ given the
conductivity field $x$ is given by (1.1) where $\Sigma_{e}=\sigma^{2}I_{n}$.

For the conductivity image prior, we adapt a Markov random field (MRF) prior
from Geman and McClure, (1987). This prior has the form

$\pi(x)\propto\exp\left\\{\beta\sum_{i\sim
j}u(x_{i}-x_{j})\right\\},\,x\in[2.5,4.5]^{m}$ (2.1)

where $\beta$ and $s$ control the regularity of the field, and $u(\cdot)$ is
the tricube function of Cleveland, (1979)

$u(d)=\left\\{\begin{array}[]{cl}\displaystyle\frac{1}{s}(1-\left|d/s\right|^{3})^{3}&\mbox{
if }-s<d<s\\\ 0&\mbox{ if }|d|\geq s\,.\end{array}\right.$

The sum is over all horizontal and vertical nearest neighbors denoted by
$i\sim j$ and given by the edges in the Markov random field graph in Figure
16.2. Hence this prior encourages neighboring $x_{i}$’s to have similar
values, but once $x_{i}$ and $x_{j}$ are more than $s$ apart, the penalty does
not grow. This allows occasional large shifts between neighboring $x_{i}$’s.
For this chapter, we fix $(\beta,s)=(.5,.3)$. A realization from this prior is
shown in the right frame of Figure 16.2. A typical prior realization shows
patches of homogeneous values, along with abrupt changes in intensity at each
patch boundary. This prior also allows an occasional, isolated, extreme single
pixel value.

Figure 16.2:  Left: the first order neighborhood MRF graph corresponding to
the prior in (2.1); each term in the sum corresponds to an edge in the MRF
graph. Right: a realization from this gray level prior.

The resulting posterior density has the form

$\pi(x|y)\propto\exp\left\\{-\frac{1}{2\sigma^{2}}(y-\eta(x))^{T}(y-\eta(x))\right\\}\times\exp\left\\{\beta\sum_{i\sim
j}u(x_{i}-x_{j})\right\\},\,x\in[2.5,4.5]^{m}.$ (2.2)

The patchiness and speckle allowed by this prior along with the rather global
nature of the likelihood make posterior exploration for this inverse problem
rather challenging, and a good test case for various MCMC schemes that have
been developed over the years. We note that nature of the posterior can be
dramatically altered by changing the prior specification for $x$. This is
discussed later in this section.

This chapter considers a number of MCMC approaches for sampling from this
posterior distribution. We start at the beginning.

####  2.1 Posterior exploration via single-site Metropolis updates

A robust and straightforward method for computing samples from the posterior
$\pi(x|y)$ is the single-site Metropolis scheme, originally carried out in
Metropolis et al., (1953) on the world’s first computer with addressable
memory, the MANIAC. A common formulation of this scheme is summarized in
Algorithm 1 using pseudo code.

Algorithm 1 Single Site Metropolis

1:

2:initialize $x$

3:for $k=1:\text{niter}$ do

4: for $i=1:m$ do

5: $x^{\prime}=x$

6: $x_{i}^{\prime}=x_{i}+z$, where $z\sim N(0,\sigma^{2}_{z})$

7: if $u<\frac{\pi(x^{\prime}|y)}{\pi(x|y)}$, where $u\sim U(0,1)$ then

8: set $x_{i}=x_{i}^{\prime}$

9: end if

10: end for

11:end for

This scheme is engineered to maintain balance with the posterior distribution
$\pi(x|y)$ – so that the nett movement between any two states $x$ and $x^{*}$
is done in proportion to the posterior density at these two points. The width
of the proposal distribution $\sigma_{z}$ should be adjusted so that
inequality in line 7 is satisfied roughly half the time (Gelman et al.,,
1996), but an acceptance rate between 70% and 30% does nearly as well for
single-site updates. The single-site Metropolis in Algorithm 1 scans through
each of the parameter elements in a fixed, deterministic order (for loop,
steps 3–9) to define one _sweep_ over parameter values. One typically records
the current value of $x$ after each sweep; we do so every 10 sweeps.

This single-site scheme was originally intended for distributions with very
local dependencies within the elements of $x$ so that the ratio in line 6
simplifies dramatically. In general, this simplification depends on the full
conditional density of $x_{i}$

$\pi(x_{i}|x_{-i},y),\mbox{ where
}x_{-i}=(x_{1},\ldots,x_{i-1},x_{i+1},\ldots,x_{n})^{T}.$

This density is determined by keeping all of the product terms in $\pi(x|y)$
that contain $x_{i}$, and ignoring the terms that don’t. Hence the ratio in
line 7 can be rewritten as

$\frac{\pi(x^{\prime}|y)}{\pi(x|y)}=\frac{\pi(x_{i}^{\prime}|x_{-i},y)}{\pi(x_{i}|x_{-i},y)}.$

In many cases this ratio becomes trivial to compute. However, in the case of
this particular inverse problem, we must still evaluate the simulator to
compute this ratio. This is exactly what makes the MCMC computation costly for
this problem.

Figure 16.3:  Five realizations from the single-site Metropolis scheme.
Realizations are separated by $1000$ sweeps over the $m$-dimensional image
parameter $x$.

Nonetheless, this straightforward sampling approach does adequately sample the
posterior, given sufficient computational effort; see Frigessi et al., (1993)
for a proof of ergodicity. Figure 16.3 shows realizations produced by the
single-site Metropolis algorithm, separated by $1000$ sweeps. Inspection of
these realizations makes it clear that posterior realizations yield a crisp
distinction between the high and low conductivity regions, as was intended by
the MRF prior for $x$. Around the boundary of the high conductivity region
there is a fair amount of uncertainty as to whether or not a given pixel has
high or low conductivity.

Figure 16.4 shows the resulting posterior mean for $x$ and shows the history
of three pixel values over the course of the single-site updating scheme. The
sampler was run until $40,000\times m$ forward simulations were carried out.
An evenly spaced sample of 4,000 values for three of the $m$ pixels is shown
in the right frame of Figure 16.4. Note that for the middle pixel (blue
circle), the marginal posterior distribution is bimodal – some realizations
have the conductivity value near 3, others near 4. Being able to move between
these modes is crucial for a well mixing chain. Getting this pixel to move
between modes is not simply a matter of getting that one pixel to move by
itself; the movement of that pixel is accomplished by getting the entire image
$x$ to move between local modes of the posterior.

computational effort (simulator evaluations $\times m$)

conductivity

posterior meanMCMC histories for 3 pixels

Figure 16.4:  Posterior mean image for $x$ and MCMC traces of three pixels:
one which is predominantly light (small conductivity); one which is
predominantly dark (high conductivity) and one which is on the edge of the
object. This MCMC run carries out $40000\times m$ forward simulator
evaluations. The value of $x_{i}$ is given every 10th iteration (i.e. every
$10\times m$ single-site updates).

This local multimodality is largely induced by our choice of prior. For
example, if we alter the prior model in (2.1) so that

$u(d)=-d^{2},$ (2.3)

we have a standard Gaussian Markov random field (GMRF) prior for $x$. If, in
addition, the simulator is a linear mapping from inputs $x$ to ouputs
$\eta(x)$, the resulting posterior is necessarily Gaussian, and hence,
unimodal. While this is not true for nonlinear forward models/simulators, the
GMRF prior still has substantial influence on the nature of the posterior;
indeed, recent theoretical results show that the posterior for EIT with GMRF
prior is typically log-concave (Nickl,, 2023), hence unimodal, confirming
these computational observations. Figure 16.5 shows two realizations and the
posterior mean resulting from such a prior with $\beta=2$. Here posterior
realizations are locally more variable – the difference between neighboring
pixels are generally larger. However the global nature of the posterior
realizations are far more controlled than those in Figure 16.3 since the GMRF
prior suppresses local modes that appear under the previous formulation. This
resulting formulation is also far easier to sample, requiring about a 10th of
the effort needed for formulation in (2.2). An alternate, controlling prior
formulation uses a process convolution prior for $x$ is given in the Appendix.
In addition to yielding a more easily sampled posterior, the prior also
represents the image $x$ with far fewer parameters than the $m$ used in the
MRF specifications.

Figure 16.5:  Two realizations and the posterior mean from the single-site
Metropolis scheme run on the posterior resulting from the GMRF prior.
Realizations are separated by $1000$ sweeps over the $m$-dimensional image
$x$.

While these alternative specifications lead to simpler posterior
distributions, they do so at the expense of overly smooth posterior
realizations. Still, such realizations may be useful for exploratory purposes,
and for initializing other samplers; we do not further pursue such
formulations here. Instead, we focus on comparison of various MCMC schemes to
sample the original gray level posterior in (2.2). We use the sample traces
from the three pixels circled in Figure 16.4 to make comparisons between a
variety of samplers which are discussed in the next sections – the movement of
these three pixels is representative of all the image pixels. In particular,
we focus on the frequency of movement between high and low conductivity at
these sites.

###  3 Multivariate updating schemes

Schemes that propose to update more than just a single component of $x$ at a
time have the potential to reduce the computational burden of producing an
MCMC sample from $\pi(x|y)$. The single-site scheme above is also applicable
when the proposal for $x^{\prime}$ changes some or all of the components of
$x$. However, producing a multivariate candidate $x^{\prime}$ that has an
appreciable chance of being accepted (i.e., satisfying the inequality in line
7 of Algorithm 1) while allowing appreciable movement, is very difficult. This
highlights a very appealing aspect of the single-site Metropolis scheme: even
fairly thoughtless 1-d proposals have an appreciable chance of being accepted
while adequately exploring the posterior.

There are clustering MCMC algorithms from statistical physics that allow for
many pixels in $x$ to be updated at once (Edwards and Sokal,, 1988). Such
methods can be adapted to this particular problem as in Higdon, (1998),
however such methods typically show decreased efficiency relative to single-
site updating when the likelihood is strong relative to the prior. This is
certainly the case with our attempts on this application whose results are not
worth discussing here. In the first edition of this handbook Higdon et al.,
(2011) we also explored the differential evolution-MCMC (DE-MCMC) sampler of
ter Braak, (2006). However, we could not find an implementation that
effectively sampled this EIT posterior so we do not further discuss DE-MCMC.
Here we’ll consider some implementations of multivariate random walk
Metropolis updating as competitors to the costly single-site Metropolis
updating for our EIT application.

Later, in Section 4 we will also explore multivariate updates generated semi-
automatically from the single-site Metropolis scheme by the MSDA algorithm,
and highlight the _adaptive_ MSDA algorithm that tunes a modified likelihood
function to improve efficiency by leveraging an otherwise inadequate
approximation to the forward map.

#### Random walk Metropolis

The multivariate random walk Metropolis (RWM) has been the focus of a number
of theoretical investigations (Tierney,, 1994; Gelman et al.,, 1996). But to
date this scheme has not been widely used in applications, and has proven
advantageous only in simple, unimodal settings. The preference for single-
site, or limited multivariate updates in practice may be attributed to how the
full conditionals often simplify computation, or may be due to the difficulty
in tuning highly multivariate proposals. In our EIT application, the
univariate full conditionals do not lead to any computational advantages. If
there is ever an application for which RWM may be preferable, this is it.
Single-site updating is very costly, and may be inefficient relative to
multivariate updating schemes for this multimodal posterior. The adaptive
Metropolis (AM) method of Haario et al., (2006) automatically tunes a
multivariate Gaussian proposal, and allows relatively automatic sampling for
EIT with GMRF prior (Taghizadeh et al.,, 2020). However, we have not found
that direct use of AM successfully samples in a complex inverse problem with
multimodal posterior (see, e.g. Cui et al.,, 2011), such as the present
posterior (2.2) that allows reconstruction of sharp conductivity boundaries in
EIT.

A multivariate Gaussian random walk Metropolis scheme for the $m$-vector $x$
is summarized in Algorithm 2 using pseudo code.

Algorithm 2 Random Walk Metropolis

1:

2:initialize $x$

3:for $k=1:\text{niter}$ do

4: $x^{\prime}=x+z$, where $z\sim N_{m}(0,\Sigma_{z})$

5: if $u<\frac{\pi(x^{\prime}|y)}{\pi(x|y)}$, where $u\sim U(0,1)$ then

6: set $x=x^{\prime}$

7: end if

8:end for

We consider three different proposals for this scheme:

$\displaystyle\Sigma_{z}$ $\displaystyle\propto$
$\displaystyle\Sigma_{1}=I_{m}$ $\displaystyle\Sigma_{z}$
$\displaystyle\propto$
$\displaystyle\Sigma_{2}=\mbox{diag}(s_{1}^{2},\ldots,s_{m}^{2})$
$\displaystyle\Sigma_{z}$ $\displaystyle\propto$
$\displaystyle\Sigma_{3}=S^{2}$

where $s_{i}^{2}$ is the posterior marginal sample variance for the
conductivity $x_{i}$, and $S^{2}$ is the $m\times m$ sample covariance matrix
– both estimated from the previously obtained single-site MCMC run. In each
case we set $\Sigma_{z}=\alpha_{i}\Sigma_{i}$ where the scalar $\alpha_{i}$ is
chosen so that the candidate $x^{\prime}$ is accepted 30% of the time, which
is close to optimal in a Gaussian setting.

RWM $\Sigma_{z}\propto I$RWM $\Sigma_{z}\propto$
diag$(s^{2}_{1},\ldots,s^{2}_{m})$RWM $\Sigma_{z}\propto S^{2}$single-site
Metropoliscomputational effort (simulator evaluations $\times m$)

conductivity

Figure 16.6:  MCMC traces of three pixels circled in Figure 16.4 under three
multivariate random walk Metropolis (RWM) schemes, and single-site Metropolis.
For each run, $40000\times m$ forward simulator evaluations are carried out.
While the RWM scheme with $\Sigma_{z}\propto S^{2}$ results in good movement
for the central pixel, the movement of the top and bottom pixels are clearly
inferior to that of single-site Metropolis.

MCMC traces for these three implementations of RWM are shown in Figure 16.6.
The traces from the single-site Metropolis scheme are also given for
comparison. Interestingly, the behavior of the traces varies with the choice
of $\Sigma_{z}$. The scheme with $\Sigma_{z}\propto S^{2}$ shows the most
movement for the central pixel, which moves between high and low conductivity
over the run. However its performance for the top, low conductivity pixel is
noticeably worse. None of the RWM schemes do as well as single-site Metropolis
when looking at the bottom, high conductivity pixel. These results suggest
that a scheme that utilizes both single-site and RWM updates with
$\Sigma_{z}\propto S^{2}$ might give slightly better posterior exploration
than single-site Metropolis alone.

###  4 Augmenting with fast, approximate simulators

In many applications, a faster, approximate simulator is available for
carrying out the forward model evaluation $\eta_{\rm approx}(x)$. There are a
number of approaches for using fast, approximate simulators to improve the
MCMC exploration of the original (fine-scale) posterior given in (2.2). In
this chapter we consider two related approaches: Metropolis-coupled MCMC
(Geyer,, 1991; Higdon et al.,, 2003; Andersen et al.,, 2003); and delayed
acceptance (DA) schemes that limit the number of calls to the expensive,
“exact” simulator (Liu,, 2001,Christen and Fox,, 2005,Cui et al.,,
2011,Banterle et al.,, 2019). We’ll also highlight the adaptive, multi-level
extension of DA given in Lykkegaard et al., (2023), demonstrating it for the
EIT application of this chapter.

We consider two different approximate forward models: 1) $\eta_{\text{m}}(x)$,
based on an incomplete multigrid solve; and 2) $\eta_{\text{c}}(x)$, based on
a coarsened representation of the conductivity image $x$. The solver used for
$\eta_{\text{m}}(x)$ is derived from the BoxMG algorithm of Dendy, (1982) and
is described in Higdon et al., (2011). Evaluating $\eta_{\text{m}}(x)$ takes
about a third of the time of evaluating the standard forward model $\eta(x)$.
It calculates rather accurate voltages; all well within $\sigma$ for each
component of $\eta(x)$. In contrast, the coarse forward model
$\eta_{\text{c}}(x)$ is well over 100 times faster that $\eta(x)$, but far
less accurate.

Figure 16.7: Fast, approximate coarse forward model $\eta_{\text{c}}(x)$. The
coarse forward model is produced by coarsening the $24\times 24$ conductivity
image (left) and then solving for the measured voltages (right). For the fine-
scale image $x$ (top left), the voltages produced by the original (fine-scale)
forward model $\eta(x)$ are given by the circles; the voltages produced by the
coarse-scale forward model $\eta_{\text{c}}(x)$ are given by the triangles.
The computed voltages differ substantially at the injector electrodes denoted
by the filled in symbols.

Its calculated voltages can be well over $100\sigma$ away from $\eta(x)$ at
the injector electrodes (see Figure 16.7). This error in $\eta_{\text{c}}(x)$
poses problems for the basic implementations of Metropolis coupling and
delayed acceptance described below, motivating an adaptive error
implementation described in Section 4.3.

####  4.1 Metropolis coupling

By augmenting the posterior of interest with auxiliary distributions one can
use Metropolis coupling (Geyer,, 1991), simulated tempering (Marinari and
Parisi,, 1992), or related schemes (Liu and Sabatti,, 1999). Here we use
Metropolis coupling (MC) and define a joint, product posterior for the pair of
$24\times 24$ conductivity images $(x,\tilde{x})$ given by

$\displaystyle\pi(x,\tilde{x}|y)$
$\displaystyle=\pi(x|y)\times\pi_{\text{m}}(\tilde{x}|y)\propto$
$\displaystyle L(y|x)\pi(x)\times L_{\text{m}}(y|\tilde{x})\pi(\tilde{x}),$

where the same priors are specified for $x$ and $\tilde{x}$, and
$L_{\text{m}}(y|\cdot)$ is the original likelihood, but with the multigrid
forward model $\eta_{\text{m}}(\cdot)$

$L_{\text{m}}(y|\tilde{x})\propto\exp\left\\{-\frac{1}{2\sigma^{2}}(y-\eta_{\text{m}}(\tilde{x}))^{T}(y-\eta_{\text{m}}(\tilde{x}))\right\\}.$

The posterior sampling proceeds by carrying out a fixed number of single site
updates for $x$ and $\tilde{x}$, using their respective posteriors, and then
proposing a swap between $x$ and $\tilde{x}$ as shown in Figure 16.8. The swap
from $(x,\tilde{x})$ to $(\tilde{x},x)$ is a deterministic Metropolis
proposal, accepted with probability
$\mbox{min}\\{1,\pi(\tilde{x},x|y)/\pi(x,\tilde{x}|y)\\}$. Since
$\eta_{\text{m}}(x)$ can be evaluated in a third of the time required for
$\eta(x)$, a parallel MC-MCMC implementation should allow three times the
number of single-site Metropolis updates on the approximate posterior relative
to the fine scale posterior before proposing a swap. This will approximately
balance the load between the two processors carrying out the single site
updates. Using the multigrid approximate forward model results in swaps that
are accepted about 80% of the time, owing to its fidelity to the fine scale
forward model. If the coarse forward model $\eta_{\text{c}}(\cdot)$ is used,
swaps are never accepted due to its large error.

Figure 16.8: A sequence from a Metropolis coupled sampler over a joint fine-
scale and approximate posterior. A scan of $3m$ single-site Metropolis updates
are carried out on the approximate posterior while a scan of $m$ updates are
carried out on the fine scale posterior. This is followed by a Metropolis
proposal that swaps the current images of the approximate and fine scale
chains. Marginally, the draws on the fine scale are from the original
posterior distribution.

####  4.2 Delayed acceptance Metropolis

The delayed acceptance approach of Christen and Fox, (2005) uses a fast,
approximate simulator to “filter” a proposal, similarly to the surrogate
transition method of Liu, (2001), for dealing with complex forward models. For
now, we define approximate posterior formulation using the multigrid simulator
$\eta_{\text{m}}(x)$:

$\displaystyle\pi_{\text{m}}(x|y)$ $\displaystyle\propto$ $\displaystyle
L_{\text{m}}(y|x)\times\pi(x).$

A simple Metropolis-based formulation of this scheme is given in Algorithm 3,

Algorithm 3 Delayed Acceptance Metropolis

1:

2:initialize $x$

3:for $k=1:\text{niter}$ do

4: for $i=1:m$ do

5: $x_{i}^{\prime}=x_{i}+z$, where $z\sim N(0,\sigma^{2}_{z})$

6: if $u_{1}<\frac{\pi_{\text{m}}(x^{\prime}|y)}{\pi_{\text{m}}(x|y)}$, where
$u_{1}\sim U(0,1)$ then

7: if
$u_{2}<\frac{\pi(x^{\prime}|y)\pi_{\text{m}}(x|y)}{\pi(x|y)\pi_{\text{m}}(x^{\prime}|y)}$,
where $u_{2}\sim U(0,1)$ then

8: set $x_{i}=x_{i}^{\prime}$

9: end if

10: end if

11: end for

12:end for

where $\pi(x|y)$ and $\pi_{\text{m}}(x|y)$ denote the posteriors using the
exact and approximate simulators respectively. Notice that the exact simulator
need only be run if the filtering condition
($u_{1}<\frac{\pi_{\text{m}}(x^{\prime}|y)}{\pi_{\text{m}}(x|y)}$) involving
the faster, approximate simulator is satisfied. Hence, if the proposal width
is chosen so that the filtering condition is satisfied only a third of the
time, the exact simulator is only run for a third of the MCMC iterations.
Using the multigrid simulator $\eta_{\text{m}}(x)$, the $40000\times m$
iterations required for our original single-site Metropolis scheme take about
66% of the computational effort using this delayed acceptance approach. Using
the coarse simulator $\eta_{\text{c}}(x)$ in this delayed acceptance algorithm
is ineffective. Any proposal accepted in the approximate test (step 5,
algorithm 3) is nearly always rejected in the fine scale test (step 6).

Figure 16.9: Two($+$) updates from a multiple step delayed acceptance (MSDA)
sampler. The current state $x$ from the fine scale is moved up
($\boldsymbol{\uparrow}$) to the approximate level. This serves as the
starting point for one (or more) single site Metropolis updates carried out
according to the approximate posterior $\pi_{\text{m}}(x|y)$ depicted by
$\stackrel{{\scriptstyle\mbox{\tiny MCMC}}}{{\longrightarrow}}$. The resulting
candidate $x^{\prime}$ is accepted on the fine scale ($\searrow$) with
probability
$\mbox{min}\\{1,\pi(x^{\prime}|y)\pi_{\text{m}}(x|y)/\pi(x|y)\pi_{\text{m}}(x^{\prime}|y)$.
If not accepted, the chain remains at $x$ for this DA step.

Multiple step delayed acceptance Algorithm 3 “filters” each proposal (step 4)
before computing an acceptance outcome (step 6). Alternatively, multiple MCMC
steps may be taken on the approximate formulation to produce a candidate
$x^{\prime}$ that differs from the current state $x$ at more than just a
single site $i$.

Algorithm 4 Multiple Step Delayed Acceptance Metropolis

1:

2:initialize $x$

3:for $k=1:\text{niter}$ do

4: set $x^{\prime}=x$

5: for $j=1:n_{\rm step}$ do

6: $x_{i}^{*}=x_{i}+z$, where $z\sim N(0,\sigma^{2}_{z})$ and $i\sim
U\\{1,\ldots,m\\}$

7: if $u_{1}<\frac{\pi_{\text{m}}(x^{*}|y)}{\pi_{\text{m}}(x^{\prime}|y)}$,
where $u_{1}\sim U(0,1)$ then

8: $x_{i}^{\prime}=x^{*}_{i}$

9: end if

10: end for

11: if
$u_{2}<\frac{\pi(x^{\prime}|y)\pi_{\text{m}}(x|y)}{\pi(x|y)\pi_{\text{m}}(x^{\prime}|y)}$,
where $u_{2}\sim U(0,1)$ then

12: set $x=x^{\prime}$

13: end if

14:end for

This multiple step delayed acceptance (MSDA) is detailed in Algorithm 4. This
algorithm produces updates that are in detailed balance with the fine scale
posterior $\pi(x|y)$ (Liu,, 2001; Lykkegaard et al.,, 2023). Figure 16.9
depicts this algorithm; now the horizontal arrows
($\stackrel{{\scriptstyle\mbox{\tiny MCMC}}}{{\longrightarrow}}$) denote
$n_{\rm step}$ single site MCMC steps using the approximate posterior
$\pi_{\text{m}}(x|y)$ as the target distribution. The collection of multiple
steps needs to define a Markov kernel that is in detailed balance with the
approximate posterior $\pi_{\text{m}}(x|y)$ for the second acceptance
probability (step 10 in Alg. 4) to be correct and give an algorithm that
produces draws from $\pi(x|y)$. Randomly choosing the site to update (step 5,
Alg. 4) achieves this, whereas the deterministic scan in Algorithm 1 does not.

This MSDA algorithm is very similar to one half of Metropolis coupling, by
following a single flow of arrows in Figure 16.8 – the acceptance of the swap
proposal in MC-MCMC is identical to acceptance of the eventual proposal
$x^{\prime}$ produced by MSDA. However, there are major differences. MC-MCMC
requires stationarity on the joint product distribution
$\pi(x|y)\times\pi_{\text{m}}(\tilde{x}|y)$. MSDA aims for stationarity only
on the fine scale posterior $\pi(x|y)$. The moves made according to the
approximate posterior $\pi_{\text{m}}(x|y)$ effectively produce a multivariate
proposal $x^{\prime}$ to be considered at the fine scale. If $n_{\rm step}$ is
small, the proposal $x^{\prime}$ is “close” to the fine scale posterior; if
$n_{\rm step}$ very large, the proposal $x^{\prime}$ is closer to the
approximate posterior. Hence, $n_{\rm step}$ could serve as an MCMC tuning
parameter, much like the proposal width used for the single site Metropolis
steps. For the example in Section 4.3, we take $n_{\rm step}=100$.

Finally, we note that Christen and Fox, (2005) and Fox et al., (2020) give a
more general formulation for the delayed acceptance sampler for which the
approximate simulator can depend on the current state $x$ of the chain; the
computed example in Christen and Fox, (2005) used a state-dependent
linearization of the forward map. While a bit more demanding computationally,
the more general algorithm can gain further efficiency by making use of local
approximations and error models. Additionally, we remark that while we used a
coarsened grid to construct the approximate simulator, any computationally
cheaper approximation can be employed. Depending on the broader context of the
problem, options include reduced-precision computing (Higham and Mary,, 2022),
early stopping of iterative solvers (Wikle et al.,, 2001), early termination
of time-dependent problems (Lykkegaard et al.,, 2023), dataset subsampling
(Quiroz et al.,, 2019), reduced grid resolution (Lykkegaard et al.,, 2023),
surrogate models (Laloy et al.,, 2013; Seelinger et al.,, 2024), and
partially-computed likelihoods (Banterle et al.,, 2019).

####  4.3 Adaptive, multiple step delayed acceptance

So far, the coarse $8\times 8$ solver-based forward model has been ineffective
in these MC-MCMC or MSDA implementations, in spite of its over 100-fold speed
up in evaluating $\eta_{\text{c}}(x)$. This is because the coarse scale
simulator $\eta_{\text{c}}(x)$ differs too much from the fine scale simulator
$\eta(x)$; simply swapping out the fine with the coarse simulator in
$L_{\text{c}}(y|x)$ results in a posterior that is too far away from the
target $\pi(x|y)$. Figure 16.10 shows the difference
$\eta(x)-\eta_{\text{c}}(x)$ for the $24\times 24$ conductivity field shown in
the left hand frame of Figure 16.7. The difference is largest at the injector
electrodes.

Figure 16.10: The difference between the fine scale simulator $\eta(x)$ and
the coarse simulator $\eta_{\text{c}}(x)$. Here $\eta(x)-\eta_{\text{c}}(x)$
is shown ($\circ$ plotting symbols) for the image $x$ given in the left hand
frame of Figure 16.1. The differences are consistently large and negative at
the injector electrodes marked by the solid circles. The adaptively estimated
bias term is also shown (- plotting symbols).

Accounting for bias and error due to using an approximate model in place of
the accurate, fine scale model has been a focus in Bayesian inverse problems
(Kaipio and Somersalo,, 2007) and in computer model emulation (Kennedy and
O’Hagan,, 2000). Here we follow Cui et al., (2011) and Fox et al., (2020) and
use adaptive MCMC (Roberts and Rosenthal,, 2009, 2007) to adapt a “coarse”
likelihood $L_{\text{c}}(y|x)$ so that multiple step proposals $x^{\prime}$
are more readily accepted in DA. At step $k$, $L_{\text{c}}(y|x)$ is modified
with an adaptive bias $16^{2}$-vector $b_{k}$ and a $16^{2}\times 16^{2}$
covariance matrix $\Sigma_{bk}$

$L_{\text{c}k}(y|x)\propto\exp\left\\{-\mbox{\small$\frac{1}{2}$}(y-\eta_{\text{c}}(x)-b_{k})^{T}(\sigma^{2}I+\Sigma_{bk})^{-1}(y-\eta_{\text{c}}(x)-b_{k})\right\\}.$

This gives an adaptive version of the coarse scale posterior
$\pi_{\text{c}k}(x|y)=L_{\text{c}k}(y|x)\pi(x)$. This adaptive, multiple step
delayed acceptance algorithm is described in Algorithm 5.

Algorithm 5 Adaptive, Multiple Step Delayed Acceptance Metropolis

1:

2:initialize $x$, $b_{0}=\eta(x)-\eta_{\text{c}}(x)$ and
$\Sigma_{b0}=\mathbf{0}$

3:for $k=1:\text{niter}$ do

4: set $x^{\prime}=x$

5: for $j=1:n_{\rm step}$ do

6: $x_{i}^{*}=x_{i}+z$, where $z\sim N(0,\sigma^{2}_{z})$ and $i\sim
U\\{1,\ldots,m\\}$

7: if
$u_{1}<\frac{\pi_{\text{c}k-1}(x^{*}|y)}{\pi_{\text{c}k-1}(x^{\prime}|y)}$,
where $u_{1}\sim U(0,1)$ then

8: $x_{i}^{\prime}=x^{*}_{i}$

9: end if

10: end for

11: if
$u_{2}<\frac{\pi(x^{\prime}|y)\pi_{\text{c}k-1}(x|y)}{\pi(x|y)\pi_{\text{c}k-1}(x^{\prime}|y)}$,
where $u_{2}\sim U(0,1)$ then

12: set $x=x^{\prime}$

13: end if

14: $b_{k}=\frac{1}{k}\left[(k-1)b_{k-1}+\eta(x)-\eta_{\text{c}}(x)\right]$

15:
$\Sigma_{bk}=\frac{1}{k}\left[(k-1)\Sigma_{bk-1}+(\eta(x)-\eta_{\text{c}}(x)-b_{k})(\eta(x)-\eta_{\text{c}}(x)-b_{k})^{T}\right]$

16:end for

Our implementation takes $n_{\rm step}=100$, $\sigma_{z}=0.3$, and initializes
$b_{0}=\eta(x)-\eta_{\text{c}}(x)$ and $\Sigma_{b0}=\bf{0}$. Figure 16.11
shows the MCMC histories for the same three conductivity pixels shown in
previous figures. The computational effort is equivalent to that of Figure
16.4 – the equivalent of 40,$000\times m$ forward model runs on the fine
scale. After many iterations, the adaptive components $b_{k}$ and $\Sigma_{k}$
stabilize; the resulting values for $b_{\infty}$ are given by the dash (-)
symbols in Figure 16.10.

We employed the free and open-source Python software library tinyDA
(Lykkegaard,, 2023) for the numerical experiments involving Delayed Acceptance
MCMC outlined here. In addition to standard Metropolis-Hastings, DA and MSDA
MCMC, tinyDA provides routines for Multilevel Delayed Acceptance (MLDA), a
multilevel extension to MSDA, which enables MCMC sampling using a multilevel
hierarchy of increasingly coarser models (Lykkegaard et al.,, 2023). When
carefully implemented, such a multilevel model hierarchy can provide even
higher gains to the MCMC sampling efficiency than the two-level sampler
employed in this work.

adaptive, multiple step delayed acceptancesingle-site Metropoliscomputational
effort (simulator evaluations $\times m$)

conductivity

Figure 16.11:  MCMC traces of three pixels circled in Figure 16.4 using the
adaptive, multiple step delayed acceptance (AMSDA) Metropolis scheme. For
comparison, the traces under the single site Metropolis scheme are shown. The
$x$-axis marks computational effort, which is the same for the two MCMC
efforts. The AMSDA sampler shows substantial improvement, making effective use
of the very fast, but highly biased, coarse scale model.

###  5 Discussion

For the EIT example, single-site Metropolis requires about 20 million
simulator evaluations to effectively sample this posterior distribution.
Multivariate updating schemes such as random walk Metropolis or DE-MCMC – as
we implemented them here, or in the previous version of this chapter – don’t
offer any real relief. Tempering and Metropolis-coupling schemes can help if
the approximate forward model is very close to the fine scale forward model.
Here the multigrid-based approximation is quite accurate, but is only three
times faster. Hence the speed-up with such a model will be of a similar
magnitude. The AMSDA approach using the much faster (over $100\times$) coarse
forward model seems to be an exception. The combination of the adaptive error
term and the multiple steps before testing the candidate with the expensive,
fine-scale forward model results in an MCMC scheme that is over an order of
magnitude more efficient. The success of this approach with other challenging
inverse problems (Fox et al.,, 2020; Lykkegaard et al.,, 2023) suggests this
isn’t just due to the special features of this particular EIT example. We also
note that for ease of presentation, we have described rather basic
implementations of delayed acceptance schemes. We point the reader to Fox et
al., (2020), Banterle et al., (2019) and Lykkegaard et al., (2023) for richer
implementations.

One challenging feature of this application is the multimodal nature of the
posterior which is largely induced by our choice of prior. By specifying a
more regularizing prior, such as the GMRF 2.3) or the process convolution
(6.1), the resulting posterior will more likely be unimodal, so that standard
MCMC schemes will be more efficient. Of course, the sacrifice is that one is
now less able to recover small-scale structure that may be present in the
inverse problem.

We considered two different fast, approximate simulators in this chapter.
There is a rather vast literature developing fast approximations to
computationally intensive models. This includes reduced order models (Benner
et al.,, 2017), cross approximations (Dolgov and Scheichl,, 2019), and
response surface/regression approaches such as polynomial chaos (Xiu and
Karniadakis,, 2002), Gaussian processes (Kennedy and O’Hagan,, 2001), or
additive regression trees (Pratola and Higdon,, 2016). The response surface
approaches are most useful when the effective dimension of the parameter
vector is fairly small ($<20$) and the output response to parameter changes is
smooth and predictable. The response surface is typically estimated from a
collection of model runs carried out at different input parameter settings. By
contrast, a well designed reduced order model can enable MCMC to handle some
rather complicated high-dimensional, non-linear problems – see Keating et al.,
(2010), for example. A comparison of reduced model and response surface
approaches can be found in Frangos et al., (2010).

Finally we note that the traditional way to speed up the computation required
to solve an inverse problem is to speed up the simulator $\eta(x)$. A
substantial amount of progress has been made to create simulators that run on
highly distributed computing machines. While MCMC does not lend itself to easy
parallelization, we are seeing advances in the use of parallelized MCMC for
inverse problems (Craiu et al.,, 2009; Vrugt,, 2016; Nishihara et al.,, 2014;
Conrad et al.,, 2018; Brockwell,, 2006). The integration of modern computing
architecture with MCMC methods will certainly extend the reach of MCMC based
solutions to inverse problems.

### Acknowledgments

DH was supported, in part, by the U.S. Department of Energy, Office of
Science, Office of Advanced Scientific Computing Research and Office of High
Energy Physics, Scientific Discovery through Advanced Computing (SciDAC)
program. CF is grateful to digiLab for support and a pleasant working
environment.

## References

  * Andersen et al., (2003) Andersen, K., Brooks, S., and Hansen, M. (2003).  Bayesian inversion of geoelectrical resistivity data.  Journal of the Royal Statistical Society: Series B, 65(3):619–642.
  * Banterle et al., (2019) Banterle, M., Grazian, C., Lee, A., and Robert, C. P. (2019).  Accelerating Metropolis-Hastings algorithms by delayed acceptance.  Foundations of Data Science, 1(2):103–128.
  * Benner et al., (2017) Benner, P., Cohen, A., Ohlberger, M., and Willcox, K. (2017).  Model reduction and approximation: theory and algorithms.  Computational Science and Engineering. SIAM.
  * Brockwell, (2006) Brockwell, A. E. (2006).  Parallel Markov chain Monte Carlo simulation by pre-fetching.  Journal of Computational and Graphical Statistics, 15(1):246–261.
  * Christen and Fox, (2005) Christen, J. and Fox, C. (2005).  Markov chain Monte Carlo using an approximation.  Journal of Computational & Graphical Statistics, 14(4):795–810.
  * Cleveland, (1979) Cleveland, W. S. (1979).  Robust locally weighted regression and smoothing scatterplots.  Journal of the American Statistical Association, 74:829–836.
  * Conrad et al., (2018) Conrad, P. R., Davis, A. D., Marzouk, Y. M., Pillai, N. S., and Smith, A. (2018).  Parallel local approximation MCMC for expensive models.  SIAM/ASA Journal on Uncertainty Quantification, 6(1):339–373.
  * Craiu et al., (2009) Craiu, R. V., Rosenthal, J., and Yang, C. (2009).  Learn from thy neighbor: Parallel-chain and regional adaptive MCMC.  Journal of the American Statistical Association, 104(488):1454–1466.
  * Cui et al., (2011) Cui, T., Fox, C., and O’Sullivan, M. J. (2011).  Bayesian calibration of a large-scale geothermal reservoir model by a new adaptive delayed acceptance Metropolis Hastings algorithm.  Water Resources Research, 47(10).
  * Dendy, (1982) Dendy, J. E. (1982).  Black box multigrid.  J. Comput. Phys., 48:366–386.
  * Dolgov and Scheichl, (2019) Dolgov, S. and Scheichl, R. (2019).  A hybrid alternating least squares–TT cross algorithm for parametric PDEs.  SIAM/ASA Journal on Uncertainty Quantification, 7(1):260–291.
  * Edwards and Sokal, (1988) Edwards, R. G. and Sokal, A. D. (1988).  Generalization of the Fortuin-Kasteleyn-Swendsen-Wang representation and Monte Carlo algorithm.  Physical Review Letters, 38:2009–2012.
  * Fox et al., (2020) Fox, C., Cui, T., and Neumayer, M. (2020).  Randomized reduced forward models for efficient Metropolis–Hastings MCMC, with application to subsurface fluid flow and capacitance tomography.  GEM-International Journal on Geomathematics, 11:1–38.
  * Fox and Nicholls, (1997) Fox, C. and Nicholls, G. (1997).  Sampling conductivity images via MCMC.  The art and science of Bayesian image analysis. Proceedings of the Leeds annual statistics research workshop, pages 91–100.
  * Frangos et al., (2010) Frangos, M., Marzouk, Y., Willcox, K., and van Bloemen Waanders, B. (2010).  Surrogate and reduced-order modeling: a comparison of approaches for large-scale statistical inverse problems.  Large-Scale Inverse Problems and Quantification of Uncertainty, pages 123–149.
  * Frigessi et al., (1993) Frigessi, A., Stefano, P. D., Hwang, C.-R., and Sheu, S.-J. (1993).  Convergence rates of the Gibbs sampler, the Metropolis algorithm and other single-site updating dynamics.  Journal of the Royal Statistical Society. Series B (Methodological), 55(1):205–219.
  * Gelman et al., (1996) Gelman, A., Roberts, G., and Gilks, W. (1996).  Efficient Metropolis jumping rules.  Bayesian Statistics 5: Proceedings of the Fifth Valencia International Meeting, June 5-9, 1994.
  * Geman and McClure, (1987) Geman, S. and McClure, D. (1987).  Statistical methods for tomographic image reconstruction.  Bulliten of the International Statistical Institute, 52:no.4 5–21.
  * Geyer, (1991) Geyer, C. J. (1991).  Monte Carlo maximum likelihood for dependent data.  In Keramidas, E., editor, Computer Science and Statistics: Proceedings of the 23rd Symposium Interface, pages 156–163.
  * Haario et al., (2006) Haario, H., Laine, M., Mira, A., and Saksman, E. (2006).  DRAM: efficient adaptive MCMC.  Statistics and computing, 16:339–354.
  * Higdon, (2002) Higdon, D. (2002).  Space and space-time modeling using process convolutions.  In Anderson, C., Barnett, V., Chatwin, P. C., and El-Shaarawi, A. H., editors, Quantitative Methods for Current Environmental Issues, pages 37–56, London. Springer Verlag.
  * Higdon et al., (2002) Higdon, D., Lee, H., and Bi, Z. (2002).  A Bayesian approach to characterizing uncertainty in inverse problems using coarse and fine scale information.  IEEE Transactions in Signal Processing, 50:389–399.
  * Higdon et al., (2011) Higdon, D., Reese, C. S., Moulton, J. D., Vrugt, J. A., and Fox, C. (2011).  Posterior exploration for computationally intensive forward models.  Handbook of Markov Chain Monte Carlo, pages 401–418.
  * Higdon, (1998) Higdon, D. M. (1998).  Auxiliary variable methods for Markov chain Monte Carlo with applications.  Journal of the American Statistical Association, 93:585–595.
  * Higdon et al., (2003) Higdon, D. M., Lee, H., and Holloman, C. (2003).  Markov chain Monte Carlo-based approaches for inference in computationally intensive inverse problems.  In Bernardo, J. M., Bayarri, M. J., Berger, J. O., Dawid, A. P., Heckerman, D., Smith, A. F. M., and West, M., editors, Bayesian Statistics 7. Proceedings of the Seventh Valencia International Meeting, pages 181–197. Oxford University Press.
  * Higham and Mary, (2022) Higham, N. J. and Mary, T. (2022).  Mixed precision algorithms in numerical linear algebra.  Acta Numerica.  Publisher: Cambridge University Press (CUP).
  * Jimenez et al., (2004) Jimenez, R., Verde, L., Peiris, H., and Kosowsky, A. (2004).  Fast cosmological parameter estimation from microwave background temperature and polarization power spectra.  Physical Review D, 70(2):23005.
  * Kaipio et al., (2000) Kaipio, J., Kolehmainen, V., Somersalo, E., and Vauhkonen, M. (2000).  Statistical inversion and Monte Carlo sampling methods in electrical impedance tomography.  INVERSE PROBL, 16(5):1487–1522.
  * Kaipio and Somersalo, (2007) Kaipio, J. and Somersalo, E. (2007).  Statistical inverse problems: discretization, model reduction and inverse crimes.  Journal of computational and applied mathematics, 198(2):493–504.
  * Kaipio and Somersalo, (2004) Kaipio, J. P. and Somersalo, E. (2004).  Statistical and Computational Inverse Problems.  Springer, New York.
  * Keating et al., (2010) Keating, E. H., Doherty, J., Vrugt, J. A., and Kang, Q. (2010).  Optimization and uncertainty assessment of strongly nonlinear groundwater models with high parameter dimensionality.  Water Resources Research, 46(10).
  * Kennedy and O’Hagan, (2000) Kennedy, M. and O’Hagan, A. (2000).  Predicting the output from a complex computer code when fast approximations are available.  Biometrika, 87:1–13.
  * Kennedy and O’Hagan, (2001) Kennedy, M. and O’Hagan, A. (2001).  Bayesian calibration of computer models (with discussion).  Journal of the Royal Statistical Society (Series B), 68:425–464.
  * Laloy et al., (2013) Laloy, E., Rogiers, B., Vrugt, J. A., Mallants, D., and Jacques, D. (2013).  Efficient posterior exploration of a high-dimensional groundwater model from two-stage Markov chain Monte Carlo simulation and polynomial chaos expansion.  Water Resources Research, 49(5):2664–2682.
  * Liu, (2001) Liu, J. (2001).  Monte Carlo Strategies in Scientific Computing.  Springer.
  * Liu and Sabatti, (1999) Liu, J. and Sabatti, C. (1999).  Simulated sintering: Markov chain Monte Carlo with spaces of varying dimensions.  Bayesian Statistics 6: Proceedings of the Sixth Valencia International Meeting, June 6-10, 1998.
  * Lykkegaard, (2023) Lykkegaard, M. B. (2023).  tinyDA.  https://github.com/mikkelbue/tinyDA.
  * Lykkegaard et al., (2023) Lykkegaard, M. B., Dodwell, T. J., Fox, C., Mingas, G., and Scheichl, R. (2023).  Multilevel delayed acceptance MCMC.  SIAM/ASA Journal on Uncertainty Quantification, 11(1):1–30.
  * Marinari and Parisi, (1992) Marinari, E. and Parisi, G. (1992).  Simulated tempering: a new Monte Carlo scheme.  Europhysics Letters, 19:451–458.
  * Metropolis et al., (1953) Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., and Teller, E. (1953).  Equations of state calculations by fast computing machines.  Journal of Chemical Physics, 21:1087–1091.
  * Moulton et al., (2008) Moulton, J. D., Fox, C., and Svyatskiy, D. (2008).  Multilevel approximations in sample-based inversion from the Dirichlet-to-Neumann map.  J. Phys.: Conf. Ser., 124:012035 (10pp).
  * Nickl, (2023) Nickl, R. (2023).  Bayesian Non-linear Statistical Inverse Problems.  EMS Press.
  * Nishihara et al., (2014) Nishihara, R., Murray, I., and Adams, R. P. (2014).  Parallel MCMC with generalized elliptical slice sampling.  The Journal of Machine Learning Research, 15(1):2087–2112.
  * Nissinen et al., (2008) Nissinen, A., Heikkinen, L. M., and Kaipio, J. P. (2008).  The Bayesian approximation error approach for electrical impedance tomography; experimental results.  Meas. Sci. Technol., 19(1):015501 (9pp).
  * Paciorek and Schervish, (2004) Paciorek, C. J. and Schervish, M. J. (2004).  Nonstationary covariance functions for Gaussian process regression.  In Thrun, S., Saul, L., and Schölkopf, B., editors, Advances in Neural Information Processing Systems 16. MIT Press, Cambridge, MA.
  * Pratola and Higdon, (2016) Pratola, M. T. and Higdon, D. M. (2016).  Bayesian additive regression tree calibration of complex high-dimensional computer models.  Technometrics, 58(2):166–179.
  * Quiroz et al., (2019) Quiroz, M., Kohn, R., Villani, M., and Tran, M.-N. (2019).  Speeding up MCMC by efficient data subsampling.  Journal of the American Statistical Association, 114(526):831–843.
  * Roberts and Rosenthal, (2007) Roberts, G. O. and Rosenthal, J. S. (2007).  Coupling and ergodicity of adaptive Markov chain Monte Carlo algorithms.  Journal of applied probability, 44(2):458–475.
  * Roberts and Rosenthal, (2009) Roberts, G. O. and Rosenthal, J. S. (2009).  Examples of adaptive MCMC.  Journal of computational and graphical statistics, 18(2):349–367.
  * Rue, (2001) Rue, H. (2001).  Fast sampling of Gaussian Markov random fields.  Journal of the Royal Statistical Society: Series B, 63(2):325–338.
  * Seelinger et al., (2024) Seelinger, L., Reinarz, A., Lykkegaard, M. B., Alghamdi, A. M. A., Aristoff, D., Bangerth, W., Bénézech, J., Diez, M., Frey, K., Jakeman, J. D., Jørgensen, J. S., Kim, K.-T., Martinelli, M., Parno, M., Pellegrini, R., Petra, N., Riis, N. A. B., Rosenfeld, K., Serani, A., Tamellini, L., Villa, U., Dodwell, T. J., and Scheichl, R. (2024).  Democratizing Uncertainty Quantification.  arXiv:2402.13768 [cs, stat].
  * Stenerud et al., (2008) Stenerud, V. R., Kippe, V., Lie, K. A., and Datta-Gupta, A. (2008).  Adaptive multiscale streamline simulation and inversion for high-resolution geomodels.  SPE Journal, 13(01):99–111.
  * Taghizadeh et al., (2020) Taghizadeh, L., Karimi, A., Stadlbauer, B., Weninger, W. J., Kaniusas, E., and Heitzinger, C. (2020).  Bayesian inversion for electrical-impedance tomography in medical imaging using the nonlinear Poisson–Boltzmann equation.  Computer Methods in Applied Mechanics and Engineering, 365:112959.
  * ter Braak, (2006) ter Braak, C. F. J. (2006).  A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces.  Statistics and Computing, 16(3):239–249.
  * Tierney, (1994) Tierney, L. (1994).  Markov chains for exploring posterior distributions (with discussion).  Annals of Statistics, 21:1701–1762.
  * Vrugt, (2016) Vrugt, J. A. (2016).  Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation.  Environmental Modelling & Software, 75:273–316.
  * Weir, (1997) Weir, I. (1997).  Fully Bayesian reconstructions from single photon emission computed tomography.  Journal of the American Statistical Association, 92:49–60.
  * Wikle et al., (2001) Wikle, C. K., Milliff, R. F., Nychka, D., and Berliner, L. M. (2001).  Spatiotemporal Hierarchical Bayesian Modeling: Tropical Ocean Surface Winds.  Journal of the American Statistical Association, 96(454):382–397.
  * Xiu and Karniadakis, (2002) Xiu, D. and Karniadakis, G. E. (2002).  The Wiener–Askey polynomial chaos for stochastic differential equations.  SIAM journal on Scientific Computing, 24(2):619–644.

###  6 Formulation based on a process convolution prior

An alternative to treating each pixel in the image as a parameter to be
estimated is to use a lower dimensional representation for the prior. Here we
describe a process convolution (Higdon,, 2002) prior for the underlying image
$x$.

We define $x(s),\,s\in\Omega$ to be a mean zero Gaussian process. But rather
than specify $x(s)$ through its covariance function, it is determined by a
latent process $u$ and a smoothing kernel $k(s)$. The latent process
$u=(u_{1},\ldots,u_{p})^{T}$ is located at the spatial sites
$\omega_{1},\ldots,\omega_{p}$, also in $\Omega$ (shown in Figure 12). The
$u_{j}$’s are then modeled as independent draws from a $N(0,\sigma^{2}_{u})$
distribution. The resulting continuous Gaussian process model for $x(s)$ is
then

$x(s)=\sum_{j=1}^{p}u_{j}k(s-\omega_{j})$ (6.1)

where $k(\cdot-\omega_{j})$ is a kernel centered at $\omega_{j}$. For the EIT
application, we define the smoothing kernel $k(\cdot)$ to be a radially
symmetric bivariate Gaussian density, with standard deviation
$\sigma_{u}=.11$. Figure 12 shows a prior draw from this model over the
$24\times 24$ pixel sites in $\Omega$.

Figure 12:  Left: A $10\times 10$ lattice of locations
$\omega_{1},\ldots,\omega_{p}$, for the $u_{j}$’s of the process convolution
prior; the $24\times 24$ image pixels are shown for reference. Middle: a
realization from the process convolution prior for $x(s)$. Right: posterior
mean from the single-site Metropolis scheme run on the $u$ vector controlling
that controls the image $x$.

Under this formulation, the image $x$ is controlled by $p=100$ parameters in
$u$. Thus a single-site Metropolis scan of $u$ takes less than 20% of the
computational effort required to update each pixel in $x$. In addition, this
prior enforces very smooth realizations for $x$. This makes the posterior
distribution better behaved, but may make posterior realizations of $x$
unreasonably smooth. The resulting posterior mean for $x$ is shown in Figure
12. For a more detailed look at process convolution models, see (Higdon,,
2002); Paciorek and Schervish, (2004) give non-stationary extensions of these
spatial models.