寻找近似Quantiles

给我写信
GitHub

算法

---

目录

1. 问题描述

---

有n个数据，如果这些数据全部在内存中，则可以在O(n)时间内从中找出某个quantile。

但是当数据量非常大的时，这些数据无法全部放入内存，怎么办？

当数据量很大时（例如PB级别的量），往往允许找出的quantile离真实值有一定的误差，但要在必要的误差限度内。

φ-quantile

把n个数字按照从小到大的顺序排列成一个有序序列L，rank值为φn的那个数据点就称为L的φ-quantile，这里 φ ∈ [0, 1]。

例如，L的中位数（median），就是L的0.5-quantile。

寻找φ-quantile的理想的算法应该有以下的优点：

1. 只需遍历一遍数据；
2. 是确定的(deterministic)算法；
3. 算出的结果，与真实结果之间的误差足够小；
4. 不需要事先知道关于数据集的分布特性；
5. 允许并行化运行该算法；
6. 除了找到第一个quantile之外，每计算一个新的quantile，所需的时间和空间开销应该是一个常量。

2. 思路

---

笨方法

排序！ 将所有的数据排序后，自然就得到了想要的quantiles。这对于数据量较小时很容易实现，代价也不大。但如果是10 trillion的数据，你排序试试看要多久。

巧方法

本文描述的OPAQ方法来自于论文 A One-Pass Algorithm for Accurately Estimating Quantiles for Disk-Resident Data，这篇论文描述了怎样通过将数据集划分为若干个等长子集的方式来寻找approximate quantils。

而论文 Approximating quantiles in very large dataset 则描述了如果原数据集被划分为不等长子集，怎样寻找approximate quantiles。

3. OPAQ算法描述

---

OPAQ是“One Pass Algorithm for Quantiles”的缩写。

首先定义一些概念。

Term	Description
m	size of each run
n	total number of elements
r	number of runs ($r = \frac{n}{m}$)
s	size of the sample for each run
φ	quantile fraction (φ ∈ [0, 1])
α	index (rank) of the quantile ($α = {φ}\times{n}$)
$e_α$	value of the quantile

OPAQ算法由2个阶段组成：1) 取样阶段（sampling phase）; 2) quantile计算阶段（quantile finding phase）。

在取样阶段，将全部数据划分为$r$个子集（称为$r$个runs）。对于每一个run，我们要为其生成$s$个sample points，即$S = [s_1, s_2, ..., s_s]$，且满足$s_i <= s_{i+1}$ 。随后，这$r$个sample list将被合并为一个长度为$rs$升序列表，该列表将在quantile计算阶段中被用于估计真实的φ-quantile的上界（upper bound）和下界（lower bound）。

OPAQ计算出的结果的准确度依赖于上述这两个阶段。

3.1 The Sampling Phase

为了对$e_φ$作出估计，我们要计算出它的lower bound $e_φ^l$以及它的upper bound $e_φ^u$，这里$e_φ ∈ [e_φ^l, e_φ^u]$，并确保$[e_φ^l, e_φ^u]$内的点的数量是bounded。

怎样为长度为$m$的某个$run$生成$s$个sample points呢？我们使用论文 On the Versatility of Parallel Sorting by Regular Sampling 中提到的regular sampling方法：

假定目标run的长度为$m$，其中每一个数据点的index从1开始，到$m$结束。将该run划分为$s$个等长的sub-run，为第i个sub-run找一个sample point $s_i$，使得$s_i$大于等于第i个sub-run中任何数据，且$s_{i}$小于等于第i+1个sub-run中任何数据。实际上，这个$s_i$可以是这个run中的数据，也可以不是这个run中的数据，但只要满足上面的条件，$s_i$即可代表第i个sub-run。

例如，如果两个sub-list为$l_1$ = {3, 5, 7, 7}以及$l_2$ = {8, 9, 10, 11}，那么则可以选择$s_1$ = 7, $s_2$ = 11；也可以选择$s_1$ = 7.5, $s_2$ = 11; 还可以选择$s_1$ = 8, $s_2 = 11.5$。

一个典型的做法是：算出list的中位数，然后根据这个中位数将该list划分为两个等长的sub-list（记为$l_1$和$l_2$），使得$l_1$中的数据都小于或等于$l_2$中的数据，分别取出其中的最大值作为两个sample point。对$l_1$和$l_2$继续按照这种方法划分出新的更短的sub-list，直到新的sub-list的长度达到了$\frac{m}{s}$。这样的迭代过程最多有$\log s$次。在最坏的情况下，我们可以在O($m\log s$)时间内找到$s$个sample points。

例子：

某个$run$中的数据点为： {7， 15， 2， 0， 4， 12， 13， 20， 0， 1， 5， 6， 2， 14， 16， 9}，要为其找出4个sample points，因此，m = 16, s = 4, m/s = 16/4 = 4。

最终生成的4个sub-run为{2, 0, 0, 1}, {4, 5, 6, 2}, {7, 12, 13, 9}, {15, 20, 14, 16}，对应的sample points为{$s_1$ = 2, $s_2$ = 6, $s_3$ = 13, $s_4$ = 20}。

生成sample points的过程如下：

在找出了全部$r$个sample lists之后，我们将它们合并起来，形成一个长度为$rs$的升序列表，记为List。

3.2 The Quantile Finding Phase

现在，我们将利用上面得到的List来计算$e_α^u$和$e_α^l$。对于List中的任意一个sample point $s_i$，它有以下的性质（证明过程可以参看后面的附录）：

1. 至少存在 $i\frac{m}{s}$ 个小于或者等于$s_i$的数据点；
2. 另外，至多存在 ${(r-1)}{(\frac{m}{s}-1)}$ 个小于$s_i$的数据点。

因此，最多可能存在有 $i\frac{m}{s} + {(r-1)}{(\frac{m}{s}-1)}$ 个小于$s_i$的数据点（这是当前i个sub-run中的数据均小于$s_i$时的情况）。

我们将利用上面的两条性质来定位$e_α^l$和$e_α^u$，并计算真实quantile与估计值之间的误差范围。

定位$e_α^l$以及$e_α^u$

首先，我们让$e_α^l$等于List中的第i个sample point，并希望它能满足：

$$i\frac{m}{s} + (r-1)(\frac{m}{s}-1) ≤ α < (i+1)\frac{m}{s} + (r-1)(\frac{m}{s}-1) \tag 1$$
之所以将给出上面这个条件，是因为我们希望找到的$e_α^l$尽可能地靠近List的大端。

解这个关于i的不等式，可以得到：

$$\frac{s}{m}[α-(r-1)(\frac{m}{s}-1))]-1 < i ≤ \frac{s}{m}[α-(r-1)(\frac{m}{s}-1))] \tag 2$$

即

$$i = \lfloor α\frac{s}{m} - (r-1)(1-\frac{s}{m})) \rfloor \tag 3$$

因此

$$e_α^l = List[\lfloor α\frac{s}{m} - (r-1)(1-\frac{s}{m})) \rfloor] \tag 4$$

我们再让$e_α^u$等于List中的第j个sample point，并希望它能满足:

$$(j-1)\frac{m}{s} < α ≤ j\frac{m}{s}$$
之所以给出上面的条件，是因为我们希望找出的$e_α^u$尽可能靠近List的小端。

解这个关于j的不等式，可以得到：

$$α\frac{s}{m} ≤ j < 1 + α\frac{s}{m} \tag 5$$

即

$$j = \lceil α\frac{s}{m} \rceil \tag 6$$

因此

$$e_α^u = List[\lceil α\frac{s}{m} \rceil] \tag 7$$

误差估算

现在我们要计算，$e_α^l$与$e_α$之间的最大范围，$e_α$与$e_α^u$之间的最大范围，以及$e_α^l$与$e_α^u$之间的最大范围。

Lemma 1: $e_α^l$与$e_α$之间最多存在$\frac{n}{s}$个数据点

证明
用$N_L$表示$e_α^l$与$e_α$之间的数据点数量，用$N_{min}({Cond})$表示满足条件Cond的数据点数量的下限。那么，我们有

$$N_L ≤ α - N_{min}{(X < e_α^l)}$$
这里的$X$代表原始数据全集中的任意一个数据点。

显然，

$$N_{min}{(X < e_α^l)} = i\frac{m}{s}$$

用(3)来替换这里的i，可以得到:

$$N_{min}{(X < e_α^l)} = \lfloor α\frac{s}{m} - (r-1)(1-\frac{s}{m})) \rfloor \frac{m}{s}$$

进而

$$\begin{aligned} N_L & ≤ α - \lfloor α\frac{s}{m} - (r-1)(1-\frac{s}{m})) \rfloor \frac{m}{s} \\ & ≤ α - [α\frac{s}{m} - (r-1)(1-\frac{s}{m})-1]\frac{m}{s} \\ & = \frac{n}{s} + 1 - r \\ & ≤ \frac{n}{s} \end{aligned}$$

Lemma 2: $e_α$与$e_α^u$之间最多存在$\frac{n}{s}$个数据点
同理可以证明

Lemma 3: $e_α^l$与$e_α^u$之间最多存在$2\frac{n}{s}$个数据点
非常明显

因此，只要每个run中的取出的sample points数量足够大（与原始数据集的大小n相比），quantile的估算值与真实值之间的差异就可以足够小。

3.3 复杂度

Phase	Complexity
Finding the $rs$ sample points	$O(rm\log{s})$
Mergeing the $r$ sample lists	$O(rs\log{r})$
Estimating $q$ quantiles	$O(q)$

实际上，OPAQ算法最复杂的部分在于为每一个subset计算出s个sample points，并将其合并成一个有序的列表。一旦这一步完成，我们可以立即得到想要的approximate $α$-quantile。

并且，在计算出了List之后，它可以被重复地用来计算其它任意的quantiles。

3.4 并行化的OPAQ

如果原始数据集（均匀地）分布在多个节点上，那么可以为每一个节点上的数据子集生成一个locally sorted sample list，再把他们合并成一个globally sorted sample list，然后再计算出α-quantile。

这里的关键在于，怎样高效地将多个locally sorted sample lists合并成一个globally sorted sample list？

3.5 处理增量数据

如果每一个run保留了原来的sorted sample list，那么，当有新的增量数据到达某个run时，该run只需要为增量数据计算出对应的sorted sample list，并将它与原来的sorted sample list进行合并即可。

3.6 不足之处

实际上，我在应用OPAQ算法时，对它并不怎么满意，因为它有以下不足:

1. 精确度受到单节点运算能力的限制

$e_α^l$与$e_α$之间最多存在$\frac{n}{s}$个数据点， $e_α^l$与$e_α^u$之间最多存在$2\frac{n}{s}$个数据点。

$s$是一个run中的sample points的数量，最多也就等于这个run的长度，一般一个run的数据都在一个节点上。因此，$s$会受到单节点的计算能力的限制，而$\frac{n}{s}$中并不含有变量$r$，这就意味着无法通过增加run的数量来提高结果的精度；相反，如果只增加run的数量，那么$n$会增大，误差反而会增大。

假设有10亿条数据点，想分为100个run进行计算，那么每个run中有1000万条数据，从每个run中选取出10万个sample points，那么quantile的偏差范围在 $2 \times \frac{10亿}{10万} = 2000$，并不是那么小。

2. 必须满足各个run中的数据点数量相同这个前提

这个条件比较苛刻，如果要满足这个条件，往往需要对数据进行shuffle。实际上，论文 Approximating quantiles in very large datasets 提出了另一个算法，它基于OPAQ算法，但是不要求每个run中的数据点的数量相同。