Sequential Monte Carlo for Graphical Models

Introduction

Bayesian inference in statistical models involving a large number of latent random variables is in general a difficult problem. This renders inference methods that are capable of efficiently utilizing structure important tools. Probabilistic Graphical Models (PGMs) are an intuitive and useful way to represent and make use of underlying structure in probability distributions with many interesting areas of applications .

Our main contribution is a new framework for constructing non-standard (auxiliary) target distributions of PGMs, utilizing what we call a sequential decomposition of the underlying factor graph, to be targeted by a sequential Monte Carlo (SMC) sampler. This construction enables us to make use of SMC methods developed and studied over the last $20$ years, to approximate the full joint distribution defined by the PGM. As a byproduct, the SMC algorithm provides an unbiased estimate of the partition function (normalization constant). We show how the proposed method can be used as an alternative to standard methods such as the Annealed Importance Sampling (AIS) proposed in , when estimating the partition function. We also make use of the proposed SMC algorithm to design efficient, high-dimensional MCMC kernels for the latent variables of the PGM in a particle MCMC framework. This enables inference about the latent variables as well as learning of unknown model parameters in an MCMC setting.

During the last decade there has been substantial work on how to leverage SMC algorithms to solve inference problems in PGMs. The first approaches were PAMPAS and nonparametric belief propagation by Sudderth et al. . Since then, several different variants and refinements have been proposed by e.g. Briers et al. , Ihler and Mcallester , Frank et al. . They all rely on various particle approximations of messages sent in a loopy belief propagation algorithm. This means that in general, even in the limit of Monte Carlo samples, they are approximate methods. Compared to these approaches our proposed methods are consistent and provide an unbiased estimate of the normalization constant as a by-product.

Another branch of SMC-based methods for graphical models has been suggested by Hamze and de Freitas . Their method builds on the SMC sampler by Del Moral et al. , where the initial target is a spanning tree of the original graph and subsequent steps add edges according to an annealing schedule. Everitt extends these ideas to learn parameters using particle MCMC . Yet another take is provided by Carbonetto and de Freitas , where an SMC sampler is combined with mean field approximations. Compared to these methods we can handle both non-Gaussian and/or non-discrete interactions between variables and there is no requirement to perform MCMC steps within each SMC step.

The left-right methods described by Wallach et al. and extended by Buntine to estimate the likelihood of held-out documents in topic models are somewhat related in that they are SMC-inspired. However, these are not actual SMC algorithms and they do not produce an unbiased estimate of the partition function for finite sample set. On the other hand, a particle learning based approach was recently proposed by Scott and Baldridge and it can be viewed as a special case of our method for this specific type of model.

Graphical models

A graphical model is a probabilistic model which factorizes according to the structure of an underlying graph $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$, with vertex set $\mathcal{V}$ and edge set $\mathcal{E}$. By this we mean that the joint probability density function (PDF) of the set of random variables indexed by $\mathcal{V}$, $X_{\mathcal{V}}:=\{x_{1},\,\dots,\,x_{|\mathcal{V}|}\}$, can be represented as a product of factors over the cliques of the graph:

We will frequently use the notation $X_{I}=\bigcup_{i\in I}\{x_{i}\}$ for some subset $I\subseteq\{1,\,\dots,\,|\mathcal{V}|\}$ and we write $\mathsf{X}_{I}$ for the range of $X_{I}$ (i.e., $X_{I}\in\mathsf{X}_{I})$. To make the interactions between the random variables explicit we define a factor graph $\mathcal{F}=\{\mathcal{V},\Psi,\mathcal{E}^{\prime}\}$ corresponding to $\mathcal{G}$. The factor graph consists of two types of vertices, the original set of random variables $X_{\mathcal{V}}$ and the factors $\Psi=\{\psi_{C}:C\in\mathcal{C}\}$. The edge set $\mathcal{E}^{\prime}$ consists only of edges from variables to factors. In Figure 1(a) we show a simple toy example of an undirected graphical model, and one possible corresponding factor graph, Figure 1(b), making the dependencies explicit. Both directed and undirected graphs can be represented by factor graphs.

Sequential Monte Carlo

In this section we propose a way to sequentially decompose a graphical model which we then make use of to design an SMC algorithm for the PGM.

SMC methods can be used to approximate a sequence of probability distributions on a sequence of probability spaces of increasing dimension. This is done by recursively updating a set of samples—or particles—with corresponding nonnegative importance weights. The typical scenario is that of state inference in state-space models, where the probability distributions targeted by the SMC sampler are the joint smoothing distributions of a sequence of latent states conditionally on a sequence of observations; see e.g., Doucet and Johansen for applications of this type. However, SMC is not limited to these cases and it is applicable to a much wider class of models.

To be able to use SMC for inference in PGMs we have to define a sequence of target distributions. However, these target distributions do not have to be marginal distributions under $p(X_{\mathcal{V}})$. Indeed, as long as the sequence of target distributions is constructed in such a way that, at some final iteration, we recover $p(X_{\mathcal{V}})$, all the intermediate target distributions may be chosen quite arbitrarily.

This is key to our development, since it lets us use the structure of the PGM to define a sequence of intermediate target distributions for the sampler. We do this by a so called sequential decomposition of the graphical model. This amounts to simply adding factors to the target distribution, from the product of factors in (1), at each step of the algorithm and iterate until all the factors have been added. Constructing an artificial sequence of intermediate target distributions for an SMC sampler is a simple, albeit underutilized, idea as it opens up for using SMC samplers for inference in a wide range of probabilistic models; see e.g., Bouchard-Côté et al. , Del Moral et al. for a few applications of this approach.

Given a graph $\mathcal{G}$ with cliques $\mathcal{C}$, let $\{\psi_{k}\}_{k=1}^{K}$ be a sequence of factors defined as follows $\psi_{k}(X_{\mathcal{I}_{k}})=\prod_{C\in\mathcal{C}_{k}}\psi_{C}(X_{C})$, where $\mathcal{C}_{k}\subset\mathcal{C}$ are chosen such that $\bigcup_{k=1}^{K}\mathcal{C}_{k}=\mathcal{C}$ and $\mathcal{C}_{i}\cap\mathcal{C}_{j}=\emptyset,~{}i\neq j$, and where $\mathcal{I}_{k}\subseteq\{1,\,\dots,\,|\mathcal{V}|\}$ is the index set of the variables in the domain of $\psi_{k}$, $\mathcal{I}_{k}=\bigcup_{C\in\mathcal{C}_{k}}C$. We emphasize that the cliques in $\mathcal{C}$ need not be maximal. In fact even auxiliary factors may be introduced to allow for e.g. annealing between distributions. It follows that the PDF in (1) can be written as $p(X_{\mathcal{V}})=\frac{1}{Z}\prod_{k=1}^{K}\psi_{k}(X_{\mathcal{I}_{k}})$. Principally, the choices and the ordering of the $\mathcal{C}_{k}$’s is arbitrary, but in practice it will affect the performance of the proposed sampler. However, in many common PGMs an intuitive ordering can be deduced from the structure of the model, see Section 5.

2 Sequential Monte Carlo for PGMs

At iteration $k$, the SMC sampler approximates the target distribution $\bar{\gamma}_{k}$ by a collection of weighted particles $\{X_{\mathcal{L}_{k}}^{i},w_{k}^{i}\}_{i=1}^{N}$. These samples define an empirical point-mass approximation of the target distribution. In what follows, we shall use the notation $\xi_{k}:=X_{\mathcal{I}_{k}\setminus\mathcal{L}_{k-1}}$ to refer to the collection of random variables that are in the domain of $\gamma_{k}$, but not in the domain of $\gamma_{k-1}$. This corresponds to the collection of random variables, with which the particles are augmented at each iteration.

After having performed this procedure for the $N$ ancestor indices and particles, they are assigned importance weights $w_{k}^{i}=W_{k}(X_{\mathcal{L}_{k}}^{i})$. The weight function, for $k\geq 2$, is given by

where, again, we write $\xi_{k}=X_{\mathcal{I}_{k}\setminus\mathcal{L}_{k-1}}$. We give a summary of the SMC method in Algorithm 1.

In the case that $\mathcal{I}_{k}\setminus\mathcal{L}_{k-1}=\emptyset$ for some $k$, resampling and propagation steps are superfluous. The easiest way to handle this is to simply skip these steps and directly compute importance weights. An alternative approach is to bridge the two target distributions $\bar{\gamma}_{k-1}$ and $\bar{\gamma}_{k}$ similarly to Del Moral et al. .

Since the proposed sampler for PGMs falls within a general SMC framework, standard convergence analysis applies. See e.g., Del Moral for a comprehensive collection of theoretical results on consistency, central limit theorems, and non-asymptotic bounds for SMC samplers.

The choices of proposal density and adjustment multipliers can quite significantly affect the performance of the sampler.

3 Estimating the partition function

The partition function of a graphical model is a very interesting quantity in many applications. Examples include likelihood-based learning of the parameters of the PGM, statistical mechanics where it is related to the free energy of a system of objects, and information theory where it is related to the capacity of a channel. However, as stated by Hamze and de Freitas , estimating the partition function of a loopy graphical model is a “notoriously difficult” task. Indeed, even for discrete problems simple and accurate estimators have proved to be elusive, and MCMC methods do not provide any simple way of computing the partition function.

On the contrary, SMC provides a straightforward estimator of the normalizing constant (i.e. the partition function), given as a byproduct of the sampler according to,

It may not be obvious to see why (3) is a natural estimator of the normalizing constant $Z_{k}$. However, a by now well known result is that this SMC-based estimator is unbiased. This result is due to Del Moral [21, Proposition 7.4.1] and, for the special case of inference in state-space models, it has also been established by Pitt et al. . For completeness we also offer a proof using the present notation in the appendix.

Since $Z_{K}=Z$, we thus obtain an estimator of the partition function of the PGM at iteration $K$ of the sampler. Besides from being unbiased, this estimator is also consistent and asymptotically normal; see Del Moral .

In we have studied a specific information theoretic application (computing the capacity of a two-dimensional channel) and inspired by the algorithm proposed here we were able to design a sampler with significantly improved performance compared to the previous state-of-the-art.

Particle MCMC and partial blocking

PMCMC methods enable blocking of the latent variables of the PGM in an MCMC scheme. Simulating all the latent variables $X_{\mathcal{L}_{K}}$ jointly is useful since, in general, this will reduce the autocorrelation when compared to simulating the variables $x_{j}$ one at a time . However, it is also possible to employ PMCMC to construct an algorithm in between these two extremes, a strategy that we believe will be particularly useful in the context of PGMs. Let $\{\mathcal{V}^{m},\,m\in\{1,\,\dots,\,M\}\}$ be a partition of $\mathcal{V}$. Ideally, a Gibbs sampler for the joint distribution $p(X_{\mathcal{V}})$ could then be constructed by simulating, using a systematic or a random scan, from the conditional distributions

We refer to this strategy as partial blocking, since it amounts to simulating a subset of the variables, but not necessarily all of them, jointly. Note that, if we set $M=|\mathcal{V}|$ and $\mathcal{V}^{m}=\{m\}$ for $m=1,\,\dots,\,M$, this scheme reduces to a standard Gibbs sampler. On the other extreme, with $M=1$ and $\mathcal{V}^{1}=\mathcal{V}$, we get a fully blocked sampler which targets directly the full joint distribution $p(X_{\mathcal{V}})$.

From (1) it follows that the conditional distributions (4) can be expressed as

where $\mathcal{C}^{m}=\{C\in\mathcal{C}:C\cap\mathcal{V}^{m}\neq\emptyset\}$. While it is in general not possible to sample exactly from these conditionals, we can make use of PMCMC to facilitate a partially blocked Gibbs sampler for a PGM. By letting $p(X_{\mathcal{V}^{m}}|X_{\mathcal{V}\setminus\mathcal{V}^{m}})$ be the target distribution for the SMC sampler of Algorithm 1, we can construct a PMCMC kernel $P_{N}^{m}$ that leaves the conditional distribution (5) invariant. This suggests the following approach: with $X_{\mathcal{V}}^{\prime}$ being the current state of the Markov chain, update block $m$ by sampling

Here we have indicated explicitly in the notation that the PMCMC kernel for the conditional distribution $p(X_{\mathcal{V}^{m}}|X_{\mathcal{V}\setminus\mathcal{V}^{m}})$ depends on both $X_{\mathcal{V}\setminus\mathcal{V}^{m}}^{\prime}$ (which is considered to be fixed throughout the sampling procedure) and on $X_{\mathcal{V}^{m}}^{\prime}$ (which defines the current state of the PMCMC procedure).

As mentioned above, while being generally applicable, we believe that partial blocking of PMCMC samplers will be particularly useful for PGMs. The reason is that we can choose the vertex sets $\mathcal{V}^{m}$ for $m=1,\,\dots,\,M$ in order to facilitate simple sequential decompositions of the induced subgraphs. For instance, it is always possible to choose the partition in such a way that all the induced subgraphs are chains.

Experiments

In this section we evaluate the proposed SMC sampler on three examples to illustrate the merits of our approach. Additional details and results are available in the appendix and code to reproduce results can be found in . We first consider an example from statistical mechanics, the classical XY model, to illustrate the impact of the sequential decomposition. Furthermore, we profile our algorithm with the “gold standard” AIS and Annealed Sequential Importance Resampling (ASIR ASIR is a specific instance of the SMC sampler by , corresponding to AIS with the addition of resampling steps, but to avoid confusion with the proposed method we choose to refer to it as ASIR.) . In the second example we apply the proposed method to the problem of scoring of topic models, and finally we consider a simple toy model, a Gaussian Markov random field (MRF), which illustrates that our proposed method has the potential to significantly decrease correlations between samples in an MCMC scheme. Furthermore, we provide an exact SMC-approximation of the tree-sampler by Hamze and de Freitas and thereby extend the scope of this powerful method.

The classical XY model (see e.g. ) is a member in the family of n-vector models used in statistical mechanics. It can be seen as a generalization of the well known Ising model with a two-dimensional electromagnetic spin. The spin vector is described by its angle $x\in(-\pi,\pi]$. We will consider square lattices with periodic boundary conditions. The joint PDF of the classical XY model with equal interaction is given by

where $\beta$ denotes the inverse temperature.

To evaluate the effect of different sequence orders on the accuracy of the estimates of the log-normalizing-constant $\log Z$ we ran several experiments on a $16\times 16$ XY model with $\beta=1.1$ (approximately the critical inverse temperature ). For simplicity we add one node at a time and all factors bridging this node with previously added nodes. Full adaptation in this case is possible due to the optimal proposal being a von Mises distribution. We show results for the following cases: Random neighbour (RND-N) First node selected randomly among all nodes, concurrent nodes selected randomly from the set of nodes with a neighbour in $X_{\mathcal{L}_{k-1}}$. Diagonal (DIAG) Nodes added by traversing diagonally ($45^{\circ}$ angle) from left to right. Spiral (SPIRAL) Nodes added spiralling in towards the middle from the edges. Left-Right (L-R) Nodes added by traversing the graph left to right, from top to bottom.

We also give results of AIS with single-site-Gibbs updates and $1\thinspace 000$ annealing distributions linearly spaced from zero to one, starting from a uniform distribution (geometric spacing did not yield any improvement over linear spacing for this case). The “true value” was estimated using AIS with $10\thinspace 000$ intermediate distributions and $5\thinspace 000$ importance samples. We can see from the results in Figure 4 that designing a good sequential decomposition for the SMC sampler is important. However, the intuitive and fairly simple choice L-R does give very good results comparable to that of AIS.

Furthermore, we consider a larger size of $64\times 64$ and evaluate the performance of the L-R ordering compared to AIS and the ASIR method. Figure 5 displays box-plots of $10$ independent runs. We set $N=10^{5}$ for the proposed SMC sampler and then match the computational costs of AIS and ASIR with this computational budget. A fair amount of time was spent in tuning the AIS and ASIR algorithms; $10\thinspace 000$ linear annealing distributions seemed to give best performance in these cases.

We can see that the L-R ordering gives results comparable to fairly well-tuned AIS and ASIR algorithms; the ordering of the methods depending on the temperature of the model. One option that does make the SMC algorithm interesting for these types of applications is that it can easily be parallelized over the particles, whereas AIS/ASIR has limited possibilities of parallel implementation over the (crucial) annealing steps.

2 Likelihood estimation in topic models

Topic models such as Latent Dirichlet Allocation (LDA) are popular models for reasoning about large text corpora. Model evaluation is often conducted by computing the likelihood of held-out documents w.r.t. a learnt model. However, this is a challenging problem on its own—which has received much recent interest —since it essentially corresponds to computing the partition function of a graphical model; see Figure 6. The SMC procedure of Algorithm 1 can used to solve this problem by defining a sequential decomposition of the graphical model. In particular, we consider the decomposition corresponding to first including the node $\theta$ and then, subsequently, introducing the nodes $z_{1}$ to $z_{M}$ in any order. Interestingly, if we then make use of a Rao-Blackwellization over the variable $\theta$, the SMC sampler of Algorithm 1 reduces exactly to a method that has previously been proposed for this specific problem . In , the method is derived by reformulating the model in terms of its sufficient statistics and phrasing this as a particle learning problem; here we obtain the same procedure as a special case of the general SMC algorithm operating on the original model.

We use the same data and learnt models as Wallach et al. , i.e. 20 newsgroups, and PubMed Central abstracts (PMC). We compare with the Left-Right-Sequential (LRS) sampler , which is an improvement over the method proposed by Wallach et al. . Results on simulated and real data experiments are provided in Figure 7. For the simulated example (Figure 7(a)), we use a small model with 10 words and 4 topics to be able to compute the exact log-likelihood. We keep the number of particles in the SMC algorithm equal to the number of Gibbs steps in LRS; this means LRS is about an order-of-magnitude more computationally demanding than the SMC method. Despite the fact that the SMC sampler uses only about a tenth of the computational time of the LRS sampler, it performs significantly better in terms of estimator variance.

The other two plots show results on real data with $10$ held-out documents for each dataset. For a fixed number of Gibbs steps we choose the number of particles for each document to make the computational cost approximately equal. Run #2 has twice the number of particles/samples as in run #1. We show the mean of $10$ runs and error-bars estimated using bootstrapping with $10\thinspace 000$ samples. Computing the logarithm of $\hat{Z}$ introduces a negative bias, which means larger values of $\log\hat{Z}$ typically implies more accurate results. The results on real data do not show the drastic improvement we see in the simulated example, which could be due to degeneracy problems for long documents. An interesting approach that could improve results would be to use an SMC algorithm tailored to discrete distributions, e.g. Fearnhead and Clifford .

3 Gaussian MRF

Finally, we consider a simple toy model to illustrate how the SMC sampler of Algorithm 1 can be incorporated in PMCMC sampling. We simulate data from a zero mean Gaussian $10\times 10$ lattice MRF with observation and interaction standard deviations of $\sigma_{i}=1$ and $\sigma_{ij}=0.1$ respectively. We use the proposed SMC algorithm together with the PMCMC method by Lindsten et al. . We compare this with standard Gibbs sampling and the tree sampler by Hamze and de Freitas .

We use a moderate number of $N=50$ particles in the PMCMC sampler (recall that it admits the correct invariant distribution for any $N\geq 2$). In Figure 12(b) we can see the empirical autocorrelation funtions (ACF) centered around the true posterior mean for variable $x_{82}$ (selected randomly from among $X_{\mathcal{V}}$; similar results hold for all the variables of the model). Due to the strong interaction between the latent variables, the samples generated by the standard Gibbs sampler are strongly correlated. Tree-sampling and PMCMC with partial blocking show nearly identical gains compared to Gibbs. This is interesting, since it suggest that simulating from the SMC-based PMCMC kernel can be almost as efficient as exact simulation, even using a moderate number of particles. Indeed, PMCMC with partial blocking can be viewed as an exact SMC-approximation of the tree sampler, extending the scope of tree-sampling beyond discrete and Gaussian models. The fully blocked PMCMC algorithm achieves the best ACF, dropping off to zero considerably faster than for the other methods. This is not surprising since this sampler simulates all the latent variables jointly which reduces the autocorrelation, in particular when the latent variables are strongly dependent. However, it should be noted that this method also has the highest computational cost per iteration.

Conclusion

We have proposed a new framework for inference in PGMs using SMC and illustrated it on three examples. These examples show that it can be a viable alternative to standard methods used for inference and partition function estimation problems. An interesting avenue for future work is combining our proposed methods with AIS, to see if we can improve on both.

We would like to thank Iain Murray for his kind and very prompt help in providing the data for the LDA example. This work was supported by the projects: Learning of complex dynamical systems (Contract number: 637-2014-466) and Probabilistic modeling of dynamical systems (Contract number: 621-2013-5524), both funded by the Swedish Research Council.

Appendix A Supplementary material

This appendix contains additional information on the experiments in the main paper as well as a simple and direct proof of the unbiasedness of the partition function estimator $\widehat{Z}^{N}_{k}$, stated in the main manuscript. It should be noted, however, that this result is not new. It has previously been established in a general setting by Del Moral [21, Proposition 7.4.1] and, additionally, by Pitt et al. who provide a more accessible proof for the special case of state-space models. Our proof is similar to that of Pitt et al. , but generalized to the PGM setting that we consider.

The classical XY model, see e.g. and references therein, is a member in the family of n-vector models used in statistical mechanics. It can be seen as a generalization of the well known Ising model with a two-dimensional electromagnetic spin. The spin vector is described by its angle $x\in(-\pi,\pi]$. We will consider a square lattice with periodic boundary conditions, i.e. the first and last row/columns are connected. The individual sites are described by their spin angle.

The full joint PDF of the classical XY model is given by

where $\beta$ is the inverse temperature and $H(X_{\mathcal{V}})$—the Hamiltonian—is a sum of pair-wise interaction described by

where the $J_{ij}$’s are parameters describing interactions between the different sites. For simplicity we set $J_{ij}=J=1$ and estimate the partition function for several sizes and $\beta$.

In our sequence of target distributions we add one variable at a time and all associated factors. A simple example where we alternate left-right, right-left, can be seen in Figure 11. To be specific we choose our sequence of intermediate target distributions as

where $\mathcal{N}_{k}=\{i:(k,i)\in\mathcal{E}\}\cap\mathcal{L}_{k-1}$ denotes the set of neighbours to variable $k$ in $\mathcal{L}_{k-1}$. The quantities $\mu(X_{\mathcal{L}_{k-1}})$ and $\kappa(X_{\mathcal{L}_{k-1}})$ follow from elementary trigonometric operations (sum of cosines). From the above expression we note that, conditionally on $X_{\mathcal{L}_{k-1}}$, the variable $x_{k}$ is von Mises distributed under $\bar{\gamma}_{k}$, with $X_{\mathcal{L}_{k-1}}$-dependent mean $\mu$ and dispersion $\kappa$. This implies that we can employ full adaption of the proposed SMC sampler. This is accomplished by choosing the aforementioned von Mises distribution as proposal distribution $r_{k}(x_{k}|X_{\mathcal{L}_{k-1}})$ and by choosing the corresponding normalizing constants $\nu(X_{\mathcal{L}_{k-1}})=2\pi I_{0}(\kappa(X_{\mathcal{L}_{k-1}}))$ (where $I_{0}$ is the modified Bessel function of order ) as adjustment weights. We use the fully adapted SMC sampler to estimate the partition function of the classical XY model.

We consider four different orderings of the nodes:

The first node is selected randomly among all nodes, concurrent nodes are then selected randomly from the set of nodes with a neighbour in $X_{\mathcal{L}_{k-1}}$.

The nodes are added by traversing from left to right with $45^{\circ}$, see Figure 9(a).

The nodes are added spiralling in towards the middle from the edges, see Figure 9(b).

The nodes are added by traversing the graph from left to right, from top to bottom, see Figure 9(c).

See illustrations of the node orderings displayed in Figure 9 for a $3\times 3$ example, numbers display at what iteration the node is added.

A.1.2 Evaluation of topic models

Here we present some additional results (Figure 10) on the synthetic example for various settings of the number of topics ($T$) and words ($W$). LRS $1$ and LRS $2$ has $10$ and $20$ samples, respectively. The number of particles where set to give comparable computational complexity.

A.1.3 Gaussian Markov random field

Consider a square lattice Gaussian Markov random field (MRF) of size $10\times 10$, given by the relation

with latent variables $X_{\mathcal{V}}=\{x_{1},\,\dots,\,x_{100}\}$ and measurements $Y_{\mathcal{V}}=\{y_{1},\,\dots,\,y_{100}\}$. The graphical representation of the latent variables in this model is shown in Figure 11.

The measurements $Y_{\mathcal{V}}$ where simulated from the model with $\sigma_{i}=1$ and $\sigma_{ij}=0.1$. Given these measurement, we seek the posterior distribution $p(X_{\mathcal{V}}\mid Y_{\mathcal{V}})$. We run four different MCMC samplers to simulate from this distribution; the proposed (fully blocked) PGAS, the proposed PGAS with partial blocking, a standard one-at-a-time Gibbs sampler, and the tree-sampler proposed by Hamze and de Freitas . For the PGAS algorithms we use $N=50$ particles. The tree-sampler exploits the fact that the model is Gaussian and it can thus not be used for arbitrary (non-Gaussian or non-discrete) graphs. By partitioning the graph into disjoint trees (in our case, chains) for which exact inference is possible, the tree-sampler implements an “ideal” partially blocked Gibbs sampler. PGAS with partial blocking can thus be seen as an SMC-based version of the tree-sampler. See Figure 11 for the ordering in the PGAS algorithm and the blocking used for tree-sampling and PGAS with partial blocking, corresponding to a partition of the graph into two chains. The variables are numbered $1,\ldots,100$ from top to bottom, left to right and $\mathcal{L}_{k}$ is taken as the $k$ first indices of $\mathcal{L}_{K}=\{1,\ldots,10,20,19,\ldots,11,21,22,\ldots,100,99,\ldots,91\}$. This ordering gives results very similar to that of the Left-Right ordering explained above.

In Figure 12(b) we can see the empirical autocorrelation funtions (ACFs) centered around the true posterior mean for variable $x_{82}$ (selected randomly from $X_{\mathcal{V}}$). Similar results hold for all the variables of the model. Due to the strong interaction between the latent variables, the samples generated by the standard Gibbs sampler are strongly correlated. Tree-sampling and PGAS with partial blocking show nearly identical gains compared to Gibbs. This is interesting, since it suggest that simulating from the SMC-based PGAS kernel can be almost as efficient as exact simulation, even using a moderate number of particles. We emphasize that the PGAS kernels leave their respective target distributions invariant, i.e. the limiting distributions is the same for all MCMC schemes. The fully blocked PGAS algorithm achieves the best ACF, dropping off to zero considerably faster than for the other methods. This is not surprising since this sampler simulates all the latent variables jointly which reduces the autocorrelation, in particular when the latent variables are strongly dependent.

However, this improvement in autocorrelation comes at a cost. For the fully blocked PGAS kernel, the maximal cardinality of the set $\mathcal{A}_{k}$ (see (16)) is 10 (one full row of variables). For the partially blocked PGAS kernel, on the other hand, $|A_{k}|\equiv 1$ since the variables in each block form a chain. This implies that the fully blocked PGAS sampler is an order of magnitude more computationally involved than the partially blocked PGAS sampler. This trade-off between autocorrelation and computational efficiency has to be taken into account when deciding which algorithm that is most suitable for any given problem.

A.2 Proof of unbiasedness

Recall that we use the convention $\xi_{k}=X_{\mathcal{I}_{k}\setminus\mathcal{L}_{k-1}}$. Define recursively the functions $f_{k}(X_{\mathcal{L}_{k}})\equiv 1$ and,

Using the definition of the weight function (see the main document) we have,

However, from the definition in (12) we have that

A.3 Ancestor sampling

To implement the PGAS sampling procedure, it remains to detail the ancestor sampling step. At each iteration $k\geq 2$, this step amount to generating a value for the ancestor index $a_{k}^{N}$ corresponding to the reference particle. Implicitly, this assigns an artificial history for the “remaining” part of the reference particle $X_{\mathcal{L}_{K}\setminus\mathcal{L}_{k-1}}^{\prime}$, by selecting one of the particles $\{X_{\mathcal{L}_{k-1}}^{i}\}_{i=1}^{N}$ as its ancestor. This results in a complete assignment for the collection of latent variables of the PGM

As shown by Lindsten et al. , the probability distribution from which $a_{k}^{N}$ should be sampled in order to ensure reversibility of the PGAS kernel w.r.t. $\bar{\gamma}_{K}$ is given by,

This expression can be understood as an application of Bayes’ theorem, where $w_{k-1}^{i}$ is the prior probability of particle $X_{\mathcal{L}_{k-1}}^{i}$ and the ratio between the target densities is the unnormalized likelihood of $X_{\mathcal{L}_{K}\setminus\mathcal{L}_{k-1}}^{\prime}$ conditionally on $X_{\mathcal{L}_{k-1}}^{i}$.

To derive an explicit expression for the ancestor sampling probabilities in our setting, note first that,

Now, let $\mathcal{A}_{k}$ be the index set of factors $\psi_{j}$, $j\geq k$ for which any of the variables $X_{\mathcal{L}_{k-1}}$ is in the domain of $\psi_{j}$; formally

It follows that any factor $\psi_{j}$ for which $j\not\in\mathcal{A}_{k}$ is independent of $X_{\mathcal{L}_{k-1}}$. Consequently, we can write (14) as

where $\widetilde{X}_{\mathcal{I}_{j}}^{i}$ is a subset of the variables in (13). In fact, the index set $\mathcal{A}_{k}$ corresponds exactly to the factors $\psi_{j}$ that depend, explicitly, both on the particle $X_{\mathcal{L}_{k-1}}^{i}$ and on the reference particle $X_{\mathcal{L}_{K}\setminus\mathcal{L}_{k-1}}^{\prime}$ (through some of their respective components). Indeed, it is only these factors that hold any information about the likelihood of $X_{\mathcal{L}_{K}\setminus\mathcal{L}_{k-1}}^{\prime}$ given $X_{\mathcal{L}_{k-1}}^{i}$.

The expression (17) is interesting, since it shows that the computational complexity of the ancestor sampling step will depend on the cardinality of the set $\mathcal{A}_{k}$. While this will depend both on the structure of the graph and on the ordering of the factors, it will for many models of interest be of a lower order than the cardinality of $\mathcal{V}$.

References