GPU-accelerated Gibbs Sampling: a case study of the Horseshoe Probit model

Gibbs sampling is a widely used Markov Chain Monte Carlo (MCMC) method for numerically approximating integrals of interest in Bayesian statistics and other mathematical sciences. Many implementations of MCMC methods do not extend easily to parallel computing environments, as their inherently sequential nature incurs a large synchronization cost. In the case study illustrated by this paper, we show how to do Gibbs sampling in a fully data-parallel manner on a graphics processing unit (GPU), for a large class of exchangeable models that admit latent variable representations. We demonstrate our method on a Horseshoe Probit regression model, and find that our implementation scales effectively to thousands of predictors and millions of data points simultaneously.


Introduction
The Bayesian statistical paradigm has a variety of desirable properties.It accounts for the uncertainty inherent in statistical inference by producing a posterior distribution, which fundamentally contains more information about the unknown quantities of interest than a point estimate.It also propagates this uncertainty to predictive distributions, and thus does not overfit in a way that paradigms that produce only point estimates may do.Unfortunately, the computational methods required to produce a posterior distribution tend to be expensive.In particular, Markov Chain Monte Carlo (MCMC) [11,12,24] methods -the cornerstone of modern Bayesian computation -often do not scale well either with data set size or model complexity.
In this paper, we present a case study of a way to implement MCMC for a large class of Bayesian models that admit exchangeable likelihoods with latent variable representations.We do so by performing the computations on a Graphics Processing Unit (GPU), a widely available parallel processor originally designed for 3D video use cases but well suited to a variety of other tasks.In the sections that follow, we describe GPUs, characterize models in which our approach is usable, and demonstrate our method on a Horseshoe Probit model with N = 1,000,000 and p = 1,000.Standard computation with such N may easily take O(days) -our runs in O(minutes).Our approach requires no new mathematical theory -instead, we consider the systems perspective.MCMC is widely thought to be inherently sequential and unsuited to parallel environments.Furthermore, many practitioners are not aware that substantially different approaches can be needed to parallelize algorithms for use on GPUs, rather than computer clusters, due to issues such as warp divergence that will be described later.To compare with an approach well-suited to compute clusters, see Terenin et al. [33].Our contribution here is to demonstrate that Gibbs sampling on GPUs is doable for a generic class of models, and present ideas that a practitioner would need to consider in order to implement GPU Gibbs sampling, in the context of a Horseshoe Probit regression model.

Previous Work
There are a number of approaches for using GPUs for accelerating computation in Bayesian inference -some approaches are completely model specific, others are generic.The following overview is brief, emphatically not exhaustive, and presented in no particular order.
• Bayesian Mixture Models.Suchard et al. [32] review and describe GPUs, and outline a method for performing calculations needed to fit Bayesian mixture models with MCMC on a GPU.
• Hamiltonian Monte Carlo.Beam et al. [2] outline a method for fitting a Bayesian multinomial logistic regression model on a GPU with Hamiltonian Monte Carlo.
• Parallel Tempering.Mingas and Bouganis [25] describe a method for sampling from multimodal distributions with parallel tempering, using hardware acceleration in the form of a fieldprogrammable grid array (FPGA), and compare their method with GPU implementations.
• Sequential Monte Carlo. Lee et al. [17] review the architecture, programming model, and performance of GPUs and describe methods for running importance-sampling-based algorithms such as Sequential Monte Carlo on GPUs.
• State-Augmented Marginal Estimation.Seita et al. [30] propose a way to use GPUs to accelerate an annealed variant of Gibbs sampling, used for obtaining high-quality point estimates without getting stuck in local optima in discrete state spaces.
• Deep Learning.Krizhevsky [14] outlines a way to employ GPUs to accelerate algorithms used for deep learning.Though these algorithms are not Bayesian, in the computational sense of producing a distribution function as output, they are meaningful to include because the success of deep learning illustrates the sheer complexity of models that can be fitted with GPUs.

Review and Description of GPUs
A GPU can be thought of as a massively parallel co-processor.Whereas a CPU typically has a few physical cores optimized for sequential serial processing, a modern GPU has thousands of smaller cores optimized for handling multiple tasks simultaneously.GPUs operate according to the single-instruction-multiple-thread (SIMT) [22] model, a special case of single-instruction-multiple-data (SIMD) [22]: every thread in a thread block must execute the same instructions, but with respect to a different location in memory.Multiple thread blocks execute concurrently on the same GPU.Whereas a CPU is designed to quickly perform a sequence of instructions with low latency, a GPU is designed to perform data-parallel, high-throughput computing, and functions optimally when it has at least tens of thousands of simultaneous threads.These differences between the CPU and GPU paradigms lead to sharply different programming models, and significant potential performance gains in applications that are able to take advantage of the parallelism a GPU offers.We now provide a brief overview of GPUs -for a more in-depth discussion, see Sanders and Kandrot [29].
GPUs have been employed in a variety of high-performance computing tasks.While originally designed for graphics, they have also been applied to general-purpose computing (GPGPU), and have been used in a variety of areas -Sanders and Kandrot [29] list use cases in physics, signal and image processing, computational fluid dynamics, medical imaging, environmental science, deep learning, and a variety of other areas.
There are two major general-purpose programming frameworks for GPUs, CUDA [26] and OpenCL [31].CUDA is a proprietary framework for Nvidia GPUs, whereas OpenCL is an open standard framework for Nvidia and AMD GPUs, as well as multi-core CPUs, digital signal processors (DSPs), field-programmable gate arrays (FPGAs) and other processors or hardware accelerators.The frameworks are similar in principle, but scientific computing at present generally focuses on CUDA, because it has a larger user base and better availability of software libraries for primitive tasks such as matrix operations and random number generation.We focus here on Nvidia GPUs and describe them in Nvidia-specific terminology -the structure and terminology for AMD GPUs is largely analogous.
A GPU is organized hierarchically in both hardware and software.The smallest unit that performs computations is a thread processor (or core).Each thread has access to its index, and typically uses it to find its unique data point on which to perform operations.Thread processors are grouped together into streaming multiprocessors, which execute code organized into thread blocks.A thread block executes groups of 32 threads, called warps, in parallel.Each streaming multiprocessor contains shared memory that may be accessed by different warps inside the corresponding thread block.The GPU itself is composed of multiple streaming multiprocessors, which have access to global memory.
A function in CUDA that runs directly on the GPU is called a kernel, and is launched with respect to a grid of thread blocks.
This architecture has a number of implications for performance.Each thread in a warp must perform the same instruction at every point in time.This means that within a warp, if an if/else statement is encountered, all threads will first execute the if portion, wait until it completes, and then execute the else portion.This is referred to as warp divergence, and should be minimized to maintain good performance.Global memory is much slower than other types of memory, and should be accessed sparingly -ideally, only to obtain data at the beginning of a computation and to output the results at the end.As noted above, due to both the large number of parallel threads and relatively high latency -which can be hidden by scheduling multiple computations -a GPU should perform at least tens of thousands of tasks in parallel to be effectively utilized.
Thus, for optimal performance, the computation of interest should permit a high degree of what may be termed fine-grained parallelism.Data should be moved to and from the GPU as little as possible -ideally, only at the beginning at end of the computation.Threads should execute code with as few if/else statements that evaluate to multiple values within one warp as possible.They should also rely on computations performed by other threads as little as possible, and only by utilizing shared memory within their respective thread block.
Many algorithms will not satisfy these requirements without significant work.For example, difficulties [17] can arise in Sequential Monte Carlo methods that require frequent bootstrap resampling steps.The most obvious way to parallelize these is via the Poisson bootstrap [7], which requires drawing Poisson random variables with different parameters in each thread at each iteration.Unfortunately, this presents two problems.Firstly, the Poisson implementation in cuRAND requires pre-processing on the CPU for each unique parameter value -a non-starter if our goal is high performance -and hence an alternative GPU-friendly Poisson generation routine would be required.Secondly, the Poisson bootstrap produces a random number of samples at each iteration, potentially creating situations in which too few or too many samples are generated -especially so if particle degeneracy issues are a concern.Other difficulties can be present in other algorithms.It's not clear how to parallelize the No U Turn Sampler [13], a widely-used extension of Hamiltonian Monte Carlo which requires building and traversing a binary tree, in non-trivial situations where gradient evaluations alone do not expose sufficient parallelism to saturate the GPU.These algorithms are likely amenable to GPU parallelism, but would require non-trivial modifications to fit the requirements of the GPU's architecture.
We demonstrate below that the class of Gibbs samplers we consider effectively meets all of the desirable criteria mentioned above.
The CUDA Toolkit includes a number of GPU-accelerated libraries.We use cuBLAS (for performing matrix operations), cuRAND (for parallel random number generation), and cuSOLVER (for solving linear systems).Our implementation is written in Scala, a general-purpose compiled language similar to and interoperable with Java.We interface with CUDA using JCuda, a set of Java bindings for CUDA, to call GPU code.Custom CUDA kernels are written in CUDA C.

Exchangeable Models
As noted above, to run efficiently on a GPU an algorithm must demonstrate a large amount of parallelism, with minimal interaction between its parts.We now show that data-augmentationbased Gibbs samplers arising from Bayesian statistical models with exchangeable likelihoods possess sufficient parallelism.

With vector y
be an exchangeable likelihood arising from an application of de Finetti's Theorem [10], and let π(θ) be the prior for θ.Suppose that for each y i there exists a (latent) z i such that yielding the data-augmented exchangeable likelihood, which we assume to exist: Such likelihoods arise in a large class of models.In fact, complex hierarchical Bayesian models are often built by starting with a set of exchangeability assumptions about the data, and specifying a generative process that yields a likelihood expressed in terms of latent variables as above by construction.Now consider the posterior distribution The full conditional distributions for with i = 1, .., N , and can be used to construct a suitable MCMC algorithm.Often, these full conditional distributions can be recognized from known families of distributions and used to construct a Gibbs sampler [11].
In what follows, " | − " stands for "given all other components of the model." Where this is not possible, Metropolis-Hastings [24] steps can instead be performed, yielding a Metropolis-within-Gibbs sampler.Both algorithms will in the limit yield samples from the posterior distribution f (θ, z | y), which yields the target distribution f (θ | y) marginally.Now, suppose that N is large.On a CPU, these algorithms will be slow: there are N full conditionals that each need to be updated at every step of the algorithm.However, notice that f (θ | z) does not depend on the data, and that for all i (in which z −i means all of the z values except z i ), i.e., all of the variables z i are fullconditionally independent.Hence, they can be sampled all at once in parallel on a GPU.This observation forms the basis of our method.
Algorithm 2 is a standard MCMC algorithm, so no additional convergence theory for it is needed for it to be valid.Our conditional independence observation is not new -indeed, it is rather obvious and has been known for decades.However, recent advances in GPU programming, such as implementation and wide availability of linear solvers needed for models that include matrix inversion, have made running GPU-accelerated Gibbs sampling tractable.We focus our attention in what follows on implementation and performance.

Case Study: Horseshoe Probit Regression
We illustrate GPU-accelerated Gibbs sampling with the Horseshoe Probit regression model [1,6], which has a likelihood function specified (for i = 1, . . ., N ) by with z i ⊥ ⊥ z −i given θ, and the following prior: Here C + (0, 1) is a standard Cauchy distribution truncated to R + , x i and β are (1 × p) and (p × 1) vectors of data (fixed constants) and regression parameters, respectively, and y i ∈ {0, 1}.The hierarchical nature of the likelihood explicitly specifies the data augmentation in Equation ( 3).
We base our algorithm on the Gibbs sampler for probit regression in Albert and Chib [1], combining it with the hierarchical representation of the Horseshoe and corresponding Gibbs steps in Makalic and Schmidt [19].By recognizing that we arrive at the following full conditionals, suppressing the conditioning bar | henceforth for readability: where X is the (n × p) matrix whose rows are the vectors , and TN(µ, σ 2 , y i ) is a truncated normal with location µ and squared scale σ 2 , truncated to R + if y i = 1 and R − if y i = 0.
Our algorithm performs the following updates: (i) Sample z, and λ −2 followed by τ −2 in parallel.
We now describe the large degree of parallelism exploited by this sampler, and show that every calculation here is well-suited to a GPU, in the sense that everything can be performed in parallel with minimal warp divergence.

β:
i. Precompute X T X before starting the algorithm using cublasSgemm.
v. Compute Rs by solving the triangular system R −1 v = s for v using cublasStrsv.
vi. Compute X T z using cublasSgemv.
ii. Draw z by using the following sampling routine, implemented as a custom CUDA kernel.a. Compute y * i = 2y i − 1 and µ * i = µ i y * i .Note that if y i = 1 then y * i = 1, and if y i = 0 then Otherwise, draw a new Gaussian and try again.
c. Else, use the Exponential rejection sampler described in Robert [27].
i. Calculate the proposal parameter θ i in [27].The formula for θ i is modified for non-zero mean and truncation to R + , and the acceptance probability i is modified accordingly as well.
Note that there are three possible kinds of warp divergence in this routine that need to be minimized.The first kind arises because the two rejection samplers cannot be performed in parallel -this is unavoidable unless all threads in a warp have µ * < 0.47 or µ * ≥ 0.47.The second kind arises from the iterative nature of rejection sampling: a warp will finish only when all threads have acceptedthis is kept under control because both samplers are efficient, with an overall worst-case acceptance probability of around 2/3 for µ * in a neighborhood around 0.47, and probability near 1 everywhere else.The third kind arises because some z i will be truncated to R + while others are truncated to R − -this is eliminated completely by introducing the y * i .λ −2 , ν −1 , ξ: i. Draw u, a vector of IID uniforms, using curandGenerateUniform.
ii. Calculate the parameter vector θ of the exponential distributions in question, and transform u to Exp(θ) via the inversion formula −θ −1 ln(u), implemented as a custom CUDA kernel.
i. Calculate β 2 using a custom CUDA kernel.
ii. Calculate the sum in the scale parameter by expressing it as the dot product λ −2 • β 2 using cublasSdot.
iii.Draw τ −2 ∼ G(a, b) using the rejection sampler in Cheng [8], implemented as a custom CUDA kernel.This kernel involves cooperation among threads, and is described as follows.
a.Each of the 32 threads draws two uniforms, performs an accept-reject step independently, and writes its results to shared memory within the block, after which all threads synchronize.
b. Thread 1 sequentially looks through shared memory and returns once it finds an accepted value.
c.If all threads reject, new proposals are computed.This is exceedingly rare for b large: each thread accepts with probability approximately 0.88 [8], so the combined acceptance probability for the entire kernel is about 1 − 3.4 × 10 −30 .
Our current implementation is well-optimized in some ways and suboptimal in others.The update for β given the precision matrix avoids matrix inversion in favor of solving a linear system and a triangular system via a single Cholesky decomposition.This works because This is far more efficient than the following process, which we call the naive implementation because it matches up well with mathematical formulas but is a poor choice for performance: (1) calculate and store (X T X + Σ 0 ) −1 by behind-the-scenes performing an LU decomposition, (2) calculate (X T X + Σ 0 ) −1 X T z, (3) use the R function mvrnorm in the package MASS [36] to sample from a Multivariate Gaussian by behind-the-scenes performing an eigenvalue-eigenvector decomposition of (X T X + Σ 0 ) −1 .
In addition to extra matrix decompositions that degrade performance and numerical stability, the naive implementation also calculates and stores a matrix inverse, which in general should never be done when that matrix inverse will subsequently be multiplied by another matrix or vector, as solving a linear system is mathematically equivalent and computationally cheaper.Our routine performs substantially better than the naive approach, from both a speed and stability standpoint, though our approach is also not perfect from a linear algebra standpoint.It would be more numerically stable to instead use a QR decomposition with rank-1 updates in place of the Cholesky, as this would avoid the construction of X T X -we did not pursue this route because no such routine is implemented in cuSOLVER.An even more stable, albeit slower, approach would involve replacing the QR decomposition with SVD.
Most of the steps in the above routines parallelize immediately -the cuBLAS, cuRAND, and cuSOLVER library calls are highly optimized to the GPU.Our custom CUDA kernels, on the other hand, are not.The update for z could be done more efficiently by, once after burn-in, sorting the data according to the mean vector Xβ, which would minimize warp divergence as it would ensure that most warps perform the same kind of rejection sampler.Due to the additional code required to track the order of the sorted data on the GPU, we did not pursue this idea.
The update for τ −2 could be made more efficient through the use of atomic operations or shuffle instructions in the final step.We did not pursue these because τ −2 is updated in parallel with z, and its computational burden is comparatively tiny with little impact on performance.We note, however, that computing τ −2 would likely not be faster if executed by fewer than 32 threads, as the GPU would need to use a 32-thread warp to run the CUDA kernel even if some threads are doing nothing.
Our implementation mixes the cuRAND host API, which generates fixed-size arrays of random numbers in parallel, and the cuRAND device API, which allows each thread to generate its own random numbers.We use the host API as much as possible, because it is far easier to implement and manage.The device API requires management of random number generator seeds and offsets on a per-thread basis.[28] does.We highly recommend Philox for its parallel efficiency, and because it is Crush-resistant (i.e., it passes all of the tests in TestU01 [16], a suite for examining the quality of random number generators) -see Manssen et al. [20] for a review of parallel random number generators.Philox is especially well-suited to situations such as ours where random numbers cannot be pre-generated because the amount of random numbers needed per thread changes at every iteration.
To further accelerate the algorithm, we overlap computation and data transfer (IO) steps -an example of hiding latency.For example, once β has been successfully updated, Monte Carlo output is downloaded off of the GPU at the same time as updates for other variables are computed.This is done by utilizing CUDA streams, a programmatic concept in CUDA that allows multiple kernels and memory operations to execute simultaneously in parallel.
We ran our algorithm on the Horseshoe Probit regression model (8,9) with multiple data sets, varying (N, p).Performance is summarized in the following section.
6 Performance Results CPU and GPU Run Time: 10,000 Monte Carlo iterations 6.1.Synthetic Data.We first ran our algorithm on a synthetic data set with known correct answer, generated by taking Our goal was to study the GPU implementation on a non-trivial model -see Chopin and Ridgway [9] for a recent discussion on the use of binary regression models for comparing performance of Bayesian algorithms.The GPU used was a Nvidia GeForce Titan X with about 12GB RAM and single-precision floating-point performance of about 11TFLOPs.
To compare CPU and GPU run times, we ran the algorithm on both a workstation and a laptop.The laptop ran a single-threaded R implementation on an Intel Core i5 3210M CPU, and the workstation ran a multi-threaded Scala implementation on two Intel Xeon E5-2650 CPUs for a total of 16 cores.
On the laptop, we used N up to 1,000 and p up to 100.On the workstation, we used N up to 1,000,000 and p up to 1,000.Both the single-threaded R and multi-threaded Scala implementations used the same optimized precision-matrix-based sampling scheme for β as in Section 5 and likewise for all other variables.In particular, the Scala implementation was parallelized in the same way as the GPU implementation -using exchangeability together with multithreaded BLAS/LAPACK calls and random number generation.The resulting comparison showcases expected performance for interpreted and compiled languages, respectively.This provides a reasonable comparison with performance that users might expect on their laptops, and on remotely accessible workstations.
Then, to study scaling, we selected a variety of different combinations (N, p) with N up to 1,000,000 and p up to 10,000.Due to memory limitations, we did not run N = 1,000,000 and p = 10,000 simultaneously, as this would have required over 40gb RAM.Larger data sets are possible but more difficult to accommodate, because such data sets would require either multiple GPUs or streaming implementations, and we chose not to explore that possibility here -see Section 7.
We ran the Gibbs sampler for 10,000 iterations in all cases, starting from z, β set to 0 and λ, ν, τ, ξ set to 1.The Horseshoe posterior correctly identified all non-zero values in β and shrank irrelevant coefficients to zero.Mixing leaves much to be desired for large N and p, but is data-dependent, so we focus on mixing in the context of real data -see Section 6.2.
Run time can be seen in Figure 1, which shows how different values of N and p affect GPU run time, along with the 16-threaded dual-CPU workstation run time in Scala, and single-threaded laptop CPU run time in R for comparison purposes.All times referred to in this paper are clock time.
The GPU is clearly multiple orders of magnitude faster than the laptop for all data sizes we examined -for instance, with (N, p) = (10,000, 1,000), the CPU and GPU calculations took 95 minutes and 41 seconds, respectively.The GPU is also significantly faster than the more expensive dual-CPU workstation for sufficiently large problems: for (N, p) =(1,000,000, 1,000), the CPU and GPU calculations took 71 minutes and 2 minutes, respectively.Furthermore, a 10-fold increase in data set size does not necessarily result in a 10-fold increase in GPU run time.This is because larger data set sizes expose a larger degree of parallelism that the GPU is able to take advantage of.This is evident for N = 10,000 going to N = 100,000, and for N =100,000 going to N =1,000,000 -a 10-fold increase in data set size -increases run time on the GPU by a factor ranging from 1.4 to 2.1.
The GPU is faster than two 8-core CPUs even when the data size is increased 10-fold: going from (N, p) = (100,000, 1,000) on the workstation to (1,000,000, 1,000) on GPU reduces run time by a factor of 2.2.By profiling our program, we found that it spent the vast majority of time performing matrix-vector multiplications, indicating that there is little we could have done to further speed up our implementation.
6.2.Real Data.One of us (DD) was recently involved in a study in which the goal was to predict, at time t for patients hospitalized in a general medical ward, whether or not such patients would have an unplanned transfer to the intensive care unit (an adverse outcome to be avoided if possible) in a 4-hour time window starting at time (t + 8) hours.Separate models were build for four sets of principal diagnoses at admission -gastro-intestinal bleeding, serious heart problems (heart attack, congestive heart failure), pneumonia and sepsis (a total-body infection) -plus a category consisting of all other major admission diagnoses that often lead to unplanned transfers.Focusing on patients in the admission category Other, we created a data set with N = 1,295,848 one-hour hospitalization episodes, which were rendered conditionally exchangeable by a rich set of 211 clinical predictor variables.At the one-hour episode level the adverse outcome was rare (our dichotomous outcome variable had a mean of 0.0048) but highly relevant to the appropriate care path, so accurate prediction was clinically crucial.
Preliminary descriptive analysis narrowed the available independent variables down to a set of p = 141 interesting predictors -where interesting was determined according to signal-to-noise ratios in maximum-likelihood estimation -and we then used our algorithm to fit a Horseshoe Probit regression model to the resulting data set, to see how many of the interesting predictors survived the regularization process imposed by the Horseshoe prior.
From a cold start, in which the z and β were all initialized to 0 and λ, ν, τ, ξ were started at 1, our GPU algorithm produced 10,000 iterations -2,500 for burn-in and 7,500 for monitoring -in 2.94 minutes.All of the components of the β vector reached equilibrium from the cold start in about 100-200 iterations.Their mixing was slow -a typical β i had an output trace that behaved like an AR 1 (ρ) time series with first-order autocorrelation 0.88 < ρ < 0.98.Trace plots can be seen in Figure 2. To ensure reliable inference, we re-ran the algorithm for 500,000 iterations, which yielded a near-identical distribution.The Horseshoe prior shrunk about 20% of the initially-interesting predictor coefficients sharply back toward 0 (using 90% posterior intervals), leading to a more parsimonious model with a perceptible improvement in predictive performance.

Discussion
Even though MCMC methods such as Gibbs sampling are inherently sequential, our results demonstrate that the steps needed at each iteration can be sufficiently parallelizable that Gibbs sampling may be made to run orders of magnitude faster on a GPU than a CPU.We have found this to be true for our model, and we expect it to be true for other models.Exchangeable models look particularly promising, because they yield full-conditionally-independent latent variables, which can immediately be updated as a block on a GPU.Models with more complicated forms of exchangeability parallelize in more complicated ways -for an example involving Latent Dirichlet Allocation [3], see Magnusson et al. [18] and Terenin et al. [34].Many non-exchangeable models, such as Gaussian Processes and Dirichlet Process Mixtures, can likely also benefit from GPU acceleration, if there exists a sufficient degree of available parallelism in the Gibbs steps.In our implementation, updating the β vector in the Horseshoe Probit regression model is significantly faster on a GPU, in spite of the fact that β i ⊥ ⊥ β j for i = j.
The parallelism inherent in exchangeability is present even if considering other algorithms.For instance, we could have used Hamiltonian Monte Carlo for β in our Horseshoe Probit model instead of the data augmentation of Albert and Chib [1], eliminating z but retaining all other parameters.This HMC-within-Gibbs sampler would substantially improve mixing, and gradient evaluations for β could be performed in parallel -see Beam et al. [2].Note, however, that all other steps would be identical to the ones described here.We have thus chosen to focus on Gibbs sampling due to its use as a basic building block in more complicated algorithms.
The larger versions of our synthetic data sets were large enough to take up most of the GPU's memory, so (N, p) = (1,000,000, 1,000) was the biggest problem size we explored.There are several approaches to computing with data sets that are too big to fit in a GPU's memory.Multiple GPUs can be used: this will introduce some performance penalty due to synchronization requirements, but may work quite well in problems that are largely limited by computationally-intensive data-parallel steps.Bigger data sets may also be analyzed by using clusters of GPUs with Asynchronous Gibbs sampling -see Terenin et al. [33].Alternatively, data can be streamed through the GPU as it is needed, rather than loaded into memory before beginning the computation.This approach has been used for point estimation in large-scale models in software packages such as Vowpal Wabbit [15] on CPUs and Bidmach [5] on GPUs, and in our view would likely also work well for Gibbs sampling.
The Nvidia GeForce Titan X GPU that we employed is a high-end consumer GPU designed for video gaming and lightweight scientific computing use cases.It is well-suited to MCMC methods, because they are robust to numerical precision and can successfully run in floating point -see Breyer et al. [4] and Mingas and Bouganis [25].
GPUs are becoming increasingly accessible.When this was written, servers with scientific-computinggrade GPUs can be rented, from a variety of cloud computing providers, for less than 1 U.S. dollar per hour.With the recent release of the Thunderbolt 3 specification, it is also possible, through the use of an external GPU case such as the Akitio Node, to connect a desktop-grade GPU to a laptop computer using a simple cable.
Software is currently the most significant barrier to widespread adoption of GPU acceleration for Bayesian inference.The functions provided by CUDA libraries are low-level, and present a high degree of difficulty to those not already familiar with high-performance computing libraries such as BLAS and LAPACK.New frameworks are needed to bring modern programmatic concepts into the GPU software stack.These would parallel recent advances in programmatic concepts and frameworks for compute clusters, such as Akka [35] and Apache Spark [38], which made computation in the cluster environment much more user-friendly through the implementation of modern concepts such as actor models, functional programming, and lazy evaluation.Similar work is needed in simplifying GPU debugging.We hope that in a few years it will be just as easy to write a Gibbs sampler on a GPU as it currently is on a CPU in a high-level language such as R.
Another difficulty in CUDA is lack of random number generation routines for less-common distributions.In the Horseshoe Probit regression model, we needed to implement custom CUDA kernels for sampling Exponential, Gamma, and Truncated Normal random variables.We also needed to implement a CUDA kernel for squaring arrays of numbers -while trivial to write, this resulted in yet another block of code that needed to be written, debugged, and compiled before we could run our algorithm.These issues are even more pronounced in OpenCL, where there are no widely available libraries for many matrix and vector routines.Future work will likely remove these difficulties.
Though our initial results on the Horseshoe Probit regression model are promising, further work is needed to study Bayesian computation on GPUs.Most parallelization strategies are unexplored, and the relative performance of different methods is unknown.We believe that our work in this case study demonstrates that fast generic methods for large classes of Bayesian inferential and predictive problems are possible, and that GPU-accelerated Bayesian computation is worth substantial further exploration.

Figure 1 .
Figure 1.Comparison of GPU, workstation (dual 8-core CPUs), and laptop (single-core CPU) run time for 10,000 iterations of the Horseshoe Probit Gibbs sampler with synthetic data.
The simplest approach is to have separate random number generators in each thread, seeded to ensure that random numbers do not overlap.The main drawback of this approach is that it requires each thread to store its random number generator state in global memory.Instead, our code stores only one random generator state, and works by (1) copying the random number generator state from global memory to local memory, (2) skipping forward k iterations in the random number generator, to ensure that random numbers generated by different threads do not overlap,(3)generating random numbers, and (4) once the computation is complete, writing the last thread's random number generator state back to global memory.To be efficient, this requires a random number generator capable of skipping ahead an arbitrarily large number of iterations in constant time without causing warp divergence.The cuRAND default XORWOW random number generator [21] does not have constant-time skip ahead, and neither does the parallel version of the popular Mersenne Twister [23], but Salmon et al.'s counter-based Philox random number generator