A Review of Multiple Try MCMC algorithms for Signal Processing

Introduction

Bayesian methods have become very popular in signal processing over the last years . They require the application of sophisticated Monte Carlo techniques, such as Markov chain Monte Carlo (MCMC) and particle filters, for the efficient computation of a-posteriori estimators . More specifically, the MCMC algorithms generate a Markov chain such that its stationary distribution coincides with the posterior probability density function (pdf) . Typically, the only requirement is to be able to evaluate the target function, where the knowledge of the normalizing constant is usually not needed.

The most popular MCMC method is undoubtedly the Metropolis-Hastings (MH) algorithm . The MH technique is a very simple method, easy to be applied: this is the reason of its success. In MH, at each iteration, one new candidate is generated from a proposal pdf and then is properly compared with the previous state of the chain, in order to decide the next state. However, the performance of MH are often not satisfactory. For instance, when the posterior is multimodal, or when the dimension of the space increases, the correlation among the generated samples is usually high and, as a consequence, the variance of the resulting estimators grows. To speed up the convergence and reduce the “burn-in” period of the MH chain, several extensions have been proposed in literature.

In this work, we provide an exhaustive review of more sophisticated MCMC methods that, at each iteration, consider different candidates as possible new state of the chain. More specifically, at each iteration different samples are compared by certain weights and then one of them is selected as possible future state. The main advantage of these algorithms is that they foster the exploration of a larger portion of the sample space, decreasing the correlation among the states of the generated chain. In this work, we describe different algorithms of this family, independently introduced in literature. The main contribution is to present them under the same the framework and notation, remarking differences, relationships, limitations and strengths. All the discussed techniques yield an ergodic chain converging to the posterior density of interest (in the following, referred also as target pdf).

The first scheme of this MCMC class, called Orientational Bias Monte Carlo (OBMC) [12, Chapter 13], was proposed in the context of molecular simulation. Later on a more general algorithm, called Multiple Try Metropolis (MTM), was introduced . MTM includes OBMC as a special case (see Section 4.1.1). The MTM algorithm has been extensively studied and generalized in different ways . Other techniques, alternative to the MTM schemes, are the so-called the Ensemble MCMC (EnMCMC) methods . They follow a similar approach to MTM but employ a different acceptance function for selecting the next state of the chain. With respect to (w.r.t.) a generic MTM scheme, EnMCMC does not require any generation of auxiliary samples (as in a MTM scheme employing a generic proposal pdf) and hence, in this sense, EnMCMC are less costly.

In all the previous techniques, the candidates are drawn in a batch way and compared jointly. In the Delayed Rejection Metropolis (DRM) algorithm , in case of rejection of the novel possible state, the authors suggest to perform an additional acceptance test considering a new candidate. If this candidate is again rejected, the procedure can be iterated until reaching a desired number of attempts. The main benefit of DRM is that the proposal pdf can be improved at each intermediate stage. However, the acceptance function progressively becomes more complex so that the implementation of DRM for a great number of attempts is not straightforward (compared to the implementation of a MTM scheme with a generic number of tries).

In the last years, other Monte Carlo methods which combine particle filtering and MCMC have become very popular in the signal processing community. For instance, this is the case of the Particle Metropolis Hastings (PMH) and the Particle Marginal Metropolis Hastings (PMMH) algorithms, which have been widely used in signal processing in order to make inference and smoothing about dynamical and static parameters in state space models . PMH can be interpreted as a MTM scheme where the different candidates are generated and weighted by the use of a particle filter . In this work, we present PMH and PMMH and discuss their connections and differences with the classical MTM approach. Furthermore, we describe a suitable procedure for recycling some candidates in the final Monte Carlo estimators, called Group Metropolis Sampling (GMS) . The GMS scheme can be also seen as a way of generating a chain of sets of weighted samples. Finally, note that other similar and related techniques can be found within the so-called data augmentation approach .

The remaining of the paper is organized as follows. Section 2 recalls the problem statement and some background material, introducing also the required notation. The basis of MCMC and the Metropolis-Hastings (MH) algorithm are presented in Section 3. Section 4 is the core of the work, which describes the different MCMC using multiple candidates. Section 6 provides some numerical results, applying different techniques in a hyperparameter tuning problem for a Gaussian Process regression model, and in a localization problem considering a wireless sensor network. Some conclusions are given in Section 7.

Problem statement and preliminaries

The analytical study of the posterior density $\bar{\pi}({\bm{\theta}})$ is often unfeasible and integrals involving $\bar{\pi}({\bm{\theta}})$ are typically intractable . For instance, one might be interested in the estimation of

where $f({\bm{\theta}})$ is a generic integrable function w.r.t. $\bar{\pi}$. Dynamic and static parameters. In some specific application, the variable of interest ${\bm{\theta}}$ can be split in two disjoint parts, ${\bm{\theta}}=[{\bf x},{\bm{\lambda}}]$, where one, ${\bf x}$, is involved into a dynamical system (for instance, $\bf x$ is the hidden state in a state-space model) and the other, ${\bm{\lambda}}$, is a static parameter (for instance, an unknown parameter of the model). The strategies for making inference about ${\bf x}$ and ${\bm{\lambda}}$ should take into account the different nature of the two parameters (e.g., see Section 4.2.2). The main notation and acronyms are summarized in Tables 1-2.

${\widehat{I}}_{T}$ converges in probability to $I$ due to the weak law of large numbers. The approximation above ${\widehat{I}}_{T}$ is known as a direct (or ideal) Monte Carlo estimator if the samples ${\bm{\theta}}_{t}$ are independent and identically distributed (i.i.d.) from $\bar{\pi}$. Unfortunately, in many practical applications, direct methods for drawing independent samples from $\bar{\pi}({\bm{\theta}})$ are not available. Therefore, different approaches are required, such as the Markov chain Monte Carlo (MCMC) techniques.

Markov chain Monte Carlo (MCMC) methods

A MCMC algorithm generates an ergodic Markov chain with invariant (a.k.a., stationary) density given by the posterior pdf $\bar{\pi}({\bm{\theta}})$ . Specifically, given a starting state ${\bm{\theta}}_{0}$, a sequence of correlated samples is generated, ${\bm{\theta}}_{0}\rightarrow{\bm{\theta}}_{1}\rightarrow{\bm{\theta}}_{2}\rightarrow....\rightarrow{\bm{\theta}}_{T}$. Even if the samples are now correlated, the estimator

is consistent, regardless the starting vector ${\bm{\theta}}^{(0)}$ . Recall we are assuming that the Markov chain is ergodic and hence the starting value is forgotten. With respect to the direct Monte Carlo approach using i.i.d. samples, the application of an MCMC algorithm entails a loss of efficiency of the estimator ${\widetilde{I}}_{T}$, since the samples are positively correlated, in general. In other words, to achieve a given variance obtained with the direct Monte Carlo estimator, it is necessary to generate more samples. Thus, in order to improve the performance of an MCMC technique we have to decrease the correlation among the states of the chain. For the sake of simplicity, we use all the generated states in the final estimators, without removing any burn-in period .

The algorithm returns the sequence of states $\{{\bm{\theta}}_{1},{\bm{\theta}}_{2},\ldots,{\bm{\theta}}_{t},\ldots,{\bm{\theta}}_{T}\}$ (or a subset of them removing the burn-in period if an estimation of its length is available). We can see that the next state ${\bm{\theta}}_{t}$ can be the proposed sample ${\bm{\theta}}^{\prime}$ (with probability $\alpha$) or the previous state ${\bm{\theta}}_{t-1}$ (with probability $1-\alpha$). Under some mild regularity conditions, when $t$ grows, the pdf of the current state ${\bm{\theta}}_{t}$ converges to the target density ${\bar{\pi}}({\bm{\theta}})$ . The MH algorithm satisfies the so-called detailed balance condition which is sufficient to guarantee that the output chain is ergodic and has ${\bar{\pi}}$ as stationary distribution . Note that the acceptance probability $\alpha$ can be rewritten as

where we have denoted $w({\bm{\theta}}^{\prime}|{\bm{\theta}}_{t-1})=\frac{\pi({\bm{\theta}}^{\prime})}{q({\bm{\theta}}^{\prime}|{\bm{\theta}}_{t-1})}$ and $w({\bm{\theta}}_{t-1}|{\bm{\theta}}^{\prime})=\frac{\pi({\bm{\theta}}_{t-1})}{q({\bm{\theta}}_{t-1}|{\bm{\theta}}^{\prime})}$ in a similar fashion of the importance sampling weights of ${\bm{\theta}}^{\prime}$ and ${\bm{\theta}}_{t-1}$ . If the proposal pdf is independent from the previous state, i.e., $q({\bm{\theta}}|{\bm{\theta}}_{t-1})=q({\bm{\theta}})$, the acceptance function depends on the ratio of the importance weights $w({\bm{\theta}}^{\prime})=\frac{\pi({\bm{\theta}}^{\prime})}{q({\bm{\theta}}^{\prime})}$ and $w({\bm{\theta}}_{t-1})=\frac{\pi({\bm{\theta}}_{t-1})}{q({\bm{\theta}}_{t-1})}$, as shown in Table 4. We refer to this special MH case as the Independent MH (I-MH) algorithm. It is strictly related to other techniques described in the following (e.g., see Section 4.2.1).

MCMC using multiple candidates

In the standard MH technique described above, at each iteration one new sample ${\bm{\theta}}^{\prime}$ is generated to be tested with the previous state ${\bm{\theta}}_{t-1}$ by the acceptance probability $\alpha({\bm{\theta}}_{t-1},{\bm{\theta}}^{\prime})$. Other generalized MH schemes generate several candidates at each iteration to be tested as new possible state. In all these schemes, an extended acceptance probability $\alpha$ is properly designed in order to guarantee the ergodicity of the chain. Figure 1 provides a graphical representation of the difference between MH and the techniques using several candidates.

Below, we describe the most important examples of this class of MCMC algorithms. In most of them, a single MH-type test is performed at each iteration whereas in other methods a sequence of tests is employed. Furthermore, most of these techniques use an Importance Sampling (IS) approximation of the target density in order to improve the proposal procedure employed within a MH-type algorithm. Namely, they build an IS approximation, and then draw one sample from this approximation (resampling step). Finally, the selected sample is compared with the previous state of the chain, ${\bm{\theta}}_{t-1}$, according to a suitable generalized acceptance probability $\alpha$. It can be proved that all the methodologies presented in this work yield a ergodic chain with the posterior ${\bar{\pi}}$ as invariant density.

The Multiple Try Metropolis (MTM) algorithms are examples of this class of methods, where $N$ samples ${\bm{\theta}}^{(1)},{\bm{\theta}}^{(2)},\ldots,{\bm{\theta}}^{(N)}$ (called also “tries” or “candidates”) are drawn from the proposal pdf $q({\bm{\theta}})$, at each iteration . Then, one of them is selected according to some suitable weights. Finally, the selected candidate is accepted or rejected as new state according to a generalized probability function $\alpha$.

The MTM algorithm is given in Table 5. For the sake of simplicity, we have considered the use of the importance weights $w({\bm{\theta}}|{\bm{\theta}}_{t-1})=\frac{\pi({\bm{\theta}})}{q({\bm{\theta}}|{\bm{\theta}}_{t-1})}$, but there is not a unique possibility, as also shown below . In its general form, when the proposal depends on the previous state of the chain $q({\bm{\theta}}|{\bm{\theta}}_{t-1})$, the MTM requires the generation of $N-1$ auxiliary samples, ${\bf v}^{(i)}$, which are employed in the computation of the acceptance function $\alpha$. They are needed in order to guarantee the ergodicity. Indeed, the resulting MTM kernel satisfies the detailed balance condition, so that the chain is reversible . Note that for $N=1$, we have ${\bm{\theta}}^{(j)}={\bm{\theta}}^{(1)}$, ${\bf v}^{(1)}={\bm{\theta}}_{t-1}$ and the acceptance probability of the MTM method becomes

that is the acceptance probability of the classical MH technique. Several variants have been studied, for instance, with correlated tries and considering the use of different proposal pdfs .

The MTM method in Table 5 needs at step 2d the generation of $N-1$ auxiliary samples and at step 2e the computation of their weights (and, as a consequence, $N-1$ additional evaluation of the target pdf are required), that are only employed in the computation of the acceptance function $\alpha$.

The importance weights are not the unique possible choice. It is possible to show that the MTM algorithm generates an ergodic chain with invariant density $\bar{\pi}$, if the weight function $w({\bm{\theta}}|{\bm{\theta}}_{t-1})$ is chosen with the form

For instance, choosing $\xi({\bm{\theta}}_{t-1},{\bm{\theta}})=\frac{1}{q({\bm{\theta}}|{\bm{\theta}}_{t-1})q({\bm{\theta}}_{t-1}|{\bm{\theta}})}$, we obtain the importance weights $w({\bm{\theta}}|{\bm{\theta}}_{t-1})=\frac{\pi({\bm{\theta}})}{q({\bm{\theta}}|{\bm{\theta}}_{t-1})}$ used above. If we set $\xi({\bm{\theta}}_{t-1},{\bm{\theta}})=1$, we have $w({\bm{\theta}}|{\bm{\theta}}_{t-1})=\pi({\bm{\theta}})q({\bm{\theta}}_{t-1}|{\bm{\theta}})$. Another interesting example can be employed if the proposal is symmetric, i.e., $q({\bm{\theta}}|{\bm{\theta}}_{t-1})=q({\bm{\theta}}_{t-1}|{\bm{\theta}})$. In this case, we can choose $\xi({\bm{\theta}}_{t-1},{\bm{\theta}})=\frac{1}{q({\bm{\theta}}_{t-1}|{\bm{\theta}})}$ and then $w({\bm{\theta}}|{\bm{\theta}}_{t-1})=w({\bm{\theta}})=\pi({\bm{\theta}})$, i.e., the weights only depend on the value of the target density at ${\bm{\theta}}$. Thus, MTM contains the Orientational Bias Monte Carlo (OBMC) scheme [12, Chapter 13] as a special case, when a symmetric proposal pdf is employed, and then one candidate is chosen with weights proportional to the target density, i.e., $w({\bm{\theta}}|{\bm{\theta}}_{t-1})=\pi({\bm{\theta}})$.

1.2 Independent Multiple Try Metropolis (I-MTM) schemes

The MTM method described in Table 5 requires to draw $2N-1$ samples at each iteration ($N$ candidates and $N-1$ auxiliary samples) and $N-1$ are only used in the acceptance probability function. The generation of the auxiliary points

can be avoided if the proposal pdf is independent from the previous state, i.e., $q({\bm{\theta}}|{\bm{\theta}}_{t-1})=q({\bm{\theta}})$. Indeed, in this case, we should draw $N-1$ samples again from $q({\bm{\theta}})$ at the step 2d of Table 5. Since we have already drawn $N$ samples from $q({\bm{\theta}})$ at step 2a of Table 5, we can set

without jeopardizing the ergodicity of the chain (recall that ${\bf v}^{(j)}={\bm{\theta}}_{t-1}$). Hence, we can avoid step 2d and the acceptance function can be rewritten as

The I-MTM algorithm is provided in Table 6.

An I-MTM method requires only $N$ new evaluations of the target pdf at each iteration, instead of $2N-1$ new evaluations in the generic MTM scheme in Table 5. Note that we can also write $\alpha({\bm{\theta}}_{t-1},{\bm{\theta}}^{(j)})$ as

Alternative version (I-MTM2). From the IS theory, we know that $\widehat{Z}_{1}=\frac{1}{N}\sum_{n=1}^{N}w({\bm{\theta}}^{(n)})$ is an unbiased estimator of the normalizing constant $Z$ of the target $\pi$ (a.k.a, Bayesian evidence or marginal likelihood). It suggests to replace $\widehat{Z}_{2}$ with other unbiased estimators of $Z$ (without jeopardizing the ergodicity of the chain). For instance, instead of recycling the samples generated in the same iteration as auxiliary points as in Eq. (14), we could reuse samples generated in the previous iteration $t-1$. This alternative version of I-MTM method (I-MTM2) is given in Table 7. Note that, in both cases I-MTM and I-MTM2, the selected candidate ${\bm{\theta}}^{(j)}$ is drawn from the following particle approximation of the target $\bar{\pi}$,

i.e., ${\bm{\theta}}^{(j)}\sim\widehat{\pi}({\bm{\theta}}|{\bm{\theta}}^{(1:N)})$. The acceptance probability $\alpha$ used in I-MTM2 can be also justified considering a proper IS weighting of a resampled particle and using the expression (7) related to the standard MH method, as discussed in . Figure 2 provides a graphical representation of the I-MTM schemes.

1.3 Reusing candidates in parallel I-MTM chains

Let us consider to run $C$ independent parallel chains yielded by an I-MTM scheme. In this case, we have $NC$ evaluations of the target function $\pi$ and $C$ resampling steps performed at each iteration (so that we have $NCT$ total target evaluations and $CT$ total resampling steps).

In literature, different authors have suggested to recycle the $N$ candidates, ${\bm{\theta}}^{(1)},{\bm{\theta}}^{(2)},\ldots,{\bm{\theta}}^{(N)}\sim q({\bm{\theta}})$, in order to reduce the number of evaluations of the target pdf . The idea is to performs $C$-times the resampling procedure considering the same set of candidates, ${\bm{\theta}}^{(1)},{\bm{\theta}}^{(2)},\ldots,{\bm{\theta}}^{(N)}$ (a similar approach was proposed in ). Each resampled candidate is then tested as possible future state of one chain. In this scenario, The number of target evaluations per iteration is only $N$ (hence, the total number of evaluation of $\pi$ is $NT$). However, the resulting $C$ parallel chains are no longer independent, and there is a lose of performance w.r.t. the independent chains. There exists also the possibility of reducing the total number of resampling steps, as suggested in the Block Independent MTM scheme (but the dependence among the chains grows even more).

2 Particle Metropolis-Hastings (PMH) method

Assume that the variable of interest is formed by only a dynamical variable, i.e., ${\bm{\theta}}={\bf x}=x_{1:D}=[x_{1}\dots,x_{D}]^{\top}$ (see Section 2). This is the case of inferring a hidden state in state-space model, for instance. More generally, let assume that we are able to factorize the target density as

The Particle Metropolis Hastings (PMH) method is an efficient MCMC technique, proposed independently from the MTM algorithm, specifically designed for being applied in this framework. Indeed, we can take advantage of the factorization of the target pdf and consider a proposal pdf decomposed in the same fashion

Then, as in a batch IS scheme, given an $n$-th sample ${\bf x}^{(n)}=x_{1:D}^{(n)}\sim q({\bf x})$ with $x_{d}^{(n)}\sim q_{d}(x_{d}|x_{d-1})$, we assign the importance weight

The previous expression suggests a recursive procedure for computing the importance weights: starting with $w_{1}^{(n)}=\frac{\pi(x_{1}^{(n)})}{q(x_{1}^{(n)})}$ and then

for $d=2,\ldots,D$. This method is usually referred as Sequential Importance Sampling (SIS). If resampling steps are also employed at some iteration, the method is called Sequential Importance Resampling (SIR), a.k.a., particle filtering (PF) (see Appendix B). PMH uses a SIR approach for providing the particle approximation $\widehat{\pi}({\bf x}|{\bf x}^{(1:N)})=\sum_{i=1}^{N}\bar{w}_{D}^{(i)}\delta({\bf x}-{\bf x}^{(i)})$ where $\bar{w}_{D}^{(i)}=\frac{w_{D}^{(i)}}{\sum_{n=1}^{N}w_{D}^{(n)}}$ and $w_{D}^{(i)}=w({\bf x}^{(i)})$, obtained using Eq. (26) (with a proper weighting of a resampled particle ). Then, one particle is drawn from this approximation, i.e., with a probability proportional to the corresponding normalized weight. Estimation of the marginal likelihood $Z$ in particle filtering. SIR combines the SIS approach with the application of resampling procedures. In SIR, a consistent estimator of $Z$ is given by

Due to the application of the resampling, in SIR the standard estimator

is a possible alternative only if a proper weighting of the resampled particles is applied (otherwise, it is not an estimator of $Z$). If a proper weighting of a resampled particle is employed, both $\widetilde{Z}$ and $\widehat{Z}$ are equivalent estimators of $Z$ . Without the use of resampling steps (i.e., in SIS), $\widetilde{Z}$ and $\widehat{Z}$ are always equivalent estimators . See also Appendix B. The complete description of PMH is provided in Table 8 considering the use of $\widetilde{Z}$. At each iteration, a particle filter is run in order to provide an approximation by $N$ weighted samples of the measure of the target. Then, a sample among the $N$ weighted particles is chosen by one resampling step. This selected sample is then accepted or rejected as next state of the chain according to an MH-type acceptance probability, which involves two estimators of marginal likelihood $Z$. PMH is also related to other popular method in molecular simulation called Configurational Bias Monte Carlo (CBMC) .

A simple look at I-MTM2 and PMH shows that they are strictly related . Indeed, the structure of the two algorithms coincides. The main difference lies that the candidates in PMH are generated sequentially, using a SIR scheme. If no resampling steps are applied, then I-MTM2 and PMH are exactly the same algorithm, where the candidates are drawn in a batch setting or sequential way. Hence, the application of resampling steps is the main difference between the generation procedures of PMH and I-MTM2. Owing to the use of resampling, the candidates $\{{\bf x}^{(1)},\ldots,{\bf x}^{(N)}\}$ proposed by PMH are not independent (differently from I-MTM2). As an example, Figure 3 shows $N=40$ particles (with $D=10$) generated and weighted by SIS and SIR procedures (each path is a generated particle ${\bf x}^{(i)}=x_{1:10}^{(i)}$). The generation of correlated samples can be also considered in MTM methods without jeopardizing the ergodicity of the chain, as simply shown for instance in , for instance. Another difference is the use of ${\widetilde{Z}}$ or ${\widehat{Z}}$. However, if a proper weighting of a resampled particle is employed, both estimators coincide . Furthermore, both I-MTM2 and PMH can be considered as I-MH schemes where a proper importance sampling weighting of a resampled particle is employed . Namely, I-MTM2 and PMH are equivalent to an I-MH technique using the following complete proposal pdf,

where $\widetilde{\pi}$ is given in Eq. (20), i.e., ${\bm{\theta}}^{(j)}\sim\widetilde{q}({\bm{\theta}})$, and then considering the generalized (proper) IS weighting, $w({\bm{\theta}}^{(j)})=\widehat{Z}^{*}$, $w({\bm{\theta}}_{t-1})=\widehat{Z}_{t-1}$ . For further details see Appendix A.

2.2 Particle Marginal Metropolis-Hastings (PMMH) method

Assume now that the variable of interest if formed by both dynamical and static variables, i.e., ${\bm{\theta}}=[{\bf x},{\bm{\lambda}}]^{\top}$. For instance, this is the case of inferring both, an hidden state ${\bf x}$ in state-space model, and static parameters ${\bm{\lambda}}$ of the model. The Particle Marginal Metropolis-Hastings (PMMH) technique is a extension of PMH which addresses this problem.

where $\pi({\bf x}|{\bm{\lambda}})=\gamma_{1}(x_{1}|{\bm{\lambda}})\prod_{d=2}^{D}\gamma_{d}(x_{d}|x_{1:d-1},{\bm{\lambda}})$. For a specific value of ${\bm{\lambda}}$, we can use a particle filter approach, obtaining the approximation $\widehat{\pi}({\mathbf{x}}|{\bm{\lambda}})=\sum_{n=1}^{N}{\bar{w}}_{D}^{(n)}\delta({\mathbf{x}}-{\mathbf{x}}^{(n)})$ and the estimator $\widetilde{Z}({\bm{\lambda}})$, as described above. The PMMH technique is then summarized in Table 9. The pdf $q_{\lambda}({\bm{\lambda}}|{\bm{\lambda}}_{t-1})$ denotes the proposal density for generating possible values of ${\bm{\lambda}}$. Observe that, with the specific choice $q_{\lambda}({\bm{\lambda}}|{\bm{\lambda}}_{t-1})=g_{\lambda}({\bm{\lambda}})$, then the acceptance function becomes

Note also that PMMH w.r.t. to ${\bm{\lambda}}$ can be interpreted as MH method where an the posterior cannot be evaluated point-wise. Indeed, $\widetilde{Z}({\bm{\lambda}})$ approximates the marginal likelihood $p({\bf y}|{\bm{\lambda}})$ .

3 Group Metropolis Sampling

The auxiliary weighted samples in the I-MTM schemes (i.e., the $N-1$ samples drawn at each iteration that are not selected to be compared with the previous state ${\bm{\theta}}_{t-1}$) can be recycled providing a consistent and more efficient estimators .

The so-called Group Metropolis Sampling (GMS) method is shown in Table 10. GMS yields a sequence of sets of weighted samples $\mathcal{S}_{t}=\{{\bm{\theta}}_{n,t},\rho_{n,t}\}_{n=1}^{N}$, for $t=1,\ldots,T$, where we have denoted with $\rho_{n,t}$ the importance weights assigned to the samples ${\bm{\theta}}_{n,t}$ (see Figure 4). All the samples are then employed for a joint particle approximation of the target. Alternatively, GMS can directly provide an approximation of a specific moment of the target pdf (i.e., given a particular function $f$). The estimator of this specific moment provided by GMS is

Unlike in the I-MTM schemes, no resampling steps are performed in GMS. However, we can recover an I-MTM chain from the GMS output applying one resampling step when $\mathcal{S}_{t}\neq\mathcal{S}_{t-1}$, i.e.,

for $t=1,\ldots,T$. More specifically, $\{{\bm{\theta}}_{t}\}_{t=1}^{T}$ is a Markov chain obtained by one run of an I-MTM2 technique. The consistency of the GMS estimators is discussed in Appendix C. GMS can be also interpreted as an iterative IS scheme where an IS approximation of $N$ samples is built at each iteration and compared with the previous IS approximation. This procedure is iterated $T$ times and all the accepted IS estimators $\widetilde{I}_{N}^{(t)}$ are finally combined to provide a unique global approximation of $NT$ samples. Note that the temporal combination of the IS estimators is obtained dynamically by the random repetitions due to the rejections in the acceptance test.

The complete weighting procedure in GMS can be interpreted as the composition of two weighting schemes: (a) by an IS approach building $\{\rho_{n,t}\}_{n=1}^{N}$ and (b) by the possible random repetitions due to the rejections in the acceptance test. Figure 4 depicts a graphical representation of the GMS outputs as chain of sets $\mathcal{S}_{t}=\{{\bm{\theta}}_{n,t},\rho_{n,t}\}_{n=1}^{N}$.

4 Ensemble MCMC algorithms

Another alternative procedure, often referred as Ensemble MCMC (EnMCMC) methods (a.k.a., called Locally weighted MCMC), involving several tries at each iteration . Related techniques has been proposed independently in different works . First, let us define the joint proposal density

and, considering $N+1$ possible elements, $\mathcal{S}=\{{\bm{\theta}}^{(1)},\ldots,{\bm{\theta}}^{(N)},{\bm{\theta}}^{(N+1)}\}$, we define the $D\times N$ matrix

with columns all the vectors in $\mathcal{S}$ with the exception of ${\bm{\theta}}^{(k)}$. For simplicity, in the followings we abuse of the notation writing $q({\bm{\theta}}^{(1)},\ldots,{\bm{\theta}}^{(N)}|{\bm{\theta}}_{t})=q({\bf\Theta}_{\neg N+1}|{\bm{\theta}}_{t})$, for instance. One simple example of joint proposal pdf is

i.e., considering independence among ${\bm{\theta}}^{(n)}$’s (and having the same marginal proposal pdf $q$). More sophisticated joint proposal densities can be employed. A generic EnMCMC algorithm is outlined in Table 11.

Note that with respect to the generic MTM method, EnMCMC does not require to draw auxiliary samples and weights. Therefore, in EnMCMC a smaller number of evaluation of target is required w.r.t. a generic MTM scheme.

In this section, we present an interesting special case, which employs a single proposal pdf $q({\bm{\theta}})$ independent on the previous state of the chain, i.e.,

In this case, the technique can be simplified as shown below. At each iteration, the algorithm described in Table 12 generates $N$ new samples ${\bm{\theta}}^{(1)},{\bm{\theta}}^{(2)},\ldots,{\bm{\theta}}^{(N)}$ and then resample the new state ${\bm{\theta}}_{t}$ within a set of $N+1$ samples, $\{{\bm{\theta}}^{(1)},\ldots,{\bm{\theta}}^{(N)},{\bm{\theta}}_{t+1}\}$ (which includes the previous state), according to the probabilities

where $w({\bm{\theta}})=\frac{\pi({\bm{\theta}})}{q({\bm{\theta}})}$ denotes the importance sampling weight. Note that Eq. (43) for $N=1$ becomes

that is the Barker’s acceptance function (see ).

As discussed in [42, Appendix B], [48, Appendix C], , the density of a resampled candidate becomes closer and closer to $\bar{\pi}$ as $N$ grows, i.e., $N\rightarrow\infty$. Hence, the performance of I-EnMCMC clearly improves with $N\rightarrow\infty$ (see Appendix A). The I-EnMCMC algorithm produces an ergodic chain with invariant density $\bar{\pi}$, by resampling $N+1$ samples at each iteration ($N$ new samples ${\bm{\theta}}^{(1)},\ldots,{\bm{\theta}}^{(N)}$ from $q$ and setting ${\bm{\theta}}^{(N+1)}={\bm{\theta}}_{t-1}$). Figure 5 summarizes the steps of I-EnMCMC.

5 Delayed Rejection Metropolis (DRM) Sampling

An alternative use of different candidates in one iteration of a Metropolis-type method is given in . The idea behind the proposed algorithm, called Delayed Rejection Metropolis (DRM) algorithm, is the following. As in a standard MH method, at each iteration, one sample is proposed ${\bm{\theta}}^{(1)}\sim q_{1}({\bm{\theta}}|{\bm{\theta}}_{t-1})$ and accepted with probability

If ${\bm{\theta}}^{(1)}$ is accepted then ${\bm{\theta}}_{t}={\bm{\theta}}^{(1)}$ and the chain is moved forward. If ${\bm{\theta}}^{(1)}$ is rejected, the DRM method suggests of drawing another samples ${\bm{\theta}}^{(2)}\sim q_{2}({\bm{\theta}}|{\bm{\theta}}^{(1)},{\bm{\theta}}_{t-1})$ (considering a different proposal pdf $q_{2}$ taking into account possibly the previous candidate ${\bm{\theta}}^{(1)}$) and accepted with a suitable acceptance probability

The acceptance function $\alpha_{2}({\bm{\theta}}_{t-1},{\bm{\theta}}^{(2)})$ is designed in order to ensure the ergodicity of the chain. If ${\bm{\theta}}^{(2)}$ is rejected we can set ${\bm{\theta}}_{t}={\bm{\theta}}_{t-1}$ and perform another iteration of the algorithm, or continue with this iterative strategy drawing ${\bm{\theta}}^{(3)}\sim q_{3}({\bm{\theta}}|{\bm{\theta}}^{(2)},{\bm{\theta}}^{(1)},{\bm{\theta}}_{t-1})$ and test it with a proper probability $\alpha_{3}({\bm{\theta}}_{t-1},{\bm{\theta}}^{(3)})$. The DRM algorithm with only 2 acceptance stages is outlined in Table 13 and summarized in Figure 6.

Note that the proposal pdf can be improved at each intermediate stage (${\bm{\theta}}^{(1)}\sim q_{1}({\bm{\theta}}|{\bm{\theta}}_{t-1})$, ${\bm{\theta}}^{(2)}\sim q_{2}({\bm{\theta}}|{\bm{\theta}}^{(1)},{\bm{\theta}}_{t-1})$ etc.), using the information provided by the previous generated samples and the corresponding target evaluations.

The idea behind DRM of creating a path of intermediate points, then improving the proposal pdf, and hence fostering larger jumps have been also considered in other works .

Summary: computational cost, differences and connections

The performance of the algorithms described above improves as $N$ grows, in general: the correlation among samples vanishes to zero, and the acceptance rate of new state approaches one (see Section 6.1). Generally, an acceptance rate close to $1$ is not an evidence of good performance for an MCMC algorithm. However, for the techniques tackled in this work, the situation is different: as $N$ grows, the procedure used for proposing a novel possible state (involving $N$ tries, resampling steps etc.) becomes better and better, yielding a better approximation of the target pdf. See Appendix A for further details. As $N$ increases, they become similar and similar to an exact sampler drawing independent samples directly from the target density (for MTM, PMH and EnMCMC schemes the explanation is given in Appendix A). However, this occurs at the expense of an additional computational cost.

In Table 14, we summarize the total number of target evaluations, $E$, and the total number of samples used in the final estimators, $Q$ (without considering to remove any burn-in period). The generic MTM algorithm has the greatest number of target evaluations. However, a random-walk proposal pdf can be used in a generic MTM algorithm and, in general, it fosters the exploration of the state space. In this sense, the generic EnMCMC seems to be preferable w.r.t. MTM, since $E=NT$ and the random-walk proposal can be applied. A disadvantage of the EnMCMC schemes is that their acceptance function seems worse in terms of Peskun’s ordering (see numerical results in Section 6.1). Namely, fixing the number of $N$ tries, the target $\pi$ the proposal $q$ pdfs, the MTM schemes seem to provide greater acceptance rates than the corresponding EnMCMC techniques. This is theoretically proved for $N=1$ , and the difference vanishes to zero as $N$ grows. The GMS technique, like other strategies , has been proposed to recycle samples or re-use target evaluations, in order to increase $Q$ (see also Section 4.1.3).

In PMH, the components of the different tries are drawn sequentially and they are correlated due to the application of the resampling steps. In DMR, each candidate is drawn in a batch way (all the components jointly) but the different candidates are drawn in a sequential manner (see Figure 6), ${\bm{\theta}}^{(1)}$ then ${\bm{\theta}}^{(2)}$ etc. The benefit of this strategy is that the proposal pdf can be improved considering the previous generated tries. Hence, if the proposal takes into account the previous samples, DMR generates correlated candidates as well. The main disadvantage of DRM is that the implementation for a generic $N>2$ is not straightforward.

I-MTM and I-MTM2 differs for the acceptance function employed. Furthermore, The main difference between the I-MTM2 and PMH schemes is the use of resampling steps during the generation the different tries. For this reason, the candidates of PMH are correlated (unlike in I-MTM2). I-MTM2 and PMH can be interpreted as I-MH methods using a sophisticated proposal density $\widehat{q}({\bm{\theta}})$ in Eq. (31), and an extended IS weighting procedure is employed. Note that, indeed, $\widehat{q}$ cannot be evaluated pointwise, hence a standard IS weighting strategy cannot be employed.

Numerical Experiments

We test different MCMC using multiple candidates in different numerical experiments. In the first example, an exhaustive comparison among several techniques with an independent proposal is given. We have considered different number of tries, length of the chain, parameters of the proposal pdfs and also different dimension of the inference problem. In the second numerical simulation, we compare different particle methods. The third one regards the hyperparameter selection for a Gaussian Process (GP) regression model. The last two examples are localization problems in a wireless sensor network (WSN): in the fourth one some parameters of the WSN are also tuned, whereas in last example a real data analysis is performed.

In order to compare the performance of different techniques, in this section we consider a multi-modal, multidimensional Gaussian target density. More specifically, we have

where ${\bm{\mu}}_{1}=[\mu_{1,1},\ldots,\mu_{1,D}]^{\top}$, ${\bm{\mu}}_{2}=[\mu_{2,1},\ldots,\mu_{2,D}]^{\top}$, ${\bm{\mu}}_{3}=[\mu_{3,1},\ldots,\mu_{3,D}]^{\top}$, with $\mu_{1,d}=-3$, $\mu_{2,d}=0$, $\mu_{3,d}=2$ for all $d=1,\ldots,D$. Moreover, the covariance matrices are diagonal, ${\bm{\Sigma}}_{i}=\delta_{i}{\bf I}_{D}$ (where ${\bf I}_{D}$ is the $D\times D$ identity matrix), with $\delta_{i}=0.5$ for $i=1,2,3$. Hence, given a random variable ${\bf\Theta}\sim{\bar{\pi}}({\bm{\theta}})$, we know analytically that $E[{\bm{\Theta}}]=[{\bar{\theta}}_{1},\ldots,{\bar{\theta}}_{D}]^{\top}$ with ${\bar{\theta}}_{d}=-\frac{1}{3}$ for all $d$, and $\mbox{diag}\{\mbox{Cov}[{\bm{\Theta}}]\}=[\xi_{1},\ldots,\xi_{D}]^{\top}$ with $\xi_{d}=\frac{85}{18}$ for all $d=1,\ldots,D$.

We apply I-MTM, IMTM2 and I-EnMCMC in order to estimate all the expected values and all the variances of the marginal target pdfs. Namely, for a given dimension $D$, we have to estimate all $\{{\bar{\theta}}_{d}\}_{d=1}^{D}$ and $\{\xi_{d}\}_{d=1}^{D}$, hence $2D$ values. The results are averaged over $3000$ independent runs. At each run, we compute an averaged square error obtained in the estimation of the $2D$ values and then calculate the Mean Square Error (MSE) averaged over the $3000$ runs. For all the techniques, we consider a Gaussian proposal density $q({\bm{\theta}})=\mathcal{N}({\bm{\theta}}|{\bm{\mu}},\sigma^{2}{\bf I}_{D})$ with ${\bm{\mu}}=[\mu_{1}=0,\ldots,\mu_{D}=0]^{\top}$ (independent from the previous state) and different values of $\sigma$ are considered.

We perform several experiments varying the number of tries, $N$, the length of the generated chain, $T$, the dimension of the inference problem, $D$, and the scale parameter of the proposal pdf, $\sigma$. In Figures 7(a)-(b)-(c)-(d)-(e), we show the MSE (obtained by different techniques) as function of $N$, $T$, $D$ and $\sigma$, respectively. In Figure 7(d), we only consider I-MTM with $N\in\{1,5,100,1000\}$ in order to show the effect of using different tries $N$ in different dimensions $D$. Note that I-MTM with $N=1$ coincides with I-MH.

Let us denote as $\phi(\tau)$ the auto-correlation function of the states of the generated chain. Figures 8(a)-(b)-(c) depicts the normalized auto-correlation function ${\bar{\phi}}(\tau)=\frac{\phi(\tau)}{\phi(0)}$ (recall that $\phi(0)\geq\phi(\tau)$ for $\tau\geq 0$) at different lags $\tau=1,2,3$, respectively. Furthermore, given the definition of the Effective Sample Size (ESS) [51, Chapter 4],

in Figure 8(d), we show of the ratio $\frac{ESS}{T}$ (approximated; cutting off the series in the denominator at lag $\tau=10$), as function of $N$. Since the MCMC algorithms yield positive correlated sequences of states, we have $ESS<T$ in general. Finally, in Figures 9(a)-(b), we provide the Acceptance Rate (AR) of a new state (i.e., the expected number of accepted jumps to a novel state), as function of $N$ and $D$, respectively.

1.2 Comment on the results

Figures 7(a)-(b), 8 and 9(a), clearly show that the performance improves as $N$ grows, for all the algorithms. The MSE values and the correlation decrease, and the ESS and the AR grow. I-MTM seems to provide the best performance. Recall that for $N=1$, I-EnMCMC becomes an I-MH with Baker’s acceptance function and I-MTM becomes the I-MH in Table 4 . For $N=1$, the results confirm the Peskun’s ordering about the acceptance function for a MH method . Observing the results, the Peskun’s ordering appears valid also for the multiple try case, $N>1$. I-MTM2 seems to have worse performance than I-MTM for all $N$. With respect to I-EnMCMC, I-MTM2 performs better for smaller $N$. The difference among the MSE values obtained by the samplers becomes smaller as $N$ grows, as shown in Figure 7(a)-(b) (note that in the first one $D=1$, in the other $D=10$, and the range of $N$ is different). The comparison among I-MTM, I-MTM2 and I-EnMCMC seems not to be affected by changing $T$ and $\sigma$, as depicted in Figures 7(c)-(e). Namely, the MSE values change but the ordering of the methods (e.g., best and worst) seems to depend mainly on $N$. Obviously, for greater $D$, more tries are required in order to obtain good performance (see Figures 7(d) and9(b), for instance). Note that for $N\rightarrow\infty$, I-MTM, I-MTM2 and I-EnMCMC perform similarly to an exact sampler drawing $T$ independent samples from ${\bar{\pi}}({\bm{\theta}})$: the correlation $\phi(\tau)$ among the samples approaches zero (for all $\tau$), ESS approaches $T$ and AR approaches $1$.

2 Numerical experiment comparing particle schemes

In this section, in order to clarify the differences between batch and particle schemes, we consider again a multidimensional Gaussian target density, that can be express as

We apply I-MTM, I-MTM2, PMH, and a variant of PMH, denote as var-PMH which uses the corresponding acceptance probability of I-MTM in Eq. (19) instead of the acceptance function of the classical PMH in Eq. (30). The goal is to estimate the vector ${\bm{\mu}}$. We compute the MSE in estimating the vector ${\bm{\mu}}=\mu_{1:10}$, averaging the results over $500$ independent simulations. The components of ${\bm{\mu}}$ are shown in Figure 10(d) with a dashed line.

For all the techniques, we employ a sequential construction of the $N$ candidates (using the chain rule, see below): in PMH and var-PMH the resampling is applied at each iteration whereas in I-MTM and I-MTM2 no resampling is applied. More specifically, the proposal density for all the methods is

where $q_{1}(\theta_{1})=\mathcal{N}(\theta_{1}|-2,4)$ and $q_{d}(\theta_{d}|\theta_{d-1})=\mathcal{N}(\theta_{d}|\theta_{d-1},\sigma_{p}^{2})$, but PMH and var-PMH employ resampling steps so that the generated tries are correlated (whereas in I-MTM and I-MTM2 the generated candidates are independent).

We test all the techniques considering different value of number of tries $N$ and number of iterations of the chain $T$. Figures 10(a)-(b) show the MSE as function of number of iterations $T$, keeping fixed the number of tries $N=3$. Figure 10(a) reports the results of the MTM schemes whereas Figure 10(b) reports the results of the PMH schemes. Figure 10(c) depicts the MSE as function of $N$ (with $T=2000$), for the PMH methods. Note that the use of only $N=3$ particles and the application of the resampling at each iteration is clearly a disadvantage for the PMH schemes. If the resampling is applied very often (as in this case), a greater number of $N$ is advisable (such as $N=100$ or $N=1000$). Hence, the results confirm that applying a resampling step at each iteration is not optimal and that a smaller rate of resampling steps could improve the performance . The results also confirm that the use of an acceptance probability of type in Eq. (19) provides smaller MSE, i.e., I-MTM and var-PMH perform better than I-MTM2 and PMH, respectively. This is more evident for small number of candidates $N$. When $N$ grows, the performance of PMH and var-PMH methods becomes similar, since the acceptance probability approaches $1$, in both cases. Figure 10(d) depicts 35 different states ${\bm{\theta}}_{t}=\theta_{1:10,t}$ at different iteration indices $t$, obtained with var-PMH ($N=1000$ and $T=1000$) and the values $\mu_{1:10}$ are given in dashed line.

3 Hyperparameter tuning for Gaussian Process (GP) regression models

Given these assumptions, the vector ${\bf f}=[f({\bf z}_{1}),\ldots,f({\bf z}_{P})]^{\top}$ is distributed as $p({\bf f}|{\bf Z},\delta,\kappa)=\mathcal{N}({\bf f};{\bf 0},{\bf K})$, where ${\bf 0}$ is a $P\times 1$ null vector, and ${\bf K}_{ij}:=\kappa({\bf z}_{i},{\bf z}_{j})$, for all $i,j=1,\ldots,P$, is a $P\times P$ matrix. The vector containing all the hyperparameters of the model is ${\bm{\theta}}=[\delta,\sigma]$, i.e., all the parameters of the kernel function in Eq. (52) and standard deviation $\sigma$ of the observation noise. In this experiment, we focus on the marginal posterior density of the hyperparameters, ${\bar{\pi}}({\bm{\theta}}|{\bf y},{\bf Z},\kappa)\propto\pi({\bm{\theta}}|{\bf y},{\bf Z},\kappa)=p({\bf y}|{\bm{\theta}},{\bf Z},\kappa)p({\bm{\theta}})$, which can be evaluated analytically, but we cannot compute integrals involving it . Considering a uniform prior within $^{2}$, $p({\mathbf{x}})$ and since $p({\bf y}|{\bm{\theta}},{\bf Z},\kappa)=\mathcal{N}({\bf y};{\bf 0},{\bf K}+\sigma^{2}{\bf I})$, we have

where $C>0$, and clearly ${\bf K}$ depends on $\delta$ . The moments of this marginal posterior cannot be computed analytically. Then, in order to compute the Minimum Mean Square Error (MMSE) estimator $\widehat{{\bm{\theta}}}=[\widehat{\delta},\widehat{\sigma}]$, i.e., the expected value $E[{\bm{\Theta}}]$ with ${\bm{\Theta}}\sim{\bar{\pi}}({\bm{\theta}}|{\bf y},{\bf Z},\kappa)$, we approximate $E[{\bm{\Theta}}]$ via Monte Carlo quadrature. More specifically, we apply I-MTM2, GMS, a MH scheme with a longer chain and a static IS method. For all these methodologies, we consider the same number of target evaluations, denoted as $E$, in order to provide a fair comparison.

We generated $P=200$ pairs of data, $\{y_{j},{\bf z}_{j}\}_{j=1}^{P}$, according to the GP model above setting $\delta^{*}=3$, $\sigma^{*}=10$, $L=1$, and drawing $z_{j}\sim\mathcal{U}()$. We keep fixed these data over the different runs. We computed the ground-truth ${\bm{\theta}}=[\delta=3.5200,\sigma=9.2811]$ using an exhaustive and costly grid approximation, in order to compare the different techniques. For I-MTM2, GMS, and MH schemes, we consider the same adaptive Gaussian proposal pdf $q_{t}({\bm{\theta}}|{\bm{\mu}}_{t},\lambda^{2}{\bf I})=\mathcal{N}({\bm{\theta}}|{\bm{\mu}}_{t},\lambda^{2}{\bf I})$, with $\lambda=5$ and ${\bm{\mu}}_{t}$ is adapted considering the arithmetic mean of the outputs after a training period, $t\geq 0.2T$, in the same fashion of (${\bm{\mu}}_{0}=^{\top}$). First, we test both techniques fixing $T=20$ and varying the number of tries $N$. Then, we set $N=100$ and vary the number of iterations $T$. Figure 11 (log-log plot) shows the Mean Square Error (MSE) in the approximation of $\widehat{{\bm{\theta}}}$ averaged over $10^{3}$ independent runs. Observe that GMS always outperforms the corresponding I-MTM2 scheme. These results confirm the advantage of recycling the auxiliary samples drawn at each iteration during an I-MTM2 run. In Figure 12, we show the MSE obtained by GMS keeping invariant the number of target evaluations $E=NT=10^{3}$ and varying $N\in\{1,2,10,20,50,100,250,10^{3}\}$. As a consequence, we have $T\in\{10^{3},500,100,50,20,10,4,1\}$. Note that the case $N=1$, $T=10^{3}$, corresponds to an adaptive MH (A-MH) method with a longer chain, whereas the case $N=10^{3}$, $T=1$, corresponds to a static IS scheme (both with the same posterior evaluations $E=NT=10^{3}$). We observe that the GMS always provides smaller MSE than the static IS approach. Moreover, GMS outperforms A-MH with the exception of two cases where $T\in\{1,4\}$.

4 Localization of a target in a wireless sensor network

Our goal is to compute the Minimum Mean Square Error (MMSE) estimator, i.e., the expected value of the posterior ${\bar{\pi}}({\bm{\theta}}|\textbf{Y})={\bar{\pi}}({\bf z},{\bm{\zeta}}|\textbf{Y})$ (recall that $D=8$). Since the MMSE estimator cannot be computed analytically, we apply Monte Carlo methods for approximating it. We compare GMS, the corresponding MTM scheme, the Adaptive Multiple Importance Sampling (AMIS) technique , and $N$ parallel MH chains with a random walk proposal pdf. For all of them we consider Gaussian proposal densities. For GMS and MTM, we set $q_{t}({\bm{\theta}}|{\bm{\mu}}_{n,t},\sigma^{2}{\bf I})=\mathcal{N}({\bm{\theta}}|{\bm{\mu}}_{t},\sigma^{2}{\bf I})$ which is adapted considering the empirical mean of the generated samples after a training period, $t\geq 0.2T$ , ${\bm{\mu}}_{0}\sim\mathcal{U}(^{D})$ and $\sigma=1$. For AMIS, we have $q_{t}({\bm{\theta}}|{\bm{\mu}}_{t},{\bf C}_{t})=\mathcal{N}({\bm{\theta}}|{\bm{\mu}}_{t},{\bf C}_{t})$, where ${\bm{\mu}}_{t}$ is as previously described (with ${\bm{\mu}}_{0}\sim\mathcal{U}(^{D})$) and ${\bf C}_{t}$ is also adapted using the empirical covariance matrix, starting ${\bf C}_{0}=4{\bf I}$. We also test the use of $N$ parallel Metropolis-Hastings (MH) chains (we also consider the case of $N=1$, i.e., a single chain), with a Gaussian random-walk proposal pdf, $q_{n}({\bm{\mu}}_{n,t}|{\bm{\mu}}_{n,t-1},\sigma^{2}{\bf I})=\mathcal{N}({\bm{\mu}}_{n,t}|{\bm{\mu}}_{n,t-1},\sigma^{2}{\bf I})$ with ${\bm{\mu}}_{n,0}\sim\mathcal{U}(^{D})$ for all $n$ and $\sigma=1$.

We fix the total number of evaluations of the posterior density as $E=NT=10^{4}$. Note that, generally, the evaluation of the posterior is the most costly step in MC algorithms (however, AMIS has the additional cost of re-weighting all the samples at each iteration according to the deterministic mixture procedure ). We recall that $T$ denotes the total number of iterations and $N$ the number of samples drawn from each proposal at each iteration. We consider ${\bm{\theta}}^{*}=[{\bf z}^{*},{\bm{\zeta}}^{*}]^{\top}$ as the ground-truth and compute the Mean Square Error (MSE) in the estimation obtained with the different algorithms. The results are averaged over $500$ independent runs and they are provided in Tables 15, 16, and 17 and Figure 13(b). Note that GMS outperforms AMIS for each a pair $\{N,T\}$ (keeping fixed $E=NT=10^{4}$), and GMS also provides smaller MSE values than $N$ parallel MH chains (the case $N=1$ corresponds to a unique longer chain). Figure 13(b) shows the MSE versus $N$ maintaining $E=NT=10^{4}$ for GMS and the corresponding MTM method. This figure again confirms the advantage of recycling the samples in a MTM scheme.

5 Localization with real data

In this section, we describe a numerical experiment involving real data. More specifically, we consider a localization problem . We have carried out an experiment with a network consisting of four nodes. Three of them are placed at fixed positions and play the role of sensors that measure the strength of the radio signals transmitted by the target. The other node plays the role of the target to be localized. All nodes are bluetooth devices (Conceptronic CBT200U2A) with a nominal maximum range of 200 m. We consider a square monitored area of $4\times 4$ m and place the sensors at fixed positions $\textbf{h}_{1}=[0.5,1]$, $\textbf{h}_{2}=[3.5,1]$ and $\textbf{h}_{3}=$, with all coordinates in meters. The target is located at $\textbf{z}=[z_{1}=2.5,z_{2}=2]$. The measurement provided by the $i$-th sensor is denoted as a random variable $Y_{i}$, considering the following model

where $B_{i}$ are again independent Gaussian random variables with pdfs $\mathcal{N}(b_{i};0,\zeta^{2})$, for all $i=1,2,3$. Differently from the previous section, we estimate in advance the following parameters of the model, ${\hat{\kappa}}\approx-26.58$ and ${\hat{\zeta}}\approx 4.73$, using a least square fitting. We obtain $N_{O}=5$ measurements from each sensor ($d_{Y}=3N_{O}=15$), and we consider a uniform prior on the $4\times 4$ m area. Given these measurements, we approximate the expected value $E[{\bf Z}]$ of the corresponding posterior ${\bar{\pi}}({\bf z})$ (here ${\bm{\theta}}={\bf z}$) using a thin deterministic bivariate grid, obtaining the ground truth $\approx[3.17,2.62]^{\top}$. We test an MH method and MTM scheme using with a random walk Gaussian proposal pdf, $q({\bf z}|{\b{z}}_{t-1})=\mathcal{N}({\bf z}|{\bf z}_{t-1},\sigma^{2}{\bf I}_{2})$, with $\sigma=1$, $T\in\{1000,5000\}$, and $N\in\{10,100,1000\}$ tries for MTM (clearly, $N=1$ for MH). We also test a Metropolis-adjusted Langevin algorithm (MALA), where the proposal is Gaussian random walk density with mean ${\bf z}_{t-1}+\beta\nabla[\log\pi({\bf z})]$ and $\nabla[\log\pi({\bf z})]$ denotes the gradient of $\log\pi({\bf z})$ . The covariance matrix of MALA Gaussian proposal is $\sigma^{2}{\bf I}_{2}$ (as the other techniques) and the drift parameter $\beta=\sigma^{2}/2$. We compute the MSE in estimating $E[{\bf Z}]\approx[3.17,2.62]^{\top}$ and averaged the results over $2000$ independent runs (at each run, we take the mean of the square error values of each component). The results are shown in Table 18. Recall that MALA uses the additional information of the gradient. Note that, a MALA-type proposal pdf can also be used in a MTM scheme. The use of multiple tries improves the mixing of the Markov chain and speeds up the convergence.

Conclusions

We have provided a thorough review of MCMC methods using multiple candidates in order to select the next state of the chain. We have presented and compared different Multiple Try Metropolis, Ensemble MCMC and Delayed Rejection Metropolis schemes. We have also described the Group Metropolis Sampling technique which generates a chain of set of weighted samples, so that some candidates are properly reused in the final estimators. Furthermore, we have shown how the Particle Metropolis-Hastings algorithm can be interpreted as an MTM scheme using a particle filter for generating the different weighted candidates. Several connections and differences have been pointed out. Finally, we have tested several techniques in different numerical experiments: two toy examples in order to provide an exhaustive comparison among the methods, a numerical example regarding the hyperparameter selection for a Gaussian Process (GP) regression model, and two localization problems, one of them involving a real data analysis.

Acknowledgements

This work has been supported by the European Research Council (ERC) through the ERC Consolidator Grant SEDAL ERC-2014-CoG 647423.

References

Appendix A Distribution after resampling

Let us also denote as ${\bm{\theta}}\in\{{\bm{\theta}}^{(1)}\ldots,{\bm{\theta}}^{(N)}\}$, a generic sample after applying one multinomial resampling step according to the normalized IS weights ${\bar{w}}_{n}$, $n=1,\ldots,N$. The density of ${\bm{\theta}}$ is given by

We also define also the matrix ${\bf m}_{\neg n}=[{\bm{\theta}}^{(1)},\ldots,{\bm{\theta}}^{(n)},{\bm{\theta}}^{(n+1)},\ldots,{\bm{\theta}}^{(N)}],$ containing all the samples except for the $n$-th. After some straightforward rearrangements, Eq. (55) can be rewritten as

where $\widehat{Z}=\frac{1}{N}\sum_{n=1}^{N}\frac{\pi({\bm{\theta}}^{(n)})}{q({\bm{\theta}}^{(n)})}$ that is that IS estimator of $Z$. The equation above represents the density of a resampled particle ${\bm{\theta}}\in\{{\bm{\theta}}^{(1)}\ldots,{\bm{\theta}}^{(N)}\}$. Note that if $\widehat{Z}=Z$ then $\widetilde{q}({\bm{\theta}})=\pi({\bm{\theta}})$. Clearly, for a finite value of $N$, there exists a discrepancy between $\widetilde{q}({\bm{\theta}})$ and $\bar{\pi}({\bm{\theta}})$, but this discrepancy decreases as $N$ grows.

Appendix B Particle Filtering

Given a proposal of type $q({\bf x})=q_{1}(x_{1})\prod_{d=2}^{D}q_{d}(x_{d}|x_{d-1})$, and a sample ${\bf x}^{(n)}=x_{1:D}^{(n)}\sim q({\bf x})$ with $x_{d}^{(n)}\sim q_{d}(x_{d}|x_{d-1})$, we assign the importance weight

The weight above can be compute with a recursive procedure for computing the importance weights: starting with $w_{1}^{(n)}=\frac{\pi(x_{1}^{(n)})}{q(x_{1}^{(n)})}$ and then

for $d=2,\ldots,D$. Let also define the partial target pdfs

SIR procedure. In SIR, a.k.a., standard particle filtering, resampling steps are incorporated during the recursion as shown of Table 19 . In general, the resampling steps are applied only in certain iterations in order to avoid the path degeneration, taking into account an approximation $\widehat{ESS}$ of the Effective Sampling Size (ESS) . If $\widehat{ESS}$ is smaller than a pre-established threshold, the particles are resampled. Two examples of ESS approximation are $\widehat{ESS}=\frac{1}{\sum_{n=1}^{N}(\bar{w}_{d}^{(n)})^{2}}$ and $\widehat{ESS}=\frac{1}{\max\bar{w}_{d}^{(n)}}$ where $\bar{w}_{d}^{(n)}=\frac{w_{d}^{(n)}}{\sum_{i=1}^{N}w_{d}^{(i)}}$ (note that $1\leq\widehat{ESS}\leq N$). Hence, the condition for the adaptive resampling can be expressed as $\widehat{ESS}<\eta N$ where $\eta\in$. SIS is given when $\eta=0$ and SIR for $\eta\in(0,1]$. When $\eta=1$, the resampling is applied at each iteration and in this case SIR is often called bootstrap particle filter . If $\eta=0$, no resampling steps are applied, and we have the SIS method described above.

Note that in Table 19, we have employed a proper weighting for resampling particles ,

Generally, it is remarked that $w_{d}^{(1)}=w_{d}^{(2)}=\ldots=w_{d}^{(n)}$ but a specific value is not given. If a different value $c\neq\widehat{Z}_{d}$ is employed, i.e., $w_{d}^{(1)}=\ldots=w_{d}^{(n)}=c$, the algorithm is still valid but the weight recursion loses part of the statistical meaning. This is the reason why the marginal likelihood estimator $\widehat{Z}=\widehat{Z}_{D}=\frac{1}{N}\sum\limits_{n=1}^{N}w_{D}^{(n)}$ is consistent only, if a proper weighting after resampling is used .

In SIS, both estimators are equivalent $\overline{Z}_{d}\equiv\widehat{Z}_{d}$.

Indeed, the classical IS estimator of the normalizing constant $Z_{d}$ at the $d$-th iteration is

An alternative formulation, denoted as $\overline{Z}_{d}$, is often used

where we have employed ${\bar{w}}_{j-1}^{(n)}=\frac{w_{j-1}^{(n)}}{\sum_{i=1}^{N}w_{j-1}^{(i)}}$ and $w_{j}^{(n)}=w_{j-1}^{(n)}\beta_{j}^{(n)}$ .

Furthermore, note that $\overline{Z}_{d}$ can be written in a recursive form as

B.2 Marginal likelihood estimators in SIR

If a proper weighting after resampling is applied in SIR, both formulations $\widehat{Z}_{d}$ and $\overline{Z}_{d}$ in Eqs. (66)-(67) provide consistent estimator of $Z_{d}$ and they are equivalent, $\widehat{Z}_{d}\equiv\overline{Z}_{d}$ (as in SIS).

If a proper weighting is not applied, only

is a consistent estimator of $Z_{d}$, in SIR. In this case, $\widehat{Z}_{d}=\frac{1}{N}\sum_{n=1}^{N}w_{d}^{(n)}$ is not a possible alternative (without using a proper weighting after resampling). However, considering the proper weighting of the resampled particles, then $\widehat{Z}_{d}$ is also a consistent estimator of $Z_{d}$ and it is equivalent to $\overline{Z}_{d}$. Below, we analyze three cases:

No Resampling ($\eta=0$): this scenario corresponds to SIS where $\widehat{Z}_{d}$, $\overline{Z}_{d}$ are equivalent as shown in Eq. (71).

Resampling at each iteration ($\eta=1$): using the proper weighting, $w_{d-1}^{(n)}=\widehat{Z}_{d-1}$ for all $n$ and for all $d$, and replacing in Eq. (68) we have

Since after resampling all particles have the same weight, we have ${\bar{w}}_{d-1}^{(n)}=\frac{1}{N}$ for all $n$. Replacing it in the expression of $\overline{Z}_{d}$ in (72), we obtain

that coincides with $\widehat{Z}_{d}$ in Eq. (74).

Adaptive resampling ($0<\eta<1$): for the sake of simplicity, let us start considering a unique resampling step applied at the $k$-th iteation with $k<d$. We check if both estimators are equal at $d$-th iteration of the recursion. Due to Eq. (71), we have $\overline{Z}_{k}\equiv\widehat{Z}_{k}$, We consider to compute the estimators before the resampling. since before the $k$-th iteration no resampling has been applied. With the proper weighting $w_{k}^{(n)}=\widehat{Z}_{k}$ for all $n$, at the next iteration we have

so that the estimators are equivalent also at the $(k+1)$-th iteration, $\overline{Z}_{k+1}\equiv\widehat{Z}_{k+1}$. Since we are assuming no resampling steps after the $k$-th iteration and until the $d$-th iteration, we have that $\overline{Z}_{i}\equiv\widehat{Z}_{i}$ for $i=k+2,\ldots,d$ due to we are in a SIS scenario for $i>k$ (see Eq. (71)). This reasoning can be easily extended for different number of resampling steps.

Appendix C Consistency of GMS estimators

Dynamic of GMS. We have already seen that we can recover an I-MTM chain from the GMS outputs applying one resampling step for each $t$ when $\mathcal{S}_{t}\neq\mathcal{S}_{t-1}$, i.e.,

for $t=1,\ldots,T$. The sequence $\{{\bm{\theta}}_{t}\}_{t=1}^{T}$ is a chain obtained by one run of an I-MTM2 technique. Note that (a) the sample generation, (b) the acceptance probability function and hence (c) the dynamics of GMS exactly coincide with the corresponding steps of I-MTM2 (or PMH; depending on candidate generation procedure). Hence, the ergodicity of the recovered chain is ensured. Parallel chains from GMS outputs. As described in Section 4.1.3, we can extend the consideration above for generation $C$ parallel I-MTM2 chains. Indeed, we resample $C$ times instead of only one, i.e.,

for $c=1,\ldots,C$, where the super-index denotes the $c$-th chain (similar procedures have been suggested in ). Clearly, the resulting $C$ parallel chains are not independent, and there is an evident loss of performance w.r.t. the case of independent chains. However, at each iteration, the number of target evaluations per iteration is only $N$ instead of $NC$. Note that that each chain in ergodic, so that each estimator $\widetilde{I}_{T}^{(c)}=\frac{1}{T}\sum_{t=1}^{T}g({\bm{\theta}}_{t}^{(c)})$ is consistent (i.e., convergence to the true value for $T\rightarrow\infty$). As a consequence, the arithmetic mean of consistent estimators,

is also consistent, for all values of $C\geq 1$. GMS as limit case. Let us consider the case $\mathcal{S}_{t}\neq\mathcal{S}_{t-1}$ (the other one is trivial), at some iteration $t$. In this scenario, the samples of the $C$ parallel I-MTM2 chains, ${\bm{\theta}}_{t}^{(1)}$,${\bm{\theta}}_{t}^{(2)}$,…,${\bm{\theta}}_{t}^{(C)}$, are obtained by resampled independently $C$ samples from the set $\{{\bm{\theta}}_{1,t},\ldots,{\bm{\theta}}_{N,t}\}$ according to the normalized weights $\bar{\rho}_{n,t}=\frac{\rho_{n,t}}{\sum_{i=1}^{N}\rho_{i,t}}$, for $n=1,\ldots,N$. Recall that the samples ${\bm{\theta}}_{t}^{(1)}$,${\bm{\theta}}_{t}^{(2)}$,…,${\bm{\theta}}_{t}^{(C)}$, will be used in the final estimator $\widetilde{I}_{C,T}$ in Eq. (77).

Let us denote as $\#j$ the number of times that a specific candidate ${\bm{\theta}}_{j,t}$ (contained in the set $\{{\bm{\theta}}_{n,t}\}_{n=1}^{N}$) has been selected as state of one of $C$ chains, at the $t$ iteration. As $C\rightarrow\infty$, The fraction $\frac{\#j}{C}$ approaches exactly the corresponding weights $\bar{\rho}_{j,t}$. Then, for $C\rightarrow\infty$, we have that the estimator in Eq. (77) approaches the GMS estimator, i.e.,

Since $\widetilde{I}_{C,T}$ as $T\rightarrow\infty$ is consistent for all values of $C$, then the GMS estimator is also consistent (and it can be obtained as $C\rightarrow\infty$).