Scalable Discrete Sampling as a Multi-Armed Bandit Problem

Introduction

Sampling a random variable from a discrete (conditional) distribution is one of the core operations in Monte Carlo methods. It is an ubiquitous and often necessary component for inference algorithms such as Gibbs sampling and particle filtering. Applying discrete sampling for large-scale problems has been a challenging task like other Monte Carlo algorithms due to the high computational burden. Various approaches have been proposed to address different dimensions of “large scales”. For example, distributed algorithms have been used to sample a model with a large number of discrete variables (Newman et al., 2009; Bratières et al., 2010; Wu et al., 2011), smart transition kernels were described for Markov chain Monte Carlo (MCMC) algorithms to sample efficiently a single variable with a large or even infinite state space (Li et al., 2014; Kalli et al., 2011). This paper is focused on another dimension of the “large-scales” where the variable to sample has a large degree of statistical dependency.

Consider a random variable with a finite domain $X\in{\cal X}$ and a distribution in the following form

A common approach to address the big data problem is divide-and-conquer that uses parallel or distributed computing resources to process data in parallel and then synchronize the results periodically or merely once in the end (Scott et al., 2013; Medlar et al., 2013; Xu et al., 2014).

An orthogonal approach has been studied for the Metropolis-Hastings (MH) algorithm in a general state space by running a sampler with subsampled data. This approach can be combined easily with the distributed computing idea for even better scalability (e.g. Ahn et al., 2015).

Maclaurin & Adams (2015) introduced an MH algorithm in an augmented state space that could achieve higher efficiency than the standard MH by processing only a subset of active data every iteration while still preserving the correct stationary distribution. But the introduction of auxiliary variables might also slow down the overall mixing rate in the augmented space.

Approximate MH algorithms have been proposed in the subsampling approach with high scalability. The stochastic gradient Langevin dynamics (SGLD) (Welling & Teh, 2011) and its extensions (Ahn et al., 2012; Chen et al., 2014; Ding et al., 2014) introduced efficient proposal distributions based on subsampled data. Approximate algorithms induce bias in the stationary distribution of the Markov chain. But given a fixed amount of runtime they could reduce the expected error in the Monte Carlo estimate via a proper trade-off between variance and bias by mixing faster w.r.t. the runtime. This is particularly important for large-scale learning problems when the runtime is one of the limiting factors for generalization performance (Bottou & Bousquet, 2008). However, the stochastic gradient MCMC approach usually skips the rejection step in order to obtain sublinear time complexity and the induced bias is very hard to estimate or control.

Another line of research on approximate subsampled MH algorithms does not ignore the rejection step but controls the error with an approximate rejection step based on a subset of data (Korattikara et al., 2014; Bardenet et al., 2014). The bias can thus be better controlled (Mitrophanov, 2005; Pillai & Smith, 2014). That idea has also been extended to slice sampling (DuBois et al., 2014) and Gibbs for binary variables (Korattikara et al., 2014).

In this paper we follow the last line of research and propose a novel approximate sampling algorithm to improve the scalability of sampling discrete distributions. We first reformulate the problem in Eq. 1 as a Multi-Armed Bandit (MAB) problem with a finite reward population via the Gumbel-Max trick (Papandreou & Yuille, 2011), and then propose three algorithms with theoretical guarantees on the approximation error and an upper bound of $N|{\cal X}|$ on the sample size. This is to our knowledge the first attempt to address discrete sampling problem with a large number of dependencies and our work will likely contribute to a more complete library of scalable MCMC algorithms. Moreover, the racing algorithm in Sec. 3.3 provides a unified framework for subsampling-based discrete sampling, MH (Korattikara et al., 2014; Bardenet et al., 2014) and slice sampling (DuBois et al., 2014) algorithms as discussed in Sec. 4. We also show in the experiments that our algorithm can be combined straightforwardly with stochastic gradient MCMC to achieve both high efficiency and controlled bias. Lastly, the proposed algorithms also deserve their own interest for MAB problems under this particular setting.

We first review an alternative way of drawing discrete variables and build a connection with MABs in Sec. 2, then propose three algorithms in Sec. 3. We discuss related work in Sec. 4 and evaluate the proposed algorithms on both synthetic data and real-world problems of Bayesian inference and graphical model inference in Sec. 5. Particularly, we show how our proposed sampler can be combined conveniently as a building component with other subsampling sampler for a hierarchical Bayesian model. Sec. 6 concludes the paper with a discussion.

Approximate Discrete Sampling

2 Approximate Discrete Sampling as a Multi-Armed Bandits Problem

In a Multi-Armed Bandit (MAB) problem, the $i$’th bandit is a slot machine with an arm, which when pulled generates an i.i.d. reward $l_{i}$ from a distribution associated with that arm with an unknown mean $\mu_{i}$. The optimal arm identification problem for MABs (Bechhofer, 1958; Paulson, 1964) in the fixed confidence setting is to find the arm with the highest mean reward with a confidence $1-\delta$ using as few pulls as possible.

Under the assumption of Eq. 1, the solution in Eq. 2 can be expressed as

where $\mathcal{L}_{i}\overset{\textnormal{def}}{=}\{l_{i,1},l_{i,2},\dots,l_{i,N}\}$. After drawing $D$ Gumbel variables $\varepsilon_{i}$, we turn the discrete sampling problem into the optimal arm identification problem in MABs where the reward $l_{i}$ is uniformly sampled from a finite population $\mathcal{L}_{i}$. An approximate algorithm that solves the problem with a fixed confidence may avoid drawing all the rewards from an obviously sub-optimal arm and save computations. We show the induced bias in the sample distribution as follows with the proof in Appx. A.1.

If an algorithm solves (2) exactly with a probability at least $1-\delta$ for any value of $\boldsymbol{\varepsilon}$, the total variation between the sample distribution $\hat{p}$ and the true distribution is bounded by

When applied in the MCMC framework as a transition kernel, we can apply immediately the theories in Mitrophanov (2005); Pillai & Smith (2014) to show that the approximate Markov chain satisfies uniform ergodicity under regular conditions and the analysis of convergence rate are readily available under various assumptions. So the discrete sampling problem of this paper reduces to finding a good MAB algorithm for Eq. 2 in our problem setting.

Algorithms for MABs with a Finite Population and Fixed Confidence

The key difference of our problem from the regular MABs is that our rewards are generated from a finite population while regular MABs assume i.i.d. rewards. Because one can obtain the exact mean by sampling all the $N$ values $l_{i,n}$ for arm $i$ without replacement, a good algorithm should pull no more than $N$ times for each arm regardless of the mean gap between arms. We introduce three algorithms in this section whose sample complexity is upper bounded by $O(ND)$ in the worst case and can be very efficient when the mean gap is large.

2 Adapted lil’UCB

We first study one of the state-of-the-art algorithms for fixed-confidence optimal arm identification problem and adjust it for the finite population setting. The lil’UCB algorithm (Jamieson et al., 2014) maintains an upper confidence bound (UCB) of $\mu_{i}$ that is inspired by the law of the iterated logarithm (LIL) for every arm. At each iteration, it draws a single sample from the arm with the highest bound and updates it. The algorithm terminates when some arm is sampled much more often than all the other arms. We refer readers to Fig. 1 of Jamieson et al. (2014) for details. The time complexity for $t$ iterations is $\mathcal{O}(\log(D)t)$. It was shown in Jamieson et al. (2014) that lil’UCB achieved the optimal sample complexity up to constants.

However, lil’UCB requires i.i.d. rewards for each arm $i$, that is, sampled with replacement from ${\cal L}_{i}$. Therefore, the total number of samples $t$ is unbounded and could be $\gg ND$ when the means are close to each other. We adapt lil’UCB for our problem with the following modifications:

Samples $l_{i,n}$ without replacement for each arm but keep different arms independent.

When $T_{i}^{(t)}=N$ for some arm $i$, the estimate $\hat{\mu}_{i}^{(t)}$ becomes exact. So set its UCB to $\hat{\mu}_{i}^{(t)}$.

The algorithm terminates either with the original stopping criterion or when the arm with the highest upper bound has an exact mean estimate, whichever comes first.

The adapted algorithm satisfies all the theoretical guarantees in Thm. 2 of Jamieson et al. (2014) with additional properties as shown in the following proposition with proof in Appx. A.2.

Theorem 2 of Jamieson et al. (2014) holds for the adapted lil’UCB algorithm. Moreover $T_{i}^{(t)}\leq N,\forall i,t$. Therefore, when the algorithm terminates, $t=\sum_{i\in{\cal X}}T_{i}^{(t)}\leq ND$.

Notice that Thm. 2 of Jamieson et al. (2014) shows that $t$ scales roughly as $O(1/\Delta^{2})$ with $\Delta$ being the mean gap and therefore $t\ll ND$ when the gap is large.

3 Racing Algorithm for a Finite Population

When rewards are sampled without replacement, the negative correlation between rewards would generally improve the convergence of $\hat{\mu}_{i}$. Unfortunately, the bound in lil’UCB ignores the negative correlation when $T_{i}^{(t)}<N$ even with the adaptations. We introduce a new family of racing algorithms (Maron & Moore, 1994) that takes advantage of the finite population setting as shown in Alg. 1. The choice of the uncertainty bound function $G$ differentiates specific algorithms and two examples will be discussed in the following sections.

Alg. 1 maintains a set of candidate set ${\cal D}$ initialized with all arms. At iteration $t$, a shared mini-batch of $m^{(t)}$ indices are drawn w/o replacement for all survived arms in ${\cal D}$. Then the uncertainty bound $G$ is used to eliminate sub-optimal arms with a given confidence. The algorithm stops when only one arm remains. We require for $m^{(t)}$ that the total number of sampled indices $T^{(t^{*})}=\sum_{t=1}^{t^{*}}m^{(t)}$ equals $N$ at the last iteration $t^{*}$. Particularly, we take a doubling schedule $T^{(t)}=2T^{(t-1)}$ (so $t^{*}=\lceil\log_{2}\frac{N}{m^{(1)}}\rceil+1$) and leave $m^{(1)}$ as a free parameter. We also require $G(\cdot,T,\cdot,\cdot)=0$ whenever $T=N$ so that Alg. 1 always stops within $t^{*}$ iterations. The computational complexity for $t$ iterations is $\mathcal{O}(DT^{(t)})$ with the marginal estimate $\hat{\sigma}_{i}$ and $\mathcal{O}(D^{2}T^{(t)})$ with the pairwise estimate $\hat{\sigma}_{i,j}$. The former version is more efficient than the latter when $D$ is large at the price of a looser bound.

for any $\delta\in(0,1)$ with a probability at least $1-\delta$, Alg. 1 returns the optimal arm with at most $ND$ samples.

The proof is provided in Appx. A.3. Unlike adapted lil’UCB, Racing draws a shared set of sample indices among all the arms and could provide a tighter bound with pairwise variance estimates $\hat{\sigma}_{i,j}$ when there is positive correlation, a typical case in Bayesian inference problems.

Serfling (1974) studied the concentration inequalities of sampling without replacement and obtained an improved Hoeffding bound. Bardenet & Maillard (2013) extended the work and provided an empirical Bernstein-Serfling bound that was later used for the subsampling-based MH algorithm (Bardenet et al., 2014): for any $\delta\in(0,1]$ and any $n\leq N$, with probability $1-\delta$, it holds that

where $\kappa=\frac{7}{3}+\frac{3}{\sqrt{2}}$, and \rho_{n}=\left\{\begin{array}[]{ll}1-\pi_{n-1}&\mbox{if }n\leq N/2\\ (1-\pi_{n})(1+\frac{1}{n})&\mbox{if }n>N/2\end{array}\right., with $\pi_{n}\overset{\textnormal{def}}{=}\dfrac{n}{N}$. The extra term $\rho_{n}$ that is missing in regular empirical Bernstein bounds reduces the bound significantly when $n$ is close to $N$. We set $m^{(1)}=2$ in Alg. 1 to provide a valid $\hat{\sigma}^{(t)}$ for any $t$ and set the uncertain bound $G$ with the empirical Bernstein-Serfling (EBS) bounds as

It is trivial to prove that $G_{\text{EBS}}$ satisfies the condition in Eq. 5 using a union bound over $t<t^{*}$.

3.2 Racing with a Normal Assumption for G𝐺G

The concentration bounds often give a conservative strategy as they assume an arbitrary bounded reward distribution. When the number of drawn samples is large, the central limit theorem suggests that $\hat{\mu}^{(t)}$ follows approximately a Gaussian distribution. Korattikara et al. (2014) made such an assumption and obtained a tighter bound. We first provide an immediate corollary of Prop. 2 in Appx. A of Korattikara et al. (2014).

With the normal assumption, we choose the uncertainty bound $G$ in the following form

Let $x^{*}$ be the best arm and $\Delta$ be the minimal normalized gap of means from other arms, defined as $\min_{i\neq x^{*}}\frac{\mu_{x^{*}}-\mu_{i}}{\sigma_{x^{*}}+\sigma_{i}}$ when using marginal variance estimate $\hat{\sigma}_{i}$ and $\min_{i\neq x^{*}}\frac{\mu_{x^{*}}-\mu_{i}}{\sigma_{x^{*},i}}$ when using pairwise variance estimate $\hat{\sigma}_{x,i}$. If Assump. 6 holds, with a probability at least $1-\delta$ Racing-Normal draws no more rewards than

where $\lceil n\rceil_{m}\overset{\textnormal{def}}{=}m2^{\lceil\log_{2}n/m\rceil}\wedge N\geq n,\forall n\leq N$. $D^{\prime}\overset{\textnormal{def}}{=}D$ if using $\hat{\sigma}_{i}$ and is $D-1$ if using $\hat{\sigma}_{x,i}$.

4 Variance Reduction for Random Rewards with Control Variates

The control variate method is mostly useful for Racing-Normal. For algorithms depending on a reward bound $C$ in order to get a tight bound for $l_{i,n}-h_{i,n}$ it requires a more restrictive condition for C as in Bardenet et al. (2015) and we might end up with an even more conservative strategy in general cases.

Related Work

The Gumbel-Max trick has been exploited in Kuzmin & Warmuth (2005); Papandreou & Yuille (2011); Maddison et al. (2014) for different problems. The closest work is Maddison et al. (2014) where this trick is extended to draw continuous random variables with a Gumbel process, reminiscent to adaptive rejection sampling.

Our work is closely related to the optimal arm identification problem for MABs with a fixed confidence. This is, to our knowledge, the first work to consider MABs with a finite population. The proposed algorithms tailored under this setting could be of interest beyond the discrete sampling problem. The normal assumption in Sec. 3.3.2 is similar to UCB-Normal in Auer et al. (2002) but the latter assumes a normal distribution for individual rewards and will perform poorly when it does not hold.

The bounds in Sec. 3.3 are based on subsampling-based MH algorithms in Bardenet et al. (2014); Korattikara et al. (2014). The proposed algorithm extends those ideas from MH to discrete sampling. In fact, let $x$ and $x^{\prime}$ be the current and proposed value in an MH iteration, Racing-EBS and Racing-Normal reduce to the algorithms in Bardenet et al. (2014) and Korattikara et al. (2014) respectively if we set

The variance reduction technique is similar to the proxies in Bardenet et al. (2015), but the control variate here is a function in the data space while the proxy in the latter is a function in the parameter space. We do not assume the posterior distribution is approximate Gaussian and our algorithm works with multi-modal distributions.

It is important not the confuse the focus of our algorithm for the big $N$ problem in Eq. 1 with other algorithms that address sampling for a large state space (big $D$) or similarly a high-dimensional vector of discrete variables (exponentially large $D$). The combination of these two approaches for problems with both big $N$ and big $D$ is possible but beyond the scope of this paper.

Experiments

Since this is the first work to discuss efficient discrete sampling for problem (1), we compare the adapted lil’UCB, Racing-EBS, Racing-Normal with the exact sampler only. We report the result of Racing-Normal in real data experiments only as the speed gains of the other two are marginal.

Surprisingly, Racing-Normal performs robustly regardless of reward distributions with the first mini-batch size $m^{(1)}=50$ while it was shown in Bardenet et al. (2014) that the algorithm with the same normal assumption in Korattikara et al. (2014) failed with LogNormal even when $m^{(1)}=500$. The dramatic improvement in robustness is mainly due to our doubling scheme where central limit theorem applies quickly with $m^{(t)}$ increasing exponentially. We do not claim that the single trick will solve the problem completely because there still exist cases in theory with extremely heavy-tailed reward distributions where our normal assumption does not hold and the algorithm will fail to meet the confidence level. In practice, we do not observe that pathological case in any of the experiments.

2 Bayesian ARCH Model Selection

We evaluate Racing-Normal in a Bayesian model selection problem for the auto-regressive conditional heteroskedasticity (ARCH) models. The discrete sampler is integrated in the Markov chain as a building component to sample the hierarchical model. Specifically, we consider a mixture of ARCHs for the return $r_{t}$ of stock price series with student-t innovations, each component with a different order $q$:

3 Author Coreference

We then study the performance in a large-scale graphical model inference problem. The author coreference problem for a database of scientific paper citations is to cluster the mentions of authors into real persons. Singh et al. (2012) addressed this problem with a conditional random field model with pairwise factors. The joint and conditional distributions are respectively

We run the experiment on the union of an unlabeled DBLP dataset of BibTex entries with about 5M authors and a Rexa corpus of about 11K author mentions with 3160 entries labeled. We monitor the clustering performance on the labeled subset with the $B^{3}$ F-1 score (Bagga & Baldwin, 1998). We use $\delta=0.05$ and the empirical error rate is about $0.046$. The number of candidate values $D$ varies in $2\sim 215$ and $N_{y}$ varies in $1\sim 1829$ upon convergence. Fig. 3(a) shows the F-1 score as a function of the number of factor evaluations with 7 random runs for each algorithm. Sub Gibbs converges about three times faster than exact Gibbs. Fig. 3(b) shows F-1 as a function of iterations that renders almost identical behavior for both algorithms, which suggests negligible bias in Sub Gibbs. The relative number of the evaluated factors of sub to exact Gibbs indicates about a 5-time speed up near convergence. The initial speed up is small because every cluster is initialized with a single mention, i.e. $N_{y}=1$.

Discussion

We consider the discrete sampling problem with a high degree of dependency and proposed three approximate algorithms under the framework of MABs with theoretical guarantees. The Racing algorithm provides a unifying approaches to various subsampling-based Monte Carlo algorithms and also improves the robustness of the original MH algorithm in Korattikara et al. (2014). This is also the first work to discuss MABs under the setting of a finite reward population.

Empirical evaluations show that Racing-Normal achieves a robust and the highest speed-up among all competitors. Whilst adaptive lil’UCB shows inferior empirical performance to Racing-Normal, it has a better sample complexity w.r.t. the number of arms $D$. It will be a future direction to combine the bound of Racing-Normal with other MAB algorithms including lil’UCB for a better scalability in $D$. Another important problem is on how to relax the assumptions for Racing-Normal without sacrificing the performance.

It would also be an interesting direction to extend our work to draw continuous random variables efficiently with the Gumbel process (Maddison et al., 2014). In continuous state space, there are infinitely many “arms” and a naive application of our algorithm will lead to infinitely large error bound. This problem can be alleviated with algorithms for contextual MAB problems.

Acknowledgements

We thank Matt Hoffman for helpful discussions on the connection of our work to the MAB problems. We also thank all the reviewers for their constructive comments. We acknowledge funding from the Alan Turing Institute, Google, Microsoft Research and EPSRC Grant EP/I036575/1.

References

Appendix A Proofs

For a discrete state space, the total variation is equivalent to half of $L_{1}$ distance between two probability vectors. Denote by $\hat{p}(X=i|\boldsymbol{\varepsilon})$ the distribution of the output of the approximate algorithm conditioned on the vector of Gumbel variables $\boldsymbol{\varepsilon}$, and $x(\boldsymbol{\varepsilon})$ the solution of Eq. 2 as a function of $\boldsymbol{\varepsilon}$. According to the premise of Prop. 1, $\hat{p}(X=x(\boldsymbol{\varepsilon})|\boldsymbol{\varepsilon})\geq 1-\delta,\forall\boldsymbol{\varepsilon}$. We can bound the $L_{1}$ error of the conditional probability as

where $\delta_{i,j}$ is the Kronecker delta function. Then we can show

A.2 Sketch of the proof of Prop. 2

As the proof of this proposition is almost identical to the proof of Jamieson et al. (2014), we only outlines the difference due to the adaptation. In the proof of Thm. 2 in Jamieson et al. (2014), the i.i.d. assumption for rewards from each arm was used only in Lemma 3 to provide Chernoff’s bound and Hoeffding’s bound. As noted in Sec. 6 of Hoeffding (1963) those bounds would still hold when rewards are sampled from a finite population without replacement. Therefore, when $T^{(t)}<N$ all the bounds hold for adapted lil’UCB.

When $T_{i}^{(t)}=N$, the second modification sets the upper bound of the mean estimate to $\hat{\mu}^{(t)}$. That is a valid upper bound of $\mu_{i}$, in fact much tighter than the bound in the original algorithm because $\hat{\mu}^{(t)}_{i}=\mu_{i}$ exactly when the entire population is observed.

Therefore, as long as $T_{i}^{(t)}\leq N,\forall i$, Theorem 2 in Jamieson et al. (2014) applies to adapted lil’UCB with modification 1 and 2 only.

With the third modification, $T^{(t)}$ could never be bigger than $N$ at the stopping time, which proves the second part of Prop 2. The proof can then be concluded if we can show modification 3 does not change the output of adapted lil’UCB with the first two modifications only. This is true because if we do not stop when the selected arm $i$ satisfies $T_{i}^{(t)}=N$, we do not need to update the upper bound of $i$ because the estimated mean is already exact. Since no upper bound is changed, the arm $i$ will always be chosen for now on and eventually the original stopping criterion of $T_{i}^{(t)}\geq 1+\lambda\sum_{j\neq i}T_{j}(t)$ is met and the same arm $i$ will be returned. ∎

A.3 Proof of Prop. 3

Denote by $x^{(t)}$ the arm with the highest estimated mean at iteration $t$ and $x^{*}$ the optimal arm with the highest true mean, $\mu_{x^{*}}>\mu_{i},\forall i\neq x^{*}$. If Alg. 1 does not stop in the first $t^{*}-1$ iterations, the estimated means of all the survived arms become exact at the last iteration $t^{*}$, $\hat{\mu}_{i}^{(t^{*})}=\mu_{i}$ because we require $T^{(t^{*})}=N$. Then $x^{(t^{*})}=x^{*}$. As we require $G(\delta,T=N,\hat{\sigma},C)=0,\forall\delta,\hat{\sigma},C$, all the sub-optimal arms will be eliminated by the last iteration and the algorithm always returns the correct best arm. This proves the upper bound of the sample size of $ND$.

Now to prove the confidence level, all we need to show is that with at least a probability $1-\delta$ arm $x^{*}$ survived all the iterations $t<t^{*}$.

Let us first consider the case when Alg. 1 uses the marginal variance estimate $\hat{\sigma}_{i}^{(t)}$. Let the events

Applying condition Eq. 5 and the union bound, we get $P(\cup_{i\in{\cal X}}E_{i})\leq\sum_{i\in{\cal X}}E_{i}=\delta.$ So with a probability at least $1-\delta$, none of those events will happen. In that case for any iteration $t<t^{*}$,

So arm $x^{*}$ won’t be eliminated at iteration $t$.

Similarly, for the case when Alg. 1 uses the pairwise variance estimate $\hat{\sigma}_{x,i}^{(t)}$, let the events

Applying condition Eq. 5 and the union bound, we get $P(\cup_{i\in{\cal X}\backslash\{x^{*}\}}E_{i,x})\leq\sum_{i\in{\cal X}\backslash\{x^{*}\}}E_{i,x}=\delta.$ So with a probability at least $1-\delta$ for any iteration $t<t^{*}$,

Therefore arm $x^{*}$ won’t be eliminated at iteration $t$. ∎

A.4 Proof of Prop. 7

Denote by $x^{(t)}$ the arm with the highest estimated mean at iteration $t$. First consider the case when Alg. 1 uses the marginal variance estimate $\hat{\sigma}_{i}^{(t)}$. With the condition in Eq. 5, it follows that $P(\cup_{i\in{\cal X}}E_{i})\leq\sum_{i\in{\cal X}}P(E_{i})\leq\delta$ where $E_{i}$ is defined in Eq. 15. So with a probability at least $1-\delta$,

Alg. 1 will stop by iteration $t$ if the RHS of the equation above satisfies the stopping criterion for all $i\neq x^{*}$, that is,

Solve the above inequality for $T^{(t)}$ and use the definition of the gap $\Delta$ we get

Since we use a doubling schedule $T^{(t)}=2T^{(t-1)}$ with $T^{(1)}=m^{(1)}$ and $T^{(t^{*})}=N$, Alg. 1 stops at an iteration no later than

And the total number of samples drawn by $t$ is upper bounded by $D(m^{(0)}2^{t-1}\wedge N)=T^{*}(\Delta)$.

Now consider the case when Alg. 1 uses the pairwise variance estimate $\hat{\sigma}_{x,i}^{(t)}$. With the condition in Eq. 5, it follows with the union bound that $P(\cup_{i\in{\cal X}\backslash\{x^{*}\}}E_{i})\leq\sum_{i\in{\cal X}\backslash\{x^{*}\}}P(E_{i})\leq\delta$ where $E_{i}$ is defined in Eq. 17. So with a probability at least $1-\delta$,

Now we can follow a similar argument as in the case with marginal variance estimate and prove the proposition. ∎

Appendix C Experiment Detailed Setting and Extra Results

The results with the marginal variance estimate $\hat{\sigma}_{i}$ for Racing are shown in Fig. 5. The Racing algorithms (both EBS and Normal) performs more conservatively compared to the plots when using pairwise variance estimate $\hat{\sigma}_{i,j}$ in Fig. 1, but the relative performance of all the algorithms are very similar to Fig. 1.

We also provide the results with $D=2$ and $D=100$ when Racing algorithms use pairwise variance estimate in Fig. 7 and 7 respectively. Racing-Normal performs the best in all situations and the empirical error never exceeds the provided bound $\delta$ with a statistical significance of $0.05$.

C.2 Details of the Bayesian ARCH Model Selection Experiment

The mixing rate of Carlin & Chib (1995) depends on a proper choice of the pseudoprior for $(\alpha_{i}^{(j)},\nu^{(j)})$. Ideally it should be similar to the parameter posterior when the model is chosen $p(\alpha_{i}^{(j)},\nu^{(j)})|q=j,{\bf r})$. We first reparameterize $(\alpha_{i}^{(j)},\nu^{(j)})$ with a softplus function $x=\log(1+\exp(x^{\prime}))$ to allow a full support along the real axis and then take the Laplace approximation at the MAP of transformed parameters as the pseudoprior for each model separately.

In order to avoid accessing the entire dataset each iteration, we use subsampling-based algorithms to sample all the conditionals except the pseudoprior as follows

where we sample $q$ with Racing-Normal Gibbs and sample $\boldsymbol{\alpha}^{(q)},\nu^{(q)}$ using MH with a proposal from SGLD and a rejection step provided by Racing-Normal MH. The rejection step controls the error introduced in SGLD when the step size is large.

We choose the step size separately for the exact and stochastic gradient Langevin dynamics (Welling & Teh, 2011) so that the acceptance rate is about 36%.

We apply the control variates by first segmenting the 2-D space of ${\bf z}_{j,t}\overset{\textnormal{def}}{=}(r_{t},\alpha_{0}^{(j)}+(\boldsymbol{\alpha}_{1:j}^{(j)})^{T}{\bf r}_{t-j:t-1})$, where $\boldsymbol{\alpha}^{(j)}$ takes the MAP value, equally into 100 bins according to marginal quantiles and then taking the reference points at the mean of each bin. We also notice that some data points have large residual reward $l_{i,n}-h_{i,n}$ when ${\bf z}_{j,t}$ is far from the reference point. We take 20% of the points with the largest distance in ${\bf z}$ as outliers, always compute them every iteration and apply the subsampling algorithm for the rest data.

C.3 Details of the Author Coreference Experiment

The main differences of this sampling problem from Eq. 1 are that

$|C_{y}|\neq|C_{y^{\prime}}|$ and the distribution of the cluster size follows approximately a power law with the value varying from as small as 1 to thousands. If we set $m^{(1)}=50$ as usual, we already draw about 33% of all the rewards in the first mini-batch. So we slightly abuse the Normal assumption and use a small size for $m^{(1)}=3$ and use doubling scheme for the rest with $m_{y}^{(2)}=(|C_{y}|-3)/10\wedge 1$. The experiment shows an empirical error $0.045$ of mis-identification of the best arm with the provided bound $\delta=0.05$.

The distribution of $\{f_{\theta}(x_{i},x_{j}):j\in C_{y}\}$ is independent from different clusters/arms. We exploit the independence of rewards and choose the bound

We obtained the dataset from the authors of Singh et al. (2012) but it is different from what is used in Singh et al. (2012) with more difficult citations. The best $B^{3}$ F-1 score reported in this paper is a reasonable value for this data set according to personal communications with the authors of Singh et al. (2012).