基数统计—— HyperLogLog 算法

文章主体内容来自于 神奇的 HyperLogLog 算法，原创链接貌似已失效，可参照大概是其转载内容 HyperLogLog ，本文在此基础上略有删改。

基数计数基本概念

**基数计数(cardinality counting)**通常用来统计一个集合中不重复的元素个数，例如统计某个网站的 UV，或者用户搜索网站的关键词数量。数据分析、网络监控及数据库优化等领域都会涉及到基数计数的需求。 要实现基数计数，最简单的做法是记录集合中所有不重复的元素集合 $S_u$，当新来一个元素 $x_i$ ， 若 $S_u$ 中不包含元素 $x_i$，则将 $x_i$ 加入 $S_u$，否则不加入，计数值就是 $S_u$ 的元素数量。这种做法存在两个问题：

1. 当统计的数据量变大时，相应的存储内存也会线性增长
2. 当集合 $S_u$ 变大，判断其是否包含新加入元素 $x_i$ 的成本变大

基数计数方法

B树

B 树最大的优势是插入和查找效率很高，如果用 B 树存储要统计的数据，可以快速判断新来的数据是否已经存在，并快速将元素插入 B 树。要计算基数值，只需要计算 B 树的节点个数。 将 B 树结构维护到内存中，可以快速统计和计算，但依然存在问题，B 树结构只是加快了查找和插入效率，并没有节省存储内存。例如要同时统计几万个链接的 UV，每个链接的访问量都很大，如果把这些数据都维护到内存中，实在是够呛。

bitmap

bitmap 可以理解为通过一个 bit 数组来存储特定数据的一种数据结构，每一个 bit 位都能独立包含信息，bit 是数据的最小存储单位，因此能大量节省空间，也可以将整个 bit 数据一次性加载到内存计算。 如果定义一个很大的 bit 数组，基数统计中每一个元素对应到 bit 数组的其中一位，例如 bit 数组 001101001代表实际数组 [2,3,5,8]。新加入一个元素，只需要将已有的 bit 数组和新加入的数字做按位或 (or) 计算。bitmap 中 1 的数量就是集合的基数值。

bitmap 有一个很明显的优势是可以轻松合并多个统计结果，只需要对多个结果求异或就可以。也可以大大减少存储内存，可以做个简单的计算，如果要统计 1 亿个数据的基数值，大约需要内存： $100000000/8/1024/1024\approx 12M$，如果用 32bit 的 int 代表每个统计数据，大约需要内存：
$32\ast 100000000/8/1024/1024\approx 381M$。

bitmap 对于内存的节约量是显而易见的，但还是不够。统计一个对象的基数值需要 12M，如果统计10000 个对象，就需要将近 120G 了，同样不能广泛用于大数据场景。

概率算法

实际上目前还没有发现更好的在大数据场景中准确计算基数的高效算法，因此在不追求绝对准确的情况下，使用概率算法算是一个不错的解决方案。概率算法不直接存储数据集合本身，而是通过一定的概率统计方法预估基数值，这种方法可以大大节省内存，同时保证误差控制在一定范围内。目前用于基数计数的概率算法包括:

• Linear Counting(LC)：早期的基数估计算法，LC 在空间复杂度方面并不算优秀，实际上 LC 的空间复杂度与简单 bitmap 方法是一样的（但是有个常数项级别的降低），都是 $O(N_{\max })$ ;

• LogLog Counting(LLC)：LogLog Counting 相比于 LC 更加节省内存，空间复杂度只有$O(\log (\log (N_{\max })))$

• HyperLogLog Counting(HLL)：HyperLogLog Counting 是基于 LLC 的优化和改进，在同样空间复杂度情况下，能够比 LLC 的基数估计误差更小。

几种基数计数方法，统计 1 个页面，数据量达到 1 亿条 ip 时，所占用的空间及时间复杂度对比：

	B树	Bitmap	HyperLogLog 算法
占用空间大小	32x12 ≈ 384MB	12MB	12KB
时间复杂度	$O(\log N)$	$O(N)$	$O(1)$

HLL

直观演示

HLLDEMO

HLL 的实际步骤

待阅读后续讲解内容后反过头重新阅读这里，效果更佳

• Value: 随机生成的值（介于 1 和 16777215 之间）；
• Hash Value: 利用 Jenkins 函数生成 value 的哈希值；
• Bit String: 用二进制表示的该哈希值的 24 位低位内容。如图中 bin(3868256379)=‘0b11100110100100001110100001111011’；
• Register Index: 最低 $b$ 位用于确定要更新其值的寄存器的索引($m=2^b$）
• Register Value: 从剩下的内容中判定前导零的个数（从右往左，从低位到高位数），根据 HLL 论文中，前导零的个数需要 +1，如上图中，绿色部分的最低位是 “1”，则前导零的个数为 0，再加 1 就是 1。二进制串 “01010101001000” 从最低位数有 3 个 0，则前导 0 的个数就是 3，再加 1 就是 4。换句话说就是找第一个 “1” 出现的位置；
• Register Values: $ m=64$ 个寄存器值。索引 0 在左上角，索引 63 在右下角。要更新的寄存器以红色突出显示。请注意，只有 Bit String 中的值大于当前值时，才会更新该值；
• m: 桶或寄存器值的数目（$m=2^b$）
• Actual Cardinality: The cardinality of the set of all inserted values (i.e. the count of distinct values);集合中所有插入值的基数（即不同值的计数）；
• Algorithm: 所使用的基数估计算法；
• Estimated Cardinality: 算法估计的集合的基数；
• % Error: 估计基数与实际基数之差除以实际基数值；
• Line Graph: $x$ 轴表示每个算法的实际基数，$y$ 轴表示误差百分比。给出了 HLL 的初始 ball 和 bins 估计。蓝色区域是理论 HLL 误差界限（$1.04/\sqrt{(m)}$）；Plots the actual cardinality on the x-axis and the % Error on the y-axis for the each algorithm. The initi
• Histogram: 理论和观测寄存器值的分布，以及观测分布的算术、几何和调和平均值。请注意，HLL 中使用的调和平均值为 $1/2^M$，而这里只是 $1/M$，仅用于说明目的；

算法来源（N次伯努利过程）

下面非正式的从直观角度描述 LLC 算法的思想来源。

设 $a$ 为待估集合（哈希后）中的一个元素，由上面对 H 的定义可知，$a$ 可以看做一个长度固定的比特串（也就是 $a$ 的二进制表示），设 H 哈希后的结果长度为 L 比特，我们将这 $L$ 个比特位从左到右分别编号为 1、2、…、L：

假设哈希函数后的值 $a$，其每个比特位服从如下分布且相互独立。

$$p(x=k)=\{\begin{array}{ll}0.5,& k=0\\ 0.5,& k=1\end{array}$$

通俗说就是 $a$ 的每个比特位为 0 和 1 的概率各为 0.5，且相互之间是独立的。

设 $\rho (a)$ 为 $a$ 的比特串中第一个 “1” 出现的位置，显然 $1\le \rho (a)\le L$，这里我们忽略比特串全为 0 的情况（概率为 $1/2^L$）。如果我们遍历集合中所有元素的比特串，取 ${\rho }_{\max }$ 为所有 $\rho (a)$ 的最大值。

此时我们可以将 $2^{{\rho }_{\max }}$ 作为基数的一个粗糙估计，即：

$$n=2^{{\rho }_{\max }}$$

解释

注意如下事实：

由于比特串每个比特都独立且服从 0-1 分布，因此从左到右扫描上述某个比特串寻找第一个 “1” 的过程从统计学角度看是一个伯努利过程，例如，可以等价看作不断投掷一个硬币（每次投掷正反面概率皆为0.5），直到得到一个正面的过程。对于 $n$ 次伯努利过程，我们可以得到 $n$ 个首次出现正面的投掷次数$k_1,k2,\cdots ,k_n$，其中最大值记为 $k_{\max }$，那么可以考虑如下两个问题：

1、进行 $n$ 次伯努利过程，所有投掷次数都不大于 $k_{\max }$ 的概率是多少？

2、进行 $n$ 次伯努利过程，至少有一次投掷次数等于 $k_{\max }$ 的概率是多少？

首先看第一个问题，在一次伯努利过程中，投掷次数大于 $k$ 的概率 $1/2^k$，即连续掷出 $k$ 个反面的概率。因此，在一次过程中投掷次数不大于 $k$ 的概率为 $1-2^k$。因此，$n$ 次伯努利过程投掷次数均不大于 $k$ 的概率为：

$$P_n(X\le k_{\max })=(1-1/2^{k_{\max }}{)}^n$$

显然第二个问题的答案是：

$$P_n(X\ge k_{\max })=1-(1-1/2^{k_{\max }-1}{)}^n$$

从以上分析可以看出，当 $n\ll 2^{k_{\max }}$ 时，$P_n(X\ge k_{\max })$ 的概率几乎为 0，同时，当 $n\gg 2^{k_{\max }}$ 时，$P_n(X\le k_{\max })$ 的概率也几乎为 0。用自然语言概括上述结论就是：当伯努利过程次数远远小于 $2^{k_{\max }}$ 时，至少有一次过程投掷次数等于 $k_{\max }$ 的概率几乎为 0；当伯努利过程次数远远大于 $2^{k_{\max }}$ 时，没有一次过程投掷次数大于 $k_{\max }$ 的概率也几乎为 0。

以上结论可以总结为：进行了 $n$ 次进行抛硬币实验，每次分别记录下第一次抛到正面的抛掷次数 $k$，那么可以用 $n$ 次实验中最大的抛掷次数 $k_{\max }$ 来预估实验组数量 $n:\hat{n}=2^{k_{\max }}$

LogLogCounting

均匀随机化

与 LC 一样，在使用 LLC 之前需要选取一个哈希函数 H 应用于所有元素，然后对哈希值进行基数估计。H 必须满足如下条件（定性的）：

1. H 的结果具有很好的均匀性，也就是说无论原始集合元素的值分布如何，其哈希结果的值几乎服从均匀分布（完全服从均匀分布是不可能的，D. Knuth 已经证明不可能通过一个哈希函数将一组不服从均匀分布的数据映射为绝对均匀分布，但是很多哈希函数可以生成几乎服从均匀分布的结果，这里我们忽略这种理论上的差异，认为哈希结果就是服从均匀分布）。
2. H 的碰撞几乎可以忽略不计。也就是说我们认为对于不同的原始值，其哈希结果相同的概率非常小以至于可以忽略不计。
3. H 的哈希结果是固定长度的。

以上对哈希函数的要求是随机化和后续概率分析的基础。后面的分析均认为是针对哈希后的均匀分布数据进行。

分桶平均

上述分析给出了 LLC 的基本思想，不过如果直接使用上面的单一估计量进行基数估计会由于偶然性而存在较大误差。因此，LLC 采用了分桶平均的思想来消减误差。具体来说，就是将哈希空间平均分成 m 份，每份称之为一个桶（bucket）。对于每一个元素，其哈希值的前 k 比特作为桶编号，其中 $2^k=m$，而后 $L-k$ 个比特作为真正用于基数估计的比特串。桶编号相同的元素被分配到同一个桶，在进行基数估计时，首先计算每个桶内元素最大的第一个 “1” 的位置，设为 $M[i]$，然后对这 $m$ 个值取平均后再进行估计，即：

$$\hat{n}=2^{\frac{\sum M[i]}{m}}$$

这相当于物理试验中经常使用的多次试验取平均的做法，可以有效消减因偶然性带来的误差。下面举一个例子说明分桶平均怎么做。

假设 H 的哈希长度为 16bit，分桶数 $m$ 定为32。设一个元素哈希值的比特串为 “0001001010001010”，由于 $m$ 为32，因此前 5 个 bit 为桶编号，所以这个元素应该归入 “00010” 即 2 号桶（桶编号从 0 开始，最大编号为 $m-1$），而剩下部分是 “01010001010” 且显然 $\rho (01010001010)=2$，所以桶编号为“00010” 的元素最大的 $\rho$ 即为 $M[2]$ 的值。

偏差修正

上述经过分桶平均后的估计量看似已经很不错了，不过通过数学分析可以知道这并不是基数 $n$ 的无偏估计。因此需要修正成无偏估计。这部分的具体数学分析在 “Loglog Counting of Large Cardinalities” 中，过程过于艰涩这里不再具体详述，有兴趣的朋友可以参考原论文。这里只简要提一下分析框架：

首先上文已经得出：

$$P_n(X\le k)=(1-1/2^k{)}^n$$

因此：

$$P_n(X=k)=(1-1/2^k{)}^n-(1-1/2^{k-1}{)}^n$$

这是一个未知通项公式的递推数列，研究这种问题的常用方法是使用生成函数（generating function）。通过运用指数生成函数和 poissonization 得到上述估计量的 Poisson 期望和方差为：

$$\begin{array}{l}{\mathcal{E}}_n\sim \left[{\left(\Gamma (-1/m)\frac{1-2^{1/m}}{\log 2}\right)}^m+{ϵ}_n\right]\cdot n\\ {\mathcal{V}}_n\sim \left[{\left(\Gamma (-2/m)\frac{1-2^{2/m}}{\log 2}\right)}^m-{\left(\Gamma (-1/m)\frac{1-2^{-1/m}}{\log 2}\right)}^{2m}+{\eta }_n\right]\cdot n^2\end{array}$$

其中 $|{ϵ}_n|$ 和 ${\eta }_n$ 不超过 $10^{-6}$。

最后通过 depoissonization 得到一个渐进无偏估计量：

$$\hat{n}={\alpha }_{m}m{2}^{\frac{1}{m}\sum M[i]}$$

其中：

$${\alpha }_m:={\left(\Gamma (-1/m)\frac{1-2^{1/m}}{\log 2}\right)}^{-m}$$

这里的

$$\Gamma (s):=\frac{1}{s}{\int }_0^{\infty }{e}^{-t}{t}^{s}dt$$

其中 $m$ 是分桶数。这就是 LLC 最终使用的估计量。

误差分析

不加证明给出如下结论：

$$\text{ Standard error }(\frac{\hat{n}}{n})\approx \frac{1.30}{\sqrt{m}}\text{ . }$$

算法应用

误差控制

在应用 LLC 时，主要需要考虑的是分桶数 $m$，而这个 $m$ 主要取决于误差。根据上面的误差分析，如果要将误差控制在 $ϵ$ 之内，则：

$$m>(1.30/ϵ{)}^2$$

内存使用分析

内存使用与 $m$ 的大小及哈希值得长度（或说基数上限）有关。假设 H 的值为 32bit，由于 ${\rho }_{\max }\le 32$，因此每个桶需要 5bit 空间存储这个桶的 ${\rho }_{\max }$，$m$ 个桶就是 $5\times m/8$ 字节。例如基数上限为一亿（约 $2^{2}7$），当分桶数 $m$ 为 1024 时，每个桶的基数上限约为 $2^{17}$，而 ${\log }_2({\log }_2(2^{17}))=4.09$ ，因此每个桶需要 5bit，需要字节数就是 $5\times 1024/8=640$，误差为 $1.30/\sqrt{1024}=0.040625$，也就是约为4%。

合并

与 LC 不同，LLC 的合并是以桶为单位而不是 bit 为单位，由于 LLC 只需记录桶的 ${\rho }_{\max }$，因此合并时取相同桶编号数值最大者为合并后此桶的数值即可。

HyperLogLog Counting

HyperLogLog Counting（以下简称 HLLC）的基本思想也是在 LLC 的基础上做改进，具体细节请参考“HyperLogLog: the analysis of a near-optimal cardinality estimation algorithm” 这篇论文。

基本算法

HLLC 的第一个改进是使用调和平均数替代几何平均数。注意 LLC 是对各个桶取算数平均数，而算数平均数最终被应用到 2 的指数上，所以总体来看 LLC 取得是几何平均数。由于几何平均数对于离群值（例如这里的 0）特别敏感，因此当存在离群值时，LLC 的偏差就会很大，这也从另一个角度解释了为什么 $n$ 不太大时 LLC 的效果不太好。这是因为 $n$ 较小时，可能存在较多空桶，而这些特殊的离群值强烈干扰了几何平均数的稳定性。

因此，HLLC使用调和平均数来代替几何平均数，调和平均数的定义如下：

$$H=\frac{n}{\frac{1}{x_1}+\frac{1}{x_2}+\dots +\frac{1}{x_n}}=\frac{n}{{\sum }_{i=1}^n\frac{1}{x_i}}$$

调和平均数可以有效抵抗离群值的扰动。使用调和平均数代替几何平均数后，估计公式变为如下：

$$\hat{n}=\frac{{\alpha }_m{m}^2}{\sum 2^{-M}}$$

其中：

$${\alpha }_m={\left(m{\int }_0^{\infty }{\left({\log }_2\left(\frac{2+u}{1+u}\right)\right)}^{m}du\right)}^{-1}$$

偏差分析

根据论文中的分析结论，与 LLC 一样 HLLC 是渐近无偏估计，且其渐近标准差为：

$$S{E}_{\text{hllc }}(\hat{n}/n)=1.04/\sqrt{m}$$

因此在存储空间相同的情况下，HLLC 比 LLC具有更高的精度。例如，对于分桶数 $m$为 $2^{1}3$（8k 字节）时，LLC 的标准误差为 1.4%，而 HLLC 为 1.1%。

分段偏差修正

在 HLLC 的论文中，作者在实现建议部分还给出了在 $n$ 相对于 $m$ 较小或较大时的偏差修正方案。具体来说，设 $E$ 为估计值：

当 $E\le \frac{5}{2}m$ 时，使用 LC 进行估计。

当 $\frac{5}{2}m<E\le \frac{1}{30}{2}^{32}$ 时，使用上面给出的 HLLC 公式进行估计。

当 $E>\frac{1}{30}{2}^{32}$ 时，估计公式如为 $\hat{n}=-2^{32}\log (1-E/2^{32})$。

关于分段偏差修正效果分析也可以在原论文中找到。

结论

并行化

这些基数估计算法的一个好处就是非常容易并行化。对于相同分桶数和相同哈希函数的情况，多台机器节点可以独立并行的执行这个算法；最后只要将各个节点计算的同一个桶的最大值做一个简单的合并就可以得到这个桶最终的值。而且这种并行计算的结果和单机计算结果是完全一致的，所需的额外消耗仅仅是小于 1k 的字节在不同节点间的传输。

应用场景

基数估计算法使用很少的资源给出数据集基数的一个良好估计，一般只要使用少于 1k 的空间存储状态。这个方法和数据本身的特征无关，而且可以高效的进行分布式并行计算。估计结果可以用于很多方面，例如流量监控（多少不同 IP 访问过一个服务器）以及数据库查询优化（例如我们是否需要排序和合并，或者是否需要构建哈希表）。

参考阅读

Redis new data structure: the HyperLogLog

HyperLogLog — Cornerstone of a Big Data Infrastructure