Paper deep dive
Differentiable Structure Learning for Cyclic Linear Gaussian Models with Latent Confounders
Sadegh Khorasani, Ali Najar, Saber Salehkaleybar, Negar Kiyavash
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 95%
Last extracted: 10/1/2026, 8:43:38 AM
Summary
This paper presents a differentiable structure learning method for cyclic linear Gaussian structural causal models (SCMs) with exogenous latent confounders. The authors derive the observed covariance family, introduce marginal quasi-equivalence to handle non-identifiable latent structures, and formulate structure learning as minimizing a Gaussian negative log-likelihood with a complexity penalty. They use Bernoulli gates to parameterize edge and latent variable inclusion, enabling joint optimization of discrete structure and continuous parameters. Theoretical consistency is established under algebraic faithfulness and structural minimality, and experiments show lower recovery error compared to baselines.
Entities (7)
Relation Signals (7)
Sadegh Khorasani → affiliatedwith → EPFL
confidence 99% · Sadegh Khorasani1 ... 1School of Computer and Communication Sciences, EPFL, Lausanne, Switzerland
Bernoulli Gates → parameterizes → Directed Edges
confidence 96% · We parameterize the inclusion of directed edges and candidate latent variables using Bernoulli gates
Bernoulli Gates → parameterizes → Candidate Latent Variables
confidence 96% · We parameterize the inclusion of directed edges and candidate latent variables using Bernoulli gates
Differentiable Structure Learning Method → handles → Cyclic Linear Gaussian Structural Causal Models
confidence 95% · We study causal structure learning from observational data in linear Gaussian structural causal models in the presence of directed cycles
Differentiable Structure Learning Method → handles → Exogenous Latent Confounders
confidence 95% · We study causal structure learning from observational data in linear Gaussian structural causal models in the presence of directed cycles and an unknown number of exogenous latent confounders
Structure Learning → formulatedas → Minimization of Gaussian Negative Log-Likelihood
confidence 94% · We formulate structure learning as minimization of the Gaussian negative log-likelihood with a logarithmically scaled complexity penalty
Marginal Quasi-Equivalence → characterizes → Observational Distribution Overlap
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We study causal structure learning from observational data in linear Gaussian structural causal models in the presence of directed cycles and an unknown number of exogenous latent confounders, bounded by a given maximum. We derive the covariance of the observed variables and introduce marginal quasi-equivalence, which characterizes when different causal models share a full-dimensional subset of the observational distributions they can generate. We formulate structure learning as minimization of the Gaussian negative log-likelihood with a logarithmically scaled complexity penalty that counts directed edges and latent variables. For a fixed number of observed variables and a fixed upper bound on latent variables, we establish consistency of global score minimizers up to marginal quasi-equivalence under algebraic faithfulness, structural minimality, and model-overlap assumptions. We parameterize the inclusion of directed edges and candidate latent variables using Bernoulli gates, whose continuous probabilities are optimized jointly with the structural coefficients. Averaging the penalized negative log-likelihood over these gates yields an objective with a closed-form differentiable complexity penalty. We prove that this expected objective has the same global infimum as the corresponding discrete structure-learning objective. Experimental results show that our approach achieves lower recovery error than previous methods in several experimental settings.
Tags
Links
- Source: https://arxiv.org/abs/2609.38618v1
- Canonical: https://arxiv.org/abs/2609.38618v1
Trouble viewing inline? Open PDF directly →
Full Text
173,866 characters extracted from source content.
Expand or collapse full text
Differentiable Structure Learning for Cyclic Linear Gaussian Models with Latent Confounders Sadegh Khorasani1, Ali Najar1, Saber Salehkaleybar2, Negar Kiyavash3 1School of Computer and Communication Sciences, EPFL, Lausanne, Switzerland 2Leiden Institute of Advanced Computer Science (LIACS), Leiden University, Leiden, The Netherlands 3College of Management of Technology, EPFL, Lausanne, Switzerland sadegh.khorasani@epfl.ch anajar13750@gmail.com s.salehkaleybar@liacs.leidenuniv.nl negar.kiyavash@epfl.ch Abstract We study causal structure learning from observational data in linear Gaussian structural causal models in the presence of directed cycles and an unknown number of exogenous latent confounders, bounded by a given maximum. We derive the covariance of the observed variables and introduce marginal quasi-equivalence, which characterizes when different causal models share a full-dimensional subset of the observational distributions they can generate. We formulate structure learning as minimization of the Gaussian negative log-likelihood with a logarithmically scaled complexity penalty that counts directed edges and latent variables. For a fixed number of observed variables and a fixed upper bound on latent variables, we establish consistency of global score minimizers up to marginal quasi-equivalence under algebraic faithfulness, structural minimality, and model-overlap assumptions. We parameterize the inclusion of directed edges and candidate latent variables using Bernoulli gates, whose continuous probabilities are optimized jointly with the structural coefficients. Averaging the penalized negative log-likelihood over these gates yields an objective with a closed-form differentiable complexity penalty. We prove that this expected objective has the same global infimum as the corresponding discrete structure-learning objective. Experimental results show that our approach achieves lower recovery error than previous methods in several experimental settings. 1 Introduction Causal graphical models describe causal relationships among variables. Learning their structure from observational data is a central problem in causal inference. In many systems, these relationships involve feedback loops and unobserved common causes, known as latent confounders. This motivates structure-learning methods that allow both directed cycles and latent variables. Several lines of research address these two challenges separately. Work on cyclic linear Gaussian models has developed distributional characterizations and score-based learning methods under the assumption of no latent confounders (Ghassami et al., 2020; Yi and Lee, 2024). In contrast, observational methods for latent confounding commonly retain acyclicity (Kaltenpoth and Vreeken, 2023; Ma et al., 2024). Approaches that allow both cycles and confounding use different model or data assumptions. In particular, Amendola et al. (2020) consider mixed graphs, whose bidirected edges represent possible correlations between structural noise variables. Their simple-graph restriction permits at most one edge of any type between each pair of observed variables. It therefore excludes reciprocal directed edges and pairs joined by both a directed and a bidirected edge. Other approaches use interventional data (Hyttinen et al., 2012; Sethuraman and Fekri, 2025) or conditional-independence assumptions for nonlinear cyclic models (Forré and Mooij, 2018). The distributional-equivalence and consistency results for cyclic models without latent variables do not directly extend after latent variables are marginalized, since candidate graphs may contain different numbers of latent variables and have different structural parameterizations. Likewise, results established under the simple-graph restriction do not directly apply to our setting, where these more general structures are allowed. We do not impose this restriction and allow arbitrary directed cycles among observed variables together with explicit exogenous latent confounders. We study observational linear Gaussian structural causal models with arbitrary directed cycles among observed variables and an unknown, but bounded, number of independent exogenous latent confounders. Rather than recovering a unique latent graph, which is generally not identifiable from the observed distribution, we formulate the problem directly in terms of the families of observed covariance matrices induced by candidate graphs. This perspective leads to a marginal notion of quasi-equivalence that captures when two graphs generate a full-dimensional common subset of observational distributions, while allowing their explicit latent structures to differ. Our contributions are as follows. • We derive the observed covariance family induced by a cyclic linear Gaussian SCM with exogenous latent variables and establish its basic geometric properties. We define marginal distribution equivalence and marginal quasi-equivalence in observed covariance space, extending the positive-overlap perspective of Ghassami et al. (2020, Definition 9) to candidate graphs with different latent-variable structures (Section 3.1). • We formulate structure learning as joint optimization over the discrete graph structure and its continuous parameters, using the Gaussian negative log-likelihood and a logarithmically scaled penalty that counts directed edges and latent variables. For a fixed number of observed variables and a fixed upper bound on latent variables, we show that global score minimizers asymptotically select covariance families whose closures contain the true covariance. Under additional conditions (in particular, under algebraic faithfulness, structural minimality, and full-dimensional model overlap), we obtain consistency up to marginal quasi-equivalence (Section 4.1). • We use Bernoulli gates for observed edges, latent-variable activation, and latent-to-observed edges, optimizing their probabilities jointly with the continuous model parameters. Averaging the penalized objective over these gates yields a differentiable expected objective with a closed-form complexity penalty. We prove that this expected objective has the same global infimum as the corresponding discrete structure-learning problem, so the continuous probability parameterization preserves the globally optimal objective value (Section 4.2). • We evaluate the approach on synthetic cyclic SCMs with varying graph density and latent confounding, as well as on simulated GeneNetWeaver networks. Across several settings with latent variables, the method achieves lower bidirectional covariance-family compatibility error than the considered baselines, while scaling to substantially larger observed graphs. (Section 5). Appendix E provides further comparisons with related work. 2 Problem Setting We study causal structure learning from observational data in linear Gaussian structural causal models (SCMs) with feedback loops among observed variables and exogenous latent confounders, meaning latent variables with no incoming edges. The unknown structure consists of edges among observed variables, edges from latent to observed variables, and the number of latent variables. 2.1 Notation and conventions For a positive integer d, let IdI_d denote the d×d× d identity matrix, dS^d the space of real symmetric d×d× d matrices, and ++d⊂dS_++^d ^d the set of symmetric positive-definite matrices. For a symmetric matrix A, the notation A≻0A 0 and A⪰0A 0 means positive definiteness and positive semidefiniteness, respectively. We write A⊤A for the transpose of A and, when A is square, det(A) (A), tr(A)tr(A), and ρ(A)ρ(A) for its determinant, trace, and spectral radius (the largest absolute value of its eigenvalues). For invertible A, we use A−⊤:=(A−1)⊤A^- :=(A^-1) . The operator diag(a1,…,ad)diag(a_1,…,a_d) forms a diagonal matrix, and ⊙ denotes the elementwise product. For a matrix A, |A||A| denotes entrywise absolute value, ‖A‖1=∑i,j|Aij|\|A\|_1= _i,j|A_ij| the entrywise ℓ1 _1 norm, ‖A‖2\|A\|_2 the spectral norm, and ‖A‖F\|A\|_F the Frobenius norm. For a finite set, vertical bars instead denote cardinality. Random vectors are treated as column vectors, and Cov(Y)Cov(Y) denotes the covariance matrix of a random vector Y. For subsets T⊆ST S of a Euclidean space, clS(T)cl_S(T) denotes the closure of T relative to S: all limits in S of convergent sequences from T. An algebraic set is the set of common solutions to finitely many polynomial equalities with real coefficients. It is irreducible if it is nonempty and cannot be written as the union of two strictly smaller algebraic sets. A semialgebraic set is a finite union of sets, each defined by finitely many polynomial equalities and inequalities. To say that a polynomial vanishes on a set means that it equals zero at every point in that set. For a nonempty semialgebraic set T⊆ℝmT ^m, dimT T is the largest k∈0,…,mk∈\0,…,m\ for which projecting T onto some k coordinates gives a set containing a nonempty open ball in ℝkR^k. Here projection means keeping the selected coordinates and discarding the others. We set dim∅=−1 =-1. For nonempty real algebraic sets, this dimension agrees with algebraic dimension (Coste, 2000, Section 3.3). For sets of symmetric p×p× p matrices, the coordinates are the p(p+1)/2p(p+1)/2 upper-triangular entries. For semialgebraic sets T⊆ST S, T is called full-dimensional in S when dimT=dimS T= S. Appendix A.3 gives formal definitions (including open balls), references, and examples. A candidate graph is denoted by G+=(V,E)G^+=(V,E), where V=O∪˙LV=O ∪L is partitioned into the observed vertex set O and the latent vertex set L. We write p=|O|p=|O|, ℓ=|L| =|L|, and 0≤ℓ≤Lmax0≤ ≤ L_ , where Lmax<∞L_ <∞ is a prescribed upper bound on the number of latent variables. The superscript ++ indicates that the graph includes both types of vertices. Matrices are indexed with observed variables before latent variables. We use the row-tail, column-head convention: BijB_ij is the coefficient associated with the directed edge i→ji→ j. Graph compatibility requires Bij=0B_ij=0 whenever (i,j)∉E(i,j)∉ E. Conversely, an edge in E permits, but does not require, its coefficient to be nonzero. Thus a graph specifies which coefficients are allowed to be nonzero, not which must be nonzero at every parameter value. Covariance-family notation and the two notions of graph equivalence are defined in Section 3.1. 2.2 Cyclic structural equations with exogenous confounders Let X=(XO⊤,XL⊤)⊤∈ℝp+ℓX=(X_O ,X_L ) ^p+ collect the observed and latent variables, and partition the structural noise vector N in the same order. The full zero-mean Gaussian SCM is X=B⊤X+N,N∼(0,Ω),Ω≻0 diagonal,X=B X+N, N (0, ), 0 diagonal, (1) where B∈ℝ(p+ℓ)×(p+ℓ)B ^(p+ )×(p+ ) respects the graph-imposed zero restrictions and Bii=0B_i=0. The diagonal Gaussian noise covariance makes the components of N mutually independent. We allow arbitrary directed cycles among observed variables and require only that Ip+ℓ−BI_p+ -B be invertible. For each noise realization, the unique solution is X=(Ip+ℓ−B)−⊤NX=(I_p+ -B)^- N, which defines a nonsingular Gaussian distribution. Latent vertices are exogenous: they have no parents, and their only permitted children are observed vertices. We additionally require each latent vertex to have at least two permitted edges to observed vertices. The resulting block structure is B=[BOO0Λ0],Ω=[ΩO00ΩL],B= bmatrixB_O&0\\ &0 bmatrix, = bmatrix _O&0\\ 0& _L bmatrix, (2) where BOO∈ℝp×pB_O ^p× p, Λ∈ℝℓ×p × p, and ΩO _O and ΩL _L are positive diagonal matrices of sizes p×p× p and ℓ×ℓ × , respectively. The hhth row of Λ contains the coefficients from the hhth latent variable to its observed children; these coefficients are called latent loadings. Equation 2 gives XO=BOO⊤XO+Λ⊤XL+NO,XL=NL.X_O=B_O X_O+ X_L+N_O, X_L=N_L. Thus the latent variables are mutually independent and independent of the observed structural noises NON_O. The block-triangular structure also gives det(Ip+ℓ−B)=det(Ip−BOO) (I_p+ -B)= (I_p-B_O), so the structural equations have a unique solution for every noise realization exactly when det(Ip−BOO)≠0 (I_p-B_O)≠ 0. The two-child restriction removes zero- and one-child latent vertices without excluding any observed distribution representable with such vertices. A latent variable with no child has no observational effect; the variance contribution of a latent variable with only one child can be absorbed into that child’s diagonal noise variance. This restriction is imposed on the permitted graph edges, not on the number of nonzero loadings at each parameter value; other latent variables may still be removable without changing the observed distribution. When ℓ=0 =0, the latent blocks are empty and their covariance contribution is understood to be the zero p×p× p matrix. 2.3 Observed marginal distribution We assume that only XOX_O is observed. The statistical analysis uses independent and identically distributed (i.i.d.) observations from its zero-mean marginal distribution. Eliminating XLX_L introduces the effective noise εO:=NO+Λ⊤NL,XO=BOO⊤XO+εO. _O:=N_O+ N_L, X_O=B_O X_O+ _O. Although the components of N are independent, the components of εO _O need not be independent. Lemma 1 (Observed marginal covariance). Under the preceding assumptions, XO∼(0,ΣO),ΣO=(Ip−BOO)−⊤Ψ(Ip−BOO)−1,X_O (0, _O), _O=(I_p-B_O)^- (I_p-B_O)^-1, (3) where the effective noise covariance is Ψ=Cov(εO)=ΩO+Λ⊤ΩLΛ. =Cov( _O)= _O+ _L . (4) Both Ψ and ΣO _O are positive definite. The matrix ΩO _O describes the observed structural noise covariance, whereas Λ⊤ΩLΛ _L describes the covariance contribution of the latent variables. The proof is given in Appendix A.1. 3 Marginal Equivalence and Score We compare graphs through their observed covariance families and define the penalized likelihood used for structure learning. 3.1 Marginal covariance families and equivalence For a fixed graph G+G^+, let (G+)P(G^+) denote the set of parameters θ=(BOO,Λ,ΩO,ΩL)θ=(B_O, , _O, _L) satisfying its graph-imposed zero restrictions, positive diagonal noise variances, and det(Ip−BOO)≠0 (I_p-B_O)≠ 0. We call these parameter values admissible. The graph’s observed covariance family is ℳ(G+)=ΣO(θ):θ∈(G+)⊆++p.M(G^+)=\ _O(θ):θ (G^+)\ _++^p. (5) Because the observed distributions are Gaussian with known zero mean, ℳ(G+)M(G^+) determines the entire observed distribution family. We use two closures of this covariance family. Write ℳ¯(G+)=cl++pℳ(G+) M(G^+)=cl_S_++^pM(G^+) for its relative Euclidean closure. It contains all positive-definite limits of covariances generated by admissible parameter sequences, but does not include singular covariance matrices. In contrast, the real Zariski closure V(G+)=ℳ(G+)¯Zar⊆pV(G^+)= M(G^+)^Zar ^p consists of all symmetric matrices satisfying every polynomial equality that holds for every covariance matrix in ℳ(G+)M(G^+). This definition uses only polynomial equalities. Matrices in V(G+)V(G^+) need not be positive definite or satisfy the other inequalities required to belong to ℳ(G+)M(G^+). Appendix A.4 explains the containment relations between these sets and the semialgebraicity of the relative closure. Proposition 2 (Model geometry). For each fixed graph, ℳ(G+)M(G^+) is semialgebraic, V(G+)V(G^+) is irreducible, and dimℳ(G+)=dimℳ¯(G+)=dimV(G+). (G^+)= M(G^+)= V(G^+). (6) The proof is in Appendix A. Definition 3 (Marginal equivalence). Let G1+G_1^+ and G2+G_2^+ have the same labeled observed vertex set; their latent vertex sets may differ. They are marginally distribution equivalent, written G1+≡OG2+G_1^+ _OG_2^+, if ℳ(G1+)=ℳ(G2+)M(G_1^+)=M(G_2^+). They are marginally quasi-equivalent, written G1+≅OG2+G_1^+ _OG_2^+, if dim(ℳ(G1+)∩ℳ(G2+))=dimℳ(G1+)=dimℳ(G2+). split (M(G_1^+) (G_2^+) )&= (G_1^+)= (G_2^+). split (7) Distribution equivalence implies quasi-equivalence. For quasi-equivalence, the intersection must have the dimension of each family, not necessarily p(p+1)/2p(p+1)/2, the dimension of pS^p. A lower-dimensional intersection in either family is insufficient. This definition follows the positive-measure overlap idea of Ghassami et al. (2020, Definition 9), expressed here through the dimension of the intersection of observed covariance families. Appendix A.6 explains this formulation for models with latent variables and why our quasi-equivalence guarantees compare each estimated graph directly with the true graph. 3.2 Candidate graph class and structural complexity For fixed p and Lmax<∞L_ <∞, let p,LmaxG_p,L_ be the class of admissible causal graphs with the same p labeled observed vertices and at most LmaxL_ exogenous latent vertices. These graphs have no self-loops or edges into latent vertices, may contain directed cycles among observed variables, and have at least two edges from each latent vertex to observed vertices. With a fixed convention for labeling latent vertices, only finitely many such graph structures are possible11 1 Finiteness is used in the statistical analysis; it does not require enumerating these graphs during optimization.. Relabeling latent vertices changes neither the observed covariance family nor the structural complexity. We measure structural complexity by c(G+)=|E(G+)|+ℓ(G+),c(G^+)=|E(G^+)|+ (G^+), (8) where |E(G+)||E(G^+)| counts edges among observed variables and edges from latent to observed variables, and ℓ(G+):=|L(G+)| [rgb]0,0,0 (G^+):=|L(G^+)| counts latent variables; we write ℓ when the graph is fixed. Each edge and each latent variable contributes one unit to c(G+)c(G^+). An edge is counted whenever it belongs to the graph, even if its fitted coefficient is zero. Its additive form gives the closed-form expected penalty in Section 4.2; a dimension-based penalty would require separate covariance-family dimension calculations. Appendix B.3 explains this distinction with a two-variable example and discusses the dimension-based alternative. 3.3 Observed-data objective and graph score Let XO(1),…,XO(n)X_O^(1),…,X_O^(n) be i.i.d. observations from (0,Σ∗)N(0, ^*), with Σ∗≻0 ^* 0. The mean is known to be zero, so the sample covariance is Sn=1n∑k=1nXO(k)XO(k)⊤.S_n= 1n _k=1^nX_O^(k)X_O^(k) . We consider n≥pn≥ p, so that Sn≻0S_n 0 almost surely (with probability one). For a candidate covariance Σ≻0 0, the observed-data Gaussian negative log-likelihood is ℒnfull(Σ)=n2plog(2π)+logdetΣ+tr(Σ−1Sn).L_n^full( )= n2 \p (2π)+ +tr( ^-1S_n) \. (9) The term nplog(2π)/2np (2π)/2 comes from the Gaussian density normalization and is independent of the graph (p is fixed) and its parameters. The superscript “full” indicates that this constant is included. Appendix A.2 derives this likelihood and its expression in structural parameters. For G+∈p,LmaxG^+ _p,L_ and θ∈(G+)θ (G^+), write ℒnfull(θ):=ℒnfull(ΣO(θ))L_n^full(θ):=L_n^full( _O(θ)). Our joint objective is Jn(G+,θ)=ℒnfull(θ)+λnc(G+),λn=12logn.J_n(G^+,θ)=L_n^full(θ)+ _nc(G^+), _n= 12 n. (10) We optimize JnJ_n over both G+∈p,LmaxG^+ _p,L_ and θ∈(G+)θ (G^+). The penalty balances fit against the number of edges and latent variables. Its scaling follows the Bayesian information criterion (BIC), but c(G+)c(G^+) need not equal the covariance-family dimension; we therefore do not claim a standard BIC interpretation. In contrast, Amendola et al. (2020) use a dimension-based BIC term; Appendix E gives the comparison. For theoretical analysis, define the penalized graph score by taking the infimum of the joint objective over the parameters of a fixed graph: sn(G+):=infθ∈(G+)Jn(G+,θ).s_n(G^+):= _θ (G^+)J_n(G^+,θ). (11) Minimizing this score over graphs gives the infimum of the joint objective. This theoretical definition does not require enumerating graphs during computation; Section 4.2 introduces Bernoulli gates for joint structure and parameter optimization. To express this same score in covariance space, define the Gaussian Kullback–Leibler (KL) divergence for S,Σ∈++pS, _++^p by D(S∥Σ)=12tr(Σ−1S)−logdet(Σ−1S)−p.D(S\| )= 12\tr( ^-1S)- ( ^-1S)-p\. (12) This is the KL divergence from (0,S)N(0,S) to (0,Σ)N(0, ); it is nonnegative and equals zero exactly when S=ΣS= . Let δn(G+)=infΣ∈ℳ¯(G+)D(Sn∥Σ) _n(G^+)= _ ∈ M(G^+)D(S_n\| ). Thus δn(G+) _n(G^+) measures the best fit to SnS_n within the graph’s covariance-family closure. Taking the infimum over this closure gives the same value as taking it over admissible parameters, even when no parameter value attains the infimum; Appendix A.2 proves this equality. Let Cn:=ℒnfull(Sn)C_n:=L_n^full(S_n), the negative log-likelihood evaluated at the sample covariance. Subtracting CnC_n from Equation 9 gives nD(Sn∥Σ)nD(S_n\| ). Thus the same graph score can be written as sn(G+)=Cn+nδn(G+)+λnc(G+),s_n(G^+)=C_n+n _n(G^+)+ _nc(G^+), (13) The constant CnC_n depends only on the observed data and has no effect on graph selection. This equivalent expression is used in the consistency analysis below. 4 Method and Theory We first establish statistical guarantees for global minimizers of the graph score induced by Equation 10. We then formulate optimization of the same joint objective using independent Bernoulli gates. 4.1 Closure consistency and geometric identification Let G∗+∈p,LmaxG^*+ _p,L_ be the data-generating graph and choose a minimum-complexity distribution-equivalent representative G†∈argminG+∈p,Lmaxc(G+):G+≡OG∗+.G ∈ *arg\,min_G^+ _p,L_ \c(G^+):G^+ _OG^*+\. (14) The minimum exists because p,LmaxG_p,L_ is finite and contains G∗+G^*+. Write c†:=c(G†)c :=c(G ). The reference G†G is not an input to the estimator. All such representatives have the same family ℳ(G∗+)M(G^*+), closure V(G∗+)V(G^*+), and complexity c†c , so the conditions below do not depend on the choice. Define the set of candidate graphs whose covariance families contain or can approximate the true covariance arbitrarily closely and whose complexity is at most c†c : 0(Σ∗):=G+∈p,Lmax:Σ∗∈ℳ¯(G+),c(G+)≤c†.G_0( ^*):= \G^+ _p,L_ : ^*∈ M(G^+),\ c(G^+)≤ c \. (15) Theorem 4 (Closure-selection consistency). Suppose p and LmaxL_ are fixed and the data are i.i.d. from (0,Σ∗)N(0, ^*), where Σ∗∈ℳ(G∗+) ^* (G^*+). Then PrargminG+∈p,Lmaxsn(G+)⊆0(Σ∗)⟶1. \! \ *arg\,min_G^+ _p,L_ s_n(G^+) _0( ^*) \ 1. (16) Equivalently, with probability tending to one, every global minimizer G^n+ G_n^+ satisfies Σ∗∈ℳ¯(G^n+),c(G^n+)≤c†. ^*∈ M( G_n^+), c( G_n^+)≤ c . The same conclusion holds for any deterministic penalty sequence satisfying λn→∞ _n→∞ and λn/n→0 _n/n→ 0. The proof is in Appendix A. This result concerns global score minimizers, not arbitrary numerical outputs, and requires none of the identification assumptions below. Recovery up to marginal quasi-equivalence needs an additional condition. Assumption 5 (Identification at the true covariance). Assume that every G+∈0(Σ∗)G^+ _0( ^*) satisfies G+≅OG†G^+ _OG . Corollary 6 (Consistency up to marginal quasi-equivalence). Under the conditions of Theorem 4 and Assumption 5, PrargminG+∈p,Lmaxsn(G+)⊆G+∈p,Lmax:G+≅OG∗+⟶1. \! \ *arg\,min_G^+ _p,L_ s_n(G^+) \G^+ _p,L_ :G^+ _OG^*+\ \ 1. (17) The proof is given in Appendix A.10. Assumption 5 is sufficient for this conclusion; necessity is not claimed. The following three assumptions are sufficient to establish it. They concern inclusions between Zariski closures, structural complexity, and intersections of covariance families. Assumption 7 (Relevant-candidate algebraic faithfulness). For every G+∈p,LmaxG^+ _p,L_ satisfying c(G+)≤c†c(G^+)≤ c , Σ∗∈ℳ¯(G+)⟹V(G†)⊆V(G+). ^*∈ M(G^+) V(G ) V(G^+). Assumption 8 (Structural minimality). For every G+∈p,LmaxG^+ _p,L_ , if V(G†)⊆V(G+)anddimℳ(G†)<dimℳ(G+),V(G ) V(G^+) (G )< (G^+), then c(G+)>c†c(G^+)>c . Assumption 9 (Full-dimensional overlap). If Σ∗∈ℳ¯(G+) ^*∈ M(G^+), c(G+)≤c†c(G^+)≤ c , and V(G+)=V(G†)=V(G^+)=V(G )=V, then dim(ℳ(G+)∩ℳ(G†))=dimV. (M(G^+) (G ) )= V. (18) These assumptions connect closure-selection consistency to recovery up to marginal quasi-equivalence. Lemma 17 gives the proof. Appendix B examines the scope of these assumptions. It shows that algebraic faithfulness and full-dimensional overlap hold for almost every admissible choice of the true model parameters (genericity). Structural minimality requires separate justification; the appendix gives sufficient conditions and a counterexample showing that it can fail. 4.2 Bernoulli formulation for joint structure optimization To optimize the joint objective in Equation 10 over edges among observed variables, edges from latent to observed variables, and the number of latent variables, we reserve LmaxL_ slots for potential latent variables and use mutually independent binary gates: Mij∼Bernoulli(qij)(i≠j),Mii=qii=0,Hhj∼Bernoulli(rhj),Zh∼Bernoulli(πh). gatheredM_ij (q_ij) (i≠ j), M_i=q_i=0,\\ H_hj (r_hj), Z_h ( _h). gathered (19) Here MijM_ij selects the edge i→ji→ j among observed variables, ZhZ_h activates latent slot h, and HhjH_hj selects its edge to observed variable j when that slot is active. The entries of Q=(qij)Q=(q_ij), R=(rhj)R=(r_hj), and π=(πh)π=( _h) lie in [0,1][0,1] and are the probabilities that the corresponding gates equal one. We optimize these probabilities; the sampled gates remain zero or one. Let ϑ=(W,Γ,ΩO) =(W, , _O), where W∈ℝp×pW ^p× p, Γ∈ℝLmax×p ^L_ × p, and ΩO≻0 _O 0 is diagonal. The matrix W contains coefficients for edges among observed variables, and Γ contains latent loadings with latent noise standard deviations absorbed into them. The same continuous parameters are used for all sampled configurations. Define BOO(M)=M⊙W,A(M)=Ip−M⊙W,Ψ(H,Z)=ΩO+(H⊙Γ)⊤diag(Z)(H⊙Γ). splitB_O(M)&=M W, A(M)=I_p-M W,\\ (H,Z)&= _O+(H ) diag(Z)(H ). split (20) For each configuration, let Ln(M,H,Z,ϑ)L_n(M,H,Z; ) be the Gaussian negative log-likelihood in Equation 9, evaluated at A(M)−⊤Ψ(H,Z)A(M)−1A(M)^- (H,Z)A(M)^-1 and with nplog(2π)/2np (2π)/2 removed. Set Ln=+∞L_n=+∞ if A(M)A(M) is singular. Let c(M,H,Z)c(M,H,Z) count selected edges among observed variables, active latent variables, and selected edges from active latents to observed variables. Appendix A.11 gives their explicit formulas. The exact Bernoulli objective is F¯n(Q,R,π,ϑ)=M,H,Z[Ln(M,H,Z,ϑ)]+λn[∑i≠jqij+∑h=1Lmaxπh(1+∑j=1prhj)]. split F_n(Q,R,π, )=&E_M,H,Z[L_n(M,H,Z; )]+ _n\! [ _i≠ jq_ij+ _h=1^L_ _h (1+ _j=1^pr_hj ) ]. split (21) The penalty is the exact expected complexity because [ZhHhj]=πhrhjE[Z_hH_hj]= _hr_hj. Only configurations with positive probability contribute to the expectation; a singular A(M)A(M) in any such configuration makes the objective infinite. If all qij∈(0,1)q_ij∈(0,1) for i≠ji≠ j, every binary observed-edge mask has positive probability, so A(M)A(M) must be invertible for every such mask. For fixed continuous parameters making every binary A(M)A(M) invertible, the objective is differentiable in the probabilities but need not be jointly convex. Appendix A.11 gives further details on this. Active latent slots with fewer than two children violate the definition of p,LmaxG_p,L_ . For each finite-likelihood configuration, the appendix shows how to remove them without changing the covariance or increasing the score. Proposition 10 (Global exactness of the Bernoulli formulation). Let Fnbin(M,H,Z,ϑ)=Ln(M,H,Z,ϑ)+λnc(M,H,Z) [rgb]0,0,0F_n^bin(M,H,Z; )=L_n(M,H,Z; )+ _nc(M,H,Z). Then infQ,R,π,ϑ[Fnbin(M,H,Z,ϑ)]=mininfϑM,H,ZFnbin(M,H,Z,ϑ). _Q,R,π, E[ [rgb]0,0,0F_n^bin(M,H,Z; )]= _M,H,Z _ [rgb]0,0,0F_n^bin(M,H,Z; ). (22) The common value equals minG+∈p,Lmaxsn(G+)−nplog(2π)/2 _G^+ _p,L_ s_n(G^+)-np (2π)/2. After removing inactive slots and active slots with fewer than two children, both formulations have the same minimizing graphs. If the left infimum is attained, a global minimizer with binary Q, R, and π exists. Remark 11. (Significance of global exactness). The Bernoulli formulation supports gradient-based joint optimization of structure probabilities and continuous parameters. Proposition 10 ensures that introducing these continuous probabilities preserves the global infimum of the corresponding binary objective. This exactness concerns the Bernoulli expectation. Appendix A.11 proves the result. 5 Experiments We evaluate the proposed method on controlled synthetic linear Gaussian SCMs and on simulated GeneNetWeaver (GNW) benchmarks (Schaffter et al., 2011). We compare against Ghassami et al. (Ghassami et al., 2020), DCCD-CONF (Sethuraman and Fekri, 2025), and Améndola et al. (Amendola et al., 2020). All methods evaluated at a shared setting receive the same generated data and trial seeds. 5.1 Setup and evaluation Synthetic experiments use n=10,000n=10,000 observational samples. Observed graphs are generated following the procedure used by Ghassami et al. (2020); for confounded settings, we augment them with randomly generated exogenous latent common causes (Appendix D). Our gradient-based optimization algorithm is explained in Appendix C. We evaluate performance by comparing the covariance families of the reference graph G∗+G^*+ and the recovered graph G^+ G^+. For the forward discrepancy dfwdd_fwd, we generate covariance matrices Σk∗ _k^* using different admissible parameter values of the reference graph. We then fit the recovered graph’s parameters to each covariance matrix while keeping its structure fixed. For the reverse discrepancy drevd_rev, we generate covariance matrices from the recovered graph and fit the reference graph’s parameters to each matrix. Each discrepancy is the largest fitting error across the evaluated covariance matrices. We report dfwd+drevd_fwd+d_rev as the primary bidirectional compatibility error; lower values indicate better agreement at these covariances. Appendix D.3 provides the exact procedure. SHD measures agreement with one particular graph representation, whereas our compatibility error allows structurally different graphs that are distributionally equivalent. An equivalence-aware version of SHD requires characterizing the relevant equivalent graphs or searching over them, which is computationally impractical here. 5.2 Synthetic Results (a) No latents, degree 4. (b) 20%20\% latent, degree 4. (c) 30%30\% latent, degree 4. (d) No latents, degree 2. (e) 20%20\% latent, degree 2. (f) 30%30\% latent, degree 2. Figure 1: Synthetic results across 10 trials. Compatibility error for varying number of observed variables, the maximum observed degree, and the latent-to-observed ratio. Figure 2: Wall-clock runtime. Figure 1 presents synthetic results across different numbers of observed variables, maximum observed degrees, and ratios of latent to observed variables. Our method achieves lower mean compatibility error in several of the larger latent-variable settings where baseline results are available. Each baseline is evaluated up to the largest graph size for which its computation remains feasible. Our method is additionally evaluated at p=256p=256 and p=512p=512 to demonstrate scalability to substantially larger graphs. Appendix D.6 provides more experiments over all degree and latent-to-observed ratio combinations. Figure 2 compares wall-clock runtime. For the baselines, we could not report the wall-clock beyond p=128p=128, because they did not complete within our computational budget. Our method additionally averaged 45444544 seconds at p=256p=256 and 2477024770 seconds at p=512p=512, demonstrating scalability to substantially larger graphs. 5.3 GeneNetWeaver benchmark We additionally evaluate fixed GNW networks with 100 observed variables and either 20 or 30 explicitly constructed exogenous latent variables. GNW generates observations using gene-regulatory dynamics rather than our assumed linear Gaussian model, providing a complementary test under model misspecification while retaining a known cyclic graph structure. Table 1 reports the compatibility error results. Across both latent settings, our method achieves a substantially lower total compatibility error, approximately half that of the closest baseline. It is also markedly faster than all baselines in both latent settings. See Appendix D.5 for more details. Table 1: GeneNetWeaver results with 100 observed variables. Values are mean ± SEM across 10 trials. Latents Method dfwd↓d_fwd drev↓d_rev (dfwd+drev)↓(d_fwd+d_rev) Time (s)↓ 20 Ghassami et al. (2020) 16.57±0.2316.57± 0.23 22.78±2.8722.78± 2.87 39.35±2.7139.35± 2.71 615.40±45.48615.40± 45.48 DCCD-CONF 28.01±1.6728.01± 1.67 10.03±0.4310.03± 0.43 38.05±1.9138.05± 1.91 1334.86±0.551334.86± 0.55 Amendola et al. (2020) 22.97±0.2422.97± 0.24 93.54±7.4393.54± 7.43 116.51±7.65116.51± 7.65 1464.31±166.351464.31± 166.35 Ours 16.33±0.2516.33± 0.25 3.31±0.133.31± 0.13 19.63±0.2419.63± 0.24 47.92±0.0747.92± 0.07 30 Ghassami et al. (2020) 16.28±0.5916.28± 0.59 24.75±1.5924.75± 1.59 41.03±1.8041.03± 1.80 743.09±58.81743.09± 58.81 DCCD-CONF 36.61±4.9336.61± 4.93 11.12±0.9311.12± 0.93 47.73±5.0847.73± 5.08 1352.79±8.691352.79± 8.69 Amendola et al. (2020) 25.46±1.5925.46± 1.59 98.94±2.9098.94± 2.90 124.40±3.57124.40± 3.57 1822.18±308.101822.18± 308.10 Ours 18.27±0.7418.27± 0.74 3.00±0.273.00± 0.27 21.27±0.7321.27± 0.73 122.54±1.87122.54± 1.87 6 Conclusion We developed a framework for learning cyclic linear Gaussian models with exogenous latent confounders and established consistency up to marginal quasi-equivalence under explicit identification conditions. Our Bernoulli formulation for cyclic Gaussian models with explicit latent variables preserves the corresponding binary objective’s global infimum. Our recovery guarantee assumes identification conditions and global score minimization; it does not imply unique graph or latent-count recovery. Broader structural-minimality conditions and numerical optimization guarantees remain open problems. References Akbulut and King (1981) S. Akbulut and H. C. King The topology of real algebraic sets with isolated singularities. Annals of Mathematics 113 (3), p. 425–446. External Links: Document, Link Cited by: §B.1. Amendola et al. (2020) C. Amendola, P. Dettling, M. Drton, F. Onori, and J. Wu Structure learning for cyclic linear causal models. In Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence, Proceedings of Machine Learning Research, Vol. 124, p. 999–1008. Cited by: §B.3, §B.3, §D.4, Appendix E, Appendix E, Appendix E, §1, §3.3, Table 1, Table 1, §5. Boyd and Vandenberghe (2004) S. Boyd and L. Vandenberghe Convex optimization. Cambridge University Press. External Links: Link Cited by: §A.9. Coste (2000) M. Coste An introduction to semialgebraic geometry. Istituti Editoriali e Poligrafici Internazionali, Pisa. External Links: Link Cited by: §A.3, §A.3, §A.3, §A.4, §A.4, §A.5, §A.5, §A.6, §A.6, §A.7, §B.1, §B.1, §B.3, §2.1, Definition 13, Definition 14. Forré and Mooij (2018) P. Forré and J. M. Mooij Constraint-based causal discovery for non-linear structural causal models with cycles and latent confounders. In Proceedings of the Thirty-Fourth Conference on Uncertainty in Artificial Intelligence, p. 269–278. External Links: Link Cited by: Appendix E, §1. Ghassami et al. (2020) A. Ghassami, A. Yang, N. Kiyavash, and K. Zhang Characterizing distribution equivalence and structure learning for cyclic and acyclic directed graphs. In Proceedings of the 37th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 119, p. 3494–3504. Cited by: §A.6, §D.1, §D.4, Appendix E, 1st item, §1, §3.1, §5.1, Table 1, Table 1, §5. Horn and Johnson (2013) R. A. Horn and C. R. Johnson Matrix analysis. 2 edition, Cambridge University Press. External Links: Document Cited by: §A.4. Hyttinen et al. (2012) A. Hyttinen, F. Eberhardt, and P. O. Hoyer Learning linear cyclic causal models with latent variables. Journal of Machine Learning Research 13 (109), p. 3387–3439. External Links: Link Cited by: Appendix E, §1. Kaltenpoth and Vreeken (2023) D. Kaltenpoth and J. Vreeken Causal discovery with hidden confounders using the algorithmic Markov condition. In Proceedings of the Thirty-Ninth Conference on Uncertainty in Artificial Intelligence, Proceedings of Machine Learning Research, Vol. 216, p. 1016–1026. Cited by: Appendix E, §1. Ma et al. (2024) P. Ma, R. Ding, Q. Fu, J. Zhang, S. Wang, S. Han, and D. Zhang Scalable differentiable causal discovery in the presence of latent confounders with skeleton posterior. In Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, p. 2148–2159. External Links: Document Cited by: Appendix E, §1. Maddison et al. (2017) C. J. Maddison, A. Mnih, and Y. W. Teh The concrete distribution: a continuous relaxation of discrete random variables. In International Conference on Learning Representations, Cited by: Appendix C, Appendix E. Schaffter et al. (2011) T. Schaffter, D. Marbach, and D. Floreano GeneNetWeaver: in silico benchmark generation and performance profiling of network inference methods. Bioinformatics 27 (16), p. 2263–2270. External Links: Document Cited by: §5. Sethuraman and Fekri (2025) M. G. Sethuraman and F. Fekri Differentiable cyclic causal discovery under unmeasured confounders. In Advances in Neural Information Processing Systems, Vol. 38. Cited by: §D.4, Appendix E, §1, §5. Yi and Lee (2024) S. Yi and S. Lee Filter, rank, and prune: learning linear cyclic gaussian graphical models. In Proceedings of the 27th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, Vol. 238, p. 1135–1143. Cited by: Appendix E, §1. Zheng et al. (2018) X. Zheng, B. Aragam, P. K. Ravikumar, and E. P. Xing DAGs with NO TEARS: continuous optimization for structure learning. In Advances in Neural Information Processing Systems, Vol. 31. Cited by: Appendix E. Łaba and Solymosi (2007) I. Łaba and J. Solymosi Incidence theorems for pseudoflats. Discrete & Computational Geometry 37 (2), p. 163–174. External Links: Document, Link Cited by: Definition 12, Definition 15. Appendix A Proofs A.1 Proof of Lemma 1 Proof. Substituting the block matrices in Equation 2 into the structural equations gives the two block equations XO=BOO⊤XO+Λ⊤XL+NO,XL=NL.X_O=B_O X_O+ X_L+N_O, X_L=N_L. The second identity follows because latent variables have no parents. Using XL=NLX_L=N_L in the observed equation and collecting the terms involving XOX_O on the left gives (Ip−BOO⊤)XO=Λ⊤NL+NO.(I_p-B_O )X_O= N_L+N_O. By assumption, Ip−BOOI_p-B_O is invertible. Therefore, XO=(Ip−BOO)−⊤(Λ⊤NL+NO).X_O=(I_p-B_O)^- ( N_L+N_O). (23) The vector on the right is Gaussian and has mean zero. Moreover, the block diagonality of Ω implies that NON_O and NLN_L are independent, so the two cross-covariance terms vanish. Consequently, Cov(Λ⊤NL+NO) ( N_L+N_O) =Λ⊤Cov(NL)Λ+Cov(NO) = Cov(N_L) +Cov(N_O) =Λ⊤ΩLΛ+ΩO=Ψ. = _L + _O= . Taking the covariance of both sides of Equation 23 now yields Cov(XO)=(Ip−BOO)−⊤Ψ(Ip−BOO)−1,Cov(X_O)=(I_p-B_O)^- (I_p-B_O)^-1, which is Equation 3. It remains to verify positive definiteness. Let ΩL1/2 _L^1/2 be the diagonal matrix obtained by taking the positive square root of each latent noise variance, and set C=ΩL1/2Λ∈ℝℓ×pC= _L^1/2 × p. The latent covariance contribution has the factorization Λ⊤ΩLΛ=(ΩL1/2Λ)⊤(ΩL1/2Λ)=C⊤C. _L =( _L^1/2 ) ( _L^1/2 )=C C. For every v∈ℝpv ^p, v⊤C⊤Cv=(Cv)⊤(Cv)≥0v C Cv=(Cv) (Cv)≥ 0, so this contribution is positive semidefinite. Since C has ℓ rows and p columns, rank(Λ⊤ΩLΛ)=rank(C⊤C)≤rank(C)≤min(ℓ,p).rank( _L )=rank(C C) (C)≤ ( ,p). Because ΩO≻0 _O 0, every nonzero v∈ℝpv ^p satisfies v⊤Ψv=v⊤ΩOv+(Cv)⊤(Cv)>0.v v=v _Ov+(Cv) (Cv)>0. Thus Ψ≻0 0, regardless of the rank of Λ . When ℓ=0 =0, the latent contribution and its rank are zero, and Ψ=ΩO≻0 = _O 0. Finally, write A=Ip−BOOA=I_p-B_O. Its invertibility implies A−1v≠0A^-1v≠ 0 whenever v≠0v≠ 0, and hence v⊤ΣOv=(A−1v)⊤Ψ(A−1v)>0.v _Ov=(A^-1v) (A^-1v)>0. Therefore ΣO≻0 _O 0, completing the proof. ∎ A.2 Derivation of the observed-data score We first derive the full Gaussian likelihood in Equation 9, then express it in structural parameters and establish the KL representation of the graph score in Equation 13. Let O∈ℝn×pX_O ^n× p have kkth row XO(k)⊤X_O^(k) . Throughout, the mean is known to be zero, so O⊤O=nSnX_O X_O=nS_n without sample-mean subtraction. Under the nonsingular Gaussian sampling model, OX_O has full column rank almost surely when n≥pn≥ p, and hence Sn≻0S_n 0. Gaussian likelihood. For a candidate covariance Σ≻0 0, the density of one observation x∈ℝpx ^p is fΣ(x)=(2π)−p/2(detΣ)−1/2exp(−12x⊤Σ−1x).f_ (x)=(2π)^-p/2( )^-1/2 \! (- 12x ^-1x ). Independence of the observations makes their joint density a product. Taking its negative logarithm gives ℒnfull(Σ) _n^full( ) =−∑k=1nlogfΣ(XO(k)) =- _k=1^n f_ (X_O^(k)) =np2log(2π)+n2logdetΣ+12∑k=1nXO(k)⊤Σ−1XO(k). = np2 (2π)+ n2 + 12 _k=1^nX_O^(k) ^-1X_O^(k). Each observation contributes plog(2π)/2p (2π)/2 through the Gaussian normalization factor, giving nplog(2π)/2np (2π)/2 in total. This term depends only on n and the observed dimension p, not on the graph, the latent count, or any continuous parameter. Using x⊤Σ−1x=tr(Σ−1xx⊤)x ^-1x=tr( ^-1x ) gives ℒnfull(Σ)=np2log(2π)+n2logdetΣ+12tr(Σ−1O⊤O).L_n^full( )= np2 (2π)+ n2 + 12tr( ^-1X_O X_O). Substituting O⊤O=nSnX_O X_O=nS_n yields Equation 9. As in the main text, we also write ℒnfull(θ)=ℒnfull(ΣO(θ))L_n^full(θ)=L_n^full( _O(θ)). Likelihood in structural parameters. For a fixed graph, write A=Ip−BOOA=I_p-B_O and Ψ=ΩO+Λ⊤ΩLΛ = _O+ _L . By Equation 3, ΣO(θ)=A−⊤ΨA−1,ΣO(θ)−1=AΨ−1A⊤, _O(θ)=A^- A^-1, _O(θ)^-1=A ^-1A , and logdetΣO(θ)=logdetΨ−2log|detA|. _O(θ)= -2 | A|. Cyclic invariance of the trace gives tr[ΣO(θ)−1O⊤O] \! [ _O(θ)^-1X_O X_O ] =tr[Ψ−1A⊤O⊤OA] =tr\! [ ^-1A X_O X_OA ] =tr[Ψ−1(OA)⊤(OA)]. =tr\! [ ^-1(X_OA) (X_OA) ]. Removing the Gaussian constant does not change the minimizing graphs or parameters. Define the constant-free negative log-likelihood by ℒn(θ):=ℒnfull(θ)−np2log(2π)=−nlog|detA|+n2logdetΨ+12tr[Ψ−1(OA)⊤(OA)]. splitL_n(θ)&:=L_n^full(θ)- np2 (2π)\\ &=-n | A|+ n2 + 12tr\! [ ^-1(X_OA) (X_OA) ]. split (24) The second equality follows by substituting the preceding identities into the full likelihood. Equivalently, the last term is n2tr(Ψ−1A⊤SnA) n2tr( ^-1A S_nA), the form used for the masked likelihood in Equation 31. For one observation x, the vector r=A⊤x=x−BOO⊤xr=A x=x-B_O x is its structural residual: what remains after subtracting the observed-parent contributions in the structural equations. Under the candidate SCM, A⊤XO=Λ⊤NL+NOA X_O= N_L+N_O has covariance Ψ . The rows of OAX_OA are therefore the transposes of these residuals. The trace term sums their Gaussian quadratic terms r⊤Ψ−1r/2r ^-1r/2, and n2logdetΨ n2 comes from the Gaussian normalization. The linear change of variables x↦A⊤x A x multiplies volume by |detA|| A|. Consequently, converting the residual density to the observation density contributes −log|detA|- | A| to each negative log-likelihood, giving −nlog|detA|-n | A| in total. These identities require only invertibility of A, not acyclicity. The latent covariance contribution remains Λ⊤ΩLΛ _L , with the loading restrictions specified by the graph; Ψ is not treated as an unrestricted covariance matrix. KL representation. Expanding Equation 12 gives nD(Sn∥Σ)=n2tr(Σ−1Sn)+logdetΣ−logdetSn−p.nD(S_n\| )= n2\tr( ^-1S_n)+ - S_n-p\. Evaluating the full likelihood at Σ=Sn =S_n gives the constant Cn=ℒnfull(Sn)=n2plog(2π)+logdetSn+p.C_n=L_n^full(S_n)= n2\p (2π)+ S_n+p\. Consequently, ℒnfull(Σ)=Cn+nD(Sn∥Σ).L_n^full( )=C_n+nD(S_n\| ). Here CnC_n includes the Gaussian normalization term and the value of the remaining likelihood terms at SnS_n. It is independent of the graph and its parameters, so optimizing the likelihood is equivalent to minimizing D(Sn∥Σ)D(S_n\| ), the KL divergence between the two zero-mean Gaussian distributions with covariances SnS_n and Σ . Equality of the likelihood infima. For every admissible parameter tuple, ℒn(θ)=Kn+nD(Sn∥ΣO(θ)),Kn=n2logdetSn+p.L_n(θ)=K_n+nD (S_n\| _O(θ) ), K_n= n2\ S_n+p\. By definition, ℳ(G+)M(G^+) is exactly the image of (G+)P(G^+) under the covariance map. Thus taking the infimum over admissible parameters equals taking the infimum over ℳ(G+)M(G^+). To pass to its relative closure, first note that set inclusion gives infΣ∈ℳ¯(G+)D(Sn∥Σ)≤infΣ∈ℳ(G+)D(Sn∥Σ). _ ∈ M(G^+)D(S_n\| )≤ _ (G^+)D(S_n\| ). Conversely, for any Σ¯∈ℳ¯(G+) ∈ M(G^+), there is a sequence Σk∈ℳ(G+) _k (G^+) converging to Σ¯≻0 0. Continuity of D(Sn∥⋅)D(S_n\|·) on ++pS_++^p implies infΣ∈ℳ(G+)D(Sn∥Σ)≤limk→∞D(Sn∥Σk)=D(Sn∥Σ¯). _ (G^+)D(S_n\| )≤ _k→∞D(S_n\| _k)=D(S_n\| ). Taking the infimum over Σ¯∈ℳ¯(G+) ∈ M(G^+) proves the reverse inequality. Consequently, infθ∈(G+)ℒn(θ)=Kn+ninfΣ∈ℳ¯(G+)D(Sn∥Σ)=Kn+nδn(G+). _θ (G^+)L_n(θ)=K_n+n _ ∈ M(G^+)D(S_n\| )=K_n+n _n(G^+). (25) Adding λnc(G+) _nc(G^+) and the omitted Gaussian constant gives infθ∈(G+)Jn(G+,θ)=sn(G+) _θ (G^+)J_n(G^+,θ)=s_n(G^+) as defined in Equation 11, and proves its equivalent expression in Equation 13. No attainment of the infimum by admissible SCM parameters is required. Every graph admits the covariance IpI_p, obtained by setting all edge coefficients to zero and all noise variances to one. Hence 0≤δn(G+)≤D(Sn∥Ip)<∞0≤ _n(G^+)≤ D(S_n\|I_p)<∞. Since p,LmaxG_p,L_ is finite, the graph score has at least one global minimizer even if a particular graph’s parameter infimum is not attained. The choice λn=12logn _n= 12 n gives BIC-like scaling. The consistency proof uses only λn→∞ _n→∞ and λn/n→0 _n/n→ 0: the penalty eventually distinguishes unnecessarily complex candidates, while a fixed positive population KL gap contributes an order-n likelihood disadvantage that dominates any fixed complexity advantage. This argument does not identify c(G+)c(G^+) with the covariance-family dimension or require the regularity conditions of a standard BIC derivation. A.3 Algebraic and semialgebraic sets and their dimensions This subsection gives the formal definitions underlying Section 2.1. We work with subsets of ℝmR^m and polynomials with real coefficients; ℝ[x1,…,xm]R[x_1,…,x_m] denotes the set of such polynomials. Symmetric matrices are identified with vectors of their upper-triangular entries. We first specify the topological terms used in the proofs. The algebraic and dimension definitions that follow include references to their sources. We temporarily use dimsa _sa and dimalg _alg to distinguish the two notions of dimension. Open sets, interior, and neighborhoods. For a∈ℝma ^m and ε>0 >0, the open ball of radius ε centered at a is Bε(a)=x∈ℝm:‖x−a‖2<εB_ (a)=\x ^m:\|x-a\|_2< \, where ‖x‖2=(∑ixi2)1/2\|x\|_2=( _ix_i^2)^1/2 is the Euclidean norm of a vector. A set U⊆ℝmU ^m is Euclidean open if, for every x∈Ux∈ U, there exists ε>0 >0 such that Bε(x)⊆UB_ (x) U. Thus openness requires a ball around every point of U, not merely one ball somewhere in U. For S⊆V⊆ℝmS V ^m, the interior of S relative to V is intVS:=x∈S:there exists ε>0 such that Bε(x)∩V⊆S.int_VS:=\x∈ S:there exists >0 such that B_ (x)∩ V S\. A set S⊆VS V is open relative to V if intVS=Sint_VS=S. When V=ℝmV=R^m, these definitions give the ordinary Euclidean interior and openness. A neighborhood of x∈Vx∈ V relative to V is a subset of V containing Bε(x)∩VB_ (x)∩ V for some ε>0 >0; the neighborhood itself need not be open. A set S is closed relative to V if clVS=Scl_VS=S, using the closure defined in Section 2.1. For symmetric matrices, we use Frobenius-norm balls; they define the same open sets as Euclidean balls in the upper-triangular coordinates. For example, [0,1][0,1] is not open in ℝR, but its interior is (0,1)(0,1), which contains an open ball. Also, the line segment S=(t,0):−1<t<1S=\(t,0):-1<t<1\ has empty interior in ℝ2R^2 but satisfies intVS=Sint_VS=S for V=(t,0):t∈ℝV=\(t,0):t \. Throughout the proofs, the specified space V determines which interior or neighborhood is meant. Definition 12 (Real algebraic sets and irreducibility). We use the real-set definitions given by Łaba and Solymosi (2007, Section 2). For finitely many polynomials f1,…,fr∈ℝ[x1,…,xm]f_1,…,f_r [x_1,…,x_m], define their common zero set by Z(f1,…,fr):=x∈ℝm:fi(x)=0 for every i=1,…,r.Z(f_1,…,f_r):=\x ^m:f_i(x)=0 for every i=1,…,r\. A set V⊆ℝmV ^m is real algebraic if V=Z(f1,…,fr)V=Z(f_1,…,f_r) for some such finite collection. A nonempty real algebraic set V is irreducible if, for any real algebraic subsets V1,V2⊆V_1,V_2 V, V=V1∪V2⟹V=V1 or V=V2.V=V_1∪ V_2 V=V_1 or V=V_2. Thus an irreducible set cannot be decomposed into two strictly smaller algebraic sets. A polynomial f vanishes on V if f(x)=0f(x)=0 for every x∈Vx∈ V. Definition 13 (Semialgebraic sets). We use the finite-union formulation in Coste (2000, Section 2.1.1, Proposition 2.1). A set S⊆ℝmS ^m is semialgebraic if it can be written as a finite union S=⋃a=1Nx∈ℝm:fai(x)=0(1≤i≤ra),gaj(x)>0(1≤j≤sa),S= _a=1^N\x ^m:f_ai(x)=0\ (1≤ i≤ r_a), g_aj(x)>0\ (1≤ j≤ s_a)\, where all faif_ai and gajg_aj are real polynomials. The integers N,ra,saN,r_a,s_a are nonnegative: an empty union denotes the empty set, and an empty list of conditions imposes no restriction. Weak inequalities are also permitted through the identity g≥0=g>0∪g=0\g≥ 0\=\g>0\∪\g=0\. Every algebraic set is semialgebraic. The converse need not hold, because inequalities can restrict a semialgebraic set to only part of an algebraic set. A map f:S→ℝkf:S ^k is called semialgebraic if its graph (x,f(x)):x∈S\(x,f(x)):x∈ S\ is a semialgebraic subset of ℝm+kR^m+k, where S⊆ℝmS ^m (Coste, 2000, Chapter 2). Definition 14 (Semialgebraic dimension). We use the coordinate-projection characterization of semialgebraic dimension; see Coste (2000, Section 3.3) for the standard dimension theory. For I=i1<⋯<ik⊆1,…,mI=\i_1<·s<i_k\ \1,…,m\, let πI:ℝm→ℝk _I:R^m ^k be the coordinate projection πI(x)=(xi1,…,xik) _I(x)=(x_i_1,…,x_i_k). For nonempty semialgebraic S⊆ℝmS ^m, define dimsaS:=max|I|:πI(S) contains a nonempty open ball in ℝ|I|. _saS:= \|I|: _I(S) contains a nonempty open ball in R^|I| \. For I=∅I= , the projection of a nonempty set is the singleton ℝ0R^0, so dimension zero is included. We use the convention dimsa∅=−1 _sa =-1. Thus dimension is the largest number of selected coordinates whose projection contains an open ball; it need not equal the total number of coordinates used to describe the set. Definition 15 (Algebraic dimension). For a nonempty real algebraic set V⊆ℝmV ^m, its algebraic dimension can be defined by chains of nested algebraic subsets (Łaba and Solymosi, 2007, Section 2): dimalgV:=maxd∈ℤ≥0:Z0⊊Z1⊊⋯⊊Zd⊆V,Z0,…,Zd are nonempty irreducible real algebraic sets. split _algV:= \d _≥ 0:&\;Z_0 Z_1 ·s Z_d V,\\ &\;Z_0,…,Z_d are nonempty irreducible real algebraic sets \. split The integer d counts strict inclusions, not sets. We set dimalg∅=−1 _alg =-1. If V is irreducible, a chain of maximum length ends at V; for a reducible set, the final set ZdZ_d need only be contained in V. Relation between the dimensions. For any S⊆ℝmS ^m, its real Zariski closure is the common zero set of all polynomials that are zero throughout S (Coste, 2000, Section 3.3): S¯Zar:=x∈ℝm:f(x)=0 for every real polynomial f satisfying f(s)=0 for all s∈S. S^\,Zar:=\x ^m:f(x)=0 for every real polynomial f satisfying f(s)=0 for all s∈ S\. Equivalently, it is the smallest real algebraic set containing S. For every nonempty semialgebraic set S, the standard dimension theorem gives (Coste, 2000, Theorem 3.20) dimsaS=dimalgS¯Zar. _saS= _alg S^\,Zar. In particular, if S is already algebraic, its Zariski closure equals S and the two dimensions agree. This justifies the common notation dim used throughout the paper. For semialgebraic T⊆ST S, full-dimensional in S means dimT=dimS T= S; it does not mean dimT=m T=m. For example, the open line segment S=(x,0)∈ℝ2:0<x<1S=\(x,0) ^2:0<x<1\ is semialgebraic and has dimension one: its first-coordinate projection is (0,1)(0,1), but S contains no open ball in ℝ2R^2. Its Zariski closure is the entire line V=(x,0):x∈ℝV=\(x,0):x \, which also has dimension one. Indeed, a polynomial that is zero on the segment restricts to a univariate polynomial that is zero on (0,1)(0,1), and must therefore be zero on the entire line. Thus S is full-dimensional in V even though S≠VS≠ V. A.4 Closure properties of covariance families We explain the containment and closure properties used in the model geometry proof. The first two properties follow from the definition of the real Zariski closure and continuity of polynomials; they do not require semialgebraicity or irreducibility. 1. The covariance family is contained in its Zariski closure. By definition, V(G+)V(G^+) consists of the symmetric matrices satisfying every polynomial equality that holds throughout ℳ(G+)M(G^+). Each matrix in ℳ(G+)M(G^+) satisfies all these equalities, so ℳ(G+)⊆V(G+)M(G^+) V(G^+). This inclusion is a consequence of the definition, not an additional model assumption. 2. Positive-definite limits are also in the Zariski closure. Let Σ¯∈ℳ¯(G+) ∈ M(G^+). By the definition of relative Euclidean closure, there is a sequence Σk∈ℳ(G+) _k (G^+) such that Σk→Σ¯≻0 _k→ 0. If a polynomial f in the upper-triangular matrix entries equals zero throughout ℳ(G+)M(G^+), then f(Σk)=0f( _k)=0 for every k. Since polynomials are continuous, f(Σ¯)=limk→∞f(Σk)=0.f( )= _k→∞f( _k)=0. Thus Σ¯ satisfies every polynomial equality defining V(G+)V(G^+) and belongs to that set. Also, every element of ℳ(G+)M(G^+) belongs to its relative closure, as witnessed by the constant sequence at that element. Together these facts give ℳ(G+)⊆ℳ¯(G+)⊆V(G+)∩++p.M(G^+) M(G^+) V(G^+) _++^p. The same continuity argument shows that any real algebraic set contains the limits of its convergent sequences, and is therefore Euclidean closed. Containment in V(G+)V(G^+) alone would not justify the inclusion of limits; its Euclidean closedness is the additional fact used here. 3. The relative closure is semialgebraic. Step 2 of the proof of Proposition 2 below establishes that ℳ(G+)M(G^+) is semialgebraic. Its Euclidean closure in pS^p is therefore semialgebraic by Coste (2000, Corollary 2.5), a consequence of the Tarski–Seidenberg theorem. By Sylvester’s criterion, positive definiteness can be expressed by finitely many strict polynomial inequalities in the matrix entries (Horn and Johnson, 2013, Section 7.2, Theorem 7.2.5). Hence ++pS_++^p is semialgebraic. Finally, finite intersections of semialgebraic sets are semialgebraic (Coste, 2000, Section 2.1.2). Since relative closure within ++pS_++^p keeps exactly the positive-definite matrices in the closure within pS^p, ℳ¯(G+)=(clpℳ(G+))∩++p M(G^+)= (cl_S^pM(G^+) ) _++^p is semialgebraic. A.5 Proof of Proposition 2 Proof. Fix G+G^+ throughout the proof. We first describe the parameter space and its map to covariance matrices, and then establish the three claims. 1. Free parameters and the covariance map. Impose the graph-required zeros on BOOB_O and Λ . List the remaining allowed coefficients and the diagonal entries of ΩO _O and ΩL _L in a vector θ∈ℝdθ ^d. Here, d is the number of free parameter coordinates, not necessarily the dimension of the covariance family. The parameter space ℝdR^d is distinct from the covariance space pS^p, which we identify with ℝp(p+1)/2R^p(p+1)/2 through the upper-triangular matrix entries. Each admissible parameter vector θ determines one covariance matrix; different admissible parameter vectors may determine the same one. Write A(θ)=Ip−BOO,Δ(θ)=detA(θ),Ψ(θ)=ΩO+Λ⊤ΩLΛ.A(θ)=I_p-B_O, (θ)= A(θ), (θ)= _O+ _L . In these free coordinates, (G+)P(G^+) is defined by (ΩO)ii>0,(ΩL)hh>0,Δ(θ)2>0.( _O)_i>0, ( _L)_h>0, (θ)^2>0. Each variance is a coordinate of θ, and Δ(θ) (θ) is a polynomial because the determinant is polynomial in the matrix entries. These finitely many polynomial inequalities therefore show that (G+)P(G^+) is semialgebraic. All inequalities are strict, and their left-hand sides are continuous. Hence, for every admissible θ, some ε>0 >0 satisfies Bε(θ)⊆(G+)B_ (θ) (G^+): the parameter set is Euclidean open in the free-coordinate space ℝdR^d. It contains the point θ0 _0 with BOO=0B_O=0, Λ=0 =0, and unit noise variances, and therefore contains an open ball centered at θ0 _0. Graph-required zero entries are fixed throughout; only the free coordinates vary. This uses the convention that an allowed edge coefficient may equal zero. The covariance map is ϕG+:(G+)⟶++p,ϕG+(θ)=A(θ)−⊤Ψ(θ)A(θ)−1. _G^+:P(G^+) _++^p, _G^+(θ)=A(θ)^- (θ)A(θ)^-1. Its values are positive definite by Lemma 1, and its image is exactly ℳ(G+)M(G^+). 2. The covariance family is semialgebraic. Consider the joint set of admissible parameters and their generated covariance matrices: =(θ,Σ)∈ℝd×p:θ∈(G+),A(θ)⊤ΣA(θ)=Ψ(θ).T= \(θ, ) ^d×S^p:θ (G^+),\;A(θ) A(θ)= (θ) \. The matrix equality represents finitely many scalar polynomial equalities in the entries of θ and Σ . Together with the polynomial inequalities defining (G+)P(G^+), these show that T is semialgebraic. Since A(θ)A(θ) is invertible on (G+)P(G^+), the matrix equality is equivalent to Σ=ϕG+(θ) = _G^+(θ). Therefore, ℳ(G+)=Σ∈p:there exists θ with (θ,Σ)∈.M(G^+)=\ ^p:there exists θ with (θ, ) \. This is the coordinate projection of T onto the covariance space. The Tarski–Seidenberg theorem states that a coordinate projection of a semialgebraic set is semialgebraic (Coste, 2000, Theorem 2.3 and Corollary 2.4). It therefore proves that ℳ(G+)M(G^+) is semialgebraic. The theorem is needed because the polynomial conditions defining T involve both parameters and covariances, whereas ℳ(G+)M(G^+) is a set of covariance matrices alone. No explicit elimination of the parameters is required. 3. The Zariski closure is irreducible. Write V=V(G+)V=V(G^+). Irreducibility means that V cannot be written as the union of two proper algebraic subsets. Suppose, to obtain a contradiction, that V=V1∪V2V=V_1∪ V_2, where V1≠V_1≠ V and V2≠V_2≠ V are algebraic. Choose a matrix Σ1∈V∖V1 _1∈ V V_1. Because V1V_1 is the common zero set of a collection of polynomials in the symmetric matrix entries, at least one defining polynomial, denoted by f, satisfies f(Σ1)≠0f( _1)≠ 0. Thus f is zero throughout V1V_1 but not throughout V. Similarly, choose a polynomial g that is zero throughout V2V_2 but not throughout V. These are polynomials in covariance-space coordinates, not in the structural parameters. Saying that a polynomial vanishes on a set means that it equals zero at every element of that set. Every generated covariance belongs to V1V_1 or V2V_2. Consequently, f(ϕG+(θ))g(ϕG+(θ))=0for every θ∈(G+).f ( _G^+(θ) )g ( _G^+(θ) )=0 every θ (G^+). To turn this into a polynomial identity in the free parameters, use the adjugate formula. Each entry of the adjugate J(θ)=adjA(θ)J(θ)=adjA(θ) is, up to sign, the determinant of a matrix obtained by deleting one row and one column of A(θ)A(θ). Determinants are polynomials in the matrix entries, so every entry of J(θ)J(θ) is a polynomial in θ. Since Δ(θ)=detA(θ)≠0 (θ)= A(θ)≠ 0 on the admissible parameter set, A(θ)−1=J(θ)Δ(θ),ϕG+(θ)=J(θ)⊤Ψ(θ)J(θ)Δ(θ)2.A(θ)^-1= J(θ) (θ), _G^+(θ)= J(θ) (θ)J(θ) (θ)^2. Both J(θ)J(θ) and Ψ(θ) (θ) have polynomial entries, so the numerator entries are polynomials in θ. Each covariance entry is therefore a ratio of polynomials, called a rational function. The polynomials f and g are still polynomials in covariance entries; substituting these ratios makes their values rational functions of θ. Putting each expression over a common denominator gives f(ϕG+(θ))=Pf(θ)Δ(θ)2r,g(ϕG+(θ))=Pg(θ)Δ(θ)2s,f ( _G^+(θ) )= P_f(θ) (θ)^2r, g ( _G^+(θ) )= P_g(θ) (θ)^2s, where Pf,PgP_f,P_g are polynomials in the free parameters and r=degfr= f, s=deggs= g. Since Δ(θ)≠0 (θ)≠ 0 on the admissible parameter set, clearing denominators yields Pf(θ)Pg(θ)=0for every θ∈(G+).P_f(θ)P_g(θ)=0 every θ (G^+). A real polynomial that equals zero on an open ball in ℝdR^d is identically zero, meaning that all its coefficients are zero. To see why, the ball contains a product of nonempty open intervals. Fixing all coordinates but one gives a polynomial in one variable that is zero throughout an interval. A nonzero polynomial in one variable has only finitely many roots, so all its coefficients must be zero. Applying this argument successively to the remaining coordinates shows that every coefficient of the original polynomial is zero. Step 1 therefore implies PfPg≡0P_fP_g≡ 0 on ℝdR^d. A product of two nonzero real polynomials cannot be the zero polynomial. Thus Pf≡0P_f≡ 0 or Pg≡0P_g≡ 0. It follows that one fixed polynomial, either f or g, is zero on every generated covariance. By the definition of Zariski closure, any polynomial that is zero on ℳ(G+)M(G^+) is also zero on V(G+)V(G^+). This contradicts the choice of f and g. Hence V(G+)V(G^+) is irreducible. The argument does not require a covariance matrix to determine a unique parameter vector. 4. The dimension identities. Step 2 and Appendix A.4 show that both ℳ(G+)M(G^+) and ℳ¯(G+) M(G^+) are semialgebraic, and that ℳ(G+)⊆ℳ¯(G+)⊆V(G+)M(G^+) M(G^+) V(G^+). For a semialgebraic set, Coste (2000, Theorem 3.20) identifies its semialgebraic dimension with the algebraic dimension of its real Zariski closure. Thus dimℳ(G+)=dimV(G+) (G^+)= V(G^+). The same theorem identifies the two dimensions of the algebraic set V(G+)V(G^+). Dimension is monotone under inclusion: any coordinate projection of a smaller set is contained in the same projection of the larger set, so an open ball in the former also belongs to the latter. Consequently, dimℳ(G+)≤dimℳ¯(G+)≤dimV(G+)=dimℳ(G+). (G^+)≤ M(G^+)≤ V(G^+)= (G^+). The dimensions must therefore all be equal, proving Equation 6. ∎ A.6 Rationale and scope of marginal quasi-equivalence Relation to the definition of Ghassami et al. Ghassami et al. (2020, Section 5, Definition 9) study linear Gaussian models without latent confounders. They define quasi-equivalence through a positive-measure set of distributions that both graphs can generate, with measure expressed in parameters describing those distributions. Their phrase “linearly independent parameters” refers to free distributional coordinates, not statistically independent random variables. Their definition is therefore phrased in distribution coordinates, rather than as an intersection of the graphs’ full structural parameter spaces. The motivation is that two graphs may share a substantial set of distributions without generating exactly the same distribution family. This positive-overlap motivation is also the basis of our definition. Why we formulate the definition in observed covariance space. In our setting, candidates may have different numbers of latent variables and different structural parameter spaces. Moreover, different structural parameter values may generate the same observed covariance. To compare graphs, we therefore use the sets of observed covariances themselves, ℳ(G1+)M(G_1^+) and ℳ(G2+)M(G_2^+), in the common space ++pS_++^p. Their dimensions count the independent covariance coordinates, rather than the number of structural parameters. Equation 7 makes the positive-overlap idea precise without choosing separate distribution parameterizations for the marginal models. This is a formulation of that idea for our setting, not an assumption that the characterizations or recovery results for models without latent variables automatically extend to models with latent variables. A precision matrix is the inverse of a covariance matrix. Using covariance rather than precision matrices does not change the dimension criterion: on ++pS_++^p, matrix inversion is one-to-one and is its own inverse. Its entries are rational functions on this semialgebraic domain, so it preserves semialgebraic dimension (Coste, 2000, Theorem 3.18). It also sends the intersection of two covariance families to the intersection of their corresponding precision families. Thus both the family dimensions and the intersection dimension are unchanged. Dimension and positive-measure overlap. Let S=ℳ(G1+)∩ℳ(G2+)S=M(G_1^+) (G_2^+), and suppose both families have dimension d. Proposition 2 ensures that the families, and hence their intersection, are semialgebraic. By Definition 14, dimS=d S=d means that some projection of S onto d selected covariance coordinates contains a nonempty open set in ℝdR^d. Such a set has positive d-dimensional Lebesgue measure, meaning positive ordinary volume in those coordinates. Conversely, a semialgebraic subset of ℝdR^d with positive Lebesgue measure contains a nonempty open set (Coste, 2000, Section 3.3). Thus the criterion requires the shared covariances to have a coordinate projection containing an open ball in ℝdR^d, rather than only a subset of dimension below d. We require the intersection to have the dimension of both families: a lower-dimensional family contained in a larger one does not qualify. The criterion does not require equality of the complete families; the parts that they do not share may also have dimension d. Why comparisons must be made directly. For arbitrary semialgebraic sets, full-dimensional overlap need not pass through an intermediate set. For example, consider the one-dimensional intervals A=(0,2),B=(1,4),C=(3,5).A=(0,2), B=(1,4), C=(3,5). Here A∩B=(1,2)A∩ B=(1,2) and B∩C=(3,4)B∩ C=(3,4) both have dimension one, but A∩C=∅A∩ C= . Thus overlap of A with B, and of B with C, does not imply overlap of A with C. This example explains why transitivity, meaning this implication between three sets, cannot be deduced from semialgebraicity alone. It is not a construction of three covariance families in our graph class. Our proofs do not need such an implication: they establish quasi-equivalence between each selected graph and the data-generating graph directly. In Lemma 17, replacing G†G by the true graph is justified by the equality ℳ(G†)=ℳ(G∗+)M(G )=M(G^*+), not by transitivity of quasi-equivalence. Why equal Zariski closures are insufficient. The same intervals also illustrate a separate distinction. All three have Zariski closure ℝR, since a polynomial that equals zero on a nonempty interval equals zero on the whole line. Nevertheless, A and C are disjoint. Consequently, equality of Zariski closures alone does not guarantee overlap for semialgebraic sets: it establishes that the sets satisfy the same polynomial equalities, but does not determine which points satisfy their inequalities. This is why Assumption 9 explicitly requires full-dimensional intersection of the covariance families, even after their Zariski closures have been shown to be equal. A.7 Sufficient conditions for geometric identification We first state the standard algebraic dimension fact used below; see Coste (2000, Section 3.3, p. 59). We include its proof to make the strict-dimension step explicit. Lemma 16 (Strict inclusion and algebraic dimension). Let V1,V2⊆ℝmV_1,V_2 ^m be nonempty irreducible real algebraic sets. If V1⊊V2V_1 V_2, then dimV1<dimV2 V_1< V_2. Proof. By Definition 15, algebraic dimension is the largest number of strict inclusions in a chain of nonempty irreducible algebraic subsets. Let d=dimV1d= V_1. Since V1V_1 is irreducible, a chain of maximal length in V1V_1 ends at V1V_1; otherwise, appending V1V_1 would give a longer chain. Thus there exist nonempty irreducible algebraic sets with Z0⊊Z1⊊⋯⊊Zd=V1.Z_0 Z_1 ·s Z_d=V_1. The assumed strict inclusion and the irreducibility of V2V_2 allow us to append one more set: Z0⊊Z1⊊⋯⊊Zd=V1⊊V2.Z_0 Z_1 ·s Z_d=V_1 V_2. This is a chain of d+1d+1 strict inclusions in V2V_2, so dimV2≥d+1>dimV1 V_2≥ d+1> V_1. The same argument applies to any nonempty proper algebraic subset W⊊V2W V_2, even if W is not irreducible. A chain of maximal length d=dimWd= W consists of irreducible sets and ends at some Zd⊆WZ_d W. Since Zd⊊V2Z_d V_2, appending V2V_2 again gives dimV2≥d+1 V_2≥ d+1. Thus every nonempty proper algebraic subset of an irreducible algebraic set has strictly smaller dimension. ∎ Lemma 17 (Geometric identification). Under Assumptions 7–9, Assumption 5 holds. Equivalently, for every G+∈p,LmaxG^+ _p,L_ , Σ∗∈ℳ¯(G+),c(G+)≤c†⟹G+≅OG∗+. ^*∈ M(G^+), c(G^+)≤ c G^+ _OG^*+. (26) Proof. Fix a candidate G+G^+ satisfying the two premises in Equation 26. Relevant-candidate algebraic faithfulness gives V(G†)⊆V(G+).V(G ) V(G^+). (27) Suppose that the inclusion is strict. Proposition 2 shows that both real algebraic sets are irreducible. Applying Lemma 16 with V1=V(G†)V_1=V(G ) and V2=V(G+)V_2=V(G^+) gives dimV(G†)<dimV(G+). V(G )< V(G^+). Using the dimension identities in Proposition 2, this is equivalent to dimℳ(G†)<dimℳ(G+). (G )< (G^+). Structural minimality would therefore imply c(G+)>c†c(G^+)>c , contradicting the premise c(G+)≤c†c(G^+)≤ c . Thus V(G+)=V(G†)=:V.V(G^+)=V(G )=:V. Equality of the Zariski closures alone is insufficient because the two semialgebraic covariance families may be disjoint subsets of V. Full-dimensional overlap supplies dim(ℳ(G+)∩ℳ(G†))=dimV. (M(G^+) (G ) )= V. (28) Proposition 2 also gives dimℳ(G+)=dimV=dimℳ(G†). (G^+)= V= (G ). Consequently, dim(ℳ(G+)∩ℳ(G†))=dimℳ(G+)=dimℳ(G†), split (M(G^+) (G ) )&= (G^+)\\ &= (G ), split which is exactly G+≅OG†G^+ _OG . Finally, ℳ(G†)=ℳ(G∗+)M(G )=M(G^*+) by construction. Replacing ℳ(G†)M(G ) by this equal set in the preceding dimension identity proves G+≅OG∗+G^+ _OG^*+. The argument applies to every member of 0(Σ∗)G_0( ^*), so Assumption 5 holds. ∎ A.8 KL separation and convergence of the minimum KL divergence In the finite-dimensional space pS^p, a set is compact if and only if it is closed and bounded in that space. A continuous real-valued function on a nonempty compact set attains its minimum. We first establish a compactness result used in both lemmas below. It applies to a fixed first covariance and, more generally, to first covariances varying over a compact set. Lemma 18 (Compactness for bounded Gaussian KL divergence). Let K⊆++pK _++^p be nonempty and compact, and let 0≤t<∞0≤ t<∞. Then the set Kt:=Σ∈++p:infS∈KD(S∥Σ)≤tK_t:= \ _++^p: _S∈ KD(S\| )≤ t \ is compact in pS^p and consists of positive-definite matrices. In particular, taking K=Σ∗K=\ ^*\ shows that Σ∈++p:D(Σ∗∥Σ)≤t\ _++^p:D( ^*\| )≤ t\ is compact. Proof. We first bound the eigenvalues of every Σ∈Kt ∈ K_t away from zero and above by a finite constant. Because K is compact and its matrices are positive definite, continuity of their smallest and largest eigenvalues gives constants 0<m≤M<∞0<m≤ M<∞ such that mIp⪯S⪯MIp(S∈K).mI_p S MI_p (S∈ K). Write u1,…,up>0u_1,…,u_p>0 for the eigenvalues of Σ . Expanding Equation 12, and using tr(Σ−1S)≥mtr(Σ−1)tr( ^-1S)≥ mtr( ^-1) and logdetS≤plogM S≤ p M, gives D(S∥Σ)≥12∑i=1p(mui+logui)−plogM−p(S∈K).D(S\| )≥ 12 \ _i=1^p ( mu_i+ u_i )-p M-p \ (S∈ K). For u>0u>0, let h(u)=m/u+loguh(u)=m/u+ u. Its derivative is h′(u)=(u−m)/u2h (u)=(u-m)/u^2, so its minimum is h(m)=1+logmh(m)=1+ m. Also, h(u)→+∞h(u)→+∞ as u↓0u 0 or u→+∞u→+∞. The preceding lower bound is independent of S. Thus, if Σ∈Kt ∈ K_t, ∑i=1ph(ui)≤2t+plogM+p. _i=1^ph(u_i)≤ 2t+p M+p. Each of the other p−1p-1 terms is at least h(m)h(m), so each eigenvalue satisfies h(ui)≤2t+plogM+p−(p−1)h(m).h(u_i)≤ 2t+p M+p-(p-1)h(m). Since the right-hand side is finite and h tends to infinity at both ends of (0,∞)(0,∞), there are constants 0<a<b<∞0<a<b<∞, depending only on K and t, such that a≤ui≤ba≤ u_i≤ b for every i. Consequently, Kt⊆B:=Σ∈p:aIp⪯Σ⪯bIp.K_t B:=\ ^p:aI_p bI_p\. The set B is closed and bounded, hence compact in the finite-dimensional space pS^p, and a>0a>0 ensures that B⊆++pB _++^p. It remains to show that KtK_t is closed in B. Define f(Σ)=minS∈KD(S∥Σ)f( )= _S∈ KD(S\| ); this minimum exists by continuity and compactness of K. For Σ,Σ~∈B , ∈ B, |f(Σ)−f(Σ~)|≤supS∈K|D(S∥Σ)−D(S∥Σ~)|.|f( )-f( )|≤ _S∈ K|D(S\| )-D(S\| )|. Joint continuity of D on the compact product K×BK× B implies uniform continuity there, so the right-hand side tends to zero as Σ~→Σ → . Hence f is continuous on B. It follows that Kt=Σ∈B:f(Σ)≤tK_t=\ ∈ B:f( )≤ t\ is a closed subset of the compact set B, and is therefore compact. ∎ Lemma 19 (Positive KL gap). Let Σ∗∈++p ^* _++^p. If C⊆++pC _++^p is closed relative to ++pS_++^p and Σ∗∉C ^*∉ C, then infΣ∈CD(Σ∗∥Σ)>0 _ ∈ CD( ^*\| )>0. Proof. If C is empty, the infimum is +∞+∞ and the result is immediate. Assume henceforth that C is nonempty. Although every Σ∈C ∈ C differs from Σ∗ ^* and therefore has positive KL divergence, this alone does not rule out an infimum of zero. We first prove that the infimum is attained; closedness of C alone does not guarantee this. Restrict the minimization without changing its infimum. Choose any Σ0∈C _0∈ C and set t=D(Σ∗∥Σ0)<∞t=D( ^*\| _0)<∞. This candidate gives infΣ∈CD(Σ∗∥Σ)≤t _ ∈ CD( ^*\| )≤ t. Define Kt=Σ∈++p:D(Σ∗∥Σ)≤t.K_t=\ _++^p:D( ^*\| )≤ t\. The intersection C∩KtC∩ K_t contains Σ0 _0, whose divergence is t. Every covariance in C∖KtC K_t has divergence greater than t, so excluding these covariances does not change the infimum: infΣ∈CD(Σ∗∥Σ)=infΣ∈C∩KtD(Σ∗∥Σ). _ ∈ CD( ^*\| )= _ ∈ C∩ K_tD( ^*\| ). Lemma 18, with K=Σ∗K=\ ^*\, proves that KtK_t is compact. Since C is closed relative to ++pS_++^p and Kt⊆++pK_t _++^p, the intersection C∩KtC∩ K_t is closed in KtK_t. It is therefore a nonempty compact set. The continuous function Σ↦D(Σ∗∥Σ) D( ^*\| ) attains its minimum on this set at some Σ^∈C∩Kt ∈ C∩ K_t. Combining this fact with Step 1 gives infΣ∈CD(Σ∗∥Σ)=minΣ∈C∩KtD(Σ∗∥Σ)=D(Σ∗∥Σ^). _ ∈ CD( ^*\| )= _ ∈ C∩ K_tD( ^*\| )=D( ^*\| ). Suppose, for a contradiction, that this minimum equals zero. Gaussian KL divergence is zero only when its two covariances are equal, so D(Σ∗∥Σ^)=0D( ^*\| )=0 implies Σ^=Σ∗ = ^*. But Σ^∈C ∈ C, whereas Σ∗∉C ^*∉ C by assumption. This is a contradiction. The attained minimum is nonnegative and cannot be zero, so it is strictly positive, as required. ∎ Lemma 20 (Convergence of the minimum empirical KL divergence). For fixed G+G^+, let Δ(G+)=infΣ∈ℳ¯(G+)D(Σ∗∥Σ). (G^+)= _ ∈ M(G^+)D( ^*\| ). (29) Then δn(G+)→Δ(G+) _n(G^+)→ (G^+) almost surely. If Σ∗∉ℳ¯(G+) ^*∉ M(G^+), then Δ(G+)>0 (G^+)>0. Here, almost-sure convergence means that the convergence event has probability one: Pr(limn→∞δn(G+)=Δ(G+))=1. \! ( _n→∞ _n(G^+)= (G^+) )=1. The probability is taken over the infinite sequence of i.i.d. observations. Proof. Because the observations have finite second moments, the strong law applied entrywise to XO(k)XO(k)⊤X_O^(k)X_O^(k) gives Sn⟶Σ∗almost surely.S_n ^* surely. Fix an outcome in this probability-one event for the remainder of the proof. Choose ε>0 >0 smaller than the smallest eigenvalue of Σ∗ ^* and let K=S∈p:‖S−Σ∗‖F≤ε.K=\S ^p:\|S- ^*\|_F≤ \. This closed ball is compact (closed and bounded in the finite-dimensional space pS^p), lies in ++pS_++^p, and contains SnS_n for all sufficiently large n. The model closure ℳ¯(G+) M(G^+) is nonempty because the admissible parameter space is nonempty. Fix any Σ0∈ℳ¯(G+) _0∈ M(G^+). Joint continuity of D(S∥Σ)D(S\| ) on ++p×++pS_++^p×S_++^p and compactness of K imply R0:=supS∈KD(S∥Σ0)<∞.R_0:= _S∈ KD(S\| _0)<∞. Thus infΣ∈ℳ¯(G+)D(S∥Σ)≤R0 _ ∈ M(G^+)D(S\| )≤ R_0 uniformly over S∈KS∈ K. We next show that the infimum over Σ can be restricted to a single compact set for every S∈KS∈ K. By Lemma 18, applied with this K and t=R0+1t=R_0+1, the set K0=Σ∈++p:infS∈KD(S∥Σ)≤R0+1K_0=\ _++^p: _S∈ KD(S\| )≤ R_0+1\ is compact and contained in ++pS_++^p. For every S∈KS∈ K, the infimum over Σ can now be restricted to ℳ¯(G+)∩K0 M(G^+)∩ K_0. Indeed, if Σ∉K0 ∉ K_0, then D(S∥Σ)>R0+1D(S\| )>R_0+1, while the fixed feasible covariance Σ0 _0 has D(S∥Σ0)≤R0D(S\| _0)≤ R_0. Thus points outside K0K_0 cannot improve the infimum. The restricted feasible set ℳ¯(G+)∩K0 M(G^+)∩ K_0 is nonempty (it contains Σ0 _0) and compact. Thus the minimum KL divergence over the covariance-family closure is attained, although admissible SCM parameters need not attain the same value. The function D is uniformly continuous on the compact product K×K0K× K_0. Since Sn→Σ∗S_n→ ^* and both matrices eventually belong to K, it follows that supΣ∈K0|D(Sn∥Σ)−D(Σ∗∥Σ)|→0 _ ∈ K_0|D(S_n\| )-D( ^*\| )|→ 0 on the fixed probability-one event. Finally, for any two real-valued functions fnf_n and f with finite infima on a common nonempty set C, |infCfn−infCf|≤supC|fn−f|. | _Cf_n- _Cf |≤ _C|f_n-f|. Applying this inequality with C=ℳ¯(G+)∩K0C= M(G^+)∩ K_0 proves δn(G+)⟶Δ(G+) _n(G^+) (G^+) almost surely. If Σ∗∉ℳ¯(G+) ^*∉ M(G^+), the set ℳ¯(G+) M(G^+) is closed relative to ++pS_++^p by definition. Lemma 19, applied with C=ℳ¯(G+)C= M(G^+), then gives Δ(G+)>0 (G^+)>0. ∎ A.9 Proof of Theorem 4 We use the following probability notation for real random variables YnY_n and positive deterministic numbers ana_n. The notation Yn=Op(an)Y_n=O_p(a_n) means that, for every ε>0 >0, there is a finite M>0M>0 such that Pr(|Yn|>Man)<ε (|Y_n|>Ma_n)< for all sufficiently large n. The notation Yn=op(an)Y_n=o_p(a_n) means that Pr(|Yn|>εan)→0 (|Y_n|> a_n)→ 0 for every ε>0 >0. We write Yn→yY_n py for convergence in probability to y, meaning Yn−y=op(1)Y_n-y=o_p(1), and Yn→+∞Y_n p+∞ if Pr(Yn>T)→1 (Y_n>T)→ 1 for every fixed real T. For deterministic sequences, an=o(n)a_n=o(n) means an/n→0a_n/n→ 0. Proof. By definition, G†G is distribution-equivalent to the true graph G∗+G^*+. Hence Σ∗∈ℳ(G∗+)=ℳ(G†)⊆ℳ¯(G†). ^* (G^*+)=M(G ) M(G ). Recall that δn(G†) _n(G ) is the infimum of D(Sn∥Σ)D(S_n\| ) over Σ∈ℳ¯(G†) ∈ M(G ). Every divergence in this infimum is nonnegative, giving the lower bound below. Since Σ∗ ^* belongs to the set over which the infimum is taken, the infimum is no greater than the divergence evaluated at Σ∗ ^*. Thus 0≤δn(G†)=infΣ∈ℳ¯(G†)D(Sn∥Σ)≤D(Sn∥Σ∗).0≤ _n(G )= _ ∈ M(G )D(S_n\| )≤ D(S_n\| ^*). This upper bound does not require the infimum to be attained or assert that Σ∗ ^* minimizes the empirical divergence. It uses only that Σ∗ ^* is one of the covariances included in the infimum. We next bound D(Sn∥Σ∗)D(S_n\| ^*) to show that nδn(G†)n _n(G ) is bounded in probability. Taylor’s theorem first reduces this task to bounding the squared covariance error ‖Sn−Σ∗‖F2\|S_n- ^*\|_F^2. Define f(S)=D(S∥Σ∗)f(S)=D(S\| ^*) for S∈++pS _++^p. The Gaussian KL formula in Equation 12 shows that f is twice continuously differentiable on ++pS_++^p. Since f(Σ∗)=0f( ^*)=0 is its minimum and Σ∗ ^* is an interior point of ++pS_++^p in pS^p, the first derivative there is zero. Choose r>0r>0 so that the closed Frobenius ball S∈p:‖S−Σ∗‖F≤r\S ^p:\|S- ^*\|_F≤ r\ lies in ++pS_++^p. The continuous second derivatives are bounded on this compact ball. The standard quadratic Taylor bound (Boyd and Vandenberghe, 2004, Section 9.1.2, Equation (9.13), p. 461), applied at Σ∗ ^*, therefore gives a finite constant C>0C>0 such that 0≤D(S∥Σ∗)≤C∥S−Σ∗∥F2whenever ∥S−Σ∗∥F≤r.0≤ D(S\| ^*)≤ C\|S- ^*\|_F^2 \|S- ^*\|_F≤ r. Both r and C are fixed independently of n. To obtain the covariance-error rate, we use Markov’s inequality: for any nonnegative random variable Z with finite expectation and any a>0a>0, Pr(Z>a)≤Za. (Z>a)≤ EZa. We will apply this inequality to Z=‖Sn−Σ∗‖F2Z=\|S_n- ^*\|_F^2, with threshold a=M2/na=M^2/n for M>0M>0. We therefore first compute ‖Sn−Σ∗‖F2E\|S_n- ^*\|_F^2. For this calculation, fix i,ji,j and set Yij,k=XO,i(k)XO,j(k)Y_ij,k=X_O,i^(k)X_O,j^(k). Then (Sn)ij=n−1∑k=1nYij,k(S_n)_ij=n^-1 _k=1^nY_ij,k. The observations are independent and identically distributed across k, so these products are also independent and identically distributed across k; independence of the two coordinates within an observation is not required. The known zero mean gives Yij,k=Σij∗EY_ij,k= ^*_ij. Their variance is finite because Gaussian coordinates have finite fourth moments and the Cauchy–Schwarz inequality gives Yij,12≤[(XO,i(1))4][(XO,j(1))4]1/2<∞.EY_ij,1^2≤ \E [(X_O,i^(1))^4 ]E [(X_O,j^(1))^4 ] \^1/2<∞. Consequently, each sample-covariance entry is unbiased. Independence across observations implies that the variance of the average is the sum of the individual variances divided by n2n^2: [((Sn)ij−Σij∗)2] [((S_n)_ij- ^*_ij)^2 ] =Var((Sn)ij) =Var((S_n)_ij) =1n2∑k=1nVar(Yij,k)=Var(Yij,1)n. = 1n^2 _k=1^nVar(Y_ij,k)= Var(Y_ij,1)n. The squared Frobenius norm is the sum of these squared entry errors. Summing their expectations therefore gives ‖Sn−Σ∗‖F2=Kn,K:=∑i,j=1pVar(Yij,1)<∞.E\|S_n- ^*\|_F^2= Kn, K:= _i,j=1^pVar(Y_ij,1)<∞. Since p and the true distribution are fixed, K does not depend on n. Substituting this expectation into Markov’s inequality gives, for any M>0M>0, Pr(‖Sn−Σ∗‖F>Mn) \! (\|S_n- ^*\|_F> M n ) =Pr(‖Sn−Σ∗‖F2>M2n) = \! (\|S_n- ^*\|_F^2> M^2n ) ≤‖Sn−Σ∗‖F2M2/n=K/nM2/n=KM2. ≤ E\|S_n- ^*\|_F^2M^2/n= K/nM^2/n= KM^2. For any ε>0 >0, choosing M>K/εM> K/ makes this probability less than ε for every n. By the definition of OpO_p stated above, this proves ∥Sn−Σ∗∥F=Op(n−1/2).\|S_n- ^*\|_F=O_p(n^-1/2). Finally, combine the local Taylor bound with the expectation bound. Outside the ball where the Taylor bound applies, Markov’s inequality gives Pr(‖Sn−Σ∗‖F>r)≤K/(nr2) (\|S_n- ^*\|_F>r)≤ K/(nr^2). Inside that ball, nD(Sn∥Σ∗)nD(S_n\| ^*) is at most nC‖Sn−Σ∗‖F2nC\|S_n- ^*\|_F^2. Hence, for every T>0T>0, Pr(nD(Sn∥Σ∗)>T)≤Knr2+CKT. \! (nD(S_n\| ^*)>T )≤ Knr^2+ CKT. For any prescribed probability bound, the second term can be made arbitrarily small by choosing T large, and the first by taking n large. Thus nD(Sn∥Σ∗)=Op(1)nD(S_n\| ^*)=O_p(1), or equivalently D(Sn∥Σ∗)=Op(n−1)D(S_n\| ^*)=O_p(n^-1). The preceding inequality 0≤nδn(G†)≤nD(Sn∥Σ∗)0≤ n _n(G )≤ nD(S_n\| ^*) now gives nδn(G†)=Op(1).n _n(G )=O_p(1). (30) Complexity of every global minimizer. Let G^n+∈argminG+∈p,Lmaxsn(G+) G_n^+∈ *arg\,min_G^+ _p,L_ s_n(G^+) be any global minimizer. Score optimality and cancellation of the common constant CnC_n give nδn(G^n+)+λnc(G^n+)≤nδn(G†)+λnc†.n _n( G_n^+)+ _nc( G_n^+)≤ n _n(G )+ _nc . Dropping the nonnegative first term on the left and dividing by λn>0 _n>0 yields c(G^n+)≤c†+nδn(G†)λn.c( G_n^+)≤ c + n _n(G ) _n. Equation 30 states that the numerator nδn(G†)n _n(G ) is bounded in probability. This does not imply that the numerator converges to zero. However, dividing it by the deterministic sequence λn→∞ _n→∞ makes the ratio converge to zero in probability. To verify this, fix ε>0 >0 and η>0η>0. Boundedness in probability provides a finite M>0M>0 such that Pr(nδn(G†)>M)<η (n _n(G )>M)<η for all sufficiently large n. For such n large enough that ελn>M _n>M, we obtain Pr(nδn(G†)λn>ε)=Pr(nδn(G†)>ελn)≤Pr(nδn(G†)>M)<η. \! ( n _n(G ) _n> )= \! (n _n(G )> _n )≤ \! (n _n(G )>M )<η. Since η can be arbitrarily small, the ratio converges to zero in probability, which is precisely nδn(G†)λn=op(1). n _n(G ) _n=o_p(1). Thus the complexity bound above can be written as c(G^n+)≤c†+op(1)c( G_n^+)≤ c +o_p(1). Because complexity is integer-valued, c(G^n+)>c†c( G_n^+)>c implies c(G^n+)−c†≥1c( G_n^+)-c ≥ 1. The bound therefore gives Pr(c(G^n+)>c†)≤Pr(nδn(G†)λn≥1)⟶0. \! (c( G_n^+)>c )≤ \! ( n _n(G ) _n≥ 1 ) 0. This conclusion holds simultaneously for every global minimizer because the upper bound does not depend on which minimizer is chosen. Excluding families whose closure does not contain the true covariance. Fix a candidate G+∈p,LmaxG^+ _p,L_ such that Σ∗∉ℳ¯(G+) ^*∉ M(G^+). Lemma 19 gives Δ(G+)=infΣ∈ℳ¯(G+)D(Σ∗∥Σ)>0. (G^+)= _ ∈ M(G^+)D( ^*\| )>0. Lemma 20 gives δn(G+)→Δ(G+) _n(G^+)→ (G^+) almost surely. Almost-sure convergence implies convergence in probability, so, for every ε>0 >0, Pr(|δn(G+)−Δ(G+)|>ε)⟶0. \! (| _n(G^+)- (G^+)|> ) 0. By the definition of op(1)o_p(1) stated above, this is equivalent to δn(G+)=Δ(G+)+op(1). _n(G^+)= (G^+)+o_p(1). Equation 13 writes each graph score as the sum of the common data-dependent term CnC_n, the empirical KL-divergence infimum multiplied by n, and the structural penalty. Applying this formula to both graphs, with c(G†)=c†c(G )=c , gives sn(G+)−sn(G†) s_n(G^+)-s_n(G ) =[Cn+nδn(G+)+λnc(G+)] = [C_n+n _n(G^+)+ _nc(G^+) ] −[Cn+nδn(G†)+λnc†] - [C_n+n _n(G )+ _nc ] =nδn(G+)−nδn(G†)+λnc(G+)−c†. =n _n(G^+)-n _n(G )+ _n\c(G^+)-c \. The two occurrences of CnC_n cancel because both scores use the same data; this term does not depend on the candidate graph. The convergence above gives δn(G+)−Δ(G+)=op(1) _n(G^+)- (G^+)=o_p(1). Multiplication by n yields nδn(G+)=nΔ(G+)+nδn(G+)−Δ(G+)=nΔ(G+)+op(n).n _n(G^+)=n (G^+)+n\ _n(G^+)- (G^+)\=n (G^+)+o_p(n). Indeed, the remainder divided by n is δn(G+)−Δ(G+) _n(G^+)- (G^+), which converges to zero in probability. For the reference graph, Equation 30 gives nδn(G†)=Op(1)n _n(G )=O_p(1). Substituting these two bounds into the exact score difference therefore gives sn(G+)−sn(G†) s_n(G^+)-s_n(G ) =nΔ(G+)+op(n)−Op(1) =n (G^+)+o_p(n)-O_p(1) +λnc(G+)−c†. + _n\c(G^+)-c \. Here the subtracted Op(1)O_p(1) term is the nonnegative reference quantity nδn(G†)n _n(G ), which is bounded in probability. Both graphs are fixed, so their complexity difference is constant. Since λn/n→0 _n/n→ 0, the penalty difference divided by n tends to zero, even if the candidate has smaller complexity than the reference. Also, an Op(1)O_p(1) term divided by n converges to zero in probability. Consequently, sn(G+)−sn(G†)n=Δ(G+)+op(1)→Δ(G+)>0. s_n(G^+)-s_n(G )n= (G^+)+o_p(1) p (G^+)>0. In particular, the score difference is at least nΔ(G+)/2n (G^+)/2 with probability tending to one. This lower bound tends to infinity, so sn(G+)−sn(G†)→+∞s_n(G^+)-s_n(G ) p+∞. Thus, for this candidate, the probability of being a global score minimizer tends to zero. The class p,LmaxG_p,L_ is finite. A finite union bound therefore makes the preceding exclusion simultaneous over all candidates whose covariance-family closures do not contain Σ∗ ^*. Combining that event with the simultaneous complexity bound shows that, with probability tending to one, every global minimizer satisfies Σ∗∈ℳ¯(G^n+),c(G^n+)≤c†. ^*∈ M( G_n^+), c( G_n^+)≤ c . Equivalently, argminG+∈p,Lmaxsn(G+)⊆0(Σ∗), *arg\,min_G^+ _p,L_ s_n(G^+) _0( ^*), which proves Equation 16. The proof used only λn→∞ _n→∞ for the complexity comparison and λn=o(n) _n=o(n) for exclusion of misspecified covariance models. Hence it also proves the stated extension to general deterministic penalty sequences. ∎ A.10 Proof of Corollary 6 Proof. Let n:=argminG+∈p,Lmaxsn(G+)A_n:= *arg\,min_G^+ _p,L_ s_n(G^+) be the set of all globally minimizing graphs. This set is nonempty because the candidate class is finite and every graph has a finite score, as shown in Appendix A.2. Define the event En:=n⊆0(Σ∗).E_n:=\A_n _0( ^*)\. Theorem 4 gives Pr(En)→1 (E_n)→ 1. Fix a data realization in EnE_n and take any G^n+∈n G_n^+ _n. By the definition of 0(Σ∗)G_0( ^*), Σ∗∈ℳ¯(G^n+),c(G^n+)≤c†. ^*∈ M( G_n^+), c( G_n^+)≤ c . Assumption 5 therefore gives G^n+≅OG† G_n^+ _OG . By the definition of marginal quasi-equivalence, this means dim(ℳ(G^n+)∩ℳ(G†))=dimℳ(G^n+)=dimℳ(G†). (M( G_n^+) (G ) )= ( G_n^+)= (G ). The reference graph G†G is distribution-equivalent to the true graph, so ℳ(G†)=ℳ(G∗+)M(G )=M(G^*+). Substituting this equality into the preceding identity gives dim(ℳ(G^n+)∩ℳ(G∗+))=dimℳ(G^n+)=dimℳ(G∗+). (M( G_n^+) (G^*+) )= ( G_n^+)= (G^*+). Thus G^n+≅OG∗+ G_n^+ _OG^*+. This step uses equality of covariance families, not transitivity of quasi-equivalence. Because the minimizer was arbitrary, the conclusion holds simultaneously for every graph in nA_n, including when the score has ties. Hence En⊆n⊆G+∈p,Lmax:G+≅OG∗+.E_n \A_n \G^+ _p,L_ :G^+ _OG^*+\ \. It follows that 1≥Prn⊆G+∈p,Lmax:G+≅OG∗+≥Pr(En)⟶1,1≥ \! \A_n \G^+ _p,L_ :G^+ _OG^*+\ \≥ (E_n) 1, which proves Equation 17. ∎ A.11 Details and proof of the Bernoulli formulation This subsection supplies the parameterization, likelihood, and differentiability details used in Section 4.2, followed by the proof of Proposition 10. Effective latent loadings. For a graph with ℓ latent variables, place ΩL1/2Λ _L^1/2 in ℓ rows of Γ and fill the remaining Lmax−ℓL_ - rows with zeros. The corresponding gates select the graph’s observed edges, active latent slots, and latent-to-observed edges. Absorbing latent noise standard deviations preserves the covariance family: conversely, any effective loadings can be represented using unit latent variances. All configurations share the same continuous parameters ϑ=(W,Γ,ΩO) =(W, , _O). For every binary (H,Z)(H,Z), the latent contribution in Equation 20 is positive semidefinite, so Ψ(H,Z)≻0 (H,Z) 0 because ΩO≻0 _O 0. Sampled likelihood and complexity. For nonsingular A(M)A(M), evaluate Equation 9 at the masked covariance and remove the Gaussian constant nplog(2π)/2np (2π)/2. Using the identities derived in Appendix A.2, this gives Ln(M,H,Z,ϑ)=−nlog|detA(M)|+n2logdetΨ(H,Z)+n2tr[Ψ(H,Z)−1A(M)⊤SnA(M)]. splitL_n(M,H,Z; )=&-n | A(M)|+ n2 (H,Z)\\ &+ n2tr\! [ (H,Z)^-1A(M) S_nA(M) ]. split (31) As in the main text, Ln=+∞L_n=+∞ when A(M)A(M) is singular. The sampled structural complexity is c(M,H,Z)=∥M∥1+∑h=1LmaxZh(1+∥Hh⋅∥1),c(M,H,Z)=\|M\|_1+ _h=1^L_ Z_h (1+\|H_h·\|_1 ), (32) where Hh⋅H_h· denotes row h of H. The first term counts selected edges among observed variables. Each active slot contributes one latent variable and one unit for each of its selected edges to observed variables. An inactive slot contributes neither covariance nor complexity, regardless of its edge gates. The active latent count is ∑hZh _hZ_h, with expectation ∑hπh _h _h. Independence of the gates gives [c(M,H,Z)]=∑i≠jqij+∑h=1Lmaxπh(1+∑j=1prhj).E[c(M,H,Z)]= _i≠ jq_ij+ _h=1^L_ _h (1+ _j=1^pr_hj ). Thus Equation 21 is exactly [Ln+λnc]E[L_n+ _nc], not an approximation obtained by replacing the binary masks in the likelihood with their probabilities. Differentiability and singular configurations. Fix ϑ and suppose every binary configuration has finite likelihood. Its probability is a product of gate probabilities or their complements, with at most one factor involving each scalar probability. The expectation is a finite sum of these products multiplied by the configuration scores. It is therefore a polynomial that is affine (linear plus a constant) in each scalar probability when all others are fixed. In particular, it is differentiable in the probabilities, but need not be convex in all of them jointly. The same finite-sum representation also gives joint differentiability in the probabilities and continuous parameters. At a parameter vector ϑ for which every binary configuration has finite likelihood, the noise variances remain positive and all required determinants remain nonzero on a sufficiently small open neighborhood of ϑ . Each configuration’s likelihood in Equation 31 is continuously differentiable there because the matrix inverses and logarithms in that expression are continuously differentiable on the stated domains. Adding each configuration’s constant complexity penalty preserves this property. The finite sum of configuration scores weighted by polynomial configuration probabilities is therefore jointly continuously differentiable in the gate probabilities and ϑ on this domain. If some configurations have singular A(M)A(M), any positive probability assigned to one of them makes the objective infinite. A configuration of zero probability is excluded from the expectation; no product of zero and infinity is evaluated. Differentiability at a finite objective value does not automatically extend to probability changes that activate a singular configuration. Instead, if some probabilities are fixed at zero or one so that every singular configuration is excluded throughout the remaining probability domain, the same polynomial argument applies to the free probabilities on that domain. Probabilities equal to zero or one are permitted; when all are binary, the distribution selects a single configuration deterministically. Simultaneous invertibility. When every observed-edge probability lies strictly between zero and one, all binary observed-edge masks have positive probability. Finite expected likelihood then requires det(Ip−M⊙W)≠0for every binary observed-edge mask M. (I_p-M W)≠ 0 every binary observed-edge mask M. For each fixed mask, this determinant is a polynomial in the off-diagonal entries of W, and it equals one at W=0W=0. It is therefore not the zero polynomial. The excluded weights form a finite union of zero sets of nonzero polynomials. Each such zero set has Lebesgue measure zero, as explained in Appendix B.1, so their union also has measure zero. The permitted weights form a Euclidean open set because there are only finitely many determinant conditions and each determinant is continuous. Thus this requirement does not force a norm bound or an acyclic graph, although approaching a singular configuration can cause numerical difficulties. Proof of Proposition 10. Proof. Let K denote the finite set of all binary triples k=(M,H,Z)k=(M,H,Z). The optimization has a finite feasible value: take W=Γ=0W= =0 and ΩO=Ip _O=I_p. Moreover, every finite likelihood is bounded below by KnK_n from Equation 25, and the complexity penalty is nonnegative. Hence the infima below are finite. For fixed gate probabilities (Q,R,π)(Q,R,π), write wk(Q,R,π)w_k(Q,R,π) for the product Bernoulli probability of configuration k. Then wk≥0w_k≥ 0, ∑k∈wk=1 _k w_k=1, and, for every ϑ for which the expected objective is finite, [Fnbin(M,H,Z,ϑ)]=∑k∈wk(Q,R,π)Fnbin(k,ϑ).E[ [rgb]0,0,0F_n^bin(M,H,Z; )]= _k w_k(Q,R,π) [rgb]0,0,0F_n^bin(k; ). Terms with zero probability make no contribution. The remaining terms form a convex combination: a weighted average with nonnegative weights summing to one. At least one score in this average is no larger than the average itself. Hence there is a configuration k0=(M0,H0,Z0)k_0=(M_0,H_0,Z_0) in the support, meaning wk0>0w_k_0>0, such that Fnbin(k0,ϑ)≤[Fnbin(M,H,Z,ϑ)]. [rgb]0,0,0F_n^bin(k_0; ) [ [rgb]0,0,0F_n^bin(M,H,Z; )]. (33) For any such tuple (Q,R,π,ϑ)(Q,R,π, ), [Fnbin(M,H,Z,ϑ)] [ [rgb]0,0,0F_n^bin(M,H,Z; )] ≥Fnbin(k0,ϑ) ≥ [rgb]0,0,0F_n^bin(k_0; ) ≥infϑ′Fnbin(k0,ϑ′) ≥ _ [rgb]0,0,0F_n^bin(k_0; ) ≥mininfϑ′k∈Fnbin(k,ϑ′). ≥ _k _ [rgb]0,0,0F_n^bin(k; ). Taking the infimum over (Q,R,π,ϑ)(Q,R,π, ) on the left proves infQ,R,π,ϑ[Fnbin]≥mininfϑk∈Fnbin(k,ϑ). _Q,R,π, E[ [rgb]0,0,0F_n^bin]≥ _k _ [rgb]0,0,0F_n^bin(k; ). For the reverse inequality, fix an arbitrary binary configuration k=(M,H,Z)k=(M,H,Z). Choosing the entries of Q, R, and π equal to the corresponding zeros and ones in k makes the product Bernoulli distribution a point mass at k, meaning wk=1w_k=1 and all other configuration probabilities are zero. The expectation then equals Fnbin(k,ϑ) [rgb]0,0,0F_n^bin(k; ) for every ϑ . Taking the infimum over ϑ and then choosing the best binary configuration gives infQ,R,π,ϑ[Fnbin]≤mininfϑk∈Fnbin(k,ϑ). _Q,R,π, E[ [rgb]0,0,0F_n^bin]≤ _k _ [rgb]0,0,0F_n^bin(k; ). The two inequalities prove Equation 22. If the infimum of the expected objective is attained at a possibly fractional tuple, apply Equation 33 to that minimizer. Binary gate probabilities concentrated on k0k_0 give a feasible point whose value is no larger than the attained global minimum. Its value must therefore be equal to the minimum, proving the existence of a binary global minimizer. It remains to show that optimizing over all binary masks has the same optimum as optimizing over the candidate class p,LmaxG_p,L_ , which requires every active latent to have at least two observed children. First, a slot with Zh=0Z_h=0 has no contribution to either Ψ or the complexity in Equation 32; it can simply be omitted. Next, suppose Zh=1Z_h=1 but ∥Hh⋅∥1=0\|H_h·\|_1=0. The hhth latent contribution to Ψ is zero. Setting Zh=0Z_h=0 therefore preserves the likelihood and strictly lowers the structural complexity. Finally, suppose an active slot has exactly one child, indexed by j. Let eje_j be the jjth standard basis vector in ℝpR^p. Its contribution to the effective noise covariance is (Hh⋅⊙Γh⋅)⊤(Hh⋅⊙Γh⋅)=Γhj2ejej⊤.(H_h· _h·) (H_h· _h·)= _hj^2e_je_j . Define a new diagonal observed-noise covariance by ΩO′=ΩO+Γhj2ejej⊤ _O= _O+ _hj^2e_je_j and deactivate slot h. The matrix ΩO′ _O remains diagonal and positive definite, and the total effective covariance Ψ is unchanged. The observed coefficient matrix A(M)A(M) is also unchanged, so every term in the likelihood is identical, while the structural complexity decreases by the latent node and its single outgoing edge. Repeating this operation removes all active slots with fewer than two children without increasing the score. This argument applies separately to a fixed binary configuration. Because the adjustment of ΩO _O depends on that configuration, it does not define a simultaneous transformation of a Bernoulli distribution whose configurations share the same continuous parameters. After inactive slots are deleted and the remaining slots are relabeled, the resulting mask represents a graph in p,LmaxG_p,L_ . Conversely, every graph in p,LmaxG_p,L_ can be represented by some binary (M,H,Z)(M,H,Z). For a fixed mask representing G+G^+, the infimum of the binary-mask objective over the continuous parameters is therefore Kn+nδn(G+)+λnc(G+)K_n+n _n(G^+)+ _nc(G^+), by Equation 25. Since Cn=Kn+nplog(2π)/2C_n=K_n+np (2π)/2, this value equals sn(G+)−nplog(2π)/2s_n(G^+)-np (2π)/2. The conversion from arbitrary masks to admissible graphs cannot increase the objective, and every admissible graph has a binary representation. Therefore the common infimum is minG+∈p,Lmaxsn(G+)−nplog(2π)/2 _G^+ _p,L_ s_n(G^+)-np (2π)/2, and the two formulations have the same minimizing graphs after the stated conversion. ∎ Scope of the exactness result. The proposition does not assert that all globally minimizing probabilities are binary. At a finite global minimum, a fractional probability vector can also be optimal if every configuration assigned positive probability attains that same minimum score at the shared continuous parameters. The proof establishes the existence of a deterministic minimizer only when the infimum is attained. Direct evaluation of the expectation sums over exponentially many configurations. The closed-form complexity expectation does not remove this cost for the likelihood expectation. Equality of global infima neither bounds the variance of sampled gradient estimates nor guarantees that a numerical optimizer reaches either infimum. Appendix B Scope and genericity of the identification assumptions We examine the three sufficient identification assumptions in Section 4.1. Subsection B.1 establishes genericity of algebraic faithfulness and full-dimensional overlap. Structural minimality concerns the target covariance family, the entire candidate graph class, and the chosen structural complexity, not the particular value of Σ∗ ^*. It is not implied by the genericity result. Subsection B.2 shows that the penalty counting edges and latent variables can violate structural minimality. Subsection B.3 presents a weaker condition at the true covariance and sufficient conditions for a restricted class of covariance families. B.1 Genericity of algebraic faithfulness and full-dimensional overlap Fix a true graph and a finite candidate class. We show that algebraic faithfulness and full-dimensional overlap hold for almost every admissible choice of the true model parameters. Here, “almost every” means that the exceptional set has Lebesgue measure zero, or zero volume, in the free parameter coordinates. This is the meaning of generic in this result; it does not establish these assumptions at every true covariance. Lemma 21 (Genericity of faithfulness and model overlap). Assume Proposition 2 holds for every graph in the finite candidate class. Let d be the number of free parameter coordinates of the true graph, and write (G∗+)⊆ℝdP(G^*+) ^d. Set V∗:=V(G∗+)=V(G†),ϕ∗:(G∗+)⟶V∗,ϕ∗(θ)=ΣO(θ).V^*:=V(G^*+)=V(G ), φ^*:P(G^*+) V^*, φ^*(θ)= _O(θ). Thus ϕ∗φ^* sends each admissible parameter vector to its observed covariance. Outside a set of Lebesgue measure zero in ℝdR^d, every admissible true parameter vector satisfies the following statements: 1. Assumption 7 holds simultaneously for all candidates; 2. Assumption 9 holds simultaneously for all candidates satisfying its premises. Moreover, for the second conclusion, the full-dimensional common subset of the two covariance models can be chosen to contain Σ∗ ^*. Proof. We first explain how an exceptional algebraic set of covariances gives a measure-zero set of true parameters. We then apply this argument to faithfulness and overlap. Exceptional covariance sets and parameter values. Let E⊊V∗E V^* be an algebraic set strictly smaller than V∗V^*. There is a polynomial f in the covariance entries that equals zero on E but not throughout V∗V^*: choose a point of V∗∖EV^* E and a defining equation of E that does not hold at that point. Write Δ(θ)=det(Ip−BOO) (θ)= (I_p-B_O). As shown in the proof of Proposition 2, each entry of ϕ∗(θ)φ^*(θ) is a ratio of polynomials in the free parameters, with denominator Δ(θ)2 (θ)^2. Substituting these entries into f and putting the terms over a common denominator gives f(ϕ∗(θ))=P(θ)Δ(θ)kf(φ^*(θ))= P(θ) (θ)^k for a polynomial P and a nonnegative integer k. The polynomial P cannot be identically zero. Otherwise, f would equal zero on every generated covariance and therefore on its Zariski closure V∗V^*, contradicting the choice of f. Since Δ(θ)≠0 (θ)≠ 0 for admissible parameters, θ∈(G∗+):ϕ∗(θ)∈E⊆θ∈ℝd:P(θ)=0.\θ (G^*+):φ^*(θ)∈ E\ \θ ^d:P(θ)=0\. The zero set of a nonzero polynomial has Lebesgue measure zero. This follows by induction on the number of variables: a nonzero polynomial in one variable has finitely many roots. In several variables, view the polynomial as a polynomial in the last variable. At least one coefficient polynomial is nonzero, so the set where all coefficient polynomials are zero has measure zero by induction. Outside that set, the polynomial has finitely many roots in the last variable; integrating over the other coordinates gives the claim. Thus the displayed set of exceptional parameter values has measure zero. This reasoning is needed because taking the inverse image of a measure-zero set does not, in general, preserve measure zero. Relevant-candidate algebraic faithfulness. Fix a candidate graph G+G^+ such that V∗⊈V(G+).V^* V(G^+). The intersection E=V∗∩V(G+)E=V^*∩ V(G^+) is an algebraic set strictly smaller than V∗V^*. The preceding argument shows that the true parameters whose covariances belong to E form a measure-zero set. Therefore, for this candidate and outside that exceptional set, Σ∗=ϕ∗(θ∗)∉V(G+). ^*=φ^*(θ^*)∉ V(G^+). The candidate class is finite, so the union of these exceptional parameter sets over all candidates with V∗⊈V(G+)V^* V(G^+) still has measure zero. Outside this union, any candidate whose Zariski closure contains Σ∗ ^* must therefore contain all of V∗=V(G†)V^*=V(G ): Σ∗∈V(G+)⟹V(G†)⊆V(G+) ^*∈ V(G^+) V(G ) V(G^+) simultaneously for every candidate. In particular, because ℳ¯(G+)⊆V(G+)∩++p M(G^+) V(G^+) _++^p (Appendix A.4), it implies Assumption 7, whose premise is restricted to candidates with c(G+)≤c†c(G^+)≤ c and Σ∗∈ℳ¯(G+) ^*∈ M(G^+). Full-dimensional model overlap. Now fix a candidate with V(G+)=V∗=:V.V(G^+)=V^*=:V. Using the interior relative to V defined in Appendix A.3, a matrix Σ∈S⊆V ∈ S V belongs to intVSint_VS exactly when some ε>0 >0 satisfies Σ′∈V:‖Σ′−Σ‖F<ε⊆S.\ ∈ V:\| - \|_F< \ S. The intersection of the ball with V, rather than the full ball in pS^p, must be contained in S. Write IG:=intVℳ(G+),I†:=intVℳ(G†)I_G:=int_VM(G^+), I_ :=int_VM(G ) for the interiors of the two covariance families relative to V. Both covariance families have dimension dimV V by Proposition 2. For a semialgebraic S⊆VS V, its boundary relative to V is clVS∖intVScl_VS _VS. A boundary point either lies outside S but is a limit of points in S, or lies in S and is a limit of points in V∖SV S. Hence clVS∖intVS⊆(clVS∖S)∪(clV(V∖S)∖(V∖S)).cl_VS _VS (cl_VS S)∪ (cl_V(V S) (V S) ) [rgb]0,0,0. (34) For any nonempty semialgebraic set T, the added limit points clpT∖Tcl_S^pT T have dimension strictly smaller than dimT T (Coste, 2000, Proposition 3.16(2)). Here V is Euclidean closed, so the closure of a subset of V within V agrees with its closure in the full symmetric matrix space. Apply this result to T=ST=S and T=V∖ST=V S; an empty set contributes no points. Since both sets have dimension at most dimV V, the displayed boundary has dimension strictly smaller than dimV V. These sets are semialgebraic by the closure and set-operation properties in Coste (2000, Section 2.1.2 and Corollary 2.5). We also need a local-coordinate property of the irreducible real algebraic set V (Akbulut and King, 1981, Section 1, p. 426–427). Write r=dimVr= V. There is an algebraic subset EV⊊VE_V V, with dimEV<r E_V<r, such that near each point of V∖EVV E_V, the set V can be described as follows: r selected covariance coordinates vary over a nonempty open subset of ℝrR^r, and the remaining coordinates are continuously differentiable functions of those selected coordinates. The choice of coordinates may depend on the point. Consequently, intersecting V with a sufficiently small open ball around such a point gives a semialgebraic set of dimension r. The exclusion of EVE_V concerns this local description of the set V, not matrix invertibility; membership in EVE_V is not defined by a zero covariance determinant. Combining this property with the boundary argument shows that each of the semialgebraic sets clVℳ(G+)∖IG,clVℳ(G†)∖I†,EVcl_VM(G^+) I_G, _VM(G ) I_ , E_V has dimension strictly smaller than dimV V. Denote their union by BGB_G. A finite union has the maximum dimension of its members. Moreover, Coste (2000, Theorem 3.20) gives dimBG¯Zar=dimBG<dimV. B_G^\,Zar= B_G< V. Since BG⊆VB_G V and V is algebraic, this Zariski closure is an algebraic subset strictly smaller than V. Apply the exceptional-parameter argument at the start of the proof with E=BG¯ZarE= B_G^\,Zar. The admissible true parameters whose covariances belong to BGB_G form a measure-zero set. The union of these sets over the finitely many candidates with V(G+)=V∗V(G^+)=V^* also has measure zero. Take a true parameter outside this union and suppose the premises of Assumption 9 hold for G+G^+. Since Σ∗∈ℳ¯(G+)⊆clVℳ(G+) ^*∈ M(G^+) _VM(G^+) and Σ∗∉BG ^*∉ B_G, we have Σ∗∈IG ^*∈ I_G. Also, Σ∗∈ℳ(G†) ^* (G ) because G†G and the true graph have the same covariance family. Exclusion from BGB_G therefore gives Σ∗∈I† ^*∈ I_ and Σ∗∉EV ^*∉ E_V. Consider the common set UG+:=IG∩I†∩(V∖EV)∩++p.U_G^+:=I_G∩ I_ ∩(V E_V) _++^p. It contains Σ∗ ^*, is semialgebraic, and is open relative to V. Indeed, IGI_G and I†I_ are interiors relative to V, V∖EVV E_V is open relative to V because EVE_V is Euclidean closed, and ++pS_++^p is open in pS^p by the strict inequalities in Sylvester’s criterion. The local-coordinate description at Σ∗ ^* therefore supplies an open neighborhood of Σ∗ ^* relative to V, contained in UG+U_G^+, whose projection contains a nonempty open subset of ℝrR^r. Our dimension definition gives dimUG+≥r U_G^+≥ r; the reverse inequality follows from UG+⊆VU_G^+ V. Thus dimUG+=dimV. U_G^+= V. Since UG+⊆ℳ(G+)∩ℳ(G†)U_G^+ (G^+) (G ), dim(ℳ(G+)∩ℳ(G†))≥dimV. (M(G^+) (G ))≥ V. On the other hand, the intersection is contained in ℳ(G+)M(G^+), whose dimension is dimV V by Proposition 2. The reverse inequality follows, and hence dim(ℳ(G+)∩ℳ(G†))=dimV. (M(G^+) (G ))= V. This is Assumption 9, and the construction shows that the full-dimensional common subset may be chosen to contain Σ∗ ^*. Finally, the union of the exceptional parameter sets for faithfulness and overlap has measure zero, so both conclusions hold outside a single such set. ∎ Lemma 21 does not prove structural minimality. Assumption 8 compares strict inclusion of Zariski closures with the discrete complexity c(G+)c(G^+) and must be verified or imposed for the chosen candidate class. Consequently, if structural minimality holds, then relevant-candidate algebraic faithfulness and full-dimensional overlap hold for generic true parameters by Lemma 21. Lemma 17 then implies that identification at the true covariance also holds generically, and Corollary 6 yields consistency up to marginal quasi-equivalence. This generic conclusion can fail at exceptional parameter values: structural minimality does not by itself establish pointwise identification at every covariance. B.2 Failure of structural minimality The following example shows why structural minimality must be checked for the chosen complexity. It does not contradict the conditional recovery result, but rules out an unconditional quasi-equivalence guarantee over the full candidate class. Take p=4p=4 and Lmax≥4L_ ≥ 4. Let G□+G_ ^+ have no observed edges and four exogenous latent variables, with respective child sets 1,2,2,3,3,4,4,1.\1,2\, \2,3\, \3,4\, \4,1\. Its complexity is c(G□+)=8+4=12c(G_ ^+)=8+4=12, and V□:=V(G□+)=Σ∈4:Σ13=Σ24=0,dimV□=8.V_ :=V(G_ ^+)=\ ^4: _13= _24=0\, V_ =8. To verify the claimed Zariski closure, fix the latent variances to one and take nonzero loadings on each permitted latent edge. Each of the four allowed off-diagonal covariances is the product of the two loadings of its own latent, so these four covariances can be varied independently to first order. The four observed noise variances independently vary the diagonal entries and compensate for changes in their latent contributions. The Jacobian of the covariance map, whose entries are partial derivatives of covariance coordinates with respect to the free parameters, therefore has rank eight at such an admissible parameter point with strictly positive observed noise variances. Selecting eight parameter coordinates and keeping the others fixed gives an invertible Jacobian for these eight covariance entries. The inverse function theorem therefore shows that their projection contains a nonempty open subset of ℝ8R^8. Definition 14 and containment in the displayed eight-dimensional linear space therefore give covariance-family dimension eight. Proposition 2 then implies that its Zariski closure has dimension eight; Lemma 16 then shows that it is the entire linear space. We next show that every candidate G+G^+ with V(G+)=V□V(G^+)=V_ has complexity at least twelve. For every candidate graph in p,LmaxG_p,L_ , an observed covariance Σij _ij is identically zero over its admissible parameter space if and only if i and j have no common ancestor, with each vertex counted as its own ancestor. Indeed, a noise source can affect an observed variable only if there is a directed path from the source to that variable, so absence of a common ancestor means that no structural noise contributes to both variables, so their covariance is zero. Conversely, if a common ancestor exists, choose all permitted coefficients strictly positive and sufficiently small that the spectral radius of the full coefficient matrix is less than one. This choice supplies an admissible parameter point; it imposes no additional assumption on the candidate class. The convergent series (I−B)−1=∑k=0∞Bk(I-B)^-1= _k=0^∞B^k then gives a strictly positive coefficient from each noise source to every vertex reachable from it. Positive noise variances make the common ancestor’s contribution to Σij _ij strictly positive, and all other contributions are nonnegative. Hence Σij _ij is not identically zero. This argument also applies to graphs containing observed directed cycles. In the following argument, adjacency refers to the undirected cycle with pairs 1,2,2,3,3,4,4,1\1,2\,\2,3\,\3,4\,\4,1\. These are exactly the pairs whose covariances are not identically zero. This cycle describes covariance pairs, not directed edges among observed variables. The opposite pairs are 1,3\1,3\ and 2,4\2,4\. An observed edge i→ji→ j in G+G^+ must join adjacent vertices of this cycle, since the opposite-pair covariances vanish identically. Let k be the other neighbor of i in the four-cycle; it is opposite j. Since Σik _ik is not identically zero, i and k have a common ancestor. The edge i→ji→ j makes this source an ancestor of j as well, contradicting the identity Σjk=0 _jk=0. Thus G+G^+ has no observed edges. Any two children of the same latent must consequently be adjacent in the cycle, since their covariance is not identically zero. Because no three vertices of the cycle are pairwise adjacent and the candidate class requires at least two children, each latent has exactly two children. Each of the four adjacent pairs needs its own latent source. There are therefore at least four latents and eight latent-to-observed edges, which proves c(G+)≥12c(G^+)≥ 12. In particular, when G∗+=G□+G^*+=G_ ^+, every distribution-equivalent representative has complexity at least twelve, and c†=12c =12. A complete directed acyclic graph on the four observed vertices, denoted Gsat+G_sat^+, has only six edges and generates every positive-definite covariance through sequential linear Gaussian regressions. Thus V□⊊V(Gsat+)=4,c(Gsat+)=6<c†,V_ V(G_sat^+)=S^4, c(G_sat^+)=6<c , which violates Assumption 8. The score comparison is stronger than this failure of a sufficient condition. Every graph quasi-equivalent to G□+G_ ^+ also has Zariski closure V□V_ . To verify this, let U=ℳ(G+)∩ℳ(G□+)U=M(G^+) (G_ ^+) for such a graph. Quasi-equivalence and Proposition 2 give dimU=dimV(G+)=dimV□=8,U⊆V(G+)∩V□. U= V(G^+)= V_ =8, U V(G^+)∩ V_ . If the two Zariski closures were different, their intersection would be a nonempty proper algebraic subset of at least one of them. Both closures are irreducible, so the proper-subset argument in the proof of Lemma 16 would give intersection dimension strictly below eight. This contradicts its containment of U, which has dimension eight. Therefore V(G+)=V□V(G^+)=V_ , and the preceding complexity bound gives c(G+)≥12c(G^+)≥ 12. For any Sn≻0S_n 0 and λn>0 _n>0, sn(Gsat+)=Cn+6λn<Cn+12λn≤sn(G+)whenever G+≅OG□+,s_n(G_sat^+)=C_n+6 _n<C_n+12 _n≤ s_n(G^+) G^+ _OG_ ^+, because δn(Gsat+)=0 _n(G_sat^+)=0 and δn(G+)≥0 _n(G^+)≥ 0 for every graph. Thus the original structural penalty can exclude every quasi-equivalent graph, even with exact global optimization. The quasi-equivalence guarantee therefore cannot extend to all graphs in the candidate class without conditions excluding this obstruction. A dimension-based penalty would avoid this particular strict-inclusion obstruction, but is not the penalty analyzed here. B.3 Weaker and sufficient conditions for structural minimality This subsection discusses a weaker condition at the true covariance and a restricted class for which structural minimality can be verified. The positive result assumes that the target covariance family admits a latent-free representation without reciprocal edge pairs, meaning that i→ji→ j and j→ij→ i never both occur. It does not establish structural minimality for families with no latent-free representation. A weaker condition at the true covariance. Structural minimality is not needed for Theorem 4; it is one sufficient route to the identification condition used in Corollary 6. The main-text assumption compares the target family with every candidate graph. For identification at a fixed Σ∗ ^*, it suffices to make this comparison only for candidates whose covariance-family closures contain Σ∗ ^*. Specifically, Assumption 8 can be replaced in Lemma 17 by Σ∗∈ℳ¯(G+),V(G†)⊊V(G+)⟹c(G+)>c†(G+∈p,Lmax). ^*∈ M(G^+), V(G ) V(G^+) c(G^+)>c (G^+ _p,L_ ). Indeed, every candidate considered in that proof already satisfies Σ∗∈ℳ¯(G+) ^*∈ M(G^+) and c(G+)≤c†c(G^+)≤ c . Algebraic faithfulness gives inclusion of Zariski closures; the displayed condition excludes strict inclusion, and full-dimensional overlap then gives quasi-equivalence. This weaker condition imposes no requirement on candidates whose covariance-family closures do not contain Σ∗ ^*. Restricting the original condition only to c(G+)≤c†c(G^+)≤ c , however, would not weaken it: every violation already has that complexity bound. We retain the global version in the main text because it is a property of the target family and competing graphs, independent of the realized covariance. Neither version can be discarded in general without an alternative identification condition, as Subsection B.2 demonstrates. The role of the complexity penalty. Proposition 2 establishes irreducibility and identifies covariance-family dimensions with those of their Zariski closures. Lemma 16 then shows that strict inclusion V(G†)⊊V(G+)V(G ) V(G^+) is equivalent to inclusion with strictly larger covariance-family dimension. If complexity were instead defined by cdim(G+)=dimℳ(G+)c_dim(G^+)= (G^+), the corresponding structural-minimality condition would hold automatically. This is not the edge-plus-latent complexity used in our score, and a standard BIC interpretation of a dimension-based penalty would still require additional regularity. The simple mixed graphs of Amendola et al. (2020) have at most one edge of any type between each pair of observed vertices. Let k count a graph’s directed edges and its bidirected edges, which permit correlations between structural noises. Their Theorem 3.1 gives covariance-family dimension p+kp+k. With p fixed, a larger dimension therefore requires more edges. Consequently, the analogous structural-minimality condition holds automatically if all candidates are simple mixed graphs and complexity counts their edges. This is a consequence of their dimension theorem, not a statistical recovery guarantee. In our explicit-latent class, the edge-plus-latent count need not increase with covariance-family dimension, so structural minimality requires a separate justification. Structural complexity and covariance dimension. Our penalty counts the edges and latent vertices in a graph; it does not count the independently varying covariance entries. For example, consider two graphs on two observed variables: one has a single latent parent of both variables and no observed edge; the other has only the observed edge 1→21→ 2. Both generate all of ++2S_++^2, a family of dimension three, but their structural complexities are three and one, respectively. The construction is given at the end of this subsection. This distinction is deliberate. The structural penalty expresses a preference for fewer edges and latent vertices, and its additive form allows the expected penalty under Bernoulli gates to be evaluated without computing model dimensions for the sampled graphs. It also introduces a limitation: greater covariance-family dimension need not imply greater structural complexity, as Subsection B.2 demonstrates. Replacing this penalty by model dimension would change the objective and its preference among graphs, rather than merely remove a redundant parameterization. A dimension criterion. The original free parameter coordinates consist of one coefficient per permitted edge and p+ℓ(G+)p+ (G^+) noise variances, giving p+c(G+)p+c(G^+) coordinates. A sharper dimension bound follows by absorbing latent variances into the loadings. Set Γ=ΩL1/2Λ = _L^1/2 , where the diagonal entries of ΩL1/2 _L^1/2 are the positive square roots of the latent variances. Then Λ⊤ΩLΛ=Γ⊤Γ _L = , and Γ has the same graph-required zeros as Λ . Conversely, any such Γ is realized by Λ=Γ = and ΩL=Iℓ _L=I_ . Thus the covariance family is also the image of (BOO,Γ,ΩO)⟼(Ip−BOO)−⊤(ΩO+Γ⊤Γ)(Ip−BOO)−1.(B_O, , _O) (I_p-B_O)^- ( _O+ )(I_p-B_O)^-1. Its domain has p+|E(G+)|p+|E(G^+)| free coordinates: one per permitted edge and one per observed noise variance. Positive observed variances and det(Ip−BOO)≠0 (I_p-B_O)≠ 0 define a nonempty Euclidean open semialgebraic subset of this coordinate space. The displayed map is rational and semialgebraic on that domain. The dimension of a semialgebraic image cannot exceed the dimension of its domain (Coste, 2000, Theorem 3.18), so dimℳ(G+)≤p+|E(G+)|=p+c(G+)−ℓ(G+). (G^+)≤ p+|E(G^+)|=p+c(G^+)- (G^+). (35) This is an upper bound, not a claim that the remaining parameters are identifiable or that their number equals model dimension. The original free coordinates remain valid for the geometry and genericity proofs. Consequently, if the target family satisfies dimℳ(G∗+)=p+c† (G^*+)=p+c , any candidate of strictly larger model dimension must have c(G+)>c†c(G^+)>c . This verifies structural minimality even without requiring inclusion of Zariski closures. The following class satisfies this criterion. Its dimension calculation is the diagonal-noise directed special case of the simple cyclic mixed-graph dimension result of Amendola et al. (2020, Theorem 3.1 and Lemma 3.3). We give a direct proof and combine this dimension equality with Equation 35 to verify structural minimality against our full candidate class, including graphs with latent variables or reciprocal observed edges. The competitors need not themselves admit a latent-free representation. Proposition 22 (A class with a latent-free representative). Suppose ℳ(G∗+)=ℳ(Gdir+)M(G^*+)=M(G_dir^+) for a latent-free graph Gdir+∈p,LmaxG_dir^+ _p,L_ with e observed edges and at most one of i→ji→ j and j→ij→ i for each distinct pair. Directed cycles are allowed. Then dimℳ(G∗+)=p+e,c†=e, (G^*+)=p+e, c =e, and Assumption 8 holds against every candidate in p,LmaxG_p,L_ , including graphs with latent variables or reciprocal edges. Under the sampling conditions of Theorem 4, global score minimizers therefore recover marginal quasi-equivalence for generic admissible parameters of G∗+G^*+. Proof. We first determine the covariance-family dimension and then compare structural complexities. Dimension of the latent-free covariance family. List the observed edges of Gdir+G_dir^+ as i1→j1,…,ie→jei_1→ j_1,…,i_e→ j_e. Write their coefficients as bk=(BOO)ikjkb_k=(B_O)_i_kj_k and the observed noise variances as ωi=(ΩO)ii _i=( _O)_i. The free parameter vector is θ=(ω1,…,ωp,b1,…,be)∈ℝp+e.θ=( _1,…, _p,b_1,…,b_e) ^p+e. Let U be its admissible domain, defined by ωi>0 _i>0 for all i and det(Ip−BOO)≠0 (I_p-B_O)≠ 0. These strict conditions make U open in ℝp+eR^p+e. The covariance map ϕ:U⟶++p,ϕ(θ)=(Ip−BOO)−⊤ΩO(Ip−BOO)−1φ:U _++^p, φ(θ)=(I_p-B_O)^- _O(I_p-B_O)^-1 is rational, continuously differentiable, and semialgebraic on U, with image ℳ(Gdir+)M(G_dir^+). The parameter point θ0=(1,…,1,0,…,0) _0=(1,…,1,0,…,0) corresponds to (BOO,ΩO)=(0,Ip)(B_O, _O)=(0,I_p) and lies in U. This point is admissible because an allowed edge coefficient may equal zero. For a parameter direction θ˙ θ, write Dϕ(θ0)[θ˙]:=dtϕ(θ0+tθ˙)|t=0Dφ( _0)[ θ]:= . ddtφ( _0+t θ) |_t=0. Differentiating the covariance map gives Dϕ(θ0)[θ˙]=Ω˙+B˙⊤+B˙,Dφ( _0)[ θ]= + B + B, where B˙ B and Ω˙ denote changes in the edge coefficients and diagonal noise variances, respectively. Indeed, the derivative of (Ip−tB˙)−1(I_p-t B)^-1 at t=0t=0 is B˙ B, and the product rule gives the three displayed terms. At θ0 _0, the partial derivative with respect to ωi _i is therefore EiiE_i, while that with respect to bkb_k is Eikjk+EjkikE_i_kj_k+E_j_ki_k; here EijE_ij has one unit entry in position (i,j)(i,j) and zeros elsewhere. To connect this derivative calculation to our definition of dimension, select the following covariance coordinates: π(Σ)=(Σ11,…,Σpp,Σi1j1,…,Σieje)∈ℝp+e.π( )=( _11,…, _p, _i_1j_1,…, _i_ej_e) ^p+e. Symmetry allows each off-diagonal entry to be identified with its upper-triangular coordinate. Since no reciprocal edge pair is present, these e off-diagonal coordinates are distinct. Thus π is a coordinate projection up to a reordering of its output coordinates, which preserves open sets. Define F=π∘ϕF=π φ, which maps U into ℝp+eR^p+e. In the parameter order specified above, its Jacobian is DF(θ0)=Ip+e.DF( _0)=I_p+e. Thus the inverse function theorem applies: there are open neighborhoods U0⊆U_0 U of θ0 _0 and W0⊆ℝp+eW_0 ^p+e of F(θ0)F( _0) such that F maps U0U_0 one-to-one onto W0W_0, with a continuously differentiable inverse. In particular, W0=F(U0)⊆π(ℳ(Gdir+)).W_0=F(U_0) π (M(G_dir^+) ). The projected covariance family therefore contains a nonempty open ball in ℝp+eR^p+e. Definition 14 gives dimℳ(Gdir+)≥p+e (G_dir^+)≥ p+e. Conversely, the semialgebraic-image dimension bound stated above gives dimℳ(Gdir+)=dimϕ(U)≤dimU=p+e. (G_dir^+)= φ(U)≤ U=p+e. Combining the two inequalities, and using the assumed equality of covariance families, yields dimℳ(G∗+)=dimℳ(Gdir+)=p+e. (G^*+)= (G_dir^+)=p+e. The Jacobian calculation thus establishes the dimension using the existing definition; it does not introduce a different definition or assume that the Jacobian has the same rank at every parameter value. Minimum complexity and structural minimality. For any distribution-equivalent candidate G+G^+, the universal dimension bound gives p+e=dimℳ(G+)≤p+c(G+)−ℓ(G+).p+e= (G^+)≤ [rgb]0,0,0p+c(G^+)- (G^+). Thus c(G+)≥e+ℓ(G+)≥e [rgb]0,0,0c(G^+)≥ e+ (G^+)≥ e, while Gdir+G_dir^+ itself has complexity e. It follows that c†=ec =e. For every candidate of strictly larger dimension, p+c(G+)−ℓ(G+)≥dimℳ(G+)>dimℳ(G∗+)=p+e, [rgb]0,0,0p+c(G^+)- (G^+)≥ (G^+)> (G^*+)=p+e, so c(G+)>e=c†c(G^+)>e=c . This proves global structural minimality without restricting the competing graphs. Lemmas 21 and 17 then establish the identification condition for generic true parameters, and Corollary 6 gives the final assertion. ∎ A cyclic example with covariance-family dimension eight. Take four observed vertices, no latents, and the directed cycle 1⟶2⟶3⟶4⟶1.1 2 3 4 1. Proposition 22 gives c†=4c =4 and dimℳ(G∗+)=8<10=dim4 (G^*+)=8<10= ^4. Its covariance family therefore has smaller dimension than the space of all positive-definite covariances. Any strictly larger Zariski closure must have complexity at least five, even when latent competitors are allowed. For generic true parameters, global score minimizers recover marginal quasi-equivalence. This conclusion does not identify a unique orientation or give a guarantee for a local optimizer. A latent example with an observed representative. A single exogenous latent with exactly two observed children and no observed edge generates every positive-definite covariance on that pair. To see this, write the desired covariance as [abbd]≻0 [ smallmatrixa&b\\ b&d smallmatrix ] 0. Positive definiteness gives d>0d>0 and ad>b2ad>b^2, so b2/d<ab^2/d<a. For b≠0b≠ 0, choose u2∈(b2/d,a)u^2∈(b^2/d,a) and set v=b/uv=b/u. Unit latent variance, loadings (u,v)(u,v), and positive observed noise variances (a−u2,d−v2)(a-u^2,d-v^2) give the desired covariance. For b=0b=0, set both loadings to zero, which is admissible because permitted coefficients may be zero, and use observed noise variances (a,d)(a,d). The latent-free graph 1→21→ 2 generates the same covariance using edge coefficient b/ab/a and noise variances a and d−b2/ad-b^2/a, which are positive. Both families therefore equal ++2S_++^2 and have dimension three. The latent graph has five original parameter coordinates and four after absorbing its latent variance, so even the normalized parameter count exceeds its model dimension. The proposition applies with c†=1c =1, although the latent graph has complexity three. Disjoint unions of such pairs, isolated vertices, and the directed graphs in the proposition also admit representatives of the same type. Thus the data-generating graph may contain both directed cycles and latent variables in separate components. The conclusion does not establish structural minimality for families that require latent variables in every representation, or identify the number of latent variables in the data-generating graph. Appendix C Practical optimization This appendix describes an optional approximation for numerical optimization of the exact Bernoulli objective in Equation 21. For latent-variable models, we use the Binary Concrete construction of Maddison et al. (2017) to obtain a differentiable surrogate for the discrete gates. This approximation is not required for the differentiability of the exact expectation or for Proposition 10. Let σ(t)=(1+exp(−t))−1σ(t)=(1+ (-t))^-1 and parameterize qij=σ(αij)q_ij=σ( _ij), rhj=σ(βhj)r_hj=σ( _hj), and πh=σ(γh) _h=σ( _h). If G=logU−log(1−U)G= U- (1-U) with U∼Uniform(0,1)U (0,1), the corresponding Binary Concrete gate is K~=σ(a+Gτ), K=σ\! ( a+Gτ ), (36) where a is the relevant logit: a probability q∈(0,1)q∈(0,1) is represented by a=log(q/(1−q))a= (q/(1-q)). The temperature τ>0τ>0 controls how closely the continuous gate approximates a binary gate. Applying Equation 36 to each nonconstant gate, with independent noise draws and M~ii=0 M_i=0, gives M~,H~,Z~ M, H, Z. These gates are differentiable functions of the logits for each fixed noise draw; the noise distribution does not depend on the parameters. For each realization with invertible Ip−M~⊙WI_p- M W, we can therefore differentiate the sampled likelihood while holding the noise draws fixed. The approximate objective at temperature τ is Fn,τ=[Ln(M~,H~,Z~,ϑ)]+λn[∑i≠jqij+∑hπh(1+∑jrhj)]. splitF_n,τ=&E [L_n( M, H, Z; ) ]\\ &+ _n\! [ _i≠ jq_ij+ _h _h (1+ _jr_hj ) ]. split (37) Invertibility of binary and continuous masks. Invertibility for every binary mask does not imply invertibility for every continuous mask. For example, with two observed variables and W=(0220),det(I2−M⊙W)=1−4M12M21,W= pmatrix0&2\\ 2&0 pmatrix, (I_2-M W)=1-4M_12M_21, every binary mask gives determinant one or minus three. In contrast, the continuous mask with M~12=M~21=1/2 M_12= M_21=1/2 gives det(I2−M~⊙W)=0 (I_2- M W)=0. Binary Concrete therefore does not guarantee invertibility for all continuous masks. Using sampled derivatives as an unbiased gradient estimator for Fn,τF_n,τ additionally requires justification for interchanging differentiation and expectation, for example through an integrable bound on the derivatives that is uniform over a parameter neighborhood. Invertibility of Ip−M~⊙WI_p- M W with probability one is not sufficient: the absolute derivatives can have infinite expectation when this matrix approaches a noninvertible matrix. Thus differentiability for individual noise realizations does not by itself justify differentiating the expectation. We deliberately retain the expected discrete complexity in closed form: at positive temperature, [M~ij]≠qijE[ M_ij]≠ q_ij in general. Replacing the closed-form penalty by a sum of sampled continuous gate values would therefore change the complexity penalty and introduce additional sampling variation. For numerical optimization, a positive lower bound on each observed-noise variance can be imposed through (ΩO)jj=ε+softplus(dj)( _O)_j= +softplus(d_j) with ε>0 >0, where softplus(t)=log(1+et)softplus(t)= (1+e^t) and dj∈ℝd_j is optimized. To assess proximity to a noninvertible matrix, one can monitor the smallest singular value σmin(Ip−M~⊙W) _ (I_p- M W), where σmin(A)=min‖x‖2=1‖Ax‖2 _ (A)= _\|x\|_2=1\|Ax\|_2. Monitoring alone does not justify interchanging differentiation and expectation. The optional restriction ‖|W|‖2≤1−εA\|\,|W|\,\|_2≤ 1- _A, with 0<εA<10< _A<1, gives ‖M~⊙W‖2≤‖|W|‖2\| M W\|_2≤\|\,|W|\,\|_2 and hence σmin(Ip−M~⊙W)≥εA _ (I_p- M W)≥ _A for every binary mask or continuous mask with entries in [0,1][0,1]. Together with the lower bound on noise variances and bounds on the remaining parameters in a neighborhood of their current values, this supplies derivative bounds on that neighborhood at fixed τ>0τ>0. Both restrictions are stronger than the statistical model and can exclude admissible covariances or well-defined cyclic systems. Numerical implementation used in the experiments. The implementation in Section 5 differs between the observed-only and latent-variable settings. When no latent variables are present, the latent activation and latent-edge gates disappear, and we optimize the exact Bernoulli formulation in Equation 21 using 256 Monte Carlo samples of hard observed-edge masks per update. Gradients of the observed-edge logits are estimated using a leave-one-out score-function estimator. We additionally use a discrete compound support-refinement step: when the score plateaus, deletion and reversal-plus-deletion candidates are evaluated, accepted deletions are permanently removed from the search space, and optimization continues from the best improving support. When latent variables are present, we use a straight-through Binary Concrete estimator. The forward likelihood is evaluated using hard binary gates, while the backward pass differentiates through their soft Binary Concrete counterparts. Thus the experimental implementation does not evaluate the forward likelihood on fractional graph masks. The value of τ is also fixed throughout training. This estimator is a numerical approximation and should be distinguished from Proposition 10, which concerns the exact Bernoulli expectation. Discrete support selection and refitting. After either optimization procedure, the learned gate probabilities and continuous parameters are converted to a binary support (M,H,Z)(M,H,Z). Inactive latent slots and active slots with fewer than two selected children are removed using the conversion in Appendix A.11. For each resulting fixed graph, we refit (W,Γ,ΩO)(W, , _O) using the original Gaussian negative log-likelihood. To target the discrete score in Equation 13, this refit optimizes over positive diagonal ΩO _O and invertible A(M)=Ip−M⊙WA(M)=I_p-M W without imposing the optional norm restriction on W. Candidate graphs are then compared using the refitted discrete penalized score. Exact thresholds, initialization, optimization budgets, and refitting settings are reported in Appendix D. These numerical procedures do not change the scope of the theoretical guarantees. Theorem 4 applies when the resulting graph is a global minimizer of the discrete score, and Corollary 6 additionally requires Assumption 5. Guarantees for Monte Carlo optimization, straight-through gradients, discrete support refinement, and finite numerical refitting require additional analysis. Appendix D Additional experimental details and sensitivity analyses D.1 Synthetic data generation Synthetic observations follow the linear Gaussian SCM in Section 2 with n=10,000n=10,000 training samples per trial. We use the observed-graph generation procedure from the public implementation of Ghassami et al. (2020). The parameter called maximum degree is a hard bound on total observed degree, i.e., observed in-degree plus observed out-degree; self-loops are excluded. For every selected observed edge, the coefficient sign is sampled uniformly and its magnitude is sampled uniformly from [0.2,0.8][0.2,0.8]. Coefficient draws are accepted only when the observed coefficient matrix satisfies spectral radius below 0.990.99 and cond(I−BOO)≤20cond(I-B_O)≤ 20. Observed noise variances are drawn independently from Uniform[1,3]Uniform[1,3]. Each exogenous latent variable independently selects an integer number of distinct observed children uniformly from 2,3,4\2,3,4\ and then selects that many observed nodes uniformly without replacement; child sets of different latent variables may overlap. Nonzero latent loadings use the same random-sign and [0.2,0.8][0.2,0.8] magnitude distribution as observed coefficients. Latent variances are normalized to one, with their scale absorbed into the latent loading rows. Latent percentages specify the nominal ratio ℓ/p /p. We evaluate ℓ/p∈0.1,0.2,0.3 /p∈\0.1,0.2,0.3\ and resolve the integer count as ⌊p(ℓ/p)+0.5⌋ p( /p)+0.5 , with at least one latent whenever the requested ratio is positive. The estimator receives Lmax=2⌊p(ℓ/p)+0.5⌋+1L_ =2 p( /p)+0.5 +1 latent slots for latent experiments and Lmax=0L_ =0 in the no-latent experiments. We uses 10 trials with different seeds. D.2 Optimization and evaluation settings The numerical procedures for the exact-Bernoulli and straight-through Binary Concrete searches are described in Appendix C. Table 2 summarizes the principal implementation, graph-extraction, refitting, and evaluation settings used in the experiments. During optimization, Adam is run for 10001000, 15001500, 30003000, 40004000, 60006000, 80008000, 1000010000, and 1000010000 updates for p=4,8,16,32,64,128,256p=4,8,16,32,64,128,256, and 512512, respectively. For the GNW experiments with p=100p=100, we use 1000010000 updates. Table 2: Principal numerical settings used in the experiments. Setting Value Structure optimization Mask samples per update 256256 Adam learning rate 10−210^-2 Binary Concrete temperature τ=1τ=1 (fixed) Structure-search initializations 11 per trial Initial observed-edge probability 0.010.01 Initial latent-edge probability 0.50.5 Initial latent activation probability 0.990.99 Initial coefficient scale 0.50.5 Candidate latent bound Lmax=2⌊p(ℓ/p)+0.5⌋+1L_ =2 p( /p)+0.5 +1 Graph extraction and refitting Gate probability threshold 0.50.5 Coefficient threshold 10−410^-4 Compatibility evaluation Covariance targets per direction Keval=20K_eval=20 Compatibility-fit starts 55 Compatibility maximum iterations 10001000 Compatibility tolerance 10−810^-8 Observed-only support refinement Plateau window 5050 updates Relative score tolerance 10−410^-4 Edges deleted per branch 11 Deletion branches 1010 Reversal-plus-deletion branches 1010 After graph extraction, inactive latent slots and active slots with fewer than two retained children are removed. Continuous parameters are then refit on each fixed support using the unrelaxed Gaussian likelihood before the discrete score and covariance-family compatibility are evaluated. D.3 Bidirectional compatibility evaluation For evaluation we use the symmetric Gaussian divergence J(S,Σ)=12D(S∥Σ)+D(Σ∥S),J(S, )= 12 \D(S\| )+D( \|S) \, where D is defined in Equation 12. Let G∗+G^*+ denote the reference graph and G^+ G^+ the recovered graph. For generic admissible parameter values θ1∗,…,θK∗∈(G∗+) _1^*,…, _K^* (G^*+), define Σk∗=ΣO(θk∗,G∗+) _k^*= _O( _k^*;G^*+) and dk∗→G^:=infθ^∈(G^+)J(Σk∗,ΣO(θ^,G^+)).d_k^*→ G:= _ θ ( G^+)J\! ( _k^*, _O( θ; G^+) ). (38) The forward discrepancy is dfwd:=max1≤k≤Kdk∗→G^.d_fwd:= _1≤ k≤ Kd_k^*→ G. (39) For the reverse direction, choose generic admissible parameter values θ~1,…,θ~K∈(G^+) θ_1,…, θ_K ( G^+) and let Σ~k=ΣO(θ~k,G^+) _k= _O( θ_k; G^+). Define dkG^→∗:=infθ∈(G∗+)J(Σ~k,ΣO(θ;G∗+)),d_k G→*:= _θ (G^*+)J\! ( _k, _O(θ;G^*+) ), (40) and drev:=max1≤k≤KdkG^→∗.d_rev:= _1≤ k≤ Kd_k G→*. (41) We use K=20K=20 in the reported experiments and summarize bidirectional compatibility by dfwd+drevd_fwd+d_rev. The continuous parameters of the fitted graph are reoptimized independently for each target covariance while its graph support remains fixed. D.4 Baseline implementations and runtime protocol Ghassami et al. (Ghassami et al., 2020) is evaluated using the public dglearn implementation with its original configuration: hill-climbing search, one initialization, at most 1000 iterations, maximum cycle/SCC parameter five, virtual refinement, and support reduction. DCCD-CONF (Sethuraman and Fekri, 2025) is evaluated using its public implementation with its provided hyperparameters: 500 epochs, batch size 512, learning rate 10−210^-2, λc=ρ=10−2 _c=ρ=10^-2, Lipschitz constant 0.90.9, linear activation, adjacency threshold 0.80.8, and confounding-covariance threshold 0.010.01. These baseline configurations are default values used in their corresponding repository. For DCCD-CONF, the observational data are supplied as a single regime with no intervention targets so that no additional interventional information is provided. The Amendola et al. (2020) baseline follows the greedy cyclic mixed-graph search with add, remove, and edge-reversal moves; to keep its neighborhood search computationally feasible, at most 64 candidate branches are evaluated at each step. All methods evaluated at a shared setting receive identical training data and trial seeds. All experiments were run on compute nodes equipped with AMD EPYC 9334 32-Core processors and 8 GB of memory. DCCD-CONF and Améndola et al. were evaluated up to p=64p=64, while dglearn was evaluated up to p=128p=128, beyond which their computation became infeasible for the experiments. D.5 GeneNetWeaver details We construct fixed GNW reference networks by first generating a weakly connected directed acyclic graph over the observed variables with a bounded total degree. Directed cycles are then introduced by adding back edges while preserving the prescribed degree constraint. Exogenous latent variables are subsequently added before simulation. Each latent variable has no parents and is connected to a randomly selected subset of observed variables, allowing different latent variables to share observed children. Data are generated using the GNW multifactorial simulator with all perturbation inputs set to zero, so that no observed variable is directly intervened upon. Each fixed network yields 20,000 observations. In each trial, 16,000 observations are used for structure learning and 4,000 are held out. The network topology and simulated dataset remain fixed across trials, while the train–test partition, stochastic optimization, mask sampling, and compatibility-evaluation draws vary. D.6 Sensitivity to graph density and latent confounding We stress test recovery over all nine combinations of maximum observed degree 1,2,4\1,2,4\ and latent ratio ℓ/p∈0.1,0.2,0.3 /p∈\0.1,0.2,0.3\. These experiments vary the structure of the data-generating graphs and therefore evaluate sensitivity to graph density and latent confounding. (a) 10%10\% latent, degree 1 (b) 10%10\% latent, degree 2 (c) 10%10\% latent, degree 4 (d) 20%20\% latent, degree 1 (e) 20%20\% latent, degree 2 (f) 20%20\% latent, degree 4 (g) 30%30\% latent, degree 1 (h) 30%30\% latent, degree 2 (i) 30%30\% latent, degree 4 Figure 3: Sensitivity to graph density and latent confounding. Rows correspond to latent ratios ℓ/p /p of 10%10\%, 20%20\%, and 30%30\%, and columns to maximum observed degrees of 11, 22, and 44. Points show mean dfwd+drevd_fwd+d_rev with SEM across trials; lower is better. Effect of density and latent confounding. Increasing the maximum observed degree permits denser cyclic structure, while increasing ℓ/p /p introduces more exogenous common causes and enriches the disturbance term Λ⊤ΩLΛ _L . The results show that recovery depends jointly on graph size, density, and latent structure rather than on any single factor alone. Because confounding also depends on the sampled child sets and loading magnitudes, ℓ/p /p should be interpreted as the number of latent common causes rather than a direct measure of confounding strength. Each baseline is evaluated up to the largest graph size for which its computation remains feasible. Appendix E Related work Distributional comparisons for cyclic Gaussian models. Ghassami et al. (2020) study cyclic linear Gaussian models without latent confounders. They compare the distribution families induced by compatible parameters, define quasi-equivalence through positive-measure overlap, and prove consistency of global penalized-likelihood minimizers under generalized faithfulness and an additional structural condition (Ghassami et al., 2020, Definition 9, Assumption 1, Theorem 3). We extend the overlap-based perspective to observed covariance families after marginalizing latent variables, requiring an intersection with the dimension of both families; Appendix A.6 explains the distinction. In the same no-latent setting, Yi and Lee (2024) propose Filter, Rank, and Prune. Under their assumptions, filtering retains every edge of at least one objective-minimizing graph with high probability for sufficiently large samples (Yi and Lee, 2024, Theorem 3.4). This is a filtering guarantee, not a quasi-equivalence consistency result for the final graph. Discovery with cycles and latent confounding. Hyttinen et al. (2012, Theorem 12) use randomized interventions to identify observed coefficients and structural-noise covariances in linear cyclic models under intervention and solvability conditions. Correlated noises represent confounding without specifying individual latent variables. Forré and Mooij (2018) study nonlinear cyclic models with latent confounders using observational and interventional conditional-independence information. They assume that these independences agree exactly with σ-separation, their graphical criterion for cyclic models; linear models need not satisfy this assumption (Forré and Mooij, 2018, Definition 2.15, Remark 2.16). We instead use observational data and compare complete Gaussian covariance families. Comparison with Améndola et al. Amendola et al. (2020, Section 2) study cyclic linear Gaussian mixed graphs. Bidirected edges permit correlated structural noises. Their simple-graph restriction allows at most one edge of any type per vertex pair, excluding reciprocal directed edges and directed–bidirected pairs. Our explicit latent model allows both reciprocal edges and a directed edge between variables sharing a latent parent, but the different representations of confounding need not define the same covariance families. Their Theorem 3.1 gives dimension p+kp+k for every simple mixed graph, where k counts directed and bidirected edges. Their Theorem 4.1 gives a graphical sufficient condition for two such graphs to have equal Euclidean covariance-family closures. Neither theorem establishes statistical consistency of their score or search. Our results instead concern selection by global score minimizers as the sample size increases, with recovery up to marginal quasi-equivalence under identification conditions. Their proposed penalty adds an edge-count term to dimension-based BIC (Amendola et al., 2020, Equations (15)–(16)); ours counts directed edges and latent variables. On latent-free directed graphs without reciprocal edges, ordinary BIC and our score select the same globally optimal graphs; this equivalence need not hold with their additional edge-count penalty. Appendix B.3 relates their dimension theorem to our structural-minimality condition for targets admitting such a latent-free representative, allowing all our latent and cyclic competitors. It does not verify structural minimality for general confounded targets. Their algorithm uses greedy edge additions, deletions, and reversals with parameter fitting (Amendola et al., 2020, Section 5). We optimize Bernoulli inclusion probabilities, although our implementation also uses discrete refinement in no-latent experiments and fixed-graph refitting (Appendix C). Latent confounding in acyclic models. Kaltenpoth and Vreeken (2023) learn independent exogenous latent confounders in acyclic linear Gaussian models. Their identification results restrict the sizes and incoming edges of groups sharing a latent confounder; identifying observed directed coefficients additionally requires equal observed-noise variances (Kaltenpoth and Vreeken, 2023, Sections 2.1–2.3, Theorem 2). SPOT uses mixed graphs without directed cycles and represents confounding through bidirected edges, rather than optimizing an explicit latent count (Ma et al., 2024, Sections 2–3). Our model allows cycles and unequal noise variances, with recovery up to marginal quasi-equivalence under identification conditions rather than unique identification. Continuous structure optimization. NOTEARS uses a smooth equality constraint characterizing acyclic coefficient matrices (Zheng et al., 2018, Theorem 1). We require invertibility of Ip−BOOI_p-B_O, not acyclicity. DCCD-CONF learns nonlinear cyclic models from interventional data, with correlated Gaussian structural noises and Bernoulli adjacency gates (Sethuraman and Fekri, 2025, Sections 2–3). Its straight-through estimator evaluates binary masks but computes derivatives through continuous gates. Its theoretical result uses a population score, with likelihood averaged under the data-generating distributions. Under its assumptions and a suitable penalty, graphs selected by exact maximizers admit the same family of distributions as the true graph for the supplied interventions (Sethuraman and Fekri, 2025, Theorem 2, Definitions A.8 and A.10). Our guarantee instead uses observational data and concerns increasing sample sizes. Bernoulli graph gates are therefore not novel here; our formulation additionally gates explicit latent activation and connections and proves equality of global infima with the corresponding discrete objective. The Concrete distribution provides continuous approximations to discrete variables (Maddison et al., 2017, Section 3.2, Appendix B). Appendix C distinguishes fully continuous masks from straight-through optimization and describes our implementation. Proposition 10 concerns the exact Bernoulli expectation, not either approximation, and does not guarantee numerical convergence to a global minimum.