 # Interdependent Gibbs Samplers

Gibbs sampling, as a model learning method, is known to produce the most accurate results available in a variety of domains, and is a de facto standard in these domains. Yet, it is also well known that Gibbs random walks usually have bottlenecks, sometimes termed "local maxima", and thus samplers often return suboptimal solutions. In this paper we introduce a variation of the Gibbs sampler which yields high likelihood solutions significantly more often than the regular Gibbs sampler. Specifically, we show that combining multiple samplers, with certain dependence (coupling) between them, results in higher likelihood solutions. This side-steps the well known issue of identifiability, which has been the obstacle to combining samplers in previous work. We evaluate the approach on a Latent Dirichlet Allocation model, and also on HMM's, where precise computation of likelihoods and comparisons to the standard EM algorithm are possible.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Gibbs sampling is a standard model learning method in Bayesian Statistics, and in particular in the field of Graphical Models,

[Gelman et al., 2014]

. In the Machine Learning community, it is commonly applied in situations where non sample based algorithms, such as gradient descent and EM are not feasible. Perhaps the most prominent application example is the Latent Dirichlet Allocation (LDA) topic model

[Blei et al., 2002], a model that is routinely used in text analysis, recommender systems, and a variety of other domains. For LDA, Gibbs sampling is considered to be the standard and the most accurate approach (see, for instance, [Smola and Narayanamurthy, 2010], [Papanikolaou et al., 2017]).

### 1.1 Local Maxima

The standard asymptotic theory of random walks guarantees that Gibbs, and more generally, MCMC samplers, eventually produce a sample from the target distribution (see [Gelman et al., 2014]). The number of steps of the sampler that is required to approximate well a sample from the target distribution is known as the mixing time of the sampler. However, it is well known to practitioners that the true mixing times of Gibbs samplers in realistic situations are unpractically high. The typical observed behavior of the random walk is not that one obtains a fresh sample from the target distribution every given number of steps, but rather that the walk converges to a certain configuration of the parameters, and stays in the vicinity of these parameters for extremely long times, which would not have had happened if mixing times where short. This is illustrated in Figure 1, and a related discussion is given in Section 4.1.

In view of this behavior, it is natural to ask how one can improve the quality of the solutions obtained as the local maxima. The simplest classical receipt is to run the sampler multiple times independently, and to choose the best solution. Note that it is usually not clear that running a sampler a few or even a few hundred times will significantly improve the quality of a solution. Moreover, in some situations, such as LDA, even the task of computing the quality of the solution (the loglikelihood in this case) is hard and can not be performed reliably.

Another approach would be to attempt to combine multiple solutions, to produce a single improved solution. However, while this general idea was discussed already in [Griffiths and Steyvers, 2004], at present there is no known way to carry this out, due to an issue known as the identifiability problem. In the next section we introduce the approach taken in this paper, and its relation to identifiability.

### 1.2 Mulitpaths and Identifiability

The approach of this paper, which we call the multipath sampler, is to use multiple, interdependent Gibbs samplers.

Consider a Bayesian setting where we a have generative model for the observed data and a set of latent variables . The model depends on a set of parameters , and this set is endowed with a prior, . We call the set of latent variables a path. For LDA, corresponds to topics, with a Dirichlet prior, and for every token , , is the topic assigned to that token. In a standard (partially collapsed) Gibbs sampler one first samples the parameters from , and then samples topic assignments , from , for each . Here are the values of all the assignments in the path except .

In the multipath sampler, instead of a single set of latent variables we maintain sets,

. We define an appropriate joint distribution on them as described in Section

3. Under this joint distribution, the marginal distribution for each given the data , coincides with the distribution of latent variables in the regular Gibbs sampler.

Next, to sample the latent variables , we first sample , and then use the usual Gibbs sampler independently for each , given . The dependence between the paths is expressed in sampling from , where depends on all the together. This step can be viewed as model averaging. The details are given in Section 3. In particular, we will show that the target distribution of the multipath sampler emphasizes the solutions that have high likelihood for the regular sampler.

As noted earlier, the idea of model averaging in the context Gibbs sampling has been discussed in the literature at least since the paper [Griffiths and Steyvers, 2004]. In particular, it is well known that it is generally impossible to average the results of independent samplers, due to the problem of identifiability (see, for instance, [Griffiths and Steyvers, 2004], [Steyvers and Griffiths, 2007], [Gelman et al., 2014]), which refers to the issue that the model parameters are only defined by the model up to a permutation. For instance, in the case of LDA, suppose and

are two sets of topics estimated by two independent Gibbs samplers. One then might attempt to combine

and , for instance by defining , to obtain a better estimate of the topic. However, even if and represent nearly the same set of topics, in does not necessarily correspond to , but may rather correspond to some other for some . Thus, to perform direct averaging, one would have to find a proper permutation expressing the correspondence, before averaging. Moreover, the situation is further complicated by the fact that and do not necessarily represent similar sets of topics, and thus for many topics in , there simply will be no “similar” topic in to average with.

With the multipath approach, we show for the first time how the identifiability problem can be avoided. Indeed, instead of running the samplers completely independently and attempting to only combine the final results, we allow the samplers to share the parameters during the optimization. This forces the samplers to find a common version of each topic, while still maintaining different and somewhat independent empirical estimates of it.

### 1.3 Experiments

We demonstrate the methods introduced in this paper on HMM and LDA Gibbs samplers, and in both cases the multipath sampler improves on the results of a regular sampler. We now motivate and briefly describe the experiments.

As mentioned earlier, in a typical problem where a Gibbs sampler is used, given the model parameters found by a sampler, there is no precise way to evaluate the likelihood of the data with respect to . See [Wallach et al., 2009], where an equivalent notion of perplexity is discussed. As a result, it is not usually easy to compare two different solutions, and , possibly found by two different methods, and perplexity, if used at all, is never used as a sole method of evaluation.

Therefore, to evaluate the performance of the multipath sampler we first compare it to the regular Gibbs sampler on a model where explicit computation of data likelihoods is possible. Specifically, we use Hidden Markov Models (HMMs), with synthetic data, where the parameters

to be estimated are the emission and transition probabilities. We compare the solutions form the Gibbs samplers to those found by the EM-type Baum-Welch algorithm, and to the ground truth.

Next, we compare the Gibbs and the multipath Gibbs sampler LDA model on synthetic data and measure the closeness of the recovered topics to the ground truth.

In addition, on the well known State of The Union corpus, [Wang and McCallum, 2006], we compare the time concentration of the topics as an additional quality indicator. Specifically, this dataset consists of speeches given between the years 1790 to 2013, and each speech is further divided into multiple documents. When such time information is available, it is natural to ask whether certain topics appear in a specific time period, or their occurrence is spread relatively uniformly over time. We call a topic time-concentrated if most of its occurrences are contained in a short time period. Consider Figure 5 for an example of a relatively well time-concentrated topic.

Next, note that the LDA model itself has no access to timestamps, and that if token assignments were given at random, then the topics would be approximately uniformly spread over time. Therefore, time concentration of the topics may be observed only if the model fits the data well, and better time concentration implies a better fit to the data. In Section 4.2 we show that topics and topic assignments found by the multipath sampler have significantly better time-concentration properties compared to topics from the regular sampler.

We note that time concentration for topics is a well studied subject. For instance, the well known Topics over Time (ToT) model, [Wang and McCallum, 2006], and its numerous extensions, incorporate the timestamp data into the model via appropriate probabilistic constraints. These additional constraints may be viewed as an expression of a prior belief about the data, in particular that a topic’s usage pattern has certain temporal properties. Using such models, one indeed obtains topics that are better concentrated in time compared to standard LDA. However, note that by extending the model, one effectively forces a bias towards learning concentrated topics. In contrast, here we show that one can obtain improved concentration by an improved learning algorithm of LDA itself, without additional assumptions. This demonstrates that time-concentration of the topics is a signal contained in the data, rather than an artifact of the extended model.

The rest of the paper is organized as follows: In Section 2 we review the literature. In Section 3 we formally define the mulipath sampler and study its relation to the regular sampler. Section 4 contains the experiments, and concluding remarks and future work are discussed in Section 5.

## 2 Literature

A few days after the first posting of this paper, we have been informed that a similar approach has already been considered previously. While the experiments and the point of view taken in the discussion in this paper are different, the core idea, Algorithm 1, is identical to the algorithm introduced in [Doucet et al., 2002].

As detailed below, there is a large body of work on performance improvements for Gibbs samplers, and in particular for the Gibbs samplers of the LDA model.

The uncollapsed Gibbs sampler for LDA was introduced in [Pritchard et al., 2000] where it was applied to genetic data. The collapsed Gibbs sampler was introduced in [Griffiths and Steyvers, 2004].

Parallelization was studied, for instance, in [Newman et al., 2007], [Smola and Narayanamurthy, 2010]. In [Neiswanger et al., 2014] parallelization methods with asymptotical exactness guarantees were given. Note that while parallelizing the fully collapsed Gibbs sampler presents several significant issues which are addressed in the above work, for LDA there is a partially collapsed sampler, studied in [Asuncion, 2011], which is straightforward to parallelize. As observed in [Magnusson et al., 2017], the partially collapsed samplers may have a performance, in terms of quality of solutions, comparable to the performance of the fully collapsed sampler.

Methods that exploit sparseness of distributions to speed up the iteration of the sampler were studied, among other work, in [Yao et al., 2009], [Li et al., 2014], and further optimized in [Chen et al., 2016]. These methods perform the same random walk as the collapsed sampler, but the sampling step itself is implemented more efficiently.

In a recent work [Papanikolaou et al., 2017], parameter inference by averaging several solutions was discussed for the case of LDA. It is important to note here that this refers to a very restricted form of averaging. Indeed, in that work, one first obtains a solution from a collapsed Gibbs sampler, and a postprocessing step considers perturbations in this solution at a single coordinate, to obtain some variability in the assignment. Thus all the perturbed solutions differ from the original solution at a single coordinate at most. It is then argued that such perturbations can be averaged, since changing only one assignment coordinate in a large corpus is unlikely to cause the identifiability issues associated with averaging. Clearly, however, due to the same reasons, such perturbations can not change the particular local maxima of the original solution and can only provide small variations around it. Thus, while this approach is useful as a local denoising method, especially in computing topic assignment distributions (rather than topic themselves), and is suggested as such in [Papanikolaou et al., 2017], the overall quality of the solution is the same as the one returned by the original sampler. As discussed in the Introduction, our approach completely avoids the identifiability issues, and allows to combine information from significantly differing topic assignments.

Finally, we note that all of the above mentioned methods can be combined with the methods introduced in this paper. In particular both multipath samplers may be parallelized and may incorporate sparsity in exactly the same way this is done for the regular sampler.

## 3 Multiple Paths

We will assume that the data is generated in a standard Bayesian setting: A set of model parameters is denoted by , endowed with a prior distribution . Given , the model specifies a distribution of latent variables , . Each takes values in the set . By analogy to HMMs, we refer to the latent variables as a path, and denote the set of all possible values of a path sequence by .

Finally, given the parameters and the path, the model specifies the distribution of the observed variables , where . We also refer to as the data variables. Each takes values in the alphabet .

In the case of HMMs, the parameters are be the transition and emission distributions. Specifically, we consider HMMs with states, and with emissions given by discrete distributions on the alphabet . Then , where

is the initial distribution on the Markov chain, for each

, is the transition distribution given the chain is at state , and are emission distributions for state . The prior is given by , , and . The path is a state sequence sampled from the Markov chain, and are the emissions given the state sequence, so that the distribution of given is .

For LDA, the parameters are the topics, , with a prior , for some fixed . The path is the set of topic assignments for each token in the corpus, and are the tokens, where is sampled from the topic .

Next, the following definitions will be useful. We refer to the quantity

 P(w,p|ϕ) (1)

as the path likelihood, and to

 P(w|ϕ)=∑p∈PP(w,p|ϕ) (2)

the data likelihood.

Given the data , we interested in the maximum likelihood parameter values, ,

 ϕmax=\operatornamewithlimitsargmaxϕP(w|ϕ), (3)

and in a sample of given and . A standard approach to obtaining an approximation of is to use a Gibbs sampler to sample from .

Now, on the set of parameters and paths, , define the function

 f(ϕ,p)=P(w,p,ϕ)=P(w,p|ϕ)⋅P(ϕ).

Denote by

 Cf=∑p∫Φf(ϕ,p)dϕ=P(w),

where is the Lebesgue measure on (that is, is a density with respect to ).

Then the probability density is, by definition, the density of the posterior distribution .

Observe that the posterior density of is proportional to the data likelihood of :

 ^f(ϕ)=∑p^f(ϕ,p)=1Cf∑pP(ϕ,p,D)= (4) =1Cf⋅P(ϕ)⋅P(w|ϕ).

We are now ready to introduce the multiple paths model. The generative model consist of simply sampling the parameters from , and then sampling independent copies of the paths, , and the data, . Here the data sequence is sampled from the original generative model given and . Therefore, the model generates paths, and versions of the data. In what follows we will condition this model on the event that all the versions of the data coincide and are equal to a given data sequence .

Specifically, given a data sequence , we seek the posterior distribution of in the path model, given the data , such that for all paths .

For , define on the function

 fm(ϕ,p1,…,pm)=P(ϕ)⋅∏j≤mP(pj,w|ϕ). (5)

Denote by the corresponding normalization constant,

 Cfm=∫Φ∑p1,…,pmfm(ϕ,p1,…,pm)dϕ.

Then gives the distribution of with respect to the multipath generative model, conditioned on the event for all .

The following observation is the key fact of this section, and we record it as a lemma.

###### Lemma 3.1.

Let be a sample from . Then the marginal density of is

 Pm(ϕ|w)=1CfmP(ϕ)P(w|ϕ)m. (6)
###### Proof.

Indeed,

 Cfm⋅Pm(ϕ|w)=∑p1,…,pmfm(ϕ,p1,…,pm)= P(ϕ)∑p1,…,pm∏j≤mP(pj,w|ϕ)= P(ϕ)∏j≤m∑pjP(pj,w|ϕ)= P(ϕ)P(w|ϕ)m.

Note that this generalizes (4). The crucial detail about (6) is that it contains the -th power of the data likelihood. Let , defined in (3), be a maximum likelihood solution given the data, and let be any other parameter. Then,

 (P(w|ϕmax)/P(w|ϕ))m= (P(ϕmax|w)/P(ϕ|w))m.

In words, we are exponentially more likely to obtain as a sample form than from (compared to any other ). Thus, as grows, any sample form should yield with data likelihood close to that of .

The discussion above shows why sampling from yields better solutions than sampling from . We now discuss the Gibbs sampler for , which is a straightforward extension of the sampler for . The sampler is given schematically in Algorithm 1. For , this is the standard Gibbs sampler.

The notation has the standard meaning of without the -th component . For the -path model, our path is the collection of all paths . However, when sampling and individual component from , since the density (5) factorizes over , we obtain:

 P(pji=s|pj−i,p1,…,pj−1,pj+1,…,pm,w,ϕ)= P(pji=s|pj−i,w,ϕ)

for every . Thus, given the parameters , one samples each path independently from the others.

For the special cases of HMM and LDA, we write out Algorithm 1 in full detail in Supplementary Material Section B. In particular, we detail the equations which describe how steps 4 and 7 of Algorithm 1 are performed for these models. As mentioned earlier, these equations are straightforward modifications of the corresponding equations for the standard sampler. Also, note that for LDA, Algorithm 1 corresponds to the so called partially collapsed Gibbs sampler, rather than to the more common fully collapsed version. However, a fully collapsed version of a sampler from is equally easy to derive, and is also presented in Supplementary Material Section B.

As discussed above, the advantage of the multipath sampling comes from the fact that it implicitly samples from a distribution proportional to . Note that while is not typically computable, the path likelihood is computable. It is then natural to ask what will happen if the Gibbs sampler is modified to sample the pair from a distribution proportional to . This can easily be accomplished, using a single path. Such a sampler will spend much more time around the pairs for which is high, compared to the regular Gibbs sampler. However, even if we were given a pair that maximizes globally over all , it is straightforward to verify that would not usually correspond to useful parameters. On the other hand, if the data is generated by a model with parameters , then consistency results (such as, for instance, [Baum and Petrie, 1966] for the case of HMMs) ensure that maximum data likelihood solutions, as in (3), reconstruct . The important feature of multipath sampler is therefore that it samples from , even though is not computable.

Finally, it is worthwhile to remark that the multipath approach bears some resemblance to the EM algorithm. Indeed, in EM algorithm, for each , one computes the distribution of given all the data and parameters . Then, given these distributions, one recomputes the parameters. One can view the multipath scheme as an approximation of the distribution with samples, instead of a computation of it in a closed form as in EM, or an approximation with one sample, as in regular Gibbs. A substantial difference between EM and the sampler is that our samples of are conditioned on the data and the rest of the assignments, , while in EM the distribution is conditioned only on the data. However, interestingly, note that if the assignments are independent given the parameters, as in the case of mixture models, than a multipath sampler directly approximates EM as grows.

## 4 Experiments

In this Section we describe the HMM and LDA experiments. In addition, in Section 4.1, after describing the HMM experiments we discuss Figure 1, referred to in Section 1, which demonstrates the local maxima behavior of the Gibbs sampler.

### 4.1 HMMs

The HMM experiments use synthetic data. The data was generated by an HMM with two states, such that the transition probability between the states is , and the emissions are distributions on symbols. The emissions are show in Figure 3 as the ground truth. All the experiments were conducted using the same fixed random sample of size from this HMM.

These particular parameters were chosen in order to make the HMM hard to learn. Otherwise, all methods would succeed all the time, and would be difficult to compare. This corresponds to realistic situations where different emissions are hard to distinguish.

Note that while may appear as a high number of samples, in this case, due to transitions probabilities which are close to , this high number of samples is in fact necessary. As we discuss later, even with this number of samples, both the EM algorithm and some runs of the Gibbs sampler manage to slightly overfit the solution.

To put the HMM learning problem in the Bayesian setting, we endow the set of parameters with a prior where transition and emission distributions are sampled from a uniform Dirchlet, . In other words, transition distributions are sampled uniformly from a 3-simplex, and emission distributions from an 11-simplex. The data is the given emissions of the HMM, and the path corresponds to the latent sequence of the states of the HMM.

Figure 1 demonstrates the local maxima behavior. It shows the data log-likelihoods of a solution along a run of a collapsed Gibbs sampler, for 16 independent runs of a collapsed Gibbs sampler.

Each run had 200000 iterations, with data log-likelihoods computed every 100 iterations. As discussed in the Introduction, each run converges fairly quickly to an area of fixed likelihood, and stays there for extremely long time intervals. Figure 1: Collapsed Gibbs HMM Data Likelihoods, illustrating the Local Maxima. Measurement is performed every 100 iterations. Total of 200000 iterations.

Figure 2 shows the data log-likelihoods of the parameters for multiple runs of different algorithms after 200000 iterations. Each algorithm was run independently 16 times (all on the same data), and the resulting data log-likelihoods are plotted in ascending order. The algorithms included are the Collapsed Gibbs sampler for 1 and five paths (C 1, 5) and a non-collapsed sampler (PC 1, 5).

For reference, the constant line (BW) corresponds to a solution found by the Baum-Welch algorithm, and Ground line corresponds to the data likelihood for the ground truth parameters. Note that the Ground line is slightly lower than the highest likelihood solutions, even with the relatively high number of samples .

Clearly most of the algorithms improve on the base-line – the Collapsed Gibbs sampler C1. The the non- collapsed sampler, PC1, which is also a base-line, provides results similar to those of C1. The best results are obtained by C 5. The emissions found by C5 run with the highest data likelihood are shown in Figure 3. Note that the imperfect approximation of the emissions may be explained by a mild overfitting due to the barely sufficient amount of data.

### 4.2 Lda

In this Section we describe the LDA experiments.

#### 4.2.1 Synthetic Data

The topics for the synthetic data are shown in Figure 4. They were constructed as follows: Choose 10 overlapping bands of equal length (except the boundary ones) on a dictionary of size , centered at equally spaced intervals over . Each topic is a linear combination of a band, with weight

and a uniform distribution on

, with weight . Most of the points in the dictionary belong to two topics (if one disregards the uniform component present in all topics).

Each document is generated by sampling the topic distribution , with , and sampling 10 tokens by sampling a topic from , and a token from that topic.

We have generated a varying number of documents from this model, and have run the collapsed LDA Gibbs sampler with a varying number of paths. Given the topics found by the sampler, to compute a distance to ground truth,, for each ground truth topic , we find a topic closest to in total variation distance, and take the mean of these distances:

 disc=110∑iminj∥∥βgi−βj∥∥1

In addition, each was averaged over 10 independent runs of the sampler (with same data and same number of paths). Samplers with different number of sources, but same number of documents, were run on the set of documents. All samplers were run for 10000 iterations, and in all cases already after 3000 the topics stabilized and practically did not change further. The results are shown in Table 1. Clearly, increasing the number of paths yields better reconstruction. Increasing the number of paths to more than 5 improved the performance further in some, but not all cases.

#### 4.2.2 SOTU Dataset

Our State of the Union speeches corpus contains speeches from years 1790 to 2013. The preprocessing included stop words removal, stemming, removal of 40 most common words, and words that appeared only once in the corpus. Each speech was divided into paragraphs, and each paragraph was treated as a document. After the preprocessing, the corpus contains around 767000 tokens, in about 21000 documents. There are on average 50-80 documents corresponding to a single year. All the models were fit with topics and prior parameters and .

As discussed in Section 1.3, we compare the time- concentration properties of LDA models on the State of the Union corpus, found by the collapsed Gibbs sampler (referred to as LDA in what follows), and by multipath collapsed Gibbs sampler with

paths (to which we refer as mpLDA). We quantify the time concentration of topics in two different ways, via the computation of yearly entropy, and via the notion of quantile buckets, as discussed below.

Given the topic assignments returned by the sampler, we compute several quantities as discussed below. For each document , let

be the topic distribution ( a probability distribution on the set

) of the document, as computed by the sampler. Specifically, if are the indexes of the tokens in a document , with the size of the document, then

 θd=1k∑l≤kδpidk.

For the multipath sampler, we take the assignments from the first path , for convenience. Typically, after convergence, for most tokens, all paths have an identical assignment for a given token.

Next, for each year, we compute the topic distribution: Let be the set of documents from year . Then the topic distribution of the year is

 θy=1∣∣Yy∣∣∑d∈Yyθd.

We refer to as the topic weight of topic for year . The Figure 5 shows in blue the topic weights for all years of a particular topic in a model found by the collapsed Gibbs sampler. The 9 words with highest probabilities in this topic are indicated in the figure. In addition, we have found the topic in the model returned by the Mulipath Gibbs, and for the aternating sampler which are closest to in total variation distance (as distributions over words). The topic weights of and are shown in orange and green.

As Figure 5 shows, the orange curve has more weight on the years where both blue and orange curves have high values, and it has less weight in other areas. Note that the second statement here (less weight) does not follow from the first (more weight) – the weight of the topic over the years does not have to sum to 1 (the sum over all topics for a fixed year is 1). We conclude that the orange topic is more time concentrated than the blue. Similar conclusion also holds for the green topic.

We now proceed to show that the situation depicted in Figure 5 for a fixed topic, is in fact typical for most topics, and overall the topics of mpLDA are better concentrated than those of LDA. First, consider the entropy of the topic distribution for every year, for all models, as shown in Figure 7 (right). The entropies of mpLDA are similar consistently lower than those of LDA for every year. This means that each year in the new models is composed of slightly fewer topics, with higher weights.

Second, given the topic weights for a fixed topic, we quantify its time concentration via the notion of quantile bucket length. A -quantile bucket of the weights is a time interval which contains a proportion of the total weight of the topic. For instance, Figure 6 show quantile buckets for topics from Figure 5. Shorter buckets mean that a fixed proportion of the weight is contained in a shorter time interval. Note that the short buckets of the orange curve in Figure 6 are shorter than than those of the blue curve.

In Figure 7 (left), the distribution of the bucket lengths, with , for all the buckets from all the topics is shown, for all models. Since the sum of each 20 buckets contributed by a single topic is 224, the total number of years, the expectation of each histogram must be

. We observe that the blue histogram is centered around this value, indicating more or less uniform (over time) weight distributions for most of the topics of LDA. On the other hand, the orange histogram is clearly skewed to the left, which implies that there are more short intervals containing a fixed weight in mpLDA.

Some additional measurements related to this experiment may be found in Supplementary Material Section A.

## 5 Conclusion and Future Work

In this work we have introduced a modification of the Gibbs sampler, which combines samples from multiple sets of latent variables and yet avoids the identifiability problem. We have shown that asymptotically this sampler is exponentially more likely, compared to the regular Gibbs sampler, to return samples which have high data likelihood in the original generative model. We have also shown empirically that the quality of the solutions indeed improves, as the number of paths grows.

In addition, we found that time concentration properties of topic models may improve when using the multipath sampler, compared to regular Gibbs. This is one instance of a phenomenon where a better optimization of a model yields better insight into the data. However, full implications of this fact, in particular in the context of trend analysis, would require a more extensive investigation.

While we have tested our methods on the HMM and LDA models, our methods were formulated in a fairly generic latent variables settings, and could be extended even further. Thus, a natural direction for future work would be to understand the impact of the methods presented here on learning in other models, such as Markov Random Fields or extended HMMs with large state spaces, for which EM algorithms are not feasible.

## References

• [Asuncion, 2011] Asuncion, Jr., A. U. (2011). Distributed and Accelerated Inference Algorithms for Probabilistic Graphical Models. PhD thesis.
• [Baum and Petrie, 1966] Baum, L. E. and Petrie, T. (1966). Statistical inference for probabilistic functions of finite state markov chains. Ann. Math. Statist., 37(6).
• [Blei et al., 2002] Blei, D. M., Ng, A. Y., and Jordan, M. I. (2002). Latent dirichlet allocation. Journal of Machine Learning Research, 3:2003.
• [Chen et al., 2016] Chen, J., Li, K., Zhu, J., and Chen, W. (2016). WarpLDA: A cache efficient o(1) algorithm for latent dirichlet allocation. Proc. VLDB Endow.
• [Doucet et al., 2002] Doucet, A., Godsill, S. J., and Robert, C. P. (2002).

Marginal maximum a posteriori estimation using markov chain monte carlo.

Statistics and Computing.
• [Gelman et al., 2014] Gelman, A., Carlin, J., Stern, H., Dunson, D., Vehtari, A., and Rubin, D. (2014). Bayesian Data Analysis. Third edition.
• [Griffiths and Steyvers, 2004] Griffiths, T. L. and Steyvers, M. (2004). Finding scientific topics. Proceedings of the National Academy of Sciences of the United States of America, 101.
• [Li et al., 2014] Li, A. Q., Ahmed, A., Ravi, S., and Smola, A. J. (2014). Reducing the sampling complexity of topic models. KDD.
• [Magnusson et al., 2017] Magnusson, M., Jonsson, L., Villani, M., and Broman, D. (2017). Sparse partially collapsed mcmc for parallel inference in topic models. Journal of Computational and Graphical Statistics.
• [Neiswanger et al., 2014] Neiswanger, W., Wang, C., and Xing, E. P. (2014). Asymptotically exact, embarrassingly parallel mcmc. UAI.
• [Newman et al., 2007] Newman, D., Asuncion, A., Smyth, P., and Welling, M. (2007). Distributed inference for latent dirichlet allocation. NIPS.
• [Papanikolaou et al., 2017] Papanikolaou, Y., Foulds, J. R., Rubin, T. N., and Tsoumakas, G. (2017). Dense distributions from sparse samples: Improved gibbs sampling parameter estimators for lda. Journal of Machine Learning Research, 18(62).
• [Pritchard et al., 2000] Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics.
• [Smola and Narayanamurthy, 2010] Smola, A. and Narayanamurthy, S. (2010). An architecture for parallel topic models. Proc. VLDB Endow.
• [Steyvers and Griffiths, 2007] Steyvers, M. and Griffiths, T. (2007). Probabilistic Topic Models.
• [Wallach et al., 2009] Wallach, H. M., Murray, I., Salakhutdinov, R., and Mimno, D. (2009). Evaluation methods for topic models. ICML.
• [Wang and McCallum, 2006] Wang, X. and McCallum, A. (2006). Topics over time: a non-markov continuous-time model of topical trends. In Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 424–433. ACM.
• [Yao et al., 2009] Yao, L., Mimno, D., and McCallum, A. (2009). Efficient methods for topic model inference on streaming document collections. KDD.

## Appendix A Additional Considerations Regarding Time Concentration

In Section 4.2.2 we have found that for topics and topic assignments from the multipath Gibbs sampler, the distribution of bucket lengths is more skewed to the left.

However, note that there might have existed degenerate topic assignments, which could have produced the concentration results of Figure 7. Indeed, consider an assignment where there is a single topic that is assigned to a major portion of every document, with the rest of the topics distributed in a concentrated way over the years. Such an assignment would lower the yearly entropy and produce a left skew in bucket lengths, but, arguably, should not be considered a better description of the data. To ensure that this is not the case, for each topic , consider the total topic weight

 wt=∑yθy(t). (7)

In Figure 8, we show the sorted topic weights for both models (right) and the histogram of bucket lengths (left), where each bucket length is weighted proportionally to the topic weight of its topic (rather then by 1, as in computation for Figure 7). The topic weight distribution is practically the same for both models, and the lengths histograms also remain practically unchanged. This shows that the mpLDA genuinely redistributes the same topic weights in a more time concentrated fashion, rather than produces anomalous assignments.

## Appendix B Explicit Implementations

In this Section we describe a fully detailed implementation of a multipath Gibbs sampler Algorithm 1 for LDA. As discussed in Section 3, in order to fully specify Algorithm 1 for a particular generative model, one has to implement the sampling in lines 4 and 7 of the Algorithm. Line 4 corresponds to line 6 in Algorithm 3, and line 7 to lines 12-19.

In addition, we describe the multipath analog of the collapsed Gibbs sampler. Schematically, the multipath collapsed Gibbs sampler is given in Algorithm 2, and Algorithm 4 is a full implementation.

The notation for the algorithms is summarized in Table 2, and is analogous to the one commonly used in the literature.

Note that for both Algorithm 3 and 4 coincide with their standard counterparts.

All computations in the derivation of the algorithms are standard Dirichlet conjugacy type computations. We omit the computations, since they are completely analogous to the computations for the derivation of the regular LDA Gibbs sampler. The algorithms for the HMM model are also derived using similar (slightly simpler) computations.

Please see the next page for the notation and the algorithms.