Alias Method——高效的离散分布采样算法

文章目录

一、采样方法
二、轮盘赌算法
三、Alias Method
四、伪代码
参考文献

一、采样方法

采样的本质是随机现象的模型，根据给定的概率分布,来模拟产生一些随机时间。另一方面,采样得到的样本集也可以看作是是一种非参数模型。即用较少量的样本点(经验分布)来近似总体分布,并刻画总体分布中的不确定性。从这个角度来说,采样其实也是一种信息降维,可以起到简化问题的作用。

对当前的数据集进行重采样,可以充分利用已有的数据集,挖掘更多信息，如Bootstrap法和jackknife法,通过对严格不能多次重采样来估计统计量的偏差、方差等。另外,利用重采样,可以在保持特定的信息下(目标信息不丢失),有意识地改变样本的分布,以更适应后续的模型训练和学习,例如利用重采样类处理分类模型的训练样本不均衡。

此外,许多模型由于结构复杂、含有隐变量等原因,导致对应的求解公式比较复杂,没有显式解析解,难以进行精确求解或推理。在这种情况下,可以利用采样方法进行随机模拟,从而对这些复杂模型进行近似求解或推理。一般会转化为某些在特定分布下的积分或期望,或者是求某些随机变量或者参数在给定数据下的后验分布等。

在一些基于随机游走的算法中，随机采样方法无论从功能上还是效率上都扮演了不可或缺的角色。首先对于node2vec图表示学习方法自身就是基于随机游走的算法，随机游走的每一步都是一个对于多项（Multinoulli）分布的采样，采样的次数是O(Nml)O(\text{Nml})O(Nml)，其中N是图中的节点数，l是随机游走路径的长度，m是随机游走算法迭代的次数。对于本课题中使用的数据接地节点数约10000，游走序列的长度和迭代次数分别为80和10，对于这样的参数设置，共需要运行8000000次离散分布采样算法，因此该算法的效率至关重要。另外在node2vec训练语料生成的过程中需要降采样算法提升模型效果，在训练过程中需要负采样算法优化训练速度。

对于离散分布的采样问题，并不是很难的问题，但是在本课题中由于大量的执行了离散采样操作，因此要考虑其运行效率。实验证明，使用Alias Method对于基于图嵌入的用户表示学习的效率提升显著。

采样算法分为连续分布采样算法和离散分布采样算法。连续分布采样算法中最常见的均匀分布采样算法是所有的采样算法——包括连续分布采样算法和离散分布采样算法——的基础。而对于均匀分布的采样也非常简单，通常建立在一个良好的硬件的随机数种子发生器之上，然后并使用该种子初始化一块指定位数的内存即可，由于浮点数机器的精度问题，不可能真正的模拟连续分布，但是通常在数值计算的实践中，该问题可以忽略不计。另外最常见的分布是高斯分布和截断离散分布，高斯分布可以使用较为通用的拒绝采样法或者逆变换采样，出于效率问题，在真正的实现中通常使用Box-Muller方法生成符合标准高斯分布的随机变量，然后再根据参数进行和平移变换。对于截断高斯分布其实就是对于同参数的高斯分布进行拒绝采样，简单的抛弃在截断边界之外的随机值并重新生成。除此之外的其他分布通常可以基于均匀分布和高斯分布的随机变量构建而成。另外，不规则的分布可以使用拒绝采样及其高效的变种，重要性采样通常用于计算期望或者逼近积分。在机器学习中，尤其是生成式模型中，以马尔科夫链蒙特卡洛采样算法为代表的一系列蒙特卡洛方法可以使用重复采样的方法逼近一些复杂或难以计算的参数并为其提供一个合理的近似值。

二、轮盘赌算法

而离散分布的采样算法较为简单，大致的接口如上图所示，通常可以抽象地认为Multinoulli分布是一个具有多面的有偏骰子模型。通常离散分布的采样算法需要两个阶段来完成，第一个阶段是根据指定分布的参数——通常是分布列——来进行准备步骤或初始化步骤，比如建立查询表，同一个分布只需要一次准备步骤。第二个阶段是正式的采样阶段，该步骤通常需要进行一次或者多次符合离散均匀分布的随机变量——通常是正整数或自然数。传统的对于任意离散分布的采样算法使用查表法或者轮盘赌算法（Roulette Wheel Selection Algorithm）。前者单次的生成时间复杂度为O(1)O(1)O(1)，因为只需要一次均匀离散分布的生成与一次查表。但是准备时间较长，因为需要建立查询表。在最好情况下，准备步骤的时间复杂度是O(n)O(n)O(n)。但是在最坏情况下，准备步骤的时间复杂度是O(∏i=0ndi)O(\prod_{i= 0}^{n}d_{i})O(∏i=0ndi),其中did_idi是离散分布中各个取值的分数表示的概率的分母。需要说明的是，这种最坏情况是指d两两互质的情况。虽然这种情况在实际情况下难以发生，但是通常来说，准备步骤的时间复杂度仍然远大于O(n)O(n)O(n)，比如分布列为 [11511056][\frac{1}{15}\frac{1}{10}\frac{5}{6}][15110165]的简单离散分布，其查询表的大小为30。而另外一个方法，轮盘赌法则在整体上表现较好，而且在各个具体应用中使用的也更加广泛，轮盘赌算法的本质思想是使用多个Bernoulli分布作为一个单独的Multinoulli分布的代替。

上图是轮盘赌法思想的简单示意图，分布遍历各个离散变量的各个输出值仅进行一次伯努利分布的采样，表示当前值是否会成为最终的输出值。如果将整个离散分布的采样看做是一个掷一个有偏多面骰子的过程，那么在轮盘赌法内部实际上是在用多次的掷硬币的操作来模拟该过程，而通过名字也可以大致反映出该思想——转一个转轮，而小钢球依次滚过每一个位置并随时可能停下。唯一不同的是实际的每次内部的伯努利分布的参数并不是分布列中对应取值的概率，而是会每次调整以保证最多一次遍历一定得到结果。该方法并不是真正的轮盘赌法，可以称之为“链式法”。在实际的实现中，轮盘赌算法不会真的遍历整个分布列，而是累积分布函数的单调递增性质，对累积概率表进行二分查找。因此，轮盘赌算法的初始化时间复杂度是O(n)O(n)O(n)，空间复杂度是O(n)O(n)O(n)，而每次的采样操作的时间复杂度仅需要O(log⁡n)O(\log n)O(logn)。另外，由于平衡树只能保证子树的高度平衡而不能保证各个节点被访问的概率平衡，因此可以根据节点被访问的概率构造最优平衡树，这使得我们在生成步骤可以有更小的底数的O(log⁡n)O(\log n)O(logn)复杂度的表现，但是我们付出的代价是 O(n2)O(n^{2})O(n2) 的初始化的复杂度。

三、Alias Method

这二者的效率均不能满足在基于游走的图表示学习方法中的优化需求,实际的随机游走算法中使用的是Alias Method。Alias Method结合了上述两种方法的各自优势。最终的效果是，使用一次离散均匀分布采样级联一次伯努利分布采样来等效进行一次多项（Multinoulli）分布的采样。基本思路是，首先进行一次均匀多项分布采样决定本次测试取值，由于各个可能的取值等概率，因此这一操作的时间复杂度是常数级。然后进行一个有偏的伯努利分布采样来决定被测试的取值是否生成。

该步操作的时间复杂度同样是常数级。如上图所示，图（a）是Alias Method的基本思想，根据图示在此实际举例说明生成过程。第一步在 O(1)O\left(1\right)O(1)的时间复杂度内“掷出骰子”随机选择了B，然后“翻转不规则硬币”得到反面，则意味着本次生成失败，沿图中“If Get Negative Bernoulli”的路径返回，重复上述过程。假设本次“掷骰子”选择了C并且“扔硬币”步骤得到正面，那么生成成功，则另随机变量取值C，结束算法。注意，直观的观察图（a）中的概率表则发现大部分区域为代表Bernoulli分布Nagetive的区域，这会导致采样失败的概率很高。注意到每次Bernoulli采样中的参数是相对值，不一定要等于原始多项分布的概率，因此可以对其进行缩放，如图（b）所示。最佳的缩放因子是多项分布中的最大的概率值，本例以0.4作为缩放因子。改进后成功概率大幅提高，对于可能取值越多的多项分布越明显，然而，依然存在持续失败的可能行。我们可以计算其成功的概率。

∑i=0n−1(1n∗pipmax⁡)=1n∑i=0n−1pipmax⁡=1n∗pmax⁡∑i=0n−1pi=1n∗pmax⁡\sum_{i = 0}^{n - 1}\left( \frac{1}{n}*\frac{p_{i}}{p_{\max}}\ \right)\ = \ \frac{1}{n}\sum_{i = 0}^{n - 1}\frac{p_{i}}{p_{\max}}\ = \ \frac{1}{n*p_{\max}}\sum_{i = 0}^{n - 1}p_{i}\ = \ \frac{1}{n*p_{\max}} i=0∑n−1(n1∗pmaxpi ) = n1i=0∑n−1pmaxpi = n∗pmax1i=0∑n−1pi = n∗pmax1

因此，期望的内部运行次数是 n∗pmax⁡n*p_{\max}n∗pmax。这意味着最大的概率越大，越容易失败，而当 pmax⁡p_{\max}pmax 等于 1n\frac{1}{n}n1 的时候，运行次数为1，此时该方法实际特殊化为普通的离散均匀分布采样算法。因此可知，该方法期望运行时间复杂度为O(n)。

实际上，考虑失败情况的根本原因，我们将查询表看成一个矩形，将先后两次的简单采样分别看作是在横轴和纵轴上的随机取值，此时的概率模型实际上为几何概型，而失败的概率就是所有失败区域的面积占矩形总面积的比例，如图中的图（a）所示。前述采用的重缩放的方法本质上是在随机数1（Bernoulli分布采样）所在的方向上尽可能去掉所有的可能取值的区域内皆失败的部分。然而我们考虑失败区域产生的原因，本质上是由于矩形总面积大于成功区域总面积导致的，那么接下来的改进就顺理成章了，只需要继续降低矩形的高度直到矩形总面积等于与成功区域的面积相等即可，然后对于高度大于矩形的部分重新分配，使得成功区域的面积恰好填充矩形面积即可。可以简单的计算得到，调整后的高度为 p‾i=1n{\overset{\overline{}}{p}}_{i} = \frac{1}{n}pi=n1，即各个可能取值的概率的平均数。或者也可以保持矩形高度为1不动，将所有可能的取值的对应概率乘以缩放因子n。

如上图所示，图（b）表示调整后的新的矩形，由于高度小于pmax⁡\text{~~}p_{\max} pmax,因此一定会有一部分成功区域的面积在矩形之外，此时我们需要对面积进行切分与调整。因为精心调整后的矩形面积等于成功区域的面积，因此一定可以存在一种方式调整成功。在本例子中，将取值C对应的多余面积移动至A的缺失部分，将取值B的多余面积移动到D的缺失部分，得到结果在图（c）中展示。对于每个列，仍然是进行一次伯努利采样，此时可知采样一定成功。移动后的面积所代表的采样结果被称为“Alias”，而原分布的采样概率通过 pip‾\frac{p_{i}}{\overset{\overline{}}{p}}ppi来计算并存储在一个表中。与之前的算法不同的是，我们需要一个额外的表来存档我们移动的路径，即每个位置上多出来的部分是来自于原本的哪个取值。两个表的具体状态如图（d）所示。在具体的采样步骤，同样使用一个离散均匀分布生成1到N之间的一个正整数然后根据概率表中的概率为参数进行一次伯努利采样，如果成功，返回对应的取值，如果失败，返回Alias表中的值。

算法	初始化时间复杂度	采样时间复杂度	空间复杂度
查表法	O(∏i=0ndi)O(\ \prod_{i = 0}^{n}d_{i}\ )O( ∏i=0ndi )	O(1)	O(∏i=0ndi)O(\ \prod_{i = 0}^{n}d_{i}\ )O( ∏i=0ndi )
链式法	O(n)	O(n)	O(n)
轮盘赌法	O(n)	O(logn)	O(n)
最优轮盘赌法	O(n2)O(n^{2})O(n2)	O(logn)	O(n)
Alias Method	O(n)	O(1)	O(n)

四、伪代码

参考文献

Alias method 别名算法
alias method
Darts, Dice, and Coins: Sampling from a Discrete Distribution
Alias Method离散分布随机取样（带python实现）