基数统计:从Linear Counting到Hyper LogLog

欢迎来我的小站转转~~

应用场景

基数统计(Cardinality Counting)指计算一个数据集中不同元素的数量，在很多场景都需要这样的功能：电商场景中的独立用户(UV)数量统计、数据库中快速计算字段取值数量以优化query、计算与某个站点相关的不同链接数量等。

为实现这样的功能，最直观的方法就是使用 Set 记录每一个元素，在数量巨大的情况下，会造成大量的内存占用。以存储 $2^{30}$ (约10亿)个长度为64bit的url为例，需要占用 $64\ast 2^{30}=64GB$ 的空间。Set方案的空间占用与元素的数量成正比，空间复杂度为 $O(N)$

如果使用bitmap来存储，则需要将不同的url映射到不同的数字。即使最理想的情况下，url之间没有冲突，且从1到 $2^{30}$ 都没有空闲，bitmap也需要占用 $2^{30}÷8=2^{27}=128MB$ 的空间。bitmap方案占用空间只与最大取值有关系，与元素数量无关。极端情况下，即使只有一个url，在这个例子中也需要占用128MB的空间，空间复杂度为 $O(N_{max})$

以上两种方式都对基数进行精确统计，但是很多时候只需要一个大概的数量，只要误差在一定范围内即可。接下来要讨论的几种算法可以在牺牲一定精确性的情况下，大幅度降低空间占用。

Linear Counting

诞生于1990年的论文“A linear-time probabilistic counting algorithm for database applications“，相比bitmap有常数项的空间复杂度优化，但整体还是 $O(N_{max})$ 。

主要步骤

1. 数据哈希：假设原始数据的基数为n，使用一组哈希空间为m的哈希函数H，将原始数据均匀分组。
2. 分桶数据统计：构建一个长度为m的bitmap，其中每一个bit对应哈希空间的一个值。生成哈希数组的值如果存在，则把相应的bit设置为1。
3. 数量估计：当所有值设置完成后，统计bitmap中为0的bit数为u。则估计的基数为 $\hat{n}=-m\log \frac{u}{m}$。

基数公式推导

记事件 $A_j$ 为“经过n个数字哈希之后，桶j仍然为空”，则 $$P(A_j)=(1-\frac{1}{m}{)}^n$$

记 $U_n$ 为经过n个数字哈希之后，仍然为空的桶数量。因为各个桶之间互相独立，所以 $E(U_n)=m(1-\frac{1}{m}{)}^n=m[(1+\frac{1}{-m}{)}^{-m}{]}^{-\frac{n}{m}}$

当m趋向于无穷大时，${\lim }_{m\to \infty }E(U_n)={\lim }_{m\to \infty }m[(1+\frac{1}{-m}{)}^{-m}{]}^{-\frac{n}{m}}=m{e}^{-\frac{n}{m}}$

得到 $n=-m\log \frac{E(U_n)}{m}$

由于每个桶的取值都服从参数相同的0-1分布，所以$U_n$符合二项分布。当n很大的时候，可以用正态分布逼近二项分布，$U_n$ 的概率密度函数为：

$f(x)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}$

因此我们实际观察到的空桶数量 $u$ 可以作为 $E(U_n)$ 的无偏估计，由于 $=-m\log \frac{x}{m}$ 是可逆函数，可以把 $\hat{n}$ 作为理论 $n$ 的无偏估计。

误差估计

偏差

记装载因子 $t=\frac{n}{m}$ , 空桶比例 $V_n=\frac{U_n}{m}$ , 空桶数量期望 $p=E(V_n)=\frac{E(U_n)}{m}=e^{-\frac{n}{m}}=e^{-t}$

记$-\ln (V_n)=f(V_n)$，则$\hat{n}=m\ast (-\ln \frac{u}{m})=mf(V_n)$

在 $V_n$ 处进行泰勒展式:

$$\hat{n}=mf(V_n)=m(f(p)+(V_n-p)f^{\prime }(p)+\frac{1}{2}(V_n-p{)}^2{f}^{\prime \prime }(p)+\frac{1}{6}(V_n-p{)}^3{f}^{\prime \prime \prime }(p)....)$$

$$=m(t-\frac{V_n-p}{p}+\frac{(V_n-p{)}^2}{2p^2}-\frac{(V_n-p{)}^3}{3p^3}+\frac{(V_n-p{)}^4}{4p^4}....)$$

取前三项:

$$E(\hat{n})=mt-\frac{mE(V_n-p)}{p}+\frac{E(V_n-p{)}^2}{2p^2}$$

由于 $p=E(V_n)$ , $E(V_n)-p=0$

所以可以化简为
$$E(\hat{n})=mt+\frac{m}{2p^2}E(V_n-p{)}^2$$

其中
$$E(V_n-p{)}^2=Var(V_n)=(\frac{1}{m})e^{-t}(1-(1+t)e^{-t})$$

$Var(V_n)$ 通过 $Var(X)=E(X^2)-E(X{)}^2$ 得到，计算的具体过程可以见论文Appendix A。

最终 $$E(\hat{n})=n+\frac{e^t-t-1}{2}$$

$$Bias(\frac{\hat{n}}{n})=E(\frac{\hat{n}}{n})-1=\frac{e^t-t-1}{2n}$$

从中得到一个显然的结论：m越大，偏差越小。

标准差

取泰勒展式的前两项， $\hat{n}=m(t-\frac{V_n-p}{p})$

$$Var(\frac{\hat{n}}{n})=\frac{1}{n^2}(\frac{m^2}{p^2}Var(V_n-p))=\frac{1}{n^2}(\frac{m^2}{p^2}Var(V_n))=\frac{1}{n^2}(\frac{m^2}{p^2}\frac{1}{m})p(1-(1+t)p)=\frac{m(e^t-t-1)}{n^2}$$

具体过程见论文的Appendix C，得 $StdError(\frac{\hat{n}}{n})=\frac{\sqrt{m(e^t-t-1)}}{n}$

其中 $t=\frac{n}{m}$ 。对m求导得到 $\frac{1}{m}(e^{\frac{n}{m}}(m-n)-1)$，由于m小于n，这是一个对m的递减函数，因此分桶数量越大，标准差越小。

长度m的选择

误差控制

如果我们需要将误差控制在 $ϵ$ 以下，即
$$\frac{\sqrt{m(e^t-t-1)}}{n}=\frac{\sqrt{\frac{e^t-t-1}{t}}}{t}<ϵ$$
得到m的最小数量为
$$m>\frac{e^t-t-1}{(ϵt{)}^2}$$

满桶控制

注意到，如果每个桶的都置为1，即 $u=0$ ，那么 $\hat{n}=-m\log \frac{u}{m}$ 就会趋向于无穷大，失去意义。但只要m<n，就总有概率导致 $u=0$，我们追求的是这个概率足够小

$n\to +\infty$，则 $U_n$ 的二项分布可近似为正态分布

可以通过让 $U_n$ 的期望相比0的偏移大于若干倍的 $U_n$ 标准差来控制：
$$E(U_n)-a\times StdDev(U_n)>0$$
代入 $E(U_n)=m{e}^{-t}$ 和 $StdDev(U_n)=(m{e}^{-t}(1-(1+t)e^{-t}){)}^{\frac{1}{2}}$，得到

$$m>a^2(e^t-t-1)$$

结合两个限制条件

结合误差控制和满桶控制，最终m在期望误差内的最小长度为:
$$m>\beta (e^t-t-1),\beta =max(a^2,\frac{1}{(ϵt{)}^2})$$

确定a

u满足二项分布，当n非常大而落进一个桶的概率 $p=\frac{1}{m}$ 非常小时，可以用泊松分布进行近似
$$\underset{n,m\to \infty }{\lim }Pr(U_n=k)=\frac{{\lambda }^k}{k!}{e}^{-\lambda }$$
因此，满桶的概率为
$$Pr(U_n=0)=e^{-\lambda }$$
如果要满桶的概率小于1%，只要 $e^{-\lambda }<0.01$，得到 $\lambda >5$，所以只要 $U_n$ 的期望与0偏离为 $\sqrt{5}$ 即可保证满桶的概率小于 1%。同理，可以将概率设置为其他值来求出对应的a，满桶概率越小需要的a越大。

空间复杂度

在 $a=\sqrt{5}$ 的情况下，作者对不同规模n下m的最小长度进行了计算，随着m和n的增大，m大约为n的十分之一。因此LC所需要的空间只有传统的bitmap直接映射方法的1/10，但整体空间复杂度仍为 $O(N_{max})$。结果如下表

仍以 $n=2^{30}$ 为例，通过计算机求解得到 $ϵ=0.01$时, $m=75402422$，因此空间占用为 $m/8/1024/1024=9$ MB，约为bitmap方案的 $9MB/128MB=\frac{1}{14}$ ，与表中的比例数量级相符。

并集与交集

只要哈希函数和m相同，Linear Counting支持对多个记录进行合并计算。
如果要计算两个集合的并集 $set\_size(R_1\cup R_2)$，可以简单进行 OR 运算得到新的bitmap，然后计算出估计结果。
如果要计算两个集合的交集，可以通过两集合和并集间接算出来： $set\_size(R_1\cap R_2)=set\_size(R_1)+set\_size(R_2)-set\_size(R_1\cup R_2)$

问题

Linear Counting在n过大时，会导致bitmap满桶，此时误差就会变为无穷大。

LogLog Counting

Linear Counting虽然有了常数级的优化，但空间复杂度仍然是 $O(N)$ 级别，而LogLog Counting却如其名只有 $O({\log }_2({\log }_2(N))$ 的空间复杂度。目前在处理大数据的基数问题时，采用的方法基本都是LogLog Counting的几个变种。

主要步骤

与Linear Counting一样，LogLog Counting 需要选取一个哈希函数H应用于所有元素，然后对哈希值进行基数估计。H必须满足如下条件：

1、具有很好的均匀性：无论原始值如何分布，其哈希结果的值都服从均匀分布（完全服从均匀分布是不可能的，但是很多哈希函数可以生成几乎服从均匀分布的结果）。

2、H的碰撞几乎可以忽略不计：对于不同的原始值，其哈希结果相同的概率非常小。

3、H的哈希结果是固定长度的。

设a是经过H哈希后的一个值，长度固定为L。因为经过H的哈希结果服从均匀分布，所以a中每个比特的取值为1的概率都是0.5，且相互独立。记 $\rho (a)$ 为a的比特串中第一个“1”出现的位置， ${\rho }_{max}$ 为所有元素哈希后 $\rho (a)$ 的最大值。

则对于总量的估计为 $\hat{n}=2^{{\rho }_{max}}$

粗糙估计

为什么可以这样估计呢？

由于比特串每个比特都独立且服从0-1分布，因此从左到右扫描上述某个比特串寻找第一个“1”的过程从统计学角度看是一个伯努利过程，例如，可以等价看作不断抛硬币直到正面的过程。这样的过程中投掷一次就得到正面的概率为$\frac{1}{2}$，投掷两次得到正面的概率是$\frac{1}{2^2}$，…，投掷k次才得到第一个正面的概率为 $\frac{1}{2^k}$。

现在如果进行n次伯努利过程，所有投掷次数都不大于k的概率和至少有一次投掷次数大于k的概率分别为:

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

$$P_n(X\ge k)=1-(1-\frac{1}{2^k}{)}^n$$

由此可见，当 $n\ll 2^k$ 时，$P_n(X\ge k)$ 几乎为零； $n\gg 2^k$ 时，$P_n(X\le k)$ 也几乎为零。从这里就可以衍生出一个结论：如果对n个元素进行哈希后，得到 ${\rho }_{m}ax$ ，则 $P(2^{\rho }\ll n)$ 和 $P(2^{\rho }\gg n)$ 都几乎为零， $2^{\rho }$ 可以作为n的一个粗糙估计。

分桶平均

但是上述的估计实在太粗糙了，因此 LogLog Counting使用分桶平均的方式，根据前k个bit，将哈希空间平分为 $m=2^k$ 份，每个桶都单独维护一个 ${\rho }_{max}$，形成一个数组 $M$。在进行哈希的时候，根据一个元素的前k个bit确定桶编号i，剩余的bit按照上述的步骤更新 $M[i]$ 的 ${\rho }_{max}$，最终估计时取所有桶的平均值再进行估计。

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

误差修正

分桶平均只是减少了误差，但是仍然基于粗糙的估计，因此论文对其进行了修正以达到无偏估计。

这部分需要大量的数学分析，只提一下分析的框架：

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

构建指数生成函数，通过poissonization得到估计量的 $Poisson$ 期望和方差后，通过de-poissonization 得到一个渐近无偏估计量：

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

其中

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

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

误差分析

不加证明得给出结论：

当 $m>64$ 时$StdError(\frac{\hat{n}}{n})\approx \frac{1.3}{\sqrt{m}}$

如果要将误差控制在 $ϵ$ 之内，则

$$m>(\frac{1.30}{ϵ}{)}^2$$

空间复杂度

内存使用与m的大小及哈希值长度有关，如果统计的基数大小为n，则${\rho }_{max}\ge O({\log }_2\frac{n}{m})$，对应保存 $M$ 的空间为 $O(m \times \log_2{\rho_{max}}) $，在m固定时，空间复杂度为 $O(lo{g}_2(lo{g}_{2}n))$
假设哈希后长度为64bit，则 ${\rho }_{max}\le 64$，因此每个桶需要 ${\log }_{2}64=6$ bit的空间来存储 ${\rho }_{max}$，分成m个桶即 $\frac{6\times m}{8}$ 字节。当分桶数m为1024时，误差为 $1.3÷\sqrt{1024}\approx 4\%$，需要的字节数就是 $6\times 1024÷8=768B$ ，比起Linear Counting的 9MB 有了5个数量级的下降。

并集与交集

LogLog Counting也支持计算两个集合的并集基数，区别在于LogLog Counting以桶为单位，而Linear Counting以bit为单位。由于LogLog Counting只需要记录桶的 ${\rho }_{max}$，所以合并的时候取相同桶编号中较大的 ${\rho }_{max}$即可。

问题

LogLog Counting比Linear Counting的节约了指数级的空间，但是它误差控制在 $\frac{1.3}{\sqrt{m}}$的前提是n远大于m。在n比较小的时候，它的估计误差就会很大。

Adaptive Counting

主要改进

Linear Counting和LogLog Counting的存储结构其实是兼容的，只要统计LogLog Counting的${\rho }_{max}$ 是否为0就可以退化为Linear Counting。
既然两者的存储结构兼容，而又各自擅长n较小和n较大的情况，Adaptive Counting就提出通过寻找一个分段点，让两者的优点互补。
Linear Counting和LogLog Counting的标准误差为
$$StdErro{r}_{linea{r}_{c}ounint}(\frac{\hat{n}}{n})=\frac{\sqrt{e^t-t-1}}{t\sqrt{m}}$$
$$StdErro{r}_{loglogcounting}(\frac{\hat{n}}{n})=\frac{1.3}{\sqrt{m}}$$

两式联立
$$\frac{\sqrt{e^t-t-1}}{t\sqrt{m}}=\frac{1.3}{\sqrt{m}}$$

即可得到相交点 $t\approx 2.89$，其中 $t=\frac{n}{m}$，此时Linear Counting的空桶率 $\beta =e^{-t}\approx 0.051$
所以Adaptive Counting中，估算公式为：

$$\hat{n}=\{\begin{array}{l}{\alpha }_{m}m{2}^{\frac{\Sigma M}{m}}\text{if}0\le \beta \le 0.051\\ -m\ln (\beta )\text{if}0.051\le \beta \le 1\end{array}$$

因为AC只是LC和LLC的简单组合，所以误差分析可以依照LC和LLC进行。值得注意的是，当β<0.051时，LLC最大的偏差不超过0.17%，因此可以近似认为是无偏的。

HyperLogLog Counting

主要改进

调和平均数以减少离群值影响

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

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

$$\hat{n}=\frac{{\alpha }_m{m}^2}{\Sigma 2^{-M}},{\alpha }_m=(m{\int }_0^{\infty }({\log }_2(\frac{2+u}{1+u}){)}^{m}du{)}^{-1}$$

HLLC通过这样改进将误差从$\frac{1.3}{\sqrt{m}}$ 优化到了 $\frac{1.04}{\sqrt{m}}$，而空间复杂度与LLC相同。

分段偏差修正

给出了在n较小或较大时的偏差修正方案，具体来说
$$估算方法=\{\begin{array}{l}\begin{array}{rl}LinearCounting& ifE\le \frac{5m}{2}\\ HyperLogLogCounting& if\frac{5m}{2}<E\le \frac{2^{32}}{30}\\ -2^{32}\log (1-\frac{E}{2^{32}})& ifE>\frac{2^{32}}{30}\end{array}\end{array}$$

参考资料

《A Linear-Time Probabilistic Counting Algorithm for Database Applications》
《LogLog Counting of Large Cardinalities (Extended Abstract)》
《HyperLogLog: the analysis of a near-optimal cardinality estimation algorithm》
《Fast and Accurate Traffic Matrix Measurement Using Adaptive Cardinality Counting》
博客：http://blog.codinglabs.org/articles/algorithms-for-cardinality-estimation-part-i.html