1 Introduction
Probabilistic programming systems (PPS) allow probabilistic models to be represented in the form of a generative model and statements for conditioning on data (Carpenter et al., 2015; Goodman et al., 2008; Goodman and Stuhlmüller, 2014; Mansinghka et al., 2014; Minka et al., 2010; Wood et al., 2014). Their core philosophy is to decouple model specification and inference, the former corresponding to the userspecified program code and the latter to an inference engine capable of operating on arbitrary programs. Removing the need for users to write inference algorithms significantly reduces the burden of developing new models and makes effective statistical methods accessible to nonexperts.
Although significant progress has been made on the problem of general purpose inference of program variables, less attention has been given to their optimization
. Optimization is an essential tool for effective machine learning, necessary when the user requires a single estimate. It also often forms a tractable alternative when full inference is infeasible
(Murphy, 2012). Moreover, coincident optimization and inference is often required, corresponding to a marginal maximum a posteriori (MMAP) setting where one wishes to maximize some variables, while marginalizing out others. Examples of MMAP problems include hyperparameter optimization, expectation maximization, and policy search
(van de Meent et al., 2016).In this paper we develop the first system that extends probabilistic programming (PP) to this more general MMAP framework, wherein the user specifies a model in the same manner as existing systems, but then selects some subset of the sampled variables in the program to be optimized, with the rest marginalized out using existing inference algorithms. The optimization query we introduce can be implemented and utilized in any PPS that supports an inference method returning a marginal likelihood estimate. This framework increases the scope of models that can be expressed in PPS and gives additional flexibility in the outputs a user can request from the program.
MMAP estimation is difficult as it corresponds to the optimization of an intractable integral, such that the optimization target is expensive to evaluate and gives noisy results. Current PPS inference engines are typically unsuited to such settings. We therefore introduce BOPP^{1}^{1}1Code available at http://www.github.com/probprog/bopp/ (Bayesian optimization for probabilistic programs) which couples existing inference algorithms from PPS, like Anglican (Wood et al., 2014), with a new Gaussian process (GP) (Rasmussen and Williams, 2006) based Bayesian optimization (BO) package (Gutmann and Corander, 2016; Jones et al., 1998; Osborne et al., 2009; Shahriari et al., 2016a).
To demonstrate the functionality provided by BOPP, we consider an example application of engineering design. Engineering design relies extensively on simulations which typically have two things in common: the desire of the user to find a single best design and an uncertainty in the environment in which the designed component will live. Even when these simulations are deterministic, this is an approximation to a truly stochastic world. By expressing the utility of a particular designenvironment combination using an approximate Bayesian computation (ABC) likelihood (Csilléry et al., 2010), one can pose this as a MMAP problem, optimizing the design while marginalizing out the environmental uncertainty.
). The bottom right plot shows convergence of the evidence of the respective model, giving the median and 25/75% quartiles.
Figure 1 illustrates how BOPP can be applied to engineering design, taking the example of optimizing the distribution of power between radiators in a house so as to homogenize the temperature, while marginalizing out possible weather conditions and subject to a total energy budget. The probabilistic program shown in Figure 2 allows us to define a prior over the uncertain weather, while conditioning on the output of a deterministic simulator (here Energy2D (Xie, 2012)a finite element package for heat transfer) using an ABC likelihood. BOPP now allows the required coincident inference and optimization to be carried out automatically, directly returning increasingly optimal configurations.
BO is an attractive choice for the required optimization in MMAP as it is typically efficient in the number of target evaluations, operates on nondifferentiable targets, and incorporates noise in the target function evaluations. However, applying BO to probabilistic programs presents challenges, such as the need to give robust performance on a wide range of problems with varying scaling and potentially unbounded support. Furthermore, the target program may contain unknown constraints, implicitly defined by the generative model, and variables whose type is unknown (i.e. they may be continuous or discrete).
On the other hand, the availability of the target source code in a PPS presents opportunities to overcome these issues and go beyond what can be done with existing BO packages. BOPP exploits the source code in a number of ways, such as optimizing the acquisition function using the original generative model to ensure the solution satisfies the implicit constaints, performing adaptive domain scaling to ensure that GP kernel hyperparameters can be set according to problemindependent hyperpriors, and defining an adaptive nonstationary mean function to support unbounded BO.
Together, these innovations mean that BOPP can be run in a manner that is fully blackbox from the user’s perspective, requiring only the identification of the target variables relative to current syntax for operating on arbitrary programs. We further show that BOPP is competitive with existing BO engines for direct optimization on common benchmarks problems that do not require marginalization.
2 Background
2.1 Probabilistic Programming
Probabilistic programming systems allow users to define probabilistic models using a domainspecific programming language. A probabilistic program implicitly defines a distribution on random variables, whilst the system backend implements generalpurpose inference methods.
PPS such as Infer.Net (Minka et al., 2010) and Stan (Carpenter et al., 2015) can be thought of as defining graphical models or factor graphs. Our focus will instead be on systems such as Church (Goodman et al., 2008), Venture (Mansinghka et al., 2014), WebPPL (Goodman and Stuhlmüller, 2014), and Anglican (Wood et al., 2014), which employ a generalpurpose programming language for model specification. In these systems, the set of random variables is dynamically typed, such that it is possible to write programs in which this set differs from execution to execution. This allows an unspecified number of random variables and incorporation of arbitrary black box deterministic functions, such as was exploited by the simulate function in Figure 2. The price for this expressivity is that inference methods must be formulated in such a manner that they are applicable to models where the density function is intractable and can only be evaluated during forwards simulation of the program.
One such general purpose system, Anglican, will be used as a reference in this paper. In Anglican, models are defined using the inference macro defquery. These models, which we refer to as queries (Goodman et al., 2008)
, specify a joint distribution
over data and variables . Inference on the model is performed using the macro doquery, which produces a sequence of approximate samples from the conditional distribution and, for importance sampling based inference algorithms (e.g. sequential Monte Carlo), a marginal likelihood estimate .Random variables in an Anglican program are specified using sample statements, which can be thought of as terms in the prior. Conditioning is specified using observe statements which can be thought of as likelihood terms. Outputs of the program, taking the form of posterior samples, are indicated by the return values. There is a finite set of sample and observe statements in a program source code, but the number of times each statement is called can vary between executions. We refer the reader to http://www.robots.ox.ac.uk/~fwood/anglican/ for more details.
2.2 Bayesian Optimization
Consider an arbitrary blackbox target function that can be evaluated for an arbitrary point to produce, potentially noisy, outputs . BO (Jones et al., 1998; Osborne et al., 2009) aims to find the global maximum
(1) 
The key idea of BO is to place a prior on that expresses belief about the space of functions within which might live. When the function is evaluated, the resultant information is incorporated by conditioning upon the observed data to give a posterior over functions. This allows estimation of the expected value and uncertainty in for all . From this, an acquisition function is defined, which assigns an expected utility to evaluating at particular , based on the tradeoff between exploration and exploitation in finding the maximum. When direct evaluation of is expensive, the acquisition function constitutes a cheaper to evaluate substitute, which is optimized to ascertain the next point at which the target function should be evaluated in a sequential fashion. By interleaving optimization of the acquisition function, evaluating
at the suggested point, and updating the surrogate, BO forms a global optimization algorithm that is typically very efficient in the required number of function evaluations, whilst naturally dealing with noise in the outputs. Although alternatives such as random forests
(Bergstra et al., 2011; Hutter et al., 2011)(Snoek et al., 2015) exist, the most common prior used for is a GP (Rasmussen and Williams, 2006). For further information on BO we refer the reader to the recent review by Shahriari et al Shahriari et al. (2016b).2.3 Gaussian Processes
Informally one can think of a Gaussian Process (GP) (Rasmussen and Williams, 2006) as being a nonparametric distribution over functions which is fully specified by a mean function and covariance function , the latter of which must be a bounded and reproducing kernel. We can describe a function as being distributed according to a GP:
(2) 
which by definition means that the functional evaluations realized at any finite number of sample points is distributed according to a multivariate Gaussian. Note that the inputs to and need not be numeric and as such a GP can be defined over anything for which kernel can be defined.
An important property of a GP is that it is conjugate with a Gaussian likelihood. Consider pairs of inputoutput data points , , and the separable likelihood function
(3) 
where is an observation noise. Using a GP prior leads to an analytic GP posterior
(4)  
(5) 
and Gaussian predictive distribution
(6) 
where we have used the shorthand and similarly for , and .
3 Problem Formulation
Given a program defining the joint density with fixed , our aim is to optimize with respect to a subset of the variables whilst marginalizing out latent variables
(7) 
To provide syntax to differentiate between and , we introduce a new query macro defopt. The syntax of defopt is identical to defquery except that it has an additional input identifying the variables to be optimized. To allow for the interleaving of inference and optimization required in MMAP estimation, we further introduce doopt, which, analogous to doquery, returns a lazy sequence where are the program outputs associated with and each is an estimate of the corresponding log marginal (see Section 4.2). The sequence is defined such that, at any time, corresponds to the point expected to be most optimal of those evaluated so far and allows both inference and optimization to be carried out online.
Although no restrictions are placed on , it is necessary to place some restrictions on how programs use the optimization variables specified by the optimization argument list of defopt. First, each optimization variable must be bound to a value directly by a sample statement with fixed measuretype distribution argument. This avoids change of variable complications arising from nonlinear deterministic mappings. Second, in order for the optimization to be well defined, the program must be written such that any possible execution trace binds each optimization variable exactly once. Finally, although any may be lexically multiply bound, it must have the same base measure in all possible execution traces, because, for instance, if the base measure of a were to change from Lebesgue to counting, the notion of optimality would no longer admit a conventional interpretation. Note that although the transformation implementations shown in Figure 3 do not contain runtime exception generators that disallow continued execution of programs that violate these constraints, those actually implemented in the BOPP system do.
4 Bayesian Program Optimization
In addition to the syntax introduced in the previous section, there are five main components to BOPP:

A program transformation, qqmarg, allowing estimation of the evidence at a fixed .

A highperformance, GP based, BO implementation for actively sampling .

A program transformation, qqprior, used for automatic and adaptive domain scaling, such that a problemindependent hyperprior can be placed over the GP hyperparameters.

An adaptive nonstationary mean function to support unbounded optimization.

A program transformation, qqacq, and annealing maximum likelihood estimation method to optimize the acquisition function subject the implicit constraints imposed by the generative model.
Together these allow BOPP to perform online MMAP estimation for arbitrary programs in a manner that is blackbox from the user’s perspective  requiring only the definition of the target program in the same way as existing PPS and identifying which variables to optimize. The BO component of BOPP is both probabilistic programming and language independent, and is provided as a standalone package.^{2}^{2}2Code available at http://www.github.com/probprog/deodorant/ It requires as input only a target function, a sampler to establish rough input scaling, and a problem specific optimizer for the acquisition function that imposes the problem constraints.
Figure 3 provides a high level overview of the algorithm invoked when doopt is called on a query q that defines a distribution . We wish to optimize whilst marginalizing out and , as indicated by the the second input to q. In summary, BOPP performs iterative optimization in 5 steps

Step 1 (blue arrows) generates unweighted samples from the transformed prior program qprior (top center), constructed by removing all conditioning. This initializes the domain scaling for .

Step 2 (red arrows) evaluates the marginal at a small number of the generated by performing inference on the marginal program qmarg (middle centre), which returns samples from the distribution along with an estimate of . The evaluated points (middle right) provide an initial domain scaling of the outputs and starting points for the BO surrogate.

Step 3 (black arrow) fits a mixture of GPs posterior Rasmussen and Williams (2006) to the scaled data (bottom centre) using a problem independent hyperprior. The solid blue line and shaded area show the posterior mean and standard deviations respectively. The new estimate of the optimum is the value for which the mean estimate is largest, with equal to the corresponding mean value.

Step 4 (purple arrows) constructs an acquisition function (bottom left) using the GP posterior. This is optimized, giving the next point to evaluate , by performing annealed importance sampling on a transformed program qacq (middle left) in which all observe statements are removed and replaced with a single observe assigning probability to the execution.

Step 5 (green arrow) evaluates using qmarg and continues to step 3.
4.1 Program Transformation to Generate the Target
Consider the defopt query q in Figure 3, the body of which defines the joint distribution . Calculating (7) (defining ) using a standard optimization scheme presents two issues:
is a random variable within the program rather than something we control and its probability distribution is only defined conditioned on
.We deal with both these issues simultaneously using a program transformation similar to the disintegration transformation in Hakaru (Zinkov and Shan, 2016). Our marginal transformation returns a new query object, qmarg as shown in Figure 3, that defines the same joint distribution on program variables and inputs, but now accepts the value for as an input. This is done by replacing all sample statements associated with with equivalent observe< statements, taking as the observed value, where observe< is identical to observe except that it returns the observed value. As both sample and observe operate on the same variable type  a distribution object  this transformation can always be made, while the identical returns of sample and observe< trivially ensures validity of the transformed program.
4.2 Bayesian Optimization of the Marginal
The target function for our BO scheme is , noting for any . The log is taken because GPs have unbounded support, while is always positive, and because we expect variations over many orders of magnitude. PPS with importance sampling based inference engines, e.g. sequential Monte Carlo (Wood et al., 2014) or the particle cascade (Paige et al., 2014), can return noisy estimates of this target given the transformed program qmarg.
Our BO scheme uses a GP prior and a Gaussian likelihood. Though the rationale for the latter is predominantly computational, giving an analytic posterior, there are also theoretical results suggesting that this choice is appropriate (Bérard et al., 2014). We use as a default covariance function a combination of a Matérn3/2 and Matérn5/2 kernel. Specifically, let be the dimensionality of and define
(8a)  
(8b) 
where indexes a dimension of and and are dimension specific length scale hyperparameters. Our prior covariance function is now given by
(9) 
where and represent signal standard deviations for the two respective kernels. The full set of GP hyperparameters is defined by . A key feature of this kernel is that it is only once differentiable and therefore makes relatively weak assumptions about the smoothness of . The ability to include branching in a probabilistic program means that, in some cases, an even less smooth kernel than (9) might be preferable. However, there is clear a tradeoff between generality of the associated reproducing kernel Hilbert space and modelling power.
As noted by (Snoek et al., 2012), the performance of BO using a single GP posterior is heavily influenced by the choice of these hyperparameters. We therefore exploit the automated domain scaling introduced in Section 4.3 to define a problem independent hyperprior and perform inference to give a mixture of GPs posterior. Details on this hyperprior are given in Appendix B.
Inference over is performed using Hamiltonian Monte Carlo (HMC) (Duane et al., 1987), giving an unweighted mixture of GPs. Each term in this mixture has an analytic distribution fully specified by its mean function and covariance function , where indexes the BO iteration and the hyperparameter sample. HMC was chosen because of the availability of analytic derivatives of the GP log marginal likelihoods. As we found that the performance of HMC was often poor unless a good initialization point was used, BOPP runs a small number of independent chains and allocates part of the computational budget to their initialization using a LBFGS optimizer (Broyden, 1970).
The inferred posterior is first used to estimate which of the previously evaluated is the most optimal, by taking the point with highest expected value , . This completes the definition of the output sequence returned by the doopt macro. Note that as the posterior updates globally with each new observation, the relative estimated optimality of previously evaluated points changes at each iteration. Secondly it is used to define the acquisition function , for which we take the expected improvement (Snoek et al., 2012), defining and ,
(10) 
where and
represent the pdf and cdf of a unit normal distribution respectively. We note that more powerful, but more involved, acquisition functions, e.g.
(HernándezLobato et al., 2014), could be used instead.4.3 Automatic and Adaptive Domain Scaling
Domain scaling, by mapping to a common space, is crucial for BOPP to operate in the required blackbox fashion as it allows a general purpose and problem independent hyperprior to be placed on the GP hyperparameters. BOPP therefore employs an affine scaling to a hypercube for both the inputs and outputs of the GP. To initialize scaling for the input variables, we sample directly from the generative model defined by the program. This is achieved using a second transformed program, qprior, which removes all conditioning, i.e. observe statements, and returns . This transformation also introduces code to terminate execution of the query once all are sampled, in order to avoid unnecessary computation. As observe statements return nil, this transformation trivially preserves the generative model of the program, but the probability of the execution changes. Simulating from the generative model does not require inference or calling potentially expensive likelihood functions and is therefore computationally inexpensive. By running inference on qmarg given a small number of these samples as arguments, a rough initial characterization of output scaling can also be achieved. If points are observed that fall outside the hypercube under the initial scaling, the domain scaling is appropriately updated^{3}^{3}3An important exception is that the output mapping to the bottom of the hypercube remains fixed such that low likelihood new points are not incorporated. This ensures stability when considering unbounded problems. so that the target for the GP remains the hypercube.
4.4 Unbounded Bayesian Optimization via NonStationary Mean Function Adaptation
Unlike standard BO implementations, BOPP is not provided with external constraints and we therefore develop a scheme for operating on targets with potentially unbounded support. Our method exploits the knowledge that the target function is a probability density, implying that the area that must be searched in practice to find the optimum is finite, by defining a nonstationary prior mean function. This takes the form of a bump function that is constant within a region of interest, but decays rapidly outside. Specifically we define this bump function in the transformed space as
(11) 
where is the radius from the origin, is the maximum radius of any point generated in the initial scaling or subsequent evaluations, and is a parameter set to by default. Consequently, the acquisition function also decays and new points are never suggested arbitrarily far away. Adaptation of the scaling will automatically update this mean function appropriately, learning a region of interest that matches that of the true problem, without complicating the optimization by overextending this region. We note that our method shares similarity with the recent work of Shahriari et al (Shahriari et al., 2016a), but overcomes the sensitivity of their method upon a userspecified bounding box representing soft constraints, by initializing automatically and adapting as more data is observed.
4.5 Optimizing the Acquisition Function
Optimizing the acquisition function for BOPP presents the issue that the query contains implicit constraints that are unknown to the surrogate function. The problem of unknown constraints has been previously covered in the literature (Gardner et al., 2014; HernándezLobato et al., 2016) by assuming that constraints take the form of a blackbox function which is modeled with a second surrogate function and must be evaluated in guessandcheck strategy to establish whether a point is valid. Along with the potentially significant expense such a method incurs, this approach is inappropriate for equality constraints or when the target variables are potentially discrete. For example, the Dirichlet distribution in Figure 2 introduces an equality constraint on powers, namely that its components must sum to .
We therefore take an alternative approach based on directly using the program to optimize the acquisition function. To do so we consider a transformed program qacq that is identical to qprior (see Section 4.3), but adds an additional observe statement that assigns a weight to the execution. By setting to the acquisition function, the maximum likelihood corresponds to the optimum of the acquisition function subject to the implicit program constraints. We obtain a maximum likelihood estimate for qacq using a variant of annealed importance sampling (Neal, 2001) in which lightweight Metropolis Hastings (LMH) (Wingate et al., 2011) with local randomwalk moves is used as the base transition kernel.
5 Experiments
We first demonstrate the ability of BOPP to carry out unbounded optimization using a 1D problem with a significant priorposterior mismatch as shown in Figure 4. It shows BOPP adapting to the target and effectively establishing a maxima in the presence of multiple modes. After 20 evaluations the acquisitions begin to explore the right mode, after 50 both modes have been fully uncovered.
5.1 Classic Optimization Benchmarks
Next we compare BOPP to the prominent BO packages SMAC Hutter et al. (2011), Spearmint Snoek et al. (2012) and TPE Bergstra et al. (2011) on a number of classical benchmarks as shown in Figure 5. These results demonstrate that BOPP provides substantial advantages over these systems when used simply as an optimizer on both continuous and discrete optimization problems. In particular, it offers a large advantage over SMAC and TPE on the continuous problems (Branin and Hartmann), due to using a more powerful surrogate, and over Spearmint on the others due to not needing to make approximations to deal with discrete problems.
5.2 Marginal Maximum a Posteriori Estimation Problems
We now demonstrate application of BOPP on a number of MMAP problems. Comparisons here are more difficult due to the dearth of existing alternatives for PPS. In particular, simply running inference on the original query does not return estimates for . We consider the possible alternative of using our conditional code transformation to design a particle marginal Metropolis Hastings (PMMH, Andrieu et al. (2010)) sampler which operates in a similar fashion to BOPP except that new are chosen using a MH step instead of actively sampling with BO. For these MH steps we consider both LMH (Wingate et al., 2011) with proposals from the prior and the randomwalk MH (RMH) variant introduced in Section 4.5.
5.2.1 Hyperparameter Optimization for Gaussian Mixture Model
We start with an illustrative case study of optimizing the hyperparameters in a multivariate Gaussian mixture model. We consider a Bayesian formulation with a symmetric Dirichlet prior on the mixture weights and a GaussianinverseWishart prior on the likelihood parameters:
(12)  
(13)  
(14)  
(15) 
Anglican code for this model is shown in Figure 4. Anglican provides stateful objects, which are referred to as random processes, to represent the predictive distributions for the cluster assignments and the observations assigned to each cluster
(16)  
(17) 
In this collapsed representation marginalization over the model parameters , , and is performed analytically. Using the Iris dataset, a standard benchmark for mixture models that contains 150 labeled examples with 4 realvalued features, we optimize the marginal with respect to the subset of the parameters and under uniform priors over a fixed interval. For this model, BOPP aims to maximize
(18) 
Figure 7 shows GP regressions on the evidence after different numbers of the SMC evaluations have been performed on the model. This demonstrates how the GP surrogate used by BO builds up a model of the target, used to both estimate the expected value of for a particular and actively sample the at which to undertake inference.
5.2.2 Extended Kalman Filter for the Pickover Chaotic Attractor
). Solid line shows median over 100 runs, whilst the shaded region the 25/75% quantiles.
We next consider the case of learning the dynamics parameters of a chaotic attractor. Chaotic attractors present an interesting case for tracking problems as, although their underlying dynamics are strictly deterministic with bounded trajectories, neighbouring trajectories diverge exponentially^{4}^{4}4It is beyond the scope of this paper to properly introduce chaotic systems. We refer the reader to Devaney et al. (1989) for an introduction.
. Therefore regardless of the available precision, a trajectory cannot be indefinitely extrapolated to within a given accuracy and probabilistic methods such as the extended Kalman filter must be incorporated
(Fujii, 2013; Ruan et al., 2003). From an empirical perspective, this forms a challenging optimization problem as the target transpires to be multimodal, has variations at different length scales, and has local minima close to the global maximum.Suppose we observe a noisy signal in some dimensional observation space were each observation has a lower dimensional latent parameter whose dynamics correspond to a chaotic attractor of known type, but with unknown parameters. Our aim will be to find the MMAP values for the dynamics parameters , marginalizing out the latent states. The established parameters can then be used for forward simulation or tracking.
To carry out the required MMAP estimation, we apply BOPP to the extended Kalman smoother
(19)  
(20)  
(21) 
where
is the identity matrix,
is a known matrix, is the expected starting position, and and are all scalars which are assumed to be known. The transition function(22a)  
(22b)  
(22c) 
corresponds to a Pickover attractor (Pickover, 1995) with unknown parameters which we wish to optimize. Note that and will give the same behaviour.
Synthetic data was generated for time steps using the parameters of , , , , a fixed matrix where and each column was randomly drawn from a symmetric Dirichlet distribution with parameter , and ground truth transition parameters of and (note that the true global optimum for finite data need not be exactly equal to this).
MMAP estimation was performed on this data using the same model and parameters, with the exceptions of , and . The prior on was set to a uniform in over a bounded region such that
(23) 
The changes and were further made to reflect the starting point of the latent state being unknown. For this problem, BOPP aims to maximize
(24) 
Inference on the transformed marginal query was carried out using SMC with 500 particles. Convergence results are given in Figure 8 showing that BOPP comfortably outperforms the PMMH variants, while Figure 9 shows the simulated attractors generated from the dynamics parameters output by various iterations of a particular run of BOPP.
5.2.3 Hidden Markov Model with Unknown Number of States
We finally consider a hidden Markov model (HMM) with an unknown number of states. This example demonstrates how BOPP can be applied to models which conceptually have an unknown number of variables, by generating all possible variables that might be needed, but then leaving some variables unused for some execution traces. This avoids problems of varying base measures so that the MMAP problem is well defined and provides a function with a fixed number of inputs as required by the BO scheme. From the BO perspective, the target function is simply constant for variations in an unused variable.
HMMs are Markovian state space models with discrete latent variables. Each latent state is defined conditionally on through a set of discrete transition probabilities, whilst each output is considered to be generated i.i.d. given . We consider the following HMM, in which the number of states , is also a random variable:
(25)  
(26)  
(27)  
(28)  
(29)  
(30)  
(31)  
(32) 
Our experiment is based on applying BOPP to the above model to do MMAP estimation with a single synthetic dataset, generated using and .
We use BOPP to optimize both the number of states and the stickbreaking parameters , with full inference performed on the other parameters. BOPP therefore aims to maximize
(33) 
As with the chaotic Kalman filter example, we compare to two PMMH variants using the same code transformations. The results, given in Figure 10, again show that BOPP outperforms these PMMH alternatives.
6 Discussion and Future Work
We have introduced a new method for carrying out MMAP estimation of probabilistic program variables using Bayesian optimization, representing the first unified framework for optimization and inference of probabilistic programs. By using a series of code transformations, our method allows an arbitrary program to be optimized with respect to a defined subset of its variables, whilst marginalizing out the rest. To carry out the required optimization, we introduce a new GPbased BO package that exploits the availability of the target source code to provide a number of novel features, such as automatic domain scaling and constraint satisfaction.
The concepts we introduce lead directly to a number of extensions of interest, including but not restricted to smart initialization of inference algorithms, adaptive proposals, and nested optimization. Further work might consider maximum marginal likelihood estimation and risk minimization. Though only requiring minor algorithmic changes, these cases require distinct theoretical considerations.
Appendix A Program Transformations in Detail
In this section we give a more detailed and language specific description of our program transformations, code for which can be found at http://www.github.com/probprog/bopp.
a.1 Anglican
Anglican is a probabilistic programming language integrated into Clojure (a dialect of Lisp) and inherits most of the corresponding syntax. Anglican extends Clojure with the special forms sample and observe (Tolpin et al., 2015). Each random draw in an Anglican program corresponds to a sample call, which can be thought of as a term in the prior. Each observe statement applies weighting to a program trace and thus constitutes a term in the likelihood. Compilation of an Anglican program, performed by the macro query, corresponds to transforming the code into a variant of continuationpassing style (CPS) code, which results in a function that can be executed using a particular inference algorithm.
Anglican program code is represented by a nested list of expressions, symbols, nonliterals for contructing data structures (e.g. [...]
for vectors), and command dependent literals (e.g.
[...] as a second argument of a let statement which is used for binding pairs). In order to perform program transformations, we can recursively traverse this nested list which can be thought of as an abstract syntax tree of the program.Our program transformations also make use of the Anglican forms store and retrieve. These allow storing any variable in the probabilistic program’s execution trace in a state which is passed around during execution and from which we can retrieve these stored values. The core use for this is to allow the outer query to return variables which are only locally scoped.
To allow for the early termination that will be introduced in Section A.5, it was necessary to add a mechanism for nonlocal returns to Anglican. Clojure supports nonlocal returns only through Java exception handling, via the keywords try throw, catch and finally. Unfortunately, these are not currently supported by Anglican and their behaviour is far from ideal for our purposes. In particular, for programs containing nested try statements, throwing to a particular try in the stack, as opposed to the most recently invoked, is cumbersome and error prone.
We have instead, therefore, added to Anglican a nonlocal return mechanism based on the Common Lisp control form catch/throw. This uses a catch tag to link each throw to a particular catch. For example
is equivalent to (max a 0). More precisely, throw has syntax (throw tag value) and will cause the catch block with the corresponding tag to exit, returning value. If a throw goes uncaught, i.e. it is not contained within a catch block with a matching tag, a custom Clojure exception is thrown.
a.2 Representations in the Main Paper
In the main paper we presented the code transformations as static transformations as shown in Figure 3. Although for simple programs, such as the given example, these transformations can be easily expressed as static transformations, for more complicated programs it would be difficult to actually implement these as purely static generic transformations in a higherorder language. Therefore, even though all the transformations dynamically execute as shown at runtime, in truth, the generated source code for the prior and acquisition transformations varies from what is shown and has been presented this way in the interest of exposition. Our true transformations exploit store, retrieve, catch and throw to generate programs that dynamically execute in the same way at run time as the static examples shown, but whose actual source code varies significantly.
a.3 Prior Transformation
The prior transformation recursively traverses the program tree and applies two local transformations. Firstly it replaces all observe statements by nil. As observe statements return nil, this trivially preserves the generative model of the program, but the probability of the execution changes. Secondly, it inspects the binding variables of let forms in order to modify the binding expressions for the optimization variables, as specified by the second input of defopt, asserting that these are directly bound to a sample statement of the form (sample dist). The transformation then replaces this expression by one that stores the result of this sample in Anglican’s store before returning it. Specifically, if the binding variable in question is phik, then the original binding expression (sample dist) is transformed into
After all these local transformation have been made, we wrap the resulting query block in a do form and append an expression extracting the optimization variables using Anglican’s retrieve. This makes the optimization variables the output of the query. Denoting the list of optimization variable symbols from defopt as optimargs and the query body after applying all the above location transformations as …, the prior query becomes
Note that the difference in syntax from Figure 3 is because defquery is in truth a syntactic sugar allowing users to bind query to a variable. As previously stated, query is macro that compiles an Anglican program to its CPS transformation. An important subtlety here is that the order of the returned samples is dictated by optimargs and is thus independent of the order in which the variables were actually sampled, ensuring consistent inputs for the BO package.
We additionally add a check (not shown) to ensure that all the optimization variables have been added to the store, and thus sampled during the execution, before returning. This ensures that our assumption that each optimization variable is assigned for each execution trace is satisfied.
a.4 Acquisition Transformation
The acquisition transformation is the same as the prior transformation except we append the acquisition function, ACQF, to the inputs and then observe its application to the optimization variables before returning. The acquisition query is thus
a.5 Early Termination
To ensure that qprior and qacq are cheap to evaluate and that the latter does not include unnecessary terms which complicate the optimization, we wish to avoid executing code that is not required for generating the optimization variables. Ideally we would like to directly remove all such redundant code during the transformations. However, doing so in a generic way applicable to all possible programs in a higher order language represents a significant challenge. Therefore, we instead transform to programs with additional early termination statements, triggered when all the optimization variables have been sampled. Provided one is careful to define the optimization variables as early as possible in the program (in most applications, e.g. hyperparameter optimization, they naturally occur at the start of the program), this is typically sufficient to ensure that the minimum possible code is run in practise.
To carry out this early termination, we first wrap the query in a catch block with a uniquely generated tag. We then augment the transformation of an optimization variable’s binding described in Section A.3 to check if all optimization variables are already stored, and invoke a throw statement with the corresponding tag if so. Specifically we replace relevant binding expressions (sample dist) with
where prologuecode refers to one of the following expressions depending on whether it is used for a prior or an acquisition transformation
We note that valid programs for both qprior and qacq should always terminate via one of these early stopping criteria and therefore never actually reach the appending statements in the query blocks shown in Sections A.3 and A.4. As such, these are, in practise, only for exposition and error catching.
a.6 Marginal/MMAP Transformation
The marginal transformation inspects all let binding pairs and if a binding variable phik is one of the optimization variables, the binding expression (sample dist) is transformed to the following
corresponding to the observe< form used in the main paper.
a.7 Error Handling
During program transformation stage, we provide three errorhandling mechanisms to enforce the restrictions on the probabilistic programs described in Section 3.

We inspect let binding pairs and throw an error if an optimization variable is bound to anything other than a sample statement.

We add code that throws a runtime error if any optimization variable is assigned more than once or not at all.

We recursively traverse the code and throw a compilation error if sample statements of different base measures are assigned to any optimization variable. At present, we also throw an error if the base measure assigned to an optimization variable is unknown, e.g. because the distribution object is from a user defined defdist where the user does not provide the required measure type metainformation.
Appendix B Problem Independent Gaussian Process Hyperprior
Remembering that the domain scaling introduced in Section 4.3 means that both the input and outputs of the GP are taken to vary between , we define the problem independent GP hyperprior as where
(34a)  
(34b)  
(34c)  
(34d)  
(34e) 
The rationale of this hyperprior is that the smoother Matérn 5/2 kernel should be the dominant effect and model the higher length scale variations. The Matérn 3/2 kernel is included in case the evidence suggests that the target is less smooth than can be modelled with the Matérn 5/2 kernel and to provide modelling of smaller scale variations around the optimum.
Appendix C Full Details for House Heating Experiment
In this case study, illustrated in Figure 1, we optimize the parameters of a stochastic engineering simulation. We use the Energy2D system from Xie (2012) to perform finitedifference numerical simulation of the heat equation and NavierStokes equations in a userdefined geometry.
In our setup, we designed a 2dimensional representation of a house with 4 interconnected rooms using the GUI provided by Energy2D. The left side of the house receives morning sun, modelled at a constant incident angle of . We assume a randomly distributed solar intensity and simulate the heating of a cold house in the morning by 4 radiators, one in each of the rooms. The radiators are given a fixed budget of total power density
. The optimization problem is to distribute this power budget across radiators in a manner that minimizes the variance in temperatures across 8 locations in the house.
Energy2D is written in Java, which allows the simulation to be integrated directly into an Anglican program that defines a prior on model parameters and an ABC likelihood for evaluating the utility of the simulation outputs. Figure 2 shows the corresponding program query. In this, we define a Clojure function simulate that accepts a solar power intensity and power densities for the radiators , returning the thermometer temperature readings . We place a symmetric Dirichlet prior on and a gamma prior on , where and are constants. This gives the generative model:
(35)  
(36)  
(37)  
(38) 
After using these to call simulate, the standard deviations of the returned temperatures is calculated for each time point,
(39) 
and used in the ABC likelihood abclikelihood to weight the execution trace using a multivariate Gaussian:
where is the identity matrix and is the observation standard deviation.
Figure 1 demonstrates the improvement in homogeneity of temperatures as a function of total number of simulation evaluations. Visual inspection of the heat distributions also shown in Figure 1 confirms this result, which serves as an exemplar of how BOPP can be used to estimate marginally optimal simulation parameters.
Acknowledgements
Tom Rainforth is supported by a BP industrial grant. Tuan Anh Le is supported by a Google studentship, project code DF6700. Frank Wood is supported under DARPA PPAML through the U.S. AFRL under Cooperative Agreement FA87501420006, Sub Award number 61160290111668.
References

Andrieu et al. (2010)
Christophe Andrieu, Arnaud Doucet, and Roman Holenstein.
Particle Markov chain Monte Carlo methods.
J Royal Stat. Soc.: Series B (Stat. Methodol.), 72(3):269–342, 2010. 
Bérard et al. (2014)
Jean Bérard, Pierre Del Moral, Arnaud Doucet, et al.
A lognormal central limit theorem for particle approximations of normalizing constants.
Electronic Journal of Probability, 19(94):1–28, 2014.  Bergstra et al. (2011) James S Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl. Algorithms for hyperparameter optimization. In NIPS, pages 2546–2554, 2011.
 Broyden (1970) Charles George Broyden. The convergence of a class of doublerank minimization algorithms 1. general considerations. IMA Journal of Applied Mathematics, 6(1):76–90, 1970.
 Carpenter et al. (2015) B Carpenter, A Gelman, M Hoffman, D Lee, B Goodrich, M Betancourt, M A Brubaker, J Guo, P Li, and A Riddell. Stan: a probabilistic programming language. Journal of Statistical Software, 2015.
 Csilléry et al. (2010) Katalin Csilléry, Michael GB Blum, Oscar E Gaggiotti, and Olivier François. Approximate Bayesian Computation (ABC) in practice. Trends in Ecology & Evolution, 25(7):410–418, 2010.
 Devaney et al. (1989) Robert L Devaney, Luke Devaney, and Luke Devaney. An introduction to chaotic dynamical systems, volume 13046. AddisonWesley Reading, 1989.
 Duane et al. (1987) Simon Duane, Anthony D Kennedy, Brian J Pendleton, and Duncan Roweth. Hybrid Monte Carlo. Physics letters B, 1987.
 Eggensperger et al. (2013) Katharina Eggensperger, Matthias Feurer, Frank Hutter, James Bergstra, Jasper Snoek, Holger Hoos, and Kevin LeytonBrown. Towards an empirical foundation for assessing Bayesian optimization of hyperparameters. In NIPS workshop on Bayesian Optimization in Theory and Practice, pages 1–5, 2013.
 Fujii (2013) Keisuke Fujii. Extended Kalman filter. Refernce Manual, 2013.
 Gardner et al. (2014) Jacob R Gardner, Matt J Kusner, Zhixiang Eddie Xu, Kilian Q Weinberger, and John Cunningham. Bayesian optimization with inequality constraints. In ICML, pages 937–945, 2014.
 Goodman et al. (2008) N Goodman, V Mansinghka, D M Roy, K Bonawitz, and J B Tenenbaum. Church: a language for generative models. In UAI, pages 220–229, 2008.
 Goodman and Stuhlmüller (2014) Noah D Goodman and Andreas Stuhlmüller. The Design and Implementation of Probabilistic Programming Languages. 2014.
 Gutmann and Corander (2016) Michael U Gutmann and Jukka Corander. Bayesian optimization for likelihoodfree inference of simulatorbased statistical models. JMLR, 17:1–47, 2016.
 HernándezLobato et al. (2014) José Miguel HernándezLobato, Matthew W Hoffman, and Zoubin Ghahramani. Predictive entropy search for efficient global optimization of blackbox functions. In NIPS, pages 918–926, 2014.
 HernándezLobato et al. (2016) José Miguel HernándezLobato, Michael A. Gelbart, Ryan P. Adams, Matthew W. Hoffman, and Zoubin Ghahramani. A general framework for constrained Bayesian optimization using informationbased search. JMLR, 17:1–53, 2016. URL http://jmlr.org/papers/v17/15616.html.
 Hutter et al. (2011) Frank Hutter, Holger H Hoos, and Kevin LeytonBrown. Sequential modelbased optimization for general algorithm configuration. In Learn. Intell. Optim., pages 507–523. Springer, 2011.
 Jones et al. (1998) Donald R Jones, Matthias Schonlau, and William J Welch. Efficient global optimization of expensive blackbox functions. J Global Optim, 13(4):455–492, 1998.
 Mansinghka et al. (2014) Vikash Mansinghka, Daniel Selsam, and Yura Perov. Venture: a higherorder probabilistic programming platform with programmable inference. arXiv preprint arXiv:1404.0099, 2014.
 Minka et al. (2010) T Minka, J Winn, J Guiver, and D Knowles. Infer .NET 2.4, Microsoft Research Cambridge, 2010.
 Murphy (2012) Kevin P Murphy. Machine learning: a probabilistic perspective. MIT press, 2012.
 Neal (2001) Radford M Neal. Annealed importance sampling. Statistics and Computing, 11(2):125–139, 2001.
 Osborne et al. (2009) Michael A Osborne, Roman Garnett, and Stephen J Roberts. Gaussian processes for global optimization. In 3rd international conference on learning and intelligent optimization (LION3), pages 1–15, 2009.
 Paige et al. (2014) Brooks Paige, Frank Wood, Arnaud Doucet, and Yee Whye Teh. Asynchronous anytime sequential monte carlo. In NIPS, pages 3410–3418, 2014.
 Pickover (1995) Clifford A Pickover. The pattern book: Fractals, art, and nature. World Scientific, 1995.
 Rainforth et al. (2016) Tom Rainforth, TuanAnh Le, JanWillem van de Meent, Michael A Osborne, and Frank Wood. Bayesian optimization for probabilistic programs. In Advances in Neural Information Processing Systems, pages 280–288, 2016.
 Rasmussen and Williams (2006) Carl Rasmussen and Chris Williams. Gaussian Processes for Machine Learning. MIT Press, 2006.
 Ruan et al. (2003) Huawei Ruan, Tongyan Zhai, and Edwin Engin Yaz. A chaotic secure communication scheme with extended Kalman filter based parameter estimation. In Control Applications, 2003. CCA 2003. Proceedings of 2003 IEEE Conference on, volume 1, pages 404–408. IEEE, 2003.
 Shahriari et al. (2016a) Bobak Shahriari, Alexandre BouchardCôté, and Nando de Freitas. Unbounded Bayesian optimization via regularization. AISTATS, 2016a.
 Shahriari et al. (2016b) Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P Adams, and Nando de Freitas. Taking the human out of the loop: A review of Bayesian optimization. Proceedings of the IEEE, 104(1):148–175, 2016b.
 Snoek et al. (2012) Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical Bayesian optimization of machine learning algorithms. In NIPS, pages 2951–2959, 2012.
 Snoek et al. (2015) Jasper Snoek, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Mostofa Patwary, Mostofa Ali, Ryan P Adams, et al. Scalable Bayesian optimization using deep neural networks. In ICML, 2015.
 Tolpin et al. (2015) David Tolpin, JanWillem van de Meent, and Frank Wood. Probabilistic programming in Anglican. Springer International Publishing, 2015. URL http://dx.doi.org/10.1007/9783319234618_36.
 van de Meent et al. (2016) JanWillem van de Meent, Brooks Paige, David Tolpin, and Frank Wood. Blackbox policy search with probabilistic programs. In AISTATS, pages 1195–1204, 2016.
 Wingate et al. (2011) David Wingate, Andreas Stuhlmueller, and Noah D Goodman. Lightweight implementations of probabilistic programming languages via transformational compilation. In AISTATS, pages 770–778, 2011.
 Wood et al. (2014) Frank Wood, Jan Willem van de Meent, and Vikash Mansinghka. A new approach to probabilistic programming inference. In AISTATS, pages 2–46, 2014.
 Xie (2012) Charles Xie. Interactive heat transfer simulations for everyone. The Physics Teacher, 50(4):237–240, 2012.
 Zinkov and Shan (2016) Robert Zinkov and ChungChieh Shan. Composing inference algorithms as program transformations. arXiv preprint arXiv:1603.01882, 2016.
Comments
There are no comments yet.