Score-Based Causal Discovery of Latent Variable Causal Models

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

• The paper introduces SALAD, a score-based framework that minimizes model degrees of freedom to identify linear Gaussian causal structures containing causally related latent variables.
• Under a one-factor model, SALAD achieves Markov-equivalence guarantees and near-perfect synthetic performance, reaching skeleton F1 scores of 0.99–1.00 at larger sample sizes.
• The method outperforms constraint-based baselines in finite-sample tests but faces scalability limits, an unproven hierarchical degrees-of-freedom conjecture, and unresolved theoretical justification for BIC.

Overview and problem setting

The paper develops SALAD (Score-bAsed Latent cAusal Discovery), a score-based framework for identifying linear latent variable causal models in which latent variables may themselves be causally related. The data-generating process is a linear Gaussian SEM over measured variables $X$ and latent variables $L$:

$L = CL + E_L, \qquad X = BL + E_X,$

with jointly independent Gaussian noise, and the goal is to recover the DAG $\mathcal{G}$ (defined by the supports of $B$ and $C$) up to a suitable notion of model equivalence. The authors position this as the first score-based method that identifies causal structures containing causally-related latent variables with identifiability guarantees.

The motivation is empirical as well as theoretical: existing methods for causally-related latents—FOFC [kummerfeld2016causal], HUANG [huang2022latent], GIN [xie2020generalized]—are constraint-based, relying on tetrad or rank-deficiency tests applied iteratively. Such procedures suffer from testing-order dependency, error propagation, and sensitivity to significance levels. Score-based methods such as GES avoid these issues in the causally sufficient case; extending them to latent variable models requires solving three problems: characterizing degrees of freedom of the marginal distribution over observed variables, formulating a scoring function, and structuring the search.

Likelihood formulation and parameter indeterminacy

A key technical observation is an indeterminacy lemma: for any parameters generating $\Sigma_X$, there exist alternative parameters with the same support but with $\Omega_L = I$ that generate the same marginal covariance. Consequently, $\Omega_L$ cannot be estimated from $\Sigma_X$ alone, and it is without loss of generality to fix $L$0 during estimation. This yields a closed-form Gaussian negative log-likelihood over $L$1, which underpins both the scoring function and the BIC criterion used in practice.

Scoring function, score equivalence, and algebraic equivalence

The core scoring function assigns structure $L$2 its degrees of freedom $L$3 if $L$4 can generate the empirical covariance matrix, and infinity otherwise. The intuition is parsimony: among structures consistent with the data, prefer the one with fewest free parameters. Two properties are established:

• Score equivalence: Markov equivalent DAGs (up to relabeling of latents) receive identical scores, so search can operate over Markov equivalence classes rather than individual DAGs.
• Consistency up to algebraic equivalence: under a generalized faithfulness assumption (an equality constraint holds only if entailed by the true structure), the global minimizer of the score satisfies $L$5, i.e., the estimated and true structures impose exactly the same equality constraints (conditional independences, rank deficiencies, possibly Verma constraints) on the observed marginal. The proof is short: any feasible estimator must generate $L$6; faithfulness gives one inclusion of constraint sets, and minimality of dimension rules out strict inclusion.

Algebraic equivalence is weaker than Markov equivalence in general—the paper gives an explicit example of two non-Markov-equivalent structures sharing the same vanishing ideal—and the authors are candid that translating algebraic equivalence into finer-grained guarantees requires additional structural assumptions. This motivates the two model classes studied next.

Linear 1-factor models

Under the structural assumption of Silva et al.—each measured variable has a single latent parent, and each latent has at least three measured children—the paper proves that $L$7, i.e., degrees of freedom equal edges plus measured variables, via generic finite-to-oneness of the polynomial parameter map and full column rank of its Jacobian. Combined with the general consistency result, this yields Markov equivalence (up to latent relabeling) between the score minimizer and the truth in the large-sample limit—a strictly stronger guarantee than algebraic equivalence.

An exact search procedure enumerates all candidate structures: latent MECs of size $L$8, ordered partitions of measured variables into child sets of size at least three, deduplicated by Markov equivalence. For computational efficiency, the authors also derive a continuous formulation (SALAD-CS) using binary mask matrices $L$9 constrained to encode the structural assumption, relaxed via Gumbel-Softmax/Gumbel-Sigmoid and solved with an augmented Lagrangian method. They acknowledge that this nonconvex program admits local optima, mitigated by running from ten random initializations.

Latent hierarchical structures

For the more general identifiable linear latent hierarchical graphs of Huang et al., where measured variables may have multiple latent parents and some latents have no measured children, the degrees-of-freedom story changes qualitatively. A proposition shows that when $L = CL + E_L, \qquad X = BL + E_X,$0 latent variables share the same parent and children sets (with enough parents or children), the degrees of freedom drop by $L = CL + E_L, \qquad X = BL + E_X,$1 relative to $L = CL + E_L, \qquad X = BL + E_X,$2: certain edges become redundant because orthogonal transformations of $L = CL + E_L, \qquad X = BL + E_X,$3 and $L = CL + E_L, \qquad X = BL + E_X,$4 (via QR-type decompositions) preserve the generated covariance while zeroing out those edges. An algorithm computes an upper bound on $L = CL + E_L, \qquad X = BL + E_X,$5 by detecting maximal groups of latents with identical adjacency to the rest of the graph.

Notably, the tightness of this upper bound is stated as a conjecture, supported by simulations over 10,000 examples (Jacobian rank computations) rather than proof; the authors note that a formal argument would require tools from algebraic statistics, citing the technical difficulty even for sparse factor analysis with independent latents. Consistency is then established only up to composition of graph operators ($L = CL + E_L, \qquad X = BL + E_X,$6, $L = CL + E_L, \qquad X = BL + E_X,$7, $L = CL + E_L, \qquad X = BL + E_X,$8) that preserve rank constraints—matching the best achievable resolution of the constraint-based method it unifies, since these operators provably do not alter the observable rank constraints. No continuous search variant is provided here, as the structural assumption does not reduce to inequality constraints.

On the BIC score

Because the pure dimension-based score is justified only asymptotically, experiments use BIC. The authors flag a genuine theoretical gap: latent-variable models form stratified exponential families, so standard BIC justifications (Schwarz; Haughton's curved-exponential-family result) do not directly apply. Empirically, however, BIC performs well—F1 scores approach one at large sample sizes—which they interpret as evidence that BIC is valid under these structural assumptions, while explicitly leaving its theoretical justification as open work.

Experiments

On synthetic 1-factor models and hierarchical structures (10–11 measured variables, three random trials per ground truth, sample sizes 100 to 10,000), SALAD dominates constraint-based baselines. Representative skeleton F1 scores:

Model	$L = CL + E_L, \qquad X = BL + E_X,$9	SALAD	SALAD-CS	HUANG	FOFC	GIN
1-factor	100	0.99	0.97	0.50	0.90 (8/12 valid runs)	0.35
1-factor	3000	1.00	0.99	0.93	0.99	0.35
Hierarchical	100	0.92	N/A	0.57	N/A (0 valid runs)	0.48
Hierarchical	10000	0.98	N/A	0.88	N/A	0.47

SHD results show the same pattern: SALAD attains near-zero SHD on 1-factor models at $\mathcal{G}$0, versus 1.67 for HUANG and roughly 15 for GIN throughout. The advantage is largest at small sample sizes, consistent with the hypothesis that constraint-based baselines degrade through error propagation. FOFC failed entirely on hierarchical structures and produced errors in several 1-factor runs regardless of significance level. The cost is runtime: SALAD takes roughly 9–45 minutes per instance versus under a minute for GIN, HUANG, and FOFC, though comparable to RCD. The authors present the exact search as a proof of concept suited to small numbers of variables.

Limitations and open questions

Several limitations are conceded explicitly. First, exact search scales poorly, and no greedy analogue of GES exists yet for latent variable models. Second, the degrees-of-freedom characterization for hierarchical structures is an unproven conjecture, albeit empirically well-supported. Third, the use of BIC lacks theoretical justification in stratified exponential families. Fourth, identifiability rests on generalized faithfulness, whose violation set has measure zero but which is still an assumption on the true parameters; the authors note that replacing it with a sparsest-Markov-representation-style condition is possible but less informative in this setting. Finally, the continuous method is nonconvex and heuristic, and the hierarchical setting admits identification only up to rank-preserving graph operators, not full Markov equivalence.

Conclusion

This paper establishes that properly formulated scores—degrees-of-freedom minimization or BIC—are score equivalent and consistent for linear latent variable causal models, delivering Markov-equivalence-level guarantees under the 1-factor assumption and operator-level guarantees under the hierarchical assumption. It thereby provides a unified score-based perspective on tetrad- and rank-based constraint methods, with substantially better finite-sample performance in experiments. The main outstanding questions are scalable greedy search, a proof of the degrees-of-freedom conjecture for hierarchical structures, and theoretical grounding for BIC in stratified exponential families.