Firefly Monte Carlo: Exact MCMC with Subsets of Data

1 Introduction

The Bayesian approach to probabilistic modeling is appealing for several reasons: the generative framework allows one to separate out modeling assumptions from inference procedures, outputs include estimates of uncertainty that can be used for decision making and prediction, and it provides clear ways to perform model selection and complexity control. Unfortunately, the fully-Bayesian approach to modeling is often very computationally challenging. It is unusual for non-trivial models of real data to have closed-form posterior distributions. Instead, one uses approximate inference via Monte Carlo, variational approximations, Laplace approximations, or other tools.

Bayesian inference can be costly. In particular, evaluating a new hypothesis usually requires examination of an entire data. For example, when running Metropolis-Hastings (MH) it is necessary to evaluate the target posterior density for each proposed parameter update, and this posterior will usually contain a factor for each datum. Similarly, typical variational Bayesian procedures need to build local approximations for each of the data in order to update the approximation to any global parameters. In both cases, these data-intensive computations are performed many times as part of an iterative procedure.

Ryan P. Adams

School of Engineering and Applied Sciences Harvard University Cambridge, MA 02138

Recent methods have been proposed to partially overcome these difficulties. Stochastic and online variational approximation procedures [Hoffman et al., 2010, 2013] can use subsets of the data to make approximations to global parameters. As these procedures are optimizations, it is possible to build on convergence results from the stochastic optimization literature and achieve guarantees on the resulting approximation. For Markov chain Monte Carlo, the situation is somewhat murkier. Recent work has shown that approximate transition operators based on subsets of data can be used for predictive prefetching to help parallelize MCMC [Angelino et al., 2014]. Other work uses approximate transition operators directly, for Metropolis-Hastings (MH) and related algorithms [Welling and Teh, 2011]. Korattikara et al. [2014] and Bardenet et al. [2014] have shown that such approximate MH moves can lead to stationary distributions which are approximate but that have bounded error, albeit under strong conditions of rapid mixing.

In this paper we present a Markov chain Monte Carlo algorithm, Firefly Monte Carlo (FlyMC), that is in line with these latter efforts to exploit subsets of data. Our transition operator, however, is exact: it leaves the true posterior distribution invariant. FlyMC is a latent variable model which introduces a collection of Bernoulli variables - one for each datum - with conditional distributions chosen so that they effectively turn on and off data points in the posterior, hence "firefly". The introduction of these latent variables does not alter the marginal distribution of the parameters of interest. Our only requirement is that it be possible to provide a "collapsible" lower bound for each likelihood term. FlyMC can lead to dramatic performance improvements in MCMC, as measured in wallclock time.

2 FireflyMonte Carlo

The Firefly Monte Carlo algorithm tackles the problem of sampling from the posterior distribution of a probabilistic model. We denote the parameters of interest as θ and assume that they have prior p(θ). We assume that N data have been observed {xn}Nn=1 and that these data are conditionally independent given θ under a likelihood p(xn | θ). Our target distribution is therefore

We will write the nth likelihood term as a function of θ as

An MCMC sampler makes transitions from a given θ to a new θ′ such that posterior distribution remains invariant. Conventional algorithms, such as Metropolis-Hastings, require evaluation of the unnormalized posterior in full at every iteration. When the data set is large, evaluating all N likelihoods is a computational bottleneck. This is the problem that we seek to solve with FlyMC.

For each data point, n, we introduce a binary auxiliary variable, zn ∈ {0, 1}, and a function Bn(θ) which is a strictly positive lower bound on the nth likelihood: 0 < Bn(θ) ≤ Ln(θ). Each zn has the following Bernoulli distribution conditioned on the parameters:

.

We now augment the posterior distribution with these N variables:

This joint distribution has a remarkable property: to evaluate the probability density over θ, conditioned on {zn}Nn=1 , it is only necessary to evaluate those likelihood terms for which zn = 1. Consider factor n from the product above.

The likelihood term Ln(θ) only appears in those factors for which zn = 1 and we can think of these data as forming a minibatch subsample of the full set. If zn = 0 for most n, transition updates for the parameters θ will be much cheaper, as these are applied to p(θ | {xn, zn}Nn=1 ).

Of course, we do have to evaluate all N bounds Bn(θ) at each iteration. We seem to have just shifted the computational burden from evaluating the Ln(θ) to evaluating the Bn(θ). However, if we choose Bn(θ) to have a convenient form such as a scaled Gaussian, then the full product ∏Nn=1 Bn(θ) can be computed for each new θ in O(1) time using the sufficient statistics of the distribution, which only need to be computed once. To make this clearer, we can rearrange the joint distribution in terms of a "pseudo-prior," ˜ p(θ) and "pseudolikelihood," ˜Ln(θ) as follows:

where the product only runs over those n for which zn = 1, and we have defined

We can generate a Markov chain for the joint distribution in Equation (2) by alternating between updates of θ conditional on {zn}Nn=1 , which can be done with any conventional MCMC algorithm, and updates of {zn}Nn=1 conditional on θ for which we discuss efficient methods in Section 3.2. We emphasize that the marginal distribution over θ is still the correct posterior distribution given in Equation (1).

At a given iteration, the zn = 0 data points are "dark": we simulate the Markov chain without computing their likelihoods. Upon a Markov transition in the space of {zn}Nn=1 , a smattering of these dark data points become "bright" with their zn = 1, and we include their likelihoods in subsequent iterations. The evolution of the chain evokes an image of fireflies, as the individual data blink on and off.

The details of choosing a lower bound and efficiently sampling the {zn} are treated in the proceeding sections, but the high-level picture is now complete. Figure 1 illustrates the augmented space, and a simple version of the algorithm is shown in Algorithm 1. Figure 2 shows several steps of Firefly Monte Carlo on a toy logistic regression model.

3 Implementation Considerations

In this section we discuss two practical matters for implementing an effective FlyMC algorithm: how to choose and compute lower bounds and how to sample the brightness variables zn. We assume a data set with N data points and a parameter vector, θ, of dimension D≪N. We also assume that it takes at least O(ND) time to evaluate the likelihoods at some θ for the whole data set and that evaluating this set of likelihoods at each iteration is the computational bottleneck for MCMC.

The goal of an effective implementation of FlyMC is to construct a Markov chain with similar convergence and mixing properties to that of regular MCMC, while only evaluating a subset of the data points on average at each iteration. If the average number of "bright" data points is M, we would like this to achieve a computational speedup of nearly N/M over regular MCMC.

3.1 Choosing a lower bound

The lower bounds, Bn(θ) of each data point's likelihood Ln(θ) should satisfy two properties. They should be relatively tight, and it should be possible to efficiently summarize a product of lower bounds ∏n Bn(θ) in time independent of N (after setup).

The tightness of the bounds is important because it determines the number of bright data points at each iteration, which determines the time it takes to evaluate the joint posterior. For a burned-in chain, the average number of bright data points, M, will be:

Therefore it is important that the bounds are tight at values of θ where the posterior puts the bulk of its mass.

The second required property is that the product of the lower bounds must be easy to compute and represent. This emerges naturally if we use scaled exponential-family lower bounds so that their product can be summarized with a set of sufficient statistics. The individual bounds Bn(θ) should be easy to compute themselves, since these are computed alongside Ln(θ) for all the bright points at each iteration.

Let's consider a concrete example. The logistic regression likelihood is

where xn ∈ RD is the set of features for the nth data point and tn ∈ {-1, 1} is its class. The logistic function has a family of scaled Gaussian lower bounds, described in Jaakkola and Jordan [1997], parameterized by ξ, the location at which the bound is tight:

where:

This is the bound shown in Fig. 1. The product of these bounds can be computed for a given θ in O(D2 ) time, provided we have precomputed the moments of the data, at a one-time setup cost of O(ND2 ).

This bound can be quite tight. For example, if we choose ξ = 1.5 the probability of a data point being bright is less than 0.02 in the region where 0.1 < Ln(θ) < 0.9. With a bit of up-front work, we can do even better than this by choosing bounds that are tight in the right places. For example, we can perform a quick optimization to find an approximate maximum a posteriori (MAP) value of θ and construct the bounds to be tight there. We explore this idea further in Section 4.

3.2 Sampling and handling the auxiliary brightness variables

The resampling of the zn variables, as shown in lines 3 to 5 of Algorithm 1, takes a step by explicitly sampling zn from its conditional distribution for a random fixed-size subset of the data. We call this approach explicit resampling and it has a clear drawback: if the fixed fraction is α (shown as RESAMPLEFRACTION in Algorithm 1), then the chain cannot have a mixing time faster than 1/α, as each data point is only visited a fraction of the time.

Nevertheless, explicit resampling works well in practice since the bottleneck for mixing is usually the exploration of θ,

1: θ0 ∼ INITIALDIST ▷ Initialize the Markov chain state. 2: for i ← 1 . . . ITERS do ▷ Iterate the Markov chain. 3: for j ← 1 . . . ⌈N × RESAMPLEFRACTION⌉ do 4: n ∼ RandInteger(1, N) ▷ Select a random data point. 5: zn ∼ Bernoulli(1 - Bn(θi-1)/Ln(θi-1)) ▷ Biased coin-flip to determine whether n is bright or dark. 6: θ′ ← θi-1 + η where η ∼ Normal(0, ϵ2 ID) ▷ Make a random walk proposal with step size ϵ. 7: u ∼ Uniform(0, 1) ▷ Draw the MH threshold. 8: if JOINTPOSTERIOR(θ′ ; {zn}Nn=1 ) JOINTPOSTERIOR(θ ; {zn}Nn=1 ) > u then ▷ Evaluate MH ratio conditioned on auxiliary variables. 9: θi ← θ′ ▷ Accept proposal. 10: else 11: θi ← θi-1 ▷ Reject proposal and keep current state. 12: function JOINTPOSTERIOR(θ ; {zn}Nn=1 ) ▷ Modified posterior that conditions on auxiliary variables. 13: P ← p(θ) × ∏Nn=1 Bn(θ) ▷ Evaluate prior and bounds. Collapse of bound product not shown. 14: for each n for which zn = 1 do ▷ Loop over bright data only. 15: P ← P × (Ln(θ)/Bn(θ) - 1) ▷ Include bound-corrected factor. 16: return P

not zn. Explicit resampling is a simple, low-overhead algorithm that is easy to vectorize for speed. The variant shown in Algorithm 1 is the simplest: data points are chosen at random, with replacement. Another variant is to deterministically choose a subset from which to Gibbs sample at each iteration. Such an approach may be appropriate for data sets which are too large to fit into memory, since we would no longer need random access to all data points.

Explicitly sampling a subset of the zn is wasteful if M ≪N, since most updates to zn will leave it unchanged. We can do better by drawing each update for zn from a pair of tunable Bernoulli proposal distributions q(z′n = 1 | zn = 0) = qd→b and q(z′n = 0 | zn = 1) = qb→d, and performing a Metropolis-Hastings accept/reject step with the true auxiliary probability p(zn | xn, θ). It is only necessary to evaluate p(zn | xn, θ) - and therefore the likelihood function - for the subset of data points which are proposed to change state. That is, if a sample from the proposal distribution sends zn = 0 to zn = 0 then it doesn't matter whether we accept or reject. If we use samples from a geometric distribution to choose the data points, it is not even necessary to explicitly sample all of the N proposals.

4 Experiments

For FlyMC to be a useful algorithm it must be able to produce effectively independent samples from posterior distributions more quickly than regular MCMC. We certainly expect it to iterate more quickly than regular MCMC since it evaluates fewer likelihoods per iteration. But we might also expect it to mix more slowly, since it has extra auxiliary variables. To see whether this trade-off works out in FlyMC's favor we need to know how much faster it iterates and how much slower it mixes.

We conducted three experiments, each with a different data set, model, and parameter-update algorithm, to give an impression of how well FlyMC can be expected to perform. In each experiment we compared FlyMC, with two choices of bound selection, to regular full-posterior MCMC. We looked at the average number of likelihoods queried at each iteration and the number of effective samples generated per iteration, accounting for autocorrelation. The results are summarized in Figure 3 and Table 1. The broad conclusion is that FlyMC offers a speedup of at least one order of magnitude compared with regular MCMC if the bounds are tuned according to a MAP-estimate of θ. In the following subsections we describe the experiments in detail.

4.1 Logistic regression

We applied FlyMC to the logistic regression task described in Welling and Teh [2011] using the Jaakkola-Jordan bounds described earlier. The task is to classify MNIST 7s and 9s, using the first 50 principal components (and one bias) as features. We used a Gaussian prior over the weights and chose the scale of that prior by evaluating performance on a held-out test set. To sample over θ, we used symmetric MetropolisHasting proposals, with step size chosen to yield an acceptance rate of 0.234 [Roberts et al., 1997], optimized for each algorithm separately. We sampled the zn using the implicit Metropolis-Hastings sampling algorithm.

We compared three different algorithms: regular MCMC, untuned FlyMC, and MAP-tuned FlyMC. For untuned FlyMC, we chose ξ = 1.5 for all data points. To compute the bounds for the MAP-tuned algorithm, we performed stochastic gradient descent optimization to find a set of weights close the the MAP value and gave each data point its own ξ to make the bounds tight at the MAP parameters: Ln(θMAP) = Bn(θMAP) for all n. For untuned FlyMC, and MAP-tuned FlyMC we used qd→b = 0.1 and qd→b = 0.01 respectively, chosen to be similar to the typical fraction of bright data points in each case.

The results are shown in Figure 3a and summarized in Table 1. On a per-iteration basis, the FlyMC algorithms mix and burn-in more slowly than regular MCMC by around a factor of two, as illustrated by the autocorrelation plots. Even on a per-likelihood basis, the na¨ıve FlyMC algorithm, with a fixed ξ, performs worse than regular MCMC, by a factor of 0.7, despite needing fewer likelihood evaluations per iteration. The MAP-tuned algorithm was much more impressive: after burn-in, it queried only 207 of the 12,2214 likelihoods per iteration on average, giving a speedup of more than 20, even taking into account the slower per-iteration mixing time. We initialized all chains with draws from the prior. Notice that the MAP-tuned algorithm performs poorly during burnin, since the bounds are less tight during this time, whereas the reverse is true for the untuned algorithm.

4.2 Softmax classification

Logistic regression can be generalized to multi-class classification problems by softmax classification. The softmax likelihood of a data point belonging to class k of K classes is

where θ is now a K × D matrix. The Jaakkola-Jordan bound does not apply to this softmax likelihood, but we can use a related bound, due to B¨ohning [1992], whose log matches the value and gradient of the log of the softmax likelihood at some particular θ, but has a tighter curvature. Murphy [2012] has the result in full in the chapter on variational inference.

We applied softmax classification to a three-class version of CIFAR-10 (airplane, automobile and bird) using 256 binary features discovered by Krizhevsky [2009] using a deep autoencoder. Once again, we used a Gaussian prior on the weights, chosen to maximize out-of-sample performance.

This time we used the Metropolis-adjusted Langevin algorithm (MALA, Roberts and Tweedie [1996]) for our parameter updates. We chose the step sizes to yield acceptance rates close to the optimal 0.57 [Roberts and Rosenthal, 1998]. Other parameters were tuned as in the logistic regression experiment.

The softmax experiment gave qualitatively similar results to the logistic regression experiment, as seen in Figure 3b and Table 1. Again, the MAP-tuned FlyMC algorithm dramatically outperformed both the lackluster untuned FlyMC and regular MCMC, offering an 11-fold speedup over the latter.

4.3 Robust sparse linear regression

Linear regression with Gaussian likelihoods yields a closedform expression for the posterior. Non-Gaussian likelihoods, however, like heavy-tailed distributions used in so-called "robust regression" do not. Our final experiment was to perform inference over robust regression weights for a very large dataset of molecular features and computed electronic properties. The data set, described by Hachmann et al. [2011, 2014] consists of 1.8 million molecules, with 57 cheminformatic features each [Olivares-Amaya et al., 2011; Amador-Bedolla et al., 2013]. The task was to predict the HOMO-LUMO energy gap, which is useful for predicting photovoltaic efficiency.

We used a student-t distribution with ν = 4 for the likelihood function and we computed a Gaussian lower bound to this by matching the value and gradient of the t distribution probability density function value at some ξ (ξ = 0 for the untuned case, ξ = θTMAP x for the MAP-tuned case). We used a sparsity-inducing Laplace prior on the weights. As before, we chose the scales of the prior and the likelihood to optimize out-of sample performance.

We performed parameter updates using slice sampling [Neal, 2003]. Note that slice sampling results in a variable number of likelihood evaluations per iteration, even for the regular MCMC algorithm. Again, we found that MAP-tuned FlyMC substantially outperformed regular MCMC, as shown in Figure 3c and Table 1.

5 Discussion

In this paper, we have presented Firefly Monte Carlo, a Markov chain Monte Carlo algorithm using subsets (minibatches) of data. FlyMC is exact in the sense that it has the true full-data posterior as its target distribution. This is achieved by introducing binary latent variables whose states represent whether a given datum is bright (used to compute the posterior) or dark (not used in posterior updates). By carefully choosing the conditional distributions of these latent variables, the true posterior is left intact under marginalization. The primary requirement for this to be efficient is that the likelihoods termmust have lower bounds that collapse in an efficient way.

There are several points that warrant additional discussion and future work. First, we recognize that useful lower bounds can be difficult to obtain for many problems. It would be helpful to produce such bounds automatically for a wider class of problems. As variational inference procedures are most often framed in terms of lower bounds on the marginal likelihood, we expect that Firefly Monte Carlo will benefit from developments in so-called "black box" variational methods [Wang and Blei, 2013; Ranganath et al., 2014]. Second, we believe we have only scratched the surface of what is possible with efficient data structures and latent-variable update schemes. For example, the MH proposals we consider here for zn have a fixed global qd→b, but clearly such a proposal should vary for each datum. Third, it is often the case that larger state spaces lead to slower MCMC mixing. We have shown empirically that the slower mixing can be more than offset by the faster per-transition computational time. In future work we hope to show that fast-mixing Markov chains on the parameter space will continue to mix fast in the Firefly auxiliary variable representation.

Firefly Monte Carlo is closely related to recent ideas in using pseudo-marginal MCMC [Andrieu and Roberts, 2009] for sampling from challenging target distributions. If we sampled each of the variables {zn} as a Bernoulli random variable with success probability 0.5, then the joint posterior we have been using becomes an unbiased estimator of the original posterior over θ, up to normalization. Running pseudomarginal MCMC using this unbiased estimator would be a special case of FlyMC: namely FlyMC with z and θ updated simultaneously with Metropolis-Hastings updates.

Acknowledgements

We thank Andrew Miller, Jasper Snoek, Michael Gelbart, Brenton Partridge, Oren Rippel and Elaine Angelino for helpful discussions. We also thank Alan Aspuru-Guzik, Johannes Hachmann, and Roberto Olivares-Amaya for the use of the OPV data set and introduction to the cheminformatic feature set. Partial funding was provided by Analog Devices (Lyric Labs), NSF IIS-1421780, and Samsung Advanced Institute of Technology.