A Statistical Perspective on Coreset Density Estimation

We formally define a coreset as follows. Throughout this work $m=o(n)$ denotes the cardinality of the coreset. Let $S=S(y|x)$ denote a conditional probability measure on $\binom{[n]}{m}$, where $x\in\mathbb{R}^{d\times n}$. In information theoretic language, $S$ is a channel from $\mathbb{R}^{d\times n}$ to subsets of cardinality $m$. We refer to the channel $S$ as a coreset scheme because it designates a data-driven method of choosing a subset of data points. In what follows, we abuse notation and let $S=S(x)$ denote an instantiation of a sample from the measure $S(y|x)$ for $x\in\mathbb{R}^{d\times n}$. A coreset $X_{S}$ is then defined to be the projection of the dataset $X=(X_{1},\ldots,X_{n})$ onto the subset indicated by $S(X)$: $X_{S}:=\{X_{i}\}_{i\in S(X)}$.

The first family of estimators that we investigate is quite general and allows the statistician to select a coreset and then employ an estimator that only manipulates data points in the coreset to estimate an unknown density. To study coresets, it is convenient to make the dependence of estimators on observations more explicit than in the traditional literature. More specifically, a density estimator $\hat{f}$ based on $n$ observations $X_{1},\ldots,X_{n}\in\mathbb{R}^{d}$ is a function $\hat{f}:\mathbb{R}^{d\times n}\to L_{2}(\mathbb{R}^{d})$ denoted by $\hat{f}[X_{1},\ldots,X_{n}](\cdot)$. Similarly, a coreset-based estimator $\hat{f}_{S}$ is constructed from a coreset scheme $S$ of size $m$ and an estimator (measurable function) $\hat{f}:\mathbb{R}^{d\times m}\to L_{2}(\mathbb{R}^{d})$ on $m$ observations. We enforce the additional restriction on $\hat{f}$ that for all $y_{1},\ldots,y_{m}\in\mathbb{R}^{d}$ and for all bijections $\pi:[m]\to[m]$, it holds that $\hat{f}[y_{1},\ldots,y_{m}](\cdot)=\hat{f}[y_{\pi(1)},\ldots,y_{\pi(m)}](\cdot)$. Given $S$ and $\hat{f}$ as above, we define the coreset-based estimator $\hat{f}_{S}:\mathbb{R}^{d\times n}\to L_{2}(\mathbb{R}^{d})$ to be the function $\hat{f}_{S}[X](\cdot):=\hat{f}[X_{S}](\cdot):\mathbb{R}^{d}\to\mathbb{R}$. We evaluate the performance of coreset-based estimators in Section 2 by characterizing their rate of estimation over Hölder classes.¹¹1Our notion of coreset-based estimators bares conceptual similarity to various notions of compression schemes as studied in the literature, e.g. Littlestone and Warmuth (1986); Ashtiani et al. (2020); Hanneke et al. (2019).

The symmetry restriction on $\hat{f}$ prevents the user from exploiting information about the ordering of data points to their advantage: the only information that can be used by the estimator $\hat{f}$ is contained in the unordered collection of distinct vectors given by the coreset $X_{S}$.

As evident from the the results in Section 2, the information-theoretically optimal coreset estimator does not resemble coreset estimators employed in practice. To remedy this limitation, we also study weighted coreset kernel density estimators (KDEs) in Section 3. Here the statistician selects a kernel $k$, bandwidth parameter $h$, and a coreset $X_{S}$ of cardinality $m$ as defined above and then employs the estimator

$$\hat{f}_{S}(y)=\sum_{j\in S}\lambda_{j}h^{-d}k\left(\frac{X_{j}-y}{h}\right),$$

where the weights $\{\lambda_{j}\}_{j\in S}$ are nonnegative, sum to one and are allowed to depend on the full dataset.

In the case of uniform weights where $\lambda_{j}=\frac{1}{m}$ for all $j\in S$, coreset KDEs are well-studied (see e.g. Bach et al., 2012; Harvey and Samadi, 2014; Phillips and Tai, 2018a, b; Karnin and Liberty, 2019). Interestingly, our results show that allowing flexibility in the weights gives a definitive advantage for the task of density estimation. By Theorems 2 and 5, the uniformly weighted coreset KDEs require a much larger coreset than that of weighted coreset KDEs to attain the minimax rate of estimation over univariate Lipschitz densities.

We reserve the notation $\lVert\cdot\rVert_{2}$ for the $L_{2}$ norm and $\left|\cdot\right|_{p}$ for the $\ell_{p}$-norm. The constants $c,c_{\beta,d},c_{L},$ etc. vary from line to line and the subscripts indicate parameter dependences.

Fix an integer $d\geq 1$. For any multi-index $s=(s_{1},\ldots,s_{d})\in\mathbb{Z}_{\geq 0}^{d}$ and $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$, define $s!=s_{1}!\cdots s_{d}!$, $x^{s}=x_{1}^{s_{1}}\cdots x_{d}^{s_{d}}$ and let $D^{s}$ denote the differential operator defined by

$$D^{s}=\frac{\partial^{|s|_{1}}}{\partial x_{1}^{s_{1}}\cdots\partial x_{d}^{s_{d}}}\,.$$

Fix a positive real number $\beta,$ and let $\lfloor\beta\rfloor$ denote the maximal integer strictly less than $\beta$. Given a multi-index $s$, the notation $\left|s\right|$ signifies the coordinate-wise application of $\left|\cdot\right|$ to $s$.

Given $L>0$ we let $\mathcal{H}(\beta,L)$ denote the space of Hölder functions $f:\mathbb{R}^{d}\to\mathbb{R}$ that are supported on the cube $[-1/2,1/2]^{d}$, are $\lfloor\beta\rfloor$ times differentiable, and satisfy

$$|D^{s}f(x)-D^{s}f(y)|\leq L\left|x-y\right|_{2}^{\beta-\lfloor\beta\rfloor}\,,\quad$$

for all $x,y\in\mathbb{R}^{d}$ and for all multi-indices $s$ such that $\left|s\right|_{1}=\lfloor\beta\rfloor$.

Let $\mathcal{P}_{\mathcal{H}}(\beta,L)$ denote the set of probability density functions contained in $\mathcal{H}(\beta,L)$. For $f\in\mathcal{P}_{\mathcal{H}}(\beta,L)$, let $\mathbb{P}_{f}$ (resp. $\mathbb{E}_{f}$) denote the probability distribution (resp. expectation) associated to $f$.

For $d\geq 1$ and $\gamma\in\mathbb{Z}_{\geq 0}$, we also define the Sobolev functions $\mathcal{S}(\gamma,L^{\prime})$ that consist of all $f:\mathbb{R}^{d}\to\mathbb{R}$ that are $\gamma$ times differentiable and satisfy

$$\lVert D^{\alpha}f\rVert_{2}\leq L^{\prime}$$

for all multi-indices $\alpha$ such that $\left|\alpha\right|_{1}=\gamma$.

Fix $\varepsilon=c^{*}(m\log n)^{-\frac{\beta}{d}}$ for $c^{*}$ to be determined and let $\mathcal{N}_{\varepsilon}$ denote an $\varepsilon$-net of $\mathcal{P}_{\mathcal{H}}(\beta,L)$ with respect to the $L_{2}([-\frac{1}{2},\frac{1}{2}]^{d})$ norm. It follows from the classical Kolmogorov-Tikhomorov bound (see, e.g., Theorem XIV of Tikhomirov, 1993) that there exists a constant $C_{\mathsf{KT}}(\beta,d,L)>0$ such that we can choose $\mathcal{N}_{\varepsilon}$ with $\log|\mathcal{N}_{\varepsilon}|\leq C_{\mathsf{KT}}(\beta,d,L)\,\varepsilon^{-d/\beta}$. In particular, there exists $\mathsf{f}\in\mathcal{N}_{\varepsilon}$ such that $\|\hat{f}_{n}-\mathsf{f}\|_{L^{2}([-1/2,1/2]^{d})}\leq\varepsilon$ where $\hat{f}_{n}$ is the minimax optimal kernel density estimator defined in (2).

We now develop our encoding procedure for $\mathsf{f}$. To that end, fix an integer $K\geq m$ such that $\binom{K}{m}\geq|\mathcal{N}_{\varepsilon}|$ and let $\phi:\binom{[K]}{m}\to\mathcal{N}_{\varepsilon}$ be any surjective map. Our procedure only looks at the first coordinates of the sample $X=\{X_{1},\ldots,X_{n}\}$. Denote these coordinates by $x=\{x_{1},\ldots,x_{n}\}$ and note that these $n$ numbers are almost surely distinct. Let $A$ denote a parameter to be determined, and define the intervals

$$B_{ik}=[(i-1)K^{-1}A+(k-1)A,\\ \quad(i-1)K^{-1}A+(k-1)A+K^{-1}A].$$

For $i=1,\ldots,K$, define

$$B_{i}=\bigcup_{k=1}^{1/A}B_{ik}.$$

The next lemma, whose proof is in the Appendix, ensures that with high probability every bin $B_{i}$ contains the first coordinate $x_{i}$ of at least one data point.

In the high-probability event $\mathcal{E}$ that every bin $B_{i}$ contains the first coordinate of some data point, choose a unique representative $x_{j}^{\circ}\in x$ such that $x_{j}^{\circ}\in B_{j}$ and pick any $T_{\mathsf{f}}\in\phi^{-1}(\mathsf{f})$. Then define $S=\{i\,:\,x_{i}=x^{\circ}_{j},j\in T_{\mathsf{f}}\}$. If there exists a bin with no observation, then let $X_{S}$ consist of two data points lying in the same bin and $m-2$ random data points. Then set $\hat{f}_{S}\equiv 0$.

Note that $\hat{f}_{S}$ is indeed a coreset-based estimator. The function $\hat{f}$ such that $\hat{f}_{S}=\hat{f}[X_{S}]$ looks at the $m$ data points in the coreset, and if their first coordinates lie in distinct bins, then $X_{S}$ is decoded as above to output the corresponding element $\mathsf{f}$ of the net $\mathcal{N}_{\varepsilon}$. Otherwise, $\hat{f}\equiv 0$.

Next, it suffices to show the upper bound of Theorem 1 in the case when $m\leq cn^{d/(2\beta+d)}$ for $c$ a sufficiently small absolute constant. For $c^{*}=c^{*}_{\beta,d,L}$ sufficiently large, by Stirling’s formula and our choice of $K$ it holds that

$$\log\binom{K}{m}\geq C_{\mathsf{KT}}(\beta,d,L)\,\left(\frac{1}{\varepsilon}\right)^{\frac{d}{\beta}}\geq\log|\mathcal{N}_{\varepsilon}|.$$

Hence, the surjection $\phi$ and our encoding estimator $\hat{f}_{S}$ are well-defined.

Next we have

$$\mathbb{E}_{f}\lVert\hat{f}_{S}-f\rVert_{2}=\mathbb{E}_{f}\big{[}\lVert\mathsf{f}-f\rVert_{2}{\rm 1}\kern-2.40005pt{\rm I}_{\mathcal{E}}\big{]}+\mathbb{E}_{f}\big{[}\lVert 0-f\rVert_{2}{\rm 1}\kern-2.40005pt{\rm I}_{\mathcal{E}^{c}}\big{]}.$$

We control the first term as follows using (1) and the fact that $\lVert\mathsf{f}-\hat{f}_{n}\rVert_{2}\leq\varepsilon$ on $\mathcal{E}$:

	$\displaystyle\mathbb{E}_{f}\big{[}\lVert\mathsf{f}-f\rVert_{2}{\rm 1}\kern-2.40005pt{\rm I}_{\mathcal{E}}\big{]}$	$\displaystyle\leq\mathbb{E}_{f}\lVert\hat{f}_{n}-f\rVert_{2}+\mathbb{E}_{f}\lVert\mathsf{f}-\hat{f}_{n}\rVert_{2}$
		$\displaystyle\leq c_{\beta,d,L}\,\big{(}n^{\frac{-\beta}{2\beta+d}}+(m\log n)^{-\frac{\beta}{d}}\big{)}.$

By the Cauchy-Schwarz inequality,

	$\displaystyle\mathbb{E}_{f}\big{[}\lVert 0-f\rVert_{2}{\rm 1}\kern-2.40005pt{\rm I}_{\mathcal{E}^{c}}\big{]}$	$\displaystyle\leq\big{(}\mathbb{E}_{f}\lVert f\rVert_{2}^{2}\,\mathbb{P}(\mathcal{E}^{c})\big{)}^{1/2}$
		$\displaystyle\leq c_{\beta,d,L}\,n^{-1}\,.$

Put together, the previous three displays yield the upper bound of Theorem 1.

Given a KDE with uniform weights and bandwidth $h$ defined by

$$\hat{f}(y)=\frac{1}{n}\sum_{j=1}^{n}k_{h}(X_{j}-y),$$

on a sample $X_{1},\ldots,X_{n}$, we define a coreset KDE $\hat{g}_{S}$ as follows in terms of a cutoff frequency $T>0$. Define $A=\{\omega\in\frac{\pi}{2}\mathbb{Z}^{d}:\,\left|\omega\right|_{\infty}\leq T\}$. Consider the complex vectors $(e^{i\langle X_{j},\omega\rangle})_{\omega\in A}$. By Carathéodory’s theorem (Carathéodory, 1907), there exists a subset $S\subset[n]$ of cardinality at most $2(1+\frac{4T}{\pi})^{d}+1$ and nonnegative weights $\{\lambda_{j}\}_{j\in S}$ with $\sum_{j\in S}\lambda_{j}=1$ such that

$$\frac{1}{n}\sum_{j=1}^{n}(e^{i\langle X_{j},\omega\rangle})_{\omega\in A}=\sum_{j\in S}\lambda_{j}(e^{i\langle X_{j},\omega\rangle})_{\omega\in A}.$$

(4)

Then $\hat{g}_{S}(y)$ is defined to be

$$\hat{g}_{S}(y)=\sum_{j\in S}\lambda_{j}k_{h}(X_{j}-y).$$

Proposition 1 is key to our results and specifies conditions on the kernel guaranteeing that the Carathéodory method yields an accurate estimator.

There exists a kernel $k_{s}\in\mathcal{C}^{\infty}$ that satisfies the conditions above for all $\beta$ and $\gamma$. We sketch the details here and postpone the full argument to the Proof of Theorem 2 in the Appendix. Let $\psi:[-1,1]\to[0,1]$ denote a cutoff function that has the following properties: $\psi\in\mathcal{C}^{\infty}$, $\psi\big{|}_{[-1,1]}\equiv 1$, and $\psi$ is compactly supported on $[-2,2]$. Define $\kappa_{S}(x)=\mathcal{F}[\psi](x)$, and let $k_{s}(x)=\prod_{i=1}^{d}\kappa_{S}(x_{i})$ denote the resulting kernel. Observe that for all $\beta>0$, the kernel $k_{s}$ satisfies

$$\text{ess sup}_{\omega\neq 0}\frac{1-\mathcal{F}[k_{s}](\omega)}{\left|\omega\right|^{\alpha}}\leq 1,\quad\forall\alpha\preceq\beta.$$

Using standard results from (Tsybakov, 2009), this implies that the resulting KDE $\hat{f}_{s}$ satisfies (5). Since $\psi=\mathcal{F}^{-1}[k_{s}]\in\mathcal{C}^{\infty}$, the Riemann–Lebesgue lemma guarantees that $\left|\kappa_{s}(x)\right|\leq c_{\beta,d}\left|x\right|^{\nu}$ is satisfied for $\nu=\lceil\beta+d\rceil$. Since $\psi$ is compactly supported, an application of Parseval’s identity yields $\kappa_{s}\in\mathcal{S}(\gamma,c_{\gamma})$. Applying Proposition 1 to $k_{s}$, we conclude that for the task of density estimation, weighted KDEs built on coresets are nearly as powerful as the coreset-based estimators studied in Section 2.

Theorem 2 shows that the Carathéodory coreset estimator achieves the minimax rate of estimation with near-optimal coreset size. In fact, a small modification yields a near-optimal rate of convergence for any coreset size as in Theorem 1.

Next we apply Proposition 1 to the popular Gaussian kernel $\phi(x)=(2\pi)^{-d/2}\exp(-\frac{1}{2}\left|x\right|_{2}^{2})$. This kernel has rapid decay in the real domain and Fourier space, and is thus amenable to our techniques. Moreover, $k_{\phi}$ is a kernel of order $\ell=1$, (Tsybakov, 2009, Definition 1.3 and Theorem 1.2) and so the standard KDE $\hat{f}_{\phi}$ on the full dataset attains the minimax rate of estimation $c_{d,L}n^{1/(2+d)}$ over the Lipschitz densities $\mathcal{P}_{\mathcal{H}}(1,L)$.

In addition, we have a nearly matching lower bound to Theorem 2 for coreset KDEs. In fact, our lower bound applies to a generalization of coreset KDEs where the vector of weights $(\lambda_{j})_{j}$ is not constrained to be in the simplex but can range within a hypercube of width that may grow polynomially with $n$: $\max_{j\in S}\left|\lambda_{j}\right|\leq n^{B}$.

This result is essentially a consequence of the lower bound in Theorem 1 because, in an appropriate sense, coreset KDEs with bounded weights are well-approximated by coreset-based estimators. Hence, in the case of bounded weights, allowing these weights to be measurable functions of the entire dataset rather than just the coreset, as would be required in Section 2, does not make a significant difference for the purpose of estimation. The full details of Theorem 4 are postponed to the Appendix.

Here we sketch the proof of Proposition 1, our main tool in constructing effective coreset KDEs. Full details of the argument may be found in the Appendix.

Let $k(x)=\prod_{i=1}^{d}\kappa(x_{i})$ denote a kernel, and suppose that $\hat{f}(y)=\frac{1}{n}\sum_{i=1}^{n}k_{h}(X_{i}-y)$ is a good estimator for an unknown density $f$ in that

$$\lVert f-\hat{f}\rVert_{2}\leq\varepsilon:=c_{\beta,d}n^{-\frac{\beta}{2\beta+d}}$$

on setting $h=n^{-1/(2\beta+d)}$. Our goal is to find a subset $S\subset[n]$ and weights $\{\lambda_{j}\}_{j\in S}$ such that

$$\frac{1}{n}\sum_{i=1}^{n}k_{h}(X_{i}-y)\approx\sum_{j\in S}\lambda_{j}k_{h}(X_{j}-y).$$

Suppose for simplicity that $\kappa$ is compactly supported on $[-1/2,1/2]$. By hypothesis and Parseval’s theorem $\kappa\in\mathcal{S}(\gamma,L^{\prime})$, and we can further show that $k\in\mathcal{S}(\gamma,c_{d,L^{\prime}})$ and $k_{h}\in\mathcal{S}(\gamma,c_{d,L^{\prime}}h^{-d/2-\gamma})$. Let $\mathcal{\bar{F}}$ denote the Fourier transform on the interval $[-1,1]$. Using the Fourier decay of $k_{h}$, we have

$$\lVert k_{h}(x)-\sum_{\left|\omega\right|_{\infty}<T}\mathcal{\bar{F}}[k_{h}](\omega)e^{i\langle x,\omega\rangle}\rVert_{2}\leq\varepsilon$$

(6)

when $T=(\frac{c_{d,\gamma,L^{\prime}}h^{-\frac{d}{2}-\gamma}}{\varepsilon})^{1/\gamma}=c_{d,\gamma,L^{\prime}}\,n^{\frac{d/2+\beta+\gamma}{\gamma(2\beta+d)}}$. Observe that this matches the setting of $T$ in Proposition 1.

The approximation (6) implies that for $X_{i}\in[-1/2,1/2]^{d}$,

$$\hat{f}(y)\approx\sum_{\left|\omega\right|_{\infty}<T}\mathcal{\bar{F}}[k_{h}](\omega)\left(\frac{1}{n}\sum_{i=1}^{n}e^{i\langle X_{i},\omega\rangle}\right)e^{i\langle y,\omega\rangle}.$$

Using the Carathéodory coreset and weights $\{\lambda_{j}\}_{j\in S}$ constructed in Section 3.1, it follows that

$$\sum_{\left|\omega\right|_{\infty}<T}\mathcal{\bar{F}}[k_{h}](\omega)\left(\frac{1}{n}\sum_{i=1}^{n}e^{i\langle X_{i},\omega\rangle}\right)e^{i\langle y,\omega\rangle}=\\ \sum_{\left|\omega\right|_{\infty}<T}\mathcal{\bar{F}}[k_{h}](\omega)\left(\sum_{i=1}^{n}\lambda_{j}e^{i\langle X_{i},\omega\rangle}\right)e^{i\langle y,\omega\rangle}.$$

Applying (6) again, we see that the right-hand-side is approximately equal to $\hat{g}_{S}(y)$, the estimator produced in Section (3.1). By the triangle inequality, we conclude that $\lVert\hat{g}_{S}(y)-f\rVert_{2}\leq c_{\beta,d}\,\varepsilon$, as desired.