Tensor-based second-order causal discovery
arXiv:2606.18074v1 [stat.ML] 16 Jun 2026
Nathan Ouyang, Kexin Wang and Anna Seigal Abstract. Causal discovery seeks to uncover the causal dependencies among variables. For this purpose, we propose an algorithm called Tensor-based Second-order Causal Discovery (TSCD). Its input is a tensor obtained from the covariance matrices of observational and interventional data. Assuming the causal dependencies follow a linear structural equation model on a directed acyclic graph (DAG), TSCD outputs the DAG and the functions on its edges, requiring only that the noise variables are uncorrelated. We also implement a version of the approach for nonlinear models. Our focus on second-order statistics (via the covariance matrices) is motivated by their statistical and computational efficiency relative to higher-order moments, their identifiability relative to first-order statistics, and that they work regardless of whether the variables are Gaussian. We show that TSCD has identifiable causal order and parameters from a number of interventions that is logarithmic in the number of variables. Experiments show that TSCD is robust to noise, competitive with existing methods, and scales to hundreds of variables.
1. Introduction Causal discovery is a fundamental problem in biology [FLNP00, SPP+ 05], neuroscience [SKPF22], economics [Imb20], public health [GPR99] and machine learning [SLB+ 21]. Its goal is to determine the causal dependencies among a set of variables and to quantify the strength of the causal effects. With this information, one seeks to understand the system mechanisms, make predictions and answer counterfactual questions. Structural equation models (SEMs) [Pea09, Wri34, PJS17] are a framework for describing causal relationships. In a SEM, each variable is a function of its direct causes, together with a noise variable. The causal structure is often represented by a directed acyclic graph (DAG) whose edges encode direct causal effects. Linear SEMs (LSEMs) are SEMs in which the functions on the edges are linear. We study both linear and nonlinear SEMs. Interventions are crucial in causal discovery because observational data alone in general identifies only a set of DAGs that induce the same conditional independence relations, known as a Markov equivalence class [VP90]. Perturbing variables alters their dependence on their parents, breaking observational symmetries. An example is Perturb-seq [DPL+ 16], where genetic perturbations are applied and their effects on gene expression are measured. In this paper, we propose a new algorithm for causal discovery, which we call Tensorbased Second-order Causal Discovery (TSCD). It recovers an unknown causal order, DAG, and model parameters, from the covariance matrices of observational and interventional data, assuming that the data follows an LSEM on an unknown DAG. The only distributional assumption required is that the noise variables are uncorrelated. We also propose an extension of the method to nonlinear SEMs, which we call TSCD-nonlinear. 1
2
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
TSCD is second-order because it takes as input only the covariance matrices of the data. A method that uses only second-order statistics is necessary for a system with multiple Gaussian noise variables. It is practical when noise variables are close to Gaussian, making it difficult to estimate higher-order statistics. Second-order statistics are practical in applications where per-context sample sizes are too small to reliably estimate full distributions. We see that such scenarios limit the accuracy of approaches that use non-Gaussianity assumptions, such as Linear Non-Gaussian Acyclic Model (LiNGAM) [SHHK06], or distributions for conditional independence testing, such as the Peter-Clark (PC) algorithm [SGS00]. Second-order statistics are also useful summaries for situations where raw data are unavailable due to privacy, storage, or communication constraints. They are more informative than first-order statistics such as mean vectors, from which the parameters in our model are not identifiable. Several works have used second-order statistics for causal discovery and parameter estimation in an LSEM. From observational data alone, the problem is well-known to be non-identifiable [PJS17]. To get around this, approaches include restricting the noise variances [PB14], applying nonlinear transformations [SB23], parameter identifiability for a fixed graph [FDD12, DFS11], and outputting an equivalence class of DAGs [Chi02]. By comparison, TSCD makes no assumptions on the graph, beyond that it is a DAG, and makes no extra assumptions on parameters, beyond those of a usual LSEM. Its only distributional assumption on the noise variables is that they are uncorrelated. To identify the causal order and model parameters, TSCD builds on a line of work to use the compatibility of causal structure across multiple contexts (also called views or environments) to recover parameters from a collection of covariance matrices [HCA+ 25, RHPM15, PBM16, SSBU23]. It studies observational data together with data generated via perfect interventions. TSCD is tensor-based because it works by constructing a tensor from the covariance matrices, whose decomposition encodes the causal relationships and intervention patterns. For p variables in k contexts, it is a p × p × k tensor that is symmetric in its first two indices; i.e., a stack of k many p × p symmetric matrices, see Figure 1. The stability of second-order statistics across contexts is a necessary condition to be a root [PJS17]. We use the tensor decomposition to devise a subspace-membership condition that is both necessary and sufficient. Subspace membership has well-studied connections to tensor decomposition, see [DL06, JLV23, KP25, Lan11, RSW26, WS23, WPKS25]. TSCD identifies candidate root nodes using a subspace-proximity criterion, refined using pairwise interventionasymmetric correlation tests. The tensor formulation yields identifiability guarantees, on the causal order and consequently on the parameters of the model. We show that ⌈log2 p⌉ interventional environments, and one observational environment, are sufficient and, in the worst case, necessary to recover the order and parameters. The paper is organized as follows. In Section 2, we give background, identifiability results, and related work. In Section 3, we show how second-order statistics across observational and intervention contexts are combined into a tensor, which we call the precision tensor. In Section 4 we prove our identifiability results. In Section 5, we explain how roots can be recovered. In Section 6, we present the TSCD algorithm. Section 7 evaluates TSCD on linear and nonlinear models and real datasets.
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
3
2. Setup and main results 2.1. Directed acyclic graphs and structural equation models. Causal relationships among a set of variables can be represented by a directed graph, where nodes are variables and directed edges encode the direct causal relationships between them. Definition 2.1. A directed acyclic graph (DAG) is a graph G = (V, E) with nodes V and directed edges E ⊆ V × V , such that G contains no directed cycles. If there is an edge j → i, we say j is a parent of i. We denote the set of parents of node i by pa(i). A root is a node with no parents. If there is a directed path j → ⋯ → i, we say that j is an ancestor of i. A structural equation model writes each variable as a function of its parents in the graph together with an exogenous noise variable. We consider DAGs with vertex set [p] = {1, . . . , p}. Definition 2.2. Given a DAG with vertex set [p], where node i represents variable xi , a structural equation model (SEM) writes xi = fi (xpa(i) , ϵi ),
for all i = 1, . . . , p,
where fi is a function, xpa(i) is the set of variables with direct edges to xi , and ϵi is an exogenous noise variable. A linear structural equation model (LSEM) is an SEM in which each function fi is linear: p
xi = ∑ Λij xj + ϵi ,
i = 1, . . . , p,
j=1
with Λij the weight of the dependency xj → xi (or 0 if there is no such edge). We say that an LSEM follows the structure of a DAG if Λij ≠ 0 if and only if j ∈ pa(i). Writing x = (x1 , . . . , xp )T and ϵ = (ϵ1 , . . . , ϵp )T , the LSEM can be written as x = Λx + ϵ or x = (I − Λ)−1 ϵ, where Λ ∈ Rp×p is the matrix of edges weights. Denote the covariance of exogenous noise variables by Ω. Let x have covariance matrix Σ and precision matrix Θ ∶= Σ−1 . Then (1)
Σ = (I − Λ)−1 Ω(I − Λ)−T
and Θ = (I − Λ)T Ω−1 (I − Λ),
see e.g. [BP02, Sul23]. The noise variables being uncorrelated implies that Ω is diagonal. The fact that the LSEM follows the structure of a DAG means there is a permutation of indices that makes Λ strictly lower triangular. This ensures that (I −Λ) is invertible. Finding the reordering that makes Λ lower triangular is the causal discovery problem. The uniqueness of causal order and recovered parameters in the model is important for interpretability and downstream analysis. Given a known causal order, the parameters in the model are identifiable from the covariance Σ, e.g. using the LDL decomposition. However, the matrix Λ in an LSEM on an unknown causal order is not identifiable from Σ. This can be seen by observing that an LDL decomposition of P ΣP T exists for every permutation matrix P . Recovery of the causal order is sufficient for recovery of the DAG and model parameters: subsequent steps, e.g. via regression, can detect the weights on the edges in the total DAG on a causal order, and hence whether edges are present or absent. Example 2.3. Consider the DAG x1 → x2 with SEM x1 = ϵ1 ,
x2 = 2x1 + ϵ2 ,
1 0 where Ω = ( ). 0 1
4
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
Now consider the reversed DAG x2 → x1 , with SEM x2 = ϵ̃2 ,
x1 = 52 x2 + ϵ̃1 ,
1 0 where Ω̃ = ( 5 ) . 0 5
1 2 Both models have covariance matrix Σ = ( ). 2 5 2.2. Interventions. Interventions can recover the parameters in LSEMs from covariance matrices identifiably. We consider perfect interventions [EGS05], which remove the influence of parents on the intervened variables. Definition 2.4 (Perfect intervention). Fix an LSEM x = (I − Λ)−1 ϵ, where ϵ is a vector of noise variables and Λij is the weight on the edge xj → xi . Let I ⊆ [p] be a set of intervened nodes. A perfect intervention on I removes all incoming edges to nodes in I, i.e., sets Λij = 0
for all
i∈I
and j ∈ pa(i),
and changes ϵi to ϵ′i for i ∈ I, such that the p noise variables remain mutually uncorrelated. Hence the intervention at I can change Ω at positions (i, i) for i ∈ I. We perform multiple interventions on a single LSEM, and call each interventional setting an intervention context. The set of intervened nodes in each intervention context is assumed to be known. 2.3. Identifiability. We present the identifiability results for our model. The input consists of the covariance matrices in the observational context and multiple intervention contexts. Assumption 2.5. The vector of variables x = (x1 , . . . , xp )T follows an unknown LSEM x = Λx + ϵ, on an unknown DAG G, with noise variables ϵ = (ϵ1 , . . . , ϵp )T mutually uncorrelated with nonzero variances, where Λij is the edge weight on the edge xj → xi . In the intervention contexts, perfect interventions are performed at a known nonempty subset of variables. We have access to the covariance matrix Σ1 of data in the observational context and covariance matrices Σ2 , . . . , Σk of data under intervention contexts. Under Assumption 2.5, we show that second-order statistics suffice to recover the causal structure from one observational and ⌈log2 p⌉ interventional contexts. Theorem 2.6 (Identifiability from intervention patterns). Consider an unknown LSEM x = Λx + ϵ, under Assumption 2.5. For each xi , define its intervention-pattern vector bi ∈ {0, 1}k by 0 (bi )j = { 1
xi intervened in context j otherwise.
If the vectors b1 , . . . , bp are nonzero and pairwise distinct, then the causal order and adjacency matrix Λ are identifiable. We prove Theorem 2.6 by forming a tensor from the covariance matrices under observational and intervention contexts. Its decomposition encodes the rows of I − Λ and the intervention-pattern vectors. Identifiability follows from the uniqueness of its decomposition. The proof is in Section 5.
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
5
Corollary 2.7. Given an observational context, ⌈log2 p⌉ perfect-intervention contexts are sufficient and in the worst case necessary to recover the DAG and adjacency matrix Λ. Corollary 2.7 follows by observing that ⌈log2 p⌉ intervention contexts suffice to assign distinct intervention-pattern vectors to all nodes. With fewer contexts, there exist two variables with the same intervention-pattern vectors, making the edge between them nonidentifiable. See Section 5. 2.4. Algorithmic overview. The proof of Theorem 2.6 suggests an algorithm for recovering the causal order and parameters of the model: stack the covariance matrices across the contexts into a tensor and decompose it. TSCD performs this tensor decomposition, with two modifications. First, modified precision (inverse covariance) matrices are stacked across contexts, rather than the covariance matrices, since they yield a lower rank decomposition that encodes the DAG and intervention patterns. Second, the algorithm does not compute a full decomposition, since this is less numerically stable. Instead, it uses the sparsity of Λ and the known intervention-pattern vectors to build the decomposition step by step. TSCD uses a subspace-proximity score that, in the population setting, takes the value 1 only at roots of the DAG. In finite samples, it keeps a few high-scoring candidates and uses correlations to choose among them. After selecting a root, it removes it and repeats. 2.5. Related Work. The PC algorithm [SGS00] is a constraint-based method that recovers the Markov equivalence class of a DAG using conditional independence tests. It iteratively removes edges based on conditional independence relations and orients edges using logical rules. Greedy equivalence search (GES) [Chi02] learns a DAG by optimizing a score (e.g., the Bayesian Information Criterion) computed over the space of equivalence classes using a greedy forward-backward procedure. Its extension to interventional data, greedy interventional equivalence search (GIES) [HB12], incorporates interventions into the scoring framework to improve identifiability. Interventional greedy sparse permutation (IGSP) [WSYU17] learns a causal DAG from observational and interventional data with unknown intervention targets by searching over permutations and selecting the sparsest graph consistent with conditional independence and invariance constraints across environments. It is nonparametric but relies on accurate conditional independence testing. The linear non-Gaussian acyclic model (LiNGAM) via Independent Component Analysis (ICA) [SHHK06] identifies a linear SEM by assuming non-Gaussian noise. It estimates the mixing matrix via ICA and recovers the causal order by permuting and scaling components to match a DAG structure. Under the assumptions that the noise variables are independent and non-Gaussian, ICA-LiNGAM is identifiable. It breaks down when some noise variables are Gaussian or dependent. DirectLiNGAM [SIS+ 11] improves upon ICA-LiNGAM by estimating a causal order through iterative regression: at each step, it identifies an exogenous variable by testing independence between residuals and candidate predictors, removes its effect, and repeats. Though more stable than ICA-based methods, it still relies on independence testing and non-Gaussianity.
6
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
NOTEARS [ZARX18] learns a DAG by solving a continuous optimization with a differentiable acyclicity constraint. It assumes access to individual samples rather than only second-order statistics and relies on optimization over a nonconvex objective. Invariant causal prediction (ICP) [PBM16] identifies causal parents of a variable by searching over subsets of variables and selecting those for which the conditional distribution of the target is invariant across environments. This exhaustive subset search is computationally expensive and requires that there is at least one intervention on each parent. We do not compare to this method in our synthetic experiments because existing implementations are not publicly available in Python and the method is designed for parent identification of a single target variable rather than full DAG recovery. The paper [RSW21] shows that commonly used data-generating processes can introduce artifacts that make the causal learning problem artificially easy. Heuristics such as sorting variables by their marginal variances can recover the causal order well under synthetic settings, not because they capture causal structure, but because the data generation encodes ordering information. To account for this, we include a sort-and-regress baseline (SortRegress) in the experiments, which orders variables by their marginal variances and then performs regression. This checks that our method does not rely on such artifacts. Our approach differs from the above methods in several ways. Unlike constraint-based methods such as PC and IGSP, it does not require access to individual samples and does not use conditional independence testing. Unlike ICA-LiNGAM and DirectLiNGAM, we only require that the noise variables are uncorrelated and do not assume independent or nonGaussian noise variables. Our method is more computationally efficient than GES/GIES for many variables, as the search space becomes large for a score-based greedy search. It is also more computationally efficient than NOTEARS for many variables, since the constrained optimization calculates an expensive matrix exponential. Unlike invariance-based approaches such as ICP, our method does not require exhaustive subset search or interventions that target each parent individually. 3. The precision tensor We construct the tensor that is the core of our identifiability proof and algorithm. It is obtained by combining precision matrices across contexts. For k contexts and p variables, it is a tensor T ∈ Rp×p×k that is symmetric in its first two factors. Its decomposition encodes the DAG and the intervention-pattern vectors. For related ideas based on stacking covariance matrices rather than precision matrices, see [WBP+ 26]. We stack precision matrices instead of covariance matrices because they lead to a decomposition involving the rows of I − Λ rather than (I − Λ)−1 . The former stays the same on the rows indexed by variables that are not intervened, but the latter does not. Definition 3.1 (Outer product). For u ∈ Rp , the outer product u ⊗ u is the matrix uuT , which has (i, j) entry ui uj . For u ∈ Rp and v ∈ Rk , the outer product u ⊗ u ⊗ v ∈ Rp×p×k is (u ⊗ u ⊗ v)ijℓ = ui uj vℓ . Given a tensor T ∈ Rp×p×k , we denote its entry at position (i, j, ℓ) by Tijℓ or T (i, j, ℓ). We say that a tensor T is symmetric under swapping the first two indices, if Tijℓ = Tjiℓ . for all i, j ∈ [p] and ℓ ∈ [k]. For example, this holds for the tensor u ⊗ u ⊗ v. Definition 3.2 (Tensor slices). For a tensor T ∈ Rp×p×k , fixing one index produces a matrix called a slice of the tensor. We denote them by T (i, ∶, ∶) ∈ Rp×k , T (∶, j, ∶) ∈ Rp×k ,
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
7
and T (∶, ∶, ℓ) ∈ Rp×p . For example, T (i, ∶, ∶) is the matrix with (j, ℓ) entry T (i, j, ℓ). If T is symmetric under swapping the first two indices, then T (i, ∶, ∶) and T (∶, i, ∶) are the same. For more on tensors and their slices, see e.g. [KB09]. We stack the decompositions in (1) together across contexts. To do so, we begin with the following observation, recorded here as a proposition for clarity, which restates the expression for the precision matrix. Let ei ∈ Rp denote the i-th standard basis vector. Proposition 3.3. Let Θ be the precision matrix from the LSEM x = (I − Λ)−1 ϵ on a DAG, where the noise covariance matrix Ω is diagonal with non-zero entries. Let vi denote the i-th row of (I − Λ), and let ωi be the i-th diagonal entry of Ω−1 . Then p
Θ = ∑ ωi vi⊗2 .
(2)
i=1
Node i is a root of the DAG if and only if vi = ei . Proof. The last statement follows because Λij = 0 for all j ∈ [p] exactly when i is a root. The matrix Θ has decomposition (I − Λ)T Ω−1 (I − Λ), see (1). Its (j, ℓ) entry is p □ ∑i=1 ωi (I − Λ)ij (I − Λ)iℓ . The product (I − Λ)ij (I − Λ)iℓ is the (j, ℓ) entry of vi⊗2 . We derive a related expression for the precision matrix under an intervention context. The nodes that are intervened become roots, as the effect of their parents is set to zero. This zeros out the corresponding rows of Λ. The rows of Λ that are not intervened stay the same. Proposition 3.4. Under the same setting as Proposition 3.3, let Θj be the precision matrix under a perfect intervention on Ij ⊆ [p]. Then Θj = ∑ ωi vi⊗2 + ∑ ωi′ e⊗2 i , i∈Ijc
i∈Ij
where ωi′ are new noise parameters. Proof. The precision matrix can be decomposed as p
−1 ⊗2 Θj = (I − Λj )T Ω−1 j (I − Λj ) = ∑(Ωj )ii (I − Λj )i,∶ , i=1
as in Proposition 3.3. It remains to relate the matrices Λj and Ωj to Λ and Ω respectively. Under a perfect intervention on node ℓ, all incoming edges to ℓ are removed, i.e., the entries of the ℓ-th row of Λ are zeroed out. Thus Λj is zero on rows i ∈ Ij and equals Λ on rows ′ i ∉ Ij . Finally, (Ω−1 j )ii is ωi for i ∈ Ij and ωi for i ∉ Ij , since the intervention can alter the variance of the exogenous noise variables at intervened variables. □ Though the above matrix decompositions are not identifiable, the recovery of certain parameters is unique: the coefficients of the vectors ei . Proposition 3.5 (Rank reduction). Let intervention context j ≥ 2 be a perfect intervention on a nonempty set Ij ⊆ [p]. For each i ∈ Ij , there is a unique scalar µ such that rank (Θj − µ e⊗2 i ) = rank (Θj ) − 1. Moreover, we have µ = (Σ1j )ii .
8
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
Proof. The matrix Θj − µei eTi drops rank by one if and only if µ=
1 eTi Θ−1 j ei
,
by Wedderburn rank reduction [Wed34]. The denominator of µ is eTi Σj ei = (Σj )ii .
□
Corollary 3.6. Under Assumption 2.5, we can recover the set of matrices ⎧ ⎫ ⎪ ⎪ ⎪ ⎪ ⊗2 M = ⎨ ∑ ωi vi ∶ j = 1, . . . , k ⎬ . ⎪ ⎪ c ⎪ ⎪ ⎩i∈Ij ⎭ Proof. The precision matrix Θj is a sum of two types of term: those of the form ωi vi⊗2 , ⊗2 for i ∉ Ij , and those of the form µe⊗2 i where µ is unique to (j, i). The coefficients of ei can be recovered by Proposition 3.5, thus the terms of the second type can be removed. □ Definition 3.7 (Precision tensor). We stack the matrices in M into a tensor T ∈ Rp×p×k with slices T (∶, ∶, j) = ∑i∈Ijc ωi vi⊗2 . We call T the precision tensor. It has decomposition p
(3)
T = ∑ ωi vi ⊗ vi ⊗ bi , i=1
where bi ∈ Rk has (bi )j = 1 if node i is not intervened in context j, and 0 otherwise. We call bi the intervention-pattern vector of the i-th variable. We collect the intervention-pattern vectors bi together to form the intervention-pattern matrix B ∈ Rp×k with (i, j) entry (bi )j , which is 0 if variable i is intervened in context j and 1 otherwise. 4. Identifiability of the precision tensor decomposition The decompositions of the individual precision matrices (2) are not unique: multiple nonequivalent choices of vectors vi and scalars ωi yield the same matrix Θ, cf. Example 2.3. In comparison, the decomposition of the precision tensor (3) is unique, under mild conditions. This typifies the benefit of combining second-order statistics across contexts. The decomposition of the precision tensor T is a CP decomposition: it writes a tensor as a sum of outer products. The decomposition is said to be unique if any such expression has summands that coincide with vi ⊗ vi ⊗ bi , up to scale and permutation [CO12]. Once the terms vi ⊗ vi ⊗ bi are recovered, the vectors vi and bi are obtained uniquely, up to scale. Lemma 4.1. The decomposition of T is unique when b1 , . . . , bp are pairwise distinct. Proof. The decomposition of T is unique when no pair of vectors bi , bj are collinear and v1 , . . . , vp are linearly independent, by the Kruskal criterion [Kru77]. Linear independence of {vi }pi=1 follows from the invertibility of I −Λ. Since each bi is binary and nonzero, collinearity only occurs if the vectors coincide. □ Theorem 2.6 follows by combining Lemma 4.1 with the structure of I − Λ. Proof of Theorem 2.6. Consider the tensor decomposition (3). The vectors vi are linearly independent since they are rows of I − Λ, which is invertible. By Lemma 4.1, the decomposition of T is identifiable. So, I −Λ is identifiable up to row scaling and permutation. The ambiguity is resolved since each vi is the i-th row of I − Λ, whose i-th entry is one. □
B 1 1 1 1 1 1 1 1 1 1
0 0 1 1 1 1 1 1 1 1
1 1 0 0 0 0 1 1 1 1
1 1 0 0 1 1 0 0 1 1
0 1 0 1 0 1 0 1 0 1
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
9
Figure 1. Constructing the precision tensor for 10 variables and five contexts (one observational). Top row: Precision matrices for each context. Middle row: The matrices after rank reduction. Bottom row: The matrices are stacked to form the precision tensor. The boxed entries are the only ones that change in the rank reduction step; they correspond to the zeros of the interventionpattern matrix B ∈ {0, 1}10×5 with rows bTi . For example, the boxed entries in the fifth matrix reflect that column 5 of B has zeros in entries i = 1, 3, 5, 7, 9. The resolution of the scaling and permutation ambiguity in the proof of Theorem 2.6 also appears in LiNGAM [SHHK06]. Corollary 2.7 quantifies how many intervention contexts are sufficient and in the worst case necessary for unique decomposition. Proof of Corollary 2.7. With one observational context, each intervention-pattern vector bi ∈ {0, 1}k has first entry equal to 1. Hence, with k − 1 additional intervention contexts, there are at most 2k−1 distinct binary patterns available. To ensure that all p nodes have distinct intervention patterns, it suffices that 2k−1 ≥ p, which holds when k − 1 ≥ ⌈log2 p⌉. In this case, we can assign distinct binary vectors bi to each node, and define the intervention sets Ij by intervening on node i in context j whenever (bi )j = 0. By Theorem 2.6, this ensures identifiability of Λ. For necessity, assume we have fewer than ⌈log2 p⌉ additional interventions, then there are two nodes, say 1 and 2, that share the same intervention pattern. Their contributions to each slice T (∶, ∶, i) of the tensor T in (3) appear only through M = ω1 v1⊗2 + ω2 v2⊗2 .
10
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
We show that this decomposition is not unique. Suppose M is only nonzero in the top-left 2 × 2 submatrix a b M[2] = ( ). b c Then there exist distinct pairs (v1 , v2 ) that produce the same M[2] . For example, v1 = (0, 1, 0, . . . , 0), v2 = (a, b, 0, . . . , 0),
and v1 = (1, 0, 0, . . . , 0), v2 = (b, c, 0, . . . , 0).
Hence, the edge direction between the first two nodes is not identifiable.
□
5. Tensor-based root selection The previous section demonstrates the identifiability of the causal order, and parameters, from the precision tensor. In principle, the parameters of the model can be recovered via usual tensor decomposition of (3), e.g. using [WPKS25, SVBDL15]. However, in practice such tensor decomposition suffers from numerical instability, worsened by the fact that the intervention-pattern vectors are often close to parallel, as most variables are not intervened. To avoid this, TSCD does not compute a full decomposition. Instead, it uses the sparsity of Λ and the known intervention-pattern vectors to build the decomposition. We show how root nodes can be identified from the precision tensor. Recall that we denote by T (i, ∶, ∶) ∈ Rp×k the slice of T obtained by fixing the first index to take the value i. The slices T (1, ∶, ∶), . . . , T (p, ∶, ∶) span a linear space, which we denote by A. The form of the decomposition in (3) implies that A = span{vi ⊗ bi ∶ i = 1, . . . , p} ⊆ Rp×k . The connection between uniqueness of tensor decomposition and of the rank one matrices or tensors in a linear space is well-studied, see e.g. [Lan11]. In our case, it is as follows. Lemma 5.1. The only rank-one matrices in A are vi ⊗ bi up to scaling for i = 1, . . . , p. Proof. Assume for contradiction that there exists a rank-one matrix a ⊗ b in A that is not collinear to any vi ⊗ bi . Then a ⊗ b = ∑pj=1 µj vj ⊗ bj , for some µj ∈ R. Without loss of generality, we assume µ1 ≠ 0. Hence we can build a new decomposition of T that uses a ⊗ b instead of v1 ⊗ b1 . That is, there exists ci ∈ Rp for i = 1, . . . , p such that T = c1 ⊗a⊗b+∑pj=2 cj ⊗vj ⊗bj . The uniqueness result in Lemma 4.1 was stated for a decomposition in which the two vectors of length p are the same (a partially symmetric decomposition), but the decomposition is still unique if all vectors are allowed to differ. The new decomposition of T is thus a contradiction. □ Lemma 5.2. Node i is a root if and only if ei ⊗ bi ∈ A. Proof. Node i is a root if and only if it has no incoming edges. This implies vi = ei , so ei ⊗ bi ∈ A. Conversely, suppose ei ⊗ bi ∈ A. Then there exist scalars {αj }pj=1 such that p
(4)
ei ⊗ bi = ∑ αj vj ⊗ bj . j=1
The vectors {vj }pj=1 are linearly independent, since they are rows of the invertible matrix I − Λ. Thus, there exists a dual set of vectors {wℓ }pℓ=1 such that ⟨wℓ , vj ⟩ = δℓj . Multiplying both sides of (4) by wℓT yields ⟨wℓ , ei ⟩bi = αℓ bℓ .
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
11
For ℓ ≠ i, the vectors bℓ and bi are not collinear, so the above equality implies ⟨wℓ , ei ⟩ = 0 and αℓ = 0. Thus only the term ℓ = i remains, and we obtain ei ⊗ bi = αi vi ⊗ bi . Since bi ≠ 0, this implies ei = αi vi . Comparing the i-th entries (both equal to 1) gives αi = 1, hence vi = ei . Therefore, node i is a root. □ The subspace-membership condition ei ⊗ bi ∈ A is equivalent to ei ⊗ ei ⊗ bi being a rankone summand of T . This says the coefficient of ei ⊗ ei in each precision matrix Θj is the same across all contexts j in which i is not intervened, cf. Proposition 3.5. Thus, ei ⊗ bi ∈ A implies that the variance of node i is stable across the contexts in which it is not intervened. We compare our approach to the simpler idea of checking the consistency of the variance of a variable across contexts. Doing so, one recovers the roots of the DAG together with variables that look like roots based on the pattern of interventions, as follows. Lemma 5.3. Under Assumption 2.5, fix i ∈ [p]. Consider an intervention context in which neither i nor any of its ancestors are intervened on. Then the variance of xi is the same as in the observational context. Proof. If none of the ancestors of i are intervened on, then the structural equations for i and all of its ancestors are unchanged, and so the joint distribution of i and its ancestors is unchanged. Since xi is a function only of its ancestors and its own noise variable, and since i is not intervened on, the distribution of xi is unchanged. Hence its variance is unchanged. □ We compare Lemmas 5.2 and 5.3. Stable variance is only a necessary condition for a node to be a root, since non-root nodes whose ancestors are never intervened also satisfy this property. In comparison, the condition ei ⊗ bi ∈ A is necessary and sufficient for i to be a root, because it also incorporates information from contexts where i is intervened. 6. Algorithm We recover a causal order of the nodes in LSEMs with perfect interventions. This is achieved in two steps: a root candidate selection procedure (Section 6.2) followed by a pairwise refinement step (Section 6.3). The two-step algorithm design is validated by the ablation experiment in Appendix B.7. Estimation of the matrix Λ given the ordering is a standard problem and is deferred to the appendix. First we explain how to construct the sample precision tensor. 6.1. The sample precision tensor. So far, we have assumed access to the true population covariance matrices, and used these to construct the precision tensor. In practice, we recover a causal order from finite samples. Let X1 ∈ Rn1 ×p be the observational dataset and X2 ∈ Rn2 ×p , . . . , Xk ∈ Rnk ×p be the interventional datasets. We construct the sample precision tensor in Algorithm 1, see Figure 1 for an illustration. The sample precision tensor only requires access to the sample covariance matrices Cov(X1 ), . . . , Cov(Xk ) ∈ Rp×p , not to individual data samples. Later for root selection we need the number of samples in each context, but again not the individual samples.
12
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
Algorithm 1: Sample precision tensor Input: Sample covariance matrices for contexts 1, . . . , k; intervention-pattern matrix B ∈ Rp×k . Output: The sample precision tensor. for j = 1, . . . , k do Mj ← Cov(Xj )−1 ; for i = 1, . . . , p do if Bij = 0 then 1 ei eTi ; Mj ← Mj − Cov(Xj )i,i end end end Form the tensor T = [M1 ∣⋯∣Mk ] ∈ Rp×p×k ; return T ;
6.2. Greedy root candidate selection. We use a greedy procedure to identify root nodes from the precision tensor. Recall that root nodes are those satisfying ei ⊗ bi ∈ A, see Lemma 5.2. In practice, to quantify how close a node is to satisfying this condition, we compute a subspace-proximity score from projection onto the subspace A: αi =
(5)
∥PA (ei ⊗ bi )∥F ∈ [0, 1], ∥bi ∥
where PA denotes orthogonal projection onto A. Lemma 6.1. Node i is a root if and only if αi = 1. Proof. Node i is a root if and only if ei ⊗ bi ∈ A, by Lemma 5.2. This is equivalent to ∥PA (ei ⊗ bi )∥F = ∥ei ⊗ bi ∥F . i )∥F Since ∥ei ∥ = 1, we have ∥ei ⊗ bi ∥F = ∥bi ∥, so this is equivalent to ∥PA (e∥bi ⊗b = 1. i∥
□
In the population setting, the top scores equal 1 and exactly identify the root nodes. The following result shows that if a node’s projection score given by (5) is high, then the node is close to being a root, in the sense that vectors ei and vi have similar directions. Lemma 6.2. Let A = span{vj ⊗ bj ∶ j = 1, . . . , p}, where vj , bj are as in Proposition 3.3. Let √ w1 , . . . , wp satisfy ⟨wℓ , vj ⟩ = δℓj . Let αi be the projection score from (5) and define δi = 1 − αi2 . Assume the bi are separated, in that for every j ≠ i, (6)
inf ∥bi − cbj ∥ ≥ γi ∥bi ∥ c
for some
γi > 0.
Then ∥wj ∥∥vj ∥ ) δi . γi j≠i
∥ei − ⟨wi , ei ⟩vi ∥ ≤ (∑
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
13
Proof. When αi is close to 1, the matrix ei ⊗ bi is well approximated by a linear combination of the matrices vj ⊗ bj . There exist scalars c1 , . . . , cp such that p
∥ei ⊗ bi − ∑ cj vj ⊗ bj ∥ ≤ ∥bi ∥δi ,
(7)
j=1
since
F
√ dist(ei ⊗ bi , A) = ∥ei ⊗ bi ∥ 1 − αi2 = ∥bi ∥δi .
Multiplying the matrix in (7) by wℓ , we obtain p T ∥⟨wℓ , ei ⟩bi − cℓ bℓ ∥ = ∥wℓ (ei ⊗ bi − ∑ cj vj ⊗ bj )∥ j=1 p
≤ ∥wℓ ∥ ∥ei ⊗ bi − ∑ cj vj ⊗ bj ∥ j=1
2
p
≤ ∥wℓ ∥ ∥ei ⊗ bi − ∑ cj vj ⊗ bj ∥ ≤ ∥wℓ ∥∥bi ∥δi ,
(8)
j=1
F
where the equality follows from ⟨wℓ , vj ⟩ = δℓj . By the separation assumption (6), (9)
inf ∥⟨wℓ , ei ⟩bi − cbℓ ∥ = ∣⟨wℓ , ei ⟩∣ inf ∥bi − cbℓ ∥ ≥ ∣⟨wℓ , ei ⟩∣γi ∥bi ∥. c
c
Combining (8) and (9) gives ∣⟨wℓ , ei ⟩∣γi ∥bi ∥ ≤ ∥wℓ ∥∥bi ∥δi which implies ∣⟨wℓ , ei ⟩∣ ≤
∥wℓ ∥ δi . γi
Finally, using the dual-basis expansion ei = ∑pℓ=1 ⟨wℓ , ei ⟩vℓ , we get ei − ⟨wi , ei ⟩vi = ∑⟨wℓ , ei ⟩vℓ . ℓ≠i
Taking norms and applying the bound above gives ∥wℓ ∥∥vℓ ∥ ) δi . γi ℓ≠i
∥ei − ⟨wi , ei ⟩vi ∥ ≤ ∑ ∣⟨wℓ , ei ⟩∣ ∥vℓ ∥ ≤ (∑ ℓ≠i
□
Searching over the space of rank one matrices, or tensors, to find those that lie close to a subspace is behind tensor decomposition algorithms including [KP25, WPKS25]. In TSCD, this idea is adapted to testing the membership of specific rank one matrices. 6.3. Pairwise intervention-asymmetric correlation tests. In practice, roots can have scores less than one and non-root nodes may have high scores. See Appendix B.6 for empirical behavior of the scores at roots. To improve robustness, we retain a small set of top candidates with high αi scores. We then use a refinement step that resolves ambiguity among candidate root nodes using asymmetry in correlations across contexts. Consider two nodes i and j that are correlated in the observational context. If intervening on i leaves the correlation with j nonzero, then the dependence cannot be explained by paths ending at i; it must flow from i to j, so i is an ancestor of j. If instead the correlation disappears, the dependence must arise from a path into i (either j → i or a latent confounder), and thus i is not a root. This is formalized in the following lemma.
14
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
Algorithm 2: Greedy causal order recovery Input: Sample precision tensor T ∈ Rp×p×k ; intervention-pattern matrix B ∈ Rp×k ; number of candidates n. Output: A causal order π = (π1 , . . . , πp ). Initialize T , S ← {1, . . . , p}, and π ← ( ); for t = 1, . . . , p do Compute the projection operator PA onto the subspace spanned by the slices T (1, ∶, ∶), . . . , T (p, ∶, ∶) of the current tensor T ; foreach i ∈ S do ∥PA (ei ⊗ bi )∥F αi ← . ∥bi ∥ end Let Ct ⊆ S be the set of nodes with the top n values of αi ; πt ← RootSelection(Ct , Cov(Xj )}kj=1 , B); // Algorithm 3 Append πt to π; Remove πt from S; Delete the πt -th row and column from each slice T (∶, ∶, 1), . . . , T (∶, ∶, k) of T ; end return π; Lemma 6.3. Let i and j be two nodes in a DAG. Assume corr(xi , xj ) ≠ 0 in the observational context. If there is a context where i is intervened on and j is not and in which corr(xi , xj ) ≠ 0, then i is an ancestor of j. Otherwise, either j is an ancestor of i or i and j share an ancestor. i intervened, j not
i
j
dependence remains
intervening i removes incoming edge
i
j
dependence broken
confounding path into i removed
i
l
j
dependence broken
Figure 2. Intervention-asymmetric correlation test. If intervening on i preserves correlation with j, then the dependence must flow from i to j, so i is an ancestor of j. If the correlation disappears, the observational dependence may instead come from j → i or from a latent common cause. 6.4. Root selection algorithm. We use pairwise intervention-asymmetric tests to select a root from a candidate set C. For each pair i, j ∈ C, we compare how their correlation behaves across intervention contexts. Correlations are classified as NONZERO, ZERO, or INCONCLUSIVE using hypothesis tests (Appendix A.1). If the correlation between i and j is NONZERO when intervening on i but not when intervening on j, this is evidence that i is upstream of j. Conversely, ZERO correlation weakens this claim. We aggregate these pairwise comparisons into node scores and select the node with the highest score as the next root in the causal order. See Algorithm 3 for details.
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
15
Algorithm 3: RootSelection by pairwise correlation tests Input: Candidate set C; sample covariance matrices for contexts ℓ = 1, . . . , k; intervention-pattern matrix B; sample sizes {Nℓ }. Output: Selected root r ∈ C. Initialize score(i) ← 0 for all i ∈ C; foreach pair i ≠ j in C do for ℓ = 1, . . . , k do (ℓ) Compute the sample correlation ρij in each context.; (ℓ)
Classify ρij as NONZERO, ZERO, or INCONCLUSIVE for all ℓ.; end Let Cobs , Ci , Cj be the collections of contexts in which neither, only i, or only j are intervened, respectively; Initialize si ← 0 and sj ← 0; if Cobs is NONZERO then if NONZERO appears in Ci then (si , sj ) ← (si + 1, sj − 1); end else si ← si − 1; end if NONZERO appears in Cj then (si , sj ) ← (si − 1, sj + 1); end else sj ← sj − 1; end end if si > sj then score(i) ← score(i) + 1; score(j) ← score(j) − 1; end else if sj > si then score(j) ← score(j) + 1; score(i) ← score(i) − 1; end end r ← arg maxi∈C score(i); return r;
One may wonder whether pairwise correlation tests alone are sufficient for causal order. While such tests exploit intervention asymmetry, they rely only on local, pairwise information and can be unreliable in finite samples or when intervention patterns are limited, e.g., when there are few or no contexts in which one node is intervened on while the other is not. The
16
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
importance of combining both global projection norm-based candidate selection and local pairwise correlation testing steps is shown in Appendix B.7. 7. Experiments We evaluate TSCD on synthetic linear and nonlinear structural equation models, comparing against causal discovery baselines. Experiment details are in Appendix B. Experiments were conducted on a MacBook Pro with an Apple M2 chip and on a university computing cluster. See this repository for the code. 7.1. Linear SEM. We generate data from random LSEMs with p = 10 nodes and edge probability 0.6, under one observational and four intervention contexts. We consider three noise regimes (Gaussian, heavy-tailed, and mixed). We compare three classes of methods: (1) Interventional methods (ours, GIES, IGSP), which use interventional data and intervention patterns; (2) Observational methods with adjacency matrix output (LiNGAM, SortRegress, NOTEARS), applied separately to each context and combined; (3) Observational methods with graph output (GES, PC), applied to five times as many observational samples to match the total sample size used by the other methods. Performance is evaluated using the relative Frobenius error, ∥Λest − Λ∥F , ∥Λ∥F for methods estimating Λ, the F1 score for edge recovery, 2 precision ⋅ recall F1 = , precision + recall and runtime. Results are shown in Figure 3. Our method achieves high accuracy and computational efficiency, performing comparably to LiNGAM in non-Gaussian settings and outperforming baselines when there is Gaussian noise. 7.2. Scalability. We consider graphs with up to p = 400 nodes under sparse LSEMs with perfect interventions. Figure 4 reports F1 score and runtime as functions of p. The F1 score remains stable above 0.94 across the range, with no degradation as the problem size increases. Runtime grows from 0.13s at p = 50 to 214s at p = 400. The experiment was conducted on a MacBook Pro with an Apple M2 chip. 7.3. TSCD-nonlinear. We consider nonlinear SEMs with perfect interventions. The goal is to test whether we can recover the causal order and the DAG based on stability of each node when the structural functions are nonlinear. The study of nonlinear SEMs is widespread, also appearing in e.g. [ZDA+ 20, MZH20, HJM+ 08]. TSCD-nonlinear constructs the causal order greedily based on stability across contexts. The initial root is selected using marginal variance stability. At each step, we predict each candidate node from previously selected variables using a neural network trained on pooled non-intervened contexts, and select the node whose prediction error is most stable across contexts. Parent sets are estimated via a gated network for feature selection. We compare TSCD-nonlinear against a baseline that orders nodes by marginal variance. The experiment is conducted on a random DAG with p = 10 nodes and edge probability 0.8.
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
100 10−1 10
2
3
10 10 Sample size per context
1 0.5 0 102
4
10
3
10 10 Sample size per context
103 104 Sample size per context
102
10
10 10 Sample size per context
103 104 Sample size per context
Relative Frob. Error
F1 Score
4
Ours (TSCD) SortRegress
1
1
F1 Score
100
3
0.5
0
1.5 Relative Frob. Error
Runtime (sec)
F1 Score
0.5
Runtime (sec) Gaussian ratio = 1
F1 Score 1
1
102
4
2
2
103 104 Sample size per context
Relative Frob. Error
100
2
0.5
0 102
103 104 Sample size per context
1.5
2
Relative Frob. Error
Runtime (sec)
Gaussian ratio = 0.5
Runtime (sec)
10
1 F1 Score
101
10
F1 Score
Relative Frob. Error 1.5 Relative Frob. Error
10 Runtime (sec)
Gaussian ratio = 0
Runtime (sec) 2
17
0.5 0 102
103 104 Sample size per context
GES NOTEARS
GIES DirectLiNGAM
0.5
0 102
103 104 Sample size per context
IGSP ICA-LiNGAM
PC
Figure 3. Performance comparison across different noise settings, with Gaussian ratios 0, 0.5, and 1. Ordering accuracy is measured by the number of parent-child errors, edges where the parent appears after the child in the recovered order. Our method recovers the correct causal order, while the baseline produces 13 parent-child errors. A representative example is in Figure 5. Our recovered graph has an edge error rate of 0.34. The runtime is 134 seconds using 5 parallel workers. 7.4. Light Causal Chamber. We evaluate TSCD and TSCD-nonlinear on a real-world light-chamber experiment from the Causal Chamber [GPB25]. The chamber is a controlled physical system in which light sources and sensor readings interact through the device. We use the light-tunnel intervention dataset lt interventions standard v1. To focus on a setting where the intervention targets can be treated as perfect interventions, we restrict to a partial ground-truth graph over nine variables: R, G, B, I˜1 , Ṽ1 , L11 , L12 , D1V , T1V .
18
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL F1 score vs. number of nodes
Runtime vs. number of nodes 150
1.02
runtime (s)
F1 score
1 0.98 0.96
100
50
0.94 0 100
200 p
300
400
100
200 p
300
400
Figure 4. Scalability of TSCD. Left: F1 score versus number of nodes. Right: runtime versus number of nodes. Each point is averaged across 10 random sparse DAGs, each with sample size p2 . The error bars are given by the standard deviation across trials.
Figure 5. Ground-truth DAG (left) and recovered graph (right) for a nonlinear SEM. The recovered structure matches the true causal order, with an edge error rate of 0.34.
Here R, G, B denote the brightness of the red, green, and blue LEDs respectively on the main light source. The variables I˜1 and Ṽ1 are the uncalibrated infrared and visible-light intensity measurements from the first light sensor, which is placed before the polarizers. The variables L11 and L12 are the brightness settings of the two auxiliary LEDs placed by this first light sensor. Finally, D1V and T1V are visible-light sensor parameters for the first sensor: the selected visible photodiode and the photodiode exposure time. We compare against several standard causal discovery baselines. TSCD-nonlinear achieves high accuracy, recovering all ground-truth edges in the selected graph with precision 0.80 and recall 1.00. TSCD recovers the linear part of the graph, the edges from R, G, B to the two
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
Method
19
Precision Recall
TSCD 0.857 0.500 TSCD-nonlinear 0.800 1.000 SortRegress 0.038 0.500 LiNGAM (ICA) 0.028 0.500 LiNGAM (Direct) 0.019 0.500 GIES 0.833 0.417 IGSP 0.000 0.167 Table 1. Causal discovery performance on the selected Causal Chamber light-tunnel graph. sink variables I˜1 and Ṽ1 , yielding precision 0.857 and recall 0.500. In contrast, the competing methods have lower recall or precision on this benchmark, see Table 1. 7.5. Flow-cytometry Dataset. We consider the flow-cytometry dataset [SPP+ 05], which records single-cell protein phosphorylation for p = 11 proteins across nine experimental conditions. We pool the two unperturbed baselines into one observational environment, and use the five kinase-inhibitor conditions as intervention contexts. We exclude the two activator conditions (PMA and β2-cAMP), since inhibitor interventions are closer to perfect interventions. We compare the graph recovered by our method in Table 2 with two reference DAGs from [SPP+ 05] and graphs recovered by ICP and hiddenICP [PBM16, MHM+ 16]. TSCD recovers eight edges, seven of which appear in the DAG from [SPP+ 05] or recovered by hiddenICP. Most missing edges are either among {Raf, PLCg, Erk, PKA, p38, JNK}, which are not intervened in any context so are not identifiable, or involve PKA, PKC, where feedback loops and nonlinear effects are expected [IOS+ 10]. 8. Conclusion We propose a tensor-based method for causal discovery from covariance matrices, across multiple intervention contexts, to recover both a causal graph and edge weights. Our approach is identifiable under a logarithmic number of intervention contexts and performs well in linear and nonlinear settings, scaling to graphs with hundreds of nodes. Our method relies on assumptions such as perfect interventions and uncorrelated noise, which may be violated in real-world data. Future directions include extending the framework to handle other types of interventions and latent confounding. These could include noiseshift interventions [RHPM15], which also admit a tensor decomposition structure, and soft interventions [MGS05] which modify rather than remove edge weights. TSCD relies only on second-order statistics through covariance and precision matrices. Incorporating first or higher-order statistics is another interesting direction for future work. Acknowledgments. NO received support from the Harvard Office of Undergraduate Research and Fellowships. AS was partially supported by NSF DMS 2608217 and an Alfred P. Sloan research fellowship.
References [BP02]
Carlos Brito and Judea Pearl. A new identification condition for recursive models with correlated errors. Structural Equation Modeling, 9(4):459–474, 2002.
20
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
Edge
[SPP+ 05]a [SPP+ 05]b [PBM16]ICP [PBM16]hiddenICP TSCD (ours)
Raf → Mek Mek → Raf Mek → Erk PLCg → PIP2 PLCg → PIP3 PLCg → PKC PIP2 → PLCg PIP2 → PKC PIP3 → PLCg PIP3 → PIP2 PIP3 → Akt Akt → Erk Akt → PKA Erk → Akt PKA → Raf PKA → Mek PKA → Erk PKA → Akt PKA → p38 PKA → JNK PKC → Raf PKC → Mek PKC → PKA PKC → p38 PKC → JNK p38 → JNK p38 → PKC JNK → PKC JNK → p38 Total
✓
✓
✓ ✓
✓ ✓ ✓
✓ ✓ ✓
✓
✓
✓ ✓
✓
✓
✓
✓
✓
✓ ✓
✓
✓
✓ ✓ ✓ ✓ ✓ ✓
✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓
18
✓
✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓
17
✓ ✓
✓
7
✓ ✓ ✓ ✓ ✓ ✓
✓ ✓
13
8
Table 2. Edges comparison across two reference graphs, graph recovered by ICP, hiddenICP and TSCD. TSCD recovers a sparse set of edges that is consistent with the other graphs. Red entries are non-identifiable edges, between {Raf, PLCg, Erk, PKA, p38, JNK}, which TSCD cannot identify.
[Chi02] [CO12] [DFS11] [DL06]
[DPL+ 16]
David M. Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3:507–554, 2002. Luca Chiantini and Giorgio Ottaviani. On generic identifiability of 3-tensors of small rank. SIAM Journal on Matrix Analysis and Applications, 33(3):1018–1037, 2012. Mathias Drton, Rina Foygel, and Seth Sullivant. Global identifiability of linear structural equation models. The Annals of Statistics, 39(2):865–886, 2011. Lieven De Lathauwer. A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization. SIAM journal on Matrix Analysis and Applications, 28(3):642–666, 2006. Atray Dixit, Oren Parnas, Biyu Li, Jenny Chen, Charles P Fulco, Livnat Jerby-Arnon, Nemanja D Marjanovic, Danielle Dionne, Tyler Burks, Raktima Raychowdhury, et al. Perturb-seq:
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
21
dissecting molecular circuits with scalable single-cell RNA profiling of pooled genetic screens. cell, 167(7):1853–1866, 2016. [EGS05] Frederick Eberhardt, Clark Glymour, and Richard Scheines. On the number of experiments sufficient and in the worst case necessary to identify all causal relations among n variables. In Proceedings of the Twenty-First Conference on Uncertainty in Artificial Intelligence, pages 178–184, 2005. [FDD12] Rina Foygel, Jan Draisma, and Mathias Drton. Half-trek criterion for generic identifiability of linear structural equation models. The Annals of Statistics, pages 1682–1713, 2012. [Fis15] Ronald A Fisher. Frequency distribution of the values of the correlation coefficient in samples from an indefinitely large population. Biometrika, 10(4):507–521, 1915. [FLNP00] Nir Friedman, Michal Linial, Iftach Nachman, and Dana Pe’er. Using Bayesian networks to analyze expression data. In Proceedings of the fourth annual international conference on Computational molecular biology, pages 127–135, 2000. [GPB25] Juan L Gamella, Jonas Peters, and Peter Bühlmann. Causal chambers as a real-world physical testbed for AI methodology. Nature Machine Intelligence, 7(1):107–118, 2025. [GPR99] Sander Greenland, Judea Pearl, and James M Robins. Causal diagrams for epidemiologic research. Epidemiology, 10(1):37–48, 1999. [HB12] Alain Hauser and Peter Bühlmann. Characterization and greedy learning of interventional Markov equivalence classes of directed acyclic graphs. Journal of Machine Learning Research, 13:2409–2464, 2012. [HCA+ 25] Ambroise Heurtebise, Omar Chehab, Pierre Ablin, Alexandre Gramfort, and Aapo Hyvärinen. Identifiable multi-view causal discovery without non-gaussianity. arXiv e-prints, pages arXiv– 2502, 2025. [HJM+ 08] Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. Advances in neural information processing systems, 21, 2008. [Imb20] Guido W Imbens. Potential outcome and directed acyclic graph approaches to causality: Relevance for empirical practice in economics. Journal of Economic Literature, 58(4):1129–1179, 2020. [IOS+ 10] Sleiman Itani, Mesrob Ohannessian, Karen Sachs, Garry P Nolan, and Munther A Dahleh. Structure learning in causal cyclic networks. In Causality: Objectives and assessment, pages 165–176. PMLR, 2010. [JLV23] Nathaniel Johnston, Benjamin Lovitz, and Aravindan Vijayaraghavan. Computing linear sections of varieties: quantum entanglement, tensor decompositions and beyond. In 2023 IEEE 64th Annual Symposium on Foundations of Computer Science (FOCS), pages 1316–1336. IEEE, 2023. [KB09] Tamara G Kolda and Brett W Bader. Tensor decompositions and applications. SIAM review, 51(3):455–500, 2009. [KP25] Joe Kileel and João M Pereira. Subspace power method for symmetric tensor decomposition. Numerical Algorithms, pages 1–38, 2025. [Kru77] Joseph B Kruskal. Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics. Linear algebra and its applications, 18(2):95– 138, 1977. [Lan11] Joseph M Landsberg. Tensors: geometry and applications: geometry and applications, volume 128. American Mathematical Soc., 2011. [MGS05] Florian Markowetz, Steffen Grossmann, and Rainer Spang. Probabilistic soft interventions in conditional Gaussian networks. In International Workshop on Artificial Intelligence and Statistics, pages 214–221. PMLR, 2005. [MHM+ 16] Nicolai Meinshausen, Alain Hauser, Joris M Mooij, Jonas Peters, Philip Versteeg, and Peter Bühlmann. Methods for causal inference from gene perturbation experiments and validation. Proceedings of the National Academy of Sciences, 113(27):7361–7368, 2016. [MZH20] Ricardo Pio Monti, Kun Zhang, and Aapo Hyvärinen. Causal discovery with general non-linear relationships using non-linear ica. In Uncertainty in artificial intelligence, pages 186–195. PMLR, 2020.
22
[PB14]
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
Jonas Peters and Peter Bühlmann. Identifiability of gaussian structural equation models with equal error variances. Biometrika, 101(1):219–228, 2014. [PBM16] Jonas Peters, Peter Bühlmann, and Nicolai Meinshausen. Causal inference using invariant prediction: Identification and confidence intervals. Journal of the Royal Statistical Society: Series B, 78(5):947–1012, 2016. [Pea09] Judea Pearl. Causality. Cambridge University Press, 2 edition, 2009. [PJS17] Jonas Peters, Dominik Janzing, and Bernhard Scholkopf. Elements of causal inference: foundations and learning algorithms. MIT press, 2017. [RHPM15] Dominik Rothenhäusler, Christina Heinze, Jonas Peters, and Nicolai Meinshausen. Backshift: Learning causal cyclic graphs from unknown shift interventions. Advances in neural information processing systems, 28, 2015. [RSW21] Alexander Reisach, Christof Seiler, and Sebastian Weichwald. Beware of the simulated DAG! causal discovery benchmarks may be easy to game. Advances in Neural Information Processing Systems, 34:27772–27784, 2021. [RSW26] Kristian Ranestad, Anna Seigal, and Kexin Wang. A real generalized trisecant trichotomy. Journal of Algebra, 2026. [SB23] Christoph Schultheiss and Peter Bühlmann. Ancestor regression in linear structural equation models. Biometrika, 110(4):1117–1124, 2023. [SGS00] Peter Spirtes, Clark N Glymour, and Richard Scheines. Causation, prediction, and search. MIT press, 2000. [SHHK06] Shohei Shimizu, Patrik O. Hoyer, Aapo Hyvärinen, and Antti Kerminen. A linear non-Gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7:2003–2030, 2006. [SIS+ 11] Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvärinen, Yoshinobu Kawahara, Takashi Washio, Patrik O. Hoyer, and Kenneth Bollen. Directlingam: A direct method for learning a linear non-Gaussian structural equation model. Journal of Machine Learning Research, 12:1225–1248, 2011. [SKPF22] Shan H Siddiqi, Konrad P Kording, Josef Parvizi, and Michael D Fox. Causal mapping of human brain function. Nature reviews neuroscience, 23(6):361–375, 2022. [SLB+ 21] Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio. Toward causal representation learning. Proceedings of the IEEE, 109(5):612–634, 2021. [SPP+ 05] Karen Sachs, Omar Perez, Dana Pe’er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308(5721):523– 529, 2005. [SSBU23] Chandler Squires, Anna Seigal, Salil S Bhate, and Caroline Uhler. Linear causal disentanglement via interventions. In International conference on machine learning, pages 32540–32560. PMLR, 2023. [Stu08] Student. The probable error of a mean. Biometrika, pages 1–25, 1908. [Sul23] Seth Sullivant. Algebraic statistics, volume 194. American Mathematical Society, 2023. [SVBDL15] Laurent Sorber, Marc Van Barel, and Lieven De Lathauwer. Structured data fusion. IEEE journal of selected topics in signal processing, 9(4):586–600, 2015. [VP90] Thomas Verma and Judea Pearl. Equivalence and synthesis of causal models. In Proceedings of the Sixth Annual Conference on Uncertainty in Artificial Intelligence, pages 255–270, 1990. [WBP+ 26] Kexin Wang, Salil Bhate, João M Pereira, Joe Kileel, Matylda Figlerowicz, and Anna Seigal. Multi-context principal component analysis. arXiv preprint arXiv:2601.15239, 2026. [Wed34] Joseph Henry Maclagan Wedderburn. Lectures on matrices, volume 17. American Mathematical Soc., 1934. [WPKS25] Kexin Wang, João M Pereira, Joe Kileel, and Anna Seigal. Multi-subspace power method for decomposing all tensors. arXiv preprint arXiv:2510.18627, 2025. [Wri34] Sewall Wright. The method of path coefficients. The annals of mathematical statistics, 5(3):161– 215, 1934. [WS23] Kexin Wang and Anna Seigal. Lower bounds on the rank and symmetric rank of real tensors. Journal of Symbolic Computation, 118:69–92, 2023.
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
23
[WSYU17] Yuhao Wang, Liam Solus, Karren Yang, and Caroline Uhler. Permutation-based causal inference algorithms with interventions. Advances in neural information processing systems, 30, 2017. [ZARX18] Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. DAGs with NO TEARS: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems, 2018. [ZDA+ 20] Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. Learning sparse nonparametric dags. In International conference on artificial intelligence and statistics, pages 3414–3425. Pmlr, 2020.
Appendix A. Algorithm Details A.1. Hypothesis testing for zero correlation. The pairwise intervention-asymmetric tests whether a sample correlation is effectively zero. Finite sample can make small correlations appear nonzero (and vice versa), so we use a three-way classification: NONZERO,
ZERO,
INCONCLUSIVE.
This allows us to avoid overconfident conclusions when the data are ambiguous. Testing for nonzero correlation. Given a sample correlation ρ̂ from N samples, we test H0 ∶ ρ = 0
against
We use the standard t-statistic
√ t = ρ̂
H1 ∶ ρ ≠ 0.
N −2 , 1 − ρ̂2
which follows a t-distribution with N − 2 degrees of freedom under H0 [Stu08]. If the two-sided p-value is less than a significance level α, we classify the correlation as NONZERO. Testing for effective zero correlation. If the null is not rejected, we next test the correlation using an equivalence test. Fix a tolerance ϵ > 0, and test H0 ∶ ∣ρ∣ ≥ ϵ
against
H1 ∶ ∣ρ∣ < ϵ.
We implement this using two one-sided tests. Using Fisher’s z-transform [Fis15], z = arctanh(ρ̂),
z± = arctanh(±ϵ),
with standard error 1 SE = √ , N −3 we compute
z − z− z − z+ , Z+ = . SE SE If Z− > z1−α and Z+ < −z1−α , we reject H0 and classify the correlation as ZERO. Inconclusive cases. If neither test is decisive, we classify the correlation as INCONCLUSIVE. This occurs when the data cannot distinguish between zero and nonzero correlation. Summary. We classify each sample correlation as follows: ● NONZERO: reject H0 ∶ ρ = 0, ● ZERO: accept equivalence ∣ρ∣ < ϵ, ● INCONCLUSIVE: neither test is decisive. This classification is used in the pairwise intervention-asymmetric test to infer causal order. Z− =
24
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
A.2. From causal order to adjacency matrix. Given a causal order π = (π1 , . . . , πp ), recovering the adjacency matrix Λ reduces to estimating the coefficients of each node from its ancestors in the order. For each node πt , its parents lie in {π1 , . . . , πt−1 }, so estimation can be restricted to this set. We consider the following approaches. Lasso regression. For each node πt , we regress xπt onto its ancestors {xπ1 , . . . , xπt−1 } using Lasso. The ℓ1 -penalty encourages sparse solutions. After identifying the support, one can optionally refit an ordinary least squares regression restricted to the parents to obtain less biased estimates. The resulting coefficients define the corresponding row of Λ. Regression. For each πt , we regress xπt onto its ancestors. The coefficients define the corresponding row of Λ. We set entries with magnitude below a threshold to zero. Cholesky decomposition. The causal order induces a permutation under which the matrix (I − Λ) is lower triangular. Applying a Cholesky decomposition to the permuted precision matrix yields a dense estimate of the structural coefficients. To obtain a sparse adjacency matrix, we set entries with magnitude below a threshold to zero. The Lasso approach requires access to individual samples; the other two methods do not. Appendix B. Experiment details B.1. Linear SEM. For each experiment instance, we generate a random DAG with p = 10 nodes. We first sample a weighted adjacency matrix Λ ∈ Rp×p , where each edge is included independently with probability 0.6. Nonzero edge weights are drawn uniformly from [0.4, 1] with random signs. To enforce acyclicity, we construct Λ to be lower-triangular and then apply a random permutation. Noise variables ϵi are independent across nodes, with variances sampled i.i.d. from U[0.1, 1]. We consider three noise settings: (i) fully Gaussian noise, (ii) half of the nodes have Gaussian noise and half follow a Student-t5 distribution, and (iii) fully Student-t5 noise. We use the intervention-pattern matrix ⎛1 ⎜1 ⎜ ⎜1 ⎜ ⎜1 ⎜ ⎜1 ⎜ B=⎜ ⎜1 ⎜ ⎜1 ⎜ ⎜1 ⎜ ⎜1 ⎜ ⎝1
1 1 1 1 1 1 1 1 0 0
1 1 1 1 0 0 0 0 1 1
1 1 0 0 1 1 0 0 1 1
1⎞ 0⎟ ⎟ 1⎟ ⎟ 0⎟ ⎟ 1⎟ ⎟ ⎟ ∈ R10×5 . 0⎟ ⎟ 1⎟ ⎟ 0⎟ ⎟ 1⎟ ⎟ 0⎠
Each row gives the intervention-pattern vector across the 5 contexts. For each context j, we generate samples from the corresponding interventional LSEM. Under perfect interventions, all incoming edges to intervened nodes are removed and their structural equations are replaced with exogenous noise. Samples are generated by solving the linear system (I − Λ(j) )X (j) = E (j) , where Λ(j) is the modified adjacency matrix for context j, and E (j) contains independent noise samples. Each context contains the same number of samples, ranging from 100 to 10,000. Each experimental configuration is repeated over 30 independent trials. The root is
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
25
selected using a candidate set of size 3 (Algorithm 3). From the causal order, Λ is estimated using Lasso regression. B.2. Scalability. We vary the number of nodes p ∈ {50, 60, . . . , 400}. For each p, we generate a random DAG with edge probability 6/(p − 1), yielding an expected in-degree of around 3. Edge weights are drawn from Unif([−1, −0.1] ∪ [0.1, 1]), and noise variances from Unif(0.01, 0.05), with Gaussian noise. We use a binary-code design with L = ⌈log2 p⌉ + 1 contexts, including one observational context. Each context contains p2 samples. The root is selected using a candidate set of size 3 (Algorithm 3). The experimental configuration was repeated over 10 independent trials. Given the recovered causal order, we estimate the adjacency matrix via regression in the observational context with threshold 0.1. Performance is evaluated using the F1 score on the recovered adjacency support. B.3. Nonlinear SEM. We consider nonlinear SEMs with perfect interventions on random DAGs with p = 10 nodes and edge probability 0.8. Data are generated by xi = fi (xpa(i) ) + ϵi , where fi is a randomly generated nonlinear function and ϵi is Gaussian noise. Each fi is implemented as a one-hidden-layer neural network with randomly initialized weights and a combination of activation functions, including tanh, ReLU, and square activations. Under interventions, if a node appears as a parent in a downstream structural equation, its input is masked to zero in that context, removing its influence. The structural equation of an intervened node is replaced by an exogenous noise term. We construct the causal order greedily. The initial root is selected based on stability of marginal variances across contexts where the node is not intervened on. At each subsequent step, for each candidate node, we train a feedforward neural network to predict the node from previously selected variables, using data pooled across contexts where the node is not intervened on. If a parent variable is intervened on in a context, its input is set to zero. To estimate the parent set, we use a gated neural network, where gates act as feature-selection weights. Variables with gate logits above a threshold are retained as parents. The model is then optionally refit using only the selected parents; if this refit improves validation loss, the reduced parent set is accepted. Otherwise, the full previously selected set is kept. Node selection is based on stability of prediction error across contexts. For each candidate node, we compute the mean squared error (MSE) in each valid context and select the node with the most stable normalized MSE, measured by the variance across contexts with a small penalty on the mean. B.4. Flow-cytometry Dataset. We evaluate our method on the flow-cytometry dataset [SPP+ 05], which records single-cell protein phosphorylation. for p = 11 signaling proteins across nine experimental conditions. We (i) pool the two unperturbed baselines (antiCD3/CD28 with and without ICAM-2 co-stimulation) into a single observational environment, and (ii) restrict the intervention contexts to the five kinase-inhibitor conditions, each targeting a single protein (Akt, PKC, PIP2, Mek, PIP3 via PI3K). We exclude the two pharmacological-activator conditions (PMA and β2-cAMP), because inhibitors more closely approximate a perfect intervention while for activators, the target remains conditionally dependent on its parents. This yields six environments and 5,846 cells in total. We compare the graph recovered by our method (Table 2) with two reference DAGs from [SPP+ 05] and the recovered graphs of ICP and hiddenICP [PBM16, MHM+ 16].
26
NATHAN OUYANG, KEXIN WANG AND ANNA SEIGAL
B.5. Light Causal Chamber. We use the Causal Chamber light-tunnel dataset with single node interventions lt interventions standard v1 from [GPB25]. The selected variables and dataset columns are {red, green, blue, ir 1, vis 1, l 11, l 12, diode vis 1, t vis 1}, corresponding to R, G, B, I˜1 , Ṽ1 , L11 , L12 , D1V , T1V . The variables R, G, B control the main red, green, and blue LEDs, while L11 , L12 control the two auxiliary LEDs placed by the first lightintensity sensor. The variables D1V and T1V are visible-light sensor parameters for the first sensor, controlling the selected visible photodiode and the photodiode exposure time. The variables I˜1 and Ṽ1 are sensor readings: the infrared and visible-light intensity measurements from the first sensor. The selected ground-truth graph contains twelve directed edges: R, G, B, L11 , L12 → I˜1 , R, G, B, L11 , L12 , D1V , T1V → Ṽ1 . We use the reference experiment together with interventions on the selected actuator and sensor-parameter variables: uniform reference, uniform red strong, uniform green strong, uniform blue strong, uniform t vis 1 strong, uniform l 11 mid, uniform l 12 mid, uniform diode vis 1 mid. TSCD recovers the six edges from R, G, B to I˜1 and Ṽ1 , corresponding to the strongest approximately linear effects in the subsystem. For TSCD-nonlinear, we set gate threshold to 2. It recovers all twelve ground-truth edges and adds three false positives, giving precision 12/15 = 0.800 and recall 12/12 = 1.000. For the other methods, we add a small independent Gaussian noise 10−5 ϵ, with ϵ ∼ N (0, 1), to each data entry before running the algorithms. This is necessary because some variables have zero empirical variance in certain intervention context, which causes numerical failures for several baselines. The perturbation is negligible relative to the scale of the measurements and is used only to help the other methods. B.6. Root Projection Scores. We study the projection score of roots, specifically, where the true roots appear in the ranking of nodes by projection score. In the population setting, the roots would all have projection score equal to 1 and appear at the top, but in practice they can be lower down the list. How far down they may appear guides the number of candidates that should be used for candidate root selection in Section 6.2. The top scoring root has a better rank when the graph is more sparse, since there are more roots. B.7. Ablation. We perform an ablation study using the LSEM data generation setup and hyperparameters in Section B.1. At each step, TSCD selects a candidate set of possible root nodes, then tests pairwise correlations to make a choice of root. We compare it to two ablated versions. The first only uses projection norm as a score to select the root at each step, i.e. there is only one root candidate. The second only uses pairwise correlation testing to determine the root at each step, i.e. all nodes are root candidates. Figure 7 reports these results. TSCD has improved performance compared to the ablated algorithms at low to moderate sample sizes. This justifies the importance of including both steps: selecting root candidates using projection norm and using pairwise correlations to choose the root from the candidates.
TENSOR-BASED SECOND-ORDER CAUSAL DISCOVERY
Rank
(a) e = 0.6
(b) e = 0.8
5
5
10
10
15
15 5
10
15
20
27
25
30
5
10
15
p
20
25
30
p
Figure 6. The distribution of best ranks for root nodes in projection score, as p (the number of nodes) varies, for two values of e (the edge probability), cf. Section B.1. All experiments were conducted with 1000 samples, half of the nodes given Gaussian noise and the rest given Student-t5 distribution noise.
103 104 Sample size per context
0.4 0.2
102
0.7 103 104 Sample size per context
102
103 104 Sample size per context
103 104 Sample size per context
F1 Score 1
0.6 F1 Score
Relative Frob. Error
Runtime (sec)
0.4 0.2
0.9 0.8 0.7
0 102
Runtime (sec)
102
103 104 Sample size per context
103 104 Sample size per context
F1 Score
Runtime (sec)
Relative Frob. Error
Relative Frob. Error
10−1 103 104 Sample size per context
1 0.6 F1 Score
Gaussian ratio = 0.5
0.8
Relative Frob. Error
10−1
102
0.9
0
Runtime (sec)
102
1
0.6 F1 Score
Relative Frob. Error
10−1 102
Gaussian ratio = 1
F1 Score
Relative Frob. Error 0.8
Runtime (sec)
Gaussian ratio = 0
Runtime (sec)
0.4 0.2
0.9 0.8 0.7
0 102
Ours (TSCD)
103 104 Sample size per context
102
Pairwise correlations only
Proj. norm only
103 104 Sample size per context
Figure 7. TSCD performance comparison with ablated versions across different noise settings, with Gaussian ratios 0, 0.5, and 1. Gaussian ratio denotes the ratio of nodes given Gaussian distributed noise, with the rest given Student-t5 distributions.