A* Sampling

Introduction

Drawing samples from arbitrary probability distributions is a core problem in statistics and machine learning. Sampling methods are used widely when training, evaluating, and predicting with probabilistic models. In this work, we introduce a generic sampling algorithm that returns exact independent samples from a distribution of interest. This line of work is important as we seek to include probabilistic models as subcomponents in larger systems, and as we seek to build probabilistic modelling tools that are usable by non-experts; in these cases, guaranteeing the quality of inference is highly desirable. There are a range of existing approaches for exact sampling. Some are specialized to specific distributions , but exact generic methods are based either on (adaptive) rejection sampling or Markov Chain Monte Carlo (MCMC) methods where convergence to the stationary distribution can be guaranteed .

This work approaches the problem from a different perspective. Specifically, it is inspired by an algorithm for sampling from a discrete distribution that is known as the Gumbel-Max trick. The algorithm works by adding independent Gumbel perturbations to each configuration of a discrete negative energy function and returning the argmax configuration of the perturbed negative energy function. The result is an exact sample from the corresponding Gibbs distribution. Previous work has used this property to motivate samplers based on optimizing random energy functions but has been forced to resort to approximate sampling due to the fact that in structured output spaces, exact sampling appears to require instantiating exponentially many Gumbel perturbations.

Our first key observation is that we can apply the Gumbel-Max trick without instantiating all of the (possibly exponentially many) Gumbel perturbations. The same basic idea then allows us to extend the Gumbel-Max trick to continuous spaces where there will be infinitely many independent perturbations. Intuitively, for any given random energy function, there are many perturbation values that are irrelevant to determining the argmax so long as we have an upper bound on their values. We will show how to instantiate the relevant ones and bound the irrelevant ones, allowing us to find the argmax — and thus an exact sample.

There are a number of challenges that must be overcome along the way, which are addressed in this work. First, what does it mean to independently perturb space in a way analogous to perturbations in the Gumbel-Max trick? We introduce the Gumbel process, a special case of a stochastic process recently defined in mathematical statistics , which generalizes the notion of perturbation over space. Second, we need a method for working with a Gumbel process that does not require instantiating infinitely many random variables. This leads to our novel construction of the Gumbel process, which draws perturbations according to a top-down ordering of their values. Just as the stick breaking construction of the Dirichlet process gives insight into algorithms for the Dirichlet process, our construction gives insight into algorithms for the Gumbel process. We demonstrate this by developing A$\!{}^{*}$ sampling, which leverages the construction to draw samples from arbitrary continuous distributions. We study the relationship between A$\!{}^{*}$ sampling and adaptive rejection sampling-based methods and identify a key difference that leads to more efficient use of bound and likelihood computations. We investigate the behaviour of A$\!{}^{*}$ sampling on a variety of illustrative and challenging problems.

The Gumbel Process

The Gumbel-Max trick is an algorithm for sampling from a categorical distribution over classes $i\in\{1,\;\ldots,\;n\}$ with probability proportional to $\exp(\phi(i))$. The algorithm proceeds by adding independent Gumbel-distributed noise to the log-unnormalized mass $\phi(i)$ and returns the optimal class of the perturbed distribution. In more detail, $G\sim\operatorname*{Gumbel}(m)$ is a Gumbel with location $m$ if $\operatorname*{P}(G\leq g)=\exp(-\exp(-g+m))$. The Gumbel-Max trick follows from the structure of Gumbel distributions and basic properties of order statistics; if $G(i)$ are i.i.d. $\operatorname*{Gumbel}(0)$, then $\operatorname*{argmax}_{i}\left\{G(i)+\phi(i)\right\}\sim\exp(\phi(i))/\sum_{i}\exp(\phi(i))$. Further, for any $B\subseteq\{1,\;\ldots,\;n\}$

Eq. 1 is known as max-stability—the highest order statistic of a sample of independent Gumbels also has a Gumbel distribution with a location that is the log partition function . Eq. 2 is a consequence of the fact that Gumbels satisfy Luce’s choice axiom . Moreover, the max and argmax are independent random variables, see Appendix for proofs.

Adapted from . Let $\mu(B)$ be a sigma-finite measure on sample space $\Omega$, $B\subseteq\Omega$ measurable, and $G_{\mu}(B)$ a random variable. $\mathcal{G}_{\mu}=\{G_{\mu}(B)\,|\,B\subseteq\Omega\}$ is a Gumbel process, if

(marginal distributions) $G_{\mu}(B)\sim\operatorname*{Gumbel}\left(\log\mu(B)\right).$

(independence of disjoint sets) $G_{\mu}(B)\;\perp\;G_{\mu}(B^{c}).$

(consistency constraints) for measurable $A,B\subseteq\Omega$, then

Top-Down Construction for the Gumbel Process

While defines and constructs a general class of stochastic processes that include the Gumbel process, the construction that proves their existence gives little insight into how to execute a continuous version of the Gumbel-Max trick. Here we give an alternative algorithmic construction that will form the foundation of our practical sampling algorithm. In this section we assume $\log\mu(\Omega)$ can be computed tractably; this assumption will be lifted in Section 4. To explain the construction, we consider the discrete case as an introductory example.

Suppose $G_{\mu}(i)\sim\operatorname*{Gumbel}(\phi(i))$ is a set of independent Gumbel random variables for $i\in\{1,\;\ldots,\;n\}$. It would be straightforward to sample the variables then build a heap of the $G_{\mu}(i)$ values and also have heap nodes store the index $i$ associated with their value. Let $B_{i}$ be the set of indices that appear in the subtree rooted at the node with index $i$. A property of the heap is that the root $(G_{\mu}(i),i)$ pair is the max and argmax of the set of Gumbels with index in $B_{i}$. The key idea of our construction is to sample the independent set of random variables by instantiating this heap from root to leaves. That is, we will first sample the root node, which is the global max and argmax, then we will recurse, sampling the root’s two children conditional upon the root. At the end, we will have sampled a heap full of values and indices; reading off the value associated with each index will yield a draw of independent Gumbels from the target distribution.

We sketch an inductive argument. For the base case, sample the max and its index $i^{*}$ using their distributions that we know from Eq. 1 and Eq. 2. Note the max and argmax are independent. Also let $B_{i^{*}}=\{0,\ldots,n-1\}$ be the set of all indices. Now, inductively, suppose have sampled a partial heap and would like to recurse downward starting at $(G_{\mu}(p),p)$. Partition the remaining indices to be sampled $B_{p}-\{p\}$ into two subsets $L$ and $R$ and let $l\in L$ be the left argmax and $r\in R$ be the right argmax. Let $[\geq\!\!p]$ be the indices that have been sampled already. Then

where $\mathcal{L}(k)$ and $\mathcal{R}(k)$ denote the left and right children of $k$ and the constraints should only be applied amongst nodes $[\geq\!\!p]\cup\{l,r\}$. This implies

Eq. 4 is the joint density of two independent Gumbels truncated at $G_{\mu}(p)$. We could sample the children maxes and argmaxes by sampling the independent Gumbels in $L$ and $R$ respectively and computing their maxes, rejecting those that exceed the known value of $G_{\mu}(p)$. Better, the truncated Gumbel distributions can be sampled efficiently via CDF inversion$G\sim\operatorname*{TruncGumbel}(\phi,b)$ if $G$ has CDF $\exp(-\exp(-\min(g,b)+\phi))/\exp(-\exp(-b+\phi))$. To sample efficiently, return $G=-\log(\exp(-b-\gamma+\phi)-\log(U))-\gamma+\phi$ where $U\sim\operatorname*{uniform}$., and the independent argmaxes within $L$ and $R$ can be sampled using Eq. 2. Note that any choice of partitioning strategy for $L$ and $R$ leads to the same distribution over the set of Gumbel values.

The basic structure of this top-down sampling procedure allows us to deal with infinite spaces; we can still generate an infinite descending heap of Gumbels and locations as if we had made a heap from an infinite list. The algorithm (which appears as Algorithm 1) begins by sampling the optimal value $G_{1}\sim\operatorname*{Gumbel}(\log\mu(\Omega))$ over sample space $\Omega$ and its location $X_{1}\sim\bar{\mu}(\cdot\,|\,\Omega)$. $X_{1}$ is removed from the sample space and the remaining sample space is partitioned into $L$ and $R$. The optimal Gumbel values for $L$ and $R$ are sampled from a Gumbel with location log measure of their respective sets, but truncated at $G_{1}$. The locations are sampled independently from their sets, and the procedure recurses. As in the discrete case, this yields a stream of $(G_{k},X_{k})$ pairs, which we can think of as being nodes in a heap of the $G_{k}$’s.

If $G_{\mu}(x)$ is the value of the perturbed negative energy at $x$, then Algorithm 1 instantiates this function at countably many points by setting $G_{\mu}(X_{k})=G_{k}$. In the discrete case we eventually sample the complete perturbed density, but in the continuous case we simply generate an infinite stream of locations and values. The sense in which Algorithm 1 constructs a Gumbel process is that the collection $\{\max\{G_{k}\,|\,X_{k}\in B\}\,|\,B\subseteq\Omega\}$ satisfies Definition 1. The intuition should be provided by the introductory argument; a full proof appears in the Appendix. An important note is that because $G_{k}$’s are sampled in descending order along a path in the tree, when the first $X_{k}$ lands in set $B$, the value of $\max\{G_{k}\,|\,X_{k}\in B\}$ will not change as the algorithm continues.

A∗\!{}^{*} Sampling

The Top-Down construction is not executable in general, because it assumes $\log\mu(\Omega)$ can be computed efficiently. A$\!{}^{*}$ sampling is an algorithm that executes the Gumbel-Max trick without this assumption by exploiting properties of the Gumbel process. Henceforth A$\!{}^{*}$ sampling refers exclusively to the continuous version.

This “linearity” allows us to decompose a target log density function into a tractable $i(x)$ and boundable $o(x)$. The tractable component is analogous to the proposal distribution in a rejection sampler. A$\!{}^{*}$ sampling searches for $\operatorname*{argmax}\{G_{k}+o(X_{k})\}$ within the heap of $(G_{k},X_{k})$ pairs from the Top-Down construction of $\mathcal{G}_{\nu}$. The search is an A$\!{}^{*}$ procedure: nodes in the search tree correspond to increasingly refined regions in space, and the search is guided by upper and lower bounds that are computed for each region. Lower bounds for region $B$ come from drawing the max $G_{k}$ and argmax $X_{k}$ of $\mathcal{G}_{\nu}$ within $B$ and evaluating $G_{k}+o(X_{k})$. Upper bounds come from the fact that

where $M(B)$ is a bounding function for a region, $M(B)\geq o(x)$ for all $x\in B$. $M(B)$ is not random and can be implemented using methods from e.g., convex duality or interval analysis. The first term on the RHS is the $G_{k}$ value used in the lower bound.

Comparison to Rejection Samplers

When bounds are refined, A$\!{}^{*}$ sampling bears similarity to adaptive rejection sampling-based algorithms. In particular, while it appears only to have been applied in discrete domains, OS∗ is a general class of adaptive rejection sampling methods that maintain piecewise bounds on the target distribution. If piecewise constant bounds are used (henceforth we assume OS∗ uses only constant bounds) the procedure can be described as follows: at each step, (1) a region $B$ with bound $M(B)$ is sampled with probability proportional to $\nu(B)\exp(M(B))$, (2) a point is drawn from the proposal distribution restricted to the chosen region; (3) standard accept/rejection computations are performed using the regional bound, and (4) if the point is rejected, a region is chosen to be split into two, and new bounds are computed for the two regions that were created by the split. This process repeats until a point is accepted.

Steps (2) and (4) are performed identically in A$\!{}^{*}$ when sampling argmax Gumbel locations and when splitting a parent node. A key difference is how regions are chosen in step (1). In OS∗, a region is drawn according to volume of the region under the proposal. Note that piece selection could be implemented using the Gumbel-Max trick, in which case we would choose the piece with maximum $G_{B}+M(B)$ where $G_{B}\sim\operatorname*{Gumbel}(\log\nu(B))$. In A$\!{}^{*}$ sampling the region with highest upper bound is chosen, where the upper bound is $G_{B}+M(B)$. The difference is that $G_{B}$ values are reset after each rejection in OS∗, while they persist in A$\!{}^{*}$ sampling until a sample is returned.

The effect of the difference is that A$\!{}^{*}$ sampling more tightly couples together where the accepted sample will be and which regions are refined. Unlike OS∗, it can go so far as to prune a region from the search, meaning there is zero probability that the returned sample will be from that region, and that region will never be refined further. OS∗, on the other hand, is blind towards where the sample that will eventually be accepted comes from and will on average waste more computation refining regions that ultimately are not useful in drawing the sample. In experiments, we will see that A$\!{}^{*}$ consistently dominates OS∗, refining the function less while also using fewer likelihood evaluations. This is possible because the persistence inside A$\!{}^{*}$ sampling focuses the refinement on the regions that are important for accepting the current sample.

Experiments

There are three main aims in this section. First, understand the empirical behavior of A$\!{}^{*}$ sampling as parameters of the inference problem and $o(x)$ bounds vary. Second, demonstrate generality by showing that A$\!{}^{*}$ sampling algorithms can be instantiated in just a few lines of model-specific code by expressing $o(x)$ symbolically, and then using a branch and bound library to automatically compute bounds. Finally, compare to OS∗ and an MCMC method (slice sampling). In all experiments, regions in the search trees are hyper rectangles (possibly with infinite extent); to split a region $A$, choose the dimension with the largest side length and split the dimension at the sampled $X_{k}$ point.

In the first experiment, we sample from $p(x)=\exp(-x)/(1+x)^{a}$ for $x>0,a>0$ using $\exp(-x)$ as the proposal distribution. In this case, $o(x)=-a\log(1+x)$ which is unimodal, so the drill down variant of A$\!{}^{*}$ sampling can be used. As $a$ grows, the function becomes peakier; while this presents significant difficulty for vanilla rejection sampling, the cost to A$\!{}^{*}$ is just the cost of locating the peak, which is essentially binary search. Results averaged over 1000 runs appear in Fig. 2 (a).

In the second experiment, we run A$\!{}^{*}$ sampling on the clutter problem , which estimates the mean of a fixed covariance isotropic Gaussian under the assumption that some points are outliers. We put a Gaussian prior on the inlier mean and set $i(x)$ to be equal to the prior, so $o(x)$ contains just the likelihood terms. To compute bounds on the total log likelihood, we compute upper bounds on the log likelihood of each point independently then sum up these bounds. We will refer to these as “constant” bounds. In $D$ dimensions, we generated 20 data points with half within $^{D}$ and half within $^{D}$, which ensures that the posterior is sharply bimodal, making vanilla MCMC quickly inappropriate as $D$ grows. The cost of drawing an exact sample as a function of $D$ (averaged over 100 runs) grows exponentially in $D$, but the problem remains reasonably tractable as $D$ grows ($D=3$ requires 900 likelihood evaluations, $D=4$ requires 4000). The analogous OS∗ algorithm run on the same set of problems requires $16\%$ to $40\%$ more computation on average over the runs.

2 Bounding Strategies

Here we investigate alternative strategies for bounding $o(x)$ in the case where $o(x)$ is a sum of per-instance log likelihoods. To allow easy implementation of a variety of bounding strategies, we choose the simple problem of estimating the mean of a 1D Gaussian given $N$ observations. We use three types of bounds: constant bounds as in the clutter problem; linear bounds, where we compute linear upper bounds on each term of the sum, then sum the linear functions and take the max over the region; and quadratic bounds, which are the same as linear except quadratic bounds are computed on each term. In this problem, quadratic bounds are tight. We evaluate A$\!{}^{*}$ sampling using each of the bounding strategies, varying $N$. See Fig. 2 (b) for results.

For $N=1$, all bound types are equivalent when each expands around the same point. For larger $N$, the looseness of each per-point bound becomes important. The figure shows that, for large $N$, using linear bounds multiplies the number of evaluations by 3, compared to tight bounds. Using constant bounds multiplies the number of evaluations by $O(\sqrt{N})$. The Appendix explains why this happens and shows that this behavior is expected for any estimation problem where the width of the posterior shrinks with $N$.

3 Using Generic Interval Bounds

Here we study the use of bounds that are derived automatically by means of interval methods . This suggests how A$\!{}^{*}$ sampling (or OS∗) could be used within a more general purpose probabilistic programming setting. We chose a number of nonlinear regression models inspired by problems in physics, computational ecology, and biology. For each, we use FuncDesigner to symbolically construct $o(x)$ and automatically compute the bounds needed by the samplers.

Several expressions for $y=f(x)$ appear in the legend of Fig. 2 (c), where letters $a$ through $f$ denote parameters that we wish to sample. The model in all cases is $y_{n}=f(x_{n})+\epsilon_{n}$ where $n$ is the data point index and $\epsilon_{n}$ is Gaussian noise. We set uniform priors from a reasonable range for all parameters (see Appendix) and generated a small (N=3) set of training data from the model so that posteriors are multimodal. The peakiness of the posterior can be controlled by the magnitude of the observation noise; we varied this from large to small to produce problems over a range of difficulties. We use A$\!{}^{*}$ sampling to sample from the posterior five times for each model and noise setting and report the average number of likelihood evaluations needed in Fig. 2 (c) (y-axis). To establish the difficulty of the problems, we estimate the expected number of likelihood evaluations needed by a rejection sampler to accept a sample. The savings over rejection sampling is often exponentially large, but it varies per problem and is not necessarily tied to the dimension. In the example where savings are minimal, there are many symmetries in the model, which leads to uninformative bounds. We also compared to OS∗ on the same class of problems. Here we generated 20 random instances with a fixed intermediate observation noise value for each problem and drew 50 samples, resetting the bounds after each sample. The average cost (heuristically set to # likelihood evaluations plus 2 $\times$ # bound evaluations) of OS∗ for the five models in Fig. 2 (c) respectively was 21%, 30%, 11%, 21%, and 27% greater than for A$\!{}^{*}$.

4 Robust Bayesian Regression

Here our aim is to do Bayesian inference in a robust linear regression model $y_{n}=\boldsymbol{w}^{\textsf{T}}\boldsymbol{x}_{n}+\epsilon_{n}$ where noise $\epsilon_{n}$ is distributed as standard Cauchy and $\boldsymbol{w}$ has an isotropic Gaussian prior. Given a dataset $\mathcal{D}=\{\boldsymbol{x}_{n},y_{n}\}_{n=1}^{N}$ our goal is to draw samples from the posterior $\operatorname*{P}(\boldsymbol{w}\,|\,\mathcal{D})$. This is a challenging problem because the heavy-tailed noise model can lead to multimodality in the posterior over $\boldsymbol{w}$. The log likelihood is $\mathcal{L}(\boldsymbol{w})=\sum_{n}\log(1+(\boldsymbol{w}^{\textsf{T}}\boldsymbol{x}_{n}-y_{n})^{2})$. We generated $N$ data points with input dimension $D$ in such a way that the posterior is bimodal and symmetric by setting $\boldsymbol{w}^{*}=[2,...,2]^{\textsf{T}}$, generating $X^{\prime}\sim\hbox{randn}(N/2,D)$ and $y^{\prime}\sim X^{\prime}\boldsymbol{w}^{*}+.1\times\hbox{randn}(N/2)$, then setting $X=\left[X^{\prime};X^{\prime}\right]$ and $y=\left[y^{\prime};-y^{\prime}\right]$. There are then equally-sized modes near $\boldsymbol{w}^{*}$ and $-\boldsymbol{w}^{*}$. We decompose the posterior into a uniform $i(\cdot)$ within the interval $^{D}$ and put all of the prior and likelihood terms into $o(\cdot)$. Bounds are computed per point; in some regions the per point bounds are linear, and in others they are quadratic. Details appear in the Appendix.

We compare to OS∗, using two refinement strategies that are discussed in . The first is directly analogous to A$\!{}^{*}$ sampling and is the method we have used in the earlier OS∗ comparisons. When a point is rejected, refine the piece that was proposed from at the sampled point, and split the dimension with largest side length. The second method splits the region with largest probability under the proposal.

We ran experiments on several random draws of the data and report performance along the two axes that are the dominant costs: how many bound computations were used, and how many likelihood evaluations were used. To weigh the tradeoff between the two, we did a rough asymptotic calculation of the costs of bounds versus likelihood computations and set the cost of a bound computation to be $D+1$ times the cost of a likelihood computation.

In the first experiment, we ask each algorithm to draw a single exact sample from the posterior. Here, we also report results for the variants of A$\!{}^{*}$ sampling and OS∗ that trade off likelihood computations for bound computations as discussed in Section 4. A representative result appears in Fig. 3 (left). Across operating points, A$\!{}^{*}$ consistently uses fewer bound evaluations and fewer likelihood evaluations than both OS∗ refinement strategies.

In the second experiment, we ask each algorithm to draw 200 samples from the posterior and experiment with the variants that reuse bound information across samples. A representative result appears in Fig. 3 (right). Here we see that the extra refinement done by OS∗ early on allows it to use fewer likelihood evaluations at the expense of more bound computations, but A$\!{}^{*}$ sampling operates at a point that is not achievable by OS∗. For all of these problems, we ran a random direction slice sampler that was given 10 times the computational budget that A$\!{}^{*}$ sampling used to draw 200 samples. The slice sampler had trouble mixing when $D>1$. Across the five runs for $D=2$, the sampler switched modes once, and it did not ever switch modes when $D>2$.

Discussion

This work answers a natural question: is there a Gumbel-Max trick for continuous spaces, and can it be leveraged to develop tractable algorithms for sampling from continuous distributions?

In the discrete case, recent work on “Perturb and MAP” (P&M) methods that draw samples as the argmaxes of random energy functions has shown value in developing approximate, correlated perturbations. It is natural to think about continuous analogs in which exactness is abandoned in favor of more efficient computation. A question is if the approximations can be developed in a principled way, like how showed a particular form of correlated discrete perturbation gives rise to bounds on the log partition function. Can analogous rigorous approximations be established in the continuous case? We hope this work is a starting point for exploring that question.

We do not solve the problem of high dimensions. There are simple examples where bounds become uninformative in high dimensions, such as when sampling a density that is uniform over a hypersphere when using hyperrectangular search regions. In this case, little is gained over vanilla rejection sampling. An open question is if the split between $i(\cdot)$ and $o(\cdot)$ can be adapted to be node-specific during the search. An adaptive rejection sampler would be able to do this, which would allow leveraging parameter-varying bounds in the proposal distributions. This might be an important degree of freedom to exercise, particularly when scaling up to higher dimensions.

There are several possible follow-ons including the discrete version of A$\!{}^{*}$ sampling and evaluating A$\!{}^{*}$ sampling as an estimator of the log partition function. In future work, we would like to explore taking advantage of conditional independence structure to perform more intelligent search, hopefully helping the method scale to larger dimensions. Example starting points might be ideas from AND/OR search or branch and bound algorithms that only branch on a subset of dimensions .

Acknowledgments

This research was supported by NSERC. We thank James Martens and Radford Neal for helpful discussions, Elad Mezuman for help developing early ideas related to this work, and Roger Grosse for suggestions that greatly improved this work.

References

Appendix for “A∗\!{}^{*} Sampling”

In this appendix we prove the main theoretical results of the paper and provide additional experimental details. First, define the following shorthand

Thus $F_{\phi}(g)$ is the CDF and $f_{\phi}(g)$ the PDF of a $\operatorname*{Gumbel}(\phi)$. The following identities are easy to verify and will be reused throughout the appendix.

Joint Distribution of Gumbel Max and Argmax

Suppose $G(i)\sim\operatorname*{TruncGumbel}(\phi(i),b)$ are $n$ independent truncated Gumbels and $Z=\sum_{i=1}^{n}\exp(\phi(i))$, then we are interested in deriving the joint distribution of $i^{*}=\operatorname*{argmax}_{i=1}^{n}G(i)$ and $G(i^{*})=\max_{i=1}^{n}G(i)$.

This is the Gibbs distribution and the density of a $\operatorname*{TruncGumbel}(\log Z,b)$. Thus, for any $B\subseteq\{1,\;\ldots,\;n\}$

These results are well-known. The fact that max Gumbel value has a location that is the log partition function means we can use samples of it as an estimator of log partition functions with known variance $\pi^{2}/6N$ for $N$ samples . Eq. 2 shows that Gumbels satisfy Luce’s choice axiom . In fact, it is a well-known result in random choice theory that the only distribution satisfying Eq. 2 is Gumbel. Notice that the argmax is also independent of the bound, and $b=\infty$ is a valid choice.

Analysis of Top-Down Construction

The goal of this section is to prove that the Top-Down Construction constructs the Gumbel process. In particular, we will argue if we run Algorithm 1 with $\mu$ on $\Omega$, then the collection

is a Gumbel process $\mathcal{G}^{\prime}_{\mu}\overset{d}{=}\mathcal{G}_{\mu}$. In order to do this we consider a special case of Algorithm 1 in which space is not subdivided, $partition(B)=(B,\emptyset)$. In this case the construction takes on a particular simple form, since no queue is needed, see Algorithm 3. We call this special case the In-Order Construction, because is produces the Gumbel values in non-increasing order.

We proceed by arguing that subdividing space has no effect on the distribution of the top $n$ Gumbels. This means that it would be impossible to distinguish a run of Algorithm 3 from a run of Algorithm 1 with the Gumbel values sorted. This allows us to use any choice of $partition$ with a run of Algorithm 1 to analyze the distribution of $\max\{G_{k}\,|\,X_{k}\in B\}$. More precisely

We argue that the top $n$ Gumbels of Algorithm 1 are distributed as in Algorithm 3 regardless of $partition$. That is, if $[i]$ is the index of the $i$th largest Gumbel, then for $1\leq n\leq|\Omega|$

This implies that the distribution over $\{\max\{G_{k}\,|\,X_{k}\in B\}\,|\,B\subseteq\Omega\}$ is invariant under the choice of $partition$ function in Algorithm 1.

We derive the following for a specific choice of $partition$

By the previous result, this is the distribution for any choice of $partition$ (provided it doesn’t produce immeasurable sets) giving us conditions 1. and 2. of Definition 1. Condition 3. is easily satisfied.

This proves the existence of the Gumbel process.

We will proceed to show that the distribution over $\{\max\{G_{k}\,|\,X_{k}\in B\}\,|\,B\subseteq\Omega\}$ is invariant under the choice of $partition$ function. To do so we argue that the top $n$ Gumbels from Algorithm 1 all have the distribution from Algorithm 3. That is, if $[k]$ is the index of the $k$th smallest Gumbel in the tree from Algorithm 1 and $\Omega_{k}=\Omega-\cup_{i=1}^{k-1}\{X_{[k]}\}$. Then for all $n\leq|\Omega|$

Notice that whenever $\mu(\Omega_{k})=\mu(\Omega_{k+1})$ we can omit the removal of $X_{k}$ and still have the same distribution. In the case of continuous $\mu$ we can completely omit all removals and set $\Omega_{k}=\Omega$.

We proceed by induction. For $n=1$, clearly

Now for $1<n\leq|\Omega|$, consider the the top $n$ nodes from a single realization of the process. Let $[<\!\!k]=\{,,\;\ldots,\;[k-1]\}$, the indices of the first $k-1$ Gumbels. By the induction hypothesis we know their distribution and they form a partial tree of the completely realized tree. Our goal is to show that

The boundary of the max partial tree are the nodes $i$ that are on the Queue and have not been expanded. We know that conditioned on $[<\!\!n+1]$ that $G_{[n+1]}=\max_{i\notin[<n+1]}G_{i}$ will come from this boundary, i.e. $G_{[n+1]}=\max_{i\in boundary}G_{i}$. The first step is to realize that the sets $B_{i}$ on the boundary of the max partial tree form a partition of $\Omega_{n+1}$. If $g_{[n+1]}=\max_{i\in boundary}g_{i}$ and $p_{i}$ is the parent of node $i$, then

Because products of indicator functions are like intersections

In other words, the boundary Gumbels are independent and $G_{i}\sim\operatorname*{TruncGumbel}(\log\mu(B_{i}),G_{[n]})$. Notice that the subsets of the boundary form a complete partition of $\Omega_{n+1}$, thus we get

The location $X_{[n+1]}$ has the following distribution:

Again, because the $B_{i}$ is a partition of $\Omega_{n+1}$, this is a mixture distribution in which subsets $B_{[n+1]}$ are sampled with probability $\mu(B_{[n+1]})/\mu(\Omega_{n+1})$ and then $X_{[n+1]}$ is sampled from $\boldsymbol{1}\!\left(x\in B_{[n+1]}\right)\exp(\phi(x))/\mu(B_{[n+1]})$. Thus,

and by the independence of the max and argmax we get that $X_{[n+1]}$ is independent of $G_{[n+1]}$. ∎

Because the joint distribution over the entire collection $\{\max\{G_{k}\,|\,X_{k}\in B\}\,|\,B\subseteq\Omega\}$ is the same regardless of $partition$, this implies that the joint of $\max\{G_{k}\,|\,X_{k}\in B\}$ and $\max\{G_{k}\,|\,X_{k}\in B^{c}\}$ for any specific choice of partition is indeed the joint marginal for any $partition$. In particular we show

Consider the $partition$ that first partitions $\Omega$ into $B$ and $B^{c}$. In this case we consider the distribution over $G_{B}=\max\{G_{k}\,|\,X_{k}\in B\}$ and $G_{B^{c}}=\max\{G_{k}\,|\,X_{k}\in B^{c}\}$ in Algorithm 1. If $X_{1}\in B$, then $G_{B}=G_{1}$ and $G_{B^{c}}=G_{3}$. Otherwise $G_{B^{c}}=G_{1}$ and $G_{B}=G_{2}$. Thus, $G_{B}>G_{B^{c}}$ iff $X_{1}\in B$. Using this knowledge we can split the distribution over $G_{B}$ and $G_{B^{c}}$ into two events.

This is the density of two independent Gumbels with locations $\log\mu(B)$ and $\log\mu(B^{c})$. This proves our result. ∎

Analysis of A∗\!{}^{*} Sampling

This section deals with the correctness and termination of A$\!{}^{*}$ sampling. We exclusively analyze the continuous version of A$\!{}^{*}$ sampling. Recall we have two continuous measures $\mu(B)=\int_{x\in B}\exp(\phi(x))$ and $\nu(B)=\int_{x\in B}\exp(i(x))$ such that we can decompose $\phi(x)$ in a tractable $i(x)$ and intractable but boundable component $o(x)$.

In this section we argue that A$\!{}^{*}$ sampling terminates with probability one by bounding it with the runtime of global-bound A$\!{}^{*}$ sampling. We analyze global-bound A$\!{}^{*}$ more closely.

Consider running A$\!{}^{*}$ with two different sets of bounds on the same realization of the Gumbel process. The returned sample, the final lower bound, and the split chosen for any region will be the same. The only thing that changes is the set of nodes in the tree that are explored. Let $U_{A^{*}}(B)=G_{k}+M(B)$ be the upper bound at node $B$ for A$\!{}^{*}$ and $U(B)=G_{k}+M$ be the upper bound at node $B$ for global-bound A$\!{}^{*}$. Because these algorithms are searching on the same realization we assume that $U(B)\geq U_{A^{*}}(B)$. Let $LB$ be the final lower bound—the optimal node. Because $U(B)\geq U_{A^{*}}(B)$, we know that global-bound A$\!{}^{*}$ visits at least the nodes for which

Finally, A$\!{}^{*}$ never visits nodes for which

So, A$\!{}^{*}$ cannot visit a node that global-bound A$\!{}^{*}$ never visits. Thus, if global-bound A$\!{}^{*}$ terminates with probability one, then so does A$\!{}^{*}$. We now analyze the run time of global-bound A$\!{}^{*}$ more closely and discover a parallel with rejection sampling.

If a global bound $M\geq o(x)$ is reused at every node in A$\!{}^{*}$ sampling, then it takes on a particularly simple form, Algorithm 4; no queue is needed, and it simplifies to a search over the stream of $(G_{k},X_{k})$ values from the In-Order construction (Alg. 3). Global-bound A$\!{}^{*}$ sampling is equivalent to rejection sampling. In particular, both rejection and A$\!{}^{*}$ sampling with constant bounds terminate after $k$ iterations with probability

Rejection terminates with this probability, because the termination condition is independent for each iteration. Thus, the distribution over the number of iterations is a geometric with probability:

That global-bound A$\!{}^{*}$ terminates with this probability is interesting, because the termination condition is not independent from the history of the Gumbel values. Nonetheless, the distribution over iterations is memoryless. Consider the stream of values $(G_{k},X_{k})$,

Global-bound A$\!{}^{*}$ terminates when $\max_{1\leq i\leq k}\{G_{i}+o(X_{i})\}\geq G_{k+1}+M$. In order to show that the distribution of $k$ is geometric we need simply to show that

First we show that the finite differences $D_{k}=G_{k+1}-G_{k}$ are mutual independent and

Inspecting the joint pdf for $G_{i}$ with $1\leq i\leq k+1$

We proceed with an inductive argument. First the base, with $k=1$, so for $d_{1}<0$

Now we do a sequence of tricky substitutions, $r_{i}=d_{i}+\sum_{j=i+1}^{k}d_{j}$ from $i=1$ to $k$, and find that this integral equals

Notice that this is basically an infinite triangle over a function that is symmetric in the $r_{i}$. Since it’s multiplied by $k!$ it is equal to the sum over all the permutations of the $r_{i}$, which ends up giving us an infinite cube:

All that remains is to evaluate $\int_{r=-\infty}^{0}H\left(r\right)\exp\left(r\right)$:

Partial Correctness of A∗\!{}^{*} Sampling

In this section we show that given termination the distribution returned by A$\!{}^{*}$ is correct. This depends only on the linearity result about the Gumbel process. We prove the linearity by using the auxiliary $\mathcal{G}_{\nu}$ to measure the bounded difference $o(x)$.

At termination A$\!{}^{*}$ returns $LB$ and $X^{*}$:

The correctness of the construction for $\mathcal{G}_{\nu}$ implies the consistency and independence requirements. Thus we need only to verify the marginal of $\max\{G_{k}+o(X_{k})\,|\,X_{k}\in B\}$ is $\operatorname*{Gumbel}(\log\mu(B))$. To see why this is the case, consider a partition $p_{1},\;\ldots,\;p_{n}$ of the range of $o(x)$ and let

$\{G_{k}\,|\,X_{k}\in B\cap P_{j}\}\sim\operatorname*{Gumbel}(\log\nu(B\cap P_{j}))$, because $G_{k}$ and $X_{k}$ are samples from the process $G_{i}(B)$. Thus,

as the partition gets finer. Since $\int_{x\in B}\exp(i(x))\exp(o(x))=\mu(B)$ we get that

Thus the distribution of $\max\{G_{k}+o(X_{k})\,|\,X_{k}\in B\}$ must be $\operatorname*{Gumbel}(\log\mu(B))$ and we’re done. ∎

Explanation of Results from Section 6.2

Consider running A$\!{}^{*}$ with two different sets of bounds on the same realization of the Gumbel process. The returned sample, the final lower bound, and the split chosen for any region will be the same. The only thing that changes is the set of regions that are explored. Let $U_{1}(B)$ be an optimal bound and $U_{2}(B)$ a suboptimal bound of region $R$. The question is how many more regions are explored by using $U_{2}$ instead of $U_{1}$. Let $R^{*}$ be the region producing the returned sample. Once this region is explored, the lower bound reaches its final value and A$\!{}^{*}$ will only explore regions with $U_{2}(B)>LB$. These regions will be new, have been unexplored by $U_{1}$, if $U_{1}(B)<LB$. Before $R^{*}$ is explored, A$\!{}^{*}$ will only explore regions with $U_{2}(B)>U_{2}(R^{*})$. Since $U_{2}(R^{*})>LB$, the condition that determines whether a region $R$ is new is $U_{2}(B)>LB>U_{1}(B)$.

Suppose $R$ is explored using $U_{1}$, but its descendants are not. How many of its descendants are explored using $U_{2}$? For the bounds we consider in Section 6.2, the suboptimality of $U_{2}$ is proportional to the region width. Let $B_{d}$ be the deepest descendant explored using $U_{2}$, and let $n$ be the width of $R$ divided by $B_{d}$. By the assumption on the suboptimality of bounds being proportional to region width, we have $(U_{2}(B)-U_{1}(B))/n>U_{2}(B_{d})-U_{1}(B_{d})$, and thus $(U_{2}(B)-U_{1}(B))/n+U_{1}(B_{d})>LB$. This implies that $n$ is bounded by a linear function of the suboptimality $U_{2}(B)-U_{1}(B)$. The balanced nature of the splitting process suggests that $B_{d}$ has depth $\log_{2}(n)$ with high probability, so $n$ also bounds the total number of explored descendants. Therefore the total number of additional regions explored by $U_{2}$ is linear in the total suboptimality of the bounds of regions explored by $U_{1}$.

When the log-likelihood is a sum of $n$ terms and we apply a constant bound to each term, the suboptimality of the total bound grows linearly with $n$. However, under the conditions of the Bernstein-von Mises theorem, the posterior will concentrate around a peak of width $O(n^{-1/2})$. This shrinks the width of the significant regions, reducing the suboptimality of an explored region to $O(\sqrt{n})$. This is the trend we see in the plot. If we apply a linear bound to each term, the suboptimality is $O(({\rm width})^{2})$ per term by the Taylor remainder theorem, and overall $O(n({\rm width})^{2})$ which is constant in $n$ for regions around the peak. Therefore we expect linear termwise bounds to explore a constant multiple of the number of regions explored by optimal bounds.

Bounding the Cauchy Log Likelihood

To perform inference in the Bayesian Robust Regression experiment, we need to upper bound $\max_{\boldsymbol{w}\in[\underline{\boldsymbol{w}},\overline{\boldsymbol{w}}]}\mathcal{L}(\boldsymbol{w})$ for each $[\underline{\boldsymbol{w}},\overline{\boldsymbol{w}}]$ encountered in the search tree. For each region, we compute the bound in two steps. First, for each $n$, compute the minimum and maximum possible values of $d_{n}=\boldsymbol{w}^{\textsf{T}}\boldsymbol{x}_{n}-y_{n}$ using interval arithmetic , yielding $\underline{d_{n}}$ and $\overline{d_{n}}$. For each $n$, we then construct a quadratic bound on the Cauchy likelihood term $C(d)=-\log(1+d^{2})$ that is guaranteed to be an upper bound so long as $d_{n}\in[\underline{d_{n}},\overline{d_{n}}]$. The bound is referred to as $B_{n}(d_{n})$ and takes the form $B_{n}(d_{n})=a_{n}d_{n}^{2}+b_{n}d_{n}+c_{n}$.

The second derivative of the Cauchy likelihood $C^{\prime\prime}(d)=\frac{4d^{2}}{(d^{2}+1)^{2}}-\frac{2}{d+1}$ changes sign only twice, at $-1$ and $1$. Outside of $$,$C^{\prime\prime}(d)$is positive (i.e.,$C$is convex), and inside it is negative (i.e.,$C$is concave). If an interval$[\underline{d_{n}},\overline{d_{n}}]$is fully in a convex region, then we use a simple linear bound of the line that passes through$(\underline{d_{n}},C(\underline{d_{n}}))$and$(\overline{d_{n}},C(\overline{d_{n}}))$. If an interval$[\underline{d_{n}},\overline{d_{n}}]$is fully in a concave region and does not contain , then we use a linear bound that is tangent to$C$at the midpoint of$[\underline{d_{n}},\overline{d_{n}}]$. If the interval contains any of$\{-1,0,1\}$then we use a quadratic bound. If the interval contains one of$-1$or$1$, then we expand the interval to include 0 and then proceed. To compute the bound, we fix the bound function to have$B_{n}(0)=C(0)$and$B^{\prime}_{n}(0)=C^{\prime}(0)$. Since the bound is quadratic, its second derivative$2a_{n}$is constant. We set this constant to be the most negative value that ensures that$B_{n}(d)$is a valid bound over all of$[\underline{d},\overline{d}]$. Concretely, the quadratic bound$B(d)=ad^{2}+bd+c$is constructed as follows. For each endpoint$d_{end}\in\{\underline{d_{n}},\overline{d_{n}}\}$ that is not equal to 0, compute

and choose the largest computed $a$. Finally, solve for $b$ and $c$ to ensure that the derivative and value of the bound match $C$ at 0:

Regression Experiment Priors

All parameters were given uniform priors. The ranges are as follows:

$y=x\frac{\cos(a)(x\sin(a)+\sqrt{x^{2}\sin(a)^{2}+2bc})}{b}$