QWO: Speeding Up Permutation-Based Causal Discovery in LiGAMs
Authors: Mohammad Shahverdikondori
Ehsan Mokhtarian
Negar Kiyavash
Venue: NeurIPS 2024 poster (2024)
Keywords: causal discovery
permutation-based
linear gaussian acyclic model
DAG learning
Primary area: causal_inference
TL;DR: We propose a method that enhances the time complexity of permutation-based causal discovery approaches in linear gaussian acyclic models.

Abstract: Causal discovery is essential for understanding relationships among variables of interest in many scientific domains. In this paper, we focus on permutation-based methods for learning causal graphs in Linear Gaussian Acyclic Models (LiGAMs), where the permutation encodes a causal ordering of the variables. Existing methods in this setting are not scalable due to their high computational complexity. These methods are comprised of two main components: (i) constructing a specific DAG, $\mathcal{G}^\pi$, for a given permutation $\pi$, which represents the best structure that can be learned from the available data while adhering to $\pi$, and (ii) searching over the space of permutations (i.e., causal orders) to minimize the number of edges in $\mathcal{G}^\pi$. We introduce QWO, a novel approach that significantly enhances the efficiency of computing $\mathcal{G}^\pi$ for a given permutation $\pi$. QWO has a speed-up of $O(n^2)$ ($n$ is the number of variables) compared to the state-of-the-art BIC-based method, making it highly scalable. We show that our method is theoretically sound and can be integrated into existing search strategies such as GRASP and hill-climbing-based methods to improve their performance.

QWO: Speeding Up Permutation-Based Causal Discovery in LiGAMs
Mohammad Shahverdikondori
College of Management of Technology
EPFL, Lausanne, Switzerland
mohammad chabvveddikkoooorr1@@pl.....
Ehsan Mokhtarian
Negar Kiyavash
School of Computer and Communication Sciences
College of Management of Technology
EPFL, Lausamne, Switzerland
EPFL, Lausanne, Switzerland negar. www.kkaaaaaspfl...cmm
ehsan. Chsan.mokkttrriaa@@ff..cc
Abstract
Causal discovery is essential for understanding relationships among variables of interest many scientific domains. In this paper, we focus on permutation- bassed methods for learning causal graphs in Linear Gaussian Acyelic Models (LiGAMs), where the permutation encodes a causal ordering of the variables. Existing meth- ods in this setting do not scale due their high computational complexity. These methods are comprised of two main components: (i) constructing a specific DAG, G". for given permutation स,, which represents the best structure that can be Iearned from the available data while adhering to and (ii) searching over the space of permutations æs (i.e., causal orders) to minimize the number of edges GT. We introduce QW-Orthogonality (QWO), novel approach that significantly enhances the efficiency of computing Gd for given permutation л. QWO has a speed-up of Onn²) (n is the number of variables) compared to the state-of-the BIC-based method, making highly scalable, We show that our method is theoreti- cally sound and can be integrated into existing search strategies such as GRaSP and hill-climbing--aal methods to improve their performance. The implementation is publicly available at https: h //gittttbbbb comm homm/ann-eee ban-eeff//WW.
Introduction
Causal discovery is fundamental understanding and modeling the relationships between variables in various scientific domains [Pea09, SSSSH00.. A causal graph, typically represented by directed acyelic graph (DAG), is a graphical model that represents how variables within system influence one another [Pe00)] The problem of causal discovery refers to learning the causal graph from available data [SGSH00, MT99, Chi02, FZ13. ZARXI&, BSSU20a. MAGK21, LAR22 MAJ+22 22 MEAK24]. which has broad applications across numerous fields such as economics Hec08], genet- ISMD 033. and social sciences [MW15]]
Score-based methods are an important class of approaches for causal discovery [US10, KGG +2 HS13. HZL 18, Sch78, Bun91]. They involve defining score function over the space of DAGs and secking the graph that maximizes this score. However, the computational complexity of exploring the space of all possible DAGs. which is in the order of 2(("") where n is the number of variables, poses significant challenge [Rob73]. Ordering- based methods, introduced by [TKO5), significantly reduce
38th Conference on Neural Information Processing Systems (NeurIPS 2024).

---

Table 1: Time complexity comparison of various methods for computing and updating 9 in LiG- AMs. Here, n is the number of variables, d is the length of the updated block of the permutation, and is the number of folds considered in k-fold cross- noss-validation.
Method
QWO
BIC
BDeu
CV General
Initial Complexity
O(n")
O(n²)
Q(n N log(N))
O
Update Complexity
d(
(nn'd)
di log(N))
O RdN
the search space of score-hased methods to (0 by considering only topological orderings of DAGs, rather than the DAGs themselves [FK03, ZNC20, MKEK23, (SLOT23.. Building on this premise, various permutation-basee methods have been proposed that further refine the search for causal structures [SdCCzIS, LAR22, SWU21 These methods have two main ingredients:
1. Constructing G*: A module that for any given permutation x constructs specific DAG, denoted by G We will formally define GT in Definition 3.I. but roughly speaking, this DAG represents the best structure that can be inferred from the data while conforming to the given permutation и.
Search over T: The work by [RU18] (which proposed the Sparsest Permutation algorithm) showed that the correct permutation minimizes the number of edges in G" Accordingly, permutation-basae methods are equipped with a search strategy that, using the module mentioned above for computing GT, searches through the space of permutations minimize the number of edges in G".
In recent years, various search methods have been developed enhance the accuracy, robustness, and scalability of the algorithm for the second part, i.e., searching over the permutations [LAR22 TBA06, MKEK23 TK0SJ. These methods typically involve traversing over the space of permutations by iteratively updating a permutation to reduce the number of edges in G using the first module to construct Gd multiple times during the algorithm's 's execution. To construct G or update G" after permutation update, most methods utilize decomposable score function and identify the parent set that maximizes the score of each variable. Table L presents time complexity comparison of different methods for computing and updating G We will delve into more detailed discussion of these methods and their computational complexity in Section S.
In this paper, we focus on linear Gaussian acyclic models (LiGAMs), an important class of continuous causal models. We propose novel method called QW- Orthoooognnaltt (QWO) for computing G for a given permutation rT, which significantly enhances the computational efficiency of this module. Specifically, as shown Table I QWO's complexity is independent of the number of data points N. Moreover, it speeds up computing/updatinn G by O(n²) compared to the state-of- -o-thh--artttt BIC-based method. Some key advantages of QWO are as follows:
1. Soundness: QWO is theoretically sound for LiGAM models. That is, it is guaranteed to learn Gm accurately for any given permutation æ when sufficient samples are available.
Scalability: The proposed method is scalable to large graphs. Its time complexity is independent of the number of samples, and its dependence on the number of variables n. is Onn)) faster than the wtate-of-the-ad BIC-based method.
3. Compatibility with existing search methods: After updating a permutation, QWO efficiently updates GD to improve its compatibility with established search methods. In our experiments. we combine QWO with both GRaSP [LAR22] and hill-climbing-basse search technique [BAA0 SG061, showing its superior performance in terms of time complexity.
2 Notations
Throughout this paper, matrices are denoted by capital letters, and sets or vectors of random variables are denoted by bold capital letters. We use the terms "variable" and "vertex" interchangeably in graphs. [r]] denotes the set of natural numbers from to n, and [[(x)) denotes the set of all n!
permutations over [n). For permutation æ [n). Pr denotes the permutation matrix, such that a permutes the rows and columns of A according to the permutation и. The identity matrix is denoted by T, and its dimension is implied by the context.
for any n x n matrix A, the product P, AP

---

where N (0,0 Σ) is a Gaussian random vector with covariance matrix E. Equivalently, we have ||Al|o denotes the number of non-zero entries in matrix A. The set of n x n diagonal matrices is represented by D,, (a. bj denotes the inner product of vectors a and b. For a set X {X)..... Xx, and S _ [n], we define Xs- = (X:li e S). diag( A) denotes a diagonal matrix with the same diagonal elements as matrix A.
In directed graph (DG), edges are directed and may form cycles. A DG without cycles is called a directed acyclic graph (DAG). If there directed edge from X to Y. then X is parent of Y. X is an ancestor of Y if there is a directed path from X to Y. In a DAG G. for any vertex X, e(g). Pag (X), and (X) denote the number of edges, the parent set of X, and the ancestor set of X,
respectively. Suppose π is a permutation over the vertices of a DAG G = (X, E). IT is topological order for G, or equivalently G is compatible with म,, if for any edge Xx()) + XnO) E, j.
Definition 2.1 (LiGM, LiGAM, G(B)). Suppose X = =[X1,..., X,)² is random vector, B is an n X n. matrix such that 1 - B is invertible, and Σ E Dnn Pair M = (B,2) is called a linear Gaussian model (LiGM) that generates X if
X= (I-
(I)
X = BX+N.
(2)
The causal graph of M, denoted by G(B), is DG with the adjacency matrix corresponding the support of BT A LiGM is a linear Gaussian acyclic model (LiGAM) if its causal graph is a DAG.
The covariance matrix of X in a LiGM M = (B, Σ2 is given by
(3) For distinct indices < i S and S S X. we use X; TL X,|Xs to denote that X; and X; are independent conditioned on Xs. The notion d-separation defined over DAGs is a graphical criterion to encode conditional independence (CI) within graph. We similarly use X; JL X,5Xx, to denote that X; and X, are d-separated given Xs in a DAG. For a formal definition of d-separation, see [[ea88]].
Cov(X) BJ- - B)-
Definition 2.2 (19). Two DAGs are Markov equivalent if they impose the same d-separations. For a DAG G, the set of its Markov equivalent DAGs are represented by [g].
In LiGAM, the Markov property states that if two variables are d-separated in the DAG, then they are conditionally independent in the corresponding probability distribution. This property holds for structural equation models (SEMs), including LiGAMs; see Theorem 1.2.5 in [Pea09]] for more details. Conversely, the faithfulneess assumption posits that if two variables are conditionally inde- pendent in the distribution, then they are separated in the DAG. Together, the Markov property and faithfulness establish one-to-one correspondence between the graphical d-separations and the conditional independencies in the distribution.
While the faithfulness assumption provides a strong correspondence between the distribution and the graph, it can be restrictive in practical applications. Recognizing this limitation, weaker versions of faithfulness have been proposed. One such alternative is the sparsest Markov representation (SMR) assumption introduced by [RU18]. Formally, for DAG G and a distribution P defined over the vertices of G. (G, P)) satisfies the SMR assumption if (G, P) satisfies the Markov property, and for every DAG g' such that (G', PP) also satisfies the Markov property and G' & [9]. it holds that [g'| > |G]. Here, |g] denotes the number of edges in G.
3 Permutation-Based Causal Discovery in LiGAMS
In this section, we discuss permutation-basee methods for causal discovery in LiGAMs. We begin by formalizing our assumptions, which will be referenced throughout the remainder of the paper.
(B", Σ*) be LiGAM that generates the random vector X (X1,X2.- .2.x2.,........ ..). Let D be N data mairix, where euch row is an i.id. sumple from X. We assume that the pair (g", Pr) satisfies the (SMR) assumption, where G" = =GBB) is the causal graph of M and pa is the joint distribution of X.
Assumption 1. Let M
3

---

Under Assumption | our goal is to learn the causal graph G" from the observational data D. How- ever, using mere observational data-even when the number of samples N approaçhes infinity- G can only be identified up to its Markov equivalence class (MEC) [SGSH00, Pea09]. Therefore, the best we can aim for in the causal discovery of LiGAMs to identify (9") a problem known to be computationally NP-hard [CHM04].
As briefly discussed in the introduction, score-based methods aim minimize score function in the space of DAGs, whereas permutation-based methods restrict this space to topological orderings of DAGs. Thus, they seek a permutation и such that ı is topological order of DAG in [g*].
Definition 3.1 (G*), Under Assumption I for any arbitrary permutation E [1([n]), we denote by 9" = (X, E²) the unique DAG constructed from P as follows: For distinct indices S i, S n,
there is an edge from Xx() to X xx(s) in E if and only if
i i<j and Xx(3) V.. (1)
Remark 3.1. For any n, 9 is compatible with x. Furthermore, if DAG G E [g*) is compatible with nI, then G" G.
[RU18] showed that a permutation π minimizes the number of edges in DAG G" if and only if is the topological order of DAG in (S*). Therefore, permutation- based methods typically formulate causal discovery as follows.
arg min E
TEU(N)) (x))
Various algorithms have been proposed for solving s LAR22 SWU21
TK05,
FK03
As discussed
in the introduction, these methods include two components:
1. Constructing G A module that for given permutation constructs G Search over x: search strategy over the space of permutations to solve
Given that minimization in b), search strategies typically traverse the space of permutations and greedily update the permutation. To maintain computational efficiency, these methods use module to update 9d incrementally rather than recomputing it from scratch. Moreover, good search strategies in permutation- based causal discovery. A consistent search method guarantees to solve b as long as it is equipped with module that correctly computes G In Appendix A methods [TBA06, SGOGJ. GRaSP is consistent, ensuring that it finds the correct permutation that minimizes the number of
are ddeaalyycconistte
we review two search strategies: GRaSP [LAR22] and hill-climmingg-aare
edges in GT. contrast, hill-climbing methods do not provide consistency guarantees but are shown to be efficient in practice.
Existing methods for computing Gd are primarily score-based and involve the following optimization:
arg max DAG g
D,
s.t. G compatible with n.
where typically decomposable score function, summing individual scores for each variable given its parents within the graph:
D,=) = Σ Pux))
1=1
To find the parent set of each variable in GT, the following optimization is performed for all 1 Sii n:
Page (Xm(i))
arg max
S(X=(), D,
(8)
UC(N=()....X
To solve various approaches have been proposed. Note that the complexity of brute-force search over all subsets U exponential, which makes it impractical. Instead, the state-of-the-at search methods apply the grow-shrink (GS) algorithm [ARSR + 23] on the candidate sets U to find the parent set of each variable. These methods require computing the score function S, O(x²) times.
It is noteworthy that the chosen score function must ensure the solution to this optimization problem g". Examples of such scores include BIC [Sch78], BDeu [Bun91], and CV Gen- eral HZL 18). Table compares the time complexity of our proposed method (QWO) against

---

these methods. The Bayesian Information Criterion (BIC) is a well-known score in the literature. particularly for LiGAMs. BIC's complexity for calculating the initial graph is O(n²); the update is Onn²d). The Bayesian Dirichlet equivalent uniform (BDeu) score applies uniform prior over the set of Bayesian networks and uses this prior to evaluating the model's accuracy. While BDeu was originally designed for discrete data, it has been adapted for continuous data by dividing the real numbers into intervals and assigning a constant value for each interval. The time complexity of BDeu depends on the number of distinct values each variable can take (i.e., intervals). Its initial and update complexities are lower-bounded by $2(n" N (V)g(N) and (2(nddNlog(N))
respectively. The Gener- alized Cross-Vate Validated Likelihood (CV General) score involves splitting the dataset into training and test sets multiple times. The final score for a variable given its parents is the average log- ligeleidlooo evaluated on the test sets using the regression functions learned from the training data. The CV General method has an initial complexity of O O(" (²² and an update complexity of O( O(Td>") k² ). This method typically uses small constant k for k-fold cross-validation, which significantly increases the computational complexity. Among the aforementioned three strategies, BIC is the fastest. Our proposed method, QWO, attains speed-up of On²²) compared to BIC.
k22
4 QWO
In this section, we present a novel approach for computing 9d for a permutation т,, with improved computational complexity, which can be easily integrated into existing search methods. Our proposed method proceeds with an alternative formulation for causal discovery in LiGAMs, but first, we need a definition.
Definition 4.1 (B()). For a random vector X, we denote by B(X) the set of coefficient matrices of LiGMs that generale X, i.e.
( {B E Dz: Cov(X) = - () 3)-
Note that this definition, the causal graphs corresponding to the elements of B(X) can be eyclic, but we restrict our attention to a subset of (() with acyelic corresponding graphs. We reformulate causal discovery in LiGAMs as follows:
arg min BEB(X)
||B||o-
s.t. G(B) is a DAG
has been shown that for any solution of @.. G(B) belongs to [g+). Furthermore, for any graph
G E [g*]. there exists a solution B to (9) such that G(B) G [Pea09]. Therefore, solving (9) is
equivalent to performing causal discovery in LiGAMs.
In the following, we establish the relationship between B(X) and GT.
Theorem 4.2. Under Assumption 1, for any permutation π E [I([n]), there exists uniquue E B(X) such that G(B) is compatible with π. Furthermore, for this B, G(B) GT.
All proofs are provided in Appendix C Theorem 4.2 implics that compute GT, can directly find the unique B E (X) whose corresponding graph G is compatible with x. Next, we propose a characterization for B(X) using whitening transformation [Fuk90, HO0d, KLS18].
Definition 4.3 (Whitening matrix W). Let Cov(X) = USUT be the singular value decomposition (SVD) of Cou(X), where S is a diagonal matrix including the singular values of Cov(X) its diagonal and U is an orthonormal matrix (i.e., UU = 1). The whitening matrix W is defined as
W = US İUT
(10)
We note that the 'W' in QWO corresponds to the whitening matrix W. Whitening is a linear transfor- mation N, == W/X that transforms the Gaussian random vector X to another Gaussian random vector N,, where (N) = T. Furthermore, for an arbitrary orthogonal matrix Q (i.e., QQ" E D,,), we define No = QN, QWX, then it straightforward to show that Nę is also a Gaussian random vector with mean zero and a diagonal covariance matrix. Therefore, (I QW, No) is a LiGM that generates X since
x= -I QW)X + No-
In the following, we show that such LiGMs create all possible LiGMs that generate X.
5

---

Algorithm The QWO module for computing and updating 9
1: Function QT) W. [Optional] idx,, [Optional] idx,, [Optional] Q)
2: Default values for optional arguments: idx; 1, idx, n, Q = any arbitrary n x n matrix
3: Denote the columns of Q" by {q: Pa=1 and the columns of W by {ws}i=1
4: for from idx, idx; do
5: r3 + [x(( 2i-i+1 (we())-0))))
1 9m())
6: (m(1) (1i:)*())
7: end for
8: G G(1 QWW
9: Return 0",Q
Theorem 4.4 (Characterizing B(X)). For any B E B(X), there exists unique orthogonal matrix Q such that B I QW, and vice versa. That is,
B(X) -I- QWI QQ DD
(11)
By combining Theorems 4.2 and 4.4. to learn G", sufficient to identify the unique orthogonal matrix Q such that G(I QW) is compatible with æ. To verify this compatibility, we impose the following two constraints:
P&QW PF is upper triangular,
di =I.
(12)
4.1 The QWO Method
In this subsection, we introduce OW-Orthogonality (WW),, our proposed approach for computing G in LiGAMs. The function QWO in Algorithm | presents the pseudocode of our method. The function takes a permutation π and a whitening matrix W as inputs, with optional arguments idxl, idxr, and Q. First, consider the case where these optional arguments are initialized as follows: idx; = 1. idx, = n, and Q to be an arbitrary n x n matrix. The goal of this function is create a matrix Q that is the unique solution of (12) and subsequently compute g"
not provided, and they
Denote the columns of QT by {a 1 and the columns of W by [w;} Since W = AT. all the eigenvalues dimensional space. Furthermore, because Q is an orthogonal matrix for each pair (i,j), S i + j S n, 4s 1 QT. The condition P2QWP* is upper triangular in (12) implies that (a(() should be orthogonal (x()) if i. Therefore, we need find (q: such that this condition holds. Note that each zero in | QW corresponds an orthogonality between { and {wi)
are positive, and thus {w:);=1 are n linearly independent vectors in an n-
With this intuition, we propose our approach for constructing Q. which maximizes the number of aforementioned orthogonality. To do so, we apply the
algorithm| to the vectors ) in the following order: (n), (n - 1). .... (1). In Algorithm ! we iteratively for each i, compute Tú, the residual of projecting (x()) the span of (Wx(I-1))--- (n)n)) This is equivalent to projecting (x(i) on the span of {9=(++1)..... (x(n) J. which are orthogonal to each other. (x(s) is set to normalized residual ri, which ensures that the product Irin and Wz(i) is 1, and consequently,
diag(A P= RW p = diag(QW)= = I. Finally, we use Theorem to construct G" = G(T QW) and return 9D and Q.
Theorem 4.5 (Soundness Algorithm [). Under Assumption and given the correct whitening matrix W as input, matrix Q, the output of Algorithum 1 is the unique solution (12). Consequently.
the returned graph corresponds to the true G defined in Desinition 3.1
The optional arguments idx;, idx,, and Q allow for incremental efficient updates to π.
Lemma 4.6. If the block between the idx th and -th positions of π is modified, the vectors (=(k) for k < idx or k >> idx, remain unchanged.
A consequence of Lemma 4.6 is that (k) for k < idx; or k > idx, from the previous computed Q can be reused and sufficed to merely recompute the vectors within the updated block by iterating through the for loop in lines 4-7.

---

Algorithm 2 Integrating QWO into a simple search method for causal discovery
1: Input: Data matrix D, Search method F.
2: Estimate Cov(X) using data matrix D
3: U,S t Compute the SVD decomposition of Cov(X), such that Cov(X) USU
4: W <- -S-
% Whitening matrix
5: π +- An initial permutation
6: G Q =QWO(T,, W)
7: while has not stopped do
8: n² +Update according to F, and let idx; and idx, be the leftmost and rightmost index of
T that have been updated, respectively
9:
G², "=''Q'- D'<=WO(m: W, idx,, idx,. Q)
10:
if 9T has less number of edges than G then
11: T, T.Q.9" c 0'
end if
13: end while
Return G
Theorem 4.7 (Time complexity of Algorithm D QWO algorithm implemented in Algorithum has the following time complexities:
( for initially computing G without optional arguments.
(n²² when called with optional arguments to update G*, where d = idx, idx.
4.2 Integrating QWO in Existing Search Methods
Algorithm E demonstrates how QWO can be integrated into existing search methods such as GRaSP [LAR2]] and Hill-Climbing (TBA06, SG06). This algorithm integrates QWO into sim- ple permutation- based search method for causal discovery that iteratively updates the permutation to minimize the number of edges in g" Initially, the algorithm estimates the covariance matrix of and applies SVD to compute the whitening matrix. Starting from an initial permutation, calls function QWO compute Q and g Subsequently, the algorithm iteratively updates the permuta- tion using the given search method A For each updated permutation, it calls function QWO with optional arguments compute the Q and g graph. Next, the algorithm checks if the per- mutation is better than the previous one by checking whether the new graph has fewer edges. the new permutation is better, it updates the permutation and proceeds the next iteration. The search algorithm stops when its stopping criterion is satisfied and returns the final æ and GT.
5 Experiments
In this section, present a comprehensive evaluation of QWO, which is designed for the module that computes G" in permutation-ban based causal discovery methods in LiGAMs. As discussed earlier, a permutation- method consists of two components: search method and module for computing g². Our comparison focuses on the module for computing G*, where we benchmark QWO against existing methods, namely BIC [Sch78], BDeu [Bun91]] and CV General [HZL+ 18]. For the search method, we utilized two specific algorithms: GRaSP [LAR22] and Hill-Climbing (HC) [TBA06 SG06J. Please refer Appendix B for the implementation details of these methods and additional results.
We generated random graphs according to Erdos-Renyi model with average degree of for each node, denoted by ERi. We did not impose any constraints on the maximum degree of the nodes. To generate the data matrix D using LiGAM (B", Σ²), sampled the entries of B from [-2, 0.55 U [0.5, 2] and the noise variances uniformly from [1,2]. Each reported number on the plots an average of 30 random graphs.
SSeacch method F could be any existing approach such as GRaSP LAR221 or Hill- Hill-Climbing TBA06 SG06].

---

Table Results for small graphs using GRaSP and Hill- Hill-Climbing search methods.
Graph
CANCER
SURVEY
ASIA
SACHS
ER2
Number of Nodes
5
8
11
5
Average Degree
1.6
2
2
3.09
2
SKFI
1
0.87
0.88
0.85
GRaSP
PSHD
0
0.87
0.63
0.2
Time (s)
0.003
0.004
0.01
0.04
0.002
QWO
SKFI
0.88
0.92
0.94
0.88
0.97
HC
PSHD
0.6
0.5
0.625
0.63
0.4
Time (s)
0.01
0.01
0.01
0.05
0.003
SKFI
1
1
0.93
GRaSP
PSHD
0
0
0.81
0
Time (s)
0.06
0.07
0.08
0.39
0.02
BIC
SKFI
0.88
0.61
0.94
0.93
0.97
HC
PSHD
0.8
1,33
0.625
0.54
0.4
Time (s)
0.01
0.05
0.08
0.83
0.15
SKFI
0.25
0.36
0.26
0.29
0.12
GRaSP
PSHD
1.4
1.5
1.5
1.72
1.4
Time (s)
44.1
72.1
116.6
211.6
41.1
BDeu
SKF1
0.5
0.54
0.4
0.37
0.6
HC
PSHD
1.2
1.16
1.374
1.81
1.2
Time (s)
64.4
95.1
224.4
496.6
78.2
SKFI
1
0.57
0.76
0.87
0.85
GRaSP
PSHID
0
1,66
1.12
0.90
0.2
Time (s)
109.6
296.3
641.1
867.5
101.6
CV General
SKFI
1
0.5
0.88
0.87
0.9
HC
PSHD
0
1.33
1.12
0.72
0.6
Time (s)
176.7
358.6
1169.4
1325.7
213.3
Two metrics were used to evaluate the performance the methods:
Skeleton F1 Score (SKF1): The F1 score between the skeleton of the predicted graph and the skeleton of the true graph, which ranges from to where indicates perfect accuracy.
Complete PDAG SHD (PSHD): For the true DAG G and predicted DAG G, PSHD calculates the complete PDAGs of [g] and [9] and evaluates the number of changes (edge addition, removal, and type change) needed to obtain one PDAG from the other, normalized by the number nodes. This metric shows the distance between the Markov equivalence classes of the predicted graph and the true graph.
We carried out experiments in four settings:
Low-Dimensional Data: We evaluated the performance of QWO and other methoods small real-world graphs, namely ASIA [LS88], CANCER [KN10]. SACHS [SPP-05, 05]. and SURVEY [SD21]]. as well as ER2 graphs with 5 nodes. The number of data samples for this part is set to 500. The results different methods using both search strategies are presented in Table As shown in the table, QWO demonstrates comparable accuracy the other approaches for different graph structures while being significantly faster.
High-Dimensions Data: In this setting, we assessed the performance of QWO on larger graphs with 10,000 data samples. Due to the exceedingly long runtime of BDeu and CV General, which are orders of magnitude slower this setting, they were not included this experiment. Instead, we compared QWO with the BIC score using both GRaSP and Hill- Hill-Climming search methods. As illustrated in Figure i while maintaining high accuracy, QWO demonstrates significant speedup over BIC for both search methods.

---

1.0
2.5
100
80
SKF]
0.6
85HD
Time(s)
1.0
0.4
HC-QWO
grasp grasp-QWO
20
grasp-BIC
0.0
150
Number of variables
Number of variables
Number variables
Results on R2 graphs.
1.0
100-
0.8
80
SKF1
PSHD
1.5
Time(s)
60
0.4
1.0
40
HC-8IC
gresp-OIC
0.0
01
150
50
100
150
Number of variables
Number of variables
Number of variables
Results on ER3 graphs.
3.0
120
2.5
0.8
SKF1
PSHD
1.5
Timo(s)
60
1.0
0.4
40
HC-QWO
0.2
HC-BIC
0.0
grasp-QUP 2WP
grasp BIC
0.5
50
100
150
100
150
0
100
150
Number of variables
Number of variables
Number of variables
Results ER4 graphs.
Figure 1: Results of performing QWO and BIC on both search methods GRaSP and HC on data generated by LiGAMs for different Erdos-Renyi graphs (ER2, ER3, ER4).
Non-Gaussian Noise Experiments: To test the robustness of method, conducted experi- ments where the main assumption of Gaussian noise in the data-generating process was violated. We evaluated QWO on linear models with non-Gaussian noise distributions (specifically expo- nential and Gumbel distributions). The results of these experiments appear in Appendix B Al- though QWO is designed for linear models with Gaussian noise, these experiments show that QWO achieves almost similar accuracy to LiGAMs on models with exponential and Gumbel noise dis- tributions.
Oracle Inverse Covariance Experiments: In practice, errors in learning the graph may arise from two sources: the error in calculating the inverse of the covariance matrix and the error in the optimization problem. To study the effect of each error separately, we designed an experiment to eliminate the first source of error by providing the correct inverse covariance matrix (similar the approach [DUF+23). We then compared the performance of the algorithms under this condition. The results of these experiments appear in Appendix B The plots show that the accuracy

---

of our method is now significantly high, indicating that the main source of error lies in the initial step of estimating the covariance matrix.
6 Conclusion, Limitations, and Future Work
In this paper, we proposed QWO, a module to accelerate permutation- based causal discovery in LiGAMs. Our method reduces the computational complexity of constructing and updating the graph 9 for a given permutation by O(n²), resulting in a significant speed-up compared to the state-of- the-art BIC-based method, as demonstrated both theoretically and through extensive experiments. Furthermore, QWO seamlessly integrates into existing search strategies, enhancing their scalability without compromising accuracy.
While our method offers substantial improvements, it is important to acknowledge its limitations. The primary assumptions underpinning QWO are that the underlying causal model is LiGAM- that is, it assumes linear relationships among variables, additive Gaussian noise, and an acyclic causal graph structure. These assumptions, though common in many applications, restrict the applicability of our method to scenarios where these conditions hold. In real-world datasets, relationships may be nonlinear, noise may not be Gaussian, the causal graph may contain cycles due to feedback loops or reciprocal relationships among variables. Notably, our experiments indicate that QWO maintains competitive performance even when the Gaussian noise assumption is violated, achieving similar accuracy on models with non-Gaussian noise distributions such as exponential and Gumbel (see Appendix B). However, the acyclicity assumption can be particularly limiting in domains where feedback mechanisms are inherent, such as economics, biology, and control systems, where causal graphs are cyclic and methods designed for acyclic graphs are not directly applicable.
Despite this limitation, our work opens avenues for extending permutation based causal discovery to more general settings. Notably, the set B(X) is defined over all LiGMs, encompassing both eyclic and acyclic models. Therefore, our characterization of ((X) in Theorem 4.4 holds for LiGMs as well. Combining this result with our reformulation of causal discovery in LiGAMs in (). we can generalize the formulation to LiGMs as follows:
arg min Q:QQT ED.
-
(13)
This suggests that causal discovery in LiGMs can be approached by finding an orthogonal matrix Q that sparsifies I - QW. Some recent works have exploredd causal discovery in cyclic models using orthogonal transformations similar to ours [GYKZ20]. However, solving the optimization problem in Equation (13) is challenging due to its non-convexity and the orthogonality constraint on Q. Developing efficient algorithms to solve this problem is a promising direction for future work.
Acknowledgments
This research was in part supported by the Swiss National Science Foundation under NCCR Automa- tion, grant agreement HINF40_180544 and Swiss SNF project (000221204335/111 /1.
References
[ARSR+23) Bryan Andrews, Joseph Ramsey, Ruben Sanchez Romero, Jazmin Camchong, and Erich Kummerfeld. Fast scalable and accurate discovery of dags using the best or- der score search and grow shrink trees. Advances in Neural Information Processing Systems, 36, 2023.
[Ber09] Dennis S Bernstein. Matrix mathematics: theory, facts, and formulas. Princeton uni- versity press, 2009.
[BSSU20a] Daniel Bernstein, Basil Saeed, Chandler Squires, and Caroline Uhler. Ordering-based causal structure learning in the presence of latent variables. In International Conference on Artificial Intelligence and Statistics, pages 4098-4108. PMLR, 2020.
[BSSU20b] Daniel Bernstein, Basil Saeed, Chandler Squires, and Caroline Uhler. Ordering-based causal structure learning in the presence of latent variables. In International conference artificial intelligence and statistics, pages 4098-41000 PMLR, 2020.
10

---

[Bun91] Wray Buntine. Theory refinement on bayesian networks. In Uncertainty proceedings 1991, pages 52-60. Elsevier. 1991.
[Chi02] David Maxwell Chickering. Optimal structure identification with greedy search. Jour- nal of machine learning research, 3NNov):507-5555 2002.
[CHM04] Max Chickering, David Heckerman, and Chris Meek. Large-sample learning of bayesian networks np-hard. Journal of Machine Learning Research, 5:1287-1330, 2004.
[DUF+23] Shuyu Dong, Kento Uemura, Akito Fujüi, Shuang Chang, Yusuke Koyanagi. Koji Maruhashi, and Michèle Sebag. Learning large causal structures from inverse covari- ance matrix via matrix decomposition. arXiv preprint arXiv:221 1.12211, 2023.
[FK03] Nir Friedman and Daphne Koller. Being bayesian about network structure. a bayesian approach to structure discovery in bayesian networks. Machine learning, 50:95-125, 2003.
[Fuk90] Keinosuke Fukunaga. Introduction to statistical pattern recognition. Elsevier, 1990.
[FZ13] Fei Fu and Qing Zhou. Learning sparse causal Gaussian networks with experimental intervention: regularization and coordinate descent. Journal of the American Statistical Association, 108(501):288 -000000, 2013.
[GL13] Gene H. Golub and Charles F. Van Loan. Matrix Computations. Johns Hopkins Uni- versity Press, 2013.
[GYKZ20] AmirEmad Ghassami, Alan Yang, Negar Kiyavash, and Kun Zhang. Characterizing distribution equivalence and structure learning for cyclic and acyclic directed graphs. International Conference on Machine Learning, pages 3494-3504. PMLR, 2020.
[Hec08] James J Heckman. Econometric causality. International Statistical Review, 76(1):1-21 2008.
[HO00] Aapo Hyvärinen and Erkki Oja. Independent component analysis: algorithms and applications. Neural networks, (3(4-5):411-4 1-300, 2000.
[HS13] Aapo Hyvärinen and Stephen M Smith. Pairwise likelihood ratios for estimation of non-gaussian structural equation models. The Journal of Machine Learning Research, 14()):111-12, 2013.
[HZL*18] Biwei Huang, Kun Zhang, Yizhu Lin, Bernhard Schölkopf, and Clark Glymour. Gener- alized score functions for causal discovery. In Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining, pages 1551-1560, 2018.
[JS10] Dominik Janzing and Bernhard Schölkopf. Causal inference using the algorithmic markov condition. IEEE Trunsactions on Information Theory, 56(10):5168-5) 2010.
[KF09] Daphne Koller and Nir Friedman. Probabilistic graphical models: principles and tech- niques. MIT press, 2009.
[KGG+22] Diviyan Kalainathan, Olivier Goudet, Isabelle Guyon, David Lopez-Paz, and Michèle Sebag. Structural agnostic modeling: Adversarial learning of causal graphs. Journal of Machine Learning Research, 23(219):1-62, 2022.
[KLS18] Agnan Kessy, Alex Lewin, and Korbinian Strimmer. Optimal whitening and decorre- lation. The American Statistician, 72(4):309-314, 2018
[KN10] Kevin B Korb and Ann E Nicholson. Bayesian artificial intelligence. CRC press, 2010.
[LAR22] Wai-Yin Lam, Bryan Andrews, and Joseph Ramsey. Greedy relaxations of the spars- est permutation algorithm. In Uncertainty Artificial Intelligence, pages 1052-1062. PMLR, 2022.

---

[LS88] SL Lanritzen and David J Spiegelhalter. Local computations with probabilities on graphical structures and their applications expert systems (with discussion). JR Star. Soe.(Se 1988.
[MAGK21] Ehsan Mokhtarian, Sina Akbari, AmirEmad Ghassami, and Negar Kiyavash. A recur- sive Markov boundary-basee approach to causal structure learning. In The KDD'21 Workshop on Causal Discovery, pages 26-54. 54. PMLR, 2021.
[MA]+22] Ehsan Mokhtarian, Sina Akbari, Fateme Jamshidi, Jalal Etesami, and Negar Kiyavash. Learning Bayesian networks in the presence of structural side information. Proceeding of AAAI-22, the Thirty-Sixthl Sixth AAAI Conference on Artificial Intelligence, 2022.
[MEAK24] Ehsan Mokhtarian, Sepehr Elahi, Sina Akbari, and Negar Kiyavash. Recursive causal discovery. arXiv preprint urXiv:2403.09300 2024.
[MKEK23] Ehsan Mokhtarian, Mohmmadsadegh Khorasani, Jalal Etesami, and Negar Kiyavash. Novel ordering-based approaches for causal structure learning in the presence of un- observed variables. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, pages 12260-12268, 2023.
[MT99] Dimitris Margaritis and Sebastian Thrun. Bayesian network induction via local neigh- borhoods. Advances in Neural Information Processing Systems, 12:505-511, 1999.
[MW15] Stephen I. Morgan and Christopher Winship. Counterfactuals and Causal Inference: Methods and Principles for Social Research. Cambridge University Press, 2015.
[Pea88] Judea Pearl. Probabilistic reasoning in intelligent systems: Networks of plausible in- ference. Morgan Kaufmann, 1988.
[Pea09] Judea Pearl. Causality. Cambridge university press, 2009.
[Rob73] Robert W Robinson. Counting labeled acyclic digraphs. New directions in the theory of graphs, pages 239-273, 1973.
[RU18] Garvesh Raskutti and Caroline Uhler. Learning directed acyclic graph models based on sparsest permutations. Stat, 7(1), 2018
[Sch78] Gideon Schwarz. Estimating the dimension of a model. The annals of statistics, pages 461-464, 1978.
[SD21] Marco Scutari and Jean- Man-BBatiste Denis. Bayesian networks: with examples in CRC press, 2021.
[SdCCZ15] Mauro Scanagatta, Cassio P de Campos, Giorgio Corani, and Marco Zaffalon. Learn- ing bayesian networks with thousands of variables. Advances neural information processing systems, 28, 2015.
[SG06] Bart Selman and Carla P Gomes. Hill-climbing search. Encyclopedia of cognitive science, 81:82, 2006.
[SGSH00] Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
[SLOT23] Pedro Sanchez, Xiao Liu, Alison Q O'Neil, and Sotirios Tsaftaris. Diffusion models for causal discovery via topological ordering. In The Eleventh International Conference on Learning Representations, 2023.
[SMD+03] Eric Schadt, Stephanie A Monks, Thomas A Drake, Aldons J Lusis, Nelson Che, Victor Colinayo, ... and Stephen H Friend. Genetics of gene expression surveyed in maize, mouse and man. Nature, 422(6929):297-302, 2003.
[SP Karen Sachs, Omar Perez, Dana Pe'er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Sci- ence, 308(522)):5222222 2005.
12

---

[SWU21] Liam Solus, Yuhao Wang. and Caroline Uhler. Consistency guarantees for greedy permutation-base causal inference algorithms. Biometrika, 108(4):795-814 2021.
[TBA06] Ioannis Tsamardinos, Laura Brown, and Constantin Aliferis. The max-min hill- climbing bayesian network structure learning algorithm. Machine learning, 65:31-78, 2006.
[TK05] Marc Teyssier and Daphne Koller. Ordering-basee search: a simple and effective algo- rithm for learning bayesian networks. In Proceedings of the Twenty-First Conference on Uncertainty in Artificial Intelligence, pages 584-590, 2005.
[ZARXI8] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with no tears: Continuous optimization for structure learning. Advances in neural information processing systems, 31, 2018.
[ZHC+24] Yujia Zheng. Biwei Huang, Wei Chen, Joseph Ramsey, Mingming Gong, Ruichu Cai, Shohei Shimizu, Peter Spirtes, and Kun Zhang. Causal-learn: Causal discovery in python. Journal of Machine Learning Research, 25(60))1-8, 2024.
[ZNC20] Shengyu Zhu, Ignavier Ng, and Zhitang Chen. Causal discovery with reinforcement learning. In International Conference on Learning Representations, 2020.
13

---

Appendix
The appendix is organized as follows:
Appendix covery.
we review some of the existing search methods for permutation- causal dis-
In Appendix
we provide additional experimental results and discuss implementation details.
In Appendix we present formal proofs for the claims made in the main text.
A
Existing Search Methods
There is a vast literature on permutation-based methods that propose, different search strategies in the space of permutations [LAR22, TBA06, SWU21, BSSU20b, TK05]. In this section, we describe two widely used methods: GRaSP and Hiil-Climbinggberd based method.
A.1 GRaSP
This method was first introduced in [LAR22]. To present GRaSP, need two definitions:
I1)) and any i. [m], where precedes in т. Write v as 1 8, 83), where each 6, is a subsequence of m. Let n and be the subsequences
Definition A.1 (Tuck). Consider any permutation
(I E de Anco (i)) and (I $ 82 Anco(j)), respectively. Then define:
(61,7,5 2,,,,, 85) ifi Ancgn (j) otherwise
luck(a, (.
Definition A.2 (Covered Edge). In a DAG G, a directed edge i + between two nodes is called covered if Pag (i) (i)
[AR22] proved that starting from any arbitrary permutation there always exists sequence of permutations (H, mT, T2, , Ta) such that (i) Gmm E [9]. (ii) for each i, we can reach T from t, by tuck operation, and (iii) E(g*++1) C E(G (see Theorem 4.5 in [[AR22)). Subsequently, they propose greedy algorithm by applying the DFS algorithm on the following graph: A graph with n!! vertices, one for each permutation in II(n), and draw directed edge from node to л if is possible to find a covered edge + in G and ma obtained by performing a tuck on the pair (i,j) in TT. Finally, they show that under the faithfulness assumption, GRaSP is sound.
A.2 Hill-Climbing
There are various search methods based on the Hill-Climbing idea, which involve searching over the space of permutations through local changes [TBA06, SG06). A common way of updating the permutation among these methods is as follows: At each step on a permutation æ, among all pairs (i,5) with j j- & (where k is constant parameter), then calculated for the new permutation Feow obtained by swapping T (t)) and (j). the number of edges in new is less than g the process continues from nnew- The process stops when better permutation is found.
pair is chosen randomly. The score is
B Additional Experiments
In this section, we provide implementation details and additional experimental results for the oracle inverse covariance and non-Gaussian noise settings.
B.1 Implementation Details
We used the implementations provided the causal-learn library [ZC+24] for JIC.BBDeu. and CV General methods. 'here are three different versions of GRaSP introduced in [[AR22]] For our experiments, we used GRaSP, The depth of the DFS algorithm for GRaSP was set to and the value of (the maximum distance of swapped indices) in the HC algorithm was set to 5.
14

---

1.0
1.0
100
0.8
80
SKF1
0.6
R5HD
0.6
Time(s))
0.4
0.4
HC-QWO
0.2
HC-Oracle
grasp grasp-QWO
20
grasp-Dracle
0.0
0.0
100 Number variables
150
Number of variables
Number of variables
Results on R2 graphs.
1.0
100
0.8
80-
SKF1
6SHD
Time(s)
60
0.4
HC-Dracle
QWO
grasp-Oracle
0.00
0.0
0
150
100
150
100
150
Number of variables
Number of variables
Number of variables
Results on E ER graphs.
1.0
0.8
80
SKF1
PSHD
Time(s)
60
0.4
40
HC-QWO
0.2
HC-Dracle
paap-QWO
20
Oraple
0
50
100
150
50
100
150
50
100
150
Number of variables
Number of variables
Number variables
Results ER4 graphs.
Figure 2: Results for our original method compared to having the true inverse covariance matrix.
For the initial permutation of the search methods, considered the initial permutation based the size of Markov boundaries of the variables, which the case of linear Gaussian models are equal the number of non-zero elements in each of the inverse covariance matrix [KF09].
To generate the data matrix D using linear model (B*, Σ*), we sampled the entries of uni- formly from [-2, -0.5] U [0.5, 2] and the noise variances uniformly from [1,2] for all of the noise distributions i.c., Gaussian, Exponential, and Gumbel.
B.2 Oracle Inverse Covariance
Figure 2 shows the comparison between our original method, which calculates the inverse covariance from data, and using the oracle inverse covariance matrix. The plots demonstrate that the accuracy of our method is significantly high when using the oracle inverse covariance, indicating that the primary source of error resides in the initial step of estimating the covariance matrix.
"For the definition of Markov boundary, PPa099].
15

---

1.0
100
1.3
80
SKF]
0.6
PSHD
Time(s)
1.0
0.4
HC-QWO
0.5
grasp grasp-QWO
20
grasp-BIC
0.0
150
Number of variables
Number of variables
Number variables
Results on R2 graphs.
1.0
3.0
100
0.8
80
Time(s)
SKF1
1.5
60
0.4
1.0
40
0.5
HC-8IC
gresp-OIC
0.0
150
50
100
150
Number of variables
Number of variables
Number of variables
Results on ER3 graphs.
120
0.8
SKF1
PSHD
Timo(s)
60
1.0
0.4
HC-QWO
0.5
40
0.2
HC-BIC
grasp-QUP PWo
0.0
grasp BIC
0.5
50
100
150
100
150
0
100
150
Number of variables
Number of variables
Number of variables
Results on ER4 graphs.
Figure 3: Results for our method compared to BIC on linear models with exponential noise.
B.3 Non-Gaussian Noise
Herein, we present the results of rumning QWO and the BIC on two linear non-Gaussian models: Exponential and Gumbel noise. Figures| and 4 show the time complexity and accuracy in terms of two metrics for both methods both search strategies. Although QWO is designed for linear models with Gaussian noise, these experiments show that QWO achieves almost similar accuracy to LiGAMs models with exponential and Gumbel noise distributions.
16

---

2.5
120
2.0
0.8
1.5
0.6
Timels)
SKF1
1.0
60
HC-QWO
0.5
40
0.2
MC-BIC
grasp-OWO
0.0
20
grasp-BI
50
100
150
50
100
150
50
100
150
Number of variables
Number of variables
Number variables
Results ER2 graphs.
120
1.0
2.5
100
2.0
80
SKF1
0.6
Time((s)
40-
HC-QWO
0.5
0.2
grasp-QWO
0.0
20
grasp-BIC
0.0
100
150
50
100
150
Number of variables
Number of variables
Number variables
Results on ER3 graphs.
120
1.0
2.5
100
SKF]
0.6
PSHD
Time(s)
1.0
0.4
HI-QWO
0.0
gnas--qW grasp BIC
-0.5
0.0
Number of variables
Number of variables
Number of variables
Results on E R4 graphs.
Figure 4: Results for our method compared BIC on linear models with Gumbel noise.
17

---

C Formal Proofs
Theorem 4.4 (Characterizing B(X)). For any B B(X), there exists unique orthogonal matrix Q such that B = 1 - -QW, and vice versa. That is,
((X) - - QWI QQ" E DD).
(11)
Proof. Recall that
( [B E D Cov(X) (I -I- B) B)
Cov(X) = USU and W US- Therefore, we have
W = US.
(14)
First, suppose B B(X), i..., E D Cov ( = (1 B) -T- B)-T. We need to show that there exists an orthogonal matrix Q such that B - QW, or equivalently, QW 1 Since W is invertible, that is equivalent to showing that A = (I B)W-1 is orthogonal. We have
AAT ((T B) ((T B)W 1)² =(I B--B)W-11g Now-ww-2e T( -
B)".
(15)
Furthermore, Equation (14) implies that
W W- = USA = USU Cov(() = (T- B) ET B) B)T
(16)
By substituting the expression for W W-W--1 from Equation (6) into Equation (15). we have
AA = (I - B) B)" E E D..
Hence, A is orthogonal. The uniqueness of Q for a given B follows from the equation Q = (I B)W
Now suppose B I QW, where QQT Dn. We need to show that B E B()). To this end, suffices to show that (I (TT)) E Dr. Since I-B = QW and Cov (X) W WW-W-T we have
(I - (I - B) QWCOr(X(QW)) = WW Q =QQ E DD.
This completes the proof.
Theorem 4.5 (Soundness of Algorithm Ih. Under Assumption I and given the correct whitening matrix W as imput, matrix Q, the output of Algorithm 1 is the unique solution to (1)). Consequently, the returned graph corresponds to the true G defined in Definition 3.1
Proof. To establish this theorem, we first demonstrate the existence of a unique matrix Q that satisfies both constraints: diag (WW I and P.QW pf is upper triangular. Subsequently, we prove that the output of Algorithm is equal to this matrix.
Suppose that Q satisfies both conditions. We proceed by induction on to show that the i-th row of the matrix PaQ, denoted by ( is unique. The constraint that P.QW pf is upper triangular implies that vi < j (T(()) I Wx(j))
Induction base, i.e., when tors in (wn(2), (3)::..... (m)m)) which determines the unique direction of [n(1). Moreover, diag( (Pag(PPPPPTT
1: The vector (x(1) is orthogonal to all of the vec-
implies ((x(1) wa()) 1. These two conditions show that (x(1) is unique.
Induction step: suppose {(x(1), (m(2)..... (x(i-1)} (.())))...0(()))) are unique and we want to show that (mai) is unique. Note that (n()) is orthogonal to vectors in both sets (1x(1)) (a(2). (x(()))))) and [Wn((+1)) (=(i=))..... (n(n)) }, and vectors in the second set orthogonal to the first set and also each other. This shows that all of the vectors in both sets are linearly independent, and the direction of (m(s) is unique. Therefore, diag(P. Pf)
implies the uniqueness of (x())-
Now, to prove the theorem, it suffices to show that the matrix Q, the output of Algorithm sat- isfies diag(P, (a,,(, RW P!) I and P,QW pf is upper-triangular. For the first one, note that in
18

---

proju(u). For matrix line 6 of the algorithm, vectors (m(() are normalized such that (qm(1), Wa(1)) = 1 which implies diag( Ps QW PT) I. For the second condition, note that for each 1 S i n, vector (=(a) is the nor- malized residual of the projection of Wn(s)) into the space of (Wr(i+1)..... Wa(n)) }. which shows the orthogonality of (m(i) to wx()) with j > i. This property implies Pa P&QWPT is upper- apprr-triangllr. M
Theorem 4.2. Under Assumption D. for any permutation π E [[(())), there exists unique B E B(X)
such that (BB) is compatible with и. Furthermore, for this B, G(B) = g"
we showed that (12) has a unique solution and then proved that the out- put of Algorithm is equal to this solution. To complete the proof, it suffices to prove that graph G(P, C(P&QWPP P equal G where Q the output of Algorithm I To show this, we need to prove for each i <
Proof. In Theorem 4.5
x(() II x (x()) -U
((x()) (x(s) =
(17)
To prove (17). we first need a few notations. For a matrix A, Ali 5] denotes the matrix consisting of rows + ........ ... of matrix A. Similarly, A ii] denotes the i-th row of A. For vector and a set of vectors u = fu1, U2,.... " (roja(v) and (veu(v) denote the projection vector e on the span of u, and the residual of this projection, respectively, ie., res)(v) A, res A (v) denotes the projection of vector on the space spanned by the rows of A.
Now, let fix 1 < i < < In the rest of the proof, will prove Equation (7 for i, Recall that Coo () UST denotes the SVD decomposition
covariance matrix X. Define
A [11 j], B= Ul 1: n], C= W(1 s). D wi + :l:
Note that AA and BB" In-ji where In denotes the k k identity matrix. Furthermore,
C = T=AA-uU" D BS ttr
(18)
We further define E be the matrix constructed as follows:
Elk] res (CTA)T)) L1 vi
(19)
To prove (17). we present the following two lemmas.
Lemma C.1. (Block Matrix Inversion Lemma [BerOg]) Suppose T is invertible matrix, which is in the following form.
Ti
T12
T
T2i
T22
In this case, T-1 can be computed using the Schur complement of T as follows:
M
TT'TiNN
T
-N
N-1
where N = T22 - T2 Ti2 is the Schur complement of T11 in T and is assumed to be invertible,
and M = (T11 Ti2T2a T21)-1
Lemma C.2. For matrix E defined in the following holds:
EET Cov(X
(20)
Recall that X = [X,,..,,, .,)T².
Proof. We will prove that (EFTT)-1 = Cov(X We define
P DT (DD")-1D,
which is the projection matrix that projects any vector to the space of rows of D. Therefore, we have
== C- (PC")
C(1 (()-)
(21)
19

---

Next, we calculate DD", P. and EE" in terms of A, B, and S. First, (88 implies the following.
BS BS -
=S U (BS
B" BS-BBB
(22)
Using
22 we have
P= 1=
(23)
Using
21
we have
FF PCC (PT[T-P))) C)-- P) C(T PJC
(24)
=C - CPC
Note that in (24), we used [ () P P. which holds true because P is projection
matrix. further refine (44, we apply (23 to calculate CC and CPC the following.
1 CPC AS 00' BB b0)) (A)
CC AS U (AS - =
AS
AS Bd (BS B BS AMT
Applying the last equations in (24). we have
EET AS-1A AT AS (BS-
(25)
Next, we present Cov(X) in terms of A, B.
A
AS-1 AT
AS-1BT
=
B
S
B
=
BS-1A"
BS-
(26)
Applying Lemma C.1 to Cov(X)-1 with (6).. we have
- - BS-1 -N- ATA-- -AAS-1A") AS-
M
Cov(X)=
N-1
where
M -AS AS- - -B
(27)
and N is matrix that we do not need to calculate. Note that M = Cov(X since M is an j x matrix. Furthermore, (27 and (25) imply that M (EE?)-1 Therefore, (EET)-1 = Cov(X )
which concludes the proof.
Now to prove the theorem we use classic result on the linear Gaussian data. Consider O is the inverse of the covariance matrix for some variables with joint Gaussian distribution, then two variables X; and Xj are independent given all the other variables if and only if O;, = 0 [KF09]. Considering this result, to show (17), it is sufficient to show the following:
(grl)), wx(s)) + (Cov(Xi))~ =0
To show the above equation, let D denotes the set (wij+11) Wx()+2) Wn(n) }. then (r() resp- (wx))) Thus, (x(j) is orthogonal to the span of the vectors in D", then we have
(Wx()) (x()) (Wa()), reSD= (wx())) = (resD= (wa()) (wa(g) (w
(respa (w=(())(=())))))) resp= (wx()))
the value of (resp (wa(s) resps (w+())) is equal to (Cov(X and this result completes the proof of the theorem.
Based on lemma C.2.
20

---

Lemma 4.6. If the block between the idxi-th and idx,-th positions of π is modified, the vectors (=(M) for k < idx; or k idx, remain unchanged.
Proof. Let x' be the new permutation obtained after modifying the permutation и. First let k > idx,- this case, foreach j > idx,, Wn(j) wx'()) thus, v remains the same for both n and Therefore, dr(k) (x'()))-
Now let idx;- As we discussed in the main text, due the construction matrix Q in Al- gorithm ! for any 1 I i I n, the span of vectors {Tm(t) (=(++1): (x(n) is equal to the span of vectors (t), (x(++1), -.((++))...).. ... (x(n)) ). Furthermore, since we just modified the block between the idx and itx positions of π to obtain n², the two sets (=(k + 1), alk + .),.... T(n)) and (r'(k +1, x'kk)))((k +),+ 2), ....: I r'(n)} } are equal. Therefore, the sets (W=(k-1) 1)- (*(k++)) (x(n)) and (w+()++) (k+(k))): (()) } are also equal. Hence, the residual of vector (a(kk) these two sets is equal, which implies that (m(k) = (=((K)- This completes the proof.
Theorem 4.7 (Time complexity of Algorithm D QWO algorithm as implemented in Algorithm has the following time complexities:
On²)) "in")for for initially computing 9 without optional arguments.
(n? d) when called with optional arguments to update G", where d = idx, idx.
Proof. In Algorithm D the steps outside the for loop include initialization of variables and computing G(T QW) are executed in ((n²). ). Inside the for loop, the time complexity for calculating each Yy is (n((n - i)) because each term in the summation is done in Onn). As iterates from to n,
this complexity is O(n²) and for the entire loop it takes O(idxx, idx;) )²²) Therefore, the time
complexity for the initial step is (n²)), ), and for the update steps (n?dd)
21

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The main claims made in abstract and introduction include introducing a novel with theoretical guarantee for causal discovery in linear Gaussian models with good practical results, which all of these have addressed in the main text, sections 4, 5, 6
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize to other settings.
Is fine include aspirational goals as motivation as long it clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
The assumptions that are needed for our method to work is presented, for example linear Gaussian model, faithfulness assumption, etc.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "L_imitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The au- thors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to ad- dress problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers díscover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proof's
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
22

---

Answer: [Yes]
Justification: For each theory result in the paper, the assumptions are clear and the full proof is provided either in the main text appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: The implementation details described in the experiments section of main text or appendix of the paper and are enough reproduce the experimental results.
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of model checkpoint, other means that appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to re- produce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case thors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
23

---

Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: The full code of experiments along with an instruction to reproduce the results is provided.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (httpss//niis..c/
pubbic//uuides//ooeuuunssss wtt.oooewwwsiiiiiiooooooooo..... for more details.
While we encourage the release of code and data, understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
/ wttp:/ccppccccccdeeeeeeepeeeeeee /tteeeee/ooeemiiseceeea asionollcyy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g.. data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Details about the parameters and experimental setup are provided in the ex- periments section.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined other appropri- ate information about the statistical significance of the experiments?
Answer: [Yes]
Justification:
The confidence intervals are shown in all of the figures including experimental
results.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
24

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. The authors should prefer- ably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetrie error bars that would yield results that are out of range (e.g. negative error rates).
If error bars reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification; We have reported the execution time in all of our experiments, reported both in Section 5 and Appendix B. We have further provided the details required reproduce the experiments.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform. in every respect, with the NeurIPS Code of Ethics httpp////oooooooooccccoo//c/tttt
Answer: [Yes]
Justification:
We have followed the NeurIPS code Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative so- cietal impacts of the work performed?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that there is socictal impact of the work performed.
25

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact spe- cific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it legitimate to point out that an improvemnnt in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitiga- tion strategies (e.g.. gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines restrictions to access the model or implement- ing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do require this, but encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (c.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: The only used existing asset is the codes in the public python library causal- learn which is cited in the min text of the paper.
Guidelines:
The answer NA means that the paper does
existing assets. The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, possible, include a URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
26

---

If assets are released, the license, copyright information. and terms of use in the pack- age should provided. For popular datasets, pwwerswithcoseoo thoode cow /ataattsett has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is available online, the authors encouraged reach out the asset's creators.
13. New Assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: The only new asset is the codes of experiments which are provided alongside the paper.
Guidelines:
The answer NA means that the paper does release assets.
Researchers should communicate the details of the dataset/code/moodel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but the main contri- bution of the paper involves human subjects, then much detail possible should be included the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
27

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
28


Reviews:
[{"review_id": "OoBKCJcRgX", "reviewer": "Reviewer_9HgY", "summary": "This paper considers the problem of speeding up permutation-based causal discovery in linear Gaussian acyclic models. A typical permutation-based causal discovery algorithm includes two components: 1) constructing a DAG permitting a given topological ordering, and 2) a search strategy over the space of permutations. While most existing work focuses on component 2), this work focuses on 1). Specifically, the authors first characterize the adjacency matrices' equivalence class using whitening transformation and orthogonal rotation. Then, with this intuition, an algorithm QW-Orthogonality (QWO) based on Gram-Schmidt algorithm is proposed to construct the DAG under a given ordering. This QWO algorithm has a time complexity of O(n3) for an ordering pi without side information.", "questions": "As in \"Weaknesses\".", "limitations": "N/A", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "+ It is novel to study the speeding up of component 2 (constructing a DAG permitting a given topological ordering) when most of the existing work considers the search strategy over permutations.\n\n+ The proposed algorithm can be integrated into various existing permutation-based algorithms and enjoys a lower time complexity in updating steps.\n\n+ The review on existing work regarding permutation-based algorithms and the overall writing is clear.", "weaknesses": "+ **Assumptions needed is not spelled out:** What is the exact assumptions needed for this work? Is it normal faithfulness, or sparsest Markovian representation as in [RU18], or something in between as in [LAR22]? Since a primary focus of permutation-based algorithms is to relax assumptions needed, I would suggest authors put assumptions explicitly (and check if they are sufficient, or sufficient and necessary). Since different lemmas/theorems may need different assumptions, it is better to specify assumptions for each of them separately.\n\n+ **The significance of the proposed method may need a better justification:** The authors defined the whole equivalence class $\\mathcal{B}$ that involves also cyclic graphs and those even with self-loops. The introduction of orthogonal transformation also reflects this. However, in this work actually only acyclic graphs are considered (as seen in the upper triangular constraints in (13)). In this case, why do the authors take a (seemingly) detour for orthogonal transformation, instead of directly apply Cholesky/LDL transformation on the permuted covariance matrix? The time complexity of this is also O(n3). With some small modification to (constrained) Cholesky decomposition, I guess an updating complexity of O(n2d) can also be achieved.\n\n+ **A more clearer justification for the necessity of speeding up Gπ construction is needed:** For example, in line 138, \"to solve (9), ..., Note that the complexity of a brute-force search over all subsets U is exponential.\" I don't quite get this -- why do we need to try all subsets as parents candidates, given that by definition (5), the existence of each edge can directly be seen without any sense of subset traversing -- though the time complexity of directly using (5) is generally O(n5) (O(n2 for testing each edge) and O(n3) for each Pearson correlation calculation).\n\n+ **Some missing references:** For example, the idea to use orthogonal rotation to characterize the equivalence class is very similar to https://arxiv.org/abs/1910.12993. The similarities and differences should be discussed.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their detailed and thoughtful comments.\n\n---\n> Assumptions needed is not spelled out\n\nWe agree with the reviewer that permutation-based algorithms aim to relax necessary assumptions, such as faithfulness. However, for our method, the sparsest Markovian representation assumption in [RU18] is indeed sufficient. We did not mention this in the text to keep it simple and avoid exceeding the page limit. Since we have an extra page for the revised version, we will explicitly mention this lesser version of faithfulness. Additionally, for the sake of completeness, we will refer to the exact assumptions in the main results of the paper.\n\n---\n> The significance of the proposed method may need a better justification\n\nWe acknowledge the reviewer's insights on the matter. While there may be alternative methods based on Cholesky decomposition as the reviewer suggested, our method offers some advantages. For instance, the update step in our method is straightforward and can be easily integrated into existing search methods. Furthermore, we propose a characterization over $\\mathcal{B}(X)$ in Theorem 4.3, which is defined for cyclic models as well, paving the way for learning cyclic models. For example, by removing the upper-triangularity constraint in Equation 13, it would be a valid optimization for learning a cyclic graph (though we might need a stronger version of faithfulness to extend this to cyclic models). However, solving the optimization for cyclic models could be a challenging problem, but it is a promising direction for future work.\n\n\n---\n> A more clearer justification for the necessity of speeding up $\\mathcal{G}^{\\pi}$ construction is needed\n\nWe believe there has been a misunderstanding regarding lines 138-141, which we include here for reference:\n\n> To solve (9), various approaches have been proposed. Note that the complexity of a brute-force search over all subsets $\\mathbf{U}$ is exponential. Instead, the state-of-the-art search methods apply the grow-shrink (GS) algorithm [ARSR+23] on the candidate sets $\\mathbf{U}$ to find the parent set of each variable. These methods require computing the score function $S$, $O(n^2)$ times.\n\nHerein, we provide a literature review of the existing methods for constructing $\\mathcal{G}^{\\pi}$. We start by pointing out that brute force is not a feasible method. Therefore, we mention alternative approaches that involve computing the score function $S$ only $O(n^2)$ number of times. We then introduce different score functions for $S$ in the following paragraph. This gives us the computation complexity of the existing work for constructing $\\mathcal{G}^{\\pi}$ (please refer to Table 1). Finally, we justify our proposed method by the fact that it gains a speed-up of $O(n^2)$ in comparison to the aforementioned class of approaches.\n\n---\n> Some missing references\n\nWe appreciate the reviewer for bringing up that paper. It looks very interesting and relevant to our problem. In the revised version, we will include a brief discussion on the similarities and differences between our work and this paper.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "vTPjdVteFQ", "reviewer": "Reviewer_D3Yo", "summary": "Authors propose a new causal discovery algorithm on permutation-based methods in the context of Linear Gaussian Acyclic Models (LiGAMs). Specifically, authors focus on the computation complexity of existing solutions and propose a novel QW-Orthogonality (QWO) that improve the efficiency of computing a new graph $\\mathcal{G}$ for a given permutation $\\pi$. The computational complexity of QWO is $\\mathcal{O}(n^2)$, that is, significantly better than BIC-based alternatives.", "questions": "- Why comparing PDAGs instead of CPDAGs when computing evaluation metrics?", "limitations": "- Authors do not discuss explicitly the limitations of the proposed method, they refer to the theoretical assumptions instead.\n- The assumptions are the usual ones that are present in every causal discovery algorithm.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- The outline of the proposed solution is clear: the computational complexity of alternative solutions is a strong limitation for the applicability of causal discovery methods. Recasting the optimization problem is the key of the contribution and it's a original contribution.\n- The quality of the paper is relevant: methods are self-contained, they are clear and easy to follow.\n- The significance of the proposed solution is evident from the experimental setup, where existing solutions exceed the running time caps.", "weaknesses": "- The proposed solution is tested in combination of just two search-based procedures, which is a rather limited evaluation, especially if we observe that the gap between HC-BIC and GRASP-BIC is significant.\n- GRASP-BIC usually achieves better results then the proposed solution in random graphs.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's comments and are pleased that they found our approach to be a clear and original contribution.\n\n---\n> Why comparing PDAGs instead of CPDAGs when computing evaluation metrics?\n\nWe thank the reviewer for pointing out this issue/typo. In our experiments, we have indeed compared the CPDAGs (complete PDAG achieved after applying Meek rules). We will correct this typo in the revised version.\n\n---\n> GRASP-BIC usually achieves better results then the proposed solution in random graphs.\n\nWe acknowledge the reviewer's observation that GRASP-BIC usually archives slightly higher accuracy in random graphs.\nHowever, we note that (i) the difference in accuracy is almost negligible, and (ii) the primary goal of our method is to reduce computational complexity, which we have demonstrated both theoretically (by a factor of $O(n^2)$) and empirically. Figure 1 shows that GRASP-BIC was scalable to nearly 50 variables and only in sparse graphs, while our method was easily scalable to over 150 variables, even in denser graphs.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "QkWG9RShXw", "reviewer": "Reviewer_rQrH", "summary": "The authors present an efficient method for evaluating a score for score-based causal discovery over LiGAMs. Their method uses the whitening matrix $W$, derived from the observed covariance matrix, as a summary statistic. Their method is $\\mathcal{O}(n^2)$ faster than the classical *BIC* method, where $n$ is the number of observed variables. \n\nBecause $W$ only needs to be calculated once, the authors' score can be evaluated equally efficiently for any number of samples $N$. This presents an advantage against state-of-the-art methods such as *BDeu* and *CV General*, which become difficult for large $N$ and are primarily intended for nonlinear models. \n\nTheir score for each topological ordering is the fewest edges a DAG can have while being consistent with at least one LiGAM that produces the observed covariance matrix.", "questions": "Can QWO be extended to incorporate finite-sample uncertainty while retaining (some of) the advantage in the scalability of the method? For instance, could a confidence set over scores be output with the same scalability? \n\nPerhaps not and QWO is intended primarily for the setting with large $N$ where other methods are too computationally expensive?", "limitations": "Yes, the required assumptions are stated.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The authors present an elegant and intuitive approach with a well-motivated score function. They write with clarity and demonstrate the clear advantage of their method over the state-of-the-art in simulated and real-world experiments.", "weaknesses": "The authors' method does not address finite-sample uncertainty and is best thought of as providing a (discrete) point estimate of an underlying population score that is a function of the population whitening matrix $W^*$. It would have been impressive for the authors to handle finite-sample uncertainty in $W^*$ as a part of their approach. \n\nIn *BIC*, for example, the likelihood of the data given a graph is explicitly part of the score function. As the number of samples $N$ increases, *BIC* will provide increasingly lower scores to DAGs outside of the Markov equivalence class of the true DAG $G^*$. This does not happen for QWO since $W$ is treated as a population whitening matrix either way.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their interesting comment regarding incorporating finite-sample uncertainty. We also appreciate the positive feedback on our method's clarity, intuitive design, and superiority in experiments.\n\n---\n> Can QWO be extended to incorporate finite-sample uncertainty while retaining (some of) the advantage in the scalability of the method? For instance, could a confidence set over scores be output with the same scalability?\n\nBelow, we drive an uncertainty analysis for each edge in $\\mathcal{G}^{\\pi}$.\nIn our experiments, our method incorporates a point estimate of $W$ using a finite set of samples. This estimation can indeed be noisy, which is the primary source of error for the final estimation. After constructing $q_1, q_2, ..., q_n$, we then proceed to check the orthogonality of vectors  $w_1, ..., w_n$ and $q_1, ..., q_n$.\nWe can show that for $i<j$, the dot product of $w_{\\pi(i)}$ and $q_{\\pi(j)}$ is an estimation for partial correlation $\\rho_{X_{\\pi(i)}, X_{\\pi(j)} |{X_{\\{\\pi(1), \\pi(2), \\ldots, \\pi(j-1)\\}\\backslash\\{\\pi(i)\\}}}}.$\nTo check whether this dot product is zero, we then apply a z-test for this partial correlation.\nUsing results in [1], we can accordingly derive the following confidence interval:\n$$ P(|\\text{Error of estimation}|>\\epsilon) \\leq \\frac{1}{ (\\text{NumOfSamples} -j -2) \\epsilon^2},\n$$\nwhere `Error of estimation' is the distance between our estimated dot product and the actual dot product.\n\n[1] Drton, Mathias, and Michael D. Perlman. \"Multiple testing and error control in Gaussian graphical model selection.\" (2007): 430-449.\n\n---\n> QWO is intended primarily for the setting with large $N$ where other methods are too computationally expensive?\n\nOur method has the advantage of being fast with large sample sizes, but it can also be applied in situations with small sample sizes. It is worth noting that there is extensive literature on computing the inverse covariance matrix with a small number of data points. For example, see [2, 3].\n\n[2] Ravikumar, Pradeep, et al. \"High-dimensional covariance estimation by minimizing l1-penalized log-determinant divergence.\" (2011): 935-980.\n\n[3] Friedman, Jerome, Trevor Hastie, and Robert Tibshirani. \"Sparse inverse covariance estimation with the graphical lasso.\" Biostatistics 9.3 (2008): 432-441.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "BaAAdIA8H0", "reviewer": "Reviewer_Znq2", "summary": "This paper focuses on permutation-based methods for causal discovery in Linear Gaussian Acyclic Models (LiGAMs). A new method called QWO is proposed to improve the efficiency of computing a causal graph given a permutation. Compared with baselines, QWO achieves superior performance.", "questions": "(See above)", "limitations": "Yes. The assumptions are restricted, such as the linear Gaussian model and faithfulness assumption.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "1.\tAccording to the theoretical analysis, the proposed method is guaranteed to learn the true graph for a given permutation when there are sufficient samples.\n2.\tThe computational complexity is quadratic, much lower than other methods, such as the BIC-based method.\n3.\tQWO can be applied to some other methods as an accerating module.\n4.\tThe experiments are sufficient to support their claims.", "weaknesses": "1.\tThe writing is poor. The format for definitions is not consistent throughout the paper. Some definitions are directly given in the content body (Lines 64-72), while some others are put in a ‘Definition’ framework (Definition 2.1, 3.1, 3.2). Some contents are repetitive. Line 89 and Line 94 both mention that “faithfulness holds”.\n2.\tSome definitions are not properly expressed. In Definition 3.3., how come “there is an edge from $X_{\\pi(i)}$ to $X_{\\pi(j)}$ if and only if $X_{\\pi(i)}$ and $X_{\\pi(j)}$ are conditionally independent”? \n3.\tAssumptions are not clearly demonstrated. Despite faithfulness, Markov assumption is also required for inferring causal graphs.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for acknowledging our method's theoretical guarantees, efficient computational complexity, and of our experimental results.\n\n---\n> The format for definitions is not consistent throughout the paper. Some definitions are directly given in the content body (Lines 64-72), while some others are put in a ‘Definition’ framework (Definition 2.1, 3.1, 3.2).\n\nWe used the ‘Definition’ framework for the following cases:\n\n- Definition 2.1 $G(B)$\n\n- Definition 3.1 $[\\mathcal{G}]$\n\n- Definition 3.2 $\\mathcal{B}(X)$\n\n- Definition 3.3 $\\mathcal{G}^{\\pi}$\n\n- Definition 4.2 Whitening matrix $W$\n\nThese notations are either new and crucial for our presentation, or there is no consensus for them in the literature. On the other hand, the notations and definitions provided in lines 64-72 (which is the notations section) are commonly used in various fields, and we did not see the necessity of assigning a separate box for them.\n\n---\n> Line 89 and Line 94 both mention that “faithfulness holds”.\n\nWe will edit line 89.\n\n---\n>  In Definition 3.3., how come “there is an edge from $X_{\\pi(i)}$ to $X_{\\pi(j)}$ if and only if $X_{\\pi(i)}$ and $X_{\\pi(j)}$ are conditionally independent”?\n\nWe thank the reviewer for pointing out this typo. The independency ($\\perp \\mathrel{\\mkern-9mu} \\perp$) in Equation (5) should be changed to dependency.\n\n---\n> Despite faithfulness, Markov assumption is also required for inferring causal graphs.\n\nWhen the underlying model is a structural equation model (SEM), which is a more general setting of our problem, the Markov property holds (please refer to Theorem 1.2.5 in [1]). Therefore, Markov property is not technically an assumption. To avoid any confusion, we will mention this in the revised version.\n\n[1] Pearl, Judea. Causality. Cambridge university press, 2009.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: All reviewers agree on the merits of this manuscript, that framing learning a CPDAG for Linear Gaussian models as an optimization problem is interesting and worthwhile by itself, as well as that it may lead to novel approaches to adjacent problems. There is also agreement that the write-up needs improvement, such as clearly stating the assumptions under which guarantees can be given. Similarly, the authors should more carefully explain what can be done for cyclic resp. acyclic models, and support this with empirical evidence. Last, the authors are strongly recommended to clearly discuss the limitations of the method beyond merely stating the assumptions. I would recommend them to also evaluate and discuss how well the method works when the assumptions are not met (as they probably wont hold in practice anyway). All in all, the general sentiment for this paper is positive and hence an accept is in order. Please use the additional page to address the further improve the manuscript, especially the three points above.
Decision: Accept (poster)