Efficient Sampling for k-Determinantal Point Processes

Introduction

Subset selection problems lie at the heart of many applications where a small subset of items must be selected to represent a larger population. Typically, the selected subsets are expected to fulfill various criteria such as sparsity, grouping, or diversity. Our focus is on diversity, a criterion that plays a key role in a variety of applications, such as gene network subsampling , document summarization , video summarization , content driven search , recommender systems , sensor placement , among many others .

Diverse subset selection amounts to sampling from the set of all subsets of a ground set according to a measure that places more mass on subsets with qualitatively different items. An elegant realization of this idea is given by Determinantal Point Processes (Dpps), which are probabilistic models that capture diversity by assigning subset probabilities proportional to (sub)determinants of a kernel matrix.

Dpps enjoy rising interest in machine learning ; a part of their appeal can be attributed to computational tractability of basic tasks such as computing partition functions, sampling, and extracting marginals . But despite being polynomial-time, these tasks remain infeasible for large data sets. Dpp sampling, for example, relies on an eigendecomposition of the Dpp kernel, whose cubic complexity is a huge impediment. Cubic preprocessing costs also impede wider use of the cardinality constrained variant $k$-Dpp .

These drawbacks have triggered work on approximate sampling methods. Much work has been devoted to approximately sample from a Dpp by first approximating its kernel via algorithms such as the Nyström method , Random Kitchen Sinks , or matrix ridge approximations , and then sampling based on this approximation. However, these methods are somewhat inappropriate for sampling because they aim to project the Dpp kernel onto a lower dimensional space while minimizing a matrix norm, rather than minimizing an error measure sensitive to determinants. Alternative methods use a dual formulation , which however presupposes a decomposition $L=XX^{\top}$ of the DPP kernel, which may be unavailable and inefficient to compute in practice. Finally, MCMC offers a potentially attractive avenue different from the above approaches that all rely on the same spectral technique.

We pursue a yet different approach. While being similar to matrix approximation methods in exploiting redundancy in the data, in sharp contrast to methods that minimize matrix norms, we focus on minimizing the total variation distance between the original Dpp and our approximation. As a result, our approximation models the true Dpp probability distribution more faithfully, while permitting faster sampling. We make the following key contributions:

An algorithm that constructs coresets for approximating a $k$-Dpp by exploiting latent structure in the data. The construction, aimed at minimizing the total variation distance, takes $O(NM^{3})$ time; linear in the number $N$ of data points. The construction works as the overhead of sampling algorithm and is much faster than standard cubic-time overhead that exploits eigendecomposition of kernel matrices. We also investigate conditions under which such an approximation is good.

A sampling procedure that yields approximate $k$-Dpp subsets using the constructed coresets. While most other sampling methods sample diverse subsets in $O(k^{2}N)$ time, the sampling time for our coreset-based algorithm is $O(k^{2}M)$, where $M\ll N$ is a user-specified parameter independent of $N$.

Our experiments indicate that our construction works well for a wide range of datasets, delivers more accurate approximations than the state-of-the-art, and is more efficient, especially when multiple samples are required.

Our sampling procedure runs in two stages. Its first stage constructs an approximate probability distribution close in total variation distance to the true $k$-Dpp. The next stage efficiently samples from this approximate distribution.

Our approximation is motivated by the diversity sampling nature of Dpps: in a Dpp most of the probability mass will be assigned to diverse subsets. This leaves room for exploiting redundancy. In particular, if the data possesses latent grouping structure, certain subsets will be much more likely to be sampled than others. For instance, if the data are tightly clustered, then any sample that draws two points from the same cluster will be very unlikely.

The key idea is to reduce the effective size of the ground set. We do this via the idea of coresets , small subsets of the data that capture function values of interest almost as well as the full dataset. Here, the function of interest is a $k$-Dpp distribution. Once a coreset is constructed, we can sample a subset of core points, and then, based on this subset, sample a subset of the ground set. For a coreset of size $M$, our sampling time is $O(k^{2}M)$, which is independent of $N$ since we are using $k$-Dpps .

Related work.

Dpps have been studied in statistical physics and probability ; they have witnessed rising interest in machine learning . Cardinality-conditioned Dpp sampling is also referred to as “volume sampling”, which has been used for matrix approximations . Several works address faster Dpp sampling via matrix approximations or MCMC . Except for MCMC, even if we exclude preprocessing, known sampling methods still require $O(k^{2}N)$ time for a single sample; we reduce this to $O(k^{2}M)$. Finally, different lines of work address learning DPPs and MAP estimation .

Coresets have been applied to large-scale clustering , PCA and CCA , and segmentation of streaming data .

Setup and basic definitions

where $e_{k}(L)$ is the $k$-th coefficient of the characteristic polynomial $\det(\lambda I-L)=\sum_{k=0}^{N}(-1)^{k}e_{k}(L)\lambda^{N-k}$. We assume that $P_{L,k}(Y)>0$ for all subsets $Y\subseteq\mathcal{Y}$ of cardinality $k$. To simplify notation, we also write $P_{k}\triangleq P_{L,k}$.

Our goal is to construct an approximation $\widehat{P}_{k}$ to $P_{k}$ that is close in total variation distance

and permits faster sampling than $P_{k}$. Broadly, we proceed as follows. First, we define a partition $\Pi=\{\mathcal{Y}_{1},\ldots,\mathcal{Y}_{M}\}$ of $\mathcal{Y}$ and extract a subset $\mathcal{C}\subset\mathcal{Y}$ of $M$ core points, containing one point from each part. Then, for the set $\mathcal{C}$ we construct a special kernel $\widetilde{L}$ (as described in Section 3). When sampling, we first sample a set $\widetilde{Y}\sim\textsc{Dpp}_{k}(\widetilde{L})$ and then, for each $c\in\widetilde{Y}$ we uniformly sample one of its assigned points $y\in\mathcal{Y}_{c}$. These second-stage points $y$ form our final sample. We denote the resulting distribution by $\widehat{P}_{k}=P_{\mathcal{C},k}$. Algorithm 1 formalizes the sampling procedure, which, after one eigendecomposition of the small matrix $\widetilde{L}$.

We begin by analyzing the effect of the partition on the approximation error, and then devise an algorithm to approximately minimize the error. We empirically evaluate our approach in Section 6.

Coreset sampling

Let $\Pi=\{\mathcal{Y}_{1},\ldots,\mathcal{Y}_{M}\}$ be a partition of $\mathcal{Y}$, i.e., $\cup_{i=1}^{M}\mathcal{Y}_{i}=\mathcal{Y}$ and $\mathcal{Y}_{i}\cap\mathcal{Y}_{j}=\emptyset$ for $i\neq j$. We call $\mathcal{C}\subseteq\mathcal{Y}$ a coreset with respect to a partition $\Pi$ if $|\mathcal{C}\cap\mathcal{Y}_{i}|=1$ for $i\in[M]$. With a slight abuse of notation, we index each part $\mathcal{Y}_{c}\in\Pi$ by its core $c\in\mathcal{C}\cap\mathcal{Y}_{c}$. Based on the partition $\Pi$, we call a set $Y\subseteq\mathcal{Y}$ singularIn combinatorial language, $Y$ is an independent set in the partition matroid defined by $\Pi$. with respect to $\Pi^{\prime}\subseteq\Pi$, if for $\mathcal{Y}_{i}\in\Pi^{\prime}$ we have $|Y\cap\mathcal{Y}_{i}|\leq 1$ and for $\mathcal{Y}_{j}\in\Pi\backslash\Pi^{\prime}$ we have $|Y\cap\mathcal{Y}_{j}|=0$. We say $Y$ is $k$-singular if $Y$ is singular and $|Y|=k$.

The following lemma shows that CoreDpp is equivalent to sampling from a $k$-Dpp where we replace each point in $\mathcal{Y}$ by its corresponding core point, and sample with the resulting induced kernel $L_{\mathcal{C}(\mathcal{Y})}$.

CoreDpp is equivalent to sampling from $\textsc{Dpp}_{k}(L_{\mathcal{C}(\mathcal{Y})})$, where in $L_{\mathcal{C}(\mathcal{Y})}$ we replace each element in $\mathcal{Y}_{c}$ by $c$, for all $c\in\mathcal{C}$.

We denote the distribution induced by Algo. 1 by $P_{\mathcal{C},k}$ and that induced by $\textsc{Dpp}_{k}(L_{\mathcal{C}(\mathcal{Y})})$ by $P^{\prime}_{k}$.

First we claim that both sampling algorithms can only sample $k$-singular subsets. By construction, $P_{\mathcal{C},k}$ picks one or zero elements from each $\mathcal{Y}_{c}$. For $P^{\prime}_{k}$, if $Y$ is $k$-nonsingular, then there would be identical rows in $(L_{\mathcal{C}(\mathcal{Y})})_{Y}=L_{\mathcal{C}(Y)}$, resulting in $\det(L_{\mathcal{C}(Y)})=0$. Hence both $P_{\mathcal{C},k}$ and $P^{\prime}_{k}$ only assign nonzero probability to $k$-singular sets $Y$. As a result, we have

For any $Y=\{y_{1},\ldots,y_{k}\}\subseteq\mathcal{Y}$ that is $k$-singular, we have

which shows that these two distributions are identical, i.e., sampling from $\textsc{Dpp}_{k}(\widetilde{L})$ followed by uniform sampling is equivalent to directly sampling from $\textsc{Dpp}_{k}(L_{\mathcal{C}(\mathcal{Y})})$. ∎

Partition, distortion and approximation error

Let us provide some insight on quantities that affect the distance $\|P_{\mathcal{C},k}-P_{k}\|_{\text{\rm tv}}$ when sampling with Algo. 1. In a nutshell, this distance depends on three key quantities (defined below): the probability of nonsingularity $\delta_{\Pi}$, the distortion factor $1+\varepsilon_{\Pi}$, and the normalization factor.

For a partition $\Pi$ we define the nonsingularity probability $\delta_{\Pi}$ as the probability that a draw $Y\sim\textsc{Dpp}_{k}(L)$ is not singular with respect to any $\Pi^{\prime}\subseteq\Pi$.

Given a coreset $\mathcal{C}$, we define the distortion factor $1+\varepsilon_{\Pi}$ (for $\varepsilon_{\Pi}\geq 0$) as a partition-dependent quantity, so that for any $c\in\mathcal{C}$, for all $u,v\in\mathcal{Y}_{c}$, and for any $(k-1)$-singular set $S$ with respect to $\Pi\setminus\mathcal{Y}_{c}$ the following bound holds:

If $\phi$ is the feature map corresponding to the kernel $L$, then geometrically, the numerator of (4.1) is the length of the projection of $\phi(u)$ onto the orthogonal complement of $\operatorname{span}\{\phi(s)\mid s\in S\}$.

The normalization factor for a $k$-Dpp ($L$) is simply $e_{k}(L)$.

Given $\Pi$, $\mathcal{C}$ and the corresponding nonsingularity probability and distortion factors, we have the following bound:

Let $Y\sim\textsc{Dpp}_{k}(L)$ and $\mathcal{C}(Y)$ be the set where we replace each $y\in Y$ by its core $c\in\mathcal{C}$, i.e., $y\in\mathcal{Y}_{c}$. With probability $1-\delta_{\Pi}$, it holds that

Let $c\in\mathcal{C}$ and consider any $(k-1)$-singular set $S$ with respect to $\Pi\setminus\mathcal{Y}_{c}$. Then, for any $v\in\mathcal{Y}_{c}$, using Schur complements and by the definition of $\varepsilon_{\Pi}$ we see that

Here, $Q_{S^{\perp}}$ is the projection onto the orthogonal complement of $\operatorname{span}\{\phi(s)\mid s\in S\}$, and $\phi$ the feature map corresponding to the kernel $L$.

With a minor abuse of notation, we denote by $\mathcal{C}(y)=c$ the core point corresponding to $y$, i.e., $y\in\mathcal{Y}_{c}$. For any $Y=\{y_{1},\ldots,y_{k}\}$, we then define the sets $Y_{i}=\{\mathcal{C}(y_{1}),\ldots,\mathcal{C}(y_{i}),y_{i+1},\ldots,y_{k}\}$, where we gradually replace each point by its core point, with $Y_{0}=Y$. If $Y$ is $k$-singular, then $\mathcal{C}(y_{i})\neq\mathcal{C}(y_{j})$ whenever $i\neq j$, and, for any $0\leq i\leq k-1$, it holds that

This bound holds when $Y$ is $k$-singular, and, by definition of $\delta_{\Pi}$, this happens with probability $1-\delta_{\Pi}$. ∎

Assuming $\varepsilon_{\Pi}$ is small, Lemma 2 states that if replacing a single element in a given subset with another one in the same part does not cause much distortion, then replacing all elements in the subset with their corresponding cores will cause little distortion. This observation is key to our approximation: if we can construct such a partition and coreset, we can safely replace all elements with core points and then approximately sample with little distortion. More precisely, we then obtain the following result that bounds the variational error. Our subsequent construction aims to minimize this bound.

Let $P_{k}=\textsc{Dpp}_{k}(L)$ and let $P_{\mathcal{C},k}$ be the distribution induced by Algo. 1. With the normalization factors $Z=e_{k}(L)$ and $Z_{\mathcal{C}}=e_{k}(\widetilde{L})$, the total variation distance between $P_{k}$ and $P_{\mathcal{C},k}$ is bounded by

From the definition of $Z$ and $Z_{\mathcal{C}}$ we know that $Z=\sum_{|Y|=k}\det(L_{Y})$ and

The last equality follows since, as argued above, $\det(L_{\mathcal{C}(Y)})=0$ for nonsingular $Y$. It follows that

where the first inequality uses the triangle inequality and the second inequality relies on Lemma 2. For the second term in (4.3), we use that, by definition of $\delta_{\Pi}$,

Thus the total variation difference is bounded as

In essence, if the probability of nonsingularity and the distortion factor are low, then it is possible to obtain a good coreset approximation. This holds, for example, if the data has intrinsic (grouping) structure. In the next subsection we provide further intuition on when we can achieve low error.

Theorem 3 depends on the data and the partition $\Pi$. Here, we aim to obtain some further intuition on the properties of $\Pi$ that govern the bound. At the same time, these properties suggest sufficient conditions for a “good” coreset $\mathcal{C}$. For each $\mathcal{Y}_{c}$, we define the diameter

Next, define the minimum distance of any point $u\in\mathcal{Y}_{c}$ to the subspace spanned by the feature vectors of points in a “complementary” set $S$ that is singular with respect to $\Pi\setminus\mathcal{Y}_{c}$:

Lemma 4 connects these quantities with $\varepsilon_{\Pi}$; it essentially poses a separability condition on $\Pi$ (i.e., $\Pi$ needs to be “aligned” with the data) so that the bound on $\varepsilon_{\Pi}$ holds.

If $d_{c}>\rho_{c}$ for all $c\in\mathcal{C}$, then

For any $c\in\mathcal{C}$ and any $u,v\in\mathcal{Y}_{c}$ and $S$ $(k-1)$-singular with respect to $\Pi\backslash\mathcal{Y}_{c}$, we have

Without loss of generality, we assume $\det(L_{S\cup\{u\}})\geq\det(L_{S\cup\{v\}})$. By definition of $\rho_{c}$ we know that

Since $0<\|Q_{S^{\perp}}\phi(v)\|\leq\|Q_{S^{\perp}}\phi(u)\|\leq\|\phi(u)\|$ by assumption, we have

Then, by definition of $\varepsilon_{\Pi}$, we have

Efficient construction

Thm. 3 states an upper bound on the error induced by CoreDpp and relates the total variation distance to $\Pi$ and $\mathcal{C}$. Next, we explore how to efficiently construct $\Pi$ and $\mathcal{C}$ that approximately minimize the upper bound.

Any set $Y$ sampled via CoreDpp is, by construction, singular with respect to $\Pi$. In other words, CoreDpp assigns zero mass to any nonsingular set. Hence, we wish to construct a partition $\Pi$ such that its nonsingular sets have low probability under $\textsc{Dpp}_{k}(L)$. The optimal such partition minimizes the probability $\delta_{\Pi}$ of nonsingularity. A small $\delta_{\Pi}$ value also means that the parts of $\Pi$ are dense and compact, i.e., the diameter $\rho_{c}$ in Equation (4.4) is small.

Finding such a partition optimally is hard, so we resort to local search. Starting with a current partition $\Pi$, we re-assign each $y$ to a part $\mathcal{Y}_{c}$ to minimize $\delta_{\Pi}$. If we assign $y$ to $\mathcal{Y}_{c}$, then the probability of sampling a set $Y$ that is singular with respect to the new partition $\Pi$ is

where s_{k}^{\Pi}(L):=\sum_{Y\text{k-sing.}}\det(L_{Y}). The matrix $L_{\backprime c}$ denotes $L$ with rows $\mathcal{Y}_{c}$ and columns $\mathcal{Y}_{c}$ deleted, and $L^{y}=L-L_{\mathcal{Y},y}L_{y,\mathcal{Y}}$. For local search, we would hence compute $L_{yy}s_{k-1}^{\Pi}(L_{{\backprime}c}^{y})$ for each point $y$ and core $c$, assign $y$ to the highest-scoring $c$. Since this testing is still expensive, we introduce further speedups in Section 5.3.

2 Constructing 𝒞𝒞\mathcal{C}

When constructing $\mathcal{C}$, we aim to minimize the upper bound on the total variation distance between $P_{k}$ and $P_{\mathcal{C},k}$ stated in Theorem 3. Since $\delta_{\Pi}$ and $\varepsilon_{\Pi}$ only depend on $\Pi$ and not on $\mathcal{C}$, we here focus on minimizing $|1-\frac{Z_{\mathcal{C}}}{Z}|$, i.e., bringing $Z_{\mathcal{C}}$ as close to $Z$ as possible. To do so, we again employ local search and subsequently swap each $c\in\mathcal{C}$ with its best replacement $v\in\mathcal{Y}_{c}$. Let $\mathcal{C}^{c,v}$ be $\mathcal{C}$ with $c$ replaced by $v$. We aim to find the best swap

Computing $Z$ requires computing the coefficients $e_{k}(L)$, which takes a total of $O(N^{3})$ timeIn theory, this can be computed in $O(N^{\omega}\log(N))$ time , but the eigendecompositions and dynamic programming used in practice typically take cubic time.. In the next section, we therefore consider a fast approximation.

3 Faster constructions and further speedups

Local search procedures for optimizing $\Pi$ and $\mathcal{C}$ can be further accelerated by a sequence of relaxations that we found to work well in practice (see Section 6). We begin with the quantity $s_{k-1}^{\Pi}(L_{{\backprime}c}^{y})$ that involves summing over sub-determinants of the large matrix $L$. Assuming the initialization is not too bad, we can use the current $\mathcal{C}$ to approximate $\mathcal{Y}$. In particular, when re-assigning $y$, we substitute all other elements with their corresponding cores, resulting in the kernel $\widehat{L}=L_{\mathcal{C}(\mathcal{Y})}$. This changes our objective to finding the $c\in\mathcal{C}$ that maximizes $s_{k-1}^{\Pi}(\widehat{L}_{\backprime c}^{y})$. Key to a fast approximation is now Lemma 5, which follows from Lemma 1.

the last equality was shown in the proof of Thm. 3. ∎

Second, when constructing $\mathcal{C}$, we observed that $Z_{\mathcal{C}}$ is commonly much smaller than $Z$. Hence, a fast approximation merely greedily increases $Z_{\mathcal{C}}$ without computing $Z$.

Third, we can be lazy in a number of updates: for example, we only consider changing cores for the part that changes. When a part $\mathcal{Y}_{c}$ receives a new member, we check whether to switch the current core $c$ to the new member. This reduction keeps the core adjustment at time $O(M^{3})$. Moreover, when re-assigning an element $y$ to a different part $\mathcal{Y}_{c}$, it is usually sufficient to only check a few, say, $\nu$ parts with cores closest to $y$, and not all parts. The resulting time complexity for each element is $O(M^{3})$.

With this collection of speedups, the approximate construction of $\Pi$ and $\mathcal{C}$ takes $O(NM^{3})$ for each iteration, which is linear in $N$, and hence a huge speedup over direct methods that require $O(N^{3})$ preprocessing. The iterative algorithm is shown in Algorithm 2. The initialization also affects the algorithm performance, and in practice we find that kmeans++ as an initialization works well. Thus we use CoreDpp to refer to the algorithm that is initialized with kmeans++ and uses all the above accelerations. In practice, the algorithm converges very quickly, and most of the progress occurs in the first pass through the data. Hence, if desired, one can even use early stopping.

Experiments

We next evaluate CoreDpp, and compare its efficiency and effectiveness against three competing approaches:

Partitioning using $k$-means (with kmeans++ initialization ), with $\mathcal{C}$ chosen as the centers of the clusters; referred to as K++ in the results.

The adaptive, stochastic Nyström sampler of (NysStoch). We used $M$ dimensions for NysStoch, to use the same dimensionality as CoreDpp.

The Metropolis-Hastings DPP sampler MCDPP . We use the well-known Gelman and Rubin multiple sequence diagnostic to empirically judge mixing.

In addition, we show results using different variants of CoreDpp: CoreDpp-z described in Sec. 5.3 and variants that are initialized either randomly (CoreDpp-r) or via kmeans++ (CoreDpp).

We first explore the effect of our fast approximate sampling on controllable synthetic data. The experiments here compare the accuracy of the faster CoreDpp from Section 5.3 to CoreDpp-z, CoreDpp-r and K++.

Fig. 1 shows the total variation distance $\|\widehat{P}_{k}-P_{k}\|_{\text{\rm tv}}$ defined in Equation (2.1) for the partition and cores generated by K++, CoreDpp, CoreDpp-r and CoreDpp-z as nClust and the length vary. We see that in general, most approximations improve as $\varepsilon_{\Pi}$ and $\delta_{\Pi}$ shrink. Remarkably, the CoreDpp variants achieve much lower error than K++. Moreover, the results suggest that the relaxations from Section 5.3 do not noticeably increase the error in practice. Also, CoreDpp-r performs comparable with CoreDpp, indicating that our algorithm is robust against initialization. Since, in addition, the CoreDpp construction makes most progress in the first pass through the data, and the kmeans++ initialization yields the best performance, we use only one pass of CoreDpp initialized with kmeans++ in the subsequent experiments.

2 Real Data

We apply CoreDpp to two larger real data sets:

MNIST . MNIST consists of images of hand-written digits, each of dimension $28\times 28$.

GENES . This dataset consists of different genes. Each sample in GENES corresponds to a gene, and the features are shortest path distances to 330 different hubs in the BioGRID gene interaction network.

For our first set of experiments on both datasets, we use a subset of 2000 data points and an RBF kernel to construct $L$. To evaluate the effect of model parameters on performance, we vary $M$ from 20 to 100 and $k$ from 2 to 8 and fix $\nu=2$ (see Section 6.1 for an explanation of the parameters). Larger-scale experiments on these datasets are reported in Section 6.3.

On these larger data sets, it becomes impossible to compute the total variation distance exactly. We therefore approximate it by uniform sampling and computing an empirical estimate.

The results in Figure 2 and Figure 3 indicate that the approximations improve as the number of parts $M$ increases and $k$ decreases. This is because increasing $M$ increases the models’ approximation power, and decreasing $k$ leads to a simpler target probability distribution to approximate. In general, CoreDpp always achieves lower error than K++, and NysStoch performs poorly in terms of total variation distance to the original distribution. This phenomenon is perhaps not so surprising when recalling that the Nyström approximation minimizes a different type of error, a distance between the kernel matrices. These observations suggest to be careful when using matrix approximations to approximate $L$.

For an intuitive illustration, Figure 4 shows a core $\mathcal{C}$ constructed by CoreDpp, and the elements of one part $\mathcal{Y}_{c}$.

3 Running Time on Large Datasets

Lastly, we address running times for CoreDpp, NysStoch and the Markov chain $k$-DPP (MCDPP ). For the latter, we evaluate convergence via the Gelman and Rubin multiple sequence diagnostic ; we run 10 chains simultaneously and use the CODA package to calculate the potential scale reduction factor (PSRF), and set the number of iterations to the point when PSRF drops below 1.1. Finally we run MCDPP again for this specific number of iterations.

For overhead time, i.e., time to set up the sampler that is spent once in the beginning, we compare against NysStoch: CoreDpp constructs the partition and $\widetilde{L}$, while NysStoch selects landmarks and constructs an approximation to the data. For sampling time, we compare against both NysStoch and MCDPP: CoreDpp uses Algo. 1, and NysStoch uses the dual form of $k$-Dpp sampling . We did not include the time for convergence diagnostics into the running time of MCDPP, giving it an advantage in terms of running time.

Fig. 5 shows the overhead times as a function of $N$. For MNIST we vary $N$ from 6,000 to 20,000 and for GENES we vary $N$ from 6,000 to 10,000. These values of $N$ are already quite large, given that the Dpp kernel is a dense RBF kernel matrix; this leads to increased running time for all compared methods. The construction time for NysStoch and CoreDpp is comparable for small-sized data, but NysStoch quickly becomes less competitive as the data gets larger. The construction time for CoreDpp is linear in $N$, with a mild slope. If multiple samples are sought, this construction can be performed offline as preprocessing as it is needed only once.

Sampling.

Fig. 6 shows the time to draw one sample as a function of $N$, comparing CoreDpp against NysStoch and MCDPP. CoreDpp yields samples in time independent of $N$ and is extremely efficient – it is orders of magnitude faster than NysStoch and MCDPP.

We also consider the time taken to sample a large number of subsets, and compare against both NysStoch and MCDPP—the sampling times for drawing approximately independent samples with MCDPP add up. Fig. 7 shows the results. As more samples are required, CoreDpp becomes increasingly efficient relative to the other methods.

Conclusion

In this paper, we proposed a fast, two-stage sampling method for sampling diverse subsets with $k$-Dpps. As opposed to other approaches, our algorithm directly aims at minimizing the total variation distance between the approximate and original probability distributions. Our experiments demonstrate the effectiveness and efficiency of our approach: not only does our construction have lower error in total variation distance compared with other methods, it also produces these more accurate samples efficiently, at comparable or faster speed than other methods.

This research was partially supposed by an NSF CAREER award 1553284 and a Google Research Award.

References