Learning DAGs from Data with Few Root Causes
Authors: Panagiotis Misiakos
Chris Wendler
Markus Püschel
Venue: NeurIPS 2023 poster (2023)
Keywords: directed acyclic graph
few root causes
structural equation models
linear SEMs
additive noise

Abstract: We present a novel perspective and algorithm for learning directed acyclic graphs (DAGs) from data generated by a linear structural equation model (SEM). First, we show that a linear SEM can be viewed as a linear transform that, in prior work, computes the data from a dense input vector of random valued root causes (as we will call them) associated with the nodes. Instead, we consider the case of (approximately) few root causes and also introduce noise in the measurement of the data. Intuitively, this means that the DAG data is produced by few data generating events whose effect percolates through the DAG. We prove identifiability in this new setting and show that the true DAG is the global minimizer of the $L^0$-norm of the vector of root causes. For data satisfying the few root causes assumption, we show superior performance compared to prior DAG learning methods.

Learning DAGs from Data with Few Root Causes
Panagiotis Misiakos, Chris Wendler and Markus Püschel Department of Computer Science, ETH Zurich {pmisiakos, wendlerc, markusp)@ethz. ch
Abstract
We present a novel perspective and algorithm for learning directed acyclic graphs (DAGs) from data generated by linear structural equation model (SEM). First, we show that a linear SEM can be viewed as a linear transform that, in prior work, computes the data from dense input vector of random valued root causes (as we will call them) associated with the nodes. Instead, we consider the case of (approximately) few root causes and also introduce noise in the measurement of the data. Intuitively, this means that the DAG data is produced by few data generating events whose effect percolates through the DAG. We prove identifiability in this new setting and show that the true DAG is the global minimizer of the L -norm of the vector of root causes. For data satisfying the few root causes assumption, we show superior performance compared to prior DAG learning methods.
1 Introduction
We consider the problem of learning the edges of an unknown directed acyclic graph (DAG) given data indexed by its nodes. DAGs represent causal dependencies (edges) between events (nodes) in the sense that an event only depends on its predecessors. Thus, DAG learning has applications in causal discovery, which, however, is more demanding problem that do consider here, since it requires further concepts of causality analysis, in particular interventions [Peters et al.) 2017). However, DAG learning is still NP-hard in general [Chickering et al. 2004). Hence, in practice, one has to make assumptions on the data generating process, to infer information about the underlying DAG in polynomial time. Recent work on DAG learning has focused on identifiable classes of causal data generating processes. A frequent assumption is that the data follow a structural equation model (SEM) [Shimizu al. 2006 Zheng et al., 20188 Gao et al., 202T], meaning that the value of every node is computed as function of the values of its direct parents plus noise. Zheng et al. (2018) introduced NOTEARS, prominent example of DAG learning, which considers the class of linear SEMs and translates the acyclicity constraint into a continuous form for easier optimization. It inspired subsequent works to use continuous optimization schemes for both linear [Ng et al.) 2020) and nonlinear SEMs [Lachapelle et al.) 2019 Zheng et al. 2020 A more expansive discussion of related work provided in Section 4
In this paper we also focus on linear SEMs but change the data generation process. We first translate the common representation of linear SEM as recurrence into an equivalent closed form. this form, prior data generation can be viewed as linearly transforming an i.i.d. random, dense vector of root causes (as we will call them) associated with the DAG nodes as input, into the actual data on the DAG nodes as output. Then we impose sparsity in the imput (few root causes) and introduce measurement noise in the output. Intuitively, this assumption captures the reasonable situation that DAG data may be mainly determined by few data generation events of predecessor nodes that percolate through the DAG as defined by the linear SEM. Note that our use the term root causes related to but different from the one in the root cause analysis by Ikram al. 2022).
Contributions. We provide a novel DAG learning method designed for DAG data generated by linear SEMs with the novel assumption of few root causes. Our specific contributions include the following:
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

We provide an equivalent, closed form of the linear SEM equation showing that it can be viewed linear transform obtained by the felexive-transitive closure the DAG's adjacency matrix. In this form, prior data generation assumed a dense, random vector as input, which we call root causes.
In contrast, we assume a sparse input, or few root causes (with noise) that percolate through the DAG produce the output, whose measurement is also subject to noise. Interestingly, Seifert et al. 2023] identify the assumption of few root causes as a form of Fourier-sparsity.
We prove identifiability of our proposed setting under weak assumptions, including zero noise in the measurements. We also show that, given enough data, the original DAG the unique minimizer of the associated optimization problem in the case of absent noise.
We propose novel and practical DAG learning algorithm, called SparseRC, for setting, based on the minimization of the I--noorm of the approximated root causes.
We benchmark SparseRC against prior DAG learning algorithms, showing significant improvements for synthetic data with few root causes and scalability thousands of nodes. SparseRC also performs among the best on real data from gene regulatory network, demonstrating that the assumption few root causes can be relevant in practice.
2 Linear SEMs and Root Causes
We first provide background on prior data generation for directed acyclic graphs (DAGs) via linear SEMs. Then we present a different viewpoint on linear SEMs based on the concept of root causes that we introduce. Based on it, we argue for data generation with few root causes, present it mathematically and motivate it, including with real-world example.
DAG. We consider DAGs 9 = (V, E) with V vertices, E the set of directed edges, no cycles including no self-loops. We assume the vertices to be sorted topologically and set accordingly V {1, 1,,,,..,dd d}. Further, dij E IR is the weight of edge (i,j) E E, and
{ aij, if (i.j) E E,
A = (aij)isev
0, else.
is the weighted adjacency matrix. A is upper triangular with zeros on the diagonal and thus A
0.
Linear SEM. Linear SEMs [Peters et al. 2017 formulate a data generating process for DAGs g. First, the values at the sources of a DAG G are initialized with random noise. Then, the remaining nodes are processed in topological order: the value vy node j is assigned the linear combination of its parents' values (called the causes of it)) plus independent noise. Mathematically. data vector
x == (Ii,..., v.) E IRd on G follows a linear SEM if x,- = xA., + n2, where the subscript j denotes column j. and n, are i.i.d. random noise variables [Shimizu et al. 20066 (Zheng et al.) 2018 Ng al. 2020). Assuming n such data vectors collected as the rows of matrix X E R=xd and the noise variables in the matrix N, the linear SEM is typically written as
X = XA+N.
2.1 Transitive closure and root causes
Equation 2 can be viewed as recurrence for computing the data values X from N. Here, we interpret linear SEMs differently by first solving this recurrence into closed form. We define
A- A A +.... A d-1, which is the Floyd- Warshall (FW) transitive closure of A over the ring (R. + Lehmann 1977]. and + A the associated -raaniii-r closure of A. Since A have (I A) + A) = I and thus can isolate X in
Lemma 2.1. The linear SEM computes data X
X N
In words, the data values in X are computed as linear combinations of the noise values N of all predecessor nodes with weights given by the reflexive- closure I A.

---

Since X is uniquely determined by N, we call the latter the roor causes of X. The root causes must be confused with the root nodes (sources) of the DAG, which are the nodes without parents. In particular, N can be non-zero in nodes that are not sources.
Few root causes. In measured output at each node. With this viewpoint, we argue that it makes sense to consider data generation process that differs two ways from (3) and thus 22) First, we assume that only few nodes produce a relevant input that we call C, up to low magnitude noise Ne. Second, we assume that the measurement of X is subject to noise N,. Formally, this yields the closed form
view N as the input to the linear SEM at each node and X as the
X = (C + No) (I A) + N
(4)
Multiplying both sides by (I A) -1 = (I- - A) yields the (standard) form as recurrence
x= XA (C+ N.) + N (I A),
(5)
i.e., N in (2) is replaced by C + N. + N₂ (I A), which in general is not i.i.d. Note that (4) generalizes (3) obtained by assuming zero root causes C C 0 and zero measurement noise N 0.
We consider C not as an additional noise variable but as the actual information, i.c., the relevant input data at each node, which then percolates through the DAG as determined by the SEM to produce the final output data X, whose measurement, as usual in practice, is subject to noise. Few root causes mean that only few nodes input relevant information in one dataset. This assumption can be stated as:
Definition 2.2 (Few Root Causes assumption). We assume that the input data (C N.) (I- () to linear SEM are approximately sparse. Formally this holds if:
|Cllo//nd < c (few root causes),
|| + N (I - A)ll /nd
<6 (negligible noise),
(6)
|C
where I|-Ilo counts the nonzero entries of C, ||-1|1 is elementwise, and e, ō are small constants.
For example, in Appendix A we show that (6) holds (in expectation) for c = = 0.1 N.. N₂ are normally distributed with zero mean and standard deviation 0.01 and the ground truth DAG g has uniformly random weights in [-1,1] 1] and average degree 4. The latter two conditions bound the amplification of N in (6) when multiplied by A.
The term root cause in the root cause analysis by [kram et al. [2022] refers to the ancestor that was the only cause of some effect(s). If there is only one root cause in our sense, characterized by non-zero value, then determining its location is the problem considered by Ikram et al. 2002
2.2 Example: Pollution model
Linear SEMs and associated algorithms have been extensively studied in the literature (Loh and Bühlmann 2014 Peters and Bühlmann| [2014 Ghoshal and Honorio 2017 Aragam and Zhou| 2015). However, we are not aware of any real-world example given in a publication. The motivation for using linear SEMs includes the following two reasons: d-variate Gaussian distributions can always be expressed as linear SEMs [Aragam and Zhou 2015) and linearity is often workable assumption when approximating non-lineas systems.
Our aim is to explain why we propose the assumption of few root causes. The high-level intuition is that it can be reasonable to assume that, with the view point of (3), the relevant DAG data is triggered by sparse events on the input size and not by random noise. To illustrate this, we present a linear SEM that describes the propagation of pollution in a river network.
DAG. We assume a DAG describing a river network. The acyclicity is guaranteed since flows only occur downstream. The nodes i E V represent geographical points of interest, e.g., cities, and edges are rivers connecting them. We assume that the cities can pollute the rivers. An edge weight agj [0, 1]. (i. E, captures what fraction of pollutant inserted at reaches the neighbour j. An example DAG with six nodes is depicted in Fig. Ta
Transitive closure. Fig. 1b shows the transitive closure of the DAG in (a). The (i, entry of the transitive closure is denoted with aij and represents the total fraction of pollutant at that reaches via all connecting paths.
3

---

1.5
4.66
0.6
0.8
0.21
10.50
0.5
6.63
0.67
Original DAG
(b) Transitive closure
(c) Measurements vector
(d) Root
Figure 1: (a) A DAG for a river network. The weights capture fractions of pollution transported between adjacent nodes. (b) The transitive closure. The weights are fractions of pollution transported between all pairs of connected nodes. (c) A possible vector measuring pollution, and (d) the root causes of the pollution, sparse in this case.
Data and root causes. Fig. Ic shows possible data vector the DAG, for example, the pollution measurement at each node done once day (and without noise in this case). The measurement the accumulated pollution from all upstream nodes. Within the model, the associated root causes in Fig. then show the origin of the pollution, two in this case. Sparsity in means that each day only small number of cities pollute. Negligible pollution from other sources is captured by noise N. and N₂ models the noise in the pollution measurements (both are assumed to be 0 in Fig. T.
We View the pollution model as abstraction that can be transferred to other real-world scenarios. For example, in gene networks that measure gene expression. few root causes would mean that few genes activated in considered dataset. In a citation network where measures the impact of keywords/ideas, few root causes would correspond to the few origins of them.
Learning the DAG
In this section present approach for recovering the DAG adjacency matrix A from given data X under the assumption of few root causes, i.e., 4 (or (5)) and (6). We first show that under no measurement noise, our proposed setting identifiable. Then generating model in the absence of noise and prove that the true DAG adjacency A is the global minimizer of the - of the root causes. Finally, to learn the DAG in practice, we perform continuous relaxation to obtain an optimization problem that is solvable with differentiation.
theoretically analyze our data
3.1 Identifiability
The of the data X generated by linear SEM commonly interpreted observations of
random row vector X = (X1, ..., Xa) [Peters and Bühlmann 2014 Zheng et al. 2018). i.e., is random noise variables, or, equivalently, as X XA + N. Given distribution Pv of the noise vector, the DAG with graph adjacency matrix A is then called identifiable if it uniquely determined by the distribution Px of X. (Shimizu et al. (2006) state that A is identifiable if all N, are non-Gaussian noise variables. This applies directly to our setting and yields the following result.
written analogous to 22)
X, XA NN; where N; L.I.d.
Theorem 3.1. Assume data generation X (C+ (I+ A), where C and N are independent. Let p E [0, 1) and assume the C₁ are independent random variables taking uniform values from [0, 1] with probability p. and are = O with probability 1 1-p. p. The noise vector No is defined as before. Then the DAG given by A is identifiable.
Proof. Using 5). the data generation equation can be viewed as a linear SEM (in standard recursive form) with noise variable C+NN which is non-Gaussian due Lévy-Cramér decomposition theorem
Lévy 1935 Cramér 1936). because C is non-Gaussian. The statement then follows from LiNGAM Shimizu al.] 2006] given that the noise variables of C' and N. are independent. Note that this independency would be violated if we assumed non-zero measurement noise NN, as in
In Theorem 3.1 C yields sparse observations whenever p is close to zero (namely dp root causes in expectation). However, identifiability follows for all p due non-Gaussianity.

---

3.2 L minimization problem and global minimizer
Suppose that the data X are generated via the noise-free version of 4
X = (I+A)
(7)
We assume C to be generated in Theorem 3.1 i.e., with randomly uniform values from (0, 1] with probability p and O with probability 1 p. where p is small such that sparse. Given the data X we aim to find A by enforcing maximal sparsity on the associated C, i... by minimizing ||
min AcRad
x (I+AA)
s.t. A is acyclic.
8
The following Theorem 3.2 states that given enough data, the true A is the unique global minimizer of (8) This means that that case the optimization problem (8) correctly specifies the true DAG. One can view the result as a form of concrete, non-probabilistic identifiability. Note that Theorem 3.2 not a sample-complexity result, but a statement of exact reconstruction in the absence of noise. The sample complexity of our algorithm is empirically evaluated later.
Theorem 3.2. Consider a DAG with weighted adjacency matrix A with d nodes. Given exponential (in d) number n of samples X the matrix A is, with high probability, the global minimizer of the optimization problem
Proof. See Appendix
3.3 Continuous relaxation
In practice the optimization problem 8 is too expensive to solve due its combinatorial nature, and, of course, the noise-free assumption rarely holds in real-world data. We consider again our general data generation in (4) assuming sparse root causes C and noises N:. N₂ satisfying the criterion (6). To solve it consider continuous relaxation of the optimization objective. Typically (e.g., (Zheng
jet al. 2018 Lee et al., 2019 Bello et al. 2022)) ),, continuous optimization DAG learning approaches
have the following general formulation:
min AERaeI
+ A)
s.t.
h (A)= =0
where / (A, X) is the loss function corresponding to matching the data, R (AA) ) is a regularizer that promotes sparsity in the adjacency matrix, usually equal to || and h (A) is a continuous constraint enforcing acyclicity.
In our case, following common practice in the literature, we substitute the -orm from 8 with its convex relaxation [Ramirez et al. 2013), the L -norm. Doing so allows for some robustmess (low magnitude) noise N. and N m- We capture the acyclicity with the continuous constraint h (A) = Ir (cAa used by /Zhnn et al. (2018). but could also use the form of Yu et al. [2019] As sparsity regularizer for the adjacency matrix and we choose (A) = A |||||| Hence, our final continuous optimization problem is
((I "I,, +AJ|A|I s.t. (AA) = O.
min
(10)
AE
2n
We call this method SparseRC (sparse root causes).
Related work
The approaches solving the DAG learning problem fall into two categories: combinatorial search or continuous relaxation. In general, DAG learning is an NP-hard problem since the combinatorial space of possible DAGs super-exponential in the number of vertices [Chickering et al., 2004). Thus, methods that search the space of possible DAGs apply combinatorial strategies to find an approximation of the ground truth DAG [Ramsey et al. 2017 Chickering 2002 I'samardinos et 2006).
5

---

Continuous optimization. Lately, with the advances of deep learning, researchers have been focusing on continuous optimization methods [Vowels et al. 202T] modelling the data generation process using SEMs. Among the first methods to utilize SEMs were CAM (Bühlmann et al. 2014) and LiNGAM [Shimizu et al., [0006] which specializes to linear SEMs with non-Gaussian noise. NOTEARS by Zheng et al. [2018] formulates the combinatorial constraint of acyclicity as a continuous one, which enables the use of standard optimization algorithms. Despite concerns, such as lack of scale-invariance [Kaiser and Sipos 2022 Reisach et al. 2021]. it has inspired many subsequent DAG learning methods. The current state-of the-art of DAG learning methods for linear
SEMs include DAGMA by |Bello et al. [2022], which introduces log det acyclicity constraint, and
GOLEM by Ng et al. 22020), which studies the role of the weighted adjacency matrix sparsity, the
acyclicity constraint, and proposes to directly minimize the data likelihood. ZZhnng et al. [2020] extended NOTEARS to apply in non-linear SEMs. Other nonlinear methods for DAG learning include DAG-GNN by Yu al. 2019]. in which also a more efficient acyclicity constraint than the one in NOTEARS is proposed, and DAG-GAN by Gao et al.) [002]. DAG-NoCur! by Yu et al.
(202T) proposes learning the DAG the equivalent space of weighted gradients of graph potential functions. Pamfil et al. (2000) implemented DYNOTEARS, a variation of NOTEARS, compatible with time series data. A recent line of works considers permutation-base
methods parameterize
the search space of the DAG [Charpentier et al. 2022 Zantedeschi et al., 2022]. The work in this
paper considers DAG learning under the new assumption of few root causes, In Misiakos et al. 2024
we further build on this work by learning graphs from time series data.
Fourier analysis. Our work is also related to the recently proposed form of Fourier analysis on DAGS from Seifert et al. argue that (what we call) the root causes can be viewed as form of spectrum of the DAG data. This means the assumption of few root causes is equivalent to Fourier-sparsity, a frequent assumption for common forms of Fourier transform [Hassanich] 2018 (Stobbe and Krause 2012 Amrollahi et al.) 2019]. Experiments in [Seifert et al. 2023 Misiakos et al. 2024 reconstruct data from samples under this assumption.
2023 2022], where equation (3) appears as special case. The authors
5 Experiments
We experimentally evaluate our DAG learning method SparseRC with both synthetically generated data with few root causes and real data from gene regulatory networks (Sachs et al., 2005). We
also mention that our method was among the winning solutions [Misiakos et al. 2023 Chevalley
al. [023) of a competition of DAG learning methods to infer causal relationships between genes
[Chevalley et al. 2022). The evaluation was done on non-public data by the organizing institution and is thus not included here.
Benchmarks. We compare against prior DAG learning methods suitable for data generated by linear SEMs with additive noise. In particular, we consider the prior GOLEM (Ng et al. [0220] NOTEARS (Zheng et af.]
2018), DAGMA (Bello et al.] 2022], DirectLiNGAM [Shimizu et al.
2011). PC (Spirtes et al. 2000 and the greedy equivalence search (GES) [Chickering] [0002].
We also compared Sparse RC against LiNGAM [Shimizu et al., 2006). causal additive models (CAM) [Bhhmaann et al. 2014). DAG-NoCurl (Yu et al. 2021]. J. fast greedy equivalence search
(FGES) Ramsey et al., 2017]. the recent simple baseline sortnregress [Reisach et al.) 2021] and max-min hill-climbing (MMHC) (Tsamardinos et al. 2006). but they where not competitive and thus
we only include the results in Appendix B
Metrics. To evaluate the found approximation À of the true adjacency matrix A, we use common performance metrics as in (Ng et 2020 RReisch al. [022]. In particular, the structural Hamming distance (SHD), which is the number of edge insertions, deletions, or reverses needed to convert Á to A, and the structural intervention distance (SID), introduced by Peters and Bühlmann [2015]. which the number of falsely estimated intervention distributions. SHD gives a good estimate of a method's 's performance when it is close to O. In contrast, when it is high it can either be that the approximated DAG has no edges in common with the ground truth, but it can also be that the approximated DAG contains all the edges of the original one together with some false positives. Thus, in such cases, we also consider the true positive rate (TPR) of the discovered edges. We also report the total number of nonzero edges discovered for the real dataset [Sachs et al. 2005) and, for the synthetic data, provide results reporting the TPR, SID, and normalized mean square error (NMSE) metrics in Appendix B For the methods that can successfully Iearn the underlying graph we further

---

Table I: SHD metric (lower is better) for learning DAGs with 100 nodes and 400 edges. Each row is an experiment. The first row is the default, whose settings are in the blue column. In each other row, exactly one default parameter is changed (Change column). The last six columns correspond to prior algorithms. The best results are shown in bold. Entries with SHD > 400 are reported as failure or shown in green if the TPR is > 0.8.
Hyperparameter
Defsult
Change
Varsoet.
SpaneRC
GOLEM
NOTEARS
DAGMA
DirecLINGAM
PC
GES
Default setting
0:95
0.0
62
50
200
262
failure
Graph type
Endos R.
Scale free
0.00
2:2 1 1.5
34 0.00
9.5
205
158
275
16
375 1.
N., N, distrbuttion
Gaussian
Gumbel
0.07
1.4 1.0
87 07:100 44
50:1 17
218 1:5-:77 7.4
267 = 05
250
1: 17
306 : S3
Fdlges / Vettios
10
0.80
46 = 7.5
212 70
263 + 20
896 +
1078 + 106
073 + 17
fsilure
Sundaniizaiion
No
fallure
Fallure
278
Larper weighes
(0.1,6.0)
(0.5.2)
1.00
fallere
96 20:2:25
92:1
200 4.1
840= 121
600-2:0.8
failure
N2 N2 deviation
0.01
a 0.1
0.97
504 19
98 + 14
170 12
239 + 13
538 45
255 + 11
failure
Dense comses C
P-0.5 0.5
0.08
1221 1 55
20 1. 2.5
126 1:32
63 8.0
257 14
24
frilure
Samples
0:97
fallure
326
10
Fixed suppert
No
Yes
0.89
failure
fallure
failure
fallure
fualure
379 37
fallure
proceed on approximating the true root causes C that generated the data X via A. With regard to this aspect, we count the number of correctly detected root causes C TPR, their false positive rate C FPR and the weighted approximation C NMSE. We also report the runtime for all methods in seconds. For each performance metric, we compute the average and standard deviation over five repetitions of the same experiment.
Our implementation. To solve (10) in practice, we implemented Py Torch model with a trainable parameter representing the weighted adjacency matrix A. This allows us utilize GPUs for the acceleration of the execution of our algorithm (GOLEM uses GPUs, NOTEARS does not). Then we use the standard Adam (Kingma and Ba 2014] optimizer to minimize the loss defined in 10
5.1 Evaluation on data with few root causes
Data generating process and defaults. In the second, blue column of Table we report the default settings for our experiment. We generate a random Erdös-Renyi Renyi graph with d = 100 nodes and assign edge directions to make it DAG as in [Zheng et al. 2018). The ratio of edges to vertices is set to 4, so the number of edges is 400, The entries of the weighted adjacency matrix are sampled uniformly at random from b. a) U (a, b), where a = 0.1 and b = 0.9. As (Zheng et al. 2018 Ng et al. 2020 Bello et al. 2022 the resulting adjacency matrix is post-processed with thresholding. In particular, the edges with absolute weight less than the threshold = 0.09 are discarded. Next, the root causes C are instantiated by setting each entry either to some random uniform value from (0. 1) with probability p = 0.1 or to with probability - p = 0.9 (thus, as in Theorem 3.1] the location of the root causes will vary). The data matrix X is computed according to 4 using Gaussian noise Ne; N₂ of standard deviation o = 0.01 to enforce (6). Finally, we do not standardize (scale for variance 1) the data, and X contains n 1000 samples (number of rows).
Experiment 1: Different application scenarios. Table cabpprcompares SparseRC six prior algorithms using the SHD metric. Every row corresponds to a different experiment that alters one particular hyperparameter of the default setting. which the first row with values of the blue column as explained above. For example, the second row only changes the graph type from Erdös-Renyi to scale-free, while keeping all other settings. Note that in the last row, fixed support means that the location of the root causes is fixed for every sample in the data.
Since the graphs have 400 edges, an SHD & 400 can be considered as good and beyond 400 can be considered failure. In row 3 this threshold is analogously set to 1000. The failure cases are indicated as such in Table or the SHD shown green if they still achieve a good TPR 0.8. In Appendix B.T we report the TPR, SID, and the runtimes of all methods.
In the first three rows, examine scenarios that perfectly match the condition 6 of few root causes. For the default settings (row 1), scale-free graphs (row 2), and different noise distribution (row 3) SparseRC performs best and almost perfectly detects all edges. The next experiment, alter parameters that deteriorate the few root causes assumption. Row 4 considers DAGs with average degree 10, i.e., about 1000 edges. High degree can amplify the measurement noise and hinder the root causes assumption (6) However, our method still performs best, but even larger degrees decay the performance of all methods shown in experiment 5 below. The standardization of data (row
'Our code is publicly available at htttg:///ghhbb.oo///mmaae//////s
7

---

Table 2: Runtime [seconds] report of the top-performing methods in Table
Hyperparameter
Defaul
Change
SparseRC (ours)
GOLEM
NOTEARS
1. Default settings
10 + 1.8
529 210
796 185
2. Graph type
Redys-Renyi Gaussian
Scale-free
11 +
460
184
180 +: 7.2
3. No. N, distribution
Gumbel
8.2 + 0.7
349 +: 125
251 = 48
4. Edges / Vertices
14 +
121
5. Samples
1000
100
13 + 0.7
194 9.6
679 + 72
6. Standardization
No
Yes
11
326 145
781 - 76
7. Larger weights in
(0.1. 0.9)
(0.5.
0.6
431 177
2834 228
8. Ne, N₂ deviation
0.01
o= 0.1
8.7
309
+ 63
433 - 53
9. Dense root causes C
p = 0.1
p= 0.5
9.1
331 121
127 = 35
10. Fixed support
Yes
15
15
360 142
660
5) is generally known to negatively affect algorithms with continuous objectives (Reisach al. 2021] as is the case here. Also, standardization changes the relative scale of the root causes and as a consequence affects (6). Row 6 considers edge weights >> which amplifies the measurement noise in (6) and our method fails. Higher standard deviation in the root causes noise (row 7) or explicitly higher density in C (row 8) decreases performance overall. However, In both cases SparseRC still discovers a significant amount of the unknown edges. For a small number of samples (row 9) most methods fail. Ours achieves high TPR but requires more data converge to the solution as we will see in experiment 2 below. Row 10 is out of the scope of our method and practically all methods fail.
Overall, the performance of SparseRC depends heavily the degree to which the assumption of few root causes 66) is fulfilled. The parameter choices in the table cover a representative set of possibilities from almost perfect recovery (row 1) to complete failure (row 10).
SparseRC is also significantly faster than the best competitors GOLEM and NOTEARS with typical speedups in the range of 10-50x as shown in Table 2
It is worth mentioning that even though LiNGAM provides the identifiability Theorem T. of our proposed setting, is not able to recover the true DAG. While both DirectLiNGAM and LiNGAM come with theoretical guarantees for their convergence, these require conditions, such as infinite amount of data, which in practice are not met. We include a more extensive reasoning for their subpar performance in Appendix B. together particularly designed experiment for LiNGAM and DirectLiNGAM.
Varsortability. In Table I we also include the varsortability for each experimental setting. Our measurements for Erdös-Renyi reported in (Reisach et al. [021 Appendix G.1] for linear SEMs, but still high in general. However, the trivial variance-sorting method sortnregress, included in Appendix [ fails overall. Note again that for fixed sparsity support (last row), all methods fail and varsortability is lower. Therefore, in this scenario, our data generating process poses a hard problem for DAG learning.
Renyi graphs (all rows except row are typically 1-2% lower than those
Experiment 2: Varying number of nodes or samples. In this experiment we first benchmark SparseRC with varying number of nodes (and thus number of edges) in the ground truth DAG, while keeping the data samples in X equal to 1000. Second, we vary the number of samples, while keeping the number of nodes fixed to d 100. All other parameters are set to default. The results are shown in Fig. reporting SHD, SID, and the runtime.
When varying the number of nodes (first row in Fig. 2 SparseRC, GOLEM, and NOTEARS achieve very good performance whereas the other methods perform significantly worse. As the number of nodes increases, SparseRC: performs best both in terms of SHD and SID while being significantly faster than GOLEM and NOTEARS. When increasing the number of samples (second row) the overall performances improve. For low number of samples SparseRC, GOLEM and NOTEARS fail. Their performance significantly improves after 500 samples where SparseRC overall achieves the best result. The rest of the methods have worse performance. SparseRC is again significantly faster than GOLEM and NOTEARS in this case.
Learning the root causes. Where method succeeds causes that generate the data. Namely, if we recover a very good estimate of the true adjacency matrix via 1OP. we may compute an approximation Ć of the root causes C, up to noise, by solving 4
Experiment 3:
also recover the root
+ (I A) X (I+
(11)
In the last row of Fig. 2 we evaluate the top-performing methods on the recovery of the root causes (and the associated values) with respect to detecting their locations (C TPR and FPR) and recovering
8

---

700
500
600
7000
6000
402-
400
303
100
Time
3000
2003
100
100
Number nodes
Number nodes
Number of nodes
8000
SparseKC
400
400-
GOLEM
6003
300
303
Time
DAGMA
200
202
DGAM
100
109
GES
100200
500
800
1000
100200
50e
aão
1030
100200
500
1000
Number samples
Number of samples
Number of samples
(a) SHD (L)
(b) SID (L)
(c) Runtime [s]
FPR
SparselCC
GOLEM
685
0.01
6.0
0.50
0.5
20
100
100
20
40
100
Number nodes
Number of nodes
Number of nodes
(d) Root causes TPR (1)
(e) Root causes FPR (L)
(f) Root causes NMSE(L)
Figure 2: Performance report the default settings. (a.b.c) illustrate SHD, SID and runtime (lower
better) while varying the number nodes with 1000 samples (first row) or varying the number of samples with 100 nodes (second row). (d, e) illustrate TPR and FPR of the estimated support of C, and (f) reports the accuracy of estimating as NMSE.
Table 3: Performance of the top-performing methods on larger DAGs.
Nodes samples n
SparseRC
NOTEARS
GOLEM
SparseRC
NOTEARS
GOLEM
d 200, a 500
155
281
1061
600
d 500, n = 1000
27
245
574
204
5789
5428
d 1000, 2 5000
26
282
699
1085
10104
37290
d 2000,
10000
50
489
time-out
7141
46213
time-out
d 3000, 10000
134
time-out
time-out
19660
time-out
time-out
(a) SHD (L)
(b) Runtime (s)
their numerical values using the C NMSE metric. The support is chosen as the entries (i, of C that are in magnitude within 10% the largest value in (i.e., relatively large): Cij > 0.1 Cmax- We consider the default experimental settings.
Experiment Larger DAGs. Table 3 we investigate the scalability of the best methods SparseRC, GOLEM, and NOTEARS from the previous experiments. We consider five different 3000 nodes and = 10000 samples. We use the prior default settings except for
settings up the sparsity of the root causes, which is now set to 5%. The metrics here are SHD and the runtime The results show that SparseRC excels in the very sparse setting with almost perfect reconstruction. far outperforming the others including in runtime. SparseRC even recovers the edge weights with an average absolute error of about 5% for each case in Table as shown in Appendix B.3
Experiment 5: Varying average degree. We evaluate the top-scoring methods on DAGs with higher average degree (dense DAGs). Note that this violates the typical assumption of sparse DAGs [Zheng et al. 2018 8) In Fig. 3 we consider default settings with DAGs of average degree to 20. SID was not computed due to the presence of cycles in the output of the algorithms and TPR is reported instead. We conclude that for dense DAGs the performance decays, expected.

---

1.0-
250
SparsaRC
2003
09
GOLEM
NOTEARS
SHD
150
1002
Time
05.
0-
0+
20
Aug. degree
degree
degree
(a) SHD (
(b) TPR (t)
(c) Runtime [s]
Figure 3: Evaluation of the oo-perffrming methods on denser DAGs.
Table 4: Performance of the top-performing methods on the dataset by Sachs et al. 2005
SHD 4
SID
Total edges
SparseRC NOTEARS
15
45
16
11
44
15
GOLEM
43
19
DAGMA
14
46
11
5.2 Evaluation on real dataset
We apply SparseRC on the causal protein-signaling network data provided by [achs et al. 2005 The dataset consists of 7466 samples (we use the first 853) for network with 11 nodes that represent proteins, and 17 edges representing their interactions. It is small, but the task of learning it has been considered a difficult benchmark [Ng al. 2020 Gao et al., 202T Yu al. 2018). It is not known whether the assumption of few root causes holds in this case. We report the performance metries SHD and SID for the most successful methods Table [ The number of total edges is used to ensure that the output of the methods are non-trivial (e.g. empty graph).
2019 [Zheng et al..
The best SID is achieved by GOLEM (equal variance), which, however, has higher SHD. NOTEARS has the best SHD equal to 11. Overall, SparseRC performs well, achieving the closest number of edges to the real one with 16 and a competitive SHD and SID. We also mention again the good performance of SparseRC on non-public biological dataset [Misiakos et al. 2023 Chevalley et al. 2023) as part of the competition by Chevalley al. (2022).
Broader Impact and Limitations
Our method inherits the broader impact of prior DAG learning methods including the important caveat for practical use that the learned DAG may not represent causal relations, whose discovery requires interventions. Further limitations that we share with prior work include (a) Learning DAGs beyond 10000 nodes are out of reach, (b) there is no theoretical convergence guarantee for the case that includes noise, (c) empirically, the performance drops in low varsortability, (d) our method is designed for linear SEMs like most of the considered benchmarks.
A specific limitation of our contribution is that it works well only for few root causes of varying location in the dataset. This in addition implies that the noise must have low deviation and the graph to have low average degree and weights of magnitude less than one. Also, in our experiments, we only consider root causes with support that follows multivariate Bernoulli distribution.
7 Conclusion
We presented a new perspective on linear SEMs by introducing the notion of root causes. Mathemati- cally, this perspective translates (or solves) the recurrence describing the SEM into an invertible linear transformation that takes as input DAG data, which we call root causes, to produce the observed data as output. Prior data generation for linear SEMs assumed dense, random valued input vector. In this paper, we motivated and studied the novel scenario of data generation and DAG learning for few root causes, i.e., a sparse input with noise, and noise the measurement data. Our solution this setting performs significantly better than prior algorithms, in particular for high sparsity where it can even recover the edge weights, is scalable to thousands of nodes, and thus expands the set of DAGs that can be learned in real-world scenarios where current methods fail.
10

---

References
A. Amrollahi, A. Zandieh, M. Kapralov, and A. Krause. Efficiently learning Fourier sparse functions. Advances Neural Information Processing Systems, 2019.
B. Aragam and Q. Zhou. Concave penalized estimation of sparse Gaussian Bayesian networks. The Journal of Machine Learning Research, 16(1):2273-2328, 2015.
K. Bello. Aragam, and P. Ravikumar. DAGMA: Learning DAGs via M-matrices and Log- Determinant Acyclicity Characterization. Advances in Neural Information Processing Systems, 35:8226- 8239, 2022.
Bühlmann, Peters, and J. Ernest. CAM: Causal additive models, high-dimensional order search and
penalized regression. The Annals of Statistics, 2266)22526,2556, 2014.
B. Charpentier, S. Kibler, and S. Günnemann, Differentiable DAG Sampling. In International Conference on Learning Representations, 2022.
M. Chevalley. Y. Roohani, A. Mehrjou, J. Leskovec, and P. Schwab. CausalBench: A sarge-scale Benchmark for Network Inference from Single-cell Perturbation Data. arXiv preprint arXiv:2210. 221.22280. 2022.
M. Chevalley, J. Sackett-Sanders, Y. Roohani, P. Notin, A. Bakulin, D. Brzezinski, K. Deng, Y. Guan, J. Hong. M. Ibrahim, et al. The CausalBench challenge: A machine learning contest for gene network inference from single-cell perturbation data. arXiv preprint arXiv:2308 15395, 2023.
D. M. Chickering. Optimal structure identification with greedy search. Journal of machine learning research. (Nov):507-554, 2002.
M. Chickering, D. Heckerman, and C. Meek. Large-sample learning of Bayesian networks is NP-hard, Journal of Machine Learning Research, 5:1287-1330, 2004.
H. Cramér. Über eine eigenschaft der normalen verteilungsfunktion. Mathematische Zeitschrift, 41(1):405-414, 1936.
Gao, L. Shen, S.-T. Xia. DAG-GAN: Causal structure learning with generative adversarial nets. In International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3320-3324. IEEE, 2021.
A. Ghoshal and J. Honorio. Learning identifiable Gaussian Bayesian networks in polynomial time and sample complexity. Advances in Neural Information Processing Systems, 30, 2017.
H. Hassanieh. The Sparse Fourier Transform: Theory and Practice, volume 19, Association for Computing Machinery and Morgan and Claypool, 2018.
A. Ikram, S. Chakraborty, S. Mitra, S. Saini, S. Bagchi, and M. Kocaoglu. Root cause analysis of failures in microservices through causal discovery. Advances Neural Information Processing Systems, 35:31 158- 31170, 2022.
M. Kaiser and M. Sipos. Unsuitability of NOTEARS for Causal Graph Discovery when dealing with Dimensional Quantitics. Neural Processing Letters, pages 1-9, 2022.
Kalainathan and Goudet. Causal discovery toolbox: Uncover causal relationships in python. arXiv preprint arXiv:1903.02277 022778 2019.
D. P. Kingma and Ba. Adam: A method for stochastic optimization. arXiv preprint a 2014.
Lachapelle, P. Brouillard, T. Deleu, and Lacoste-Julien. Gradient-based neural dag learning. International Conference on Learning Representations, 2019.
H.-C Lee, M. Danieletto, R. Miotto, Chemg, and J. T. Dudley, Scaling structural learning with NO-BEARS to infer causal transcriptome networks. In Pacific Symposium on Biocomputing, pages 391-402, World Scientific, 2019.
D. J. Lehmann. Algebraic structures for transitive closure. Theoretical Computer Science, 4(1):59-777) -6, 1977.
P. Lévy. Propriétés asymptotiques des sommes de variables aléatoires enchainées. Bull. Sci. Math, 59(84-96): 109-128. 1935.
P.-.. Loh and P. Bühlmann. High- diad-nddienionl learning of linear causal networks via inverse covariance estimation. The Journal of Machine Learning Research, (5())3305-3100 2014.

---

P. Misiakos, C. Wendler, and M. Püschel. Learning Gene Regulatory Networks under Few Root Causes assumption. OpenReview, 2023. URL http://gwggppppppeeeeas//////aaa
P. Misiakos, V. Mihal, and M. Püschel. Learning signals and graphs from time-series graph data with few causes. In International Conference in Acoustics, Speech and Signal Processing (ICASSP). IEFE, 2024. To appear.
L. Ng, A. Ghassami, and K. Zhang. On the role of sparsity and DAG constraints for learning linear DAGs. Advances Neural Information Processing Systems, 33:1 7994--17954 2020.
R. Pamfil, N. Stiwattanaworachai, S. Desai, P. Pilgerstorfer, K. Georgatzis, P. Beaumont, and B. Aragam. Dynotears: Structure learning from time-series data. In International Conference on Artificial Intelligence and Statistics, pages 1595-1605. PMLR, 2020.
Peters and Bühlmann. Identifiability Gaussian structural equation models with equal error variances. Biometrika, 101(1):219 228, 2014.
Peters and Bühhmann. Structural intervention distance for evaluating causal graphs. Neural computation, 27 (3):771-799, 2015.
Peters, D. Janzing, and Schölkopf. Elements of causal inference: foundations and learning algorithms. The MIT Press, 2017.
C. Ramirez, Kreinovich, and M. Argacz. Why & is a good approximation to to: A Geometric Explanation. Journal of Uncertain Systems, 7(3):203-207, 2013.
Ramsey, M. Glymour, R. Sanchez-Romera and C. Glymour. million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images. International journal of data science and analytics, 3(2):121-12,, 2017.
A. Reisach, C. Seiler, and S. Weichwald. Beware of the simulated DAG! causal discovery benchmarks may be easy game. Advances Neural Information Processing Systems, 34:27772-27784 2021.
K. Sachs, O. Perez, D. Pe' D. A. Lauffenburger, and G. P. Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308(572)) 529, 2005.
B. Seifert, C. Wendler, and M. Püschel. Learning Pourier-Sparse Functions DAGs. In ICLR2022 Workshop on the Elements of Reasoning: Objects, Structure and Causality. 2022.
B. Seifert, C. Wendler, and M. Püschel.. Causal Fourier Analysis Directed Acyclic Graphs and Posets. IEEE Trans. Signal Process., 71:3805-3820, 2023. doi: 10.1 11099TSP.2023333249999
Shimizu, P. O. Hoyer, A. Hyvärinen, and A. Kerminen. Linear Non-Gaussiaa Acyelic Model for Causal Discovery. Journal of Machine Learning Research, 7(72):2003-2030 2006. URL htttt://jjjmr..oo swperss/w7/shhiizu0.. html
Shimizu, Inazumi, Y. Sogawa, A. Hyvarinen, Kawahara, T. Washio, P. O. Hoyer, K. Bollen, and Hoyer. Directlingam: A direct method for learning linear non-gaussian structural equation model. Journal of Machine Learning Research-JMLR 12(Apr):1225-1244) 2011.
P. Spirtes, C. N. Glymour, R. Scheines, and D. Heckerman. Causation, prediction, and search. MIT press, 2000.
P. Stobbe and A. Krause. Learning Fourier sparse set functions. In Artificial Intellgence and Statistics, pages 1125-1132 133. PMIR, 2012.
M. Tsagris, C. Beneki, and H. Hassani. On the folded normal distribution. Mathematics, 2(1):12-28, 2014. ISSN 2227-7390. doi: 10.3390/aath2100012
httlhhtts://www mdpi com/ cwm.c227-ceo/ 2227 /390/2/1/13
L. Tsamardinos. L. E. Brown, and C. F. Aliferis. The max-min hill-climbing Bayesian network structure learning algorithm. Machine learning, (5(1):31-78, 2006.
M. J. Vowels, N. C. Camgoz, and R. Bowden. D'ya like DAGs? A survey on structure learning and causal discovery. ACM Computing Surveys (CSUR), 2021.
Y. Yu, J. Chen, T. Gao, and M. Yu. DAG-GNN: DAG structure learning with graph neural networks. In International Conference on Machine Learning, pages 7154-7163, PMLR, 2019.
Y. Yu. T. Gao, N. Yin, and Q. Ji. DAGs with no curl: An efficient DAG structure learning approach. In International Conference on Machine Learning pages 12156 12166.1166 PMLR, 2021.
12

---

V. Zantedeschi, J. Kaddour, L. Franceschi, M. Kusner, and V. Niculae. DAG Learning the Permutahedron. In ICLR2022 Workshop on the Elements of Reasoning: Objects, Structure and Causality, 2022.
X. Zheng, B. Aragam, P. K. Ravikumar, and E. P. Xing. Dags with tears: Continuous optimization for structure learning. Advances in Neural Information Processing Systems, 31, 2018.
X. Zheng, C. Dan, B. Aragam, P. Ravikumar, and E. Xing. Learning Sparse Nonparametric DAGs. In International Conference on Artificial Intelligence and Statistics, pages 3414-3425. PMLR, 2020.
13

---

A Few Root Causes Assumption
Lemma A.I. Consider random row vector gxxd of root causes that are generated as in Section 5.1 with
probability = 0.1 being non-zero and value taken uniformly at random from (0, 1). Also, let the noise vectors Ne, N2 1 be Gaussian Ne, N₂ N (0, aI) with 0.01 and the DAG has average degree 8(C) and weights aiy [-1,1]. Then show that for = 0.1:
S c (few rool causes),
|I + No (I A)
<
(negligible noise),
(12)
[[]]]]
Proof. From the Bernuulii distribution for the support of the root causes we get
B [|[||||] = pd d
(13)
Also the root causes take uniform values in (0, 1) and have expected value 0.5. Therefore,
[[[[]] :] = 0.5d
(14)
Finally, we have that:
E [ N₂ (I AIIL
R [[ + (T A) 16 S
(15)
E [e]b] + (8(C) + I)E [I]
6do
(16)
computation of the expected value follows from mean the Folded normal distribution 2014 Dividing the last two relations gives the required relation
sagris
[| + N. (I A)]l]
2
0.1
(17)
[|
51
B Additional experimental results
B.1 Different application scenarios
Experiment 1: Different application scenarios. We present additional experimental results include metrics baselines. results support our observations and conclusions from the main text. First, we expand reepotnng the SHD metric and provide the computation of SHD, TPR, SID and runtime in the tables below. Those further include the methods LiNGAM (Shimizu et al., 2006), CAM [Bühlmann et al. 2014). DAG-NoCurl (Yu al. 2021 (GES [Ramsey et al. 2017 sortnregress [Reisach et al. 2021 and MMHC
T'samardinos al. 2006
14

---

B.2 Varying number of nodes or samples
Experiment 2: Varying number of nodes or samples. In addition to Fig. 2 we include the plots of Fig. 4 for the experiments that vary the number of nodes of the ground truth DAG the number of samples in the data. The plots include the methods LiNGAM, CAM, DAG-NoCurl, fGES, sortnregress and MMHC and additionally contain the metrics TPR and total number of edges regarding the unweighted adjacency matrix and NMSE with respect the weighted approximation of the adjacency matrix.
B.3 Larger DAGs
Experiment 3: Larger DAGs. In the large scale experiment we don't report SID as is computationally too expensive for a large number of nodes. Since, our method achieved perfect reconstruction, we further investigate whether the true weights were recovered. Table reportss metrics that evaluate the weighted approximation of the true adjacency matrix. Denoting with [ the number of edges of the ground truth DAG (V, E). compute:
||AAA|
the average L loss
[E]
the Max-L loss maxis
Aij
Aus.
||A-Â]
the average L loss
and
|
the NMSE
A
B.4 LINGAM's performance
The subpar performance of DirecILiNGAM arises as a contradiction to the fact that (Shimizu et al. 2006 provide the identifiability result for our setting. A possible reasoning to this can be that the corresponding linear noise variances distribution. Fig we evaluate LiNGAM and DirectLiNGAM with respect ours, using default settings, with zero measurement noise N = O, which translates 44 to a linear SEM with i.idd noise. We conclude that the methods cannot recover the true graph even in these settings. The reasoning we give for their failure is the following: First, LiNGAM can stay a local optimum due to badly chosen initial state in the ICA step of the algorithm, even compute wrong ordering of the variables (Shimizu al. 2011 Secondly, DirM guaranteed to converge to the truc solution only the conditions are striclty met, among which the infinite number of available data, which is not the case in our experiments.
(5) are non i.i.d., where both LiNGAM and DirectLiNGAM consider i.i.d. noise
B.5
Real Data
the real dataset from (Sachs et al. 2005) we follow Lachapelle al. (2099) and only use first 853 samples from the dataset. Differences may occur between the reported results and the literature due to different choice of hyperparameters and the use of 853 samples, where others might utilize the full dataset.
Root causes in real data. To question whether the few root causes assumption is valid in the real dataset of [Sachs et al. 2005 perform the following experiment. Using the ground truth matrix in the (Sachs et al.)
Table SHD metric (lower is better) for learning DAGs with 100 nodes and 400 edges. Each row is an experiment. The first row is the default, whose settings are in the blue column. In each other row, exactly one default parameter changed (change). The best results are shown in bold. Entries with SHD higher that 400 are reported as failure and in such a case the corresponding TPR is reported if higher than 0.8.
Hyperparameter
Default
Change
LINGAM
CAM
DAG-NoCuri
IGES
sortnregress
MMHC
1. Default settings
285 1:
325 28
737 1: 94
547 86
800 136
failure
Graph type
Erdös-Renyi Gaussion
Scale-free
failure
45
36
failure
3. N., N, distribution
Gumbel
478 + 342
350 + 18
642 + 55
407 +
796 0 86
failure
4. Eidges/ / Vertices
10
1102 108
961 : 20
1455 277
failure
1950 97
failure
Samples
100
emor
out
failure
enor
emerr
failure
6. Stundardization
No
Yes
363 =: 38
329 + 34
failure
370 109
failure
failure
7. Larger weights in A
0.9)
(0.5.2)
800 = 129
failure
135 + 5G
failure
1275 # 133
failure
N., N₂ deviation
0.01
failure
29
failure
failure
53
failure
Dense causes C
p= 0.1
0.5
faiture
271 +: 21
586 58
528 72
746 134
failure
10. Fixed support
No
Yes
failure
390 1: 51
failure
failure
failure
failure
15

---

Table 6: TPR (higher is better). Entries with TPR lower that 0.5 are reported as failure.
(a)
Hyperparaneter
Detault
Change
SpansRCC
GOLEM
NOTEARS
DAGMA
DirecLNGAM
PC
GEES
Default eettings
1.00 + 0.00
0.82 0:00
0.92 :: 0:00
failere
0.59 0:00
fuilure
0.72 :t 0.04
2. Grophh type
Fodis--Renyi
Seale-fiee
0.99 + 0.00
0.96 000++0001
0.96 + +.002
0.50 + +.001
Falure
0.84 +
N., N, distribution
Caussian
Gumbel
1.00 0.00
0.82 0.111
0.92 0.01
failure
0.99 1.000
failure
0.81 0.022
Edges/ Veruces
0.00
0
0.50 0.01
failure
0.00
fallure
failure
5. Samples
1000
n- 100
0.90 + 0.02
fullure
0.76 + 0.02
failure
enar
failure
0.61 + 0.06
Stumdaddization
No
Yes
0.84 0.01
failure
failkrre
failere
0.90 1: 0.01
failure 0.78 1: 0.00
Larger weights
(0.1,0.9)
(0,5,2)
failure
0.83 0.04
0.01
0.53
0.00
fallure
0.70 0.00
8. N., N₂ devintion
0.01
0.1
0.88 0.01
0.85+00000
0.55 + 0.01
fuilure
0.97 +: 0.01
failure
0.70 + 0.04
Dense causes C
0.1
0.5
0.88 0.02
0.94 0:01
0.93 + 0:00
0.81 0.00
0.99 0.00
failure 0.50 + 0.05
Fixed support
failure
0.02
fallure
0.05
(b)
Hyperparaeeter
Default
Change
LINGAM
CAM
DAG-NoCCrr
(GBS
sortmieggeess
MMHC
1. Default settings
0,99 0.00
failure
0.84 0.02
0.80 1: 0.04
0.57 4: 01.02
0.58 0.03
Graph type
Rends
Scale free
0.99 0.01
failure
0.92 d: 0.02
0.76 +: 0.06
0.96 4: 0.01
fuilure
3. Ng.. N2 distribution
Gaussian
Gumbel
0.94 0.10
failure
0.82 + 0.02
0.82 + 0.03
0.88 + 0.02
0.56 156++00 0.04
Edges Vertiees
failure
0.05
0.54
0.90 0.01
failure
5. Samples
n = 1000
n = 100
etror
time-out
0.76 = 0.03
crror
error
failure
Standardization
Yes
failure
failure
0.80
failure
0.05
7. Larger weights in A
(0.1.0.9)
(0.5.2)
1.00 = 0.00
failure
0.92 + 0.04
0.74 0.06
0.91 + (1.01
failure
N.,. N, deviation
0.01
0.12
failure
0.78
0.76
0.85 0.01
0.02
Dense mont canses C
0.1
p 0.5
0.70 0.03
0.51 0.03
0.84 + 0.04
0.83 + 0.03
0.88 + 0.02
0.57 + 0.04
10. Fixed support
Yes
0.52 :: 0.04
failure
failure
0.54 - 0.03
0.68 + 0.04
failure
Table 7: SID (lower is better).
Hyperpanaettrr
Defaul:
Change
SparseKC (ours)
GOLEM
NOTEARS
DAGMA
DireILINGAM PC
GES
Defaaih vattings
51 =b 68
4005 1441
2720 + 7.00
0089 + 265
257 9-
7948 115
0:00 + 561
Gooph type
Fudii-Rayi
Seale-fie
72 42
442 +:2++11
425 50
1681 + +1172
98
=
2570 + 314
1052 +-11
N., distriborion
Gaussian
Gumbel
125 1 91
4218 10:111110
2945 1 60)
7007 1: 220
153
7966 1. 405
6967 377
Edges Vertioes
10
7174
156
Samples
1000
a- 100
SID timo-nut
7410 320
SID timovuul
timo-cul
SID timo-nst
SID timo-ul
Standedication
No
Yes
SID lims--uut
9881 :: 147
8887 + 201
88G1 + +185
208 +: 174
7801 +: 368
5268 :: 1078
Larger weights in
(..1.0.9)
(0.5. 2)
GUS8 30G
3165 549
2179 2: 2012
0624 4: 151
17 2: 34
SID time out
SID time
N. N, deviation
0.01
7140 411
3808 635
6388 255
7227 + 573
882 : 157
T811 305
6310 + 702
0.1
SID time-nut
200
102
SID time-nst
SID
10. Finod suppoct
Yes
7000
418
7707
77:00 d:
6800 601
d: 314
(b)
Hyperparameter
Default
Change
LINGAM
CAM
DAG-NoCurl
RGES
sortinegrees
MMHC
1. Default scttings
595 4: 212
818- : 476
5800 594
7706 4: 456
1118 1: 161
SID time-oul
2. Graph
Endris-Renyi
Scale-free
90 + 45
SID timc-oul
8U7 219
7412 + 493
427 + 83
SID time-out
3. N.. N, distribution
Gaussian
Gumbel
1796 :: 2754
8916 : 317
5723 575
9010 :: 170
1522 797
SID time-oul
Edges / Vertices
10
1516 + 234
95:00 + 30
G200 620
9057 + 172
4967 507
SID time-out
5. Samples
n= 1000
100
error
time-oot
4190 + 507
emot
arror
SID time-oul
Stanfardization
No
Yes
2847 + 577
8707 + 287
8065 68
7872 + 600
7044 124
SID time-out
7. Larger weights in A
(0.1.0.9)
(0.5.2)
16 - 33
81/7 :: 231
1660 781
7137 1 196
968 181
SID time-oul
N₂, N₂ deviation
0.01
U.1
7403 + 862
8804 + 168
6006 506
7842 +: 022
4715 01::::61
SID time-out
Dense canses
p-
p- 0.5
6944
8034 + 278
5008
7510 + 187
4072
time-out
10. Fixed support
No
Yes
6002 +: 619
7176 +: 578
7350 0:50++00 569
8632 353
5477 + 279
SID time-out
dataset we perform an approximation of the underlying root causes. As the known matrix is unweighted we choose to assign them the weights that would minimize the sparsity of the root causes [[C]] The root causes
2005
C that are computed based on this optimization are depicted in Fig. [aa] In Fig. 6b We see the distribution of the
values of root causes: there many values with low magnitude close
and few with significantly high
magnitude. Also many at roughly the same location (Fig. 6a
B.6 Implementation details
Training details. In the implementation of SparceRC we initiate PyTorch Module with single d x d weighted matrix represented as a linear layer. The model trained so that the final result will have a weight matrix which approximates the original DAG adjacency matrix. To train our model we use the Adam optimizer [Kingma and Ba 2014] with learning rate 10 10"3, We choose À = 10 as the coefficient for the L regularizer of optimization problem (00 For NOTEARS also choose À 1 The rest of the methods are used with
16

---

Table 8: Runtime [s].
(a)
Hyperparameter
Default
Change
SparsoRC (ours)
NOTEARS
DAGMA
DiReLLNGAM
Default settings
10 1.8
= 210
196 = 185
320 =: 94
103 ::
22
7.1
Graph type
Erdox Renyi
Scale-free
L.1
181
180 L.
258 49
108 1.08
115 58
5.2
distribution
4. Edges Vertices
10
141:10 1.0
317 121
471 L. 82
217 11
1691121
5,6 1.1
0..
Samples Stundardization
1000
13:
26
222
3.9+ 0.5
No
Yes
11 1.9
326 145
781 76
157
107 0.6
174.44
5.3
Larger weights in A
(..1.
(0.5,22
8.4 0.6
431 177
2834 = 228
493 = 46
117
2.5
21
deviation
Dense not couses
0.1
0.5
9.1
331 121
(27 L. 15
133 20
142 1.36
1.41
4.9 1.2
Fixed support
2.0
501
106
53=332
6.2 2.1
Hyperparanetor
Default
Change
LINGAM
CAM
DAG-NoCut
IGES
sorturegress
MMHC
1. Default settings
2.4 0.D
372 7.4
15
7.9
2.6 0.2
3.5 0.3
Gruph type
Endös Renyi
Scale
2.7
3. N.. Nr distribution
Gausssan
Gumbel
2.5
312 1.35
1.110
6.2
2.1 0.3
34 0.6
4. Edges Vertices
4
10
3.1 0.5
341 9.7
21
573 167
2.1 0.0
3.1
Samples
1000
100
error
error
6. Standardization
No
Yes
2.4 0.1
304 +
13+15
7.2
2.0
3.7
7. Larger weights in
(0.1.0.)
(0.5,2)
3.0
319 10.0
106
118 46
2.3 0.0
2-5-02
N., deviation
2.28
Dease root causes C
0.1
p= 0.5
6.8 0.7
306 7.4
111
6.8
2.0
311
10. Fixed support
No
Yes
5.5 0.6
342
52
6.1 6.145.1
2.0
0.6
Table 9: SparseRC weight reconstruction performance on larger DAGs.
Nodes d, samples
Avg. I loss
Max-T,
loss
Avg. I.2 loss
NMSE
200, 500
0.071
0.317
0.003
0.087
500, n 1000
0.066
0.275
0.002
0.079
d 1000, 5000
0.050
0.287
0.001
0.060
2000, " 10000
0.050
0.388
0.001
0.062
3000, 10000
0.054
0.399
0.001
0.067
default parameters, except for the post-processing threshold which is always to 0.09, except to the case of larger weights (Table row 6) where is set slightly below weights 0.4.
Resources. Our experiments were run on a single laptop machine core CPU with 32GB RAM and an NVIDIA Gelorce RTX 3080 GPU.
Licenses. We use code repositories that are open-source and publicly available on github. All repos-
itories licensed under the Apache 2.0 or MIT license. In particular we use the github reposito-
ries of DAGMA [Bello al. 2022
GOLEM (Ng et al. 2020
https///iiuub.comiinaaiiiei
NOTEARS [Zheng al. [2018] https://github.com/// htggghuubb.oom//uunnenngtnttttnt
DAG-NoCurl Yu et al., 2021] https://gihub..om/iihhh...
234/DAG-NoCut LiNGAM [Shimizu
et al., 2006) https://ittub... com/edt 5//inggm the implementation of sortnregress RRisach et al. 2021
http://guuu..co/Srriddtt Varsortability and the causal discovery toolbox (Kalainathan and Goudet
2019
https://ittub.coo http:/ggthub.com/Fnnn (eccisollliinn//Caaaaaassooooo
Our code is licensed under the MIT license and
available publicly in nighhbb
C
Global Minimizer
We provide the proof of Theorem 3.2 of the main paper. The data X assumed he generated via noise-less root causes described. We denote by A the optimal solution of the optimization problem S We want to prove the following result.
Theorem C.1. Consider DAG with weighted adjacency matrix A with d nodes. Given exponential (in d) nomber samples X the matrix A with high probability, the global minimizer of the optimization problem
In particular what will show the following result, from which theorem follows directly.
17

---

500
802
403
SHD
600
Time (6)
300
400
200
100
Number nodees
Number nodes
Number nodes
500
400
0000
400
SparseRC
300
Time
300
NOTEARS
200
203
DAGMA
100
2003
100
DirectLINGAM
0 1000200
500
800
1000
100200
500
500
1000
102200
500
500
1000
GES
Numbar samples
Number of samples
Numbar of samples
LINGAM
CAM
(a) SHD (1)
(b) SID (L)
(c) Runtime [s]
DAG-NoCurl
1.0
1300
IGES
sortnregress
1000
MMHC
800
NM55
20-
600
0.4
400
300
02
160
20
40
60
Number nodes
Number nodes
Number of nodes
1.0-
1000
0.8
2500
g
0.7
1100
1000
0.3
150
sčo
10200
1000
150
850
samples
Number samples
samples
(d) TPR (t)
(e) Total edges
(f) NMSE (L)
Figure 4: Plots illustrating performance metrics (a) SHD (lower is better), (b) SID (lower is better), (c) Time [seconds], (d) TPR (higher is better), (e) Total number of proposed edges and (f) NMSE (lower better). Each metric is evaluated in two experimental scenarios: when varying the number of rows (upper figure) and when varying the number of samples (lower figure).
400
2000
Sporwit:
200
1000
200
Drecdt NGAM
100
0-
20
40
60
80
100
20
40
60
80
100
40
62
80
100
Number of nodes
Number of nodes
Number nades
(a) SHID (1)
(b) SID (J)
(c) Runtime [s]
Figure 5: Evaluation of LiNGAM's performance when N = 0.
Theorem C.2. Consider a DAG with weighted adjacency matrix A. Given that the number of data satisfies
233
2 (1- (1-p)a-
(18)
where k = [dp] and
2
In
2 In
(19)
p+(1
p)
18

---

Columns (nodes)
Count values
2500-
2000
1500
1000
500-
Rows (samples)
500 1000 1502000250000033000
|cyl
(a) Illustration of C
(b) Histogram of values of C
Figure 6: Root causes C in Sachs et al. 2005 dataset. (a) Root causes illustration: darker color means higher value. (b) Histogram of values of C.
Then with probability (1 c) the solution the optimization problem coincides with the true DAG, namely A A.
Remark C.3. The reason we choose & = [dp] is that because of the Bernoulli trials for the non-zero root causes, the expected value of the cardinality of the support will be exactly [dp] Thus expect have concentration that value & which something we desire, since in proof we use all possible patterns with support equal to k.
Proof sketch. Given large (exponential), but finite number of data increases), pairs of sub-matrices of rool causes matrices C (true) and C (corresponding the solution of optimization), with the same support and covering all possibilities of support with & non-zero elements, due to the assumption randomly varying support. Exploiting the support the pairs iteratively prove that the entries of the matrices A (true) and À (solution of the optimization) equal. The complete proof provided next.
find, with high likelihood
(approaching
We begin with some important observations and definitions.
Lemma C.4. If A
then =
Proof. We have that
I+Ā=++ (I+ = =[1+ (I+ (I+ - -I- -1- A
(20)
Definition C.5. We define Č= X XA the root causes corresponding the optimal adjacency matrix.
Definition C.6. Let S C (1,2,3, ..., df be set of indices. We say that root causes vector (Ca, ... Ca) support if
for E S.
Definition C.7. For given support S we consider the set R C [n] of the rows of C that have support S. Then, CR, ÖR denote the submatrices consisting of the rows with indices in R.
Lemma C.8. For any row subset R C [r4] we have that runk (c)) rank (Cx)
Proof. We have that
è (I+ = -X -
(21)
Therefore, since both A, Å are acyclic and I A, + Ā are invertible we have that
CR (I+ A A) = Cr (I A) + rank (CR) = rank (Cn)
(22)
Lemma C.9. For any row subset R. C [n] of root causes C with the same support S such that |S] and [R] 2 k the nonzero columns of CR are linearly independent with probability and therefore
rank (Cn) = k.
19

---

Proef. Assume Ct, .t..... the non-zero columns of CR. Then each is vector of dimension whose each
entry is sampled uniformly random in the range (0. 1]. Given any k 1 vectors out of Ci ...: C2 their linear span can form a subspace of (0, 1]" dimension at most & However, since every e, is sampled uniformly at random from (0, 1J" and r 2 > - 1 the probability that c, lics in the linear span of the other vectors has measure O. Therefore, the required result holds with probability 1.
Lemma C.10. With probability (1 =) for every different support with |S] there are rows R such that CR have support and
[R]
2d 23d-22((d 1)
(23)
(1-6)' Also set N be the random variable representing the number of root causes with ssupp(c)] [c) = & and N; the number repetitions of the i-th l: kusupport pattern, à 1,... ( We first use conditional probability.
Proof. Set ( 2dd- 1) and K
IP (N 2 K N.N, 2 1, vi) = P (NN: 2 4, WIN 2 PN 2 K)
(24)
Now we will show that P (NV 2 K) 2 (1 <) and P (N: 2 1, W 2 A) 2 (1 - c) using Chernoff bounds.
From Hoeffding inequality we get:
(N - (1 87₂) Se nn
(25)
where =
=mpp(1 p111-)p-=((0 is the expected value of N, 8 2 p1 " V In (2) and n 2
K
so
11 (1-p))a-k
S (1 5)m) S - pn (N
(26)
Now given that we have least K root causes with support exactly k, will show that, exactly K such root causes, each different pattern will appear least | times with probability (1 - c). From the union bound
P (UM
EP(NS
(27)
So we need to show that TP (N; (N) S We use again the Hoeffding inequality.
P(M S 6))) Se - ///2
(28)
where the expected value is =
K
I-8 since now the probability is uniform over all possible s-spsist
patterns. Given 8 2
/ In
we derive
PNN. 5 = P(N (1 Se 2/
(29)
The result follows.
Lemma C.11. Assume that the data number is large enough gh that for every different support S we have the rows R such that On have support S and [R] 2dI where
1> 2d d(d- 1) =) 2d Σ k(k 11)
(30)
Then there exist rows Rt such that [R] > k. Cg has support and C has support S with I$ k. Moreover. both the non-zero columns of Cg. CR form two linearly independent set of vectors. In words, this means that the data are many enough so we can always find pairs of corresponding submatrices of root causes with each consisting of k non-zero columns.
Proof. According to the assumption there exists row-submatrix CR where all rows have support R and [R] > 241, Group the root causes C₁ into all different supports, which are 2d in number. Take any such group R' C I of rows of Ò that all have the same support S'. If [R'| then from Lemma rank (Öx²
rank (Cr) & so CR will have least non-zero columns, and since the support is fixed, it will have least k non-zero elements each which means at least many non-zeros as Cg. Therefore CR can only
20

---

have less non-2ero elements if [R'| < k. and in that case Cx' has most k(A - 1) more elements, If we count for all & = I.... d all different supports of Czy for all possible supports of Cn this gives that C can have at most Ei= kk) - 1) less non-zero elements compared to C.
Due to the pigeonhole principle, there exists ft R and [it] with Cg having the same support S, not necessarily equal to S. According to our previous explanation we need to have at least & non-zero columns in Cg. we had k columns then this would give more non-zero elements, but
d=2 d
1> 2dd d(d =
d(d
2a
k(k- -1 (31)
k 2
k=00
So then |/ > ||C|| which is a contradiction due to the optimality of A. Therefore, C2 has exactly l: non-zero columns which necessarily need to be linearly independent in order have rank C rank (C k as necessary. The linear independence the columns & follows from Lemma
C.9 since k.
Definition C.12. A pair CR, CR constructed according to Lemma C.1 that have fixed support S, S respectively with cardinality & each and |tt large enough so that rank (Cn)) raan (Cn) = k, will be called a k-pair of submatrices.
Remark C.13. For notational simplicity will drop the index R whenever the choice of the rows according the sparsity pattern S is defined in context.
The complete Theorem 3.2 follows after combining the following two propositions C. and C.18 Before, the proof of the first proposition need the following lemmas.
Lemma C.14. A is (strictly) upper triangular if and only if A is (strictly) upper triangular.
Proef. If A is (strictly) upper triangular then it is straightforward to show the same for A = A A +... Ad-1, d-1 For the other way round, remember that the inverse of upper triangular matrix (if exists) matrix. However, the inverse of I + A is I A which means, that if A is strictly upper triangular, then I A and so is I A are upper triangular which in turn reduces to A being strictly upper triangular.
upper triangular
Lemma C.15. Consider C, C a l papii.. Also let X be the corresponding submatrix of data. If X = 0 then
C = 0 and āz = 0 Vj E supp
(32)
Proof.
0 0= x
Σ az C.u CC.
Σ
C.tC.. + C..i
(33)
jEsupp(C)
C. * O then E supp and the expression above constitutes a linear combination of the support columns with not all coefficients non-zero. This contradicts Lemma C.9 Thus C., = and ā5. = 0 Vj
supp
We are now ready to prove Prop. C. 16
Proposition C.16. If the data X are indexed such that A is (strictly) upper triangular, then so is A. Remark C717. Based on lemma C14 in proof of Prop. C.16 next may only about the transitive
weights, namely we will only need to show that ājz- = 0 for all <.j.
Proof. We first choose a p-pa C, Ò such that the support C concetrated to the last columns, namely C. = 0 for i 1,.... d k. Then the corresponding data submatrix X will necessarily have the same sparsity pattern since the values are computed according to predecessors, which in our case lie in smaller indices. Thus
X 0 Vi 1,..., d & and necessarily C will follow the same sparsity pattern. Moreover, according to
Lemma C.15 we get that
āyz = 0, for all S d l: and d k++ 1
(34)
We notice that the desired condition is fullfilled for = - k. i.c. it has zero influence from a node with larger index. We now prove the same sequentially for = d - k - 1, d k - 2, ..., 1 by moving each time (one step) left the leftmost index of the support S. We prove the following statement by induction:
P(I) a..- 0 for
(35)
21

---

We know that P(d - k) is true and P(1) gives the required relation for all indices i = 1,..,, d k. Now
assume that P()) holds. If we pick h- p-pair C, C such that C has support S = {1, d + 2, d - k + 3..., d).
Then ; 0 for all i < which with Lemma C.15 gives C. 0 for I which means supp C
Note, that by the induction hypothesis we have that ر,ة = 0 for all < j. However it is true that X C,, and also
[1, + 1, 1.... d) and ā5z 0 for j E supp C))
X Σ
C., + +.,, C.
(36)
jEsppp C))
Therefore I supp (c) and thus âr all which combined with P()) gives
نز 0 for < J. j
(37)
which is exactly P(I - 1) and the induction is complete.
Now remains show that āj: = for + 1 1 i and < j. We will again proceed constructively using induction. This time we will sequentially choose support S that is concentrated on the last k: + columns and at cach step we will move the zero column one index to the right. For / = d - k + 1,.., d let's define:
Q(1) ā,z for < and ررو Jee for d k S <1
(38)
First we show the base case Q - 1: 1). For this we choose k pair C, C such that C has support S = {d - k, k 2, d k + 3.... d}. It is true that X.,. 0 for i = 1,... d - ko 1 hence C.,, for doesn't have any non-zero parents, since also (s(a-)) for +
i S k and therefore node
< j from previous claim. Therefore C. = 8 = C.d-E. Also, for |
X C. C.
(39)
The equation from Å gives:
X
Σ C.C. + C. (āa dd) C+
Σ C.C. CC, j=d-k+1
(40)
j=d-k
From the linear Independence Lemma C.9 of the support of Ć necessarily need to have C. 0, ad-Ri ād-k, and ā= = 0 for I < j which gives the base case.
For the rest of the induction we proceed similar manner. We assume with strong induction that all Q(d - k (()..)))) are true and proceed to prove Q(I 1). Given these assumptions we have that
dz. J. for d - S d -ikkiis and az. f < j. d kiis (d k.d k 1,... + 2..... d) (the
(41)
Consider root
support S
(1 + 1)-th column
0) for the k-pair C, C. Then
have the equations:
Č.
X d4 Od- C.a-k =
C. = C.d-a-1 C.d-141
X.1 = Σ C., + CC =
+c,
+ č.,
j-d-
j=d-k
Where used the linear independence lemma and sequentially proved that the root causes columns up / are equal. The equation for the (1 1) column becomes:
X
Σ C. + C. Σ j=d-k
T + C.
(42)
S=d-k
ā,1+1 = 0 for 1+1 <<j
(2 - Čng + C.441 0
4,141 = J,,141 for j < < ++1
(43)
C =
where the last set of equalities follows from linear independence. This concludes the induction and the proof.
L
22

---

To complete the proof of Theorem 3.2 it remains to show the following proposition. Proposition C.18. !f both A, A are upper triangular then À = A. For our proof we use the following definition.
Definition C.19. We denote by Pa the set of all l- k-pairs Cn: On for all possible support patterns.
Now proceed to the proof.
Proof. We will show equivalently that A = A using two inductions. First we show for I 1,,,, & the following statement.
P(() For all p-pairs C, Ò in Pk the first I non-zero columns C..,, Catgri.. Cii,, and C.G,, C. C,ui are in the same positions, i.e. in ig and
either they are respectively equal C..,, = Č C. is in the last possible index, namely i
(1 (1-1)
For the base case P(1), consider a k-pair C, and let ij be the position of the first non-zero root causes column of C. Then x = 0 for i < in and therefore from Lemma C.15 C ; for < i Hence
X
Σ
ā5a C. + C.,, C.
(44)
J-31
Therefore C,, C.,51 and we proved P(I). by satisfying both the positioning and the first requirement.
Assuming now that P(() holds, will show P(l + 1). Take any k-pair of Pa which denote C, Č Then, if C..,, is the last possible position, then necessarily Cr is in the last possible position. Moreover, from the induction hypothesis the first root causes columns are the same positions. Therefore in the same manner, C,,,, is in the last position and C. too. This case fulfills the desired statement.
If C is not in the last position, then from induction hypothesis, the first root causes columns equal. If C, is in the last position and the same holds C. the requirement satisfied. Otherwise 441 < i+t1 and the equation for 4441 gives:
X X Σ Cagin C.,,+ Oimm
Σ anyiinn C,,
Σ (anseen ayyhoo)) C.G, + C
According to linear independence Lemma C.9 we necessarily derive C.iy = 0, absurd. Thus itf1 4441 and the induction statement is fullfilled in that case.
It remains to consider the case where Crice is not in the last position. Since, 4441 is not the last position there exists k pair C², O' such that the column i441 is zero and the (1 + 1)- root causes column of C' lies at ب i441. The equation at 1 for X' gives:
Σ Tajan Ciay = Dagan Ci., + + C
(45)
For the pair C', C' the induction hypothesis holds, thus we derive:
Σ Taiyteen ag.te1)) C'e,, + C =0
(46)
Therefore, Lemma
aly.te for all j = 1, 1,..., I and returning back to the equation for C.9 gives a 3e41
X X. we derive:
X
X
T d1y Cuay, + C.in =
11iie C.,, +
C
(47)
P(U)
C.e,,
C
(48)
j=l
C Ci
C
(49)
23

---

which completes the induction step. Notice that P(k) gives that for all io B--pair, the & root causes columns are in the same position. Given this fact we will now show that A A. We will prove by induction that for 1 == k 1,..., 1.0
Q(1 ن Tiy for Si < d-1
(50)
To prove Q(k 1) choose all the k pairs C, C, such that C..b is in the last possible position, 12 = dd +
Then for i < d - & the columns i. C..i, lie at the same position and are equal. Choosing 41 i and computing the equation for X. where j
k 1 gives:
X J 3. C = ā, = ā. C..
(51)
Therefore م- = asy for all 1 < jiijj < d (k - 1) and Q - 1) is satisfied. Next, assume that Q(A - I)
truc. We want to show Q(k ((k-1-1)) 1). Similarly to the base case we consider all k-pairs such that C.,42 lies in its last possible position 1r42 d + + 2, and it-1 d le + Since the + 1) column is not in the last position, from the previous induction have that:
C C..
C C..i2
(52)
Ch = C..411
The equation for + 1+1 gives:
1+1
X a--t-11
az d k+1+1 Cauy
(53)
1+1
Bzy d ++++1 Ca.y
(54)
Ci ++++++1100 Ci.y
(55)
0z ktie1 d
(56)
By choosing all such possihle k pairs the indices i, span all possibilities <i<d this with Qik - 1) we get that Q(k: - 1) is true and the desired result follows.
d d + 1+ 1. Combining
24


Reviews:
[{"review_id": "bScx9WfsAU", "reviewer": "Reviewer_NZKD", "summary": "This paper presents a new formulation of linear SEM by specifying a structure on the noise variables, which seems to impose zero-inflated distrbutions to achieve the \"few roots\" modelling goal. Identifiablity is given and guarantee for $L^0$ minimization estimator is provided for a special noise-free case. The minimization problem is further formulated into continuous optimization and numeric experiments are conducted.", "questions": "- Thm 3.1 is proved by transforming the proposed model into a non-Gaussian DAG model. Just want to clarify: is the proposed formulation covered by the orginal DAG model? If so, are all other DAG learning approaches applicable? In that case, how do the authors justify the superiority of the new formulation beyond the empirical experiments?\n- In Figure 2(a) and (b), the SHD and SID for mobius are larger than others at the beginning but drop in the end, any explanation on that?", "limitations": "Lmitations are discussed in Section 6.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The idea of the new formulation is interesting, based on the illustrative river pollution example.\n- The experiments consider many different setups, also including real datasets. The proposed algorithm shows comparable performance in the simulation study, though not by much.", "weaknesses": "- The motivation of the proposed formulation needs more elaboration, see questions. And the current definitions of \"few roots\" and \"negaligible noise\" in (6) are not formal or clear.\n- Thm 3.2 only works for noise-free setting, thus there is no guarantee for consistency of $L^0$ minimization estimator for general model (4) beyond noise-free.\n- Neither the related work and experiments discuss or compare with constraint-based methods, even the basic PC algorithm. Is it becaus we lose Markov property and faithfulness in this new formulation?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Response to reviewer NZKD\n\n## Strengths\n\n**Improvement against other methods.** Our algorithm offers significant improvements over the baselines in most scenarios of the simulation experiment (Table 1 and Fig. 2) and especially when the number of nodes in the graph is scaled up (Table 3).\n\n## Weaknesses\n\n**Motivation of our setting.** Please read our general reply for a better explanation of the motivation and also our success in a causal discovery competition on real-world data.\n\n**Assumption on the noise and root causes.** Our assumptions are formally defined in Eq. (6). In words, the root causes are approximately (as captured by $\\mathbf{N}_c$) sparse and the data is (as always) subject to measurement noise $\\left(\\mathbf{N}_x\\right)$, where both noise terms are of low magnitude as described in (6). In the concrete experiments we consider concrete numbers for the error energies. We will clarify.\n\n**Theorem 3.2 and the presence of noise.** Indeed, we have included this statement as limitation (b) in section 6.\n\n**Constraint-based methods.** We did not include constraint-based methods since prior work already has shown to outperform them. So we focused on the best baselines. However, we include a comparison against the PC algorithm in the attached pdf, Fig. 1. The PC algorithm performs inferior.\n\n## Questions\n\n**Superiority of our formulation.** Yes, our model still assumes a linear SEM via Eq. (5). Thus, it differs only in the assumption on the distribution and structure of what was viewed before as noise, which we view now as the input provided by each node to produce the DAG data as output (Eq. (6)). In this new formulation we assume (approximately) sparse input and the optimization problem we solve captures that. Thus the proposed input-output viewpoint of the linear SEM and motivates the assumption of having sparse root causes.\n\n**Low data regime.** Yes, it appears that in the low data regime our method performs worse. This is likely due to the sparsity of the root causes (which have varying support by assumption): to identify all edges the root cause supports of the data vectors may need to cover a significant subset of the nodes.\n"}, {"review_id": "ROsg7cJSEf", "reviewer": "Reviewer_STdi", "summary": "This paper considers learning of linear SEMs (weight matrix) under a data generation process that differs from the common formulation. It is assumed that each sample is generated from only few number of non-zero noise variables in which the set of noise variables is stochastic. The main theoretical result is identifiability of the weight matrix via an L0 norm objective. A relaxation of the objective with an L1 norm loss is used for extensive synthetic experiments and shown to be effective for learning the weight matrix under the proposed data generation mechanism.", "questions": "\n- The authors can elaborate on the Weaknesses item 1.\n\n- An additional question for the presentation. In L103-105,  “The high-level intuition is that it can be reasonable to assume that, with the viewpoint of (3), the relevant DAG data is triggered by sparse events on the input size and not by random noise.” \nAt first, I thought that you only consider a fixed subset of the nodes that have “random noise” with large magnitudes at each sample. If that was the case, the data is not i.i.d, the events are sparse on the input size, but it doesn’t necessarily occur on the same set of variables and I was puzzled with the importance of the proposed mechanism. Then I realized in Theorem 3.1 that at every sample, the non-zero noise variables are chosen randomly and it made sense (please correct me if I misunderstood anything). For a paper that proposes a new data generation model, the presentation should have been cleaner.\n\n- Theorem 3.2.. Given a large enough but finite number of samples, but not a sample complexity result. This is a rather weird statement.\n\n- Note: figure 2 has Möbius as the algo name. Perhaps it’s forgotten in the main text.", "limitations": "The limitations of this work are discussed very clearly in Section 6. I thank the authors for that paragraph.\n", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- Sparse root causes are an interesting concept to explore.\n- Experiments on synthetic data demonstrate strong performance. Specifically, the proposed algorithm is fast and seems to scale up well, and when the assumptions are violated slightly, the algorithm still yields reasonable results.", "weaknesses": "\n- While the sparse root causes assumption is interesting, it is not well-motivated. For the pollution example, if I understand correctly, the top left and bottom left nodes do not really play “causal” roles in the sense that they are deterministic mediators in the system, and the whole causal system can be represented without using those nodes. In this sense, “few root nodes” become the effective causal nodes. Am I missing something? The numerical evaluations using real data also do not provide much help for motivation. \n- NNZ is not a strong metric when SHD is already provided, and the only edge of the proposed algorithm in the real data experiments is this metric.\n- Presentation can be improved. For instance, Theorem 2.1 formulation of the linear SEM is trivial. Similarly, having $N_c$ and $N_x$ separately is superfluous; having a single small variance, not necessarily isotropic noise suffices for the description.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Response to reviewer STdi\n\n## Weaknesses\n\n**Motivation for root causes.** Please see our general reply for a better motivation and also our success in a causal discovery competition with real-world data.\n\n**Root causes support.** Note that we assume varying support of the root causes in the dataset $\\mathbf{X}$. The example in Fig. 1 shows only one data vector with a specific root cause support of size 2. Because of the varying support, the system cannot be simplified and all nodes are significant, possibly in distinct samples (data vectors).\n\n**NNZ metric.** We will adjust the text to only refer to SHD, SID as indicators for the performance of our algorithm. However, NNZ is valuable as an indicator that the output of the algorithms is not trivial. For example, the empty DAG would result in an SHD of 17 in Table 4 (since the true DAG only has 17 edges). We will clarify this in the text.\n\n**Anisotropic noise.** Yes, one can fuse the noise variables $\\mathbf{N}_c$ and  $\\mathbf{N}_x$ into a single (anisotropic) one:\n$$\\mathbf{N} =\\mathbf{N}_c + \\mathbf{N}_x\\left(\\mathbf{I} - \\mathbf{A}\\right).$$\nHowever, we prefer to keep them separate since they have distinct and intuitive meanings in a real-world setting: $\\mathbf{N}_c$ captures approximately (not exact) sparse root causes, $\\mathbf{N}_x$ is the measurement noise always present when obtaining real data $\\mathbf{X}$. Also, this allows us to manipulate them separately later in our experiments (Table 1, rows 9, 10).\n\n## Questions\n\n**Fixed support of root causes.** Yes, we have to clarify earlier that the root causes have varying support. We only consider fixed support in one of the synthetic experiments (Table 1, row 13).\n\n**Sample complexity.** Indeed not well written, we will rewrite. In the supplementary material we provide a bound on the number of samples in the noise-free case. Empirically, our algorithm works with much less samples and in the presence of noise, which shows that this bound is very loose.\n\n**Möbius.** It is a mistake, we'll fix.\n"}, {"review_id": "vRiq0CWyev", "reviewer": "Reviewer_LucQ", "summary": "This paper considers a new setting of linear DAG learning problem. Based on a linear transform of linear SEM, authors propose to study a new setting where there are few \"root causes\", with potential measurement noise in the data. Identifiability is proved and the true DAG is shown to be the global minimizer of the L0-norm of the vector of \"root causes\", under a specific distribution on the \"root cause\" variables. ", "questions": "- My first concern is about the new setting of learning linear DAGs; it is not clear whether the new setting is indeed meaningful in practice. In the pollution model example, it is stated that \"the relevant DAG data is triggered by sparse events on the input size and not by random noise \", and \"We assume a DAG describing a river network. The acyclicity is guaranteed since flows only occur downstream. .. We assume that the cities can pollute the rivers.\" In this example, why do we need to learn DAGs? The graph structure can be more accurately obtained by getting the information of flows. As such, I suggest authors give more practical examples in the context of DAG learning, to make the new setting indeed meaningful. \n\n- root causes: in the DAG learning literature, \"root causes\" generally refer to the source nodes of DAGs. Not sure if it is suitable to use a (somewhat) conventional name to refer to something new in the same context. \n\n* Theorem 2.1 is not new and may be not stated as a theorem. \n* Regarding Thm 3.1: similarly, the result simple follows from the LiNGAM result, by assuming a specific distribution on the \"root causes\", so maybe consider put it as a lemma or proposition. Besides, in the experiments in the supplementary material, I can see LiNGAM failed. Can you explain why? After all, the linear SEM falls exactly into the setting of LiNGAM if there no measurement noise. \n* after Eq. 8, \"Among all possible DAG matrices, the solution of the optimization problem (8) is the one that minimizes the number of the root causes X\": can you give more details about this claim? \n* This may be a bit picky, but only sparse graphs (with edge/node=2 and 3) are considered. Please try other degrees of graphs. (But this is not very important and may be added after the rebuttal.) \n* please use \\cite, \\citet properly; e.g., line 131 line 181-182\n\nOverall, I like the new setting of learning linear DAGs, but every new setting should be validated with more examples/details. I look forward to author response.", "limitations": "Authors discussed a limitation that the proposed method only works for few root causes in the paper. To me, another important limitation is the specific distribution assumption on the \"root causes\", as the proof heavily depends on this assumption. ", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "- a new setting for the linear DAG learning problem\n- useful identification result (Thm. 3.2) with a complete proof\n\n", "weaknesses": "- the new setting and its motivating example are not sufficiently convincing. \n- authors only consider specific distribution on the \"root causes\" variables, making  theoretic result somewhat limited\n- some results are trivial from the literature (e.g., Thm 2.1) ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Response to reviewer LucQ\n\n## Weaknesses\n\n**Motivation for root causes.** Please see our general reply for a better motivation and also our success in a causal discovery competition with real-world data.\n\n**Assumptions on the root causes.** Yes, within the scope of sparse root causes we only consider the special case that the support of the root causes is a multivariate Bernoulli variable (since the values are either zero or nonzero) and we should mention it in the limitations. We considered it a reasonable assumption to execute our ideas. Identifiability will hold for larger classes of distributions but we did not run experiments.\n\n**Theorem 2.1.** We will call it a lemma. It is indeed an easy consequence of basic linear algebra (and thus not novel of course) but usually not explicitly stated in papers on linear SEMs. We do so here since this input/output view of a linear SEM is fundamental to our work.\n\n## Questions\n\n**Meaningfulness in practice.** Please see our general reply for a better motivation and also our success in a causal discovery competition with real-world data.\n\n**Term: Root causes.** Source nodes are sometimes called root nodes but the term root causes seems not common. In any event, we will better clarify better its use in our work to avoid ambiguity.\n\n**Theorem 2.1.** See reply above.\n\n**Theorem 3.1 and LiNGAM's performance.** Yes, Theorem 3.1. is a consequence of the input being non-Gaussian, which we also state. Better to call it a corollary to the LiNGAM result than a lemma.\n\nThe subpar performance of LinGAM is an interesting question and a possible reasoning can be found in [Shimizu et.al, 2006], Section 5. Theoretically, the LiNGAM algorithm is guaranteed to find the true DAG. However, in practice, there exist small estimation errors. Moreover, exhaustive search over all possible permutations to compute the true causal order is infeasible for a large number of nodes. Thus, it turns out that the algorithm approximates with zeros the smallest elements and then computes the corresponding permutation to make it upper triangular. This makes the algorithm approximate and it is unclear whether it is guaranteed to find the true matrix in practice. For example, as we see in our experiments, it can fail.\n\n**Optimization Objective, after Eq. (8).**\nFor a DAG matrix $\\widehat{\\mathbf{A}}$ the quantity $\\mathbf{X}\\left(\\mathbf{I} + \\overline{\\widehat{\\mathbf{A}}}\\right)^{-1}$ is equal to the root causes $\\widehat{\\mathbf{C}}$ that would generate the data $\\mathbf{X}$ via equation (7). Therefore, the optimization objective in (8) minimizes the overall number of root causes (number of nonzero entries of $\\mathbf{C}$) for the entire data matrix $\\mathbf{X}$. We will write this better.\n\n**Edge density.** We chose to follow the experimental settings of prior work, which mostly generates sparse DAGs. For completeness, we conducted an experiment with varying average degree included it in the attached pdf, Fig. 2. For degree up to roughly half the maximal possible degree, we still perform best.\n\n**Citations.** Thanks, we will fix it.\n\n## Limitations\n\n**Root causes distribution assumption.** Yes, we will mention in the limitations that we only considered multivariate Bernoulli. Identifiability will generalize to other distributions (due to the non-Gaussianity needed in Theorem 3.1) but we did not run any experiments with any. We did mention in the limitations that fixed support is not supported (and it seems also not identifiable).\n"}, {"review_id": "TIRBN1qP2a", "reviewer": "Reviewer_Aevs", "summary": "The paper studies a new causal discovery method, in which it assumes that the DAG data is produced by few data-generating events whose effect percolates through the DAG. They propose a simple but effective method to learn the true DAG based on the few roots assumption. The proposed method outperforms baselines in various settings.", "questions": "Please reply to the questions in the weaknesses. I am willing to raise the score if all the concerns are solved.", "limitations": "Yes", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The few roots assumption is reasonable. And the paper motivates it well.\n2. The proposed solution is simple and effective.\n3. The paper conducts experiments on both synthetic and real-world datasets, indicating the effectiveness of the proposed method.\n", "weaknesses": "1. It is not clear why the objective function Eq.(10) contains noise. It needs more clear explanation and derivation here.\n2. I wonder whether the proposed method could find the root nodes at the same time rather than just learning the DAG.\n3. The proposed method could not achieve the best results on real-world datasets. Hence, I doubt the few roots assumptions satisfy the real-world scenarios. BTW, it is better to say the network is a protein network rather than a gene network.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Response to reviewer Aevs\n\n## Weaknesses\n\n**Optimization robustness to noise.** The optimization problem doesn't contain noise explicitly. It is the convex  $L^1$ relaxation of the noise-free version of the optimization problem. Doing this relaxation allows some robustness to (low magnitude) noise as we later confirm in our experiments (Table 1, row 7).\n\nThe objective in Eq. (10) only contains noise implicitly, after rewriting it with Eq. (4).\nIf we have found the true adjacency matrix $\\mathbf{A}$, the objective becomes  \n$$\\left\\|\\left\\|\\mathbf{X}\\left(\\mathbf{I} + \\overline{\\mathbf{A}}\\right)^{-1}\\right\\|\\right\\| =\\left\\|\\left\\|\\mathbf{C} + \\mathbf{N}_c + \\mathbf{N}_x\\left(\\mathbf{I} + \\overline{\\mathbf{A}}\\right)^{-1}\\right\\|\\right\\| = \\left\\|\\left\\|\\mathbf{C} + \\mathbf{N}_c + \\mathbf{N}_x\\left(\\mathbf{I} -\\mathbf{A}\\right)\\right\\|\\right\\|$$\nThis means that we approximate the root causes up to noise (also in next answer).\n\nWe will explain better.\n\n**Learning the root nodes.** Yes, we can do this in the following way. If we recover the true adjacency matrix $\\mathbf{A}$ via optimization Eq. (10),\nwe may compute an approximation $\\widehat{\\mathbf{C}}$ of the root causes $\\mathbf{C}$ up to noise, by solving Eq. (4): $$\\widehat{\\mathbf{C}}= \\mathbf{C} + \\mathbf{N}_c + \\mathbf{N}_x\\left(\\mathbf{I} -\\mathbf{A}\\right) = \\mathbf{X}\\left(\\mathbf{I} + \\overline{\\mathbf{A}}\\right)^{-1}.$$\nThis obviously requires a very good (weighted) estimate of the original adjacency matrix. We evaluate the top-performing methods on the recovery of the root causes (and the associated values) with an additional experiment in Fig. 3 of the attached pdf.\n\n**Real-world dataset performance.** In Fig. 4 of the attached pdf we include an experiment that illustrates the sparsity in the root causes for the particular dataset [Sachs et. al., 2005]. Please see our general reply for further explanation and also our success in a causal discovery competition with real-world data.\n"}]

Author rebuttal:
# General Comments
We thank all the reviewers for their kind reviews and their effort and interest in understanding and commenting on our work. We will incorporate these comments in an improved revision. Here, we would like to address two main points that arose from the comments.

**Motivation for few root causes**
We view the pollution example as a metaphor that can be moved to various DAG scenarios as roughly sketched next.

For example, as pointed out by reviewer LucQ, the data could measure (amount of) water flow in which case the root causes would capture cities with major input (e.g., through rainfalls).

In gene networks that measure gene expression, few root causes would mean that few genes are activated in a considered dataset.

In a citation network where one measures the impact of keywords/ideas, few root causes would correspond to the few origins of them.

[Peters et. al., 2017] mention in pp. 19-21 the *Principle of independent mechanisms* which in essence makes the assumption that any causal data generation relies on independent mechanisms that turn the input into output. This corresponds to our viewpoint on the linear SEM with the independent root causes as input.

As a thought: with the input/output view of linear SEMs in Eq. (3) one can question why all real world DAG data should be generated from i.i.d. noise $\mathbf{N}$ as input.

**Application in real-world data**
Our method doesn't achieve the best performance on the dataset from [Sachs et. al., 2005], but is reasonably competitive, which is all we wanted to show.

*More importantly*, after the NeurIPS submission our method was among the three winning entries in a causal discovery competition, run by a major pharma company, with the results presented in a (non-archival) report and at a major ML conference (no specifics due to double-blind review). One reviewer commented "... the few root causes assumption may have biological relevance worthy of further investigation." Our method was the only entry that offered some theoretical guarantees and performed well even with a linearity assumption.

Also, motivated by a comment of Reviewer Aevs we did an experiment to estimate the root causes of the dataset from [Sachs et. al., 2005]. In the pdf, Fig. 4 we show the result. The dataset comes with an unweighted adjacency matrix. We chose weights to enforce sparsity on the root causes. The sparsity roughly holds (Fig. 4, right). Interestingly (Fig. 4, left), the root causes tend to have fixed support, which could explain why we are unable to perfectly recover the DAG.


Meta-review: After discussion, there was a weak consensus to accept this paper. The authors provide a novel perspective on identifying DAGs via root causes, and while the results have some limitations (outlined in the reviews), proposing new directions for identifying causal graphs is an important direction for the community. The inclusion of additional experiments during the rebuttal was convincing for several reviewers. For this reason, I recommend accept.

Nonetheless, I hope the authors will take the author feedback into account carefully. In particular, the following:
- There is a lack of rigor in some places that can be tightened up
- Even though previous work may have shown that the PC algorithm has limitations, the authors should replicate this behavior on their own, and include the results in the appendix (I would consider including GES as well)
Decision: Accept (poster)