应用场景
基数统计(Cardinality Counting)指计算一个数据集中不同元素的数量，在很多场景都需要这样的功能：电商场景中的独立用户(UV)数量统计、数据库中快速计算字段取值数量以优化query、计算与某个站点相关的不同链接数量等。
为实现这样的功能，最直观的方法就是使用 Set 记录每一个元素，在数量巨大的情况下，会造成大量的内存占用。以存储 $2^{30}$ (约10亿)个长度为64bit的url为例，需要占用 $64 \times 2^{30} = 64GB$ 的空间。Set方案的空间占用与元素的数量成正比，空间复杂度为 $O(N)$
如果使用bitmap来存储，则需要将不同的url映射到不同的数字。即使最理想的情况下，url之间没有冲突，且从1到 $2^{30}$ 都没有空闲，bitmap也需要占用 $2^{30}\div 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})$ 。
主要步骤
- 数据哈希：假设原始数据的基数为n，使用一组哈希空间为m的哈希函数H，将原始数据均匀分组。
- 分桶数据统计：构建一个长度为m的bitmap，其中每一个bit对应哈希空间的一个值。生成哈希数组的值如果存在，则把相应的bit设置为1。
- 数量估计：当所有值设置完成后，统计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}} = me^{-\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*(-\ln{\frac{u}{m}})=mf(V_n)$
在 $V_n$ 处进行泰勒展式:
$$ \hat{n} = mf(V_n) = m(f(p)+(V_n-p)f’(p)+\frac{1}{2}(V_n-p)^2f’'(p)+\frac{1}{6}(V_n-p)^3f’''(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的选择
误差控制
如果我们需要将误差控制在 $\epsilon$ 以下，即 $$\frac{\sqrt{m(e^t-t-1)}}{n} = \frac{\sqrt{\frac{e^t-t-1}{t}}}{t} \lt \epsilon$$ 得到m的最小数量为 $$m \gt \frac{e^t-t-1}{(\epsilon t)^2}$$
满桶控制
可以通过让 $U_n$ 的期望相比0的偏移大于若干倍的 $U_n$ 标准差来控制： $$E(U_n)-a\times StdDev(U_n) >0$$ 代入 $E(U_n)=me^{-t}$ 和 $StdDev(U_n)=(me^{-t}(1-(1+t)e^{-t}))^{\frac{1}{2}}$，得到
$$m>a^2(e^t-t-1)$$
结合两个限制条件
结合误差控制和满桶控制，最终m在期望误差内的最小长度为: $$m>\beta (e^t-t-1), \quad \beta=max(a^2, \frac{1}{(\epsilon t)^2})$$
确定a
u满足二项分布，当n非常大而落进一个桶的概率 $p=\frac{1}{m}$ 非常小时，可以用泊松分布进行近似 $$\lim_{n,m \to \infty}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})$。结果如下表
并集与交集
只要哈希函数和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_max$ ，则 $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_m2^{\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 \gt 64$ 时 $StdError(\frac{\hat{n}}{n}) \approx \frac{1.3}{\sqrt{m}}$
如果要将误差控制在 $\epsilon$ 之内，则
$$m \gt (\frac{1.30}{\epsilon})^2$$
空间复杂度
内存使用与m的大小及哈希值长度有关，如果统计的基数大小为n，则$\rho_{max} \ge O(\log_2{\frac{n}{m}})$，对应保存 $M$ 的空间为 $O(m \times \log_2{\rho_{max}}) $，在m固定时，空间复杂度为 $O(log_2(log_2n))$ 假设哈希后长度为64bit，则 $\rho_{max} \le 64$，因此每个桶需要 $\log_2{64} = 6$ bit的空间来存储 $\rho_{max}$，分成m个桶即 $\frac{6\times m}{8}$ 字节。当分桶数m为1024时，误差为 $1.3 \div \sqrt{1024} \approx 4%$，需要的字节数就是 $6\times 1024 \div 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的标准误差为 $$StdError_{linear_counint}(\frac{\hat{n}}{n}) = \frac{\sqrt{e^t-t-1}}{t\sqrt{m}}$$ $$StdError_{loglog counting}(\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{cases} \alpha_m{m2^{\frac{\Sigma{M}}{m}}} , \text{if} , 0 \le \beta \le 0.051 \\ -m\ln(\beta) , \text{if} , 0.051 \le \beta \le 1 \end{cases}$$
因为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}}}, \quad \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{cases} \begin{aligned} Linear Counting, &if , E \le \frac{5m}{2} \\ HyperLogLog Counting , &if , \frac{5m}{2} \lt E \le \frac{2^{32}}{30} \\ -2^{32}\log(1-\frac{E}{2^{32}}) , &if , E \gt \frac{2^{32}}{30} \end{aligned}\end{cases}$$
参考资料
《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
loglog demo: http://content.research.neustar.biz/blog/hll.html
