Efficient Differentiable Discovery of Causal Order

## Abstract

Abstract In the algorithm Intersort, Chevalley et al. ( 2024 ) proposed a score-based method to discover the causal order of variables in a Directed Acyclic Graph (DAG) model, leveraging interventional data to outperform existing methods. However, as a score-based method over the permutahedron, Intersort is computationally expensive and non-differentiable, limiting its ability to be utilised in problems involving large-scale datasets, such as those in genomics and climate models, or to be integrated into end-to-end gradient-based learning frameworks. We address this limitation by reformulating Intersort using differentiable sorting and ranking techniques. Our approach enables scalable and differentiable optimization of causal orderings, allowing the continuous score function to be incorporated as a regularizer in downstream tasks. Empirical results demonstrate that causal discovery algorithms benefit significantly from regularizing on the causal order, underscoring the effectiveness of our method. Our work opens the door to efficiently incorporating regularization for causal order into the training of differentiable models and thereby addresses a long-standing limitation of purely associational supervised learning.

## 1 Introduction

Causal discovery is fundamental for understanding complex systems by identifying underlying causal relationships from data. It has significant applications across various fields, including biology (Meinshausen et al., ; Chevalley et al., , ), medicine (Feuerriegel et al., ), and social sciences (Imbens & Rubin, ), where causal insights inform decision-making and advance scientific knowledge. Traditionally, causal discovery has relied heavily on observational data due to the practical challenges and costs associated with conducting large-scale interventional experiments. However, observational data alone often necessitates strong assumptions about the data distribution to ensure identifiability beyond the Markov equivalence class (Spirtes et al., ; Shimizu et al., ; Hoyer et al., ).

The emergence of large-scale interventional datasets, particularly in domains like single-cell genomics (Replogle et al., ; Datlinger et al., ; Dixit et al., ), introduces new opportunities and challenges for causal discovery. Interventional data, obtained through targeted manipulations of variables, offers a unique perspective by revealing causal mechanisms that may be obscured in purely observational studies. Recently (Chevalley et al., ) introduced Intersort, which uses the notion of *interventional faithfulness* and a score on causal orders, enabling the inference of causal orderings by comparing marginal distributions across observational and interventional settings.

Despite its promising performance, Intersort faces significant limitations. Specifically, it lacks differentiability, which hinders its integration into gradient-based learning frameworks commonly used in modern machine learning. Additionally, scalability remains a challenge, making this method impractical for applications involving large numbers of variables—a common scenario in fields such as genomics (Replogle et al., ; Chevalley et al., ) and neuroscience that involve up to tens of thousands of variables.

In this work, we address these challenges by extending Intersort to a continuous realm to make it both differentiable and scalable. We achieve this by expressing the score in terms of a potential function and utilizing differentiable sorting and ranking techniques, including the Sinkhorn operator (Cuturi, ). This novel formulation allows us to incorporate a useful inductive bias into downstream tasks as a differentiable constraint or regularizer, enabling seamless integration into gradient-based optimization frameworks.

We develop a causal discovery algorithm that integrates differentiable Intersort score into its objective function. Our empirical evaluations on diverse simulated datasets—including linear, random Fourier features, gene regulatory networks (GRNs) and neural network models—demonstrate that the proposed regularized algorithm significantly outperforms baseline methods such as GIES (Hauser & Bühlmann, ) and DCDI (Brouillard et al., ) on RFF and GRN data. Moreover, we demonstrate that our approach exhibits robustness across different data distributions and noise types. The algorithm efficiently scales with large datasets, maintaining consistent performance regardless of data size.

Our contributions pave the way for fully leveraging interventional data in various causal tasks in a scalable and differentiable manner. By addressing the limitations of existing methods, we enable the application of interventional faithfulness in large-scale settings and facilitate its integration into modern causal machine learning pipelines.

## 2 Related work

Causal discovery is a fundamental problem in machine learning and statistics, aiming to uncover causal relationships from data. Traditional methods can be broadly categorized into constraint-based, score-based, and functional approaches.
Algorithms like PC (Spirtes et al., ) and FCI (Spirtes, ) use conditional independence tests to infer causal structures. These methods often require faithfulness and causal sufficiency assumptions, which may not hold in real-world scenarios.
Score-based approaches, such as GES (Chickering, ) and GIES (Hauser & Bühlmann, ), search over possible graph structures to maximize a goodness-of-fit score. More recently, Zheng et al. () introduced NOTEARS, a differentiable approach that formulates causal discovery as a continuous optimization problem using an acyclicity constraint. Most recently proposed models follow the idea of NOTEARS and are thus continuously differentiable, e.g. using neural networks to model the functional relationships (Brouillard et al., ; Lachapelle et al., ).

Chevalley et al. () recently introduced Intersort, a score-based algorithm to derive the causal order of the variables when many single-variable interventions are available. Intersort relies on a light assumption on the changes in marginal distributions across observational and interventional settings, as measure by a statistical distance (see definition 3 subsequently). Intersort is a two-step algorithm, where the first step, SORTRANKING, finds an initial ordering by taking into account the scale of the measured distances, and the second step, LOCALSEARCH, search in a close neighbourhood in permutation space to improve the score. Their algorithm suffers from scalability issues, as the second step, LOCALSEARCH, becomes computationally intractable for large datasets. Moreover, it is unclear how this inductive bias can be applied to downstream causal task, as the score is not differentiable with respect to a causal order and thus it cannot be integrated into continuously differentiable models. Previous to Intersort, most methods to infer the causal order of a system worked primarily on observational data. For example, EASE (Gnecco et al., ) leverages extreme values to identify the causal order. Reisach et al. () observed that the performance of continuously differentiable causal discovery models on synthetic datasets strongly correlates with *varsortability*, which measures to what extent the marginal variance corresponds to the true causal order. Similar ideas were consequently developed to derive the causal order based on score matching (Rolland et al., ; Montagna et al., , ).

## 3 Method

### 3.1 Definitions and Assumptions

In this section, we introduce notations and definitions that are used throughout the paper inspired by (Pearl, ; Peters et al., ).

Let $(\mathcal{M},d)$ be a metric space, and let $\mathcal{P}(\mathcal{M})$ denote the set of probability measures over $\mathcal{M}$. We define $D$ to be a statistical distance function $D:\mathcal{P}(\mathcal{M})\times\mathcal{P}(\mathcal{M})\rightarrow[0,\infty)$ that measures the divergence between probability distributions on $\mathcal{M}$.
Consider a set of $d$ random variables $\mathbf{X}=(X_{1},X_{2},\dots,X_{d})$ indexed by $V=\{1,2,\dots,d\}$, with joint distribution $P_{\mathbf{X}}$. We denote the marginal distribution of each variable as $P_{X_{i}}$ for $i\in V$.
A causal graph is a tuple $\mathcal{G}=(V,E)$ of nodes and edges that form a Directed Acyclic Graph (DAG), where $V$ is the set of nodes (variables), and $E\subseteq V\times V$ is the set of directed edges representing causal relationships. An edge $(i,j)\in E$ indicates that variable $X_{i}$ is a direct cause of variable $X_{j}$.
Let $\mathbf{A}^{\mathcal{G}}$ be the adjacency matrix of $\mathcal{G}$, where $\mathbf{A}^{\mathcal{G}}_{ij}=1$ if $(i,j)\in E$, and $\mathbf{A}^{\mathcal{G}}_{ij}=0$ otherwise. For each node $j\in V$, the set of parents $\mathrm{Pa}(j)$ consists of all nodes with edges pointing to $j$, i.e., $\mathrm{Pa}(j)=\{i\in V\mid(i,j)\in E\}$.
We denote the set of descendants of node $i$ as $\mathrm{De}_{\mathcal{G}}(i)$, which includes all nodes reachable from $i$ via directed paths. Similarly, the set of ancestors of $i$ is denoted as $\mathrm{An}_{\mathcal{G}}(i)$.

An SCM $\mathcal{C}=(\mathbf{S},P_{N})$ consists of a set of structural assignments $\mathbf{S}$ and a joint distribution over exogenous noise variables $P_{N}$. Each variable $X_{j}$ is assigned via a structural equation:

$$ $X_{j}=f_{j}\left(\mathbf{X}_{\mathrm{Pa}(j)},N_{j}\right),$ $$

where $N_{j}$ is an exogenous noise variable, and $\mathbf{X}_{\mathrm{Pa}(j)}$ are the parent variables of $X_{j}$.

In our work, we focus on interventions that modify the structural assignments of certain variables. Specifically, we consider interventions where the structural assignment of a variable $X_{k}$ is replaced by a new exogenous variable $\tilde{N}_{k}$, independent of its parents
$X_{k}=\tilde{N}_{k}.$

###### Definition 1 .

A *causal order* is a permutation $\pi:V\rightarrow\{1,2,\dots,d\}$ such that for any edge $(i,j)\in E$, we have $\pi(i)<\pi(j)$. This ensures that causes precede their effects in the ordering (Peters et al., ).

Since $\mathcal{G}$ is acyclic, at least one causal order exists, though it may not be unique. We denote the set of all causal orders consistent with $\mathcal{G}$ as $\Pi^{*}$.

###### Definition 2 .

To measure the discrepancy between a proposed permutation $\pi$ and the true causal graph $\mathcal{G}$, we use the *top order divergence* (Rolland et al., ), defined as:

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

This divergence counts the number of edges that are inconsistent with the ordering $\pi$, i.e., edges where the cause appears after the effect in the proposed ordering. For any causal order $\pi^{*}\in\Pi^{*}$, we have $D_{\text{top}}(\mathcal{G},\pi^{*})=0$.

###### Assumption 1 (Interventional Faithfulness) .

Interventional faithfulness (Chevalley et al., ) assumes that all directed paths in the causal graph manifest as significant changes in the distribution under interventions as measured by a statistical distance.
Specifically, if intervening on variable $X_{i}$ leads to a detectable change in the distribution of variable $X_{j}$, then there must be a directed path from $X_{i}$ to $X_{j}$ in the causal graph $\mathcal{G}$. Conversely, if there is no directed path from $X_{i}$ to $X_{j}$, then intervening on $X_{i}$ does not affect the distribution of $X_{j}$ beyond a significance threshold $\epsilon$.

Interventional faithfulness allows us to use statistical divergences between observational and interventional distributions to infer the causal ordering of variables. By assuming interventional faithfulness, we can relate changes observed under interventions to the underlying causal structure. More formally, it is defined as:

###### Definition 3 ( Chevalley et al. ( 2024 ) ) .

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}$.

###### Definition 1 .

###### Definition 2 .

###### Assumption 1 (Interventional Faithfulness) .

###### Definition 3 ( Chevalley et al. ( 2024 ) ) .

### 3.2 Differentiable score

While Intersort demonstrates cutting-edge results in discerning causal order among variables, its primary drawback is the substantial computational cost, which restricts its application to small-scale problems. The authors of the original paper acknowledged this limitation, confining their evaluation to a mere 30 nodes Chevalley et al. (). A covariate set of this size is prohibitively small for many real-world problems, such as those in genomics and climate change, where tens of thousands of variables are considered. We aim to enhance the scalability of Intersort through a differentiable objective function. This not only facilitates scaling to a considerably larger number of variables but also enables the integration of this algorithm in end-to-end gradient-based model training. In the subsequent sections, we initially revisit the fundamental score that underpins Intersort. Following this, we proceed to present a differentiable formulation, DiffIntersort, that addresses these shortcomings.

Intersort score– 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$, Chevalley et al. () 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\geq 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) $$

We parameterize an ordering of the variable as determined by a permutation of the variables $\pi$ via a potential ${\bm{p}}\in\mathbb{R}^{d}$ such that $\pi(i)<\pi(j)\iff p_{i}>p_{j}$. We write the permutation matrix associated to ${\bm{p}}$ as ${\bm{\sigma}}({\bm{p}})$, which is a $d\times d$ binary matrix. Let $s:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ be the sorting function, which sorts a vector of $d$ numbers in descending order. The Jacobian of $s({\bm{p}})$ with respect to ${\bm{p}}$ is equal to the permutation matrix ${\bm{\sigma}}({\bm{p}})$. We define $(\operatorname*{grad}({\bm{p}}))_{ij}=p_{i}-p_{j}$, which is non-negative if and only if $\pi(i)<\pi(j)$ in the associated topological order. Applying the element-wise $\operatorname*{Step}$ function produces $(\operatorname*{Step}(\operatorname*{grad}({\bm{p}})))_{ij}=\mathrm{1}_{p_{i}-
p_{j}>0}$ which is a matrix of the possible edges according to the potential ${\bm{p}}$.

We aim to rewrite the score such that it is parameterized by a potential ${\bm{p}}$. By building the matrix ${\mathbf{D}}\in\mathbb{R}^{d\times d}$ as

$$ ${\mathbf{D}}_{ij}=\left\{\begin{array}[]{lcl}D\left(\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}&\mbox{if}& i\in\mathcal{I}\\ 0&\mbox{if}&i\notin\mathcal{I}\end{array}\right.$ (2) $$

we can write the score in terms of the potential instead of permutation as follows:

$$ $S({\bm{p}},\epsilon,D,\mathcal{I},P_{X}^{\mathcal{C},(\emptyset)},\mathcal{P}_ {int},c)=\sum_{i,j}\left({\mathbf{D}}\odot\operatorname*{Step}(\operatorname*{ grad}({\bm{p}}))\right)_{ij}.$ (3) $$

The relationship between the potential and permutation is clarified through the following theoretical result.

###### Theorem 1 .

Let ${\mathbb{P}}=\operatorname*{arg\,max}_{{\bm{p}}}S({\bm{p}},\epsilon,D,\mathcal
{I},P_{X}^{\mathcal{C},(\emptyset)},\mathcal{P}_{int},c)\text{ s.t. }{\bm{p}}_
{i}\neq{\bm{p}}_{j}\forall i,j\in\{1,\cdots,d\}$ be the set of potentials that maximize the score, such that no two potentials are equal. $\Pi=\operatorname*{arg\,max}_{\pi}S(\pi,\epsilon,D,\mathcal{I},P_{X}^{\mathcal
{C},(\emptyset)},\mathcal{P}_{int},c)$ be the set of permutations that maximize the Intersort score. For all $\pi\in\Pi$, there is a ${\bm{p}}\in{\mathbb{P}}$ such that $\pi(i)<\pi(j)\iff p_{i}>p_{j}$.

The proof can be found in the appendix in section 6.1. This score is still not practically useful as it provides non-informative gradients for ${\bm{p}}$. To remedy this, inspired by Annadani et al. () we define ${\bm{L}}\in\{0,1\}^{d\times d}$ as a matrix with upper triangular part to be $1$, and vector ${\bm{o}}=[1,\ldots,d]$. They propose the formulation

$$ $\displaystyle\operatorname*{Step}(\operatorname*{grad}({\bm{p}}))={\bm{\sigma} }({\bm{p}}){\bm{L}}{\bm{\sigma}}({\bm{p}})^{T}\$ $\displaystyle\text{where}\;{\bm{\sigma}}({\bm{p}})=\operatorname*{arg\,max}_{{ \bm{\sigma}}^{\prime}\in\bm{\Sigma}_{d}}{\bm{p}}^{T}({\bm{\sigma}}^{\prime}{ \bm{o}})$ (4) $$

where $\bm{\Sigma}_{d}$ represents the space of all $d$ dimensional permutation matrices. This formulation opens the door to use tools from the differentiable permutation optimization literature. More specifically, we need to build a smooth approximation to the argmax operator in the definition of ${\bm{\sigma}}({\bm{p}})$. Fortunately, theoretical results for this transformation is already known in the optimization literature which is built upon the concept of Sinkhorn Operator that we briefly discuss in the following. We refer readers to the original paper (Sinkhorn, ) and further applications (Adams & Zemel, ) of Sinkhorn operator for detailed presentation.

The Sinkhorn operator, ${\mathcal{S}}({\bm{M}})$, on a matrix ${\bm{M}}$, involves a sequence of alternating row and column normalizations, known as Sinkhorn iterations. Mena et al. () demonstrated that the non-differentiable $\operatorname*{arg\,max}$ problem

$$ ${\bm{\sigma}}=\operatorname*{arg\,max}_{{\bm{\sigma}}^{\prime}\in\Sigma_{d}} \left\langle{\bm{\sigma}}^{\prime},{\bm{M}}\right\rangle_{F}$ (5) $$

can be relaxed using an entropy regularizer, where the solution is given by ${\mathcal{S}}({\bm{M}}/t)$. Specifically, they showed that ${\mathcal{S}}({\bm{M}}/t)=\operatorname*{arg\,max}_{{\bm{\sigma}}^{\prime}\in{
\mathcal{B}_{d}}}\left\langle{\bm{\sigma}}^{\prime},{\bm{M}}\right\rangle+th({
\bm{\sigma}}^{\prime})$, where $h(\cdot)$ denotes the entropy function. This regularized solution converges to the solution of eq. 5 as $t\rightarrow 0$, shown by $\lim_{t\rightarrow 0}{\mathcal{S}}({\bm{M}}/t)$. The implication for our setting is that we can write eq. 4 as

$$ $\operatorname*{arg\,max}_{{\bm{\sigma}}^{\prime}\in\bm{\Sigma}_{d}}{\bm{p}}^{T }({\bm{\sigma}}^{\prime}{\bm{o}})=\operatorname*{arg\,max}_{{\bm{\sigma}}^{ \prime}\in\bm{\Sigma}_{d}}\langle{\bm{\sigma}}^{\prime},{\bm{p}}{\bm{o}}^{T} \rangle_{F}=\lim_{t\rightarrow 0}{\mathcal{S}}\left(\frac{{\bm{p}}{\bm{o}}^{T} }{t}\right).$ (6) $$

In practice, we approximate the limit with a value of $t>0$ and a certain number of iterations $T$, which results in a differentiable and doubly stochastic matrix in the $d$-dimensional
Birkhoff polytope ${\mathcal{B}_{d}}$. The parameter $t>0$ acts as a temperature controlling the smoothness of the approximation. In our experiments, we use $t=0.05$ and $T=500$. We then apply the Hungarian algorithm Kuhn () to obtain a binary matrix, and use the straight-through estimator in the backward pass. The resulting binary matrix is denoted as ${\mathcal{S}}^{T}_{\rm bin}({\bm{p}}{\bm{o}}^{T}/t)$ with "bin" emphasizing a binary-valued matrix. As a result, the score becomes differentiable and can be differentiated through the iterations of the Sinkhorn operator. By replacing the non-differentiable part of eq. 2 with this matrix, the complete form of the differentiable score (we call it DiffIntersort) is derived as

$$ $S({\bm{p}},\epsilon,D,\mathcal{I},P_{X}^{\mathcal{C},(\emptyset)},\mathcal{P}_ {int},t,T)=\sum_{i,j}\left({\mathbf{D}}\odot\left({\mathcal{S}}^{T}_{\textrm{ bin}}\left(\frac{{\bm{p}}{\bm{o}}^{T}}{t}\right){\bm{L}}{\mathcal{S}}^{T}_{ \textrm{bin}}\left(\frac{{\bm{p}}{\bm{o}}^{T}}{t}\right)^{T}\right)\right)_{ij}.$ (7) $$

For the rest of the paper, we drop in subscript "bin" and use $S({\bm{p}})$ for conciseness.The maximizers of the DiffIntersort score and the Intersort score are equal for $t\rightarrow 0$ and $T\rightarrow\infty$ (Theorem 1). The DiffIntersort score $S(\mathbf{p})$ can be maximized with respect to the potential vector $\mathbf{p}$ using gradient descent algorithms. This allows us to find the ordering of variables that best aligns with the interventional data, according to the statistical distances captured in $\mathbf{D}$.

###### Theorem 1 .

### 3.3 Causal Discovery algorithm

After deriving the differentiable score in the previous section, we now proceed to use the score in a causal discovery algorithm.
Let’s consider a dataset $\mathbf{X}\in\mathbb{R}^{n\times d}$ consisting of $n$ observations of $d$ variables $\{X_{1},X_{2},\dots,X_{n}\}$. We assume the data to be generated from an unknown Structural Equation Model (SEM), which can be described by a Directed Acyclic Graph (DAG) $\mathcal{G}$ representing the causal relationships among variables. Our goal is to recover the causal structure and ordering of the variables from both observational and interventional data. We introduce a potential vector $\mathbf{p}\in\mathbb{R}^{d}$ that induces a permutation representing the causal ordering of variables. Let $S(\mathbf{p})$ be the DiffIntersort score, which measures the consistency of the ordering induced by $\mathbf{p}$ with the interventional data. The causal discovery can then be formulated as a constrained optimization problem

$$ $\min_{\theta,\mathbf{p}}\ \mathcal{L}_{\text{fit}}(\theta,\mathbf{p})\quad \text{subject to}\quad\mathbf{p}\in\arg\min_{\mathbf{p}^{\prime}}S(\mathbf{p}^ {\prime})$ (8) $$

where $\theta$ represents the parameters of the causal mechanisms (e.g., weight matrices in linear models), and $\mathcal{L}_{\text{fit}}(\theta,\mathbf{p})$ is the fitting loss that measures how well the model with parameters $\theta$ explains the observed data. The constraint ensures that the potential vector $\mathbf{p}$ minimizes the DiffIntersort score, thus enforcing a causal ordering consistent with the interventional data. We transform the constrained optimization into an unconstrained penalized problem

$$ $\min_{\theta,\mathbf{p}}\quad\mathcal{L}_{\text{fit}}(\theta,\mathbf{p})+ \lambda S(\mathbf{p}),$ (9) $$

where $\lambda>0$ is a regularization parameter controlling the trade-off between fitting the data and enforcing the causal ordering through the DiffIntersort score.

As an example, a linear causal model can be constructed as

$$ $X_{j}=\sum_{i=1}^{d}W_{ji}X_{i}+b_{j}+N_{j},$ (10) $$

where $W_{ji}$ are the entries of the weight matrix $\mathbf{W}\in\mathbb{R}^{d\times d}$, $b_{j}$ is the bias term, and $N_{j}$ is a noise term. To enforce the causal ordering induced by $\mathbf{p}$, we use the permuted upper-triangular matrix ${\bm{M}}_{\mathbf{p}}={\mathcal{S}}^{T}_{\textrm{bin}}({\bm{p}}{\bm{o}}^{T}/t)
{\bm{L}}{\mathcal{S}}^{T}_{\textrm{bin}}({\bm{p}}{\bm{o}}^{T}/t)^{T}$, which is a $d\times d$ matrix with $d(d-1)/2$ entries equal to $1$. The matrix represents the possible locations of edges in the graph according to the causal ordering $\mathbf{p}$. By element-wise multiplication $\tilde{\mathbf{W}}=\mathbf{W}\circ{\bm{M}}_{\mathbf{p}}^{T}$, matrix ${\bm{M}}_{\mathbf{p}}$ acts as a mask to ensure that variable $X_{j}$ can only depend on variables preceding it in the causal ordering. The predicted values can be written in terms of the entries of $\tilde{\mathbf{W}}$ as $\hat{X_{j}}=\sum_{i=1}^{d}\tilde{\mathbf{W}}_{ji}X_{i}+b_{j}$.

Inspired by the fitting loss in Shen et al. (), we define the fitting loss $\mathcal{L}_{\text{fit}}(\theta)$ as:

$$ $\mathcal{L}_{\text{fit}}(\theta,\mathbf{p})=\frac{1}{n^{0}}\sum_{i=1}^{n^{0}} \ell(\mathbf{x}_{i},\hat{\mathbf{x}}_{i};\theta,\mathbf{p})+\gamma\sum_{e\in \mathcal{E}}\omega^{e}\left(\frac{1}{n^{e}}\sum_{i=1}^{n^{e}}\ell^{e}(\mathbf{ x}_{i},\hat{\mathbf{x}}_{i};\theta,\mathbf{p})-\frac{1}{n^{0}}\sum_{i=1}^{n^{0 }}\ell^{0}(\mathbf{x}_{i},\hat{\mathbf{x}}_{i};\theta,\mathbf{p})\right),$ (11) $$

where $\ell(\mathbf{x}_{i},\hat{\mathbf{x}}_{i};\theta)$ is the mean absolute error (MAE) loss function for observational sample $i$, $\ell^{e}(\mathbf{x}_{i},\hat{\mathbf{x}}_{i};\theta)$ is the loss for samples in environment $e$. In our case, an environment corresponds to an intervention on one variable. $\gamma\geq 0$ is a parameter controlling the emphasis on invariance across environments. We use $\gamma=0.5$. $\omega^{e}$ are weights for each environment. We set $\omega^{e}=\frac{1}{|\mathcal{E}|}$. $\mathcal{E}$ is the set of environments, with $0\in\mathcal{E}$ denoting the reference observational environment. $n^{e}$ is the number of samples in environment $e\in\mathcal{E}$. The loss encourages the model to fit the data in the reference environment while penalizing deviations in performance across different environments, promoting robustness to interventions. This should encourage the weights to corresponds to the true causal weights, as the equivalence between robustness and causality is well established (Meinshausen, ).

Combining the fitting loss and the regularization terms, the final loss function is:

$$ $\begin{split}\mathcal{L}(\theta,\mathbf{p})=&\frac{1}{n^{0}}\sum_{i=1}^{n^{0}} \ell(\mathbf{x}_{i},\hat{\mathbf{x}}_{i};\theta,\mathbf{p})\\ &+\gamma\sum_{e\in\mathcal{E}}\omega^{e}\left(\frac{1}{n^{e}}\sum_{i=1}^{n^{e} }\ell^{e}(\mathbf{x}_{i},\hat{\mathbf{x}}_{i};\theta,\mathbf{p})-\frac{1}{n^{0 }}\sum_{i=1}^{n^{0}}\ell^{0}(\mathbf{x}_{i},\hat{\mathbf{x}}_{i};\theta, \mathbf{p})\right)\\ &+\lambda_{1}\left\|\mathbf{W}\right\|_{1}+\lambda_{2}S(\mathbf{p}).\end{split}$ (12) $$

This loss function includes all the components:
(1) *Data Fitting Loss*Measures how well the model predicts the observed data, adjusted for interventions,
(2) *Environment Invariance Penalty*Encourages the model to have consistent performance across different environments,
(3) *$L_{1}$ Regularization*: Promotes sparsity in the weight matrix $\mathbf{W}$,
(4) *DiffIntersort Regularization*: Incorporates interventional faithfulness by penalizing with the DiffIntersort score $S(\mathbf{p})$ eq. 7
. We also note that no acyclicity constraint is needed as the weight matrix is enforced to be acyclic through the masking based on the causal order $\mathbf{p}$.

## 4 Empirical Results

We next evaluate the proposed DiffIntersort differentiable score both in its effectiveness in deriving the causal order of a system, as well as it usefulness as a differentiable constraint in a causal discovery model.

We first evaluate the DiffIntersort score in it ability to recover the causal order in simulated graphs and distance matrices. We here reproduce the experiment of (Chevalley et al., ). We compare the top order divergence of DiffIntersort to SORTRANKING, and to Intersort for $5$ and $30$ variables, and the upper-bounds of Thm 2 and Thm 4 derived in (Chevalley et al., ). Intersort does not scale beyond $100$ variables. The upper-bounds serve as a sanity check to assess how far to the true optimum of the score the approximate solution is. We evaluate on both Erdős-Rényi distribution (Erdős et al., ) and scale-free network modelled by the Barabasi-Albert distribution Albert & Barabási (), with varying edge densities and intervention coverage. The results are reported in fig. 1 for $2000$ variables and in figs. 5 and 6 for $5$, $30$, $100$ and $1000$ variables. It is crucial that our score is optimizable up to at least $2000$ variables as it is a common scale in real world datasets such as single-cell transcriptomics (Replogle et al., ). As observed, DiffIntersort fulfills the upper-bounds for all settings, even at large scale. At large scale, it also outperforms SORTRANKING in almost all settings. Those results validate our proposed approach of solving the Intersort problem in a continuous and differentiable framework, and guarantees that it is not limited by scale.

Figure: (a) Simulation ER with 2000 variables
Refer to caption: extracted/5919639/figures/violin_bound_2000.png

We now evaluate our method, DiffIntersort, on simulated data and compare its performance to various baseline methods. We follow the experimental setup of Chevalley et al. () to ensure a fair and consistent evaluation across different domains. Specifically, we generate graphs from an Erdős-Rényi distribution (Erdős et al., ) with an expected number of edges per variable $c\in\{1,2\}$. Data is simulated using both linear relationships and random Fourier features (RFF) additive functions to capture non-linear dependencies. In addition to these synthetic datasets, we apply our models to simulated single-cell RNA sequencing data generated using the SERGIO tool (Dibaeinia & Sinha, ), utilizing the code provided by Lorch et al. () (MIT License, v1.0.5). We also test our method on neural network functional data following the setup of Brouillard et al. (), using the implementation from Nazaret et al. () (MIT License, v0.1.0). To assess the impact of interventions, we vary the ratio of intervened variables in the set ${25\%,50\%,75\%,100\%}$. All datasets are standardized based on the mean and variance of the observational data to eliminate the Varsortability artifact identified by Reisach et al. (). For the linear and RFF domains, the noise distribution is chosen uniformly at random from the following options: uniform Gaussian (noise scale independent of the parents), heteroscedastic Gaussian (noise scale functionally dependent on the parents), and Laplace distribution. In the neural network domain, the noise distribution is Gaussian with a fixed variance. We conduct experiments on 10 simulated datasets for each domain and each ratio of intervened variables. The observational datasets contain 5,000 samples, and each intervention dataset comprises 100 samples, mirroring the sample sizes typically found in real single-cell transcriptomics studies (Replogle et al., ).

We compare the performance of DiffIntersort and SORTRANKING (Chevalley et al., ) as measured by the top order divergence $D_{top}$ on $100$ variables in fig. 2. For the DiffIntersort score and the Intersort score, we use the same parameters as in Chevalley et al. (): $\epsilon=0.3$ for linear, RFF and NN data, and $\epsilon=0.5$ for GRN data, and $c=0.5$. We use the Wasserstein distance (Villani et al., ) for the statistical metric. Results for $10$ and $30$ variables, additionally compared to Intersort, can be found in the appendix in fig. 7. As can be observed, the performance of the two algorithms is close. This demonstrates that the optimizing DiffIntersort can be solved at scale using continuously differentiable optimization also on realistic synthetic data.

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

We now evaluate our causal discovery method on synthetic datasets generated using linear structural equation models (SEMs), gene regulatory network (GRN) models, random Fourier features (RFF) models, and neural network (NN) models as presented previously. For each model type, we consider variable sizes of $10$, $30$, and $100$ to assess scalability and performance across different problem dimensions. We use two evaluation metrics: Structural Hamming Distance (SHD) (Tsamardinos et al., ) and Structural Intervention Distance (SID) (Peters & Bühlmann, ) to compare the inferred graphs with the true causal graphs. We compared to two baselines, namely GIES (Hauser & Bühlmann, ) and DCDI (Brouillard et al., ). We note that those two baselines do not scale to $100$ variables. For our model, we compare the performance of our proposed causal discovery model with and without the DiffIntersort constraint. We present the results for the SHD metrics at $30$ variables in fig. 3. The results for $10$ and $30$ variables for SHD can be found in the appendix in fig. 9. The results for SID can be found in fig. 12(l) in the appendix.

As can be seen, the DiffIntersort constraint is consistently beneficial in terms of performance on both metrics, for all types of data and at all considered scales. This comparison validates the usefulness of inducing the interventional faithfulness inductive bias to a causal models via the DiffIntersort score. We expect that this approach may be applicable to other causal tasks of interest, in settings where a large set of single variable interventions are available. Compared to baselines, our model outperforms on the GRN and RFF data. GIES is the best model on linear data, and DCDI has a slightly better performance on NN data. GIES and DCDI do not scale to $100$ variables but we would expect the results to be the same, as our algorithm has an F1 score that is almost unaffected by the number of variables (see fig. 4). The results on the F1 score also shows the robustness of our causal discovery model with the DiffIntersort constraint to the number of variables.

Figure: (a) GRN, 30 vars, SHD
Refer to caption: extracted/5919639/figures/results_test_const_gene_30_SHD_violin_plot.png

Figure: (a) Linear
Refer to caption: x1.png

## 5 Conclusion

In this work, we addressed the scalability and differentiability limitations of Intersort, a score-based method for discovering causal orderings using interventional data. By reformulating the Intersort score through differentiable sorting and ranking techniques—specifically utilizing the Sinkhorn operator—we enabled the scalable and differentiable optimization of causal orderings. This reformulation allows the Intersort score to be integrated as a continuous regularizer in gradient-based learning frameworks, facilitating its use in downstream causal discovery tasks.
Our proposed approach not only preserves the theoretical advantages of Intersort but also significantly improves its practical applicability to large-scale problems. Empirical evaluations demonstrate that incorporating the differentiable Intersort score into a causal discovery algorithm leads to superior performance compared to existing methods, particularly in complex settings involving non-linear relationships and large numbers of variables. The algorithm exhibits robustness across various data distributions and noise types, effectively scaling with increasing data size without compromising performance.
By bridging the gap between interventional faithfulness and differentiable optimization, our work opens new avenues for integrating interventional data into modern causal machine learning pipelines. This advancement holds promise for a wide range of applications dealing with large covariate sets where understanding causal relationships is crucial, such as genomics, neuroscience, environmental, and social sciences.

While our approach enhances the scalability and differentiability of causal discovery using interventional data, several avenues for future research remain. One potential direction is to integrate our differentiable Intersort score into more complex models, such as deep neural networks, to further improve causal discovery in high-dimensional and highly non-linear settings.
Another promising area is the application of the differentiable Intersort approach to real-world datasets in domains like genomics or healthcare, where large-scale interventional data are increasingly available.

## Acknowledgements and disclosure

MC, AM and PS are employees and shareholders of GSK plc.

## 6 Appendix

### 6.1 Proofs

###### Proof of Theorem 1 .

First, let us recall that we have ${\bm{p}}\in\mathbb{R}^{d}$, and $\pi\in\{0,1,\dots,d\}^{d}$, where $\forall i,j\in\{0,1,\dots,d\},\pi_{i}\neq\pi_{j}$. We thus trivially have that any permutation $\pi$ can be represented by a potential ${\bm{p}}$, by ${\bm{p}}_{i}=-\pi_{i}\forall i\in\{0,1,\dots,d\}$.
We now have to prove that if $\pi\in\Pi$, then the corresponding potential ${\bm{p}}_{\pi}\in{\mathbb{P}}$. Let $s=\max_{\pi}S(\pi,\epsilon,D,\mathcal{I},P_{X}^{\mathcal{C},(\emptyset)},
\mathcal{P}_{int},c)$ be the maximum achievable score. The sum of the score is over the elements of ${\mathbf{D}}_{ij}$ where $\pi_{i}<\pi_{j}$. For all these pairs of indices, we also have that $p_{\pi_{i}}>p_{\pi_{j}}$, and thus for all those pairs, we also have $(\operatorname*{Step}(\operatorname*{grad}({\bm{p}}_{\pi})))_{ij}=1$. This exactly corresponds to the elements that are non-zero and thus contribute to the sum in eq. 3. Thus we have that $S({\bm{p}}_{\pi})=s$, and as such ${\bm{p}}_{\pi}\in{\mathbb{P}}$, which concludes the proof.
∎

###### Proof of Theorem 1 .

### 6.2 Additional experiments

Figure: (a) Simulation ER with 5 variables
Refer to caption: extracted/5919639/figures/violin_bound_5.png

Figure: (a) Simulation SF with 5 variables
Refer to caption: extracted/5919639/figures/violin_bound_sf_5.png

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

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

Figure: (a) GRN, 10 vars, SHD
Refer to caption: extracted/5919639/figures/results_test_const_gene_10_SHD_violin_plot.png

Figure: (a) GRN, 10 vars, SHD
Refer to caption: extracted/5919639/figures/results_test_const_sf_gene_10_SHD_violin_plot.png

Figure: (a) GRN, 10 vars, SID
Refer to caption: extracted/5919639/figures/results_test_const_gene_10_SID_violin_plot.png

Figure: (a) GRN, 10 vars, SID
Refer to caption: extracted/5919639/figures/results_test_const_sf_gene_10_SID_violin_plot.png

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 5
presentation: 2
questions: - As mentioned in weaknesses, constrained optimization only appears in Eq. (8), and the experiments and loss function proposed do not solve the constrained problem. Is there a particular reason to consider this constrained problem, instead of focusing instead on regularization from the start?
- What is the role of regularization strength? How is $\lambda_2$ determined, and what is the effect of varying it?
- I also am not sure I quite understand the final model used; a linear model is given in Equation (10), but loss functions are specified in terms of $\ell(\mathbf{x}_i, \hat{\mathbf{x}}_i; \theta)$; is the linear model used in all experiments, with $\theta = \\\{ \mathbf{W}, \mathbf{b} \\\} $?

## References

Reisach, A., Seiler, C., & Weichwald, S. (2021). Beware of the simulated DAG! Causal discovery benchmarks may be easy to game. Advances in Neural Information Processing Systems, 34, 27772-27784.

Reisach, A., Tami, M., Seiler, C., Chambaz, A., & Weichwald, S. (2024). A scale-invariant sorting criterion to find a causal order in additive noise models. Advances in Neural Information Processing Systems, 36.

Chevalley, M., Schwab, P., & Mehrjou, A. (2024). Deriving Causal Order from Single-Variable Interventions: Guarantees & Algorithm. arXiv preprint arXiv:2405.18314.
soundness: 2
strengths: - I find the idea of regularizing based on topological ordering very interesting, especially because this can cause many downstream errors. Using a differentiable version of Intersort as a regularizer seems a fitting and effective way to realize this idea. It is also clever to derive a score function which uses learned causal orders to avoid acyclicity constraints.
- Empirically and theoretically, the results of DiffIntersort seem close to Intersort, which is reassuring for the differentiable relaxation.
summary: This paper focuses on a differentiable version of the Intersort score, a score for identifying the topological ordering of a causal graph via interventional data (Chevalley et al. 2024). This proceeds by rewriting the Intersort score to be parameterized by a potential, and applying tools used for gradient-based optimization of permutations, resulting in DiffIntersort. A differentiable score for topological ordering allows for a novel regularization scheme in score-based causal discovery, and a novel method CausalDisco is presented, which optionally includes DiffIntersort as a regularization term, and is tested empirically against DCDI and GIES. It is shown empirically that DiffIntersort provides scores generally consistent with Intersort with significantly less computational cost. CausalDisco provides mixed results, outperforming benchmark methods in two of four tested settings, but regularization with DiffIntersort is shown to consistently benefit CausalDisco.
weaknesses: **On the Presentation of the Main Ideas** 

I find the novelty in Intersort -> DiffIntersort somewhat limited; I recognize the use of a potential is nontrivial, but this is still mostly just combining Chevalley et al. (2024) and Annadani et al. (2023) in a straightforward way. I find this particularly unfortunate because I think regularization based on ordering to be a much more novel and interesting idea. While subjective, I think the paper would benefit from focusing on order-based regularization, with DiffIntersort as a way to achieve that (rather than focusing on DiffIntersort, with order-based regularization an application thereof).

Note that the effects of this propagate to the experimental section: for example, I do not see any experiments considering the effect of varying $\lambda_2$, which is an interesting ablation from the perspective of regularization. (Indeed, I cannot seem to find the actual value of $\lambda_2$ used -- please let me know if I am missing it). Causal discovery with interventional data has also been an increasingly popular topic, and I imagine DiffIntersort could be used to regularize a more popular objective function.

**On Reproducibility** 

Code is provided to compute DiffIntersort scores, which is great. However, it doesn't seem like there is code available for CausalDisco, or for the experiments included in the paper. Especially since the details in the manuscript regarding CausalDisco are somewhat sparse (e.g.., hyperparmeters used for regaularization), code would be very helpful in ensuring the reproducibility of results.

**On the "DiffIntersort Constraint"** 

I have doubts about the presentation of "CausalDisco with constraint." A constraint only appears in Eq. (8), and it is likely more appropriate to present as a bilevel optimization problem. Moreover, the actual CausalDisco loss is presented as an unconstrained regularization problem, solutions of which are (a priori) not generally equivalent to the corresponding bilevel problem.

**Minor Note on Naming** 

The closest ideas to the proposed CausalDisco are likely those of Reisach et al. (2021) and Reisach et al. (2023). The library developed for their work is called "CausalDisco", which can be a source of possible confusion.

## Minor Typos
- Line (220) Kuhn (1955) should use `\citep{}`
- (Line 316) Albert & Barabasi (2002) should use `\citep{}`
- (Line 408) I think the reference should be to simply Figure 9.

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: See above.
soundness: 2
strengths: 1. The experimental results is relatively comprehensive.

2. The introduction of differentiable operator is novel to me and may be beneficial to other permutation based methods.
summary: The paper considers the problem of identifying the causal ordering from interventional data. The proposed method is based on an existing score-based algorithm named "Intersort". Different from Intersort that requires score evaluation over permutations, this work uses differentiable sorting and ranking techniques, namely the Sinkhorn operator, to enable gradient-based optimization for causal discovery. The proposed approach is claimed to be efficient and scalable.
weaknesses: 1. The motivation of the work is unclear:
   - What is being discovered? In section 3.3 the authors say "Our goal is to recover the causal structure and ordering of the variables from both observational and interventional data", but in the title and elsewhere just the causal ordering is emphasized. So which is exactly the objective?
   - Note that with the interventional data available the whole graph's equivalence class can be identified and the complexity for such whole structure learning is almost the same -- given a causal ordering, the exact DAG from its CPDAG (or smaller equivalence class due to interventions) can be identified trivially, using the local Markov property. Therefore, I doubt that the causal ordering learning in this setting is  not a necessary problem.

2. The identifiable guarantee of the work needs more justification:
   - What are exactly all the assumptions needed in this work? I only see Interventional Faithfulness. Does it mean that this work, together with the score used and the distance function used, is fully nonparametric? No assumptions on e.g., the functional model, the noise distributions (which other differential methods usually require) are required? And what are the hidden assumptions under the distance measure, e.g., is it meaningful to compare the distances?
   - For the identifiability of the causal order itself, what is needed? The authors says they "derive the causal order of the variables when many single-variable interventions are available", but how `many' is many? That a causal ordering can be fixed yields that the whole DAG structure can be identified (from its CPDAG or smaller equivalence class). Then do the authors implicitly assume this?
   - Overall, since Intersort may not appear as a well-known method for now, some background over that algorithm should be given, before providing any details of implementation or efficiency improvements in the paper.

3. The presentation can be improved. E.g., the definition 2 (top order divergence) serves only as a metric of evaluation, and should not be put at the beginning of the method part to distract readers.

4. More detailed literature review on causal discovery from interventional data, and on the differentiable ordering learning is expected.

---

confidence: 3
contribution: 1
initial_rating: 3
presentation: 3
questions: - How does the computational complexity of DiffIntersort compare to existing differentiable approaches for interventional data such as DP-DAG?
- How do you handle $p_i=p_j$? Does that connect the nodes bidirectionally or does it remove the edge altogether?
- What is your wall time when you claim that the baselines do not scale over 100 variables?
soundness: 2
strengths: The development of scalable approaches for causal discovery is valuable and has meaningful applications in important contexts, such as genomics. The proposed algorithm performs better than the other baselines on some scenarios (gene-regulatory and random Fourier features) and on par with them on the remaining (linear SEMs and neural network simulated data).
summary: The authors introduce DiffIntersort as a differentiable approximation of the existing Intersort algorithm (Chevalley et al., 2024). By minimizing a score-function, Intersort first fits from data a complete order between variables that is then refined by exploring the space of permutations. In this paper, the authors define the score as a differentiable function using the Sinkhorn operator and then propose to jointly minimize this score and a loss function on the observations. DiffIntersort can handle interventional samples, and it is thus compared with GIES and DCDI on synthetic data and simulated but plausible gene regulatory networks.
weaknesses: The paper does not reference in the related works significant works of recent literature. First, the decomposition of the matrix to ensure acyclicity without constraints (Eq. 3) strongly resembles a parameterization introduced by NO-CURL using a ReLU function (Yu et al. 2021) and by COSMO using a smooth acyclic orientation (Massidda et al. 2024), both in the observational context. Overall, the algorithm has even more similarities with BayesDAG (Annadani et al. 2023) and DP-DAG (Charpentier et al. 2022). In particular, DP-DAG also supports interventional data and uses the Gumbel-Sinkhorn operator and the Hungarian algorithm to get a differentiable representation of a permutation. Overall, the paper should have a more comprehensive related work section and compare experimentally with other differentiable approaches for interventional data.

The authors motivate their differentiable extension of Intersort to avoid its substantial computational cost. It would aid the presentation if the authors could report it and, more importantly, it could improve the contribution if the authors actually reported what the computational complexity of DiffIntersort is. On this matter, it's worth recalling that DP-DAG has quadratic complexity on the number of nodes. Furthermore, I would expect in a paper focused on the *efficiency* of the algorithm to have at least one comparison of the execution time of the algorithm and its competitors.

Also, please note that your solution is incorrectly, I guess, referred to as "causaldisco" in Figure 3.

*References* in alphabetical order.
- Annadani, Yashas, et al. "Bayesdag: Gradient-based posterior inference for causal discovery." Advances in Neural Information Processing Systems 36 (2023): 1738-1763.
- Charpentier, Bertrand, Simon Kibler, and Stephan Günnemann. "Differentiable DAG Sampling." International Conference on Learning Representations.
- Massidda, Riccardo, et al. "Constraint-Free Structure Learning with Smooth Acyclic Orientations." The Twelfth International Conference on Learning Representations.
- Yu, Yue, et al. "DAGs with no curl: An efficient DAG structure learning approach." International Conference on Machine Learning. Pmlr, 2021.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 1
questions: ## Main ##
As recently discussed by Ng et al., 2024 (https://proceedings.mlr.press/v236/ng24a/ng24a.pdf), non-convexity of the loss might constitute the biggest challenge to differentiable causal discovery methods. What are the authors' thoughts on this? This should be a subject of discussion in the paper.

## Others ##
1. Do the authors assume causal sufficiency? Do they assume independence of the noise terms? 
2. The experiments of Figure 1 are carried out on very sparse graphs, in the ER setting: specifically, authors sample graphs with 2000 nodes, and the expected number of edges {200, 100, ~40} corresponding to $p_e$ = {1e-4, 5e-5, 2e-5}. I would ask to consider experiments in much denser settings.
soundness: 2
strengths: The problem of differentiable causal discovery is of great interest given that, as mentioned by the authors, causal discovery algorithms that can be trained by backpropagation can be used as components of deep learning pipelines. The main contribution of this work is the definition of an alternative score (alternative to that of Intersort) that is differentiable and ensures recovery of an optimal causal order when the score is maximized. Empirically, Diffintersort appears to be scaling to large graphs, with impressive performance when the number of interventions increases. I believe this algorithm can make a good contribution to the causal discovery community
summary: In this paper, the authors propose a modification of the Intersort algorithm from Chevalley et al., 2024. Intersort adopts a scoring function, that, under interventional faithfulness assumption and access to single variable interventions for each node, if maximized yields the correct causal order. The scoring function measures the fit of the casual order based on its compatibility with interventional data. The downside of this method is that it does not scale to graphs with many nodes (up to 100, according to the authors): this work introduces Diffintersort as an alternative that is differentiable and proposes to scale to thousands of nodes.
weaknesses: ## Main points

1. One main issue I found is that even though *efficiency* is in the title, its empirical and theoretical analysis is not treated in the paper. The novel contribution of this work consists of the algorithmic part, which is a variation of Intersort (as the theory heavily builds on Chevalley, 2024, and Annadani, 2023): however, for this to be considered a sufficiently relevant contribution, simply stating that *Intersort does not scale beyond 100 variables* (L314), whereas DiffIntersort can scale to 2000 variables, is not enough. This must be theoretically and empirically analyzed. In detail:
    1. An analysis of the algorithmic complexity of DiffIntersort is missing, and this should be done in comparison to an analysis of the algorithmic complexity of Intersort
    2. There is no trace of the execution times of DiffIntersort in the experiments, under different conditions (especially varying the number of nodes, as I believe this is what the authors mean by *efficiency*), and in comparison with the execution times of Intersort. 
2. The experimental section has some obscure points and limitations in its reproducibility
    1. It is not mentioned what optimization algorithm is used for the experiments. 
    2. It is unclear how hyperparameters are set. The authors simply say that they choose $\epsilon = 0.3$, $\gamma = 0.5$, there is no mention of how $\lambda_1, \lambda_2$ parameters of the loss are set, and the same goes for $c$.
    3. Some of the proposed mechanisms are unclear or need more discussion. The authors say that they use neural networks to generate nonlinear mechanisms as in Brouillard et al., 2020: in that paper, the authors define different classes of structural causal models (with and without additive noise) using neural networks. Which one are they referring to? Moreover, the discussion of random Fourier features generation is overlooked: I believe it would be beneficial to provide some reference to the topic and ideally some discussion in the appendix. Finally, there is no mention of how weights are sampled for the linear mechanisms. In general, please discuss in much greater detail the data generation.

## Other points ##
1. In the nonlinear case (potentially with non-additive noise) it is not discussed how $\hat X_i$ in the loss formulation is computed. The authors detail computations for the linear case as an example, using masking on the weight matrices that ensures acyclicity, yet it is unclear how this is achieved in the nonlinear case.
2. The authors claim that they prove DiffIntersort inference abilities on a variety of noise terms, yet they only adopt Gaussian and Laplacian distribution. This is quite a limited choice: I suggest adding more distributions, or you could also consider densities as nonlinear transformations parametrized by neural networks as in Montagna et al., 2023 (https://arxiv.org/abs/2310.13387)

Meta-review: This paper introduces a differentiable version of the Intersort algorithm, and proposes CausalDisco, a causal discovery method using DiffIntersort for regularization. Reviewers found the scalability and differentiability promising but raised concerns about limited theoretical and empirical analysis, unclear objectives, incomplete implementation details, and missing comparisons with relevant methods. An ethics concern was also flagged regarding potential improper paraphrasing. While promising, the paper requires significant revisions to clarify its contributions and address these issues.

Decision: Reject