Bit Blasting Probabilistic Programs

新的概率编程语言HyBit

https://github.com/Tractables/Dice.jl/tree/hybit

https://arxiv.org/pdf/2312.05706v2

摘要

概率编程语言（Probabilistic Programming Languages, PPLs）是一种用于创建和推理概率模型的富有表现力的工具。不幸的是，当今的PPL对混合概率程序（即同时涉及连续和离散结构的概率程序）支持不够良好。在本文中，我们开发了一种新的近似推理算法，用于混合概率程序，该算法首先对连续分布进行离散化，然后在得到的程序上执行离散推理。我们的核心创新是一种称为“位爆破”（bit blasting）的离散化形式，它使用数字的二进制表示，使得一个包含2^个离散点的域可以简洁地表示为一个由poly()个布尔随机变量构成的离散概率程序。令人惊讶的是，我们证明了许多常见的连续分布可以通过“位爆破”进行离散化，而不会相对于显式离散化造成任何精度损失，并且还支持高效的概率推理。我们构建了一个名为HyBit的概率编程系统，专门用于处理混合程序，该系统采用位爆破方法，随后进行离散概率推理。我们通过实验验证了我们的方法相比于现有的基于采样的和符号推理方法的优势。

CCS概念：• 计算数学 → 概率表示；概率推理问题
附加关键词和短语：离散化、位爆破、概率推理
ACM参考格式：
Poorva Garg, Steven Holtzen, Guy Van den Broeck 和 Todd Millstein. 2024. Bit Blasting Probabilistic Programs. Proc. ACM Program. Lang. 8, PLDI, Article 182 (June 2024), 24 pages. https://doi.org/10.1145/3656412

1 引言
概率编程语言（PPLs）是一种用于创建和推理概率模型的富有表现力的工具。许多这样的模型本质上是混合的，即同时涉及连续结构（例如高斯分布）和离散结构（例如伯努利随机变量、if语句及其他控制流）。例如，在医疗诊断、基因表达和信息物理系统等应用中就会出现混合模型 [Chen et al. 2020; Lee and Seshia 2017]。

然而，今天的PPL并不能很好地支持混合程序。概率编程语言中的主要分析任务是概率推理，即计算某个事件根据程序定义的分布发生的概率。现有的推理算法通常采用采样方法来进行近似推理。一些方法，特别是哈密顿蒙特卡洛方法（Hamiltonian Monte Carlo），被用于Pyro和Stan等PPL中 [Bingham et al. 2019; Gorinova et al. 2021]，但它们不支持离散随机变量，而是要求将这些变量（手动或自动）边缘化出去。然而，这种方法在许多情况下会导致指数级爆炸，尤其是在离散变量数量较多时。其他基于采样的方法具有通用性，因此可以处理离散性，如重要性采样、马尔可夫链蒙特卡洛和序列蒙特卡洛等 [Koller and Friedman 2009]。然而，这些算法在面对由于离散性产生的多峰分布 [Yao et al. 2021] 以及那些以低概率事件为条件的程序时，常常表现不佳。

在本文中，我们提出了一种通过离散化来处理混合概率程序的新推理算法：我们将混合程序中的连续分布转换为离散分布。这产生了一个完全离散的概率程序，可以在其上使用现有的离散推理算法。离散化将连续分布近似为一系列区间，每个区间对应值落在该区间内的概率。先前的研究 [Albarghouthi et al. 2017; Beutner et al. 2022; Claret et al. 2013; Huang et al. 2021] 中已经使用过各种形式的离散化方法，但它们的规模都与区间的数量呈线性关系。这就带来了一个明显的权衡：需要足够多的小区间以避免过多的精度损失，但随着区间数量的增加，推理的成本迅速变得不可承受。

我们引入了一种新的离散化方法，称之为“位爆破”（bit blasting），这一术语借鉴自验证领域中的同名技术 [Bruttomesso and Sharygina 2009]。位爆破离散化的一个关键特性是，它仅需使用 poly(log) 个布尔随机变量即可表示 个区间的离散化。这是通过使用数字的二进制表示，并将离散化表示为基于这种二进制表示的离散概率程序来实现的。乍一看，这种简洁的表示方式似乎会丢失太多精度，从而无法成为一个可行的策略，但我们通过理论和实验证明事实并非如此。

首先，我们证明了一大类常见的连续密度函数可以被“位爆破”安全地处理，也就是说，其精度与朴素离散化相比没有任何损失。表1列出了属于这一类的示例分布；我们将整个这类分布称为混合伽马分布（mixed-gamma distributions）。

例如，考虑将0到1之间的连续均匀分布进行离散化。朴素地将其离散化为2³²个区间需要枚举2³²个值。相反，这个分布可以以二进制形式表示为一个由32个flip(0.5)形式的伯努利随机变量组成的元组，即掷硬币，其结果等概率为0或1。这一观察并不新颖，类似的结果也适用于指数分布 [Marsaglia 1971]。然而，其他混合伽马密度函数的位爆破离散化是新颖的。此外，与均匀分布和指数分布不同，这些离散化并不是简单地定义为一组独立的伯努利随机变量，而是需要基于这些变量构建完整的离散概率程序。

一种简洁的表示方式并不一定意味着推理效率高，因为一般情况下推理本身就很困难。作为我们的第二个贡献，我们证明了位爆破后的混合伽马分布在不仅是正确且简洁的前提下，还支持在位精度数量上的多项式时间推理。具体来说，我们证明了使用知识编译（knowledge compilation）方法进行离散概率推理 [Chavira and Darwiche 2008; Chavira et al. 2006; De Raedt et al. 2007; Fierens et al. 2015; Holtzen et al. 2020]，该方法将推理归约为布尔公式上的加权模型计数（weighted model counting），对于位爆破后的混合伽马分布具有上述性质。因此，对混合伽马分布进行位爆破后，再通过知识编译进行离散推理的过程，其运行时间被保证为离散化所用位宽的多项式级别。

第三，我们利用上述理论成果设计了一个新的概率编程系统HyBit，它通过位爆破实现对混合概率程序的非随机近似推理。HyBit是一种离散概率语言，支持定点二进制数及其相关的算术运算，并作为嵌入式的领域特定语言在Julia中实现 [Bezanson et al. 2017]。HyBit的API允许用户生成任意连续分布的位爆破定点近似。混合伽马分布可以通过其正确的位爆破离散分布来表示。对于其他分布，API允许用户使用分段离散近似，其中每一段本身就是一个位爆破后的混合伽马分布。

HyBit利用知识编译对给定的离散概率程序执行精确推理。在最坏情况下，通过对混合概率程序进行位爆破推理的时间复杂度可能是位数的指数级，但我们通过实验验证了我们方法的优势。我们展示了HyBit在一套全面的混合概率程序基准测试中，表现优于现有的基于采样和符号推理的方法。

总体而言，本文提出了以下贡献：

  1. 我们在第2节中阐述了混合概率程序推理中的挑战。

  2. 我们在第3节中提出了一种新的离散化方法，称为位爆破，其特点是离散区间数量的简洁性。

  3. 我们提出了一类连续分布，即混合伽马分布，在这类分布中存在正确的位爆破表示。我们在第3节中形式化了这种构造并证明了其性质。我们进一步证明，基于知识编译的推理在这些分布上能够随位宽呈多项式级扩展。

  4. 我们在第4节中描述了HyBit概率编程语言及其通过位爆破实现的新推理算法。

  5. 在第5节中，我们将HyBit与其他PPL在来自现有文献的基准测试中进行了实证比较。我们还分析了HyBit在其超参数（即位数和分段数）下的行为表现。

HyBit的项目地址为：https://github.com/Tractables/Dice.jl/tree/hybit 。
本文的完整版本可在arXiv上获取，编号为[Garg et al. 2023]，其中包含完整的证明过程。

2 动机示例

本节通过三个示例说明对混合概率程序进行推理的挑战。首先，我们展示一个来自计算生物学的例子，其中包含固有的逻辑结构。接下来，我们展示文献中的一个例子，由于离散控制流的存在，其后验分布呈现多模态。最后，我们通过共轭高斯分布展示低概率观测值的一个例子。随后，我们研究了包括 HyBit 在内的各种推理算法的表现。这些示例展示了 HyBit 相对于其他近似推理算法的优势。

我们在第5节中提供了与多个近似和精确基线方法的详细比较。

2.1 逻辑结构

我们展示了一个来自计算生物学的简化示例，该示例将基因表达与血糖水平联系起来。图1展示了对应的概率程序，其中的任务是：在已知患者的血糖水平的前提下，更新关于该患者体内某个基因出现概率的信念。

图1中概率程序的前四行使用beta分布作为一般人群中个基因出现的先验概率。语法flip()表示一个成功概率为的伯努利随机变量。在第5行，程序使用reduce(|, gene)语法来表示表达式⨀ₜᵢ₌₁ gene[i]，换句话说，如果至少有一个基因被表达，则认为该患者患有糖尿病。

接下来是关于患者血糖水平的多次测量读数。对于每一次读数，程序首先根据患者是否患有糖尿病定义了一个随机变量来表示血糖水平。然后我们使用observe(y, v)语法来对随机变量取值为进行条件约束——在程序中，这用于对患者实际的血糖读数进行建模。最后，在第12行，程序查询了第一个基因出现的后验分布的期望值。

图2展示了使用不同的推理算法对该程序进行推理的结果，在20分钟的超时限制下，随着基因数量（T）的增加，各算法的表现情况。

Stan 使用的是哈密顿蒙特卡洛方法（Hamiltonian Monte Carlo），这种方法不直接支持离散随机变量。取而代之的是，需要通过手动或自动的方式将这些变量边缘化出去，例如使用变量消除（variable elimination）[Gorinova 等 2020]。如图所示，当基因数量超过15个时，Stan 就会超时——因为其计算复杂度随 T 呈指数级增长。

同样的指数爆炸问题也困扰着 GuBPI [Beutner 等 2022]，该方法采用符号求值与离散化相结合的方式来计算上下界。正如图中所示，通用采样方法如带有 Metropolis Hastings 核的马尔可夫链蒙特卡洛（WebPPL MH）和序列蒙特卡洛（WebPPL SMC）能够扩展规模并提供合理的准确性。

精确推理策略 Psi [Gehr 等 2016] 在程序设计避免产生大规模离散状态空间的前提下，也能很好地应对离散随机变量数量的增长。AQUA 对混合程序进行离散化处理 [Huang 等 2021]，但无法支持本例中的这个程序。

我们的系统与方法 HyBit 能够扩展到 50 个基因，并且在所有近似推理算法中具有最小的绝对误差。当该程序以 HyBit 编写时，它通过其离散的位级抽象进行表示。用户在编写程序时，可以使用 HyBit API 将所有的连续分布（具体来说是图1中的第2、6、9行）替换为其位爆破后的离散近似版本。

结果是，我们现在得到了一个离散程序，其中包含对布尔值和定点数的分布。请注意，现在第8行的 observe 条件是基于这样一个事实：blood_sugar1 的离散分布取值落在对应数值79的区间内。

在本次实验中，我们使用了 25位 的位宽——每个连续分布被离散化为一个由25位组成的元组所表示的程序，这些位被解释为一个定点数。此外，我们使用了 4096个分段 来近似这些连续分布，每个分段本身是一个位爆破后的指数分布。

如果采用朴素的离散化方式，使用25位将是非常缓慢以至于不可行的，因为它会产生 2²⁵ 个区间，即约1.34亿（134M）个区间。然而，我们通过位爆破后的程序仅使用了 53K 次掷硬币 （即布尔随机变量）来表示它们。

此外，基于知识编译的推理方法 [Fierens 等 2015; Holtzen 等 2020] 自动识别并利用了程序逻辑结构中的条件独立性，从而帮助推理过程实现良好的扩展性。

有关此实验的更多细节，请参见附录。

2.2 处理多模态性

本节展示了一个多模态分布的例子，用以突出混合概率程序推理中的另一个挑战。多模态分布在多个高峰之间被低概率区域分隔，这类分布在诸如传感器网络定位、宇宙学等多个应用领域中非常常见 [Shaw 等 2007; Tak 等 2018]。我们改编了文献中的一个例子 [Yao 等 2021]，如图3所示。

图3中所示的概率程序对现有的概率推理方法来说具有很大挑战性。第3行的datapts包含九个数据项，其中三分之二为5，其余为−5。这导致了(₁, ₂)的后验分布围绕(5, −5)和(−5, 5)呈现双峰分布。

随着数据点数量按相同比例增加，(₁, ₂)的后验会收敛到(5, −5)。然而，在数据点数量有限的情况下，₁的后验在5和−5附近呈现双峰分布。

多个模态的存在给基于采样的算法带来了挑战，因为它们容易陷入某一个模态中无法跳脱。具体来说，使用Metropolis Hastings核的WebPPL MCMC算法，以及使用HMC的Stan和WebPPL，最终都会任意地停留在其中一个模态中，并未能探索另一个模态。

图4a和图4b分别展示了使用WebPPL MCMC和Stan HMC得到的结果，其中两次不同的运行分别被困在了不同的模态中。另一方面，HyBit对其离散抽象执行精确推理，因此能够全局探索整个分布，从而识别出两个模态。序列蒙特卡洛（SMC）（图4c）、Psi 和 GuBPI 同样能够进行全局探索，因此也能找到这两个模态。最后，为了应对直接离散化带来的计算挑战，AQUA通过调整其离散化区间，专注于高概率区域，结果只识别出了更高概率的那个模态（图4d）。

2.3 处理低概率观测

本节展示了一个具有低概率观测值的共轭高斯示例（见图5）。在图5中，在第2行和第3行对低概率数据进行条件化之后，程序查询了随机变量mu的后验分布。

为什么低概率观测值难以处理？直观上来说，通用的采样算法从先验分布开始采样，但在寻找具有显著权重的样本时会遇到困难。只有在进行大量采样之后，这些算法才能逐渐接近真实的后验分布。

另一方面，HyBit 通过对位爆破抽象执行精确推理，能够全局探索后验分布的定义域，因此不受这一问题的影响。

图6绘制了真实的先验分布和后验分布，以及来自不同推理算法的结果。

对于基于采样的算法——使用Metropolis Hastings核的MCMC 和 SMC，我们运行了相应的WebPPL程序 [Goodman and Stuhlmüller 2014]，并在之后获得了1000个样本并进行了绘图。重要性采样（importance sampling）算法未能为此程序获得任何具有非可忽略权重的样本。

这些采样器正在将后验分布向真实后验靠近，但需要更多的样本才能实现这一点。即使分别采样了约1600万次和65,000次之后，MCMC和SMC所获得样本的期望值仍然存在绝对误差 分别为 0.549798 和 1.520776 。

GuBPI 报告了每个区间上概率的上下界，并产生了 2.33 的绝对误差。

另一方面，HyBit 所得到的后验分布与真实分布完全重合。

Stan HMC 在处理低概率观测方面表现良好，并获得了很高的精度。

Psi 同样得到了后验分布的精确符号表达式。

最后，AQUA 所报告的后验均值的误差为 5.66 ，因为它未能对的先验做出任何更新。

3 位爆破：核心见解

为了在混合概率程序的离散结构上实现推理的扩展性，我们需要一种将离散性作为头等公民对待，并能有效消除连续结构的算法。本节定义了“位爆破”（bit blasting）的语义概念，并将其设定为具有优良性质的离散化的一种特殊情况。接着，我们为常见类别的连续分布提供了位爆破函数。我们所提供的离散化技术具备正确性 （精度可达位）、简洁性 ，并适用于高效的推理。

3.1 离散化与位爆破

在概率论的标准术语中 [Rosenthal 2006]，一个概率空间 (Ω, Σ, ) 由样本空间 Ω、Ω 上的 -代数 Σ，以及定义在 Σ 上的概率测度 组成。

广义而言，一个离散化函数 的输入是一个这样的概率空间 (Ω, Σ, )，输出则是一个离散概率空间(Ω, Σ, )，其中 Ω 是一个可数集合。

我们将研究一种更具体的离散化概念：该离散化以有限区间上的连续分布 为输入，并输出一个在 2ᵇ 个点上的离散分布 ，其中 为位数。

形式上，设 [, ) 是一个区间，其中 , ∈ ℝ。对于输入概率空间，我们使用 B([, )) 来表示区间 [, ) 的子集所构成的 Borel -代数 [Rosenthal 2006]。

对于输出概率空间，我们用 P() 表示集合 的幂集（即其对应的 -代数）。

此外，我们假设样本空间的离散化方式如下所示。

在我们定义一个 位的位爆破函数 之前，我们需要先确定一个通用的离散概率分布表示方式。为此，我们定义了一个称为离散概率闭包 （discrete probabilistic closure）的概念，类似于概率图灵机 [Arora and Barak 2006]。

每个概率闭包是一个从一组有偏硬币掷（biased coin flips）到某个离散集合的确定性函数 。该函数通过与硬币掷结果相关的概率，在其输出上诱导出一个概率分布。此外，它还包含一个接受布尔公式 （accepting Boolean formula），用于处理观测，并限制输入硬币掷所能取的值的集合。

下面我们将给出其形式化定义：

请注意，正如示例3所解释的那样，任何在 2ᵇ 个值上的离散分布 都可以表示为一个离散概率闭包 (, , )，其中 || = 2ᵇ − 1。

Dice [Holtzen 等 2020] 和 Problog [Fierens 等 2015] 就是直接符合离散概率闭包这一范式的概率编程语言（PPL）的例子。

我们希望离散概率闭包能够更加简洁 ——其大小应为精度位数 的多项式级别。为此，我们定义了位爆破函数 。

定义5（位位爆破函数） ：一个 位位爆破函数 [.]ᵇ 是一种 位的离散化函数，它输出一个离散概率闭包 (, , )，该闭包所使用的布尔随机变量数量是位数 的多项式级别，即 || ∈ (poly())。

根据定义3、定义4和定义5可知，对于任意整数 > 0，一个 位位爆破函数对于给定的概率空间 ([, ), B([, )), ) 是正确 （sound）的，当且仅当：

我们想指出的是，bitblast_unif 仅使用了 2次硬币掷 就表示了一个在 4个值 上的分布。它可以推广为使用 次硬币掷 来表示一个在 2ᵇ 个值 上的均匀分布。

相比之下，naïve_unif 使用了 3次硬币掷 ，并会推广为使用 2ᵇ −1 次硬币掷 。

因此，bitblast_unif 是一个适用于 位位爆破函数（-bit blasting function）的有效输出。

在下一节中，我们将详细介绍指数分布和混合伽马分布的 位位爆破函数的具体实现。

3.2 具体的位爆破函数：预备知识

接下来，我们的目标是为混合伽马分布 （mixed-gamma distributions）提供一个具体的、满足正确性要求的位爆破函数实例 。混合伽马分布的概率密度函数定义如下。

为了记号上的方便，我们将连续分布限制在单位区间 内，从而得到定义在 位单位区间 上的离散分布。随后我们会将我们的方法推广到任意有限区间，以构建基于位爆破的概率编程系统 HyBit 。

为了描述我们对位爆破函数的构造，我们使用了 Dice [Holtzen 等 2020]。Dice 已经能将其程序编译为带权布尔公式 （通过 ⇝ 判断），这些公式符合离散概率闭包 的定义。³
这使得我们只需定义一个从概率密度函数到 Dice 程序的 ⇝ᵇ 判断，即可指定一个位爆破函数。

Dice 还定义了一个分布语义函数 J.Kᴰ：p → V → [0, 1]，它以一个 Dice 程序 p 作为输入，并输出一个归一化的概率分布（表示为从值集合 V 到概率的一个函数）。我们在后续使用函数 J.Kᴰ 来论证我们构造的正确性。

有关 Dice 的语法和语义的更多细节，请参见附录。

混合伽马密度函数的编译判断形式为 Υ ⇝ᵇ p ，其中 p 是 Dice 程序。

我们进一步给出以下定义：

  • 密度函数与 Dice 程序之间的 等价性

    （-equivalence）

  • Dice 程序的 简洁性

    （-succinctness）

我们定义 -简洁性 （-succinctness），要求 Dice 程序 p 中使用的硬币掷次数与位数 成线性关系。请注意，-简洁性 所施加的条件比 位位爆破函数 所需的条件更为严格（后者只要求使用 poly() 次硬币掷）。这意味着，如果我们为某个混合伽马密度函数建立了一个满足 -简洁性 的判断，那么我们就可以为该分布构造一个合法的

3.3.1 指数分布 ₀,

让我们首先考虑均匀分布 （₀,₀），它是指数分布的一个特例。

如果我们使用 位 对均匀分布进行位爆破，并将其划分为 2ᵇ 个区间，那么我们会得到一个在 [0, 1]ᵇ 上定义的离散分布 ᵇ，它包含 2ᵇ 个离散点，每个点的概率为 1 / 2ᵇ。

一种直接的离散化策略是枚举这 2ᵇ 个值，这需要使用 2ᵇ −1 次硬币掷。

但我们可以用更高效的方式实现同样的效果：使用一个由 位 构成的元组，其中每一位都是一个无偏的硬币掷 flip(0.5)。

之所以可以这样对均匀分布进行位爆破，是因为二进制数字之间具有独立性 。这一策略同样可以推广到一般的指数分布。这一观点早在统计学的经典论文中就被提出 [Marsaglia 1971]。我们通过以下规则形式化这一思想。

我们在附录中提供了上述引理的详细证明。在本文的其余部分，我们将指数分布视为构建其他分布的基本元素，因为只有它们具备二进制位之间相互独立 这一特性。然而，正如我们接下来将展示的那样，对于其他分布，构造正确的位爆破函数 仍然是可能的。

3.3.2 伽马分布 ₁,

为了为 ₁, 构造一个正确的位爆破函数，我们提出一个关键的数学洞察。

考虑图8a中的程序。连续随机变量 X 和 Y 分别服从均匀分布（₀,₀）和指数分布（₀,）。该程序返回在条件 Y < X 下 X 的新分布。

结果表明，这个后验分布是一个特定的伽马分布 ₁,。

下面我们展示相应的推导过程，其中 pdf 表示概率密度函数。

如果我们使用 位 对图8a中的程序进行离散化，会发生什么？

我们会得到图8b中的程序，其中每个连续随机变量都被其对应的位爆破版本 所替代（例如 X 被替换为 Xᵇ，依此类推）。我们已经知道，对于均匀分布 和指数分布 ，存在满足 等价性 的 Dice 程序。但问题是：其他构造 （如条件语句）会怎样呢？

正如图8d所示，observe(Yᵇ < Xᵇ) 相较于其连续版本 observe(Y < X) 会产生误差。好消息是，我们可以通过以下公式来量化并处理这种误差 ：

规则 Expo1 和 Trans-expo1zero （见附录）捕捉了上述直觉。

其中，unifObs(y, b) = p 是一个辅助判断，程序 p 构造了一个均匀分布 ，并通过 observe 对其进行条件约束，使其取值小于 y 。

3.3.3 广义伽马分布 ₐ,

前一小节展示了通过对不等式 (Y < X) 进行条件约束，如何在 ₀, 的基础上引入一个线性因子，从而得到 ₁,。请注意，无论随机变量 X 的初始概率密度函数是什么形式，对 (Y < X) 进行条件约束都会引入一个线性因子 。也就是说，如果 X 的初始概率密度为 ()，那么如下概率程序的输出密度将是 ()。这意味着：如果我们能够对 () 进行位爆破，我们也就能够对 () 进行位爆破。但我们仍需考虑由 observe(Yᵇ < Xᵇ) 所带来的误差，即 Pr(ᵇ | ᵇ == ᵇ, < )。结果表明，这个修正项是一个伽马分布的混合分布 ，而这类分布也可以被正确地进行位爆破 。

我们在附录中提供了相应的判断规则以及以下引理的证明。

3.3.4 伽马分布的混合形式 ∑ᵢ ᵢₐᵢ,ᵢ

由于广义伽马密度函数 ₐ, 可以被进行位爆破，因此混合伽马密度函数 也可以被正确地进行位爆破。

我们对每一个单独的广义伽马密度进行位爆破，然后使用 if-then-else 结构来构建它们的混合形式，如下所示。

之前的研究 [Holtzen 等 2020] 定义了判断 ⇝，它以 Dice 程序 p 作为输入，并输出一个带权布尔公式 (, , )，该公式与离散概率闭包 （定义4）的定义一致。

并且由于根据定理8 ，对于所有混合伽马密度函数，判断 ⇝ᵇ 是 -简洁的 （-succinct），因此 ⇝ 总是输出一个使用 poly() 次硬币掷的 。

因此，复合判断 ⇝ ◦ ⇝ᵇ 构成了一个合法的 位位爆破函数 。

此外，之前的研究 [Holtzen 等 2020] 还证明了将程序编译为带权布尔公式的正确性，即其结果与 Dice 程序的语义是一致的。

结合这一结论和定理7 ，我们可以得出：复合判断 ⇝ ◦ ⇝ᵇ 是一个正确的 位离散化函数 。

详细的证明可以在附录中找到。

3.3.5 示例：拉普拉斯分布

前面的章节描述了当混合伽马密度被限制在单位区间时，如何对其进行位爆破。那么，对于那些被平移 或缩放 到其他有限区间上的分布，又该如何进行位爆破呢？

我们通过拉普拉斯分布 （Laplace distribution）来解释这一点。

拉普拉斯分布有两个参数：

  • 位置参数（location）

  • 尺度参数（scale）

其概率密度函数如下所示，其中 ∈ ℝ：

我们考虑在区间 [ − , + ) 上被截断的拉普拉斯分布。我们假设 是一个合适的 2 的幂次 ，这样与 的乘法运算（记作 × p）可以简化为小数点移位操作 ，同时假设 是一个可以用 位 精确表示的数，以便对 p 进行精确移位。

首先，我们生成在宽度为 而非 1 的区间上的缩放后的指数分布：

3.4 位爆破如何助力推理？

我们已经展示了针对混合伽马分布 的正确位爆破函数 。但这种位爆破方式如何帮助对包含这些分布的概率程序进行推理呢？

为了回答这个问题，我们聚焦于一种特定的推理策略——知识编译 （knowledge compilation）。

我们首先介绍有关知识编译的一些必要预备知识，然后论证通过 ⇝ᵇ 所获得的程序为何在知识编译中具有高效性。

基于知识编译的方法 [Fierens 等 2015; Holtzen 等 2020] 用于精确的离散概率推理 ，其核心思想是将离散概率程序编译为带权布尔公式 ，并使用有序二元决策图 （Ordered Binary Decision Diagrams, OBDDs）来表示这些公式。

  • 当程序返回一个 单个布尔随机变量 时，对应的 OBDD 是 单根 的；

  • 当程序返回一个 布尔随机变量的元组 时，对应的 OBDD 是 多根 的。

通过为程序中的所有硬币掷（即带有权重的有偏硬币掷）设定一个取值，并沿着这些取值在 OBDD 中进行遍历（实线表示 true，虚线表示 false），我们可以到达对应每一位输出值的终端节点。

加权模型计数 （Weighted Model Counting）操作可以计算每个输出位到达“1 终端”的概率。这是一个动态规划算法，其运行时间与 OBDD 的大小成线性关系。

一个布尔公式 对应的 OBDD 大小，记作 OBDD()，是指该 OBDD 中的节点数量。

因此，如果我们能为某个分布构造出更小的 OBDD 表示，就可以更高效地将其与其他结构组合进一个离散概率程序中。

接下来我们正式讨论并证明：通过判断 ⇝ᵇ 所获得的每一个程序都可以被编译为一个带权布尔公式，进而被编译为一个多根 OBDD ，其大小随着位数呈线性增长 ，而不是最坏情况下的指数增长。

回顾一下，Dice 程序会被编译为一个带权布尔公式 (, , )：

  • 是一个（或一组）布尔公式，对应程序的返回值；

  • 是一个接受布尔公式，用于编码观测条件；

  • 是一个权重函数，表示各个硬币掷的概率。

定理10 ：对于所有 Υ, p, , , ，存在某个常数 ，使得对任意位数 ，如果 Υ ⇝ᵇ p 且 p ⇝ (, , )，那么存在一个布尔随机变量在 中的变量顺序 Π，使得：

OBDD() + OBDD() ≤

我们现在为上述定理的证明提供直觉解释。

请注意，在通过判断 ⇝ᵇ 所获得的程序中，只有两个构造依赖于位数 ：

  1. 指数分布 ₀, 的构造
  2. 通过对不等式的条件约束

    （即通过 unifObs(y, b) 对指数分布与均匀分布之间的关系进行建模）

接下来我们将说明这些构造所对应的 OBDD 大小如何随着位数 线性增长 。

3.4.1 指数分布

在图9中，我们对一个指数分布进行了 3位爆破 ，也就是说，我们在8个值上得到了一个离散的指数分布：
{0, 0.125, 0.25, ..., 0.875}。

图中展示了一个三根节点的 OBDD ，其中每个根节点分别标记为 ₁、₂ 和 ₃，代表返回的3位值中的每一位。

假设我们将对应节点 ₁ 的硬币掷设为 true，那么对于 ₁ 来说，我们会到达终端 1，并将其值设为 1。

加权模型计数 （WMC）操作会计算 ₁ 取值为 1 的概率为 6.14×10⁻⁶，因为节点 ₁ 为真的概率是 6.14×10⁻⁶。

同样地，WMC 也可以用于该 OBDD 的其他根节点。

由于每一位只需要一个 OBDD 节点，因此整体的 OBDD 大小随位数 线性增长 。

又因为 WMC 的运行时间与 OBDD 的大小成线性关系，所以对指数分布的推理时间复杂度为 O()。

附录中还给出了一个关于均匀分布的 OBDD 示例。

对于所有通过判断 ⇝ᵇ 得到的程序 p，如果 p ⇝ (, , )，那么：

  • 是表示程序返回值的一组布尔公式。

我们论证说，p 的返回值始终是一个指数分布的混合体 ，因此其对应的 OBDD 大小也与位数呈线性关系。

3.4.2 对指数分布与均匀分布之间的不等式进行条件约束

判断 ⇝ᵇ 的规则使用了辅助判断 unifObs(y, b) = p，用于对均匀分布与指数分布的二进制表示之间的不等式 进行条件约束。

由于唯一依赖于位数的构造（即指数分布和不等式）都随着位数呈线性增长，因此定理10在直觉上是成立的。我们在附录中提供了形式化的证明。

4 HyBit：一个概率编程系统

前一节介绍了如何对混合伽马分布 进行位爆破 （bit blasting）。我们进一步利用这一技术构建了一个用于混合概率程序 的概率编程系统 HyBit 。

本节将描述其语法与实现 ，并详细阐述两个重要方面：

  1. 连续分布的分段近似方法
  2. 采用二进制表示的优势
4.1 HyBit —— 语法与实现细节

我们构建了一个概率编程系统 HyBit ，它围绕混合伽马密度函数的正确位爆破 （sound bit blasting）以及其他连续分布的近似位爆破 展开。

HyBit 是作为一门浅嵌入式领域特定语言 （shallow embedded DSL）在 Julia [Bezanson 等 2017] 中实现的。

图10给出了 HyBit 表达式的核心语法。它支持以下功能：

  • 对布尔值的概率分布建模： flip （即伯努利分布）

  • 对定点数的概率分布建模： general_gamma 和 bitblast

  • 布尔运算：¬, ∧, ∨

  • 算术运算：+, −, *, /, %, <, ==

  • 概率条件建模中的硬性观测： observe

对于图10中列出的所有构造，HyBit 会执行一种非标准解释执行 （non-standard execution），并将其编译为 OBDD （有序二元决策图），以进行概率推理。

由于 HyBit 是作为 Julia 的一个库实现的，因此程序员还可以使用 Julia 提供的语言结构，例如（有界）循环、元组和函数。

举个例子，可以在循环体中使用 HyBit 构造配合 Julia 的 for 循环来构建概率模型。

HyBit 作为一个开源项目提供，并附带了丰富的示例集。

图11详细展示了 HyBit 的 API。

  • DistFix{W, F} 表示一个定点数类型，总位宽为 W，其中小数点后占 F 位。

  • 函数 general_gamma 可对指定的广义伽马密度函数进行正确的位爆破 ，将其映射到给定的定点数精度 W 和 F 上。

  • 通过在广义伽马密度上使用 if-then-else 结构，可以实现对混合伽马密度函数 的正确位爆破。

  • 函数 bitblast 用于对任意连续分布进行分段近似 的位爆破，用户可通过参数 W、F 和 pieces 来指定使用的位数和分段数量。

  • API 还允许用户选择用于分段近似的离散分布类型：线性分布 或 指数分布 。

  • 线性分段的参数（斜率）和指数分段的参数（）会自动选择，使得第一个区间与最后一个区间的概率比值 与朴素离散化保持一致。

  • 最后，API 还提供了用于查询随机变量的概率分布 、期望值 和 方差 的函数。

接下来的小节将更详细地描述分段近似方法 ，以及期望值与方差的计算方式 。

4.2 分段近似

尽管混合伽马分布能够涵盖许多常见的自然分布,但仍有一些常见分布不在其中,例如高斯分布 。

目前仍是一个开放问题：高斯分布是否可以被正确地进行位爆破（更不用说编译成紧凑的 OBDD）。对于这类分布，我们可以使用分段近似 （piece-wise approximation）方法来处理。

设 是定义在区间 [, ) 上的一个任意连续概率分布。

要使用具有 个分段的分布对 进行位爆破，我们通过一个由 个离散概率分布 组成的混合分布来近似 ，这些分布分别对应互不重叠的子区间。

对于每一个分段，我们构造一个经过平移和缩放的位爆破混合伽马密度函数 ，然后将它们组合成一个混合分布。

请注意，由于每个分段使用了 O() 次硬币掷，因此具有 个分段的分段分布总共使用 O() 次硬币掷。

第5节展示了这种策略在实验上的优势。

这种使用线性分段 或指数分段 的近似方法，可以很方便地通过 HyBit 中提供的 bitblast API 实现（见图11）。

图12 展示了使用 2、4、8 和 16 个分段 对高斯分布进行位爆破的结果，其中每个分段都是一个位爆破后的指数分布。这为用户提供了在精度与性能之间的传统权衡 。

我们将在第5.2节中对此进行更详细的说明。

4.3 二进制表示的优势

除了位爆破所提供的简洁性 之外，二进制表示 在概率推理中还具有重要的额外优势。

首先，许多混合概率模型涉及对连续随机变量进行算术运算 。由于我们使用的是定点数的二进制表示，像 +、*、/、< 这样的算术操作会被编译为作用于二进制数上的布尔公式（类似于计算机体系结构中的 ALU 电路）。

这种表示方式使得概率推理（特别是我们采用的知识编译方法 ）能够识别并利用算术运算中存在的结构信息，例如在一个计算过程中各个输出位之间的条件独立性 。

最近的研究 [Cao et al. 2023] 描述了这类编译过程，并通过实验证明了其在整数运算中的优势；而 HyBit 则将这些优势扩展到了定点数 的计算上。

此外，二进制表示 也使期望值和方差的计算更加高效。

对一个定义在 2ᵇ 个值上的分布来说，朴素地计算其期望值和方差需要分别计算这 2ᵇ 个值的概率。

而使用按位表示（bitwise representation），我们只需计算每个比特位的概率，从而实现了指数级的效率提升 。

需要注意的是，在最坏情况下，对于任意的混合概率程序，获得其对应的二进制表示的 OBDD 本身可能是关于位数呈指数增长 的。

但对于混合伽马分布类 ，这种转换的复杂度是线性的 （见定理10）。

我们在以下两个定理中形式化了期望值与方差的计算方式，并在附录中提供了相应的证明。

定理12 ：设 D 是定义在区间 [0, 2ⁿ) 上的一个离散概率分布，该分布以 位二进制形式表示为 (ₙ, ₙ₋₁, ..., ₁)，则 D 的期望值可以利用期望的线性性质如下计算:

5 实验评估

我们评估了在真实概率程序中使用位爆破方法的实用性。为了实现这一目标，我们进行了相关实验，以探讨以下问题：

  • Q1

    ：HyBit 相较于现有的推理算法表现如何？见第5.1节

  • Q2

    ：分段近似的效果如何？见第5.2节

5.1 与现有推理算法的比较 5.1.1 近似推理算法

我们将 HyBit 与两类近似推理算法进行了对比评估。

采样方法
我们对比了以下代表性算法：

  • WebPPL 的拒绝采样（rejection sampling）

  • 使用 Metropolis Hastings 核的 MCMC 采样

  • SMC 采样

  • Stan 的 HMC 采样

离散化方法
我们还与以下两类基于离散化的推理算法进行了比较：

  • AQUA

  • GuBPI [Beutner 等 2022; Huang 等 2021]

比较不同概率编程系统的性能是一项具有挑战性的任务，因为性能直接受到程序结构的影响。

为公平起见，我们在每个系统中编写了等效的基准程序，并尽最大努力对其进行了优化。

本节及其后续章节中的表格报告了随机算法在10次运行中的平均绝对误差 ；对于其他推理算法，则报告单次运行的结果。

所有实验均采用单线程 方式执行，运行环境为配备 2.4 GHz CPU 和 512 GB RAM 的服务器。

表2报告了 HyBit 与其他近似推理算法在性能评估中的对比结果。

我们选取了所有 Psi [Gehr 等 2016] 曾评估过的混合与连续基准程序 ，并额外补充了一些来自现有研究 [Huang 等 2021] 的相关基准程序。

我们尽最大努力通过解析方法 或使用计算机代数系统 为这些基准程序计算出真实值 （ground truth）。

在本评估中，我们仅包含那些我们可以可靠地获得真实值 的基准程序。

对于所有基准程序，我们都报告了各算法相对于真实值的绝对误差 。对于返回非布尔值的基准程序，我们计算了每种方法所得到的期望值的绝对误差 。

我们报告的是各推理算法在20分钟超时限制内 所能达到的最小误差 。

在所有基准测试中：

  • HyBit 将混合伽马分布替换为其 正确的位爆破分布 ；

  • 其他分布则被替换为 线性分段近似 版本，即 ₁,₀ 分布。

每个基准程序所使用的位数 （bits）和分段数 （pieces）也在表2中列出。

为了在这些基准程序上运行 Stan，我们使用了 SlicStan [Gorinova 等 2020] 来生成经过离散随机变量边缘化处理 后的 Stan 程序。

对于所有 WebPPL 基线方法，我们在20分钟内对所有采样算法均使用默认设置 ，并尽可能多地进行采样。

如表2所示，使用位爆破的 HyBit 在所有基准测试中与现有方法表现相当，甚至在 19 个中的 11 个 基准上表现更优。对于其余 8 个基准，HyBit 的准确率也非常接近。

  • AQUA

    仅在 4 个基准上表现更好；

  • GuBPI

    未能获得良好的准确率，这主要是因为其采用的枚举式离散化方法在高精度下扩展性较差。

WebPPL 和 Stan （通过 SlicStan 实现自动边缘化）支持大多数基准程序，但在限定时间内未能达到较高的准确率。这是因为基于采样的算法具有随机性，在有限时间内无法从真实的后验分布中获取足够多的样本。

5.1.2 精确推理算法

表3比较了 HyBit 与一个使用代数方法进行精确推理 的概率编程系统 Psi [Gehr 等 2016]。

我们尽最大努力对基准程序进行了优化翻译以适配 Psi 的运行环境。Psi 通常会输出一个符号表达式 ，我们需要将其输入 Mathematica 进行进一步简化。

然而，这些代数表达式的计算和简化并非易事：

  • Psi 在其中 6 个基准上超时；

  • Mathematica 也无法简化其中 4 个基准的结果。

相比之下，HyBit 能够处理全部 19 个基准程序 ，因为它将计算归约为对布尔随机变量的离散推理，并对推理查询进行近似处理。

5.2 分段近似的效果如何？

我们分析了在使用不同数量的分段对连续分布进行近似时，性能与准确率之间的权衡 。

图14 展示了在 四个基准程序 上，随着分段数的增加，运行时间与准确率的变化趋势 ，并且展示了在不同位宽（bitwidth）下的表现。

  • 随着线性分段数的增加，运行时间 先下降后上升 。

  • 准确率则如下面四个子图所示， 随着分段数的增加而提升 。这是因为当我们增加分段数量时，连续分布被更精确的位爆破分布所替代。

  • 然而，在达到某个“最佳点”之后，准确率的提升是以运行时间增加为代价的。

附录中提供了更多实验，进一步证明了相较于基于中心极限定理 的方法，使用分段近似 更具优势。

6 相关工作

概率编程一直是语义和推理两个研究方向上的热门领域 [Dahlqvist 等 2023; Milch 等 2005]。

本节将 HyBit 与相关工作进行对比。

总体而言，HyBit 的关键创新在于提出了一种位爆破方法 ，用于对混合概率程序进行简洁的离散化表示 。这一方法显著区别于现有工作。

离散化方法

一些早期的方法通过对连续或混合概率程序进行离散化，并通过枚举所有离散化值 来估计后验分布 [Huang et al. 2021]，但这种方法在很多情况下无法扩展到足够高的精度。

一种早期的离散化技术也使用了二进制表示 [Claret et al. 2013]，但其方法并不是我们所说的“位爆破”（bit blasting），因为它并不简洁，且最终生成的表示形式仍与离散点的数量成正比。

最近的一项工作则使用离散化来为概率程序的后验分布提供上下界估计 [Beutner et al. 2022]。

混合概率程序的推理算法

还有一些研究专门针对混合概率程序：

  • Leios

    [Laurel and Misailovic 2020] 通过对混合程序进行 连续化处理 ，以利用现有连续推理算法的能力。

  • HyBit

    则采取相反的策略，对混合程序进行 离散化处理 ，这有助于在面对具有复杂离散结构的混合程序时提升推理的可扩展性。

  • SPPL

    通过将混合程序转换为特定的表示形式来进行推理 [Saad et al. 2021]，但这种表示方式限制了所能支持的混合程序类型。例如，SPPL 不支持对连续随机变量进行算术运算，而 HyBit 可以做到这一点。

  • 此外，概率逻辑编程语言也被扩展用于支持混合模型，例如通过使用 区间轨迹 （interval traces） [Gutmann et al. 2011]。

代数方法

一些 PPL 的推理算法通过生成闭式代数表达式 来编码概率分布，并使用符号计算技术执行精确推理 [Gehr et al. 2016; Hur et al. 2014; Narayanan et al. 2016]。

然而，这些系统在表达能力和可处理程序范围上都存在一定的限制，如表3所示。

基于路径的推理算法

PPL 中常见的一类推理算法是操作语义式的 （operational）：它们通过使用随机变量的具体取值来记录程序的执行路径。

这类方法包括：

  • 各种采样算法（如拒绝采样、MCMC）

  • 变分近似方法 [Bingham et al. 2019; Carpenter et al. 2017; Chaganty et al. 2013; Dillon et al. 2017; Goodman et al. 2008; Hur et al. 2015; Kucukelbir et al. 2015; Mansinghka et al. 2013, 2018; Minka et al. 2014; Pfeffer 2007; Saad and Mansinghka 2016; Tristan et al. 2014; van de Meent et al. 2015; Wingate and Weber 2013; Wood et al. 2014]

像拒绝采样和 MCMC 这样的采样算法虽然通用，但也存在已知的局限性，例如难以处理多模态分布和低概率证据（见第2节）。

更高级的技术如哈密顿蒙特卡洛（HMC）和变分近似方法可以缓解这些问题，但它们要求函数具有连续性和几乎处处可微性，因此必须通过边缘化消除所有离散结构。

使用二进制表示

“位爆破”是一种广泛应用于软件验证 领域的技术，常被约束求解器用来基于二进制表示进行算术推理 [Bruttomesso 和 Sharygina 2009]。

最近也有研究在整数上的概率程序推理中采用数字的二进制表示 [Cao et al. 2023]，以便利用该表示中的条件独立性。

HyBit 中的位爆破受到这些技术的启发，但其目标不同，因而采用了截然不同的技术路线：即开发简洁且在许多情况下可证明正确 的连续概率分布近似方法。

7 结论与未来工作

在本研究中，我们阐明了为混合概率程序 开发新型推理方法的必要性。

我们提出了位爆破 （bit blasting）方法：通过将混合概率程序进行简洁的离散化表示 ，然后使用离散推理算法 对其进行分析。

我们定义了一类连续密度函数——混合伽马密度 （mixed-gamma densities），并证明对于这类分布，位爆破不仅具有简洁性 ，而且相较于显式的离散化方法是正确的 （sound），并且在推理上是可证明高效 的。

在此基础上，我们提出了一种新的概率编程语言 HyBit ，它基于位爆破技术，采用了一种全新的针对混合程序的推理算法。我们通过实验展示了 HyBit 在性能上优于现有的近似推理算法。在未来的工作中，我们希望进一步扩展能够被正确位爆破的分布类别 。我们计划研究如何将 HyBit 扩展以支持分层贝叶斯模型 （hierarchical Bayesian models）。此外，我们还计划提升其易用性，不再要求用户为每个概率程序手动指定超参数 。我们也感兴趣于探索将 HyBit 与其他推理方法进行集成，以结合它们各自的优势，从而支持更广泛的混合概率程序。

原文链接：https://github.com/Tractables/Dice.jl/tree/hybit