Deriving Causal Order from Single-Variable Interventions: Guarantees & Algorithm

## Abstract

Abstract Targeted and uniform interventions to a system are crucial for unveiling causal relationships. While several methods have been developed to leverage interventional data for causal structure learning, their practical application in real-world scenarios often remains challenging. Recent benchmark studies have highlighted these difficulties, even when large numbers of single-variable intervention samples are available. In this work, we demonstrate, both theoretically and empirically, that such datasets contain a wealth of causal information that can be effectively extracted under realistic assumptions about the data distribution. More specifically, we introduce the notion of interventional faithfulness , which relies on comparisons between the marginal distributions of each variable across observational and interventional settings, and we introduce a score on causal orders . Under this assumption, we are able to prove strong theoretical guarantees on the optimum of our score that also hold for large-scale settings. To empirically verify our theory, we introduce Intersort , an algorithm designed to infer the causal order from datasets containing large numbers of single-variable interventions by approximately optimizing our score. Intersort outperforms baselines (GIES, DCDI, PC and EASE) on almost all simulated data settings replicating common benchmarks in the field. Our proposed novel approach to modeling interventional datasets thus offers a promising avenue for advancing causal inference, highlighting significant potential for further enhancements under realistic assumptions.

## 1 Introduction

Causal structure learning is pivotal for understanding complex systems, aiming to discover causal relationships from data. This field spans disciplines such as biology (Meinshausen et al., ; Yu et al., ), medicine (Farmer et al., ; Joffe et al., ; Feuerriegel et al., ), and social sciences Baum-Snow & Ferreira (); Imbens & Rubin (), where causal insights drive informed decision-making and deepen our understanding of underlying processes. Traditionally, causal discovery has relied heavily on observational data, primarily because expansive interventional experiments are often impractical or costly in many application domains. The inherent limitations of observational data necessitate assumptions about data distribution to ensure identifiability beyond the Markov equivalence class (Spirtes et al., ; Shimizu et al., ; Hoyer et al., ). However, the emergence of large-scale interventional data, especially in single-cell transcriptomics data (Replogle et al., ; Datlinger et al., , ; Dixit et al., ), introduces both new opportunities and challenges.

Interventional data, characterized by targeted alterations to a system, offers a unique perspective for causal discovery. These controlled interventions reveal causal mechanisms often obscured in observational studies. As such, many causal discovery methods using interventional data have been proposed (Brouillard et al., ; Lorch et al., ; Lopez et al., ; Hauser & Bühlmann, ). However, their practicality remains challenging. For example, recent benchmarks in the gene network domain (Chevalley et al., , ) demonstrate that simpler models based on straightforward statistical comparisons between observational and interventional distribution outperforms established causal discovery methods, particularly when many single-variable interventions are available. The motivation for this research thus stemmed from these recent advances in applying causal discovery techniques to gene network inference. Our methodology, including the use of the Wasserstein distance, aligns with evaluation metrics for predicted graphs as established in the benchmark by Chevalley et al. (). Additionally, a component of our algorithm builds upon and formalizes approaches from leading models in gene network inference, such as those highlighted in the CausalBench challenge (Kowiel et al., ; Nazaret & Hong, ; Deng & Guan, ; Chevalley et al., ). This work not only addresses a critical domain of application but also enhances the theoretical framework surrounding contemporary advancements in the field (Yu et al., ; Chai et al., ; Akers & Murali, ; Hu et al., ).

Here, we study further the potential of single-variable interventional data for causal structure discovery, which we reduce to finding the *causal order* (topological order) of the variables. Even though the causal order represents a superset of the true causal relationships, it is useful *per se*. Intuitively and informally, correlations are more likely to be causations if they follow the causal order. More practically, it can for example be used to guide experimental design, as it greatly reduces the hypothesis space. As a specific example, in biology, a task of interest is to select pairs of genes to intervene on, which is a space much larger than single gene interventions. For a target gene that let us assume lies in the middle of the causal order, the number of candidate double interventions goes from ${d\choose 2}$ down to ${\frac{d}{2}\choose 2}$, where $d$ is the number of variables.

To find the causal order, we introduce a data-based score to rank candidate causal orders, and we establish theoretical guarantees on its optimum, providing upper-bounds on the expected error depending on the graph density and on the probability for a variable to be intervened on, that also hold for large scale settings. Crucially, our approach makes light assumptions on the data distributions, requiring only what we call *$\epsilon$-interventional faithfulness*, which characterizes the strength of changes in marginal distributions between observational and interventional distributions. We then introduce a novel algorithm called Intersort that leverages this assumption to derive the causal order.

Our empirical evaluations on diverse simulated datasets (linear, random Fourier features (Lorch et al., ), neural network (Nazaret et al., ; Brouillard et al., ) and single cell (Dibaeinia & Sinha, ) with various noise distributions), which are close replications of real-world systems and follow common practice in the field, confirm our theoretical results. The efficacy of our approach in accurately determining the causal order from interventional data signifies a notable advancement in causal inference. Intersort outperforms the three baselines, namely PC (Spirtes et al., ), GIES (Hauser & Bühlmann, ), DCDI (Brouillard et al., ) and EASE (Gnecco et al., ) on almost all settings. It is also robust to data normalization (Reisach et al., ), which is a failure mode of many recently proposed continuous optimization causal discovery methods. Moreover, our empirical results show that Intersort makes more efficient use of interventional information compared to existing approaches which is a critical advantage especially in domains where interventions are costly and sometimes impossible to perform. This suggests that the task of identifying causal structures with many single-variable interventions may be more feasible than previously believed. Moreover, our method demonstrates similar performance and robustness across different types of simulated data, emphasizing its reliability and versatility.

Our method relies on *$\epsilon$-interventional faithfulness* which is a lighter version of many of the existing assumptions in causal discovery literature, while still providing theoretical guarantees. Realistic assumptions on the data distributions are crucial, as in many application domains such as gene expression data, the validity of common assumptions in the field may be unverifiable or known to not hold. Furthermore, many methods catastrophically fail when their assumptions do not hold, rendering them inapplicable (Montagna et al., ; Heinze-Deml et al., ).

## 2 Related work

##### Causal ordering

Even though the causal order does not contain the full causal information as it does not uniquely identify the causal graph, it can subsequently be used with for example penalized regression techniques to recover the edges (Bühlmann et al., ; Shimizu et al., ). Also, the correct causal order is useful in itself as a fully connected graph can be constructed from it which describes the interventional distributions (Peters & Bühlmann, ; Bühlmann et al., ). More recently, following the discovery of Reisach et al. () that the causal order can be recovered by sorting the variable by variance in simulated datasets, many algorithms to recover the causal order from observational data have been proposed, for example via score matching (Rolland et al., ; Montagna et al., , ). To our knowledge, we are the first to propose an algorithm to infer the causal order from interventional data.

## 3 Causal Order: Definitions

We here introduce our framework for causality following Pearl () and we follow the notations of Peters et al. (). We also present the definition of a causal order as well as a divergence to measure how faithfully a causal order recalls the causal information.

##### Statistical metric

Let $(M,d)$ be a metric space, and let $\mathcal{P}(M)$ be the set of probability measures over $M$. We define $D:\mathcal{P}(M)\times\mathcal{P}(M)\rightarrow[0,\infty)$ as a distance between probability measures on $M$.

##### Distribution

We consider a set of $d$ random variables $X=(X_{1},...,X_{d})$ indexed by $V=\{1,...d\}$, with associated joint distribution $P_{X}$. We denote the marginal distribution of each random variable as $P_{X_{i}}$, $i\in V$.

##### Causal Graph

We write $\mathcal{G}=(V,E)$ a DAG between variables $v\in V$, $E\subseteq V^{2}$, such that $(v,v)\notin E$ for all $v\in V$. We write $A^{\mathcal{G}}$ the $d\times d$ adjacency matrix, where $A^{\mathcal{G}}(i,j)=1\iff(i,j)\in E$. We note $\textbf{Pa}_{j}$ the parents of $j$, where $i$ is a parent if $(i,j)\in E$. We note $\textbf{DE}_{\mathcal{G}}^{i}$ the descendants of $i$ in graph $\mathcal{G}$. $j$ is a descendant of $i$ if there is a directed path from $i$ to $j$. We note $\textbf{AN}_{\mathcal{G}}^{i}$ the ancestors of $i$ in graph $\mathcal{G}$. $j$ is an ancestor of $i$ if there is a directed path from $j$ to $i$.

##### Causal Order

Given a causal graph $\mathcal{G}$, we call a permutation of the variables:

$$ $\pi:\{1,...d\}\rightarrow\{1,...d\},$ $$

a causal ordering of the variables in $\mathcal{G}$ if it satisfies $\pi(i)<\pi(j)$ if $j\in\textbf{DE}_{\mathcal{G}}^{i}$ (Peters et al., ). A permutation is a bijective mapping of the indices to new indices.

For a causal graph, a causal ordering always exists, but may not be unique. We note the set of causal orderings $\Pi^{*}$, and any member of $\Pi^{*}$ as $\pi^{*}$. We can rearrange the variables according to a permutation $\pi$, and we write the new associate adjacency matrix $A^{\mathcal{G}}_{\pi}$, where $A^{\mathcal{G}}_{\pi}$ is upper triangular if $\pi\in\Pi^{*}$. Lastly, we define $\overline{A^{\mathcal{G}}_{\pi}}=A^{\mathcal{G}}_{\pi}+{A^{\mathcal{G}}_{\pi}}
^{2}+...+{A^{\mathcal{G}}_{\pi}}^{d-1}$ the transitive closure of the adjacency matrix. It is straightforward to show that $\overline{A^{\mathcal{G}}_{\pi^{*}}}(i,j)>0$ if and only if there is a directed path from $i$ to $j$ in $\mathcal{G}$. Furthermore, there may be a total causal effect(^1^11There is a total causal effect of $i$ on $j$ if there exists some $\tilde{N}$ such that $X_{i}\not\!\perp\!\!\!\perp X_{j}$ in $P_{X}^{C,do(X_{k}:=\tilde{N}))}$ (Definition 6.12 in Peters et al. ()).) of $i$ on $j$ only if $\overline{A^{\mathcal{G}}_{\pi^{*}}}(i,j)>0$. $\overline{A^{\mathcal{G}}_{\pi^{*}}}$ is also upper-triangular.

##### Top Order Divergence D t ⁢ o ⁢ p subscript 𝐷 𝑡 𝑜 𝑝 {D_{top}} italic_D start_POSTSUBSCRIPT italic_t italic_o italic_p end_POSTSUBSCRIPT

To measure the discrepancy between a permutation $\pi$ and a graph $\mathcal{G}$, we use the top order divergence (Rolland et al., ). The top order divergence $D_{top}:(\{1,...d\},\{1,...d\}\times\{1,...d\})\times(\{1,...d\}\rightarrow\{1
,...d\})\rightarrow\mathbb{N}_{0}$ is defined as:

$$ $D_{top}(\mathcal{G},\pi)=\sum_{\pi(i)>\pi(j)}A^{\mathcal{G}}(i,j).$ $$

$D_{top}$ counts the number of edges that cannot be recovered given a topological order $\pi$, or said differently, the number of false negative edges of the fully connected DAG corresponding to $\pi$. As such, it represents a lower bound on the structural hamming distance (SHD) achievable when constrained to a particular topological order.

###### Observation 1 .

The divergence $D_{top}$ can also be written as $D_{top}(\mathcal{G},\pi)=\sum_{i<j}A^{\mathcal{G}}_{\pi}(i,j)$.

###### Lemma 1 .

Let $\pi^{*}\in\Pi^{*}$ be a causal ordering for $\mathcal{G}$. Then $D_{top}(\mathcal{G},\pi^{*})=0$. (Rolland et al., )

###### Observation 1 .

###### Lemma 1 .

##### SCMs and interventions

A structural causal model (SCM) $\mathcal{C}=(\textbf{S},P_{N})$ is determined by a collection of $d$ assignments $s\in\textbf{S}$

$$ $X_{j}:=f_{j}(X_{pa_{j}},N_{j}),$ $$

where $\textbf{Pa}_{j}\in V\setminus j$ and $P_{N}$ is a joint distribution over the noise variables that is jointly independent (Peters et al., ).

##### Distribution with Intervention P X C , d o ( X k := N ~ k ) ) P_{X}^{C,do(X_{k}:=\tilde{N}_{k}))} italic_P start_POSTSUBSCRIPT italic_X end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_C , italic_d italic_o ( italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT := over~ start_ARG italic_N end_ARG start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ) ) end_POSTSUPERSCRIPT

An intervention is then defined as a replacement of a subset of the structural assignments of an SCM $\mathcal{C}$, such that no cycle is created. In this work, we consider intervention on a single variable, where the structural assignment is replaced by a new random variable independent of the parents. We denote the new distribution entailed by an intervention $P_{X}^{\mathcal{C},do(X_{k}:=\tilde{N}_{k})}$. We denote $\mathcal{I}\subseteq\{1,...d\}$ the set of nodes that are intervened on. We also assume that for each intervention, we have access to the interventional distribution $P_{X}^{\mathcal{C},do(X_{k}:=\tilde{N}_{k})},k\in\mathcal{I}$. The set of interventional random variables is denoted as $\tilde{N}=\bigcup_{k\in\mathcal{I}}\{\tilde{N}_{k}\}$.

## 4 ϵ italic-ϵ \epsilon italic_ϵ -Interventional faithfulness

We here introduce the main assumption of our work that we name *interventional faithfulness*. Akin to the common assumption of faithfulness of a joint distribution with respect to the true graph, which assumes that all d-separations in the graph are reflected by conditional independences in the entailed distribution, $\epsilon$-interventional faithfulness assumes that all paths in the graphs are revealed by changes in distribution under intervention above a significance threshold $\epsilon$. This assumption is also similar to how Mooij et al. () define causality as the existence of interventions inducing a change in distribution. The main difference with comparable assumptions made in the literature is that this definition depends on a statistical distance and significance threshold, instead of an oracle of the statistical dependence. We formalize our assumption as follows:

###### Definition 1 .

Given the distributions $P_{X}^{\mathcal{C},(\emptyset)}$ and $P_{X}^{\mathcal{C},do(X_{k}:=\tilde{N}_{k})},\forall k\in\mathcal{I}$, we say that the tuple $(\tilde{N},\mathcal{C})$ is $\epsilon$*-interventionally faithful* to the graph $\mathcal{G}$ associated to $\mathcal{C}$ if for all $i\neq j,i\in\mathcal{I},j\in V$, $D\left(P_{X_{j}}^{\mathcal{C},(\emptyset)},P_{X_{j}}^{\mathcal{C},do(X_{i}:=
\tilde{N}_{i})}\right)>\epsilon$ if and only if there is a directed path from $i$ to $j$ in $\mathcal{G}$.

The concept of interventional faithfulness, where interventions on a node cause significant deviations in the distribution of downstream nodes, is discussed in various contexts. It aligns closely with the notion of ”c-faithfulness” in causal graphs, where interventional distributions respect the causal structure and result in significant changes in downstream nodes (Shanmugam et al., ; Squires et al., ). Assumption 4.4 in Yang et al. () explicitly states that interventions on upstream nodes should affect downstream nodes. The causal Markov condition and faithfulness for both observational and interventional distributions also support this idea (Addanki & Kasiviswanathan, ). The use of interventional data to improve the identifiability of causal models is another related concept, emphasizing how interventions can reveal causal relationships not identifiable from observational data alone (Hauser & Bühlmann, ).

###### Lemma 2 .

If all interventions in $\tilde{N}$ have full support and all directed paths in $\mathcal{G}$ imply a total causal effect, then $(\tilde{N},\mathcal{C})$ is $\epsilon$-interventionally faithful for $\epsilon=0$. (Follows from Proposition 6.13 in Peters et al. ().)

###### Lemma 3 .

Consider a linear structural causal model $\mathcal{C}$ defined by the equations $X_{j}=\sum_{i\in\text{Pa}_{j}}\beta_{ij}X_{i}+N_{j}$, where $\beta_{ij}$ are the weights representing causal effects, drawn independently from a continuous distribution, and $N_{j}$ are noise variables that are independent of each other and of the parent variables, each with a distribution having full support. Let the $\tilde{N}_{i}$ be distributions with full support. Then, $(\tilde{N},\mathcal{C})$ is almost surely $\epsilon$-interventionally faithful for $\epsilon=0$, as the total causal effect along any directed path is non-zero with probability one.

###### Remark .

If $(\tilde{N},\mathcal{C})$ is $\epsilon$-interventionally faithful for $\epsilon=0$, then it is also $\epsilon$-interventionally faithful for all $0<\epsilon<\min\{D\left(P_{X_{j}}^{\mathcal{C},(\emptyset)},P_{X_{j}}^{
\mathcal{C},do(X_{i}:=\tilde{N}_{i})}\right),\forall(i,j)\in E:D\left(P_{X_{j}
}^{\mathcal{C},(\emptyset)},P_{X_{j}}^{\mathcal{C},do(X_{i}:=\tilde{N}_{i})}
\right)>0\}$

Lemma 3 shows that the $\epsilon$-interventionally faithful assumption covers a large class of SCMs. We leave further theoretical work to characterize how large the class of $\epsilon$-interventionally faithful $(\tilde{N},\mathcal{C})$ tuples is as future work.

###### Definition 1 .

###### Lemma 2 .

###### Lemma 3 .

###### Remark .

## 5 A score on causal orders

Given an observational distribution $P_{X}^{\mathcal{C},(\emptyset)}$ and a set of interventional distributions $\mathcal{P}_{int}=\{P_{X}^{\mathcal{C},do(X_{k}:=\tilde{N}_{k})},k\in\mathcal{
I}\}$, $\mathcal{I}\subseteq V$, we define the following score for a permutation $\pi$, for some statistical distance $D:\mathcal{P}(M)\times\mathcal{P}(M)\rightarrow[0,\infty)$, $\epsilon>0$, $c>0$:

$$ $\begin{split}S(\pi,\epsilon,D,\mathcal{I},P_{X}^{\mathcal{C},(\emptyset)}, \mathcal{P}_{int},c)=&\sum_{\pi(i)<\pi(j),i\in\mathcal{I},j\in V}\left(D\left( P_{X_{j}}^{\mathcal{C},(\emptyset)},P_{X_{j}}^{\mathcal{C},do(X_{i}:=\tilde{N} _{i})}\right)-\epsilon\right)\\ &+c\cdot d\cdot\bm{1}_{D\left(P_{X_{j}}^{\mathcal{C},(\emptyset)},P_{X_{j}}^{ \mathcal{C},do(X_{i}:=\tilde{N}_{i})}\right)>\epsilon}\end{split}$ (1) $$

Usually, we can set the $\epsilon$ of eq. 1 to the same value of $\epsilon$-interventional faithfulness as the input dataset, unless it is equal to zero. However, the remark in remark Remark ensures the existence of a suitable $\epsilon$ in that case, as it can be chosen in a data-driven way by taking a non-zero value for $\epsilon$ smaller than the smallest non-negative distance.

###### Theorem 1 .

Assume that we are given $P_{X}^{\mathcal{C},(\emptyset)}$ and $P_{X}^{\mathcal{C},do(X_{k}:=\tilde{N}_{k})},\forall k\in\mathcal{I}$ and $\mathcal{I}=V$, such that $(\tilde{N},\mathcal{C})$ is $\epsilon$-interventionally faithful for some $\epsilon>0$, and let $\pi_{opt}\in\operatorname*{arg\,max}_{\pi}\mathcal{S}(\pi)$. Then $D_{top}(\mathcal{G},\pi_{opt})=D_{top}(\mathcal{G},\pi^{*})=0$.

The above theorem shows that in the case where we observe interventions on all the variables, we are guaranteed to find a valid topological order by maximizing our score when the data is $\epsilon$-interventionally faithful. Similar results were proven for this case where all the variables are intervened (Hauser & Bühlmann, ; Tian & Pearl, ). We now turn to the case where variables have a probability of being intervened smaller than one. In that case, the error will not be zero in expectation, and we thus aim to upper-bound the expected error as measured by $D_{top}$. We first prove an important property of the optimum to our proposed score. Lemma 4 shows that for a given edge, a single intervention among a large set of candidate variables is sufficient to correctly order the two variables of an edge.

###### Lemma 4 .

Assume that we are given $P_{X}^{\mathcal{C},(\emptyset)}$ and $P_{X}^{\mathcal{C},do(X_{k}:=\tilde{N}_{k})},\forall k\in\mathcal{I}$, such that $(\tilde{N},\mathcal{C})$ is $\epsilon$-interventionally faithful for some $\epsilon>0$, and let $\pi_{opt}\in\operatorname*{arg\,max}_{\pi}\mathcal{S}(\pi)$. Let $(i,j)\in E$, then if $j\in\mathcal{I}$ or for some $k\in\textbf{AN}_{j}^{\mathcal{G}}\setminus\textbf{AN}_{i}^{\mathcal{G}},k\in
\mathcal{I}$, then $\pi_{opt}(i)<\pi_{opt}(j)$.

We now use this result to derive a graph dependent upper-bound to the expected error, given a uniform probability for a variable to be intervened. We leave theoretical derivations for other intervention distributions, for example data dependent intervention selection, as future work.

###### Theorem 2 .

Assume that we are given $P_{X}^{\mathcal{C},(\emptyset)}$ and $P_{X}^{\mathcal{C},do(X_{k}:=\tilde{N}_{k})},\forall k\in\mathcal{I}$, such that $(\tilde{N},\mathcal{C})$ is $\epsilon$-interventionally faithful for some $\epsilon>0$, and let $\mathcal{I}$ be chosen uniformly at random, where $p_{int}:=P(i\in\mathcal{I})\forall i\in V,0<p_{int}<1$, then $\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]\leq\sum_{(i,j)\in\mathcal{G}}(1-p_{
int})^{|\textbf{AN}_{j}^{\mathcal{G}}\cup\{j\}\setminus\textbf{AN}_{i}^{
\mathcal{G}}|}$.

###### Theorem 1 .

###### Lemma 4 .

###### Theorem 2 .

## 6 Relaxation of ϵ italic-ϵ \epsilon italic_ϵ -Interventional Faithfulness

We now look at the case where $\epsilon$-interventional faithfulness may not be fully fulfilled and we derive some bounds for that setting. We consider the worst case to be the one where interventions only reveal direct causal relationships, i.e., only children are affected by a node intervention such that the distance is larger than $\epsilon$. We can hypothesize that real distributions lie on a continuum between those two extremes, that is, some interventions have an $\epsilon$-strong effect only on children and other interventions have an effect on more downstream variables. This restricted setting also covers the case when multiple paths between two variables $i$ and $j$ ”cancel out” the causal effect of $i$ on $j$. In that restricted setting, we can rewrite lemma 4 as follows:

###### Lemma 5 .

Let $(i,j)\in E$, then if $j\in\mathcal{I}$ or for some $k\in\textbf{Pa}_{j}^{\mathcal{G}}\setminus\textbf{Pa}_{i}^{\mathcal{G}},k\in
\mathcal{I}$, then $\pi_{opt}(i)<\pi_{opt}(j)$.

It is then straightforward that the bound on $D_{top}$ for $\pi_{opt}$ in that case is:

###### Theorem 3 .

Let $\mathcal{I}$ be chosen uniformly at random, where $p_{int}:=P(i\in\mathcal{I})\forall i\in V,0<p_{int}<1$, then $\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]\leq\sum_{(i,j)\in\mathcal{G}}(1-p_{
int})^{|\textbf{Pa}_{j}^{\mathcal{G}}\cup\{j\}\setminus\textbf{Pa}_{i}^{
\mathcal{G}}|}$.

Now that the bound only depends on the parent of the variables, it becomes possible to find a closed form for the expectation of $D_{top}$ for a particular graph distribution.

###### Theorem 4 .

Let $\mathcal{I}$ be chosen uniformly at random, $\mathcal{G}$ be a random Erdos-Renyi directed acyclic graph with edge probability $p_{e}$, where $p_{int}:=P(i\in\mathcal{I})\forall i\in V,0<p_{int}$, then $\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]\leq\frac{(1-p_{int})^{2}}{p_{int}}
\left[d-(1-p_{int}p_{e})\frac{1-(1-p_{int}p_{e})^{d}}{p_{int}p_{e}}\right]$. We also have a looser, but independent of $p_{e}$, bound $\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]\leq\frac{(1-p_{int})^{2}}{p_{int}}d$.

The natural next step is to look at how the bound behaves as $d$ goes to infinity. More precisely, let’s look at the normalized error $\frac{\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]}{d}$. If $p_{e}$ is fixed, we simply get $\lim_{d\rightarrow\infty}\frac{\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]}{d}
\leq\frac{(1-p_{int})^{2}}{p_{int}}$. However, in realistic scenarios, we can expect the edge probability to depend on the number of variables. A common assumption is that the expected number of edges per variable is constant.

###### Lemma 6 .

Let $\mathbb{E}|E|=\frac{d(d-1)}{2}p_{e}=\frac{c}{2}(d-1)$, where $c>0,c\in\mathbb{R}$ is a constant, then $\lim_{d\rightarrow\infty}\frac{\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]}{d}
\leq\frac{(1-p_{int})^{2}}{p_{int}}\left[1-\frac{1}{p_{int}\cdot c}\left(1-e^{
-p_{int}\cdot c}\right)\right]$.

The above lemma shows that the expected error as $d$ goes to infinity is $\mathcal{O}(d)$, providing a strong guarantee on $\pi_{opt}$ even in a large scale setting.

###### Lemma 5 .

###### Theorem 3 .

###### Theorem 4 .

###### Lemma 6 .

## 7 Intersort

Given that our score for causal order eq. 1 is optimized over a set of permutations, finding the optimum by exhaustively trying all possible solutions becomes intractable as the number of variables grows larger. We thus develop an approximation algorithm to tractably find a solution that optimizes for our proposed score in eq. 1. This algorithm consists of two steps. The first step uses a different score to find an initial solution. In the second step, this solution is refined by performing a local search using our original score eq. 1.

### 7.1 Intersort step 1: Finding an initial solution

To find an initial solution, we design a new score, similar to eq. 1, but this time on graphs. The intuition is that we want to score edges directly, and edges $(i,j)$ that have a distance $D\left(P_{X_{j}}^{\mathcal{C},(\emptyset)},P_{X_{j}}^{\mathcal{C},do(X_{i}:=
\tilde{N}_{i})}\right)$ greater than $\epsilon$ should be on the upper-triangular part. We thus write this new score $\mathcal{S}_{init}$ as follows:

$$ $\mathcal{S}_{init}\left(G,\epsilon,D,\mathcal{I},P_{X}^{\mathcal{C},(\emptyset )},\mathcal{P}_{int}\right)=\sum_{(i,j)\in G,i\in\mathcal{I}}D\left(P_{X_{j}}^ {\mathcal{C},(\emptyset)},P_{X_{j}}^{\mathcal{C},do(X_{i}:=\tilde{N}_{i})} \right)-\epsilon|G|$ (2) $$

Then, we aim to find a DAG $G$ that maximizes this score.

$$ $G_{opt}=\operatorname*{arg\,max}_{G\in\mathcal{G}_{DAG}}\mathcal{S}_{init} \left(G,\epsilon,D,\mathcal{I},P_{X}^{\mathcal{C},(\emptyset)},\mathcal{P}_{ int}\right)$ (3) $$

Once again, the search over all possible DAGs is computationally expensive. We thus design an approximation algorithm based on the following heuristic: the highest values $D\left(P_{X_{j}}^{\mathcal{C},(\emptyset)},P_{X_{j}}^{\mathcal{C},do(X_{i}:=
\tilde{N}_{i})}\right)$, which can be seen as a score on the corresponding edge $(i,j)$, should be edges in the solution DAG to maximize the score $\mathcal{S}_{init}$. To do so, we sort the edges based on this ”edge score” from highest to lowest. We then construct our solution by adding the edges following the sorted order, adding an edge only if it does not break the acyclity of the solution. We then stop adding edges when those have a score $D_{ij}=D\left(P_{X_{j}}^{\mathcal{C},(\emptyset)},P_{X_{j}}^{\mathcal{C},do(X_
{i}:=\tilde{N}_{i})}\right)$ lower than $\epsilon$ in the matrix $\mathbf{D}=(D_{ij})$. The topological order of this solution $G_{opt}$ for $\mathcal{S}_{init}$ is then used as initial condition for our original objective of maximizing $\mathcal{S}$. The runtime of this algorithm is dominated by the sorting step, which is $\mathcal{O}(d\cdot|\mathcal{I}|\log(d\cdot|\mathcal{I}|))$. We call this procedure sortranking and illustrate it in algorithm 2 in the appendix.

### 7.2 Intersort step 2: A local search for π o ⁢ p ⁢ t subscript 𝜋 𝑜 𝑝 𝑡 \pi_{opt} italic_π start_POSTSUBSCRIPT italic_o italic_p italic_t end_POSTSUBSCRIPT

Figure: Algorithm 1 Complete algorithm

Starting from the initial solution found by maximizing eq. 2, we then perform a greedy local search to optimize our original objective eq. 1. At each step, given a current solution, we search in a close neighbourhood set for a candidate permutation $\pi$ with a higher score than our current solution. To derive the neighborhood set, we introduce an operator $\Sigma:\{\{1,...d\}\rightarrow\{1,...d\}\}\rightarrow\{\{1,...d\}\rightarrow\{
1,...d\}\}$, defined such that for a given set of permutations $\Pi$, $\Sigma(\Pi)$ returns a new set of permutations $\Pi^{\prime}$. Each permutation $\pi\in\Pi$ is modified by applying $\sigma:\{1,...d\}\rightarrow\{1,...d\}\rightarrow\{\{1,...d\}\rightarrow\{1,..
.d\}\}$, which for each variable $i\in V$, returns all the permutations such that $\pi(i)$ takes any of the other $d-1$ possible values. $\sigma$ takes one permutation as input and return a set of permutations which is $\mathcal{O}(d^{2})$. The total number of applications of $\Sigma$, denoted as $k$, directly defines the neighborhood size which is $\mathcal{O}(d^{2k})$. This definition allows us to explore varying depths of the search space, which is $\mathcal{O}(d!)$, by controlling $k$, thereby balancing between the breadth of the search and the computational resources required. A smaller $k$ value implies a tighter search around the initial solution and encode our confidence in the initial solution, optimizing runtime at the potential cost of missing broader, potentially better solutions. In our experiments, we use $k=1$.
We then stop once an iteration of the algorithm brings no improvement to the score. This procedure, that we name localsearch, is described in algorithm 3 in the appendix. Our complete algorithm for finding the causal order is called Intersort (see algorithm 1).

Now arises the question of how close the solutions from the approximation algorithm intersort are to the optimum of eq. 1. We evaluate this empirically for graphs with $5$ variables, where computing the exact solution is still viable. The results are presented in fig. 1(a). As can be observed, the mean $D_{top}$ scores of the two algorithms are very close and have overlapping confidence intervals (at level $95\%$).

Figure: (a) Simulation with 5 variables
Refer to caption: extracted/5916887/figures/violin_bound_intersort_5.png

We note that more refined algorithms could be developed to solve this optimization over permutations. For example, for large $d$, the local search may become computationally expensive. However, for the purpose of this study, we find the proposed algorithm to have satisfying properties, as it scales to settings with $30$ variables and finds a solution close to the optimum. We also report results for a scale-free network modelled with the Barabasi-Albert distribution Albert & Barabási () in the appendix in fig. 3, which supports the generality of our approach. We also report the performance of SORTRANKING in large scale settings in figs. 5 and 4. The results demonstrate that the our score can be optimized at scale, even though we can observe room for improvement, as the error is above the upper-bounds for many settings. As such, we leave further algorithmic development for large scale settings as potential future work.

## 8 Estimating the distances from i.i.d samples

We now turn to the question of computing the statistical distances between observational and interventional distributions from real samples. Indeed, until now, we assumed access to population distributions, such that $D(P,Q)=0$ if and only if $P=Q$. However, in a real world setting, we may only assume that we have access to $i.i.d$ samples from the distributions, that is $\{X^{i}\}_{i=1}^{n_{0}}\sim P_{X}^{\mathcal{C},(\emptyset)}$ and for $j\in\{1,\ldots,d\},\{X^{i}\}_{i=1}^{n_{j}}\sim P_{X}^{\mathcal{C},do(X_{j}:=
\tilde{N}_{j})}$, where $n_{0},\ldots,n_{d}\in\mathbb{N}$ are the sample sizes of the observational set and of the interventional sets. Given this setting, we require a statistical distance with the following properties: it should be applicable to sample distributions, the distance between $i.i.d$ samples should converge to the distance between the population distribution as $n$ goes to infinity, with some guarantees on the convergence rate. A distance that fulfills all those criteria is the Wasserstein distance Villani et al. (). More specifically, it converges to the population distance in $\mathcal{O}\left(\frac{\sqrt{\ln{n}}}{\sqrt{n}}\right)$, where $n$ is the number of samples. We recall the definition of this distance in the appendix.

## 9 Empirical results

Figure: (a) Linear 30 variables
Refer to caption: extracted/5916887/figures/sweep_new_intersort_lin_30_violin_plot.png

We now evaluate Intersort on simulated empirical data and compare its performance to various baselines. The graphs are simulated from a Erdos-Renyi distribution (Erdős et al., ), with an expected number of edges per variable $c\in\{1,2\}$. We follow a setup close to (Lorch et al., ) to simulate data from linear and random Fourier features (RFF) additive functional relationships. We also apply the models to simulated single-cell data from the SERGIO (Dibaeinia & Sinha, ) model, using the code from Lorch et al. () (MIT License, v1.0.5). Lastly, we apply on neural network functional data following the setup of Brouillard et al. () and using the implementation of Nazaret et al. () (MIT License, v0.1.0). We run the models on varying levels of intervention, where the ratio of intervened variables is in $\{25\%,50\%,75\%,100\%\}$. The data is always standardized based on the mean and variance of the observational dataset, which removes the Varsortability (Reisach et al., ) artefact in the data. For the RFF and linear domain, the distribution of the noise is chosen uniformly at random from uniform Gaussian (noise scale independent from the parents), heteroscedastic Gaussian (noise scale functionally dependent on the parents), and Laplace. For the neural network domain, the noise distribution is Gaussian with a fixed variance. We run on $10$ simulated datasets for each domain and each ratio of intervened variables. We simulate $5000$ samples for the observational datasets and $100$ samples for each interventions, which is a setting similar to real single cell transcriptomics datasets (Replogle et al., ).

Choosing appropriate baselines to compare to is not trivial, given that, to our knowledge, our setting of predicting the causal order from interventional data has not been considered in the literature, as the focus in prior literature has been on observational data (Reisach et al., ; Bühlmann et al., ; Rolland et al., ; Shimizu et al., ). As such, we construct causal orders from the CPDAG predicted by two standard methods, namely PC (Spirtes et al., ) and GIES (Hauser & Bühlmann, ), by creating a DAG from the oriented edges and then computing the topological order of this DAG. PC is a constrained based method relying on conditional independence tests to estimate the structure of the graph. We run the implementation of Zheng et al. () (MIT License, causal-learn v0.1.3). GIES is an extension to interventional data of the GES model (Chickering, ), which is a score based method that greedily adds and removes edges to the estimated CPDAG. We use the Gaussian BIC score, and run the package implementation of Gamella () (BSD 3-Clause License, v0.0.1). Additionally, we compare to EASE (Gnecco et al., ), which is a method that aims to learn the causal order from observational data, leveraging the insight that extreme tail-values reveal the causal order. Given that interventions can lead to extreme values not present in the observational data, we also compare to this method. For EASE, we run our own Python implementation of the algorithm. Lastly, we compare to DCDI (Brouillard et al., ), which is a continuously differentiable causal discovery model that can leverage interventional data. For Intersort, we use $\epsilon=0.3$ for the linear, RFF and neural network domain, and $\epsilon=0.5$ for the single-cell domain. We set $c=0.5$. We compute the Wassertein distances with the SciPy python package (Virtanen et al., ).

The results for each method on the different types of data considered are displayed in fig. 2, where the distribution of the $D_{top}$ score of each model is plotted against the ratios of intervened variables considered. As can be observed, Intersort performs well on all data domains, and shows decreasing error as more interventions are available, exhibiting the model’s capability to capitalize on the interventional information to recover the causal order across diverse settings. Compared to the baselines, only GIES in the linear domain and DCDI in the neural network domain at $25\%$ and $50\%$ perform better. These experiments demonstrate that $\epsilon$-interventional faithfulness is fulfilled by a diverse set of data types, and that this property can be robustly exploited to recover causal information. Only when strong assumptions about the data distributions are fulfilled, such as linear data for GIES, can better results be obtained in low intervention settings.

## 10 Discussion and conclusion

We note that Intersort’s performance could still be improved in various ways. For example, a more scalable or closer to the optimum approximation algorithm could be designed. Also, the optimal parameter $\epsilon$ could be chosen in a data dependent way instead of using a fixed value. Lastly, other statistical distances, especially ones that converge faster to 0 0 in terms of number of samples when the two distributions are equal, could lead to better results. We leave such potential improvements as future work. Another question that may need further investigation is the $\epsilon$-interventional faithfulness assumption. We argue that it is not any stronger than common assumption in the literature such as faithfulness, and lemma 3 shows that it already covers a sizeable set of SCMs. Furthermore, our empirical results demonstrate that it probably holds for a rather large class of SCMs. We also derive theoretical results for a relaxation of this assumption in section 6. An interesting direction of theoretical work would be to further analyse the intricate relationship between the statistical distance $D$, the number of variables $d$, the parameter $\epsilon$, the interventional random variables $\tilde{N}$, and the sample size for non-asymptotic settings. For example, out framework makes it amenable to study the strength of the interventions needed to reveal the causal descendants above a statistically significant threshold.

In this study, we introduced Intersort, an innovative algorithm designed to uncover the causal order from interventional data by optimizing for a new proposed score on causal orders. We introduced the $\epsilon$-interventional faithfulness assumption and proved that interventional datasets fulfilling it have strong guarantees in terms of upper-bounds on the expected error of the optimum of our score. The performance of Intersort and the validity of our theoretical findings was extensively demonstrated on a diverse set of simulated datasets, across functional relationships, noise types and graph densities. Intersort has demonstrated superior performance when compared to established benchmarks such as PC, GIES, and EASE. The robustness and versatility of Intersort underscore the potential of the $\epsilon$-interventional faithfulness assumption to reshape causal inference methodologies. We envision that this work should spearhead new model development based on our proposed perspective on dataset with high numbers of single variable interventions. This includes downstream tasks such as causal discovery, for example by reconstructing the graph given the estimated causal order. It could inform active intervention selection as well, as designing a policy to choose which variables to target next may lead to better guarantees than our results based on interventions selected uniformly at random. All these directions and method developments can have real world impact on critical application domains such as biology, from which this work and its assumptions were inspired.

## Acknowledgements and disclosure

The authors thank Djordje Miladinovic, Aleksei Triastcyn and Lachlan Stuart for feedback and edits on the manuscript. The authors also thank Prof. Nicolai Meinshausen for valuable suggestions on the theoretical part, as well as Prof. Stefan Bauer for discussions about related work and relevance of the project. MC, AM and PS are employees and shareholders of GSK plc.

## Appendix A Proofs of Theorems

###### Proof of Lemma 3 .

Recall that for a linear SCM, the total causal effect between any pair $(i,j)$ is equal to the sum of the products of the weights along each paths from $i$ to $j$. Given that the weights are continuously random and independent, the probability that the sum is exactly 0 0 is 0 0. As such, the tuple $(\tilde{N},\mathcal{C})$ is $\epsilon$-interventionally faithful for $\epsilon=0$ almost surely.
∎

###### Proof of Theorem 1 .

Let $\mathbf{D}$ be a $d\times d$ matrix containing all the distances, where $\mathbf{D}(i,j)=D\left(P_{X_{j}}^{\mathcal{C},(\emptyset)},P_{X_{j}}^{\mathcal
{C},do(X_{i}:=\tilde{N}_{i})}\right)$. For the case of every variables being intervened on, we can set $c=0$.

We are going to prove that $\overline{A}_{\pi_{opt}}$ is upper-triangular, and from that it follows that $D_{top}(\mathcal{G},\pi_{opt})=0$. We know that such a permutation exists, as for any $\pi^{*}$, $\overline{A}_{\pi^{*}}$ is upper-triangular.

Second, we note that given that $(\tilde{N},\mathcal{C})$ is interventionally faithful, $\overline{A}$ has the non zero entries at the same position than the matrix $D$. We thus need to show that for any $D_{\pi}$ that is not upper-triangular, then $\pi\notin\operatorname*{arg\,max}_{\pi}\mathcal{S}(\pi)$. Then we can compute:

$$ $\displaystyle\mathcal{S}(\pi^{*})$ $\displaystyle=\sum_{\pi^{*}(i)<\pi^{*}(j),}\left(\mathbf{D}(i,j)-\epsilon\right)$ $\displaystyle=\sum_{\pi^{*}(i)<\pi^{*}(j)}\mathbf{D}\left(i,j\right)-\frac{d \cdot(d-1)}{2}\epsilon$ $\displaystyle=\sum_{i,j}\mathbf{D}\left(i,j\right)-\frac{d\cdot(d-1)}{2}\epsilon$ $$

Then, let $\pi$ be such that $D_{\pi}$ has some non-zero entries in the lower-triangular part. Then:

$$ $\displaystyle\mathcal{S}(\pi^{*})-\mathcal{S}(\pi)$ $\displaystyle=\sum_{\pi(i)>\pi(j)}\mathbf{D}(i,j)>0,$ $$

which implies that $\pi\notin\operatorname*{arg\,max}_{\pi}\mathcal{S}(\pi)$.
∎

###### Proof of lemma 4 .

Let $k$ be a variable that satisfies the condition, that is, $k\in\mathcal{I}$ and $k\in\textbf{AN}_{j}^{\mathcal{G}}\setminus\textbf{AN}_{i}^{\mathcal{G}}$. Then, for this $k$, the hypothesis implies that $D\left(P_{X_{j}}^{\mathcal{C},(\emptyset)},P_{X_{j}}^{\mathcal{C},do(X_{k}:=
\tilde{N}_{k})}\right)>\epsilon$ and $D\left(P_{X_{i}}^{\mathcal{C},(\emptyset)},P_{X_{i}}^{\mathcal{C},do(X_{k}:=
\tilde{N}_{k})}\right)=0$. We thus need to show that given those distances, $\pi_{opt}$ is such that $\pi_{opt}(i)<\pi_{opt}(j)$.

Let’s note that only distances that are larger than $\epsilon$ contribute positively to the score. All the distances that are equal to 0 0 contribute $-\epsilon$ to the sum. As such, a necessary condition for an optimal permutation is that all distances $D\left(P_{X_{i}}^{\mathcal{C},(\emptyset)},P_{X_{i}}^{\mathcal{C},do(X_{j}:=
\tilde{N}_{k})}\right)>\epsilon$ contribute to the sum, which happens only if $\pi_{opt}(i)<\pi_{opt}(j)$. Given that $(\tilde{N},\mathcal{C})$ is $\epsilon$-interventionally faithful, we have that any $\pi\in\Pi^{*}$ are such that all distances larger than $\epsilon$ contribute to $\mathcal{S}(\pi)$, which shows existence of such a permutation. $\pi_{opt}$ thus must have that $\pi_{opt}(k)<\pi_{opt}(j)$. If $i\in\mathcal{I}$, then with the same argument we have $\pi_{opt}(i)<\pi_{opt}(j)$.

It remains to show that $\pi_{opt}(i)<\pi_{opt}(k)$ if $i\notin\mathcal{I}$. We have $D\left(P_{X_{i}}^{\mathcal{C},(\emptyset)},P_{X_{i}}^{\mathcal{C},do(X_{k}:=
\tilde{N}_{k})}\right)=0$, which contributes negatively to the score if $\pi_{opt}(k)<\pi_{opt}(i)$. So changing to $\pi_{opt}(i)<\pi_{opt}(k)$ improves the score by $\epsilon$, as $i\notin\mathcal{I}$. This implies that $\pi_{opt}(i)$ must be smaller than $\pi_{opt}(k)$.

Lastly, if $j\in\mathcal{I}$, but not $i$ and no $k\in\textbf{AN}_{j}^{\mathcal{G}}\setminus\textbf{AN}_{i}^{\mathcal{G}}$, we need to show that we still have $\pi_{opt}(i)<\pi_{opt}(j)$. With the same argument as before, having $\pi_{opt}(i)<\pi_{opt}(j)$ can only improves the score as $i\notin\mathcal{I}$. This concludes the overall proof.
∎

###### Proof of theorem 2 .

$$ $\displaystyle\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]$ $\displaystyle=\mathbb{E}\sum_{(i,j)\in\mathcal{G}}\mathrm{1}_{\pi_{opt}(i)>\pi _{opt}(j)}$ $\displaystyle=\sum_{(i,j)\in\mathcal{G}}\mathbb{E}[\mathrm{1}_{\pi_{opt}(i)> \pi_{opt}(j)}]$ $\displaystyle=\sum_{(i,j)\in\mathcal{G}}P({\pi_{opt}(i)>\pi_{opt}(j)})$ $$

Now let’s look at $P({\pi_{opt}(i)>\pi_{opt}(j)})$, where there is an edge from $i$ to $j$ in the graph. Let’s note the condition of lemma 4 as $C$.

$$ $\displaystyle P({\pi_{opt}(i)>\pi_{opt}(j)})$ $\displaystyle=P({\pi_{opt}(i)>\pi_{opt}(j)}|\neg C)P(\neg C)+P({\pi_{opt}(i)> \pi_{opt}(j)}|C)P(C)$ $\displaystyle=P({\pi_{opt}(i)>\pi_{opt}(j)}|\neg C)P(\neg C)$ $\displaystyle=P({\pi_{opt}(i)>\pi_{opt}(j)}|\neg C)P(\forall k\in\textbf{AN}_{ j}^{\mathcal{G}}\cup\{j\}\setminus\textbf{AN}_{i}^{\mathcal{G}}:k\notin \mathcal{I})$ $\displaystyle=P({\pi_{opt}(i)>\pi_{opt}(j)}|\neg C)(1-p_{int})^{|\textbf{AN}_{ j}^{\mathcal{G}}\cup\{j\}\setminus\textbf{AN}_{i}^{\mathcal{G}}|}$ $\displaystyle\leq(1-p_{int})^{|\textbf{AN}_{j}^{\mathcal{G}}\cup\{j\}\setminus \textbf{AN}_{i}^{\mathcal{G}}|}$ $$

We can then take the sum over all the edges to conclude the proof.

∎

###### Proof of theorem 4 .

Without loss of generality, we assume that the variable are already in the causal order. We then write the error as a sum on all the upper-triangular elements that contain an edge and for which the variable are misplaced by $\pi_{opt}$.

$$ $\displaystyle\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]$ $\displaystyle=\mathbb{E}\sum_{i=1}^{d-1}\sum_{j=i+1}^{d}\mathrm{1}_{(i,j)\in E ,\pi_{opt}(i)>\pi_{opt}(j)}$ $\displaystyle=\sum_{i=1}^{d-1}\sum_{j=i+1}^{d}P((i,j)\in E,\pi_{opt}(i)>\pi_{ opt}(j))$ $\displaystyle=\sum_{i=1}^{d-1}\sum_{j=i+1}^{d}p_{e}\cdot P(\pi_{opt}(i)>\pi_{ opt}(j)|(i,j)\in E)$ $$

We now upper-bound the probability of error given the existence of an edge:

$$ $\displaystyle P(\pi_{opt}(i)>\pi_{opt}(j)|(i,j)\in E)$ $\displaystyle\leq P(j\notin\mathcal{I}\wedge\forall k<j,k\in\textbf{Pa}_{j}^{ \mathcal{G}}\wedge k\notin\textbf{Pa}_{i}^{\mathcal{G}}\implies k\notin \mathcal{I})$ $\displaystyle=(1-p_{int})\prod_{k=1}^{i}P(k\in\textbf{Pa}_{j}^{\mathcal{G}} \wedge k\notin\textbf{Pa}_{i}^{\mathcal{G}}\implies k\notin\mathcal{I})\prod_{ k=i+1}^{j-1}P(k\in\textbf{Pa}_{j}^{\mathcal{G}}\implies k\notin\mathcal{I})$ $\displaystyle\leq(1-p_{int})^{2}\prod_{k=i+1}^{j-1}P(k\in\textbf{Pa}_{j}^{ \mathcal{G}}\implies k\notin\mathcal{I})$ $\displaystyle=(1-p_{int})^{2}\prod_{k=i+1}^{j-1}(1-P(k\in\textbf{Pa}_{j}^{ \mathcal{G}}\wedge k\in\mathcal{I}))$ $\displaystyle=(1-p_{int})^{2}(1-p_{int}p_{e})^{j-i-1}$ $$

Plugging this back into the original sum and working out the geometric sums:

$$ $\displaystyle\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]$ $\displaystyle\leq{p_{e}(1-p_{int})^{2}}\sum_{i=1}^{d-1}\sum_{j=i+1}^{d}(1-p_{ int}p_{e})^{j-i-1}$ $\displaystyle={p_{e}(1-p_{int})^{2}}\sum_{i=1}^{d-1}\sum_{j=0}^{d-i-1}(1-p_{ int}p_{e})^{j}$ $\displaystyle={p_{e}(1-p_{int})^{2}}\sum_{i=1}^{d-1}\frac{1-(1-p_{int}p_{e})^{ d-i}}{p_{int}p_{e}}$ $\displaystyle=\frac{p_{e}(1-p_{int})^{2}}{p_{int}p_{e}}\left[d-\sum_{i=0}^{d-1 }(1-p_{int}p_{e})^{d-i}\right]$ $\displaystyle=\frac{(1-p_{int})^{2}}{p_{int}}\left[d-\sum_{i=1}^{d}(1-p_{int}p _{e})^{i}\right]$ $\displaystyle=\frac{(1-p_{int})^{2}}{p_{int}}\left[d-(1-p_{int}p_{e})\frac{1-( 1-p_{int}p_{e})^{d}}{p_{int}p_{e}}\right]$ $\displaystyle=\frac{(1-p_{int})^{2}}{p_{int}^{2}p_{e}}\left[p_{int}p_{e}\cdot d -(1-p_{int}p_{e})+(1-p_{int}p_{e})^{d+1}\right]$ $\displaystyle\leq\frac{(1-p_{int})^{2}}{p_{int}}d$ $$

∎

###### Proof of lemma 6 .

First, we can rewrite the equation on $p_{e}$ to get $p_{e}=\frac{c}{d}$. Then, plugging it into the bound we get:

$$ $\displaystyle\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]$ $\displaystyle\leq\frac{(1-p_{int})^{2}}{p_{int}}\left[d-(1-p_{int}p_{e})\frac{ 1-(1-p_{int}p_{e})^{d}}{p_{int}p_{e}}\right]$ $\displaystyle=\frac{(1-p_{int})^{2}}{p_{int}}\left[d-(1-p_{int}\frac{c}{d}) \frac{d}{c}\frac{1-(1-p_{int}\frac{c}{d})^{d}}{p_{int}}\right]$ $$

Finally, we normalize by $d$ and take the limit:

$$ $\displaystyle\lim_{d\rightarrow\infty}\frac{\mathbb{E}[D_{top}(\mathcal{G},\pi _{opt})]}{d}$ $\displaystyle\leq\lim_{d\rightarrow\infty}\frac{(1-p_{int})^{2}}{p_{int}}\left [1-(1-p_{int}\frac{c}{d})\frac{1-(1-p_{int}\frac{c}{d})^{d}}{c\cdot p_{int}}\right]$ $\displaystyle=\frac{(1-p_{int})^{2}}{p_{int}}\left[1-\frac{1-e^{-c\cdot p_{int }}}{c\cdot p_{int}}\right]$ $$

∎

###### Proof of Lemma 3 .

###### Proof of Theorem 1 .

###### Proof of lemma 4 .

###### Proof of theorem 2 .

###### Proof of theorem 4 .

###### Proof of lemma 6 .

## Appendix B Pen and paper example

We here spell out all the possible configurations for a graph with $3$ variables with the associated expected error of $\pi_{opt}$ in tables 1 and 2. Readers may find those example useful to build intuition on how the algorithm works.

**Table 1: Expected top-order divergence for graphs of size $3$, with $2$ out of the three variables observed under intervention.**
| Number edges | Adjacency Matrix | Transitive Closure | $\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]$ |
| --- | --- | --- | --- |
| 0 | $\begin{pmatrix}0&0&0\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | $\begin{pmatrix}0&0&0\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | 0 |
| 1 | $\begin{pmatrix}0&1&0\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | $\begin{pmatrix}0&1&0\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | 0 |
| 2 | $\begin{pmatrix}0&1&0\\<br>0&0&1\\<br>0&0&0\end{pmatrix}$ | $\begin{pmatrix}0&1&1\\<br>0&0&1\\<br>0&0&0\end{pmatrix}$ | 0 |
| 2 | $\begin{pmatrix}0&1&1\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | $\begin{pmatrix}0&1&1\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | 0 |
| 3 | $\begin{pmatrix}0&1&1\\<br>0&0&1\\<br>0&0&0\end{pmatrix}$ | $\begin{pmatrix}0&1&1\\<br>0&0&1\\<br>0&0&0\end{pmatrix}$ | 0 |

**Table 2: Expected top-order divergence for graphs of size $3$, with $1$ out of the three variables observed under intervention.**
| Number edges | Adjacency Matrix | Transitive Closure | $\mathbb{E}[D_{top}(\mathcal{G},\pi_{opt})]$ |
| --- | --- | --- | --- |
| 0 | $\begin{pmatrix}0&0&0\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | $\begin{pmatrix}0&0&0\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | 0 |
| 1 | $\begin{pmatrix}0&1&0\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | $\begin{pmatrix}0&1&0\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | $\frac{1}{3}\cdot\frac{1}{2}=\frac{1}{6}$ |
| 2 | $\begin{pmatrix}0&1&0\\<br>0&0&1\\<br>0&0&0\end{pmatrix}$ | $\begin{pmatrix}0&1&1\\<br>0&0&1\\<br>0&0&0\end{pmatrix}$ | $\frac{1}{3}\cdot\frac{1}{2}+\frac{1}{3}\cdot\frac{1}{2}=\frac{1}{3}$ |
| 2 | $\begin{pmatrix}0&1&1\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | $\begin{pmatrix}0&1&1\\<br>0&0&0\\<br>0&0&0\end{pmatrix}$ | $\frac{1}{3}\cdot\frac{1}{2}+\frac{1}{3}\cdot\frac{1}{2}=\frac{1}{3}$ |
| 3 | $\begin{pmatrix}0&1&1\\<br>0&0&1\\<br>0&0&0\end{pmatrix}$ | $\begin{pmatrix}0&1&1\\<br>0&0&1\\<br>0&0&0\end{pmatrix}$ | $\frac{1}{3}\cdot\frac{1}{2}+\frac{1}{3}\cdot\frac{1}{2}=\frac{1}{3}$ |

## Appendix C Pseudocodes

Figure: Algorithm 2 Finding an initial solution

Figure: Algorithm 3 Local search optimization

## Appendix D Wasserstein Distance

The $\mathcal{W}_{p}$ $p$-Wasserstein distance between two distributions is defined as follows (Villani et al., ):

###### Definition 2 .

The $p$-Wasserstein distance $\mathcal{W}_{p}$ between two probability measures $P$ and $Q$ on $\mathbb{R}^{d}$ with finite $p$-moments, $p\in[1,\infty)$,

$$ $W_{p}(P,Q)=\inf_{\gamma\in\Gamma(P,Q)}\left(\int_{X\times X}\|x-y\|^{p}d\gamma (x,y)\right)^{\frac{1}{p}}$ (4) $$

where $\Gamma(P,Q)$ is the set of joint probability measure on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginal $P$ and $Q$.

In this work, we use $p=1$, and we have $d=1$ as we only compute the distance between the marginal distribution of each variable. Let $P_{n}$ and $Q_{n}$ be distributions of samples from $P$ and $Q$ of size $n$. We have that the distance $W_{p}(P_{n},Q_{n})$ tends to $W_{p}(P,Q)$ as $n\rightarrow+\infty$. Furthermore, for $d\leq 2$, Ajtai et al. () show that:

$$ $\mathbb{E}|W_{p}(P_{n},Q_{n})-W_{p}(P,Q)|=\mathcal{O}\left(\frac{\sqrt{\ln{n}} }{\sqrt{n}}\right).$ (5) $$

###### Definition 2 .

## Appendix E Details of empirical evaluation

### E.1 Linear and RFF domain

Each causal variable $x_{j}$ is defined with respect to its parents $x_{\text{pa}(j)}$ through the equation:

$$ $x_{j}\leftarrow f_{j}(x_{\textbf{Pa}_{j}^{\mathcal{G}}},\epsilon_{j})=f_{j}(x_ {\textbf{Pa}_{j}^{\mathcal{G}}})+h_{j}(x_{\textbf{Pa}_{j}^{\mathcal{G}}}) \epsilon_{j}$ (6) $$

Here, $\epsilon_{j}$ represents additive, potentially heteroscedastic noise, where the noise scale $h_{j}(x_{\text{pa}(j)})$ is modeled as:

$$ $h_{j}(x)=\log(1+\exp(g_{j}(x)))$ (7) $$

with $g_{j}(x)$ being a nonlinear function. For the heteroscedastic case, we use random Fourier features with length scale $10.0$ and output scale of $2.0$.

To simulate interventions, the value of the intervened variable is set to fixed constant that is drawn from a signed Uniform distribution between $1.0$ and $5.0$.

#### E.1.1 Domain-Specific Modelling

- •
LINEAR Domain:
The causal functions $f_{j}$ are linear transformations:
$f_{j}(x_{\textbf{Pa}_{j}^{\mathcal{G}}})=w_{j}^{\top}x_{\textbf{Pa}_{j}^{
\mathcal{G}}}+b_{j}$
(8)
where $w_{j}$ and $b_{j}$ are independently sampled for each transformation. As in , we sample $w_{j}$ from a signed Uniform distribution between $1$ and $3$, and $b_{j}$ from a Uniform distribution between $-3$ and $3$.
- •
RFF Domain:
Each $f_{j}$ is drawn from a Gaussian Process (GP) approximated using random Fourier features:
$f_{j}(\mathbf{x}_{\textbf{Pa}_{j}^{\mathcal{G}}})=b_{j}+c_{j}\sqrt{\frac{2}{M}
}\sum_{m=1}^{M}\alpha^{(m)}\cos\left(\frac{1}{\ell_{j}}\omega^{(m)}\cdot
\mathbf{x}_{\textbf{Pa}_{j}^{\mathcal{G}}}+\delta^{(m)}\right)$
(9)
where $\alpha^{(m)}\sim\mathcal{N}(0,1)$, $\omega^{(m)}\sim\mathcal{N}(0,\mathbf{I})$, and $\delta^{(m)}\sim\text{Uniform}(0,2\pi)$, and the parameters $b_{j}$, $c_{j}$, and the length scale $\ell_{j}$ are sampled independently. The parameter $M$ is set to 100. The length scale $\ell_{j}$ is drawn from a Uniform distribution between $7.0$ and $10.0$, the output scale $c_{j}$ from a Uniform distribution between $10.0$ and $20.0$ and the bias $b_{j}$ from a Uniform distribution between $-3$ and $3$.

### E.2 Simulation of Single-Cell Gene Expression Data

As in Lorch et al. (), we use the SERGIO simulator by Dibaeinia & Sinha () to generate realistic high-throughput scRNA-seq data, implemented in the AVICI package by Lorch et al. ().

#### E.2.1 Simulation Process with SERGIO

SERGIO generates gene expression data as snapshots from the steady state of a dynamical system governed by the chemical Langevin equation. The system’s evolution is influenced by master regulators (MRs) with fixed rates and non-linear interactions among genes defined by the causal graph $G$. Cell types emerge from variations in MR rates, contributing to biological system noise and type-specific expression profiles.

#### E.2.2 Parameters for Simulation

To simulate $c$ cell types ($5$) across $d$ genes, SERGIO requires parameters such as interaction strengths $k$ (Uniform between $1.0$ and $5.0$), MR production rates $b$ (Uniform between $1.0$ and $3.0$), Hill coefficients $\gamma$ ($2.0$), decay rates $\lambda$ ($0.8$), and stochastic noise scales $\zeta$ ($1.0$). These are selected within recommended ranges to ensure realistic simulation. The values in parenthesis represents the parameter distributions we used. We did not simulate technical noise.

Intervention on a variable corresponds to knocking out a gene, that is, the expression of the gene is kept at 0 0.

### E.3 Simulation of Neural Network-Based Data for Causal Discovery

In the context of simulating data for causal discovery, we construct a domain where the data is generated by a set of random fully connected neural networks. These networks serve as the backbone for defining the observational conditional distributions within our simulated environment.

#### E.3.1 Neural Network Specification

Each neural network is a multilayer perceptron (MLP) with a single hidden layer consisting of 10 neurons. Formally, the MLPs are functions $\text{MLP}:\mathbb{R}^{j}\to\mathbb{R}^{1}$, equipped with ReLU activation functions, mapping a $j$-dimensional input to a scalar output. The output of each MLP represents the mean $\mu$ of the conditional distribution for the corresponding variable:

$$ $p_{j}(x_{j}|x_{\textbf{Pa}_{j}^{\mathcal{G}}})\sim\mathcal{N}(\mu=\text{MLP}(x _{\textbf{Pa}_{j}^{\mathcal{G}}}),\sigma=1.0).$ (10) $$

#### E.3.2 Interventional Data Generation

For the simulation of interventional data, the procedure modifies the distribution of intervened nodes. Instead of being determined by the MLP, the distribution of an intervened node is fixed to a marginal normal distribution:

$$ $p_{j}(x_{j}|\text{do}(x_{j}))\sim\mathcal{N}(2,1.0).$ (11) $$

This alteration simulates the effect of an intervention in the data generation process.

## Appendix F Additional experiments

### F.1 Simulations

Figure: (a) Simulation with 5 variables
Refer to caption: extracted/5916887/figures/violin_bound_intersort_sf_5.png

Figure: (a) Simulation with 1000 variables
Refer to caption: extracted/5916887/figures/violin_bound_intersort_sf_1000.png

Figure: (a) Simulation with 1000 variables
Refer to caption: extracted/5916887/figures/violin_bound_intersort_1000.png

### F.2 Empirical data

Figure: (a) Linear 10 variables
Refer to caption: extracted/5916887/figures/sweep_new_intersort_lin_10_violin_plot.png

Figure: (a) Linear 10 variables
Refer to caption: extracted/5916887/figures/results_sf_lin_10_violin_plot.png

Figure: (a) Linear 100 variables
Refer to caption: extracted/5916887/figures/results_test_const_lin_100_violin_plot.png

## Appendix G Computational resources

The simulations of fig. 1 were run on an Apple M1 Macbook pro. The experiments for figs. 2 and 6 were run on a cluster with 20 CPUs, 16Gb of memory per CPU.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: 1. What are the main differences between the current work and (Tian and Pearl, 2013) in terms of faithfulness assumption and methodology?

2. In lines 65-68, the authors are giving the advantage of knowing the order. Did they mean that the number of candidates is divided by 4 if we know that the target gene is in the middle of the causal order?

3. In eq. (1), in the second term, why is there a coefficient of $d$?

4. In Lemma 4, based on the chosen $p_e$, it seems that the graph is disconnected with high probability. Can the authors elaborate more on this?

5. In line 283, the authors argued that the expected error is growing with the order $O(d)$ and they mentioned that it is a strong guarantee. I think the upper bound on the expected error is $d^2$. Therefore, I am not sure this is indeed a strong result.

6. The pseudocode in Algorithm 1 is somehow useless. The most important parts (SORTRANKING and LOCACLSEARCH) are not described there.

7. Evaluating empirically does not imply the quality of the output of the algorithm in theory. I suggest removing the term "approximation algorithm" in line 359.

8. What is the computational complexity of the proposed algorithm? It is also good to compare the algorithms in terms of runtime empirically.

9. In Fig. 2, in some settings (such as NN 30 variables or GRN with a fraction of 23.33%), other methods have better performance. I think it is good to elaborate more on these cases in the experiment section.
soundness: 2
strengths: Originality/Quality: The authors asserted at the end of the Related Work section that they are proposing the first algorithm to infer the causal order from the interventional data. I do not think this is true. First, (Tian and Pearl, 2013) is one of the earliest works using the changes in marginal distributions due to intervention in inferring some orders among the variables in the system (See Section 4 there). Second, there are extensive works on recovering an equivalence class of models from the observational and interventional data (whether the intervention locations are known or not). A few works are cited in the paper (such as GIES and DCDI). These works can also provide some information about the causal orders that are encoded in the recovered equivalence class. 

Clarity: The paper is generally well-written.

Significant: Based on what was mentioned, I think that the authors should carefully compare their methods with previous work (such as (Tian and Pearl, 2013) when $\epsilon=0$). Moreover, from the experimental results, it is not clear that the proposed algorithm indeed improves SOTA.
summary: The authors considered the problem of recovering a causal order in structural causal models (SCM) under the causal sufficiency assumption (it seems that is the case but it is not clearly mentioned in the paper). They assumed that they have access to observational data and some single-variable interventional ones and used the changes observed in the marginal distributions of variables to recover a causal order. The authors introduced a notion of faithfulness (interventional faithfulness) and showed that a causal order can be recovered by maximizing a score (defined in (1)) if we have single-variable interventions on all the variables in the system. They proposed a heuristic algorithm (called INTERSORT) aiming to improve the score. The experimental results showed that in some specific settings, the proposed algorithm has a better performance compared to some of the previous work.
weaknesses: Comparison with previous work: I think there is no clear comparison with previous work (especially with (Tian and Pearl, 2013)) and discussion about the advantages of the current approach. 

Theoretical result in a very limited setting: I think the assumption of having single-variable intervention on all the variables is very restrictive (what can we say in theory about the recovered causal order if a portion of variables are intervened on?). Moreover, the proposed algorithm is designed under the causal sufficiency assumption (which as far as I checked, is not clearly mentioned in the paper). 

The notion of interventional faithfulness: I think this assumption is an extension of "influentially" in (Tian and Pearl, 2013). The connection to that definition is not discussed in the paper.

Theoretical guarantee about INTERSORT: Although the authors used the term "approximation algorithm" in line 359, there is no theoretical guarantee about the quality of the output of INTERSORT.

---

confidence: 4
contribution: 4
initial_rating: 8
presentation: 4
questions: 1. Are there scenarios where the proposed method cannot be applied? In other words, are there real-world systems that do not satisfy e-interventional faithfulness or its relaxations?
2. How could INTERSORT be extended to handle a large number of nodes? I presume for calculating distance, one could replace Wasserstein with e.g. MMD, but what would be the main bottleneck of the algorithm aside from the distance metric?
3. Could the authors elaborate on how the proposed framework could be extended to derive causal order over a subset of variables?
4. Do the authors have any intuition on how this framework might be used to develop better causal discovery algorithms?
soundness: 3
strengths: * The paper introduces a new definition of faithfulness with both theoretical guarantees and empirical results.
* It is well-written and clear.
* Numerical experiments are designed in a sensible manner that adequately supports the claims.
* This is an important and highly relevant contribution to the community, with developed theory that has the potential to further advance causality.
summary: The paper introduces the concept of "interventional faithfulness," which relies on comparisons between the marginal distributions of each causal variable across observational and interventional settings, and also develops a scoring function for ranking causal orders of variables. Based on this concept, the authors 1) provide strong theoretical guarantees on the optimum of the proposed score and 2) propose INTERSORT to infer causal order from datasets containing large amounts of single-variable interventions. INTERSORT outperforms existing methods in various simulations, demonstrating the potential of the new theory for advancing causal inference in domains like biology.
weaknesses: * It’s unclear how the findings would generalize to different settings, such as different distributions over random interventional variables, or in the case of having discrete causal variables.
* Empirical experiments are conducted on a limited set of underlying models.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: Are there useful heuristics for choosing which variables to intervene on ($\mathcal{I}$)?
soundness: 4
strengths: This paper presents an original solution to a significant problem: a method for directly learning the causal order using interventional data.
 * Although technically dense, the paper is surprisingly easy to read and well-organized.
 * The theoretical guarantees for optimality under a reasonable faithfulness assumption appear solid.
 * The approximation algorithm is computationally tractable and validated with a variety of benchmarks.
summary: Learning the topological ordering of the causal graph behind a set of variables is useful in a variety of scientific domains. For instance, when gene expressions are correlated, the causal order can help discern the structure of the regulatory network. This paper proposes a novel approach for ordering variables using many single-variable interventions. Strong theoretical guarantees and diverse empirical evaluations are presented.
weaknesses: * The majority of the analysis assumes access to oracle statistical distances between interventional distributions. Little attention is paid to the estimation of these distances using samples, and how they affect the sorting algorithm. 
 * The main score objective (Equation 1) is difficult to understand and not explained much.

Minor comments
 * Please define all the terms on like 152, like the noise variable $N_j$.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: 1. If setting the threshold following the Remark in Line 203, then it means every pair with distribution change after perturbation will be considered. Why not directly set it to 0, which means detecting distribution change directly?
2. The selection processes also have a large distance, however, this is not a causal relationship. Moreover, selection is common in gene regulatory networks. In your method, how to guarantee the ones you find are causal relations?
3. In eq. (1), what does the second part mean? I scratch my head to understand it.
4. In Algorithm 1, how many kinds of permutations will be considered. Usually, it will be non-traversable with an increase in the number of variables.
5. In the part of Distribution with intervention, ' the structural assignment is replaced by a new random variable independent of the parents' means the hard intervention right? Moreover, the author also mentioned that distribution after intervention can be accessed, do you have assumptions about the distribution of intervention? and why you need the distribution needs to be accessed.
soundness: 2
strengths: The score of causal order is designed to identify causal relationships with the help of distribution change engendered by intervention. Theoretically, the upper bound of the expected error is provided for the optimal causal order.
summary: This paper proposed INTERSORT, an algorithm for causal discovery with interventional data, to find a casual relationship with distribution change after intervention. Theoretically, the upper bound of the error is provided for the optimal causal order. Experimental results on simulated and real-world datasets verify the effectiveness of the proposed method.
weaknesses: 1. If setting the threshold following the Remark in Line 203, then it means every pair with distribution change after perturbation will be considered. Why not directly set it to 0, which means detecting distribution change directly?
2. The selection processes also have a large distance, however, this is not a causal relationship. Moreover, selection is common in gene regulatory networks. In your method, how to guarantee the ones you find are causal relations?
3. In eq. (1), what does the second part mean? I scratch my head to understand it. Could you provide descriptions of the symbols used in the equation?
4. In Algorithm 1, how many kinds of permutations will be considered. Usually, it will be non-traversable with an increase in the number of variables.
5. In the part of Distribution with intervention, ' the structural assignment is replaced by a new random variable independent of the parents' means the hard intervention right? Moreover, the author also mentioned that distribution after intervention can be accessed, do you have assumptions about the distribution of intervention? and why you need the distribution needs to be accessed.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: - What would be an application in which one would only be interested in inferring the causal order and not the full graph/Markov equivalence class?
- How exactly is the algorithm implemented? E.g., which p-Wasserstein distance is chosen? Which algorithm is used to compute the distance?
soundness: 3
strengths: - The proposed method is the first to estimate the causal order from interventional data
- The paper contains several novel and theoretically grounded contributions, including the causal score as well as the algorithm to approximate it
- Promising empirical results on simulated data
summary: The paper proposes a new method to infer the causal order of an (unknown) causal graph from interventional data. For this purpose, the authors introduce a new assumption (interventional faithfulness), a theoretically grounded score on causal orders, and a practical algorithm called intersort approximating the proposed score. Intersort is evaluated on simulated data.
weaknesses: - As common in causal discovery, the method relies on generally untestable assumptions (e.g., epsilon interventional faithfulness)
- The method relies on estimating distributional distances (e.g., Wasserstein distance) which can be statistically challenging
- No experiments using real-world data are provided. I understand that benchmarking methods are challenging due to unknown causal ground truth, however, I think it would be nice to sketch an application of the method on real-world data and potentially obtain insights
- The related work on causal discovery could be expanded

Minor 
- Appendix B seems to be missing
- There is some weird formatting on page 16 (Appendix C and D intersect)
- There is no text in Appendix F.1/ F.2 so it is hard to understand which figures belong to which section. Also, Appendix G intersects with F2

Meta-review: The paper introduces an approach for learning the order of a causal DAG model from single variable interventions. Reviewers found it useful, and I'm generally favourable to accept based on the framing of the algorithms being novel and for providing results others may want to build upon. But I believe scholarship could be improved, and the claims of originality better contextualized: for avoidance of doubt, the paper is not the first to address causal ordering learning from interventional data, a problem which broadly speaking is not hard under some typical assumptions. Some of my points below are mentioned by reviewer R6EY, and I'm disappointed that they are not addressed by the authors.

* I'm somewhat confused why the authors mention Tian and Pearl repeatedly as a key citation for their contribution. First, please correct the reference: Tian and Pearl's is a 2001 (!) paper from UAI. I see limited connections between the motivation of that paper and the presented paper. Tian and Pearl's main point was how to make use of natural experiments, while my interpretation is that the present contribution is motivated by experimental design. If the authors want to refer to a historical paper on using interventional data, it would make far more sense to refer to Cooper and Yoo's "Causal discovery from a mixture of experimental and observational data" (UAI, 1999). Tian and Pearl explain very clearly the connection between their paper and (last century's) literature on causal discovery from interventional data, please refer to them for details.

* As correctly mentioned in the paper, the notion of learning from pairwise orderings appears in LiNGAM-like algorithms, including DirectLiNGAM (Shimizu et al., 2011). Readers would benefit from an explicit discussion between methods at the more abstract level of being told pairwise directions, regardless of the source of data (observational vs experimental) and parametric family (linear or not). But claims such as "Reisach et al. (2021) that the causal order can be recovered by sorting the variable by variance in simulated datasets" are strangely phrased: causal order can *not* be recovered by variance in simulated datasets in general (!), only those which are either badly designed or, for some hopefully well-justified reason that is nothing but a very special case, follow a physical process where variance increases with depth of causal connection (e.g., some classes of additive error models with constant error term variance).

* Like reviewer R6EY, I'm surprised that papers like Eberhardt, Glymour and Scheines "On the Number of Experiments Sufficient and in the Worst Case Necessary to Identify All Causal Relations Among N Variables" (UAI, 2005) and the older "N-1 Experiments Suffice to Determine the Causal Relations Among N Variables" (published as a tech report in 2004) are not cited. Please do so, and explain the connections. Those papers solve a closely related (and to some extent, harder) problem than the one posed in this paper, using more classical faithfulness. In fact, a paper such as Saengkyongam and Silva "Learning Joint Nonlinear Effects from Single-variable Interventions in the Presence of Hidden Confounders" (UAI, 2020) jumps straight into the problem of causal effect identification with a given structure because it recognizes that ordering can be readily learned from the same data using methods such as simpler versions of Eberhardt et al. (as only ordering is necessary, instead of a full causal graph).

* broadly speaking, "interventional faithfulness" is not a novel concept as mentioned in the abstract, as numerous variants of it appear in papers where discovery with a mixture of interventional datasets take place, including some cited (any method which claims correctness of causal discovery with interventional data needs some notion of interventional faithfulness). The $\epsilon$ variant given here is novel to the best of my understanding, but it is worthwhile pointing to the reader that there are several $epsilon$-like variants of faithfulness, for instance the one found in Zhang and Spirtes "Strong Faithfulness and Uniform Consistency in Causal Inference" (UAI, 2003). The notion introduced in the paper is novel and worth discussing as far as I know, but it should be noticed that much work has been done in relaxing faithfulness assumptions (for instance, from the citations of Zhang and Spirtes). Tian and Pearl's influentially is a type of marginal faithfulness and has no extra degree of freedom.

If this paper is to be published at ICLR, I hope all of the above is addressed in the final version -- which I think it's doable and I'll assume good faith from the authors.

Decision: Accept (Poster)