Score-Based Causal Discovery of Latent Variable Causal Models

# Score-Based Causal Discovery of Latent Variable Causal Models
Source: [https://arxiv.org/html/2605.20396](https://arxiv.org/html/2605.20396)
###### Abstract

Identifying latent variables and the causal structure involving them is essential across various scientific fields\. While many existing works fall under the category of constraint\-based methods \(with e\.g\. conditional independence or rank deficiency tests\), they may face empirical challenges such as testing\-order dependency, error propagation, and choosing an appropriate significance level\. These issues can potentially be mitigated by properly designed score\-based methods, such as Greedy Equivalence Search \(GES\)\(Chickering,[2002](https://arxiv.org/html/2605.20396#bib.bib19)\)in the specific setting without latent variables\. Yet, formulating score\-based methods with latent variables is highly challenging\. In this work, we develop score\-based methods that are capable of identifying causal structures containing causally\-related latent variables with identifiability guarantees\. Specifically, we show that a properly formulated scoring function can achieve score equivalence and consistency for structure learning of latent variable causal models\. We further provide a characterization of the degrees of freedom for the marginal over the observed variables under multiple structural assumptions considered in the literature, and accordingly develop both exact and continuous score\-based methods\. This offers a unified view of several existing constraint\-based methods with different structural assumptions\. Experimental results validate the effectiveness of the proposed methods\.

Machine Learning, ICML

## 1Introduction

At the core of understanding complex systems lies causal discovery, the identification of causal relations from observational data\(Spirtes et al\.,[2001](https://arxiv.org/html/2605.20396#bib.bib72); Pearl,[2009](https://arxiv.org/html/2605.20396#bib.bib54)\)\. One common assumption in causal discovery algorithms is the absence of latent confounders, known ascausal sufficiency, positing that the observed correlations stem either from true causation or can be sufficiently explained by other observed variables\. Yet, real\-world scenarios often defy this assumption\. For instance, in psychological studies, the measured questionnaires are indirect proxies of latent mental factors\. In unstructured data like images and texts, the observed pixels and words are confounded by latent semantic variables\. Directly applying causal discovery methods without considering these latent variables can lead to false discoveries, as latent variables may introduce spurious correlations among observed ones that cannot be attributed to true causation\.

Notable efforts have thus been made to identify the true causal relations in the presence of latent variables\. Earliest attempts include Fast Causal Inference \(FCI\)\(Spirtes et al\.,[2001](https://arxiv.org/html/2605.20396#bib.bib72); Zhang,[2008](https://arxiv.org/html/2605.20396#bib.bib86)\)and its variants\(Colombo et al\.,[2012](https://arxiv.org/html/2605.20396#bib.bib23); Spirtes et al\.,[2013](https://arxiv.org/html/2605.20396#bib.bib73); Claassen et al\.,[2013](https://arxiv.org/html/2605.20396#bib.bib22); Akbari et al\.,[2021](https://arxiv.org/html/2605.20396#bib.bib2)\)that exploit conditional independence information\. There are two main limitations of FCI: First, the results, presented by partial ancestral graphs \(PAG\)\(Richardson,[1996](https://arxiv.org/html/2605.20396#bib.bib58)\), tend to be overgeneralized – e\.g\., whenever two observed variables may be confounded, it indicates so\. Second, it focuses solely on causal relations among observed variables and does not provide information about those among latent variables\. In short, FCI does not require specific assumptions about the latent structure, at the cost of having a less informative output\. In contrast, one may often be interested in identifying the causal relations among latent variables \(e\.g\., the latent mental and semantic variables in the above examples\)\.

Hence, another line of work has been developed to discover the causal structure also among latent variables\. For the identifiability conditions, these methods typically introduce additional parametric assumptions to mitigate the large model indeterminacies faced by FCI\. This includes rank or tetrad condition\-based methods with linearity assumption\(Silva et al\.,[2003](https://arxiv.org/html/2605.20396#bib.bib67),[2006](https://arxiv.org/html/2605.20396#bib.bib68); Silva & Scheines,[2005](https://arxiv.org/html/2605.20396#bib.bib66); Choi et al\.,[2011](https://arxiv.org/html/2605.20396#bib.bib20); Kummerfeld & Ramsey,[2016](https://arxiv.org/html/2605.20396#bib.bib44); Huang et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib38); Dong et al\.,[2023](https://arxiv.org/html/2605.20396#bib.bib28)\), high\-order moments\-based methods\(Shimizu et al\.,[2009](https://arxiv.org/html/2605.20396#bib.bib64); Zhang et al\.,[2018](https://arxiv.org/html/2605.20396#bib.bib87); Cai et al\.,[2019](https://arxiv.org/html/2605.20396#bib.bib17); Salehkaleybar et al\.,[2020](https://arxiv.org/html/2605.20396#bib.bib59); Xie et al\.,[2020](https://arxiv.org/html/2605.20396#bib.bib83); Adams et al\.,[2021](https://arxiv.org/html/2605.20396#bib.bib1); Dai et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib26); Chen et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib18); Améndola et al\.,[2023](https://arxiv.org/html/2605.20396#bib.bib4); Wang & Drton,[2023](https://arxiv.org/html/2605.20396#bib.bib82)\), matrix decomposition\-based methods\(Anandkumar et al\.,[2013](https://arxiv.org/html/2605.20396#bib.bib5)\), copula model\-based methods\(Cui et al\.,[2018](https://arxiv.org/html/2605.20396#bib.bib25)\), mixture oracles\-based methods\(Kivva et al\.,[2021](https://arxiv.org/html/2605.20396#bib.bib42)\), and multiple domains\-based methods\(Zeng et al\.,[2021](https://arxiv.org/html/2605.20396#bib.bib85); Sturma et al\.,[2023](https://arxiv.org/html/2605.20396#bib.bib75)\)\. For the algorithmic procedures, these methods generally fall under the category of constraint\-based methods, by matching the statistical properties to possible structural patterns and constructing the whole causal structure iteratively\. A typical constraint\-based method in the causally sufficient case is PC\(Spirtes & Glymour,[1991](https://arxiv.org/html/2605.20396#bib.bib71)\)\. Despite the asymptotic consistency, the empirical reliability of constraint\-based methods may be limited due totesting\-order dependencyanderror propagation\(Spirtes,[2010](https://arxiv.org/html/2605.20396#bib.bib70); Colombo et al\.,[2012](https://arxiv.org/html/2605.20396#bib.bib23)\), especially when the number of variables is large\.

To address such empirical issues of constraint\-based methods, score\-based causal discovery methods have been introduced, and may be more favored in practical applications\(Nandy et al\.,[2018](https://arxiv.org/html/2605.20396#bib.bib47); Ramsey et al\.,[2017](https://arxiv.org/html/2605.20396#bib.bib55)\)\. Unlike the iterative construction of a single causal graph by constraint\-based methods, score\-based methods assign a score to each potential graph reflecting how well it explains the observed data and generally search over the graph space to find the optimal graph\. In the causally sufficient case, one typical score\-based method is the Greedy Equivalence Search \(GES\)\(Chickering,[2002](https://arxiv.org/html/2605.20396#bib.bib19)\)\. There also exists several score\-based methods that can handle latent variables\(Shpitser et al\.,[2012](https://arxiv.org/html/2605.20396#bib.bib65); Triantafillou & Tsamardinos,[2016](https://arxiv.org/html/2605.20396#bib.bib77); Nowzohour et al\.,[2017](https://arxiv.org/html/2605.20396#bib.bib52); Bhattacharya et al\.,[2021](https://arxiv.org/html/2605.20396#bib.bib11); Shahin & Chechik,[2020](https://arxiv.org/html/2605.20396#bib.bib63); Bernstein et al\.,[2020](https://arxiv.org/html/2605.20396#bib.bib8); Bellot & van der Schaar,[2021](https://arxiv.org/html/2605.20396#bib.bib6); Claassen & Bucur,[2022](https://arxiv.org/html/2605.20396#bib.bib21)\)\. Similar to FCI, most of them do not discover the causal relations among latent variables, except the method byZhang \([2004](https://arxiv.org/html/2605.20396#bib.bib89)\)without identifiability guarantee\. When latent variables are introduced and relations among them are further allowed in the causal structures, challenges arise in characterizing the degrees of freedom\(Geiger et al\.,[1996](https://arxiv.org/html/2605.20396#bib.bib32),[2001](https://arxiv.org/html/2605.20396#bib.bib33)\), formulating a scoring function, and structuring the search procedure\. We tackle these challenges in this paper, and to the best of our knowledge, this is the first score\-based method that identifies causal structures containing causally\-related latent variables with identifiability guarantees\.

Contributions\.We develop score\-based methods, called SALAD \(which stands for Score\-bAsed Latent cAusal Discovery\), for causal discovery of latent variable causal models, providing a unified view for several existing constraint\-based methods\(Silva et al\.,[2003](https://arxiv.org/html/2605.20396#bib.bib67); Huang et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib38)\)\. Our contributions can be summarized as follows:

- •We develop a formulation of scoring function for identifying linear latent variable causal models\. We show \(1\) that it is score equivalent and \(2\) that minimizing it yields a structure that is algebraic equivalent to the true structure\. The latter implies that both structures have the same equality constraints \(on the marginal over the observed variables\), including conditional independence and rank deficiency constraints\.
- •We provide a characterization of the degrees of freedom for the marginal over the observed variables under the structural assumptions considered bySilva et al\. \([2003](https://arxiv.org/html/2605.20396#bib.bib67)\); Huang et al\. \([2022](https://arxiv.org/html/2605.20396#bib.bib38)\)\.
- •We develop exact score\-based methods for estimating the causal structure, and show that they can asymptotically identify the true equivalence class of the whole structure\. We also provide continuous score\-based methods in some of the settings to improve the computational efficiency\.
- •We demonstrate that the proposed score\-based methods achieve improved performance over existing constraint\-based methods for estimating the structures of latent variable causal models, which further validate the effectiveness of score\-based methods\.

Notations\.For a matrixMM, we define its support set assupp⁡\(M\)≔\{\(i,j\):Mi,j≠0\}\\operatorname\{supp\}\(M\)\\coloneqq\\\{\(i,j\):M\_\{i,j\}\\neq 0\\\}\. We denote byM𝐒,:M\_\{\\mathbf\{S\},:\}the rows inMMindexed by set𝐒\\mathbf\{S\}, and similarly byM:,𝐒M\_\{:,\\mathbf\{S\}\}for the columns\. For a directed acyclic graph \(DAG\)𝒢\\mathcal\{G\}, we denote by\|𝒢\|\|\\mathcal\{G\}\|the number of edges in𝒢\\mathcal\{G\}\. Also, letdiag⁡\(ℝ\>0m\)\\operatorname\{diag\}\(\\mathbb\{R\}\_\{\>0\}^\{m\}\)be the set ofm×mm\\times mdiagonal matrices with positive diagonal entries,𝕌m\\mathbb\{U\}^\{m\}be the set ofm×mm\\times mstrictly upper triangular matrices, and𝔾m\\mathbb\{G\}^\{m\}be the set of graphs withmmmeasured variables that follow[Equation1](https://arxiv.org/html/2605.20396#S2.E1)\. For set𝐒\\mathbf\{S\}, we define itskk\-partition as a partition of its elements intokknon\-empty subsets\.

## 2Latent Variable Causal Models

In this section, we discuss several aspects of latent variable causal models\. Specifically, we describe the preliminaries and problem setting in[Section2\.1](https://arxiv.org/html/2605.20396#S2.SS1), as well as the formulation of likelihood function in[Section2\.2](https://arxiv.org/html/2605.20396#S2.SS2)\. We provide a discussion of latent variable causal models in[SectionA\.1](https://arxiv.org/html/2605.20396#A1.SS1)\.

### 2\.1Preliminaries and Problem Setting

We consider a linear latent variable causal model with DAG𝒢\\mathcal\{G\}, in which the measured variablesX=\(X1,…,Xm\)X=\(X\_\{1\},\\dots,X\_\{m\}\)and latent \(unmeasured\) variablesL=\(L1,…,Ln\)L=\(L\_\{1\},\\dots,L\_\{n\}\)follow the data generating procedure:

L=C​L\+ELandX=B​L\+EX,L=CL\+E\_\{L\}\\quad\\text\{and\}\\quad X=BL\+E\_\{X\},\(1\)whereEXE\_\{X\}andELE\_\{L\}are jointly independent noise terms that follow Gaussian distributions\. The structure of DAG𝒢\\mathcal\{G\}is defined by the support of matricesBBandCC, i\.e\.,Lj→LiL\_\{j\}\\rightarrow L\_\{i\}is an edge in𝒢\\mathcal\{G\}ifCi,j≠0C\_\{i,j\}\\neq 0andLj→XiL\_\{j\}\\rightarrow X\_\{i\}is an edge in𝒢\\mathcal\{G\}ifBi,j≠0B\_\{i,j\}\\neq 0\. For DAG𝒢\\mathcal\{G\}, we denote byB𝒢∈\{0,1\}m×nB\_\{\\mathcal\{G\}\}\\in\\\{0,1\\\}^\{m\\times n\}the binary adjacency matrix that represent the edges from latent variablesLLto measured variablesXX, and byC𝒢∈\{0,1\}n×nC\_\{\\mathcal\{G\}\}\\in\\\{0,1\\\}^\{n\\times n\}the binary adjacency matrix that represent the edges among latent variablesLL\. Without loss of generality, we assume that matricesCCandC𝒢C\_\{\\mathcal\{G\}\}are strictly upper triangular\.

LetΣX\\Sigma\_\{X\}andΣL\\Sigma\_\{L\}be the population covariance matrices of measured variablesXXand latent variablesLLrespectively\. Also letΩX\\Omega\_\{X\}andΩL\\Omega\_\{L\}be the \(diagonal\) covariance matrices of noise termsEXE\_\{X\}andELE\_\{L\}respectively\.ΣL\\Sigma\_\{L\}can be written as

ΣL=\(I−C\)−1​ΩL​\(I−C\)−⊤\.\\Sigma\_\{L\}=\(I\-C\)^\{\-1\}\\Omega\_\{L\}\(I\-C\)^\{\-\\top\}\.ByΣX=B​ΣL​B⊤\+ΩX\\Sigma\_\{X\}=B\\Sigma\_\{L\}B^\{\\top\}\+\\Omega\_\{X\}, we then have

ΣX=B​\(I−C\)−1​ΩL​\(I−C\)−⊤​B⊤\+ΩX\.\\Sigma\_\{X\}=B\(I\-C\)^\{\-1\}\\Omega\_\{L\}\(I\-C\)^\{\-\\top\}B^\{\\top\}\+\\Omega\_\{X\}\.\(2\)We say that a DAG𝒢\\mathcal\{G\}can generate a covariance matrix if there exists a parameterization of𝒢\\mathcal\{G\}such that[Equation2](https://arxiv.org/html/2605.20396#S2.E2)holds\. Furthermore, since the labeling of latent variables in general cannot be identified, we say that two DAGs are Markov equivalent if they are Markov equivalent after relabeling of latent variables\. GivenTTi\.i\.d\. samples of variablesXX, denoted as𝐃\\mathbf\{D\}with empirical covariance matrixSS, the goal is to estimate the structure𝒢\\mathcal\{G\}up to certain type of model equivalence \(specified in[Sections4](https://arxiv.org/html/2605.20396#S4)and[5](https://arxiv.org/html/2605.20396#S5)\)\.

### 2\.2Formulation of Likelihood Function

We first discuss about the indeterminacy of parameterΩL\\Omega\_\{L\}via the following lemma, since it affects how we formulate the likelihood\. The proof is given in[SectionB\.1](https://arxiv.org/html/2605.20396#A2.SS1)\.

###### Lemma 1\(Indeterminacy ofΩL\\Omega\_\{L\}\)\.

For any parametersB,C,ΩX,ΩLB,C,\\Omega\_\{X\},\\Omega\_\{L\}, andΣX\\Sigma\_\{X\}that follow[Equation2](https://arxiv.org/html/2605.20396#S2.E2), there exist parametersB~\\tilde\{B\}andC~\\tilde\{C\}withsupp⁡\(B\)=supp⁡\(B~\)\\operatorname\{supp\}\(B\)=\\operatorname\{supp\}\(\\tilde\{B\}\)andsupp⁡\(C\)=supp⁡\(C~\)\\operatorname\{supp\}\(C\)=\\operatorname\{supp\}\(\\tilde\{C\}\)such that

ΣX=B~​\(I−C~\)−1​\(I−C~\)−⊤​B~⊤\+ΩX\.\\Sigma\_\{X\}=\\tilde\{B\}\(I\-\\tilde\{C\}\)^\{\-1\}\(I\-\\tilde\{C\}\)^\{\-\\top\}\\tilde\{B\}^\{\\top\}\+\\Omega\_\{X\}\.

In other words, any covariance matrixΣX\\Sigma\_\{X\}resulting from DAG𝒢\\mathcal\{G\}and arbitraryΩL\\Omega\_\{L\}can be generated by alternative parameters from the same DAG withΩ~L=I\\tilde\{\\Omega\}\_\{L\}=I\. This implies that the parameterΩL\\Omega\_\{L\}cannot be estimated fromΣX\\Sigma\_\{X\}without additional information and further assumption\. Furthermore, since the goal is to estimate the structure𝒢\\mathcal\{G\}, this suggests that one may assumeΩL\\Omega\_\{L\}to be an identity matrix during estimation without loss of generality\. It is worth noting that such indeterminacy ofΩL\\Omega\_\{L\}has been discussed in various existing works\(Squires et al\.,[2023](https://arxiv.org/html/2605.20396#bib.bib74)\), which we make precise here, as it is crucial for formulating the likelihood\.

We now provide the likelihood formulation for the linear latent variable causal model in[Equation1](https://arxiv.org/html/2605.20396#S2.E1)\. As suggested by[Lemma1](https://arxiv.org/html/2605.20396#Thmlemma1), we setΩL=I\\Omega\_\{L\}=Iin the likelihood\. Given the empirical covariance matrixSSobtained fromTTsamples, the negative log\-likelihood is given up to additive constant by

ℒ​\(B,C,ΩX;𝐃\)\\displaystyle\\mathcal\{L\}\\left\(B,C,\\Omega\_\{X\};\\mathbf\{D\}\\right\)=T2​tr⁡\(S​\(B​\(I−C\)−1​\(I−C\)−⊤​B⊤\+ΩX\)−1\)\\displaystyle\\quad=\\frac\{T\}\{2\}\\operatorname\{tr\}\\left\(S\\left\(B\(I\-C\)^\{\-1\}\(I\-C\)^\{\-\\top\}B^\{\\top\}\+\\Omega\_\{X\}\\right\)^\{\-1\}\\right\)\+T2​log⁡det⁡\(B​\(I−C\)−1​\(I−C\)−⊤​B⊤\+ΩX\)\.\\displaystyle\\quad\\quad\\quad\+\\frac\{T\}\{2\}\\log\\operatorname\{det\}\\left\(B\(I\-C\)^\{\-1\}\(I\-C\)^\{\-\\top\}B^\{\\top\}\+\\Omega\_\{X\}\\right\)\.

## 3Score\-Based Identification of Latent Variable Causal Models

In this section, we discuss how to learn linear latent variable causal models with scoring function\. First, we introduce the notion of distribution sets and equality constraints in[Section3\.1](https://arxiv.org/html/2605.20396#S3.SS1)\. We formulate the scoring function in[Section3\.2](https://arxiv.org/html/2605.20396#S3.SS2), and show how it enables structure identification up to algebraic equivalence in[Section3\.3](https://arxiv.org/html/2605.20396#S3.SS3)\. We then discuss about the BIC score in[Section3\.4](https://arxiv.org/html/2605.20396#S3.SS4)\.

### 3\.1Distribution Sets and Equality Constraints

We describe the notion of distribution set that is a key ingredient of our score\-based search procedure\. It refers to the set of marginal distributions generated by a specific structure\.

###### Definition 1\(Distribution set\)\.

The distribution set of DAG𝒢\\mathcal\{G\}, denoted byℳ​\(𝒢\)\\mathcal\{M\}\(\\mathcal\{G\}\), is defined as

ℳ\(𝒢\)≔\{\\displaystyle\\mathcal\{M\}\(\\mathcal\{G\}\)\\coloneqq\\\{B​\(I−C\)−1​ΩL​\(I−C\)−⊤​B⊤\+ΩX:\\displaystyle B\(I\-C\)^\{\-1\}\\Omega\_\{L\}\(I\-C\)^\{\-\\top\}B^\{\\top\}\+\\Omega\_\{X\}:supp⁡\(B\)⊆supp⁡\(B𝒢\),supp⁡\(C\)⊆supp⁡\(C𝒢\),\\displaystyle\\,\\operatorname\{supp\}\(B\)\\subseteq\\operatorname\{supp\}\(B\_\{\\mathcal\{G\}\}\),\\operatorname\{supp\}\(C\)\\subseteq\\operatorname\{supp\}\(C\_\{\\mathcal\{G\}\}\),ΩX∈diag\(ℝ\>0m\),ΩL∈diag\(ℝ\>0n\)\}\.\\displaystyle\\,\\Omega\_\{X\}\\in\\operatorname\{diag\}\(\\mathbb\{R\}\_\{\>0\}^\{m\}\),\\Omega\_\{L\}\\in\\operatorname\{diag\}\(\\mathbb\{R\}\_\{\>0\}^\{n\}\)\\\}\.

Specifically,ℳ​\(𝒢\)\\mathcal\{M\}\(\\mathcal\{G\}\)is the set of covariances matricesΣX\\Sigma\_\{X\}that can be generated by DAG𝒢\\mathcal\{G\}by varying the parameters in matricesBB,CC,ΩX\\Omega\_\{X\}, andΩL\\Omega\_\{L\}\. Moreover, sinceC𝒢C\_\{\\mathcal\{G\}\}is acyclic by assumption, we have\(I−C\)−1=∑k=0n−1Ck\(I\-C\)^\{\-1\}=\\sum\_\{k=0\}^\{n\-1\}C^\{k\}\. It follows that[Equation2](https://arxiv.org/html/2605.20396#S2.E2)is a polynomial map, and thus the distribution setℳ​\(𝒢\)\\mathcal\{M\}\(\\mathcal\{G\}\)is, by Tarski–Seidenberg theorem\(Benedetti & Risler,[1990](https://arxiv.org/html/2605.20396#bib.bib7)\), a semialgebraic set\. Note that a set is said to be*semialgebraic*if it can be equivalently represented by a finite number of polynomial equalities and inequalities\(Benedetti & Risler,[1990](https://arxiv.org/html/2605.20396#bib.bib7)\)\.

Structure𝒢\\mathcal\{G\}imposes various types of equality \(i\.e\., algebraic\) constraints on the covariance matrices, such as conditional independence \(i\.e\., vanishing partial correlation\) constraints\(Spirtes et al\.,[2001](https://arxiv.org/html/2605.20396#bib.bib72)\), rank deficiency \(i\.e\., vanishing determinant\) constraints\(Spirtes et al\.,[2001](https://arxiv.org/html/2605.20396#bib.bib72); Sullivant et al\.,[2010](https://arxiv.org/html/2605.20396#bib.bib76)\), and possibly Verma constraints\(Verma & Pearl,[1991](https://arxiv.org/html/2605.20396#bib.bib79)\)\. We refer the readers toDrton \([2018](https://arxiv.org/html/2605.20396#bib.bib29)\)for an overview\. LetH​\(𝒢\)H\(\\mathcal\{G\}\)be the set of equality constraints imposed by structure𝒢\\mathcal\{G\}on the distribution setℳ​\(𝒢\)\\mathcal\{M\}\(\\mathcal\{G\}\), andℍm≔⋃𝒢∈𝔾mH​\(𝒢\)\\mathbb\{H\}^\{m\}\\coloneqq\\bigcup\_\{\\mathcal\{G\}\\in\\mathbb\{G\}^\{m\}\}H\(\\mathcal\{G\}\)be the set of possible equality constraints imposed by any structure𝒢\\mathcal\{G\}\(withmmmeasured variables\)\. Two structures𝒢1\\mathcal\{G\}\_\{1\}and𝒢2\\mathcal\{G\}\_\{2\}are said to be*algebraic equivalent*if they lead to the same equality constraints, i\.e\.,H​\(𝒢1\)=H​\(𝒢2\)H\(\\mathcal\{G\}\_\{1\}\)=H\(\\mathcal\{G\}\_\{2\}\)\(van Ommen & Mooij,[2017](https://arxiv.org/html/2605.20396#bib.bib78)\)\.111In the terminology of algebraic geometry,ℳ​\(𝒢1\)\\mathcal\{M\}\(\\mathcal\{G\}\_\{1\}\)andℳ​\(𝒢2\)\\mathcal\{M\}\(\\mathcal\{G\}\_\{2\}\)share the same vanishing ideal or Zariski closure\(Cox et al\.,[2015](https://arxiv.org/html/2605.20396#bib.bib24)\)\.

Furthermore, letdim\(𝒢\)\\dim\(\\mathcal\{G\}\)denote the model dimension or degrees of freedom of DAG𝒢\\mathcal\{G\}for the marginal over the observed variables, which can be viewed as the number of free parameters for the distribution setℳ​\(𝒢\)\\mathcal\{M\}\(\\mathcal\{G\}\)\. In general, the degrees of freedom are not necessarily equal to the number of parameters \(i\.e\., sum of number of edges and measured variables\) in the presence of latent variables\(Geiger et al\.,[1996](https://arxiv.org/html/2605.20396#bib.bib32),[2001](https://arxiv.org/html/2605.20396#bib.bib33)\)\. In[Sections4](https://arxiv.org/html/2605.20396#S4)and[5](https://arxiv.org/html/2605.20396#S5), we further characterize the degrees of freedom under specific structural assumptions\.

### 3\.2Formulation of Scoring Function

We now provide the formulation of the scoring function for identifying linear latent variable causal models\. Specifically, our score\-based method involves searching for the structure with the smallest degrees of freedom that can generate the covariance matrix\. Given structure𝒢\\mathcal\{G\}with samples𝐃\\mathbf\{D\}and empirical covariance matrixSS, the scoring function is

scoredim⁡\(𝒢,𝐃\)≔\{dim\(𝒢\)if​𝒢​can generate​S,∞otherwise\.\\operatorname\{score\}\_\{\\textrm\{dim\}\}\(\\mathcal\{G\},\\mathbf\{D\}\)\\coloneqq\\begin\{cases\}\\dim\(\\mathcal\{G\}\)&\\text\{ if \}\\mathcal\{G\}\\text\{ can generate \}S,\\\\ \\infty&\\text\{ otherwise\.\}\\end\{cases\}To determine whether structure𝒢\\mathcal\{G\}can generateSS, one may minimize the squared errors betweenSSand the covariance matrix parameterized by𝒢\\mathcal\{G\}, or compare the maximum likelihood w\.r\.t\.𝒢\\mathcal\{G\}\(e\.g\., see[Equation3](https://arxiv.org/html/2605.20396#S3.E3)\) to the likelihood ofSS\. Moreover, we show in[Section3\.4](https://arxiv.org/html/2605.20396#S3.SS4)that the scoring function satisfies the property of score equivalence\.

Similar scoring function has been discussed byRaskutti & Uhler \([2014](https://arxiv.org/html/2605.20396#bib.bib56)\)\. Roughly speaking, there may exist multiple structures that can generate the same distribution; the scoring function identifies one of them with the smallest degrees of freedom\. In[Sections3\.3](https://arxiv.org/html/2605.20396#S3.SS3),[4](https://arxiv.org/html/2605.20396#S4)and[5](https://arxiv.org/html/2605.20396#S5), we discuss how this scoring function identifies structures up to different types of model equivalence in the large sample limit\. Specifically, we show in[Section3\.3](https://arxiv.org/html/2605.20396#S3.SS3)that it yields a structure that is algebraic equivalent to the ground truth\.

Since the key idea is to identify the structure that generates the covariance matrix \(in the large sample limit\) with the smallest degrees of freedom, different types of scoring functions that can achieve so can also be used\. In[Section3\.4](https://arxiv.org/html/2605.20396#S3.SS4), we further discuss the use of the BIC score\.

### 3\.3Identifying Structures up to Algebraic Equivalence

Having formulated the scoring function in[Section3\.2](https://arxiv.org/html/2605.20396#S3.SS2), the question remains as how to leverage it to identify the underlying structure𝒢\\mathcal\{G\}\. To do so, a key ingredient is to establish the correspondence between the covariance matrix and the structure𝒢\\mathcal\{G\}\. As discussed in[Section3\.1](https://arxiv.org/html/2605.20396#S3.SS1), the structure𝒢\\mathcal\{G\}imposes different types of constraints on the entries of covariance matrices, including equality and inequality constraints\. Here, we adopt the following assumption which requires that the equality constraints are imposed by the structure𝒢\\mathcal\{G\}\.

###### Assumption 1\(Generalized faithfulness\(Ghassami et al\.,[2020](https://arxiv.org/html/2605.20396#bib.bib34)\)\)\.

A distributionΣX\\Sigma\_\{X\}is said to be generalized faithful to DAG𝒢\\mathcal\{G\}if the entries ofΣX\\Sigma\_\{X\}satisfy an equality constraintκ∈ℍm\\kappa\\in\\mathbb\{H\}^\{m\}only ifκ∈H​\(𝒢\)\\kappa\\in H\(\\mathcal\{G\}\)\.

It is worth noting that different types of faithfulness assumptions have been adopted in causal discovery\(Spirtes et al\.,[2001](https://arxiv.org/html/2605.20396#bib.bib72); Ghassami et al\.,[2020](https://arxiv.org/html/2605.20396#bib.bib34); Huang et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib38)\)to relate the constraints of the distributions \(e\.g\., conditional independence and rank deficiency constraints\) to the underlying structure\. This is often motivated by the fact that the set of parameters violating these assumptions has Lebesgue measure zero \(see, e\.g\.,Ghassami et al\. \([2020](https://arxiv.org/html/2605.20396#bib.bib34), Proposition 8\)\)\.

We then present the following result that describes the notion of equivalence achieved by minimizing the scoring function\. The proof is provided in[SectionB\.3](https://arxiv.org/html/2605.20396#A2.SS3), which is partly inspired by the proof ofGhassami et al\. \([2020](https://arxiv.org/html/2605.20396#bib.bib34), Theorem 3\)\.

###### Theorem 1\(Algebraic equivalence\)\.

Suppose the true DAG𝒢∗\\mathcal\{G\}^\{\*\}and the distributionΣX\\Sigma\_\{X\}satisfy the generalized faithfulness assumption\. Let𝒢^∈arg​min𝒢∈𝔾m⁡scoredim⁡\(𝒢,𝐃\)\\hat\{\\mathcal\{G\}\}\\in\\operatorname\*\{arg\\,min\}\_\{\\mathcal\{G\}\\in\\mathbb\{G\}^\{m\}\}\\operatorname\{score\}\_\{\\textrm\{dim\}\}\(\\mathcal\{G\},\\mathbf\{D\}\)\. Then,𝒢^\\hat\{\\mathcal\{G\}\}and𝒢∗\\mathcal\{G\}^\{\*\}are algebraic equivalent, i\.e\.,H​\(𝒢^\)=H​\(𝒢∗\)H\(\\hat\{\\mathcal\{G\}\}\)=H\(\\mathcal\{G\}^\{\*\}\), in the large sample limit\.

In general, relating the estimated structure to the true one, which are algebraic equivalent, can be challenging without any restrictions on the structures\. In[Sections4](https://arxiv.org/html/2605.20396#S4)and[5](https://arxiv.org/html/2605.20396#S5), we show that, under specific structural assumptions,[Theorem1](https://arxiv.org/html/2605.20396#Thmtheorem1)helps achieve notions of model equivalence that are more fine\-grained than algebraic equivalence \(including Markov equivalence in[Section4](https://arxiv.org/html/2605.20396#S4)\)\. Therefore, a general recipe may involve identifying suitable structural assumptions that allow algebraic equivalence to translate into more fine\-grained notions of model equivalence\. This enables the application of the score based procedure in[Theorem1](https://arxiv.org/html/2605.20396#Thmtheorem1), given an appropriate characterization of the degrees of freedom\. We give a further discussion of generalized faithfulness and algebraic equivalence in[SectionsA\.2](https://arxiv.org/html/2605.20396#A1.SS2)and[A\.3](https://arxiv.org/html/2605.20396#A1.SS3), respectively\.

### 3\.4Remark on the BIC Score

The scoring function discussed in[Section3\.2](https://arxiv.org/html/2605.20396#S3.SS2)is justified in the large sample limit and may not perform well for finite\-sample cases\. We consider the BIC score\(Schwarz,[1978](https://arxiv.org/html/2605.20396#bib.bib62); Chickering,[2002](https://arxiv.org/html/2605.20396#bib.bib19)\)that maximizes the likelihood while penalizing the degrees of freedom of structure𝒢\\mathcal\{G\}:

scoreBIC⁡\(𝒢,𝐃\)≔scoreℒ⁡\(𝒢,𝐃\)\+log⁡T2​dim\(𝒢\)\.\\operatorname\{score\}\_\{\\textrm\{BIC\}\}\(\\mathcal\{G\},\\mathbf\{D\}\)\\coloneqq\\operatorname\{score\}\_\{\\mathcal\{L\}\}\(\\mathcal\{G\},\\mathbf\{D\}\)\+\\frac\{\\log T\}\{2\}\\dim\(\\mathcal\{G\}\)\.wherescoreℒ⁡\(𝒢,𝐃\)\\operatorname\{score\}\_\{\\mathcal\{L\}\}\(\\mathcal\{G\},\\mathbf\{D\}\)denotes the optimal negative log\-likelihood w\.r\.t\. structure𝒢\\mathcal\{G\}, given by

scoreℒ⁡\(𝒢,𝐃\)≔min\(B,C,ΩX\):supp⁡\(B\)⊆supp⁡\(B𝒢\),supp⁡\(C\)⊆supp⁡\(C𝒢\),ΩX∈diag⁡\(ℝ\>0m\)⁡ℒ​\(B,C,ΩX;𝐃\)\.\\operatorname\{score\}\_\{\\mathcal\{L\}\}\(\\mathcal\{G\},\\mathbf\{D\}\)\\coloneqq\\min\_\{\\begin\{subarray\}\{c\}\(B,C,\\Omega\_\{X\}\):\\\\ \\operatorname\{supp\}\(B\)\\subseteq\\operatorname\{supp\}\(B\_\{\\mathcal\{G\}\}\),\\\\ \\operatorname\{supp\}\(C\)\\subseteq\\operatorname\{supp\}\(C\_\{\\mathcal\{G\}\}\),\\\\ \\Omega\_\{X\}\\in\\operatorname\{diag\}\(\\mathbb\{R\}\_\{\>0\}^\{m\}\)\\end\{subarray\}\}\\mathcal\{L\}\\left\(B,C,\\Omega\_\{X\};\\mathbf\{D\}\\right\)\.\(3\)Since it may not be straightforward to derive a closed\-form solution, various numerical solvers or continuous optimization methods, such as L\-BFGS\(Byrd et al\.,[1995](https://arxiv.org/html/2605.20396#bib.bib15)\)and gradient descent, as well as the expectation\-maximization algorithm\(Dempster et al\.,[1977](https://arxiv.org/html/2605.20396#bib.bib27)\), can be used to compute the maximum likelihood above\.

It is worth noting that the BIC score has been widely adopted in score\-based causal discovery\(Chickering,[2002](https://arxiv.org/html/2605.20396#bib.bib19)\)\.Haughton \([1988](https://arxiv.org/html/2605.20396#bib.bib36)\)showed that it is an asymptotic approximation for the log marginal likelihood of*curved*exponential families, which include Gaussian DAG models without latent variables\(Geiger et al\.,[2001](https://arxiv.org/html/2605.20396#bib.bib33); Richardson & Spirtes,[2002](https://arxiv.org/html/2605.20396#bib.bib57)\)\. In the presence of latent variables, the models are*stratified*exponential families, and complications arise in using BIC for model selection\. Although the typical theoretical justifications of using BIC\(Schwarz,[1978](https://arxiv.org/html/2605.20396#bib.bib62); Haughton,[1988](https://arxiv.org/html/2605.20396#bib.bib36)\)may not apply for identifying latent variable causal models in our setting, we apply it in place ofscoredim⁡\(𝒢,𝐃\)\\operatorname\{score\}\_\{\\textrm\{dim\}\}\(\\mathcal\{G\},\\mathbf\{D\}\)in our experiments, since the latter is justified in the large sample limit and may not perform well for finite\-sample cases\. Surprisingly, using the BIC score leads to a superior empirical performance, specifically under the structural assumptions described in[Sections4](https://arxiv.org/html/2605.20396#S4)and[5](https://arxiv.org/html/2605.20396#S5)\. This suggests that BIC may be a valid scoring criterion in these cases\. Therefore, future works involve studying the theoretical justifications of using BIC score under these structural assumptions\.

Recall that a scoring function is*score equivalent*if every pair of Markov equivalent structures have the same score\(Chickering,[2002](https://arxiv.org/html/2605.20396#bib.bib19)\)\. This is a desirable property for score\-based procedure as it implies that we can search in the space of Markov equivalence classes \(MECs\) instead of DAGs\. That is, one does not have to compute the score multiple times for the DAGs in the same MEC, which may help improve the runtime\. We show that our scoring functions satisfy such a property, with a proof given in[SectionB\.2](https://arxiv.org/html/2605.20396#A2.SS2)\.

###### Proposition 1\(Score equivalence\)\.

Suppose that DAGs𝒢1\\mathcal\{G\}\_\{1\}and𝒢2\\mathcal\{G\}\_\{2\}are Markov equivalent\. Then, we havescoredim⁡\(𝒢1,𝐃\)=scoredim⁡\(𝒢2,𝐃\)\\operatorname\{score\}\_\{\\textrm\{dim\}\}\(\\mathcal\{G\}\_\{1\},\\mathbf\{D\}\)=\\operatorname\{score\}\_\{\\textrm\{dim\}\}\(\\mathcal\{G\}\_\{2\},\\mathbf\{D\}\)andscoreBIC⁡\(𝒢1,𝐃\)=scoreBIC⁡\(𝒢2,𝐃\)\\operatorname\{score\}\_\{\\textrm\{BIC\}\}\(\\mathcal\{G\}\_\{1\},\\mathbf\{D\}\)=\\operatorname\{score\}\_\{\\textrm\{BIC\}\}\(\\mathcal\{G\}\_\{2\},\\mathbf\{D\}\)\.

## 4Linear 1\-Factor Latent Variable Models

In the previous section, we show that the scoring function can produce a structure algebraic equivalent to the ground truth\. We now discuss how such result helps estimate a structure up to Markov equivalence\. In this section, we focus on the structural assumption bySilva et al\. \([2003](https://arxiv.org/html/2605.20396#bib.bib67),[2006](https://arxiv.org/html/2605.20396#bib.bib68)\)\.222In our setting, it may be sufficient to require that each latent variable has at least two measured variables as children\.

###### Assumption 2\(Silva et al\. \([2003](https://arxiv.org/html/2605.20396#bib.bib67)\)\)\.

Each measured variable has a single latent parent, and each latent variable has at least three measured variables as children\.

An example illustrating the above assumption is provided in[Figure1](https://arxiv.org/html/2605.20396#S4.F1)\.Silva et al\. \([2003](https://arxiv.org/html/2605.20396#bib.bib67)\)proposed a search procedure based on statistical tests of tetrad constraints that can identify structures under this assumption\. In this section, we develop a score\-based method based on this structural assumption\. We first characterize the degrees of freedom of the structure in[Section4\.1](https://arxiv.org/html/2605.20396#S4.SS1), as required by the scoring function\. We then establish the consistency and provide an exact score\-based search procedure in[Section4\.2](https://arxiv.org/html/2605.20396#S4.SS2)\. We also develop a continuous search procedure in[Section4\.3](https://arxiv.org/html/2605.20396#S4.SS3)that may be more computationally efficient\.

![Refer to caption](https://arxiv.org/html/2605.20396v1/x1.png)Figure 1:Example of 1\-factor latent variable model\.### 4\.1Degrees of Freedom

The scoring function requires a proper specification of the degrees of freedom during the search procedure\. For the structural assumption in[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2), the degrees of freedom, as one may expect, equals the number of edges in DAG𝒢\\mathcal\{G\}plus the number of measured variables\. Here, the number of edges include those among the latent variables and those from the latent variables to the measured ones\. The proof follows straightforwardly from parameter identifiability under[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)\(Bollen,[1989](https://arxiv.org/html/2605.20396#bib.bib12)\), which is provided in[SectionB\.4](https://arxiv.org/html/2605.20396#A2.SS4)for completeness\.

###### Proposition 2\(Degrees of freedom\)\.

Suppose that DAG𝒢\\mathcal\{G\}satisfies[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)\. Then,dim⁡\(𝒢\)=\|𝒢\|\+m\\operatorname\{dim\}\(\\mathcal\{G\}\)=\|\\mathcal\{G\}\|\+m\.

To illustrate, the degrees of freedom of the example in[Figure1](https://arxiv.org/html/2605.20396#S4.F1)are simply equal to2424\. The above property holds in many other settings such as the typical setting without latent confounders\(Chickering,[2002](https://arxiv.org/html/2605.20396#bib.bib19)\), as well as those with bow\-free acyclic mixed graphs\(Brito & Pearl,[2002](https://arxiv.org/html/2605.20396#bib.bib13)\)and cycles \(excluding 2\-cycles\)\(Amendola et al\.,[2020](https://arxiv.org/html/2605.20396#bib.bib3)\)\. Note that such property, while desirable, does not hold in general for structures with latent variables\(Geiger et al\.,[1996](https://arxiv.org/html/2605.20396#bib.bib32)\)\. For instance, the degrees of freedom of the structures that we consider in[Section5](https://arxiv.org/html/2605.20396#S5)are generally not equal to\|𝒢\|\+m\|\\mathcal\{G\}\|\+m\.

### 4\.2Consistency and Exact Score\-Based Search

Having characterized the degrees of freedom, we now establish the correctness of score\-based approach under[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)and accordingly develop an exact search procedure\. Specifically, under[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)and the generalized faithfulness assumption, we show that the structure with the optimal score is Markov equivalent to the true structure\.

###### Theorem 2\(Correctness\)\.

Suppose that the true DAG𝒢∗\\mathcal\{G\}^\{\*\}and the distributionΣX\\Sigma\_\{X\}satisfy the generalized faithfulness assumption, and that𝒢∗\\mathcal\{G\}^\{\*\}satisfies[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)\. Let𝒢^\\hat\{\\mathcal\{G\}\}be a global minimizer of the following optimization problem:

min𝒢∈𝔾m\\displaystyle\\min\_\{\\mathcal\{G\}\\in\\mathbb\{G\}^\{m\}\}scoredim⁡\(𝒢,𝐃\)\\displaystyle\\operatorname\{score\}\_\{\\textrm\{dim\}\}\(\\mathcal\{G\},\\mathbf\{D\}\)\(4\)subject​to\\displaystyle\\mathrm\{subject\\ to\}𝒢​satisfies[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2),\\displaystyle\\mathcal\{G\}\\text\{ satisfies \\lx@cref\{creftype~refnum\}\{assumption:silva\_graphical\_criterion\}\},wheredim\(𝒢\)=\|𝒢\|\+m\\dim\(\\mathcal\{G\}\)=\|\\mathcal\{G\}\|\+m\. Then,𝒢^\\hat\{\\mathcal\{G\}\}and𝒢∗\\mathcal\{G\}^\{\*\}are Markov equivalent in the large sample limit\.

The proof can be found in[SectionB\.5](https://arxiv.org/html/2605.20396#A2.SS5), which leverages[Theorem1](https://arxiv.org/html/2605.20396#Thmtheorem1)that shows how the scoring function produces a structure that is algebraic equivalent to the true structure\. Moreover, the theorem above indicates that one could perform exact search for all structures under[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)\. A naive approach is to iterate over all possible structures and check if each of them satisfy[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)\. This may be computationally infeasible because much of the time may be spent on structures that do not fall within the model class\.

The question is then how to efficiently enumerate and perform exact search for these structures\. We provide an algorithm to do so in[Algorithm1](https://arxiv.org/html/2605.20396#alg1)\. Leveraging the score equivalence property in[Proposition1](https://arxiv.org/html/2605.20396#Thmproposition1), we consider only structures that are not Markov equivalent to one another, since they are indistinguishable based on[Theorem2](https://arxiv.org/html/2605.20396#Thmtheorem2)and give rise to the same score\. First, we generate the possible structuresC𝒢C\_\{\\mathcal\{G\}\}among the latent variables that are not Markov equivalent to one another\. To construct the structureB𝒢B\_\{\\mathcal\{G\}\}from latent variables to measured variables, we then find all ordered partitions of measured variables and add each subset from the partition to be the children of each latent variable\. We compute the score for each structure enumerated by[Algorithm1](https://arxiv.org/html/2605.20396#alg1), and find the structure with the optimal score\. Under[Theorem2](https://arxiv.org/html/2605.20396#Thmtheorem2), such an exact search procedure will output a structure Markov equivalent to the true one\.

Algorithm 1Enumerating structures under[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)Input:Measured variables

X1,…,XmX\_\{1\},\\dots,X\_\{m\}
Output:Set of DAGs

𝐀\\mathbf\{A\}
Initialize

𝐀\\mathbf\{A\}as an empty set;

for

n=1n=1to

⌊m/3⌋\\lfloor m/3\\rfloordo

foreachlatent MEC with

nnvariablesdo

Generate a latent DAG

C𝒢C\_\{\\mathcal\{G\}\}from the latent MEC;

foreachordered

nn\-partition

\(𝐏j\)j=1n\(\\mathbf\{P\}\_\{j\}\)\_\{j=1\}^\{n\}of

\{Xi\}i=1m\\\{X\_\{i\}\\\}\_\{i=1\}^\{m\}do

if

\|𝐏j\|≥3\|\\mathbf\{P\}\_\{j\}\|\\geq 3for

j=1,…,nj=1,\\dots,nthen

Construct DAG

𝒢\\mathcal\{G\}with latent DAG

C𝒢C\_\{\\mathcal\{G\}\}and each latent

LjL\_\{j\}pointing to the variables in

𝐏j\\mathbf\{P\}\_\{j\};

if

𝒢\\mathcal\{G\}is not Markov equivalent to all DAGs in

𝐀\\mathbf\{A\}then

𝐀←𝐀∪\{𝒢\}\\mathbf\{A\}\\leftarrow\\mathbf\{A\}\\cup\\\{\\mathcal\{G\}\\\};

returnset

𝐀\\mathbf\{A\}

### 4\.3Continuous Search

The exact search procedure presented in the previous section requires computing the score for each structure satisfying the structural assumption, which can be computationally intensive when there is a large number of variables\. For instance, when using the BIC score, each computation involves solving a continuous optimization problem in[Equation3](https://arxiv.org/html/2605.20396#S3.E3); the same applies toscoredim⁡\(𝒢,𝐃\)\\operatorname\{score\}\_\{\\textrm\{dim\}\}\(\\mathcal\{G\},\\mathbf\{D\}\)\. This is inherent to discrete score\-based search that assigns a score to each structure\. A key question naturally arises: how do we unify the structure search part and likelihood computation into a single continuous optimization problem? Such a unified procedure helps reduce the computational burden of separately computing the score for each structure in a discrete search\. Furthermore, this aligns with recent studies in continuous optimization for causal discovery\(Zheng et al\.,[2018](https://arxiv.org/html/2605.20396#bib.bib90); Ng et al\.,[2020](https://arxiv.org/html/2605.20396#bib.bib48); Vowels et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib81)\)\.

We first provide a reformulation of[Equation4](https://arxiv.org/html/2605.20396#S4.E4)with the BIC score that is more amenable to continuous optimization\. The key lies in characterizing the penalty term\|𝒢\|\|\\mathcal\{G\}\|and the constraint involving[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)\. Specifically, we solve the following constrained optimization problem:

minMB∈\{0,1\}m×n¯,MC∈\{0,1\}n¯×n¯,B∈ℝm×n¯,C∈𝕌n¯,ΩX∈diag⁡\(ℝ\>0m\)\\displaystyle\\min\_\{\\begin\{subarray\}\{c\}M\_\{B\}\\in\\\{0,1\\\}^\{m\\times\\bar\{n\}\},\\\\ M\_\{C\}\\in\\\{0,1\\\}^\{\\bar\{n\}\\times\\bar\{n\}\},\\\\ B\\in\\mathbb\{R\}^\{m\\times\\bar\{n\}\},C\\in\\mathbb\{U\}^\{\\bar\{n\}\},\\\\ \\Omega\_\{X\}\\in\\operatorname\{diag\}\(\\mathbb\{R\}\_\{\>0\}^\{m\}\)\\end\{subarray\}\}\\quad\(1Tℒ\(MB⊙B,MC⊙C,ΩX;𝐃\)\\displaystyle\\Big\(\\frac\{1\}\{T\}\\mathcal\{L\}\\left\(M\_\{B\}\\odot B,M\_\{C\}\\odot C,\\Omega\_\{X\};\\mathbf\{D\}\\right\)\+λ∥MB∥1\+λ∥MC∥1\)\\displaystyle\\qquad\+\\lambda\\\|M\_\{B\}\\\|\_\{1\}\+\\lambda\\\|M\_\{C\}\\\|\_\{1\}\\Big\)subject​to\\displaystyle\\mathrm\{subject\\ to\}\\quad∑k=1n¯\(MB\)i,k−1=0,i∈\[m\],\\displaystyle\\sum\_\{k=1\}^\{\\bar\{n\}\}\(M\_\{B\}\)\_\{i,k\}\-1=0,\\;i\\in\[m\],\(5\)\(\(∑k=1m\(MB\)k,j\+∑k=1n¯\(MC\)k,j\)\\displaystyle\\Bigg\(\\bigg\(\\sum\_\{k=1\}^\{m\}\(M\_\{B\}\)\_\{k,j\}\+\\sum\_\{k=1\}^\{\\bar\{n\}\}\(M\_\{C\}\)\_\{k,j\}\\bigg\)\(∑k=1m\(MB\)k,j−3\)\)≥0,j∈\[n¯\],\\displaystyle\\qquad\\bigg\(\\sum\_\{k=1\}^\{m\}\(M\_\{B\}\)\_\{k,j\}\-3\\bigg\)\\Bigg\)\\geq 0,\\;j\\in\[\\bar\{n\}\],wheren¯=⌊m/3⌋\\bar\{n\}=\\lfloor m/3\\rflooris an upper bound of the number of latent variables andλ=log⁡T/2​T\\lambda=\\log T/2T\. In the formulation above, the matricesMBM\_\{B\}andMCM\_\{C\}can be viewed as the support matrices ofBBandCC; they act as binary masks which indicate which edges are present in the structure\. Furthermore, the two constraints in[Equation5](https://arxiv.org/html/2605.20396#S4.E5)serve as a characterization of[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)using the support matricesMBM\_\{B\}andMCM\_\{C\}\. Specifically, the first constraint requires each row ofMBM\_\{B\}to have one nonzero entry \(i\.e\., each measured variable has one a single latent parent\)\. The second constraint requires each column ofMBM\_\{B\}andMCM\_\{C\}to satisfy the following: either the column ofMBM\_\{B\}has at least three nonzero entries, or the column ofMBM\_\{B\}andMCM\_\{C\}have no nonzero entries \(i\.e\., each latent variable has at least33measured variables as children, or no child at all\)\.

We now discuss how to solve[Equation5](https://arxiv.org/html/2605.20396#S4.E5)using continuous constrained optimization procedure\. We first introduce slack variableti≥0t\_\{i\}\\geq 0and convert the inequality constraints into equality constraints\. To estimate the binary matricesMBM\_\{B\}andMCM\_\{C\}, we apply the Gumbel\-Softmax technique\(Maddison et al\.,[2017](https://arxiv.org/html/2605.20396#bib.bib46); Jang et al\.,[2017](https://arxiv.org/html/2605.20396#bib.bib40)\)that is widely used to sample and approximate samples from a categorical distribution, which has also been adopted in continuous optimization approaches for causal discovery\(Ng et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib49); Brouillard et al\.,[2020](https://arxiv.org/html/2605.20396#bib.bib14)\)\. Specifiaclly, we apply Gumbel\-Sigmoid for each entry ofMCM\_\{C\}, and Gumbel\-Softmax withn¯\\bar\{n\}categories for each row ofMBM\_\{B\}; the latter also incorporates the first constraint of[Equation5](https://arxiv.org/html/2605.20396#S4.E5)that requires each row ofMBM\_\{B\}to have one nonzero entry\. The resulting continuous constrained optimization problem can then be solved using standard methods such as augmented Lagrangian method and quadratic penalty method\(Bertsekas,[1982](https://arxiv.org/html/2605.20396#bib.bib9),[1999](https://arxiv.org/html/2605.20396#bib.bib10); Nocedal & Wright,[2006](https://arxiv.org/html/2605.20396#bib.bib51)\)\. These methods transform the constrained problem into a series of unconstrained problems, each of which can be solved via continuous optimization methods such as gradient descent or L\-BFGS\(Byrd et al\.,[1995](https://arxiv.org/html/2605.20396#bib.bib15)\)\. In this work, we adopt augmented Lagragian method that is commonly used in causal discovery with continuous optimization\(Zheng et al\.,[2018](https://arxiv.org/html/2605.20396#bib.bib90); Vowels et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib81)\)\.

## 5Linear Latent Hierarchical Structures

The structural assumption in[Section4](https://arxiv.org/html/2605.20396#S4)requires \(i\) each measured variable to have only one latent parent and \(ii\) each latent variable to have measured children\. In real\-world cases, the structure may be more complex – the measurement model may not be a tree and some latent variables may not have measured children\. Thus, we also consider a more general assumption formulated byHuang et al\. \([2022](https://arxiv.org/html/2605.20396#bib.bib38)\)\.

We use𝐋\\mathbf\{L\}to denote a set of latent variables and writeCh𝒢​\(𝐋\)≔⋃Li∈𝐋Ch𝒢​\(Li\)\\text\{Ch\}\_\{\\mathcal\{G\}\}\(\\mathbf\{L\}\)\\coloneqq\\bigcup\_\{L\_\{i\}\\in\\mathbf\{L\}\}\\text\{Ch\}\_\{\\mathcal\{G\}\}\(L\_\{i\}\),Pa𝒢​\(𝐋\)≔⋃Li∈𝐋Pa𝒢​\(Li\)\\text\{Pa\}\_\{\\mathcal\{G\}\}\(\\mathbf\{L\}\)\\coloneqq\\bigcup\_\{L\_\{i\}\\in\\mathbf\{L\}\}\\text\{Pa\}\_\{\\mathcal\{G\}\}\(L\_\{i\}\), andDe𝒢​\(𝐋\)≔⋃Li∈𝐋De𝒢​\(Li\)\\text\{De\}\_\{\\mathcal\{G\}\}\(\\mathbf\{L\}\)\\coloneqq\\bigcup\_\{L\_\{i\}\\in\\mathbf\{L\}\}\\text\{De\}\_\{\\mathcal\{G\}\}\(L\_\{i\}\), as their parents, children, and descendants, respectively\. We now explain the notions of pure children and latent atomic cover, which serve as the fundamental building blocks of the whole structure\.

###### Definition 2\(Pure children set\(Huang et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib38)\)\)\.

LetPCh𝒢​\(𝐋\)\\text\{PCh\}\_\{\\mathcal\{G\}\}\(\\mathbf\{L\}\)denote the collection of pure\-child sets of𝐋\\mathbf\{L\}\. We say that𝐕\\mathbf\{V\}is a pure\-child set of𝐋\\mathbf\{L\}, written𝐕∈PCh𝒢​\(𝐋\)\\mathbf\{V\}\\in\\text\{PCh\}\_\{\\mathcal\{G\}\}\(\\mathbf\{L\}\), if \(i\)𝐕∩𝐋=∅\\mathbf\{V\}\\cap\\mathbf\{L\}=\\emptyset, \(ii\)𝐕⊆Ch𝒢​\(𝐋\)\\mathbf\{V\}\\subseteq\\text\{Ch\}\_\{\\mathcal\{G\}\}\(\\mathbf\{L\}\), \(iii\)Pa𝒢​\(𝐕\)=𝐋\\text\{Pa\}\_\{\\mathcal\{G\}\}\(\\mathbf\{V\}\)=\\mathbf\{L\}, and \(iv\)De𝒢​\(𝐕\)∩𝐋=∅\\text\{De\}\_\{\\mathcal\{G\}\}\(\\mathbf\{V\}\)\\cap\\mathbf\{L\}=\\emptyset\.

###### Definition 3\(Latent atomic cover\(Huang et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib38)\)\)\.

Let𝐋\\mathbf\{L\}be a set of latent variables in𝒢\\mathcal\{G\}with\|𝐋\|=k\|\\mathbf\{L\}\|=k\.𝐋\\mathbf\{L\}is a latent atomic cover if the following conditions hold:

- \(i\)There exists a set of variables𝐂\\mathbf\{C\}such that𝐂∈PCh𝒢​\(𝐋\)\\mathbf\{C\}\\in\\text\{PCh\}\_\{\\mathcal\{G\}\}\(\\mathbf\{L\}\)and\|𝐂\|≥k\+1\|\\mathbf\{C\}\|\\geq k\+1\.
- \(ii\)IfV∈Pa𝒢​\(𝐋\)\{V\}\\in\\text\{Pa\}\_\{\\mathcal\{G\}\}\(\\mathbf\{L\}\), thenV\{V\}is also the parent of all elements in𝐋\\mathbf\{L\}\. Plus, there exists a set𝐍\\mathbf\{N\}with\|𝐍\|≥k\+1\|\\mathbf\{N\}\|\\geq k\+1such that𝐍∩𝐋=∅\\mathbf\{N\}\\cap\\mathbf\{L\}=\\emptyset, every element in𝐍\\mathbf\{N\}is a neighbor of every element in𝐋\\mathbf\{L\}, and𝐋\\mathbf\{L\}d\-separates𝐍\\mathbf\{N\}and𝐂\\mathbf\{C\}\.
- \(iii\)There does not exist a partition of𝐋\\mathbf\{L\}such that all elements in the partition are latent atomic covers\.

Having introduced the required notions, we now provide the structural assumptions considered byHuang et al\. \([2022](https://arxiv.org/html/2605.20396#bib.bib38)\)\.

###### Assumption 3\(Identifiable linear latent hierarchical graph\(Huang et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib38)\)\)\.

A graph𝒢\\mathcal\{G\}is an identifiable linear latent hierarchical graph if \(i\) every latent variableLiL\_\{i\}belongs to at least one latent atomic cover and there is no triangle structure in the graph, and \(ii\) if there exists a set of variables𝐕\\mathbf\{V\}such that every variable in𝐕\\mathbf\{V\}is a collider of two latent atomic covers𝐋1\\mathbf\{L\}\_\{1\},𝐋2\\mathbf\{L\}\_\{2\}, and denote by𝐓\\mathbf\{T\}the minimal set of variables that d\-separates𝐋𝟏\\mathbf\{L\_\{1\}\}from𝐋𝟐\\mathbf\{L\_\{2\}\}, then we must have\|𝐕\|\+\|𝐓\|≥\|𝐋1\|\+\|𝐋2\|\|\\mathbf\{V\}\|\+\|\\mathbf\{T\}\|\\geq\|\\mathbf\{L\}\_\{1\}\|\+\|\\mathbf\{L\}\_\{2\}\|\.

The assumption above requires that each latent variable belongs to at least one latent atomic cover, since a latent atomic cover is the minimal identifiable substructure in a graph using rank constraints of covariance over observed variables\. Also, the assumption requires certain graphical patterns that are related to the common children across different latent atomic covers for the identifiability of the whole latent structure\. The assumption above may be more general than Assumption[2](https://arxiv.org/html/2605.20396#Thmassumption2)as it allows each measured variable to have multiple latent variables together as parents, and also allow some latent variables to not have any measured child at all; see[SectionA\.4](https://arxiv.org/html/2605.20396#A1.SS4)for more details\. An example is given in[Figure2](https://arxiv.org/html/2605.20396#S5.F2)\. Under the assumption above,Huang et al\. \([2022](https://arxiv.org/html/2605.20396#bib.bib38)\)developed a constraint\-based method based on rank deficiency test to estimate the equivalence class of the true structure\.

In this section, we develop a score\-based method under[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)\. We characterize the degrees of freedom in[Section5\.1](https://arxiv.org/html/2605.20396#S5.SS1)and provide an exact search method in[Section5\.2](https://arxiv.org/html/2605.20396#S5.SS2)\. We do not provide a continuous search method for this structural assumption, since it cannot be straightforwardly formulated as inequality constraints, similar to[Equation5](https://arxiv.org/html/2605.20396#S4.E5)\.

![Refer to caption](https://arxiv.org/html/2605.20396v1/x2.png)Figure 2:Example of latent hierarchical structure\.### 5\.1Degrees of Freedom

As noted in[Sections3\.2](https://arxiv.org/html/2605.20396#S3.SS2)and[4](https://arxiv.org/html/2605.20396#S4), a key ingredient of the score\-based method is the specification of the degrees of freedom\. It may be natural to expect that the degrees of freedom are equal to the number of edges and measured variables, similar to the structural assumption considered in[Section4](https://arxiv.org/html/2605.20396#S4)and the standard setting without latent variables\(Chickering,[2002](https://arxiv.org/html/2605.20396#bib.bib19)\)\. However, this property does not hold for latent hierarchical structures, as illustrated by the following lemma\.

###### Proposition 3\.

Suppose that DAG𝒢\\mathcal\{G\}follows the linear latent variable causal model in[Equation1](https://arxiv.org/html/2605.20396#S2.E1)\. Suppose also that there existk≥2k\\geq 2latent variables in𝒢\\mathcal\{G\}with the same set of parents and children, where either the number of parents or children is at leastkk\. Then,dim\(𝒢\)≤\|𝒢\|\+m−k​\(k−1\)/2\\dim\(\\mathcal\{G\}\)\\leq\|\\mathcal\{G\}\|\+m\-k\(k\-1\)/2\.

The proof is given in[SectionB\.6](https://arxiv.org/html/2605.20396#A2.SS6)\. As shown in the proof, under these circumstances, there exists an alternative structure𝒢~\\tilde\{\\mathcal\{G\}\}obtained by removingk​\(k−1\)/2k\(k\-1\)/2edges \(corresponding to edges involving the parents or children\) from𝒢\\mathcal\{G\}such that𝒢~\\tilde\{\\mathcal\{G\}\}leads to the same distribution set as𝒢\\mathcal\{G\}, i\.e\.,ℳ​\(𝒢~\)=ℳ​\(𝒢\)\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\}\)=\\mathcal\{M\}\(\\mathcal\{G\}\)\. Thus, the degrees of freedom in this scenario are reduced compared to the number of edges and measured variables\. This reduction can be intuitively explained by the redundancy of certain edges in such situations\. For instance, the degrees of freedom for the structure are4444instead of4646\(i\.e\., the sum of number of edges and measured variables\), because the variables\{L2,L3\}\\\{L\_\{2\},L\_\{3\}\\\}and\{L4,L5\}\\\{L\_\{4\},L\_\{5\}\\\}\(i\.e\., in the same atomic covers\) have the same parents and children\.

Algorithm 2Degrees of freedom under[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)Input:Structure

𝒢\\mathcal\{G\}
Output:Degrees of freedom

dd
Initialize degrees of freedom

d←\|𝒢\|\+md\\leftarrow\|\\mathcal\{G\}\|\+m;

foreachsubset

𝐋\\mathbf\{L\}of latent variables where

\|𝐋\|≥2\|\\mathbf\{L\}\|\\geq 2do

ifvariables in

𝐋\\mathbf\{L\}have the same parents and childrenandno proper superset of

𝐋\\mathbf\{L\}satisfies the previous conditionthen

d←d−\|𝐋\|​\(\|𝐋\|−1\)/2d\\leftarrow d\-\|\\mathbf\{L\}\|\(\|\\mathbf\{L\}\|\-1\)/2;

returndegrees of freedom

dd

Building on the result above, we develop a procedure in[Algorithm2](https://arxiv.org/html/2605.20396#alg2)to calculate the degrees of freedom under[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)\. Specifically, the algorithm iterates over all subsets of latent variables and calculate the degrees of freedom that can be reduced\. The following proposition shows that the algorithm outputs the upper bound of the degrees of freedom, with a proof provided in[SectionB\.7](https://arxiv.org/html/2605.20396#A2.SS7)\.

###### Proposition 4\(Degrees of freedom\)\.

Suppose that DAG𝒢\\mathcal\{G\}satisfies[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)\. Then,[Algorithm2](https://arxiv.org/html/2605.20396#alg2)outputs the upper bound ofdim⁡\(𝒢\)\\operatorname\{dim\}\(\\mathcal\{G\}\)\.

We conjecture, supported by simulations over10,00010,000examples \(by computing the rank of Jacobian matrices\(Geiger et al\.,[2001](https://arxiv.org/html/2605.20396#bib.bib33)\)\) and the experiments in[Section6](https://arxiv.org/html/2605.20396#S6), that the upper bound provided by this algorithm is tight under operator𝒪min​\(𝒪skeleton​\(⋅\)\)\\mathcal\{O\}\_\{\\textrm\{min\}\}\(\\mathcal\{O\}\_\{\\textrm\{skeleton\}\}\(\\cdot\)\), although a proof seems to involve tools from algebraic statistics and is not straightforward\. For instance,Drton et al\. \([2023](https://arxiv.org/html/2605.20396#bib.bib30)\)analyzed the degrees of freedom for sparse factor analysis, which is technically complex even with independent latent variables\.

### 5\.2Consistency and Exact Score\-Based Search

With the algorithm to compute the degrees of freedom, we now develop a score\-based method to estimate latent hierarchical structures\. We first establish the correctness of the score\-based method\. Under[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)and the generalized faithfulness assumption, we prove that the structure with the optimal score is Markov equivalent to the true hierarcahical structure, up to certain rank equivalent graph operators\. The proof and definition of the operators together with illustrative examples can be found in[SectionB\.8](https://arxiv.org/html/2605.20396#A2.SS8)\.

###### Theorem 3\(Correctness\)\.

Suppose that the true DAG𝒢∗\\mathcal\{G\}^\{\*\}and the distributionΣX\\Sigma\_\{X\}satisfy the generalized faithfulness assumption, and that𝒢∗\\mathcal\{G\}^\{\*\}satisfies[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)\. Let𝒢^\\hat\{\\mathcal\{G\}\}be a global minimizer of the following optimization problem:

min𝒢∈𝔾m\\displaystyle\\min\_\{\\mathcal\{G\}\\in\\mathbb\{G\}^\{m\}\}scoredim⁡\(𝒢,𝐃\)\\displaystyle\\quad\\operatorname\{score\}\_\{\\textrm\{dim\}\}\(\\mathcal\{G\},\\mathbf\{D\}\)subject​to\\displaystyle\\mathrm\{subject\\ to\}𝒢​satisfies[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3),\\displaystyle\\quad\\mathcal\{G\}\\text\{ satisfies \\lx@cref\{creftype~refnum\}\{assumption:llh\_graphical\_criterion\}\},𝒢=𝒪min​\(𝒪skeleton​\(𝒢\)\)\.\\displaystyle\\quad\\mathcal\{G\}=\\mathcal\{O\}\_\{\\textrm\{min\}\}\(\\mathcal\{O\}\_\{\\textrm\{skeleton\}\}\(\\mathcal\{G\}\)\)\.Then,𝒪atomic​\(𝒢^\)\\mathcal\{O\}\_\{\\textrm\{atomic\}\}\(\\hat\{\\mathcal\{G\}\}\)and𝒪atomic​\(𝒪min​\(𝒪skeleton​\(𝒢∗\)\)\)\\mathcal\{O\}\_\{\\textrm\{atomic\}\}\(\\mathcal\{O\}\_\{\\textrm\{min\}\}\(\\mathcal\{O\}\_\{\\textrm\{skeleton\}\}\(\\mathcal\{G\}^\{\*\}\)\)\)are Markov equivalent in the large sample limit\.

Algorithm 3Enumerating structures under[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)Input:Measured variables

X1,…,XmX\_\{1\},\\dots,X\_\{m\}
Output:Set of DAGs

𝐀\\mathbf\{A\}
Initialize

𝐀\\mathbf\{A\}as an empty set;

for

n=1n=1to

n¯\\bar\{n\}do

forpartition

\{𝐂i\}i=1l\\\{\\mathbf\{C\}\_\{i\}\\\}\_\{i=1\}^\{l\}of

\{Li\}i=1n\\\{L\_\{i\}\\\}\_\{i=1\}^\{n\}as atomic coversdo

forDAG

𝒢𝐂→𝐂\\mathcal\{G\}\_\{\\mathbf\{C\}\\rightarrow\\mathbf\{C\}\}among

\{𝐂i\}i=1l\\\{\\mathbf\{C\}\_\{i\}\\\}\_\{i=1\}^\{l\}do

forDAG

𝒢𝐂→X\\mathcal\{G\}\_\{\\mathbf\{C\}\\rightarrow X\}from

\{𝐂i\}i=1l\\\{\\mathbf\{C\}\_\{i\}\\\}\_\{i=1\}^\{l\}to

\{Xi\}i=1m\\\{X\_\{i\}\\\}\_\{i=1\}^\{m\}do

Construct DAG

𝒢\\mathcal\{G\}by combining

𝒢𝐂→𝐂\\mathcal\{G\}\_\{\\mathbf\{C\}\\rightarrow\\mathbf\{C\}\}and

𝒢𝐂→X\\mathcal\{G\}\_\{\\mathbf\{C\}\\rightarrow X\};

if

𝒢\\mathcal\{G\}satisfies[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)and is not Markov equivalent to all DAGs in

𝐀\\mathbf\{A\}then

𝐀←𝐀∪\{𝒢\}\\mathbf\{A\}\\leftarrow\\mathbf\{A\}\\cup\\\{\\mathcal\{G\}\\\};

returnset

𝐀\\mathbf\{A\}

Table 1:F1 scores of skeletons across various structural assumptions and sample sizes\. For each setting, the top two methods are in bold\. For FOFC, the number within the brackets indicates the number of valid runs \(for which an error did not occur\)\.Based on the theorem above, we develop an exact search procedure for structures under[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)\. We introduce a procedure in[Algorithm3](https://arxiv.org/html/2605.20396#alg3)for enumeration of these structures, wheren¯\\bar\{n\}is a hyperparameter indicating the maximal number of latent variables\. Note that a possible upper bound forn¯\\bar\{n\}is3​m3m\. Similar to the algorithm developed in[Algorithm1](https://arxiv.org/html/2605.20396#alg1), we enumerate only structures that are not Markov equivalent to one another, leveraging the score equivalence property\. The whole procedure of[Algorithm3](https://arxiv.org/html/2605.20396#alg3)is roughly as follows\. We maintain a set of DAGs,𝐀\\mathbf\{A\}\. Given the number of observed variables, we first decide the possible number of latent variables, and then enumerate all possible combinations of atomic covers\. For each combination of atomic covers, we enumerate all possible DAGs among atomic covers and all possible DAGs from atomic covers to observed variables\. Finally, we combine both enumerated DAGs to get a possible graph𝒢\\mathcal\{G\}; if𝒢\\mathcal\{G\}satisfies[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)and is not Markov equivalent to all structures in𝐀\\mathbf\{A\}, we add it to𝐀\\mathbf\{A\}\. Once the search space is constructed, the algorithm identifies the structure with the optimal score, with the degrees of freedom for each structure computed using[Algorithm2](https://arxiv.org/html/2605.20396#alg2)\.

## 6Experiments

We conduct experiments to validate our score\-based methods, by comparing them to existing methods that support causally\-related latent variables, such as FOFC\(Kummerfeld & Ramsey,[2016](https://arxiv.org/html/2605.20396#bib.bib44)\), HUANG\(Huang et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib38)\), and GIN\(Xie et al\.,[2020](https://arxiv.org/html/2605.20396#bib.bib83)\)\. We do not include the comparison with FCI because, even when working perfectly, it will output complete PAGs over the observed variables for most ground truths considered here, which do not have any information of the orientation and are not informative\. Moreover, we denote our exact search method by SALAD and continuous one by SALAD\-CS, and adopt the BIC score here\.

For the ground truths, we consider the 1\-factor models and hierarchical structures provided in[Figures4](https://arxiv.org/html/2605.20396#A4.F4)and[5](https://arxiv.org/html/2605.20396#A4.F5)in[AppendixC](https://arxiv.org/html/2605.20396#A3)\. For each structure, the nonzero elements of matricesBBandCCare generated uniformly at random from the interval\[−2,−0\.5\]∪\[0\.5,2\.0\]\[\-2,\-0\.5\]\\cup\[0\.5,2\.0\]\. For GIN, the noise termsEXE\_\{X\}andELE\_\{L\}followUniform⁡\[−α,α\]\\operatorname\{Uniform\}\[\-\\alpha,\\alpha\], whereα\\alphais sampled uniformly from\[6,15\]\[\\sqrt\{6\},\\sqrt\{15\}\]\. For the other methods, the noise terms follow Gaussian distributions with variances sampled uniformly from interval\[2,5\]\[2,5\]\. We consider sample sizeT∈\{100,300,1000,3000,10000\}T\\in\\\{100,300,1000,3000,10000\\\}\. We evaluate the estimated structures using F1 scores calculated over the skeleton and structural Hamming distance \(SHD\) over the MEC\. We run three random trials for each ground truth, and report the mean and standard devation for each metric\. Further details about the metrics and baselines can be found in[AppendixC](https://arxiv.org/html/2605.20396#A3)\.

The F1 scores of skeletons are reported in[Table1](https://arxiv.org/html/2605.20396#S5.T1), while the SHDs of MECs are given in[Table2](https://arxiv.org/html/2605.20396#A4.T2)in the supplementary material\. One observes that our methods achieve much better F1 scores and SHDs as compared to the other baselines, especially for small sample sizes\. For instance, for100100samples, our SALAD method achieves average F1 scores of0\.990\.99and0\.920\.92for 1\-factor model and hierarchical structures, respectively, while the second best baseline achieves F1 scores of0\.900\.90and0\.570\.57, respectively\. Note that although FOFC achieves an F1 score of0\.900\.90for 1\-factor model in this case, four of the runs are not valid \(i\.e\., an error occurred\)\. A possible reason of the improvement is that the existing constraint\-based baselines, as discussed in[Section1](https://arxiv.org/html/2605.20396#S1), may be prone to the issue of error propagation during the estimation procedure, while our score\-based method is not susceptible to such an issue\. Furthermore, the F1 scores of our method are close to one for both structural assumptions when the sample size is large, which suggest that BIC may be a valid scoring function in our setting and help verify the correctness established in[Theorems2](https://arxiv.org/html/2605.20396#Thmtheorem2)and[3](https://arxiv.org/html/2605.20396#Thmtheorem3)\. The runtime and computational efficiency are discussed in[AppendixD](https://arxiv.org/html/2605.20396#A4)\.

## 7Conclusion and Discussion

In this work, we propose SALAD, a score\-based causal discovery method capable of identifying causal relations among latent variables\. Achieving score equivalence and consistency, along with degrees of freedom characterization and exact and continuous score\-based methods, our work provides a unified view on multiple existing constraint\-based methods with latent variables, and further validates the effectiveness of score\-based methods\. We hope that this work could spur future research on developing score\-based methods for latent variable causal models\.

Indeed, our exact methods require a relatively long runtime, similar to the exact score\-based methods even without latent variables\(Singh & Moore,[2005](https://arxiv.org/html/2605.20396#bib.bib69); Yuan & Malone,[2013](https://arxiv.org/html/2605.20396#bib.bib84)\)\. Future works include developing greedy approaches similar to GES to make the search procedure more efficient and scalable, and studying the theoretical justifications of using BIC score under the structural assumptions considered\.

## Acknowledgements

The authors would like to thank the anonymous reviewers for helpful comments and suggestions\. The authors would also like to acknowledge the support from NSF Grant 2229881, the National Institutes of Health \(NIH\) under Contract R01HL159805, and grants from Apple Inc\., KDDI Research Inc\., Quris AI, and Florin Court Capital\.

## Impact Statement

This paper presents work whose goal is to advance the field of Machine Learning\. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here\.

## References

- Adams et al\. \(2021\)Adams, J\., Hansen, N\., and Zhang, K\.Identification of partially observed linear causal models: Graphical conditions for the non\-gaussian and heterogeneous cases\.*Advances in Neural Information Processing Systems*, 34:22822–22833, 2021\.
- Akbari et al\. \(2021\)Akbari, S\., Mokhtarian, E\., Ghassami, A\., and Kiyavash, N\.Recursive causal structure learning in the presence of latent variables and selection bias\.*Advances in Neural Information Processing Systems*, 34:10119–10130, 2021\.
- Amendola et al\. \(2020\)Amendola, C\., Dettling, P\., Drton, M\., Onori, F\., and Wu, J\.Structure learning for cyclic linear causal models\.In*Conference on Uncertainty in Artificial Intelligence*, 2020\.
- Améndola et al\. \(2023\)Améndola, C\., Drton, M\., Grosdos, A\., Homs, R\., and Robeva, E\.Third\-order moment varieties of linear non\-gaussian graphical models\.*Information and Inference: A Journal of the IMA*, 12\(3\):iaad007, 2023\.
- Anandkumar et al\. \(2013\)Anandkumar, A\., Hsu, D\., Javanmard, A\., and Kakade, S\.Learning linear bayesian networks with latent variables\.In*International Conference on Machine Learning*, pp\. 249–257\. PMLR, 2013\.
- Bellot & van der Schaar \(2021\)Bellot, A\. and van der Schaar, M\.Deconfounded score method: Scoring dags with dense unobserved confounding\.*arXiv preprint arXiv:2103\.15106*, 2021\.
- Benedetti & Risler \(1990\)Benedetti, R\. and Risler, J\.\-J\.*Real algebraic and semi\-algebraic sets*\.Actualités mathématiques\. Hermann, Paris, 1990\.
- Bernstein et al\. \(2020\)Bernstein, D\., Saeed, B\., Squires, C\., and Uhler, C\.Ordering\-based causal structure learning in the presence of latent variables\.In*International Conference on Artificial Intelligence and Statistics*, pp\. 4098–4108\. PMLR, 2020\.
- Bertsekas \(1982\)Bertsekas, D\. P\.*Constrained Optimization and Lagrange Multiplier Methods*\.Academic Press, 1982\.
- Bertsekas \(1999\)Bertsekas, D\. P\.*Nonlinear Programming*\.Athena Scientific, 2nd edition, 1999\.
- Bhattacharya et al\. \(2021\)Bhattacharya, R\., Nagarajan, T\., Malinsky, D\., and Shpitser, I\.Differentiable causal discovery under unmeasured confounding\.In*International Conference on Artificial Intelligence and Statistics*, 2021\.
- Bollen \(1989\)Bollen, K\. A\.*The General Model, Part I: Latent Variable and Measurement Models Combined*, chapter Eight, pp\. 319–394\.John Wiley & Sons, Ltd, 1989\.ISBN 9781118619179\.
- Brito & Pearl \(2002\)Brito, C\. and Pearl, J\.A new identification condition for recursive models with correlated errors\.*Structural Equation Modeling: A Multidisciplinary Journal*, 9\(4\):459–474, 2002\.doi:10\.1207/S15328007SEM0904“˙1\.
- Brouillard et al\. \(2020\)Brouillard, P\., Lachapelle, S\., Lacoste, A\., Lacoste\-Julien, S\., and Drouin, A\.Differentiable causal discovery from interventional data\.In*Advances in Neural Information Processing Systems*, 2020\.
- Byrd et al\. \(1995\)Byrd, R\. H\., Lu, P\., Nocedal, J\., and Zhu, C\.A limited memory algorithm for bound constrained optimization\.*SIAM Journal on Scientific Computing*, 16\(5\):1190–1208, 1995\.
- Byrne \(2001\)Byrne, B\.*Structural Equation Modeling With AMOS: Basic Concepts, Applications, and Programming*\.Multivariate Applications Series\. Taylor & Francis, 2001\.
- Cai et al\. \(2019\)Cai, R\., Xie, F\., Glymour, C\., Hao, Z\., and Zhang, K\.Triad constraints for learning causal structure of latent variables\.*Advances in neural information processing systems*, 32, 2019\.
- Chen et al\. \(2022\)Chen, Z\., Xie, F\., Qiao, J\., Hao, Z\., Zhang, K\., and Cai, R\.Identification of linear latent variable model with arbitrary distribution\.In*Proceedings of the AAAI Conference on Artificial Intelligence*, volume 36, pp\. 6350–6357, 2022\.
- Chickering \(2002\)Chickering, D\. M\.Optimal structure identification with greedy search\.*Journal of Machine Learning Research*, 3\(Nov\):507–554, 2002\.
- Choi et al\. \(2011\)Choi, M\. J\., Tan, V\. Y\., Anandkumar, A\., and Willsky, A\. S\.Learning latent tree graphical models\.*Journal of Machine Learning Research*, 12:1771–1812, 2011\.
- Claassen & Bucur \(2022\)Claassen, T\. and Bucur, I\. G\.Greedy equivalence search in the presence of latent confounders\.In*Conference on Uncertainty in Artificial Intelligence*, 2022\.
- Claassen et al\. \(2013\)Claassen, T\., Mooij, J\., and Heskes, T\.Learning sparse causal models is not np\-hard\.*arXiv preprint arXiv:1309\.6824*, 2013\.
- Colombo et al\. \(2012\)Colombo, D\., Maathuis, M\. H\., Kalisch, M\., and Richardson, T\. S\.Learning high\-dimensional directed acyclic graphs with latent and selection variables\.*The Annals of Statistics*, pp\. 294–321, 2012\.
- Cox et al\. \(2015\)Cox, D\. A\., Little, J\., and O’Shea, D\.*Ideals, Varieties, and Algorithms*\.Springer, New York, fourth edition, 2015\.
- Cui et al\. \(2018\)Cui, R\., Groot, P\., Schauer, M\., and Heskes, T\.Learning the causal structure of copula models with latent variables\.2018\.
- Dai et al\. \(2022\)Dai, H\., Spirtes, P\., and Zhang, K\.Independence testing\-based approach to causal discovery under measurement error and linear non\-gaussian models\.*Advances in Neural Information Processing Systems*, 35:27524–27536, 2022\.
- Dempster et al\. \(1977\)Dempster, A\. P\., Laird, N\. M\., and Rubin, D\. B\.Maximum likelihood from incomplete data via the EM algorithm\.*Journal of the Royal Statistical Society: Series B*, 39:1–38, 1977\.
- Dong et al\. \(2023\)Dong, X\., Huang, B\., Ng, I\., Song, X\., Zheng, Y\., Jin, S\., Legaspi, R\., Spirtes, P\., and Zhang, K\.A versatile causal discovery framework to allow causally\-related hidden variables\.*arXiv preprint arXiv:2312\.11001*, 2023\.
- Drton \(2018\)Drton, M\.Algebraic problems in structural equation modeling\.In*Advanced Studies in Pure Mathematics*, pp\. 35–86\. Mathematical Society of Japan, 2018\.
- Drton et al\. \(2023\)Drton, M\., Grosdos, A\., Portakal, I\., and Sturma, N\.Algebraic sparse factor analysis\.*arXiv preprint arXiv:2312\.14762*, 2023\.
- Forster et al\. \(2017\)Forster, M\., Raskutti, G\., Stern, R\., and Weinberger, N\.The frugal inference of causal relations\.*The British Journal for the Philosophy of Science*, 69, 04 2017\.
- Geiger et al\. \(1996\)Geiger, D\., Heckerman, D\. E\., and Meek, C\.Asymptotic model selection for directed networks with hidden variables\.In*Conference on Uncertainty in Artificial Intelligence*, 1996\.
- Geiger et al\. \(2001\)Geiger, D\., Heckerman, D\., King, H\., and Meek, C\.Stratified exponential families: Graphical models and model selection\.*The Annals of Statistics*, 29\(2\):505–529, 2001\.
- Ghassami et al\. \(2020\)Ghassami, A\., Yang, A\., Kiyavash, N\., and Zhang, K\.Characterizing distribution equivalence and structure learning for cyclic and acyclic directed graphs\.In*International Conference on Machine Learning*, 2020\.
- Goldberg \(1992\)Goldberg, L\.The development of markers for the big five factor structure\.*Psychological Assessment*, 4:26–42, 03 1992\.
- Haughton \(1988\)Haughton, D\. M\. A\.On the choice of a model to fit data from an exponential family\.*The Annals of Statistics*, 16\(1\):342–355, 1988\.
- Himi et al\. \(2019\)Himi, S\. A\., Buehner, M\., Schwaighofer, M\., Klapetek, A\., and Hilbert, S\.Multitasking behavior and its related constructs: Executive functions, working memory capacity, relational integration, and divided attention\.*Cognition*, 189:275–298, 08 2019\.
- Huang et al\. \(2022\)Huang, B\., Low, C\., Xie, F\., Glymour, C\., and Zhang, K\.Latent hierarchical causal structure discovery with rank constraints\.In*Advances in Neural Information Processing Systems*, 2022\.
- Hyvärinen et al\. \(2023\)Hyvärinen, A\., Khemakhem, I\., and Monti, R\.Identifiability of latent\-variable and structural\-equation models: from linear to nonlinear\.*Annals of the Institute of Statistical Mathematics*, 2023\.
- Jang et al\. \(2017\)Jang, E\., Gu, S\., and Poole, B\.Categorical reparameterization with gumbel\-softmax\.In*International Conference on Learning Representations*, 2017\.
- Kingma & Ba \(2014\)Kingma, D\. and Ba, J\.Adam: A method for stochastic optimization\.In*International Conference on Learning Representations*, 2014\.
- Kivva et al\. \(2021\)Kivva, B\., Rajendran, G\., Ravikumar, P\., and Aragam, B\.Learning latent causal graphs via mixture oracles\.*Advances in Neural Information Processing Systems*, 34:18087–18101, 2021\.
- Koller & Friedman \(2009\)Koller, D\. and Friedman, N\.*Probabilistic Graphical Models: Principles and Techniques*\.MIT Press, Cambridge, MA, 2009\.
- Kummerfeld & Ramsey \(2016\)Kummerfeld, E\. and Ramsey, J\.Causal clustering for 1\-factor measurement models\.In*Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining*, pp\. 1655–1664, 2016\.
- Leung et al\. \(2015\)Leung, D\., Drton, M\., and Hara, H\.Identifiability of directed Gaussian graphical models with one latent source\.*Electronic Journal of Statistics*, 10, 05 2015\.
- Maddison et al\. \(2017\)Maddison, C\. J\., Mnih, A\., and Teh, Y\. W\.The concrete distribution: A continuous relaxation of discrete random variables\.In*International Conference on Learning Representations*, 2017\.
- Nandy et al\. \(2018\)Nandy, P\., Hauser, A\., and Maathuis, M\. H\.High\-dimensional consistency in score\-based and hybrid structure learning\.*The Annals of Statistics*, 46\(6A\):3151–3183, 2018\.
- Ng et al\. \(2020\)Ng, I\., Ghassami, A\., and Zhang, K\.On the role of sparsity and DAG constraints for learning linear DAGs\.In*Advances in Neural Information Processing Systems*, 2020\.
- Ng et al\. \(2022\)Ng, I\., Zhu, S\., Fang, Z\., Li, H\., Chen, Z\., and Wang, J\.Masked gradient\-based causal structure learning\.In*SIAM International Conference on Data Mining*, 2022\.
- Ng et al\. \(2024\)Ng, I\., Huang, B\., and Zhang, K\.Structure learning with continuous optimization: A sober look and beyond\.In*Proceedings of the Third Conference on Causal Learning and Reasoning*, 2024\.
- Nocedal & Wright \(2006\)Nocedal, J\. and Wright, S\. J\.*Numerical optimization*\.Springer series in operations research and financial engineering\. Springer, 2nd edition, 2006\.
- Nowzohour et al\. \(2017\)Nowzohour, C\., Maathuis, M\. H\., Evans, R\. J\., and Bühlmann, P\.Distributional equivalence and structure learning for bow\-free acyclic path diagrams\.2017\.
- Paszke et al\. \(2019\)Paszke, A\., Gross, S\., Massa, F\., Lerer, A\., Bradbury, J\., Chanan, G\., Killeen, T\., Lin, Z\., Gimelshein, N\., Antiga, L\., Desmaison, A\., Kopf, A\., Yang, E\., DeVito, Z\., Raison, M\., Tejani, A\., Chilamkurthy, S\., Steiner, B\., Fang, L\., Bai, J\., and Chintala, S\.PyTorch: An imperative style, high\-performance deep learning library\.In*Advances in Neural Information Processing Systems*, 2019\.
- Pearl \(2009\)Pearl, J\.*Causality*\.Cambridge university press, 2009\.
- Ramsey et al\. \(2017\)Ramsey, J\., Glymour, M\., Sanchez\-Romero, R\., and Glymour, C\.A million variables and more: the fast greedy equivalence search algorithm for learning high\-dimensional graphical causal models, with an application to functional magnetic resonance images\.*International Journal of Data Science and Analytics*, 3\(2\):121–129, 2017\.
- Raskutti & Uhler \(2014\)Raskutti, G\. and Uhler, C\.Learning directed acyclic graphs based on sparsest permutations\.*arXiv preprint arXiv:1307\.0366v3*, 2014\.
- Richardson & Spirtes \(2002\)Richardson, T\. and Spirtes, P\.Ancestral graph Markov models\.*The Annals of Statistics*, 30\(4\):962–1030, 2002\.
- Richardson \(1996\)Richardson, T\. S\.*Models of feedback: interpretation and discovery*\.PhD thesis, Carnegie\-Mellon University, 1996\.
- Salehkaleybar et al\. \(2020\)Salehkaleybar, S\., Ghassami, A\., Kiyavash, N\., and Zhang, K\.Learning linear non\-gaussian causal models in the presence of latent variables\.*The Journal of Machine Learning Research*, 21\(1\):1436–1459, 2020\.
- Scheines et al\. \(1998\)Scheines, R\., Spirtes, P\., Glymour, C\., Meek, C\., and Richardson, T\.The TETRAD project: Constraint based aids to causal model specification\.*Multivariate Behavioral Research*, 33:65–117, 1998\.
- Schölkopf et al\. \(2021\)Schölkopf, B\., Locatello, F\., Bauer, S\., Ke, N\. R\., Kalchbrenner, N\., Goyal, A\., and Bengio, Y\.Towards causal representation learning\.*Proceedings of the IEEE*, 109\(5\):612–634, 2021\.
- Schwarz \(1978\)Schwarz, G\.Estimating the dimension of a model\.*The Annals of Statistics*, 6\(2\):461–464, 1978\.
- Shahin & Chechik \(2020\)Shahin, R\. and Chechik, M\.Automatic and efficient variability\-aware lifting of functional programs\.*Proceedings of the ACM on Programming Languages*, 4\(OOPSLA\):1–27, 2020\.
- Shimizu et al\. \(2009\)Shimizu, S\., Hoyer, P\. O\., and Hyvärinen, A\.Estimation of linear non\-gaussian acyclic models for latent factors\.*Neurocomputing*, 72\(7\-9\):2024–2027, 2009\.
- Shpitser et al\. \(2012\)Shpitser, I\., Richardson, T\. S\., Robins, J\. M\., and Evans, R\.Parameter and structure learning in nested Markov models\.*arXiv preprint arXiv:1207\.5058*, 2012\.
- Silva & Scheines \(2005\)Silva, R\. and Scheines, R\.Generalized measurement models\.Technical report, Carnegie\-Mellon Univ Pittsburgh PA School of Computer Science, 2005\.
- Silva et al\. \(2003\)Silva, R\., Scheines, R\., Glymour, C\., and Spirtes, P\.Learning measurement models for unobserved variables\.In*Conference on Uncertainty in Artificial Intelligence*, 2003\.
- Silva et al\. \(2006\)Silva, R\., Scheines, R\., Glymour, C\., and Spirtes, P\.Learning the structure of linear latent variable models\.*Journal of Machine Learning Research*, 7\(8\):191–246, 2006\.URL[http://jmlr\.org/papers/v7/silva06a\.html](http://jmlr.org/papers/v7/silva06a.html)\.
- Singh & Moore \(2005\)Singh, A\. P\. and Moore, A\. W\.Finding optimal Bayesian networks by dynamic programming\.Technical report, Carnegie Mellon University, 2005\.
- Spirtes \(2010\)Spirtes, P\.Introduction to causal inference\.*Journal of Machine Learning Research*, 11\(5\), 2010\.
- Spirtes & Glymour \(1991\)Spirtes, P\. and Glymour, C\.An algorithm for fast recovery of sparse causal graphs\.*Social Science Computer Review*, 9:62–72, 1991\.
- Spirtes et al\. \(2001\)Spirtes, P\., Glymour, C\., and Scheines, R\.*Causation, Prediction, and Search*\.MIT press, 2nd edition, 2001\.
- Spirtes et al\. \(2013\)Spirtes, P\. L\., Meek, C\., and Richardson, T\. S\.Causal inference in the presence of latent variables and selection bias\.*arXiv preprint arXiv:1302\.4983*, 2013\.
- Squires et al\. \(2023\)Squires, C\., Seigal, A\., Bhate, S\., and Uhler, C\.Linear causal disentanglement via interventions\.In*International Conference on Machine Learning*, 2023\.
- Sturma et al\. \(2023\)Sturma, N\., Squires, C\., Drton, M\., and Uhler, C\.Unpaired multi\-domain causal representation learning\.*arXiv preprint arXiv:2302\.00993*, 2023\.
- Sullivant et al\. \(2010\)Sullivant, S\., Talaska, K\., and Draisma, J\.Trek separation for gaussian graphical models\.*The Annals of Statistics*, 38\(3\):1665–1685, 2010\.
- Triantafillou & Tsamardinos \(2016\)Triantafillou, S\. and Tsamardinos, I\.Score\-based vs constraint\-based causal learning in the presence of confounders\.In*Cfa@ uai*, pp\. 59–67, 2016\.
- van Ommen & Mooij \(2017\)van Ommen, T\. and Mooij, J\. M\.Algebraic equivalence of linear structural equation models\.In*Conference on Uncertainty in Artificial Intelligence*, 2017\.
- Verma & Pearl \(1991\)Verma, T\. and Pearl, J\.Equivalence and synthesis of causal models\.In*Conference on Uncertainty in Artificial Intelligence*, 1991\.
- Virtanen et al\. \(2020\)Virtanen, P\., Gommers, R\., Oliphant, T\. E\., Haberland, M\., Reddy, T\., Cournapeau, D\., Burovski, E\., Peterson, P\., Weckesser, W\., Bright, J\., van der Walt, S\. J\., Brett, M\., Wilson, J\., Jarrod Millman, K\., Mayorov, N\., Nelson, A\. R\. J\., Jones, E\., Kern, R\., Larson, E\., Carey, C\., Polat, İ\., Feng, Y\., Moore, E\. W\., Vand erPlas, J\., Laxalde, D\., Perktold, J\., Cimrman, R\., Henriksen, I\., Quintero, E\. A\., Harris, C\. R\., Archibald, A\. M\., Ribeiro, A\. H\., Pedregosa, F\., van Mulbregt, P\., and Contributors, S\. \. \.SciPy 1\.0: Fundamental algorithms for scientific computing in Python\.*Nature Methods*, 17:261–272, 2020\.
- Vowels et al\. \(2022\)Vowels, M\. J\., Camgoz, N\. C\., and Bowden, R\.D’ya like DAGs? a survey on structure learning and causal discovery\.*ACM Computing Surveys*, 55\(4\), nov 2022\.ISSN 0360\-0300\.
- Wang & Drton \(2023\)Wang, Y\. S\. and Drton, M\.Causal discovery with unobserved confounding and non\-Gaussian data\.*Journal of Machine Learning Research*, 24\(271\):1–61, 2023\.
- Xie et al\. \(2020\)Xie, F\., Cai, R\., Huang, B\., Glymour, C\., Hao, Z\., and Zhang, K\.Generalized independent noise condition for estimating latent variable causal graphs\.*Advances in neural information processing systems*, 33:14891–14902, 2020\.
- Yuan & Malone \(2013\)Yuan, C\. and Malone, B\.Learning optimal Bayesian networks: A shortest path perspective\.*Journal of Artificial Intelligence Research*, 48\(1\):23–65, 2013\.
- Zeng et al\. \(2021\)Zeng, Y\., Shimizu, S\., Cai, R\., Xie, F\., Yamamoto, M\., and Hao, Z\.Causal discovery with multi\-domain LiNGAM for latent factors\.In*Causal Analysis Workshop Series*, pp\. 1–4\. PMLR, 2021\.
- Zhang \(2008\)Zhang, J\.On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias\.*Artificial Intelligence*, 172\(16\-17\):1873–1896, 2008\.
- Zhang et al\. \(2018\)Zhang, K\., Gong, M\., Ramsey, J\., Batmanghelich, K\., Spirtes, P\., and Glymour, C\.Causal discovery with linear non\-Gaussian models under measurement error: Structural identifiability results\.In*Conference on Uncertainty in Artificial Intelligence*, 2018\.
- Zhang et al\. \(2024\)Zhang, K\., Xie, S\., Ng, I\., and Zheng, Y\.Causal representation learning from multiple distributions: A general setting\.*arXiv preprint arXiv:2402\.05052*, 2024\.
- Zhang \(2004\)Zhang, N\. L\.Hierarchical latent class models for cluster analysis\.*Journal of Machine Learning Research*, 5:697–723, dec 2004\.
- Zheng et al\. \(2018\)Zheng, X\., Aragam, B\., Ravikumar, P\., and Xing, E\. P\.DAGs with NO TEARS: Continuous optimization for structure learning\.In*Advances in Neural Information Processing Systems*, 2018\.

Supplementary Material

## Appendix AFurther Discussions

We provide supplementary discussions below as complements to various sections in the main paper\.

### A\.1Latent Variable Causal Models

In real\-world scenarios, one may often encounter the latent variable causal model in[Equation1](https://arxiv.org/html/2605.20396#S2.E1), where the measured variables do not influence each other and are effects of latent variables\. Thus, there have been many works that aim to estimate this type of linear latent variable causal models\(Silva et al\.,[2003](https://arxiv.org/html/2605.20396#bib.bib67),[2006](https://arxiv.org/html/2605.20396#bib.bib68); Silva & Scheines,[2005](https://arxiv.org/html/2605.20396#bib.bib66); Zhang,[2004](https://arxiv.org/html/2605.20396#bib.bib89); Choi et al\.,[2011](https://arxiv.org/html/2605.20396#bib.bib20); Kummerfeld & Ramsey,[2016](https://arxiv.org/html/2605.20396#bib.bib44); Cai et al\.,[2019](https://arxiv.org/html/2605.20396#bib.bib17); Xie et al\.,[2020](https://arxiv.org/html/2605.20396#bib.bib83); Dai et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib26); Huang et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib38); Chen et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib18)\)\.

To provide some examples, in psychometrics, multiple questions are often used as indirect proxies for each latent personality dimension \(e\.g\., openness, extraversion, self\-esteem\)\(Goldberg,[1992](https://arxiv.org/html/2605.20396#bib.bib35); Byrne,[2001](https://arxiv.org/html/2605.20396#bib.bib16); Himi et al\.,[2019](https://arxiv.org/html/2605.20396#bib.bib37)\), forming a latent variable causal model\. When analyzing fMRI data, a large number of voxels are measured, which do not necessarily have clear semantic meanings\. A hierarchical structure can then be used to model functionally meaningful brain regions at different levels\(Huang et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib38)\)\. In representation learning, recent works typically assumed that measured variables \(e\.g\., image pixels\) are effects of latent variables and that there are no direct causal influences among the measured variables\(Schölkopf et al\.,[2021](https://arxiv.org/html/2605.20396#bib.bib61); Hyvärinen et al\.,[2023](https://arxiv.org/html/2605.20396#bib.bib39); Zhang et al\.,[2024](https://arxiv.org/html/2605.20396#bib.bib88)\)\.

### A\.2Generalized Faithfulness Assumption

We discuss the necessity of the generalized faithfulness assumption adopted in our results\. One of the advantages of score\-based causal discovery is that it typically relies on the sparsest Markov representation \(SMR\) assumption \(or unique frugality assumption\)\(Forster et al\.,[2017](https://arxiv.org/html/2605.20396#bib.bib31); Raskutti & Uhler,[2014](https://arxiv.org/html/2605.20396#bib.bib56)\), which is strictly weaker than the faithfulness assumption\(Spirtes et al\.,[2001](https://arxiv.org/html/2605.20396#bib.bib72)\)in the setting without latent variables\. In our setting with latent variables, it is possible to modify Theorems[1](https://arxiv.org/html/2605.20396#Thmtheorem1),[2](https://arxiv.org/html/2605.20396#Thmtheorem2), and[3](https://arxiv.org/html/2605.20396#Thmtheorem3)to replace generalized faithfulness with a formulation similar to the SMR assumption\. However, doing so may not be informative because \(1\) SMR may not be strictly weaker than the faithfulness assumption in our setting, and \(2\) simply assuming SMR in our setting may not provide insights into what structural assumptions the true structure should obey\. Thus, we adopt the generalized faithfulness assumption and various structural assumptions to make the results more informative\.

### A\.3Algebraic Equivalence

We provide a further discussion of algebraic equivalence\(van Ommen & Mooij,[2017](https://arxiv.org/html/2605.20396#bib.bib78)\)as a complement to[Section3\.3](https://arxiv.org/html/2605.20396#S3.SS3)\. First note that two algebraic equivalent structures are not necessarily Markov equivalent\. The reason is that, without any restriction on the structures, one may construct different structures that entail the same equality constraints\. For instance, consider the structure in[Figure5](https://arxiv.org/html/2605.20396#A4.F5)\(a\), denoted as𝒢1\\mathcal\{G\}\_\{1\}, and another structure𝒢2\\mathcal\{G\}\_\{2\}that is identical to𝒢1\\mathcal\{G\}\_\{1\}, except that the edgesL1→X1L\_\{1\}\\rightarrow X\_\{1\},L2→X1L\_\{2\}\\rightarrow X\_\{1\}, andL1→X2L\_\{1\}\\rightarrow X\_\{2\}are removed in𝒢2\\mathcal\{G\}\_\{2\}\. One can show𝒢1\\mathcal\{G\}\_\{1\}and𝒢2\\mathcal\{G\}\_\{2\}are algebraic equivalent, but clearly they are not Markov equivalent\.

Nonetheless, algebraic equivalence may be a reasonable way for estimating linear latent variable causal models, because equality constraints \(of the covariance matrices\) are some of the major footprints in the data that one could leverage \(without considering higher\-order statistics\) to identify the underlying structures\. This can be done by relating these constraints to the structures via the generalized faithfulness assumption\.

### A\.4Structural Assumptions

We discuss the similarities and differences between the structural assumptions considered in our work\. First, both[Assumptions2](https://arxiv.org/html/2605.20396#Thmassumption2)and[3](https://arxiv.org/html/2605.20396#Thmassumption3)require that the observed variables are leaf nodes, and that there are no direct causal influences among observed variables\. The key differences between them are as follows\. \(1\)[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)requires that each latent variable has at least three measured variables as its children, while[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)allows latent variables to form a hierarchical structure \- some latent variables may only have latent variables as their children\. \(2\)[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)requires each observed variable to be caused by a single latent variable, while[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)allows an observed variable to be caused by a group of latent variables\. Since[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)does not require each latent variable to have measured variables as children, the structure among latent variables cannot be arbitrary, and thus there is a tradeoff between[Assumptions2](https://arxiv.org/html/2605.20396#Thmassumption2)and[3](https://arxiv.org/html/2605.20396#Thmassumption3)\.

## Appendix BProofs

### B\.1Proof of[Lemma1](https://arxiv.org/html/2605.20396#Thmlemma1)

See[1](https://arxiv.org/html/2605.20396#Thmlemma1)

###### Proof\.

LetB~≔B​ΩL12\\tilde\{B\}\\coloneqq B\\Omega\_\{L\}^\{\\frac\{1\}\{2\}\}andC~≔ΩL−12​C​ΩL12\\tilde\{C\}\\coloneqq\\Omega\_\{L\}^\{\-\\frac\{1\}\{2\}\}C\\Omega\_\{L\}^\{\\frac\{1\}\{2\}\}\. We havesupp⁡\(B\)=supp⁡\(B~\)\\operatorname\{supp\}\(B\)=\\operatorname\{supp\}\(\\tilde\{B\}\),supp⁡\(C\)=supp⁡\(C~\)\\operatorname\{supp\}\(C\)=\\operatorname\{supp\}\(\\tilde\{C\}\), and

ΣX\\displaystyle\\Sigma\_\{X\}=B​\(I−C\)−1​ΩL​\(I−C\)−⊤​B⊤\+ΩX\\displaystyle=B\(I\-C\)^\{\-1\}\\Omega\_\{L\}\(I\-C\)^\{\-\\top\}B^\{\\top\}\+\\Omega\_\{X\}=B​ΩL12​ΩL−12​\(I−C\)−1​ΩL12​ΩL12​\(I−C\)−⊤​ΩL−12​ΩL12​B⊤\+ΩX\\displaystyle=B\\Omega\_\{L\}^\{\\frac\{1\}\{2\}\}\\Omega\_\{L\}^\{\-\\frac\{1\}\{2\}\}\(I\-C\)^\{\-1\}\\Omega\_\{L\}^\{\\frac\{1\}\{2\}\}\\Omega\_\{L\}^\{\\frac\{1\}\{2\}\}\(I\-C\)^\{\-\\top\}\\Omega\_\{L\}^\{\-\\frac\{1\}\{2\}\}\\Omega\_\{L\}^\{\\frac\{1\}\{2\}\}B^\{\\top\}\+\\Omega\_\{X\}=\(B​ΩL12\)​\(I−ΩL−12​C​ΩL12\)−1​\(I−ΩL−12​C​ΩL12\)−⊤​\(B​ΩL12\)⊤\+ΩX\\displaystyle=\(B\\Omega\_\{L\}^\{\\frac\{1\}\{2\}\}\)\(I\-\\Omega\_\{L\}^\{\-\\frac\{1\}\{2\}\}C\\Omega\_\{L\}^\{\\frac\{1\}\{2\}\}\)^\{\-1\}\(I\-\\Omega\_\{L\}^\{\-\\frac\{1\}\{2\}\}C\\Omega\_\{L\}^\{\\frac\{1\}\{2\}\}\)^\{\-\\top\}\(B\\Omega\_\{L\}^\{\\frac\{1\}\{2\}\}\)^\{\\top\}\+\\Omega\_\{X\}=B~​\(I−C~\)−1​\(I−C~\)−⊤​B~⊤\+ΩX∎\\displaystyle=\\tilde\{B\}\(I\-\\tilde\{C\}\)^\{\-1\}\(I\-\\tilde\{C\}\)^\{\-\\top\}\\tilde\{B\}^\{\\top\}\+\\Omega\_\{X\}\\qed

### B\.2Proof of[Proposition1](https://arxiv.org/html/2605.20396#Thmproposition1)

The following proof is partly inspired by that of the score equivalence property in the setting without latent variables\(Koller & Friedman,[2009](https://arxiv.org/html/2605.20396#bib.bib43)\)\.

See[1](https://arxiv.org/html/2605.20396#Thmproposition1)

###### Proof\.

Because structures𝒢1\\mathcal\{G\}\_\{1\}and𝒢2\\mathcal\{G\}\_\{2\}are Markov equivalent, they can generate the same set of covariance matrices over variablesXXandLL\. Thus, for any parametersB,C,ΩXB,C,\\Omega\_\{X\}of𝒢1\\mathcal\{G\}\_\{1\}withΩL=I\\Omega\_\{L\}=I, there exists parametersB′,C′,ΩX′,ΩL′B^\{\\prime\},C^\{\\prime\},\\Omega\_\{X\}^\{\\prime\},\\Omega\_\{L\}^\{\\prime\}of𝒢2\\mathcal\{G\}\_\{2\}that can generate the same covariance matrix overXXandLL, which imply thatB′,C′,ΩX′,ΩL′B^\{\\prime\},C^\{\\prime\},\\Omega\_\{X\}^\{\\prime\},\\Omega\_\{L\}^\{\\prime\}can generate the same covariance matrix overXX\. By[Lemma1](https://arxiv.org/html/2605.20396#Thmlemma1), there exists parameters of𝒢2\\mathcal\{G\}\_\{2\}, denoted asB~,C~,Ω~X=ΩX′\\tilde\{B\},\\tilde\{C\},\\tilde\{\\Omega\}\_\{X\}=\\Omega\_\{X\}^\{\\prime\}andΩ~L=I\\tilde\{\\Omega\}\_\{L\}=I, that can generate the covariance matrix\. Note that the likelihood function depends only on the covariance matrix, which indicates

scoreℒ⁡\(𝒢1,𝐃\)=ℒ​\(B^,C^,Ω^X;𝐃\)=ℒ​\(B~,C~,Ω~X;𝐃\)≥scoreℒ⁡\(𝒢2,𝐃\),\\operatorname\{score\}\_\{\\mathcal\{L\}\}\(\\mathcal\{G\}\_\{1\},\\mathbf\{D\}\)=\\mathcal\{L\}\(\\hat\{B\},\\hat\{C\},\\hat\{\\Omega\}\_\{X\};\\mathbf\{D\}\)=\\mathcal\{L\}\(\\tilde\{B\},\\tilde\{C\},\\tilde\{\\Omega\}\_\{X\};\\mathbf\{D\}\)\\geq\\operatorname\{score\}\_\{\\mathcal\{L\}\}\(\\mathcal\{G\}\_\{2\},\\mathbf\{D\}\),whereB^,C^,Ω^X\\hat\{B\},\\hat\{C\},\\hat\{\\Omega\}\_\{X\}are the solutions of the optimization problem in[Equation3](https://arxiv.org/html/2605.20396#S3.E3)for𝒢=𝒢1\\mathcal\{G\}=\\mathcal\{G\}\_\{1\}, and, as described above,B~,C~,Ω~X\\tilde\{B\},\\tilde\{C\},\\tilde\{\\Omega\}\_\{X\}are the corresponding parameters of structure𝒢2\\mathcal\{G\}\_\{2\}\.

Similarly, the same reasoning impliesscoreℒ⁡\(𝒢2,𝐃\)≥scoreℒ⁡\(𝒢1,𝐃\)\\operatorname\{score\}\_\{\\mathcal\{L\}\}\(\\mathcal\{G\}\_\{2\},\\mathbf\{D\}\)\\geq\\operatorname\{score\}\_\{\\mathcal\{L\}\}\(\\mathcal\{G\}\_\{1\},\\mathbf\{D\}\)\. Combining both cases, we havescoreℒ⁡\(𝒢1,𝐃\)=scoreℒ⁡\(𝒢2,𝐃\)\\operatorname\{score\}\_\{\\mathcal\{L\}\}\(\\mathcal\{G\}\_\{1\},\\mathbf\{D\}\)=\\operatorname\{score\}\_\{\\mathcal\{L\}\}\(\\mathcal\{G\}\_\{2\},\\mathbf\{D\}\)\. Furthermore, since𝒢1\\mathcal\{G\}\_\{1\}and𝒢2\\mathcal\{G\}\_\{2\}can generate the same set of covariance matrices over variablesXXandLL, they can generate the same set of covariance matrices over variablesXX\. This impliesdim\(𝒢1\)=dim\(𝒢2\)\\dim\(\\mathcal\{G\}\_\{1\}\)=\\dim\(\\mathcal\{G\}\_\{2\}\)\. Therefore, we havescoreBIC⁡\(𝒢1,𝐃\)=scoreBIC⁡\(𝒢2,𝐃\)\\operatorname\{score\}\_\{\\textrm\{BIC\}\}\(\\mathcal\{G\}\_\{1\},\\mathbf\{D\}\)=\\operatorname\{score\}\_\{\\textrm\{BIC\}\}\(\\mathcal\{G\}\_\{2\},\\mathbf\{D\}\)andscoredim⁡\(𝒢1,𝐃\)=scoredim⁡\(𝒢2,𝐃\)\\operatorname\{score\}\_\{\\textrm\{dim\}\}\(\\mathcal\{G\}\_\{1\},\\mathbf\{D\}\)=\\operatorname\{score\}\_\{\\textrm\{dim\}\}\(\\mathcal\{G\}\_\{2\},\\mathbf\{D\}\)\. ∎

### B\.3Proof of[Theorem1](https://arxiv.org/html/2605.20396#Thmtheorem1)

The overall proof strategy below is partly inspired by the proof ofGhassami et al\. \([2020](https://arxiv.org/html/2605.20396#bib.bib34), Theorem 3\)\. See[1](https://arxiv.org/html/2605.20396#Thmtheorem1)

###### Proof\.

Since the search space contains the true DAG𝒢∗\\mathcal\{G\}^\{\*\}that can generateΣX\\Sigma\_\{X\}in the large sample limit, the estimated DAG𝒢^\\hat\{\\mathcal\{G\}\}can also generateΣX\\Sigma\_\{X\}, because otherwise its score will be infinity and will not be a solution of the optimization problem\. Therefore,ΣX\\Sigma\_\{X\}belongs to the distribution set of𝒢^\\hat\{\\mathcal\{G\}\}, i\.e\.,ΣX∈ℳ​\(𝒢^\)\\Sigma\_\{X\}\\in\\mathcal\{M\}\(\\hat\{\\mathcal\{G\}\}\), which implies thatΣX\\Sigma\_\{X\}contains all the equality and inequality constraints of𝒢^\\hat\{\\mathcal\{G\}\}\. Under the generalized faithfulness assumption, we have

H​\(𝒢^\)⊆H​\(𝒢∗\)\.H\(\\hat\{\\mathcal\{G\}\}\)\\subseteq H\(\\mathcal\{G\}^\{\*\}\)\.\(6\)Now suppose by contradiction thatH​\(𝒢^\)⊊H​\(𝒢∗\)H\(\\hat\{\\mathcal\{G\}\}\)\\subsetneq H\(\\mathcal\{G\}^\{\*\}\)\. This impliesdim\(𝒢^\)\>dim\(𝒢∗\)\\dim\(\\hat\{\\mathcal\{G\}\}\)\>\\dim\(\\mathcal\{G\}^\{\*\}\), which is a contradiction because the objective function impliesdim\(𝒢^\)≤dim\(𝒢∗\)\\dim\(\\hat\{\\mathcal\{G\}\}\)\\leq\\dim\(\\mathcal\{G\}^\{\*\}\)\. Thus, we obtain

H​\(𝒢^\)⊊̸H​\(𝒢∗\)\.H\(\\hat\{\\mathcal\{G\}\}\)\\not\\subsetneq H\(\\mathcal\{G\}^\{\*\}\)\.\(7\)By[Equations6](https://arxiv.org/html/2605.20396#A2.E6)and[7](https://arxiv.org/html/2605.20396#A2.E7), we haveH​\(𝒢^\)=H​\(𝒢∗\)H\(\\hat\{\\mathcal\{G\}\}\)=H\(\\mathcal\{G\}^\{\*\}\)\. ∎

### B\.4Proof of[Proposition2](https://arxiv.org/html/2605.20396#Thmproposition2)

We first state the following lemma adapted fromLeung et al\. \([2015](https://arxiv.org/html/2605.20396#bib.bib45)\)that relates the parameter identifiability from a given structure to the underlying degrees of freedom\.

###### Lemma 2\(Leung et al\. \([2015](https://arxiv.org/html/2605.20396#bib.bib45)\)\)\.

Supposef:𝐒→ℝdf:\\mathbf\{S\}\\rightarrow\\mathbb\{R\}^\{d\}is a polynomial map defined on an open set𝐒⊆ℝp\\mathbf\{S\}\\subseteq\\mathbb\{R\}^\{p\}\. The following statements are equivalent:

1. \(i\)ffis generically finite\-to\-one\.
2. \(ii\)The Jacobian matrix offfis generically of full column rank\.

We now provide the proof of the following proposition\. See[2](https://arxiv.org/html/2605.20396#Thmproposition2)

###### Proof\.

By[Corollary1](https://arxiv.org/html/2605.20396#Thmcorollary1), it suffices to consider the case whereΩL=I\\Omega\_\{L\}=I\. Since the structure𝒢\\mathcal\{G\}satisfies[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2), byBollen \([1989](https://arxiv.org/html/2605.20396#bib.bib12)\), the corresponding parametersB,CB,CandΩX\\Omega\_\{X\}of𝒢\\mathcal\{G\}are identifiable fromΣX\\Sigma\_\{X\}up to certain indeterminacy\. Specifically,BBis identifiable up to column permutations and sign changes,CCis identifiable up to equal row and column permutations, andΩX\\Omega\_\{X\}is identifiable\. Therefore, the map fromB,CB,C, andΩX\\Omega\_\{X\}toΣX\\Sigma\_\{X\}is generically finite\-to\-one\.

The map fromB,CB,C, andΩX\\Omega\_\{X\}\(withΩL=I\\Omega\_\{L\}=I\) toΣX\\Sigma\_\{X\}is a polynomial map\. By[Lemma2](https://arxiv.org/html/2605.20396#Thmlemma2), the Jacobian matrix of such map is generically of full column rank\. Note that the degrees of freedom \(or dimension\) of a polynomial map are equal to the maximal rank of the corresponding Jacobian matrix\(Geiger et al\.,[2001](https://arxiv.org/html/2605.20396#bib.bib33), Theorem 10\)\. Therefore, the degrees of freedom are equal to the number of parameters inB,CB,C, andΩX\\Omega\_\{X\}, i\.e\.,dim⁡\(𝒢\)=‖B𝒢‖0\+‖C𝒢‖0\+m=\|𝒢\|\+m\\operatorname\{dim\}\(\\mathcal\{G\}\)=\\\|B\_\{\\mathcal\{G\}\}\\\|\_\{0\}\+\\\|C\_\{\\mathcal\{G\}\}\\\|\_\{0\}\+m=\|\\mathcal\{G\}\|\+m\. ∎

### B\.5Proof of[Theorem2](https://arxiv.org/html/2605.20396#Thmtheorem2)

See[2](https://arxiv.org/html/2605.20396#Thmtheorem2)

###### Proof\.

Since the search space contains the true DAG𝒢∗\\mathcal\{G\}^\{\*\}that can generateΣX\\Sigma\_\{X\}in the large sample limit, the estimated DAG𝒢^\\hat\{\\mathcal\{G\}\}can also generateΣX\\Sigma\_\{X\}, because otherwise its score will be infinity and will not be a solution of the optimization problem\. Because𝒢∗\\mathcal\{G\}^\{\*\}andΣX\\Sigma\_\{X\}satisfy the generalized faithfulness assumption, we haveH​\(𝒢^\)=H​\(𝒢∗\)H\(\\hat\{\\mathcal\{G\}\}\)=H\(\\mathcal\{G\}^\{\*\}\)in the large sample limit by[Proposition2](https://arxiv.org/html/2605.20396#Thmproposition2)and restricting the set of structures to those satisfying[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)in[Theorem1](https://arxiv.org/html/2605.20396#Thmtheorem1)\. This indicates that𝒢^\\hat\{\\mathcal\{G\}\}andΣX\\Sigma\_\{X\}also satisfy the generalized faithfulness assumption\.

Since𝒢^\\hat\{\\mathcal\{G\}\}and𝒢∗\\mathcal\{G\}^\{\*\}satisfy[Assumption2](https://arxiv.org/html/2605.20396#Thmassumption2)and both can faithfully generate the covariance matrixΣX\\Sigma\_\{X\}, we have:

- •BySilva et al\. \([2003](https://arxiv.org/html/2605.20396#bib.bib67), Corollary 1\), the measurement models of𝒢^\\hat\{\\mathcal\{G\}\}and𝒢∗\\mathcal\{G\}^\{\*\}are identical \(up to relabeling of latent variables\)\. In other words, the columns ofB𝒢^B\_\{\\hat\{\\mathcal\{G\}\}\}are a permutation of the columns ofB𝒢∗B\_\{\\mathcal\{G\}^\{\*\}\}\.
- •With a correct measurement model, by leveraging the transitivity of Markov equivalence, it follows straightforwardly fromSilva et al\. \([2006](https://arxiv.org/html/2605.20396#bib.bib68), Theorems 20\)that the structural models \(i\.e\., the subgraphs over all and only the latent variables\) of𝒢^\\hat\{\\mathcal\{G\}\}and𝒢∗\\mathcal\{G\}^\{\*\}are Markov equivalent \(up to relabeling of latent variables\)\.

Combining the reasoning for both measurement model and structural model,𝒢^\\hat\{\\mathcal\{G\}\}and𝒢∗\\mathcal\{G\}^\{\*\}are Markov equivalent \(up to relabeling of latent variables\)\. ∎

### B\.6Proof of[Proposition3](https://arxiv.org/html/2605.20396#Thmproposition3)

For structure𝒢\\mathcal\{G\}, we define the following distribution set with the constraintΩL=I\\Omega\_\{L\}=I:

ℳ​\(𝒢;ΩL=I\)≔\{B​\(I−C\)−1​\(I−C\)−⊤​B⊤\+ΩX:supp⁡\(B\)⊆supp⁡\(B𝒢\),supp⁡\(C\)⊆supp⁡\(C𝒢\),ΩX∈diag⁡\(ℝ\>0m\)\}\.\\mathcal\{M\}\(\\mathcal\{G\};\\Omega\_\{L\}=I\)\\coloneqq\\\{B\(I\-C\)^\{\-1\}\(I\-C\)^\{\-\\top\}B^\{\\top\}\+\\Omega\_\{X\}:\\operatorname\{supp\}\(B\)\\subseteq\\operatorname\{supp\}\(B\_\{\\mathcal\{G\}\}\),\\operatorname\{supp\}\(C\)\\subseteq\\operatorname\{supp\}\(C\_\{\\mathcal\{G\}\}\),\\Omega\_\{X\}\\in\\operatorname\{diag\}\(\\mathbb\{R\}\_\{\>0\}^\{m\}\)\\\}\.Note that the dimension of the domain and the image space are upper bounds for the dimension ofℳ​\(𝒢;ΩL=I\)\\mathcal\{M\}\(\\mathcal\{G\};\\Omega\_\{L\}=I\)\. By[Lemma1](https://arxiv.org/html/2605.20396#Thmlemma1), it is straightforward to obtain the following corollary\.

###### Corollary 1\.

For any structure𝒢\\mathcal\{G\}satisfying[Equation1](https://arxiv.org/html/2605.20396#S2.E1), we have

ℳ​\(𝒢\)=ℳ​\(𝒢;ΩL=I\)​and​dim\(𝒢\)≤min⁡\(\|𝒢\|\+m,12​m​\(m\+1\)\)\.\\mathcal\{M\}\(\\mathcal\{G\}\)=\\mathcal\{M\}\(\\mathcal\{G\};\\Omega\_\{L\}=I\)\\text\{\\qquad and\\qquad\}\\dim\(\\mathcal\{G\}\)\\leq\\min\\left\(\|\\mathcal\{G\}\|\+m,\\frac\{1\}\{2\}m\(m\+1\)\\right\)\.

Therefore, it suffices to analyze the degrees of freedom forℳ​\(𝒢;ΩL=I\)\\mathcal\{M\}\(\\mathcal\{G\};\\Omega\_\{L\}=I\)instead ofℳ​\(𝒢\)\\mathcal\{M\}\(\\mathcal\{G\}\)\.

We first provide the following lemma which shows that an appropriate orthogonal transformation ofBBandCCcan generate the same covariance matrix\.

###### Lemma 3\(Orthogonal transformation\)\.

Consider any set of parametersB,C,ΩXB,C,\\Omega\_\{X\}, andΣX\\Sigma\_\{X\}that satisfy

ΣX=B​\(I−C\)−1​\(I−C\)−⊤​B⊤\+ΩX\.\\Sigma\_\{X\}=B\(I\-C\)^\{\-1\}\(I\-C\)^\{\-\\top\}B^\{\\top\}\+\\Omega\_\{X\}\.For any orthogonal matrixQQ, i\.e\.,Q​Q⊤=IQQ^\{\\top\}=I\. the parametersB~=B​Q\\tilde\{B\}=BQandC~=Q⊤​C​Q\\tilde\{C\}=Q^\{\\top\}CQalso satisfy

ΣX=B~​\(I−C~\)−1​\(I−C~\)−⊤​B~⊤\+ΩX\.\\Sigma\_\{X\}=\\tilde\{B\}\(I\-\\tilde\{C\}\)^\{\-1\}\(I\-\\tilde\{C\}\)^\{\-\\top\}\\tilde\{B\}^\{\\top\}\+\\Omega\_\{X\}\.

###### Proof\.

The proof follows from straightforward algebraic manipulations:

ΣX\\displaystyle\\Sigma\_\{X\}=B​\(I−C\)−1​\(I−C\)−⊤​B⊤\+ΩX\\displaystyle=B\(I\-C\)^\{\-1\}\(I\-C\)^\{\-\\top\}B^\{\\top\}\+\\Omega\_\{X\}=B​Q​Q⊤​\(I−C\)−1​Q​Q⊤​\(I−C\)−⊤​Q​Q⊤​B⊤\+ΩX\\displaystyle=BQQ^\{\\top\}\(I\-C\)^\{\-1\}QQ^\{\\top\}\(I\-C\)^\{\-\\top\}QQ^\{\\top\}B^\{\\top\}\+\\Omega\_\{X\}=\(B​Q\)​\(I−Q⊤​C​Q\)−1​\(I−Q⊤​C​Q\)−⊤​\(B​Q\)⊤\+ΩX\\displaystyle=\(BQ\)\(I\-Q^\{\\top\}CQ\)^\{\-1\}\(I\-Q^\{\\top\}CQ\)^\{\-\\top\}\(BQ\)^\{\\top\}\+\\Omega\_\{X\}=B~​\(I−C~\)−1​\(I−C~\)−⊤​B~⊤\+ΩX\.∎\\displaystyle=\\tilde\{B\}\(I\-\\tilde\{C\}\)^\{\-1\}\(I\-\\tilde\{C\}\)^\{\-\\top\}\\tilde\{B\}^\{\\top\}\+\\Omega\_\{X\}\.\\qed

The following result shows that, in specific cases, some of the edges can be removed from the structure while still leading to the same distribution set\.

###### Lemma 4\.

Suppose that DAG𝒢\\mathcal\{G\}follows the linear latent variable causal model in[Equation1](https://arxiv.org/html/2605.20396#S2.E1)\. Suppose also that there existk≥2k\\geq 2latent variables𝐋\\mathbf\{L\}in𝒢\\mathcal\{G\}with the same set of parents and children, where the number of children \(parents\) is at leastkk\. Then, there exists a structure, denoted by𝒢~\\tilde\{\\mathcal\{G\}\}, such that: \(1\)𝒢~\\tilde\{\\mathcal\{G\}\}is identical to𝒢\\mathcal\{G\}, except thatk​\(k−1\)/2k\(k\-1\)/2edges among those latent variables𝐋\\mathbf\{L\}and their children \(parents\) are removed in𝒢~\\tilde\{\\mathcal\{G\}\}, and \(2\)ℳ​\(𝒢~\)=ℳ​\(𝒢\)\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\}\)=\\mathcal\{M\}\(\\mathcal\{G\}\)\.

###### Proof\.

Consider any set of parametersB,C,ΩXB,C,\\Omega\_\{X\}, andΣX\\Sigma\_\{X\}that satisfy

ΣX=B​\(I−C\)−1​\(I−C\)−⊤​B⊤\+ΩX\.\\Sigma\_\{X\}=B\(I\-C\)^\{\-1\}\(I\-C\)^\{\-\\top\}B^\{\\top\}\+\\Omega\_\{X\}\.\(8\)We first consider the case where the number of children is at leastkk\. Denote by𝐒\\mathbf\{S\}the set of indices of the latent variables in𝐋\\mathbf\{L\}\. Since the latent variables𝐋\\mathbf\{L\}have the same set of children, the rows that correspond to the nonzero entries in each column ofB:,𝐒B\_\{:,\\mathbf\{S\}\}andC:,𝐒C\_\{:,\\mathbf\{S\}\}are the same, which we denote by𝐑1\\mathbf\{R\}\_\{1\}and𝐑2\\mathbf\{R\}\_\{2\}forB:,𝐒B\_\{:,\\mathbf\{S\}\}andC:,𝐒C\_\{:,\\mathbf\{S\}\}, respectively\. LetD=\(B𝐑1,𝐒,C𝐑2,𝐒\)D=\(B\_\{\\mathbf\{R\}\_\{1\},\\mathbf\{S\}\},C\_\{\\mathbf\{R\}\_\{2\},\\mathbf\{S\}\}\)be a matrix by concatenating the rows ofB𝐑1,𝐒B\_\{\\mathbf\{R\}\_\{1\},\\mathbf\{S\}\}andC𝐑2,𝐒C\_\{\\mathbf\{R\}\_\{2\},\\mathbf\{S\}\}; that is,DDis a matrix of dimension\(\|𝐑1\|\+\|𝐑2\|\)×k\(\|\\mathbf\{R\}\_\{1\}\|\+\|\\mathbf\{R\}\_\{2\}\|\)\\times k, where\|𝐑1\|\+\|𝐑2\|≥k\|\\mathbf\{R\}\_\{1\}\|\+\|\\mathbf\{R\}\_\{2\}\|\\geq k\. Applying orthogonal transformation as in the QR\-decomposition,DDcan be written asD=D~​QD=\\tilde\{D\}Q, whereD~\\tilde\{D\}is a lower\-triangular matrix andQQis an orthogonal matrix\. We rewrite the equation asD~=D​Q−1\\tilde\{D\}=DQ^\{\-1\}whereQ−1Q^\{\-1\}is also an orthogonal matrix\.

Consider the reversed mapping of the indices𝐀1≔\{1,2,…,\|𝐑1\|\}\\mathbf\{A\}\_\{1\}\\coloneqq\\\{1,2,\\dots,\|\\mathbf\{R\}\_\{1\}\|\\\}and𝐀2≔\{\|𝐑1\|\+1,\|𝐑1\|\+2,…,\|𝐑1\|\+\|𝐑2\|\}\\mathbf\{A\}\_\{2\}\\coloneqq\\\{\|\\mathbf\{R\}\_\{1\}\|\+1,\|\\mathbf\{R\}\_\{1\}\|\+2,\\dots,\|\\mathbf\{R\}\_\{1\}\|\+\|\\mathbf\{R\}\_\{2\}\|\\\}\. We now construct ann×nn\\times northogonal matrixUUas follows: \(1\)U𝐒,𝐒=Q−1U\_\{\\mathbf\{S\},\\mathbf\{S\}\}=Q^\{\-1\}, \(2\) the other non\-diagonal entries are zero, and \(3\) the other diagonal entries are one\. LetB~=B​U\\tilde\{B\}=BUandC~=U⊤​C​U\\tilde\{C\}=U^\{\\top\}CU\. Clearly, the entries inB~\\tilde\{B\}are the same asBB, except thatB𝐑1,𝐒B\_\{\\mathbf\{R\}\_\{1\},\\mathbf\{S\}\}is replaced withB𝐑1,𝐒​U𝐒,𝐒=B𝐑1,𝐒​Q−1=\(D​Q−1\)𝐀1,:=D~𝐀1,:B\_\{\\mathbf\{R\}\_\{1\},\\mathbf\{S\}\}U\_\{\\mathbf\{S\},\\mathbf\{S\}\}=B\_\{\\mathbf\{R\}\_\{1\},\\mathbf\{S\}\}Q^\{\-1\}=\(DQ^\{\-1\}\)\_\{\\mathbf\{A\}\_\{1\},:\}=\\tilde\{D\}\_\{\\mathbf\{A\}\_\{1\},:\}\. Similarly, the entries inC~\\tilde\{C\}are the same asCC, except that \(1\)C𝐑2,𝐒C\_\{\\mathbf\{R\}\_\{2\},\\mathbf\{S\}\}is replaced withC𝐑2,𝐒​U𝐒,𝐒=C𝐑2,𝐒​Q−1=\(D​Q−1\)𝐀2,:=D~𝐀2,:C\_\{\\mathbf\{R\}\_\{2\},\\mathbf\{S\}\}U\_\{\\mathbf\{S\},\\mathbf\{S\}\}=C\_\{\\mathbf\{R\}\_\{2\},\\mathbf\{S\}\}Q^\{\-1\}=\(DQ^\{\-1\}\)\_\{\\mathbf\{A\}\_\{2\},:\}=\\tilde\{D\}\_\{\\mathbf\{A\}\_\{2\},:\}, and \(2\)C𝐒,:C\_\{\\mathbf\{S\},:\}is replaced withU𝐒,𝐒⊤​C𝐒,:U\_\{\\mathbf\{S\},\\mathbf\{S\}\}^\{\\top\}C\_\{\\mathbf\{S\},:\}\. Since the latent variables𝐋\\mathbf\{L\}have the same set of parents, we havesupp⁡\(U𝐒,𝐒⊤​C𝐒,:\)⊆supp⁡\(C𝐒,:\)\\operatorname\{supp\}\(U\_\{\\mathbf\{S\},\\mathbf\{S\}\}^\{\\top\}C\_\{\\mathbf\{S\},:\}\)\\subseteq\\operatorname\{supp\}\(C\_\{\\mathbf\{S\},:\}\)\. This implies thatB~,C~\\tilde\{B\},\\tilde\{C\}, andΩX\\Omega\_\{X\}are parameterization of, e\.g\., structure𝒢~\\tilde\{\\mathcal\{G\}\}, where𝒢~\\tilde\{\\mathcal\{G\}\}has the same edges as𝒢\\mathcal\{G\}, except thatk​\(k−1\)/2k\(k\-1\)/2of the edges from𝒢\\mathcal\{G\}are removed in𝒢~\\tilde\{\\mathcal\{G\}\}\(which correspond to the “non\-lower\-triangular” entries fromDDthat become zero inD~\\tilde\{D\}after QR\-decomposition\)\. Clearly,𝒢~\\tilde\{\\mathcal\{G\}\}is identical to𝒢\\mathcal\{G\}, except thatk​\(k−1\)/2k\(k\-1\)/2edges among those latent variables𝐋\\mathbf\{L\}and their children are removed in𝒢~\\tilde\{\\mathcal\{G\}\}\. Furthermore, by[Lemma3](https://arxiv.org/html/2605.20396#Thmlemma3), the parametersB~\\tilde\{B\},C~\\tilde\{C\}, andΩX\\Omega\_\{X\}can generate the same covariance matrixΣX\\Sigma\_\{X\}\.

Since we are able to construct the same structure𝒢~\\tilde\{\\mathcal\{G\}\}using the above procedure for every parameterizationBB,CC, andΩX\\Omega\_\{X\}of𝒢\\mathcal\{G\}in[Equation8](https://arxiv.org/html/2605.20396#A2.E8), we haveℳ​\(𝒢~;ΩL=I\)=ℳ​\(𝒢;ΩL=I\)\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\};\\Omega\_\{L\}=I\)=\\mathcal\{M\}\(\\mathcal\{G\};\\Omega\_\{L\}=I\), which, by[Corollary1](https://arxiv.org/html/2605.20396#Thmcorollary1), impliesℳ​\(𝒢~\)=ℳ​\(𝒢\)\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\}\)=\\mathcal\{M\}\(\\mathcal\{G\}\)\. The same reasoning also applies when the number of parents is at leastkk, i\.e\., such a structure𝒢~\\tilde\{\\mathcal\{G\}\}can also be constructed\. ∎

We now provide the proof of the following proposition\. See[3](https://arxiv.org/html/2605.20396#Thmproposition3)

###### Proof\.

By[Lemma4](https://arxiv.org/html/2605.20396#Thmlemma4), there exists a structure, denoted by𝒢~\\tilde\{\\mathcal\{G\}\}, such that: \(1\)𝒢~\\tilde\{\\mathcal\{G\}\}is identical to𝒢\\mathcal\{G\}, except thatk​\(k−1\)/2k\(k\-1\)/2edges are removed in𝒢~\\tilde\{\\mathcal\{G\}\}, and \(2\)ℳ​\(𝒢~\)=ℳ​\(𝒢\)\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\}\)=\\mathcal\{M\}\(\\mathcal\{G\}\)\. By[Corollary1](https://arxiv.org/html/2605.20396#Thmcorollary1), this implies

dim\(𝒢\)=dim\(𝒢~\)≤\|𝒢~\|\+m=\|𝒢\|\+m−12​k​\(k−1\)\.\\displaystyle\\dim\(\\mathcal\{G\}\)=\\dim\(\\tilde\{\\mathcal\{G\}\}\)\\leq\|\\tilde\{\\mathcal\{G\}\}\|\+m=\|\\mathcal\{G\}\|\+m\-\\frac\{1\}\{2\}k\(k\-1\)\.∎

### B\.7Proof of[Proposition4](https://arxiv.org/html/2605.20396#Thmproposition4)

See[4](https://arxiv.org/html/2605.20396#Thmproposition4)

###### Proof\.

Let𝐋1,𝐋2,…,𝐋p\\mathbf\{L\}\_\{1\},\\mathbf\{L\}\_\{2\},\\dots,\\mathbf\{L\}\_\{p\}be the pairwise disjoint sets of latent variables in structure𝒢\\mathcal\{G\}such that \(1\) each𝐋i\\mathbf\{L\}\_\{i\}has at least two latent variables, \(2\) the variables of each𝐋i\\mathbf\{L\}\_\{i\}have the same set of parents and children in𝒢\\mathcal\{G\}, \(3\) no proper superset of each𝐋i\\mathbf\{L\}\_\{i\}has the same set of parents and children𝒢\\mathcal\{G\}\. Since the structure𝒢\\mathcal\{G\}is a DAG, we assume without loss of generality that𝐋1,𝐋2,…,𝐋p\\mathbf\{L\}\_\{1\},\\mathbf\{L\}\_\{2\},\\dots,\\mathbf\{L\}\_\{p\}are sorted based on the reversed causal ordering in𝒢\\mathcal\{G\}\. That is, variables𝐋i1\\mathbf\{L\}\_\{i\_\{1\}\}cannot be the ancestors of variables𝐋i2\\mathbf\{L\}\_\{i\_\{2\}\}in structure𝒢\\mathcal\{G\}fori1<i2i\_\{1\}<i\_\{2\}\.

Although[Algorithm2](https://arxiv.org/html/2605.20396#alg2)does not impose any specific order on the sets of latent variables, we suppose that the algorithm computes the dimension based on𝐋1,𝐋2,…,𝐋p\\mathbf\{L\}\_\{1\},\\mathbf\{L\}\_\{2\},\\dots,\\mathbf\{L\}\_\{p\}\(sorted according to reversed causal ordering\), which does not affect the computed degrees of freedom\. Denote bydjd\_\{j\}the output of the algorithm in thejj\-th iteration wherej∈\[p\]j\\in\[p\]\. The final output is thendpd\_\{p\}\. It suffices to showdim\(𝒢\)≤dj\\dim\(\\mathcal\{G\}\)\\leq d\_\{j\}in each iteration\.

We provide a proof by induction\. Specifically, we show that, for thejj\-th iteration, there exists a structure𝒢~j\\tilde\{\\mathcal\{G\}\}\_\{j\}such that:

1. 1\.𝒢~j\\tilde\{\\mathcal\{G\}\}\_\{j\}is identical to𝒢\\mathcal\{G\}, except that\|𝒢\|\+m−dj\|\\mathcal\{G\}\|\+m\-d\_\{j\}edges among the variables⋃i=1j𝐋i\\bigcup\_\{i=1\}^\{j\}\\mathbf\{L\}\_\{i\}and their children are removed in𝒢~j\\tilde\{\\mathcal\{G\}\}\_\{j\}\.
2. 2\.ℳ​\(𝒢~j\)=ℳ​\(𝒢\)\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\}\_\{j\}\)=\\mathcal\{M\}\(\\mathcal\{G\}\)\.

By[Corollary1](https://arxiv.org/html/2605.20396#Thmcorollary1), this implies the desired outcome

dim\(𝒢\)=dim\(𝒢~j\)≤\|𝒢~j\|\+m=\|𝒢\|\+m−\(\|𝒢\|\+m−dj\)=dj\.\\dim\(\\mathcal\{G\}\)=\\dim\(\\tilde\{\\mathcal\{G\}\}\_\{j\}\)\\leq\|\\tilde\{\\mathcal\{G\}\}\_\{j\}\|\+m=\|\\mathcal\{G\}\|\+m\-\(\|\\mathcal\{G\}\|\+m\-d\_\{j\}\)=d\_\{j\}\.
For induction, we first consider the base casej=1j=1\. By assumption, the variables𝐋1\\mathbf\{L\}\_\{1\}have the same set of parents and children in𝒢\\mathcal\{G\}, where, under[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3), the number of children is at least\|𝐋1\|\|\\mathbf\{L\}\_\{1\}\|\. By[Lemma4](https://arxiv.org/html/2605.20396#Thmlemma4), there exists a structure, denoted by𝒢~1\\tilde\{\\mathcal\{G\}\}\_\{1\}, such that: \(1\)𝒢~1\\tilde\{\\mathcal\{G\}\}\_\{1\}is identical to𝒢\\mathcal\{G\}, except that\|𝐋1\|​\(\|𝐋1\|−1\)/2=\|𝒢\|\+m−d1\|\\mathbf\{L\}\_\{1\}\|\(\|\\mathbf\{L\}\_\{1\}\|\-1\)/2=\|\\mathcal\{G\}\|\+m\-d\_\{1\}edges among the variables𝐋1\\mathbf\{L\}\_\{1\}and their children are removed in𝒢~1\\tilde\{\\mathcal\{G\}\}\_\{1\}, and \(2\)ℳ​\(𝒢~1\)=ℳ​\(𝒢\)\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\}\_\{1\}\)=\\mathcal\{M\}\(\\mathcal\{G\}\)\. Therefore, the base case is done\.

Suppose that the statements hold forj=tj=t, i\.e\., there exists a structure𝒢~t\\tilde\{\\mathcal\{G\}\}\_\{t\}such that: \(1\)𝒢~t\\tilde\{\\mathcal\{G\}\}\_\{t\}is identical to𝒢\\mathcal\{G\}, except that\|𝒢\|\+m−dt\|\\mathcal\{G\}\|\+m\-d\_\{t\}edges among the variables⋃i=1t𝐋i\\bigcup\_\{i=1\}^\{t\}\\mathbf\{L\}\_\{i\}and their children are removed in𝒢~t\\tilde\{\\mathcal\{G\}\}\_\{t\}, and \(2\)ℳ​\(𝒢~t\)=ℳ​\(𝒢\)\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\}\_\{t\}\)=\\mathcal\{M\}\(\\mathcal\{G\}\)\.

Now considerj=t\+1j=t\+1\. Note that only the edges among⋃i=1t𝐋i\\bigcup\_\{i=1\}^\{t\}\\mathbf\{L\}\_\{i\}and their children are removed in𝒢~t\\tilde\{\\mathcal\{G\}\}\_\{t\}\(as compared to𝒢\\mathcal\{G\}\); by assumption, the variables⋃i=1t𝐋i\\bigcup\_\{i=1\}^\{t\}\\mathbf\{L\}\_\{i\}are not ancestors of the variables𝐋t\+1\\mathbf\{L\}\_\{t\+1\}\. Therefore, the incoming and outgoing edges of variables𝐋t\+1\\mathbf\{L\}\_\{t\+1\}in𝒢~t\\tilde\{\\mathcal\{G\}\}\_\{t\}are the same as those in𝒢\\mathcal\{G\}\. This implies that𝐋t\+1\\mathbf\{L\}\_\{t\+1\}have the same set of parents and children in𝒢~t\\tilde\{\\mathcal\{G\}\}\_\{t\}, because they have the same set of parents and children in𝒢\\mathcal\{G\}\. Furthermore, under[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3), the number of children is at least\|𝐋t\+1\|\|\\mathbf\{L\}\_\{t\+1\}\|\. By[Lemma4](https://arxiv.org/html/2605.20396#Thmlemma4), there exists a structure, denoted by𝒢~t\+1\\tilde\{\\mathcal\{G\}\}\_\{t\+1\}, such that: \(1\)𝒢~t\+1\\tilde\{\\mathcal\{G\}\}\_\{t\+1\}is identical to𝒢~t\\tilde\{\\mathcal\{G\}\}\_\{t\}, except that\|𝐋t\+1\|​\(\|𝐋t\+1\|−1\)/2\|\\mathbf\{L\}\_\{t\+1\}\|\(\|\\mathbf\{L\}\_\{t\+1\}\|\-1\)/2edges among the variables𝐋t\+1\\mathbf\{L\}\_\{t\+1\}and their children are removed in𝒢~t\+1\\tilde\{\\mathcal\{G\}\}\_\{t\+1\}, and \(2\)ℳ​\(𝒢~t\+1\)=ℳ​\(𝒢~t\)\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\}\_\{t\+1\}\)=\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\}\_\{t\}\)\.

By the induction hypothesis, we haveℳ​\(𝒢~t\+1\)=ℳ​\(𝒢~t\)=ℳ​\(𝒢\)\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\}\_\{t\+1\}\)=\\mathcal\{M\}\(\\tilde\{\\mathcal\{G\}\}\_\{t\}\)=\\mathcal\{M\}\(\\mathcal\{G\}\)\. Also, it is clear that𝒢~t\+1\\tilde\{\\mathcal\{G\}\}\_\{t\+1\}is identical to𝒢\\mathcal\{G\}, except that

\|𝒢\|\+m−dt\+12​\|𝐋t\+1\|​\(\|𝐋t\+1\|−1\)=\|𝒢\|\+m−dt\+1\|\\mathcal\{G\}\|\+m\-d\_\{t\}\+\\frac\{1\}\{2\}\|\\mathbf\{L\}\_\{t\+1\}\|\(\|\\mathbf\{L\}\_\{t\+1\}\|\-1\)=\|\\mathcal\{G\}\|\+m\-d\_\{t\+1\}edges among the variables⋃i=1t\+1𝐋i\\bigcup\_\{i=1\}^\{t\+1\}\\mathbf\{L\}\_\{i\}and their children are removed in𝒢~t\+1\\tilde\{\\mathcal\{G\}\}\_\{t\+1\}\. Therefore, the induction step is done\. ∎

![Refer to caption](https://arxiv.org/html/2605.20396v1/figures/operation1.png)\(a\)Graph𝒢\\mathcal\{G\}\.
![Refer to caption](https://arxiv.org/html/2605.20396v1/figures/operation2.png)\(b\)𝒪skeleton​\(𝒢\)\\mathcal\{O\}\_\{\\text\{skeleton\}\}\(\\mathcal\{G\}\)\.
![Refer to caption](https://arxiv.org/html/2605.20396v1/figures/operation3.png)\(c\)𝒪min​\(𝒪skeleton​\(𝒢\)\)\\mathcal\{O\}\_\{\\text\{min\}\}\(\\mathcal\{O\}\_\{\\text\{skeleton\}\}\(\\mathcal\{G\}\)\)\.
![Refer to caption](https://arxiv.org/html/2605.20396v1/figures/operation4.png)\(d\)𝒪atomic​\(𝒪min​\(𝒪skeleton​\(𝒢\)\)\)\\mathcal\{O\}\_\{\\text\{atomic\}\}\(\\mathcal\{O\}\_\{\\text\{min\}\}\(\\mathcal\{O\}\_\{\\text\{skeleton\}\}\(\\mathcal\{G\}\)\)\)\.

Figure 3:Example to illustrate graph operators𝒪atomic\\mathcal\{O\}\_\{\\text\{atomic\}\},𝒪min\\mathcal\{O\}\_\{\\text\{min\}\}, and𝒪skeleton\\mathcal\{O\}\_\{\\text\{skeleton\}\}\.
### B\.8Definition of Graph Operators and Proof of[Theorem3](https://arxiv.org/html/2605.20396#Thmtheorem3)

We provide the definition of structure operations𝒪atomic\\mathcal\{O\}\_\{\\text\{atomic\}\},𝒪min\\mathcal\{O\}\_\{\\text\{min\}\}, and𝒪skeleton\\mathcal\{O\}\_\{\\text\{skeleton\}\}below, with an example in[Figure3](https://arxiv.org/html/2605.20396#A2.F3)\.

###### Definition 4\(Minimal\-graph operator\(Huang et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib38); Dong et al\.,[2023](https://arxiv.org/html/2605.20396#bib.bib28)\)\)\.

For every two atomic covers𝐋\\mathbf\{L\}and𝐏\\mathbf\{P\}in structure𝒢\\mathcal\{G\}, we merge𝐋\\mathbf\{L\}to𝐏\\mathbf\{P\}if the following conditions hold: \(i\)𝐋\\mathbf\{L\}is the pure children of𝐏\\mathbf\{P\}, \(ii\) all elements of𝐋\\mathbf\{L\}and𝐏\\mathbf\{P\}are latent and\|𝐋\|=\|𝐏\|\|\\mathbf\{L\}\|=\|\\mathbf\{P\}\|, and \(iii\) the pure children of𝐋\\mathbf\{L\}form a single atomic cover, or the siblings of𝐋\\mathbf\{L\}form a single atomic cover\. We denote such an operator as minimal\-graph operator𝒪min​\(𝒢\)\\mathcal\{O\}\_\{\\text\{min\}\}\(\\mathcal\{G\}\)\.

###### Definition 5\(Skeleton operator\(Huang et al\.,[2022](https://arxiv.org/html/2605.20396#bib.bib38); Dong et al\.,[2023](https://arxiv.org/html/2605.20396#bib.bib28)\)\)\.

Given an atomic cover𝐋\\mathbf\{L\}in structure𝒢\\mathcal\{G\}\. Consider𝒮\\mathcal\{S\}as the set of atomic covers such that for all𝐒∈𝒮\\mathbf\{S\}\\in\\mathcal\{S\}, we have𝐒⊆𝐋\\mathbf\{S\}\\subseteq\\mathbf\{L\}\. Let𝐂=PCh𝒢\(𝐋\)\\∪𝐒∈𝒮PCh𝒢\(𝐒\)\\mathbf\{C\}=\\text\{PCh\}\_\{\\mathcal\{G\}\}\(\\mathbf\{L\}\)\\backslash\\cup\_\{\\mathbf\{S\}\\in\\mathcal\{S\}\}\\text\{PCh\}\_\{\\mathcal\{G\}\}\(\\mathbf\{S\}\)\. We add edges from elements in𝐋\\mathbf\{L\}to elements in𝐂\\mathbf\{C\}, and we denote such an operator as skeleton operator𝒪skeleton​\(𝒢\)\\mathcal\{O\}\_\{\\text\{skeleton\}\}\(\\mathcal\{G\}\)\.

###### Definition 6\(Intra atomic operator\)\.

For every atomic cover𝐋\\mathbf\{L\}in structure𝒢\\mathcal\{G\}, if\|𝐋\|≥2\|\\mathbf\{L\}\|\\geq 2, then we add edges between elements in𝐋\\mathbf\{L\}such that𝐋\\mathbf\{L\}form a fully connected DAG\. We denote such an operator as intra atomic operator𝒪atomic​\(𝒢\)\\mathcal\{O\}\_\{\\text\{atomic\}\}\(\\mathcal\{G\}\)\.

###### Example 1\(Example for graph operations\)\.

Let the graph in[Figure3](https://arxiv.org/html/2605.20396#A2.F3)\(a\) be𝒢\\mathcal\{G\}\. By the skeleton operator, we add edges fromL2L\_\{2\}andL3L\_\{3\}toX7X\_\{7\}, and we arrive at𝒪skeleton​\(𝒢\)\\mathcal\{O\}\_\{\\text\{skeleton\}\}\(\\mathcal\{G\}\), which is shown in[Figure3](https://arxiv.org/html/2605.20396#A2.F3)\(b\)\. By the minimal graph operator, we deleteL5L\_\{5\}and directly linkL1L\_\{1\}toX8,X9,X10X\_\{8\},X\_\{9\},X\_\{10\}, and arrive at𝒪min​\(𝒪skeleton​\(𝒢\)\)\\mathcal\{O\}\_\{\\text\{min\}\}\(\\mathcal\{O\}\_\{\\text\{skeleton\}\}\(\\mathcal\{G\}\)\), which is shown in[Figure3](https://arxiv.org/html/2605.20396#A2.F3)\(c\)\. Finally, by the intra atomic operator, we add edges amongL2,L3,L4L\_\{2\},L\_\{3\},L\_\{4\}such that they are fully connected, and arrive at𝒪atomic​\(𝒪min​\(𝒪skeleton​\(𝒢\)\)\)\\mathcal\{O\}\_\{\\text\{atomic\}\}\(\\mathcal\{O\}\_\{\\text\{min\}\}\(\\mathcal\{O\}\_\{\\text\{skeleton\}\}\(\\mathcal\{G\}\)\)\), which is shown in[Figure3](https://arxiv.org/html/2605.20396#A2.F3)\(d\)\.

We now provide the proof of the following result\. See[3](https://arxiv.org/html/2605.20396#Thmtheorem3)

###### Proof\.

Since the search space contains the DAG𝒪min​\(𝒪skeleton​\(𝒢∗\)\)\\mathcal\{O\}\_\{\\textrm\{min\}\}\(\\mathcal\{O\}\_\{\\textrm\{skeleton\}\}\(\\mathcal\{G\}^\{\*\}\)\)that can generateΣX\\Sigma\_\{X\}in the large sample limit, the estimated DAG𝒢^\\hat\{\\mathcal\{G\}\}can also generateΣX\\Sigma\_\{X\}, because otherwise its score will be infinity and will not be a solution of the optimization problem\. Because𝒢∗\\mathcal\{G\}^\{\*\}andΣX\\Sigma\_\{X\}satisfy the generalized faithfulness assumption, we haveH​\(𝒢^\)=H​\(𝒢∗\)H\(\\hat\{\\mathcal\{G\}\}\)=H\(\\mathcal\{G\}^\{\*\}\)in the large sample limit by restricting the set of structures to those satisfying[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)and𝒢=𝒪min​\(𝒪skeleton​\(𝒢\)\)\\mathcal\{G\}=\\mathcal\{O\}\_\{\\textrm\{min\}\}\(\\mathcal\{O\}\_\{\\textrm\{skeleton\}\}\(\\mathcal\{G\}\)\)in[Theorem1](https://arxiv.org/html/2605.20396#Thmtheorem1)\. This indicates that𝒢^\\hat\{\\mathcal\{G\}\}andΣX\\Sigma\_\{X\}also satisfy the generalized faithfulness assumption\.

Let𝒢′\\mathcal\{G\}^\{\\prime\}be the structure estimated by Algorithm 1 inHuang et al\. \([2022](https://arxiv.org/html/2605.20396#bib.bib38)\)based on the covariance matrixΣX\\Sigma\_\{X\}\. Since𝒢^\\hat\{\\mathcal\{G\}\}satisfies[Assumption3](https://arxiv.org/html/2605.20396#Thmassumption3)and can faithfully generateΣX\\Sigma\_\{X\},Huang et al\. \([2022](https://arxiv.org/html/2605.20396#bib.bib38), Theorem 10\)andDong et al\. \([2023](https://arxiv.org/html/2605.20396#bib.bib28), Theorem 13\)imply that𝒪atomic​\(𝒢′\)\\mathcal\{O\}\_\{\\textrm\{atomic\}\}\(\\mathcal\{G\}^\{\\prime\}\)and𝒪atomic​\(𝒪min​\(𝒪skeleton​\(𝒢^\)\)\)\\mathcal\{O\}\_\{\\textrm\{atomic\}\}\(\\mathcal\{O\}\_\{\\textrm\{min\}\}\(\\mathcal\{O\}\_\{\\textrm\{skeleton\}\}\(\\hat\{\\mathcal\{G\}\}\)\)\)are Markov equivalent\. With similar reasoning, we can show that𝒪atomic​\(𝒢′\)\\mathcal\{O\}\_\{\\textrm\{atomic\}\}\(\\mathcal\{G\}^\{\\prime\}\)and𝒪atomic​\(𝒪min​\(𝒪skeleton​\(𝒢∗\)\)\)\\mathcal\{O\}\_\{\\textrm\{atomic\}\}\(\\mathcal\{O\}\_\{\\textrm\{min\}\}\(\\mathcal\{O\}\_\{\\textrm\{skeleton\}\}\(\\mathcal\{G\}^\{\*\}\)\)\)are Markov equivalent\. Therefore, by the transitivity of Markov equivalence,𝒪atomic​\(𝒪min​\(𝒪skeleton​\(𝒢^\)\)\)\\mathcal\{O\}\_\{\\textrm\{atomic\}\}\(\\mathcal\{O\}\_\{\\textrm\{min\}\}\(\\mathcal\{O\}\_\{\\textrm\{skeleton\}\}\(\\hat\{\\mathcal\{G\}\}\)\)\)and𝒪atomic​\(𝒪min​\(𝒪skeleton​\(𝒢∗\)\)\)\\mathcal\{O\}\_\{\\textrm\{atomic\}\}\(\\mathcal\{O\}\_\{\\textrm\{min\}\}\(\\mathcal\{O\}\_\{\\textrm\{skeleton\}\}\(\\mathcal\{G\}^\{\*\}\)\)\)are Markov equivalent\.

Recall that𝒢^=𝒪min​\(𝒪skeleton​\(𝒢^\)\)\\hat\{\\mathcal\{G\}\}=\\mathcal\{O\}\_\{\\textrm\{min\}\}\(\\mathcal\{O\}\_\{\\textrm\{skeleton\}\}\(\\hat\{\\mathcal\{G\}\}\)\)\. This implies that𝒪atomic​\(𝒢^\)\\mathcal\{O\}\_\{\\textrm\{atomic\}\}\(\\hat\{\\mathcal\{G\}\}\)and𝒪atomic​\(𝒪min​\(𝒪skeleton​\(𝒢∗\)\)\)\\mathcal\{O\}\_\{\\textrm\{atomic\}\}\(\\mathcal\{O\}\_\{\\textrm\{min\}\}\(\\mathcal\{O\}\_\{\\textrm\{skeleton\}\}\(\\mathcal\{G\}^\{\*\}\)\)\)are Markov equivalent\. ∎

## Appendix CSupplementary Experiment Details

Implementation details\.To improve the efficiency of[Algorithm2](https://arxiv.org/html/2605.20396#alg2), we iterate over the latent atomic covers to identify latent variables with the same set of parents and children\. For our exact search method, we use L\-BFGS\(Byrd et al\.,[1995](https://arxiv.org/html/2605.20396#bib.bib15)\)implemented throughSciPy\(Virtanen et al\.,[2020](https://arxiv.org/html/2605.20396#bib.bib80)\)andPyTorch\(Paszke et al\.,[2019](https://arxiv.org/html/2605.20396#bib.bib53)\)packages \(with the default hyperparameters\) to solve the optimization problem in[Equation3](https://arxiv.org/html/2605.20396#S3.E3)when computing BIC\. The experiments for the exact search method are conducted on1616CPUs in parallel\. For the continuous search method, i\.e\., SALAD\-CS, we use augmented Lagrangian method\(Bertsekas,[1982](https://arxiv.org/html/2605.20396#bib.bib9),[1999](https://arxiv.org/html/2605.20396#bib.bib10); Nocedal & Wright,[2006](https://arxiv.org/html/2605.20396#bib.bib51)\)to solve the continuous constrained optimization problem, in which each subproblem is solved using the Adam optimizer\(Kingma & Ba,[2014](https://arxiv.org/html/2605.20396#bib.bib41)\)with30003000iterations\. Furthermore,[Equation5](https://arxiv.org/html/2605.20396#S4.E5)involves a nonconvex optimization problem; similar to continuous optimization methods for causal discovery\(Ng et al\.,[2024](https://arxiv.org/html/2605.20396#bib.bib50)\), this procedure may yield suboptimal local solutions\. Thus, we run the SALAD\-CS method from1010random initializations, and select the final solution with the best score\.

For HUANG, GIN, and RCD, we use the publicly available implementations with default hyperparameters\. For FOFC, we use the implementation through thepy\-causalpackage\(Scheines et al\.,[1998](https://arxiv.org/html/2605.20396#bib.bib60)\)with Wishart test and significance level of0\.0010\.001\. Note that we also experimented with significance level of0\.010\.01,0\.050\.05and0\.10\.1, for which many of the runs are invalid \(because an error occurred\)\.

Metrics\.Since the goal is to recover the structure up to Markov equivalence, we compute the SHDs over the MECs\. Note that the labeling of the latent variables is not important; therefore, we calculate the SHDs of the estimated MECs over all possible permutations of the latent variables, and select the smallest SHD\. Similarly, we also compute the F1 scores of the estimated skeletons over all permutations of latent variables, and select the highest F1 score\.

For FOFC, an error occurred in some of the experimental runs\. Therefore, we additionally report the number of valid runs \(for which an error did not occur\)\.

## Appendix DRuntime and Computational Efficiency

In this section, we report the runtime for different methods considered\. For the 1\-factor models, our SALAD method has a runtime of8\.77±0\.738\.77\\pm 0\.73and44\.88±8\.0244\.88\\pm 8\.02minutes for1010and1111measured variables, respectively, while for hierarchical structures, it takes16\.11±2\.0116\.11\\pm 2\.01minutes\. For the SALAD\-CS method, each optimization run takes14\.17±0\.6914\.17\\pm 0\.69and14\.73±1\.9114\.73\\pm 1\.91minutes for1010and1111measured variables, respectively\. For the baselines, GIN, HUANG, and FOFC generally finish within one minute\. For the 1\-factor models, RCD requires17\.23±34\.8417\.23\\pm 34\.84and13\.06±24\.3313\.06\\pm 24\.33minutes for1010and1111measured variables, respectively, while for hierarchical structures, it has a runtime of7\.76±15\.337\.76\\pm 15\.33minutes\.

Our methods have a comparable runtime as RCD, but achieve better performance\. Although the runtime of our methods exceeds that of GIN, HUANG, and FOFC, the improvement in the causal discovery performance is significant\. As discussed in[AppendixC](https://arxiv.org/html/2605.20396#A3), our experiments are conducted on CPUs\. It is worth noting that the runtime may be further decreased by \(i\) conducting experiments with GPU acceleration \(specifically when using gradient\-based optimization to solve[Equations3](https://arxiv.org/html/2605.20396#S3.E3)and[5](https://arxiv.org/html/2605.20396#S4.E5)\), or \(ii\) performing more score computations of different structures \(specifically for exact search\) concurrently on different CPUs\.

Indeed, the relatively long runtime of our methods may be unsurprising because, even without latent variables, exact score\-based methods\(Singh & Moore,[2005](https://arxiv.org/html/2605.20396#bib.bib69); Yuan & Malone,[2013](https://arxiv.org/html/2605.20396#bib.bib84)\)are known to require a long runtime\. The search procedure developed in our work serves as a proof of concept, tailored for scenarios involving a relatively small number of variables\. Nonetheless, the empirical performance validates the effectiveness of score\-based methods for estimating latent variable causal models\. Future works include developing greedy approaches similar to GES to make the search procedure more efficient and scalable\.

Table 2:SHDs of MECs across various structural assumptions and sample sizes\. For each setting, the top two methods are in bold\. For FOFC, the number within the brackets indicates the number of valid runs \(for which an error did not occur\)\.![Refer to caption](https://arxiv.org/html/2605.20396v1/x3.png)\(a\)Example 1\.
![Refer to caption](https://arxiv.org/html/2605.20396v1/x4.png)\(b\)Example 2\.
![Refer to caption](https://arxiv.org/html/2605.20396v1/x5.png)\(c\)Example 3\.
![Refer to caption](https://arxiv.org/html/2605.20396v1/x6.png)\(d\)Example 4\.

Figure 4:Ground truths for 1\-factor latent variable models\.![Refer to caption](https://arxiv.org/html/2605.20396v1/x7.png)\(a\)Example 1\.
![Refer to caption](https://arxiv.org/html/2605.20396v1/x8.png)\(b\)Example 2\.
![Refer to caption](https://arxiv.org/html/2605.20396v1/x9.png)\(c\)Example 3\.
![Refer to caption](https://arxiv.org/html/2605.20396v1/x10.png)\(d\)Example 4\.

Figure 5:Ground truths for latent hierarchical structures\.