Augur: a Modeling Language for Data-Parallel Probabilistic Inference

2013·Arxiv

Abstract

Abstract

It is time-consuming and error-prone to implement inference procedures for each new probabilistic model. Probabilistic programming addresses this problem by allowing a user to specify the model and having a compiler automatically generate an inference procedure for it. For this approach to be practical, it is important to generate inference code that has reasonable performance. In this paper, we present a probabilistic programming language and compiler for Bayesian networks designed to make effective use of data-parallel architectures such as GPUs. We show that the compiler can generate data-parallel inference code scalable to thousands of GPU cores by making use of the conditional independence relationships in the Bayesian network.

1. Introduction

Machine learning, and especially probabilistic modeling, can be difficult to apply. A user needs to not only design the model, but also implement the right inference procedure. There are many different inference algorithms, most

Copyright 2014 by the author(s).

of which are conceptually complicated and difficult to implement at scale. Despite the enthusiasm that many people who practice data analysis have for machine learning, this complexity is a barrier to deployment. Any effort to simplify the use of machine learning would thus be very useful.

Probabilistic programming (?), as introduced in the BUGS project (?), is a way to simplify the application of machine learning based on Bayesian inference. The key feature of probabilistic programming is separation of concerns: the user specifies what needs to be learned by describing a probabilistic model, while the compiler automatically generates the how, that is, the inference procedure. In particular, the programmer writes code that describes a probability distribution. Using a compiler-generated inference algorithm, the programmer then samples from this distribution.

However, doing inference on probabilistic programs is computationally intensive and challenging. As a result, developing algorithms to perform inference is an active area of research. These include deterministic approximations (such as variational methods) and Monte Carlo approximations (such as MCMC algorithms). The problem is that most of these algorithms are conceptually complicated, and it is not clear, especially for non-experts, which one would work best for a given model.

To address the performance issues, our work has been driven by two observations. The first observation is that good performance starts with an appropriate inference algorithm, and selecting the right algorithm is often the hardest problem. For example, if our compiler only emits Metropolis-Hastings inference, there are models for which our programming language will be of no use, even given large amounts of computational power. We must design the compiler in such a way that we can include the latest research on inference while reusing pre-existing analyses and optimizations, or even mix inference techniques. Consequently, we have designed our compiler as a modular framework where one can add a new inference algorithm while reusing already implemented analyses and optimizations. For that purpose, our compiler uses an intermediate representation (IR) for probability distributions that serves as a target for modeling languages and as a basis for inference algorithms. We will show this IR is key to scaling the compiler and the inference to very large networks.

The second observation is if we wish to continue to bene-fit from advances in hardware we must focus on producing highly parallel inference algorithms. We claim that many MCMC inference algorithms are highly data-parallel (??) within a single Markov Chain if we take advantage of the conditional independence relationships of the input model (e.g. the assumption of i.i.d. data makes the likelihood independent across data points). Moreover, we can automatically generate good data-parallel inference with a compiler. Such inference will run very efficiently on highly parallel architectures such as Graphics Processing Units (GPUs). It is important to note that parallelism brings an interesting trade-off for performance since some inference techniques can result in less parallelism and will not scale as well.

In this paper, we present our compilation framework, named Augur. To start (Section 2), we demonstrate how to specify two popular models in our language and use them to do learning and prediction. Then, (Section 3), we describe how we support different modeling languages by embedding them into Scala using Scala’s macro system, which provides type checking and IDE support, and we describe the probabilistic IR. Next, (Section 4), we describe data-parallel versions of Metropolis-Hastings and Gibbs sampling that scale on the GPU and speed up sampling. Our compiler also includes a Metropolis-Within-Gibbs sampler but we do not detail these in this paper. Then (Section 5), we present the results of some benchmarks, which include comparisons against other implementations of inference for a regression, a Gaussian Mixture Model, and LDA. Finally (Section 6), we review the literature in probabilistic programming.

Our main results are: first, not only are some inference algorithms highly data-parallel and amenable to GPU execution, but a compiler can automatically generate such GPU implementations effectively; second, for the compiler to be able to handle large model specifications (such as LDA) it is key to use a symbolic representation of the distribution

Figure 1. Specification of the latent Dirichlet allocation model in Augur. The model specifies the probability distribution w). The keyword bayes introduces the modeling language for Bayesian networks.

rather than constructing the graphical model.

2. The Augur Language

As examples, we first present the specification of two models in Augur, Latent Dirichlet Allocation (LDA) (?) and multivariate regression. Then we show how the LDA model is used to learn the topics present in a set of documents. The supplementary material contains five examples of probabilistic models in Augur including a polynomial regression, a categorical and a Gaussian Mixture Model (GMM), a Naive Bayes classifier, and a Hidden Markov Model (HMM).

2.1. Specifying the Models

The LDA model specification is shown in Figure 1. The probability distribution is defined as a Scala object (object LDA) and is composed of two declarations. First, we declare the support of the probability distribution as a class that must be named sig. Here the support is composed of four arrays, with one each for the distribution of topics per document (theta), the distribution of words per topic (phi), the topics assigned to the words (z), and the words in the corpus (w). The support is used to store the inferred model parameters. These last two arrays are flat representations of ragged arrays, and so we do not require the documents to be of equal length.

The second declaration specifies the Bayesian network associated with LDA and makes use of our domain specific language for Bayesian networks. The DSL is marked by the bayes keyword and delimited by the following enclosing

Figure 2. Specification of a multivariate regression in Augur.

brackets. The model first declares the parameters of the model: K for the number of topics, V for the vocabulary size, M for the number of documents, and N for the array that associates each document with its size.

In the model itself, we define the hyper-parameters (values alpha and beta) for the Dirichlet distributions and draw K Dirichlet samples of dimension V for the distribution of words per topic (phi) and M Dirichlet samples of dimension K for the distribution of topics per document (theta). Then, for each word in each document, we draw a topic z from theta, and finally a word from phi based on the topic we drew for z.

The regression model in Figure 2 is defined in the same way and uses similar language features. In this example the support comprises the (x,y) data points, the weights w, the bias b, and the noise tau. The model uses an additional sum function to sum across the feature vector.

2.2. Using the model

Once a model is specified, it can be used as any other Scala object by writing standard Scala code. For instance, one may want to use the LDA model with a training corpus to learn a distribution of words per topic and then use it to learn the per-document topic distribution of a test corpus. An implementation is presented in Figure 3. First the programmer must allocate the parameter arrays which contain the inferred values. Then the signature of the model is constructed which encapsulates the parameters. The LDA.model.map command returns the MAP estimate of the parameters given the observed words.

To test the model, a new signature is constructed containing the test documents, and the previously inferred phi values. Then LDA.model.map is called again, but with both the phis and the words observed (by supplying Set("phi")). The inferred thetas for the test documents are stored in stest.theta.

3. A Modular Compilation Framework

Before we detail the architecture of our compiler, it is useful to understand how a model goes from a specification down to CUDA code running on the GPU. There are two distinct compilation phases. The first happens when the programmer compiles the program with scalac (assuming that the code from Figure 1 is in a file named LDA.scala)

The file augur.jar is the package containing our compiler. The first phase of compilation happens statically, during normal scalac compilation. In this phase, the block of code following the bayes keyword is transformed into our intermediate representation for probability distributions. The second compilation phase happens at runtime, when the programmer calls the LDA.model.map method. At that point, the IR is transformed, analyzed, and optimized, and finally, CUDA code is emitted and run.

Our framework is therefore composed of two distinct components that communicate through the IR: the front end, where domain specific languages are converted into the IR, and the back end, where the IR can be compiled down to various inference algorithms (currently MetropolisHastings, Gibbs sampling, and Metropolis-Within-Gibbs). To define a modeling language in the front end, we make use of the Scala macro system. The macro system allows us to define a set of functions (called “macros”) that will be executed by the Scala compiler on the code enclosed by the macro. We are currently focusing on Bayesian networks, but other DSLs (e.g., Markov random fields) could be added without modifications to the back end. The implementation of the macros to define the Bayesian network language is conceptually uninteresting so we omit further details. Our Bayesian network language is fairly standard, with the notable exception that it is implicitly parallel.

Separating the compilation into two distinct phases gives us many advantages. As our language is implemented using Scala’s macro system it provides automatic syntax highlighting, method name completion and code refactoring in any IDE which supports Scala. This greatly improves the usability of the DSL as no special tools need to be developed to support it. This macro system allows Augur to use Scala’s parser, semantic analyzer (e.g., to check that variables have been defined), and type checker. Also we bene-fit from the Scala compiler’s optimizations such as constant folding and dead code elimination.

Then, because the IR is compiled to CUDA code at runtime, we know the values of all the hyper-parameters and

Figure 3. Example use of the LDA. Function LDA.model.map returns a maximum a posteriori estimation. It takes as arguments the set of variables to observe (on top of the ones declared as observed in the model specification), the hyperparameters, the initial parameters, the output parameters, the number of iterations and the inference to use. The parameters are stored in LDA.sig.

the size of the dataset. This enables better optimization strategies, and also gives us key insights into how to extract parallelism (Section 4.2). For example, when compiling LDA, we know that the number of topics is much smaller than the number of documents and thus parallelizing over documents will produce more parallelism than parallelizing over topics.

Finally, we also provide a library which defines standard distributions such as Gaussian, Dirichlet, etc. In addition to these standard distributions, each model denotes its own user-defined distribution. All of these distributions are subtypes of the Dist supertype. Currently, the Dist interface provides two methods: map, which implements maximum a posteriori estimation, and sample, which returns a sequence of samples.

4. Generation of Data-Parallel Inference

When an inference procedure is invoked on a model (e.g. LDA.model.map), the IR is compiled down to CUDA inference code for that model. Informally, our IR expressions are generated from this Backus-Naur form grammar:

The goal of the IR is to make the sources of parallelism in the model more explicit and to support analysis of the probability distributions present in the model. For example, a indicates that each sub-term can be evaluated in parallel.

The use of such a symbolic representation for the model is key to scale to large networks. Indeed, as we will show in the experimental evaluation (Section 5), popular probabilistic programming language implementations such as JAGS or Stan reify the graphical model, resulting in unreasonable memory consumption for models such as LDA. A consequence of our symbolic representation is that it be- comes more difficult to discover conjugacy relationships, a

point we will come back to.

In the rest of this section, we explain how the compiler generates data-parallel samplers that exploit the conditional independence structure of the model. We will use our two examples to explain how the compiler analyzes the model and generates the inference code.

4.1. Generating data-parallel MH samplers

If the user wants to use Metropolis-Hastings inference on a model, the compiler needs to emit code for a function f that is proportional to the distribution the user wants to sample from. This function is then linked with our library implementation of Metropolis-Hastings. The function f is composed of the product of the prior and the likelihood of the model and is extracted automatically from the model specification. For example, applied to our regression example, f is defined as

which is equal to (and represented in our IR as)

In this form, the compiler knows that the distribution factorizes into a large number of terms that can be evaluated in parallel and then efficiently multiplied together; more specifically, it knows that the data is i.i.d. and that it can optimize accordingly. In this case, each (x, y) contributes to the likelihood independently, and they can be evaluated in parallel. In practice, we work in log-space, so we perform summations. The compiler can then generate the CUDA code for the evaluation of f from the IR representation. This code generation step is conceptually simple and we will not explain it further.

It is interesting to note that despite the simplicity of this parallelization the code scales reasonably well: there is a large amount of parallelism because it is roughly proportional to the number of data points; uncovering the parallelism in the code does not increase the overall quantity of computation that has to be performed; and the ratio of computation to global memory accesses is high enough to hide memory latency bottlenecks.

4.2. Generating data-parallel Gibbs samplers

Alternatively, and more interestingly, the compiler can generate a Gibbs sampler for some models. For instance, we would like to generate a Gibbs sampler for LDA, as a simple Metropolis-Hastings sampler will have a very low acceptance ratio. Currently we cannot generate a collapsed or blocked sampler, but there is interesting work related to dynamically collapsing or blocking variables (?), and we leave it to future work to extend our compiler with this capability.

To generate a Gibbs sampler, the compiler needs to Figure out how to sample from each univariate distribution. As an example, to draw as part of the th sample, the compiler needs to generate code that samples from the following distribution

As we previously explained, our compiler uses a symbolic representation of the model: the upside is that it makes it possible to scale to large networks, but the downside is that it becomes more challenging to uncover conjugacy relations and independence between variables. To accomplish this, the compiler implements an algebraic rewrite system that attempts to rewrite the above expression in terms of expressions it knows (i.e., the joint distribution of the entire model). We show a few selected rules below to give a flavor of the rewrite system.

(a)

Rule (a) states that like terms can be canceled. Rule (b) says that terms that do not depend on the variable of integration can be pulled out of the integral. Rule (c) says that we can partition a product over N-terms into two products, one where a predicate q is satisfied on the indexing variable and one where it is not. Rule (d) is a combination of the product and sum rule. Currently, the rewrite system is just comprised of rules we found useful in practice, and it is easy to extend the system to add more rewrite rules.

Going back to our example, the compiler rewrites the de-

sired expression into the one below:

In this form, it is clear that each is independent of the others after conditioning on the other random variables. As a result, they may all be sampled in parallel.

At each step, the compiler can test for a conjugacy relation. In the above form, the compiler recognizes that the are drawn from a categorical distribution and is drawn from a Dirichlet, and can exploit the fact that these are conjugate distributions. The posterior distribution for is:

where is a vector whose kth entry is the number of z associated with document m that were assigned topic k. Importantly, the compiler now knows that the drawing of each z must include a counting phase.

The case of the variables is more interesting. In this case, we want to sample from

After the application of the rewrite system to this expression, the compiler discovers that this is equal to

The key observation that the compiler takes advantage of to reach this conclusion is the fact that the z are distributed according to a categorical distribution and are used to index into the array. Therefore, they partition the set of words w into K disjoint sets , one for each topic. More concretely, the probability of words drawn from topic k can be rewritten in partitioned form using rule (c) as

This expresses the intuition that once a word’s topic is fixed, the word depends on only one of the distributions. In this form, the compiler recognizes that it should draw from

where is the count of words assigned to topic k. In general, the compiler detects patterns like the above when it notices that samples drawn from categorical distributions are being used to index into arrays.

Finally, the compiler turns to analyzing the . In this case, it will again detect that they can be sampled in parallel but it will not be able to detect a conjugacy relationship. It will then detect that the are drawn from discrete distributions, so that the univariate distribution can be calculated exactly and sampled from. In cases where the distributions are continuous, it can try to use another approximate sampling method as a subroutine for drawing that variable.

One concern with such a rewrite system is that it may fail to find a conjugacy relation if the model has a complicated structure. So far we have found our rewrite system to be robust and it can find all the usual conjugacy relations for models such as LDA, Gaussian Mixture Models or Hidden Markov Models, but it suffers from the same shortcomings as implementations of BUGS when deeper mathematics are required to discover a conjugacy relation (as would be the case for instance for a non-linear regression). In the cases where a conjugacy relation cannot be found, the compiler will (like BUGS) resort to using MetropolisHastings and therefore exploit the inherent parallelism of these algorithms.

Finally, note that the rewrite rules are applied deterministically and the process will always terminate and produce the same result. Overall, the cost of analysis is negligible compared to the sampling time for large data sets. Although the rewrite system is simple, it enables us to use a concise symbolic representation for the model and thereby scale to large networks.

4.3. Data-parallel Operations on Distributions

To produce efficient parallel code, the compiler needs to uncover parallelism, but we also need to rely on a good library of data-parallel operations for distributions. For instance, in the case of LDA, there are two steps in which we need to draw from many Dirichlet distributions in parallel. In the first case, when drawing the topic distributions for the documents, each thread can draw one of the by generating K Gamma variates and normalizing them (?). Since the number of documents is usually very large, this produces enough parallelism to make full use of the GPU’s cores.

However, this will not produce sufficient parallelism when drawing the , because the number of topics is usually small compared to the number of cores. Consequently, we use a different procedure (Algorithm 1). To generate K Dirichlet variates over V categories with concentration parameters , we first need to generate a matrix A where and then normalize each row of this matrix. For sampling the , we were effectively

launching a thread for each row. Now that the number of columns is much larger than the number of rows, we launch a thread to generate the gamma variates for each column, and then separately compute a normalizing constant for each row by multipying the matrix by an all-ones vector using CUBLAS. This is an instance where the twostage compilation procedure (Section 3) is useful, because the compiler is able to use information about the relative sizes of K and V to decide that Algorithm 1 will be more efficient than the simple scheme.

This sort of optimization is not unique to the Dirichlet distribution. For example, when generating a large number of multinormal variates by applying a linear transformation to a vector of normal variates, the strategy for extracting parallelism may change based on the number of variates to generate, the dimension of the multinormal, and the number of GPU cores. We found that to use the GPU effectively we had to develop the language in concert with the creation of a library of data-parallel operations on distributions.

4.4. Parallelism & Inference Tradeoffs

It is difficult to give a cost model for Augur programs. Traditional approaches are not necessarily appropriate for probabilistic programs because there are tradeoffs between faster sampling times and convergence which are not easy to characterize. In particular, different inference methods may affect the amount of parallelism that can be exploited in a model. For example, in the case of multivariate regression, we can use the Metropolis-Hastings sampler presented above, which lets us sample from all the weights in parallel. However, we may be better off generating a Metropolis-Within-Gibbs sampler where the weights are sampled one at a time. This reduces the amount of exploitable parallelism, but it may converge faster, and there may still be enough parallelism in each step of MetropolisHastings by evaluating the likelihood in parallel.

In the Hidden-Markov model, once again, one may try to sample the state of the Markov chain in parallel using a Metropolis-Hastings sampler just for these variables. If the HMM is small this may be a good way to make use of a GPU. Of course, for a large HMM it will be more effective to sample the states of the Markov chain in sequence, and in this case there is less parallelism to exploit.

In each of these cases, it is not clear which of these alternatives is better, and different models may perform best with different settings. Despite these difficulties, because Augur tries to parallelize operations over arrays, users can maximize the amount of parallelism in their models by structuring them so that data and parameters are stored in large, flattened arrays. In addition, as more options and inference strategies are added to Augur, users will be able to experiment with the tradeoffs of different inference methods in a way that would be too time-consuming to do manually.

5. Experimental Study

To assess the runtime performance of the inference code generated by Augur, we present the results of benchmarks for the two examples presented throughout the paper and for a Gaussian Mixture Model (GMM). More detailed information on the experiments can be found in the supplementary material.

For the multivariate regression and GMM, we compare Augur’s performance to those of two popular languages for statistical modeling, JAGS (?) and Stan (?). JAGS is an implementation of the BUGS language, and performs inference using Gibbs sampling, adaptative MetropolisHastings, and slice sampling. Stan uses No-U-Turn sampling, a variant of Hamiltonian Monte Carlo sampling. In the regression experiment, we configured Augur to use Metropolis-Hastings1, while for the GMM experiments Augur generated a Gibbs sampler.

In addition to JAGS and Stan, for the LDA benchmarks we also compare Augur to a handwritten CUDA implementation of a Gibbs sampler and a Scala implementation of a collapsed Gibbs sampler (?) from the Factorie library (?). The former gives us a reference comparison for what might be possible for a manually optimized GPU implementation, while the latter gives a baseline for a Scala implementation that does not use GPUs.

5.1. Experimental Setup

For the linear regression experiment, we used data sets from the UCI regression repository (?). The Gaussian Mixture Model experiments used two synthetic data sets, one generated from 3 clusters, the other from 4 clusters. For the LDA benchmark, we used a corpus extracted from the simple English variant of Wikipedia, with standard stopwords removed. This gives a corpus with 48556 documents, a vocabulary size of 37276 words, and approximately 3.3 million tokens. From that we sampled 1000 documents to use as a test set, removing words which appear only in the test set. To evaluate the model fit we use the log predictive probability measure (?) on the test set.

All experiments were run on a single workstation with an Intel Core i7 3770 @3.4GHz CPU, 32 GB RAM, a Samsung 840 SSD, and an NVIDIA Geforce Titan. The Titan runs on the Kepler architecture. All probability values are calculated in double precision. The CPU performance results using Factorie are calculated using a single thread, as the multi-threaded samplers are neither stable nor performant in the tested release. The GPU results use all 896 double-precision ALU cores available in the Titan2.

5.2. Results

In general, our results show that once the problem is large enough we can amortize Augur’s startup cost of model compilation to CUDA, nvcc compilation to a GPU binary, and copying the data to and from the GPU. This cost is approximately 10 seconds on average across all our experiments. After this point Augur scales to larger numbers of samples in shorter runtimes than comparable systems. More experimental detail is in the supplementary material.

Our experiments with regression show that Augur’s inference is similar to JAGS in runtime and performance, and better than Stan. This is not a surprise because the regressions have few random variables and the data sets are relatively small, and so making use of the GPU is not justified (except maybe for much larger data sets). However, the results are very different for models with latent variables where the number of variables grows with the data set.

For instance, using the GMM example, we show (Figure 8) that Augur scales better than JAGS and Stan. For a hundred thousand data points, Augur draws a thousand samples in about 3 minutes whereas JAGS needs more than 21 minutes and Stan requires more than 6 hours. Each system found the correct means and variances for the clusters, the aim here is to measure the scaling of runtime with problem size.

Results from the LDA experiment are presented in Figure 5 and use predictive probability to compare convergenceover time. We compute the predictive probability and record the time after drawing samples, for i ranging from 0 to 11 inclusive. The time is reported in seconds. It takes Augur 8.1 seconds to draw its first sample for LDA. The inference for

100

Figure 4. Evolution of runtime to draw a thousand samples for varying data set sizes for a Gaussian Mixture Model. Stan’s last data point is cropped, it took 380 minutes.

Augur is very close to that of the hand-written CUDA implementation, and much faster than the Factorie collapsed Gibbs sampler. Indeed, it takes 6.7 more hours for the collapsed LDA implementation to draw 2048 samples than it does for Augur.

We also implemented LDA in the JAGS and Stan systems but they run into scalability issues. The Stan version of LDA uses 55 gigabytes of RAM but failed to draw a second sample given a week of computation time. Unfortunately, we could not apply JAGS because it requires more than 128 gigabytes of memory. In comparison, Augur uses less than a gigabyte of memory for this experiment.

6. Related Work

To our knowledge, BUGS (?) was the first probabilistic programming language. Interestingly, most of the key concepts of probabilistic programming already appeared in the first paper to introduce BUGS (?). Since then, research in probabilistic programming languages has been focused in two directions: improving performance and scalability through better inference generation; and, increasing expressiveness and building the foundations of a universal probabilistic programming language. These two directions are useful criteria to compare probabilistic programming languages.

In terms of language expressiveness, Augur is currently limited to the specification of Bayesian networks. It is possible to extend this language (e.g., non-parametric models) or to add new modeling languages (e.g., Markov random field), but our current focus is on improving the inference generation. That is in contrast with languages like Han-

Figure 5. Evolution of the predictive probability over time for up to 2048 samples and for three implementations of LDA inference: Augur, hand-written CUDA, Factorie’s Collapsed Gibbs.

sei (?), Odds (?), Stochastic Lisp (?) and Ibal (?) which focus on increasing expressiveness, at the expense of performance. However, as Augur is embedded in the Scala programming language, we have access to the wide variety of libraries on the JVM platform and benefit from Scala tools. Augur, like Stan (?) and BUGS (??) is a domain specific probabilistic language for Bayesian networks, but it is embedded in such a way that it has a very good integration with the rest of Scala, which is crucial to software projects where data analysis is only one component of a larger artifact.

Augur is not the only system designed for scalability and performance. It is also the case of Dimple (?), Factorie (?), Infer.net (?) and Figaro (??), and the latest versions of Church (?). Dimple focuses on performance using specialized inference hardware, though it does provide an interface for CPU code. Factorie mainly focuses on undirected networks, and is a Scala library rather than a DSL (unlike all the other systems mentioned). It has multiple inference backends, and aims to be a general purpose machine learning package. Infer.net is the system most similar to Augur, in that it has a two phase compilation approach, though it is based around variational methods. A block Gibbs sampler exists but is only functional on a subset of the models. Figaro focuses on a different set of inference techniques, including techniques which use exact inference in discrete spaces (they also have Metropolis-Hastings inference). Church provides the ability to mix different inference algorithms and has some parallel capability, but it is focused on task-parallelism for multicores rather than on data-parallelism for parallel architectures. GraphLab (?) is another framework for parallel machine learning which is more focused on multiprocessor and distributed computing than on the kind of data-parallelism available on GPUs. The key difference between Augur and these other languages is the systematic generation of data-parallel algorithms for large numbers of cores (i.e., thousands) on generally available GPU hardware, and the use of a symbolic representation of the model in the compiler.

7. Conclusion

We find that it is possible to automatically generate parallel MCMC-based inference algorithms, and it is also possible to extract sufficient parallelism to saturate a modern GPU with thousands of cores. Our compiler achieves this with no extra information beyond that which is normally encoded in a graphical model description and uses a symbolic representation that allows scaling to large models (particularly for latent variable models such as LDA). It also makes it easy to run different inference algorithms and evaluate the tradeoffs between convergence and sampling time. The generated inference code is competitive in terms of sample quality with other probabilistic programming systems, and for large problems generates samples much more quickly.

A. Examples of Model Specification

We present a few examples of model specifications in Augur, covering three important topics in machine learning: regression (A.1), clustering (A.2, A.3,A.5), and classifica-tion (A.4). Our goal is to show how a few popular models can be programmed in Augur. For each of these examples, we first describe the support of the model, and then sketch the generative process, relating the most complex parts of the program to their usual mathematical notation.

A.1. Univariate polynomial regression

Our first example model is for univariate polynomial regression (Figure 6). The model’s support is composed of the array w for the weights of each mononomial, x for the domain data points and y for their image. The parameters of the model are: N, the dataset size and M, the order of the polynomial. For simplicity, this example assumes that the domain of x ranges from 0 to 2.

The generative process is: We first independently draw each of the M weights, , then draw (x, y) as follows:

For simplicity, the model is presented with many “hardwired” parameters, but it is possible to parameterize the model to control the noise level, or the domain of x.

A.2. Categorical mixture

The third example is a categorical mixture model (Figure 7). The model’s support is composed of an array z for the cluster selection, x for the data points that we draw, theta for the priors of the categorical that represents the data, and phi for the prior of the indicator variable. The parameters of the model are: N data size, K number of clusters, and V for the vocabulary size.

The generative process is: For each of the N data points we want to draw, we select a cluster z according to their distribution phi and then draw from the categorical with distribution given by theta(z).

A.3. Gaussian Mixture Model

The fourth example is a univariate Gaussian mixture model (Figure 8). The model’s support is composed of an array z for the cluster selection, x for the data points that we draw, mu for the priors over the cluster means, sigma for the priors of the cluster variances, and phi for the prior of the indicator variable. The parameters of the model are: N data size, K number of clusters.

The generative process is: For each of the N data points we want to draw, we select a cluster z according to their distribution pi and then draw from the Gaussian centered at mu(z) and of standard deviation sigma(z).

A.4. Naive bayes classifier

The fifth example is a binary naive Bayes classifier (Figure 9). The support is composed of an array c for the class and an array f for the features, pC the prior on the positive class, and pFgivenC an array for the probability of each binary feature given the class. The hyperparameters of the model are: N the number of data points, K the number of features and. The features form a 2-dimensional matrix but again the user has to “flatten” the matrix into an array.

The generative process is: First we draw the probability of an event being in one class or the other as pC. We use pC has the parameter to decide for each event in which class it falls (c). Then, for each feature, we draw the probability of the feature occurring, pFgivenC, depending on whether the event is in the class or not. Finally, we draw the features f for each event.

A.5. Hidden Markov Model

The sixth example is a hidden Markov model (Figure 10) where the observation are the result of coin flips. The support is composed of the result of the coin flips flips, the priors for each of the coins bias, the transition matrix to decide how to change coin transitionmatrix, and the states of the Markov chain that indicates which coin is

Figure 6. Specification of a univariate polynomial regression

Figure 7. Specification of a categorical mixture model

Figure 8. Specification of a Gaussian mixture model

Figure 9. Specification of a naive Bayes classifier

Figure 10. Specification of a Hidden Markov Model

being used for the flip MCstates. The two parameters of the model are the size of the data N, and the number of coins being used numberstates.

The generative process is: draw a transition matrix for the Markov chain, a bias for each of the coins, decide what coin is to be used in each state, using the transition matrix, and based on this, flip the coin that should be used for each state.

B. Simple Data-Parallel Sampling from M Dirichlet Distributions

The algorithm in 2 presents a simple way to draw from a number of Dirichlet distributions in parallel on a GPU. It works well if the number M is very large. On the contrary, it is a bottleneck if M is small or much lesser than the dimension of the Dirichlet distributions.

C. Experimental study

This section contains additional data from the benchmarks.

C.1. Multivariate Regression

In our regression experiment, we compare Augur against two other models, one implemented in Jags (11) and one in Stan (12). These models are both based upon the BMLR code developed by Kruschke (?). Each system uses the same priors and hyperparameters.

The regression experimental protocol was as follows: each dataset had 10 90%/10% train/test splits generated, and each dataset was tested using 10 different random initialisations across each of the train/test splits. Then the number of samples was varied between 100, 200, 500, 1000, 2000, 50003. This gives a total of 600 runs of each system on

Figure 13. Result of Multivariate Regression on the Concrete Compressive Strength data set.

each dataset. The presented figures average across both the random seeds and the train/test splits to produce one point per number of samples. We then plot average RMSE on the test sets against average runtime.

In figures 13, 14, 15 and 16 we present results on the Concrete compressive, winequality-red, winequality-white and Yacht Hydrodynamics datasets from the UCI repository (?). In Concrete and Yacht we present results from JAGS, Stan and Augur. On the winequality experiments we only present results from JAGS and Augur due to machine time constraints (Stan’s runtime was too high to perform suffi-ciently many experiments on larger datasets). JAGS is using a Gibbs sampler for the weights and the bias, and uses a slice sampler for the variance of the noise. Augur uses random walk Metropolis-Hastings, and Stan is using the No-U-Turn variant of Hamiltonian Monte Carlo. We can see that Augur has a startup cost of about 10 seconds, and Stan has a startup cost of about 20 seconds. After that point Augur can draw samples more quickly than both Stan and JAGS, though due to JAGS’s low startup time (¡¡1 second) it is only on large datasets with many samples that Augur provides a speedup.

The RMSEs of JAGS and Augur converge to approximately similar values, though Augur takes longer to converge (in terms of the number of samples, and total runtime) as Metropolis-Hastings is a less efficient inference algorithm for regression than a tuned Gibbs sampler. As mentioned in section 5 of the paper JAGS has a special case for working with linear regression models which alters the sampling procedure, and this feature is not currently available in Augur.

Figure 11. Multivariate Regression in Jags

Figure 12. Multivariate Regression in Stan

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 0

500 1000 20005000150 300 750 1500 3000 7500

Figure 14. Result of Multivariate Regression on the winequality-red data set.

Figure 15. Result of Multivariate Regression on the winequality-white data set.

We find that the regression results show that Augur is competitive with other systems, though linear regression models tend not to be large enough to properly exploit all the computation available in the GPU.

C.2. Gaussian Mixture Model

The Gaussian Mixture Model results presented in section 5 of the paper show how each of the three systems scale as the dataset size is increased. We sampled 100, 000 datapoints from two different mixture distributions: one with 4 gaussians centered at {-5,-1,1,5} with standard deviation {1,0.1,2,1}, and one with 3 gaussians centered at {-5, 0, 5} with standard deviations {0.1,0.1,0.1}. Each dataset had a flat mixing distribution, that is draws from each gaus- sian were equiprobable. From each dataset we subsampled

Figure 16. Result of Multivariate Regression on the Yacht Hydrodynamics data set.

Figure 17. GMM in Jags

smaller datasets using 100, 1000 and 10, 000 datapoints.

We used the GMM presented in the paper for Augur, for Stan we used the GMM listed in the modelling handbook, and for JAGS we wrote a standard GMM (shown in figure 17, based upon Augur’s. Each model used the same prior distributions and hyperparameters.

Figure 4 in the paper is from the dataset with 4 centres. In figure 18 we show the runtime of the remaining dataset with 3 centres. For computational reasons we stopped Stan’s fi-nal run after 3 hours (Stan took approx. 6 hours to complete on the first dataset). Here we can see that Augur’s runtime scales much more slowly as the dataset size is increased. JAGS remains reasonably competitive until 100, 000 data points, at which point Augur is faster by a factor of 7. Stan is also relatively competitive but scales extremely poorly as the number of datapoints is increased.

C.3. LDA

In an attempt to confirm the result presented in the paper, we present another result (figure 19) measuring the predic-

100

Figure 18. Evolution of runtime to draw a thousand samples for varying data set sizes for a Gaussian Mixture Model. Stan’s 100, 000 data point was not generated.

tive probability averaging across multiple runs using different train/test splits. In this experiment, we averaged across 10 runs with different train/test splits and present the timings with error bars. We also ran an experimment across 10 different random initializations and seeds, and all algorithms again showed robustness to the variation. We reduced the maximum number of samples to 512 as generating results for the Collapsed Gibbs sampler was proving prohibitive in terms of runtime for repeated experiments.

A third experiment (figure 20) reports on the natural logarithm of run time in milliseconds to draw 512 samples as the number of topics varies. The sparse implementation’s running time does not increase as quickly as Augur’s as the number of topics increases. As a result, it runs faster when the number of topics is large. This is because Augur’s Gibbs sampler is linear in the number of topics during the step of sampling each of the . The collapsed Gibbs sampler’s performance worsen when the number of topics is increased, as seen in our results and in the experiments in (?). Again, Augur’s generated code is on par with the hand-written CUDA implementation.

We experimented with the SparseLDA implementation which forms Factorie’s standard LDA model, but this implementation proved to be unreliable. The predictive probability measure actually decreased as more samples were drawn using the SparseLDA implementation. We are working with the developers of Factorie to investigate this problem. The SparseLDA implementation is interesting as it uses a set of LDA specific assumptions to generate a highly optimised Gibbs sampler. We found Augur to be competitive in terms of runtime when drawing more than 256 sam- ples. With smaller sample sizes there is insufficient com-

Figure 19. Average over 10 runs of the evolution of the predictive probability over time.

Figure 20. Comparison of scalability of Augur, hand-written CUDA, and Factorie’s collapsed Gibbs w.r.t the number of topics.

putation to amortize the compilation costs.

designed for accessibility and to further open science