Sampling Sketches for Concave Sublinear Functions of Frequencies

TL;DR

This paper introduces efficient, scalable sampling sketches for estimating statistics involving concave sublinear functions of key frequencies in massive distributed datasets, with strong theoretical guarantees and practical effectiveness.

Contribution

We develop composable sampling sketches tailored to any concave sublinear frequency function, achieving near-optimal size and high-quality statistical estimates.

Findings

Sketches are close in size to the target sample

Samples provide guarantees near those of ideal samples

Experimental results show simplicity and effectiveness

Abstract

We consider massive distributed datasets that consist of elements modeled as key-value pairs and the task of computing statistics or aggregates where the contribution of each key is weighted by a function of its frequency (sum of values of its elements). This fundamental problem has a wealth of applications in data analytics and machine learning, in particular, with concave sublinear functions of the frequencies that mitigate the disproportionate effect of keys with high frequency. The family of concave sublinear functions includes low frequency moments ($p \leq 1$), capping, logarithms, and their compositions. A common approach is to sample keys, ideally, proportionally to their contributions and estimate statistics from the sample. A simple but costly way to do this is by aggregating the data to produce a table of keys and their frequencies, apply our function to the frequency values, and then apply a weighted sampling scheme. Our main contribution is the design of composable sampling sketches that can be tailored to any concave sublinear function of the frequencies. Our sketch structure size is very close to the desired sample size and our samples provide statistical guarantees on the estimation quality that are very close to that of an ideal sample of the same size computed over aggregated data. Finally, we demonstrate experimentally the simplicity and effectiveness of our methods.

Tables

Table 1. Table 1 : Experimental Results: f ​ ( ν ) = ν 0.5 𝑓 𝜈 superscript 𝜈 0.5 f({\nu})={\nu}^{0.5} , 200 200 200 rep.

$k$	NRMSE		Benchmark		max #keys		max #elem
	bound	actual	ppswor	Pri.	ave	max	ave	max
Dataset: abcnews ($7.07\times {10}^6$ elements, $91.7\times {10}^3$ keys)
$25$	$0.834$	$0.213$	$0.213$	$0.217$	$31.7$	$37$	$50.9$	$76$
$50$	$0.577$	$0.142$	$0.128$	$0.137$	$58.5$	$66$	$95.1$	$136$
$75$	$0.468$	$0.120$	$0.111$	$0.110$	$85.4$	$94$	$134.8$	$181$
$100$	$0.404$	$0.105$	$0.098$	$0.103$	$111.2$	$120$	$171.1$	$256$
Dataset: flickr ($7.64\times {10}^6$ elements, $572.4\times {10}^3$ keys)
$25$	$0.834$	$0.200$	$0.190$	$0.208$	$31.2$	$37$	$53.1$	$77$
$50$	$0.577$	$0.144$	$0.147$	$0.142$	$57.8$	$64$	$94.6$	$130$
$75$	$0.468$	$0.123$	$0.114$	$0.110$	$83.7$	$91$	$131.7$	$175$
$100$	$0.404$	$0.115$	$0.095$	$0.099$	$108.9$	$116$	$173.4$	$223$
Dataset: zipf1.1 ($2.00\times {10}^6$ elements, $652.2\times {10}^3$ keys)
$25$	$0.834$	$0.215$	$0.198$	$0.217$	$31.8$	$39$	$52.5$	$75$
$50$	$0.577$	$0.123$	$0.137$	$0.131$	$58.7$	$66$	$95.0$	$130$
$75$	$0.468$	$0.109$	$0.115$	$0.114$	$84.7$	$91$	$135.2$	$186$
$100$	$0.404$	$0.106$	$0.103$	$0.097$	$111.2$	$119$	$176.3$	$221$
Dataset: zipf1.2 ($2.00\times {10}^6$ elements, $237.3\times {10}^3$ keys)
$25$	$0.834$	$0.199$	$0.208$	$0.214$	$31.1$	$38$	$53.2$	$83$
$50$	$0.577$	$0.144$	$0.138$	$0.145$	$57.9$	$65$	$98.4$	$139$
$75$	$0.468$	$0.122$	$0.116$	$0.124$	$83.9$	$90$	$138.2$	$173$
$100$	$0.404$	$0.098$	$0.109$	$0.096$	$109.6$	$115$	$179.2$	$227$
Dataset: zipf1.5 ($2.00\times {10}^6$ elements, $22.3\times {10}^3$ keys)
$25$	$0.834$	$0.201$	$0.207$	$0.194$	$30.1$	$35$	$53.4$	$74$
$50$	$0.577$	$0.152$	$0.139$	$0.142$	$56.1$	$60$	$101.5$	$136$
$75$	$0.468$	$0.115$	$0.115$	$0.112$	$81.6$	$86$	$151.8$	$199$
$100$	$0.404$	$0.098$	$0.094$	$0.086$	$107.1$	$113$	$196.3$	$248$

Table 2. Table 2 : Experimental Results: f ​ ( ν ) = ln ⁡ ( 1 + ν ) 𝑓 𝜈 1 𝜈 f({\nu})=\ln(1+{\nu}) , 200 200 200 rep.

$k$	NRMSE		Benchmark		max #keys		max #elem
	bound	actual	ppswor	Pri.	ave	max	ave	max
Dataset: abcnews ($7.07\times {10}^6$ elements, $91.7\times {10}^3$ keys)
$25$	$0.834$	$0.208$	$0.217$	$0.194$	$29.5$	$34$	$49.1$	$71$
$50$	$0.577$	$0.138$	$0.136$	$0.142$	$54.9$	$60$	$80.9$	$110$
$75$	$0.468$	$0.130$	$0.099$	$0.117$	$80.0$	$85$	$111.1$	$152$
$100$	$0.404$	$0.102$	$0.115$	$0.103$	$104.9$	$109$	$140.7$	$184$
Dataset: flickr ($7.64\times {10}^6$ elements, $572.4\times {10}^3$ keys)
$25$	$0.834$	$0.227$	$0.199$	$0.180$	$28.0$	$31$	$41.4$	$69$
$50$	$0.577$	$0.144$	$0.151$	$0.129$	$53.3$	$59$	$72.2$	$101$
$75$	$0.468$	$0.119$	$0.121$	$0.109$	$78.2$	$83$	$99.8$	$135$
$100$	$0.404$	$0.097$	$0.104$	$0.095$	$102.7$	$106$	$130.3$	$166$
Dataset: zipf1.1 ($2.00\times {10}^6$ elements, $652.2\times {10}^3$ keys)
$25$	$0.834$	$0.201$	$0.204$	$0.234$	$29.2$	$34$	$48.8$	$71$
$50$	$0.577$	$0.127$	$0.132$	$0.129$	$54.4$	$58$	$80.4$	$119$
$75$	$0.468$	$0.116$	$0.122$	$0.110$	$79.6$	$84$	$110.9$	$142$
$100$	$0.404$	$0.107$	$0.106$	$0.104$	$104.5$	$109$	$139.8$	$165$
Dataset: zipf1.2 ($2.00\times {10}^6$ elements, $237.3\times {10}^3$ keys)
$25$	$0.834$	$0.209$	$0.195$	$0.218$	$28.5$	$33$	$48.0$	$72$
$50$	$0.577$	$0.147$	$0.144$	$0.139$	$53.7$	$57$	$80.5$	$113$
$75$	$0.468$	$0.120$	$0.111$	$0.113$	$78.8$	$84$	$111.4$	$143$
$100$	$0.404$	$0.098$	$0.106$	$0.102$	$103.9$	$108$	$140.3$	$173$
Dataset: zipf1.5 ($2.00\times {10}^6$ elements, $22.3\times {10}^3$ keys)
$25$	$0.834$	$0.210$	$0.197$	$0.226$	$27.2$	$30$	$45.2$	$66$
$50$	$0.577$	$0.141$	$0.146$	$0.149$	$52.1$	$55$	$78.9$	$104$
$75$	$0.468$	$0.124$	$0.112$	$0.106$	$76.9$	$79$	$110.5$	$146$
$100$	$0.404$	$0.100$	$0.101$	$0.099$	$101.9$	$104$	$139.1$	$173$

Equations

$$\sum _{x\in \mathcal{X}}{L}_{x}f({\nu }_x)$$

$$\widehat{\sum _{x\in \mathcal{X}}{L}_x{f}_x}:=\sum _{x\in S}{L}_x\widehat{f_x}$$

$$\text{E}\left[\widehat{\sum _{x\in \mathcal{X}}{L}_x{f}_x}\right]=\text{E}\left[\sum _{x\in S}{L}_x\widehat{f_x}\right]$$

$$=\text{E}\left[\sum _{x\in \mathcal{X}}{L}_x\widehat{f_x}\right]=\sum _{x\in \mathcal{X}}{L}_x\text{E}\left[\widehat{f_x}\right]=\sum _{x\in \mathcal{X}}{L}_x{f}_x.$$

$$\widehat{f_x}=\{\begin{array}{ll}\frac{f_x}{{\mathrm{Pr}}_{\text{seed}(x)\sim {\text{SeedDist}}_x}[\text{seed}(x)<\tau ]}& x\in S\\ 0& x\notin S\end{array}.$$

$${\tau }_x:=\{\text{seed}(z)\mid z\in \mathcal{X}\setminus \{x\}{\}}_{(k-1)},$$

$$\text{Var}\left(\widehat{f_x}\right)\le \frac{1}{k-2}{f}_x\sum _{z\in \mathcal{X}}{f}_z.$$

$$\text{Var}\left(\widehat{f_x}\right)={\text{E}}_{{\tau }_x}\left[\text{Var}\left(\widehat{f_x}\mid {\tau }_x\right)\right]={\text{E}}_{{\tau }_x}\left[f_x^2\left(\frac{1}{\mathrm{Pr}[\text{seed}(x)<{\tau }_x]}-1\right)\right].$$

$$\text{Var}\left(\widehat{f_x}\mid {\tau }_x\right)=f_x^2\left(\frac{1}{\mathrm{Pr}[\text{seed}(x)<{\tau }_x]}-1\right)\le \frac{f_x}{{\tau }_x}.$$

$$\text{Var}\left[\widehat{f_x}\right]$$

$$\le {\text{E}}_{t\sim B}\left[f_x^2\left(\frac{1}{\mathrm{Pr}[\text{seed}(x)<t]}-1\right)\right].$$

$$f(\nu )={\mathcal{L}}^{\mathrm{c}}[a](\nu ):={\int }_0^{\infty }a(t)(1-e^{-\nu t})dt.$$

$$a(t)=({\mathcal{L}}^{\mathrm{c}}{)}^{-1}[f](t).$$

$${\mathcal{L}}^{\mathrm{c}}[a](\nu {)}_{\alpha }^{\beta }:={\int }_{\alpha }^{\beta }a(t)(1-e^{-\nu t})dt.$$

$$f(\nu )={\int }_0^{\infty }a(t)\min \{1,\nu t\}dt$$

$$(1-1/e)f(\nu )\le \tilde{f}(\nu )\le f(\nu ).$$

$$\text{Var}\left[{\hat{v}}_x\right]\le \frac{4v_{x}V}{k-2}.$$

$${\mathrm{SumMax}}_D(x):=\sum _{z\mid z.p=x}{\mathrm{Max}}_D(z)$$

$$\text{ElementScore}(e):=h(e.key)/e.val$$

$${\text{seed}}_D(x):=\underset{e\in D\mid e.key.p=x}{\min }\text{ElementScore}(e).$$

$$f({\nu }_x)={\mathcal{L}}^{\mathrm{c}}[a]({\nu }_x{)}_0^{\gamma }+{\mathcal{L}}^{\mathrm{c}}[a]({\nu }_x{)}_{\gamma }^{\infty },$$

$$\text{Var}\left[\widehat{f({\nu }_x)}\right]\le \frac{4f({\nu }_x){\sum }_{z\in \mathcal{X}}f({\nu }_z)}{(1-\epsilon {)}^2(k-2)}.$$

$$A(\gamma )$$

$$B(\gamma )$$

$$(1-\epsilon )\left({\int }_0^{\gamma }a(t)tdt\right){\nu }_x\le {\int }_0^{\gamma }a(t)(1-e^{-{\nu }_{x}t})dt\le \left({\int }_0^{\gamma }a(t)tdt\right){\nu }_x.$$

$${\int }_0^{\gamma }a(t)(1-e^{-{\nu }_{x}t})dt\le {\int }_0^{\gamma }a(t){\nu }_{x}tdt=\left({\int }_0^{\gamma }a(t)tdt\right){\nu }_x.$$

$${\int }_0^{\gamma }a(t)(1-e^{-{\nu }_{x}t})dt\ge {\int }_0^{\gamma }a(t)\left({\nu }_{x}t-\frac{({\nu }_{x}t{)}^2}{2}\right)dt.$$

$${\int }_0^{\gamma }a(t)(1-e^{-{\nu }_{x}t})dt\ge (1-\epsilon )\left({\int }_0^{\gamma }a(t)tdt\right){\nu }_x.$$

$$V_x=\frac{1}{r}{\mathrm{SumMax}}_E(x)+{\nu }_x{\int }_0^{\gamma }ta(t)dt.$$

$${\mathrm{Max}}_E((x,i))=\underset{e\in E\mid e.key=(x,i)}{\max }e.val=A(\max \{\gamma ,M_{x,i}\}).$$

$$\text{E}[V]\ge \frac{A(\gamma )}{r}\cdot k.$$

$$\text{E}[{\mathrm{Max}}_E((x,i))]={\mathcal{L}}^{\mathrm{c}}[a]({\nu }_x{)}_{\gamma }^{\infty }={\int }_{\gamma }^{\infty }a(t)(1-e^{-{\nu }_{x}t})dt.$$

Full text

1 Introduction

We consider massive distributed datasets that consist of elements that are key-value pairs $e=(e.\mathit{key},e.\mathit{val})$ with $e.\mathit{val}>0$. The elements are generated or stored on a large number of servers or devices. A key $x$ may repeat in multiple elements, and we define its frequency ${\nu}_{x}$ to be the sum of values of the elements with that key, i.e., ${\nu}_{x}:=\sum_{e\mid e.\mathit{key}=x}{e.\mathit{val}}$. For example, the keys can be search queries, videos, terms, users, or tuples of entities (such as video co-watches or term co-occurrences) and each data element can correspond to an occurrence or an interaction involving this key: the search query was issued, the video was watched, or two terms co-occurred in a typed sentence. An instructive common special case is when all elements have the same value $1$ and the frequency ${\nu}_{x}$ of each key $x$ in the dataset is simply the number of elements with key $x$.

A common task is to compute statistics or aggregates, which are sums over key contributions. The contribution of each key $x$ is weighted by a function of its frequency ${\nu}_{x}$. One example of such sum aggregates are queries of domain statistics $\sum_{x\in H}{{\nu}_{x}}$ for some domain (subset of keys) $H$. The domains of interest are often overlapping and specified at query time. Sum aggregates also arise as components of a larger pipeline, such as the training of a machine learning model with parameters $\boldsymbol{\theta}$, labeled examples $x\in\mathcal{X}$ with frequencies ${\nu}_{x}$, and a loss objective of the form $\ell(\mathcal{X};\boldsymbol{\theta})=\sum_{x}f({\nu}_{x})L(x;\boldsymbol{\theta})$. The function $f$ that is applied to the frequencies can be any concave sublinear function. Concave sublinear functions, which we discuss further below, are used in applications to mitigate the disproportionate effect of keys with very high frequencies. The training of the model typically involves repeated evaluation of the loss function (or of its gradient that also has a sum form) for different values of $\theta$. We would like to compute these aggregates on demand, without needing to go over the data many times.

When the number of keys is very large it is often helpful to compute a smaller random sample $S\subseteq\mathcal{X}$ of the keys from which aggregates can be efficiently estimated. In some applications, obtaining a sample can be the end goal. For example, when the aggregate is a gradient, we can use the sample itself as a stochastic gradient. To provide statistical guarantees on our estimate quality, the sampling needs to be weighted (importance sampling), with heavier keys sampled with higher probability, ideally, proportional to their contribution ($f({\nu}_{x})$). When the weights of the keys are known, there are classic sampling schemes that provide estimators with tight worst-case variance bounds [33, 17, 8, 13, 37, 38].

The datasets we consider here are presented in an unaggregated form: each key can appear multiple times in different locations. The focus of this work is designing composable sketch structures (formally defined below) that allow to compute a sample over unaggregated data with respect to the weights $f({\nu}_{x})$. One approach to compute a sample from unaggregated data is to first aggregate the data to produce a table of key-frequency pairs $(x,{\nu}_{x})$, compute the weights $f({\nu}_{x})$, and apply a weighted sampling scheme. This aggregation can be performed using composable structures that are essentially a table with an entry for each distinct key that occurred in the data. The number of distinct keys, however, and hence the size of that sketch, can be huge. For our sampling application, we would hope to use sketches of size that is proportional to the desired sample size, which is generally much smaller than the number of unique keys, and still provide statistical guarantees on the estimate quality that are close to that of a weighted sample computed according to $f({\nu}_{x})$.

Concave Sublinear Functions.

Typical datasets have a skewed frequency distribution, where a small fraction of the keys have very large frequencies and we can get better results or learn a better model of the data by suppressing their effect. The practice is to apply a concave sublinear function $f$ to the frequency, so that the importance weight of the key is $f({\nu}_{x})$ instead of simply its frequency ${\nu}_{x}$. This family of functions includes the frequency moments ${\nu}_{x}^{p}$ for $p\leq 1$, $\ln(1+{\nu}_{x})$, $\textsf{cap}_{T}({\nu}_{x})=\min\{T,{\nu}_{x}\}$ for a fixed $T\geq 0$, their compositions, and more. A formal definition appears in Section 2.4.

Two hugely popular methods for producing word embeddings from word co-occurrences use this form of mitigation: word2vec [30] uses $f({\nu})={\nu}^{0.5}$ and $f({\nu})={\nu}^{0.75}$ for positive and negative examples, respectively, and GloVe [35] uses $f({\nu})=\min\{T,{\nu}^{0.75}\}$ to mitigate co-occurrence frequencies. When the data is highly distributed, for example, when it originates or resides at millions of mobile devices (as in federated learning [29]), it is useful to estimate the loss or compute a stochastic gradient update efficiently via a weighted sample.

The suppression of higher frequencies may also directly arise in applications. One example is campaign planning for online advertising, where the value of showing an ad to a user diminishes with the number of views. Platforms allow an advertiser to specify a cap value $T$ on the number of times the same ad can be presented to a user [23, 34]. In this case, the number of opportunities to display an ad to a user $x$ is a cap function of the frequency of the user $f({\nu}_{x})=\min\{T,{\nu}_{x}\}$, and the number for a segment of users $H$ is the statistics $\sum_{x\in H}f({\nu}_{x})$. When planning a campaign, we need to quickly estimate the statistics for different segments, and this can be done from a sample that ideally is weighted by $f({\nu}_{x})$.

Our Contribution.

In this work, we design composable sketches that can be tailored to any concave sublinear function $f$, and allow us to compute a weighted sample over unaggregated data with respect to the weights $f({\nu}_{x})$. Using the sample, we will be able to compute unbiased estimators for the aggregates mentioned above. In order to compute the estimators, we need to make a second pass over the data: In the first pass, we compute the set of sampled keys, and in the second pass we compute their frequencies. Both passes can be done in a distributed manner.

A sketch $S(D)$ is a data structure that summarizes a set $D$ of data elements, so that the output of interest for $D$ (in our case, a sample of keys) can be recovered from the sketch $S(D)$. A sketch structure is composable if we can obtain a sketch $S(D_{1}\cup D_{2})$ of two sets of elements $D_{1}$ and $D_{2}$ from the sketches $S(D_{1})$ and $S(D_{2})$ of the sets. This property alone gives us full flexibility to parallelize or distribute the computation. The size of the sketch determines the communication and storage needs of the computation.

We provide theoretical guarantees on the quality (variance) of the estimators. The baseline for our analysis is the bounds on the variance that are guaranteed by PPSWOR on aggregated data. PPSWOR [37, 38] is a sampling scheme with tight worst-case variance bounds. The estimators provided by our sketch have variance at most $4/((1-\varepsilon)^{2})$ times the variance bound for PPSWOR. The parameter $\varepsilon\leq 1/2$ mostly affects the run time of processing a data element, which grows near-linearly in $1/\varepsilon$. Thus, our sketch allows us to get approximately optimal guarantees on the variance while avoiding the costly aggregation of the data.

We remark that these guarantees are for soft concave sublinear functions. This family approximates any concave sublinear function up to a multiplicative factor of $1-1/e$. As a result, our sketch can be used with any (non-soft) concave sublinear function while incurring another factor of $\left(1+\frac{1}{e-1}\right)^{2}$ in the variance.

The space required by our sketch significantly improves upon the previous methods (which all require aggregating the data). In particular, if the desired sample size is $k$, we show that the space required by the sketch at any given time is $O(k)$ in expectation. We additionally show that, with probability at least $1-\delta$, the space will not exceed $O\left(k+\min\{\log{m},\log{\log{\left(\frac{\operatorname{\textsf{Sum}}_{D}}{\operatorname{\textsf{Min}}(D)}\right)}}\}+\log{\left(\frac{1}{\delta}\right)}\right)$ at any time while processing the dataset $D$, where $m$ is the number of elements, $\operatorname{\textsf{Sum}}_{D}$ the sum of weights of all elements, and $\operatorname{\textsf{Min}}(D)$ is the minimum value of an element in $D$. In the common case where all elements have weight $1$, this means that for any $\delta$, the needed space is at most $O\left(k+\log{\log{m}}+\log{\left(\frac{1}{\delta}\right)}\right)$ with probability at least $1-\delta$.111For streaming algorithms we are typically interested in deterministic worst-case bounds on the space, but streaming algorithms with randomized space have also been considered in some cases, in particular when studying the sliding window model [3, 7].

We complement our work with a small-scale experimental study. We use a simple implementation of our sampling sketch to study the actual performance in terms of estimate quality and sketch size. In particular, we show that the estimate quality is even better than the (already adequate) guarantees provided by our worst-case bounds. We additionally compare the estimate quality to that of two popular sampling schemes for aggregated data, PPSWOR [37, 38] and priority (sequential Poisson) sampling [33, 17]. In the experiments, we see that the estimate quality of our sketch is close to what achieved by PPSWOR and priority sampling, while our sketch uses much less space by eliminating the need for aggregation.

The paper is organized as follows. The preliminaries are presented in Section 2. We provide an overview of PPSWOR and the statistical guarantees it provides for estimation. Then, we formally define the family of concave sublinear functions. Our sketch uses two building blocks. The first building block, which can be of independent interest, is the analysis of a stochastic PPSWOR sample. Typically, when computing a sample, the data from which we sample is a deterministic part of the input. In our construction, we needed to analyze the variance bounds for PPSWOR sampling that is computed over data elements with randomized weights (under certain assumptions). We provide this analysis in Section 3. The second building block is the $\operatorname{\textsf{SumMax}}$ sampling sketch, which is discussed in Section 4. This is an auxiliary sketch structure that supports datasets with a certain type of structured keys. We put it all together and describe our main result in Section 5. The experiments are discussed in Section 6.

Related Work.

There are multiple classic composable weighted sampling schemes for aggregated datasets (where keys are unique to elements). Schemes that provide estimators with tight worst-case variance bounds include priority (sequential Poisson) sampling [33, 17] and VarOpt sampling [8, 13]. We focus here on PPSWOR [37, 38] as our base scheme because it extends to unaggregated datasets, where multiple elements can additively contribute to the frequency/weight of each key.

There is a highly prolific line of research on developing sketch structures for different task over streamed or distributed unaggregated data with applications in multiple domains. Some early examples are frequent elements [31] and distinct counting [20], and the seminal work of [1] providing a theoretical model for frequency moments. Composable sampling sketches for unaggregated datasets were also studied for decades. The goal is to meet the quality of samples computed on aggregated frequencies while using a sketch structure that can only hold a final-sample-size number of distinct keys. These include a folklore sketch for distinct sampling ($f({\nu})=1$ when ${\nu}>0$) [27, 39] and sketch structures for sum sampling ($f({\nu})={\nu}$) [12]. The latter generalizes the discrete sample and hold scheme [22, 18, 14] and PPSWOR. Sampling sketches for cap functions ($f({\nu})=\min\{T,{\nu}\}$) were provided in [11] and have a slight overhead over the aggregated baseline. The latter work also provided multi-objective/universal samples that with a logarithmic overhead simultaneously provide statistical guarantees for all concave sublinear $f$. In the current work we propose sampling sketches that can be tailored to any concave sublinear function and only have a small constant overhead.

An important line of work uses sketches based on random linear projections to estimate frequency statistics and to sample. In particular, $\ell_{p}$ sampling sketches [21, 32, 2, 26, 25] sample (roughly) according to $f({\nu})={\nu}^{p}$ for $p\in[0,2]$. These sketches have a higher logarithmic overhead on the space compared to sample-based sketches, and do not support all concave sublinear functions of the frequencies (for example, $f({\nu})=\ln{(1+{\nu})}$). In some respects they are more limited in their application – for example, they are not designed to produce a sample that includes raw keys. Their advantage is that they can be used with super-linear ($p\in(1,2]$) functions of frequencies and can also support signed element values (the turnstile model). For the more basic problem of sketches that estimate frequency statistics over the full data, a complete characterization of the frequency functions for which the statistics can be estimated via polylogarithmic-size sketches is provided in [6, 4]. Universal sketches for estimating $\ell_{p}$ norms of subsets were recently considered in [5]. The seminal work of Alon et al. [1] established that for some functions of frequencies (moments with $p>2$), statistics estimation requires polynomial-size sketches. A double logarithmic size sketch, extending [19] for distinct counting, that computes statistics over the entire dataset for all soft concave sublinear functions is provided in [10]. Our design builds on components of that sketch.

2 Preliminaries

Consider a set $D$ of data elements of the form $e=(e.\mathit{key},e.\mathit{val})$ where $e.\mathit{val}>0$. We denote the set of possible keys by $\mathcal{X}$. For a key $z\in\mathcal{X}$, we let $\operatorname{\textsf{Max}}_{D}(z):=\max_{e\in D\mid e.\mathit{key}=z}e.\mathit{val}$ and $\operatorname{\textsf{Sum}}_{D}(z):=\sum_{e\in D\mid e.\mathit{key}=z}e.\mathit{val}$ denote the maximum value of a data element in $D$ with key $z$ and the sum of values of data elements in $D$ with key $z$, respectively. Each key $z\in\mathcal{X}$ that appears in $D$ is called active. If there is no element $e\in D$ with $e.\mathit{key}=z$, we say that $z$ is inactive and define $\operatorname{\textsf{Max}}_{D}(z):=0$ and $\operatorname{\textsf{Sum}}_{D}(z):=0$. When $D$ is clear from context, it is omitted. For a key $z$, we use the shorthand ${\nu}_{z}:=\operatorname{\textsf{Sum}}_{D}(z)$ and refer to it as the frequency of $z$. The sum and the max-distinct statistics of $D$ are defined, respectively, as $\operatorname{\textsf{Sum}}_{D}:=\sum_{e\in D}e.\mathit{val}$ and $\operatorname{\textsf{MxDistinct}}_{D}:=\sum_{z\in\mathcal{X}}\operatorname{\textsf{Max}}_{D}(z)$. For a function $f$, $f_{D}:=\sum_{z\in\mathcal{X}}f(\operatorname{\textsf{Sum}}_{D}(z))=\sum_{z\in\mathcal{X}}f({\nu}_{z})$ is the $f$-frequency statistics of $D$.

For a set $A\subseteq\mathbb{R}$, we use $A_{(i)}$ to denote the $i$-th order statistic of $A$, that is, the $i$-th lowest element in $A$.

2.1 The Composable Bottom-$k$ Structure

In this work, we will use composable sketch structures in order to efficiently summarize streamed or distributed data elements. A composable sketch structure is specified by three operations: The initialization of an empty sketch structure $s$, the processing of a data element $e$ into a structure $s$, and the merging of two sketch structures $s_{1}$ and $s_{2}$. To sketch a stream of elements, we start with an empty structure and sequentially process data elements while storing only the sketch structure. The merge operation is useful with distributed or parallel computation and allows us to compute the sketch of a large set $D=\bigcup_{i}D_{i}$ of data elements by merging the sketches of the parts $D_{i}$.

In particular, one of the main building blocks that we use is the bottom-$k$ structure [15], specified in Algorithm 1. The structure maintains $k$ data elements: For each key, consider only the element with that key that has the minimum value. Of these elements, the structure keeps the $k$ elements that have the lowest values.

2.2 The PPSWOR Sampling Sketch

In this subsection, we describe a scheme to produce a sample of $k$ keys, where at each step the probability that a key is selected is proportional to its weight. That is, the sample we produce will be equivalent to performing the following $k$ steps. At each step we select one key and add it to the sample. At the first step, each key $x\in\mathcal{X}$ (with weight $w_{x}$) is selected with probability $w_{x}/\sum_{y}w_{y}$. At each subsequent step, we choose one of the remaining keys, again with probability proportional to its weight. Since the total weight of the remaining keys is lower, the probability that a key is selected in a subsequent step (provided it was not selected earlier) is only higher. This process is called probability proportional to size and without replacement (PPSWOR) sampling.

A classic method for PPSWOR sampling is the following scheme [37, 38]. For each key $x$ with weight $w_{x}$, we independently draw $\textnormal{{seed}}(x)\sim\textsf{Exp}(w_{x})$. Outputting the sample that includes the $k$ keys with smallest $\textnormal{{seed}}(x)$ is equivalent to PPSWOR sampling as described above. This method together with a bottom-$k$ structure can be used to implement PPSWOR sampling over a set of data elements $D$ according to ${\nu}_{x}=\operatorname{\textsf{Sum}}_{D}(x)$. This sampling sketch is presented here as Algorithm 2. The sketch is due to [12] (based on [22, 18, 14]).

Proposition 2.1**.**

Algorithm 2 maintains a composable bottom-$k$ structure such that for each key $x$, the lowest value of an element with key $x$ (denoted by $\textnormal{{seed}}(x)$) is drawn independently from $\textsf{Exp}({\nu}_{x})$. Hence, it is a PPSWOR sample according to the weights ${\nu}_{x}$.

The proof is provided in Appendix A.

2.3 Estimation Using Bottom-$k$ Samples

PPSWOR sampling (Algorithm 2) is a special case of bottom-$k$ sampling [38, 15, 16].

Definition 2.2**.**

Let $k\geq 2$. A bottom-$k$ sample over keys $\mathcal{X}$ is obtained by drawing independently for each active key $x$ a random variable $\textnormal{{seed}}(x)\sim\textsf{SeedDist}_{x}$. The $k-1$ keys with lowest $\textnormal{{seed}}(x)$ values are considered to be included in the sample $S$, and the $k$-th lowest value $\tau:=\{\textnormal{{seed}}(x)\mid x\in\mathcal{X}\}_{(k)}$ is the inclusion threshold.

The distributions $\textsf{SeedDist}_{x}$ are such that for all $t>0$, $\Pr[\textnormal{{seed}}(x)<t]>0$, that is, there is positive probability to be below any positive $t$. Typically the distributions $\textsf{SeedDist}_{x}$ come from a family of distributions that is parameterized by the frequency ${\nu}_{x}$ of keys. The frequency is positive for active keys and [math] otherwise. In the special case of PPSWOR sampling by frequency, $\textsf{SeedDist}_{x}$ is $\textsf{Exp}({\nu}_{x})$.

We review here how a bottom-$k$ sample is used to estimate domain statistics of the form $\sum_{x\in H}f({\nu}_{x})$ for $H\subseteq\mathcal{X}$. More generally, we will show how to estimate aggregates of the form

[TABLE]

for any set of fixed values $L_{x}$. Note that we can represent domain statistics in this form by setting $L_{x}=1$ for $x\in H$ and $L_{x}=0$ for $x\notin H$. For the sake of this discussion, we treat $f_{x}:=f({\nu}_{x})$ as simply a set of weights associated with keys, assuming we can have $f_{x}>0$ only for active keys.

In order to estimate statistics of the form (1), we will define an estimator $\widehat{f_{x}}$ for each $f_{x}$. The estimator $\widehat{f_{x}}$ will be non-negative ($\widehat{f_{x}}\geq 0$), unbiased ($\textsf{E}\left[\widehat{f_{x}}\right]=f_{x}$), and such that $\widehat{f_{x}}=0$ when the key $x$ is not included in the bottom-$k$ sample ($x\not\in S$ using the terms of Definition 2.2).

As a general convention, we will use the notation $\widehat{z}$ to denote an estimator of any quantity $z$. We define the sum estimator of the statistics $\sum_{x\in\mathcal{X}}{L_{x}f_{x}}$ to be

[TABLE]

We also note that since $\widehat{f_{x}}=0$ for $x\notin S$, computing the sum only over $x\in S$ is the same as computing the sum over all $x\in\mathcal{X}$, that is, $\sum_{x\in S}{L_{x}\widehat{f_{x}}}=\sum_{x\in\mathcal{X}}{L_{x}\widehat{f_{x}}}$.

Note that the sum estimator can be computed as long as the fixed values $L_{x}$ and the per-key estimates $\widehat{f_{x}}$ for $x\in S$ are available. From linearity of expectation, we get that the sum estimate is unbiased:

[TABLE]

We now define the per-key estimators $\widehat{f_{x}}$. The following is a conditioned variant of the Horvitz-Thompson estimator [24].

Definition 2.3**.**

Let $k\geq 2$ and consider a bottom-$k$ sample, where $S$ is the set of $k-1$ keys in the sample and $\tau$ is the inclusion threshold. For any $x\in\mathcal{X}$, the inverse-probability estimator of $f_{x}$ is

[TABLE]

In order to compute these estimates, we need to know the weights $f_{x}$ and distributions $\textsf{SeedDist}_{x}$ for the sampled keys $x\in S$. In particular, in our applications when $f_{x}=f({\nu}_{x})$ is a function of frequency and the seed distribution is parameterized by frequency, then it suffices to know the frequencies of sampled keys.222In the general case, we assume these functions of the frequency are computationally tractable or can be easily approximated up to a small constant factor. For our application, the discussion will follow in Section 5.4.

Claim 2.4**.**

The inverse-probability estimator is unbiased, that is, $\textsf{E}\left[\widehat{f_{x}}\right]=f_{x}$.

Proof.

We first consider $\widehat{f_{x}}$ when conditioned on the seed values of all other keys $\mathcal{X}\setminus\{x\}$ and in particular on

[TABLE]

which is the $k-1$ smallest seed on $\mathcal{X}\setminus\{x\}$. Under this conditioning, a key $x$ is included in $S$ with probability $\Pr_{\textnormal{{seed}}(x)\sim\textsf{SeedDist}_{x}}[\textnormal{{seed}}(x)<\tau_{x}]$. When $x\not\in S$, the estimate is [math]. When $x\in S$, we have that $\tau_{x}=\tau$ and the estimate is the ratio of $f_{x}$ and the inclusion probability. So our estimator is a plain inverse probability estimator and thus $\textsf{E}\left[\widehat{f_{x}}\mid\tau_{x}\right]=f_{x}$.

Finally, from the fact that the estimator is unbiased when conditioned on $\tau_{x}$, we also get that it is unconditionally unbiased: $\textsf{E}\left[\widehat{f_{x}}\right]=\textsf{E}_{\tau_{x}}\left[\textsf{E}\left[\widehat{f_{x}}\mid\tau_{x}\right]\right]=\textsf{E}_{\tau_{x}}\left[f_{x}\right]=f_{x}$. ∎

We now turn to analyze the variance of the estimators. The guarantees we can obtain on the quality of the sum estimates depend on how well the distributions $\textsf{SeedDist}_{x}$ are tailored to the values $f_{x}$, where ideally, keys should be sampled with probabilities proportional to $f_{x}$. PPSWOR, where $\textnormal{{seed}}(x)\sim\textsf{Exp}(f_{x})$, is such a “gold standard” sampling scheme that provides us with strong guarantees: For domain statistics $\sum_{x\in H}f_{x}$, we get a tight worst-case bound on the coefficient of variation333Defined as the ratio of the standard deviation to mean. For our unbiased estimators it is equal to the relative root mean squared error. of $1/\sqrt{q(k-2)}$, where $q=\sum_{x\in H}f_{x}/\sum_{x\in\mathcal{X}}f_{x}$ is the fraction of the statistics that is due to the domain $H$. Moreover, the estimates are concentrated in a Chernoff bounds sense. For objectives of the form (1), we obtain additive Hoeffding-style bounds that depend only on sample size and the range of $L_{x}$.

When we cannot implement “gold-standard” sampling via small composable sampling structures, we seek guarantees that are close to that. Conveniently, in the analysis it suffices to bound the variance of the per-key estimators [9, 11]: A key property of bottom-$k$ estimators is that $\forall x,z,\ \textsf{cov}(\widehat{f_{x}},\widehat{f_{z}})\leq 0$ (equality holds for $k\geq 3$) [11]. Therefore, the variance of the sum estimator can be bounded by the sum of bounds on the per-key variance. This allows us to only analyze the per-key variance $\textsf{Var}\left(\widehat{f_{x}}\right)$. To achieve the guarantees of the “gold standard” sampling, the desired bound on the per-key variance for a sample of size $k-1$ (a bottom-$k$ sample where the $k$-th lowest seed is the inclusion threshold) is

[TABLE]

So our goal is to establish upper bounds on the per-key variance that are within a small constant of (2). We refer to this value as the overhead. The overhead factor in the per-key bounds carries over to the sum estimates.

We next review the methodology for deriving per-key variance bounds. The starting point is to first bound the per-key variance of $\widehat{f_{x}}$ conditioned on $\tau_{x}$.

Claim 2.5**.**

With the inverse probability estimator we have

[TABLE]

Proof.

Follows from the law of total variance and the unbiasedness of the conditional estimates for any fixed value of $\tau_{x}$, $\textsf{E}\left[\widehat{f_{x}}\mid\tau_{x}\right]=f_{x}$. ∎

For the “gold standard” PPSWOR sample, we have $\Pr[\textnormal{{seed}}(x)<t]=1-\exp(-f_{x}t)$ and using $\frac{e^{-x}}{1-e^{-x}}\leq\frac{1}{x}$, we get

[TABLE]

In order to bound the unconditional per-key variance, we use the following notion of stochastic dominance.

Definition 2.6**.**

Consider two density functions $a$ and $b$ both with support on the nonnegative reals. We say that $a$ is dominated by $b$ ($a\preceq b$) if for all $z\geq 0$, $\int_{0}^{z}a(y)dy\leq\int_{0}^{z}b(y)dy$.

Note that since they are both density functions, it implies that the CDF of $a$ is pointwise at most the CDF of $b$. In particular, the probability of being below some value $y$ under $b$ is at least that of $a$.444We note in our applications, lower values mean higher inclusion probabilities. In most applications, higher values are associated with better results, and accordingly, first-order stochastic dominance is usually defined as the reverse.

When bounding the variance, we use a distribution $B$ that dominates the distribution of $\tau_{x}$ and is easier to work with and then compute the upper bound

[TABLE]

With PPSWOR, the distribution of $\tau_{x}$ is dominated by the distribution $\textsf{Erlang}[\sum_{z\in\mathcal{X}}f_{z},k-1]$, where $\textsf{Erlang}[V,k]$ is the distribution of the sum of $k$ independent exponential random variables with parameter $V$. The density function of $\textsf{Erlang}[V,k]$ is $B_{V,k}(t)=\frac{V^{k}t^{k-1}e^{-Vt}}{(k-1)!}$. Choosing $B$ to be $\textsf{Erlang}[\sum_{z\in\mathcal{X}}f({\nu}_{z}),k-1]$ in (4) and using (3), we get the bound in (2).

Note that if we have an estimator that gives a weaker bound of $c\cdot\frac{f_{x}}{\tau_{x}}$ on the conditional variance and the distribution of $\tau_{x}$ is similarly dominated by $\textsf{Erlang}[\sum_{z\in\mathcal{X}}f_{z},k-1]$, we will obtain a corresponding bound on the unconditional variance with overhead $c$.

2.4 Concave Sublinear Functions

A function $f:[0,\infty)\rightarrow[0,\infty)$ is soft concave sublinear if for some $a(t)\geq 0$ it can be expressed as555The definition also allows $a(t)$ to have discrete mass at points (that is, we can add a component of the form $\sum_{i}{a(t_{i})(1-e^{-{\nu}t_{i}})}$). We generally ignore this component for the sake of presentation, but one way to model this is using Dirac delta.

[TABLE]

$\operatorname{\mathcal{L}^{c}}[a]({\nu})$ is called the complement Laplace transform of $a$ at ${\nu}$. The function $a(t)$ is the inverse Laplacec (complement Laplace) transform of $f$:

[TABLE]

A table with the inverse Laplacec transform of several common functions (in particular, the moments ${\nu}^{p}$ for $p\in(0,1)$ and $\ln{(1+{\nu})}$) appears in [10]. We additionally use the notation

[TABLE]

The sampling schemes we present in this work will be defined for soft concave sublinear functions of the frequencies. However, this will allow us to estimate well any function that is within a small multiplicative constant of a soft concave sublinear function. In particular, we can estimate concave sublinear functions. These functions can be expressed as

[TABLE]

for $a(t)\geq 0$.666Here we also allow $a(t)$ to have discrete mass using Dirac delta. For the sake of presentation, we also assume bounded ${\nu}$ – otherwise we need to add a linear component $A_{\infty}{\nu}$ for some $A_{\infty}\geq 0$. The component $A_{\infty}{\nu}$ can easily be added to the final sketch presented in Section 5, for example, by taking the minimum with another independent PPSWOR sketch. The concave sublinear family includes all functions such that $f(0)=0$, $f$ is monotonically non-decreasing, $\partial_{+}f(0)<\infty$, and $\partial^{2}f\leq 0$.

Any concave sublinear function $f$ can be approximated by a soft concave sublinear function as follows. Consider the corresponding soft concave sublinear function $\tilde{f}$ using the same coefficients $a(t)$. The function $\tilde{f}$ closely approximates $f$ pointwise [10]:

[TABLE]

Our weighted sample for $\tilde{f}$ will respectively approximate a weighted sample for $f$ (later explained in Remark 5.11).

3 Stochastic PPSWOR Sampling

In this section, we provide an analysis of PPSWOR for a case that will appear later in our main sketch. The case we consider is the following. In the PPSWOR sampling scheme described in Section 2.2, the weights $w_{x}$ of the keys were part of the deterministic input to the algorithm. In this section, we consider PPSWOR sampling when the weights are random variables. We will show that under certain assumptions, PPSWOR sampling according to randomized inputs is close to sampling according to the expected values of these random inputs.

Formally, let $\mathcal{X}$ be a set of keys. Each key $x\in\mathcal{X}$ is associated with $r_{x}\geq 0$ independent random variables $S_{x,1},\ldots,S_{x,r_{x}}$ in the range $[0,T]$ (for some constant $T>0$). The weight of key $x$ is the random variable $S_{x}:=\sum_{i=1}^{r_{x}}{S_{x,i}}$. We additionally denote its expected weight by $v_{x}:=\textsf{E}[S_{x}]$, and the expected sum statistics by $V:=\sum_{x}v_{x}$.

A stochastic PPSWOR sample is a PPSWOR sample computed for the key-value pairs $(x,S_{x})$. That is, we draw the random variables $S_{x}$, then we draw for each $x$ a random variable $\textnormal{{seed}}(x)\sim\textsf{Exp}[S_{x}]$, and take the $k$ keys with lowest seed values.

We prove two results that relate stochastic PPSWOR sampling to a PPSWOR sample according to the expected values $v_{x}$. The first result bounds the variance of estimating $v_{x}$ using a stochastic PPSWOR sample. We consider the conditional inverse-probability estimator of $v_{x}$ (Definition 2.3). Note that even though the PPSWOR sample was computed using the random weight $S_{x}$, the estimator $\widehat{v_{x}}$ is computed using $v_{x}$ and will be $\frac{v_{x}}{\Pr[\textnormal{{seed}}(x)<\tau]}$ for keys $x$ in the sample. Based on the discussion in Section 2.3, it suffices to bound the per-key variance and relate it to the per-key variance bound for a PPSWOR sample computed directly for $v_{x}$. We show that when $V\geq Tk$, the overhead due to the stochastic sample is at most $4$ (that is, the variance grows by a multiplicative factor of $4$). The proof details would also reveal that when $V\gg Tk$, the worst-case bound on the overhead is actually closer to $2$.

Theorem 3.1**.**

Let $k\geq 3$. In a stochastic PPSWOR sample, if $V\geq Tk$, then for every key $x\in\mathcal{X}$, the variance $\textsf{Var}\left[\hat{v}_{x}\right]$ of the bottom-$k$ inverse probability estimator of $v_{x}$ is bounded by

[TABLE]

Note that in order to compute these estimates, we need to be able to compute the values $v_{x}=\textsf{E}[S_{x}]$ and $\Pr[\textnormal{{seed}}(x)<\tau]$ for sampled keys. With stochastic sampling, the precise distribution $\textsf{SeedDist}_{x}$ depends on the distributions of the random variables $S_{x,i}$. For now, however, we assume that $\textsf{SeedDist}_{x}$ and $v_{x}$ are available to us with the sample. In Section 5, when we use stochastic sampling, we will also show how to compute $\textsf{SeedDist}_{x}$.

The second result in this section provides a lower bound on the probability that a key $x$ is included in the stochastic PPSWOR sample of size $k=1$. We show that when $V\geq\frac{1}{\varepsilon}\ln{\left(\frac{1}{\varepsilon}\right)}T$, the probability key $x$ is included in the sample is at least $1-2\varepsilon$ times the probability it is included in a PPSWOR sample according to the expected weights.

Theorem 3.2**.**

Let $\varepsilon\leq\frac{1}{2}$. Consider a stochastic PPSWOR sample of size $k=1$. If $V\geq\frac{1}{\varepsilon}\ln{\left(\frac{1}{\varepsilon}\right)}T$, the probability that any key $x\in\mathcal{X}$ is included in the sample is at least $(1-2\varepsilon)\frac{v}{V}$.

The proofs of the two theorems are deferred to Appendix B.

4 $\operatorname{\textsf{SumMax}}$ Sampling

Sketch

In this section, we present an auxiliary sampling sketch which will be used in Section 5. The sketch processes elements $e=(e.\mathit{key},e.\mathit{val})$ with keys $e.\mathit{key}=(e.\mathit{key}.p,e.\mathit{key}.s)$ that are structured to have a primary key $e.\mathit{key}.p$ and a secondary key $e.\mathit{key}.s$. For each primary key $x$, we define

[TABLE]

where $\operatorname{\textsf{Max}}$ is as defined in Section 2. If there are no elements $e\in D$ such that $e.\mathit{key}.p=x$, then by definition $\operatorname{\textsf{Max}}_{D}(z)=0$ for all $z$ with $z.p=x$ (as there are no elements in $D$ with key $z$) and therefore $\operatorname{\textsf{SumMax}}_{D}(x)=0$. Our goal in this section is to design a sketch that produces a PPSWOR sample of primary keys $x$ according to weights $\operatorname{\textsf{SumMax}}_{D}(x)$. Note that while the key space of the input elements contains structured keys of the form $e.\mathit{key}=(e.\mathit{key}.p,e.\mathit{key}.s)$, the key space for the output sample will be the space of primary keys only. Our sampling sketch is described in Algorithm 3.

The sketch structure consists of a bottom-$k$ structure and a hash function $h$. We assume we have a perfectly random hash function $h$ such that for every key $z=(z.p,z.s)$, $h(z)\sim\textsf{Exp}[1]$ independently (in practice, we assume that the hash function is provided by the platform on which we run). We process an input element $e$ by generating a new data element with key $e.\mathit{key}.p$ (the primary key of the key of the input element) and value

[TABLE]

and then processing that element by our bottom-$k$ structure. The bottom-$k$ structure holds our current sample of primary keys.

By definition, the bottom-$k$ structure retains the $k$ primary keys $x$ with minimum

[TABLE]

To establish that this is a PPSWOR sample according to $\operatorname{\textsf{SumMax}}_{D}(x)$, we study the distribution of $\textnormal{{seed}}_{D}(x)$.

Lemma 4.1**.**

For all primary keys $x$ that appear in elements of $D$, $\textnormal{{seed}}_{D}(x)\sim\textsf{Exp}[\operatorname{\textsf{SumMax}}_{D}(x))]$. The random variables $\textnormal{{seed}}_{D}(x)$ are independent.

The proof is deferred to Appendix C.

Note that the distribution of $\textnormal{{seed}}_{D}(x)$, which is $\textsf{Exp}[\operatorname{\textsf{SumMax}}_{D}(x)]$, does not depend on the particular structure of $D$ or the order in which elements are processed, but only on the parameter $\operatorname{\textsf{SumMax}}_{D}(x)$. The bottom-$k$ sketch structure maintains the $k$ primary keys with smallest $\textnormal{{seed}}_{D}(x)$ values. We therefore get the following corollary.

Corollary 4.2**.**

Given a stream or distributed set of elements $D$, the sampling sketch in Algorithm 3 produces a PPSWOR sample according to the weights $\operatorname{\textsf{SumMax}}_{D}(x)$.

5 Sampling Sketch for Functions of Frequencies

In this section, we are given a set $D$ of elements $e=(e.\mathit{key},e.\mathit{val})$ and we wish to maintain a sample of $k$ keys, that will be close to PPSWOR according to a soft concave sublinear function of their frequencies $f({\nu}_{x})$. At a high level, our sampling sketch is guided by the sketch for estimating the statistics $f_{D}$ due to Cohen [10]. Our sketch uses a parameter $\varepsilon$ that will tradeoff the running time of processing an element with the bound on the variance of the inverse-probability estimator.

Recall that a soft concave sublinear function $f$ can be represented as $f(w)=\operatorname{\mathcal{L}^{c}}[a](w)_{0}^{\infty}=\int_{0}^{\infty}a(t)(1-e^{-wt})dt$ for $a(t)\geq 0$. Using this representation, we express $f({\nu}_{x})$ as a sum of two contributions for each key $x$:

[TABLE]

where $\gamma$ is a value we will set adaptively while processing the elements. Our sampling sketch is described in Algorithm 4. It maintains a separate sampling sketch for each set of contributions. The sketch for $\operatorname{\mathcal{L}^{c}}[a]({\nu}_{x})_{0}^{\gamma}$ is discussed in Section 5.1, and the sketch for $\operatorname{\mathcal{L}^{c}}[a]({\nu}_{x})_{\gamma}^{\infty}$ is discussed in Section 5.2. In order to produce a sample from the sketch, these separate sketches need to be combined. Algorithm 5 describes how to produce a final sample from the sketch. This is discussed further in Section 5.3. Finally, we discuss the computation of the inverse-probability estimators $\widehat{f({\nu}_{x})}$ for the sampled keys in Section 5.4. In particular, in order to compute the estimator, we need to know the values $f({\nu}_{x})$ for the keys in the sample, which will require a second pass over the data. The analysis will result in the following main theorem.

Theorem 5.1**.**

Let $k\geq 3$, $0<\varepsilon\leq\frac{1}{2}$, and $f$ be a soft concave sublinear function. Algorithms 4 and 5 produce a stochastic PPSWOR sample of size $k-1$, where each key $x$ has weight $V_{x}$ that satisfies $f({\nu}_{x})\leq\textsf{E}[V_{x}]\leq\frac{1}{(1-\varepsilon)}f({\nu}_{x})$. The per-key inverse-probability estimator of $f({\nu}_{x})$ is unbiased and has variance

[TABLE]

The space required by the sketch at any given time is $O(k)$ in expectation. Additionally, with probability at least $1-\delta$, the space will not exceed $O\left(k+\min\{\log{m},\log{\log{\left(\frac{\operatorname{\textsf{Sum}}_{D}}{\operatorname{\textsf{Min}}(D)}\right)}}\}+\log{\left(\frac{1}{\delta}\right)}\right)$ at any time while processing $D$, where $m$ is the number of elements in $D$, $\operatorname{\textsf{Min}}(D)$ is the minimum value of an element in $D$, and $\operatorname{\textsf{Sum}}_{D}$ is the sum of frequencies of all keys.

Remark 5.2**.**

The parameter $\varepsilon$ mainly affects the run time of processing an element. For each element processed by the stream, we generate $r=\frac{k}{\varepsilon}$ output elements that are then further processed by the sketch. Hence, the run time of processing an element grows with $\frac{1}{\varepsilon}$. The space is affected by $\varepsilon$ when considering worst case over the randomness. The total number of possible keys for output elements is $r$ times the number of active keys, and in the worst case (over the randomness), we may store all of them in Sideline.

The sketch and estimator specification use the following functions in a black-box fashion

[TABLE]

where $a(t)$ is the inverse complement Laplace transform of $f$, as specified in Equation (6) (Section 2.4). Closed expressions for $A(t)$ and $B(t)$ for some common concave sublinear functions $f$ are provided in [10]. These functions are well-defined for any soft concave sublinear $f$. Also note that it suffices to approximate the values of $A$ and $B$ within a small multiplicative error (which will carry over to the variance, see Remark 5.11), so one can also use a table of values to compute the function.

While processing the stream, we will keep track of the sum of values of all elements $\operatorname{\textsf{Sum}}_{D}=\sum_{x\in\mathcal{X}}{{\nu}_{x}}$. We will then use $\operatorname{\textsf{Sum}}_{D}$ to set $\gamma$ adaptively to be $\frac{2\varepsilon}{\operatorname{\textsf{Sum}}_{D}}$. Thus, this is a running “candidate” value that can only decrease over time. The final value of $\gamma$ will be set when we produce a sample from the sketch in Algorithm 5. In the discussion below, we will show that setting $\gamma=\frac{2\varepsilon}{\operatorname{\textsf{Sum}}_{D}}$ satisfies the conditions needed for each of the sketches for $\operatorname{\mathcal{L}^{c}}[a]({\nu}_{x})_{0}^{\gamma}$ and $\operatorname{\mathcal{L}^{c}}[a]({\nu}_{x})_{\gamma}^{\infty}$.

5.1 The Sketch for $\operatorname{\mathcal{L}^{c}}[a]({\nu}_{x})_{0}^{\gamma}$

For the contributions $\operatorname{\mathcal{L}^{c}}[a]({\nu}_{x})_{0}^{\gamma}=\int_{0}^{\gamma}a(t)(1-e^{-{\nu}_{x}t})dt$, we will see that as long as we choose a small enough $\gamma$, $\int_{0}^{\gamma}a(t)(1-e^{-{\nu}_{x}t})dt$ will be approximately $\left(\int_{0}^{\gamma}a(t)tdt\right){\nu}_{x}$, up to a multiplicative $1-\varepsilon$ factor. Note that $\left(\int_{0}^{\gamma}a(t)tdt\right){\nu}_{x}$ is simply the frequency ${\nu}_{x}$ scaled by $B(\gamma)=\int_{0}^{\gamma}a(t)tdt$. A PPSWOR sample is invariant to the scaling, so we can simply use a PPSWOR sampling sketch according to the frequencies ${\nu}_{x}$ (Algorithm 2). The scaling only needs to be considered in a final step when the samples of the two sets of contributions are combined to produce a single sample.777In our implementation (Section 6) we incorporated an optimization where we only keep in the PPSWOR sample elements that may contribute to the final sample.

Lemma 5.3**.**

Let $\varepsilon>0$ and $\gamma\leq\frac{2\varepsilon}{\max_{x}{{\nu}_{x}}}$. Then, for any key $x$,

[TABLE]

Proof.

Consider a key $x$ with frequency ${\nu}_{x}$. Using $1-e^{-z}\leq z$, we get

[TABLE]

Now, using $1-e^{-z}\geq z-\frac{z^{2}}{2}$ for $z\geq 0$,

[TABLE]

Note that $\gamma\leq\frac{2\varepsilon}{\max_{y}{{\nu}_{y}}}\leq\frac{2\varepsilon}{{\nu}_{x}}$. Hence, for every $0\leq t\leq\gamma$, ${\nu}_{x}t\leq 2\varepsilon$, and ${\nu}_{x}t-\frac{({\nu}_{x}t)^{2}}{2}\geq(1-\varepsilon){\nu}_{x}t$. As a result, we get that

[TABLE]

∎

Note that our choice of $\gamma=\frac{2\varepsilon}{\operatorname{\textsf{Sum}}_{D}}$ satisfies the condition of Lemma 5.3.

5.2 The Sketch for $\operatorname{\mathcal{L}^{c}}[a]({\nu}_{x})_{\gamma}^{\infty}$

The sketch for $\operatorname{\mathcal{L}^{c}}[a]({\nu}_{x})_{\gamma}^{\infty}=\int_{\gamma}^{\infty}a(t)(1-e^{-{\nu}_{x}t})dt$ processes elements in the following way. We map each input element $e=(e.\mathit{key},e.\mathit{val})$ into $r=\frac{k}{\varepsilon}$ output elements with keys $(e.\mathit{key},1)$ through $(e.\mathit{key},r)$ and values $Y_{i}\sim\textsf{Exp}[e.\mathit{val}]$ drawn independently. Each of these elements is then processed separately.

The main component of the sketch is a $\operatorname{\textsf{SumMax}}$ sampling sketch of size $k$. Our goal is that for each generated output element $((e.\mathit{key},i),Y_{i})$, the $\operatorname{\textsf{SumMax}}$ sketch will process an element $((e.\mathit{key},i),A(\max\{Y_{i},\gamma\}))$. However, since $\gamma$ decreases over time and we do not know its final value, we only process the elements with $Y_{i}\geq\gamma$ into the $\operatorname{\textsf{SumMax}}$ sketch. We keep the rest of the elements in an auxiliary structure (implemented as a maximum priority queue) that we call Sideline. Every time we update the value of $\gamma$, we remove the elements with $Y_{i}\geq\gamma$ and process them into the $\operatorname{\textsf{SumMax}}$ sketch. Thus, at any time the Sideline structure only contains elements with value less than $\gamma$.888In our implementation (see Section 6) we only keep in Sideline elements that have the potential to modify the SumMax sketch when inserted.

For any active input key $x$ and $i\in[r]$, let $M_{x,i}$ denote the minimum value $Y_{i}$ that was generated with key $(x,i)$. We have the following invariants that we will use in our analysis:

Either the element $((x,i),M_{x,i})$ is in Sideline or an element $((x,i),A(M_{x,i}))$ was processed by the $\operatorname{\textsf{SumMax}}$ sketch. 2. 2.

Since $\gamma$ is decreasing over time, all elements ejected from the Sideline have value that is at least $\gamma$.

We also will use the following property of the sketch.

Lemma 5.4**.**

For any key $x$ that was active in the input and $i\in[r]$, $M_{x,i}\sim\textsf{Exp}[{\nu}_{x}]$ and these random variables are independent for different pairs $(x,i)$.

Proof.

$M_{x,i}$ by definition is the minimum of independent exponential random variables with sum of parameters ${\nu}_{x}$. ∎

In the following lemma, we bound the size of Sideline (and as a result, the entire sketch) with probability $1-\delta$ for $0<\delta<1$ of our choice.

Lemma 5.5**.**

For a set of elements $D$, denote by $m$ the number of elements in $D$, and let $\operatorname{\textsf{Min}}(D)$ denote the minimum value of any element in $D$. The expected number of elements in Sideline at any given time is $O(k)$, and with probability at least $1-\delta$, the number of elements in Sideline will not exceed $O\left(k+\min\{\log{m},\log{\log{\left(\frac{\operatorname{\textsf{Sum}}_{D}}{\operatorname{\textsf{Min}}(D)}\right)}}\}+\log{\left(\frac{1}{\delta}\right)}\right)$ at any time while processing $D$.

The proof is deferred to Appendix D.

5.3 Generating a Sample from the Sketch

The final sample returned by Algorithm 5 is the merge of two samples:

The PPSWOR sample for $\operatorname{\mathcal{L}^{c}}[a]({\nu}_{x})_{0}^{\gamma}$ with frequencies scaled by $B(\gamma)=\int_{0}^{\gamma}ta(t)dt$. 2. 2.

The $\operatorname{\textsf{SumMax}}$ sample for $\operatorname{\mathcal{L}^{c}}[a]({\nu}_{x})_{\gamma}^{\infty}$ with weights scaled by $\frac{1}{r}$. Before the scaling, the $\operatorname{\textsf{SumMax}}$ sample processes an element $(e.\mathit{key},A(\gamma))$ for each remaining $e\in\textnormal{{Sideline}}$.

The scaling is performed using a property of exponential random variables and is formalized in the following lemma.

Lemma 5.6**.**

Given a PPSWOR sample where each key $x$ has frequency ${\nu}_{x}$, we can obtain a PPSWOR sample for the weights $c\cdot{\nu}_{x}$ by returning the original sample of keys but dividing the seed value of each key by $c$.

Proof.

A property of exponential random variables is that if $Y\sim\textsf{Exp}[w]$, then for any constant $c>0$, $y/c\sim\textsf{Exp}[cw]$. Consider the set of seed values $\{\textnormal{{seed}}(x)\mid x\in\mathcal{X}\}$ computed for the original PPSWOR sample according to the frequencies ${\nu}_{x}$. If we divided each seed value by $c$, the seed of key $x$ would come from the distribution $\textsf{Exp}(c{\nu}_{x})$. Hence, a PPSWOR sample according to the weights $c{\nu}_{x}$ would contain the $k$ keys with lowest seed values after dividing by $c$, and these $k$ keys are the same keys that have lowest seed values before dividing by $c$. ∎

Denote by $E$ the set of all elements that are passed on to the $\operatorname{\textsf{SumMax}}$ sketch, either during the processing of the set of elements $D$ or in the final phase.

Lemma 5.7**.**

The final sample computed by Algorithm 5 is a PPSWOR sample with respect to weights

[TABLE]

Proof.

The sample $s.\textit{ppswor}$ in Algorithm 5 is a PPSWOR sample with respect to frequencies ${\nu}_{x}$, which is then scaled by $\int_{0}^{\gamma}ta(t)dt$ to get a PPSWOR sample according to the weights ${\nu}_{x}\int_{0}^{\gamma}ta(t)dt$. The sample $s.\operatorname{\textsf{SumMax}}$ is a $\operatorname{\textsf{SumMax}}$ sample, which by Corollary 4.2 is a PPSWOR sample according to the weights $\operatorname{\textsf{SumMax}}_{E}(x)$. This sample is scaled to be a PPSWOR sample according to the weights $\frac{1}{r}\operatorname{\textsf{SumMax}}_{E}(x)$.

Note that these samples are independent. When we perform a bottom-$k$ merge of the two samples, the seed of key $x$ is then the minimum of two independent exponential random variables with parameters $\int_{0}^{\gamma}ta(t)dt$ and $\frac{1}{r}\operatorname{\textsf{SumMax}}_{E}(x)$. Therefore, the distribution of $\textnormal{{seed}}(x)$ in the merged sample is $\textsf{Exp}(\int_{0}^{\gamma}ta(t)dt+\frac{1}{r}\operatorname{\textsf{SumMax}}_{E}(x))$, which means that the sample is a PPSWOR sample according to those weights, as desired. ∎

We next interpret $\frac{1}{r}\operatorname{\textsf{SumMax}}_{E}(x)$. From the invariants listed in Section 5.2 and the description of Algorithm 5, we have that for any active input key $x$ and $i\in[r]$, the element $((x,i),A(\max\{\gamma,M_{x,i}\}))$ was processed by the $\operatorname{\textsf{SumMax}}$ sketch. Because $A$ is monotonically non-increasing,

[TABLE]

Now, $\frac{1}{r}\operatorname{\textsf{SumMax}}_{E}(x)=\sum_{i=1}^{r}{\frac{1}{r}\operatorname{\textsf{Max}}_{E}((x,i))}$. By Lemma 5.4, the summands $\frac{1}{r}\operatorname{\textsf{Max}}_{E}((x,i))$ are independent random variables for every $x$ and $i$. By the monotonicity of the function $A$, each summand $\frac{1}{r}\operatorname{\textsf{Max}}_{E}((x,i))$ is in the range $\left[0,\frac{A(\gamma)}{r}\right]$.

The final sample returned by Algorithm 5 is then a stochastic PPSWOR sample as defined in Section 3. The weight of key $x$ also includes the deterministic component ${\nu}_{x}\int_{0}^{\gamma}ta(t)dt$, which can be larger than $\frac{A(\gamma)}{r}$. However, since this is summand is deterministic, we can break it into smaller deterministic parts, each of which will be at most $\frac{A(\gamma)}{r}$. This way, the sample still satisfies the condition that the weight of every key is the sum of independent random variables in $\left[0,\frac{A(\gamma)}{r}\right]$. The next step is to show that it satisfies the conditions of Theorem 3.1.

Lemma 5.8**.**

For a key $x$, define $V_{x}=\frac{1}{r}\operatorname{\textsf{SumMax}}_{E}(x)+{\nu}_{x}\int_{0}^{\gamma}ta(t)dt$. Let $V=\sum_{x\in\mathcal{X}}{V_{x}}$. Then, for any $0<\varepsilon\leq\frac{1}{2}$ and $r\geq\frac{k}{\varepsilon}$,

[TABLE]

The proof, which is deferred to Appendix D, uses the following lemma which is due to Cohen [10].

Lemma 5.9**.**

[10]** For every input key $x$ and $i=1,\ldots,r$,

[TABLE]

The following theorem combines the previous lemmas to show how to estimate $f({\nu}_{x})$ using the sample and bound the variance. Note that we need to specify how to compute the estimator. For getting $f({\nu}_{x})$ we make another pass, and the computation of the conditioned inclusion probability is described in the next subsection.

Theorem 5.10**.**

The sample returned by Algorithm 5 is a stochastic PPSWOR sample, where each key $x$ has weight $V_{x}$ that satisfies $f({\nu}_{x})\leq\textsf{E}[V_{x}]\leq\frac{1}{(1-\varepsilon)}f({\nu}_{x})$. The per-key inverse-probability estimator according to the weights $f({\nu}_{x})$,

[TABLE]

is unbiased and has variance

[TABLE]

where $V=\sum_{x\in\mathcal{X}}{V_{x}}$.

Proof.

We first prove that $f({\nu}_{x})\leq\textsf{E}[V_{x}]\leq\frac{1}{(1-\varepsilon)}f({\nu}_{x})$. The randomized weight $V_{x}$ is the sum of two terms:

[TABLE]

As in the proof of Lemma 5.8, using Lemma 5.9,

[TABLE]

The quantity ${\nu}_{x}\int_{0}^{\gamma}ta(t)dt$ is deterministic, and since $\gamma\leq\frac{2\varepsilon}{\max_{x}{{\nu}_{x}}}$, by Lemma 5.3,

[TABLE]

Recall that $f({\nu}_{x})=\int_{0}^{\infty}a(t)(1-e^{-{\nu}_{x}t})dt=\int_{0}^{\gamma}a(t)(1-e^{-{\nu}_{x}t})dt+\int_{\gamma}^{\infty}a(t)(1-e^{-{\nu}_{x}t})dt$. Combining Equations (8) and (9), we get that

[TABLE]

Now consider the estimator $\widehat{f({\nu}_{x})}$ for $f({\nu}_{x})$. To show that the estimator is unbiased, that is, $\textsf{E}\left[\widehat{f({\nu}_{x})}\right]=f({\nu}_{x})$, we can follow the proof of Claim 2.4 exactly as written earlier. It is left to bound the variance. For the sake of the analysis, consider the following estimator for $\textsf{E}[V_{x}]$:

[TABLE]

This is again the inverse-probability estimator from Definition 2.3. Our sample is a stochastic PPSWOR sample according to the weights $V_{x}$, where each one of $V_{x}$ is a sum of independent random variables in $\left[0,\frac{A(\gamma)}{r}\right]$ (recall that the deterministic part can also be expressed as a sum with each summand in $\left[0,\frac{A(\gamma)}{r}\right]$). Lemma 5.8 shows that $\textsf{E}[V]\geq\frac{A(\gamma)}{r}\cdot k$. Hence, we satisfy the conditions of Theorem 3.1, which in turn shows that

[TABLE]

Finally, note that $\widehat{f({\nu}_{x})}=\frac{f({\nu}_{x})}{\textsf{E}[V_{x}]}\cdot\widehat{\textsf{E}[V_{x}]}$. We established above that $\frac{f({\nu}_{x})}{\textsf{E}[V_{x}]}\leq 1$ and $\textsf{E}[V_{x}]\leq\frac{1}{(1-\varepsilon)}f({\nu}_{x})$. We conclude that

[TABLE]

∎

Remark 5.11**.**

The theorem establishes a bound on the variance of the estimator for $\textsf{E}[V_{x}]$, and then uses it to bound the variance of the estimator for $f({\nu}_{x})$, which is possible since $\textsf{E}[V_{x}]$ approximates $f({\nu}_{x})$. This results in increasing the variance by an $\varepsilon$-dependent factor. Similarly, if we wish to estimate $f({\nu}_{x})$ for a concave sublinear function $f$ (and not a soft concave sublinear function), we can use the same idea and lose another constant factor in the variance.

5.4 Expressing Conditioned Inclusion Probabilities

To compute the estimator $\widehat{f({\nu}_{x})}$, we need to know both $f({\nu}_{x})$ and the precise conditioned inclusion probability $\Pr[\textnormal{{seed}}(x)<\tau]$. In order to get $f({\nu}_{x})$, we perform a second pass over the data elements to obtain the exact frequencies ${\nu}_{x}$ for $x\in S$. This can be done via a simple composable sketch that collects and sums the values of data elements with keys that occur in the sample $S$.

We next consider computing the conditioned inclusion probabilities. The following lemma considers the seed distributions of keys in the final sample. It shows that the distributions are parameterized by ${\nu}_{x}$ and describes their CDF.

Lemma 5.12**.**

Algorithms 4 and 5 describe a bottom-$k$ sampling scheme, where in the output sample the seed of each key $x$ is drawn from a distribution $\textsf{SeedDist}^{(F)}[{\nu}_{x}]$. The distribution $\textsf{SeedDist}^{(F)}[w]$ has the following cumulative distribution function:

[TABLE]

where

[TABLE]

The proof is deferred to Appendix D.

6 Experiments

We implemented our sampling sketch and report here the results of experiments on real and synthetic datasets. Our experiments are small-scale and aimed to demonstrate the simplicity and practicality of our sketch design and to understand the actual space and error bounds (that can be significantly better than our worst-case bounds).

6.1 Implementation

Our Python 2.7 implementation follows the pseudocode of the sampling sketch (Algorithm 4), the PPSWOR (Algorithm 2) and $\operatorname{\textsf{SumMax}}$ (Algorithm 3) substructures, the sample production from the sketch (Algorithm 5), and the estimator (that evaluates the conditioned inclusion probabilities, see Section 5.4). We incorporated two practical optimizations that are not shown in the pseudocode. These optimizations do not affect the outcome of the computation or the worst-case analysis, but reduce the sketch size in practice.

Removing redundant keys from the PPSWOR subsketch

The pseudocode (Algorithm 4) maintains two samples of size $k$, the PPSWOR and the $\operatorname{\textsf{SumMax}}$ samples. The final sample of size $k$ is obtained by merging these two samples. Our implementation instead maintains a truncated PPSWOR sketch that removes elements that are already redundant (do not have a potential to be included in the merged sample). We keep an element in the PPSWOR sketch only when the seed value is lower than $rB(\gamma)$ times the current threshold $\tau$ of the $\operatorname{\textsf{SumMax}}$ sample. This means that the “effective” inclusion threshold we use for the PPSWOR sketch is the minimum of the $k$th largest (the threshold of the PPSWOR sketch) and $rB(\gamma)\tau$. To establish that elements that do not satisfy this condition are indeed redundant, recall that when we later merge the PPSWOR and the $\operatorname{\textsf{SumMax}}$ samples, the value of $B(\gamma)$ can only become lower and the $\operatorname{\textsf{SumMax}}$ threshold can only be lower, making inclusion more restrictive. This optimization may result in maintaining much fewer than $k$ elements and possibly an empty PPSWOR sketch. The benefit is larger for functions when $A(t)$ is bounded (as $t$ approaches [math]). In particular, when $a(t)=0$ for $t\leq\gamma$ we get $B(\gamma)=0$ and the truncation will result in an empty sample.

Removing redundant elements from Sideline

The pseudocode may place elements in Sideline that have no future potential of modifying the $\operatorname{\textsf{SumMax}}$ sketch. In our implementation, we place and keep an element $((e.\mathit{key},i),Y)$ in Sideline only as long as the following condition holds: If $((e.\mathit{key},i),A(Y))$ is processed by the current $\operatorname{\textsf{SumMax}}$ sketch, it would modify the sketch. To establish redundancy of discarded elements, note that when an element is eventually processed, the value it is processed with is at most $A(Y)$ (can be $A(\gamma)$ for $\gamma\geq Y$) and also at that point the $\operatorname{\textsf{SumMax}}$ sketch threshold can only be more restrictive.

6.2 Datasets and Experimental Results

We used the following datasets for the experiments:

• •

abcnews [28]: News headlines published by the Australian Broadcasting Corp. For each word, we created an element with value $1$.

• •

flicker [36]: Tags used by Flickr users to annotate images. The key of each element is a tag, and the value is the number of times it appeared in a certain folder.

• •

Three synthetic generated datasets that contain $2\times 10^{6}$ data elements. Each element has value $1$, and the key was chosen according to the Zipf distribution (numpy.random.zipf), with Zipf parameter values $\alpha\in\{1.1,1.2,1.5\}$. The Zipf family in this range is often a good model to real-world frequency distributions.

We applied our sampling sketch with sample size parameter values $k\in\{25,50,75,100\}$ and set the parameter $\varepsilon=0.5$ in all experiments. We sampled according to two concave sublinear functions: the frequency moment $f({\nu})={\nu}^{0.5}$ and $f({\nu})=\ln(1+{\nu})$.

Tables 1 and 2 report aggregated results of 200 repetitions for each combination of dataset, $k$, and $f$ values. In each repetition, we were using the final sample to estimate the sum $\sum_{x\in\mathcal{X}}{f({\nu}_{x})}$ over all keys. For error bounds, we list the worst-case bound on the CV (which depends only on $k$ and $\varepsilon$ and is proportional to $1/\sqrt{k}$) and report the actual normalized root of the average squared error (NRMSE). In addition, we report the NRMSE that we got from 200 repetitions of estimating the same statistics using two common sampling schemes for aggregated data, PPSWOR and priority sampling, which we use as benchmarks.

Also, we consider the size of the sketch after processing each element. Since the representation of each key can be explicit and require a lot of space, we separately consider the number of distinct keys and the number of elements stored in the sketch. We report the maximum number of distinct keys stored in the sketch at any point (the average and the maximum over the 200 repetitions) and the respective maximum number of elements stored in the sketch at any point during the computations (again, the average and the maximum over the 200 repetitions). The size of the sketch is measured at end of the processing of each input element – during the processing we may store one more distinct key and temporarily store up to $r=k/\varepsilon$ additional elements in the Sideline.

We can see that the actual error reported is significantly lower than the worst-case bound. Furthermore, the error that our sketch gets is close to the error achieved by the two benchmark sampling schemes. We can also see that the maximum number of distinct keys stored in the sketch at any time is relatively close to the specified sample size of $k$ and that the total sketch size in terms of elements rarely exceeded $3k$, with the relative excess seeming to decrease with $k$. In comparison, the benchmark schemes require space that is the number of distinct keys (for the aggregation), which is significantly higher than the space required by our sketch.

7 Conclusion

We presented composable sampling sketches for weighted sampling of unaggregated data tailored to a concave sublinear function of the frequencies of keys. We experimentally demonstrated the simplicity and efficacy of our design: Our sketch size is nearly optimal in that it is not much larger than the final sample size, and the estimate quality is close to that provided by a weighted sample computed directly over the aggregated data.

Acknowledgments

Ofir Geri was supported by NSF grant CCF-1617577, a Simons Investigator Award for Moses Charikar, and the Google Graduate Fellowship in Computer Science in the School of Engineering at Stanford University. The computing for this project was performed on the Sherlock cluster. We would like to thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that contributed to these research results.

Appendix A Proofs Deferred from Section 2

Proof of Proposition 2.1.

Each data element $e=(e.\mathit{key},e.\mathit{val})$ is processed by giving it a score $\textnormal{{ElementScore}}(e)\sim\textsf{Exp}(e.\mathit{val})$ and then processing the element $(e.\mathit{key},\textnormal{{ElementScore}}(e))$ by the bottom-$k$ structure. For a key $x$, we define

[TABLE]

to be the smallest score of an element with key $x$.

Since $\textnormal{{seed}}(x)$ is the minimum of independent exponential random variables, its distribution is $\textsf{Exp}(w_{x})$. After processing the elements in $D$, the bottom-$k$ structure contains the $k$ pairs $(x,\textnormal{{seed}}(x))$ with smallest $\textnormal{{seed}}(x)$ values, and hence obtains the respective PPSWOR sample.

Consider now the sketch resulting from merging two sampling sketches computed for $D_{1}$ and $D_{2}$. For each key $x$, denote by $\textnormal{{seed}}_{1}(x)$ and $\textnormal{{seed}}_{2}(x)$ the values of $x$ in the sketches for $D_{1}$ and $D_{2}$, respectively. Then, the merged sketch contains the $k$ pairs $(x,\textnormal{{seed}}(x))$ with smallest $\textnormal{{seed}}(x):=\min\{seed_{1}(x),seed_{2}(x)\}$ values. As $\textnormal{{seed}}(x)$ is the minimum of two independent exponential random variables with parameters $\operatorname{\textsf{Sum}}_{D_{1}}(x)$ and $\operatorname{\textsf{Sum}}_{D_{2}}(x)$, we get that $\textnormal{{seed}}(x)\sim\textsf{Exp}(\operatorname{\textsf{Sum}}_{D_{1}\cup D_{2}}(x))$, as desired. ∎

Appendix B Proofs Deferred from Section 3

B.1 Threshold Distribution and Fixed-Threshold Inclusion Probability

Before proceeding, we establish two technical lemmas that will be useful later. The first lemma shows that the distribution of the $k$-th lowest seed is dominated by the Erlang distribution which takes the sum of the expected weights $V$ as a parameter. The lemma will be useful later when we consider the inclusion threshold $\tau_{x}$.

Lemma B.1**.**

Consider a set of keys $\mathcal{X}$, such that the weight of each $x\in\mathcal{X}$ is a random variable $S_{x}\geq 0$. Let $V=\sum_{x\in\mathcal{X}}{\textsf{E}[S_{x}]}$, and for each $x$ with $S_{x}>0$, draw $\textnormal{{seed}}(x)\sim\textsf{Exp}[S_{x}]$ independently. Then, the distribution of the $k$-th lowest seed, $\{\textnormal{{seed}}(x)\mid x\in X\}_{(k)}$, is dominated by $\textsf{Erlang}[V,k]$.

We first establish the dominance relation $\textsf{Exp}[S]\preceq\textsf{Exp}[\textsf{E}[S]]$ for any nonnegative random variable $S$.

Lemma B.2**.**

Let $S\geq 0$ be a random variable. Let $X$ be a random variable such that $X\sim\textsf{Exp}[S]$ when $S>0$ and $X=\infty$ otherwise.999$X$ takes values in $\mathbb{R}$ and cannot be $\infty$. However, the random variable $X$ represents the minimum seed value (or the inclusion threshold $\tau_{x}$ as in Section 2.3), and the event $X\leq\tau$ represents whether the inclusion threshold for key $x$ is at most $\tau$. The case $S=0$ corresponds to no elements generated with keys in $\mathcal{X}\setminus\{x\}$, so we can say that the inclusion threshold for $x$ is $\infty$ (the event we care about is whether $x$ enters the sample or not). Here we are trying to show a distribution that dominates the distribution of the inclusion threshold, and for that purpose, any threshold $\tau>0$ is more restrictive than $\infty$. From a technical perspective, when $S=0$, we can still use the CDF of $\textsf{Exp}[S]$ since $\Pr[X\leq\tau]=0=1-e^{-S\tau}$. Later, when we consider the $k$-th lowest seed, we will similarly allow it to be $\infty$ when less than $k$ keys are active. Then, the distribution of $X$ is dominated by $\textsf{Exp}[\textsf{E}[S]]$, that is, $\forall\tau,\Pr_{X\sim\textsf{Exp}[S]}[X\leq\tau]\leq 1-e^{-E[S]\tau}$.

Proof.

Follows from Jensen’s inequality:

[TABLE]

∎

Therefore, for any key $x$, the $\textnormal{{seed}}(x)$ distribution with stochastic weights is dominated by $\textsf{Exp}[\textsf{E}[S_{x}]]$, which is the distribution used by PPSWOR according to the expected weights. We now consider the distribution of $\{\textnormal{{seed}}(x)\mid x\in X\}_{(k)}$, which is the $k$-th lowest seed value. We show that the distribution of the $k$-th lowest seed is dominated by $\textsf{Erlang}[V,k]$ (recall that $V=\sum_{x}{\textsf{E}[S_{x}]}$). In the proof, we will use the following property of dominance (the proof of the following property is standard and included here for completeness).

Claim B.3**.**

Let $X_{1},\ldots,X_{k},Y_{1},\ldots,Y_{k}$ be independent random variables such that the distribution of $X_{i}$ is dominated by that of $Y_{i}$. Then the distribution of $X_{1}+\ldots+X_{k}$ is dominated by that of $Y_{1}+\ldots+Y_{k}$.

Proof.

We prove for $k=2$ (the proof for $k>2$ follows from a simple induction argument). Denote by $f_{i}$ and $F_{i}$ the PDF and CDF functions of $X_{i}$, respectively, and by $g_{i}$ and $G_{i}$ the PDF and CDF of $Y_{i}$. From the dominance assumption, we know that $F_{i}(t)\leq G_{i}(t)$ for all $i$. Now,

[TABLE]

∎

We are now ready to prove Lemma B.1.

Proof of Lemma B.1.

Conditioned on the values of $S_{x}$, the distribution of $\{\textnormal{{seed}}(x)\mid x\in X\}_{(k)}$ is dominated by $\textsf{Erlang}[\sum_{x\in\mathcal{X}}{S_{x}},k]$ (for a proof, see Appendix C in [11]). The distribution of $\{\textnormal{{seed}}(x)\mid x\in X\}_{(k)}$ (unconditioned on the values of $S_{x}$) is a linear combination of distributions, which are each dominated by the respective Erlang distribution. Using the definition of dominance (Definition 2.6) and the law of total probability, we get that the distribution of $\{\textnormal{{seed}}(x)\mid x\in X\}_{(k)}$ (unconditioned on $S_{x}$) is dominated by $\textsf{Erlang}[\sum_{x\in\mathcal{X}}{S_{x}},k]$ (unconditioned on $S_{x}$). A random variable drawn from $\textsf{Erlang}[\sum_{x\in\mathcal{X}}{S_{x}},k]$ has the same distribution as the sum of $k$ independent random variables drawn from $\textsf{Exp}(\sum_{x\in\mathcal{X}}{S_{x}})$. The distribution of each of these $k$ exponential random variables is dominated by $\textsf{Exp}(V)$ (by Lemma B.2). Using Claim B.3, we get that $\textsf{Erlang}[\sum_{x\in\mathcal{X}}{S_{x}},k]\preceq\textsf{Erlang}[V,k]$. The assertion of the lemma then follows from the transitivity of dominance. ∎

The second lemma provides lower bounds on the CDF of $\textsf{Exp}(S)$ under certain conditions.

Lemma B.4**.**

Let the random variable $S=\sum_{i=1}^{r}{S_{i}}$ be a sum of $r$ independent random variables in the range $[0,T]$. Let $v=\textsf{E}[S]$. Then,

[TABLE]

In the regime $\tau T<1$, we get that the probability of being less than $\tau$ is close to that of $\textsf{Exp}(\textsf{E}[S])$.

Lemma B.5**.**

Let $S$ be a random variables in $[0,T]$ with expectation $\textsf{E}[S]=v$. Then for all $\tau$,

[TABLE]

Proof.

Denote the probability density function of $S$ by $p_{S}$. Conditioned on the value of $S$, the probability of $X\sim\textsf{Exp}[S]$ being below $\tau$ is

[TABLE]

It follows that

[TABLE]

Consider the function $f(x)=1-e^{-x\tau}$ for a fixed $\tau\geq 0$. Since $f$ is concave, for every $x\in[0,T]$,

[TABLE]

By monotonicity,

[TABLE]

and finally,

[TABLE]

∎

Lemma B.6**.**

Let the random variable $S=\sum_{i=1}^{r}{S_{i}}$ be a sum of $r$ independent random variables in the range $[0,T]$. Let $v=\textsf{E}[S]$. Then,

[TABLE]

Proof.

Let $X\sim\textsf{Exp}[S]$. Since $S=\sum_{i=1}^{r}{S_{i}}$, we could define $r$ independent exponential random variables $X_{i}\sim\textsf{Exp}[S_{i}]$. $X$ has the same distribution as $\min_{1\leq i\leq r}{X_{i}}$. Hence,

[TABLE]

where the last inequality follows from the arithmetic mean-geometric mean inequality. Now, using the inequality $1-x\leq e^{-x}$ (for any $x\in\mathbb{R}$), and the fact that $f(x)=x^{r}$ is non-decreasing for $x,r\geq 0$, we get that

[TABLE]

Consequently,

[TABLE]

∎

Proof of Lemma B.4.

Follows from Lemma B.6 using the inequality $1-e^{-x}\geq x-x^{2}/2$ for $x\geq 0$. ∎

B.2 Variance Bounds for the Inverse-Probability Estimator

Proof of Theorem 3.1.

We start bounding the per-key variance as in Claim 2.5:

[TABLE]

By Lemma B.1, we know that the distribution of $\tau_{x}$ (the $k-1$ lowest seed of the keys in $\mathcal{X}\setminus\{x\}$) is dominated by $\textsf{Erlang}[V,k-1]$, hence

[TABLE]

To bound the first summand, since $t\leq\frac{1}{T}$, we get from Lemma B.4 (applied to $\textnormal{{seed}}(x)$) that $\Pr[\textnormal{{seed}}(x)<t]\geq 1-e^{-v_{x}t\left(1-\frac{tT}{2}\right)}\geq 1-e^{-v_{x}t/2}$. It follows that

[TABLE]

To bound the second summand, since $t>\frac{1}{T}$, $\Pr[\textnormal{{seed}}(x)<t]\geq\Pr[\textnormal{{seed}}(x)<1/T]\geq 1-e^{-v_{x}/2T}$. Subsequently,

[TABLE]

Combining, we get that

[TABLE]

∎

B.3 Inclusion Probability in a Stochastic Sample

Proof of Theorem 3.2.

We first separately deal with the case where there is only one key, which we denote $x$. In this case, $V=v_{x}$, and if $S_{x}>0$, then $x$ is included in the sample. Otherwise, the sample is empty. In the proof of Lemma B.4, when $S=0$, we used $\Pr[X\leq\tau]=1-e^{-s\tau}=0$ and the event $X\leq\tau$ does not happen. Hence, we can use Lemma B.4 to bound $\Pr[S_{x}>0]\geq\Pr[\textnormal{{seed}}(x)\leq\tau]$ for any $\tau>0$. We pick $\tau=\frac{2\varepsilon}{T}$ and get that $x$ is included in the sample with probability

[TABLE]

using $V\geq\frac{1}{\varepsilon}\ln{\left(\frac{1}{\varepsilon}\right)}T$ and $2(1-\varepsilon)\geq 1$.

If there is more than one key, a key $x$ is included in the sample if $\textnormal{{seed}}(x)$ is smaller than the seed of all other keys. The distribution of $\min_{z\neq x}{\textnormal{{seed}}(z)}$ is $\textsf{Exp}\left(\sum_{z\neq x}{S_{z}}\right)$, which is dominated by $\textsf{Exp}\left(\sum_{z\neq x}{v_{z}}\right)$ (Lemma B.2). Then,

[TABLE]

∎

Appendix C Proofs Deferred from Section 4

Proof of Lemma 4.1.

With a slight abuse of notation, for a full key $z=(z.p,z.s)$ we define

[TABLE]

Now, since we use the same value $h(z)$ for all elements with key $z$, the minimum $\textnormal{{ElementScore}}(e)$ value generated for an element $e\in D$ with key $e.\mathit{key}=z$ is $h(z)/\operatorname{\textsf{Max}}_{D}(z)$:

[TABLE]

Recall that for $X\sim\textsf{Exp}[1]$ and $a>0$, the distribution of $X/a$ is $\textsf{Exp}[a]$, and that $h(z)\sim\textsf{Exp}[1]$. Therefore, the algorithm effectively draws $\textnormal{{seed}}_{D}(z)\sim\textsf{Exp}[\operatorname{\textsf{Max}}_{D}(z)]$ for every key $z$. Moreover, from our assumption of independence of $h$, the variables $\textnormal{{seed}}_{D}(z)$ of different keys $z$ are also independent.

We now notice that for a primary key $x$,

[TABLE]

That is, $\textnormal{{seed}}_{D}(x)$ is the minimum, over all keys $z$ with primary key $z.p=x$ that appeared in at least one element of $D$, of $\textnormal{{seed}}_{D}(z)$.

The random variables $\textnormal{{seed}}_{D}(z)$ for input keys $z$ are independent and exponentially distributed with respective parameters $\operatorname{\textsf{Max}}_{D}(z)$. From properties of the exponential distribution, their minimum is also exponentially distributed with a parameter that is equal to the sum of their parameters $\operatorname{\textsf{Max}}_{D}(z)$:

[TABLE]

that is, $\textnormal{{seed}}_{D}(x)\sim\textsf{Exp}[\operatorname{\textsf{SumMax}}_{D}(x)]$. Moreover, the independence of $\textnormal{{seed}}_{D}(x)$ (for primary keys $x$) follows from the independence of $\textnormal{{seed}}_{D}(z)$ (for input keys $z$). ∎

Appendix D Proofs Deferred from Section 5

Proof of Lemma 5.5.

Consider a fixed time during the processing of $D$ by Algorithm 4 (after some but potentially not all elements have been processed). For each key $x$, let $v_{x}$ be the sum of values of elements with key $x$ that have been processed so far.

For any $t>0$, key $x\in\mathcal{X}$, and $i\in[r]$, we define an indicator random variable $I^{t}_{x,i}$ for the event that an element with key $(x,i)$ was generated with value less than $t$. In particular, the number of elements in Sideline is $\sum_{x\in\mathcal{X}}{\sum_{i=1}^{r}{I^{\gamma}_{x,i}}}$. The event $I^{t}_{x,i}=1$ is the event that the minimum value of the elements generated with key $(x,i)$ is at most $t$. The distribution of the minimum value of these elements is $\textsf{Exp}(v_{x})$, and it follows that

[TABLE]

In particular, when $t=\gamma=\frac{2\varepsilon}{\sum_{z\in\mathcal{X}}{v_{z}}}$ and $r=\frac{k}{\varepsilon}$, we get

[TABLE]

From Chernoff bounds,

[TABLE]

Applying this each time an element is processed and taking union bound, we get that the size of Sideline increases beyond $4k+3\ln{m}+3\ln{\left(\frac{1}{\delta}\right)}$ at any time with probability at most $\delta$.

We now improve the bound to use $\log{\log{\left(\frac{\operatorname{\textsf{Sum}}_{D}}{\operatorname{\textsf{Min}}(D)}\right)}}$ instead of $\log{m}$. Let $t>0$. Consider all the times where the value of $\gamma$ is in the interval $\left[\frac{1}{2}t,t\right]$, and for every $\gamma^{\prime}$ in that interval, let $v_{x}(\gamma^{\prime})$ denote the frequency of key $x$ at the time where $\gamma=\gamma^{\prime}$. Since $\gamma$ decreases over time as elements are processed, any generated element stored in Sideline when $\gamma\in\left[\frac{1}{2}t,t\right]$ must have value at most $t$. Since only more elements are generated as $\gamma$ decreases, we can look all the elements that have been generated until $\gamma$ reached $\frac{1}{2}t$.101010It may be the case that $\gamma$ is never $\frac{1}{2}t$, but in that case we consider the minimum value of $\gamma$ that is at least $\frac{1}{2}t$.

We bound the number of elements with value at most $t$ that have been generated until the time where $\gamma=\frac{1}{2}t$. From the way we set $\gamma$ in Algorithm 4, we get as long as $\gamma\geq\frac{1}{2}t$, $\sum_{x\in\mathcal{X}}{v_{x}(\gamma)}t\leq 4\varepsilon$. Now, consider the indicator $I^{t}_{x,i}$ as defined above for the time where $\gamma=\frac{1}{2}t$. The number of elements stored in Sideline at any time when $\gamma\in\left[\frac{1}{2}t,t\right]$ is at most $\sum_{x\in\mathcal{X}}{\sum_{i=1}^{r}{I^{t}_{x,i}}}$. We get that

[TABLE]

and using Chernoff bounds,

[TABLE]

Finally, the minimum value $\gamma$ can get is $\frac{2\varepsilon}{\operatorname{\textsf{Min}}(D)}$, and the maximum value is $\frac{2\varepsilon}{\operatorname{\textsf{Sum}}_{D}}$. Hence, we can divide the interval of possible values for $\gamma$ into $\lceil\log{\left(\frac{\operatorname{\textsf{Sum}}_{D}}{\operatorname{\textsf{Min}}(D)}\right)}\rceil$ intervals of the form $\left[\frac{1}{2}t,t\right]$, and apply the bound in Equation (10) to each one of them. By the union bound, we get that the probability that the size of Sideline exceeds $8k+3\ln{\lceil\log{\left(\frac{\operatorname{\textsf{Sum}}_{D}}{\operatorname{\textsf{Min}}(D)}\right)}\rceil}+3\ln{\left(\frac{1}{\delta}\right)}$ at any time during the processing of $D$ is at most $\delta$. ∎

Proof of Lemma 5.8.

Using Lemma 5.9, for every key $x$,

[TABLE]

Recall that $\gamma=\frac{2\varepsilon}{\operatorname{\textsf{Sum}}_{D}}$. Then, using $1-e^{-x}\geq\frac{x}{2}$ for $0\leq x\leq 1$,

[TABLE]

Since $r\geq\frac{k}{\varepsilon}$, we conclude that

[TABLE]

∎

Proof of Lemma 5.12.

Consider a key $x$. The seed $\textnormal{{seed}}^{(F)}(x)$ in the output sample is the minimum of $\textnormal{{seed}}^{(1)}(x)$ and $\textnormal{{seed}}^{(2)}$, which are the seed values obtained by the scaled PPSWOR and the $\operatorname{\textsf{SumMax}}$ samples, respectively.

The scaled PPSWOR sample is computed with respect to the weights ${\nu}_{x}B(\gamma)$, and thus $\textnormal{{seed}}^{(1)}(x)\sim\textsf{Exp}[{\nu}_{x}B(\gamma)]$. Therefore using the density function of $\textsf{Exp}[{\nu}_{x}B(\gamma)]$, we get that for all $t>0$,

[TABLE]

The scaled $\operatorname{\textsf{SumMax}}$ sample is a PPSWOR sample with respect to weights $\frac{1}{r}\operatorname{\textsf{SumMax}}_{E}(x)$. Therefore, $\textnormal{{seed}}^{(2)}(x)\sim\textsf{Exp}[\frac{1}{r}\operatorname{\textsf{SumMax}}_{E}(x)]$. Note however that $\frac{1}{r}\operatorname{\textsf{SumMax}}_{E}(x)$ is itself a random variable and in particular, the value $\operatorname{\textsf{SumMax}}_{E}(x)$ is not available to us with the sample. We recall that $\operatorname{\textsf{SumMax}}_{E}(x)=\sum_{i=1}^{r}\operatorname{\textsf{Max}}_{E}((x,i))$ where $\operatorname{\textsf{Max}}_{E}((x,i))$ are i.i.d. random variables. Using properties of the exponential distribution, we know that $\textsf{Exp}[\frac{1}{r}\operatorname{\textsf{SumMax}}_{E}(x)]$ is the same distribution as the minimum of $r$ independent random variables drawn from $\textsf{Exp}[\frac{1}{r}\operatorname{\textsf{Max}}_{E}((x,1))],\ldots,\textsf{Exp}[\frac{1}{r}\operatorname{\textsf{Max}}_{E}((x,r))]$. Therefore, for $t>0$,

[TABLE]

We now express $\Pr[\textsf{Exp}[\frac{1}{r}\operatorname{\textsf{Max}}_{E}((x,i))]>t]$ using the fact that $\operatorname{\textsf{Max}}_{E}((x,i))=A(\max\{y,\gamma\})$ for $y\sim\textsf{Exp}[{\nu}_{x}]$:

[TABLE]

Since $\Pr[\textnormal{{seed}}^{(2)}(x)>t]=p_{2}^{r}$ and using the fact that $\textnormal{{seed}}^{(1)}(x)$ and $\textnormal{{seed}}^{(2)}(x)$ are independent, we conclude that

[TABLE]

∎