AN ADAPTIVE INDEPENDENCE SAMPLER MCMC ALGORITHM
FOR INFINITE-DIMENSIONAL BAYESIAN INFERENCES∗

arXiv:1508.03283v2 [math.NA] 10 Apr 2016

ZHE FENG† AND JINGLAI LI‡
Abstract. Many scientific and engineering problems require to perform Bayesian inferences in
function spaces, in which the unknowns are of infinite dimension. In such problems, many standard
Markov Chain Monte Carlo (MCMC) algorithms become arbitrary slow under the mesh refinement,
which is referred to as being dimension dependent. In this work we develop an independence sampler
based MCMC method for the infinite dimensional Bayesian inferences. We represent the proposal
distribution as a mixture of a finite number of specially parametrized Gaussian measures. We show
that under the chosen parametrization, the resulting MCMC algorithm is dimension independent.
We also design an efficient adaptive algorithm to adjust the parameter values of the mixtures from
the previous samples. Finally we provide numerical examples to demonstrate the efficiency and
robustness of the proposed method, even for problems with multimodal posterior distributions.
Key words. adaptive Markov chain Monte Carlo, Bayesian inference, Gaussian mixture, independence sampler, inverse problem

1. Introduction. Nonparametric Bayesian inferences have applications in many
scientific problems, ranging from regression [15] to inverse problems [17, 34]. In those
problems the unknown that we want to infer is of infinite-dimension, for example, a
function of space or time. In many practical problems, the posterior distributions do
not admit a closed form and need to be computed numerically. Specifically one first
represents the unknown function with a finite-dimensional parametrization, for example, by discretizing the function on a pre-determined mesh grid, and then solve the
resulting finite dimensional inference problem with the Markov Chain Monte Carlo
(MCMC) simulations. It has been known that standard MCMC algorithms, such as
the random walk Metropolis-Hastings (RWMH), can become arbitrarily slow as the
discretization mesh of the unknown is refined [31, 33, 6, 26]. That is, the mixing time of
an algorithm can increase to infinity as the dimension of the discretized parameter approaches to infinity, and in this case the algorithm is said to be dimension-dependent.
To this end, a very interesting line of research is to develop dimension-independent
MCMC algorithms by requiring the algorithms to be well-defined in the function
spaces. In particular, a family of dimension-independent MCMC algorithms were
presented in [8] by constructing a preconditioned Crank-Nicolson (pCN) discretization of a stochastic partial differential equation (SPDE) that preserves the reference
measure.
Just like its finite dimensional counterparts, the sampling efficiency of the infinite
dimensional MCMC can be improved by incorporating the data information in the
proposal design. One way of doing so is to guide the proposal with the local derivative information of the likelihood function. Methods in this category include: the
stochastic Newton MCMC [25, 28], the operator-weighted proposal method [20], the
infinite-dimensional Metropolis-adjusted Langevin algorithm (MALA) [7, 5], and the
dimension-independent likelihood-informed (DILI) MCMC [9], just to name a few. An
alternative type of methods to improve the efficiency with the data information is the
∗ The

work was supported by the NSFC, under grant 11301337.
of Mathematics, Zhiyuan College, Shanghai Jiao Tong University, 800 Dongchuan
Rd, Shanghai 200240, China, (sjtufz@sjtu.edu.cn).
‡ Corresponding Author, Institute of Natural Sciences, Department of Mathematics, and the
MOE Key Laboratory of Scientific and Engineering Computing, Shanghai Jiao Tong University, 800
Dongchuan Rd, Shanghai 200240, China, (jinglaili@sjtu.edu.cn).
† Department

1

2

Z. FENG AND J. LI

adaptive MCMC (c.f. [1, 2, 32] and the references therein), which automatically adjusts the proposal as the algorithm proceeds. While the first type of approaches utilize
the gradient or the Hessian of the likelihood function to accelerate the computation,
the adaptive methods do not require such information, which makes it particularly
convenient for problems with black-box models.
In this paper we propose an adaptive MCMC algorithm with independence sampler (IS) [35] for infinite dimensional problems. IS, also known as the independent
Metropolis-Hastings (MH) [16], or the Metropolized independent sampling [23], is an
alternative to the popular RWMH algorithm, which proposes from a stationary distribution, i.e., one that is independent of the present position. The design principle
for the independence sampler method is rather straightforward: loosely speaking, one
should choose the proposal distribution to be as close to the target distribution as
possible. The basic idea here is to represent the proposal distribution with a mixture
of a finite number of parametrized Gaussian measures, and optimize the parameters
as the algorithm proceeds. Our specific parametrization ensures the algorithm to be
well defined in function spaces and therefore dimension independent. As is mentioned
earlier, a major advantage of the proposed method is that it can propose efficiently
without using the derivative information of the likelihood function. Moreover as is
demonstrated by our numerical examples in Section 5, our method performs well for
multimodal posterior distributions which can be challenging for many existing algorithms.
The rest of the paper is organized as the following. In section 2 we introduce
the basic setup of the infinite dimensional Bayesian inference problem. In section 3
we present the Gaussian mixture based independence sampler and show that it is
well-defined in the function space. Section 4 is devoted to a detailed description of
the complete algorithm and and section 5 provides several numerical examples of the
proposed method.
2. Problem setup. We consider a separable Hilbert space X with inner product
h·, ·iX . Our goal is to estimate the unknown u ∈ X from data y ∈ Y where Y is the
data space and y is related to u via the likelihood function
L(u, y) =

1
exp(−Φy (u)),
Z

where Z is a normalization constant. In what follows, without causing any ambiguity,
we shall drop the superscript y in Φy for simplicity. In this work we require the
functional Φ satisfies the Assumptions (6.1) in [8], i.e.,
(a) there exists q > 0, Q > 0 such that, for all u ∈ X,
0 ≤ Φ(u) ≤ Q(1 + kukqX );
(b) for every r > 0 there is Qr > 0 such that, for all u, v ∈ X with
max{kukX , kvkX } < r,
|Φ(u) − Φ(v)| ≤ Qr ku − vkX .
We do not have any restrictions on the space Y .
In the Bayesian inference we assume that the prior µ0 of u, is a (without loss of
generality) zero-mean Gaussian measure defined on X with covariance operator C0 ,
i.e. µ0 = N (0, C0 ). Note that C0 is symmetric positive and of trace class. The range

ADAPTIVE MCMC FOR INFINITE DIMENSIONAL INFERENCE

3

1

of C02 ,
1

E = {u = C02 x | x ∈ X} ⊂ X,
which is a Hilbert space equipped with inner product [10],
−1

−1

h·, ·iE = hC0 2 ·, C0 2 ·iX ,
is called the Cameron-Martin space of measure µ0 . In this setting, the posterior
measure µy of u conditional on data y is provided by the Radon-Nikodym derivative:
dµy
1
(u) = exp(−Φ(u)),
dµ0
Z

(2.1)

which can be interpreted as the Bayes’ rule in the infinite dimensional setting. Our
goal is to draw samples from the posterior µy with MCMC algorithms.
Note that the definition of the maximum a posteriori (MAP) estimator in finite
dimensional spaces does not apply here, as the measures µy and µ0 are not absolutely
continuously with respect to the Lebesgue measure; instead, the MAP estimator in
X is defined as the minimizer of the Onsager-Machlup functional (OMF) [11, 21]:
1
(2.2)
I(u) := Φ(u) + kuk2E ,
2
over the Cameron-Martin space E. In Section 5, we shall use OMF as an indicating
quantity to compare the performance of various MCMC algorithms. Finally we quote
the following lemma ([10], Chapter 1), which will be useful in next section:
Lemma 2.1. There exists a complete orthonormal basis {ek }k∈N on
X and a seP∞
quence of non-negative numbers {αk }k∈N such that C0 ek = αk ek and k=1 αk < ∞,
i.e., {ek }k∈N and {αk }k∈N being the eigenfunctions and eigenvalues of C0 respectively.
Without loss of generality, we assume that the eigenvalues {αk }∞
k=1 are in a descending order.
3. Gaussian mixture based independence sampler. In this section, we
present our Gaussian mixture based independence sampler and show that it is welldefined in the function space.
3.1. Independence sampler MCMC. We start by briefly reviewing the independence sampler MCMC algorithm. Given a proposal distribution µ, we define
measures
ν(du, du0 ) = µ(du0 )µy (du),
ν † (du, du0 ) = µ(du)µy (du0 ),
on the product space X × X. When ν † is absolute continuous with respect to ν, we
can define the acceptance probability [36]


dν †
0
0
A(u, u ) = min 1,
(u, u ) ,
(3.1)
dν
where
dν †
dµy 0 dµ
(u, u0 ) =
(u ) y (u).
dν
dµ
dµ
The IS MCMC in a function space proceeds as follows in each iteration:

(3.2)

4

Z. FENG AND J. LI

1. Draw a sample uproposed from the proposal µ.
2. Let unext = uproposed with probability A(ucurrent , uproposed ) and
unext = ucurrent with probability 1 − A(ucurrent , uproposed ).
We reinstate that the function space IS algorithm is well-defined if and only
if ν † is absolutely continuous with respect to ν, which requires that µ and µy are
equivalent to each other. Since µy and µ0 are equivalent, it suffices to require µ and
µ0 to be equivalent. Interestingly, the pCN scheme with a specific choice of parameter
values yields a dimension-independent IS whose proposal distribution is simply the
prior. Despite its dimension-independence property, simply proposing according to
the prior is inefficient when the data is highly informative, i.e., the posterior being far
apart from the prior. Next we shall introduce a more efficient proposal measure than
the prior that is to be used in IS MCMC algorithms.
3.2. Gaussian mixture proposals. In finite dimensional Bayesian inference
problems, Gaussian mixture (GM) distributions [27] are often used as the IS proposal
distributions for their flexibility and convenience to draw samples from. We now
extend the use of GM to the infinite dimensional setting. Let {µj }Jj=1 be a set of
Gaussian measures on X with µj = N (mj , Cj ) for j = 1...J, and we define the
Gaussian mixture proposal as
µ(dx) =

J
X

wj µj (dx)

(3.3)

j=1

PJ
where {wj }Jj=1 are the mixing weights with j=1 wj = 1. It should be clear that µ
is equivalent to µ0 as long as each µj is equivalent to µ0 , and moreover the RadonNikodym derivative of µ to µ0 is
J

X dµj
dµ
(u) =
wj
(u).
dµ0
dµ0
j=1

(3.4)

Next we discuss our parametrization of each µi . First recall that, according to
Lemma 2.1, {ek }k∈N form a complete basis set of X. Our parametrization of µi
is in the form of:
mj =

∞
X

xj,k αk ek ,

(3.5a)

k=1

Cj−1 = C0−1 + Hj

(3.5b)

where each Hj is defined as
Hj · =

∞
X

hj,k hek , ·iek

(3.5c)

k=1

and xj,k and hj,k are coefficients. The following Theorem provides a sufficient condition for µj = N (mj , Cj ) to be a well defined Gaussian measure on X and equivalent
to µ0 .
Theorem 3.1. If xj , hj ∈ l2 , and hj,k > − α1k for all k ∈ N, µj = N (mj , Cj ) is
a Gaussian measure on X that is equivalent to µ0 .

ADAPTIVE MCMC FOR INFINITE DIMENSIONAL INFERENCE

5

Proof. We let {βj,k }k∈N be the eigenvalues of Cj , i.e, Cj ek = βj,k ek for all k ∈ N.
And it is easy to see that,
βj,k = (αk−1 + hj,k )−1 =
As xj , hj ∈ l2 ,

1
1+αk hj,k

is bounded and thus

αk
.
1 + αk hj,k

P∞

k=1

(3.6)

βj,k < ∞. It follows that Cj ∈

L+
1 (X)

and µj = N (mj , Cj ) defines a Gaussian measure on X.
We now show that µj is equivalent to µ0 . First we introduce µ0j = N (0, Cj ).
Using Eq. (3.6) and hj,k > − α1k for all k ∈ N, we can get
∞
X
(βj,k − αk )2
k=1

(βj,k + αk )2

=

∞
X
k=1

∞
X
αk2 h2j,k
≤
αk2 h2j,k < ∞,
(2 + αk hj,k )2
k=1

as limk→∞ αk = 0 and hj ∈ l2 . By the Feldman-Hajek theorem [10], we have µ0j
1

1

is equivalent to µ0 . Now recall that mj ∈ E = C02 (X) = Cj2 (X), and so we have
µ0j = N (0, Cj ) and µj = N (mj , Cj ) are equivalent, which completes the proof.
Let us assume for now that the conditions in Theorem 3.1 is satisfied and we shall
verify this assumption later. It is easy to show that
1
dµj
|C0 |1/2
1 −1/2
mj k2X + hu, Cj−1 mj iX − hu, Hj uiX )
(u) =
exp(− kCj
1/2
dµ0
2
2
|Cj |


∞ r
∞  2
Y
X
αk
αk 2
1
2αk
=
exp −
xj,k + hj,k u2k −
xj,k uk , (3.7)
βj,k
2
βj,k
βj,k
k=1

k=1

where uk = hu, ek i is the projection of u onto ek . Note that the density dµj /dµ0
actually depends on mj and hj , and thus for convenience’s sake, we define a function
f (·, ·, ·) such that
f (u, xj , hj ) =

dµj
(u),
dµ0

and we then can derive from Eq (3.4) that
J

X
dµy
1
(u) = exp(−Φ(u))/(
wj f (u, xj , hj )),
dµ
Z
j=1
and the density dµ/dµy can be computed accordingly.
3.3. Minimizing the Kullback-Leibler divergence. Now recall that for the
algorithm to be efficient we need the proposal µ to be close to µy and a natural choice
is to determine µ by minimizing the Kullback-Leibler divergence (KLD) between µy
and µ:
Z
dµy
DKL (µy ||µ) = log
(u)µy (du),
(3.8)
dµ
where µ is parametrized with Eq. (3.5). Note that xj and hj are set to be of infinite
dimensions in the formulation above. In numerical simulations, however, xj and hj
must be truncated at some finite number K. Such a truncation is also reasonable

6

Z. FENG AND J. LI

from a practical point of view. In fact, one often can realistically assume that the
data is only informative on a finite number of directions [9, 8] in X, and under this
assumption, we only need to keep a finite number of components of each xj and hj .
We emphasize that K which represents the number of dimensions that are informed
by the data (i.e., the so-called intrinsic dimensionality), should not be confused with
the discretization dimensionality of the problem, i.e., the number of mesh points used
to represent the unknown. Determining the value of K is an important task for our
algorithm and here we choose K with a heuristic approach:
K = min{k ∈ N |

αk
< },
α1

where  is a prescribed threshold. In what follows, we shall adopt this finite, Kdimensional formulation, and thus we have the following optimization problem:
min

{xj , hj ∈RK , wj ∈[0,1]}J
j=1

DKL (µy ||µ),

(3.9)

PJ
subject to j=1 wj = 1. By some elementary calculations, we can show that Eq. (3.9)
is equivalent to
Z
min

{xj , hj ∈RK , wj ∈[0,1]}J
j=1

−

J
X
wj f (u, xj , hj )]µy (du),
log[

(3.10)

j=1

PJ
subject to j=1 wj = 1. We now show that the proposal µ constructed this way is
well-defined in function space, and to this end we have the following corollary.
Corollary 3.2. If {xj , hj , wj }Jj=1 is a solution of Eq. (3.10) , the resulting µ is
equivalent to µ0 .
Proof. It is obvious that if {xj , hj , wj }Jj=1 is a solution of Eq. (3.10), xj , hj ∈ l2 .
Taking partial derivative of the objective function in Eq. (3.10) with respect to hj,k
and setting it to be zero yields the following equation:
Z

αk
dµ
=
PJ
1
+
αk hj,k
l=1 wl f (u, xl , hl )
wj f (u, xj , hj )

y

Z

wj f (u, xj , hj )(αk xj,k − uk )2
)dµy .
PJ
l=1 wl f (u, xl , hl )

As the following two integrals are obviously positive:
Z

wj f (u, xj , hj )
PJ

l=1

wl f (u, xl , hl )

dµy > 0,

Z
and

wj f (u, xj , hj )(αk xj,k − uk )2
)dµy > 0,
PJ
l=1 wl f (u, xl , hl )

we have 1 + αk hj,k > 0. Thus all the conditions of Theorem 3.1 are satisfied and the
corollary follows immediately.
Finally we note that, in the special case where J=1, namely the proposal being simply a Gaussian distribution, our parametrization is similar to the finite rank
representation used in [29, 30]. In fact, the aforementioned works also proposed to approximate the posterior with a Gaussian distribution by minimizing the KLD between
the two distributions. The major difference is the KLD (recall that it is asymmetric)
formulation: the authors of [29, 30] compute the divergence from the Gaussian approximation to the true posterior, while here we compute the divergence the other way
around. An advantage of the present formulation is that the solution to Eq. (3.10)

ADAPTIVE MCMC FOR INFINITE DIMENSIONAL INFERENCE

7

can be explicitly obtained:
Z
1
uk dµy ,
αk
1
1
−
hk = R
,
αk
(αk xk − uk )2 dµy
xk =

(3.11a)
(3.11b)

for k = 1...K, while in their formulation the resulting optimization problem has to be
solved with a stochastic optimization algorithm. The explicit solutions (3.11) are of
essential importance in our adaptive algorithm.
4. The adaptive algorithm. In this section we discuss the algorithm to implement the IS method proposed in Section 3, starting with an introduction to the
adaptive MCMC.
4.1. Adaptive MCMC. The basic idea of the adaptive MCMC is to repeatedly
adjust the proposal parameters using the information in the previous samples. Here
we are focused on the adaptive algorithms with IS [16, 13, 19, 12], while noting
that other types of adaptive algorithms include the adaptive MH [14], the adaptive
MALA [3, 24], and the adaptive Metropolis-within-Gibbs (MwG) [32]. Specifically
our adaptive algorithm has the following three key ingredients. First, to enforce
the asymptotic ergodicity, we terminate the adaptation in a finite number of steps.
Secondly we use a tempered pre-run to obtain the initial parameter values for the
iteration. Simply speaking the technique of tempering is to construct a sequence
of intermediate distributions that converge to the true posterior µy and use these
intermediate distributions to guide the MCMC samples to the true posterior. This
strategy is particularly useful for multimodal posterior distributions. Without loss of
generality, we assume that the tempering distributions are augmented by a tempering
parameter λ:
dµy,λ
∝ exp(−λΦ(u)),
dµ0
and clearly µy,λ = µy when λ = 1 and the tempering distribution is “wider” than the
true posterior for 0 ≤ λ < 1. In practice we can choose a finite number of tempering
Itemp
parameters {λi }i=1
where 0 ≤ λ1 < λ2 < ... < λItemp = 1. We also note that for
problems where the posterior is not too far apart from the prior, tempering may not
be necessary. Finally we estimate and update the proposal parameters after every
fixed number of iterations. The adaptive scheme is summarized as the following:
• Initialization: the total number of iterations Itol , the number of adapted
iterations Iadp , the number of pre-run (tempering) iterations Itemp , a set of
Itemp
tempering parameters {λi }i=1
, the number of samples used in each tempered
iteration Ntemp , the number of samples in each iteration NS .
• Pre-run (optional): let µ0(0) = µ0 ; for i = 1 : Itemp perform:
1. Run MCMC with proposal µ0(i−1) to draw a set of Ntemp samples from
µy,λi , denoted by Si0 .
2. Update the parameter values with samples Si0 obtaining proposal µ0(i) ;
• Iteration: let S = ∅ and µ(0) = µ0(Itemp ) ; for i= 1 to Itol perform:
1. Run MCMC with proposal µ(i−1) to draw a set of NS samples from µy ,
denoted by Si . Let S = S ∪ Si .

8

Z. FENG AND J. LI

2. If i < Iadp , update the parameter values with samples S obtaining proposal µ(i) ; otherwise, let µ(i) = µ(i−1) .
The adaptive algorithm presented above is rather simple; we note, however, that
our method is rather flexible and one can pair it with any desired adaptive IS algorithm. A key step in the adaptive algorithm is to estimate the parameters from the
samples, which is done by solving the sample average estimator of the optimization
problem (3.10):
N
X

max

{xj ,hj ,wj }J
j=1

subject to

PJ

j=1

n=1

J
X
log[
wj f (un , xj , hj )],

(4.1)

j=1

wj = 1. Next we discuss two methods to solve Eq. (4.1).

4.2. Expectation Maximization algorithm. The expectation maximization
(EM) is one of the most popular methods to determine the parameters in mixture
models [27]. Simply put, the EM algorithm iteratively updates the parameter values
in a way that the function value is always increased until convergence is achieved.
Each iteration consists of an Expectation-step and a Maximization-step. It should be
noted that, the EM algorithm, is not guaranteed to converge to the optimal solutions
in general [37]. The theory and implementation details of the EM algorithm and its
application to mixture models can be found in the aforementioned references , and
we shall not repeat them here. When applied to our problem, the update formula
in each iteration can be explicitly obtained. In the Expectation-Step, the membership probability qjn , namely the probability that a sample un is in the mixture j, is
computed,
wj f (un , hj , mj )
,
qjn = PJ
n
j=1 wj f (u , hj , mj )

(4.2)

for each j = 1...J and n = 1...N ; in the Maximization-Step, the parameter values are
updated using the following equations:
N
1 X i
q ,
wj =
N i=1 j

xj,k =

hj,k =

(4.3a)

N
X
1
q n un
N αk wj n=1 j k
N
X
n=1

qjn (

N
X

n=1

qjn (αk xj,k − unk )2 )−1 −

(4.3b)
1
,
αk

(4.3c)

where unk = hun , ek i. The EM algorithm is arguably the most common method to estimate the parameters of mixtures. However, our numerical tests indicate that in some
practical problems the EM algorithm is not sufficiently reliable especially when the
sample set only contains a small number of accepted draws. Moreover, our algorithm
frequently updates the proposal parameters, which makes the computationally intensive EM algorithms less attractive from an efficiency perspective. For these reasons,
we propose an alternative method to EM, which estimates the mixture parameters
using clustering.

ADAPTIVE MCMC FOR INFINITE DIMENSIONAL INFERENCE

9

4.3. Estimating parameters with clustering. Our estimation method with
clustering is largely based on the finite dimensional method developed in [13]. The
idea is rather simple: one first partitions the samples into several clusters and then fit
each cluster with a Gaussian distribution. A difficulty here is that our MCMC samples
are of infinite dimension, which makes clustering challenging. To solve the problem,
we first project the samples onto the K eigenfunctions of the covariance operator and
n
n
then cluster the resulting K dimensional data {(un1 , ..., unK )}N
n=1 and uk = hu , ek i.
Specifically we use the k-means algorithm to cluster the data, and the number of
clusters J is determined with the Bayesian information criteria (BIC) method [27]. In
fact we have found in our numerical tests that the algorithm is rather robust against
the number of clusters. We then use the Gaussian distribution parametrized in the
form of Eq (3.5) to fit each cluster, and thanks to Eq. (3.11), the parameters values
can be estimated explicitly as,
xj,k =

X
1
unk ,
Nj αk n

(4.4a)

u ∈Θj

hj,k =

1
Nj

1
1
−
,
n
2
2
α
k
un ∈Θj (uk ) − mj,k

P

(4.4b)

where Θj is the j-th cluster of samples, Nj is the sample size of Θj , for j = 1...J and
k = 1...K. The mixture weights are simply determined by the fraction of samples
in each cluster. We note that the clustering based method does not generally yield
a solution to Eq. (4.1) and thus we regard it as an approximate method to estimate
the parameters. We conclude the section with a pseudo code (Algorithm 1) of our
algorithm, and interested readers can use it as a basis for their own implementation.
5. Numerical examples.
5.1. An ordinary differential equation example. Our first example is a
simple inverse problem where the forward model is governed by an ordinary differential
equation (ODE):
dx(t)
= −u(t)x(t)
dt

(5.1)

with a prescribed initial condition. We assume that the solution x(t) is observed at
several times in the interval [0, T ] and we want to infer the unknown coefficient u(t)
for t ∈ [0, T ].
In our numerical experiments, we let the initial condition be η(0) = 1 and T = 1.
Now suppose that the solution is measured every T /20 time unit from 0 to T and
the error in each measurement is assumed to be an independent zero-mean Gaussian
random variable with variance 0.052 . In the computation, 100 equally spaced grid
points are used to represent the unknown. Moreover, we assume that the state space
for u is X = L2 ([0, T ]) and the prior is a zero-mean Gaussian measure in X with an
exponential covariance function:
C(t, t0 ) = exp(−

|t − t0 |
).
2

(5.2)

The true coefficient u(t) is a realization from the prior (shown in Fig. 5.1) and the
data is simulated accordingly.

10

Z. FENG AND J. LI

I

temp
, Ntemp , Ntol , Nmax , Nadp .
input : Itemp , {λi }i=1
tol
output: Ntol samples drawn from µy : {un }N
n=1 .
µ ← µ0 ;
for i ← 1 to Itemp do
draw u0 ∼ µ;
for n ← 1 to Ntemp do
draw u0 ∼ µ;
draw a ∼ U [0, 1] and compute

A ← min{1,

dµy,λi 0 dµ
(u ) y,λi (un−1 )};
dµ
dµ

if A > a then un ← u0 ; else un ← un−1 ;
end
cluster {u0 , ..., uNtemp } into J subsets: Θ1 , ..., ΘJ ;
for j ← 1 to J do
Nj ← sample size of Θj , wj ← Nj /Ntemp ;
compute parameters xj and hj using Eqs. (4.4);
compute µj using Eq. (3.7);
end
PJ
µ ← j=1 wj µj ;
end
draw u0 ∼ µ;
for n ← 1 to N do
draw u0 ∼ µ;
draw a ∼ U [0, 1] and compute
A ← min{1,

dµy 0 dµ n−1
(u ) y (u
)};
dµ
dµ

if A > a then un ← u0 ; else un ← un−1 ;
if (n < Nmax )&(n mod Nadp = 0) then
cluster {u0 , ..., un } into J subsets: Θ1 , ..., ΘJ ;
for j ← 1 to J do
Nj ← sample size of Θj , wj ← Nj /n;
compute parameters xj and hj using Eqs. (4.4);
compute µj using Eq. (3.7);
end
PJ
µ ← j=1 wj µj ;
end
end
Algorithm 1: The complete algorithm for the adaptive IS with GM. Itemp is the
Itemp
number of tempered iterations. {λi }i=1
are the tempering parameters. Ntemp is
the number of samples used in each tempered iteration. Ntol is the total number
of samples drawn by the algorithm. Nadp is the number of samples drawn between
two consecutive parameter updates. Nmax is the maximum length of chain before
the adaptation is terminated.

ADAPTIVE MCMC FOR INFINITE DIMENSIONAL INFERENCE

11

1.2
1
0.8

u(t)

0.6
0.4
0.2
Truth
Prior
Mixture
Gaussian
pCN

0
-0.2
-0.4
-0.6

0

0.2

0.4

0.6

0.8

1

t
Fig. 5.1. (for example 1) The posterior mean computed with the four different MCMC schemes.
The truth is also plotted for comparison.

We now draw samples from the posterior of u(t) with four different MCMC
schemes: prior-based IS, adaptive IS with Gaussian approximation, adaptive IS with
Gaussian mixtures, and the random walk pCN (RW-pCN). In each MCMC scheme,
3 × 105 draws are generated. In the prior based IS, one simply proposes according to
the prior distribution, and no adaptation is used. In the adaptive IS with Gaussian
approximation, the proposal is restricted to be a single Gaussian (i.e. J = 1), and in
this case clustering is not needed. In both of the adaptive IS methods, the parameters
are updated after every 1000 draws, and the parameter adaptation is terminated in
the last 105 iterations. We do not use tempering in this example. The RW-pCN
algorithm used in this work p
iterates as follows
1. propose uproposed = 1 − β 2 ucurrent + βw, where w ∼ µ0
2. Let unext = uproposed with probability
a = min{1, exp(Φ(uproposed ) − Φ(ucurrent ))},
and let unext = ucurrent with probability 1 − a.
In this example we use β = 0.1. Note that, in all the numerical examples, we choose
the stepsize β so that the resulting acceptance probability is in the range 20% − 30%,
as is recommended in [33].
In Fig. 5.1, we show the posterior mean computed by the four MCMC schemes,
while the truth is also shown for comparison purpose. One can see that the results
of the four algorithms are nearly identical, suggesting that all the algorithms can
estimate the posterior mean to a similar level of accuracy. We then use the OMF as
an indicative parameter and show the trace plots of it in Fig. 5.2. We see from the
plots that the two adaptive IS algorithms achieve much faster mixing rate than the
other two methods. To further compare the efficiency of the methods, we compute
the autocorrelation functions (ACF) of various quantities with the samples drawn by
the four methods, and plot the ACF results in Fig. 5.3. In particular, we plot the
ACF of the OMF as a function of lag in Fig. 5.3 (left) and show the lag 1 ACF for the

12

Z. FENG AND J. LI

Prior

100
90

90

80

80

70

70

60

60

50

50

40

40

30

0

1

2

3
×10 5

Gaussian

100

30

90

80

80

70

70

60

60

50

50

40

40
0

1

0

1

2

3

×10

5

30

2

3
×10 5

2

3
×10 5

pCN

100

90

30

Mixture

100

0

1

Fig. 5.2. (for example 1) The trace plots of the OMF for the four different MCMC schemes.

unknown u at each grid point in in Fig. 5.3 (right). It can be seen from the figure that,
our adaptive algorithms with single Gaussian proposal and with mixtures both result
in much lower ACF values than the other two methods. When comparing the two
adaptive algorithms, the mixture proposal outperforms the single Gaussian. For the
IS algorithms, the acceptance probability is also a useful performance indicator, where
higher acceptance rates are usually preferred, while it is not the case for random walk
algorithms [33]. In Fig 5.4 (left) we plot the acceptance probability as a function of

13

ADAPTIVE MCMC FOR INFINITE DIMENSIONAL INFERENCE

1

1

0.8
0.6

ACF (lag 1)

0.9

ACF

Prior
Mixture
Gaussian
pCN

0.4

0.7

0.2
0

Prior
Mixture
Gaussian
pCN

0.8

0

50

100
lag

150

0.6

200

0

0.2

0.4

0.6

0.8

1

t

Fig. 5.3. (for example 1) Autocorrelation functions (ACF) for the four different MCMC methods. Left: ACF of the OMF plotted as a function of lags. Right: the lag 1 ACF for u at each grid
point.
4
3 ×10

0.4

2.5

0.3
0.25

2
ESS

Acceptance rate

0.35

0.2
Prior
Mixture
Gaussian
pCN

0.15
0.1

1

0.05
0
0

50

100
150
200
Iterations (× 10 3)

Prior
Mixture
Gaussian
pCN

1.5

250

300

0.5

0

0.4

0.8

1.2

1.6

2.0

t

Fig. 5.4. (for example 1) Left: the acceptance rate of the four MCMC schemes. Right: the
ESS at each grid point.

iterations for all the methods. For the three IS algorithms, one can see that the two
adaptive algorithms have significantly higher acceptance probability than the prior
based method. Meanwhile, the acceptance probability of IS with mixtures is higher
than that of the one with the single Gaussian. The effective sample size (ESS) is
another common measure of the sampling efficiency of MCMC [18]. ESS is computed
by
ESS =

N
,
1 + 2τ

where τ is the integrated autocorrelation time and N is the total sample size, and it
gives an estimate of the number of effectively independent draws in the chain. We
computed the ESS of the unknown u at each grid point and show the results in Fig. 5.4
(right). Once again, the plots indicate that the adaptive algorithms produce much
more effectively independent samples than the prior based IS and the RW-pCN, while
the mixture proposal outperforms the single Gaussian one in most of the dimensions.
In summary, in this simple nonlinear inverse problem, we show that our adaptive
algorithms are significantly more efficient than the prior based IS and the RW-pCN.
Meanwhile, the mixture proposal outperforms the single Gaussian one, indicating that
the more flexible mixture representation does improve the efficiency.

14

Z. FENG AND J. LI

Prior

1.5
1

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1

-1.5

0

0.2

0.4

0.6

0.8

1

Gaussian

1.5

-1.5

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1
0

0.2

0.4

0.6

0

0.2

0.4

0.8

1

-1.5

0.6

0.8

1

0.6

0.8

1

pCN

1.5

1

-1.5

Mixture

1.5

0

0.2

0.4

Fig. 5.5. (for example 2) 100 samples randomly selected from the chain drawn by each method.

5.2. A bimodal likelihood function example. Our second example is an
artificially constructed bimodal problem. Once again we assume the unknown u ∈
X = L2 ([0, 1]) and the prior is a zero mean Gaussian measure with the same covariance
function Eq. (5.2) as the first example. We consider a bimodal likelihood function,
given by,
1
1
exp(−Φ(u)) ∝ exp(− ku − sin(2πt)k22 ) + exp(− ku + sin(2πt)k22 ),
2
2
and it can be verified that the Φ(·) chosen this way satisfies the Assumptions (6.1)
in [8]. It is easy to see that the posterior distribution should have two modes: one is
close to sin(2πt) and the other is close to − sin(2πt).
We draw samples from the posterior of u(t) with the same four MCMC schemes
used in the first example, and in each MCMC scheme, 5 × 105 draws are generated. In
both of the adaptive IS methods, the parameters are updated after every 1000 draws,
and the adaptation is terminated in the last 105 iterations, with no tempering used.
In the RW-pCN, we choose β = 0.5. In all the computations, 100 grid points are used
to represent the unknown function u.
As has been mentioned, the posterior distribution has two modes and we shall
exam if the algorithms can capture both of them. In this respect, we randomly select
100 samples from the chain generated by each algorithm and plot them in Fig. 5.5.
We can see that the results of each algorithm can capture the two models of the
posterior. Next we shall compare the efficiency of the four algorithms. As before,
we first show the trace plots of the OMF for the four algorithms in Fig 5.6 and one

ADAPTIVE MCMC FOR INFINITE DIMENSIONAL INFERENCE

Prior

100

80

60

60

40

40

20

0

1

2

3

4

5
×10 5

Gaussian

100

20

80

60

60

40

40

0

1

2

3

0

1

2

4

5

×10

5

20

3

4

5
×10 5

3

4

5
×10 5

pCN

100

80

20

Mixture

100

80

15

0

1

2

Fig. 5.6. (for example 2) The trace plots of the OMF for the four different MCMC schemes.

can see that the results of the two adaptive methods and pCN all obtain fairly good
mixing results, while the prior based IS seems to have a much slower mixing rate than
the other three. Fig. 5.7 (left) plots the ACF of the OMF as a function of lag and
Fig. 5.7 (right) shows the lag 1 ACF for the unknown at each grid point. Both figures
indicate that the adaptive IS with mixtures has the best performance in terms of
ACF values. Fig. 5.8 (left) plots the acceptance rate against the number of iterations,
which shows that the three IS algorithms perform very differently: the prior based

16

Z. FENG AND J. LI

1
Prior
Mixture
Gaussian
pCN

0.8

1

0.8
ACF (lag 1)

ACF

0.6
0.4
0.2

0.6
Prior
Mixture
Gaussian
pCN

0.4

0
0

50

100
lag

150

0.2
0

200

0.2

0.4

0.6

0.8

1

t

Fig. 5.7. (for example 2) Autocorrelation functions (ACF) for the four different MCMC methods. Left: ACF of the OMF plotted as a function of lags. Right: the lag 1 ACF for u at each grid
point.

4

10

x 10

8

0.6
Prior
Mixture
Gaussian
pCN

0.4

6
ESS

acceptance rate

0.8

4

Prior
Mixture
Gaussian
pCN

0.2
2
0
0

100

200
300
Iterations (× 103)

400

500

0
0

0.2

0.4

0.6

0.8

1

t

Fig. 5.8. (for example 2) Left: the acceptance rate of the four MCMC schemes. Right: the
ESS at each grid point.

Mixture

The mean value of Different proposal

pCN

1.5

3

3

1

2

2

0.5

1

1

0

0

−0.5

−1

−1

−1

−2

−2

Mixture
pCN

0

−1.5
0

0.2

0.4

0.6

0.8

1

−3
0

0.5

1

−3
0

0.5

1

Fig. 5.9. (for example 2) Left: the sample mean of the mixture-based IS method (solid) and
that of the pCN method (dashed). Right: samples drawn by the mixture IS method and by the pCN
method.

ADAPTIVE MCMC FOR INFINITE DIMENSIONAL INFERENCE

17

IS results in an acceptance rate less than 1%, the adaptive IS with one Gaussian
results in a rate up to 17%, and that of the adaptive IS with mixtures rises to around
80% as the iteration proceeds. We compute the ESS of each dimension and show the
results in Fig. 5.8 (right), and we see that the ESS of the adaptive IS with mixtures is
significantly higher than that of the other three methods, indicating that the adaptive
IS with mixtures has a substantial advantage in this multimodal problem.
Finally to understand the limitation of the proposed method, we test it on another
bimodal likelihood function:
1
1
exp(−Φ(u)) ∝ exp(− ku − 2 sin(2πt)k22 ) + exp(− ku + 2 sin(2πt)k22 ).
2
2
We drew 5 × 105 samples with the mixture based IS algorithm and with the pCN. We
plot the mean of the samples drawn by both method in 5.9 (left), and in 5.9 (right),
we plot 100 samples drawn by each algorithms. It can be seen from the figures that,
both methods can only capture one mode of the posterior distribution, indicating
that, the problem becomes challenging for our method and the pCN when the modes
of the target distribution are far apart.
5.3. Inverse heat conduction under model uncertainty. Our last example
is the inverse heat conduction (IHC) problems, which consist of estimating temperature or heat flux density on an inaccessible boundary from a measured temperature
history inside a solid. These problems have been studied over several decades due
to their importance in a variety of scientific and engineering applications [4]. The
IHC problems become nonlinear if the thermal properties are temperature dependent, where the inversion is significantly more difficult than the linear ones. In this
example we consider a one-dimensional heat conduction equation


∂
∂u
∂u
=
c(u)
,
(5.3)
∂t
∂x
∂x
with initial u(x, 0) = uo (x). Here x and t are the spatial and temporal variable, u(x, t)
is the temperature, and c(u) is the temperature dependent thermal conductivity, and
the length of the medium is L, all in dimensionless units. We now assume that a heat
flux is injected through the left boundary (x = 0), yielding a Neumann boundary
condition:
∂
u(0, t) = q(t).
∂x
The boundary condition (BC) at x = L is subject to uncertainty: with probability
0.8 it is
∂
u(L, t) = 0,
∂x

(5.4a)

∂
u(L, t) = −u.
∂x

(5.4b)

and with probability 0.2 it is

The interpretation is that, the system has two possible states: one with a perfectly
insulted boundary at x = L, and the other has heat diffusion at x = L.
Suppose that we place a temperature sensor in the medium (x = xs ) and the goal
is to infer the heat flux q(t) for t ∈ [0, T ], from the temperature history measured by

18

Z. FENG AND J. LI

L

q(t)

sensor

Fig. 5.10. Schematic diagram of the IHC problem.

the sensor in the time interval. The schematic of this problem is shown in Fig. 5.10.
A similar problem without model uncertainty has been studied in [22].
In the simulation, we let L = 1, T = 2, c(u) = u2 + 1, xs = 0.9, and the initial
condition be uo (x) = 0. The temperature is measured 50 times (equally spaced) and
the error in each measurement is assumed to be an independent zero-mean Gaussian
random variable with variance 0.12 . We assume the prior on q(t) is a stationary
zero-mean Gaussian process with a squared exponential covariance function:
C(t, t0 ) = exp(−

|t − t0 |2
),
2d2

(5.5)

where d = 0.3. The ”truth flux” q(t) is a realization of the prior (shown in Fig. 5.12)
and the data is simulated with the generated flux q(t) and the boundary condition (5.4b). In this problem the likelihood function becomes:
dµy
= 0.8 exp(−Φ1 (u)) + 0.2 exp(−Φ2 (u)),
dµ0
where Φ1 (u) corresponds to Eq. (5.3) with BC (5.4a) and Φ2 (u) corresponds to
Eq. (5.3) with BC (5.4b).
We draw samples from the posterior of u(t) with the four MCMC schemes used
in the previous examples. In each MCMC scheme, 1.5 × 105 draws are generated. In
both of the adaptive IS methods, the parameters are updated after every 500 draws,
and the adaptation is terminated after 105 draws. To accelerate the convergence,
we use tempering in the first 11 iterations (5,500 draws) with tempering parameter
λ = (i − 1)/10 for i = 1...11 . In the RW-pCN, we choose β = 0.1.
We first show the trace plot of the OMF in Fig. (5.11), and it is quite clear that
the results of the two adaptive methods are better than those of the prior based IS
and the pCN. Because of the multimodality of the likelihood function, the posterior
may have multiple modes, and to verify this, we apply the K-means method described
in Section 4 to cluster the samples drawn by the four methods. The samples of the
adaptive IS with mixtures can be successfully classified into two groups and we plotted
the mean of each group in Fig. 5.12 (left), compared against the true heat flux. The
K-means method, however, fails to separate the samples drawn by the other three
methods, likely because the chains have not reached the target posterior distribution
yet. We plot the means of the samples of the three methods in Fig. 5.12 (right). Like
the previous examples we show the ACF results of the four methods in Figs. 5.13, and
the acceptance rates and the ESS in Figs 5.14. In all the plots, the adaptive IS with
mixtures exhibits the best performance, followed by the IS with a single Gaussian.
6. Conclusions. In conclusion, we have presented an adaptive IS algorithm for
infinite dimensional Bayesian inference. Namely we choose a Gaussian mixture with
a particular parametrization as our proposal, and adaptively adjust the parameter
values using sample history. We prove that the proposed algorithm is well-defined

19

ADAPTIVE MCMC FOR INFINITE DIMENSIONAL INFERENCE

Prior

Mixture

40

40

30

30

20

20

10

10

0

5

10

0

15

5

10

4

x 10

Gaussian

pCN

40

40

30

30

20

20

10

10

0

5

15
4

x 10

10

0

15

5

4

10

x 10

15
4

x 10

1

1

0.5

0.5

u(t)

u(t)

Fig. 5.11. (for example 3) The trace plots of the OMF for the four different MCMC schemes.

0

0
Truth
Prior
Gaussian
pCN

−0.5

−0.5
truth
Mixture
−1
0

0.5

1
t

1.5

2

−1
0

0.5

1
t

1.5

2

Fig. 5.12. (for example 3) Left: the means of the samples in each cluster of the chain drawn
by the IS with mixtures and the true flux. Right: the means of the samples drawn by the adaptive
IS with a single Gaussian, the prior based IS and the RW-pCN.

in function space and thus is dimension independent. We also develop an efficient
algorithm based on clustering to compute the parameter values in each iteration. We
demonstrate the efficiency of the proposed method with numerical examples and in
particular we show that it performs well for multimodal posteriors. We emphasize
that, the proposed method is easy to implement, treating the problem as a black box
model, and requiring no information on the mathematical structure of the forward
model.
As has been demonstrated by the numerical examples, the mixture proposals
can generally provide faster mixing rates than the single Gaussian, thanks to its
higher flexibility. On the other hand, given that the Gaussian approximation is less

20

Z. FENG AND J. LI

1

1

0.6
Prior
Mixture
Gaussian
pCN

0.4
0.2
0
0

Prior
Mixture
Gaussian
pCN

0.95
ACF (lag 1)

ACF

0.8

50

100
lag

0.9
0.85
0.8

150

200

0.75
0

0.4

0.8

1.2

1.6

2

t

Fig. 5.13. (for example 3) Autocorrelation functions (ACF) for the four different MCMC
methods. Left: ACF of the OMF plotted as a function of lags. Right: the lag 1 ACF for u at each
grid point.
4

0.25

1.8

0.2

1.6

0.15
Prior
Mixture
Gaussian
pCN

0.1
0.05
0
50
100
iterations (× 103)

150

ESS

Acceptance rate

2

x 10

1.4

Prior
Mixture
Gaussian
pCN

1.2
1
0.8
0.6

0.4

0.8

1.2

1.6

2

t

Fig. 5.14. (for example 3) Left: the acceptance rate of the four MCMC schemes. Right: the
ESS at each grid point.

complex computationally (without the clustering step), we recommend to use the
single Gaussian approximation in problems where the posterior distributions do not
deviate too much from a Gaussian measure, and to use mixtures for strongly nonGaussian posteriors.
There are number of possible extensions of the work. First in this work we
approximate the solution to the KLD minimization problem with clustering. It is
possible that, if we can modify the standard EM algorithm and use it to solve the
optimization problem directly, we may obtain better a mixture proposal in each iteration and improve the sampling efficiency. Secondly, the intrinsic dimensionality K
is of essential importance for our method, and in the present work, K is determined
rather heuristically. Thus developments of more effective and theoretically justified
methods certainly deserve further studies. Finally, the algorithm developed here is
based on independence sampler, and we are also interested in extending the ideas to
the development of adaptive infinite dimensional random walk algorithms. We plan
to investigate these problems in the future.
REFERENCES

ADAPTIVE MCMC FOR INFINITE DIMENSIONAL INFERENCE

21

[1] C. Andrieu and J. Thoms, A tutorial on adaptive MCMC, Statistics and Computing, 18
(2008), pp. 343–373.
[2] Y. Atchade, G. Fort, E. Moulines, and P. Priouret, Adaptive markov Chain Monte Carlo:
theory and methods, Preprint, (2009).
[3] Y. F. Atchade, An adaptive version for the Metropolis adjusted Langevin algorithm with a
truncated drift, Methodology and Computing in applied Probability, 8 (2006), pp. 235–254.
[4] J. Beck, C. St Clair, and B. Blackwell, Inverse heat conduction: Ill-posed problems, John
Wiley and Sons Inc., New York, NY, 1985.
[5] A. Beskos, A stable manifold MCMC method for high dimensions, Statistics & Probability
Letters, 90 (2014), pp. 46 – 52.
[6] A. Beskos, G. Roberts, A. Stuart, et al., Optimal scalings for local metropolis–hastings
chains on nonproduct targets in high dimensions, The Annals of Applied Probability, 19
(2009), pp. 863–898.
[7] A. Beskos, G. Roberts, A. Stuart, and J. Voss, Mcmc methods for diffusion bridges,
Stochastics and Dynamics, 8 (2008), pp. 319–350.
[8] S. L. Cotter, G. O. Roberts, A. Stuart, D. White, et al., MCMC methods for functions:
modifying old algorithms to make them faster, Statistical Science, 28 (2013), pp. 424–446.
[9] T. Cui, K. J. Law, and Y. M. Marzouk, Dimension-independent likelihood-informed MCMC,
arXiv preprint arXiv:1411.3688, (2014).
[10] G. Da Prato, An introduction to infinite-dimensional analysis, Springer, 2006.
[11] M. Dashti, K. J. H. Law, A. M. Stuart, and J. Voss, MAP estimators and their consistency
in Bayesian nonparametric inverse problems, Inverse Problems, 29 (2013), p. 095017.
[12] J. Gåsemyr, On an adaptive version of the Metropolis–Hastings algorithm with independent
proposal distribution, Scandinavian Journal of Statistics, 30 (2003), pp. 159–173.
[13] P. Giordani and R. Kohn, Adaptive independent Metropolis–Hastings by fast estimation
of mixtures of normals, Journal of Computational and Graphical Statistics, 19 (2010),
pp. 243–259.
[14] H. Haario, E. Saksman, and J. Tamminen, An adaptive Metropolis algorithm, Bernoulli,
(2001), pp. 223–242.
[15] N. L. Hjort, C. Holmes, P. Müller, and S. G. Walker, Bayesian nonparametrics, vol. 28,
Cambridge University Press, 2010.
[16] L. Holden, R. Hauge, and M. Holden, Adaptive independent Metropolis-Hastings, The Annals of Applied Probability, (2009), pp. 395–413.
[17] J. Kaipio and E. Somersalo, Statistical and computational inverse problems, vol. 160,
Springer, 2005.
[18] R. E. Kass, B. P. Carlin, A. Gelman, and R. M. Neal, Markov Chain Monte Carlo in
Practice: A Roundtable Discussion, The American Statistician, 52 (1998), pp. 93–100.
[19] J. M. Keith, D. P. Kroese, and G. Y. Sofronov, Adaptive independence samplers, Statistics
and Computing, 18 (2008), pp. 409–420.
[20] K. J. Law, Proposals which speed up function-space mcmc, Journal of Computational and
Applied Mathematics, 262 (2014), pp. 127–138.
[21] J. Li, A note on the Karhunene-Loeve expansions for infinite-dimensional Bayesian inverse
problems, Statistics and Probability Letters, 106 (2015), pp. 1 – 4.
[22] J. Li and Y. M. Marzouk, Adaptive construction of surrogates for the bayesian solution of
inverse problems, SIAM Journal on Scientific Computing, 36 (2014), pp. A1163–A1186.
[23] J. S. Liu, Metropolized independent sampling with comparisons to rejection sampling and importance sampling, Statistics and Computing, 6 (1996), pp. 113–119.
[24] T. Marshall and G. Roberts, An adaptive approach to Langevin MCMC, Statistics and
Computing, 22 (2012), pp. 1041–1057.
[25] 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.
[26] J. C. Mattingly, N. S. Pillai, A. M. Stuart, et al., Diffusion limits of the random walk
metropolis algorithm in high dimensions, The Annals of Applied Probability, 22 (2012),
pp. 881–930.
[27] G. McLachlan and D. Peel, Finite mixture models, John Wiley & Sons, 2004.
[28] N. Petra, J. Martin, G. Stadler, and O. Ghattas, A computational framework for infinitedimensional 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.
[29] F. Pinski, G. Simpson, A. Stuart, and H. Weber, Kullback-Leibler approximation for probability measures on infinite dimensional spaces, arXiv preprint arXiv:1310.7845, (2013).

22

Z. FENG AND J. LI

[30] F. J. Pinski, G. Simpson, A. M. Stuart, and H. Weber, Algorithms for kullback-leibler approximation of probability measures in infinite dimensions, arXiv preprint arXiv:1408.1920,
(2014).
[31] G. O. Roberts, A. Gelman, W. R. Gilks, et al., Weak convergence and optimal scaling of
random walk metropolis algorithms, The annals of applied probability, 7 (1997), pp. 110–
120.
[32] G. O. Roberts and J. S. Rosenthal, Examples of adaptive MCMC, Journal of Computational
and Graphical Statistics, 18 (2009), pp. 349–367.
[33] G. O. Roberts, J. S. Rosenthal, et al., Optimal scaling for various metropolis-hastings
algorithms, Statistical science, 16 (2001), pp. 351–367.
[34] A. M. Stuart, Inverse problems: a Bayesian perspective, Acta Numerica, 19 (2010), pp. 451–
559.
[35] L. Tierney, Markov chains for exploring posterior distributions, the Annals of Statistics,
(1994), pp. 1701–1728.
, A note on metropolis-hastings kernels for general state spaces, Annals of Applied Prob[36]
ability, (1998), pp. 1–9.
[37] C. J. Wu, On the convergence properties of the EM algorithm, The Annals of statistics, (1983),
pp. 95–103.

