\section{Introduction}
Selecting $k$ diverse items out of a larger collection of $n$ items
is a classical problem in computer science which naturally emerges in many
tasks such as summarization (select $k$ phrases) and recommendation (select
$k$ articles/ads to show to the user). An increasingly popular approach
to model and quantify diversity in this
subset selection problem is that of determinantal point processes (DPPs).
Given a set $[n] \stackrel{\text{\tiny{def}}}{=} \{1, \ldots, n\}$ of $n$ items and a target
size $k$, one can define a DPP of size $k$ (known as a $k$-DPP)
through an $n\times n$ posivite semi-definite (PSD) similarity matrix
$\mathbf{L}$ (also known as the kernel matrix).
The matrix $\mathbf{L}$ encodes the similarities between items, and the user must choose
it so that $[\mathbf{L}]_{ij}$ is larger the more the $i$-th and $j$-th items are similar.
Given $k$ and $\mathbf{L}$, we define $S \sim k$-${\mathrm{DPP}}(\mathbf{L})$
as a distribution over all ${n \choose k}$ index subsets
$S\subseteq [n]$ of size $k$, such that $\Pr(S) \propto \det(\mathbf{L}_S)$
is proportional to the determinant of the sub-matrix $\mathbf{L}_S$ induced by
the subset. DPPs have found numerous applications in machine
learning, not only for summarization \cite{dpp-summarization,dpp-video,dpp-salient-threads,pmlr-v80-celis18a}
and recommendation \cite{dpp-shopping,NIPS2018_7805}, but also in experimental
design \cite{bayesian-experimental-design,symmetric-polynomials},
stochastic optimization \cite{dpp-minibatch,randomized-newton},
Gaussian Process optimization \cite{NIPS2016_6452},
low-rank approximation
\cite{pca-volume-sampling,more-efficient-volume-sampling,nystrom-multiple-descent}, and
more (recent surveys include
\cite{dpp-ml,dpp-stats,dpp-randnla}). Note that early work on DPPs
focused on a \emph{random-size} variant, which we
denote $S\sim {\mathrm{DPP}}(\mathbf{L})$, where the subset size is allowed to take any value
between $0$ and $n$, and the role of parameter $k$ is replaced by the expected size $\mathbb E[|S|]=
d_{\textnormal{eff}}(\mathbf{L}) \stackrel{\text{\tiny{def}}}{=} \mathrm{tr}\left(\mathbf{L}(\mathbf{L} + \mathbf I)^{-1}\right)$.
The quantity $d_{\textnormal{eff}}(\mathbf{L})$ is known in randomized linear algebra
\cite{ridge-leverage-scores,dpp-randnla} and learning theory
\cite{caponnetto2007optimal} as the effective dimension. While
random-size DPPs exhibit deep connections to many scientific domains
\cite{dpp-independence}, the \emph{fixed-size} $k$-DPPs are typically
more practical from a machine learning stand-point \cite{k-dpp}.

\textbf{Sampling from a $k$-DPP.} The first $k$-DPP samplers
scaled poorly, as they all relied on an eigendecomposition \cite{k-dpp}  of
$\mathbf{L}$ taking $\mathcal{O}(n^3)$ time.
Replacing the eigendecomposition with a Cholesky factorization can increase numerical stability \cite{dpp-noeig},
and empirical performance \cite{poulson2019high-performance} thanks to
dynamically-scheduled, shared-memory parallelizations, but still
ultimately require $\mathcal{O}(n^3)$ time. A number of methods have
been proposed which use approximate eigendecomposition
\cite{dpp-coreset,dpp-nystrom} to reduce the computational cost,
however these approaches provide limited guarantees on the accuracy of sampling.

To improve scalability, several approaches
based on Monte-Carlo sampling were introduced, using rejection or Gibbs sampling.
The fastest MCMC sampler for  $k$-DPPs, to the best of our knowledge, is
by \cite{rayleigh-mcmc} and has  $n\cdot{\mathrm{poly}}(k)$
complexity, i.e., asymptotically much faster than the cost of eigendecomposition.
However these MCMC methods do \emph{not} sample exactly from the $k$-DPP distribution,
and can only guarantee that the final sample will be close in distribution
to a $k$-DPP. Moreover these guarantees only hold \emph{after mixing}, which is difficult
to verify and requires at least $\mathcal{O}(nk^2)$ time, making MCMC methods not applicable
when $n$ is large.

\begin{wrapfigure}{R}{.44\textwidth}
  \centering
\vspace{-2\baselineskip}
    \begin{tabular}{r|l}
      &Complexity
      \\
              \hline
      \cite{k-dpp,dpp-noeig,poulson2019high-performance,more-efficient-volume-sampling}
&$n^3$\\[1mm]
      \textsc{DPP-VFX}\xspace \cite{derezinski2019exact}& $n\cdot k^{10} + k^{15}$\\[1mm]
$\alpha$-DPP (this paper)& $(\beta n\cdot k^6+k^9)\sqrt k$
    \end{tabular}
    \captionof{table}{\small Runtime comparison of exact $k$-DPP
      sampling algorithms. Here, $\beta\leq 1$ is the fraction of
      items observed by $\alpha$-DPP (see Theorem~\ref{thm:informal-main}).}
    \vspace{-3mm}
    \label{tab:comp}
\end{wrapfigure}

A recent line of works \cite{derezinski2019exact,dpp-intermediate},
using the ideas from \cite{leveraged-volume-sampling,correcting-bias},
developed sampling algorithms
specially designed for a \emph{random-size} DPP (as opposed to a $k$-DPP), which avoid
expensive decomposition of the kernel while sampling exactly from $S\sim{\mathrm{DPP}}(\mathbf{L})$.
In particular, they showed that it is sufficient to first choose an intermediate subset
$\sigma \subseteq [n]$ sampled i.i.d.\,from the \emph{marginal}
distribution of the DPP,  i.e.,\xspace $\mathbf P(i \in \sigma) \approx \mathbf P(i \in S)$, and
then sample from a DPP restricted to the items indexed by
$\sigma$. Since the size of $\sigma$ is typically much less
than $n$, this leads to a more efficient algorithm.
Note that rescaling ${\mathrm{DPP}}(\mathbf{L})$ into ${\mathrm{DPP}}(\alpha\mathbf{L})$ using some constant~$\alpha$
only changes the expected size of $S$ from $d_{\textnormal{eff}}(\mathbf{L})$ to $d_{\textnormal{eff}}(\alpha\mathbf{L})$.
By accurately choosing an appropriate $\alpha_{\star}$, one can boost the probability that the random
size of $S$ is exactly $k$, and convert a DPP sampler into a $k$-DPP sampler
by repeatedly sampling $S \sim {\mathrm{DPP}}(\alpha_{\star}\mathbf{L})$ until $S$ has
size~$k$. Based on this reduction, \citet{derezinski2019exact}
gave the first algorithm (\textsc{DPP-VFX}\xspace) which is capable of exact sampling
from a $k$-DPP in time $n\cdot{\mathrm{poly}}(k)$. However,
when sampling from $k$-DPPs, the approach of
\cite{derezinski2019exact} has \emph{two major limitations}:
\begin{enumerate}[leftmargin=*]
  \item \textsc{DPP-VFX}\xspace has an $\Omega(n)$ runtime bottleneck, since it requires
    computing all $n$ marginals, one for each item, in
    order to define the i.i.d.\/ distribution of $\sigma$, which may
    be infeasible for very large $n$.
\item
  The reduction used by \cite{derezinski2019exact} to convert a DPP
  sampler into a $k$-DPP sampler
   increases the time complexity by a factor
  of at least $k^4$, resulting in a $\acO(n\cdot k^{10}+k^{15})$ runtime.
\end{enumerate}

In this paper, we address both of these limitations by introducing a
new algorithm called $\alpha$-DPP, which 1) does not need to compute all
of the marginals, and 2) uses a new efficient reduction to convert
from a random-size DPP to a fixed-size $k$-DPP (see Table
\ref{tab:comp} for comparison).

\textbf{Main contribution: uniform intermediate sampling for $k$-DPPs.}
To resolve the $\Omega(n)$ runtime bottleneck, we use an additional
intermediate sample $\rho$ based on \emph{uniform} sub-sampling.
Since uniform sampling can be implemented without looking at the actual items
in the collection, this means that we do not even have to look at any
item outside of $\rho$.
The only necessary assumption required by our approach is that the maximum
entry (i.e.,\xspace similarity) of $\mathbf{L}$ is bounded by a constant $\kappa^2$. However,
to simplify exposition we also assume w.l.o.g.~that $d_{\textnormal{eff}}(\mathbf{L}) \geq k$ (see \Cref{s:binary-search}).

In particular, we 1) sample $\rho$ uniformly out of $[n]$,
then 2) only approximate the marginal probabilities of items in $\rho$
to compute $\sigma$, and finally 3) downsample $\sigma$ into
a DPP sample $S$. To guarantee that $S$ is distributed exactly according to the DPP
it is crucial that $\rho$ is diverse enough. We show
that sampling a $k^2/d_{\textnormal{eff}}(\mathbf{L})$ fraction of $[n]$ into $\rho$ (i.e.,\xspace $|\rho| \approx k^2/d_{\textnormal{eff}}(\mathbf{L}) \cdot n$)
is enough. Since all the expensive computation is performed only
on $\rho$, this gives us a $d_{\textnormal{eff}}(\mathbf{L})/k^2$ speedup
over existing methods.
\vspace{2mm}

\begin{theorem}\label{thm:informal-main}
Given any $\mathbf{L} \succeq \mathbf{0}$ with $\max_{ij} \mathbf{L}_{ij} \leq
\kappa^2$ and $1 \leq k \leq d_{\textnormal{eff}}(\mathbf{L})$, there exists an algorithm that returns
$S \sim k$-${\mathrm{DPP}}(\mathbf{L})$, and with probability $1 - \delta$ runs in time
$$\acO\big((\beta n\cdot k^6+k^9)\sqrt k\log(1/\delta)\big),$$
where $\beta \leq\min\big\{k^2\kappa^2 /d_{\textnormal{eff}}(\mathbf{L}),1\big\}$ is the fraction
of items observed by the algorithm.
\end{theorem}
In the derivation of \Cref{thm:informal-main} we make several novel contributions.
First, we provide a DPP sampler that
given $\mathbf{L}$ and a rescaling $\alpha \leq 1$ leverages a mixture of uniform and rejection sampling
to sample from ${\mathrm{DPP}}(\alpha\mathbf{L})$ observing only an $\alpha\kappa^2k$ fraction of
the items. We then show that the optimal rescaling $\alpha_{\star}$ required
by the reduction from $k$-DPP to DPP can be bounded with $\alpha_{\star} \leq \mathcal{O}(k/d_{\textnormal{eff}}(\mathbf{L}))$,
and thus our rescaling-aware sampler can sample from $k$-DPPs looking
only at a $k^2/d_{\textnormal{eff}}(\mathbf{L})$ fraction of the items. Finally, we provide
an efficient search algorithm to find a close approximation $\hat{\alpha}$
of $\alpha_{\star}$.

\textbf{Model misspecification and computational free lunch.}
Our result can be also interpreted from a perspective of model misspecification.
Note that every time the users define a $k$-DPP they
also implicitly define a random size ${\mathrm{DPP}}(\mathbf{L})$.
Moreover, the natural expected sample size (i.e.,\xspace implicit number of unique items in $[n]$)
of ${\mathrm{DPP}}(\mathbf{L})$ is $d_{\textnormal{eff}}(\mathbf{L})$,
which does not depend on the desired size~$k$.
Therefore, if $\mathbf{L}$ is not chosen appropriately
$d_{\textnormal{eff}}(\mathbf{L})$ might be much larger than $k$, and the $k$-DPP is selecting
$k$ unique items out of a much larger implicit pool of $d_{\textnormal{eff}}(\mathbf{L}) \gg k$ unique items.
In this case, it is possible to consider only a small $k^2/d_{\textnormal{eff}}(\mathbf{L})$ fraction of the items
selected uniformly at random and still have enough unique items to sample a diverse $k$-subset.
Our result shows for the first time that it is possible to take advantage of this modeling disagreement between $k$ and $d_{\textnormal{eff}}(\mathbf{L})$ to gain computational
savings while still sampling \emph{exactly} from the DPP, i.e.,\xspace a computational free lunch.

\textbf{Binary search reduction from k-DPP to DPP.} Both our approach and the one of \citet{derezinski2019exact} rely on first implementing
an efficient random-size DPP sampler, followed by the usage of a black-box construction
based on rejection sampling to transform the DPP sampler into a $k$-DPP sampler.
However the reduction of \citet{derezinski2019exact} requires access to a
high-precision estimate of $d_{\textnormal{eff}}(\alpha\mathbf{L})$ in order to appropriately
tune $\alpha$. This makes optimizing $\alpha$ the bottleneck in the reduction from $k$-DPP to DPP,
and therefore there is a large computational gap between the two problems.
We close this gap thanks to a novel approach to find a suitable rescaling $\alpha$
based not on optimization but rather on binary search. Crucially, to find a suitable $\alpha$
this approach does not require an estimate of $d_{\textnormal{eff}}(\alpha\mathbf{L})$, but only $\mathcal{O}(\sqrt{k}\log(n))$ black-box calls
to a DPP sampler. Therefore, it can transform any random size DPP sampler
into a $k$-DPP sampler with only a $\sqrt{k}$ overhead,
and could be applied to any future improved sampler beyond this paper.

\section{Sampling from a rescaled DPP with intermediate uniform subsampling}
\label{s:alpha-sampler}
In this section we focus on a specific class of DPPs, $S \sim {\mathrm{DPP}}(\alpha\mathbf{L})$,
specified using a rescaling $\alpha \leq 1$
and a similarity matrix $\mathbf{L}$, which we refer to as rescaled DPPs. The main result
of the section is showing that a sufficiently large subset selected uniformly
at random can be used as an intermediate sample to sample from a
rescaled DPP without looking at all of the items.
The main reason to focus on rescaled DPPs is because they naturally appear when
reducing $k$-DPP sampling to DPP sampling, where rescaling is used to align
the random size of the DPP and $k$. This is going to be the focus of the next section.
However the approach proposed in this section
is not limited to rescaled DPPs, but under the right assumptions can be extended
to accelerate sampling from generic DPPs. We will discuss these extensions
at the end of the section.\vspace{-1mm}
\paragraph{Notation}
We use $[n]$ to denote the set $\{1,\dots,n\}$. For a matrix
$\mathbf B\in\mathbb R^{n\times m}$ and index sets $C$, $D$, we use $\mathbf B_{C,D}$
to denote the submatrix of $\mathbf B$ consisting of the intersection of rows
indexed by $C$ with columns indexed by $D$.  If $C=D$, we use a
shorthand $\mathbf B_{C}$ and if $D=[m]$, we may write $\mathbf B_{C,[m]}$.
Finally, we also allow $C,D$ to be multisets or
sequences, in which case each row/column is duplicated in the matrix
according to its multiplicity (and in the case of sequences, we
order the rows/columns as they appear in the sequence).
Note that with this notation if $\mathbf{L}=\mathbf B\B^\top$ then
$\mathbf{L}_{C,D}=\mathbf B_{C,[n]}\mathbf B_{D,[n]}^\top$.

\vspace{-2mm}
\subsection{Background: distortion-free intermediate sampling.}\label{subs:background}
Rather than sampling directly from the target DPP, intermediate
sampling \cite{dpp-intermediate,leveraged-volume-sampling}
first selects an intermediate subset $\sigma$ from $[n]$, and then refines it
by extracting $S$ from $\sigma$. Crucially, if $\sigma$ is selected according
to a so-called Regularized DPP (${\mathrm{R\textnormal{-}DPP}}$), this is equivalent to sampling $S$ from a DPP.
\begin{definition}\label{d:r-dpp}
  For any psd matrix $\mathbf{L} \in \mathbb R^{n\times n}$, distribution
  $p\stackrel{\text{\tiny{def}}}{=} \{p_i\}_{i=1}^n$ and $r>0$, define $\widetilde{\mathbf{L}} \in \mathbb R^{n\times n}$ with~$\widetilde{\mathbf{L}}_{i,j} \stackrel{\text{\tiny{def}}}{=} \frac{\mathbf{L}_{i,j}}{r\sqrt{p_ip_j}}$. We define an ${\mathrm{R\textnormal{-}DPP}}_p^r(\mathbf{L})$
  as distribution over events $A\subseteq\bigcup_{k=0}^\infty[n]^k$ such that\vspace{-2mm}
  \begin{align*}
    \Pr(A) \stackrel{\text{\tiny{def}}}{=}\mathbb E_{\sigma}\big[\mathbf 1_{[\sigma\in A]}\det(\mathbf I +
    \widetilde{\mathbf{L}}_\sigma)\big]/{\det(\mathbf I + \mathbf{L})}\mathbin{\raisebox{0.5ex}{,}}\quad\text{for}
    \quad\sigma=(\sigma_1,\dots,\sigma_t)\simiid p,\quad
    t\sim\mathrm{Poisson}(r).
  \end{align*}
\end{definition}
\begin{proposition}[{\citealp[Theorem 8]{dpp-intermediate}}]\label{t:composition}
  For any $\mathbf{L}$, $p$, $r,$ and $\widetilde{\mathbf{L}}$ defined as in Definition
  \ref{d:r-dpp},\vspace{-1mm}
  \begin{align*}
    \text{if}\quad \sigma\sim{\mathrm{R\textnormal{-}DPP}}_p^r(\mathbf{L})\quad\text{and}\quad
    S\sim{\mathrm{DPP}}(\widetilde{\mathbf{L}}_\sigma)\quad\text{then}\quad
    \{\sigma_i:i\!\in\! S\}\sim{\mathrm{DPP}}(\mathbf{L}).
  \end{align*}
\end{proposition}
A computationally inefficient but conceptually simple approach to rejection sample $\sigma$ is the following:\\*
1) compute all marginals $\mathbf P(i \in S) = \ell_i(\mathbf{L}) \stackrel{\text{\tiny{def}}}{=} [\mathbf{L}(\mathbf I+\mathbf{L})^{-1}]_{i}$ and sum to $\sum_{i=1}^n\ell_i(\mathbf{L}) = d_{\textnormal{eff}}(\mathbf{L})$ \cite{ridge-leverage-scores};\\*
2) sample $t \sim \text{Poisson}(c)$ and $\sigma \sim \text{Multinomial}\left(t, \frac{\ell_1(\mathbf{L})}{d_{\textnormal{eff}}(\mathbf{L})}, \ldots, \frac{\ell_n(\mathbf{L})}{d_{\textnormal{eff}}(\mathbf{L})}\right)$ for an appropriate constant $c$;\\*
3) accept $\sigma$ w.p.~$\tfrac{\det(\mathbf I + \wt{\mathbf{L}}_{\sigma})}{C\det(\mathbf I + \mathbf{L})}$, where $C$ is an appropriate constant used to make the rejection step valid.

All existing intermediate sampling algorithms \cite{dpp-intermediate,derezinski2019exact,leveraged-volume-sampling,correcting-bias} rely on this approach,
refining it to make use of efficient approximations of the marginals
$\ell_i(\mathbf{L})$ and adapting the constants $c$ and $C$ to the data. However they all share a common bottleneck:
to sample $\sigma$ i.i.d.~they need to approximate all marginals $\ell_i(\mathbf{L})$
and the normalization constant $d_{\textnormal{eff}}(\mathbf{L})$, and therefore the final
runtime scales as $n\cdot{\mathrm{poly}}(k)$.
While this is much smaller than the $\mathcal{O}(n^3)$ required by an exact sampler,
it still becomes quickly unfeasible when $n$ is very large. In what follows
we will introduce another approach to sample from an ${\mathrm{R\textnormal{-}DPP}}$ that does not
require to approximate the marginals of all items, but only the items selected
in a preliminary \emph{uniform} intermediate sample.

\subsection{Faster DPP sampling with uniform intermediate sampling}

\begin{wrapfigure}[21]{R}{.45\textwidth}
   \vspace{-2\baselineskip}
 \begin{minipage}{1.02\linewidth}
 \begin{algorithm}[H]
   \caption{\textsc{$\alpha$-DPP}\xspace sampler}
 \label{alg:alpha-dpp-sampler}
 \begin{algorithmic}[1]
 \renewcommand{\algorithmicrequire}{\textbf{Input:}}
 \REQUIRE $\alpha$, $\mathbf{L}$, $\mathcal{D}$, $\mathbf W$, $r \geq 1$
\STATE Set $\widehat{\mathbf{L}} = \mathbf W^{1/2}\mathbf{L}_{\mathcal{D},\mathcal{D}}\mathbf W^{1/2} \in \mathbb R^{m \times m}$
\REPEAT
\STATE Sample $u \sim \text{Poisson}(re^{1/r}\alpha n\kappa^2)$
\STATE Sample $\rho = \text{Uniform}(u, [n])$
\FOR{$j = \{1, \dots, u\}$}
\STATE Compute $l_{\rho_j}$ using Eq.~(\ref{eq:rls-approx})
\STATE
Sample $z_j \sim \text{Bernoulli}(l_{\rho_j}/(\alpha\kappa^2))$\label{line:reject-binomial}\hspace*{-1cm}
\ENDFOR
\STATE Set $\sigma = \{\rho_j: z_j = 1\}$, $t = |\sigma|$
\STATE Set $[\widetilde{\mathbf{L}}_{\sigma}]_{ij} = \tfrac{1}{r\sqrt{l_{\sigma_i} l_{\sigma_j}}}[\mathbf{L}]_{\sigma_i\sigma_j}$
\STATE $\textit{Acc}\sim\!
   \text{Bernoulli}\Big(\tfrac{{\mathbb E}^{\tilde{s}}\det(\mathbf I+ \alpha\widetilde{\mathbf{L}}_\sigma)}{{\mathbb E}^{t/r}\det(\mathbf I+\alpha\widehat{\mathbf{L}})}
 \Big)$ \label{line:acc}
 \UNTIL{$\textit{Acc}=\text{true}$\label{line:rep2}}
 \STATE Sample $\widetilde{S}\sim {\mathrm{DPP}}\big(\alpha\widetilde{\mathbf{L}}_\sigma\big)$\label{line:sub}
   \RETURN $S = \{\sigma_i: i\!\in\! \widetilde{S}\}$
 \end{algorithmic}
 \end{algorithm}
\end{minipage}
\end{wrapfigure}
We now introduce our novel $\alpha$-rescaled DPP sampler, called
\textsc{$\alpha$-DPP}\xspace (see \Cref{alg:alpha-dpp-sampler}).
It requires as input a rescaling $\alpha$, a similarity matrix $\mathbf{L}$ and a parameter
$r$ that will be used to tune the Poisson sampling step of \Cref{t:composition}
approach.
It also requires as input a \emph{dictionary} $\mathcal{D}$ containing $m$ elements,
and set of weights stored in a diagonal matrix $\mathbf W \in \mathbb R^{m \times m}$.
A dictionary is a subset of items $\mathcal{D} \subseteq [n]$
such that reweighting the items in $\mathcal{D}$ by $\mathbf W$ provides a good approximation
of $\mathbf{L}$, so that the approximate marginals
\begin{align}\label{eq:rls-approx}
l_i \stackrel{\text{\tiny{def}}}{=} \alpha[\mathbf{L} - \alpha\mathbf{L}_{[n],\mathcal{D}}^\top(\alpha\mathbf{L}_{\mathcal{D}} + \mathbf W^{-1})^{-1}\mathbf{L}_{[n],\mathcal{D}}]_{i}
\vspace{-.25\baselineskip}
\end{align}
computed using $\mathcal{D}$ and $\mathbf W$ are close to the true marginals $\ell_i$ (see \Cref{a:misc-proofs}).
Compared to the meta-approach of \Cref{t:composition}, the main technical difference
is that rather than sampling directly $t$ from an appropriate Poisson,
and then $\sigma$ from a Multinomial, we introduce an intermediate uniform sampling
step. In particular, we first sample a Poisson $u$, and then uniformly sample
a subset $\rho$ containing $u$ items. We then compute an approximation
$l_i$ of the marginal $\ell_i$ only for the items in $\rho$,
and downsample $\rho$ into $\sigma$ using rejection sampling (\Cref{line:reject-binomial}).
Finally, we accept or reject $\sigma$ (\Cref{line:acc})
and then downsample $\sigma$ into $S$ using a standard DPP sampler on
the smaller $\wt{\mathbf{L}}_{\sigma}$.

\Cref{alg:alpha-dpp-sampler} is not simply a different
implementation of the approach of \Cref{t:composition}, since even if Multinomial sampling is implemented with lazy evaluations of $l_i$,
we would still need to compute the \emph{normalization} constant of the Multinomial,
which strictly requires computing all $l_i$. Similarly, the rejection test of \Cref{line:acc} is also designed to accept as many candidates as possible without requiring the computation of the normalization constant as in \cite{derezinski2019exact}. Rather
our approach is a novel method to sample from an ${\mathrm{R\textnormal{-}DPP}}$
using Poisson rejection sampling. In particular, we prove not only that
it always returns an $S$ sampled according to the exact DPP distribution,
but also that if the dictionary satisfies certain conditions,
the main of which is $(\varepsilon,\alpha)$-accuracy (see \Cref{a:misc-proofs} and \cite{calandriello_disqueak_2017}),
then the algorithm will generate $S$ quickly.
\begin{theorem}\label{thm:alpha-sampling}
Given any $\mathbf{L} \succeq \mathbf{0}$, dictionary $\mathcal{D}$, $\mathbf W\succ \mathbf{0}$, $r \geq 1$ and $\alpha > 0$, \textsc{$\alpha$-DPP}\xspace returns
$S \sim {\mathrm{DPP}}(\alpha\mathbf{L})$. Moreover, if $r \geq d_{\textnormal{eff}}(\alpha\mathbf{L}) \geq 1/2$, $\mathcal{D}$ and $\mathbf W$ are $(1/d_{\textnormal{eff}}(\alpha\mathbf{L}), \alpha)$-accurate,
$\mathcal{D}$ satisfies $|\mathcal{D}| \leq 10d_{\textnormal{eff}}(\alpha\mathbf{L})$, and $d_{\textnormal{eff}}(\alpha\widehat{\mathbf{L}}) \leq 10d_{\textnormal{eff}}(\alpha\mathbf{L})$, w.p.~$1 - \delta$ \textsc{$\alpha$-DPP}\xspace runs in time\vspace{-1mm}
\begin{align*}
  \mathcal{O}\Big(\left[ \min\{\alpha\kappa^2d_{\textnormal{eff}}(\alpha\mathbf{L}), 1\}\cdot n \cdot
  d_{\textnormal{eff}}(\alpha\mathbf{L})^6\log^2(n/\delta) +
  d_{\textnormal{eff}}(\alpha\mathbf{L})^9\log^3(n/\delta)\right ] \cdot \log(1/\delta)\Big).
  \end{align*}
\end{theorem}
\vspace{-1mm}

The main implication of our result is that the intermediate distribution based
on marginals can be replaced more and more accurately with a uniform distribution
as $\alpha$ becomes smaller. This results in having to compute marginals
only for a $\min\{\alpha\kappa^2d_{\textnormal{eff}}(\alpha\mathbf{L}), 1\}$ fraction of the $n$ items.
This speedup can be significant when the rescaling $\alpha$ is very small,
as is the case when we want to sample a small number of items out of a large
collection.
Compared to other exact DPP samplers, such as \textsc{DPP-VFX}\xspace,
our \textsc{$\alpha$-DPP}\xspace is strictly faster by roughly a $1/(\alpha\kappa^2d_{\textnormal{eff}}(\alpha\mathbf{L}))$ factor
when implemented with an appropriate caching strategy for the
estimates $l_i$ (see \Cref{a:misc-proofs}).
Further, unlike MCMC samplers, \textsc{$\alpha$-DPP}\xspace is an \emph{exact} sampler. Moreover, there is no
known MCMC approach that can achieve a runtime sub-linear in $n$ when $\alpha$
is small as \textsc{$\alpha$-DPP}\xspace.

An $(\varepsilon, \alpha)$-accurate dictionary
that also satisfies the other conditions can be generated using
a slight modification of the \textsc{BLESS}\xspace algorithm \cite{NIPS2018_7810}, that we call \textsc{BLESS-I}\xspace algorithm, presented
in \Cref{a:poisson-bless}.
However, note that since the marginals $\ell_i$ are equivalent to
the ridge leverage scores \cite{ridge-leverage-scores} of item $i$,
we can replace \textsc{BLESS-I}\xspace with any present or future algorithm for leverage score sampling that can
be modified to be rescaling-aware \cite{calandriello_disqueak_2017,NIPS2018_7810}.
Moreover, note that \textsc{BLESS-I}\xspace also returns an estimate
of $d_{\textnormal{eff}}(\alpha\mathbf{L})$ that is sufficiently accurate to tune $r$ and $\varepsilon$.
At the same time, our analysis could be excessively conservative, and instead of trying
to set $r$ and $\varepsilon$ using $d_{\textnormal{eff}}(\alpha\mathbf{L})$ as suggested by \Cref{thm:alpha-sampling},
a more practical strategy is to start with a constant
$r$ and increase it slowly if the sampler is rejecting with a too low probability,
using a doubling schedule to preserve overall time complexity.

\textbf{Proof sketch.} The proof is divided in two parts, proving that \textsc{$\alpha$-DPP}\xspace is an exact sampler
(\Cref{lem:a-dpp-exact}) and that under the right conditions it is efficient (\Cref{lem:a-dpp-efficient}).

For the first part we once again rely on the approach of \Cref{t:composition},
but with the added difficulty of not being allowed to compute all the marginals.
To avoid this bottleneck, we show that:\vspace{-1mm}
\begin{enumerate}
\item[A)] sampling $t \sim \text{Poisson}(r)$ and $\sigma \sim
  \text{Multinomial}(t, \{\ell_i/d_{\textnormal{eff}}(\alpha\mathbf{L})\}_{i=1}^n)$;\quad and
\item[B)] sampling $n$ independent $s_i \sim \text{Poisson}(r'\ell_i(\alpha\mathbf{L}))$, and adding $s_i$ copies of item $i$ to $\sigma$,
\end{enumerate}\vspace{-1mm}
are equivalent for an appropriate choice of $r$ and $r'$, i.e.,\xspace
we prove that the $\sigma$ generated by both approach A and B follow the same distribution.
However, unlike approach A, approach B does not require computing
a normalization constant, i.e.,\xspace it samples from \emph{unnormalized} probabilities.
Moreover, if we know an upper bound on the marginals we can further reduce the number of marginals that need
to be computed. In our case we use
the bound $\ell_i \leq \alpha\kappa^2$, and show that
\vspace{-1mm}
\begin{enumerate}
  \item[C)] sampling $n$ Poisson independently $u_i \sim
\text{Poisson}(r'\alpha\kappa^2)$, only if $u_i > 0$ computing
$\ell_i$ and sampling $s_i \sim \text{Binomial}(u_i,
\ell_i/(\alpha\kappa^2))$, and adding $s_i$ copies of item $i$ to
$\sigma$
\end{enumerate}\vspace{-1mm}
once again generates $\sigma$ strictly equivalent to the ones of approach B and A.
The added advantage of approach C over the others is that only the marginals of items with $u_i > 0$
are actually computed, and there is no need to compute a normalization constant.
Starting from this new approach, to obtain our \textsc{$\alpha$-DPP}\xspace sampler (\Cref{alg:alpha-dpp-sampler})
we simply replace the $n$ Poisson $u_i \sim \text{Poisson}(r'\alpha\kappa^2)$ with a single
$u \sim \text{Poisson}(r'\alpha\kappa^2 n)$ followed by uniform sampling,
and replace the exact $\ell_i$ with approximate $l_i$.

For the second part we derive a lower bound on the acceptance probability
similar to the one from \citet{derezinski2019exact}. However, while they
use an $n \times d_{\textnormal{eff}}(\alpha\mathbf{L})$ Nystr\"om\xspace approximation of the matrix $\mathbf{L}$,
to avoid direct dependencies on $n$ we are forced to use a less stable
approximation $\widehat{\mathbf{L}}$. As a result, controlling $d_{\textnormal{eff}}(\alpha\widehat{\mathbf{L}})$
requires a more careful analysis.

\textbf{Beyond uniform subsampling.} One of the implications of our analysis
is that more adaptive upper bounds on the marginals $\ell_i$ could further
speedup our \textsc{$\alpha$-DPP}\xspace sampling approach.
In particular, we chose uniform sampling, i.e.,\xspace a uniform upper bound,
for its conceptual simplicity and because knowing an upper bound $\kappa^2$
on the entries of the similarity matrix usually does not require looking at
the items, e.g.,\xspace $\kappa^2$ is always equal to 1 for Gaussian similarity, Cosine similarity or other
self normalized similarities. However for other similarities, such as linear
similarity, this bound could be very loose. A simple replacement is using
the actual diagonal of $\mathbf{L}$, which requires to look at all items and $\mathcal{O}(n)$
time to compute but is usually very scalable. Ideally, one could imagine
designing a sequence of upper bounds starting from cheaper to more computationally
expensive, where more advanced techniques such as random projection are used
near the end to further filter candidate items.
 \section{Efficient reduction from k-DPP to rescaled DPP via binary search}
\label{s:binary-search}

Given our fast DPP sampler, we can see a $k$-DPP as a sampling process where we first sample
$S \sim {\mathrm{DPP}}(\alpha\mathbf{L})$, check if the sample size $|S|$ is equal to $k$, and then
accept or reject the sample accordingly. Rescaling
$\mathbf{L}$ by a constant factor $\alpha$ only changes the expected size
$d_{\textnormal{eff}}(\alpha\mathbf{L})$ (and not the $k$-DPP), with
$\alpha > 1$ increasing the expected size and $\alpha < 1$ decreasing it.
Thus, it is natural to imagine that there exists some $\alpha_\star$ for which
the acceptance probability is high. Indeed this was recently proven to be possible.
\citet{derezinski2019exact} show that if the \emph{mode}
$m_{\alpha_{\star}}$ of $S_{\alpha_\star}\sim{\mathrm{DPP}}(\alpha_\star \mathbf{L})$ is equal to
$k$, then we will accept with probability at least $\Omega(1/\sqrt{k})$.
They also provide an
algorithm to find such an $\alpha_{\star}$. However,
this algorithm has a prohibitively high computational cost,
$\acO(nk^{10} + k^{15})$, because ensuring that the mode of
${\mathrm{DPP}}(\alpha_\star\mathbf{L})$ is \emph{exactly} $k$ requires an extremely
accurate approximation of $\mathbf{L}$. Instead, our approach is to run a binary search
to find a good rescaling $\alpha$, which will terminate once the
acceptance probability is high enough, regardless of whether $k$ is
exactly the mode. Crucially, this binary search only
requires a black box ${\mathrm{DPP}}(\alpha\mathbf{L})$ sampler (such as
our \textsc{$\alpha$-DPP}\xspace), and it only queries the sampler $\acO(\sqrt k)$ many
times. To prove that the binary search finds a good $\alpha$ in a
small number of steps, we establish a new property (Lemma~\ref{l:pb-new}) of the Poisson
Binomial distribution (the distribution of the subset sizes of
${\mathrm{DPP}}(\alpha\mathbf{L})$), which should be of independent interest.

\subsection{Binary search}
Our main result in this subsection is \Cref{alg:binary-search}, which requires only oracle
access to the samples from a random-size DPP, and finds a rescaling
$\hat\alpha$ which enables efficient rejection sampling from a
$k$-DPP. Note that the provided oracle sampler does not have to be
our \textsc{$\alpha$-DPP}\xspace sampler, so the algorithm could be paired with other samplers.

\begin{lemma}[{restate=[name=restated]binarysearch}]\label{l:binary-search}
  Suppose that we are given an integer $k$, a range
  $I=[\alpha_{\min},\alpha_{\max}]$ where
  $\alpha_{\max}=\gamma\alpha_{\min}$, and access
  to an oracle which, for any
  $\alpha\in I$, returns $S\sim
  {\mathrm{DPP}}(\alpha \mathbf{L})$. If there exists $\alpha_\star\in I$ such that $k$
  is the mode of $|S|$ for $S\sim{\mathrm{DPP}}(\alpha_\star \mathbf{L})$, then using
  $O\big(\sqrt k\log^2(k\log(\gamma)/\delta)\big)$
  calls to the oracle we can find $\hat\alpha\in I$ such that with
  probability $1-\delta$ we have
  \begin{align*}
    \Pr(|S| = k) = \Omega\big(\tfrac1{\sqrt k}\big),\qquad\text{for}\quad
    S\sim{\mathrm{DPP}}(\hat\alpha \mathbf{L}).
  \end{align*}
\end{lemma}\vspace{-2mm}
The distribution of subset size $|S|$ for $S\sim{\mathrm{DPP}}(\mathbf{L})$
can be defined via the eigenvalues $\lambda_1\geq\lambda_2\geq ...$ of
$\mathbf{L}$ (see \cite{dpp-independence}): if we let $b_i\sim
\mathrm{Bernoulli}(\frac{\lambda_i}{\lambda_i+1})$ for $i\geq 1$, then
$\sum_i b_i$ is distributed identically to $|S|$. This distribution is known as the \emph{Poisson
  Binomial}, and it has been extensively studied in the probability
literature \cite{pb-review}. The recent result of
\cite{derezinski2019exact} on the probability of the mode of a Poisson
Binomial shows that it is possible to find $\hat\alpha$ satisfying the
condition of Lemma \ref{l:binary-search}.
\begin{lemma}[{restate=[name=restated]pbmode}]\label{l:pb-mode}
   There is an absolute constant $0<c<1$ such that for any Poisson
   Binomial distribution $p:\mathbb Z_{\geq 0}\rightarrow\mathbb R_{\geq 0}$,
with mode $k^*$ we have $p(k^*)\geq \frac{c}{\sqrt{k^*\!+1}}\cdot$
\end{lemma}\vspace{-2mm}
This result, however, does not provide an efficient way of
finding an $\hat\alpha$ such that the mode of the subset size distribution of
${\mathrm{DPP}}(\hat\alpha\mathbf{L})$ is $k$. We circumvent this problem by performing
a binary search (\Cref{alg:binary-search}) that looks for such an
$\hat\alpha$, but stops early when it reaches a
sufficiently good candidate, avoiding excess computations. To make this rigorous, we establish the
following new property of the Poisson Binomial distribution, which
should be of independent interest.
\begin{lemma}[{restate=[name=restated]pbnew}]\label{l:pb-new}
  Let $p:\mathbb Z_{\geq 0}\rightarrow\mathbb R_{\geq 0}$ be a Poisson Binomial distribution,
  and let $k\geq 1$ satisfy
$p(k)<\frac{c}{12\sqrt{3(k+1)}}\mathbin{\raisebox{0.5ex}{,}}$ where $c$ comes from Lemma
\ref{l:pb-mode}. Then, $P_{<k}=\sum_{i<k}p(k)$ and
$P_{>k}=\sum_{i>k}p(k)$ satisfy:
\begin{enumerate}[leftmargin=1cm]
  \vspace{-3mm}
\item if the mode of $p$ is less than $k$, then
$P_{>k}\leq \frac12 - \frac c{12}$; \vspace{-2mm}
\item if the mode of $p$ is greater than $k$, then $P_{<k}\leq \frac12
  - \frac{c}{12}\cdot$
\end{enumerate}
\end{lemma}\vspace{-1mm}
Informally, the above result states the following:
For any $k$, either its probability under the given Poisson Binomial
is at least $\Omega(\frac1{\sqrt k})$, or this $k$ splits the probability mass into two
uneven parts, with the larger one containing the mode. Thus, as long
as our candidate $\alpha$ does not yield high acceptance probability
for $k$, it is easy to make the branching decision in the binary search by estimating the
quantities $P_{>k}$ and $P_{<k}$ simply by repeated sampling from
${\mathrm{DPP}}(\alpha\mathbf{L})$. Note that if the condition on $p(k)$ is not
satisfied, then performing the branching decision could be very
expensive, but our algorithm avoids this possibility. The proof of
\Cref{l:binary-search} (\Cref{a:bin-search}) follows from Lemmas
\ref{l:pb-mode} and \ref{l:pb-new}.

\setlength{\textfloatsep}{5pt}
\begin{algorithm}[t]
  \caption{Binary search for initializing the $k$-${\mathrm{DPP}}(\mathbf{L})$ sampler}
\label{alg:binary-search}
{
\setstretch{1.4}
\begin{algorithmic}[1]
  \renewcommand{\algorithmicrequire}{\textbf{Input:}}
  \renewcommand{\algorithmicensure}{\textbf{Output:}}
\REQUIRE $0<\alpha_{\min}<\alpha_{\max}$, sampling oracle for
${\mathrm{DPP}}(\alpha \mathbf{L})$, integer $k$ and constants $C>0, \delta\in(0,1)$
\ENSURE $\hat\alpha$ such that ${\mathrm{DPP}}(\hat\alpha \mathbf{L})$ can be used
to efficiently sample $k$-${\mathrm{DPP}}(\mathbf{L})$
\FOR{$s = \{1, \ldots, \lceil\log(\gamma)\rceil\}$}
\STATE \textbf{if} {$\alpha_{\max}/\alpha_{\min} < (1+\frac1{(k+3)^2})$\label{line:if1}} \textbf{then return} $\hat\alpha = \alpha_{\min}$
\STATE $\bar\alpha \leftarrow\sqrt{\alpha_{\min}\alpha_{\max}}$
\STATE Sample $S_1,...,S_t\simiid {\mathrm{DPP}}(\bar\alpha \mathbf{L})$\quad
where\quad $t = C\sqrt k\log(s/\delta)$
\STATE $\hat P_k \leftarrow \frac1t \sum_{i=1}^t\mathbf 1_{[|S_i|=k]}$
\STATE \textbf{if} {$\hat P_k\geq \frac12\cdot \frac{c}{12\sqrt{3(k+1)}}$\label{line:if2}} \textbf{then return}
$\hat\alpha = \bar\alpha$
\STATE  $(\hat P_{<k}, \hat P_{>k}) \leftarrow \big(\frac1t
\sum_{i=1}^t\mathbf 1_{[|S_i|<k]},\ \frac1t
\sum_{i=1}^t\mathbf 1_{[|S_i|>k]}\big)$

\vspace{.15cm}
\STATE \textbf{if} $\hat P_{<k}> \hat P_{>k}$ \textbf{then} $(\alpha_{\min},\alpha_{\max}) = (\bar\alpha,\alpha_{\max})$ \textbf{else} $(\alpha_{\min},\alpha_{\max}) = (\alpha_{\min},\bar\alpha)$
\label{line:branch}
\ENDFOR
\end{algorithmic}
}
\end{algorithm}

\subsection{Constructing the initial interval}

To initiate our binary search, we must first find a range of values
$[\alpha_{\min},\alpha_{\max}]$, which contains the desired $\alpha_\star$, and also
construct a sampling oracle for ${\mathrm{DPP}}(\alpha\mathbf{L})$.
The binary search procedure is deliberately presented in a way that is agnostic to
how these two steps are accomplished, because a number of existing DPP
samplers could be adapted to take advantage of
\Cref{alg:binary-search}, including
\cite{dpp-noeig,poulson2019high-performance,dpp-intermediate,derezinski2019exact}.
Our implementation of these two steps is different than these previous
approaches in that it takes advantage of the structure of the kernel
so that it only has to look at a potentially small fraction of the
data points. We achieve this with a modified version
of the BLESS algorithm \cite{NIPS2018_7810}.

\begin{lemma}[{restate=[name=restated]algovalidinterval}]\label{lem:algo-valid-interval}
  W.p.~$1-\delta$ \textsc{BLESS-I}\xspace runs in time
 $\acO\left(\min\{\alpha_{\max}\kappa^2, 1\}nk^6 + k^9\right)$
  and satisfies:
  \vspace{-.5\baselineskip}
  {
  \setlength{\itemsep}{-2pt}
  \begin{enumerate}[leftmargin=1cm]
    \item The interval $[\alpha_{\min},\alpha_{\max}]$ is bounded by $\frac14(k-1)/\mathrm{tr}(\mathbf{L}) \leq \alpha_{\min} \leq \alpha_{\max} \leq 8(k + 2)/d_{\textnormal{eff}}(\mathbf{L})$
  \item    There is $\alpha_\star\in[\alpha_{\min},\alpha_{\max}]$
    for which $k$ is the mode of $|S|$ where
    $S\sim{\mathrm{DPP}}(\alpha_\star\mathbf{L})$;
    \item The dictionary $\mathcal{D}^{\alpha_{\max}}$ satisfies the
      conditions from \Cref{thm:alpha-sampling} for any $\alpha\in[\alpha_{\min},\alpha_{\max}]$.
  \end{enumerate}
  }
\end{lemma}

The first two parts of the lemma ensure that the interval
$I=[\alpha_{\min},\alpha_{\max}]$ is a valid input for
the binary search in \Cref{alg:binary-search} and that its size
$\gamma = \alpha_{\max}/\alpha_{\min} \leq 4\mathrm{tr}(\mathbf{L})/d_{\textnormal{eff}}(\mathbf{L})$
is bounded in the log-scale. The last part implies that
$\alpha$-DPP can be used by that algorithm as the oracle sampler.

At a high level, \Cref{alg:adjusted-bless} proceeds by starting with a
small $\alpha^0$ that is guaranteed to be a valid lower bound for the interval,
and for which a dictionary $\mathcal{D}^0$ can be constructed simply via
uniform sampling. Then we repeatedly double the $\alpha$ and refine
the dictionary, until we reach $\alpha^i$ such that we can ensure that
with high probability $d_{\textnormal{eff}}(\alpha^i\mathbf{L}) \geq k + 1$ which makes it
a valid upper bound for the interval (then, this $\alpha^i$
becomes $\alpha_{\max}$).

\subsection{Overall time complexity of $k$-DPP sampling}
Putting together all the results from the previous sections, we can
finally bound the computational complexity of our $k$-DPP sampler,
which first uses \textsc{BLESS-I}\xspace (\Cref{alg:adjusted-bless}) to construct a dictionary and search interval,
and then applies the binary search of (\Cref{alg:binary-search}) using
our $\alpha$-DPP sampler (\Cref{alg:alpha-dpp-sampler}) as the sampling oracle.
Once again note that in the following computational analysis we will use conservative
values for many parameters, notably $r$ from \textsc{$\alpha$-DPP}\xspace and $q$ from \textsc{BLESS-I}\xspace,
as they are suggested from the theory. However in practice it is always better to
start from a more optimistic value, and keep doubling them only if the sampler
repeatedly fails to accept. Importantly, samples generated this way will still
be exactly distributed according to the DPP, as all the approximations used
in our approach only influence the runtime of our algorithm, and not the correctness of
its acceptance, which always holds.

By \Cref{lem:algo-valid-interval}, the preprocessing step of running
\textsc{BLESS-I}\xspace takes $\acO\left(\min\{\alpha_{\max}\kappa^2,1\}nk^6 + k^9\right)$
and generates a dictionary $\mathcal{D}$ with size $\acO(k^3)$.
Since $d_{\textnormal{eff}}(\alpha\mathbf{L}) \leq d_{\textnormal{eff}}(\alpha_{\max}\mathbf{L}) \leq \mathcal{O}(k)$ for all
$\alpha$ in the search interval, each call to the \textsc{$\alpha$-DPP}\xspace sampler also requires
at most $\acO\left(\min\{\alpha_{\max}\kappa^2k, 1\}nk^6 + k^9\right)$.
Finally, the binary search invokes $\alpha$-DPP at most $\acO(\sqrt{k})$ times
so the overall runtime is
$\acO\big((\min\{\alpha_{\max}\kappa^2k, 1\}nk^{6} +
    k^{9})\cdot\sqrt{k}\big)$.
We now provide a bound on $\alpha_{\max}$.
\begin{lemma}[{restate=[name=restated]lambdalowerbound}]\label{lem:lambda-lower-bound}
For any matrix $\mathbf{L}$ and $0 < \alpha \leq 1$, we have
$d_{\textnormal{eff}}(\alpha\mathbf{L})/d_{\textnormal{eff}}(\mathbf{L}) \geq \alpha \geq d_{\textnormal{eff}}(\alpha\mathbf{L})/\mathrm{tr}(\mathbf{L})$.
\end{lemma}
Applied to $\alpha_{\max}$, we obtain $\alpha_{\max} \leq d_{\textnormal{eff}}(\alpha_{\max}\mathbf{L})/d_{\textnormal{eff}}(\mathbf{L}) \leq
\mathcal{O}(k/d_{\textnormal{eff}}(\mathbf{L}))$, giving us the final runtime of
$\acO\big((\min\{k^2\kappa^2/d_{\textnormal{eff}}(\mathbf{L}), 1\}nk^{6} + k^{9})\cdot\sqrt{k}\big)$
reported in \Cref{thm:informal-main}.

\noindent\begin{minipage}{0.49\textwidth}
\centering
\includegraphics[width=.95\textwidth]{fig/fig_runtime_repro_noresample_first_k_dpp.pdf}
\vspace{-3mm}
\captionof{figure}{Small scale experiment}\label{fig:small-scale}
\vspace{-2mm}
\end{minipage}\hfill
\begin{minipage}{0.49\textwidth}
\centering
\includegraphics[width=.95\textwidth]{fig/fig_runtime_repro_k_succ_rbf_time.pdf}
\vspace{-3mm}
\captionof{figure}{Large scale experiment.}\label{fig:large-scale}
\vspace{-2mm}
\end{minipage}
\section{Experiments}
\vspace{-.75\baselineskip}
In this section, we evaluate our \textsc{$\alpha$-DPP}\xspace sampler on a
benchmark\footnote{\url{https://github.com/LCSL/dpp-vfx}}
introduced by \cite{derezinski2019exact} (see
\Cref{a:experiments}). The benchmark uses
subsets of the infinite MNIST dataset \cite{loosli-canu-bottou-2006}
with $d = 784$ and $n$ varying up to $10^6$.
All experiments are executed on a 28 core Xeon E5-2680 v4. Each experiment is repeated multiple times, and we report mean
values and a 95\% confidence interval.

\textbf{Baselines:} we compare \textsc{$\alpha$-DPP}\xspace with \textsc{DPP-VFX}\xspace \cite{derezinski2019exact}, an MCMC sampler \cite{rayleigh-mcmc}
and a sampler based on eigendecompositions \cite{dpp-independence, gillenwater2014approximate}.
All algorithms  are implemented in \texttt{python} as part
\ifisarxiv
\else
\footnote{We provide our modified version of the library which includes \textsc{$\alpha$-DPP}\xspace in the supplementary material.}
\fi
of \texttt{DPPy} \cite{gautier2018dppy}.
Due to their similar input, we use the same oversampling parameters (see \Cref{a:experiments}) for \textsc{$\alpha$-DPP}\xspace and \textsc{DPP-VFX}\xspace.
We run the MCMC sampler for $\mathcal{O}(nk)$ iterations to guarantee mixing \cite{rayleigh-mcmc}. For more details
on hyperparameter tuning we refer to \Cref{a:experiments}.

\textbf{Results} \ We begin by reporting results on a smaller subset of data (\Cref{fig:small-scale})
where even the non-efficient samplers can be run. We use an \texttt{rbf}
similarity with $\sigma = \sqrt{3d}$, and set $k=10$ to match the number of
digit classes in MNIST. Note that for $n = 70000$ \textsc{BLESS-I}\xspace estimates
$d_{\textnormal{eff}}(\mathbf{L}) \approx 300$, validating our assumption of $d_{\textnormal{eff}}(\mathbf{L}) \gg k$.
Thanks to this mismatch, we can see how \textsc{$\alpha$-DPP}\xspace maintains a \emph{constant}
runtime as $n$ grows, and increasingly matches or outpaces competing baselines as $n$ grows.
In particular, it becomes faster than the eigendecomposition based sampler (which cannot scale beyond $n = 24000$)
or the MCMC sampler. However, the gap is still sufficiently small that
\textsc{DPP-VFX}\xspace, the previously fastest $k$-DPP sampler
available, remains competitive.

For larger datasets we consider only the scalable samplers, \textsc{$\alpha$-DPP}\xspace and
\textsc{DPP-VFX}\xspace. We consider again an \texttt{rbf} similarity, but this time we choose $n$ up
to $10^6$ and $\sigma = \sqrt{10}$. This further increases the gap between $k$ and $d_{\textnormal{eff}}(\mathbf{L})$,
with \textsc{BLESS-I}\xspace estimating $d_{\textnormal{eff}}(\mathbf{L}) \approx 1000$.
We report results in \Cref{fig:large-scale}, with runtime shown in
log-scale. In this regime, the gap between \textsc{DPP-VFX}\xspace and
\textsc{$\alpha$-DPP}\xspace widens, as \textsc{DPP-VFX}\xspace
cannot use rescaling to reduce the final dictionary size from $d_{\textnormal{eff}}(\mathbf{L})$ to
$d_{\textnormal{eff}}(\hat\alpha\mathbf{L}) \approx k$, and has to compute $n$ marginal probabilities
since it does not leverage uniform intermediate subsampling.
In particular, thanks to the uniform sampling step, we see that \textsc{$\alpha$-DPP}\xspace's
runtime does not grow as $n$ grows, since all the expensive computations
are performed in the small intermediate subset which is hardly
sensitive to $n$. We note that, due to using a smaller dictionary,
\textsc{$\alpha$-DPP}\xspace requires about 2-5x more trials in the rejection sampling
step, which leads to larger variance in the runtime.

\begin{wrapfigure}[12]{R}{0.43\textwidth}
  \begin{minipage}{\linewidth}
    \vspace{-5mm}
\centering
\includegraphics[width=.95\textwidth]{fig/fig_runtime_repro_k_succ_rbf_fraction.pdf}
\vspace{-2mm}
\captionsetup{font=small}
\captionof{figure}{Fraction of items observed by \textsc{$\alpha$-DPP}\xspace.}\label{fig:rejections}
\vspace{-4mm}
    \end{minipage}
    \end{wrapfigure}

In \Cref{fig:rejections}, we report the fraction of data that
is observed by \textsc{$\alpha$-DPP}\xspace in the large scale experiment. This
quantity, denoted as $\beta$ in \Cref{thm:informal-main}, is
responsible for much of the computational gains of the algorithm over
\textsc{DPP-VFX}\xspace, reported in \Cref{fig:large-scale}. Note that the remaining
$1-\beta$ portion of the data does not ever need to be loaded into
the program's memory, which leads to a significant reduction in
memory accesses. We observe that as the
data size increases, the fraction of items observed by \textsc{$\alpha$-DPP}\xspace goes down to as
little as $1$\% for $n=10^6$, which is why the
runtime of \textsc{$\alpha$-DPP}\xspace stays roughly flat, whereas the runtime of
$\textsc{DPP-VFX}\xspace$ grows.

\newpage
\section*{Broader impact}
DPPs were discovered in the 70s by Odile Macchi to model repulsion of
particle distributions in fermions, so improvements in samplers may help
in modelling physical simulations. In bringing faster DPP samplers to
machine learning we aim to enable a
better handling of diversity through this rigorous theoretical framework.
\begin{ack}
  MD thanks the NSF for funding via the NSF TRIPODS program.
\end{ack}

\bibliographystyle{plainnat}