# Investigating the Potential of Using Large Language Models for Scheduling

Deddy Jobson<EMAIL_ADDRESS>MercariTokyoJapan and  Li Yilin
<EMAIL_ADDRESS>MercariTokyoJapan

(2024; 2024-04-24; 2024-05-10; 20 February 2007; 12 March 2009; 5 June 2009)

###### Abstract.

The inaugural ACM International Conference on AI-powered Software introduced
the AIware Challenge, prompting researchers to explore AI-driven tools for
optimizing conference programs through constrained optimization. We
investigate the use of Large Language Models (LLMs) for program scheduling,
focusing on zero-shot learning and integer programming to measure paper
similarity. Our study reveals that LLMs, even under zero-shot settings, create
reasonably good first drafts of conference schedules. When clustering papers,
using only titles as LLM inputs produces results closer to human
categorization than using titles and abstracts with TFIDF. The code has been
made publicly available111https://github.com/deddyjobson/llms-for-scheduling.

large language models, scheduling, mathematical optimization, clustering

††copyright: acmlicensed††doi: 10.1145/3664646.3665084††journalyear:
2024††submissionid: fsews24aiwaremain-p27-p††isbn:
979-8-4007-0685-1/24/07††conference: Proceedings of the 1st ACM International
Conference on AI-Powered Software; July 15–16, 2024; Porto de Galinhas,
Brazil††booktitle: Proceedings of the 1st ACM International Conference on AI-
Powered Software (AIware ’24), July 15–16, 2024, Porto de Galinhas,
Brazil††ccs: Applied computing Multi-criterion optimization and decision-
making††ccs: Information systems Data mining

##  1\. Introduction

The inaugural ACM International Conference on AI-powered Software debuts the
AIware Challenge, encouraging researchers and practitioners to implement AI-
driven tools to create an optimized conference program, a fundamental issue in
constraint optimization. Constrained optimization is a prevalent mathematical
task, needing intensive computational resources and tailored algorithms to
address the unique complexities of the decision space. Applying Large Language
Models (LLMs) has demonstrated significant success in various domains such as
mathematical problem-solving and reasoning.

This report investigates the feasibility of leveraging LLMs to address the
program scheduling challenge. Our investigation revolves around two primary
approaches: 1) zero-shot learning to generate the schedules and 2) integer
programming with LLMs to measure the similarity between papers.

##  2\. Problem Statement

Our paper focuses on automating the allocation of papers to predetermined
sessions, treating it as an allocation and constrained clustering problem. We
do not factor in the specific timing of paper presentations or the impact of
parallel tracks, earmarking these as potential areas for future investigation
on the broader goal of automating conference scheduling.

##  3\. Experiments

###  3.1. Zero-Shot Scheduling by LLMs

Our first AIware approach directly prompts a large language model to make the
schedule. While LLMs have a lot of room to improve on arithmetic reasoning
tasks, recent approaches have been proposed to improve their abilities for
math problems(Yang et al., 2023; Zhang et al., 2024; Levonian et al., 2023).
Furthermore, LLMs are expected to improve at solving LLMs at zero-shot
configurations, too, due to the scaling laws(Wei et al., 2022).

However, from our experiments, we find the LLMs are still unable to create
schedules without violating constraints. We, therefore, additionally
experiment with smaller-scale problems by downsampling the number of sessions
or papers within each session. To measure the performance of our approach, we
measure the following metrics:

  * •

Completeness and Homogeneity wrt the ground truth(Rosenberg and Hirschberg,
2007)

  * •

% of constraints broken, for various types of constraints.

We use a temperature setting of 0.8 in our experiments.

We find that GPT-4 in zero-shot settings can generate reasonably good
conference schedules. We also counted the number of violations according to
the schedule created by the GPT -4. Looking at violations on the number of
sessions and papers in the proposed schedule, we find that, on average, less
than 3 papers were missing, and two sessions were added. Furthermore, checking
if the session durations are within bounds, we find that while many sessions
proposed were too long (in some cases, more than 50%), in most cases, the
proposed session was off by less than 10% of the length of the session. This
suggests that by moving around a couple of papers between sessions, a human
can finish the work mostly done by an LLM, suggesting a collaboration with
LLMs be the best approach.

Computing the homogeneity and completeness scores with respect to the actual
schedule of MSR 2022 (Figure 1), we find that while homogeneity drops with an
increase in problem size, completeness is not as much affected.

Figure 1. Comparing the homogeneity and completeness score for different
number of papers per session. Homogeneity deteriorates when the scale of the
data is larger, but not completeness.

###  3.2. Integer Program using LLM to measure similarity

In Section 3.1, we discover that LLMs face challenges in strict adherence to
all constraints, particularly when dealing with many papers. We plan to
explore how to leverage the text-understanding capabilities of LLMs to
facilitate the problem.

We incorporate the similarity of papers as an objective in our optimization
process to build an integer programming problem and use LLM to generate the
similarity. Recognizing that similarity is a subjective concept, we have
utilized the existing session arrangement as the benchmark for comparison. To
assess the proximity of our results to the original schedule, we employ the
completeness score and homogeneity score metrics (Rosenberg and Hirschberg,
2007).

Leveraging LLMs’ text-understanding capabilities, we cluster the papers into 5
groups and contrast the outcomes with a Bag of Words approach featuring TFIDF
normalization. The comparison, outlined in Table 1, demonstrates the superior
performance of LLM when utilizing only the paper titles compared to TFIDF with
both titles and abstracts.

To further evaluate the clustering outcomes, we formulate an integer
programming problem (formula see Appendix) using the best results from LLM
with text alone and TFIDF with text and abstract as the similarity parameters.
The similarity metric indicates that papers $i$ and $j$ sharing the same
cluster exhibit a similarity value of 1, and 0 otherwise. completeness score
and homogeneity score are outlined in Table 2. Notably, the comprehensive and
homogeneity scores align LLM’s clustering result closely with TFIDF. This
alignment reflects a common scenario in optimization problems where the
parameters are derived from predicted outcomes.

Approaches |  Completeness
Score |  Homogeneity
Score
---|---|---
TFIDF (title only) |  0.33 |  0.14
TFIDF (title and abstract) |  0.37 |  0.20
LLM (title only) |  0.41 |  0.22
Table 1. clustering results applying LLM and TFIDF. This is the average results from 5 random trials. LLM with title only exhibits the highest performance in completeness score and homogeneity score when considering the current session label as the ground truth. Approaches |  Completeness
Score |  Homogeneity
Score
---|---|---
TFIDF (title and abstract) |  0.42 |  0.42
LLM (title only) |  0.41 |  0.41
Table 2. The results of integer programming optimization using similarity
measures from LLM and TFIDF. Due to the inferior performance of TFIDF with
title-only, it is omitted from this comparison. The results show that the LLM
has a similar result to the TFIDF(title and abstract).

##  4\. Conclusion

Our experimentation shows that LLMs are a promising tool in managing
optimization problems even when they lack precision when the number of
decision variables exceeds a certain threshold. To address this limitation, a
potential strategy involves involving humans in the loop, or combining the
strengths of LLMs and numerical solvers as we demonstrated through our
clustering + Integer Programming approach.

##  5\. Acknowledgements

We would like to thank Mercari Inc. for supporting this research. LLMs were
used to reword sections and to reformat tables in the paper. However, we did
not generate entire sections through LLMs.

## References

  * (1)
  * Levonian et al. (2023) Zachary Levonian, Chenglu Li, Wangda Zhu, Anoushka Gade, Owen Henkel, Millie-Ellen Postle, and Wanli Xing. 2023.  Retrieval-augmented Generation to Improve Math Question-Answering: Trade-offs Between Groundedness and Human Preference.  https://doi.org/10.48550/arXiv.2310.03184 arXiv:2310.03184 [cs].
  * Rosenberg and Hirschberg (2007) Andrew Rosenberg and Julia Hirschberg. 2007.  V-measure: A conditional entropy-based external cluster evaluation measure. In _Proceedings of the 2007 joint conference on empirical methods in natural language processing and computational natural language learning (EMNLP-CoNLL)_. 410–420.
  * Wei et al. (2022) Jason Wei, Yi Tay, Rishi Bommasani, Colin Raffel, Barret Zoph, Sebastian Borgeaud, Dani Yogatama, Maarten Bosma, Denny Zhou, Donald Metzler, Ed H. Chi, Tatsunori Hashimoto, Oriol Vinyals, Percy Liang, Jeff Dean, and William Fedus. 2022.  Emergent Abilities of Large Language Models.  http://arxiv.org/abs/2206.07682 arXiv:2206.07682 [cs].
  * Yang et al. (2023) Chengrun Yang, Xuezhi Wang, Yifeng Lu, Hanxiao Liu, Quoc V. Le, Denny Zhou, and Xinyun Chen. 2023.  Large Language Models as Optimizers.  https://doi.org/10.48550/arXiv.2309.03409 arXiv:2309.03409 [cs].
  * Zhang et al. (2024) Tianjun Zhang, Aman Madaan, Luyu Gao, Steven Zheng, Swaroop Mishra, Yiming Yang, Niket Tandon, and Uri Alon. 2024.  In-Context Principle Learning from Mistakes.  https://doi.org/10.48550/arXiv.2402.05403 arXiv:2402.05403 [cs].

##  Appendix A Prompt Design

The following is the template of the prompt we used. Note that the terms
surrounded by curly braces ({}) show where the list of papers was added.

    Program creation is the process of taking all the
    accepted papers to a conference and allocating a
    presentation slot for each paper with parallel sessions.
    The PC chairs of a conference typically do this manually.
    Assign sessions to the following papers based on
    the following constraints:
    1. The total length of all paper presentations
       within a session should be less than or equal to
       the length of the session they are in.
    2. No new sessions should be added.
    3. The "Discussions and Q/A" event should only be
       at the end of a session if present.
    4. All papers must be assigned to some session.

    The output should contain the schedule in the form
    of the csv representation of a data frame. This
    csv representation should be in three quotes (‘‘‘)
    on both sides so that I can easily extract it
    from your result and make a data frame.
    Example output format(use as many rows as the
    actual number of papers):

    ‘‘‘
    session@talk_title@duration
    231@An Empirical Study on Maintainable Method ...@7
    223@Improve Quality of Cloud Serverless ...@7
    15@Extracting corrective actions from code ...@7
    15@How to Improve Deep Learning for Software ...@7
    11@ReCover: a Curated Dataset for Regression...@4
    ‘‘‘

    The list of paper and session titles are below:

    Session Lengths:
    {sessions_df_string}

    Paper durations:
    {papers_df_string}

    Make sure, above all else, that your response
    is formatted as requested, with the proper headers.

Some points to note:

  * •

On top of giving the list of papers, we gave a list of sessions with session
length.

  * •

We shuffled the list of papers so the LLM would not copy the original
ordering.

  * •

We opted for values separated by the commercial at symbol (@) rather than
commas to make handling titles with commas easier.

  * •

The final statement was added to increase the chances of getting the output in
the right format.

##  Appendix B Integer Programming Formulation

Define binary decision variables $x_{i,j}$:

$x_{i,j}=\begin{cases}1&\text{if paper }i\text{ is scheduled in session }j,\\\
0&\text{otherwise}.\end{cases}$

We introduce a new variable to indicate whether two papers are scheduled for
the same session. Define binary decision variables $z_{i,j,m}$:

$z_{i,j,m}=\begin{cases}1&\text{if papers }i\text{ and }j\text{ are in the
same session }m,\\\ 0&\text{otherwise}.\end{cases}$

###  B.1. Objective Function

Maximize the similarity between papers scheduled in the same session (if the
paper is clustered into the same cluster, the similarity is 1 otherwise it is
0):

$\text{Maximize}\quad\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{m=1}^{M}\text{simlarity}{(i,j)}\cdot
z_{i,j,m}$

###  B.2. Constraints

####  B.2.1. Link $z$ and $x$ variables

For each pair of papers both scheduled in the same session:

$z_{i,j,m}\geq x_{i,m}+x_{j,m}-1$ $z_{i,j,m}\leq x_{i,m}$ $z_{i,j,m}\leq
x_{j,m}$

####  B.2.2. Each paper scheduled exactly once

Ensure every paper is scheduled in exactly one session:

$\sum_{j=1}^{M}x_{i,j}=1\quad\forall i$

####  B.2.3. Session length constraints

Ensure the total duration of papers in any session does not exceed the
available time:

$\sum_{i=1}^{N}x_{i,j}\cdot\text{papers\\_duration}_{i}\leq\text{session\\_lengths}_{j}\quad\forall
j$


$\displaystyle\mathcal{QLQ}=\mathcal{L}_{S}\mathcal{Q}+\mathcal{L}_{E}\mathcal{Q}+\alpha\mathcal{QL}_{I}\mathcal{Q}.$
(85d)

We then substitute these expressions into Eq. (80b) and solve for
$\mathcal{Q}\rho$ as before, obtaining

$\displaystyle\mathcal{Q}\tilde{\rho}(s)=\frac{1}{s-(\mathcal{L}_{S}+\mathcal{L}_{E}+\alpha\mathcal{QL}_{I})}\alpha\mathcal{L}_{I}\mathcal{P}\tilde{\rho}(s)$
(86)

where we have once again ignored the term $\propto\mathcal{Q}\rho(0)$, as this
vanished for factorising initial conditions. This can be substituted once more
into Eq. (80a) alongside Eqs. (85) to find

$\displaystyle s\mathcal{P}\tilde{\rho}(s)-$
$\displaystyle\mathcal{P}\rho(0)=\mathcal{P}\mathcal{L}_{S}\mathcal{P}\tilde{\rho}(s)$
(87)
$\displaystyle+\alpha^{2}\mathcal{PL}_{I}\frac{1}{s-(\mathcal{L}_{S}+\mathcal{L}_{E}+\alpha\mathcal{QL}_{I})}\mathcal{QL}_{I}\mathcal{P}\tilde{\rho}(s).$

we then take the lowest order in $\alpha$ in each the homogeneous and
inhomogeneous terms, Gardiner and Zoller (2000)

$\displaystyle s\mathcal{P}\tilde{\rho}(s)$
$\displaystyle=\mathcal{P}\rho(0)+\mathcal{P}\mathcal{L}_{S}\mathcal{P}\tilde{\rho}(s)$
(88)
$\displaystyle+\alpha^{2}\mathcal{PL}_{I}\frac{1}{s-(\mathcal{L}_{S}+\mathcal{L}_{E})}\mathcal{L}_{I}\mathcal{P}\tilde{\rho}(s).$

Finally, once again exploiting the convolution property of the Laplace
transform we may take its inverse and trace over the environment to find

$\displaystyle\frac{d}{dt}\rho_{S}(t)$
$\displaystyle=\mathcal{L}_{S}\rho_{S}(t)$
$\displaystyle+\alpha^{2}\int_{0}^{t}d\tau\operatorname{Tr}_{E}[\mathcal{PL}_{I}e^{(\mathcal{L}_{S}+\mathcal{L}_{E})\tau}\mathcal{L}_{I}\mathcal{P}\rho(t-\tau)]$
$\displaystyle=-i[H_{S},\rho_{S}(t)]$ (89)
$\displaystyle-\alpha^{2}\int_{0}^{t}d\tau
Tr_{E}[[H_{I}(t),[H_{I}(t-\tau),\rho_{S}(t-\tau)\otimes\rho_{B}]]],$

where $H_{I}(\tau)=e^{(\mathcal{L}_{S}+\mathcal{L}_{E})\tau}H_{I}$ is the
interaction Hamiltonian in the interaction picture. We thus recover the
familiar Bloch-Redfield equation from which we may straightforwardly obtain
the GKSL master equation Breuer and Petruccione (2002).


# Bayesian Intervention Optimization for Causal Discovery

Yuxuan Wang1
&Mingzhou Liu1
&Xinwei Sun2
&Wei Wang311footnotemark: 1
&Yizhou Wang1
Correspondence to<EMAIL_ADDRESS><EMAIL_ADDRESS>
(

1 School of Computer Science, Peking University
2 School of Data Science, Fudan University
3 National Key Laboratory of General Artificial Intelligence, BIGAI
)

###### Abstract

Causal discovery is crucial for understanding complex systems and informing
decisions. While observational data can uncover causal relationships under
certain assumptions, it often falls short, making active interventions
necessary. Current methods, such as Bayesian and graph-theoretical approaches,
do not prioritize decision-making and often rely on ideal conditions or
information gain, which is not directly related to hypothesis testing. We
propose a novel Bayesian optimization-based method inspired by Bayes factors
that aims to maximize the probability of obtaining decisive and correct
evidence. Our approach uses observational data to estimate causal models under
different hypotheses, evaluates potential interventions pre-experimentally,
and iteratively updates priors to refine interventions. We demonstrate the
effectiveness of our method through various experiments. Our contributions
provide a robust framework for efficient causal discovery through active
interventions, enhancing the practical application of theoretical
advancements.

##  1 Introduction

Causal discovery is fundamental for advancing Artificial General Intelligence
(AGI), enabling agents to build comprehensive causal world models essential
for understanding complex systems and testing hypotheses [19]. While
observational data is abundant and easy to collect, it often falls short in
uncovering true causal relationships due to assumptions like the Markov
property and faithfulness [21]. Active interventions are necessary to recover
structural causal models, but they can be prohibitively costly in real-world
scenarios, such as genetic knockout experiments used to discover causal links
between genes and biological traits. Therefore, efficient active intervention
algorithms that can quickly and reliably uncover causal relationships are
urgently needed.

To address these challenges, we propose a method inspired by the Bayes factor
[11], widely used in hypothesis testing across various fields including
Bayesian statistics [24], physics [1], and neuroscience [12]. The Bayes factor
provides a comparative measure of evidence for two competing hypotheses based
on observed data, making them valuable for hypothesis testing in uncertain and
data-limited environments [11]. Our approach aims to maximize the probability
of achieving decisive and correct evidence, thus directly supporting
hypothesis testing.

Our method leverages observational data to estimate causal models under
different hypotheses and computes a pre-experimental evaluation of the
potential success of interventions. By identifying the most effective
interventions and iteratively updating our priors using Bayesian methods, we
refine our interventions to increase the likelihood of obtaining strong and
truthful evidence through Bayesian optimization. This method addresses the
shortcomings of existing approaches by providing a direct measure for
acquiring decisive and correct evidence rather than relying on information
gain, which does not necessarily translate to hypothesis-testing confidence.

#### Related Works

Existing approaches to active intervention can be broadly categorized into
Bayesian and graph-theoretical methods [7]. The Bayesian approach, pioneered
by Tong [23] and Murphy [17], typically involves selecting an intervention
that maximizes a Bayesian utility function, often using mutual information
between the graph and the post-intervention samples. While these methods have
advanced to consider multiple simultaneous interventions [22], budget
constraints [2], and expert knowledge [16], they still use information gain,
which is not directly related to hypothesis testing.

Graph-theoretical methods explore selecting interventions to recover the
maximum number of edges in the graph [9, 8], minimize intervention costs [13,
15], and address hidden variables [14]. However, these methods often assume
perfect identification of the Markov equivalence class and infinite
interventional data, which are impractical in real-world scenarios.
Additionally, current active intervention methods predominantly use
interventional data without fully utilizing observational data.

Our approach addresses these limitations by integrating observational data to
model intervention distribution and focusing on identifying decisive and
correct hypotheses. By leveraging the Bayes factor, which provides a
comparative measure of evidence for choosing the null or alternative
hypotheses, we guide our interventions more effectively. Additionally, by
utilizing Bayesian optimization, our framework offers a more direct and robust
approach to causal discovery, ensuring that our interventions are grounded in
strong statistical evidence.

Our contributions can be summarized as follows:

  1. 1.

We propose an active intervention algorithm that maximizes the probability of
obtaining decisive and correct evidence in causal discovery.

  2. 2.

We develop a Monte Carlo estimation-based Bayesian optimization method for
estimating and optimizing the probability of identifying decisive and correct
evidence.

  3. 3.

We demonstrate the effectiveness of our method compared to existing approaches
using generalized linear models and diverse noise distributions.

##  2 Preliminary

#### Notation

In our study, we consider a Directed Acyclic Graph (DAG) $G=(V,E)$, where
$V=\\{1,\dots,d\\}$ represents the set of vertices, and each vertex
corresponds to a random variable. These variables, denoted as
$\mathbf{X}_{V}=\\{X_{1},\dots,X_{d}\\}\in\mathcal{X}^{d}$, are indexed by
$V$. Our initial observational dataset
$\mathbf{D}=\\{\mathbf{x}_{V}^{(i)}\\}_{i=1}^{n}$ consists of instances drawn
from the probability distribution $p_{G}(x_{1},\dots,x_{d})$, representing the
joint distribution of these variables.

#### Structural Causal Models (SCM)

The underlying data generative mechanism of a DAG $G$ can be descirbed by a
set of structural equations:

$X_{i}\coloneqq f_{i}(X_{\mathrm{pa}_{G}(i)},\varepsilon_{i})\quad\forall i\in
V$ (1)

Here, $f_{i}$’s are the causal mechanisms that remain invariant under
interventions, $\mathrm{pa}_{G}(i)$ denotes parents of variable $X_{i}$ in DAG
$G$, and $\varepsilon_{i}$’s are mutually independent exogenous noise
variables. These structural equations embody both the conditional
distributions of random variables and the potential impact of interventions.
While these functions can be nonparametric, we often consider parametric
approximations with parameters $\gamma$. A typical example of SCM is the
Additive Noise Model (ANM), where the structural equation for each variable is
a function of its parents plus a random noise independent with the parent
variables:

$X_{i}\coloneqq
f_{i}(X_{\mathrm{pa}_{G}(i)};\gamma_{i})+\varepsilon_{i}\quad\varepsilon_{i}\rotatebox[origin={c}]{90.0}{$\models$}X_{\mathrm{pa}_{G}(i)}\quad\forall
i\in V$ (2)

#### Do-calculus

Do-calculus [18] is a framework used to determine the effect of interventions
on a set of variables within a causal model. Specifically, an intervention,
denoted by $\mathrm{do}(X_{\mathbf{W}}=\mathbf{x}^{\prime}_{\mathbf{W}})$,
modifies the original distribution by setting the variables $X_{\mathbf{W}}$
to specific values $\mathbf{x}^{\prime}_{\mathbf{W}}$, effectively breaking
any causal dependencies from their parents. The resulting interventional
distribution is given by the truncated factorization formula:

$\displaystyle
P(\mathbf{X}_{V}\mid\mathrm{do}(\mathbf{X}_{W}=\mathbf{x}^{\prime}_{W})=\prod_{i\in
V\setminus W}P(X_{i}\mid
X_{\text{pa}_{G}(i)})\mathbb{I}(\mathbf{X}_{W}=\mathbf{x}^{\prime}_{W}),$

where $\mathbb{I}(\cdot)$ is an indicator function. A perfect intervention on
any variable $X_{j}$ removes all dependencies with its parents, resulting in a
modified DAG $G^{\prime}=(V,E\setminus\\{(\text{pa}_{G}(j),j)\\})$. This
approach allows us to isolate the causal effect of the intervention on the
target variables, providing a clearer understanding of the underlying causal
relationships.

#### Hypothesis testing with Bayes factors

The Bayes factor [11] is an important tool in Bayesian hypothesis testing,
providing a comparative measure of the evidence for two competing hypotheses.
Given two hypotheses $\mathbb{H}_{0}$ and $\mathbb{H}_{1}$, and observed data
$\mathbf{D}$, the Bayes factor BF is defined as the ratio of the likelihoods
of the data under each hypothesis:

$\text{BF}_{01}=\frac{P\left(\mathbf{D}\mid\mathbb{H}_{0}\right)}{P\left(\mathbf{D}\mid\mathbb{H}_{1}\right)}$
(3)

This ratio assesses the relative support for each hypothesis provided by the
data. A value of $\text{BF}_{01}>1$ suggests stronger evidence in favor of
$\mathbb{H}_{0}$. This method is particularly valuable in decision-making
processes where uncertainty and limited data are prevalent.

Table 1: Classification for the evidence levels of the Bayes factor $\textnormal{BF}_{01}$ (from [20] adapted from [10]). Bayes factor |  | Evidence Level
---|---|---
$>100$ |  | Extreme evidence for $\mathbb{H}_{0}$
$30-100$ |  | Very strong evidence for $\mathbb{H}_{0}$
$10-30$ |  | Strong evidence for $\mathbb{H}_{0}$
$3-10$ |  | Moderate evidence for $\mathbb{H}_{0}$
$1-3$ |  | Anecdotal evidence for $\mathbb{H}_{0}$
$1$ |  | No evidence
$1/3$ \- 1 |  | Anecdotal evidence for $\mathbb{H}_{1}$
$1/10$ \- 1/3 |  | Moderate evidence for $\mathbb{H}_{1}$
$1/30$ \- 1/10 |  | Strong evidence for $\mathbb{H}_{1}$
$1/100$ \- 1/30 |  | Very strong evidence for $\mathbb{H}_{1}$
$<1/100$ |  | Extreme evidence for $\mathbb{H}_{1}$

The Bayes factor can be classified into different categories of evidence
strength, which helps interpret its value in the context of supporting
$\mathbb{H}_{0}$ or $\mathbb{H}_{1}$ [6]. Table 1 summarizes these categories.
This classification provides a nuanced understanding of how compelling the
data is in favor of one hypothesis over the other. To simplify, consider the
evidence levels for $\mathbb{H}_{0}$ as described in the table above, and
apply symmetry for evidence levels supporting $\mathbb{H}_{1}$.

##  3 Methodology

In this section, we introduce the problem formulation first and detail our
methodology for uncovering causal relationships using active sampling and
Bayesian optimization within the Additive Noise Model framework. We focus on
determining whether there is a causal edge between two random variables by
formulating this problem as a hypothesis test. Specifically, the conditional
distributions after applying do-calculus for the hypotheses $\mathbb{H}_{0}$
and $\mathbb{H}_{1}$ are:

$\displaystyle\mathbb{H}_{0}:p(y\mid\mathrm{do}(X=x))=p(y)\quad\mathbb{H}_{1}:p(y\mid\mathrm{do}(X=x))=p(y\mid
x)$ (4)

When $\mathbb{H}_{1}$ holds, the structural causal model is simply $X\to Y$.
If $\mathbb{H}_{0}$ holds, we assume $X\not\to Y$, meaning there is no direct
causal influence from $X$ to $Y$. Under this hypothesis, we consider several
possible relationships between $X$ and $Y$ (1)
$X\rotatebox[origin={c}]{90.0}{$\models$}Y$: $X$ and $Y$ are independent. (2)
$X\leftarrow Y$: There is a direct causal influence from $Y$ to $X$. (3)
$X\leftarrow U\to Y$: There exists an unobservable confounder $U$ that
influences both $X$ and $Y$. Among these scenarios, the case where
$X\rotatebox[origin={c}]{90.0}{$\models$}Y$ represents a situation of
independence, which is relatively straightforward and can be easily tested
using standard independence tests. Therefore, we will not focus on this
scenario in our subsequent analysis. Instead, we will concentrate on the more
complex relationships where $X$ and $Y$ are either causally related through
$Y\to X$ or via a common confounder $C$.

###  3.1 $P_{DC}$ as optimization objective

In this work, we define a specific metric, $P_{DC}$, to evaluate the
probability of obtaining decisive and correct evidence after performing a
specific intervention $\mathrm{do}(X=x)$. Let $\mathbf{D}_{\text{int}}$ denote
the already obtained interventional data. Our $P_{DC}$ metric is formulated to
determine the probability of achieving decisive and correct evidence by
combining the new data obtained from the intervention $\mathrm{do}(X=x)$ with
the existing interventional data $\mathbf{D}_{\text{int}}$. We define the
probabilities that the experiment yields a correct and decisive result under
the hypotheses $\mathbb{H}_{0}$ and $\mathbb{H}_{1}$ hold as $P_{DC}^{0}$ and
$P_{DC}^{1}$, respectively.

$\displaystyle P_{DC}^{0}(k_{0},\mathbf{D}_{\text{int}},\mathrm{do}(X=x))$
$\displaystyle:=P\left(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\mathbf{D}_{\text{new}})>k_{0}\mid\mathbf{D}_{\text{int}},\mathrm{do}(X=x),\mathbb{H}_{0}\right)$
(5) $\displaystyle P_{DC}^{1}(k_{1},\mathbf{D}_{\text{int}},\mathrm{do}(X=x))$
$\displaystyle:=P\left(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\mathbf{D}_{\text{new}})<k_{1}\mid\mathbf{D}_{\text{int}},\mathrm{do}(X=x),\mathbb{H}_{1}\right)$

Here, $\mathbf{D}_{\text{int}}$ represents the interventional data already
collected, and $\mathbf{D}_{\text{new}}$ represents the new data obtained from
the intervention $\mathrm{do}(X=x)$. In the definitions of $P_{DC}^{0}$ and
$P_{DC}^{1}$, $\mathbf{D}_{\text{new}}$ is sampled under the assumptions that
$\mathbb{H}_{0}$ and $\mathbb{H}_{1}$ hold true, respectively. Using
theseedefinitions, the overall $P_{DC}$ metric can be defined as:

$\displaystyle P_{DC}(k_{0},k_{1},\mathbf{D}_{\text{int}},\mathrm{do}(X=x)):=$
$\displaystyle
P_{DC}^{0}(k_{0},\mathbf{D}_{\text{int}},\mathrm{do}(X=x))P(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}})$
(6)
$\displaystyle+P_{DC}^{1}(k_{1},\mathbf{D}_{\text{int}},\mathrm{do}(X=x))P(\mathbb{H}_{1}\mid\mathbf{D}_{\text{int}})$

In this formulation,
$P_{DC}^{i}(k_{0},\mathbf{D}_{\text{int}},\mathrm{do}(X=x))$ represents the
probability of obtaining decisive and correct evidence after performing this
additional intervention $\mathrm{do}(X=x)$ and combining it with the existing
data and hypothesis $\mathbb{H}_{i}$. This allows us to select the $x$ that
maximizes $P_{DC}$, thereby obtaining decisive and correct statistical
evidence. Formally, solving the optimization problem:

$\max_{x}P_{DC}(k_{0},k_{1},\mathbf{D}_{\text{int}},\mathrm{do}(X=x))$ (7)

enables us to identify the intervention that maximizes the probability of
obtaining strong and correct evidence. We will describe the detailed algorithm
in the following sections.

#### Comparison with previous works on $P_{DC}$

Previous definitions of $P_{DC}$ in the field of Bayesian sample size
determination [6], such as $P_{DC}(k_{0},k_{1},n)$, primarily focused on
estimating the required sample size $n$ before conducting experiments to
ensure a high probability of obtaining decisive and correct evidence. These
methods emphasized determining the number of observations needed, rather than
specifying any intervention strategies:

$P_{DC}(k_{0},k_{1},n)=P_{DC}^{0}(\mathbf{y}^{n},k_{0}\mid\mathbb{H}_{0})P(\mathbb{H}_{0})+P_{DC}^{1}(\mathbf{y}^{n},k_{1}\mid\mathbb{H}_{1})P(\mathbb{H}_{1})$

While these methods aimed to provide substantial evidence for the correct
hypothesis, they did not account for active intervention in experiments.
Instead, they relied on passively observing the environment’s output, focusing
solely on the sample size required. Although [3] considered active
intervention in experiments, they modeled the sampling strategy as a Gaussian
distribution without optimizing the intervention itself to expedite causal
discovery.

Our approach, in contrast, optimizes the intervention variable
$\mathrm{do}(X=x)$ to maximize the probability of obtaining decisive and
correct evidence in a single experiment. This shift in perspective allows us
to design interventions that achieve the highest possible evidential support,
offering a more targeted and efficient approach to hypothesis testing.

Throughout the following of this paper, we will use the notation
$P_{DC}(\mathbf{D}_{\text{int}},\mathrm{do}(X=x))$ to refer to our defined
metric $P_{DC}(k_{0},k_{1},\mathbf{D}_{\text{int}},\mathrm{do}(X=x))$ when
emphasizing the intervention $\mathrm{do}(X=x)$ and the existing data
$\mathbf{D}_{\text{int}}$ and $P_{DC}$ for simplification and metric notation.

###  3.2 Optimization for the $P_{DC}$

In this subsection, we discuss our approach to optimizing the $P_{DC}$ metric
using Bayesian optimization. We employ Monte Carlo (MC) sampling to estimate
$P_{DC}$ and use gradient descent to identify the optimal intervention.
Specifically, we detail the process of estimating the distributions $m_{0}$
under $\mathbb{H}_{0}$ and $m_{1}$ under $\mathbb{H}_{1}$, and subsequently
computing the Bayes (BF) as defined in Eq. (3).

#### Estimation of $P_{DC}$

To estimate and optimize $P_{DC}$, we propose a Monte Carlo (MC) estimation
method that provides a differentiable surrogate estimator as follow. The
comprehensive derivation of this method is available in Appendix A.1.

$\displaystyle P_{DC}(\mathbf{D}_{\text{int}},\mathrm{do}(X=x))$
$\displaystyle\approx\frac{1}{N}\sum_{i=1}^{N}\exp\left(-\frac{1}{\beta}\text{ReLU}(k_{0}-\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\\{y_{0i},x\\}))\right)P(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}})$
(8)
$\displaystyle+\frac{1}{N}\sum_{i=1}^{N}\exp\left(-\frac{1}{\beta}\text{ReLU}(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\\{y_{1i},x\\})-k_{1})\right)P(\mathbb{H}_{1}\mid\mathbf{D}_{\text{int}})$

where $y_{0i}\sim m_{0}(y)$ is one of $N$ samples drawn from the estimated
distribution $m_{0}(\cdot)$ and $y_{1i}\sim m_{1}(y\mid x)$ is one of $N$
samples drawn from the estimated distribution $m_{1}(\cdot\mid x)$. $\beta$ is
the scale parameter to approximate the indicator function detailed in the
Appendix A.1. Similarly, the second term in the equation can be estimated by
modifying the indicator function to reflect the corresponding event under
$\mathbb{H}_{1}$ and sampling from the distribution $m_{1}(y\mid x)$.

#### Estimation of the interventional distribution

Based on our introduction in Section 3.1, we need to estimate the
interventional distributions under each hypothesis to sample from them and
compute the Bayes factor. Although the observational data has not undergone
do-calculus, it can be used to estimate the probability density functions
under each hypothesis. Specifically:

  * •

Under $\mathbb{H}_{0}$, the hypothesis implies that
$p(y)=p(y\mid\mathrm{do}(X=x),\mathbb{H}_{0})$, and we estimate
$m_{0}(y)\approx p(y)=p(y\mid\mathrm{do}(X=x),\mathbb{H}_{0})$.

  * •

Under $\mathbb{H}_{1}$, the hypothesis implies that $p(y\mid
x)=p(y\mid\mathrm{do}(X=x),\mathbb{H}_{1})$, and we estimate $m_{1}(y\mid
x)\approx p(y\mid x)=p(y\mid\mathrm{do}(X=x),\mathbb{H}_{1})$.

When $\mathbb{H}_{0}$ holds, we can estimate $m_{0}$ by performing maximum
likelihood estimation (MLE) on the observational data $\\{y_{i}\\}_{i=1}^{N}$.
Formally, this estimation is:

$\max_{\phi}\sum_{i=1}^{N_{\text{obs}}}\log m_{0}(y_{i})$

where $\\{y_{i},x_{i}\\}\in D_{\text{obs}}$. Since the mixture of normals can
approximate any distribution [5], we parameterize $m_{0}$ using a mixture of
normals and solve the above problem.

When $\mathbb{H}_{1}$ holds, we can write down the structural equations 3.2,
and we have

$\displaystyle X=n_{X},\quad Y=f(X,\gamma)+n_{Y}$
$p(y\mid\mathrm{do}(X=x),\mathbb{H}_{1})=p(y\mid x)=p_{Y}(y-f(x,\gamma))$

We similarly parameterize $m_{1}$ using a mixture of normals to approximate
the ground truth interventional distribution as:

$m_{1}(y\mid x)=p_{\psi}(y-f(x,\theta))\approx
p(y\mid\mathrm{do}(X=x),\mathbb{H}_{1})$

where $\theta$ and $\psi$ are the parameters of the link function and the
mixture of normals for $m_{1}$, respectively. To estimate the interventional
distribution under $\mathbb{H}_{1}$, we solve the following optimization
problem:

$\displaystyle\max_{\theta,\psi}\log
p(X,Y)\Leftrightarrow\max_{\theta,\psi}\sum_{i=1}^{N_{\text{obs}}}\log
p(x_{i})+\log m_{1}(y_{i}\mid
x_{i})\Leftrightarrow\max_{\theta,\psi}\sum_{i=1}^{N_{\text{obs}}}\log
m_{1}(y_{i}\mid x_{i})$

where $\\{y_{i},x_{i}\\}\in D_{\text{obs}}$. By parameterizing $m_{0}$ and
$m_{1}$ using mixtures of normals, we effectively approximate the
interventional distributions under $\mathbb{H}_{0}$ and $\mathbb{H}_{1}$,
respectively. This allows us to sample from these estimated distributions and
compute the Bayes factor necessary for our analysis.

#### Estimation of Bayes Factor

We employ a non-parametric method to estimate the Bayes factor. Given that we
have obtained the maximum likelihood estimates (MLE) $\hat{m}_{0}$ and
$\hat{m}_{1}$ for the densities $m_{0}$ and $m_{1}$ respectively, we can
compute the likelihoods directly. Using these likelihoods, we calculate the
Bayes factor for $N_{\text{int}}$ interventional data points
$\mathbf{D}_{\text{int}}={y_{i},x_{i}}_{i=1}^{N_{\text{int}}}$ as follows:

$\displaystyle\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\\{y,x\\})$
$\displaystyle=\frac{P(\mathbf{D}_{\text{int}}\cup\\{y,x\\}\mid\hat{m}_{0},\mathbb{H}_{0})}{P(\mathbf{D}_{\text{int}}\cup\\{y,x\\}\mid\hat{m}_{1},\mathbb{H}_{1})}=\prod_{i=1}^{N_{\text{int}}}\frac{\hat{m}_{0}(y_{i})}{\hat{m}_{1}(y_{i}\mid
x_{i})}\times\frac{\hat{m}_{0}(y)}{\hat{m}_{1}(y\mid x)}$

With these estimates, we can optimize $P_{DC}$ through Bayesian optimization.
Based on the optimization results, we actively intervene in the environment to
obtain interventional samples. Once new samples are collected, we need to
update our estimates of $P_{DC}$ refining our calculations accordingly.

#### Update of the prior of the hypothesis

As we actively engage with the environment and collect interventional data, we
update the posterior for the hypothesis. This is crucial for calculating
$P_{DC}$ in Eq. (8) and evaluating the posterior of the ground truth given the
interventional data. Using Bayes’ theorem, the posterior updates are:

$\displaystyle P\left(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}}\right)$
$\displaystyle=\frac{P\left(\mathbf{D}_{\text{int}}\mid\mathbb{H}_{0}\right)P(\mathbb{H}_{0})}{\sum_{i\in[0,1]}P\left(\mathbf{D}_{\text{int}}\mid\mathbb{H}_{i}\right)P(\mathbb{H}_{i})}\quad
P\left(\mathbb{H}_{1}\mid\mathbf{D}_{\text{int}}\right)$
$\displaystyle=\frac{P\left(\mathbf{D}_{\text{int}}\mid\mathbb{H}_{1}\right)P(\mathbb{H}_{1})}{\sum_{i\in[0,1]}P\left(\mathbf{D}_{\text{int}}\mid\mathbb{H}_{i}\right)P(\mathbb{H}_{i})}$
(9)

Algorithm 1 Actively intervention for hypothesis testing

1: Input: Observational data $\mathbf{D}_{\text{obs}}$; intervention sample
size $M$

2: Output:
$\text{BF}_{01}(\mathbf{D}_{\text{int}}),P\left(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}}\right),P\left(\mathbb{H}_{1}\mid\mathbf{D}_{\text{int}}\right)$

3: $\mathbf{D}_{\text{int}}=\emptyset$

4:
$P\left(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}}\right)=P\left(\mathbb{H}_{1}\mid\mathbf{D}_{\text{int}}\right)=\frac{1}{2}$

5: Estimate the interventional distribution for the hypothese
$\hat{m}_{0},\hat{m}_{1}$ with the observational data
$\mathbf{D}_{\text{obs}}$ as described in 3.2

6: for $m=1,\dots,M$ do

7: Estimate $P_{DC}(\mathbf{D}_{\text{int}},\mathrm{do}(X=x))$ using Monte
Carlo method as described in Eq. (8)

8: Optimize $x_{\text{opt}}$ to maximize
$P_{DC}(\mathbf{D}_{\text{int}},\mathrm{do}(X=x_{\text{opt}}))$

9: Sample $y\mid\mathrm{do}(X=x_{\text{opt}})\sim
p_{G}(y\mid\mathrm{do}(X=x_{\text{opt}}))$ from the environment with
intervention

10:
$\mathbf{D}_{\text{int}}\leftarrow\mathbf{D}_{\text{int}}\cup\\{(y,x_{\text{opt}})\\}$

11: Update
$P\left(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}}\right),P\left(\mathbb{H}_{1}\mid\mathbf{D}_{\text{int}}\right)$
with Eq. (9)

12: end for

Now we can effectively calculate the Probability of Decisive and Correct
evidence ($P_{DC}$). Due to the differentiability of the entire computational
process, we can directly utilize some gradient-based methods for optimization.
This capability allows us to precisely refine our intervention strategies and
improve the efficiency and accuracy of our causal discovery process. In
summary, we propose a Bayesian optimization-based active intervention
algorithm, outlined in Algorithm 1. The algorithm can be divided into four
main steps:

  1. 1.

Data-Dependent Prior Formulation: Using the observational dataset
$\mathbf{D}_{\text{obs}}$ to inform the creation of priors.

  2. 2.

Monte Carlo Estimation: Utilizing the formulated priors in a Monte Carlo
sampling process to estimate the Probability of Decisive and Correct evidence
($P_{DC}$) for our Bayes factor analysis.

  3. 3.

Optimization of Intervention Strategy: Optimizing the intervention strategy
$\mathrm{do}(X=x)$ to maximize $P_{DC}$.

  4. 4.

Active Sampling and Updating Priors: Actively sampling from the environment
using the optimized intervention strategy and integrating the newly collected
interventional data ($\mathbf{D}_{\text{int}}$) to update the priors.

##  4 Experiment

In this section, we present a series of experiments designed to evaluate the
effectiveness of our proposed causal inference methodology. These experiments
are conducted on synthetic datasets that represent a variety of causal
structures and noise distributions. We generate 5000 observational data points
without any intervention for each experimental setup to form the basis of our
analysis.

#### Compared methods

We compare our method against the following baselines:

  * •

Random Sampling: Interventions are chosen uniformly at random.

  * •

ABCD [2]/CBED [22] Strategies: These methods maximize information gain to
determine the optimal intervention point. One can refer to Appendix A.2 for
implementation details.

  * •

Our Method: A strategy that maximizes the probability of decisive and correct
evidence ($P_{DC}$) to determine the intervention value.

#### Evaluation metrics

We evaluate the performance of each method using the following metrics:

  * •

$P_{DC}$: Probability of Decisive and Correct Evidence. A higher $P_{DC}$
indicates a greater chance that the intervention correctly and decisively
determines the causal relationship.

  * •

$\log\text{BF}_{01}$: The logarithm of the Bayes factor, representing the
likelihood ratio of the interventional data under the null hypothesis to that
under the alternative hypothesis. Higher values indicate stronger evidence for
the null hypothesis, while lower values support the alternative hypothesis.

  * •

$P(\mathbb{H}_{\text{gt}}\mid\mathbf{D}_{\text{int}})$: Posterior probability
of the ground truth hypothesis given the interventional data. Higher values
indicate a better likelihood of correctly identifying the true causal
relationship, thus are more favorable.

#### Synthetic data generation

As described in Section 3, we investigate two hypotheses regarding whether
there exists direct causal relationship from $X$ to $Y$. Previously, we
discussed the triviality of the independence scenario
($X\rotatebox[origin={c}]{90.0}{$\models$}Y$) under $\mathbb{H}_{0}$. Here, we
focus on the remaining possible relationships under $\mathbb{H}_{0}$: (1)
$X\leftarrow Y$: There is a direct causal influence from $Y$ to $X$. (2)
$X\leftarrow U\to Y$: There exists an unobservable confounding variable $U$
that influences both $X$ and $Y$.

For $\mathbb{H}_{1}$ ($X\to Y$), the structural equations are:

$\displaystyle X=n_{X},\quad Y=f(X,\gamma)+n_{Y}$

where $f$ is parameterized by $\gamma$. In our experiments, we employ a
generalized linear model with an invertible link function, specifically
$f=\tanh$. The function is parameterized as $A\cdot\tanh(BX)$, with $A=2$ and
$B=1$.

When considering $\mathbb{H}_{0}$ under the scenario $X\leftarrow Y$, the
structural equations are:

$\displaystyle Y=n_{Y},\quad X=f(Y,\gamma)+n_{X}$

using the same parameterization for $f$ as the ground truth of
$\mathbb{H}_{1}$.

For the confounding scenario $X\leftarrow U\to Y$, the structural equations
become:

$\displaystyle U=n_{U},\quad Y=f(U,\gamma)+n_{Y},\quad X=g(U,\beta)+n_{X}$

where $f$ and $g$ are parameterized similarly, with $f=A\cdot\tanh(BU)$ and
$g=C\cdot\tanh(DU)$ , ensuring the unobservable confounder $U$ influences both
$X$ and $Y$. We set $A=2,B=1,C=2,D=1$. In each ground truth scenario, we
generate 5000 observational data points based on the specified structural
equations. These observational datasets are then provided to all baseline
methods for consistent evaluation and comparison.

To represent a wide range of noise distributions in our experiments, we use a
mixture of normals due to its flexibility in approximating various
distributions [5]. For each experimental setup, we randomly generate three
distributions to ensure a diverse range of noise characteristics. To prevent
the distributions from degenerating into overly simple forms, we fix the
position and variance of each component within the distributions while
randomizing their weights. The weights for the $n$ components are generated
using a combination of uniform and softmax distributions. The mixture model is
defined as follows:

$n_{Y}=\sum_{i=1}^{k}\pi_{i}\mathcal{N}(\mu_{i},\sigma_{i}^{2}),\quad\mathbf{\pi}=\frac{1}{2n}\mathbf{1}+\frac{1}{2}\text{Softmax}(\mathbf{z})$

where $k$ is the number of components, $\pi_{i}$ are the weights, $\mu_{i}$
are the means, and $\sigma_{i}^{2}$ are the variances of the components,
$\mathbf{1}_{n}$ is a $n$-dimension vector of ones and
$\mathbf{z}=(z_{1},z_{2},\ldots,z_{n})$ are $n$ independent samples from a
standard normal distribution $\mathcal{N}(0,I_{n})$. This approach ensures
that the weights are both diverse and appropriately normalized, thus avoiding
degeneracy. This method allows us to simulate a wide variety of noise
distributions, enhancing the robustness of our experiments. For each
experimental setup, the ground truth noise distribution is randomly chosen
from the three randomly generated distributions to mitigate the impact of
randomness. We also conducted experiments to show that even when the position
and variance of each component are randomized, the results and conclusions
remain consistent. All experiments are repeated with 10 different random seeds
to ensure the robustness and reproducibility of our results.

Figure 1: Results under different ground truths with
$k_{0}=\frac{1}{k_{1}}=10$: $P_{DC}$, $\log\text{BF}_{01}$, and
$P(\mathbb{H}_{gt}\mid\mathbf{D}_{\text{int}})$. The first row corresponds to
$\mathbb{H}_{0}$ ($X\leftarrow Y$)). The second row corresponds to
$\mathbb{H}_{0}$ ($X\leftarrow U\to Y$), and the last row corresponds to
$\mathbb{H}_{1}$ ($X\to Y$).

Figure 2: Results under different ground truths with
$k_{0}=\frac{1}{k_{1}}=10$: $P_{DC}$, $\log\text{BF}_{01}$, and
$P(\mathbb{H}_{gt}\mid\mathbf{D}_{\text{int}})$. But with random mean and
random covariance compared to the settings in the Figure 1.The first row
corresponds to $\mathbb{H}_{0}$ ($X\leftarrow Y$)). The second row corresponds
to $\mathbb{H}_{0}$ ($X\leftarrow U\to Y$), and the last row corresponds to
$\mathbb{H}_{1}$ ($X\to Y$).

#### Results

We delve into the outcomes of our tests corresponding to the earlier described
experimental setups. For each hypothesis, we present graphs that illustrate
how $P_{DC}$, $\log\text{BF}_{01}$, and
$P(\mathbb{H}_{gt}\mid\mathbf{D}_{\text{int}})$ evolve with the number of
samples. Both $P_{DC}$ and $P(\mathbb{H}_{gt}\mid\mathbf{D}_{\text{int}})$ are
higher values indicating better performance, as described in Section 2. The
$\log\text{BF}_{01}$ is evaluated such that higher values are better when
$\mathbb{H}_{0}$ is the ground truth and lower values are better when
$\mathbb{H}_{1}$ is the ground truth. Each graph shows the mean and standard
deviation across all experiments and random seeds, with solid lines
representing the mean and shaded areas representing the standard deviation.

To demonstrate the effectiveness of our method in obtaining different levels
of $P_{DC}$, we present results for various thresholds $k_{0}$ and $k_{1}$,
specifically $k_{0}=\frac{1}{k_{1}}=10$. Results for
$k_{0}=\frac{1}{k_{1}}=30$ and $k_{0}=\frac{1}{k_{1}}=100$ are provided in
Appendix B. This includes results for different confounders and robustness
tests with random means, variances, and weights.

Our results in Figure 1 reveal that when $\mathbb{H}_{0}$ holds, our algorithm
consistently outperforms information gain and random sampling methods in most
scenarios and metrics. Specifically, our method demonstrates faster
improvement (or decline, depending on whether the metric is maximized or
minimized) compared to the baselines. This indicates that our approach is more
effective in conducting interventions and achieving decisive and correct
evidence. Under $\mathbb{H}_{1}$, our method performs comparably to existing
methods. This could be due to the problem’s complexity when $\mathbb{H}_{1}$
is true, where accumulating a decisive Bayes factor of the alternative
hypothesis over the null hypothesis while ensuring correctness is a slower
process, leading to a smaller performance gap.

Additionally, our method shows robustness in various randomized settings,
including those with random means, variances, and weights as shown in Figure
2. While the gap between our method and the baselines decreases in these
scenarios, our approach still maintains a competitive edge. This robustness is
critical in practical applications where noise and variability are inherent,
demonstrating that our method can adapt to different data characteristics and
still perform effectively. Our approach’s ability to handle diverse conditions
while achieving consistent results highlights its practical utility and
reliability in real-world causal discovery tasks.

##  5 Conclusion and Limitations

#### Conclusion

In this paper, we propose an active sampling and Bayesian optimization-based
methodology for hypothesis testing, optimizing the probability of decisive and
correct evidence ($P_{DC}$). Through extensive experiments on synthetic data,
we demonstrate that our method consistently outperforms baseline strategies,
including random sampling and information gain-based method in the aspect of
the probability of getting the decisive and correct evidence.

#### Limitations

A key limitation of our method is the need for prior knowledge of the
structural equation’s link function when modeling $\hat{m}_{0}(y)$ and
$\hat{m}_{1}(y\mid x)$. This requirement can be restrictive in scenarios where
such information is unavailable. Future work could explore approaches that
require less specific prior knowledge, potentially employing non-parametric
techniques.

## References

  * [1] Elcio Abdalla, Guillermo Franco Abellán, Amin Aboubrahim, Adriano Agnello, Özgür Akarsu, Yashar Akrami, George Alestas, Daniel Aloni, Luca Amendola, Luis A Anchordoqui, et al.  Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies.  Journal of High Energy Astrophysics, 34:49–211, 2022.
  * [2] Raj Agrawal, Chandler Squires, Karren Yang, Karthikeyan Shanmugam, and Caroline Uhler.  Abcd-strategy: Budgeted experimental design for targeted causal structure discovery.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 3400–3409. PMLR, 2019.
  * [3] Federico Castelletti and Guido Consonni.  Bayesian sample size determination for causal discovery.  arXiv preprint arXiv:2206.00755, 2022.
  * [4] Thomas M Cover.  Elements of information theory.  John Wiley & Sons, 1999.
  * [5] SR Dalal and WJ Hall.  Approximating priors by mixtures of natural conjugate priors.  Journal of the Royal Statistical Society: Series B (Methodological), 45(2):278–286, 1983.
  * [6] Fulvio De Santis.  Statistical evidence and sample size determination for bayesian hypothesis testing.  Journal of statistical planning and inference, 124(1):121–144, 2004.
  * [7] Juan L Gamella and Christina Heinze-Deml.  Active invariant causal prediction: Experiment selection through stability.  Advances in Neural Information Processing Systems, 33:15464–15475, 2020.
  * [8] Alain Hauser and Peter Bühlmann.  Two optimal strategies for active learning of causal models from interventional data.  International Journal of Approximate Reasoning, 55(4):926–939, 2014.
  * [9] Yang-Bo He and Zhi Geng.  Active learning of causal networks with intervention experiments and optimal designs.  Journal of Machine Learning Research, 9(Nov):2523–2547, 2008.
  * [10] Harold Jeffreys.  Theory of Probability (3rd Edition).  Oxford, University Press, 1961.
  * [11] Robert E Kass and Adrian E Raftery.  Bayes factors.  Journal of the american statistical association, 90(430):773–795, 1995.
  * [12] Christian Keysers, Valeria Gazzola, and Eric-Jan Wagenmakers.  Using bayes factor hypothesis testing in neuroscience to establish evidence of absence.  Nature neuroscience, 23(7):788–799, 2020.
  * [13] Murat Kocaoglu, Alex Dimakis, and Sriram Vishwanath.  Cost-optimal learning of causal graphs.  In International Conference on Machine Learning, pages 1875–1884. PMLR, 2017.
  * [14] Murat Kocaoglu, Karthikeyan Shanmugam, and Elias Bareinboim.  Experimental design for learning causal graphs with latent variables.  Advances in Neural Information Processing Systems, 30, 2017.
  * [15] Erik Lindgren, Murat Kocaoglu, Alexandros G Dimakis, and Sriram Vishwanath.  Experimental design for cost-aware learning of causal graphs.  Advances in Neural Information Processing Systems, 31, 2018.
  * [16] Andrés R Masegosa and Serafín Moral.  An interactive approach for bayesian network learning using domain/expert knowledge.  International Journal of Approximate Reasoning, 54(8):1168–1181, 2013.
  * [17] Kevin P Murphy.  Active learning of causal bayes net structure.  Technical report, technical report, UC Berkeley, 2001.
  * [18] Judea Pearl.  Causality.  Cambridge university press, 2009.
  * [19] Jonathan Richens and Tom Everitt.  Robust agents learn causal world models.  arXiv preprint arXiv:2402.10877, 2024.
  * [20] F. D. Schönbrodt and E.J. Wagenmakers.  Bayes factor design analysis: Planning for compelling evidence.  Psychon. B. Rev., 25(1):128–142, 2017.
  * [21] Peter Spirtes, Clark Glymour, and Richard Scheines.  Causation, prediction, and search.  MIT press, 2001.
  * [22] Panagiotis Tigas, Yashas Annadani, Andrew Jesson, Bernhard Schölkopf, Yarin Gal, and Stefan Bauer.  Interventions, where and how? experimental design for causal models at scale.  Advances in Neural Information Processing Systems, 35:24130–24143, 2022.
  * [23] Simon Tong and Daphne Koller.  Active learning for structure in bayesian networks.  In International joint conference on artificial intelligence, volume 17, pages 863–869. Citeseer, 2001.
  * [24] Rens van de Schoot, Sarah Depaoli, Ruth King, Bianca Kramer, Kaspar Märtens, Mahlet G Tadesse, Marina Vannucci, Andrew Gelman, Duco Veen, Joukje Willemsen, et al.  Bayesian statistics and modelling.  Nature Reviews Methods Primers, 1(1):1, 2021.

##  Appendix A Optimization objective and the problems computation derivation

###  A.1 Derivation for equation 8

Following the definition of $P_{DC}$ in eq 6, we have:

$\displaystyle P_{DC}$
$\displaystyle(\mathbf{D}_{\text{int}},\mathrm{do}(X=x))=P_{DC}^{0}(\mathbf{D}_{\text{int}},\mathrm{do}(X=x))P(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}})+P_{DC}^{1}(\mathbf{D}_{\text{int}},\mathrm{do}(X=x))P(\mathbb{H}_{1}\mid\mathbf{D}_{\text{int}})$
Expanding $P_{DC}^{0}$ and $P_{DC}^{1}$ using their definitions in eq 5
$\displaystyle=P\left(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\mathbf{D}_{\text{new}})>k_{0}\mid\mathbf{D}_{\text{int}},\mathrm{do}(X=x),\mathbb{H}_{0}\right)P(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}})$
$\displaystyle\quad+P\left(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\mathbf{D}_{\text{new}})<k_{1}\mid\mathbf{D}_{\text{int}},\mathrm{do}(X=x),\mathbb{H}_{1}\right)P(\mathbb{H}_{1}\mid\mathbf{D}_{\text{int}})$
Using the expectation form of the probability of an event $A$,
$P(A)=\mathbb{E}[\mathbb{I}(A)]$ $\displaystyle=\mathbb{E}_{y\sim
m_{0}(y)}\left[\mathbb{I}\left(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\\{y,x\\})>k_{0}\right)\mid\mathbf{D}_{\text{int}},\mathrm{do}(X=x),\mathbb{H}_{0}\right]P(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}})$
$\displaystyle\quad+\mathbb{E}_{y\sim m_{1}(y\mid
x)}\left[\mathbb{I}\left(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\\{y,x\\})<k_{1}\right)\mid\mathbf{D}_{\text{int}},\mathrm{do}(X=x),\mathbb{H}_{1}\right]P(\mathbb{H}_{1}\mid\mathbf{D}_{\text{int}})$

Calculating $P_{DC}$ for general data likelihoods can be extremely challenging
due to the complexity and high dimensionality of the data. To address this
difficulty, we adopt a Monte Carlo estimation approach combined with
optimization techniques. For this approach to be effective, it is crucial that
the terms within the Monte Carlo estimator are differentiable. Therefore, we
need to handle the indicator functions within the expectations. The indicator
function $\mathbb{I}(A)$ is equal to 1 if the condition $A$ is true and 0
otherwise. For our purposes, we transform the indicator function into a
Heaviside step function, $H(x)$, which is defined as:

$H(x)=\begin{cases}0&\text{if }x<0\\\ 1&\text{if }x\geq 0\end{cases}$

Using this, we can rewrite $\mathbb{I}(A)$ as $H(x-k)$ for some threshold $k$.
However, the Heaviside function is not differentiable, and its "derivative" is
the Dirac delta function, which is not a function in $L_{1}$ space but a
distribution i.e. genenralized function. To facilitate Monte Carlo estimation
and optimization, we need a smooth, differentiable approximation of the
Heaviside function.

We approximate the Heaviside function with a smooth function that transitions
smoothly around the threshold. One common approach is to use the following
exponential smoothing:

$H_{\beta}(x)=\begin{cases}\exp\left(-\frac{x}{\beta}\right)&\text{if }x<0\\\
1=\exp\left(-\frac{0}{\beta}\right)&\text{if }x\geq 0\end{cases}$

As $\beta$ approaches 0, $H_{\beta}(x)$ converges to the Heaviside function
$H(x)$. By using this smooth approximation, we can maintain differentiability.
Noting that $H_{\beta}$ is in the form of
$\exp\left(-\frac{1}{\beta}(\cdot)\right)$, we can rewrite it as
$H_{\beta}(x)=\exp\left(-\frac{1}{\beta}\text{ReLU}(x)\right)$.

With this approximation, we can rewrite the expression. Taking the first term
in Eq. A.1 as an example:

$\displaystyle\mathbb{E}_{y\sim
m_{0}(y)}\left[\mathbb{I}\left(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\\{y,x\\})>k_{0}\mid\mathbf{D}_{\text{int}},\mathrm{do}(X=x),\mathbb{H}_{0}\right)\right]P(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}})$
$\displaystyle=\mathbb{E}_{y\sim
m_{0}(y)}\left[H\left(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\\{y,x\\})-k_{0}\right)\mid\mathbf{D}_{\text{int}},\mathrm{do}(X=x),\mathbb{H}_{0}\right]P(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}})$
Take a small enought $\beta$ $\displaystyle\approx\mathbb{E}_{y\sim
m_{0}(y)}\left[H_{\beta}\left(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\\{y,x\\})-k_{0}\right)\mid\mathbf{D}_{\text{int}},\mathrm{do}(X=x),\mathbb{H}_{0}\right]P(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}})$
$\displaystyle=\mathbb{E}_{y\sim
m_{0}(y)}\left[\exp\left(-\frac{1}{\beta}\text{ReLU}\left(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\\{y,x\\})-k_{0}\right)\right)\mid\mathbf{D}_{\text{int}},\mathrm{do}(X=x),\mathbb{H}_{0}\right]P(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}})$

Thus, we replace the indicator function with this smooth approximation in our
Monte Carlo estimator. The expectation can then be approximated using Monte
Carlo sampling:

$\displaystyle P_{DC}(\mathbf{D}_{\text{int}},\mathrm{do}(X=x))$
$\displaystyle\approx\frac{1}{N}\sum_{i=1}^{N}\exp\left(-\frac{1}{\beta}\text{ReLU}(k_{0}-\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\\{y_{0i},x\\}))\right)P(\mathbb{H}_{0}\mid\mathbf{D}_{\text{int}})$
$\displaystyle+\frac{1}{N}\sum_{i=1}^{N}\exp\left(-\frac{1}{\beta}\text{ReLU}(\text{BF}_{01}(\mathbf{D}_{\text{int}}\cup\\{y_{1i},x\\})-k_{1})\right)P(\mathbb{H}_{1}\mid\mathbf{D}_{\text{int}})$

In our experiments, we set $N=4096$ and $\beta=0.2$. This choice of parameters
ensures a balance between approximation accuracy and computational efficiency.

In this form, the Monte Carlo estimator is differentiable, allowing us to
efficiently optimize $P_{DC}$ by selecting the intervention $x$ that maximizes
the probability of obtaining decisive and correct evidence.

###  A.2 Computation of Information Gain

In this section, we detail the computation of information gain for causal
discovery. Our goal is to estimate the information gain conditioned on the
existing interventional data $\mathbf{D}_{\text{int}}$ when considering the
new data $\mathbf{D}_{\text{new}}$ obtained from the intervention
$\mathrm{do}(X=x)$. Specifically, we treat the existence of a direct causal
link from $X$ to $Y$ (the null hypothesis or the alternative hypothesis) as a
random variable, and calculate the mutual information between this variable
and the union of the pre-existing interventional samples with the newly
obtained data.

When using information gain as a utility function for causal discovery, we aim
to maximize the value of the intervention point to achieve the highest
possible information gain. According to the conditional information equation
[4], we have:

$\displaystyle
I(\mathbf{D}_{\text{new}};H\mid\mathbf{D}_{\text{int}},\mathrm{do}(X=x))=I(\mathbf{D}_{\text{int}}\cup\mathbf{D}_{\text{new}};H\mid\mathrm{do}(X=x))-I(\mathbf{D}_{\text{int}};H\mid\mathrm{do}(X=x))$

Let $\mathbf{D}=\mathbf{D}_{\text{int}}\cup\mathbf{D}_{\text{new}}$. Note that
$\mathbf{D}_{\text{int}}$ is already a constant and does not contain any
randomness, thus $\mathbf{D}\mid\mathrm{do}(X=x),H\sim p(y\mid x,H)$. The
second term in Eq. (A.2) is a constant since $\mathbf{D}_{\text{int}}$ is
specified, and we can therefore neglect it.

$\displaystyle
I(\mathbf{D}_{\text{int}}\cup\mathbf{D}_{\text{new}};H\mid\mathrm{do}(X=x))$
$\displaystyle=\mathbb{E}_{\mathbf{D}\sim p(y\mid H,\mathrm{do}(X=x)),H\sim
P(H)}\log\frac{P(\mathbf{D},H\mid\mathrm{do}(X=x))}{P(\mathbf{D}\mid\mathrm{do}(X=x))P(H\mid\mathrm{do}(X=x))}$
$\displaystyle=\mathbb{E}_{\mathbf{D}\sim p(y\mid H,\mathrm{do}(X=x)),H\sim
P(H)}\log\frac{P(\mathbf{D}\mid
H,\mathbf{D}_{\text{int}},\mathrm{do}(X=x))}{P(\mathbf{D}\mid\mathrm{do}(X=x))}$
$\displaystyle=\mathbb{E}_{\mathbf{D}\sim p(y\mid H,\mathrm{do}(X=x)),H\sim
P(H)}\log\frac{P(\mathbf{D}\mid
H,\mathrm{do}(X=x))}{\sum_{j\in\\{0,1\\}}P(\mathbf{D}\mid
H,\mathrm{do}(X=x))P(H=\mathbb{H}_{j}\mid\mathrm{do}(X=x))}$
$\displaystyle=\mathbb{E}_{\mathbf{D}\sim p(y\mid H,\mathrm{do}(X=x)),H\sim
P(H)}\log\frac{P(\mathbf{D}\mid
H,\mathrm{do}(X=x))}{\sum_{j\in\\{0,1\\}}P(\mathbf{D}\mid
H,\mathrm{do}(X=x))P(H=\mathbb{H}_{j})}$ $\displaystyle=\int p(y\mid
H=\mathbb{H}_{0},\mathrm{do}(X=x))P(H=\mathbb{H}_{0})\times$
$\displaystyle\log\frac{P(\mathbf{D}\mid
H=\mathbb{H}_{0},\mathrm{do}(X=x))}{P(\mathbf{D}\mid
H=\mathbb{H}_{0},\mathrm{do}(X=x))P(H=\mathbb{H}_{0})+P(\mathbf{D}\mid
H=\mathbb{H}_{1},\mathrm{do}(X=x))P(H=\mathbb{H}_{1})}\mathrm{d}y$
$\displaystyle+\int p(y\mid
H=\mathbb{H}_{1},\mathrm{do}(X=x))P(H=\mathbb{H}_{1})\times$
$\displaystyle\log\frac{P(\mathbf{D}\mid
H=\mathbb{H}_{1},\mathrm{do}(X=x))}{P(\mathbf{D}\mid
H=\mathbb{H}_{0},\mathrm{do}(X=x))P(H=\mathbb{H}_{0})+P(\mathbf{D}\mid
H=\mathbb{H}_{1},\mathrm{do}(X=x))P(H=\mathbb{H}_{1})}\mathrm{d}y$
$\displaystyle=\int p(y\mid
H=\mathbb{H}_{0},\mathrm{do}(X=x))P(H=\mathbb{H}_{0})\log\frac{\text{BF}_{01}(\mathbf{D})}{\text{BF}_{01}(\mathbf{D})P(H=\mathbb{H}_{0})+P(H=\mathbb{H}_{1})}\mathrm{d}y$
$\displaystyle+\int p(y\mid
H=\mathbb{H}_{1},\mathrm{do}(X=x))P(H=\mathbb{H}_{1})\log\frac{1}{\text{BF}_{01}(\mathbf{D})P(H=\mathbb{H}_{0})+P(H=\mathbb{H}_{1})}\mathrm{d}y$

Given this setup, we recognize that this formulation can be expressed in terms
of the Bayes factor. Therefore, we can leverage the previous calculation
methods for Bayes factors (Section 3.2) and use a Monte Carlo approach similar
to that used for estimating and optimizing $P_{DC}$ to estimate and optimize
this objective as well.

##  Appendix B Experiment results for $k_{0}=\frac{1}{k_{1}}=30$ and
$k_{0}=\frac{1}{k_{1}}=100$

In this section, we present detailed results for the experiment setting. The
following results complement the findings discussed in the main text and
provide further insights into the performance of our method across varying
conditions.

We include results for different confounders and robustness tests with random
means, variances, and weights. The method for generating the weights is
consistent with the approach described in the main text, where weights are
normalized using a combination of uniform and softmax distributions.

For the robustness tests, The means ($\mu_{i}$) of the components are sampled
from a uniform distribution in the range $[-4,4]$. The variances
($\sigma_{i}^{2}$) are sampled from a standard chi-squared distribution with 3
degrees of freedom.

These additional experiments ensure that our method’s effectiveness and
conclusions remain consistent even under varying noise characteristics and
distribution parameters. All experiments are repeated with 10 different random
seeds to ensure the robustness and reproducibility of our results.

Figure 3: Results under different ground truths with
$k_{0}=\frac{1}{k_{1}}=10$: $P_{DC}$, $\log\text{BF}_{01}$, and
$P(\mathbb{H}_{gt}\mid\mathbf{D}_{\text{int}})$. The first row corresponds to
$\mathbb{H}_{0}$ ($X\leftarrow Y$)). The second row corresponds to
$\mathbb{H}_{0}$ ($X\leftarrow U\to Y$), and the last row corresponds to
$\mathbb{H}_{1}$ ($X\to Y$).

Figure 4: Results under different ground truths with
$k_{0}=\frac{1}{k_{1}}=10$: $P_{DC}$, $\log\text{BF}_{01}$, and
$P(\mathbb{H}_{gt}\mid\mathbf{D}_{\text{int}})$. But with random mean and
random covariance compared to the settings in the Figure 3. The first row
corresponds to $\mathbb{H}_{0}$ ($X\leftarrow Y$)). The second row corresponds
to $\mathbb{H}_{0}$ ($X\leftarrow U\to Y$), and the last row corresponds to
$\mathbb{H}_{1}$ ($X\to Y$).

Figure 5: Results under different ground truths with
$k_{0}=\frac{1}{k_{1}}=10$: $P_{DC}$, $\log\text{BF}_{01}$, and
$P(\mathbb{H}_{gt}\mid\mathbf{D}_{\text{int}})$. The first row corresponds to
$\mathbb{H}_{0}$ ($X\leftarrow Y$)). The second row corresponds to
$\mathbb{H}_{0}$ ($X\leftarrow U\to Y$), and the last row corresponds to
$\mathbb{H}_{1}$ ($X\to Y$).

Figure 6: Results under different ground truths with
$k_{0}=\frac{1}{k_{1}}=10$: $P_{DC}$, $\log\text{BF}_{01}$, and
$P(\mathbb{H}_{gt}\mid\mathbf{D}_{\text{int}})$. But with random mean and
random covariance compared to the settings in the Figure 5. The first row
corresponds to $\mathbb{H}_{0}$ ($X\leftarrow Y$)). The second row corresponds
to $\mathbb{H}_{0}$ ($X\leftarrow U\to Y$), and the last row corresponds to
$\mathbb{H}_{1}$ ($X\to Y$).