Quantiles over data streams: experimental comparisons, new analyses, and further improvements

A fundamental problem in data management and analysis is to generate descriptions of the distribution of data. It is most common to give such descriptions in terms of the cumulative distribution, which is characterized by the quantiles of the data. The design and engineering of efficient methods to find these quantiles has attracted much study, especially in the case where the data are given incrementally, and we must compute the quantiles in an online, streaming fashion. While such algorithms have proved to be extremely useful in practice, there has been limited formal comparison of the competing methods, and no comprehensive study of their performance. In this paper, we remedy this deficit by providing a taxonomy of different methods and describe efficient implementations. In doing so, we propose new variants that have not been studied before, yet which outperform existing methods. To illustrate this, we provide detailed experimental comparisons demonstrating the trade-offs between space, time, and accuracy for quantile computation.

1. Note that floating-point numbers in standard representations (e.g., IEEE 754) can be mapped to integers in a fixed universe in an order-preserving fashion.

2. http://www.minorplanetcenter.net/iau/ecs/mpcat-obs/mpcat-obs.html.

3. Right ascension is an astronomical term used to locate a point (a minor planet in this case) in the equatorial coordinate system.

4. http://www.ncfloodmaps.com.

5. It is possible for the error to be affected due to more duplicates in smaller universes, but we found this effect negligible in practice.

The first three authors are supported by HKRGC Grants GRF-621413, GRF-16211614, and GRF-16200415, and by a Microsoft Grant MRA14EG05. The work of GC is supported in part by European Research Council Grant ERC-2014-CoG 647557, the Yahoo Faculty Research and Engagement Program and a Royal Society Wolfson Research Merit Award.

Authors and Affiliations

1. Department of Compute Science and Engineering, HKUST, Kowloon, Hong Kong

   Ge Luo, Lu Wang & Ke Yi

2. Department of Computer Science, University of Warwick, Coventry, UK

   Graham Cormode

Corresponding author

Correspondence to Ge Luo.

Additional lemmas and proofs

1.1 Size of the truncated tree \(\hat{T}\)

Lemma 1

The truncated tree \(\hat{T}\) has size \(O({1 \over \varepsilon } \log u)\) in expectation.

Proof

Recall that only nodes whose estimated frequency is above \(\eta \varepsilon n\) are added to the truncated tree \(\hat{T}\), where \(\eta \) is a constant. We classify these nodes into heavy nodes and non-heavy nodes. A node is heavy if its true frequency is above \({1\over 2}\eta \varepsilon n\), and non-heavy otherwise.

We first observe that, on any level of the dyadic structure, there are at most \({n \over {1\over 2} \eta \varepsilon n} = O(1/\varepsilon )\) heavy nodes, treating \(\eta \) as a (small) constant. So even if they are all added to \(\hat{T}\), there are only \(O({1\over \varepsilon } \log u)\) of them in total.

For a non-heavy node v, the Count-Sketch may overestimate its frequency to above \(\eta \varepsilon n\) with a constant probability, say 1 / 4. Note that there can be \(\varTheta (u)\) non-heavy nodes being overestimated in expectation, but we will argue below that only \(O({1\over \varepsilon }\log u)\) will be added during the top-down construction of \(\hat{T}\).

Note that for a non-heavy node to be added, its parent must be a heavy node or another non-heavy node that has been overestimated. Thus, all the non-heavy nodes in \(\hat{T}\) make up a number of subtrees, where the root of each subtree must be a child of a heavy node. Let \({\mathbf {t}}\) be any such non-heavy subtree, and we will bound \({\mathrm {E}}[|{\mathbf {t}}|]\). Let r be the root of \({\mathbf {t}}\). For any node v below r, let \(I_v\) be an indicator variable where \(I_v=1\) if \(v\in {\mathbf {t}}\) and \(I_v=0\) otherwise. Let d(v) be the depth of v in \({\mathbf {t}}\), and we define \(d(r) = 0\). For v to be added to \({\mathbf {t}}\), all of its d(v) ancestors must have been added, so \({\mathrm {E}}[I_v] \le (1/4)^{d(v)+1}\) since each level uses an independent Count-Sketch. We then have

Finally, we observe that at most two such \({\mathbf {t}}\)’s can be attached to a heavy node, and there are only \(O({1\over \varepsilon }\log u)\) heavy nodes, so we conclude that there are \(O({1\over \varepsilon }\log u)\) non-heavy nodes in \(\hat{T}\) in expectation. \(\square \)

1.2 Constraints on the \(x^*_i\)’s

Lemma 2

Let \(x^*\) be the BLUE of (4). Let \(\lambda _v\) be any solution to (5). For any leaf w, let \(Z_w= \lambda _w \sum _{z\in \text {anc}(w)\setminus r}y_z/\sigma ^2_z\), and for any internal node v, \(Z_v= \sum _{w\prec v}\lambda _w Z_w\). For any node v except the root r, let \(F_v = \sum _{w\in \text {anc}(v)\setminus \{r\}}x^*_w/\sigma ^2_w\). Let \(\varDelta = (Z_r-y_r\pi _{r\text {'s left child}})/\lambda _r\). Then, we have

Proof

We will follow the method of Lagrange multiplier [19] to find the BLUE of (4). Since only the subtree root is known exactly, we introduce a single Lagrange multiplier \(\eta \). We set \(\sigma ^2_r=1/\eta \) instead of 0, and will later take the limit as \(\eta \) goes to \(\infty \). Denote by \({\text {diag}}(1/\sigma _v)\) the diagonal matrix with \(1/\sigma _v\) at entry (v, v). We further define \(Z={\text {diag}}(1/\sigma _v){y}\) and \(U={\text {diag}}(1/\sigma _v)A\). Then, the Lagrange function can be rewritten as \((Z-Ux^*)^T(Z-Ux^*)\). By differentiation, we can derive (i) \(y_r=x_r^*\); and (ii)

This is sufficient to define a solution; we can solve for \(x^*\) by premultiplying by \((U^TU)^{-1}\), at the cost of a computing a matrix inverse. In the following, to derive the equations stated in the lemma from (7), which leads to a much more efficient algorithm to compute \(x^*\).

Let \(\text {anc}(u,v)=\text {anc}(u)\cap \text {anc}(v)\). Let u be any node of \(T_r\). We also use \([\tau ]\) denote the \(\tau \) leaves of \(T_r\). Then, by simple calculation, we can see that \((U^TU)_{u,w}=\sum _{v\in \text {anc}(u,w)}\sigma _v^{-2}\), and \((U^TZ)_u=\sum _{v\in \text {anc}(u)}y_v/\sigma ^2_v\).

First, we take the weighted sum of corresponding rows on the LHS of (7) to obtain \(\sum _{u\prec v}\lambda _u(U^TU)_ux^*=\)

Note that in the last line of (8), the second component can be written as

We can also derive that the first component is

To see this, let us assume that this holds for any descendant of v. Then, we can derive \(\sum _{u\prec v}\sum _{w\in \text {anc}(u)\setminus \text {anc}(v)}\frac{\lambda _ux^*_w}{\sigma ^2_w}\)

Combining the above two results, we have

Secondly, we take the weighted sum of corresponding rows on the RHS of (7) to obtain

Finally, by combing (9) and (10), we derive that \(\forall v\),

Substituting v by r in (11), we can derive

We already have \(x^*_r=y_r\). Then, we can conclude that either \(\eta =+\infty \) or \(y_r\pi _r=Z_r\). As we know, \(y_r\) is given and irrelevant to \(\pi _rZ_r\) which implies \(\eta =+\infty \). To handle this infinity, we first express

where s is a child of root r. Now we take limit of \(\varDelta (\eta )\) and derive that

Finally by taking limit on (11) for any \(v\ne r\), we derive \(x^*_v=(Z_v-\lambda _v F_{\text {parent}(v)}-\lambda _v\varDelta )/\pi _v \), and the lemma is proved. \(\square \)

1.3 Covariance analysis

Lemma 3

Suppose we build a Count-Sketch with \(d=1\) row and w columns on a vector x. For any two elements \(u\ne v\), let \(y_u\) and \(y_v\) be the estimators for \(x_u\) and \(x_v\). Then, \({\mathrm {Cov}}(y_u,y_v) = x_u x_v/w\).

Proof

Recall that in the Count-Sketch, for any element v, the estimator is computed as \(y_v = g(v)\sum _{z}g(z)x_zI_v(z)\), where \(I_v(z)=1\) if \(h(v)=h(z)\) and 0 otherwise. So we have

Let \(f(i,j)={\mathrm {E}}[g(u)g(v)g(i)g(j)]x_ix_j{\mathrm {E}}[I_u(i)I_v(j)]\). We find that \(f(u,v)=x_ux_v\) and \(f(v,u)= x_u x_v/w\) (note, the last term in f(i, j) is not symmetric in i and j). If \(\{i,j\}\ne \{u,v\}\), then \(f(i,j)=0\), since \(g(\cdot )\) is 4-wise independent hash function. Therefore, we derive that \({\mathrm {Cov}}(y_u,y_v)=x_u x_v/w\). \(\square \)

On the other hand, prior analysis on the Count-Sketch shows that \({\mathrm {Var}}(y_v) = {1\over w} \sum _{i} x_i^2\). Thus, the covariance is usually order of magnitude smaller than the variance.

Reprints and permissions