干预，在何处及如何进行？大规模非线性 SCM 的主动因果发现

来自单个细胞的蛋白质信号网络的结构是什么？不同的习惯如何影响疾病的出现？这样的问题涉及由非线性的有噪声的过程控制的复杂系统中的因果效应。在大多数场合下，这样的系统的被动观测是不足以发现真实的因果效应关系的，需要开销大的实验来消除相互竞争的假设间的不明确性。因此，实验的设计是很令人感兴趣的，有效的实验方案有助于降低实验的成本，同时通过（闭环的 / 策略驱动的）科学方法帮助产出知识的过程（图 1）。

在因果的语言中，因果关系被一个有向无环图（DAG）定性地表示，图中节点对应于研究的系统中不同的变量，边表示变量间的信息流。DAG 的抽象使得我们可以表示眼前的观测的可能解释（假设）空间。以贝叶斯概率（信念）表示这样的假设允许我们将科学方法的问题形式化为一个贝叶斯推断，其中目标是要估计后验分布 $p(\text{DAGs}\,|\,\text{Observations})$。DAG 上的后验概率可以允许我们使用信息论的获得函数，该函数指导实验朝着最具信息量的变量进行，以消除互相竞争的假设间的不明确性。这样的设计过程属于用于因果发现的贝叶斯最优实验设计（BOECD）的领域。

在贝叶斯最优实验设计（BOED）的框架中，要寻找能够最大化关于某些感兴趣参数的期望信息增益的实验。在因果发现中，实验是因果干预的形式，感兴趣的参数是结构因果模型（SCM）及其相关联的 DAG。

因果模型中的干预指的是我们操作的变量（或目标）和我们要将该变量设定为的值（或强度）。因此，对于学习因果模型的设计空间就是干预目标的所有子集和所选目标的干预值的可能可数无穷集合的集合。干预值封装了在许多因果推断应用中重要的语义。例如，在药物应用中，干预可以对应于不同药物的给药，干预值以每种药物的剂量水平的形式存在。尽管该值的适当选择对于确定底层的因果模型至关重要，但现有的主动因果发现的工作只关注于选择干预目标。这样的工作中，干预值通常都是某个任意的固定值（比如 0），这是次优的（见图 2a）。因此，在非线性因果模型的一般情况下，缺乏对干预值和干预目标选择的整体处理。本文提出了一个贝叶斯实验设计方法（CBED，读音为“seabed”），通过贝叶斯优化来获得最优的干预目标和干预值。

此外，一些情境要求选择一批干预。批量干预的问题计算很昂贵，因为它需要评估所有可能的干预组合。我们将 CBED 扩展到批设定下，并提出了两种不同的批策略，用于获取可解的贝叶斯最优干预目标和值。第一个策略 Greedy-CBED 贪心地构建干预集。由于互信息的子模块性，贪心启发式算法仍然是接近最优的。第二种策略 Soft-CBED 通过从有限的候选集合中随机抽样构建一组干预，从而显著提高了计算效率，恢复 DAG 结构和 SCM 的参数，速度与贪心策略一样快。这种策略非常适合资源受限的情境。

在整个工作中，我们做了如下的因果发现的标准假设：

假设 1（因果充分性）没有隐藏的混淆因子，所有感兴趣的随机变量都是可观测的。

假设 2（有限样本）可用的观测 / 干预样本数量有限。

假设 3（有加性噪声的非线性 SCM）结构因果模型具有附加高斯噪声的非线性条件期望。

假设 4（单目标）每个干预都是原子的，并应用于 SCM 的单个目标。

此外，我们假设干预以大小为 $\mathcal{B}$ 的批来设计并执行的，全部干预的固定预算由 $\text{Number of Batches}\times\mathcal{B}$ 给出。我们也假设底层的图是稀疏的，因为所有现实世界设定中也是这样。实验设计在稀疏图设定中更合适，因为信息量大的干预目标和值的数量会比稠密图的少得多。对应于稀疏图的许多节点具有父集的概率非常小，因此使用随机策略进行实验不能最大程度地提供信息。最后，我们对于用少数量的批来恢复全图 $\mathbf{G}$ 很有兴趣。和所有的因果推断任务一样，我们上面做出的假设要被小心地对感兴趣的应用验证。

我们证明了我们的方法 Greedy-CBED 和 Soft-CBED 在线性和非线性 SCM 设定中表现得比最先进的主动因果发现基线更好。此外，我们的方法在受现实世界启发的非线性数据集 DREAM 中取得了优异的结果。

记号 令 $\mathbf{V}=\{1,\dots,d\}$ 为任一 DAG $\mathbf{g}=(\mathbf{V},E)$ 的顶点集且 $\mathbf{X_V}=\{\text{X}_1,\dots,\text{X}_d\}\subseteq\mathcal{X}$ 为由 $\mathbf{V}$ 索引的随机变量。我们有一个初始观测数据集 $\mathcal{D}=\{\mathbf{x_V}^{(i)}\}_{i=1}^n$，由实例 $\mathbf{x_V}\sim P(\text{X}_1=\text{x}_1,\dots,\text{X}_d=\text{x}_d)=p(\text{x}_1,\dots,\text{x}_d)$。

因果贝叶斯网络 因果贝叶斯网络（CBN）是 $(\mathbf{g},P)$ 对，使得对于任意 $\mathbf{W}\subset\mathbf{V}$，

其中 $\text{do}(X_\mathbf{W})$ 表示在变量 $X_\mathbf{W}$ 上的假设干预，$1(\cdot)$ 是指示函数，$\text{pa}_\mathbf{g}(i)$ 表示 DAG $\mathbf{g}$ 中变量 $X_i$ 的双亲。在任意变量 $X_j$ 上的完美干预完全移除了对于其双亲的所有依赖，即 $P(X_j|X_{\text{pa}_\mathbf{g}(j)})=P(X_j)$ 从而得到一个残缺 DAG $\mathbf{g}'=(\mathbf{V},E\setminus(\text{pa}_\mathbf{g}(j),j))$。

结构因果模型 从数据生成机制的视角，$\mathbf{X_V}$ 上的 DAG $\mathbf{g}$ 与一组结构等式相匹配：

其中 $f_i$ 是（可能非线性）因果机制，在干预任何变量 $X_j\ne X_i$ 时保持不变，$\epsilon_i$ 是分布任意的外源噪声变量，其相互独立，即 $\epsilon_i$ 独立于 $\epsilon_j$ $\forall i\ne j$。(1) 表示在一个因果贝叶斯网络中的条件分布，如果机制已知，可以进一步揭示干预的效果。这些等式共同形成了结构因果模型（SCM），与相关联的 DAG $\mathbf{g}$ 一起。尽管机制 $f$ 一般情况下可能是非参数化的，我们假设存在一个参数为 $\gamma\in\Gamma$ 的这些机制的参数化近似。在线性 SCM 的情况中，$\gamma$ 对应于 $E$ 中边的权重。在非线性的情况中，它们可以表示非线性函数的参数，其参数化了一个高斯分布的均值。

(1) 的一个常见形式对应于高斯加性噪声模型（ANM）：

一个 ANM 由一个 DAG $\mathbf{g}$ ，参数为 $\gamma=[\gamma_1,\dots,\gamma_d]$ 与方差 $\sigma^2=[\sigma^2_1,\dots,\sigma^2_d]$ 的机制 $f(\cdot;\gamma)=[f_1(\cdot;\gamma_1),\dots,f_d(\cdot;\gamma_d)]$ 完全指定。为了记号简洁，此后我们记 $\theta=(\gamma,\sigma^2)$ 且所有感兴趣的参数为 $\phi=(\mathbf{g},\theta)$。

贝叶斯因果发现 因果推断中的一个常见假设是因果关系是定性已知的且可以被一个 DAG 表示。虽然在某些情况下可以从领域知识中获得这种定性信息，但在大多数应用中是不可行的。因果发现的目标是要给定一个数据集 $\mathcal{D}$，恢复 SCM 以及相关联的 DAG。通常来说，没有对机制 $f$ 的类型进一步的假设（如线性或非线性），真的 SCM 可能不是只从观测数据中可识别的。这种不可识别性是因为可能有多个 DAG（从而有多个 $P(\mathbf{X_V})$ 的因式分解）都同样很好地解释了数据。这样的 DAG 被称为马尔可夫等价。干预可以提高可识别性。除了可识别性问题，使用有限数据估计节点间的函数关系是另一个不确定性的来源。未知 SCM 上的贝叶斯参数估计提供了一个条理化的方式来量化这些不确定性并得到给定观测数据下 SCM 上的后验分布。然后实验者可以使用后验编码的知识来设计信息量大的实验，其有效地取得干预数据来解决未知边朝向和函数不确定性。

SCM 和 DAG 的贝叶斯推断 进行 SCM 和 DAG 上联合贝叶斯推断的关键挑战是 DAG 的空间是离散的且是变量数量的超指数。然而，近期的基于变分推断的技术提供了一个可解的、可扩展的进行这些参数的后验推断的方式。给定一个近似了后验 $p(\phi|\mathcal{D})$ 的可解的分布 $q_\psi(\phi)$，变分推断最大化（对数）证据下界：

这些技术中的关键思想是用于 DAG 的变分族 $\Psi$ 被参数化的方式。变分因果网络（VCN）方法的变分族是一个邻接矩阵上的自回归伯努利分布。他们通过先验进一步强制了无环性限制。BCD-Net 通过玻尔兹曼分布考虑节点排列上的分布并用 Gumbel-Sinkhorn 算子进行推断。DiBS 考虑邻接矩阵项上的潜变量并在这些潜变量上用 SVGD 进行推断。我们证明了这个工作中使用 DiBS 模型的 BOECD 的实验结果，因为其很容易扩展到非线性 SCM 中。

贝叶斯最优实验设计 贝叶斯最优实验设计是一个信息论的对于选择最优的实验来估计任意参数 $\theta$ 的问题的方法。对于 BOED，实验 $\xi$ 的效用是观测 $\mathbf{y}$ 和 $\theta$ 之间的互信息（MI）：

BOED 的目标是要选择能最大化目标 $\xi^*=\argmax_\xi\mathrm{U_{BOED}}(\xi)$ 的实验。不幸的是，由于嵌套期望，求值与优化这个目标很有挑战，几个估计量已被引入，是 BOED 目标的下界，可以接着被并入许多优化方法来选择设计。

一个常见的设定，叫做静态、固定或批设计，是要同时优化 $\mathcal{B}$ 个设计 $\{\xi_1,\dots,\xi_\mathcal{B}\}$。这些设计接着被执行，实验结果被收集用于以一种贝叶斯风格更新模型参数。

随机变量 $\mathbf{X_V}$ 上的真的 SCM 和相关联的 DAG $\tilde{\phi}=(\tilde{\mathbf{g}},\tilde{\theta})$ 是事实，但我们对 $\tilde{\phi}$ 的信念由于多个原因是不确定的。首先，从观测数据 $\mathcal{D}$ 中最多只可能将 DAG $\tilde{g}$ 学到马尔可夫等价类（MEC）。不确定性也产生于 $\mathcal{D}$ 为有限样本，我们通过引入随机变量 $\Phi$ 来建模，其中 $\phi$ 是结果。令 $\phi\sim p(\phi|\mathcal{D})\propto p(\mathcal{D}|\phi)p(\phi)$ 为随机变量 $\Phi$ 的一个实例，其是从我们在观测了数据集 $\mathcal{D}$ 之后的 SCM 上的后验中抽样得到的。

我们想要设计一个实验来表明一个干预 $\xi\coloneqq\{(j,v)\}\coloneqq\mathrm{do}(\mathrm{X}_j=v)$，其最大化在观测到干预结果 $\mathbf{y}\sim P(\mathrm{X}_1=\mathrm{x}_1,\dots,\mathrm{X}_d=\mathrm{x}_d|\mathrm{do}(\mathrm{X}_j=v))=p(\mathbf{y}|\xi)$ 之后关于 $\Phi$ 的信息增益。这里，$\mathbf{y}$ 是随机变量 $\mathbf{Y}\sube\mathcal{X}$ 的一个实例，其是根据由在干预下的残缺真图 $\tilde{g}'$ 指定的分布而分布的。一次只看一个干预，可以将 BOECD 形式化为在观测到一次实验的结果 $\mathbf{y}$ 之后关于 $\Phi$ 的信息的增益。最大化信息增益的实验 $\xi\coloneqq\{(j,v)\}$ 是最大化 $\Phi$ 和 $\mathbf{Y}$ 之间的互信息的实验：

上面的目标不仅在干预目标 $j\in\mathbf{V}$ 的离散集合上，而且在干预值 $v\sub\mathcal{X}_j$ 的不可数集合上考虑取 $\argmax$。虽然现存的 BOECD 的工作仅仅考虑干预目标的设计来限制复杂度，我们的方法同时解决了两个问题。我们首先概述单一设计的方法，然后在 3.2 节论述如何将单一设计扩展到批设定中。

为了最大化公式 (3) 中的目标，我们需要（1）估计候选干预的 MI 以及（2）通过对每个候选干预目标在干预值的域上进行优化来最大化估计的 MI。

估计 MI 因为互信息是不可解的，有许多方式来估计互信息，取决于我们是否能够从后验中采样以及是否似然可以求值。由于我们考虑的模型允许后验采样和似然求值，其足够得到一个仅需要似然求值和期望的蒙特卡洛近似的估计量。为了这样做，我们得到一个类似于贝叶斯分歧主动学习（Bayesian Active Learning by Disagreement，BALD）的估计量，其将 MI 考虑为结果 $\mathbf{Y}$ 上的条件熵的差：

其中 $\mathrm{H}(\cdot)$ 是熵。推导见附录 B.1。上述公式的蒙特卡洛估计量可以被用作近似（附录 B.2）。公式 (4) 有一个直观的解释。它把高的互信息分配给模型关于结果最不一致的干预。我们记单一设计的 MI 为 $\mathcal{I}(\{(j,v)\})\coloneqq\mathrm{I}(\mathbf{Y};\Phi|\{j,v\},\mathcal{D})$。

选择干预值 正如 (3) 中所示，取得最大化目标不仅要选择干预目标而且要设定合适的干预值。尽管优化干预目标是可解的（从离散且有限数量个节点中选择），选择干预值通常是不可解的，因为它们是连续的。对于任意给定的目标节点 $j$，MI 是一个 $v\in\mathcal{X}_j$ 上的非线性函数（见图 2），因此用梯度上升方法求解只能得到一个局部最大值。考虑到 MI 的求值是很昂贵的，我们将一个给定的目标节点 $j$ 的 MI 视为一个黑盒函数，使用贝叶斯优化（BO）来获得它的最大值。BO 试图在整个集合 $\mathcal{X}_j$ 上用尽可能少的求值次数找到这个函数的最大值 $\max_{v\in\mathcal{X}_j}\mathcal{I}(\{(j,v)\})$。细节见附录 E。

BO 通常通过安排一个函数 $\mathcal{I}(\{j,\cdot\})$ 上的高斯过程（GP）先验，然后用查询点 $\mathbf{v}^*=\{v^{(1)*},\dots,v^{(T)*}\}$ 得到这个函数的后验来进行。令在每个优化步骤 $t$ 查询的互信息值为 $\hat{\text{U}}_j^t=\mathcal{I}(\{j,v^{(t)*}\})$。一点 $v^{(t+1)}$ 的后验预测可以以闭式形式为一个均值为 $\boldsymbol{\mu_j^{(t+1)}}(v)$ 且方差为 $\boldsymbol{\sigma_j^{(t+1)}}(v)$ 的高斯分布得到。下标 $j$ 表示我们对每个干预目标维护不同的 $\mu$ 和 $\sigma$，上标 $t$ 表示 BO 步数。查询通过有一个定义在这个后验上的采集函数，其建议下一个查询的点来进行。对于 BO，我们使用一个叫做置信上界（Upper Confidence Bound，UCB）的采集函数，其通过用一个超参数 $\beta_j^t:v_j^{(t+1)*}=\argmax_v\boldsymbol{\mu_j^t}(v)+\sqrt{\beta_j^{t+1}}\boldsymbol{\sigma_j^t}(v)$ 权衡探索和利用来建议下一个查询的点。

我们独立地在每个候选干预目标 $j=\{1,\dots,d\}$ 上通过查询在一个固定域 $[-k,k]\sub\mathbb{R}$ 内的点来运行 GP-UCB。注意这个域可以基于应用来选择，例如，如果我们必须限制剂量等级在一个固定的范围内。每一个 GP 在我们的设定中是一维的，因为一些 UCB 的求值足够得到一个很好的最大候选值了。此外，每个候选目标的 GP-UCB 是并行的，使其更有效率。我们最后选择在所有候选干预目标中有着最大 MI 值的设计。

在许多应用中，需要选择最有信息量的干预集合而非一次一个干预。举例来说，一个生物学家带着一份实验脚本进入一个潮湿的实验室。批量实验消除了直到执行下一个实验等待实验完成和分析的瓶颈。给定每个批 $\mathcal{B}$ 的预算，其表示在一个批中实验的次数，选择批的问题就变成了 $\argmax_\Xi\text{I}(\mathbf{Y};\Phi|\Xi,\mathcal{D})$，使得 $\text{cardinality}(\Xi)=\mathcal{B}$，其中 $\Xi$ 是干预集合 $\bigcup_{i=1}^\mathcal{B}(j_i,v_i)$，$\mathbf{Y}$ 表示批的干预的结果的随机变量。我们记批设计的 MI 为 $\mathcal{I}(\Xi)\coloneqq\text{I}(\mathbf{Y};\Phi|\Xi,\mathcal{D})$。

贪心算法 计算最优解 $\mathcal{I}(\Xi^*)$ 是计算不可行的。然而，由于条件互信息是子模的且非递减的（证明见附录 B.4），我们可以得到一个简单的贪心算法（算法 1），其可以得到最优解的至少 $(1-1/e)\approx0.64$ 的近似。我们记这个策略为 $\text{Greedy-CBED}$。

Soft Top-K 尽管贪心算法是可解的，其需要 $O(\mathcal{B}d)$ GP-UCB 的实例。Kirsch 等人 (2021) 证明了软 top-k 选择策略表现和贪心算法相似，将计算需要降低到 $O(d)$ 次运行 GP-UCB。为了实现这一方法，我们通过为每个节点 $j=\{1,\dots,d\}$ 保留所有 $T$ 个 GP-UCB 的求值来构建一个候选干预目标-值对的有限集合。因此，对于 $d$ 个节点，我们的候选集是由 $d\times T$ 个实验组成的。我们用 MI 估计来为候选集中的每个实验进行打分。接着我们按照 MI 分数的 softmax 的比例进行 $\mathcal{B}$ 次的不放回抽样（算法 2）。我们记这个策略为 $\text{Soft-CBED}$。

算法 1：Greedy-CBED

输入：$\mathcal{E}$ 环境，$N$ 初始观测样本，$\mathcal{B}$ 批大小，$d$ 节点数量

$\rhd$ 初始化实验集合 $\Xi$ 为空

1 $\Xi\leftarrow\emptyset$

2 for $n=1\dots\mathcal{B}$ do

3 for $j=1\dots d$ do

$\rhd$ 使用 GP-UCB 选择每个节点 $j$ 的最优干预值

4 $V_j\leftarrow\argmax_v\mathcal{I}(\Xi\cup\{(j,v)\})$

5 $U_j\leftarrow\mathcal{I}(\Xi\cup\{(j,V_j)\})$

6 $j^*\leftarrow\argmax_jU_j$

7 $v^*\leftarrow V_{j^*}$

8 $\Xi\leftarrow\Xi\cup\{(j^*,v^*)\}$

9 返回 $\Xi$

算法 2：Soft-CBED

输入：$\mathcal{E}$ 环境，$N$ 初始观测样本，$\mathcal{B}$ 批大小，$d$ 节点数量，$\zeta$ softmax 温度

1 for $j=1\dots d$ do

$\rhd$ 使用 GP-UCB 选择每个节点 $j$ 的最优干预值

2 初始化 $\mu_j^0$ 和 $\sigma_j^0$

3 for $t=1\dots T$ do

4 $V_j^t\leftarrow\argmax_v\mu_j^{t-1}(v)+\sqrt{\beta^t}\sigma_j^{t-1}(v)$

5 $\hat{U}_j^t\leftarrow\mathcal{I}(\Xi\cup\{(j,V_j^t)\})$

6 更新 GP 来得到 $\mu_j^t$ 和 $\sigma_j^t$

7 $\{(t_i,j_i)\}_{i\in\{1,\dots,\mathcal{B}\}}\leftarrow\mathcal{B}$ 个样本不放回抽样 $\propto\exp(\hat{U}_j^t/\zeta)$

8 $\Xi\leftarrow\{(j_i,V_{j_i}^{t_i})\}_{i\in\{1,\dots,\mathcal{B}\}}$

9 返回 $\Xi$

我们概述了我们的方法和两个主要的现有的主动因果发现方法的对比。

ABCD（Agrawal 等人，2019）ABCD 中使用的 MI 估计量是基于加权重要性采样的。然而，对于 ABCD 中使用的重要性采样权重的具体选择，它们的 MI 估计结果与我们的方法近似相同（见附录 B.5）。然而，ABCD 没有选择干预值，而是次优地将它们设定为一个固定值。此外，我们提出的 $\text{Soft-CBED}$ 是一个更快且更高效的批策略，尤其是当值也需要获取时。从这一点上看，我们的方法是 ABCD 的非线性假设、值获取以及软 top-k 批量策略的扩展。

AIT（Scherrer 等人，2021）AIT 是一个基于 F 分数的干预目标获取策略。尽管它不是一个 BOECD 方法，我们在这里证明了其可以被视为一个 MI 的蒙特卡洛估计量近似，当结果 $\mathbf{Y}$ 是高斯的时。然而，AIT 和 ABCD 一样没有选择干预值并且没有批量策略。

定义 3.1 令 $\mathbf{Y}$ 为一个高斯随机变量。那么 Scherrer 等人 (2021) 的差异分数是互信息（公式 (4)）的一个近似的一个蒙特卡洛估计量。证明见附录 B.6。

早期的使用贝叶斯最优实验设计用于因果发现（BOECD）的尝试可以在 Murphy，2001 和 Tong 和 Kolloer，2001 的工作中找到。然而，这些方法解决了像限制图为拓扑有序结构、序列化干预、线性模型和离散变量等简单设定下。

在 Cho 等人，2016 和 Ness 等人，2017 中，BOECD 被应用于学习生物网络结构。BOECD 也在 Greenwald 等人，2019 中被探索，在图的无向边总能形成一颗树的假设之下。更近期，ABCD 框架（Agrawal 等人，2019）扩展了 Murphy，2001 和 Tong 和 Kolloer，2001 的工作，设定是干预可以应用于连续变量的批量中。为了实现这一目标，他们（近似）解决了最大化干预（实验）、结果和观察数据之间的批量互信息的子模性问题，给定 DAG。DAG 假设是使用 DAG-bootstrap（Friedman 等人，2013）采样得到的。我们的工作在几个方面与 ABCD 不同：我们通过在 DAG 上使用最先进的后验模型（Lorch 等人，2021）来处理线性和非线性 SCM，我们使用 BO 来选择要干预的值，但我们也使用 softBALD（Kirsch 等人，2021）来准备批量，其明显快于 ABCD 方法的贪心逼近。

在 von Kugelgen 等人，2019 中，作者提出了使用高斯过程来建模 DAG 上的后验分布，并使用 BO 来识别要干预的值，然而该方法并未展示出对于比二变量图更大的图的可扩展性，由于其依赖于多维高斯过程建模条件分布。

在可微分因果发现领域出现了一个新的工作主体，在这个领域中，通常从观测数据中找到结构的问题，是用梯度上升和函数近似（如神经网络）来解决的（Zheng 等人，2018，Ke 等人，2019，Brouillard 等人，2020，Bengio 等人，2019）。在最近的工作中（Cundy 等人，2021，Lorch 等人，2021，Annadani 等人，2021），作者提出了一个 DAG 上后验的变分近似，其考虑建模一个分布而非 DAG 的一个能最好地解释观测数据 $\mathcal{D}$ 的点估计。这样的工作可以用来替代 DAG-bootstrap（Friedman 等人，2013），允许具有更大支持的后验分布的建模。

除了基于 BOECD 的方法之外，一些主动因果学习方法也被提出（He 和 Geng，2008，Gamella 和 Heinze-Deml，2020，Scherrer 等人，2021，Shanmugam 等人，2015，Squires 等人，2020，Kocaoglu 等人，2017）。主动 ICP（Gamella 和 Heinze-Deml，2020）使用 ICP（Peters 等人，2016）用于因果学习当选择一个主动的选择目标的策略时，然而，这个工作没适用于全图都需要被恢复的场景。在 Zhang 等人，2021 中，作者提出了一个主动学习方法来解决识别将动态因果网络推向理想状态的干预的问题。一些方法解决了主动获取干预数据以定向骨架图中边的问题（Shanmugam 等人，2015，Squires 等人，2020，Kocaoglu 等人，2017）。更接近我们提出的方法的是 AIT（Scherrer 等人，2021），其使用图上的一个基于神经网络的后验模型，但评估 F 值来选择干预。

我们在人造的和真实世界的因果实验设计问题和一系列基线上评估了我们的方法的性能。我们目的是在经验上研究下面几个方面：（1）Greedy-CBED 和 Soft-CBED 整体提出的策略在大规模（50 个节点）人造数据集上的竞争力；（2）基于 GP-UCB 的值采集策略的性能；（3）提出的方法在真实世界启发的数据集上的性能。

随机 随机基线随机地获取干预目标。

AIT / softAIT 主动干预目标（AIT）（Scherrer 等人，2021）使用一个基于 f 分数的采集策略来选择干预目标。更多细节见附录 B.6。因为原来提出的方法没有考虑批设定，我们引入了一个将 AIT 用如同 3.2 节中描述的 soft 批量化增强后的变种。

CBED / Greedy-CBED / Soft-CBED 这些是如同在第 3 节中描述的 MI 的蒙特卡洛估计量。CBED 选择单一的最大化 MI 的干预（目标与值）且这个干预被用在整个批中。Greedy-CBED 中，批是以贪心的方式来建立的，一次选择一个目标-值对（算法 1）。Soft-CBED 按 MI 分数的比例采样（目标，值）对来选择一个批，如同 3.2 节以及算法 2 中描述的。

固定：这个值选择策略假定将干预值设定为一个固定值。实验中，我们将值固定为 0。

Sample-Dist：这个值选择策略从观测数据的支持中采样。

GP-UCB：这个策略使用提出的 GP-UCB 贝叶斯优化策略来选择最大化 MI 的值。

人造图 在这个设定中，我们生成了大小为 20 的 Erdos-Rényi（ER）和大小为 50 的 Scale-Free（SF）图。对于线性 SCM，我们均匀随机采样边的权重 $\gamma$。对于非线性的 SCM，我们将每个变量参数化为均值为其父母的一个非线性函数的高斯分布。我们用神经网络建模非线性函数。在所有设定中，我们将噪声方差设置为 $\sigma^2=0.1$。对于两种图，我们将每个节点的期望边的数量设置为 1。我们在附录 D.1 中提供了实验的更多细节。

单细胞蛋白质信号网络 DREAM 基准族被设计用于评估一个细胞的控制网络的因果发现算法。一组 ODE 和 SDE 生成了数据集，其模拟单个细胞的逆向工程的网络。我们使用 GeneNetWeaver 来模拟稳态 wind-type 表示和单基因 knockouts。参考附录 D.2 精确的设定。

$\mathbb{E}$-SHD：定义为图上的后验模型中采样的样本和真图 $\mathbb{E}$-SHD $\coloneqq\mathbb{E}_{\mathbf{g}\sim p(\mathcal{G}|\mathcal{D})}[\text{SHD}(\mathbf{g},\tilde{\mathbf{g}})]$ 之间的期望结构汉明距离。

$\mathbb{E}$-SID：由于 SHD 对于干预的概念是不可知的，期望结构干预距离（$\mathbb{E}$-SID）量化了关于因果推断表述和干预分布的图之间的差异。

AUROC：预测所有边是否存在 / 不存在的二元分类任务的接收者操作特征曲线下面积。

AUPRC：预测所有边是否存在 / 不存在的二元分类任务的准确率-召回率曲线下面积。

对于每个采集目标和数据集，我们给出了期望结构汉明距离 $\mathbb{E}$-SHD 和 期望结构干预距离 $\mathbb{E}$-SID（Peters and Buhlmann，2015）的平均值和标准差，接收者操作特征曲线下面积 AUROC 和准确率-召回率曲线下面积 AUPRC。我们将这些指标作为获得的干预样本数量的函数进行评估（或实验），这有助于定量地比较不同的获取策略。除 $\mathbb{E}$-SID 外，我们将其他指标的结果归入附录 I 。

在合成图（图 3(a)）中，我们可以看到，对于具有 $50D$ 变量和非线性函数关系的 ER 和 SF 图，基于软 top-k 选择使用 GP-UCB 的批量方法在 $\mathbb{E}$-SID 度量方面优于所有基线。另一方面，即使在与所提出的值获取相结合后，单独的 AIT 也不能快速收敛到真实图，但当与所提出的软策略进一步增强时，softAIT 恢复真实因果图的速度可提高 4 倍，并且与 $\text{Soft-CBED}$ 相比具有竞争力。我们在其他指标上也观察到类似的表现。此外，我们发现这种趋势适用于 $20D$ 变量和线性模型。全部结果载于附录 I。

接下来，我们将研究具有主动因果发现的值选择策略的重要性。我们在公式 4 中使用 MI 估计量；此外，我们使用两种启发式方法——固定值策略和从支持中采样值来测试所提出的 GP-UCB。正如我们在图 3(b) 中所看到的，使用 GP-UCB 选择值显然有利于因果发现过程。我们期望这一发现，因为互信息相对于干预值不是恒定的。为了使这一点更清楚，我们在附录 G 中展示了一个简单的双变量图中值的影响。此外，我们注意到，从观察数据集的支持中简单地采样比将值固定为 0 表现更差。我们假设这是由于支持高密度区域的认知不确定性较低，暗示这些区域可能信息较少。

为了进一步了解软批策略与其他批选择策略的比较，我们将 $\text{Soft-CBED}$ 与 $\text{Greedy-CBED}$ 和 $\text{CBED}$ 的结果进行了比较。我们观察到（图 3(c)），Greedy-CBED 和 $\text{Soft-CBED}$ 总体上给出了非常相似的结果。虽然 $\text{Greedy-CBED}$ 在某些条件下是最优的（Kirsch 等人，2019），但 $\text{Soft-CBED}$ 仍然具有竞争力，并且具有可以一次性选择批的优势。从表 1 中这两种批处理策略的运行时性能也可以看出这一点。这两种批选择策略都比选择一个干预目标/值对并执行 $\mathcal{B}$ 次（$\text{CBED}$）要好得多。

最后，在 DREAM 任务中，我们看到我们的方法在 $\mathbb{E}$-SID 指标上优于 softAIT 和随机基线（见图 4）。在这些实验中，由于干预是模拟基因敲除设定，我们只使用固定值策略，其值为 0.0。尽管随机基线仍然是一种有竞争力的选择，但在某些情况下，$\text{Soft-CBED}$ 目标明显更好（Ecoli1，Ecoli2 数据集）。

本文研究了有效选择贝叶斯最优实验来发现因果模型的问题。我们提出的框架同时回答了在批量设定中在何处以及如何进行干预的问题。我们提出了一种贝叶斯优化策略来获取干预目标和干预值。此外，我们提出了两种不同的批量策略：一种基于贪心选择，另一种基于软 top-k 选择。所提出的在批量设定中选择干预目标-值对的方法为多达 $50D$ 变量的因果模型提供了优于最先进的的性能。我们使用合成数据集和单细胞调节网络的受现实世界启发的数据集验证了这一点，显示了对生物学和其他实验科学等领域的潜在影响。

B.1 推导结果上的互信息

在下面的引理中，我们推导在 (4) 中给出的结果上的互信息。

引理 B.1

证明

B.2 估计结果上的互信息

对于允许对实验结果密度（似然）$p(\mathbf{y}|\phi,\{(j,v)\})$ 求值的模型，我们可以使用下面的估计量来估计 $\mathrm{I}(\mathbf{Y};\Phi|\{(j,v)\},\mathcal{D})$：

定义 B.1 实验结果的边缘熵 $\text{H}(\mathbf{Y}|\{(j,v)\},\mathcal{D})$ 的蒙特卡洛估计量 $\widehat{\text{H}}(\mathbf{Y}|\{(j,v)\},\mathcal{D})$ 由下式给出：

其中 $\widehat{\mathbf{y}}_{i,k}\sim p(\mathbf{y}|\widehat{\phi}_i,\{(j,v)\})$ 是从由被干预 $\{(j,v)\}$ 增强后的 $c_o$ 个 SCM 中第 $i$ 个 $\widehat{\phi}_i\sim p(\phi|\mathcal{D})$ 参数化的密度中抽样的 $m$ 个样本之一。样本 $\widehat{\mathbf{y}}_{i,k}$ 的似然接着在被干预 $\{(j,v)\}$ 增强后的额外的 $c_{in}$ 个 SCM 中的第 $l$ 个 $\widehat{\phi}_l\sim p(\phi|\mathcal{D})$ 的参数化下求值。

$\widehat{\text{H}}(\mathbf{Y}|\{(j,v)\},\mathcal{D})$ 是一个一致但有偏的 $\text{H}(\mathbf{Y}|\{(j,v)\},\mathcal{D})$ 的估计量，由于非线性 $\log$ 函数内部的期望。或者，我们可以观察下面的 $\text{H}(\mathbf{Y}|\{(j,v)\},\mathcal{D})$ 的下界：

由琴生不等式。那么我们就可以定义一个这个下界的无偏估计量。

定义 B.2 实验结果的边缘熵的下界 $-\mathbb{E}_{p(\mathbf{y}|\{(j,v)\},\mathcal{D})}\left[\mathbb{E}_{p(\phi|\mathcal{D})}[\log p(\mathbf{y}|\phi,\{(j,v)\})]\right]$ 的无偏蒙特卡洛估计量 $\widehat{\text{H}}^*(\mathbf{Y}|\{(j,v)\},\mathcal{D})$ 由下式给出：

最终，我们定义 $\text{H}(\mathbf{Y}|\Phi,\{(j,v)\},\mathcal{D})$ 的估计量。

定义 B.3 实验结果在给定 $\Phi$ 条件下的熵 $\text{H}(\mathbf{Y}|\Phi,\{(j,v)\},\mathcal{D})$ 的蒙特卡洛估计量 $\widehat{\text{H}}(\mathbf{Y}|\Phi,\{(j,v)\},\mathcal{D})$ 由下式给出：

其中 $\widehat{\mathbf{y}}_{i,k}\sim p(\mathbf{y}|\widehat{\phi}_i,\{(j,v)\})$ 是从由被干预 $\{(j,v)\}$ 增强后的 $c_o$ 个 SCM 中第 $i$ 个 $\widehat{\phi}_i\sim p(\phi|\mathcal{D})$ 参数化的密度中抽样的 $m$ 个样本之一。

算法 3：互信息计算

输入：后验 $q(\phi|\mathcal{D}_{o}\cup\mathcal{D}_{in})$，后验采样数量 $c$，干预采样数量 $m$，干预 $\{(j,v)\}$。我们记 $Y$ 为观测数据。

$\rhd$ 从后验中采样

$\{\widehat{\phi}_i\sim q(\phi|\mathcal{D}_{o}\cup\mathcal{D}_{in})\}_{i=1}^c$

$\rhd$ 从残缺的 SCM 中采样

$\{\widehat{\mathbf{y}}_{i,j,k}\sim p(\mathbf{y}|\widehat{\phi}_i,\{(j,v)\})\}_{k=1}^m$

返回 $-\frac{1}{c\times m}\sum_{i=1}^c\sum_{k=1}^m\log(\frac{1}c{\sum_{l=1}^c}p(\widehat{\mathbf{y}}_{i,k}|\widehat{\phi}_l,\{(j,v)\}))+\frac{1}{c\times m}\sum_{i=1}^c\sum_{k=1}^m\log p(\widehat{\mathbf{y}}_{i,k}|\widehat{\phi}_i,\{(j,v)\})$

B.3 批互信息的蒙特卡洛估计量

尽管公式 (4) 适用于单一设计的 MI，我们在这里该处批设计的 MI 估计量。

B.4 互信息的子模性和单调性证明

贝叶斯优化（BO）是一种全局优化技术用于优化黑盒函数。更正式地说，对于任意定义在一个集合 $\mathcal{X}$ 上的求值很昂贵的函数 $\text{U}$ ，BO 试图在整个集合 $\mathcal{X}$ 上用尽可能少的求值次数来寻找函数的最大值。

BO 通常通过安排一个未知函数上的先验，然后用查询点 $\mathbf{x}^*=\{x^{*}_1,\dots,x^{*}_t\}$ 得到这个函数的后验来进行。一个常见的先验是高斯过程（GP），其均值为 0，协方差函数由一个核 $k(x,x')$ 定义。令 $\mathbf{U}_{\mathbf{x}^*}=[\text{U}(x_1^*),\dots,\text{U}(x_t^*)]$ 表示函数求值的向量，$\mathbf{K}$ 表示核矩阵，其中 $\mathbf{K}_{i,j}=k(x_i^*,x_j^*)$，令 $\mathbf{k}_{t+1}=[k(x_1^*,x_{t+1}),\dots,k(x_t^*,x_{t+1})]$。点 $x_{t+1}$ 的后验预测可以以闭式形式得到：

对于高斯过程，我们使用了下面的超参数。Matern 核，$\text{length scale}$ 1.0，$\text{length scale bounds}$（下=$1e-5$，上=$1e5$），$\text{Nu}$（学到的函数的平滑度）2.5。同时对核矩阵的对角线增加 $1e-6$。