Diversity Subsampling: Custom Subsamples from Large Data Sets

1. Introduction

Diversity subsampling selects a subset of points from a data set with the goal of having the selected points spread out evenly over the data space. The goal of diverse samples is ubiquitous in the field of experimental design and is a common design criterion. In the specific situation in which some supervised learning model will be fit to the selected subsample (after some follow-up experiment in which response labels are observed for each point in the sample), a uniform or space-filling subsample will typically be beneficial. For instance, training the model with a space-filling subsample, as opposed to a random subsample that follows the same distribution as the training data, provides some robustness to the covariate drift issue (i.e., the distribution of the future test data to which the supervised learning model will be applied is different from the training data distribution). As another example, sometimes the goal may be to use the fitted supervised learning algorithm to optimize the response (regarded as a function of the covariates), a more space-filling subsample seems advantageous because it reveals the behavior of the response in the entire covariate space. In the literature, diversity subsampling has been found useful in a variety of situations, for example, Yu and Kim (2010), Puzyn et al. (2011), Haussmann et al. (2020), Song et al. (2022b), and Silveira and Barbeira (2022), to name a few.

Suppose we have data set D consisting of an identically and independently distributed (i.i.d.) sample of size N of some random vector $\mathit{X}\in {\mathbb{R}}^q$; that is, $D=\{{\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_N\}$, where ${\mathit{x}}_i\in {\mathbb{R}}^q$. Let ${\mathcal{S}}_n=\{{\mathit{x}}_{j_1},\dots ,{\mathit{x}}_{j_n}\}\subset D$, for $\{j_1,\dots ,j_n\}\subset \{1,2,\dots ,N\}$ denote a subsample of size n selected from D. In the literature, loosely speaking, ${\mathcal{S}}_n$ is called a diverse subsample from D if the points it constitutes are as different from each other as possible (sometimes called the repulsiveness property; Wang et al. 2018, Biyik et al. 2019) and cover as much of the effective support of the data as possible. In some machine learning applications, each observation ${\mathit{x}}_i$ is already labeled with a response ${\mathit{y}}_i$; our paper focuses on the setting where the labels ${\{{\mathit{y}}_i\}}_{i=1}^N$ are not available, and one will observe the responses only for the n cases in ${\mathcal{S}}_n$ (e.g., by conducting some follow-up experimental or observational study, surveying the cases, etc.) which will then serve as the training data for fitting some supervised learning model. For situations in which all N response labels are available and the goal is supervised compression of the data (Joseph and Mak 2021, Chen and Zhang 2022).

Existing popular methods for diversity subsampling from given data include determinantal point processes (DPP; Biyik et al. 2019), computer-aided design of experiments (CADEX; Kennard and Stone 1969), Poisson disk sampling (PDS; Cook 1986, Yuksel 2015), and k-means clustering (MacQueen et al. 1967, Wu 2018). DPP can be used for diversity subsampling by first specifying an appropriately constructed kernel matrix so that the most diverse set of points will have the largest probability of being selected under this DPP; this set is thus the mode of the DPP. It has been shown that finding this mode is an NP-hard problem (Ko et al. 1995). Various algorithms have been proposed to approximate the mode of a DPP (Biyik et al. 2019, Han and Gillenwater 2020). CADEX selects points sequentially from data such that the distance between each newly selected point and its closest point in the previously selected point set is as large as possible. Song et al. (2022b) recently proposed an efficient heuristic for the CADEX algorithm, which is called the scSampler algorithm. PDS is a conceptually simple algorithm that selects subsequent points outside of the union of hyper-spherical neighborhoods of all selected points. The radius of these hyper-spherical neighborhoods is usually user specified and taken as an input of the algorithm, which ensures that the pairwise distances between all selected points are no smaller than twice the specified radius (Cook 1986). Existing PDS algorithms that do not require the users to prespecify a fixed radius include McCool and Fiume (1992) and Yuksel (2015). Wu (2018) proposed the RD ALR method (representativeness and diversity, active learning for regression) that uses k-means clustering (MacQueen et al. 1967) to select a diverse subsample from D; it selects each point in ${\mathcal{S}}_n$ sequentially. At iteration j, RD ALR clusters the entire data into j + 1 clusters using k-means, and the newly selected point are chosen from clusters that do not contain any previously selected points.

One common property of the previously mentioned methods is that the subsamples they select are all repulsive, in the sense that they attempt to ensure that the selected points are as far away from each other as possible. Although repulsiveness is an appealing property of a uniform or space-filling subsample, it might not be achievable when there are many replicated observations in D. In contrast to the existing subsampling algorithms, we define our diversity subsampling goal as to select a subsample that follows, as closely as possible, a mutually independent uniform distribution over the effective support of D, and we refer to our proposed algorithm as the diversity subsampling (DS) algorithm. Our DS algorithm has substantial benefits compared with the existing ones. First, aiming to select an i.i.d uniform subsample over the effective support of D allows replications to occur in the subsample. This has advantages in applications in which the response variable Y is observed with error and/or when there are additional latent variables because in these situations, it may be desirable to observe Y at multiple very similar X values. Second, our DS algorithm handles data with large numbers of replicated values along all or certain dimensions much better than existing algorithms: In this challenging situation, the subsample selected by the DS algorithm is still as evenly spread out over the data space as possible, but the subsamples selected by existing algorithms degrade to random sampling (we provide an example to illustrate this shortly). Third, our DS algorithm turns out to be much more computationally efficient, especially with larger N and n (which is increasingly common in the big data era) than existing methods (see Section 5.4 for details). Last, it is easy to generalize our DS algorithm to select subsamples following any desired distribution (not just uniform).

To illustrate the better uniformity properties of the DS subsample than existing methods when there are many replicated observations in D, consider the following simple example in which X follows a bivariate standard normal distribution with independent components (see Section 5.3 for experimental results using quantitative measures with additional examples). We set q = 2, N = 10,000, and n = 2,000 in this example. The 10,000 data points were generated using 2,000 independently drawn points from the bivariate standard normal distribution, each of which was replicated five times. Figure 1 compares heat maps of the estimated density of typical subsamples selected by three different methods: random sampling from a data set, scSampler-sp1 (Song et al. 2022b), and our DS algorithm (see Section 2 for details of the DS algorithm). The densities were estimated using a Gaussian kernel with a bandwidth of 0.1. As Figure 1 shows, the points selected by random sampling (the left plot of Figure 1) follow the original bivariate standard normal distribution; that is, they concentrate in the middle region of the plotted area, where the data points (not shown in Figure 1) are denser. The subsample selected by the scSampler-sp1 algorithm (the middle plot of Figure 1) also exhibits the same phenomenon. In contrast, the subsample selected by DS (the right plot of Figure 1) is much closer to being uniformly distributed over the effective support of the data.

The structure of this paper is as follows. Section 2 presents our DS algorithm. Section 3 discusses generalizations of the proposed DS algorithm to select subsamples following any desired distribution. Section 4 discusses theoretical convergence of the sample produced by the DS algorithm to the desired target distribution. Section 5 numerically demonstrates the effectiveness of the DS algorithm and its advantages compared with known methods using data exhibiting a variety of distributions. Section 6 concludes the paper.

2. DS Algorithm for Diversity Subsampling

The goal of our DS algorithm is to select a subsample ${\mathcal{S}}_n$ of size n (typically, n is much less than the data size N) points in $D=\{{\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_N\}$ that is distributed as closely as possible to an i.i.d. sample from the uniform distribution over S, the effective support of the data distribution. First, we discuss the case when the distribution of X is continuous, and then we extend our DS method for general data sets with discrete or mixed continuous and discrete data distributions. Let $f(\mathit{x})$ and S denote the probability density function (p.d.f) of X and its support, both of which are assumed unknown. We develop two versions of the DS algorithm, for sampling from D with or without replacement. The former is relevant in a smaller class of applications in which it is desirable to observe the response Y multiple times for the same observation X. For brevity, we focus on the sampling-without-replacement version of the DS algorithm (denoted as the DS algorithm) in this section and discuss the sampling-with-replacement version (denoted as the DS-WR algorithm) briefly. We present certain asymptotic performance results for the DS algorithm in a more general way in Section 4, where we demonstrate that under certain assumptions the uniformity and independence of points in ${\mathcal{S}}_n$ can be guaranteed asymptotically.

Our approach adapts the well-known sampling/importance resampling (SIR) algorithm (Rubin 1987, 1988), which subsamples proportional to certain weights from simulated data generated from some proposal distribution (as opposed to subsampling from a given data set). The SIR algorithm works as follows. To generate i.i.d. samples from a target distribution with p.d.f. $g(\mathit{x})$, SIR samples from some proposal p.d.f. $f_{\text{prop}}(\mathit{x})$ that is ideally similar to g and convenient to sample from. SIR first generates a sample from $f_{\text{prop}}$, say $D_{\text{SIR}}=\{{\mathit{x}}_1,\dots ,{\mathit{x}}_M\}$, and then randomly samples (with or without replacement) a smaller set ${\mathcal{S}}_{\text{SIR}}$ from $D_{\text{SIR}}$ with the probability of each point being selected proportional to $g({\mathit{x}}_i)/f_{\text{prop}}({\mathit{x}}_i),\hspace{0.17em}i=1,\dots ,M$. The selected ${\mathcal{S}}_{\text{SIR}}$ is asymptotically an i.i.d. sample from $g(\mathit{x})$ as M approaches infinity (Skare et al. 2003).

The SIR algorithm subsamples from a simulated data set generated from a known proposal distribution and is not directly applicable to subsample from a given data set with unknown distribution. We adapt it to the latter setting, as follows. If we view the (unknown) data distribution $f(\mathit{x})$ as the proposal distribution and the uniform distribution over S as the target distribution, then we can obtain an approximately i.i.d. uniform distribution over S by subsampling from D with the probability of each point being selected proportional to $(1/|S|)/f({\mathit{x}}_i)\hspace{0.17em}\propto 1/f({\mathit{x}}_i)\approx 1/{\widehat{f}}_N({\mathit{x}}_i)$, for $i=1,\dots ,N$, for some suitable estimator ${\widehat{f}}_N$ of the density f. Procedure 1 summarizes the DS algorithm for sampling without replacement. In Section 4, we present theoretical convergence results for the DS algorithm with known $f(\mathit{x})$, and in Section 5, we demonstrate numerically that it produces samples that reasonably approximate the desired target distribution. In the DS algorithm, we first estimate $f(\mathit{x})$ at each data point in D (see Appendix B for details). Denote the estimated density at ${\mathit{x}}_i$ as $\widehat{f}({\mathit{x}}_i)$, for $i=1,\dots ,N$. Then the DS algorithm selects the n points in ${\mathcal{S}}_n$ from D sequentially at each iteration k as follows.

At iteration k = 1, the DS algorithm selects the index $j_1\in \{1,\dots ,N\}$ of the first point in ${\mathcal{S}}_n$, denoted by ${\mathit{x}}_{j_1}$, with probability

At iteration $k\ge 2$, the DS algorithm selects the index $j_k\in \{1,\dots ,N\}\backslash \{j_1,\dots ,j_{k-1}\}$ of the next sampled point ${\mathit{x}}_{j_k}$ with probability

Three practical issues must be addressed when applying the previous idea to real data sets. One issue is that n might not be negligibly small compared with N. Because we are sampling without replacement, the size of the remaining data set $D\backslash {\mathcal{S}}_{k-1}$ at iteration k decreases with k; thus the distribution of the remaining data can change dramatically from the distribution of the original data D. We remedy this by updating the estimated density regularly so that it reflects the distribution of the remaining data and not the distribution of D. In particular, we update the density of remaining points in D after every $M=\max \{100,\lfloor n/\mathcal{U}\rfloor \}$ points have been selected, where U is a user-chosen integer (we use U = 10 in all examples). Here the symbol $\lfloor n/\mathcal{U}\rfloor$ denotes the largest integer not larger than $n/\mathcal{U}$. We provide a method to update density efficiently in Appendix B.

The second issue is that the density estimation method proposed in Appendix B only applies when the data distribution is continuous. In the case that the data distribution is discrete or mixed continuous and discrete, we approximate the true distribution with a continuous one by adding a small Gaussian noise (lines 3–10 of Procedure 1). For simplicity and because data are often rounded, we apply this Gaussian perturbation for all data sets regardless of the true data distribution being discrete or not. To be specific, for the density estimation step, $\forall i\in \{1,\dots ,N\}$, we replace ${\mathit{x}}_i$ with ${\mathit{x}}_i+{\sigma }_p{\mathit{Z}}_i$, where ${\{{\mathit{Z}}_i\}}_{i=1}^N$ follow i.i.d. multivariate normal distribution $\mathcal{N}({\mathbf{0}}_q,{\mathbf{1}}_{q\times q})$, where ${\mathbf{0}}_q$ is the zero vector in ${\mathcal{R}}^q$ and ${\mathbf{1}}_{q\times q}$ is the identity matrix in ${\mathcal{R}}^{q\times q}$. We choose ${\sigma }_p>0$ to be a small number compared with the pairwise distances of D. In Procedure 1, we chose σ_p to be approximately $\frac{1}{8}$ times the smallest nonzero pairwise distance among D. In the case that D consists only of discrete observations, this choice aims to ensure that a perturbed data point stays close to its original location. Precisely, for any point ${\mathit{x}}_i\in D,\hspace{0.17em}i=1,\dots ,N$, let its closest nonoverlapping neighbor in D be ${\mathit{x}}_{i_1}$ and let $\mathit{u}\in {\mathbb{R}}^{q\times 1}$ denote the normalized vector aligned with ${\mathit{x}}_{i_1}-{\mathit{x}}_i$. Regarding ${\mathit{x}}_i$ as the origin, the projected univariate noise (added to ${\mathit{x}}_i$) along the direction of u will be a normal random variable with mean ${\mathit{u}}^T{\mathbf{0}}_{q\times 1}=0$ and variance ${\mathit{u}}^T({\sigma }_p^2{\mathbf{1}}_{q\times q})\mathit{u}={\sigma }_p^2$. By choosing σ_p at least as large as 0.125$\Vert {\mathit{x}}_{i_1}-{\mathit{x}}_i\Vert$, we ensure that the probability of the projected univariate noise in the subspace of u being located closer to ${\mathit{x}}_i$ than ${\mathit{x}}_{i_1}$ is not smaller than $1-6.4\times {10}^{-5}$. In the case that D consists only of continuous observations, the added noise will be negligible. Adding Gaussian noise is actually consistent with the notion of Gaussian kernel density estimation (KDE), which can be viewed as convolving the empirical distribution of D with a Gaussian “noise” density with standard deviation related to the kernel bandwidth.

The third issue is related to the second and regards estimating the density $f({\mathit{x}}_i)$. Let ${\tilde{\mathit{x}}}_i$ denote the perturbed version of ${\mathit{x}}_i\in D,\hspace{0.17em}i=1,\dots ,N$. We can consider estimating $f({\mathit{x}}_i)$ by $\widehat{f}({\mathit{x}}_i)$ or by $\widehat{f}({\tilde{\mathit{x}}}_i)$. When the data distribution is continuous, these two options yield similar results because, by design, the noises added to continuous data are negligible. When the data distribution is discrete, the former choice may seem intuitively preferable because it ensures that whenever ${\mathit{x}}_i={\mathit{x}}_j$, we will have $\widehat{f}({\mathit{x}}_i)=\widehat{f}({\mathit{x}}_j),\hspace{0.17em}\{i,j\}\in \{1,\dots ,N\}$. However, in practice, we observed that our DS algorithm appears to tolerate the standard deviation in the estimated density better than it tolerates bias, and the latter choice provided more diverse subsamples (see Appendix A for an numerical example).

Besides the SIR algorithm (Rubin 1987, 1988), one can also consider adapting the importance support points (ISP) resampling algorithm (algorithm 3 in Huang et al. 2022) into a diversity subsampling algorithm. The ISP resampling algorithm was developed and used for the standard Markov China Monte Carlo setting of approximating a posterior distribution. We choose to adapt SIR in this paper due to computational and storage issues with the ISP resampling algorithm.

Procedure 1

(DS Algorithm: Diversity Subsampling Without Replacement)

Input: n, $D=\{{\mathit{x}}_1,\dots ,{\mathit{x}}_N\}\subset {[0,1]}^q$

• 1: $N\leftarrow |D|$

• 2: Standardize D into ${[0,1]}^q$, and store this new data set as $\tilde{D}={\{{\tilde{\mathit{x}}}_i\}}_{i=1}^N$

• 3: $n_s\leftarrow \min \{2000,\lfloor \frac{N}{4}\rfloor \}$

• 4: Randomly select $\{r_1,\dots ,r_{n_s}\}\subset \{1,\dots ,N\}$ with equal probabilities

• 5: $D_{\text{sub}}\leftarrow$ all unique points in ${\{{\tilde{\mathit{x}}}_{r_j}\}}_{j=1}^{n_s}$

• 6: Repeat lines 4 and 5 until $|D_{\text{sub}}|>1$

• 7: ${\sigma }_p\leftarrow \frac{1}{8}{\min }_{\{\tilde{\mathit{x}},\tilde{\mathit{y}}\}\subset D_{\text{sub}}}{\Vert \tilde{\mathit{x}}-\tilde{\mathit{y}}\Vert }_2$

• 8: for $i=1,\dots ,N$ do

• 9:  Independently generate ${\mathit{Z}}_i\sim \mathcal{N}({\mathbf{0}}_q,{\mathbf{1}}_{q\times q})$

• 10:  ${\tilde{\mathit{x}}}_i\leftarrow {\tilde{\mathit{x}}}_i+{\sigma }_p{\mathit{Z}}_i$

• 11: end for

• 12: Estimate density $\widehat{f}({\tilde{\mathit{x}}}_i)$, for each $i=1,2,\dots ,N$ (see Appendix B)

• 13: $\widehat{f}({\mathit{x}}_i)\leftarrow \widehat{f}({\tilde{\mathit{x}}}_i)$, for each $i=1,2,\dots ,N$

• 14: Set index array $I_1\leftarrow \{1,\dots ,N\}$

• 15: Select j₁ randomly from I₁ with probabilities given by Equation (1)

• 16: ${\mathcal{S}}_1\leftarrow \{{\mathit{x}}_{j_1}\}$

• 17: Set an estimated number of density estimation updates U (e.g., $U\leftarrow 10$)

• 18: Set the number of points to select before conducting a density update. $M\leftarrow \max \{100,\lfloor \frac{n}{U}\rfloor \}$

• 19: Set $\mathit{\text{UpdateState}}\leftarrow \{M,2M,\dots ,\lfloor \frac{n}{M}\rfloor M\}$

• 20: for $k=2,\dots ,n$ do

• 21:  $I_k\leftarrow I_{k-1}\backslash j_{k-1}$

• 22:  if $k-1\in \mathit{\text{UpdateState}}$ then

• 23:   Re-Estimate ${\{\widehat{f}({\tilde{\mathit{x}}}_i)\}}_{i\in I_k}$ by efficient updating (See Appendix B for details)

• 24:   Reset $\widehat{f}({\mathit{x}}_i)\leftarrow \widehat{f}({\tilde{\mathit{x}}}_i)$, for each $i\in I_k$

• 25:  end if

• 26:  Select j_k randomly from I_k with probabilities given by Equation (2)

• 27:  ${\mathcal{S}}_k\leftarrow {\mathcal{S}}_{k-1}\cup \{{\mathit{x}}_{j_k}\}$

• 28: end for

Output: Subsample set ${\mathcal{S}}_n=\{{\mathit{x}}_{j_1},{\mathit{x}}_{j_2},\dots ,{\mathit{x}}_{j_n}\}$

2.1. Sampling-with-Replacement Version of the DS Algorithm

In certain situations, subsampling from D with replacement may be relevant. In our primary setting where D is unlabeled data, and the intent is to observe a response Y for each point in ${\mathcal{S}}_n$, sampling with replacement is only relevant if it is possible to experiment on (i.e., observe a response y_i for) the same subject ${\mathit{x}}_i$ multiple times, each time observing a potentially different stochastic value for y_i. As an example, suppose ${\mathit{x}}_i$ is a set of measured characteristics of sample i of a chemical compound from a large set D of samples over which X varies randomly, y_i is some output property of a subsequent chemical reaction with reactants taken from sample i, and y_i can vary randomly if the reaction is repeated using reactants from sample i again. In this case, one may want to use sampling from D with replacement, so the same sample can be experimented on multiple times if needed. Conversely, suppose ${\mathit{x}}_i$ is a set of clinical, phenotypical, and demographic characteristics of a patient i having a particular condition, and y_i is the efficacy of a treatment for that condition administered to patient i. In this case, it is not meaningful to apply an experimental treatment to the same patient twice, so sampling with replacement becomes irrelevant. See Rubin (1988) for additional discussion of sampling with or without replacement.

For the sampling-with-replacement version of the DS algorithm, denoted by DS-WR, at each iteration $k\in \{1,2,\dots ,n\}$, we select the next point ${\mathit{x}}_{j_k}$ from D as follows. The index $j_k\in \{1,2,\dots ,N\}$ is selected with probability

independently of which points have been selected in earlier iterations. The DS-WR algorithm is identical to Procedure 1, except that for DS-WR, one repeats line 15 n times to produce the indices $\{j_1,j_2,\dots ,j_n\}$ of the desired subsample of size n, and all lines after line 15 are omitted.

If the DS-WR version is relevant for a particular problem, there is a potential drawback that should be considered if there are severe outliers in the data, as also noted in section 10.4 in Gelman et al. (1995). Most existing density estimation techniques tend to underestimate the density in extremely sparse tail regions of S, where there are only a few outlier points in D. This causes the estimated density at the outlier x points to have extremely low values and therefore very high probabilities of getting selected in an iteration of the DS algorithm. To provide the most diversity, it is generally desirable to select these points for the subsample, which the DS algorithm does. However, in the DS-WR algorithm, the same outlier points may get selected repeatedly for the subsample, which ends up being counterproductive for the diversity goal.

Figure 2 shows an example to illustrate this. The subsample of size n = 100 was selected by the DS-WR algorithm from a data set D with N = 3,000 and q = 2. The open circles are the selected points in ${\mathcal{S}}_n$, and the numbers beside them indicate how many times each point was selected. There are only 14 unique points in the subsample of size 100, and over half of the 100 points are repeated values of only three ${\mathit{x}}_i$ points that were repeated 38, 13, and 10 times, respectively. The DS algorithm (without replacement, see Procedure 1) would likely have selected each of the outlier points, but only once each, which seems more desirable in this case.

3. Generalization to Customized Nonuniform Subsamples

In this section, we discuss how to generalize the DS method to select subsamples following desired distributions other than uniform. We also provide an example using the generalized DS method to select a subsample with desired properties without using a well-defined target distribution. For use in various machine learning contexts, the generalized DS algorithm could be used to incorporate information on the response, if they are available. As noted in Joseph and Mak (2021), for some applications this may be more desirable than using a diverse subsample in the input space. The theoretical performance of the generalized DS algorithm is shown in Section 4.

Consider the same setting as in Section 2, and suppose we want to select an i.i.d. subsample from D following some desired distribution with specified density $g(\mathit{x})=c\tilde{g}(\mathit{x})$ having support $S_g\subset S$. Here $\tilde{g}(\mathit{x})$ is known but c (the constant necessary for $g(\mathit{x})$ to be a proper density) can be unknown. The generalized DS algorithm differs from the DS algorithm only in that each ${\mathit{x}}_i$ is selected with probability

for the generalized DS. Specifically, at iteration k = 1, the generalized DS algorithm selects the index $j_1\in \{1,\dots ,N\}$ of the first point ${\mathit{x}}_{j_1}$, with probability

At iteration $k\ge 2$, the generalized DS algorithm selects the index $j_k\in \{1,\dots ,N\}\backslash \{j_1,\dots ,j_{k-1}\}$ of the next sampled point ${\mathit{x}}_{j_k}$ with probability

Consequently, Procedure 1 applies to the generalized DS algorithm with only line 15 and line 26 modified according to Equation (4) and Equation (5), respectively. Because N is finite in practice, the performance of the generalized DS algorithm will degrade if the tail of $f(\mathit{x})$ happens to be the mode of $g(\mathit{x})$, because in D there may not be enough points to choose to form a subsample with density $g(\mathit{x})$. The same issue is of course present in the DS algorithm, because the uniform target density also will be much larger than $f(\mathit{x})$ in its tail regions. Thus, the points in D falling in the tail regions of $f(\mathit{x})$ will be depleted too quickly to allow for a uniform subsample over the tails. We consider this phenomenon in our examples in Section 5, where we also observe that this is a unavoidable problem for any existing diversity subsampling algorithm.

We now provide an example using the generalized DS algorithm in the case that $g(\mathit{x})$ is not a well-defined p.d.f. Balancing the diversity property in the input space and some criterion in the response space has been found beneficial for some applications in the literature (Joseph and Mak 2021, Ren et al. 2021). For regression problems, Joseph and Mak (2021) observed that balancing diversity in the input space and the loss function in the response space makes the selected subsample more robust with respect to different modeling strategies. For classification problems in active learning, uncertainty and diversity are two desirable properties of the selected subsample. Many works have found that a subsample incorporating both properties can be more beneficial (Ren et al. 2021). Here we provide a simple classification example illustrating how to use the generalized DS for such purposes.

Consider the following binary classification problem with q = 2. The data set D of predictor vectors is the same as in the MGM example described in Section 5.2, and we generate the binary response variable Y ($=0$ or 1) via

We use the Euclidean norm of the gradient of Equation (6) as a simple surrogate measure of classification uncertainty (a large gradient means the predicted probability $P(Y=1|\mathit{X})$ is changing rapidly in that region, which usually translates to higher classification uncertainty) and denote it as $u(\mathit{x})\stackrel{\text{def}}{=}{\Vert {\nabla }_{\mathit{X}}P(Y=1|\mathit{X})\Vert }_2$. To incorporate both uncertainty and diversity into the subsample, we let $g(\mathit{x})=u(\mathit{x})+u_{\alpha }$, where $u_{\alpha }$ is a constant (and thus represents a uniform distribution to promote the diversity, whereas $u(\mathit{x})$ considers the classification uncertainty) that we take to be the lower α quantile of set ${\{u({\mathit{x}}_i)\}}_{i=1}^N$ for some specified α. Here $g(\mathit{x})$ typically will not be a well-defined p.d.f. For the generalized DS algorithm, the sampling probability of each data point will be proportional to $g({\mathit{x}}_i)/f({\mathit{x}}_i)$. Selecting a larger α will promote more diversity and selecting a smaller α will result in more subsample points chosen in the areas with higher uncertainty. Figure 3 shows the selected subsamples for various α. The colorbar shows the values of Equation (6). We observe that when α = 0, most of the selected points locate close to the decision boundary $\mathit{x}:P(Y=1|\mathit{x})=0.5$ (Figure 3(a)), where the highest uncertainty occurs in this model. Using a larger α, such as $\alpha =50\%$ allows more subsample points in other regions in the predictor space (Figure 3(c)). When $\alpha =100\%$, the subsample appears to be quite diverse.

4. Theoretical Performance Analysis

This section presents results on the convergence rate of the generalized DS algorithm (Section 3) when the true $f(\mathit{x})$ is known (Theorem 1). For convenience, we say that a function u(N) of N is $O(N^r)$ if $\exists c>0,\hspace{0.17em}k>0$, ($c,\hspace{0.17em}k$ are both constants that do not depend on N), such that $|u(N)|\le c{N}^r,\hspace{0.17em}\forall N\ge k$. We say u(N) is $o(N^r)$ if $\forall c>0,\hspace{0.17em}\exists k(c)>0$ that does not depend on N, such that $|u(N)|<c{N}^r,\hspace{0.17em}\forall N\ge k(c)$.

When the true p.d.f of the data distribution is available, the generalized DS algorithm (Section 3) is the same as the SIR algorithm (Rubin 1987, 1988), for which Skare et al. (2003) proves convergence and which we restate as Theorem 1.

Theorem 1

(Theorem 2, Part 1 in Skare et al. 2003). Let ${\mathit{X}}^N=({\mathit{X}}_i:i=1,2,\dots ,N)$, where ${\mathit{X}}_1,{\mathit{X}}_2,\dots ,{\mathit{X}}_N\in {\mathbb{R}}^q$ are i.i.d. copies of random vector X having density function $f(\mathit{x})$ with support $S\subset {\mathbb{R}}^q$. Let $g(\mathit{x})=c\tilde{g}(\mathit{x})$ be the desired density function for the selected subsample that is known up to a constant $c\in {\mathbb{R}}^{+}$. Assume that

• A.1. The support of $g(\mathit{x}),\hspace{0.17em}{S}_g\subset S$

• A.2. The variance of $\frac{g(\mathit{X})}{f(\mathit{X})},\hspace{0.17em}\mathit{\text{Var}}[\frac{g(\mathit{X})}{f(\mathit{X})}]={\int }_S\frac{g^2(\mathit{x})}{f(\mathit{x})}\text{d}x-1={\sigma }^2<+\infty$

Suppose $n\ge 1$ is fixed. For each $N=1,2,\dots$, let $({\mathit{Z}}_N^k:k=1,2,\dots ,n)$ be drawn sequentially from ${\mathit{X}}^N$, without replacement, according to the conditional probability mass function (p.m.f.) (for k = 1)

$$P({\mathit{Z}}_N^1={\mathit{X}}_i|{\mathit{X}}^N)=\frac{\frac{\tilde{g}(\mathit{x})}{f({\mathit{X}}_i)}}{{\sum }_{j=1}^N\frac{\tilde{g}(\mathit{x})}{f({\mathit{X}}_{\mathit{j}})}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(i=1,\dots ,N),$$(7)

and p.m.f. (for $k\ge 2$)

$$P({\mathit{Z}}_N^k={\mathit{X}}_i|{\mathit{X}}^N,{\mathit{Z}}_N^1,\dots ,{\mathit{Z}}_N^{k-1})=\frac{\frac{\tilde{g}({\mathit{X}}_i)}{f({\mathit{X}}_i)}}{{\sum }_{j=1}^N\frac{\tilde{g}({\mathit{X}}_j)}{f({\mathit{X}}_j)}-{\sum }_{j=1}^{k-1}\frac{\tilde{g}({\mathit{Z}}_N^j)}{f({\mathit{Z}}_N^j)}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(i\in \{1,2,\dots ,N\}\backslash \{j_1^N,j_2^N,\dots ,j_{k-1}^N),$$(8)

where $j_1^N,j_2^N,\dots ,j_n^N$ are the indices of the sampled observations in ${\mathit{X}}^N$, that is, ${\mathit{Z}}_N^k={\mathit{X}}_{j_k^N}$, for $k=1,2,\dots ,n$. Denote the joint p.d.f. of ${\mathit{Z}}_N^1,\dots ,{\mathit{Z}}_N^n$ as $p_N^n({\mathit{z}}_1,\dots ,{\mathit{z}}_n)$. Then as $N\to +\infty$,

$$\frac{p_N^n({\mathit{z}}_1,\dots ,{\mathit{z}}_n)}{{\prod }_{k=1}^{n}g({\mathit{z}}_k)}-1=\frac{1}{N}({\displaystyle \sum _{k=1}^{n}k}(1-\frac{g({\mathit{z}}_k)}{f({\mathit{z}}_k)}+{\sigma }^2))+O\left(\frac{1}{N^2}\right)\to 0.$$(9)

Furthermore, when ${\mathit{\text{sup}}}_{\mathit{x}\in S}\frac{g(\mathit{x})}{f(\mathit{x})}<+\infty$, the convergence of $p_N^n({\mathit{z}}_1,\dots ,{\mathit{z}}_n)$ is uniform.

In practice, $f(\mathit{x})$ must be estimated, in which case the generalized DS algorithm may not have the convergence properties as SIR, due to inaccuracies in estimating $f(\mathit{x})$. It is not unreasonable to assume that if $\widehat{f}(\mathit{x})$ converges to $f(\mathit{x})$ at a sufficient rate, then the sample produced by the DS algorithm using $\widehat{f}(\mathit{x})$ will converge to the target distribution. Although we were unable to prove this under general convergence conditions that are known to hold for our estimator of $f(\mathit{x})$, preliminary results indicate that convergence of the DS algorithm holds for certain specified conditions on the density estimator. More importantly, our numerical results presented later demonstrate that the DS algorithm with estimated $f(\mathit{x})$ does indeed produce samples that effectively approximate the target distribution.

Remark 1.

Theorem 1 apply to the DS algorithm (Procedure 1) for producing a uniform sample by setting $\tilde{g}(\mathit{x})=1,\hspace{0.17em}\forall \mathit{x}\in S_g$.

5. Numerical Performance Analysis

In this section, we test the performance of the DS algorithm using data sets following various distributions and compare it with competing subsampling methods. Section 5.1 applies various subsampling algorithms to a real-world data set and visually compares the diversity of the selected subsamples. Section 5.2 discusses numerical performance evaluation criteria and the experimental settings for the comparative examples. Section 5.3 provides numerical results comparing uniformity performances of DS with competing algorithms, and Section 5.4 compares runtime of the algorithms.

5.1. Example: Bike-Sharing Data

We apply DS, scSampler-sp1 (Song et al. 2022b), and random sampling to the open-source bike-sharing data set used in Fanaee-T and Gama (2014). The data set is of size 17,379, and we included the following 6 variables in our experiments: hr, holiday, weathersit, atemp, hum, and windspeed. Here hr records at which hour of a day the number of rented bikes was counted; holiday indicates whether the day is a holiday; weathersit uses four grades to record the weather situation; atemp stands for apparent temperature; hum records humidity and windspeed stores the wind speed. A subsample of size n = 300 was selected from the data set using each method (we standardized the data into a unit hypercube for DS and scSampler-sp1), and Figure 4 show projected two-dimensional (2D) scatter plots in the hum-hr and hum-atemp subspaces. We observe from Figure 4 that the subsamples selected by DS and scSampler-sp1 appear more spread-out and diverse than the one selected by random sampling. For example, both DS and scSampler-sp1 selected points at the lowest hum level, whereas random sampling did not. Although for brevity we omit projected scatter plots for other pairs of variables, similar conclusions apply for them.

5.2. Experimental Settings

As a performance evaluation criterion, we use the sample version of energy distance (Székely 2003, Székely and Rizzo 2004) to measure the extent to which a subsample is drawn from an i.i.d. uniform distribution over some specified $\mathrm{\Omega }\subset S$, where Ω is compact. We provide further details of the choice of Ω below. Let $\mathit{X},\mathit{U}$ be mutually independent random variables in ${\mathbb{R}}^q$. The sample version of energy distance is as follows. Let ${\mathcal{A}}_n=\{{\mathit{X}}_1,\dots ,{\mathit{X}}_n\}$ and $\mathcal{U}=\{{\mathit{U}}_1,\dots ,{\mathit{U}}_{N_1}\}$ be i.i.d. samples of size n and N₁, respectively, from the same distributions followed by X and U, respectively. The energy distance between samples ${\mathcal{A}}_n$ and $\mathcal{U}$ is (Székely 2003)

The Euclidean norm will be used throughout this paper unless otherwise specified. By proposition 1 in Székely (2003), as $n\to \infty$ and $N_1\to \infty$, one can show that $E[e({\mathcal{A}}_n,\mathcal{U})]\to 0+$ if and only if $\mathit{X}\stackrel{\text{D}}{=}\mathit{U}$. Here $\stackrel{\text{D}}{=}$ denotes equality in distribution.

We use Equation (10) in the following way. We choose $\mathcal{U}$ to be a large uniform sample over Ω of size $N_1\gg n$, and we take ${\mathcal{A}}_n$ to be the subsample of size n selected from D using some subsampling algorithm. If points in subsample ${\mathcal{A}}_n$ (subsample points located outside of Ω are excluded under this criterion) follow an i.i.d. uniform distribution over Ω, $E[e({\mathcal{A}}_n,\mathcal{U})]$ should be close to zero if n is sufficiently large (we take a fixed large N₁ in this paper), and the more ${\mathcal{A}}_n$ differs from an i.i.d. uniform distribution over Ω, the larger we expect $E[e({\mathcal{A}}_n,\mathcal{U})]$ to be. For selected experiments, we also report the p values for testing the null hypothesis that ${\mathcal{A}}_n$ and $\mathcal{U}$ are generated from the same distribution using the test proposed in Székely and Rizzo (2004).

Intuitively, to produce a diverse and space-filling subsample, it is typically desirable to select all or most of the points in the regions for which $f(\mathit{x})$ has low-density, which we denote by ${\mathrm{\Omega }}^c$, where $S=\mathrm{\Omega }\cup {\mathrm{\Omega }}^c$. Because we cannot expect uniformity over both Ω and ${\mathrm{\Omega }}^c$ for large n (unless N is extremely large) as the points in ${\mathrm{\Omega }}^c$ will be quickly depleted, we use a separate measure for performances within ${\mathrm{\Omega }}^c$. Because data are sparse in ${\mathrm{\Omega }}^c$, it is generally desirable to select as high a percentage of points in ${\mathrm{\Omega }}^c$ as possible. Consequently, we define the following low-density region sampling ratio to measure this property. Suppose Ω is chosen such that $\forall \mathit{x}\in \mathrm{\Omega },\hspace{0.17em}f(\mathit{x})\ge \delta$ and $\forall \mathit{x}\in {\mathrm{\Omega }}^c,\hspace{0.17em}f(\mathit{x})<\delta$, for some specified value δ. Then for a subsample ${\mathcal{S}}_n$ of size n selected from D, the low-density region sampling ratio of ${\mathcal{S}}_n$ is defined as

where ${\mathit{I}}_{\{·\}}$ denotes the indicator function. For instance, ${\mathit{I}}_{\{{\mathit{x}}_{j_k}\in {\mathrm{\Omega }}^c\}}$ equals one if ${\mathit{x}}_{j_k}\in {\mathrm{\Omega }}^c$, and equals zero otherwise. Only unique subsample indices are used in Equation (11). In situations where a uniform subsample over S of size n will exhaust points in ${\mathrm{\Omega }}^c$, a larger $r({\mathcal{S}}_n)$ is desirable, with $r({\mathcal{S}}_n)=1$ typically being the ideal goal (in conjunction with a uniform distribution over Ω).

For the examples, we consider data sets in which each element of X follows independent Gaussian, gamma, exponential, geometric, and multivariate Gaussian mixture (MGM) distributions. We refer to these five different data distributions as the normal, gamma, exponential, geometric, and MGM. The normal example illustrates the performances of the subsampling algorithms when the data distribution is symmetric; the mode of the density is achieved in the interior region of S; S is unbounded. The gamma example illustrates the performances when the data distribution is skewed and S had boundaries but is unbounded. The exponential distribution, which is an extreme case of the gamma distribution, illustrates the performances when the data distribution is extremely skewed and the mode is located at the boundary. The MGM examples are to compare the subsampling algorithms when the data distribution is multimodal with correlated predictors. We also include geometric examples to illustrate the performances of different algorithms when the data distribution is discrete. For all examples in this section, the size of D is $N={10}^4$ for q = 2 and $N={10}^5$ for q = 10.

The data distributions of each simulation example are as follows. For the normal, gamma, and exponential examples, the elements of X are independent and the data distribution has density $f(\mathit{x})={\prod }_{j=1}^{q}f(x_j)$, where $f(x)=e^{-x^2/2}/\sqrt{2\pi },\hspace{0.17em}f(x)=x{e}^{-x/2}/(4\mathrm{\Gamma }(2))$ and $f(x)=e^{-x}$, respectively, where $\mathrm{\Gamma }()$ denotes the gamma function. For the MGM example, D was generated as a mixture of two Gaussian clusters in ${\mathbb{R}}^q$. The true density function is as follows:

wherefor i = 1, 2. That is, the two Gaussian clusters have a common covariance matrix $\mathbf{\Sigma }$ but different means ${\mathit{\mu }}_1,{\mathit{\mu }}_2$. We let ${\mathit{\mu }}_1={\mathbf{0}}_q$ and ${\mathit{\mu }}_2=5.0{({(-1)}^i)}_{i=1}^q$. For the common covariance matrix $\mathbf{\Sigma }$, we set $\mathbf{\Sigma }={\sigma }^2{\mathit{I}}_{q\times q}+\alpha \mathit{a}{\mathit{a}}^T$. Here ${\sigma }^2=4.0,\hspace{0.17em}\alpha =1.0,\hspace{0.17em}{\mathit{I}}_{q\times q}\in {\mathbb{R}}^{q\times q}$ is the identity matrix, and $\mathit{a}={(a_i)}_{i=1}^q$, for $a_i=0.2(i-2){(-1)}^i$. The data distribution of the geometric example also has independent covariates; thus the joint p.m.f. is $f(\mathit{x})={\prod }_{j=1}^{q}f(x_j)$, where $f(x)={(1-p)}^{x-1}p$. We set p = 0.5 for q = 2 and p = 0.9 for q = 10. The supports of the p.d.f. functions of the normal, gamma, exponential, MGM, and geometric examples are, respectively, $S={\mathbb{R}}^q,{\mathbb{R}}_{+}^q,{\mathbb{R}}_{+}^q,{\mathbb{R}}^q,{\mathbb{Z}}_{+}^q$.

We choose $\mathrm{\Omega }=\{\mathit{x}\in S|f(\mathit{x})\ge \delta \}$, where δ is chosen such that $99\%$ of the data points in D fall inside of Ω when q = 2 and 99.9% of the data points in D do so when q = 10. Figure 5 depicts Ω for selected examples with q = 2. To be clear, for discrete distributions, the target uniform distribution is a discrete uniform distribution over S. This distinction between continuous and discrete distributions is not used in any way in the DS algorithm, although it is used in our performance measures.

5.3. Numerical Results for Performance Comparisons

The competing algorithms that we compare include scSampler-sp1, scSampler-sp16, DPP, WSE, (weighted sample elimination algorithm, a heuristic for PDS proposed by Yuksel 2015), and simple random sampling. Here scSampler-sp1 (Song et al. 2022b) serves as an efficient heuristic for the CADEX algorithm. scSampler-sp16 (Song et al. 2022b) is included as a modification of scSampler-sp1 for better computational efficiency. We used published codes for scSampler-sp1 and scSampler-sp16 (Song et al. 2022a), DPP (Biyik 2019), and WSE (Yuksel 2016). We omit the results for RD ALR in this paper, as we found that it generally did not perform as well as other methods.

The implementation details for competing methods are as follows. For the two versions of the scSampler algorithm, no additional parameters need to be set. For DPP, we set hyperparameters $\alpha =4.0,\hspace{0.17em}\gamma =0.0$ and steps = 0 (Biyik et al. (2019) recommended setting $\alpha \ge 2.0$ and $\gamma =0.0$ in their paper to get the most diverse subsample; we chose steps = 0 for computational efficiency at the cost of the resulted subsample being deterministic). For WSE, we used the default weighting function in Yuksel (2016) and set Progressive = True to make the selected subsample sequential.

Figure 6 plots the average value of $e({\mathcal{S}}_n,\mathcal{U})$ across 200 replicates of the experiments (its standard error across the 200 replicates was negligibly small and is omitted from the figures). For DPP and WSE, only a single replicate was used, since WSE and DPP (with steps = 0) are deterministic and produce the same subsamples each replicate. The results for the low density region sampling ratio (Equation (11)) are consistent with the $e({\mathcal{S}}_n,\mathcal{U})$ results and are shown in Appendix D. The size of the large uniform sample $\mathcal{U}$ inside Ω is the same as $|D|$. As a reference point, we have included the $e({\mathcal{S}}_n,\mathcal{U})$ measures for true uniform samples (denoted as uniform sampling) inside Ω.

We assessed the performances of the subsampling algorithms for $n=20,520,\dots ,4020$ at q = 2 (DPP was only tested for n = 20, 250,…1,860 due to expensive computational cost; see Section 5.4) and for $n=200,900,\dots ,6500$ at q = 10.¹ DPP was not included in the comparisons at q = 10 as its computation costs was excessive.

Figure 6 shows that average value of $e({\mathcal{S}}_n,\mathcal{U})$ for uniform sampling converges to zero as n increases, whereas for most subsampling algorithms it first decreases and then increases instead of monotonically decreasing with n. This is because uniform sampling generates sample points from Ω directly, whereas the subsampling algorithms select points from a finite D without replacement. After points in D in the lower-density regions in Ω are depleted, the subsampling algorithms can no longer produce uniform samples. We discuss this phenomenon in more detail in Appendix C.

We also observe from Figure 6 that some algorithms have an average $e({\mathcal{S}}_n,\mathcal{U})$ measure that is actually less than uniform sampling, such as DPP and WSE for the normal example at q = 2 and $n\le 1000$ (Figure 6(a)). This is a well-known phenomenon, by which points on a uniform grid can have a lower $e({\mathcal{S}}_n,\mathcal{U})$ measure than a truly uniform sample (Mak and Joseph 2018). For instance, in separate experiments (not shown here, for brevity), we found that, a grid of size n = 100 in ${[0,1]}^2$ often has a smaller sample energy distance to $\mathcal{U}$ than a true uniformly distributed sample. Consequently, we can view a subsampling algorithm whose performance is as close to the uniform sampling as possible as the most effective method in terms of selecting i.i.d. uniformly distributed subsamples over Ω. From Figure 6, we see that overall, the average sample energy distances of DS deviates least from the uniform sampling results for most cases compared with other subsampling algorithms.

For the normal examples in 2D and 10D, we also test whether the subsamples selected by various algorithm are equal in distribution to $\mathcal{U}$ using the equal-distribution test proposed by Székely and Rizzo (2004) and implemented as the eqdist.etest function in the R package energy (Rizzo and Szekely 2022). The value of parameter R for the eqdist.etest function was set to be 200 (i.e., the precision of the estimated p value was set to $5\times {10}^{-3}$), and the resulting p values are shown by Figure 7. For computational reasons, the size of $\mathcal{U}$ was set to 10,000 for this experiment. Figure 7 appears consistent with Figure 6, (a) and (b): In Figure 7(a), DS had the largest p value for most n, and in Figure 7(b), only WSE had a larger p value for most n. However, as discussed in Section 5.4, WSE has far larger (often prohibitively large) computational expense than DS.

We also repeated the experiments in Figure 6 using data sets in which observations were replicated for the continuous distributions. Specifically, we generated a data set of size $N/5$ and then replicated this data set five times to generate the set D of size N, which we refer to as replicated data sets. The results using the replicated data sets are shown by Figure 8. Compared with Figure 6, the DS performance has barely changed, whereas the scSampler-sp1 performance has degraded significantly. We conclude that the performance of DS relative to the existing methods further improves, often substantially, for replicated data. The results in Figure 8(a) are consistent with the visualization in Figure 1: at n = 2,000, the performances of scSampler-sp1 are comparable with random sampling, whereas the performances of DS are substantially better.

5.4. Computation Time Comparisons

Table 1 provides runtimes for the DS, scSampler-sp1, scSampler-sp16, DPP, and WSE² algorithms. We implemented the DS algorithm in Python 3 and tested the runtimes of the previously mentioned algorithms on the Quest High Performance Computing Cluster of Northwestern University. The runtimes were for the normal examples (runtime is relatively invariant to the underlying distribution) with varying data set sizes. The reported runtime statistics in Table 1 are in seconds and were averaged over 100 independent replications, although there was not much replicate-to-replicate variability. NA indicates that a subsampling algorithm took longer than one hour to run (for each replicate), in which case the experiment was terminated, and the approach was considered to be computationally prohibitive for that case.

Table 1. Runtime Comparisons of Subsampling Algorithms

From Table 1, DS and scSampler-sp16 are far more efficient than other methods, often multiple orders of magnitude faster. One should keep in mind, however, that the uniformity performance of scSampler-sp16 was overall far inferior to DS in Figures 6 and 8. For n ≥ 1,000, the runtime of DS remains relatively unchanged with varying n, whereas the runtime of scSampler-sp16 (and of scSampler-sp1 and DPP) grows super linearly with n increasing. Consequently, for the largest n in Table 1, DS becomes faster than scSampler-sp16. Moreover, for larger scale examples with even larger n, we can expect DS to be substantially faster than even scSampler-sp16 (in addition to having substantially better performance). DPP and WSE were prohibitively slow for even the smallest size examples in Table 1 and were too slow to be included in the table for the larger size examples.

The DS algorithm mainly consists of three parts: estimating the density of D upfront from scratch (line 12 in Procedure 1), updating the density periodically (line 23 in Procedure 1) and selecting subsample points (line 15 and line 26 in Procedure 1). For convenience, we will refer these three parts to initial density estimation, density updating, and subsample selection, respectively. We observe that for all tested examples, initial density estimation and density Updating cost around 95% of the total computation time, and the computational costs of initial density estimation and density updating are roughly comparable (results are not shown here for brevity). Considering that the density estimation accounts for most of the computational cost, for large N the DS expense could be reduced by randomly sampling some fraction of D to use for the density estimation. This is somewhat similar to the computational strategy used by scSampler-sp16 to reduce computational expense relative to scSampler-sp1. The former randomly partitions D into 16 subsets and applies scSampler-sp1 to each subset and then combines the 16 subsamples.

6. Conclusions

This paper presents a novel diversity subsampling algorithm, the DS algorithm, that selects an i.i.d. uniform subsample from a data set over the effective support of the empirical data distribution to the largest extent possible. The asymptotic performance of the DS algorithm was presented for the case that $f(\mathit{x})$ is known, and its effectiveness and advantages over existing diversity subsampling algorithms were demonstrated numerically for the intended situation in which $f(\mathit{x})$ must be estimated: Overall, the DS algorithm selects subsamples more similar to a true i.i.d. uniform sample than existing algorithms do with much lower computational cost. We have also presented a generalized version of the DS algorithm to select subsamples following any desired target density (or nonnegative target function in general). All methods proposed in the paper are implemented in the FADS (Shang et al. 2022) Python package.