QWO: Speeding Up Permutation-Based Causal Discovery in LiGAMs

Authors: mohammad  shahverdikondori
ehsan  mokhtarian
negar  kiyavash
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 20:04:39

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.

Main paper: 
 introduction: 
 Causal discovery is fundamental to understanding and modeling the relationships between variables in various scientific domains [Pea09, SGSH00]. A causal graph, typically represented by a directed acyclic graph (DAG), is a graphical model that represents how variables within a system influence one another [Pea09] . The problem of causal discovery refers to learning the causal graph from available data [SGSH00, MT99, Chi02, FZ13, ZARX18, BSSU20a, MAGK21, LAR22, MAJ + 22, MEAK24], which has broad applications across numerous fields such as economics [Hec08] , genetics [SMD + 03], and social sciences [MW15] .Score-based methods are an important class of approaches for causal discovery [JS10, KGG + 22, HS13, HZL + 18, Sch78, Bun91] . They involve defining a score function over the space of DAGs and seeking 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 Ω(n 2 ) , where n is the number of variables, poses a significant challenge [Rob73] . Ordering-based methods, introduced by [TK12a] , significantly reduce Table 1: Time complexity comparison of various methods for computing and updating G π in LiGAMs. Here, n is the number of variables, d is the length of the updated block of the permutation, and k is the number of folds considered in k-fold cross-validation.Method QWO BIC BDeu CV General Initial Complexity O(n 3 ) O(n 5 ) Ω(n 3 N log(N )) O( n 2 N 3 k 2 ) Update Complexity O(n 2 d) O(n 4 d) Ω(n 2 dN log(N )) O( ndN 3 k 2 )the search space of score-based methods to 2 O(n log(n)) by considering only topological orderings of DAGs, rather than the DAGs themselves [FK03, ZNC19, MKEK22, SLOT22] . Building on this premise, various permutation-based methods have been proposed that further refine the search for causal structures [SdCCZ15, LAR22, SWU21] . These methods have two main ingredients:1. Constructing G π : A module that for any given permutation π constructs a specific DAG, denoted by G π . We will formally define G π in Definition 3.1, but roughly speaking, this DAG represents the best structure that can be inferred from the data while conforming to the given permutation π. 2. Search over π: The work by [RU18] (which proposed the Sparsest Permutation algorithm) showed that the correct permutation minimizes the number of edges in G π . Accordingly, permutation-based methods are equipped with a search strategy that, using the module mentioned above for computing G π , searches through the space of permutations to minimize the number of edges in G π .In recent years, various search methods have been developed to enhance the accuracy, robustness, and scalability of the algorithm for the second part, i.e., searching over the permutations [LAR22, TBA06, MKEK22, TK12b] . 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 G π multiple times during the algorithm's execution. To construct G π or update G π after a permutation update, most methods utilize a decomposable score function and identify the parent set that maximizes the score of each variable. Table 1 presents a time complexity comparison of different methods for computing and updating G π . We will delve into a more detailed discussion of these methods and their computational complexity in Section 3.In this paper, we focus on linear Gaussian acyclic models (LiGAMs), an important class of continuous causal models. We propose a novel method called QW-Orthogonality (QWO) for computing G π for a given permutation π, which significantly enhances the computational efficiency of this module. Specifically, as shown in Table 1 , QWO's complexity is independent of the number of data points N . Moreover, it speeds up computing/updating G π by O(n 2 ) compared to the state-of-the-art 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 G π accurately for any given permutation π when sufficient samples are available. 2. 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 O(n 2 ) faster than the state-of-the-art BIC-based method. 3. Compatibility with existing search methods: After updating a permutation, QWO efficiently updates G π to improve its compatibility with established search methods. In our experiments, we combine QWO with both GRaSP [LAR22] and a hill-climbing-based search technique [TBA06, SG06] , showing its superior performance in terms of time complexity. 
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.[n] denotes the set of natural numbers from 1 to n, and Π([n]) denotes the set of all n! permutations over [n] . For a permutation π ∈ [n], P π denotes the permutation matrix, such that for any n × n matrix A, the product P π AP T π permutes the rows and columns of A according to the permutation π. The identity matrix is denoted by I, and its dimension is implied by the context. ∥A∥ 0 denotes the number of non-zero entries in matrix A. The set of n × n diagonal matrices is represented by D n . ⟨a, b⟩ denotes the inner product of vectors a and b. For a set X = {X 1 , . . . , X n } and S ⊆ [n], we define X S = {X i |i ∈ S}. diag(A) denotes a diagonal matrix with the same diagonal elements as matrix A.In a directed graph (DG), edges are directed and may form cycles. A DG without cycles is called a directed acyclic graph (DAG). If there is a directed edge from X to Y , then X is a 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), Pa G (X), and Anc G (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). π is a topological order for G, or equivalently G is compatible with π, if for any edgeEQUATIONwhere N ∼ N (0, Σ) is a Gaussian random vector with covariance matrix Σ. Equivalently, we haveEQUATION)The causal graph of M, denoted by G(B), is a DG with the adjacency matrix corresponding to the support of B T . 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, Σ) is given byEQUATIONFor distinct indices 1 ≤ i, j ≤ n, and S ⊆ X, we use X i ⊥ ⊥ X j |X S to denote that X i and X j are independent conditioned on X S . The notion of d-separation defined over DAGs is a graphical criterion to encode conditional independence (CI) within a graph. We similarly use X i ⊥ ⊥ X j |X S to denote that X i and X j are d-separated given X S in a DAG. For a formal definition of d-separation, see [Pea88] . In a 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 faithfulness assumption posits that if two variables are conditionally independent in the distribution, then they are d-separated in the DAG. Together, the Markov property and faithfulness establish a 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 a 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 ′ , P ) also satisfies the Markov property andG ′ / ∈ [G], it holds that |G ′ | > |G|.Here, |G| denotes the number of edges in G.A search strategy over the space of permutations to solve (5).Given that minimization in (5), search strategies typically traverse the space of permutations and greedily update the permutation. To maintain computational efficiency, these methods use a module to update G π incrementally rather than recomputing it from scratch. Moreover, good search strategies are ideally consistent in permutation-based causal discovery. A consistent search method guarantees to solve (5) as long as it is equipped with a module that correctly computes G π . In Appendix A, we review two search strategies: GRaSP [LAR22] and hill-climbing-based methods [TBA06, SG06] . GRaSP is consistent, ensuring that it finds the correct permutation that minimizes the number of edges in G π . In contrast, hill-climbing methods do not provide consistency guarantees but are shown to be efficient in practice.Existing methods for computing G π are primarily score-based and involve the following optimization:EQUATIONwhere S is typically a decomposable score function, summing individual scores for each variable given its parents within the graph:EQUATIONTo find the parent set of each variable in G π , the following optimization is performed for all 1 ≤ i ≤ n:EQUATIONTo solve (8), various approaches have been proposed. Note that the complexity of a brute-force search over all subsets U is exponential, which makes it impractical. The final score for a variable given its parents is the average log-likelihood evaluated on the test sets using the regression functions learned from the training data. The CV General method has an initial complexity of O( n 2 N 3 k 2 ) and an update complexity of O( ndN 3 k 2 ). This method typically uses a 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 a speed-up of O(n 2 ) compared to BIC.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 a separate "Limitations" 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 authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a 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 address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover 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 important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
permutation-based causal discovery in ligams: 
 In this section, we discuss permutation-based methods for causal discovery in LiGAMs. We begin by formalizing our assumptions, which will be referenced throughout the remainder of the paper. Assumption 1. Let M * = (B * , Σ * ) be a LiGAM that generates the random vector X = [X 1 , X 2 , . . . , X n ] T . Let D be an N × n data matrix, where each row is an i.i.d. sample from X. We assume that the pair (G * , P * ) satisfies the (SMR) assumption, where G * = G(B * ) is the causal graph of M * and P * is the joint distribution of X.Under Assumption 1, our goal is to learn the causal graph G * from the observational data D. However, using mere observational data-even when the number of samples N approaches 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 is to identify [G * ], a problem known to be computationally NP-hard [CHM04] .As briefly discussed in the introduction, score-based methods aim to minimize a 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 a topological order of a DAG in [G * ]. Definition 3.1 (G π ). Under Assumption 1, for any arbitrary permutation π ∈ Π([n]), we denote by G π = (X, E π ) the unique DAG constructed from P * as follows: For distinct indices 1 ≤ i, j ≤ n, there is an edge from X π(i) to X π(j) in E π if and only if i < j and X π(i) ̸⊥ ⊥ X π(j) |X {π(1),π(2),...,π(j-1)}\{π(i)} .(4) 1. Constructing G π : A module that for a given permutation π constructs G π .Remark 3.1. For any π, G π is compatible with π. Furthermore, if a DAG G ∈ [G * ] is compatible with π, then G π = G. [RU18]Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?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 or 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 crossreferenced.• All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a 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. 
qwo: 
 In this section, we present a novel approach for computing G π 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(X)). For a random vector X, we denote by B(X) the set of coefficient matrices of LiGMs that generate X, i.e.,B(X) = {B|∃Σ ∈ D n : Cov(X) = (I -B) -1 Σ(I -B) -T }.Note that in this definition, the causal graphs corresponding to the elements of B(X) can be cyclic, but we restrict our attention to a subset of B(X) with acyclic corresponding graphs. We reformulate causal discovery in LiGAMs as follows: arg min In the following, we establish the relationship between B(X) and G π . Theorem 4.2. Under Assumption 1, for any permutation π ∈ Π([n]), there exists a unique B ∈ B(X)such that G(B) is compatible with π. Furthermore, for this B, G(B) = G π .All proofs are provided in Appendix C. Theorem 4.2 implies that to compute G π , we can directly find the unique B ∈ B(X) whose corresponding graph G(B) is compatible with π. Next, we propose a characterization for B(X) using whitening transformation [Fuk90, HO00, KLS18]. Definition 4.3 (Whitening matrix W ). Let Cov(X) = U SU T be the singular value decomposition (SVD) of Cov(X), where S is a diagonal matrix including the singular values of Cov(X) on its diagonal and U is an orthonormal matrix (i.e., U U T = I). The whitening matrix W is defined asEQUATION)We note that the 'W' in QWO corresponds to the whitening matrix W . Whitening is a linear transformation N I := W X that transforms the Gaussian random vector X to another Gaussian random vector N I , where Cov(N I ) = I. Furthermore, for an arbitrary orthogonal matrix Q (i.e., QQ T ∈ D n ), if we define N Q := QN I = QW X, then it is straightforward to show that N Q is also a Gaussian random vector with mean zero and a diagonal covariance matrix. Therefore,(I -QW, N Q ) is a LiGM that generates X since X = (I -QW )X + N Q .In the following, we show that such LiGMs create all possible LiGMs that generate X.The QWO module for computing and updating G π1: Function QWO(π, W , [Optional] idx l , [Optional] idx r , [Optional] Q) 2:Default values for optional arguments: idx l = 1, idx r = n, Q = any arbitrary n × n matrix 3: Denote the columns of Q T by {q i } n i=1 and the columns of W by {w i } n i=1 4: for i from idx r to idx l do 5:r i ← w π(i) - n j=i+1 ⟨w π(i) ,q π(j) ⟩ ∥q π(j) ∥ 2 2 q π(j) 6: q π(i) ← ri ⟨ri,w π(i) ⟩ 7: end for 8: G π ← G(I -QW ) 9: Return G π , Q Theorem 4.4 (Characterizing B(X)).For any B ∈ B(X), there exists a unique orthogonal matrix Q such that B = I -QW , and vice versa. That is,EQUATIONBy combining Theorems 4.2 and 4.4, to learn G π , it is 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:EQUATIONIn this subsection, we introduce QW-Orthogonality (QWO), our proposed approach for computing G π in LiGAMs. The function QWO in Algorithm 1 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, we consider the case where these optional arguments are not provided, and they are initialized as follows: idx l = 1, idx r = n, and Q to be an arbitrary n × n matrix. The goal of this function is to create a matrix Q that is the unique solution of (12) and subsequently to compute G π .Denote the columns of Q T by {q i } n i=1 and the columns of W by {w i } n i=1 . Since W = U S -1 2 U T , all the eigenvalues of W are positive, and thus {w i } n i=1 are n linearly independent vectors in an n-dimensional space. Furthermore, because Q is an orthogonal matrix for each pair (i, j), 1 ≤ i ̸ = j ≤ n, q i ⊥ q j . The condition P π QW P T π is upper triangular in (12) implies that q π(i) should be orthogonal to w π(j) if j > i. Therefore, we need to find {q i } n i=1 such that this condition holds. Note that each zero in ∥I -QW ∥ 0 corresponds to an orthogonality between {q i } n i=1 and {w i } n i=1 . With this intuition, we propose our approach for constructing Q, which maximizes the number of aforementioned orthogonality. To do so, we apply the Gram-Schmidt algorithm [GL13] to the vectors {w i } n i=1 in the following order: π(n), π(n -1), ..., π(1). In Algorithm 1, we iteratively for each i, compute r i , the residual of projecting w π(i) on the span of {w π(i+1) , . . . , w π(n) }. This is equivalent to projecting w π(i) on the span of {q π(i+1) , . . . , q π(n) }, which are orthogonal to each other. q π(i) is set to normalized residual r i , which ensures that the product of q π(i) and w π(i) is 1, and consequently, diag(P π QW P T π ) = diag(QW ) = I. Finally, we use Theorem 4.2 to construct G π = G(I -QW ) and return G π and Q. Theorem 4.5 (Soundness of Algorithm 1). Under Assumption 1 and given the correct whitening matrix W as input, matrix Q, the output of Algorithm 1 is the unique solution to (12). Consequently, the returned graph corresponds to the true G π defined in Definition 3.1.The optional arguments idx l , idx r , and Q allow for incremental efficient updates to π. Lemma 4.6. If the block between the idx l -th and idx r -th positions of π is modified, the vectors q π(k) for k < idx l or k > idx r remain unchanged.A consequence of Lemma 4.6 is that q π(k) for k < idx l or k > idx r from the previous computed Q can be reused and it 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 ← Compute the SVD decomposition of Cov(X), such that Cov(X) = U SU T 4: W ← U S -1 2 U T % Whitening matrix 5: π ← An initial permutation 6: G π , Q ← QWO(π, W ) 7: while F has not stopped do 8: π ′ ← Update π according to F, and let idx l and idx r be the leftmost and rightmost index of π that have been updated, respectively 9:G π ′ , Q ′ ← QWO(π, W , idx l , idx r , Q) 10: if G π ′has less number of edges than G π then 11:π, Q, G π ← π ′ , Q ′ , G π ′ 12:end if 13: end while 14: Return G π Theorem 4.7 (Time complexity of Algorithm 1). QWO algorithm as implemented in Algorithm 1 has the following time complexities:• O(n 3 ) for initially computing G π without optional arguments.• O(n 2 d) when called with optional arguments to update G π , where d = idx ridx l .Algorithm 2 demonstrates how QWO can be integrated into existing search methods such as GRaSP [LAR22] and Hill-Climbing [TBA06, SG06] . This algorithm integrates QWO into a simple 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 X and applies SVD to compute the whitening matrix. Starting from an initial permutation, it calls function QWO to compute Q and G π . Subsequently, the algorithm iteratively updates the permutation using the given search method F 1. For each updated permutation, it calls function QWO with optional arguments to compute the new Q and G π graph. Next, the algorithm checks if the new permutation is better than the previous one by checking whether the new graph has fewer edges. If the new permutation is better, it updates the permutation and proceeds to the next iteration. The search algorithm stops when its stopping criterion is satisfied and returns the final π and G π .Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it 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 is described in the experiments section of main text or appendix of the paper and are enough to 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 a 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 a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it 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 a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , 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 authors 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 to have some path to reproducing or verifying the results. 
experiments: 
 In this section, we present a comprehensive evaluation of QWO, which is designed for the module that computes G π in permutation-based causal discovery methods in LiGAMs. As discussed earlier, a permutation-based method consists of two components: a search method and a 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, SG06] . Please refer to Appendix B for the implementation details of these methods and additional results.We generated random graphs according to an Erdos-Renyi model with an average degree of i 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 a LiGAM (B * , Σ * ), we sampled the entries of B * from [-2, -0.5] ∪ [0.5, 2] and the noise variances uniformly from [1, 2]. Each reported number on the plots is an average of 30 random graphs. Two metrics were used to evaluate the performance of the methods:1Search method F could be any existing approach such as GRaSP [LAR22] or Hill-Climbing [TBA06, SG06] . • Skeleton F1 Score (SKF1): The F1 score between the skeleton of the predicted graph and the skeleton of the true graph, which ranges from 0 to 1, where 1 indicates perfect accuracy. • Complete PDAG SHD (PSHD): For the true DAG G and predicted DAG Ĝ, PSHD calculates the complete PDAGs of [G] and [ Ĝ] and evaluates the number of changes (edge addition, removal, and type change) needed to obtain one PDAG from the other, normalized by the number of nodes. This metric shows the distance between the Markov equivalence classes of the predicted graph and the true graph.We carried out our experiments in four settings:• Low-Dimensional Data: We evaluated the performance of QWO and other methods on small real-world graphs, namely ASIA [LS88] , CANCER [KN10] , SACHS [SPP + 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 of different methods using both search strategies are presented in Table 2 . As shown in the table, QWO demonstrates comparable accuracy to the other approaches for different graph structures while being significantly faster. • High-Dimensional 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 in this setting, they were not included in this experiment. Instead, we compared QWO with the BIC score using both GRaSP and Hill-Climbing search methods. As illustrated in Figure 1 , while maintaining high accuracy, QWO demonstrates a significant speedup over BIC for both search methods. • Non-Gaussian Noise Experiments: To test the robustness of our method, we conducted experiments 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 exponential and Gumbel distributions). The results of these experiments appear in Appendix B. Although 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 distributions.• 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 to the approach in [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.Question: Does the paper provide open access to the data and code, with sufficient instructions to 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. Guidelines:• The answer NA means that the paper does not 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 in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
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 G π for a given permutation by O(n 2 ), resulting in a significant speed-up compared to the state-ofthe-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 a 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, or 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 cyclic and acyclic models. Therefore, our characterization of B(X) in Theorem 4.4 holds for LiGMs as well. Combining this result with our reformulation of causal discovery in LiGAMs in (9), we can generalize the formulation to LiGMs as follows: arg minEQUATIONThis suggests that causal discovery in LiGMs can be approached by finding an orthogonal matrix Q that sparsifies I -QW . Some recent works have explored 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. Results on 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, we considered the initial permutation based on the size of Markov boundaries2 of the variables, which in the case of linear Gaussian models are equal to the number of non-zero elements in each row of the inverse covariance matrix [KF09] .To generate the data matrix D using a linear model (B * , Σ * ), we sampled the entries of B * uniformly from [-2, -0.5]∪[0.5, 2] and the noise variances uniformly from [1, 2] for all of the noise distributions i.e., Gaussian, Exponential, and Gumbel. 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: We have followed the NeurIPS code Ethics. 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA]Justification:Guidelines:• The answer NA means that there is no societal impact of the work performed.• 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 specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement 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 as 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 or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation 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). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification:Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing 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 is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 
licenses for existing assets: 
 Question: Are the creators or original owners of assets (e.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. 
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 as details about compensation (if any)? Answer: [NA] Justification: Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to 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. 
 
 Appendix: 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. Results on ER4 graphs.Figure 3 : Results for our method compared to BIC on linear models with exponential noise.Herein, we present the results of running QWO and the BIC on two linear non-Gaussian models: Exponential and Gumbel noise. Figures 3 and 4 show the time complexity and accuracy in terms of two metrics for both methods on both search strategies. Although 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 distributions. Theorem 4.4 (Characterizing B(X)). For any B ∈ B(X), there exists a unique orthogonal matrix Q such that B = I -QW , and vice versa. That is,EQUATION)Proof. Recall thatB(X) = {B|∃Σ ∈ D n : Cov(X) = (I -B) -1 Σ(I -B) -T },Cov(X) = U SU T , and W = U S -1 2 U T . Therefore, we haveEQUATION)First, suppose B ∈ B(X), i.e., ∃Σ ∈ D n : Cov(X) = (I -B) -1 Σ(I -B) -T . We need to show that there exists an orthogonal matrix Q such that B = I -QW , or equivalently, QW = I -B. Since W is invertible, that is equivalent to showing that A := (I -B)W -1 is orthogonal. We haveEQUATION)Furthermore, Equation ( 14) implies thatEQUATION)By substituting the expression for W -1 W -T from Equation ( 16) into Equation ( 15), we haveAA T = (I -B)(I -B) -1 Σ(I -B) -T (I -B) T = Σ ∈ D n .Hence, A is orthogonal. The uniqueness of Q for a given B follows from the equation Q = (I -B)W -1 . Now suppose B = I -QW , where QQ T ∈ D n . We need to show that B ∈ B(X). To this end, it suffices to show that (I -B)Cov(X)(I -B) T ∈ D n . Since I -B = QW and Cov(X) = W -1 W -T , we have(I -B)Cov(X)(I -B) T = QW Cov(X)(QW ) T = QW W -1 W -T W T Q T = QQ T ∈ D n .This completes the proof.Theorem 4.5 (Soundness of Algorithm 1). Under Assumption 1 and given the correct whitening matrix W as input, matrix Q, the output of Algorithm 1 is the unique solution to (12). 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(P π QW P T π ) = I and P π QW P T π is upper triangular. Subsequently, we prove that the output of Algorithm 1 is equal to this matrix.Suppose that Q satisfies both conditions. We proceed by induction on i to show that the i-th row of the matrix P π Q, denoted by q T π(i) , is unique. The constraint that P π QW P T π is upper triangular implies that ∀i < j : q π(i) ⊥ w π(j) . Induction base, i.e., when i = 1: The vector q π(1) is orthogonal to all of the vectors in {w π(2) , w π(3) , . . . , w π(n) } which determines the unique direction of q π(1) . Moreover, diag(P π QW P T π ) = I implies ⟨q π(1) , w π(1) ⟩ = 1. These two conditions show that q π(1) is unique.Induction step: suppose {q π(1) , q π(2) , . . . , q π(i-1) } are unique and we want to show that q π ( i) is unique. Note that q π ( i) is orthogonal to vectors in both sets {q π(1) , q π(2) , . . . , q π(i-1) } and {w π(i+1) , w π(i+2) , . . . , w π(n) }, and vectors in the second set are 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 q π(i) is unique. Therefore, diag(P π QW P T π ) = I implies the uniqueness of q π(i) . Now, to prove the theorem, it suffices to show that the matrix Q, the output of Algorithm 1 satisfies diag(P π QW P T π ) = I and P π QW P T π is upper-triangular. For the first one, note that in line 6 of the algorithm, vectors q π(i) are normalized such that ⟨q π(i) , w π(i) ⟩ = 1 which implies diag(P π QW P T π ) = I. For the second condition, note that for each 1 ≤ i ≤ n, vector q π(i) is the normalized residual of the projection of w π(i) into the space of {w π(i+1) , . . . , w π(n) }, which shows the orthogonality of q π(i) to w π(j) with j > i. This property implies P π QW P T π is upper-triangular. Theorem 4.2. Under Assumption 1, for any permutation π ∈ Π([n]), there exists a unique B ∈ B(X) such that G(B) is compatible with π. Furthermore, for this B, G(B) = G π . Proof. In Theorem 4.5, we showed that (12) has a unique solution and then proved that the output of Algorithm 1 is equal to this solution. To complete the proof, it suffices to prove that graph G(P π QW P T π ) is equal to G π , where Q is the output of Algorithm 1. To show this, we need to prove for each i < j,EQUATION)To prove (17), we first need a few notations. For a matrix A, A[i : j] denotes the matrix consisting of rows {i, i + 1, . . . , j} of matrix A. Similarly, A[i] denotes the i-th row of A. For a vector v and a set of vectors u = {u 1 , u 2 , . . . , u k }, proj u (v) and res u (v) denote the projection of vector v on the span of u, and the residual of this projection, respectively, i.e., res u (v) = v -proj u (v). For a matrix A, res A (v) denotes the projection of vector v on the space spanned by the rows of A.Now, let us fix 1 ≤ i < j ≤ n. In the rest of the proof, we will prove Equation ( 17) for i, j. Recall that Cov(X) = U SU T denotes the SVD decomposition of the covariance matrix of X. DefineA = U [1 : j], B = U [j + 1 : n], C = W [1 : j], D = W [j + 1 : n].Note that AA T = I j and BB T = I n-j , where I k denotes the k × k identity matrix. Furthermore,EQUATIONWe further define E to be the j × n matrix constructed as follows:EQUATIONTo prove (17), we present the following two lemmas. Lemma C.1. (Block Matrix Inversion Lemma [Ber09] ) Suppose T is an invertible matrix, which is in the following form.T = T 11 T 12 T 21 T 22In this case, T -1 can be computed using the Schur complement of T 11 as follows:T -1 = M -T -1 11 T 12 N -1 -N -1 T 21 T -1 11 N -1 ,where N = T 22 -T 21 T -1 11 T 12 is the Schur complement of T 11 in T and is assumed to be invertible, and M = (T 11 -T 12 T -1 22 T 21 ) -1 . Lemma C.2. For matrix E defined in (19), the following holds:EE T = Cov(X [j] ) -1 (20) Recall that X [j] = [X 1 , . . . , X j ] T .Proof. We will prove that (EE T ) -1 = Cov(X [j] ). We defineP = D T (DD T ) -1 D,which is the projection matrix that projects any vector to the space of rows of D. Therefore, we haveEQUATIONNeurIPS 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 score 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 in 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 can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction).• The answer NA means that the paper does not use 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, if 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.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets 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 not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new 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 not release new assets.• Researchers should communicate the details of the dataset/code/model 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. 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: + **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.

+ **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.

+ **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).

+ **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.
questions: questions: As in "Weaknesses".
limitations: limitations: N/A
paper_summary: paper_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.
strength_weakness: strength_weakness: - 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.
- GRASP-BIC usually achieves better results then the proposed solution in random graphs.
questions: questions: - Why comparing PDAGs instead of CPDAGs when computing evaluation metrics?
limitations: limitations: - Authors do not discuss explicitly the limitations of the proposed method, they refer to the theoretical assumptions instead.
- The assumptions are the usual ones that are present in every causal discovery algorithm.
paper_summary: paper_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. 

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

Their 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.
strength_weakness: strength_weakness: 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. 

In *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.
questions: 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? 

Perhaps not and QWO is intended primarily for the setting with large $N$ where other methods are too computationally expensive?
limitations: limitations: Yes, the required assumptions are stated.
paper_summary: paper_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.
strength_weakness: strength_weakness: 1.	The 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”.
2.	Some 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”? 
3.	Assumptions are not clearly demonstrated. Despite faithfulness, Markov assumption is also required for inferring causal graphs.
questions: questions: (See above)
limitations: limitations: Yes. The assumptions are restricted, such as the linear Gaussian model and faithfulness assumption.

Comments: comment: Many thanks for the clarifications. I maintain my score 7: Accept.
Thanks for the responses, and I will maintain my score.
Thank you taking your time to reply. I'm satisfied with the explanations provided.