\section{Introduction}










Efficiently generating exact samples from a given target distribution, known up to normalization, has been a fundamental problem since the early days of algorithm design~\citep{marsaglia1963generating, walker1974new, kronmal1979alias, bratley2011guide, knuth2014art}. It is a basic building block of randomized algorithms and simulation, and understanding its theoretical limits is of intellectual and practical merit. Formally, let $p$ be a probability distribution on the set of integers $[N]:=\{1, \ldots, N\}$, and assume we are given query access to $\tilde p := Z p$, with an unknown constant $Z$. The Alias algorithm \citep{walker1974new} takes $\mathcal{O}(N)$ preprocessing time, after which one can repeatedly sample from $p$ in constant expected time. More recently, sophisticated algorithms have been devised to allow time-varying $\tilde p$ \citep{hagerup1993optimal, matias2003dynamic}, which also require $\mathcal{O}(N)$ preprocessing time. Unsurprisingly, for arbitrary $\tilde{p}$, the $\mc O(N)$ time is the best one can hope for, as shown in \cite{bringmann2017efficient} via a reduction to searching arrays. 


A common element of the aforementioned algorithms is the powerful idea of rejection. Rejection sampling, along with Monte Carlo simulation and importance sampling, can be traced back to the work of Stan Ulam and John von Neumann \citep{vonNeumann1951, Ulam}. As the name suggests, rejection sampling is an algorithm which proposes candidate samples, which are then accepted with a probability carefully chosen to ensure that accepted samples have distribution $p$. Despite the fundamental importance of rejection sampling in the applied sciences, there is surprisingly little work exploring its theoretical limits. In this work, we adapt the minimax perspective, which has become a staple of the modern optimization~\citep{nesterov2018lectures} and statistics~\citep{tsybakov2009nonparametric} literature, and we seek to characterize the number of queries needed to obtain rejection sampling algorithms with constant acceptance probability (e.g.\ at least $1/2$), uniformly over natural classes of target distributions.

We consider various classes of shape-constrained discrete distributions that exploit the ordering of the set $[N]$ (monotone, strictly unimodal, discrete log-concave). We also consider a class of distributions on the complete binary tree of size $N$, were only a partial ordering of the alphabet is required. For each of these classes, we show that the rejection sampling complexity scales sublinearly in the alphabet size $N$, which can be compared with the literature on sublinear algorithms~\citep{Gol10, Gol17}. This body of work is largely focused on statistical questions such as estimation or testing and the present paper extends it in another statistical direction, namely sampling from a distribution known only up to normalizing constant, which is a standard step of Bayesian inference.

To illustrate the practicality of our methods, we present an application to adversarial bandits \citep{bubeckcesabianchi2012bandits} and describe a variant of the classical \EXPthree\ algorithm whose iteration complexity scales as $\mathcal{O}(\log^2(K))$ where $K$ is the number of arms.

\section{Background on rejection sampling complexity}\label{scn:rej_sampling}

\subsection{Classical setting with exact density queries}

To illustrate the idea of rejection sampling, we first consider the classical setting where we can make queries to the \emph{exact} target distribution $p$. 
Given a \emph{proposal distribution} $q$ and an upper bound $M$ on the ratio $\max_{x\in [N]} p(x)/q(x)$, rejection sampling proceeds by
drawing a sample $X \sim q$ and a uniform random variable $U \sim \unif(0,1)$. If $U \leq p(X)/(Mq(X))$, the sample $X$ is returned; otherwise, the whole process is repeated. Note that the rejection step is equivalent to flipping a biased coin: conditionally on $X$, the sample $X$ is accepted with probability $p(X)/(Mq(X))$ and rejected otherwise. We refer to this procedure as \emph{rejection sampling with acceptance probability $p/(Mq)$}.

It is easy to check that the output of this algorithm is indeed distributed according to $p$. Since $Mq$ forms an upper bound on $p$, the region $G_q = \{(x, y)\,:\,x\in [N],\,y \in [0,Mq(x)]\}$ is a superset of $G_p = \{(x,y)\,:\,x\in [N],\,y \in [0,p(x)]\}$. Then, a uniformly random point from $G_q$ conditioned on lying in $G_p$ is in turn uniform on $G_p$, and so its $x$-coordinate has distribution $p$. A good rejection sampling scheme hinges on the design of a good proposal $q$ that leads to few rejections.

If $q=p$, then the first sample $X$ is accepted. More generally, 
the number of iterations required before a variable is accepted follows a geometric distribution with parameter $1/M$ (and thus has expectation $M$). In other words, the bound $M$ characterizes the quality of the rejection sampling proposal $q$, and the task of designing an efficient rejection sampling algorithm is equivalent to determining a strategy for building the proposal $q$ which guarantees a small value of the ratio $M$ using few queries.

\subsection{Density queries up to normalization}\label{scn:density_up_to_normalization}

In this paper, we instead work in the setting where we can only query the target distribution up to normalization, which is natural for Bayesian statistics, randomized algorithms, and online learning. Formally, let $\mc P$ be a class of probability distributions over a finite alphabet $\ms X$, and consider a target distribution $p\in \mc P$. We assume that the algorithm $\mc A$ has access to an oracle which, given $x \in \ms X$, outputs the value $Zp(x)$, where $Z$ is an unknown constant. The value of $Z$ does not change between queries. Equivalently, we can think of the oracle as returning the value $p(x)/p(x_0)$, where $x_0 \in \ms X$ is a fixed point with $p(x_0) > 0$. 

To implement rejection sampling in this query model, the algorithm must construct an \emph{upper envelope} for $\tilde p$, i.e., a function $\tilde q$ satisfying $\tilde q \ge \tilde p$. We can then normalize $\tilde q$ to obtain a probability distribution $q$. To draw new samples from $p$, we first draw samples $X \sim q$, which are then accepted with probability $\tilde p(X)/\tilde q(X)$. The following theorem shows that the well-known guarantees for rejection sampling also extend to our query model. The proof is provided in Appendix~\ref{scn:proof_of_rej_sampling}.

\begin{theorem}\label{thm:rejection_sampling}
    Suppose we have query access to the unnormalized target $\tilde{p} = p Z_p$ supported on $\ms X$, and that we have an upper envelope $\tilde q \ge \tilde p$.
    Let $q$ denote the corresponding normalized probability distribution, and let $Z_q$ denote the normalizing constant, i.e., $\tilde q = q Z_q$.
   
   
   
   
   
    Then, rejection sampling with acceptance probability $\tilde{p}/\tilde q$ outputs a point distributed according to $p$, and the number of samples drawn from $q$ until a sample is accepted follows a geometric distribution with mean $Z_q/Z_p$.
\end{theorem}

After $n$ queries to the oracle for $p$ (up to normalization), the output $\mc A(n, \tilde p)$ of the algorithm is an upper envelope $\tilde q \ge \tilde p$, and in light of the above theorem it is natural to define the \emph{ratio}
\begin{align*}
    r(\mc A, n, \tilde p)
    &:= \frac{Z_q}{Z_p}
    = \frac{\sum_{x\in \ms X} \tilde q(x)}{\sum_{x\in \ms X} \tilde p(x)} \,.
\end{align*}
The ratio achieved by the algorithm determines the expected number of queries to $\tilde p$ needed to generate each new additional sample from $p$.

As discussed in the introduction, our goal when designing a rejection sampling algorithm is to minimize this ratio uniformly over the choice of target $p\in \mc P$.
We therefore define the rejection sampling complexity of the class $\mc P$ as follows.

\begin{definition}\label{defn:complexity}
   For a class of distributions $\mc P$, \emph{the rejection sampling complexity of $\mc P$} is the minimum number $n\in \N$ of queries needed, such that there exists and algorithm $\mc A$ that satisfies
   \begin{align*}
       \sup_{\tilde p \in \tilde{\mc P}} r(\mc A, n,\tilde p) \le 2\,,
   \end{align*}
   where $\tilde{\mc P} := \{ \tilde p= Zp\,:\, Z >0\}$ is the set of all positive rescalings of distributions in $\mc P$.
\end{definition}
The constant $2$ in Definition~\ref{defn:complexity} is arbitrary and could be replaced by any number strictly greater than $1$, but we fix this choice at $2$ for simplicity.
With this choice of constant, and once the upper envelope is constructed, new samples from the target can be generated with a constant ($\leq 2$) expected number of queries per sample.

Note that when the alphabet $\ms X$ is finite and of size $N$, then $N$ is a trivial upper bound for the complexity of $\mc P$, simply by querying all of the values of $\tilde p$ and then returning the exact upper envelope $\mc A(N, \tilde p) = \tilde p$. Therefore, for the discrete setting, our interest lies in exhibiting natural classes of distributions whose complexity scales \emph{sublinearly} in $N$.

In this work, we specifically focus on \emph{deterministic} algorithms $\mc A$. In fact, we believe that adding internal randomness to the algorithm does not significantly reduce the query complexity.
Using Yao's minimax principle~\citep{yao1977minimax}, it seems likely that our lower bounds can extended to hold for randomized algorithms.
We leave this extension for future work.

\section{Results for shape-constrained discrete distributions}

In order to improve on the trivial rate of $\mc O(N)$ on an alphabet of size $N$, we need to assume some structure of the target distributions. A well-known set of structural assumptions are shape constraints \citep{groeneboom2014nonparametric, silvapulle2011constrained}, which have been extensively studied in the setting of estimation and inference. When the alphabet is $[N]$, shape constraints are built on top of the linear ordering of the support. We show that such assumptions indeed significantly reduce the complexity of the restricted classes of distributions to sublinear rates. We also consider the setting where the linear ordering of the support is relaxed to a partial ordering, and show it also results in sublinear complexity

Our complexity results for various classes of discrete distributions are summarized in~\autoref{tab:discrete}. We define the various classes below, and give the sublinear complexity algorithms that construct the upper envelopes in Figure~\ref{fig:algs}.

\begin{table}
  \caption{Rejection sampling complexities for classes of discrete distributions. Here, $N$ always denotes the alphabet size, $\ms X = \{1,\dotsc,N\}$.}
  \label{tab:discrete}
  \centering
  \begin{tabular}{ccccc}
    \toprule
    Class     & Definition     & Complexity & Theorem & Algorithm \\
    \midrule
    monotone & \autoref{defn:monotone}  & $\Theta(\log N)$  & \autoref{thm:monotone} & Algorithm~\ref{algo:monotone_upper}  \\
    strictly unimodal     & \autoref{defn:unimodal} & $\Theta(\log N)$  & \autoref{thm:unimodal} & Algorithm~\ref{alg:unimodal} \\
    cliff-like     & \autoref{defn:cliff}  & $\Theta(\log \log N)$ & \autoref{thm:cliff} & Algorithm~\ref{alg:log_concave} \\
    discrete log-concave     & \autoref{defn:lc}  & $\Theta(\log \log N)$ & \autoref{thm:discrete_lc} & Algorithm~\ref{alg:log_concave}\\
    monotone on a binary tree & \autoref{defn:tree} & $\Theta(N/\log N)$ & \autoref{thm:tree} & Algorithm~\ref{alg:tree} \\
    \bottomrule
  \end{tabular}
\end{table}

\begin{figure}
  \begin{algorithm}[H]
  \caption{Construct upper envelope for monotone distributions on $[N]$}\label{algo:monotone_upper}
  \begin{algorithmic}[1]
  \State Query the values $\tilde p(2^i)$, $0 \le i\le \lceil \log_2 N \rceil - 1$.
  \State Construct the upper envelope $\tilde q$ as follows: set $\tilde q(1) := \tilde p(1)$, and
  \begin{align*}
      \tilde q(x) := \tilde p(2^i) \, , \qquad \text{ for } x\in (2^i, 2^{i+1}]\,.
  \end{align*}
  \end{algorithmic}
\end{algorithm}
  \begin{algorithm}[H]
  \caption{Construct upper envelope for strictly unimodal distributions on $[N]$}\label{alg:unimodal}
  \begin{algorithmic}[1]
  \State Use binary search to find the mode of $\tilde p$.
  \State Use Algorithm~\ref{algo:monotone_upper} to construct an upper envelope on each side of the mode.
  \end{algorithmic}
\end{algorithm}
  \begin{algorithm}[H]
  \caption{Construct upper envelope for discrete log-concave distributions on $[N]$}\label{alg:log_concave}
  \begin{algorithmic}[1]
  \State Use binary search to find the first index $1 \le i \le \lceil \log_2 N \rceil$ such that $\tilde p(2^i) \le \tilde p(1)/2$, or else determine that $i$ does not exist.
  \State If $i$ does not exist, output the constant upper envelope $\tilde q \equiv \tilde p(1)$.
  \State Otherwise, output
  \begin{align*}
      \tilde q(x)
      := \begin{dcases} \tilde p(1)\,, & x < 2^{i}\,, \\ \tilde p(2^i) \exp\bigl[- \frac{\log(\tilde p(1)/\tilde p(2^i))}{2^i - 1} \, (x - 2^i) \bigr]\,, & x \ge 2^i\,. \end{dcases}
  \end{align*}
  \end{algorithmic}
\end{algorithm}
  \begin{algorithm}[H]
  \caption{Construct upper envelope for monotone distributions on binary trees of size $[N]$}\label{alg:tree}
  \begin{algorithmic}[1]
  \State Query $\tilde p(x)$ for all vertices $x$ which are at depth at most $\ell_0 := \ell - \lfloor \log_2\ell \rfloor + 1$, where $\ell$ is the maximum depth of the tree.
  \State Output
\begin{align*}
    \tilde q(x)
    &:= \begin{cases} \tilde p(x)\,, & \text{if}~\on{depth}(x) \le \ell_0\,, \\ \tilde p(y)\,, & \text{if}~\on{depth}(x) > \ell_0\,,\; \on{depth}(y) = \ell_0\,,\; \text{and}~x~\text{is a descendant of}~y\,. \end{cases}
\end{align*}
  \end{algorithmic}
\end{algorithm}
\caption{Algorithms for constructing rejection sampling upper envelopes which attain the minimax rates described in Table~\ref{tab:discrete}.}
\label{fig:algs}
\end{figure}

\subsection{Structured distributions on a linearly ordered set}

A natural class of discrete distributions which exploits the linear ordering of the set $[N]$ is the class of monotone distributions, defined below.

\begin{definition}\label{defn:monotone}
    The class of \emph{monotone} distributions on $[N]$ is the class of probability distributions $p$ on $[N]$ with $p(1) \ge p(2) \ge p(3) \ge \cdots \ge p(N)$.
\end{definition}





We show in~\autoref{thm:monotone} that the rejection sampling complexity of the class of monotone distributions is $\Theta(\log N)$, achieved via Algorithm~\ref{algo:monotone_upper}. It is also straightforward to extend Algorithm~\ref{algo:monotone_upper} to handle the class of strictly unimodal distributions defined next (see~\autoref{thm:unimodal} and Algorithm~\ref{alg:unimodal}).

\begin{definition}\label{defn:unimodal}
    The class of \emph{strictly unimodal} 
    distributions on $[N]$ is the class of probability distributions $p$ on $[N]$ such that: there exists a point $x \in [N]$ with $p(1) < p(2) < \cdots < p(x)$ and $p(x) > p(x+1) > \cdots > p(N)$.
\end{definition}

It is natural to ask whether further structural properties can yield even faster algorithms for sampling. This is indeed the case, and we start by illustrating this on a simple toy class of distributions.

\begin{definition}\label{defn:cliff}
    The class of \emph{cliff-like} distributions on $[N]$ is the class of probability distributions $\unif([N_0])$ for $N_0 \in [N]$.
\end{definition}

Since the class of cliff-like distributions is contained in the class of monotone distributions, Algorithm~\ref{algo:monotone_upper} yields a simple upper bound of $\mc O(\log N)$ for this class. However, we can do better by observing that in order to construct a good rejection sampling upper envelope for this class, we do not need to locate the index $N_0$ of the cliff exactly; it suffices to find it approximately, which in this context means finding an index $N_0'$ such that $N_0' \le N_0 \le 2N_0'$. Since we only need to search over $\mc O(\log N)$ possible values for $N_0'$, binary search can accomplish this using only $\mc O(\log \log N)$ queries. We prove in~\autoref{thm:cliff} that this rate is tight.

\begin{remark}
The class of cliff-like distributions provides a simple example of a class for which obtaining queries to the exact distribution is \emph{not} equivalent to obtaining queries for the distribution up to a normalizing constant.
Indeed, in the former model, the value of $p(1) = 1/N_0$ reveals the distribution in one query, implying a complexity of $\Theta(1)$, whereas we prove in~\autoref{thm:cliff} that the complexity under the second model is $\Theta(\log \log N)$.
\end{remark}

Instead of formally describing the algorithm for sampling from cliff-like distributions, we generalize the algorithm to cover a larger class of structured distributions: the class of \emph{discrete log-concave distributions} \cite[see \S 4]{saumard2014log}.


\begin{definition}\label{defn:lc}
    The class of \emph{discrete log-concave} distributions on $[N]$ is the class of probability distributions $p$ on $[N]$ such that for all $x \in \{2,\dotsc,N-1\}$, we have ${p(x)}^2 \ge p(x-1) p(x+1)$. Equivalently it is the class of distributions $p$ on $[N]$ for which there exists a convex function $V : \R\to\R \cup \{\infty\}$ such that $p(x) = \exp(-V(x))$ for all $x \in [N]$. In addition, we assume that the common mode of all of the distributions is at $1$.\footnote{Without this condition, the class of discrete log-concave distributions includes the family of all Dirac measures on $[N]$, and the rejection sampling complexity is then trivially $\Theta(N)$.}
\end{definition}

We prove in~\autoref{thm:discrete_lc} that the rejection sampling complexity of discrete log-concave distributions is $\Theta(\log \log N)$, achieved by Algorithm~\ref{alg:log_concave} (note that this algorithm also applies for cliff-like distributions, since cliff-like distributions are discrete log-concave).

\begin{remark}
The class of discrete log-concave distributions is another case for which rejection sampling with exact density queries is much easier than with queries up to a normalizing constant. In the former model, \cite{devroye1987simple} requires only a single query to construct a rejection sampling upper envelope with ratio $\le 5$. In contrast, we show in~\autoref{thm:discrete_lc} that the complexity under the second model is $\Theta(\log\log N)$.
\end{remark}

\subsection{Monotone on a binary tree}

The previous examples of structured classes all rely on the linear ordering of $[N]$. We now show that it is possible to develop sublinear algorithms when the linear ordering is relaxed to a partial ordering.
Specifically, we consider a structured class of distributions on balanced binary trees (note that the previously considered distributions can be viewed as distributions on a path graph).

\begin{definition}\label{defn:tree}
   The class of \emph{monotone distributions on a binary tree} with $N$ vertices is the class of probability distributions $p$ on a binary tree with $N$ vertices, with maximum depth $\lceil \log_2(N+1)\rceil$, such that for every non-leaf vertex $x$ with children $x_1$ and $x_2$, one has $p(x) \ge p(x_1) + p(x_2)$.
\end{definition}

We prove in~\autoref{thm:tree} that the rejection sampling complexity of this class is $\Theta(N/\log N)$; the corresponding algorithm is given as Algorithm~\ref{alg:tree}.

In a sense, \autoref{defn:tree} reduces to the class of monotone distributions when the underlying graph is a path, since each vertex in the (rooted) path graph has one ``child''.The reader may wonder whether replacing the condition $p(x) \ge p(x_1) + p(x_2)$ with $p(x) \ge p(x_1) \vee p(x_2)$ is more natural. In~\autoref{thm:alternate_monotone}, we show that rejection sampling cannot achieve sublinear complexity under the latter definition.

\section{Application to bandits}

Rejection sampling does not just provide us with a method for sampling from a target distribution; it provides us with the stronger guarantee of an upper envelope $\tilde q \ge \tilde p$, with a bound on the ratio of the normalizing constants of $\tilde q$ and $\tilde p$ (see Section~\ref{scn:density_up_to_normalization}). In this section, we show how this stronger property can be used to provide a faster, approximate implementation of the anytime variant of the \EXPthree\ algorithm. We expect that rejection sampling can yield similar computational speedups while retaining performance guarantees for other randomized algorithms.

Recall the adversarial bandit problem~\cite[Ch.\ 3]{bubeckcesabianchi2012bandits}: given $K$ arms, at each step $t \in [T]$ the player chooses an arm $I_t \in [K]$ to play. Simultaneously, an adversary chooses a loss vector $\ell_t \in [0,1]^K$. The chosen arm is then played, and the player incurs a loss of $\ell_t(I_t)$. The aim of the player is to find a strategy that minimizes the pseudo-regret, defined by
\begin{equation*}
    \overline{R}_n = \mathbb{E} \sum\limits_{t=1}^T \ell_t(I_t) - \min\limits_{k \in [K]} \mathbb{E} \sum\limits_{t=1}^T \ell_t(k)\,.
\end{equation*}
See Algorithm \ref{algo:EXP3} for the strategy known as  \EXPthree, which achieves a pseudo-regret of at most $2\sqrt{TK\log K}$ \cite[Theorem 3.1]{bubeckcesabianchi2012bandits}, which is minimax optimal up to the factor of $\sqrt{\log K}$~\cite[Theorem 3.4]{bubeckcesabianchi2012bandits}. In what follows $x(i)$ denotes the $i$'th coordinate of a vector $x$, and $e_j$ denotes the $j$'th standard basis vector in $\R^K$.
\begin{wrapfigure}[14]{L}{0.49\textwidth}
\vspace{.2cm}
\begin{minipage}{0.49\textwidth}
\begin{algorithm}[H]
  \caption{The  \EXPthree\ algorithm.}\label{algo:EXP3}
  \begin{algorithmic}[1]
    \Procedure{\EXPthree}{$T$, ${(\eta_t)}_{t=1}^T$}
    \State set $L_0 := 0$ and $p_0 := \unif\{1,\dots,K\}$
    \For{$t = 1,\dotsc, T$}
\State draw and play $I_t \sim p_{t-1}$\label{EXP3:draw}
\State observe loss $\ell_t(I_t)$
\State set $L_t := L_{t-1} + e_{I_t} \ell_t(I_t) / p_{t-1}(I_t)$\label{EXP3:loss update}
\State set $p_t \propto \exp(-\eta_t L_t)$\label{EXP3:prob update}
\EndFor
    \EndProcedure
  \end{algorithmic}
\end{algorithm}
\end{minipage}
\end{wrapfigure}

The computationally intensive steps of the iteration in Algorithm \ref{algo:EXP3} are drawing the sample on line \ref{EXP3:draw} and updating the distribution on line \ref{EXP3:prob update}. For each $t$, let us write $\tilde{p}_t = \exp(-\eta_t L_t)$ for the unnormalized version of $p_t$. Note that $\tilde{p}_t$ is fully determined by $L_t$ and $\eta_t$. Thus, if we can sample from $\tilde{p}_{t-1}$ on line \ref{EXP3:draw} in $o(K)$ time, then we can improve the na\"{\i}ve iteration complexity of $\Theta(K)$ since we can just skip line \ref{EXP3:prob update}. We achieve this by constructing a specialised data structure $\mathcal{D}$ that maintains the empirical loss vector $L$ in sorted order, thereby allowing fast sampling via Algorithm~\ref{algo:monotone_upper}. We record the requirements on $\mathcal{D}$ in the lemma below. 
\begin{lemma}\label{lem:data structure}
There exists a data structure $\mathcal{D}$ that stores a length-$K$ array $L$ and supports the following operations in $\mathcal{O}(\log K)$ worst-case time:
\begin{enumerate}
\item Given $(L[i], i)$ and a number $\ell$, set $L[i] = \ell$. \label{data structure op2}
\item Given $k \in [K]$, output the $k$-th largest element of the array $(L[k],k)_{k \in [K]}$ (in the dictionary order). \label{data structure op3}
\end{enumerate}
\end{lemma}

For a proof of~\autoref{lem:data structure}, see Section~\ref{sec:data structure proof}. Let us now describe a minor modification of the  \EXPthree\ algorithm which has (virtually) identical performance guarantees with per-iteration complexity $\mathcal{O}(\log^2 K)$. First, instead of sampling from $p_{t-1}$ directly, we do so using the rejection sampling proposal $q_{t-1}$ constructed from $\tilde p_{t-1}$ via our algorithm for monotone distributions (Algorithm~\ref{algo:monotone_upper}); this is possible because~\autoref{lem:data structure} gives us query access to the sorted version of $L_t$. Second, we modify the unbiased estimator of the loss in line~\ref{EXP3:loss update} accordingly.

The new unbiased estimator of the loss is defined as follows.
Draw another independent arm $J \sim q_t$ and replace line~\ref{EXP3:loss update} with
\begin{equation*}
    L_t := L_{t-1} + \frac{e_{I_t} \ell_t(I_t)}{\tilde{p}_{t-1}(I_t)} \, \frac{\tilde{p}_{t-1}(J)}{q_{t-1}(J)}\,.
\end{equation*}
Observe that if
$\mathcal{F}_t$ denotes the $\sigma$-algebra generated by all rounds up to time $t-1$ as well as the randomness of the adversary in step $t$, then
\begin{align*}
    \mathbb{E}\Bigl[ \frac{e_{I_t} \ell_t(I_t)}{\tilde{p}_{t-1}(I_t)} \, \frac{\tilde{p}_{t-1}(J)}{q_{t-1}(J)} \Bigm\vert \mathcal{F}_t\Bigr]
    &= \mathbb{E}\Bigl[ \frac{e_{I_t} \ell_t(I_t)}{p_{t-1}(I_t)} \Bigm\vert \mathcal{F}_t\Bigr] \mathbb{E}\Bigl[ \frac{p_{t-1}(J)}{q_{t-1}(J)} \Bigm\vert \mathcal{F}_t\Bigr]
    = \mathbb{E}\Bigl[ \frac{e_{I_t} \ell_t(I_t)}{p_{t-1}(I_t)} \Bigm\vert \mathcal{F}_t\Bigr]
    = \ell_t\,.
\end{align*}
The modified algorithm is given as Algorithm~\ref{algo:EXP3 mod}.

\begin{algorithm}[H]
  \caption{The modified version of the  \EXPthree\ algorithm.}\label{algo:EXP3 mod}
  \begin{algorithmic}[1]
    \Procedure{Fast-\EXPthree}{$T$, ${(\eta_t)}_{t=0}^{T-1}$}
    \State set $L_0 := 0$ and $p_0 := \unif\{1,\dots,K\}$
    \For{$t = 1,\dotsc, T$}
    \State build rejection sampling proposal $q_{t-1}$ of $\tilde{p}_{t-1}\label{EXP3 mod:envelope} :=\exp(-\eta_{t-1}L_{t-1})$
    \State draw and play $I_t \sim p_{t-1}$ via rejection sampling using $q_{t-1}$\label{EXP3 mod:sample}
    \State draw $J \sim q_{t-1}$ independently\label{EXP3 mod: sample2}
\State observe loss $\ell_t(I_t)$
\State set $L_t := L_{t-1} + e_{I_t} \ell_t(I_t) \tilde{p}_{t-1}(J)/ (\tilde{p}_{t-1}(I_t) q_{t-1}(J))$\;\label{EXP3 mod:loss update}
\EndFor
    \EndProcedure
  \end{algorithmic}
\end{algorithm}

Let us verify the claimed $\mc O(\log^2 K)$ per-iteration complexity of Algorithm~\ref{algo:EXP3 mod}. Let $\mathcal{D}$ be an instance of the data structure described in~\autoref{lem:data structure} and initialize it with $L = L_0$. Building the rejection envelope $q_{t-1}$ on line \ref{EXP3 mod:envelope} requires $\mc O(\log K)$ calls to operation \ref{data structure op3}, for a total complexity of $\mathcal{O}(\log^2 K)$. Sampling from $q_{t-1}$ on line \ref{EXP3 mod:sample} requires $\mathcal{O}(\log K)$ time and performing a rejection step takes one call to operation \ref{data structure op3}, so the expected complexity of this step is $\mathcal{O}(\log K)$. Drawing $J$ on line \ref{EXP3 mod: sample2} requires $\mc O(\log K)$ time and one call to operation \ref{data structure op3}, and finally, line \ref{EXP3 mod:loss update} requires one call to operation \ref{data structure op2}.

The following result (proven in Appendix~\ref{scn:proof_regret}) provides a pseudo-regret guarantee.

\begin{proposition}\label{prop:EXP3 mod performance}
Algorithm \ref{algo:EXP3 mod} with step size $\eta_t := \frac 12 \sqrt{\frac{\log K}{K(t+1)}}$ satisfies
\begin{equation*}
    \overline{R}_n \leq 4\sqrt{TK\log K}\,. 
\end{equation*}
\end{proposition}
We regard the above result as a proof of concept for the use of rejection sampling to more efficiently implement subroutines in randomized algorithms.
Our proof of~\autoref{prop:EXP3 mod performance} follows well-known arguments from the bandit literature, and the key new ingredient is our strong control of the ratio between the target and proposal distribution, which allows us to bound the variance of our unbiased estimator of the loss. 


\begin{remark}\label{rem:exp3 reducing variance}
In Algorithm \ref{algo:EXP3 mod} one may replace $\tilde{p}_{t-1}(J)/q_{t-1}(J)$ by $(1/m) \sum_{i=1}^m \tilde{p}_{t-1}(J_i)/q_{t-1}(J_i)$ for i.i.d.\ $J_i \sim q_{t-1}$ in order to further reduce the variance of the estimator. 
\end{remark}


\begin{figure}
\centering
\begin{minipage}{0.47\textwidth}
  \centering
  \includegraphics[width=\linewidth]{EXP3-cliff.jpg}
  \caption{Error bars denote $4$ standard deviations over $20$ runs.}
  \label{fig:changing-cliff regret}
\end{minipage}%
\hfill
\begin{minipage}{0.47\textwidth}
  \centering
  \includegraphics[width=\linewidth]{speed-comparison.jpg}
  \caption{Comparison of iteration speed.}
  \label{fig:EXP3 speed main}
\end{minipage}
\end{figure}

We conduct a small simulation study to confirm that Algorithm \ref{algo:EXP3 mod} is competitive with \EXPthree. In ~\autoref{fig:changing-cliff regret} we plot the regret of the two algorithms over $T=20\text{k}$ steps using the theoretical step size $\eta_t = \sqrt{\log K / (K\,(t+1))}$ and $m=5$ (see~\autoref{rem:exp3 reducing variance}). We run the algorithms on a toy problem with $K=256$ arms, where $10\%$ of the arms always return a loss of $0$, and the remaining arms always return the maximal loss of $1$. In~\autoref{fig:EXP3 speed main} we compare the time it takes for the two algorithms to complete an iteration on the same problem, but with varying number of arms $K$. In the case of \EXPthree,\, when $\log_2 K > 20$ the values are not computed, and for $\log_2 K \leq 20$ they are extrapolated from $100$ iterations, as the running time becomes prohibitive for such large experiments. 

\begin{remark}
We note that for constant step size $\eta_t \equiv \eta$ (e.g.\ when the time horizon $T$ is known in advance), theoretically it is possible to implement \EXPthree\, with constant iteration cost by applying the results of \cite{matias2003dynamic}. Further, using the `doubling trick' such an algorithm can be turned into an anytime algorithm by restarting it $\approx \log T$ times at a computational cost of $\mathcal{O}(K)$ each time. However, our Algorithm \ref{algo:EXP3 mod} is the first to implement the more elegant solution of decaying step size in sublinear time per iteration.
\end{remark}



\section{Conclusion and outlook}

We studied the query complexity of rejection sampling within a minimax framework, and we showed that for various natural classes of discrete distributions, rejection sampling can obtain exact samples with an expected number of queries which is sublinear in the size of the support of the distribution. Our algorithms can also be run in sublinear time, which make them substantially faster than the baseline of multinomial sampling, as shown in our application to the  \EXPthree\ algorithm.

A natural direction for future work is to investigate the complexity of rejection sampling on other structured classes of distributions, such as distributions on graphs, or distributions on continuous spaces. In many of these other settings, the complexity of algorithms based on Markov chains has been studied extensively, but the complexity of rejection sampling remains to be understood.

\section*{Acknowledgments}

Sinho Chewi was supported by the Department of Defense (DoD) through the National Defense Science \& Engineering Graduate Fellowship (NDSEG) Program. Thibaut Le Gouic was supported by NSF award IIS-1838071. Philippe Rigollet was supported by NSF awards IIS-1838071, DMS-1712596,  and DMS-2022448.

