Large-Scale Differentiable Causal Discovery of Factor Graphs
Authors: Romain Lopez
Jan-Christian Huetter
Jonathan Pritchard
Aviv Regev
Venue: NeurIPS 2022 Accept (2022)
Keywords: Causal discovery learning
causal learning gene regulatory networks
factor graphs
graph learning from interventions

Abstract: A common theme in causal inference is learning causal relationships between observed variables, also known as causal discovery. This is usually a daunting task, given the large number of candidate causal graphs and the combinatorial nature of the search space. Perhaps for this reason, most research has so far focused on relatively small causal graphs, with up to hundreds of nodes. However, recent advances in fields like biology enable generating experimental data sets with thousands of interventions followed by rich profiling of thousands of variables, raising the opportunity and urgent need for large causal graph models.  Here, we introduce the notion of factor directed acyclic graphs ($f$-DAGs) as a way to restrict the search space to non-linear low-rank causal interaction models. Combining this novel structural assumption with recent advances that bridge the gap between causal discovery and continuous optimization, we achieve causal discovery on thousands of variables. Additionally, as a model for the impact of statistical noise on this estimation procedure, we study a model of edge perturbations of the $f$-DAG skeleton based on random graphs and quantify the effect of such perturbations on the $f$-DAG rank. This theoretical analysis suggests that the set of candidate $f$-DAGs is much smaller than the whole DAG space and thus may be more suitable as a search space in the high-dimensional regime where the underlying skeleton is hard to assess. We propose Differentiable Causal Discovery of Factor Graphs (DCD-FG), a scalable implementation of $f$-DAG constrained causal discovery for high-dimensional interventional data. DCD-FG uses a Gaussian non-linear low-rank structural equation model and shows significant improvements compared to state-of-the-art methods in both simulations as well as a recent large-scale single-cell RNA sequencing data set with hundreds of genetic interventions.

Large-Scale Differentiable Causal Discovery of Factor Graphs
Romain Lopez 1.2 Jan-Christian Hütter Jonathan K. Pritchard 2.3.t Aviv Regev
Division of Research and Early Development, Genentech {lopez. romain, huettej1, rggeve)@gene. com
Department of Geneties, Stanford University pritch@stanford.co edu
Department of Biology, Stanford University
These authors contributed equally to this work.
Abstract
A common theme causal inference learning causal relationships between observed variables, also known as causal discovery. This is usually daunting task, given the large number of candidate causal graphs and the combinatorial nature of the search space. Perhaps for this reason, most research has far focused on relatively small causal graphs, with up to hundreds of nodes. However, recent advances in fields like biology enable generating experimental data sets with thousands of interventions followed by rich profiling of thousands of variables, raising the opportunity and urgent need for large causal graph models. Here, we introduce the notion factor directed acyclic graphs (J-DAGs) as a way to restrict the search space to non-linear low-rank causal interaction models. Combining this novel structural assumption with recent advances that bridge the gap between causal discovery and continuous optimization, we achieve causal discovery on thousands of variables. Additionally, on this estimation procedure, we study a model of edge perturbations of the f-DAG skeleton based on random graphs and quantify the effect of such perturbations on the f-DAG rank. This theoretical analysis suggests that the set of candidate f-DAGs is much smaller than the whole DAG space and thus may be more suitable as search space the high-dimensional regime where the underlying skeleton is hard to assess. We propose Differentiable Causal Discovery of Factor Graphs (DCD-FG), a scalable implementation f-DAG constrained causal discovery for high-dimensiioal interventional data. DCD-FG uses Gaussian non-linear low-rank structural equation model and shows significant improvements compared fatt--oftthe-art methods in both simulations well recent large-scale single-cell RNA sequencing data set with hundreds of genetic interventions.
model for the impact of statistical noise
Introduction
Characterizing causal dependencies between variables is fundamental problem in science [1, 2, 3]. Such relationships are often described via directed acyelic graph (DAG) where nodes represent variables and directed edges encode causal links between pairs of variables. One may then consider
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

Table 1: Comparison with related work. denotes the number of features, denotes the number of learned factors. Because m remains low (e.g., 20), we have m < d. 'additive cubic model only.
non-linear link function
likelihood eval.
DAG penalty eval.
intervention model
Related work
complexity
complexity
NO-TEARS [19]
O
(
X
NO-BEARS [26]
(
O(d
X
NO-TEARS-LR [20]
X
((md)
X
DCDI [6]
(
DCD-FG (this work)
O(nd)
(mnd)
set of conditional probability distributions in addition to the DAG to define a causal graphical model (CGM) [4]. The graphical model is referred as "causal" because it provides clear semantic for defining the effect of interventions on the joint likelihood [5].
The combinatorial nature of the set of DAGs and its size make inference of causal structure from data a hard problem. Indeed, the number of DAGs (or the number of permutations) grows super exponentially in the number of nodes, limiting the practical application of most methods to tens of nodes [6, 7, 8]. This highlights a disparity between the development of the field and emerging data from real-world problems in finance [9, 10], neuroscience [11, 12] and high-throughput biology [13, 14, 15, 16]. For example, experimental advances in biology have enabled the generation of data sets where expression profiles of thousands of RNA transcripts are measured in hundreds of thousands of cells following genetic interventions in each of hundreds or thousands of different genes.
To reduce this computational burden, promising direction is to limit the complexity of the search space. This idea has become commonplace the interface of statistics and optimization, where sparsity [17] as well as low-rank assumptions [18] are often exploited in machine learning algorithms. Although sparsity is widely used causal structure learning [7, 19], the use of low-rank constraints in this specific setting remains largely under-explored, with few notable exceptions [20, 21J. This may be due to the fact that such an assumption may be hard to incorporate into common methods for permutations.
DAG learning that rely graphical constraints
Recently, the NOTEARS methodology [19] introduced a continuous relaxation of DAG learning, effectively closing the gap between causal structure learning and continuous optimization. NOTEARS facilitates incorporating more flexible structural assumptions into DAG learning, such as neural network parametrizations of the conditional probability distributions. Although NOTEARS has been the subject of several follow-up works [22, 23, 24, 25, 26], including learning from interventional data [6], several challenges remain. First, the likelihood function usually decomposes as a product of local likelihoods for each node, and existing options for those local models are either fast fit but likely underfitting the data (linear models [19]) or computationally expensive and prone to overfitting (a neural network taking as input all of the other nodes [6]). Second, most methods rely on differentiable penalty for acyclicity whose evaluation has cubic complexity in the number of nodes. As a result, these methods are impractical for graphs with more than a hundred nodes.
In this work, we investigate both challenges and propose methodology for large-scale discovery of causal structure and prediction of unseen interventions that scales to millions of samples and thousands of nodes. Our key idea is to limit the search space to what we refer to factor directed acyclic graphs (f-DAGs), a class of low-rank graphs defined in Section 1. This constraint assumes that many nodes share similar sets of parents and children, which is the case in scale-free topologies [20] and many biological systems, where genes act together in programs [3, 27, 28]. Based on this class of graphs, we design flexible model and scalable inference procedure (Table 1) that we refer to as Differentiable Causal Discovery of Factor Graphs (DCD-FG).
After introducing the necessary background (Section 2), we define the class of f-DAGs (Section 3) and draw connections to several flavors of matrix factorization. particular, we show connections between the number of factors in an f-DAG and the Boolean rank [29] of its adjacency matrix. We exploit these connections to present a scalable acyclicity score with linear complexity in the number of observed variables. Finally, we characterize the identifiability of these graphs under an Erdös- Rényi random graph model and prove the instability of the Boolean rank under edge perturbations. The latter analysis highlights that restricting inference to f-DAGs is more efficient in the high-

---

dimensional causal discovery regime. Then, we posit a flexible class of likelihood models as well as scalable inference method for our DCD-FG framework (Section 4). Finally, we present runtime experiments, simulation studies, and case study on single-cell RNA sequencing data with hundreds of genetic interventions (Section 5). In this last challenging instance of interventional data with high-dimensional measurements, we show that our framework outperforms current state-of-the-ar causal discovery approaches for predicting the effect held-out interventions.
Background
Our work builds upon continuous relaxations of the causal structure learning problem. We therefore first briefly introduce causal graphical models and then how the inference problem is solved with a gradient-based optimization framework.
2.1 Causal Graphical Models
Following the framework introduced in [5], let X = (X1,..., Xa) denote a set of random variables with joint probability distribution P. Let G E) be DAG where each vertex v; V is associated with a random variable X, and each edge (v,, Uj) represents a causal relationship from X; to Xj. A correspondence between the graph and the probability distribution /' is obtained via the factorization condition
P(X1 X.) II (x*)).
where (i) denotes the vector of random variables formed by all parents of node vj in G. We refer to pair (P. G) that satisfies the factorization condition as a causal graphical model (CGM).
Unlike classical graphical models, CGMs provide principled way to model interventions [30]. With interventional data, each observation is measured under a specific regime & with interventional joint density p(k)_ For each regime, a subset the target variables IŞ is subject to intervention. Each of these interventions affects the relationship between the target node and its parents, by altering the conditional distribution. In the case of perfect interventions [8], the dependency intervened upon nodes from their parents is removed, and the interventional joint density p(k) becomes
p(t) = P (a())
ΠΙ p()) (X,),
(2)
JEIz
jeIn
where p(k) models the effect of the perfect (stochastic) interventions on each targeted feature of IR-
2.2 Differentiable Causal Structure Learning
significant challenge is to identify and learn CGMs from data when the causal relationships between features not known beforehand. Several methods have been proposed for this problem, such as cesttraint-bssed methods (c.g., the PC algorithm [8]), score-based methods (c.g., GSP [31]), and their extensions for modeling interventions (e.g., IGSP [7)), all reviewed in [32].
The NOTEARS methodology [19] proposes to solve a continuous relaxation of the optimization problem of score-based methods, which in the case of observational data can be written as
n Σ -''''-xx' / such that Trexp W o W) =d W
min
(3)
where n denotes the total number of i.idd observations, X Rdd denotes the i-th observation of X, denotes the Hadamard product, W igd denotes the parameters of a lincar Gaussian conditional distribution, assuming equal variance for all nodes, and À > 0 denotes a tuning parameter. The search space is restricted to DAGs by enforcing that the trace of the exponential (tr-exp) of W W is equal to its dimension d. Importantly, evaluating the tr-exp of W o W as well its gradient with respect to W has ((d)) time complexity, which makes it potentially prohibitive for large-scale applications.
This fundamental work has been subject to multiple extensions, of which the most relevant ones are listed in Table 1. First, different models of relationships between variables were introduced with

---

11
U/2
D3
11
3/
fi
f2
f3
U7
f3
Ut
(a) Complete graph Gy
(b) Half Half-square []]
(c) Half square G= G3[V]
Figure I: A factor graph and its induced half-square graphs. Red circles represent variable nodes. Orange diamonds represent factor nodes.
the aim of better fitting the data, including several non-linear models, such as neural networks [6, 22, 24, 33, 23]. Because these methods use one neural network per conditional distribution (or one deep sigmoidal flow), they computationally expensive, emphasized in [6]. The additive cubic model proposed by NOBEARS [26] does not have these scalability issues, but may be of limited flexibility. To reduce the number of parameters, graph neural networks have been proposed [25], but that architecture cannot be easily applied to the case of interventional data. In addition, NOTEARS- LR [20] proposes a low-rank decomposition of the linear model from NOTEARS but did not consider that the computational complexity could become a bottleneck for large graphs
Second, studies also investigated variants of the acyclicity penalty. [25] proposes matrix power variant the trace exponential criterion, for numerical stability. NOBEARS [26] uses an algebraic characterization of acyclicity based on the spectral radius of the adjacency matrix, that can be approximated in (dd). LoRAM [21] exploits a low-rank assumption to obtain a (dd²) acyclicity score, where denotes the rank. However, the proposed framework is tailored to projections of DAG into a low-rank space but not to causal discovery learning.
Factor Directed Acyclic Graphs
Next, we introduce factor directed acyclic graphs (f-DAGs) and draw connections low-rank matrix factorizations. Then, we motivate their use in causal inference by studying the effect of edge perturbation on their rank. Additionally, we provide two differentiable acyclicity penalties computable in linear time in the number of nodes for a fixed set of factors by exploiting the low-rank structure.
3.1 Definitions and relationship to low-rank decomposition
Factor graphs, especially undirected ones, are commonly used in the graphical models literature to factorize probability distributions and describe the complexity of message passing algorithms [34]. In this work, we use them to construct causal graphs over features. In addition to the set of feature vertices V = (v1,--. (Vi--..,aa) we consider set of factor vertices F {S1,..., .......,,ee Sm ). for m N. When the edge set E links only vertices of different types, the graph G, = (V, F. E) is bipartite graph factor directed graph (f- (f-DiGraph). G, additionally does not contain cycles, we call it an f-DAG, G5 canonically induces two half-square graphs G?LV] and GGLF]. defined by drawing an edge between vertices of V (or F) of distance exactly two in Gy (see example in Figure 1). Unless otherwise mentioned, we will refer to G GJ[V] the half-square of G2 over vertices. We define the set of half half-square graphs of f-DiGraphs with d variables and m factors as
and we refer
GT = (G (V G2 (V, F, F)) is factor directed graph and |F| m)
(4)
Intriguingly, the set 9TH may be identified as the set of matrices with Boolean rank [29] of at most m. Proposition 1 (Bounded Boolean rank of half half-square adjacency matrix). For a factor graph G = (V. E), let E {0, 1)dxm (resp. V {0, 1 mxx be the binary matrix encoding the presence or absence of edges directed towards factor nodes (resp. variable nodes), according to E. The adjacency matrix A(G) of the half-square graph G may be decomposed as A(G) = = V where e [d]. Consequently, A(G) denotes the Boolean matrix product, (U V)ú V Us A Vaj. i,j je

---

has Boolean rank bounded above by m. Conversely, every adjacency matrix over the feature nodes with Boolean rank bounded by m can be written as half-square of an f f-DiGraph with m factors.
This result, proven in Appendix A. establishes a connection between inference of causal f-DAGs and Boolean matrix factorization. Additionally, it is easy to notice that the (integer-valued) matrix product UV provides a valid (weighted) adjacency matrix for G. In this work, we show that the Boolean decomposition is suited for theoretical analysis of the method, while the linear one is useful for efficient algorithm design. We further note that this proposed class of graphs is smaller than the one arising from adjacency matrices with (unconstrained) linear matrix rank bounded by m, which we detail further in Appendix A.4.
3.2 Statistical Properties of Random Causal Factor Graphs
An important theoretical question pertains to the assumptions necessary for identification of acyclic graphs in -T (which we denote by D(")) from data. Under the classical set of assumptions (causal sufficiency, causal Markov property and faithfulness), we may identify the causal DAG from ob- servational data only up to its Markov equivalence class (MEC) [35]. However, graphs in D can have many v-structures (emerging from having several feature parents), potentially making them identifiable. Indeed, we prove in Appendix that under an adapted Erdós Kénés- random graph model [36] over Dj graphs are identifiable with high probability (i.e., their MEC is reduced to one graph). We also verified this with simulations (Figure 2A).
Although valid, the previous result may be disconnected from real-world applications in high- dimensional regimes for which stronger assumptions are required (e.g., strong faithfulness [37J), but rarely hold [38]. More concretely, we expect errors in the estimated skeleton of the feature node graph. As toy model of these errors, we assume the true causal graph G is in G and apply stochastic edge perturbation operator A that randomly removes or adds one edge to G. While there are no general rules as to how identifiability is changed under these perturbations [39] (i.e., the size of the MEC could increase or decrease), we show that the Boolean rank of the graph strictly increases with high probability:
Theorem (Boolean rank instability for edge perturbation). Let G e 9 be sampled according to an Erdős-Rény random directed graph model. The probability that adding or removing an edge increases the Boolean rank is arbitrarily highfor large d:
P (Ranky (A(G)) > Rankg (G)) 2 1 aa""
(5)
where Rankg denotes the Boolean rank, and a > and q (0, depend on m and the parameters of the random graph model.
Precise definitions of the random graph model and proof appear in Appendix C. This result, analogous the upper semicontinuity of the matrix rank, suggests that in noisy settings, the skeleton of graphs in D may be more easily recoverable compared to arbitrary DAGS. We verified those results with simulations (Figure 2B) for small graphs using a mixed integer linear programming approach [40]. We performed surger-scale analysis by randomly perturbing fraction g of the edges in the graph and reporting the resulting matrix rank of the adjacency matrix for larger graphs (Figure 2C).
3.3 Characterizing Acyclicity of Factor Graphs
We start by relating the acyclicity of a f-DiGraph with the one of its induced half-squares. A simple graphical argument is enough show that acyclicity for a a f-DiGraph need only be enforced on the smaller of its half-square graphs.
Lemma 1 (Induced acyclicity).
be an f-DiGraph. Then, G is acyclic + G G2[ is acyclic G?. is acyclic.
Let G=
(6)
The proof can be found in Appendix D. We further note that the matrix UV (VU, resp.) counts the number of paths between two nodes in V (F, resp.) and is thus a valid (weighted) adjacency matrix
Using the matrix rank is only heuristic approximation of the quantitative increase in complexity of the underlying matrix, because in general provides neither upper nor lower bound on the Boolean rank [41].
5

---

MEC (m= 0) MEC (m=15)
1,0
0.05
qu0
9=0.01
Boolean Rank Increase
0.0
0.04
qw0.05
q=0.1
samples (m=5)
q-0.16
Size
10"
LogReg (mn5)
Density
0.03
qu0.2
samples (m=10)
0.4
LogRigg (m=10)
0.02
10
0.2
0.0
A.
Number of nodes
B.
Number of nodes
C.
Linear Rank (d=500; m=40)
Figure 2: Properties of simulated Erdós- Rénéí random f-DiGraphs. (A) Size of MEC of the simulated half-square graph with varying m and d. (B) Probability increase of Boolean rank after edge perturbation (Theorem 1). A single point denotes the result of experiment, the solid line is the probability estimated by logistic regression. (C) Matrix rank change after multiple random edge perturbations.
Thexp
Penalty gradlent runtime (so1)
16!
Power
Spectral Radius
Likelthordd gradient runtiime (w/1Q)
10"
NOTEARS-LR
NOBEARS
10a
Spectral Factor
Theip
Total runtime (ailt)
DCD-FG
Lineardensee LincarLR
DCD-FG (Unear)
10-*
10
MLP Factor
10%
10-3
10a
10
10
10
10"
10a
10"
10*
A.
Number nodes
B.
Number nodes
C.
Number nodes
Figure 3: Runtime analysis. Time for gradient calculation of likelihood (A), penalty (B) their (C) for different variants of optimization- based DAG inference method (on one NVIDIA Tesla T4 GPU with 15Gb of RAM).. We selected batch size of 128 datapoints and number of m 40 factors for all experiments. If value is non-reported, a memory error was raised at runtime. NOTEARS and NOTEARS-LR have almost identical runtime in this analysis.
for G (G²LF) resp.). As result, we may characterize acyclicity by applying the tr-exp penalty to the matrix VU in time O(m² + m² compared the O(d mdd) steps needed for evaluating the penalty on UV. Alternatively, we may use the spectral radius of the adjacency matrix as an acyclicity score [26] that
be approximately computed in O(Tmd) steps, with T iterations of the power method between each gradient step (details in Appendix D). The resulting computational gains are showcased in Figure 3A. For small m. both variants have similar runtime.
4 Differentiable Discovery of Causal Factor Graphs
To define a CGM, we couple the f-DiGraph with a likelihood model. We follow recent work that partitions the parameter space into conditional distribution parameters O and parameters & encoding the causal graph [6]. In particular, let us assume we have at our disposal parameterized distribution M(4) [U(#). [(Φ) over adjacency matrices of f-DiGraphs. The score function, assuming perfect interventions, is defined as:
S( θ) = EM'MM()
Σ E
p(k)
x.;) - [A]][M
(7)
jeZz
where ((t) denotes the distribution of data points X under regime k. w/ denotes a density model for data
feature X,, conditioned on all other features X-, that parents of the feature according to the

---

sampled matrix M. We optimize the score function S under an acyclicity constraint,
max S (b, ①) such that C(E[M( - 0. 0.6
(8)
where C may correspond to either the spectral radius, or the tr-exp characterization of acyclicity. Numerically, first-order optimization techniques with reparameterized gradients and the augmented Lagrangian method are used solve problem (8). We outline here some key features of DCD-FG, and provide the complete implementation details in Appendix E.
Differentiable Sampling of Factor Graphs A first important challenge specific to work constructing a density M(Φ) over f-DiGraphs. The DCDI framework [6], and a few earlier methods [23, 42], parameterize the set of adjacency matrices with entry-wise Gumbel-sigmoid [43] samples, and zeros in the diagonal entries. Naively applying this parameterization for sampling matrices M [U, V] causes the induced feature graphs to have a large number of self-loops, i.e.,
edges of the form (v, v) that we found to be detrimental to the performance of the model. To
circumvent this issue, we propose an alternative model which the matrices U and V are correlated. More precisely, for W e {0,-1, 1} dddm sampled according to a Gumbel-softmax distribution [43], the entries of U and V are constructed from W U.j (( ij 1} and V = 1[ -11 for i E [d],.] j [m]. Because the entries Uis and V5s may never be both equal to 1, there are no self-loops in the induced half-square graph.
hybrid likelihood model A second important challenge is to propose flexible density models pe that have reasonable runtime as well enough capacity for practical purposes. For this, we further exploit the semantics of the factor graph by introducing deterministic factor variables hy at each factor node f & F. These variables are calculated the output of multi-layer perceptron on the input variables of each factor defined by the matrix U, h MLP(U.,f X; ,)) for neural networks parameters y Then, the conditional distribution of each node depends linearly on its parent factors defined by the matrix V, X, ~ Normal (axJ (V.,, h) + (y: for parameters aj RT", ,3 e R. J > O. The resulting computational gains are highlighted in Figure 3B,C. Although we present here the specific case of a Gaussian likelihood model, the same strategy may be adopted for more complex distributional models such as sigmoidal flows [6].
Experiments
We tested DCD-FG both synthetic and real-world data sets with d = 100 to d = 1,000, and a large number observations (n 2 50,000). In this large-scale setting, many state-of-the-art causal discovery methods fail to terminate (DCDI and IGSP). Therefore,
compared DCD-FG to NOTEARS [19], its additive non-linear variant NOBEARS [26] and its linear low-rank variant NOTEARS-LR [20]. In order to have baseline that is external to the NOTEARS framework, applied IGSP (after feature aggregation via clustering with different resolutions when necessary for a reasonable runtime). For every model, performed hyperparameter search using goodness of fit metric on a small validation set. We provide further details on all experiments, including the grids used for hyperparameter search, as well supplementary experiments, in Appendix F.
5.1 Gaussian Structural Causal Models
We consider synthetic data sets with perfect interventions and known targets. Each data set has d= 100 nodes and n 50,000 observations, sampled from interventional distributions governed by either a linear causal mechanism [44] or a nonlinear causal mechanism with additive noise (NN) [45]. Graphs sampled from an Ardōs-Réénj random directed graph model with 10 factors. A total of 100 interventions were performed, each sampling up to 3 target nodes. Datapoints from 20 interventional regimes were held-out for model evaluation.
We assessed the performance of each method by the negative log ligeliikeihood of datapoints from unseen interventions (Interventional NLL) [46], as well as two metrics comparing the estimated graph to the ground truth graphs: the structural Hamming distance (SHD), and the precision and recall of edge detection. We report results for 10 randomly generated graphs and data sets in Figure 4. NOTEARS-LE and DCD-FG both outperformed all methods for all metrics on the linear dataset (p 0.01, Wilcoxon signed-rank test), showing that both effectively exploit the low-rank structure
7

---

Linear (d=100; m=10)
1.00
DCD-FG
NOTEARS-LR
+|I
0.75
Recall
NOTEARS
0.50
NOBEARS
0.25
0.00
0.5
0.6
0.7
500
1000
0.0
0.5
1.0
I-NIL
SHD
Precision
NN (d=100; m=10)
DCD-FG
NOTEARS-LR
Rocall
NOTEARS
NOBEARS
IGSP
0.0
600
Interventional NLL
Precision
Figure 4: Results on simulated Gaussian structural causal models (Section 5.1).
of the causal graph. On the non-linear (NN) dataset, DCD-FG outperformed all methods in terms of interventional NLL, as well as recall and F1 -score (combining precision and recall, Appendix F) (p 0.01, Wilcoxon signed-rank test). DCD-FG has high SHD but we attribute this the fact that all methods besides DCD-FG discover very sparse graphs.
5.2 Genetic Interventions and Gene Expression Data
As an application to real world data, we present an experiment focused on causal learning of gene regulatory networks from gene expression data with genetic interventions, central problem in modern molecular biology [3, 47,48) 48]. Although this particular task has been studied by computational biologists for over two decades, there have been substantial experimental advances in the last few years. In particular, method called Perturb-SS -earrb--equuue now allows us to perform interventions targeting hundreds or thousands of genes and measure the effect on full gene expression profiles hundreds of thousands of single cells using single cell RNA-seq [13]. Surprisingly, however, little to no causal learning work has focused on these advanced datasets. few notable exceptions, such as [35], focused on early data for which the benchmarked methods were tractable (d 24 genes).
We focus on a recent Berturb-CITE-see experiment [14] that contains expression profiles from 218,331 melanoma (cancer) cells, after interventions targeting each of 249 genes. Each measurement from single-cell combines the identity of the intervention (target gene) and a count vector where each entry is the expression level of each gene in the genome. Because of experimental limitations [49], we observe signal only for subset of several thousand genes (here we selected 1,000 genes) out of the approximately 20,000 genes in the genome. This dataset includes patient-derived melanoma cells with same genetic interventions but exposed three conditions: co-culture with cells derived from the patient's tumor (73, 114 cells) (which can recognize and kill melanoma cells), interferon (IFN)~> treatment (87,590 cells) and control (57,627 cells) that we treat as three separate datasets. The goal to identify gene networks in the melanoma cancer cells that either confer
of the experiment resistance or sensitivity to T cell mediated killing, to identify targets for therapeutic intervention in cancer. For every dataset, we retain cells from 20% of the interventions a test set unavailable during training.
We applied our baseline methods as well as DCD-FG to each of the three datasets. Because we do not have a ground truth causal graph, we use datapoints from held-out interventions to evaluate the models [46], reporting both the interventional NLL (I-NLL) and the mean absolute error (I- MAB) across those interventions (Figure 5). Note that accurately predicting the outcome of genetic interventions that were not measured experimentally is of high utility biologists.

---

Control
IFN
Co-culture
NOTEARS-LR
-II
NOTEARS-
1.1
1.2
0.7
0.8
1.2
1.3
0.7
0.8
1.2
1.3
0.70 0.75
I-NLL
I-NLL
I-MAE
I-NLL
I-MAE
Figure Results on the Perturb-CITE- -iT dataset [14] (Section 5.2, lower is better).
DCD-FG outperformed all variants of the NOTEARS method by a large margin, including NOTEARS- IR, and for all metrics (p 0.01, Wilcoxon signed-rank test). In order to diagnose the poor perfor- mance of the competing methods, we looked at the number of inferred edges by each method. All variants of NOTEARS identified extremely sparse graphs (less than hundred edges for NOTEARS and NOBEARS, and few thousand edges for NOTEARS-LRY which may explain their inability to predict the effect of held-out interventions. Interestingly, IGSP identified hundreds of thousands of edges, a number that was comparable to DCD-FG, but still had poor performance. This suggests that the IGSP-inferres graph did not recapitulate well the causal relationships between genes.
In particular, we carefully examined the S-DAG G obtained with the best performing model from our hyperparameter sweep the IFN-r treated cells. That model has half-square [V] has 196,303 edges. On average, each module has 194 ingoing edges and 116
20 factors, and the
outgoing edges. To facilitate visualization, we display the half-square over factors [F[ in Figure 6.
Antigan processing g and presentation
Cell eycle
Chemasins sigraing pailway
JAK2
CDKS
Imarlienn signaling
Imtoricusin 1 regulation extracol matrix
Interlesicin22 signaling pallway
feepproe elecaled sletwie cyloentc créckum
TCF-bata regulatoon of estrscelullr manis
STATI
TNF-aipha effects
actiity
CDKS Upstream gene (JNCSTAT & CDKAS)
HLA-A HLA-C
IRF1
Downstraam gone (CR pathway)
Figure 6: Half-square [F] (completed with a few genes) of an f-DAG identified by DCD-FG on IFN-y treated malignant cells with interventions. Circles: factors, colored by the gene set enrichment analysis result for its incoming and outgoing genes. Empty nodes: genes, labelled name. Red/blue: genes expected to be up or downstream (resp.) by prior biological knowledge.
[F]
To begin to assess the biological relevance of the graph, we performed two analyses. First, we tested the list of incoming and outgoing genes for each factor for enrichment in genes from known biological processes (via [50]). The top hits (node colors in Figure 6) captured many relevant processes in perturbed malignant cells treated by IFN including antigen presentation (needed for recognition by the T cells), multiple innate immune and related signaling pathways (chemokine, interferon, TNF-a, and TGF-8 signaling: all can affect the ability immune cells to target cancer cells) and stages of the cell cycle. Thus, the graph captured key regulated processes in the response. Second, and more crucially, to highlight key genes the context of the graph, we displayed them onto the half-square G2 [F] based on their strongest link (details in Appendix F). As proof of concept, we focused on two classes of genes: key known regulators we expect to be positioned upstream, such as the interferon receptors (IFNGRI and which sense the IFN signal, JAK/STAT, needed transduce the signal, and CDK4 and 6, which regulate the cell cycle but also repress antigen presentation genes; and others we expect be downstream, in particular those from an immune cancer resistance pathway we previously discovered in patient tumors (ICR [51]; downstream) in Figure 6. Excitingly, while no such information was used to constrain the model, captured these ordered relations, including a causal path from the interferon receptors to interferon signaling modules, from JAK to interferon signaling to STAT, IRFI and HIAA genes, from CDK4 the cell cycle, and from CDK6 and the IFNGRI HLA genes (antigen presentation, MHCI genes) ([14, 51] and references therin). Notably,

---

there are also some connections that may not be borne out biologically. such as the separation of CDK4 / CDK6 different pathways. Overall, DCD-FG a promising starting point for deciphering gene regulation at the scale of the whole transcriptome with Perturb-seq data, and predicting the outcome of interventions that were not tested experimentally.
Discussion
We have proposed DCD-FG, novel approach for large-scale causal discovery that restricts the search space to factor directed acyclic graphs, and efficiently exploits this structure during infer- ence. Our theoretical results suggest that this class of graphs offers statistical benefits under either the faithfulness assumption or under some stochastic edge perturbation model in random graphs. Our numerical experiments show that in important real-world examples, method outperforms NOTEARS, NOBEARS as well as low-rank variants of those methods.
Since the publication of NOTEARS [19], two manuscripts highlighted that the method's evaluation may be confounded by the design of the simulations [52, 53]. However, those studies exclusively focus on causal discovery from observational data. The recent results obtained by DCDI on interventional data with a more suitable simulation design [6], as well as the results of this manuscript on real data show the potential promise of the overall framework.
Recently, several other papers identified the acyclicity constraint as bottleneck for causal discovery learning, and propose to either use the constraint a soft penalty [54], or discard it from the objective function [55] by using a different parameterization of the causal graph learning approach, and interventional data. Although the approach from [54] is currently restricted to learning linear causal models from observational data, [55] could complement our approach by simplifying the optimization procedure but still explicitly model low-rank interactions.
We presented an application of DCD-FG to a large-scale high-throughput gene expression dataset with genetic perturbations ("Perturb-Seq"). The method had better predictive performance for held-out perturbations than state-of-the ahtt-at, and identified both well established relations and new intriguing ones, offering great utility to biologists. Notably, some of the causal relationships may not be accurate, and it is likely that several assumptions of the underlying model may be violated. The biological evaluation and validation of the method will therefore be important to more deeply assess the performance of DCD-FG.
Future work will explore the specification of the noise model, for which count distributions (potentially as part of latent variable model) may be more appropriate [56], well as the absence of confounding variables such as cell cycle [14]. Additionally, we plan to investigate extending the framework of DCD-FG to the inference of causal models with feedback loops [57] in order generate an even more exact and biologically interpretable causal graph. Finally, having a Bayesian alternative of DCD-FG (e.g., based on [58]) would allow scientists to apply those methods for automated experimental design and scientific discovery.
Acknowledgments and Disclosure of Funding
We thank Sébastien Lachapelle, Philippe Brouillard, Alexandre Drouin, Chandler Squires and Gonçalo Rui Alves Faria for insightful conversations about causal structure learning problems. We thank Natasa Tagasovska, Stephen Ra, and Kyunghyun Cho for general conversations about causal inference and biology. We also acknowledge Katie Sciger-Schuller, Chris Frangich, Taka Kudo, Josh Weinstock and Basak Eraslan for discussions about Perturb-seg and the CerturbCCITE-ee dataset. We warmly thank Geoffrey Négiar, Natasa Tagasovska, and Tara Chari for their constructive eriticisms on draft of this paper.
Disclosures: Romain Lopez and Jan-Christian Huetter are employees of Genentech. Jan-Christian Huetter has equity in Roche. Jonathan Pritchard acknowledges support from grant R0THG008 140 from the National Human Genome Research Institute. Aviv Regev is co-founder and equity holder of Celsius Therapeutics and an equity holder in Immunitas. She was an SAB member of ThermoFisher Scientific, Syros Pharmaceuticals, Neogene Therapeutics, and Asimov until July 31st, 2020; she has been an employee of Genentech since August 1st, 2020, and has equity in Roche.
10

---

References
[1] Karen Sachs, Omar Perez, Dana Pe' er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein- signaling networks derived from multiparameter single-cell data. Science, 308(5721):523 529, 2005.
[2] Bin Zhang, Chris Gaiteri, Liviu-Gabriel Bodca, Zhi Wang, Joshua McFlwee, Alexei A Podtelezhnikov, Chunsheng Zhang, Tao Xie. Linh Tran, Radu Dobrin, et al. Integrated systems approach identifies genetic nodes and networks in late-onset alzheimer's disease. Cell, 153(3)7707-700, 2013.
[3] Eran Segal, Dana Pe'er, er, Aviv Regev, Daphne Koller, and Nir Friedman. Learning module networks. Journal of Machine Learning Research, 6(4), 2005,
[4] Daphne Koller and Nir Friedman. Probabilistic Graphical Models. MIT Press, 2009.
[5] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of Causal Inference Foundations and Learning Algorithms. MIT Press, 2017.
[6] Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien. and Alexandre Drouin. Differentiable causal discovery from interventional data. In Advances Neural Information Processing Systems, 2020.
Yuhao Wang, Liam Solus, Karren Dai Yang. and Caroline Uhler. Permutation-based causal inference algorithms with interventions. In Advances in Neural Information Processing Systems, 2017.
[8] Peter Spirtes, Clark N. Glymour, and Richard Scheines. Causation, Prediction, and Search. The MIT Press, 2000.
[9] Theo Diamandis, Yonathan Murin, and Andrea Goldsmith. Ranking causal influence of financial markets via directed information graphs. In Annual Conference Information Sciences and Systems (CISS), 2018.
[10] Katerina Rigana, Erast C Wit, and Samantha Cook. Using network-based causal inference to detect the sources of contagion in the currency market. arXiv, 2021.
[11] Rakesh Malladi, Giridhar Kalamangalam, Nitin Tandon, and Behnaam Aazhang. Identifying seizure onset zone from the causal connectivity inferred using directed information. IEEF Journal of Selected Topics in Signal Processing, 1077):126--128, 2016.
[12] Joseph Ramsey, Madelyn Glymour, Ruben Rombeo, and Clark Glymour. A million variables and more: the fast greedy equivalence scarch algorithm for lcarning high-dimensional graphical causal models, with an application to functional magnetic resonance images. International Journal of Data Science and Analytics, 3(2) 2017.
[13] Atray Dixit, Oren Parnas, Biyu Li, Jenny Chen, Charles P Fulco, Livnat Jerby-Arnon, Nemanja D Marjanovic, Danielle Dionne, Tyler Burks, Raktima Raychowdhury, al. Perturb-seq: dissecting molecular circuits with scalable single-cell RNA profiling pooled genetic screens. Cell, 167(7):1853 1866, 2016.
[14] Chris J Frangich, Johannes C Melms, Pratiksha I Thakore, Kathryn R Geiger-Schuller, Patricia Ho, Adrienne M Luoma, Brian Cleary. Livnat Jerby-Aron, Shruti Malu, Michael S Cuoco, et al. Multimodal pooled Perturb-CITE-see screens in patient models define mechanisms of cancer immune evasion. Nature Genetics, 53(3):332-341) 2021.
[15] Thomas Norman, Max A Horlbeek, Joseph M Replogle, Alex Y Ge, Albert Xu, Marco Jost, Luke Gilbert, and Jonathan Weissman. Exploring genetic interaction manifolds constructed from rich single- cell phenotypes. Science, 365(6455))786-793, 2019.
[16] Joseph M. Replogle, Reuben A. Saunders, Angela Pogson, Jeffrey A. Hussmann, Alexander Lenail, Alina Guna, Lauren Mascibroda, Fric Wagner, Karen Adelman, Jessica I.. Bonnar, Marco Jost, Thomas M. Norman, and Jonathan S. Weissman. Mapping information-rich genotype-phenatype landscapes with genome-scale Perturb-seq. Cell, 2022.
[17] Trevor Hastie, Robert Tibshirani, and Martin Wainwright. Statistical Learning with Sparsity: The Lasso and Generalizations. Chapman & Hall/CRC, 2015.
[18] Yehuda Koren, Robert Bell, and Chris Volinsky. Matrix factorization techniques for recommender systems. Computer, 42(8)333337 2009.
[19] Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric Xing- DAGs with NO TEARS: continuous optimization for structure learning. Advances Neural Information Processing Systems, 2018.
[20] Zhuangyan Fang, Shengyu Zhu, Jiji Zhang, Yue Liu, Zhitang Chen, and Yangbo He. Low rank directed acyelic graphs and causal structure learning. arXiv, 2020.
[21] Shuyu Dong and Michèle Sebag. From graphs to DAGs: low-complexity model and a scalable algorithm. arXiv, 2022.
[22] Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien caadient-based neural DAG learning. In International Conference Learning Representations, 2020.
11

---

[23] Ignavier Ng. Shengyu Zhu, Zhuangyan Fang. Haoyang Li. Zhitang Chen, and Jun Wang. Masked gradient- causal structure learning. In SIAM International Conference Data Mining (SDM). 2022.
[24] Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. Learning sparse nonpara-
metric DAGs. In International Conference Artificial Intelligence and Statistics, 2020.
[25] Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN; DAG structure learning with graph neural networks.
In International Conference Machine Learning, 2019.
[26] Hao- Choh Lee. Matteo Damieletto, Riccardo Miotto, Sarah T. Cherng, and Joel I. Dudley. Scaling structural learning with NO-BEARS to infer causal transcriptome networks. In Pacifie Symposium on Biocomputing. 2019.
[27] Brian Cleary, Le Cong, Anthea Cheung, Eric S Lander, and Aviv Regev. Efficient generation of transcrip- tomic profiles by random composite measurements. Cell, 171(6):1424-146 2017.
[28] Graham Heimberg, Rajat Bhatnagar, Hana El-Samad, and Matt Thomson. Low dimensionality in gene expression data enables the accurate extraction of transcriptional programs from shallow sequencing. Cell Systems, 2(4):239-250, 2016.
[29] Pauli Miettinen and Stefan Neumann. Recent developments in Boolean matrix factorization. In Interna- tional Joint Conferences Artificial Intelligence, 2021.
[30] Frederick Eberhardt and Richard Scheines. Interventions and causal inference. Philosophy of Science, 2007.
[31] Liam Solus, Yuhao Wang, and Caroline Uhler. Consistency guarantees for greedy pernutation-based causal inference algorithms. arXiv, 2017.
[32] Clark Glymour, Kun Zhang, and Peter Spirtes. Review of causal discovery methods based on graphical models. Frontiers in Genetics, 10:524, 2019.
[33] Joris M Mooij, Sara Magliacane, and Tom Claassen. Joint causal inference from multiple contexts, Journal of Machine Learning Research, 2020.
[34] Hans-Andrea Loeliger. An introduction to factor graphs. IEEE Signal Processing Magazine, 21(1):28-41, 2004.
[35] Karren D. Yang, Abigail Katcoff, and Caroline Uhler. Characterizing and learning equivalence classes of causal DAGS under interventions. In International Conference on Machine Learning, 2018.
[36] Paul Erdōs. Alfréd Rényi, et al. On the evolution of random graphs. Publication of the Mathematical Institute of the Hungarian Academy of Sciences, 5(1):17-60, 1960.
[37] Samuel Wang and Mathias Drton. High-dimensional causal discovery under non-Gaussianity. Biometrika, 107(1):4-59)
2019.
[38] Caroline Uhler. Garvesh Raskutti, Peter Bühlmann, and Bin Yu. Geometry of the faithfulness assumption in causal inference. The Annals of Statisties, pages 436-463, 2013.
[39] Marco F Eigenmann, Prectam Nandy, and Marloes H Maathuis. Strueture learning of linear gaussian structural equation models with weak edges. Uncertainty in Artificial Intelligence (UAI), 2017.
[40] Reka A Kovacs, Oktay Gunluk, and Raphael A Hauser. Binary matrix factorisation via column generation. In AAAI Conference Artificial Intelligence, volume 35, pages 3823-3831, 2021.
[41] Derek DeSantis. Erik Skau, Due P. Truong. and Boian Alexandrov. Factorization of binary matrices: Rank relations, uniqueness and model selection of boolean decomposition. ACM Transactions on Knowledge Discovery Data, 2022.
[42] Diviyan Kalainathan, Olivier Goudet, Isabelle Guyon, David Lopez-Pau Paz, and Michèle Sebag. Structural agnostic modeling: Adversarial learning of causal graphs. arXiv, 2018
[43] Chris J. Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous relaxation of discrete random variables. In International Conference on Machine Learning, 2017.
[44] Chandler Squires, Yuhao Wang, and Caroline Uhler. Permutation based causal structure learning with unknown intervention targets. Uncertainty in Artificial Intelligence, 2020.
[45] Diviyan Kalainathan and Olivier Goudet. Causal discovery toolbox: Uncover causal relationships in python. Journal of Machine Learning Research, 21:1-5, 2020.
[46] Amanda Gentzel, Dan Garant, and Jensen. The case for evaluating causal models using interventional measures and empirical data. In Advances in Neural Information Processing Systems, pages 11722-1 1732, 2019.
[47] Nir Friedman, Michal Linial, Iftach Nachman, and Dana Pe'er. Using Bayesian networks to analyze expression data, Journal of Computational Biology, 7(3-4):601-62) 620, 2000.
12

---

[48] Dana Pe'er, Aviv Regev, Gal Elidan, and Nir Friedman. Inferring subnetworks from perturbed expression profiles. Bioinformaticx, 17, 2001.
[49] Dominic Grun, Lennart Kester, and Alexander van Oudenaarden. Validation of noise models for single-cell transcriptomics. Nature Methods, 2014.
[50] Edward Y Chen, Christopher M Tan, Yan Kou, Qiaonan Duan, Zichen Wang, Gabriela Vaz Meirelles. Neil R Clark, and Avi Ma' ayan. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformaties, (4(1):1-13 2013.
[51] Livnat Jerny-Amon, Parin Shah, Michael S Cuoco, Christopher Rodman, Mci-Ju Su, Johannes Melms, Rachel Lecson, Abhay Kanodia, Shaolin Mci, Jia-Ren Lin, ot al. A cancer cell program promotes t cell exclusion and resistance to checkpoint blockade. Cell, 175(4):9849997, 2018.
[52] Marcus Kaiser and Maksim Sipos. Unsuitability of NOTEARS for causal graph discovery. arXiv, 2021.
[53] Alexander G. Reisach, Christof Sciler, and Sebastian Weichwald. Beware of the simulated DAG! causal discovery benchmarks may be casy to game. In Advances in Neural Information Processing Systems, 2021.
[54] Ignavier Ng. AmirEmad Ghassami, and Kun Zhang. On the role sparsity and dag constraints for learning linear dags. In Advances in Neural Information Processing Systems, 2020.
[55] Phillip Lippe, Taco Cohen, and Efstratios Gavves. Efficient neural causal discovery without acyclicity constraints. International Conference on Learning Representations, 2022.
[56] Romain Lopez. Jeffrey Regier, Michael B Cole, Michael I Jordan, and Nir Yosef. Deep generative modeling for single-cell transcriptomics. Nature Methods, 15(12):0033-1088 2018.
[57] Jacob W Freimer, Oren Shaked, Sahin Naqvi, Nasa Sinnott- Arwa Kathiria, Christian M Garrido, Amy F Chen. Jessica T Cortez, William J Greenleaf, Jonathan K Pritchard, et al. Systematic discovery and perturbation of regulatory genes in human T cells reveals the architecture of immune networks. Nature Genetics, 54(8):1133-1144, 2022.
[58] Lars Lorch, Jonas Rothfuss, Bernhard Schölkopf, and Andrcas Krause, DiBS: Differentiable bayesian structure learning. Advances in Neural Information Processing Systems, 34, 2021.
[59] Tom S. Verma and Judea Pearl. Equivalence and synthesis of causal models. In Uncertainty in Artificial Intelligence, 1990.
[60] Pauli Micttinen, Taneli Mielikäinen, Aristides Gionis, Gautam Das, and Heikki Mannila. The discrete basis problem. IEFE transactions on knowledge and data engineering,
2008
[61] Eric Jang, Shixiang Gu, and Ben Poolc. Categorical reparamctrization with Gumble-Softmat In Interna- tional Conference on Learning Representations, 2017.
[62] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Andreas Müller, Joel Nothman, Gilles Louppe, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau. Matthieu Brucher, Matthieu Perrot, and Edouard Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-2830, 201
[63] Alexander Wolf, Philipp Angerer, and Fabian J Theis. SCANPY: large-scale single-cell gene expression data analysis. Genome Biology, 19(1):1-5, 2018.
[64] Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. International Conference on Artificial Intelligence and Statistics, 2010.
[65] Tijmen Tieleman and Geoffrey Hinton. Lecture 6.5- 6.5-RMSProp: Divide the gradient by running average of its recent magnitude. Coursera: Neural networks for machine learning, 2012.
13

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] (Discussion section)
(c) Did you discuss any potential negative societal impacts of your work? [Yes] (Discussion section)
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] (Appendix)
(b) Did you include complete proofs of all theoretical results? [Yes] (Appendix)
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or URL)? [Yes] (Appendix)
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] (Appendix)
(c) Did you report error bars (e.g.. with respect to the random seed after running experi- ments multiple times)? [Yes] (Box plots in main figure)
(d) Did you include the total amount of compute and the type of resources used (c.g., type of GPUs, internal cluster, cloud provider)? [Yes] (Appendix)
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] (Appendix)
(b) Did you mention the license of the assets? [Yes] (Appendix)
(c) Did you include any assets either in the supplemental material URL? [Yes]
(Appendix)
(d) Did you discuss whether and how consent was obtained from people whose data you're ror're using/curating? [N/A]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
14


Reviews:
[{"review_id": "zj_swDhfZd2", "reviewer": "Reviewer_ADLi", "summary": "The paper proposed a new DAG constraints for low rank adjacency matrices. Compared to original DAG constraints, the new constraint can scale to larger graphs. On binary graph, the proposed DAG constraints is a necessary and sufficient condition for an adjacency matrix to form a DAG. However, for a weighted adjacency matrix, it seems that without further work, the proposed DAG constraint is a sufficient but not necessary condition. In that case the proposed method actually is optimising over a space that is smaller than the true DAG space.", "questions": "See my comments above.", "limitations": "1. The main result is trivial and it is only an sufficient condition of DAG.\n2. There may be no real identifiability result.", "rating": 5, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The idea of the paper is very similar to [1], and [1] is rejected by ICLR2020 and now it is only a preprint. Unfortunately, the paper share the same fatal problem as [1]. Let us consider an adjacency matrix W = UV, if we enforce that \nWhen apply such low rank DAG constraints, it is easy to derive a necessary and sufficient condition on binary adjacency matrices. If W is a weighted matrix, then we have to first make the entries of W to be positive to apply the NOTEARS type DAG constraints. In this procedure, in order to use the low rank property, one has to apply absolute value, or square over entries of U and V. In that case, the DAG constraint becomes sufficient only.\n\nThe main result (Proposition 2) of the paper is actually trivial. The proof in the supplementary is over complicated. In fact the proof of this trivial result is very straightforward. For a matrix W \\in R_{>=0}^{n\\times n} = UV, where all entries of U and V are also no-negative, W is nilpotent if and only W is a DAG. Furthermore, VU is nilpotent if and only if UV is nilpotent, and thus VU is DAG is equivalent to UV is DAG.\n\nThere is also a minor problem. The rank of DAG graphs can be very large. For example, if we consider a chain structured DAG with n nodes, the rank of its adjacency matrix is n - 1. Thus the proposed method only suits sparse graph structures. \n\nFinally, in the identifiability section, the authors provide a theorem about the rank of the graph? How the the theorem related to the identifiability? Will this enforce the faithfulness of the discovered graph?\n[1] Fang, Zhuangyan, et al. \"Low rank directed acyclic graphs and causal structure learning.\" arXiv preprint arXiv:2006.05691 (2020).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}, {"review_id": "q_sN4pbIeW3", "reviewer": "Reviewer_AHb3", "summary": "This paper proposes an approach to causal discovery in high dimensional settings where interventional data is available. The authors use a low-rank assumption on a factor graph, and incorporate the model into a differentiable causal discovery algorithm by utilizing a likelihood model which assumes perfect interventions. Empirical results show strong performance with respect to prior art.", "questions": "* Is there a parameterization that would allow for the rank to be chosen as well? Perhaps I am missing something along these lines.\n* The work assumes perfect interventions in the likelihood. Would this line of work extend to imperfect interventions? Purely observational settings? It would be nice if the authors were able to clearly delineate within the text where the interventions are strictly required.\n* Is it possible to have a set of experiments, as detailed above, that show robustness to misspecification?", "limitations": "The authors do a nice job of describing their work within the parameters of necessary assumptions and constraints. I don't foresee clear societal impact issues ehre.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n* The introduction of a method that allows for a low rank factor assumption is very nice, and provides a valuable tool in causal discovery where sparsity can be a difficult assumption to enforce.\n* The approach is a nice and relatively simple solution to a difficult problem.\n* The parameterization of the model with the Gumbel-softmax prior is quite nice.\n\nWeaknesses:\n\n* Unclear how the number of factors should be chosen. \n* Experimental evidence which shows the robustness to misspecification (both in terms of rank and the Gaussian assumption) would be helpful.\n* (Minor) It would be good if the authors more clearly delineated that this paper deals solely with interventional data. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "4FaJo88pvT", "reviewer": "Reviewer_WAgt", "summary": "This paper aims to solve the large-scale causal discovery problem with hundreds of nodes. To this end, the authors first introduce factor directed acyclic graphs (f-DAGs) which have many promising properties by incorporating the structural assumption to restrict the search space and reduce the optimization difficulty. Then, they propose Differentiable Causal Discovery of Factor Graphs (DCD-FG) to complete large-scale causal discovery. Theoretical analyses and experimental results demonstrate the superiority of f-DAGs and DCD-FG.", "questions": "1. How to justify that the low rank assuption is commonly met in real scenarios?\n\n2. Although intuitively making sense, is there a theoretical result for the claim that \"more statistically robust in the high-dimensional regime where the underlying skeleton is hard to assess\"?\n\n3. It seems most experimental results suround the SHD, and F1 score. Since your claim is fast searching, can you also show some computational efficiency results? ", "limitations": "There are no  limitations and potential negative societal impact of their work.", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n1. This paper effectively reduce the search space in the large-scale causal discovery problem by introducing factor directed acyclic graphs. Based on which they propose Differentiable Causal Discovery of Factor Graphs to realize scalable causal discovery.\n\n2. The experiments are well designed. This paper provides runtime experiments, simulation studies, and a case study on single-cell RNA sequencing data with hundreds of genetic interventions to demonstrate the efficacy of the proposed method.\n\nWeakness:\n1. It seems quite confusing in line 41-44. In the first sentence, it says “sparsity [17] as well as low-rank assumptions [18] are often exploited in algorithms”. However, in the next sentence, it points out “the use of low-rank constraints remains largely under-explored for this purpose”.\n\n2. The theoretic property of this work is not well presented. The proposition 1 is only a very special case of lpw rank models. It would be better to justify the method in a wider scenarios.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept