gess
An MCMC algorithm for sampling continuous probability distributions in parallel
view repo
Many probabilistic models introduce strong dependencies between variables using a latent multivariate Gaussian distribution or a Gaussian process. We present a new Markov chain Monte Carlo algorithm for performing inference in models with multivariate Gaussian priors. Its key properties are: 1) it has simple, generic code applicable to many models, 2) it has no free parameters, 3) it works well for a variety of Gaussian process based models. These properties make our method ideal for use while model building, removing the need to spend time deriving and tuning updates for more complex algorithms.
READ FULL TEXT VIEW PDFAn MCMC algorithm for sampling continuous probability distributions in parallel
:exclamation: This is a read-only mirror of the CRAN R package repository. fastHorseshoe — The Elliptical Slice Sampler for Bayesian Horseshoe Regression. Homepage: http://jingyuhe.com/software.html
The multivariate Gaussian distribution is commonly used to specify a priori beliefs about dependencies between latent variables in probabilistic models. The parameters of such a Gaussian may be specified directly, as in graphical models and Markov random fields, or implicitly as the marginals of a Gaussian process (GP). Gaussian processes may be used to express concepts of spatial or temporal coherence, or may more generally be used to construct Bayesian kernel methods for nonparametric regression and classification. Rasmussen and Williams (2006) provide a recent review of GPs.
Inferences can only be calculated in closed form for the simplest Gaussian latent variable models. Recent work shows that posterior marginals can sometimes be well approximated with deterministic methods (Kuss and Rasmussen, 2005; Rue et al., 2009). Markov chain Monte Carlo (MCMC) methods represent joint posterior distributions with samples (e.g. Neal, 1993). MCMC can be slower but applies more generally.
In some circumstances MCMC provides good results with minimal model-specific implementation. Gibbs sampling, in particular, is frequently used to sample from probabilistic models in a straightforward way, updating one variable at a time. In models with strong dependencies among variables, including many with Gaussian priors, Gibbs sampling is known to perform poorly. Several authors have previously addressed the issue of sampling from models containing strongly correlated Gaussians, notably the recent work of Titsias et al. (2009). In this paper we provide a technique called elliptical slice sampling that is simpler and often faster than other methods, while also removing the need for preliminary tuning runs. Our method provides a drop-in replacement for MCMC samplers of Gaussian models that are currently using Gibbs or Metropolis–Hastings and we demonstrate empirical success against competing methods with several different GP-based likelihood models.
Our objective is to sample from a posterior distribution over latent variables that is proportional to the product of a multivariate Gaussian prior and a likelihood function that ties the latent variables to the observed data. We will use
to indicate the vector of latent variables that we wish to sample and denote a zero-mean Gaussian distribution with covariance
by(1) |
We also use to state that is drawn from
a distribution with the density in (1). Gaussians
with non-zero means can simply be shifted to have zero-mean with a
change of variables; an example will be given in
Section 3.3. We use to
denote the likelihood function so that our target
distribution for the MCMC sampler is
(2) |
where is the normalization constant, or the marginal likelihood, of the model.
Our starting point is a Metropolis–Hastings method introduced by Neal (1999). Given an initial state , a new state
(3) |
is proposed, where is a step-size parameter. The proposal is a sample from the prior for
and more conservative for values closer to zero. The move is accepted with probability
(4) |
otherwise the next state in the chain is a copy of .
reported that for some Gaussian process classifiers the Metropolis–Hastings method was many times faster than Gibbs sampling. The method is also simpler to implement and can immediately be applied to a much wider variety of models with Gaussian priors.
A drawback, identified by Neal (1999), is that the step-size needs to be chosen appropriately for the Markov chain to mix efficiently. This may require preliminary runs. Usually parameters of the covariance and likelihood function are also inferred from data. Different step-size parameters may be needed as the model parameters are updated. It would be desirable to automatically search over the step-size parameter, while maintaining a valid algorithm.
For a fixed auxiliary random draw, , the locus of possible proposals by varying in (3) is half of an ellipse.
(5) |
defining a full ellipse passing through the current state and the auxiliary draw . For a fixed there is an equivalent that gives the same proposal distribution in the original algorithm. However, if we can search over the step-size, the full ellipse gives a richer choice of updates for a given .
‘Slice sampling’ (Neal, 2003) provides a way to sample along a line with an adaptive step-size. Proposals are drawn from an interval or ‘bracket’ which, if too large, is shrunk automatically until an acceptable point is found. There are also ways to automatically enlarge small initial brackets. Naively applying these adaptive algorithms to select the value of in (3) or in (5) does not result in a Markov chain transition operator with the correct stationary distribution. The locus of states is defined using the current position , which upsets the reversibility and correctness of the update.
We would like to construct a valid Markov chain transition operator on the ellipse of states that uses slice sampling’s existing ability to adaptively pick step sizes. We will first intuitively construct a valid method by positing an augmented probabilistic model in which the step-size is a variable. Standard slice sampling algorithms then apply to that model. We will then adjust the algorithm for our particular setting to provide a second, slightly tidier algorithm.
Our augmented probabilistic model replaces the original latent variable with prior with
(6) |
The marginal distribution over the original latent variable is still , so the distribution over data is identical. However, we can now sample from the posterior over the new latent variables:
and use the values of deterministically derived from these samples. Our first approach applies two Monte Carlo transition operators that leave the new latent posterior invariant.
Operator 1: jointly resample the latents given the constraint that is unchanged. Because the effective variable of interest doesn’t change, the likelihood does not affect this conditional distribution, so the update is generic and easy to implement.
Operator 2: use a standard slice sampling algorithm to update the step-size given the other variables.
Input: current state , a routine that samples from , log-likelihood function . Output: a new state . When is drawn from , the marginal distribution of is also . |
Sample from :
return
The resulting algorithm is given in Figure 1. The auxiliary model construction makes the link to slice sampling explicit, which makes it easy to understand the validity of the approach. However, the algorithm can be neater and the range for slice sampling is unnatural on an ellipse. The algorithm that we will present in detail results from eliminating and and a different way of setting slice sampling’s initial proposal range. The precise connection will be given in Section 2.4. A more direct, technical proof that the equilibrium distribution of the Markov chain is the target distribution is presented in Section 2.3.
Elliptical slice sampling, our proposed algorithm is given in Figure 2, which includes the details of the slice sampler. An example run is illustrated in Figure 3(a–d). Even for high-dimensional problems, the states considered within one update lie in a two-dimensional plane. In high dimensions and are likely to have similar lengths and be an angle of apart. Therefore the ellipse will typically be fairly close to a circle, although this is not required for the validity of the algorithm.
As intended, our slice sampling approach selects a new location on the randomly generated ellipse in (5). There are no rejections: the new state is never equal to the current state unless that is the only state on the ellipse with non-zero likelihood. The algorithm proposes the angle from a bracket which is shrunk exponentially quickly until an acceptable state is found. Thus the step size is effectively adapted on each iteration for the current and .
Drawing costs , for -dimensional and general . The usual implementation of a Gaussian sampler would involve caching a (Cholesky) decomposition of , such that draws on subsequent iterations cost . For some problems with special structure drawing samples from the Gaussian prior can be cheaper.
In many models the Gaussian prior distribution captures dependencies: the observations are independent conditioned on . In these cases, computing will cost computation. As a result, drawing the random variate will be the dominant cost of the update in many high-dimensional problems. In these cases the extra cost of elliptical slice sampling over Neal’s Metropolis–Hastings algorithm will be small.
As a minor performance improvement, our implementation optionally accepts the log-likelihood of the initial state, if known from a previous update, so that it doesn’t need to be recomputed in step 2.
Input: current state , a routine that samples from , log-likelihood function . Output: a new state . When is drawn from , the marginal distribution of is also . |
Choose ellipse:
Log-likelihood threshold:
Draw an initial proposal, also defining a bracket:
if then:
Accept: return
else:
Shrink the bracket and try a new point:
if then: else:
GoTo 4.
Elliptical slice sampling considers settings of an angle variable, . Figure 2 presented the algorithm as it would be used: there is no need to index or remember the visited angles. For the purposes of analysis we will denote the ordered sequence of angles considered during the algorithm by with .
We first identify the joint distribution over a state drawn from the target distribution (
2) and the other random quantities generated by the algorithm:(7) |
where the vertical level was drawn uniformly in , that is, . The final term, , is a distribution over a random-sized set of angles, defined by the stopping rule of the algorithm.
Given the random variables in (
7) the algorithm deterministically computes positions, , accepting the first one that satisfies a likelihood constraint. More generally each angle specifies a rotation of the two a priori Gaussian variables:(8) |
For any choice of
this deterministic transformation has unit Jacobian. Any such rotation also leaves the joint prior probability invariant,
(9) |
for all , which can easily be verified by substituting values into the Gaussian form (1).
It is often useful to consider how an MCMC algorithm could make a reverse transition from the final state back to the initial state . The final state was the result of a rotation by in (8). Given an initial state of , the algorithm could generate in step 1. Then a rotation of would return back to the original pair. Moreover, the same ellipse of states is accessible and rotations of will reproduce any intermediate locations visited by the initial run of the algorithm.
In fact, the algorithm is reversible:
(10) |
the equilibrium probability of a forwards draw (7) is the same as the probability of starting at , drawing the same (possible because ), and
(11) |
resulting in the original state being returned. The reverse configuration corresponding to the result of a forwards run in Figure 3(d) is illustrated in Figure 3(e). Substituting (9) into (7) shows that ensuring that the forward and reverse angles are equally probable,
(12) |
results in the reversible property (10).
The algorithm does satisfy (12): The probability of the first angle is always . If more angles were considered before an acceptable state was found, these angles were drawn with probabilities . Whenever the bracket was shrunk in step 8, the side to shrink must have been chosen such that remained selectable as it was selected later. The reverse transition uses the same intermediate proposals, making the same rejections with the same likelihood threshold, . Because the algorithm explicitly includes the initial state, which in reverse is at , the reverse transition involves the same set of shrinking decisions as the forwards transitions. As the same brackets are sampled, the probabilities for drawing angles are the same for the forwards and reverse transitions.
The reversibility of the transition operator (10) implies that the target posterior distribution (2) is a stationary distribution of the Markov chain. Drawing from the stationary distribution and running the algorithm draws a sample from the joint auxiliary distribution (7). The deterministic transformations in (8) and (11) have unit Jacobian, so the probability density of obtaining a joint draw corresponding to is equal to the probability given by (7) for the original variables. The reversible property in (10) shows that this is the same probability as generating the variables by first generating from the target distribution and generating the remaining quantities using the algorithm. Therefore, the marginal probability of is given by the target posterior (2).
Given the first angle, the distribution over the first proposed move is . Therefore, there is a non-zero probability of transitioning to any region that has non-zero probability under the posterior. This is enough to ensure that, formally, the chain is irreducible and aperiodic (Tierney, 1994). Therefore, the Markov chain has a unique stationary distribution and repeated applications of elliptical slice sampling to an arbitrary starting point will asymptotically lead to points drawn from the target posterior distribution (2).
There is some amount of choice in how the slice sampler on the ellipse could be set up. Other methods for proposing angles could have been used, as long as they satisfied the reversible condition in . The particular algorithm proposed in Figure 2 is appealing because it is simple and has no free parameters.
The algorithm must choose the initial edges of the bracket randomly. It would be aesthetically pleasing to place the edges of the bracket at the opposite side of the ellipse to the current position, at . However this deterministic bracket placement would not be reversible and gives an invalid algorithm.
The edge of a randomly-chosen bracket could lie on the ‘slice’, the acceptable region of states. Our recommended elliptical slice sampling algorithm, Figure 2, would accept this point. The initially-presented algorithm, Figure 1, effectively randomly places the endpoints of the bracket but without checking this location for acceptability. Apart from this small change, it can be shown that the algorithms are equivalent.
In typical problems the slice will not cover the whole ellipse. For example, if is a representative sample from a posterior, often will not be. Increasing the probability of proposing points close to the current state may increase efficiency. One way to do this would be to shrink the bracket more aggressively (Skilling and MacKay, 2003). Another would be to derive a model from the auxiliary variable model (6
), but with a non-uniform distribution on
. Another way would be to randomly position an initial bracket of width less than — the code that we provide optionally allows this. However, as explained in section 2.2, for high-dimensional problems such tuning will often only give small improvements. For smaller problems we have seen it possible to improve the cpu-time efficiency of the algorithm by around two times.Another possible line of research is methods for biasing proposals away from the current state. For example the ‘over-relaxed’ methods discussed by Neal (2003) have a bias towards the opposite side of the slice from the current position. In our context it may be desirable to encourage moves close to , as these moves are independent of the previous position. These proposals are only likely to be useful when the likelihood terms are very weak, however. In the limit of sampling from the prior due to a constant likelihood, the algorithm already samples reasonably efficiently. To see this, consider the distribution over the outcome after iterations initialized at :
where and are values drawn at iteration . Only one angle is drawn per iteration when sampling from the prior, because the first proposal is always accepted. The only dependence on the initial state is the first term, the coefficient of which shrinks towards zero exponentially quickly.
A common modeling situation is that an unknown constant offset, , has been added to the entire latent vector . The resulting variable, , is still Gaussian distributed, with the constant added to every element of the covariance matrix. Neal (1999) identified that this sort of covariance will not tend to produce useful auxiliary draws . An iteration of the Markov chain can only add a nearly-constant shift to the current state. Indeed, covariances with large constant terms are generally problematic as they tend to be poorly conditioned. Instead, large offsets should be modeled and sampled as separate variables.
No algorithm can sample effectively from arbitrary distributions. As any distribution can be factored as in (2), there exist likelihoods for which elliptical slice sampling is not effective. Many Gaussian process applications have strong prior smoothness constraints and relatively weak likelihood constraints. This important regime is where we focus our empirical comparison.
Elliptical slice sampling builds on a Metropolis–Hastings (M–H) update proposed by Neal (1999). Neal reported that the original update performed moderately better than using a more obvious M–H proposal,
(13) |
and much better than Gibbs sampling for Gaussian process classification. Neal also proposed using Hybrid/Hamiltonian Monte Carlo (Duane et al., 1987; Neal, 1993), which can be very effective, but requires tuning and the implementation of gradients. We now consider some other alternatives that have similar requirements to elliptical slice sampling.
Elliptical slice sampling builds on the family of methods introduced by Neal (2003). Several of the existing slice sampling methods would also be easy to apply: they only require point-wise evaluation of the posterior up to a constant. These methods do have step-size parameters, but unlike simple Metropolis methods, typically the performance of slice samplers does not crucially rely on carefully setting free parameters.
The most popular generic slice samplers use simple univariate updates, although applying these directly to would suffer the same slow convergence problems as Gibbs sampling. While Agarwal and Gelfand (2005) have applied slice sampling for sampling parameters in Gaussian spatial process models, they assumed a linear-Gaussian observation model. For non-Gaussian data it was suggested that “there seems to be little role for slice sampling.”
Elliptical slice sampling changes all of the variables in at once, although there are potentially better ways of achieving this. An extensive search space of possibilities includes the suggestions for multivariate updates made by Neal (2003).
One simple possible slice sampling update performs a univariate update along a random line traced out by varying in (13). As the M–H method based on the line worked less well than that based on an ellipse, one might also expect a line-based slice sampler to perform less well. Intuitively, in high dimensions much of the mass of a Gaussian distribution is in a thin ellipsoidal shell. A straight line will more rapidly escape this shell than an ellipse passing through two points within it.
Titsias et al. (2009) introduced a sampling method inspired by sparse Gaussian process approximations. control variables are introduced such that the joint prior is Gaussian, and that still has marginal prior . For Gaussian process models a parametric family of joint covariances was defined, and the model is optimized so that the control variables are informative about the original variables: is made to be highly peaked. The optimization is a pre-processing step that occurs before sampling begins.
The idea is that the small number of control variables will be less strongly coupled than the original variables, and so can be moved individually more easily than the components of . A proposal involves resampling one control variable from the conditional prior and then resampling from . This move is accepted or rejected with the Metropolis–Hastings rule.
Although the method is inspired by an approximation used for large datasets, the accept/reject step uses the full model. After pre-processing it costs to evaluate a proposed change to the -dimensional vector . One ‘iteration’ in the paper consisted of an update for each control variable and so costs — roughly elliptical slice sampling updates. The control method uses fewer likelihood evaluations per iteration, although has some different minor costs associated with book-keeping of the control variables.
In some applications it may make sense to update only a subset of the latent variables at a time. This might help for computational reasons given the scaling for drawing samples of subsets of size . Titsias et al. (2009) also identified suitable subsets for local updates and then investigated sampling proposals from the conditional Gaussian prior.
In fact, local updates can be combined with any transition operator for models with Gaussian priors. If is a subset of variables to update and are the remaining variables, we can write the prior as:
(14) |
and the conditional prior is:
A change of variables allows us to express the conditional posterior as: . We can then apply elliptical slice sampling, or any alternative, to update (and thus ). Updating groups of variables according to their conditional distributions is a standard way of sampling from a joint distribution.
We performed an empirical comparison on three Gaussian process based probabilistic modeling tasks. Only a brief description of the models and methods can be given here. Full code to reproduce the results is provided as supplementary material.
Each of the models associates a dimension of the latent variable, , with an ‘input’ or ‘feature’ vector . The models in our experiments construct the covariance from the inputs using the most common method,
(15) |
the squared-exponential or “Gaussian” covariance. This covariance has “lengthscale” parameter
and an overall “signal variance”
. Other covariances may be more appropriate in many modeling situations, but our algorithm would apply unchanged.Gaussian regression: given observations of the latent variables with Gaussian noise of variance ,
(16) |
the posterior is Gaussian and so fully tractable. We use this as a simple test that the method is working correctly. Differences in performance on this task will also give some indication of performance with a simple log-concave likelihood function.
We generated ten synthetic datasets with input feature dimensions from one to ten. Each dataset was of size , with inputs drawn uniformly from a -dimensional unit hypercube and function values drawn from a Gaussian prior, , using covariance (15) with lengthscale and unit signal variance, . Noise with variance was added to generate the observations.
Gaussian process classification: a well-explored application of Gaussian processes with a non-Gaussian noise model is binary classification:
(17) |
where are the label data and
is a sigmoidal function:
for the logistic classifier; a cumulative Gaussian for the probit classifier.We ran tests on the USPS classification problem as set up by Kuss and Rasmussen (2005). We used , and the logistic likelihood.
Log Gaussian Cox process: an inhomogeneous Poisson process with a non-parametric rate can be constructed by using a shifted draw from a Gaussian process as the log-intensity function. Approximate inference can be performed by discretizing the space into bins and assuming that the log-intensity is uniform in each bin (Møller et al., 1998). Each bin contributes a Poisson likelihood:
(18) |
where the model explains the counts in bin
as drawn from a Poisson distribution with mean
. The offset to the log mean, , is the mean log-intensity of the Poisson process plus the log of the bin size.We perform inference for a Cox process model of the dates of mining disasters taken from a standard data set for testing point processes (Jarrett, 1979). The 191 events were placed into 811 bins of 50 days each. The Gaussian process parameters were fixed to and days (a third of the range of the dataset). The offset in (18) was set to , to match the empirical mean rate.
A trace of the samples’ log-likelihoods, Figure 4
, shows that elliptical slice sampling and control variables sampling have different behavior. The methods make different types of moves and only control variables sampling contains rejections. Using long runs of either method to estimate expectations under the target distribution is valid. However, sticking in a state due to many rejections can give a poor estimator as can always making small moves. It can be difficult to judge overall sampling quality from trace plots alone.
As a quantitative measure of quality we estimated the “effective number of samples” from log-likelihood traces using R-CODA (Cowles et al., 2006). Figure 6 shows these results along with computer time taken. The step size for Neal’s Metropolis method was chosen using a grid search to maximize performance. Times are for the provided implementations under Matlab v7.8 on a single bit, GHz Intel Xeon CPU. Comparing runtimes is always problematic, due to implementation-specific details. Our numbers of effective samples are primarily plotted for the same number of updates with the understanding that some correction based loosely on runtime should be applied.
The control variables approach was particularly recommended for Gaussian processes with low-dimensional input spaces. On our particular low-dimensional synthetic regression problems using control variables clearly outperforms all the other methods. On the model of mining disasters, control variable sampling has comparable performance to elliptical slice sampling with about 50% less run time. On higher-dimensional problems more control variables are required; then other methods cost less. Control variables failed to sample in high-dimensions (Figure 5). On the USPS classification problem control variables ran exceedingly slowly and we were unable to obtain any meaningful results.
Elliptical slice sampling obtained more effective samples than Neal’s M–H method with the best possible step size, although at the cost of increased run time. On the problems involving real data, elliptical slice sampling was better overall whereas M–H has more effective samples per unit time (in our implementation) on the synthetic problems. The performance differences aren’t huge; either method would work well enough.
Elliptical slice sampling takes less time than slice sampling along a straight line (line sampling involves additional prior evaluations) and usually performs better.
The slice samplers use many more likelihood evaluations than the other methods. This is partly by choice: our code can take a step-size parameter to reduce the number of likelihood evaluations (Section 2.4). On these problems the time for likelihood computations isn’t completely negligible: speedups of around may be possible by tuning elliptical slice sampling. Our default position is that ease-of-use and human time is important and that the advantage of having no free parameters should often be taken in exchange for a factor of two in runtime.
We fixed the parameters of and in our experiments to simplify the comparison. Fixing the model potentially favors the methods that have adjustable parameters. In problems were and change dramatically, a single step-size or optimized set of control variables could work very poorly.
Elliptical slice sampling is a simple generic algorithm with no tweak parameters. It performs similarly to the best possible performance of a related M–H scheme, and could be applied to a wide variety of applications in both low and high dimensions.
Thanks to Michalis Titsias for code, Sinead Williamson and Katherine Heller for a helpful discussion, and to Radford Neal, Sam Roweis, Christophe Andrieu and the reviewers for useful suggestions. RPA is funded by the Canadian Institute for Advanced Research.
Journal of Machine Learning Research
, 6:1679–1704, 2005.Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.
Journal of the Royal Statistical Society: Series B, 71(2):319–392, 2009.
Comments
There are no comments yet.