Using parallel computation to improve Independent Metropolis–Hastings based estimation

We recall here the basic IMH algorithm, assuming the availability of a proposal distribution that we can sample, and which probability density $\mu$ is known up to a normalization constant. The independent Metropolis–Hastings algorithm, described in Algorithm 1, generates a Markov chain with invariant density $\pi$, corresponding to the target distribution.

In the larger picture of Monte Carlo and MCMC algorithms, the IMH algorithm holds a rather special status. It has certainly been studied more often than other MCMC schemes (Robert and Casella, 2004), but it is undoubtedly a less practical solution than the more generic random walk Metropolis–Hastings algorithm. For instance, it is rather rarely used by itself because it requires the derivation of a tolerably good approximation to the true target, approximation that most often is unavailable. On the other hand, first-order approximations and Metropolis-within-Gibbs schemes are not foreign to calling for IMH local moves based on Gaussian representations of the targets. The reason theoretical studies of the IMH algorithm abound is that it has strong links with the non-Markovian simulation methods such as importance sampling. Contrary to random-walk Metropolis–Hastings schemes, IMH algorithms may enjoy very good convergence properties and may also reach acceptance probabilities that are close to one. Furthermore, the potentially large gain in variance reduction provided by the parallelization scheme developped in this paper may counteract the lesser efficiency of the original IMH compared with the random walk Metropolis–Hastings algorithm.

An important feature of the IMH algorithm, when addressing parallelism, is that it cannot work but in an iterative manner, since the outcome of step $t$, namely the value $x_{t}$, is required to compute the acceptance ratio at step $t+1$. This sequential construction is compulsory for the validation of the algorithm given the Markov property at its core (Robert and Casella, 2004). Nonetheless, given that, in the IMH algorithm, the proposed values $(y_{t})$ are generated independently from the current state of the Markov chain, $x_{t}$, it is altogether possible to envision the generation of $T$ proposed values $y_{t}$ first, along with the computation of the associated ratios $\omega_{t}=\pi(y_{t})/\mu(y_{t})$. Once this computation requirement is completed, only the acceptance steps need to be considered iteratively. This two-step perspective makes for a huge saving in computing time when the simulation of the $y_{t}$’s and the derivation of the $\omega_{t}$’s can be achieved in parallel since both the remaining computation of the ratios $\rho(x_{t-1},y_{t})$ given the $\omega_{t}$’s and their subsequent comparison with uniform draws typically are orders of magnitude faster.

In this respect the IMH algorithm strongly differs from the random walk Metropolis–Hastings (RWMH) algorithm, for which acceptance ratios cannot be processed beforehand because the proposed simulated values depend on the current value of the Markov chain. The universal availability of parallel processing schemes may thus lead to a new surge of popularity for the IMH algorithm. Indeed, when taking advantage of $p$ parallel processing units, an IMH can be run for $p$ times as many iterations as RWMH, at almost the same computing cost since RWMH cannot be directly parallelized.

In order to better describe this increased computing power, we first note that, once $T$ successive values of a Markov chain have been produced, the sequence is usually processed as a regular Monte Carlo sample to obtain an approximation of an expectation under the target distribution, $\mathbb{E}_{\pi}\left[h(X)\right]$ say, for some arbitrary functions $h$. We propose in this paper a technique that improves the precision of the estimation of this expectation by taking advantage of parallel processing units without jeopardizing the Markov property.

Before presenting our improvement scheme, we introduce the notation $\vee$ (read “or”) for the operator that represents a single step of the IMH algorithm. Using this notation, given the current value $x_{t}$ and a sequence of $p$ independent proposed values $y_{1},\ldots,y_{p}\sim\mu$, the IMH algorithm goes from step $t$ to step $t+p$ according to the diagram in Figure 1.

We propose to take full advantage of the simulated proposed values and of the computation of their corresponding $\omega$ ratios. To this effect, we introduce the block IMH algorithm, made of successive simulations of blocks of size $p\times p$. In this alternative scheme, the number of blocks $b$ is such that the number of desired iterations $T$ is equal to $b*p$, in order to keep the comparison with a standard IMH output fair. Usually $p$ needs not be calibrated since it represents the number of physical parallel processing units that can be exploited by the code. However, in principle, this number $p$ can be set arbitrarily high and based on virtual parallel processing units, the drawback being then an increase in the computing cost. (Note that the block IMH algorithm can also be implemented with no parallel abilities, still it provides a gain in variance that may counteract the increase in time.) In the following examples, we take $p$ varying from $4$ to $100$. We first explain how a block is simulated, and then how to move from one block to the next.

A $p\times p$ block consists in the generation of $p$ parallel generations of $p$ values of Markov chains, all starting at time $t$ from the current state $x_{t}$ and all based on the same proposed simulated values $y_{1},\ldots,y_{p}$. The different between the $p$ floes is the orders in which those $y_{i}$’s are included. For instance, these orders may be the $p$ circular permutations of $y_{1},\ldots,y_{p}$, or they may be instead random permutations, as discussed in detail (and compared) in Section 3. The block IMH algorithm is illustrated in Figure 2 for the circular set of permutations.

It should be clear that each of the $p$ parallel chains found in this block is a valid MCMC sequence of length $p$ when taken separately. As such, it can be processed as a regular MCMC output. In particular, if $x_{t}$ is simulated from the stationary distribution, any of the subsequent $x^{(j)}_{t+i}$ is also simulated from the stationary distribution. However, the point of the $p$ parallel flows is double:

• •

  it aims at integrating out part of the randomness resulting from the ancillary order in which the $y_{k}$’s are chosen, getting close to the conditioning on the order statistics of the $y_{k}$’s advocated by Perron (1999);

• •

  it also aims at partly integrating out the randomness resulting from the generation of uniform variables in the selection process, since the block implementation results in drawing $p^{2}$ uniform realizations instead of $p$ uniform realizations for a standard IMH setting.

Both of those points essentially amount to implementing a new Rao–Blackwellization technique (a more precise connection is drawn in Section 4). In an independent setting, each of the $y_{k}$’s occurs a number $n_{k}\geq 0$ of times across the $p$ steps of the $p$ parallel chains, i.e. for a number $p^{2}$ of realizations. Therefore, when considering the standard estimator $\hat{\tau}_{1}$ of $\mathbb{E}_{\pi}\left[h(X)\right]$, based on a single MCMC chain,

this estimator necessarily has a larger variance than the double average

where $y_{0}:=x_{t}$ and $n_{0}$ is the number of times $x_{t}$ is repeated. (The proof for the reduction of the variance from $\hat{\tau}_{1}$ to $\hat{\tau}_{2}$ easily follows from a double integration argument.) We again insist on the compelling feature that computing $\hat{\tau}_{2}$ using $p$ parallel processing units does not cost more time than computing $\hat{\tau}_{1}$ using a single processing unit.

In order to preserve its Markov validation, the algorithm must properly continue at time $t+p$. An obvious choice is to pick one of the $p$ sequences at random and to take the corresponding $x_{t+p}^{(j)}$ as the value of $x_{t+p}$, starting point of the next parallel block. This mechanism is represented in Figure 3. While valid from a Markovian perspective, since the sequences are marginally produced by a regular IMH algorithm, this means that the chain deduced from the block IMH algorithm is converging at exactly the same speed as the original IMH algorithm. An alternative choice for the starting points of the blocks takes advantage of the weights $n_{k}$ on the $y_{k}$’s that are computed via the block structure. Indeed, those weights essentially act as importance weights and they allow for a selection of any of the $p^{2}$ $x_{t+i}^{(j)}$’s as the starting point of the incoming block, which corresponds to choosing one of the proposed $y_{k}$’s with probability proportional to $n_{k}$. While this proposal does reduce the length of the resulting chain, it does not impact the estimation aspects (which still involve all of the $p^{2}$ values) and it could improve convergence, given that the weighted $y_{k}$’s behave like a discretized version of a sample from the target density $\pi$. We will not cover this alternative any further.

The original version of the block IMH algorithm is described in Algorithm 2, The algorithm is made of a loop on the $b$ blocks and an inner loop on the $p$ parallel chains of each block. The $p$ steps of this inner loop are actually meant to be implemented in parallel. The output of Algorithm 2 is double:

• •

  a standard Markov chain of length $T$, which is made of $b$ chains of length $p$, each of which is chosen among $p$ chains at line 12 of Algorithm 2,

• •

  a $p\times T$ array $(x_{t}^{k})_{t=1:T}^{k=1:p}$, on which the estimator $\hat{\tau}_{2}$ is based.

Since the point-wise evaluation of the target density $\pi(y_{k})$ is usually the most computer-intensive part of the algorithm, sampling additional uniform variables has a negligible impact here, as do further costs related to the storage of vectors larger than in the original IMH. This is particularly compelling since the multiple chains do not need to be stored further than during a single block execution time. That is why, although we sample $p$ times more uniforms in the block IMH algorithm, we still consider it to be running at roughly of the same cost as the IMH algorithm. The number of target density evaluations indeed is the same for both and most often represent the overwhelming part of the computing time in the Metropolis–Hastings algorithm. Besides, pseudo-random generation of uniforms can also benefit from parallel processing, see e.g. L’Ecuyer et al. (2001).

In the following Monte Carlo experiment, various versions of the block IMH algorithm are compared one to another, as well as to standard IMH and importance sampling. We stress that a straightforward reason for not conducting a comparison with a plain parallel algorithm based on $p$ independent parallel chains is that it does not make much sense cost-wise. Indeed, running $p$ parallel MCMC chains of the same length $T$ does cost $p$ times more in terms of target density evaluations. Obviously, if one insists on running $p$ independent chains, for instance as to initialize an MCMC algorithm from several well-dispersed starting points, each of those chains can benefit from our stabilizing method, which will improve the resulting estimation.

The method is presented here for square blocks of dimension $(p,p)$, but blocks could be rectangular as well: the algorithm is equally valid when using $r\neq p$ permutations, leading to $(r,p)$ blocks. We focus here on square blocks because when the machine at hand provides $p$ parallel processing units, then it is most efficient to simulate the proposed values and the uniforms, and to compute the target densities and the acceptance ratios at the $p$ proposed values in parallel. Once again, the block IMH algorithm with $p\times p$ square blocks has about the same cost as the original IMH algorithm, because computing target densities and acceptance ratios does more than compensate for the cost of randomly picking an index at the end of each block. This amounts to say that line 4 of Algorithm 1 and line 5 of Algorithm 2 are (by far) the most computationally demanding ones in the respective algorithms.

We now introduce a toy example that we will follow throughout the paper. The target $\pi$ is the density of the standard $\mathcal{N}(0,1)$ normal distribution and the proposal $\mu$ is the density of the $\mathcal{C}(0,1)$ Cauchy distribution. Hence, the density ratio is

We only consider the integral $\int{x\pi(\text{d}x)}$, the mean of $\pi$ equal to zero in this case. The acceptance rate of the IMH algorithm for this example is around $70\%$. (Note that IMH with higher acceptance rates are considered to be more efficient, in opposition to other Metropolis–Hastings algorithms, see Robert and Casella, 2004.)

In all results related to the toy example presented thereafter, $10,000$ independent runs are used to compute the variance of the estimates. The value of $p$ represents the number of parallel processing units that are available, ranging from $4$ for a desktop computer to $100$ for a cluster or a graphics processing unit (GPU) (this value could even be larger for computers equipped with multiple GPUs).

The results of the simulation experiments are presented in Figures 4–8 as barplots, which indicate the percentage of variance decrease associated with the estimators under comparison, the reference estimator being the standard IMH output $\hat{\tau}_{1}$. In agreement with the block sampling perspective, the same proposed values and uniform draws were used for all the estimators that are plotted on the same graph (that is, for a given value of $p$), so that the comparison is not perturbed by an additional noise associated with the simulation.

Let $\mathcal{S}$ be the set of permutations of $\{1,\ldots,p\}$. Its size is $p!$, therefore too large to allow for averaging over all permutations, although this solution would be ideal. We consider the simpler option of finding $p$ efficient permutations in $\mathcal{S}$, denoted by $(\sigma_{1},\ldots,\sigma_{p})$, the goal being a choice favoring the largest possible decrease in the variance of the estimator $\hat{\tau}_{2}$ defined in Section 2.

We compare the five described types of permutations on the toy example. Figure 4 shows the results for various values of $p$, displaying the variance reduction of $\hat{\tau}_{2}$ associated with each of the permutation orders, compared to the variance of the original IMH estimator $\hat{\tau}_{1}$. For each of the $10,000$ independent replications, the block IMH algorithm was launched on one single $p\times p$ block, e.g. with $b=1$ using the notation of Section 2, since $b$ plays no role whatsoever in this comparison.

As mentioned above, using the same order in the $y_{k}$’s for each of the $p$ parallel chains already produces a significant decrease of about $20\%$ in the variance of the estimators. This simulation experiment shows that the three random permutations (random, half-random half-reversed and stratified) are quite equivalent in terms of variance improvement and that they are significantly better than the circular permutation proposal, which only slightly improves over the “same order” scheme. Therefore, in the next Monte Carlo experiments, we will only use the random order solution, simplest of the random schemes. An amount of improvement like $35\%$ when $p\geq 32$ is quite impressive when considering that it is essentially obtained cost-free for a computer with parallel abilities (Holmes et al., 2011).

Within the standard IMH algorithm of Section 2.1, a cost-free improvement can be obtained by a straightforward Rao–Blackwellization argument. Given that at step $t+i$, $y_{i}$ is accepted with probability $\rho(x_{t+i-1},y_{i})$ and rejected with probability $1-\rho(x_{t+i-1},y_{i})$, the weight of $y_{i}$ can be updated by $\rho(x_{t+i-1},y_{i})$ and the weight of the simulated value $y_{j}$ corresponding to $x_{t+i-1}$ can be similarly updated by the probability $1-\rho(x_{t+i-1},y_{i})$. Considering next the block IMH algorithm, at the beginning of each block we can define $p$ weights, denoted by $(w_{k})_{k=1}^{p}$, initialized at $0$ and then, for the first of the $p$ parallel chains, denoting by $j$ the index such that $x^{(1)}_{t+i-1}=y_{j}$, we update these weights at each time $t+i$ as

This is obviously repeated for each of the other parallel chains, ending up with $\sum_{k}w_{k}=p^{2}$. This leads to a new estimator

This estimator still depends on all uniform generations created within the block, since those weights $w_{k}$ depend upon the acceptances and rejections of the $y_{k}$’s made during the block update. However, along the steps of the block, the $w_{k}$’s are basically updated by the expectations of the acceptance indicators conditionally upon the results of the previous iterations, whereas the $n_{k}$ of Section 2 are directly updated according to the acceptance indicators. Hence, the $w_{k}$’s have a smaller variance than the $n_{k}$’s by virtue of the Rao–Blackwell theorem, leading to $\hat{\tau}_{3}$ necessarily having a smaller variance than $\hat{\tau}_{2}$.

We now discuss a more involved Rao-Blackwellization technique first proposed by Casella and Robert (1996).

Exploiting the Rao–Blackwellization technique of Casella and Robert (1996) within each parallel chain provides via a conditioning argument an even more stable approximation of arbitrary posterior quantities. As developed in Casella and Robert (1996), for a single Markov chain $(x_{1}^{(i)},\ldots,x_{p}^{(i)})$, a Rao–Blackwell weighting scheme on the proposed values $y_{t}$, with weights $\varphi_{t}$, is given by a recursive scheme

where $(t>0)$

and

associated with the Metropolis–Hastings ratios

The cumulated computation of the $\delta_{t}$’s, of the $\rho_{tu}$’s and of the $\xi_{tu}$’s requires an $\text{O}(p^{2})$ computing time. Given that $p$ is usually not very large, this additional cost is not redhibitory as in the original proposal of Casella and Robert (1996) who were considering the application of this Rao–Blackwellization technique over the whole chain, with a cost of $\text{O}(T^{2})$ (see also Perron, 1999).

Therefore, starting from the estimator $\hat{\tau}_{2}$, the weight $n_{k}$ counting the number of occurrences of $y_{k}$ in the $p\times p$ block can be replaced with the expected number $\varphi_{k}$ of times $y_{k}$ occurs in this block (given the $p$ proposed values), which is the sum of the expected numbers of times $y_{k}$ occurs in each of the $p$ parallel chain:

Since the $p$ parallel chains incorporate the proposed values with different orders, the $\varphi$’s may differ for each chain and must therefore be computed $p$ times. Note that the cost is still in $\text{O}(p^{2})$ if this computation can be implemented in parallel. Then, by a Rao-Blackwell argument, $\hat{\tau}_{2}$ and $\tau_{3}$ are dominated by $\hat{\tau}_{4}$ defined as follows:

Therefore, this Rao–Blackwellization scheme involves no uniform generation for the computation of $\hat{\tau}_{4}$: the randomness associated with these uniforms is completely integrated out.

The four estimators defined up to now can be summarized as follows:

• •

  $\hat{\tau}_{1}$ is the basic IMH estimator of $\mathbb{E}_{\pi}\left[h(X)\right]$,

• •

  $\hat{\tau}_{2}$ improves $\hat{\tau}_{1}$ by averaging over permutations of the proposed values, and by using $p$ times more uniforms than $\hat{\tau}_{1}$,

• •

  $\hat{\tau}_{3}$ improves upon $\hat{\tau}_{2}$ by a basic Rao-Blackwell argument,

• •

  $\hat{\tau}_{4}$ improves upon the above by a further Rao-Blackwell argument, integrating out the ancillary uniform variables, but at a cost of $\text{O}(p^{2})$.

Note that these four estimators all involve the same number $p$ of target density evaluations, which again represent the overwhelming part of the computing time.

Figure 5 gives a comparison between the variances of the three improved estimators defined above and the variance of the basic IMH estimator. The permutations are random in this case. As was already apparent on Figure 4, the block estimator $\hat{\tau}_{2}$ is significantly better than $\hat{\tau}_{1}$ for any value of $p$. Moreover, both Rao-Blackwellization modifications seem to improve only very slightly the estimation when compared with $\hat{\tau}_{2}$, even though the improvement increases with $p$.

Recall that the “same order” scheme already provided a significant decrease in the variance of the estimation. In the light of our results, our interpretation is that using $p$ parallel chains with the same proposed values acts like a ”poor man” Rao–Blackwellization technique since $p$ times more uniforms are used. Specifically, each of the $p$ proposed values is proposed $p$ times instead of once, thus reducing the impact of each single uniform draw on the overall estimation.

When we use Rao–Blackwellization on top of the block IMH, in the estimators $\hat{\tau}_{3}$ and $\hat{\tau}_{4}$, we try indeed to integrate out a randomness that already is partly gone. This explains why, although Rao–Blackwellization techniques provide a significant improvement over standard IMH, the improvement is much lower and thus rather unappealing when used in the block IMH setting. This outcome was at first frustrating since Rao–Blackwellization is indeed affordable at a cost of only $\text{O}(p^{2})$. However, this shows in fine that the improvement brought by the block IMH algorithm roughly provides the same improvement as the Rao–Blackwell solution, at a much lower cost.

The proposal density $\mu$ may also be used to construct directly an importance sampling (IS) estimator

where the values $y_{t}$ are drawn from $\mu$. It therefore makes sense to compare the IMH algorithm with an IS approximation because IS is similarly easy to parallelize, and straightforward to program. Furthermore, since the IS estimator does not involve ancillary uniform variables, it is comparable to the Rao–Blackwellized version of IMH, and hence to the block IMH. Obviously, IS cannot necessarily be used in the settings when IMH algorithms are used, because the latter are also considered for approximating simulations from the target density $\pi$. In particular, when considering Metropolis-within-Gibbs algorithms, IS cannot be used in a straightforward manner, even for approximating integrals.

Before giving numerical results for a comparison run on the toy example, we now explain why in this comparison we took the number of blocks to be larger than $1$. The previous comparisons were computed with $b=1$, i.e. on a single $p\times p$ block and for a large number of independent runs. The choice of $b$ was then irrelevant since we were comparing methods that were exploiting in different ways the $p$ proposed values generated in each block. When considering the block IMH algorithm as a whole, as explained in Section 2, the end of each block sees a new starting value chosen from the current block. This ensures that the algorithm produces a valid Markov chain. However, our construction also implies that the successive blocks produced by the algorithm are correlated, which should lead to lesser performances than for the IS estimator.

In the comparison between IMH and IS, we therefore need to take into account this correlation between successive blocks. To this effect, we produce the variance reductions for several values of $b$. Those reductions are presented in Figure 6 for $p=16$ and different values of $b=1,10,100$. Once again, the permutations in the block IMH algorithm are chosen to be random.

Figure 6 shows the a priori surprising result that, when selecting $b=1$ in the experiment, the variance results are in favor of the block IMH solutions over the IS estimator, but, for any realistic application, $b$ is (much) larger than $1$. For all $b\geq 10$, the IS estimator has a smaller variance than the three alternative block IMH estimators, if only by a small margin. (Note that the variance improvement over the original MCMC estimator is slightly increasing with $b$ despite the correlation between blocks, given that the correlation between the $p^{2}$ terms involved in the block IMH estimators is lower than the correlation in the original MCMC chain.) This experiment thus shows that the block IMH solution gets very close to the IS estimator, while preserving the Markovian features of the original IMH algorithm.