从 fair coin 说起

这是一个经典老题了，如何用 fair coins 来构造一个 biased coin，以及如何用 fair coins 来构造一个 discrete uniform distribution。基于这两个问题，我又 yy 了另外几个小 follow-up，试图去分析其时间复杂度。

Case 1: \(\text{Bernoulli}(p)\)

先来一个 warm-up exercise：在只有 fair coins （即 \(\text{Bernoulli}(0.5)\)）的情况下，给定任意一个 \(0 < p < 1\)，如何从 \(\text{Bernoulli}(p)\) 中 sample 一个元素？

这个很好做：我们只要 sample \([0, 1)\) 之间的一个实数 \(x\)，并和 \(p\) 比较就好了。我们考虑逐位 sample 并比较：每次得到 \(x\) 的二进制表示第 \(i\) 位，\(x_i \sim \text{Bernoulli}(p)\)，如果 \(x_i \neq p_i\)，则我们就可以得到 \(x\) 和 \(p\) 的相对大小关系了，否则继续下一位。

def sample(p):
    while True:
        x = random.randint(2)
        p_i = int(p * 2)
        if x != p_i:
            return int(x < p_i)
        p = p * 2 - p_i

follow-up： 时间复杂度

上面这个算法有最坏时间复杂度保证吗？

显然没有。

那么存在一个最坏时间复杂度保证的算法吗？

对于几乎所有的 \(p\) 都不存在。假设存在一个算法能保证 \(n\) 步内结束，那么我们考虑用到的所有 \(n\) 个 random bit，一共有 \(2^n\) 种组合。对于每种组合我们一定会返回一个 0 或者 1，所以算法一定是从 \(\text{Bernoulli}(\frac{q}{2^n})\) 中 sample，其中 \(q\) 为整数。而几乎所有的 \(p\) 都不能表示成 \(\frac{q}{2^n}\) 这样的形式。

那上面这个算法的期望复杂度是多少？

显然，每一步有 \(\frac{1}{2}\) 的概率直接结束，所以期望就是 2 步结束。挺意外的一个结果，和 \(p\) 无关。

Case 2: \(U(n)\)

令 \(n \in \mathbb{Z}^+\)，如何用 fair coins 来得到一个 \(U(n)\) 的 sample 呢？这里 \(U(n)\) 指 \(\{0, 1, \dots, n - 1\}\) 上的 uniform distribution。

我们每次 sample \(k = \lceil \log_2 n \rceil\) 个 bit，得到一个 \(U(2^k)\) 的 sample。如果这个数小于 \(n\)，那么返回即可，否则重新开始。

def sample(n):
    k = int(math.ceil(math.log2(n)))
    while True:
        x = 0
        for i in range(k):
            x = x * 2 + random.randint(2)
        if x < n:
            return x

follow-up：时间复杂度

上面这个算法有最坏时间复杂度保证吗？

显然没有。

那存在一个有时间复杂度保证的算法吗？

和之前一模一样的 argument，没有。

那这个算法的期望时间复杂度是多少？

每一轮需要 \(k\) 个 bit，并且有 \(\frac{n}{2^k}\) 的概率结束，所以期望时间复杂度就是 \(\frac{k 2^k}{n}\)。

follow-up：更优的算法

存在期望复杂度更低的算法吗？

存在。和 warm-up exercise 类似，还是 sample 一个 \([0, 1)\) 之间的实数 \(x\)，答案即为 \(\lfloor xn \rfloor\)。我们从高位到低位枚举 \(x\) 二进制表示，然后用一个区间 \([l, r)\) 表示 \(x\) 可能取值，一旦发现 \(\lfloor ln \rfloor = \lfloor r n \rfloor\)，那么我们就已经得到了 \(\lfloor xn \rfloor\)，返回即可。

def sample(n):
    l, r = 0, 1
    while True:
        if random.randint(2):
            l = (l + r) / 2
        else:
            r = (l + r) / 2
        if int(l * n) == int(r * n):
            return int(l * n)

follow-up：避免浮点数操作

上面这个做法需要用到两个浮点数 \(l, r\)，其精度有限。能避免使用无限精度的数据类型吗？

可以。我们回到一开始的算法，它没有达到最优是因为，如果 \(n \leq x < 2^k\)，我们只会重来，白白浪费了这里面用到的 randomness。新算法如下：我们维护两个数 \(x, m\)，表示现在我们有一个 \(U(m)\) 的 sample \(x\)。如果 \(m \geq n\)，那么我们观察 \(x\)：

• 如果 \(x < n\)，那么 \(x\) 就是 \(U(n)\) 的一个 sample，直接返回即可；
• 否则我们就知道 \(x\) 是 \(U(n, m)\) 的一个 sample，故 \(x' = x - n\) 就是 \(U(m - n)\) 的一个 sample。 而每次拿到一个新 bit 后，我们可以得到一个 \(U(2m)\) 下的 sample \(x' = 2x + \text{Bernoulli}(0.5)\)，重复操作直到 \(m \geq n\) 为止。

def sample(n):
    x, m = 0, 1
    while True:
        x, m = x * 2 + random.randint(2), m * 2
        if m >= n:
            if x < n:
                return x
            x, m = x - n, m - n

follow-up：时间复杂度分析

上面这个算法期望时间复杂度，\(f(n)\)，具体为多少？

显然对 \((x, m)\) 整一个高斯消元是可行的，但是这有点复杂了。聪明一点的同学可以考虑 \(f_m\) 表示在已经有了一个 \(U(m)\) 的 sample 后得到 \(U(n)\) 的 sample 的期望复杂度，我们有： \[ f_m = \begin{cases} \frac{m-n}{m} f_{m-n}, & m \geq n \\ f_{2m} + 1 , & m < n\end{cases} \] 这样未知数会少一些……

但是 somehow 我们有更聪明的方法。令 \(Y\) 表示几轮之后算法结束，我们先来一个常用技巧： \[ \mathbb{E}[Y] = \sum_{y \geq 0}^\infty \Pr[Y \geq y]. \] 于是问题转化为计算 \(y\) 步之后算法还没有结束的概率是多少。我们把算法作一点点修改：之前是如果 \(x < n\) 那么我们返回 \(x\)，现在是如果 \(x \geq m - n\) 我们返回 \(x \bmod n\)，在 \(x \sim U(m)\) 的情况下这两者是等价的。

def sample(n):
    x, m = 0, 1
    while True:
        x, m = x * 2 + random.randint(2), m * 2
        if m >= n:
            if x >= m - n:
                return x % n
            m = m - n

我们可以观察到，第 \(y\) 轮结束之后的 \(m\) 刚好就是 \(2^y \bmod n\)，于是假设有无限精度，再改动一点点点点：

def sample(n):
    x, m = 0, 1
    while True:
        x, m = x * 2 + random.randint(2), m * 2
        if x >= m % n:
            return x % n

然后就可以很清晰的看到，\(y\) 轮之后仍然不结束的概率是 \[ \Pr[Y \geq y] = \Pr[x < (2^y \bmod n)] = \frac{2^y \bmod n}{2^y}. \] 故算法的期望时间复杂度为 \[ f(n) = \mathbb{E}[Y] = \sum_{y=0}^\infty \frac{2^y \bmod n}{2^y}. \]

虽然这个求和有无数项，但是只要注意到 \(2^y \bmod n\) 有循环节就好办了。

follow-up：时间复杂度上下界

刚才我们得到了 \(f(n)\) 的精确表达式，但是看起来不存在一个对于任何一个 \(n\) 均成立的 closed-form。那我们能否得到一个 \(f(n)\) 的上下界？

一个显然的，不需要任何推导的下界是： \[ f(n) \geq \log_2 n, \] 因为 \(U(n)\) 包含 \(\log_2 n\) bit 的信息，所以我们期望至少 \(\log_2 n\) 次才能得到一个 sample……

稍微一分析，我们就可以得到一个 \(f(n)\) 的更强的下界。令 \(k = \lceil \log_2 n\rceil\)，而 \([0, k-1]\) 内的 \(y\) 均有 \(\frac{2^y \bmod n}{2^y} = 1\)，故有： \[ f(n) \geq k \geq \log_2 n. \]

那上界呢？我能想到的最好的一个是：我们把 \(y\) 分成三组来估计 \(\frac{2^y \bmod n}{2^y}\),

• \(0 \leq y < k\)：很显然，此时 \(\frac{2^y \bmod n}{2^y} = 1\)，这一组内所有的和为 \(k\)；
• \(y = k\)：有 \(\frac{2^k \bmod n}{2^k} \leq \frac{2^k - n}{2^k} = 1 - \frac{n}{2^k}\)；
• \(y > k\)：有 \(\frac{2^y \bmod n}{2^y} < \frac{n}{2^y}\)，这一组内所有的和有上界 \(\sum_{y > k}^\infty \frac{n}{2^y} = \frac{n}{2^k}\)；

故 \(f(n)\) 有上界 \[ f(n) = \sum_{y=0}^\infty \frac{2^y \bmod n}{2^y} < k + 1 - \frac{n}{2^k} + \frac{n}{2^k} = k + 1. \] 综上，我对 \(f(n)\) 的最好的 bound 就是 \[ \log_2 n \leq \lceil \log_2 n \rceil \leq f(n) < \lceil \log_2 n \rceil + 1, \] 前面的 \(\log_2 n\) 就是信息论的 lower bound。可以看到，这个算法最多比理论下界多 2 bits，还算是不错的吧……