Exact sampling of determinantal point processes with sublinear time preprocessing

Introduction

Our method is based on a technique developed recently by and later extended by . In this approach, we carefully downsample the index set $[n]=\{1,...,n\}$ to a sample $\sigma=(\sigma_{1},...,\sigma_{t})\in[n]^{t}$ that is small but still sufficiently larger than the expected target size $k$, and then run a DPP on $\sigma$. As the downsampling distribution we use a regularized determinantal point process (R-DPP), proposed by , which (informally) samples $\sigma$ with probability $\Pr(\sigma)\sim\det(\mathbf{I}+\widetilde{\mathbf{L}}_{\sigma})$, where $\widetilde{\mathbf{L}}$ is a rescaled version of $\mathbf{L}$. Overall, the approach is described in the diagram below, where $|S|\leq t\ll n$,

Given a psd matrix $\mathbf{L}$, its $i$th $\lambda$-ridge leverage score (RLS) $\tau_{i}(\lambda)$ is the $i$th diagonal entry of $\mathbf{L}(\lambda\mathbf{I}+\mathbf{L})^{-1}$. The $\lambda$-effective dimension $d_{\textnormal{eff}}(\lambda)$ is the sum of the leverage scores, $\sum_{i}\tau_{i}(\lambda)$.

Intuitively, if the marginal probability of $i$ is high, then this index should likely make it into the intermediate sample $\sigma$. This suggests that i.i.d. sampling of the indices $\sigma_{1},...,\sigma_{t}$ proportionally to 1-ridge leverage scores, i.e. $\Pr(\sigma_{1}=i)\propto\tau_{i}(1)$, should serve as a reasonable and cheap heuristic for constructing $\sigma$. In fact, we can show that this distribution can be easily corrected by rejection sampling to become the R-DPP that we need. Computing ridge leverage scores exactly costs $\mathcal{O}(n^{3})$, so instead we compute them approximately by first constructing a Nyström approximation of $\mathbf{L}$.

Let $\mathbf{L}$ be a psd matrix and $C$ a subset of its row/column indices with size $m\triangleq|C|$. Then we define the Nyström approximation of $\mathbf{L}$ based on $C$ as the $n\times n$ matrix $\widehat{\mathbf{L}}\triangleq(\mathbf{L}_{C,\mathcal{I}})^{\scriptscriptstyle{\top}}\mathbf{L}_{C}^{+}\mathbf{L}_{C,\mathcal{I}}$.

While our main algorithm can sample only from the random-size DPP, and not from the fixed-size $k$-DPP, we present a rigorous reduction argument which lets us use our DPP algorithm to sample exactly from a $k$-DPP (for any $k$) with a small computational overhead.

Prior to our work, fast exact sampling from generic DPPs has been considered out of reach. The first procedure to sample general DPPs was given by and even most recent exact refinements , when the DPP is represented in the form of an $L$-ensemble, require preprocessing that amounts to an expensive $n\times n$ matrix diagonalization at the cost $\mathcal{O}(n^{3})$, which is shown as the first-sample complexity column in Table 1.

Nonetheless, there are well-known samplers for very specific DPPs that are both fast and exact, for instance for sampling uniform spanning trees , which leaves the possibility of a more generic fast sampler open. Since the sampling from DPPs has several practical large scale machine learning applications , there are now a number of methods known to be able to sample from a DPP approximately, outlined in the following paragraphs.

As DPPs can be specified by kernels ($L$-kernels or $K$-kernels), a natural approximation strategy is to resort to low-rank approximations . For example, provides approximate guarantee for the probability of any subset being sampled as a function of eigengaps of the $L$-kernel. Next, construct coresets approximating a given $k$-DPP and then use them for sampling. In their Section 4.1, show in which cases we can hope for a good approximation. These guarantees become tight if these approximations (Nyström subspace, coresets) are aligned with data. In our work, we aim for an adaptive approach that is able to provide a good approximation for any DPP.

The second class of approaches are based on Markov chain Monte-Carlo techniques . There are known polynomial bounds on the mixing rates of MCMC chains with arbitrary DPPs as their limiting measure. In particular, showed them for cardinality-constrained DPPs and for the general case. The two chains have mixing times which are, respectively, linear and quadratic in $n$ (see Table 1). Unfortunately, for any subsequent sample we need to wait until the chain mixes again.

Neither the known low-rank approximations or the known MCMC methods are able to provide samples that are exactly distributed (also called perfect sampling) according to a DPP. This is not surprising as having scalable and exact sampling is very challenging in general. For example, methods based on rejection sampling are always exact, but they typically do not scale with the dimension and are adversely affected by the spikes in the distribution , resulting in high rejection rate and inefficiency. Surprisingly, our method is based on both low-rank approximation (a source of inaccuracy) and rejection sampling (a common source of inefficiency). In the following section, we show how to obtain a perfect DPP sampler from a Nyström approximation of the $L$-kernel. Then, to guarantee efficiency, in Section 3 we bound the number of rejections, which is possible thanks to the use of intermediate downsampling.

Exact sampling using any Nyström approximation

As discussed in the introduction, our method relies on an intermediate downsampling distribution to reduce the size of the problem. The exactness of our sampler relies on the careful choice of that intermediate distribution. To that end, we use regularized determinantal processes, introduced by . In the below definition, we adapt them to the kernel setting.

For any $\mathbf{L}$, $p$, $r$ and $\widetilde{\mathbf{L}}$ defined as in Definition 3,

To sample from the R-DPP, DPP-VFX uses rejection sampling, where the proposal distribution is sampling i.i.d. proportionally to the approximate 1-ridge leverage scores $l_{i}\approx\tau_{i}(1)$ (see Definition 1 and the following discussion), computed using any Nyström approximation $\widehat{\mathbf{L}}$ of matrix $\mathbf{L}$. Apart from $\widehat{\mathbf{L}}$, the algorithm also requires an additional parameter $q$, which controls the size of the intermediate sample. Because of rejection sampling and Proposition 1, the correctness of the algorithm does not depend on the choice of $\widehat{\mathbf{L}}$ and $q$, as demonstrated in the following result. The key part of the proof involves showing that the acceptance probability in Line 4 is bounded by 1. Here, we obtain a considerably tighter bound than the one achieved by , which allows us to use a much smaller intermediate sample $\sigma$ (see Section 3) while maintaning the efficiency of rejection sampling.

Proof We start by showing that the Bernoulli probability in Line 4 is bounded by 1. Note that this is important not only to sample correctly, but also when we later establish the efficiency of the algorithm. If we showed a weaker upper bound, say $c>1$, we could always divide the expression by $c$ and retain the correctness, however it would also be $c$ times less likely that $\textit{Acc}=1$.

Since $\widehat{\mathbf{L}}$ is a Nyström approximation for some $C\subseteq\mathcal{I}$, it can be written as

for any $\mathbf{B}$ such that $\mathbf{L}=\mathbf{B}\mathbf{B}^{\scriptscriptstyle{\top}}$, where $\mathbf{P}\triangleq\mathbf{B}_{C,\mathcal{I}}^{\scriptscriptstyle{\top}}\mathbf{L}_{C}^{+}\mathbf{B}_{C,\mathcal{I}}$ is a projection (so that $\mathbf{P}^{2}=\mathbf{P}$). Let $\widetilde{\mathbf{L}}\triangleq\widetilde{\mathbf{B}}\widetilde{\mathbf{B}}^{\scriptscriptstyle{\top}}$, where the $i$th row of $\widetilde{\mathbf{B}}$ is the rescaled $i$th row of $\mathbf{B}$, i.e. $\widetilde{\mathbf{b}}_{i}^{\scriptscriptstyle{\top}}\triangleq\sqrt{\!\frac{s}{ql_{i}}}\,\mathbf{b}_{i}^{\scriptscriptstyle{\top}}$. Then, we have

Thus, we showed that the expression in Line 4 is valid. Let $\widetilde{\sigma}$ denote the random variable distributed as $\sigma$ is after exiting the repeat loop. It follows that

Let $C\subseteq[n]$ be a random set variable with any distribution. Suppose that $S_{1}$ and $S_{2}$ are returned by two executions of DPP-VFX, both using inputs constructed from the same $\mathbf{L}$ and $\widehat{\mathbf{L}}=\mathbf{L}_{\mathcal{I},C}\mathbf{L}_{C}^{+}\mathbf{L}_{C,\mathcal{I}}$. Then $S_{1}$ and $S_{2}$ are (unconditionally) independent.

Conditions for fast sampling

The complexity cost of DPP-VFX can be roughly summarized as follows: we pay a large one-time cost to precompute $\widehat{\mathbf{L}}$ and all its associated quantities, and then we pay a smaller cost in the rejection sampling scheme which must be multiplied by the number of times we repeat the loop until acceptance. We first show that if the sum of the approximate RLS (i.e., $\sum_{i}l_{i}$, denoted by $s$) is sufficiently close to $k$, then we will exit the loop with high probability. We then analyze how accurate the precomputing step needs to be to satisfy this condition.

If the Nyström approximation $\widehat{\mathbf{L}}$ and the intermediate sample size parameter $q$ satisfy

then $\Pr(\textit{Acc}=\text{true})\geq e^{-2}$. Therefore, with probability $1-\delta$ Algorithm 1 exits the rejection sampling loop after at most $O(\log\delta^{-1})$ iterations and, after precomputing all of the inputs, the time complexity of the rejection sampling loop is O\big{(}k^{6}\log\delta^{-1}+\log^{4}\!\delta^{-1}\big{)}.

Proof Le $\sigma$ be distributed as in Line 3. The probability of exiting the repeat loop at each iteration is

Let $\widehat{\mathbf{L}}$ be constructed by sampling $m=\mathcal{O}(k^{3}\log\frac{n}{\delta})$ columns proportionally to their RLS. Then, with probability $1-\delta$, $\widehat{\mathbf{L}}$ satisfies the assumption of Theorem 3.

There exist many algorithms to sample columns proportionally to their RLS. For example, we can take the BLESS algorithm from with the following guarantee.

There exists an algorithm that with probability $1-\delta$ samples $m$ columns proportionally to their RLS in $\mathcal{O}(nk^{2}\log^{2}\!\frac{n}{\delta}+k^{3}\log^{4}\!\frac{n}{\delta}+m)$ time.

We can now compute the remaining preprocessing costs, given a Nyström approximation $\widehat{\mathbf{L}}$.

Given $\widehat{\mathbf{L}}$ with rank $m$, we can compute $l_{i}$, $s$, $z$, and $\widetilde{\mathbf{L}}$ in $\mathcal{O}(nm^{2}+m^{3})$ time.

We are finally ready to combine these results to fully characterize the computational cost.

subset $S_{1}$ in: $\mathcal{O}(nk^{6}\log^{2}\!\frac{n}{\delta}+k^{9}\log^{3}\!\frac{n}{\delta}+k^{3}\log^{4}\!\frac{n}{\delta})$ time,

then, $S_{2}$ in: \mathcal{O}\big{(}k^{6}\log\frac{1}{\delta}+\log^{4}\!\frac{1}{\delta}\big{)} time.

Discussion. Due to the nature of rejection sampling, as long as we exit the loop, i.e., we accept the sample, the output of DPP-VFX is guaranteed to follow the DPP distribution for any value of $m$. In Theorem 1 we set $m=(k^{3}\log\frac{n}{\delta})$ to satisfy Theorem 3 and guarantee a constant acceptance probability in the rejection sampling loop, but this might not be necessary or even desirable in practice. Experimentally, much smaller values of $m$, starting from $m=\Omega(k\log\frac{n}{\delta})$ seem to be sufficient to accept the sample, while at the same time a smaller $m$ greatly reduces the preprocessing costs. In general, we recommend to separate DPP-VFX in three phases. First, compute an accurate estimate of the RLS using off-the-shelf algorithms in $\mathcal{O}(nk^{2}\log^{2}\!\frac{n}{\delta}+k^{3}\log^{4}\!\frac{n}{\delta})$ time. Then, sample a small number $m$ of columns to construct an explorative $\widehat{\mathbf{L}}$, and try to run DPP-VFX If the rejection sampling loop does not terminate sufficiently fast, then we can reuse the RLS estimates to compute a more accurate $\widehat{\mathbf{L}}$ for a larger $m$. Using a simple doubling schedule for $m$, this procedure will quickly reach a regime where DPP-VFX is guaranteed to accept w.h.p., maintaining its asymptotic complexity, while at the same time resulting in faster sampling in practice.

Reduction from DPPs to k-DPPs

We next show that with a simple extra rejection sampling step we can efficiently transform any exact DPP sampler into an exact $k$-DPP sampler.

A common heuristic to sample $S$ from a $k$-DPP is to first sample $S$ from a DPP, and then reject the sample if the size of $S$ is not exactly $k$. The success probability of this procedure can be improved by appropriately rescaling $\mathbf{L}$ by a constant factor $\alpha$,

Experiments

In this section, we experimentally evaluate the performance of DPP-VFX compared to exact sampling and MCMC-based approaches . In particular, since Section 2 proves that DPP-VFX samples exactly from the DPP, we are interested in evaluating computational performance. This will be characterized by showing how DPP-VFX and baselines scale with the size $n$ of the matrix $\mathbf{L}$ when taking a first sample, and how DPP-VFX achieves constant time when resampling.

To construct $\mathbf{L},$ we use random subsets of the infinite MNIST digits dataset , where $n$ varies up to $10^{6}$ and $d=784$. We use an RBF kernel with $\sigma=\sqrt{3d}$ to construct $\mathbf{L}$. All algorithms are implemented in python. For exact and MCMC sampling we used the DPPy library, , while for DPP-VFX we reimplemented BLESS , and used DPPy to perform exact sampling on the intermediate subset. All experiments are carried out on a 24-core CPU and fully take advantage of potential parallelization. For the Nyström approximation we set $m=10d_{\textnormal{eff}}(1)\approx 10k$. While this is much lower than the $\mathcal{O}(k^{3})$ value suggested by the theory, as we will see it is already accurate enough to result in drastic runtime improvements over exact and MCMC. For each algorithm we controlFor simplicity we do not perform the full $k$-DPP rejection step, but only adjust the expected size of the set. the size of the output set by rescaling the input matrix $\mathbf{L}$ by a constant, following the strategy of Section 4. In Figure 1 we report our results, means and 95% confidence interval over 10 runs, for subsets of MNIST that go from $n=10^{3}$ to $n=7\cdot 10^{4}$, i.e., the whole original MNIST dataset.

Exact sampling is clearly cubic in $n$, and we cannot push our sampling beyond $n=1.5\cdot 10^{4}$. For MCMC, we enforce mixing by runnning the chain for $nk$ steps, the minimum recommended by . However, for $n=7\cdot 10^{5}$ the MCMC runtime is $358$ seconds and cannot be included in the plot, while DPP-VFX completes in $38$ seconds, an order of magnitude faster. Moreover, DPP-VFX rarely rejects more than 10 times, and the mode of the rejections up to $n=7\cdot 10^{5}$ is $1$, that is we mostly accept at the first iteration. Figure 2 reports the cost of the second sample, i.e., of resampling. For exact sampling, this means that an eigendecomposition of $\mathbf{L}$ is already available, but as the plot shows the resampling process still scales with $n$. On the other hand, DPP-VFX’s complexity (after preprocessing) scales only with $k$ and remains constant regardless of $n$.

Finally, we scaled DPP-VFX to $n=10^{6}$ points, a regime where neither exact nor MCMC approaches are feasible. We report runtime and average rejections, with mean and 95% confidence interval over 5 runs. DPP-VFX draws its first sample in $68.479\pm 2.63$ seconds, with only $9\pm 4.5$ rejections.

MD thanks the NSF for funding via the NSF TRIPODS program.

References

Appendix A Omitted proofs for the main algorithm

In this section we present the proofs omitted from Sections 2 and 3, which regarded the correctness and efficiency of DPP-VFX. We start by showing that multiple samples drawn using the same Nyström approximation are independent.

Let $C\subseteq[n]$ be a random set variable with any distribution. Suppose that $S_{1}$ and $S_{2}$ are returned by two executions of DPP-VFX, both using inputs constructed from the same $\mathbf{L}$ and $\widehat{\mathbf{L}}=\mathbf{L}_{\mathcal{I},C}\mathbf{L}_{C}^{+}\mathbf{L}_{C,\mathcal{I}}$. Then $S_{1}$ and $S_{2}$ are (unconditionally) independent.

Proof Let $A$ and $B$ be two subsets of $[n]$ representing elementary events for $S_{1}$ and $C$, respectively. Theorem 2 implies that

Now, for any $A_{1},A_{2}\subseteq[n]$ representing elementary events for $S_{1}$ and $S_{2}$ we have that

Since $\sum_{B\in[n]}\Pr(C\!=\!B)=1$, we get that $S_{1}$ and $S_{2}$ are independent. We now bound the precompute cost, starting with the construction of the Nyström approximation $\widehat{\mathbf{L}}$.

Let $\widehat{\mathbf{L}}$ be constructed by sampling $m=\mathcal{O}(k^{3}\log\frac{n}{\delta})$ columns proportionally to their RLS. Then, with probability $1-\delta$, $\widehat{\mathbf{L}}$ satisfies the assumption of Theorem 3.

Proof Let $\mathbf{L}=\mathbf{B}\mathbf{B}^{\scriptscriptstyle{\top}}$ and $\widehat{\mathbf{L}}=\mathbf{B}\mathbf{P}\mathbf{B}^{\scriptscriptstyle{\top}}$ (where $\mathbf{P}$ is a projection matrix). Using algebraic manipulation, we can write

The $\mathbf{P}\mathbf{B}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{P}$ matrix in the above expression been recently analyzed by in the context of RLS sampling who gave the following result that we use in the proof.

Let the projection matrix $\mathbf{P}$ be constructed by sampling $\mathcal{O}(k\log(\frac{n}{\delta})/\varepsilon^{2})$ columns proportionally to their RLS. Then,

Tuning $\varepsilon=1/(2k+2)$ we obtain $s\leq k+1/2$ and reordering gives us the desired accuracy result. Similarly, we can invert the bound of Proposition 3 to obtain

Given $\mathbf{L}$ and an arbitrary Nyström approximation $\widehat{\mathbf{L}}$ of rank $m$, computing $l_{i}$, $s$, $z$, and $\widetilde{\mathbf{L}}$ requires $\mathcal{O}(nm^{2}+m^{3})$ time.

Appendix B Omitted proofs for the reduction to k-DPPs

In this section we present the proofs omitted from Section 4. Recall that our approach is based on the following rejection sampling strategy:

First, we show the existence of the factor $\alpha^{\star}$ for which the rejection sampling is efficient.

Our starting point is a standard Chernoff bound for $|S_{\alpha}|$.

The distribution of $|S_{\alpha}|$ is given by $\Pr(|S_{\alpha}|=i)\propto e_{i}(\alpha\mathbf{L})$, where $e_{i}(\cdot)$ is the $i$th elementary symmetric polynomial of the eigenvalues of a matrix. Denoting $\lambda_{1},\dots,\lambda_{n}$ as the eigenvalues of $\mathbf{L}$, we can express the elementary symmetric polynomials as the coefficients of the following univariate polynomial with real non-positive roots,

The non-negative coefficients of such a real-rooted polynomial form a unimodal sequence (Lemma 1.1 in ), i.e., $e_{0}(\alpha\mathbf{L})\leq\dots\leq e_{M_{\alpha}}(\alpha\mathbf{L})\geq\dots\geq e_{n}(\alpha\mathbf{L})$, with the mode (shared between no more than two positions $k,k+1$) being close to the mean $k_{\alpha}$: $|M_{\alpha}-k_{\alpha}|\leq 1$ (Theorem 2.2 in ). Moreover, it is easy to see that $M_{0}=0$ and $M_{\alpha}=n$ for large enough $\alpha$, so since the sequence is continuous w.r.t. $\alpha$, for every $k\in[n]$ there is an $\alpha^{\star}$ such that $\Pr(|S_{\alpha^{\star}}|=k)=\Pr(|S_{\alpha^{\star}}|=M_{\alpha^{\star}})$ (every $k$ can become one of the modes). In light of (2), this means that

where the last inequality holds because $|k-k_{\alpha^{\star}}|\leq 1$. Finally, we show how to find $\alpha^{\star}$ efficiently.

Given a Poisson binomial r.v. $|S_{\alpha}|$ with mean $k_{\alpha}$, let $k\triangleq\lfloor k_{\alpha}\rfloor$. The mode $M_{\alpha}$ is

Let $\widehat{\mathbf{L}}$ be constructed by sampling $m=\mathcal{O}((k_{\alpha}/\varepsilon^{2})\log(n/\delta))$ columns proportionally to their RLS. Then with probability $1-\delta$

Proof of Lemma 9 We simply apply the same reasoning of Lemma 2 on both sides. Let $(1-\varepsilon)s_{\alpha^{\star}}=k$, with $\varepsilon$ that will be tuned shortly. Then proving the first inequality to satisfy Proposition 5 is straightforward: $k=(1-\varepsilon)s_{\alpha^{\star}}\leq k_{\alpha^{\star}}$. To satisfy the other side we upper bound $k_{\alpha^{\star}}\leq(1+\varepsilon)s_{\alpha^{\star}}=(1-\varepsilon)s_{\alpha^{\star}}+2\varepsilon s_{\alpha^{\star}}.$ We must now choose $\varepsilon$ such that $2\varepsilon s_{\alpha^{\star}}=1/(k+3)<1/(k+2)$. Substituting, we obtain $\varepsilon=\frac{1}{2(k+3)s_{\alpha^{\star}}}\cdot$ Plugging this in the definition of $s_{\alpha^{\star}}$ we obtain that $\alpha^{\star}$ must be optimized to satisfy

which we plug in the definition of $\varepsilon$ obtaining our neccessary accuracy $\varepsilon=1/(2k^{2}+6k+1)$. Therefore, sampling $m=\widetilde{\mathcal{O}}(k_{\alpha^{\star}}k^{4})$ columns gives us a $s_{\alpha}$ sufficiently accurate to be optimized. However, we still need to bound $k_{\alpha^{\star}}$, which we can do as follows using Lemma 9 and $k\geq 1$

Therefore $m=\widetilde{\mathcal{O}}(k_{\alpha^{\star}}k^{4})\leq\widetilde{\mathcal{O}}(k^{5})$ suffices accuracy wise. Moreover, since $s_{\alpha}$ is parametrized only in terms of the eigenvalues of $\widehat{\mathbf{L}}$, which can be found in $\widetilde{\mathcal{O}}(nm^{2}+m^{3})$ time, we can compute an $\alpha^{\star}$ such that $s_{\alpha^{\star}}=\tfrac{2k^{2}+6k+1}{2k+6}$ in $\widetilde{\mathcal{O}}(nk^{10}+k^{15})$ time, which guarantees $k\leq k_{\alpha^{\star}}<k+\tfrac{1}{k+2}\cdot$ Finally, note that these bounds on the accuracy of $\widehat{\mathbf{L}}$ are extremely conservative. In practice, it is much faster to try to optimize $\alpha^{\star}$ on a much coarser $\widehat{\mathbf{L}}$ first, e.g., for $m=\mathcal{O}(k_{1})$, and only if this approach fails to increase the accuracy of $\widehat{\mathbf{L}}$.