为什么我的眼里常含泪水？因为我有一个算法不会。为了节约点眼泪，今天我们就来介绍著名的Flajolet-Martin算法，尽管这是一种始于上世纪80年代中期的算法，但是今天这个时代再来探讨它似乎更有意义，因为它是一种真正的“大数据”算法。

大数据时代，我们所面对的将是海量数据。所以进行任何一个看似微不足道的操作之前，都应该慎重地考察一下它的执行代价。假设你有一个很大的数据集，非常非常大，以至于不能全部存入内存。并且这个数据集中可能包含有重复的数据。如果想知道这个数据集中有多少个不同的（也就是不重复的）元素，最straightforward的方法似乎就是去‘直接数一数’，但这样做其实并不实际。

我们都知道对无序数据进行排序，时间复杂度为O(n*log(n))，然后扫描一遍这样的数据就能知道其中包含有多少个不同的元素。扫描时间复杂度为O(n)，所以总的时间复杂度为O(n*log(n))。但是别忘了，排序就意味着我们必须得找个同样大小的空间来储存排序后的结果。于是空间复杂度是O(n)的。由于数据量太大，可能内存根本装不下，所以排序是不切实际的。此外，我们还可以考虑使用哈希，将时间复杂度降到O(n)。但空间窘境仍然没有得意改善。所以，在大数据情况下，不论是基于排序还是基于哈希的方法都不具备可行性。

对于海量数据而言，其实最终我们并不需要得到一个十分精确的结果，或者说一个近似的结果也是可以接受的。那么我们能不能设计一个近似的方法，并保证算法在时间和空间上都具有较高效率呢？这就是本文要介绍的Flajolet-Martin算法。该算法最初出现在Flajolet和Martin于1985年发表的论文”Probabilistic Counting Algorithms for Data Base Applications”中。后来，Flajolet等人又陆续发表了该算法的改进版本（主要是LogLog counting of large cardinalities 和  HyperLogLog：The analysis of a near-optimal cardinality estimation algorithm这两篇文章）。

我们先来从一个intuitive的情况入手。我们怎么估计结果数据集中有多少非重复的数字呢？了解到原来的数据集是随机数，且充分分散，一个非常简单的方法是：找出最小的数字。如果最大的可能的数值是 m，最小的值是 x，我们 可以估计大概有 m/x 个非重复的数字在数据集里面。举个例子，如果我们扫描一个数字在 0 到 1 之间的数据集，发现最小的数字是 0.01。我们有理由猜想可能数据集里大概有 100 个非重复的数字。如果我们找到一个更小的最小值的话，可能包含的数据个数可能就更多了。请注意不管每个数字重复了多少次都没关系，这也是很显然的。

这个过程的优点是非常直观，但同时它也很不精确。不难举出一个反例：一个只包含少数几个非重复数字的数据集里面有一个很小的数。同样的一个含有许多非重复数字的数据集含有一个比我们想像中更大的最小值，用这种估计方法也会很不精确。最后，很少有数据充分分散充分随机的数据集（所以一个改进的灵感就是把数据转换为一个更加均匀的分散的分布）。但是这个算法原型给了我们一些灵感使得我们有可能达到我们的目的，我们需要更精致一些的算法.

Flajolet和Martin最早发现通过一个良好的哈希函数，可以将任意数据集映射成服从均匀分布的（伪）随机值序列。根据这一事实，可以将任意数据集变换为均匀分布的随机数集合，然后就可以使用上面的方法进行估计了，不过只是这样是远远不够的。

接下来，他们陆续发现一些其它的基数（非重复数个数）估计方法，而其中一些方法的效果优于之前提到的方法（也就是记录最小的哈希值的方法）。Flajolet和Martin计算了哈希值的二进制表示的0后缀（即哈希后的值的末尾的 0 字符的个数），结果发现在随机数集合中，通过计算每一个元素的二进制表示的0后缀，设k为最长的0后缀的长度，则平均来说集合中大约有2^k个不同的元素；（显然在一个随机的数据集中，平均每 2^k 个元素就出现一个长度为 k 的全为 0 的比特序列。我们要做的就是找出这些序列并记录最长的来估计非重复元素的个数。）

假定哈希函数H(e)能够把元素e映射到[0, 2^m-1]区间上；再假定函数TailZero(x)能够计算正整数x的二进制表示中末尾连续的0的个数，譬如TailZero(2) = TailZero(0010) = 1，TailZero(8) = TailZero(1000) = 3，TailZero(10) = TailZero(1010) = 1；我们对每个元素e计算TailZero(H(e))，并求出最大的TailZero(H(e))记为Max，那么对于独立元素数目的估计为2^Max。

举例来说，给定序列{e1, e2, e3, e2}，独立元素数目N = 3。假设给定哈希函数H(e)，有

• H(e1) = 2 = 0010，TailZero(H(e1)) = 1
• H(e2) = 8 = 1000，TailZero(H(e2)) = 3
• H(e3) = 10 = 1010，TailZero(H(e3)) = 1

第1步，将Max初始化为0；
第2步，对于序列中第1项e1，计算TailZero(H(e1)) = 1 > Max，更新Max = 1；
第3步，对于序列中第2项e2，计算TailZero(H(e2)) = 3 > Max，更新Max = 3；
第4步，对于序列中第3项e3，计算TailZero(H(e3)) = 1 ≤ Max，不更新Max；
第5步，对于序列中第4项e2，计算TailZero(H(e2)) = 3 ≤ Max，不更新Max；
第6步，估计独立元素数目为N’ = 2^Max = 2^3 = 8。

在这个简单例子中，实际值N = 3，估计值N’ = 8，误差比较大。此外，估计值只能取2的乘方，精度不够高。

我们可以用这个方法估计基数。但是，这仍然不是很理想的估计方法，因为和基于最小值的估计一样，这个方法的方差很大。不过另一方面，这个估计方法比较节省资源：对于32位的哈希值来说，只需要5比特去存储0后缀的长度。

在实际应用中，为了减小误差，提高精度，我们通常采用一系列的哈希函数H1(e), H2(e), H3(e)……，计算一系列的Max值Max1, Max2, Max3……，从而估算一系列的估计值2^Max1, 2^Max2, 2^Max3……，最后进行综合得到最终的估计值。具体做法是：首先设计A*B个互不相同的哈希函数，分成A组，每组B个哈希函数；然后利用每组中的B个哈希函数计算出B个估计值;接着求出B个估计值的算术平均数为该组的估计值；最后选取各组的估计值的中位数作为最终的估计值。

举例来说，对于序列S，使用3*4 = 12个互不相同的哈希函数H(e)，分成3组，每组4个哈希函数，使用12个H(e)估算出12个估计值：

第1组的4个估计值为<2, 2, 4, 4>，算术平均值为(2 + 2 + 4 + 4) / 4 = 3；
第2组的4个估计值为<8, 2, 2, 2>，算术平均值为(8 + 2 + 2 + 2) / 4 = 3.5；
第3组的4个估计值为<2, 8, 8, 2>，算术平均值为(2 + 8 + 8 + 2) / 4 = 5；

3个组的估计值分别为<3, 3.5, 5>，中位数为3.5；

因此3.5 ≈ 4即为最终的估计值。

分析到目前为止给出的FM算法的时间复杂度。假定序列长度为N，哈希函数H(e)的数目为K。初始化K个Max值的时间复杂度为O(K)；对N个元素e使用K个哈希函数H(e)计算TailZero(H(e))并更新Max值的时间复杂度为O(N*K)；综合K个Max值给出最终估计值的时间复杂度为O(K)。因此总的时间复杂度为O(N*K)。

分析FM算法的空间复杂度。该算法需要存储K个Max值，而每个元素e在进行相关计算后就可以丢掉。因此总的空间复杂度为O(K)。

综上所述，FM算法的时间复杂度为O(N*K)，空间复杂度为O(K)。一般来说K比较小，可以认为FM算法的时间复杂度为O(N)，空间复杂度为O(1)。

从实验统计上来看这给了我们一个相当好的结果，但哈希的代价的是很高的。一个更好的方式是一个叫做随机平均（stochastic averaging）的方法。这种方法不是使用多个哈希函数，而是使用一个哈希函数，但是将哈希值的区间按位切分成多个桶（bucket）。也就是把哈希的输出进行分割，然后使用它的一部分作为桶序号来放到许多桶中一个桶里去。

假设我们需要 1024 个值，我们可以使用哈希函数的前 10 个比特值作为桶的序号，然后使用剩下的哈希值来计算前导 0 比特序列。这个方法并不会损失精确度，但是节省了大量的哈希计算。把我们目前学到的应用一下，这里有一个简单的基于Python的实现。这和 Durand-Flajolet 的论文中的算法是等价的，为了实现方便和清晰所以我计算的是尾部的 0 比特序列。结果是完全等价的。

def trailing_zeroes(num):
        """Counts the number of trailing 0 bits in num."""
        if num == 0:
            return 32 # Assumes 32 bit integer inputs!
        p = 0
        while (num >> p) & 1 == 0:
            p += 1
        return p
     
def estimate_cardinality(values, k):
        """Estimates the number of unique elements in the input set values.
     
        Arguments:
            values: An iterator of hashable elements to estimate the cardinality of.
            k: The number of bits of hash to use as a bucket number; there will be 2**k buckets.
        """
        num_buckets = 2 ** k
        max_zeroes = [0] * num_buckets
        for value in values:
            h = hash(value)
            bucket = h & (num_buckets - 1) # Mask out the k least significant bits as bucket ID
            bucket_hash = h >> k
            max_zeroes[bucket] = max(max_zeroes[bucket], trailing_zeroes(bucket_hash))
        return 2 ** (float(sum(max_zeroes)) / num_buckets) * num_buckets * 0.79402

这段代码实现了上面讨论的估计算法：保持一个计算前导（或尾部）0个数的数组，然后在最后对个数求平均值，如果我们的平均值是 x，我们的估计就是 2^x 乘以桶的个数。前面没有说到的是这个魔法数 0.79402。数据统计表明我们的程序存在一个可预测的偏差，它会给出一个比实际更大的估计值。这个在 Durand-Flajolet 的论文中导出的魔法常数是用来修正这个偏差的。实际上这个数字随着使用的桶的个数（最大2^64）而发生变化，但是对于更多数目的桶数，它会收敛到我们上面用到的算法的估计数字。有兴趣深入研究的读者可以参看完整的原始论文，以了解跟多信息，包括那个魔术数是怎么导出的。

这个程序给了我们一个非常好的估计，对于 m 个桶来说，平均错误率大概在 1.3/sqrt(m) 左右。因此对于分桶数为1024的情况（所需内存1024 * 5 = 5120，即 640 字节），我们大约会有4%的平均误差。每桶5比特的存储已经足以估计227的数据集，而我们只用的不到1k的内存！

来看一下试验结果：

>>> [100000/estimate_cardinality([random.random() for i in range(100000)], 10) for j in range(10)]
    [0.9825616152548807, 0.9905752876839672, 0.979241749110407, 1.050662616357679, 
     0.937090578752079, 0.9878968276629505, 0.9812323203117748, 1.0456960262467019, 
     0.9415413413873975, 0.9608567203911741]

不错！虽然有些估计误差大于4%的平均误差，但总体来说效果良好。如果你准备自己做一下这个试验，有一点需要注意：Python内置的 hash() 方法将整数哈希为它自己。因此诸如 estimate_cardinality(range(10000), 10) 这种方式得到的结果不会很理想，因为内置 hash() 对于这种情况并不能生成很好的散列。但是像上面例子中使用随机数会好很多。

虽然我们已经得到了一个非常好的估计，但它有可能做到更好。Durand 和 Flajolet 发现极端数值会很大地影响估计结果的准确度。通过在求平均前舍弃一些最大值，准确度可以得到提高。特别地，舍弃前 30% 大的桶，仅仅计算 70% 的桶的平均值，精确度可以用 1.30/sqrt(m) 提高到 1.05/sqrt(m)! 这意味着在我们之前的例子中，用 640 字节的状态，平均错误率从 4% 变成了大约 3.2%。但并没增加空间的使用。

最后，Flajolet et al 的论文的贡献就是使用了一个不同类型的平均数。使用调和平均数而不是几何平均数。通过这么做，我们可以把错误率降到 1.04/sqrt(m)，同样不增加需要的空间。当然完整的算法要更复杂一点，因为它必须修正小的和大的基数误差。有兴趣的读者还是请参考一下原文。

参考文献：

http://blog.notdot.net/2012/09/Dam-Cool-Algorithms-Cardinality-Estimation