Monte Carlo methods, such as Markov chain Monte Carlo (MCMC) algorithms, have
become very popular in signal processing over the last years. In this work, we
introduce a novel MCMC scheme where parallel MCMC chains interact, adapting
cooperatively the parameters of their proposal functions. Furthermore, the
novel algorithm distributes the computational effort adaptively, rewarding the
chains which are providing better performance and, possibly even stopping other
ones. These extinct chains can be reactivated if the algorithm considers
necessary. Numerical simulations shows the benefits of the novel scheme

Elvira, V.

Louzada, F.

Luengo, D.

Martino, L.

English

arXiv

L. Martino

V Elvira

D. Luengo

F Louzada

Crossref

Parallel metropolis chains with cooperative adaptation

Monte Carlo methods, such as Markov chain Monte Carlo (MCMC) algorithms, have become very popular in signal processing over the last years. In this work, we introduce a novel MCMC scheme where parallel MCMC chains interact, adapting cooperatively the parameters of their proposal functions. Furthermore, the novel algorithm distributes the computational effort adaptively, rewarding the chains which are providing better performance and, possibly even stopping other ones. These extinct chains can be reactivated if the algorithm considers it necessary. Numerical simulations show the benefits of the novel scheme

Martino, Luca

Elvira Arregui, Víctor

Luengo García, David

Louzada, Francisco

Servicio de Coordinación de Bibliotecas de la Universidad Politécnica de Madrid

PARALLEL METROPOLIS CHAINS WITH COOPERATIVE ADAPTATION 
L. Martino , V. Elvira , D. Luengo , F. Louzada 
Institute of Mathematical Sciences and Computing, Universidade de Sao Paulo, Sao Carlos (Brazil). 
Dep. of Signal Theory and Communic, Universidad Carlos III de Madrid, Leganes (Spain). 
Dep. of Signal Theory and Communic, Universidad Politecnica de Madrid, Madrid (Spain). 
ABSTRACT 
Monte Carlo methods, such as Markov chain Monte Carlo 
(MCMC) algorithms, have become very popular in signal pro-
cessing over the last years. In this work, we introduce a 
novel MCMC scheme where parallel MCMC chains inter-
act, adapting cooperatively the parameters of their proposal 
functions. Furthermore, the novel algorithm distributes the 
computational effort adaptively, rewarding the chains which 
are providing better performance and, possibly even stopping 
other ones. These extinct chains can be reactivated if the al-
gorithm considers it necessary. Numerical simulations show 
the benefits of the novel scheme. 
Index Terms— Interacting Parallel MCMC, Adaptive 
MCMC, cooperative adaptation, Bayesian inference. 
1. INTRODUCTION 
Markov Chain Monte Carlo (MCMC) algorithms [1] are 
widely employed in signal processing and communications 
for Bayesian inference and optimization [2, 3, 4, 5, 6]. 
They draw random samples from a complicated multidimen-
sional target probability density function (pdf), 7r(x) with 
x G V C Rd, generating a Markov chain which converges 
to 7r(x). The performance, i.e., the speed of the convergence, 
depends strictly on the choice of a suitable proposal function 
<z(x), and more specifically, on the discrepancy between q(x.) 
and7r(x). 
Speeding up the convergence has motivated an intense re-
search activity. On the one hand, one active research line 
considers the design of an adaptive proposal density within 
MCMC techniques. Several schemes have been developed 
in order to tune online the parameters of the proposal den-
sity, learning them from the previously generated samples 
[7, 8, 5, 6]. On the other hand, the use of several parallel 
chains instead of a single long chain has been studied for 
different reasons. First of all, employing parallel chains al-
lows the use of different proposal pdfs. Moreover, another 
This work has been partially supported by the Spanish government 
through the DISSECT (TEC2012-38058-C03-01) project, the BBVA Foun-
dation through the MG-FIAR project, by the Grant 2014/23160-6 of Sao 
Paulo Research Foundation (FAPESP) and by the Grant 305361/2013-3 of 
the National Council for Scientific and Technological Development (CNPq). 
important motivation is the interest in the implementation of 
MCMC techniques within a parallel architecture [9, 10, 11]. 
Finally, a third reason is to speed up the exploration of the 
state space [12, 4, 13, 14, 11]. Several works in the literature 
are focused on producing an interaction among the different 
parallel chains [15, 7, 13, 14, 16]. The exchange of infor-
mation among the chains can yield jointly two related bene-
fits: produce a faster convergence of the chains to the target 
(e.g., reallocating "lost" chains around a mode of the target 
[11]), and help the cooperative exploration of the state space 
(for instance, generating a repulsion among the chains dur-
ing a certain number of iterations [14]). As a consequence, 
the combined use of interacting parallel chains and adaptive 
proposal pdfs is considered of great interest in the literature 
[7]. 
In this work, we propose a novel scheme, called paral-
lel adaptive independent Metropolis (PAIM), involving par-
allel MCMC chains which exchange information in order to 
adapt online their proposal densities. Namely, the interac-
tion among the chains is carried out by a cooperative adapta-
tion of the proposal functions. Each MCMC chain employs a 
proposal pdf, independent of the previous state of the chain, 
which are formed by a mixture of two densities, each one de-
termined by two parameters: a mean vector and covariance 
matrix. All the parameters are updated using empirical esti-
mators (as in [8, 5]) applied to the complete set or only to a 
subset of the previously generated states. On the one hand, 
the first component of each mixture aims to provide a global 
adaptation, the complete set of states is used. On the other 
hand, the second component of each mixture is adapted con-
sidering only a subset of the previously generated states in 
order to learn local features of the target pdf. These sub-
sets of states are built using a simple clustering-type strategy. 
This generates a cooperative adaptation, preventing that dif-
ferent proposals cover the same region, and allowing them to 
cover different modes of the target function, for instance. Fur-
thermore, the novel algorithm is able to identify the proposal 
functions better located, and to allocate more computational 
effort to the corresponding chains. Indeed, PAIM adapts the 
number of iterations of every chain, stopping those chains 
which are using a badly located proposal. PAIM is also able 
to turn back on certain chains when it is necessary. 
2. GENERAL SETUP 
In many different applications [12, 3, 4], it is necessary 
to draw samples from a complicated d-dimensional tar-
get probability density function (pdf), 7f(x) oc 7r(x), with 
x e D C l d . For this purpose, we consider the use of N 
parallel Metropolis-Hastings chains [1], each one employ-
ing a different independent proposal pdf ipn(x), defined as a 
mixture of two pdfs </iiK(x) and </2,n(x). More specifically, 
explicitly indicating the parameters of the two components, 
we have 
1 1 
VVi(x) = -</l j n(x|/LXi j n , Cltn) + T ^ , n (x|/i>2,n, ^2,n), — c 
2 ' 
—i 
2 ' 
for n = 1 , . . . , N, where /Ltjjn represents a d-dimensional 
mean vector and C j „ is a d x d covariance matrix for i G 
{ 1 , 2}. In the novel method, the N parallel chains interact 
in order to jointly adapt all the parameters, /LXJ„ and Cj „ for 
«€ {1 ,2} and Vn. The parameters /Ltijn and C i „ , of the first 
component q\n of every proposal tpn are updated to provide 
a global adaptation, whereas the parameters /Lt2,n and C2 j n of 
</2,n are adapted to extract local features of the target pdf. 
Furthermore, in the novel scheme, each chain performs 
a different number of iterations, Kn for n = 1 , . . . , N. In-
deed, the proposed algorithm is also able to determine online 
a subset of the N chains which are obtaining the best per-
formance. These chains are employed more often than the 
rest. Namely, the total number of iterations Kn performed by 
the corresponding chains is increased. Let us denote the to-
tal number of desired samples as L, chosen in advance by the 
user. Then, we have K\ + K2... + K^ = L, i.e., PAIM is 
stopped when L samples have been generated.1 
3. ADAPTIVE PARALLEL INDEPENDENT 
METROPOLIS ALGORITHM 
The parallel adaptive independent Metropolis (PAIM) algo-
rithm works on two different time scales. First of all, the in-
dex t denotes the current step of the algorithm, with 
t = 0 , . . . , Ttot, (1) 
where the total number of steps Ttot is not decided by the user, 
but automatically tuned by PAIM (we have always Ttot < L). 
Furthermore, we have have a different iteration index kn for 
each chain, such that 
kn = 0 , . . . , Kn, (2) 
for n = 1 , . . . , N. The value of each Kn is also decided 
by PAIM (recall that z^n=i K-n = T). At each step, the set 
At contains the indices corresponding to the active chains. 
At every t-th step, all the active chains perform one iteration 
(i.e., i f the j-th chain is active, then kj = kj +1), whereas the 
inactive chains remain frozen. For every t < Ttrain, where 
Ttrain is chosen by the user, all the chains are active, i.e., 
At = {1, 2, . . . , N}, t < Ttrain. 
Whereas for t > Ttrain we have At Q { 1 , 2 , . . . , N}. The in-
teracting adaptation is performed at any t such that Ttrain < 
t < Tstop. We consider the possibility of stopping the adap-
tation after Tstop steps, since the adaptation could jeopardize 
the ergodicity of the chains. However, the numerical results 
described in Section 5 show that the algorithm seems to main-
tain the correct ergodicity properties. 
The adaptation is performed as follows. During the first 
Ttrain time steps the algorithm simply assigns the new cur-
rent states, generated at the t-th step, to one chain among the 
N possible, according to the minimum Euclidean distance be-
tween them and the means, /L4 „ for n = 1, . . . , N. Thus, we 
allow the method to use a few iterations (t = 1 , . . . , Train) to 
collect information about the target, as in [8, 5]. Afterwards, 
the algorithm adapts all the parameters, /J,\
 n and C^ n for 
i G { 1 , 2}, and the set of active chains At, until t = Tstop. 
On the one hand, the parameters ii2 n and C;> n are updated 
using the empirical estimators for the means and covariances 
considering only the samples assigned to the n-th chain. On 
the other hand, the parameters /J,\
 n and C\'n are adapted 
considered all the states generated so far. The chains are 
turned on or off taking into account the number of samples 
assigned to the n-th chain. PAIM is detailed below. 
1. Initialization: 
1.1- Parameters: choose the number of chains N and the 
desired samples L. Select the positive values e, Ttrain, 
Tstop,2 the initial parameters, [i\
 n and C^ 'n for i G 
{1,2} , the initial states x n o, and the counters rnn = 1 
for n = 1 , . . . , N. Define the initial set of the active 
chains as *4o = { 1 , 2 , . . . , N}. 
1.2- Indices: set £ = 0, t = — 1 , and kn = 0 for n = 
1 , . . . , N. 
2. M H steps: 
2.1- Set t = t + 1, Z = 0 and i = 0. 
2.2- For all the active chains, i.e., for all the indices j & At: 
(a) Sample x ' from the j-th proposal pdf, 
2 
/ / (t) \ J- V ^ (t), , (t) „ ( 
x ~ ipy (xj = — y. il j (,xl As j i ^l 
1Note that we are including all the states in the “burn-in” periods of the 
N chains. However, they can be discarded if desired. 2We recall that it is necessary to set T t r a i n < T s t o p . 
(b) Accept X j fc i + i = x ' with probability 4. FURTHER CONSIDERATIONS ABOUT PAIM 
a = m m 1, 
fc) 
3 ' (3) 
Otherwise, s e t x ^ . + i = xjjkj. 
(c) Set Z = Z U { Z J + I = Xj^.-i- i}, and i = i + 1. 
(d) Set #£+i = x^fc.-i-i, kj = kj + 1 and^ = ^ + 1. 
(e) Stop condition: If £ > L then stop and return 
{ # i , . . . , 8L}. 
3. Assignment (if t < Tstop): 
3.1- For all the vectors z r e 2 , i.e., for r = 1 , . . . , \Z\ : 
(a) Find the closest mean to z$ (w.r.t. the Euclidean 
distance), i.e., find the index 
n* = argmin ||/i,2 „ — z r ||2. (4) 
n ' 
(b) Set s n * j m r i t + i = Zj and mn* = mn* + 1. 
4. Adaptation (if Ttrain <t < Tstop): 
4.1- Set At+i = 0. 
4.2- For n = 1,... ,N, update the mean vectors, 
<t+i) ~(v+-n 1 v ^ /1 
Ml n I-1 H / ®ii 
(t+1) 1 
(5) 
i = l 
update the covariance matrices, 
Ci = C"+ ' = y (Oi — fi"+ ') + eld 
£ — 1 -^^ i=i 
C ( t+i) 1 
TOn — 1 
/ (sn,i — 1^2 n ) + £ Idi 
where Id is d x d unit matrix. Moreover, if 
\ AT m « A 
a n = iv —T7 > U, 
5^7 =1 mj 
(6) 
(7) 
then set At+i = At+i U { n } (where [a] denotes the 
smallest integer larger or equal than a, with a e R). 
Otherwise, i f an = 0 the n-th chain is deactivated, i.e., 
not used at the next iteration. 
5. Repeat from step 2. 
Observe that the set Z contains all the new states generated 
in one specific step. This is refreshed at the beginning of ev-
ery step t. Since the number of active chains is variable, the 
cardinality of Z is also changing with t. Moreover, we have 
denoted with sn j the «-th state assigned to the n-th chain. 
The counter mn indicates the number of states associated to 
the n-th chain. Note also that for all t > Ttrain we have 
(*) 
Ml,n M w C\ = C"-
1
, for n = 1 , . . . , N. 
Namely, all the functions </in are updated using the empirical 
estimators of the mean and covariance of the target obtained 
from all the generated samples. 
It is important to remark that PAIM is able to distribute 
the computational efforts efficiently. Indeed, let us consider 
the initial use of a huge number N of parallel chains, with 
proposal pdfs localized randomly over the state space. A l l the 
chains using a badly located proposal pdf would be quickly 
deactivated, i.e., only the chains with proposal pdfs located 
close to high probability regions would survive. However, the 
algorithm is also able to start up again certain chains if, after 
some steps, new states have been assigned to them. Finally 
note that, in the description of the algorithm, the parameters 
are updated using a block procedure, but efficient recursive 
update formulas can be employed (e.g., see [5]), so that PAIM 
can be efficiently applied in high dimensional problems. 
5. NUMERICAL RESULTS 
We consider a bi-dimensional “banana-shaped” target distri-
bution [8], which is a benchmark function commonly used in 
the literature. Mathematically, it is given by 
7r(xi , X2) oc e x p 
1 
(4 — Bx\ — x 2 ) 
2 ^ 2ri 
where, we have set B = 10, 771 = 4, r/2 = 5, and rjs = 
5. The goal is to estimate the expected value i?[X], where 
X = [X i ,X2] ~ 7f(xi,X2). We compute the true value 
E[X\ « [—0.4845,0]T approximately by using an exhaus-
tive deterministic numerical method (with an extremely thin 
grid), in order to obtain the mean square error (MSE) of PAIM 
and the corresponding independent parallel M H chains with 
the same initial parameters but no adaptation. 
We consider N s {5,10, 50,100} chains with Gaus-
sian proposals, </jin(x|/Lx^„, C^„) for i e {1,2} and n = 
1 , . . . , N. The initial means and the initial states are chosen 
randomly at each run. More specifically, we have x n 0 ~ 
U{[—15, —15] x [—15,15]) and ii\n ~ W([—15, —15] x 
[-15,15]). The initial covariance matrices are C\l = 
[a2 0; 0 a2]T with a = 10. We consider L = 5000 total 
number of samples. For PAIM, we test Ttrain & {1,10,20}, 
10 
0 
-10 
-10 0 
(a) 
10 
10 
5 
0 
-5 
-10 
-10 -5 0 5 10 
(b) 
Fig. 1. (a) Contour plot of the banana-shaped target 7¯ , and the initial means [v2 n (circles) with N = 20. (b) We show the 
generated samples (dots), the initial (circles) and final means (x-marks) of the second component </2,n of the proposal of the 
final active chains. The final covariance ellipsoids are also depicted. The final mean of the first component of the proposals, 
Ml,n 
p ) 
V {TstoP) is shown with a square jointly with the corresponding covariance ellipsoid with dashed line. 
50  
40  
30 
20 
10 
Table 1 . Percentage of reduction in the MSE obtained using 
PAIM, with respect to IPCs. 
Fig. 2. The indices of the active chains (denoted with black 
marks) at each step t of PAIM in one specific run, setting 
N = 50, Ttrain = 2, L = 1000. In this run, the final number 
of active chains is 16 and Ttot = 38. 
setting e = 0.4 and Tstop = oo (i.e., we never stop the 
adaptation). The results are averaged over 500 independent 
simulations for each combination of parameters. PAIM al-
ways provides a smaller mean square error (MSE) in the 
estimation of E[X] (averaging the two components) with re-
spect to the corresponding independent parallel chains (IPCs). 
Table 1 shows the percentage of reduction in MSE obtained 
using PAIM with different values of N and Ttrain, which 
ranges from 22 % up to 60 % approximately. Indeed, in this 
example the minimum train period (Ttrain = 1) provides 
the best results. Figure 1(a) shows the initial means of the 
second components of the proposals in PAIM, i.e., ii2 „ , and 
J- train N = 5 N = 10 N = 50 N = 100 
1 51.34% 58.05% 63.78% 60.31% 
10 46.95% 41.73% 44.23% 35.65% 
20 29.81% 35.51% 33.29% 22.33% 
the contour plot of the target n. Figure 1(b) depict the initial 
(circles) and final (x-marks) configurations of the means of 
the second components of the final active proposals, i.e., /L4 „ 
and [i\ n°P , obtained in one specific run (setting N = 20 
and L = 5 x 104). The final mean /J,\ ^°p = fH,^Tstop) of the 
first component common to all the proposals is depicted with 
a square. The figures also show with solid line the covariance 
ellipsoids, corresponding to the « 90% of the probability 
mass, of the second components of the proposals of the final 
active chains. The covariance ellipsoid corresponding to the 
first common component of the proposal pdfs is displayed 
with a dashed line. Figure 2 shows the indices corresponding 
to the active chains, as a function of t, in one specific run 
(N = 50 and L = 103). 
6. CONCLUSIONS 
In this work, we have introduced a new interacting parallel 
MCMC scheme, where a cooperative adaptation of the pro-
posal densities is performed. Furthermore, the computational 
effort is efficiently distributed by the the novel method among 
the set of parallel chains, according to their performance. 
7. REFERENCES 
[1] C. P. Robert and G. Casella, Monte Carlo Statistical 
Methods, Springer, 2004. 
[2] A. Doucet and X. Wang, “Monte Carlo methods for 
signal processing,” IEEE Signal Processing Magazine, 
vol. 22, no. 6, pp. 152–170, Nov. 2005. 
[3] W. J. Fitzgerald, “Markov chain Monte Carlo methods 
with applications to signal processing,” Signal Process­
ing, vol. 81, no. 1, pp. 3–18, January 2001. 
[4] A. Jasra, D. A. Stephens, and C. C. Holmes, “On 
population-based simulation for static inference,” Statis­
tics and Computing, vol. 17, no. 3, pp. 263–279, 2007. 
[5] D. Luengo and L. Martino, “Fully adaptive Gaussian 
mixture Metropolis-Hastings algorithm,” in IEEE Inter­
national Conference on Acoustics, Speech, and Signal 
Processing (ICASSP), 2013, pp. 6148–6152. 
[6] L. Martino, J. Read, and D. Luengo, “Independent dou­
bly adaptive rejection Metropolis sampling within Gibbs 
sampling,” IEEE Transactions on Signal Processing, 
vol. 63, no. 12, pp. 3123–3138, 2015. 
[7] R. Craiu, J. Rosenthal, and C. Yang, “Learn from thy 
neighbor: Parallel-chain and regional adaptive MCMC,” 
Journal of the American Statistical Association, vol. 
104, no. 448, pp. 1454–1466, 2009. 
[8] Heikki Haario, Eero Saksman, and Johanna Tamminen, 
“An adaptive Metropolis algorithm,” Bernoulli, vol. 7, 
no. 2, pp. 223–242, April 2001. 
[9] B. Calderhead, “A general construction for paralleliz­
ing Metropolis-Hastings algorithms,” Proceedings of 
the National Academy of Sciences of the United States 
of America (PNAS), vol. 111, no. 49, pp. 17408–17413, 
2014. 
[10] P. Jacob, C. P. Robert, and M. H. Smith, “Using par­
allel computation to improve Independent Metropolis-
Hastings based estimation,” Journal of Computational 
and Graphical Statistics, vol. 3, no. 20, pp. 616–635, 
2011. 
[11] L. Martino, D. Luengo V. Elvira, J. Corander, and 
F. Louzada, “Orthogonal parallel MCMC methods for 
sampling and optimization,” arXiv:1507.08577, 2015. 
[12] J. Corander, M. Ekdahl, and T. Koski, “Parallel inter­
acting MCMC for learning of topologies of graphical 
models,” Data Mining and Knowledge Discovery, vol. 
17, no. 3, pp. 431–456, 2008. 
[13] L. Martino, V. Elvira, D. Luengo, A. Arte´s-Rodr´ ıguez, 
and J. Corander, “Orthogonal MCMC algorithms,” 
IEEE Workshop on Statistical Signal Processing (SSP), 
pp. 364–367, June 2014. 
[14] L. Martino, V. Elvira, D. Luengo, A. Arte´s-Rodr´ ıguez, 
and J. Corander, “Smelly parallel MCMC chains,” IEEE 
International Conference on Acoustics, Speech, and Sig­
nal Processing (ICASSP), pp. 4070–4074, 2015. 
[15] R. Casarin, R. V. Craiu, and F. Leisen, “Interacting 
multiple try algorithms with different proposal distribu­
tions,” Statistics and Computing, vol. 23, no. 2, pp. 185– 
200, 2013. 
[16] L. Martino, V. Elvira, D. Luengo, and J. Corander, “Lay­
ered adaptive importance sampling,” arXiv:1505.04732, 
May 2015. 


Archivo Digital UPM

Parallel Metropolis chains with cooperative adaptation

Parallel Metropolis chains with cooperative adaptation

Abstract

Similar works

Full text

Available Versions

Crossref

Servicio de Coordinación de Bibliotecas de la Universidad Politécnica de Madrid

Archivo Digital UPM