An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores for Complex Non-linear and Categorical Data
Authors: Nikita Lagrange
Herve Isambert
Venue: ICML 2025 poster (2025)
Keywords: causal discovery
search-and-score structure learning
latent variable
multivariate information
Primary area: probabilistic_methods->structure_learning
TL;DR: We propose an efficient search-and-score algorithm for ancestral graphs based multivariate information scores

Abstract: We propose a greedy search-and-score algorithm for ancestral graphs, which include directed as well as bidirected edges, originating from unobserved latent variables. The normalized likelihood score of ancestral graphs is estimated in terms of multivariate information over relevant "$ac$-connected subset" of vertices, $\boldsymbol{C}$, that are connected through collider paths confined to the ancestor set of $\boldsymbol{C}$. For computational efficiency, the proposed two-step algorithm relies on local information scores limited to the close surrounding vertices of each node (step 1)  and edge (step 2). This computational strategy, although restricted to information contributions from $ac$-connected subsets containing up to two-collider paths, is shown to outperform state-of-the-art causal discovery methods on challenging benchmark datasets.

An Efficient Search-and-Sowe Algorithm for Ancestral Graphs using Multivariate Information Scores for Complex Non-linear and Categorical Data
Nikita Lagrange'
Hervé Isambert
Abstract
the asymptotic probability that the N independent sam- ples, xx((), are drawn from the model distri- bution, g(x)), i.e. Lp.g = q(xx11) ,.. xx()) = II ((x:()) rather than the empirical distribution, p(z). This leads to, log LDIg = lo [0,00((x(1))) which converges towards N "x p(æ) log g(z)
We propose a greedy search-and-scon algorithm for ancestral graphs, which include directed as well as bidirected edges, originating from unob- served latent variables. The normalized likeli- hood score of ancestral graphs is estimated in terms of multivariate information over relevant "ac-connected subsets" of vertices, C, that are connected through collider paths confined the ancestor of C. For computational efficiency, the proposed two-step algorithm relies local in- formation scores limited to the close surrounding vertices of each node (step 1) and edge (step 2). This computational strategy, although restricted to information contributions from ac-connected sub- sets containing up to two-collider paths, shown to outperform state-of-the-art causal discovery methods on challenging benchmark datasets.
-N Hp. q) in the large sample size limit, N + oc, with log ZDD Ollog N).
The structural constraints of the model G translate into the factorization form of the theoretical probability distribution, (x:) (Pearl & Paz
1985 Pearl 1988 2009 Spirtes et al.) 2000 Richardson 2009). In particular, the probability dis-
tribution of Bayesian networks (BN) factorizes in terms of conditional probabilities of each variable given its parents, Jnx (z) II, ( |pax)). where pax, denote the val- ues of the parents of node X; in G. Pax.- For Bayesian net- works, the factors of the model distribution, 9(z; |pax,). ). can be directly estimated with the empirical conditional proba- bilities of each node given its parents as, (zx:paxx)) = p(r; |pax,) leading to the well known estimation of the likelihood function in terms of conditional entropies HX:: = - En p(I;, pax,) ) log p(z;)pax;;)
1. Introduction
The likelihood function plays a central role in the selec- tion of graphical model G based on observational data D. Given N independent samples from D, the likelihood LDig that they might have been generated by the graphical model is given by (Koller& Friedman, 2009).
vertices
LDIGEN ZD.
exp (-N Σ H(X,| (2)
X;EV
This paper concerns the experimental setting for which some variables of the underlying Bayesian model are not observed. This frequently occurs in practice for many applications. We derive an explicit likelihood function for the class of ances- tral graphs, which include directed well as bidirected edges, arising from the presence of unobserved latent vari- ables. Tian and Pearl 2002 showed that the probability distribution of such graphs factorizes into c-components including subsets of variables connected through bidirected paths (i.e. containing only bidirected edges). Richardson 2009 later proposed refined factorization of the model distribution of the broader class of acyclic directed mixed graphs in terms of conditional probabilities over "head" and "tail" subsets of variables within each ancestrally closed sub- sets of vertices. However, unlike with Bayesian networks, the contributions of c-components or head-and-tall factors to the likelihood function cannot simply be estimated in terms of empirical distribution p( p(x:, as shown below. This
1 ZD.
exp (-NH(..g)
CDgg
(1)
where Hp).g) 9) = Ea p(x) log (αx) is the cross-entropy between the empirical probability distribution p(æ) of the observed data D and the theoretical probability dis- tribution (xx) of the model G and ZD, data- and model-dependent factor ensuring proper normalization con- dition for finite dataset. In short, Eq 1 results from
'CNRS, UMR168. Institut Curie, PSL University, Sorbonne University, Paris, France. Correspondence to: Hervé Isambert therve.sambbert@curme
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
leaves the likelihood function of ancestral graphs difficult to estimate from empirical data, in general, although iterative methods have been developped when the data is normally distributed (Richardson & Spirtes) 2002 Drton et al. 2009
McGill
1954
Ting
1962
Han
1980
Yeung 1991> as,
T(X) = H(X)
I(X;Y) = H(X) + H(Y)) H(X,)
Evans & Richardson 2010
Triantafillou & T'amaadinos
I(V) -- --1 HS))
2016
Rantanen et al. 2021
Claassen & Bucur
2022).
sgV
The present paper provides an explicit decomposition of the likelihood function of ancestral graphs in terms of multivari- ate cosssinffrratiin over relevant 'ac-connected" subsets of variables, Fig. I which do not rely on the head-and-tail factorization but coincide with the parametrizing sets (Hu & Evans 2020) derived from the head-and-tail factorization. It suggests a natural estimation of these revelant contributions to the likelihood function in terms of empirical distribu- tion p(æ). This result extends the likelihood expression of Bayesian Networks (Eq. 2 to include the effect of unob- served latent variables and enables the implementation of a greedy search-and-scon algorithm for ancestral graphs. For computational efficiency, the proposed two-step algorithm relies on local information scores limited to the surrounding vertices of each node (step 1) and edge (step 2). This com- putational strategy shown to outperform state-of-the-art causal discovery methods on challenging benchmarks. This paper extends a manuscript opted-in at NeurIPS 2024 (La- grange & Isambert 2024) to include benchmarks on linear as well as non-linear continuous models, Figs. 2 & B.2
where the semicolon separators are needed to distinguish multipoint (cross) information from joint variables as in T([X, Z): = T(X;Y) + T(Z;Y) - T(X; Y; Z). Be- low, implicit separators between non-conditioning vari- ables in multivariate (cross) information will always cor- respond to semicolons, e.g. as in I(V) in Eq. 5 Unlike multivariate (cross) entropies, which are always positive, H .. Xx)) > O, multivariate (cross) information, [(): Xx), can be positive or negative for ls > 3. while they remain always positive for k < 3, e.g. I(X;Y) > 0.
In turn, multivariate (cross) entropies can be expressed through the Principle of expression but in terms of multivariate (cross) information,
Exclusion into the same
- Σ --V==- I(S), SCV
(6)
Conditional multivariate (cross) information I(V Z) are defined similarly as multivariate (cross) information ((V) but in terms of conditional (cross) entropies as,
IV 2)
- SCV
E(1)) H(S)Z)
2. Theoretical Results
2.1. Multivariate Cross-entropy and Cross-information
Eqs.
5
7 lead to a decomposition rule relative to a variable which can conditioned on a set of joint variables,
Z. Eq
The theoretical result of the paper (Theorem 1) is expressed in terms of multivariate cross-information derived from mul- tivariate cross-entropies through the Inclusion-Exclusion Principle. The same expressions can be written between multivariate information and multivariate entropies by sim- ply substituting q( ((zz)) with p((x;}) in the equations be- low and will be used to estimate the likelihood function of ancestral graphs (Proposition 3).
A
A.... Amm. with implicit comma separators,
I(V) = T(VI2) + /(V; Z) I(VIA) = I(VIZ, A) + ((A
(8)
Alternatively. conditional (cross) information, such as T(X; A), can be expressed in terms of non-conditional (cross) entropies using Eq. 4
As recalled above, the cross-entropy between m variables, V {X1,".. ,,m,, J, is defined as,
I(X;Y A) = H(Y))))))(
(10)
= H(X, A) + HY, A) H(X, Y, A) ()H((A)
H(V) = - p(=1.-... -Epe1 ................................ 2Tm) logg pm) (3)
which can in turn be expressed in terms of non-conditional (cross) information as,
{a.}
where p( (xxx)) is the empirical joint probability distribution of the variables {X,} and g({x;}) the joint probability dis- tribution of the model. Bayes formula. g(xx,), J. {M₂}) = q(fz;) Hyrh) 9((N,)) }), directly translates into the definition of conditional cross-entropy through the decomposition,
I(X 4X) - (-1)* I(XY - ::A.))
--.' (X;Y;A;;:> Am)
X.YSS'
Σ (-1)1s' S'CS
(11)
(((x)), (G))) = H((XH)()) + H(())) (4)
Multivariate (cross) information, I(V)= = I(X;: ... X,m)
where S = (X,Y) This corresponds to all (cross)
are defined from multivariate (cross) entropies through Inclusion-Exchin formulas over all subsets of variables
information terms including both X and Y in the expression of the multivariate (cross) entropy, H(X, Y, A), Eq.
2

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
2.2. Graphs and Connection Criteria
(through any type of edge) there is an ac-connecting path between X and Y given C.
2.2.1. DIRECTED MIXED GRAPHS, ANCESTRAL GRAPHS
Two vertices are said to be adjacent if there an edge (of any type) between them, X*-xY, where stands for any (head or tail) end mark. X and Y are said to be neighbors if X Y, parent and child if X - Y and spouses if X + Y g.
2.3. I. Likelihood Decomposition of Ancestral Graphs
Theorem 1. [likelihood of ancestral graphs] The cross- entropy H(p g) and likelihood LDig of an ancestral graph G is decomposable in terms of multivariate cross- information, I(C), summed over all ac-connected subsets of variables, C (Deftnition 3),
A path in G is a sequence of distinct ver- tices Vi,..., Vn consecutively adjacent in g, as, In particular,
ac-connected Σ CCV
Vi*+V/2** -. - collider path between Vi and Vn has the form +1++ 1/2 +
(p.
(-1))C
V V and directed path corresponds to V1 - V/2 + -....V......
ac- connecteeee
D
exp
(N
(-1) (12)
ZD.
X is called an ancestor of Y and Y a descendant of X if X = Y or there is a directed path from X to Y. X + + Y. AY) denotes the set of ancestors of Y in g. By extension, for any subset of vertices, C CVV, (C) denotes the set of ancestors for all Y e c g.
CCV
where N is the number of üid samples in the dataset DD and ZD. a data- and dependdd normalization constant.
The proof of Theorem 1 is left to Appendix B. It is based on partition of the cross-entropy (Eq. (6) into cross- information contributions from ac-connected and non-ac- connected subsets of variables, which does not rely on head-and-tail factorizations nor on imset formalism Hu and Evans proposed an equivalent result (Proposition 3.3 in (Hu & Evans 2020 with a proof using head-and-tail decomposition to define parametrizing sets, which happen to coincide with the ac-connected sets defined here (Defini- tion 3). Theorem 1 characterizes in particular the Markov equivalence class of ancestral graphs (Richardson & Spirtes
A directed mixed graph (DMG) is a vertex-edge graph g = (V, E) that can contain two types of edges: directed (-)) and bidirected (+)) edges.
directed cycle occurs if X E (Y) and X Y. An almost directed cycle occurs if X E Ang and X + Y.
Definition 1. An ancestral graph is a DMG:
i) without directed cycles;
ii) without almost directed cycles.
An ancestral graph is said to be maximal if every miss- ing edge corresponds to a structural independence. If an ancestral graph G is not maximal, there exists a unique maximal ancestral graph 9 by adding bidirected edges to g (Richardson & Spirtes 2002
2002
Richardson 2003 Ali & Richardson 2002 Ali et al.
2005
Tian 2005
Ali et al. 2009) as,
Corollary 2. Two ancestral graphs are Markov equivalent iff they have the same ac-connected subsets of vertices.
Note, in particular, that Eq. 12 holds for maximal ancestral graphs (MAG), for which all pairs of ac-connected variables are connected by an edge, and their Markov equivalent repre- sentatives, the partial ancestral graphy (PAG) (Richardson)
2.2.2. ac- LC-CONNECTING PATHS, ac-CONNECTED SUBSETS
Let us now define ancestor collider connecting paths or ac-connecting paths, which entail simpler path connecting criterion than the traditional m-connecting criterion, dis- cussed in the Appendix A. Yet, ac-connecting paths and ac-connected subsets will turn out to be directly relevant to characterize the likelihood decomposition and Markov equivalent classes of ancestral graphs.
& Spirtes 2002 1999 Zhang 2007 2008
Proposition 3. The likelihood decomposition of ancestral graphs (Eq. [2] Theorem 1) can be estimated by replac- ing the model distribution q by the empirical distribution in the retained multivariate cross-information terms ((C) corresponding to all ac-connected subsets of variables, C.
Definition 2. [ac-connecting path] An ac-connecting path between X and Y given a subset of variables C' (possibly including X and Y) is a collider path, X * Zi >> + Zg+ +*Y; with all Z, Ang (YU)).CC.. that is, with Z, in C or connected to {X, Y) UC by ancestor path, i.e. Z: ..... + T with T E {X,Y)UC
Hence, Proposition amounts to estimating all relevant cross-information terms in the likelihood function with the corresponding multivariate information terms computed from the available data, while assuming by construction that the model distribution obeys all local and global condi-
Definition 3. [ac-connected subset] A subset C' is said to be ac-connected if VX, Y E C. and Y are connected
'The genesis of Theorem 1 and Proposition 3 is fur- ther discussed in https::/ 72 22f4 Bo (dt.:nnooooel@@GGG.. wtt.eewwwee.feeiiiii....
3

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
tional independences entailed by the ancestral graph. The corresponding factorization of the model distribution can be expressed in terms of empirical distribution, assuming positive distributions, see Appendix C.
Markov equivalent, as discussed in (Ali et al. 2009). In- deed, their cross-entropy decompositions differ by two ac- connected contributing terms: three-point cross informa- tion I with a collider path not confined in C (i.e. X ~ Z em T ++ Y and corresponding ancestor path Z --+ Y) and a four-point information term [(X; Y; T) due to the two-collider path (X + Z x- + T < Y). More quantitatively, it shows that the graph of Fig. [ with two-collider path is more likely than the graph of Fig. DF whenever (X; T) - ((; Y; 7; - [(X; 22))
Fig. D Iillustrates the cross-entropy decomposition for a few graphical models in terms of coss-inforrmtton contribu- tions from their ac-connected subsets of vertices. In par- ticular, an unshielded non-collider (e.g. X + Z + W. Fig. [T], is less likely (i.e. higher (osss-enroopy) than unshielded collider or 'v-structure' (e.g. X + Z + W. Fig. [BB, ). if the corresponding three-point information term is negative, [[;; Z: W) < 0, in agreement with earlier results (Affeldt & Isambert 2015 Verny et al. 2017). How- ever, this early approach, exploiting the sign and mag- nitude of three-point information to orient v-structures, does not include higher order terms involving multiple v- structures, which can lead to orientation conflicts between unshielded triples, in practice. Resolving such orientation conflicts requires to include information contributions from higher-order ac-connected subgraphs, such as star-like ac- connected subsets including three or more parents, Fig. [ Similarly, the cross-entropies of collider paths involving several colliders also include higher-order terms, as with the simple example of two-collider path, Fig. DP. By contrast, the cross-entropy based on the head-and-tail fac- torization of the same two-collider path, i.e. ((x, 2, M, w) = q(2, y/ær w) ())))((() (Richardson 2009). is found to be equivalent to the cross-entropy of Bayesian graph without bidirected edge, Fig. D when estimated with the empirical distribution p(-), see Appendix C. This observation illus- trates the difficulty to estimate the likelihood functions of ancestral graphs using head-and-tail factorization.
T(X; V1>) - T(X; YI 7) <0. Finally, the Markov equiv- alent graphs of Fig. D and Fig. D. also due to (Ali et al.) 2009), illustrate the fact that the actual ancestor collider path between unconnected pairs does not need to unique nor conserved between Markov equivalent graphs (as long as their cross-entropies share the same multivariate cross- information decomposition).
3. Efficient Search-and-Score Causal Discovery using Local Information Scores
The likelihood estimation of ancestral graphs (Theorem I and Proposition 3) enables the implementation of a search- and-score algorithm for this broad class of graphs, which has attracted a number of contributions recently (Triantafillou & Tsamardinos 2016: Rantanen et al., 202 Claassen & Bucur 2022] Anndews et al.) 2022 Hu & Evans, 2024b a))
Our specific objective is not to develop an exact method limited to simple graphical models with few nodes and small datasets but to implement an efficient and reliable heuristic method applicable to more challenging graphical models and large datasets.
Further examples of graphical models, Figs. -I,, show the relative simplicity of the decomposition with only few (non- trivial) ac-connected contributing subsets C with |C| > 3, as compared to the much larger number of non-ac-connected non-contributing subsets, that cancel each other by con- struction due to conditional independence constraints of the underlying model. Note, in particular, that most con- tributing multivariate information (C) only concern direct connections or collider paths within a single component subgraph induced by C' (solid line edges in Fig. D. How- ever, occasionally, collider paths extending beyond C into Ang(C) \ C (marked with wiggly edges) with correspond- ing ancestor path(s) (marked with dashed edges) do occur, as shown in Fig. DG.
Indeed, search-and- to rely on heuristic rather than exhaustive search, in gen- eral, given that the number of ancestral graphs grows super- exponentially the number of vertices increases. This can be implemented for instance with Monte Carlo algorith- mic scheme with random restarts, which efficiently probes relevant graphical models. Here, we opt, instead, to use the prediction of an efficient hybrid causal discovery method, MIIC (Verny et al. [017 Cabeli et al. [02T Ribeiro- -antas et al. 2024), starting point for a subsequent search-and- score approach based on the proposed likelihood estimation of ancestral graphs (Eq. 12 and Proposition 3).
structure learning methods need
Moreover, while the likelihood decomposition of ances- tral graphs may involve extended ac-connected subsets of variables, as illustrated in Fig. T we aim to implement a computationally efficient scarch-and-score causal discov- ery method based on approximate local scores limited to the close surrounding vertices of each node and edge. Yet, while MIIC only relies on unshielded triple scores, the novel search-and-score extension, MIIC_ search&score, uses also higher-order local information scores to compare alternative
In addition, the present information-theoreii decomposi- tion of the likelihood of ancestral graphs can readily distin- guish their Markov equivalence classes according to Corol- lary 2. For instance, the ancestral graphs of Fig. MF and Fig. DG. despite sharing the same edges and the same unshielded collider (X + Z + T), turn out not to be
4

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
accommend
H(p.g)
Σ
Mxg x))
Mx xp Triples
1[x x (2 4-oples
Σ
+
Veriies
Edges
Σ(x)
Σ
Verces
Edges
Σ x)))
Σ Mxxxx)
Vertios
Edgrs
Σ I[ x
Mxx x))
Veriees
Edges
Σ(x)
Σ Mxx x))
Vartioa
Edge
E x))
Σ (X
Verices
Edges
Σ x
Σ
+
Vertice
Edge
Σ ( +
Vertioes
Edgas
Σ
E E x)) +
Verniess
Edgrs
x())
Σ MKK
Verices
Edges
Figure Cross-entropy decomposition ancestral graphs. Examples cross-entropy decomposition ancestral graphs (red edges, Ihs) in terms of relevant multivariate cross-information contributions I(C) with C V (red nodes, rhs). Simple graphs: (A) without unshielded colliders, (B) single non-overlapping unshielded colliders. (C) with overlapping unshielded colliders through three more (conditionally) independent parents or (D) through two-(or more)-collider path. (E) Bayesian graph corresponding to the head-and-tail factorization the two-collider path in (D) estimated using the empirical distribution p(.), see Appendix (F) Simple Bayesian graph not Markov equivalent to ancestral graph (G) sharing the same edges and unshielded collider (Ali et al. 2009). Solid black edges correspond to direct connections or collider paths confined to the corresponding ac-connected subset C,, while wiggly edges indicate collider paths extending beyond C yet indirectly connected C by ancestor path, marked with dashed edges, Definition 2. By contrast, graphs and illustrate the fact that collider paths may not be unique nor conserved between two Markov equivalent graphs (i.e. sharing the same cross-information terms) (Ali al. 2009)
subgraphs, as detailed below.
3.1. MIIC, an Hybrid Causal Discovery Method based on Unshielded Triple Scores
The proposed method is shown to outperform MIIC and other state-of-the-art causal discovery methods on challeng- ing datasets including latent variables.
MICC an hybrid causal discovery method combining constraint-based and informatiin--hherrtic frameworks
5

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
VVeny et al. 2017 Cabeli et al.
2020).
Unlike traditional
to the global likelihood decomposition, Eq. 12 (Step 2). The two steps are illustrated on a running example, Fig. D.1
constraint-based methods (Pearl 2009 Spirtes al., 2000), MIIC does not directly attempt to uncover conditional inde- pendences but, instead, iteratively substracts the most sig- nificant three-point (conditional) information contributions of successive contributors, A1, A2. A.. from the mutual information between each pair of variables, I(X;Y).
3.2.1. STEP 1: NODE SCORES FOR EDGE ORIENTATION PRIMING AND EDGE REMOVAL
The first step consists in minimizing a node score corre- sponding to the local normalized log likelihood of cach node w.r.t. its possible parents or spouses amongst the con- nected nodes predicted by MIIC. To this end, the node score assesses the conditional entropy of each node w.r.t. se- lection of parents, spouses or neighbors, Pa'x, C Pax U Sp x U Ne and factorized Normalized Maximum Like- lihood (fNML) regularization (Affeldt & Isambert (015) (conditional entropies can be estimated as differences of mutual information for continuous variables, Appendix D),
T(X;Y) A() T(X;Y;; A1)
... [[X:Y:A,11 (A)] = (A.).) (13)
where Au hx) > 0 is the positive informa- tion contribution from Aş to I( (Affeldt & Isam- bert, 2015] Affeldt et at.) (016). Conditional indepen- dence is eventually established when the residual condi- tional mutual information on the right hand side of Eq. [3] T(X; (A)))) becomes smaller than complexity term, L.e (x,Y(AA > ((X;Y[AA) >> O, which depen- dents on the considered variables and sample size N.
1 9ng
Score,(X)) = H(X,|Paxx)) )+
Σ log logCin Cijy
(14)
N
This leads to an undirected skeleton, which MIIC then (par- tially) orients based on the sign and amplitude of the reg- ularized conditional 3-point information terms (Affeldt & Isambert 2015; Verny al. 2017). In particular, negative conditional 3-point information terms, / [X;; Y; (A.)) O, correspond to the signature of causality in observational data (Affeldt& & Isambert, 2015) and lead to the prediction of a v-structure, X + Z < Y. if X and Y are not connected in the skeleton. By contrast, a positive conditional 3- 3-point information term, (x; Y; ({A}) > 0, implies the absence of v-structure and suggests to propagate the orientation of a previously directed edge X + Z Y as X + Z + Y.
where qo: corresponds to the combination of levels of Pa' xa while Ta, is the number of levels of Xt, and ny the number of samples corresponding to a particular combination of levels j in each summand, with Σ₁ n = N, the total number of samples. log City is the INML regulatization cost summed over all combinations of levels, qu,, (Kontkanen &
Myllymäki 2007 Roos al. 2008). see Appendix D.
This first algorithm is looped over each node, priming the orientations of their surrounding edges (as directed, bidi- rected or undirected), until convergence. Edges without orientation priming at either extremity are assumed to be false positive edges and removed at the end of Step 1.
In practice, MIIC's 's strategy to circumvent spurious condi- tional independences significantly improves recall, that is, the fraction of correctly recovered edges, compared to tradi- tional constraint- methods. Yet, MIIC only relies on unshielded triple scores to reliably uncover significant con- tributors and orient v-structures, as outlined above. MIIC has been recently improved to ensure the consistency of the separating set in terms of indirect paths in the final skeleton or (partially) oriented graphs (Li et al. 2019) RRbeiio-Dantaa jet al. 2024) and to improve the reliably of predicted orien- tations (Cabeli et al. 2021 Riheiro-Dantas et al. 2024).
3.2.2. STEP 2: EDGE ORIENTATION SCORES, AS LOCAL CONTRIBUTIONS TO THE GLOBAL LIKELIHOOD
The second step consists in orienting the edges retained after Step 1, based on the optimization of their local contribu- tions to the global likelihood score, Eq. 12 restricted to ac-connected subsets containing up to two-collider paths. This amounts to minimizing each edge orientation score w.r.t. its nodes' parents and spouses, corresponding to minus the conditional information plus fNML complexity cost, Table 1 given three sets of parents and spouses of X and Y. i.c. Pa' xx = Pa, U Sp SY, \Y, Pa yux = US,USp XX and Pa xy = Pa xxy U Pa ru with their corresponding combinations of levels, 93s- any and q.v. These orien- tation scores, Table [ include symmetrized fNML com- plexity terms to enforce Markov equivalence, if X and Y share the same parents or spouses (excluding X and Y), see Appendix D. Indeed, all three scores become equals if Pa'nx = Pa x Pa' x implying also the same combi- nations of parent and spouse levels, 9pr 4.ng Lzv'
The predictions of this recent version of MIIC, which in- clude 3 types of edges (directed, bidirected and undirected), have been used as starting point for the subsequent local search-and-score method implemented in the present paper.
3.2. New Search-and-Score Method based on Higher-Orded Local Information Scores
Starting from the structure predicted by MIIC, as detailed above, MIIC_ See_sseeeacrrrreeeeeeeeeee proceeds in two steps, first to remove likely false positive edges (Step 1) and then to orient the remaining edges based on their estimated contributions
While orientation scores cannot be summed over individual edges due to multiple countings of ac-connected contri-
6

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
Table 1. Local scores for the orientation of a single directed or bidirected edge, see Appendix D. Symmetrized fNML. complexity (Markow equivalent)
Edge
Information
x I(X:Y|Pa's
Σ
log CR. - Σ, En' log Ch₃y Σ)
log Ciny
Mote log Ci3)
NN
X
"
log Crj - Σ) log Chüy
log Cuzy
"yve log Ciay.
2N
Y
Paxx))
Ray
log Cig - Σ)
New log Cngy + Σ yuz
log Ch
Byy log Ciss
butions, score differences between alternative orientations provide an estimate of the global score change. Hence, step 2 algorithm is looped over each edge to compute orientation score decrement, given its current orientation and the orientations of surrounding edges. The orientation change corresponding to the largest global score decrement, without forming new directed or almost directed cycles, is then chosen each iteration until convergence or until limit cycle is reached. Limit cycles may originate from the local two-collider approximation of the global score.
LINEAR Comtinouss Models
NON-LINEAR comtinuous Models
Pradisio
Rocall
Precsion
Recall
10% Lv
4. Experimental Results
orientation scores (Table effectively predicts bidirected orientations on three simple ancestral models, Fig. E. in Appendix E, when the end nodes do not share the same parents (Fig. E Model 1). share some parents (Fig. E.1 Model 2) or when the bidirected edge is part of a longer than two-collider paths (Fig. E. Model 3). The prediction of the edge orien- tation scores are summarized in Table E.I and show good predictions for large enough datasets.
We first tested whether
20 LV
1,000
10.000
1,000
10,000
100 1,000 10,000
100
1,000 10,000
We then tested MIIC search&score performance on ances- tral graphs obtained by hiding up to 20% of variables in either linear or non-linear continuous Bayesian networks of 50-150 nodes and average degree 3-5. Figs. 2 & E.2 compare MIIC_search&aee used as starting point for MIIC search&score and to M3HC (Triantafillou & Tsamardinos 2016), GFCI (Ogarrio al.] 2016) and DAG-GNN (Yu et al. 2019). a DL method based on continuous constrained optimization. The parameter settings for each method are detailed in Appendix E.
Sample Size (N)
M3HC
GFCI
DAG-GNN
MIIC
MICC search&score
Figure 2. Benchmark results on continuous datasets. Bench- mark results are averaged over 30 independent ancestral graph models obtained by hiding 0%, 10% or 20% of variables in linear Gaussian models (left) in more complex models with non-linear couplings between variables (right), and including 50 nodes av- crage degree 5. results used as starting point for MIIC search&score, M3HC (Tr
performance to MIIC results
results compared to MIIC
antafillou & Tsamardinos] 2016). GFCI (Ogarrio et al.
2016) and
DAG-GN (Yu et al. 2019). See also Fig. E.2 Causal discovery performance is assessed in terms of Precision and Recall relative the theoretical PAGs, while counting as false positive all correctly predicted edges but with a different orientation as the directed or bidirected edges of the PAG. Error bars: 95% confidence interval.
While MIIC_ SUCC_ssarc&&scoree is outperformed on linear Gaussian models by causal discovery methods assum- ing linear combinations of the variables (with possible pre- or post-nonlinear transformations for DAG-GNN), sulC_scerch&scoree clearly outperforms these methods on more complex models including non-linear couplings be- tween variables, Figs. 2 & F.2 MIIC_search&sccore also outperforms MIIC in terms of edge precision with little to no decrease in edge recall, demonstrating the benefit of MMICCseearhhalll search&score's rationale to improve MIIC predic- tions by extending MIIC information scores from unshielded triples higher-order information contributions.
We also analyzed challenging benchmarks for categorical datasets from the bnlearn repository (Scutari) 2010). They concern ancestral graphs obtained by hiding up to 20% of variables in Discrete Bayesian Networks of increasing com- plexity (number of nodes and parameters), such as Alarm (37 nodes, 509 parameters), Insurance (27 nodes, 1,008 pa- rameters), Barley (48 nodes, 84 links, 1 14.005 parameters). and Mildew (35 nodes, 540, 150 parameters), Figs.
7

---

An Efficient Scarch-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
ALARM (509 parameters)
INSURANCE (1,008 parameters)
BARLEY (114,005 parameters)
MILDEW (540,150 parameters)
Precision
Recall
Precision
Recall
Precision
Recall
Precision
Recall
0% LV
10 Lu
28 % LV
100 1,000 10,000
100
1,000
10,000
100
1,000
10,000
100
1,000
10,000
100
1,000
10,000
100
1,000
10,000
100
1,000
10,000
100
1,000
10,000
Sample Size (N)
FCI
MIIC
MIIC search&score
Figure 3. Benchmark results on 'real-world' categorical datasets from the bnlearn repository. Benchmark results are averaged over 50 independent categorical datascts from ancestral graphs obtained by hiding 0%. 10% or 20% of variables in Discrete Bayesian Networks of increasing complexity (see main text): Alarm, Insurance, Barley and Mildew. MIIC search&score results are compared, as in Fig, 2 to MIIC results used as starting point for MIIC search&score and to FCI Zheng et al.
2024). Error bars: 95% confidence interval.
compare MIIC search&score performance to MIIC results and to FCI from the python causal-learn package (Zheng et al. 2024). The parameter settings for each method are detailed in Appendix E. Fig. results obtained from independent datasets for each ancestral graph and sample size, while Fig. E.3 results provide bootstrap sensitivity analysis to sampling noise for each method based on independent resamplings with replaceement of single datasets of increasing size.
Figs.
Hence, MIIC search&score is shown to outperform all other tested causal discovery methods on complex non-Gaussian datasets including non-linear couplings between variables, while GFCI and DAG-GNN are the best performers on lin- ear Gaussian datasets. These results demonstrate GFCI, M3HC and DAG-GNN's clear advantage in assuming intrin- sically linear combinations of the variables, when analyz- ing multivariate Gaussian datasets. However, these results also highlight the clear limitation of this assumption, when analyzing more complex datasets including non-linear cou- plings between variables, which goes beyond the variable transformation of pre- post-nonlinear causal models. By contrast, MIIC seccssechh&tccrree and MIIC make no particular assumption on the data distributions and achieve similar good performance across a broad range of Gaussian or non- Gaussian multimodal distributions Figs. 2 E.2 as well as on complex
As with continuous data, MiIC_search&score outperforms MIIC in terms of edge precision with little to decrease in edge recall, Figs. 3 E.3 @IlC_search&sccore is also found to outperform FCI on both precision and recall on small datasets (e.g. N < 10, 000 samples) of complex graphical models (i.e. Barley and Mildew), while reaching similar performance at larger sample sizes or for simpler graphical models (i.e. similar precision on Alarm and Insurance), as expected from the asymptotic consistency of FCI for very large datasets, Fig. We also observed that FCI had a hard time to converge on bootstrapped datasets, explaining the lack of FCI comparison with MIIC and MIIC_ search&score in Fig. E.3 for complex models at large sample sizes.
categorical datasets Figs.
Importantly, the benchmark PAGs used to discovery methods with increasing proportions of latent vari- ables (Figs. 2 & 3 and Figs. E.2 & E.3) include not only bidirected edges originating from hidden common causes but also additional directed or undirected edges arising, in
the causal
8

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
particular, from indirect effects of hidden variables with observed parents. Irrespective of their orientations, all these additional edges originating from indirect effects of hidden variables generally correspond to weaker effects (i.e. lower mutual information of indirect effects due to the Data Pro- cessing Inequality) and are more difficult to uncover than the edges of the original DAG model without hidden vari- ables. This explains the steady decrease in recall for com- plex ancestral models with higher proportions of hidden variables, while precision remains essentially unaffected. This is clearly apparent for complex random models with average degree 5 (Figs, & E.2) and for complex real-world models, Insurance, Barley and Mildew (Figs. 3 E.
detailed in section 3. This already shows significant im- provement in causal discovery performance (i.e. combining good precision and good recall on challenging benchmarks) as compared to existing oftt--o-thhearle art methods.
In addition, for ancestral graphs including only ac- connected subsets with a maximum of two-collider paths, the local consistency of the likelihood score (Eq. 12) guar- antees that the two-step algorithm gives a correct estimate of the global likelihood, based on Proposition 3, just like with the likelihood decomposition of Bayesian networks (Eq, 2 Besides, for more complex ancestral graphs includ- ing longer collider paths, an oracle version of the approach could also be proposed, as the general likelihood decom- position (Eq. 12 could in principle be used in conjunction with an exhaustive search-and- algorithm over MAGs or PAGs, which can be generated rather efficiently (Hu & Evans 2020). Alternatively, and in practice for graphs with more than about 10-15 nodes, an MCMC algorithm based on the likelihood decomposition (Eq. [2] could be used to efficiently search for high-scoring MAGs or PAGs.
MIIC_search&&scor scalability is primarily limited by the quadratic complexity of MIIC w.r.t. the number of nodes, with only small time increase when including latent vari- ables thanks to MIIC's greedy approach, see Fig. S5 in (Verny et al. 2017). By contrast, the two-step search-and- score scheme of MelC_seaacc&&scorreee is essentially linear in the numbers of nodes (step 1) and edges (step 2) for a fixed degree. In practice, Steps 1 and 2 take similar running time as for MIIC to provide a starting graph for MIlC_searc&&score. Hence, MMIC__serrc&&er ity remains competitive against recent differential causal dis- covery methods using continuous optimization techniques, whose complexity is generally cubic in the number of vari- ables, although more scalable methods have recently been developped, e.g. (Lopez et al. 2022 Montagna et al. 2023 Amin& & Wilson 2024).
Finally, another limitation of the proposed likelihood de- composition (Eq. [2) is its restriction to MAGs without undirected edge, as shown in the proof of Theorem 1, Ap- pendix B section (iii). However, it would be interesting to see how this result could be extended to MAGs including undirected edges. In addition, there are other existing scores under latent variables for specific parametric assumptions like discrete variables, exponential or stratified Gaussian families. In particular, discrete chain graph models (Drton 2009, fully bidirected graph models (Drton & Richard- son 2008), and discrete nested Markov models (Richardson et al. 2023) have been shown to be curved exponential models, which can be scored consistently using BIC scores. One extension of these BIC scores has also been recently proposed (Bellot et al., 2024) to allow for the distinction between different acyclic directed mixed graphs (ADMG) from (Verma & Pearl 1990), see Fig. 1c & d in (Bellot et al.) 2024). These ADMGs imply the same set of conditional in- dependence constraints and yet are distinguishable because they imply an equality between different functionals of the probability distribution (see Eq. in (Bellot et al., (0024). It would be interesting to further explore the relation of these other existing scores under latent variables with the present information theoretical score (Eq. 12).
scalabil-
ca- pacity to efficiently and robustly learn complex graphi- cal models from limited available data, which is a fre- quent setting for many real-world applications. In addition, which has been implemented to ana- lyze challenging categorical or complex continuous datasets including non-lineat couplings between variables, is quite unique in this regard, all other search-and-score methods for ancestral graphs (Triantafillou & Tsamardinos 2016
All in all, these results highlight
Rantanen et al. 2021 Claassen & Bucur [022]] Andrews ct al. 2022; Hu & Evans 2024a) have only been demon- strated on continuous datasets from linear Gaussian models.
5. Limitationns
The main limitation of the paper concerns the local scores used in the search-and-score algorithm, which are limited to ac-connected subsets of vertices with a maximum of two-collider paths.
A Acknovledgements
We would like to thank MICC Team members for their con- tributions to the MIIC projects and the reviewers for their constructive remarks and suggestions. NL acknowledges funding from the CNRS-Imperial College London joint PhD programme. HI acknowledges funding from CNRS, Institut Curie, ANR PEPR A14scMed (ANR- --N-22--R 22---ES---002) and
While this approach could be extended to higher-order infor- mation contributions including three-or-mont cor-mor--colldee paths, it allows for simple two-step search-and-scom scheme at the level of individual nodes (step 1) and edges (step 2), as

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
affiliation to Prairie Institute.
Chickering, D. M. A Transformational Characterization of Equiv- alent Bayesian Network Structures. In UAI '95: Proceedings of the Eleventh Annual Conference on Uncertainty in Artificial Intelligence, pp. 87-98. Morgan Kaufmann, 1995.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Claassen, T. and Bucur, 1. G. Greedy equivalence search in the presence of latent confounders. In Cussens, J. and Zhang, K. (cds.), Proceedings of the Thirty-Eighth Conference Uncer- tainty in Artificial Intelligence, volume 180 of Proceedings of Machine Learning Research, pp. 443-454422 PMLR, 01-05 Aug 2022.
References
Drton, M. Discrete chain graph models. Bernoulli, 15(3):736-733 aug 2009, doi: 10.31550088-b bej17.
Affeldt, S. and Isambert, H. Robust reconstruction of causal graphi- cal models based on conditional 2- 2-point and 3-point information. In Proceedings of the Thirty-First Conference Uncertainty in Artificial Intelligence, UAI 2015, July 12-16, 2015, Amsterdam, The Netherlands, pp. 42-51, 2015.
Drton, M. and Richardson, T. S. Graphical methods for efficient likelihood inference in gaussian covariance models. J. Mach. Learn. Res., 9:893-914, 2008.
Drton, M., Eichler, M., and Richardson, T. S. Computing maxi- mum likelihood estimates recursive linear models with cor- related errors. Journal of Machine Learning Research, 10(81): 2329-2348, 2009.
Affeldt, S., Verny. L., and Isambert, Н. 3of12: network re- construction algorithm based on 2-point and 3-point informa- tion statistics. BMC Bioinformatics, 17(S2). jan 2016. doi: 11/8///1189901---118
Evans, R. and Richardson, T. S. Maximum likelihood fitting of acyclic directed mixed graphs to binary data. In Proceedings of the 26th Conference on Uncertainty in Artificial Intelligence, UAT 10, Corvallis, OR, USA, 2010. AUAT Press.
Ali, R. A. and Richardson, T. S. Markov equivalence classes for maximal ancestral graphs. In Proceedings of the Eighteenth Conference on Uncertainry in Artificial Intelligence, 0AT'02, pp. 1-9, San Francisco, CA, USA, 2002. Morgan Kaufmann Publishers Inc. ISBN 1-55860-897-4.
Han, T. S. Multiple mutual informations and multiple interactions in frequency data. Information and Control, 46( 1):26-45, 1980.
Ali, R, A., Richardson, T. S., Spirtes, P., and Zhang. J. Towards characterizing Markow equivalence classes for directed acyelic graphs with latent variables. In Proceedings of the Fifteenth Conference Uncertainty in Artificial Intelligence. UAT'S,, San Francisco, CA, USA, 2005. Morgan Kaufmann Publishers Inc.
Hu, Z. and Evans, R. Faster algorithms for Markov equivalence. In Peters, J. and Sontag, D. (eds.), Proceedings of the 36th Confer- ence on Uncertainty in Artificial Intelligence (UAI), volume 124 of Proceedings of Machine Learning Research, pp. 739-748. PMLR, 03 06 Aug 2020.
Ali, R. A.. Richardson, T. S., and Spirtes, P. Markov equivalence for ancestral graphs. Ann. Statist., 37(5B):2808-2883 oct 2009. doi: 10.1214/08- 12214/00--00666.
Hu, Z. and Evans, R. A fast score-based search algorithm for maximal ancestral graphs using entropy. Preprint, 2024a. arxiv 2402.04777.
Amin, A. and Wilson, A. Scalable and flexible causal discovery with an efficient test for adjacency. In Proceedings of the 41 st International Conference Machine Learning, volume 235 of Proceedings of Machine Learning Research. PMLR, 2024. doi: 448.5500/rxx..2466.009
Hu. Z. and Evans, R. Towards standard imsets for maximal ances- tral graphs. Bernoulli, 30(3), 2024b. doi: 10.3 bej223bbee6633...
Kalisch, M., Mächler, M., Colombo, D., Maathuis, M. H., and
Bühlmann, P. Causal inference using graphical models with the package pealg. J. Stat. Softw., 47(11) 26, 2012.
Andrews, B., Cooper, G. F., Richardson, T. S., and Spirtes, The m-connecting imset and factorization for ADMG models. Preprint, 2022. arxiv 2207.08963.
Koller, D. and Friedman, N. Probabilistic Graphical Models: Principles and Techniques. MIT Press, 2009.
Kontkanen, P. Computationally Efficient Methods for MDL-
Bellot, A., Zhang, J., and Bareinboim, F. Scores for learning discrete causal graphs with unobserved confounders. Proceed- ings of the AAAI Conference on Artificial Intelligence, 38(10): 11043-11051, 2024. doi: 10.1609/aaai.v wtt.aaaaaw.i.i.....2.. 38.10.28980.
Optimal Density Estimation and Data Clustering. PhD thesis, 2009.
Kontkanen, P. and Myllymäki, P. lincar-time algorithm for computing the multinomial stochastic complexity. Inf. Process. Lett., 103(6) 2007.
Cabeli, V. Verny L.., Sella, N., Uguzzoni, G., Verny, M., and Isambert, H. Learning clinical networks from medical records based on information estimates in mixed-type data. PLoS Com- put. Biol., 16(5))e1007866. 2020. doi: Vo.iral 1007866.
Kontkanen, P., Buntine, W., Myllymäki, P., Rissanen, J., and Tirri, H. Efficient computation of stochastic complexity. in: C. Bishop, B. Frey (Eds.) Proceedings of the Ninth International Conference Artificial Intelligence and Statistics, Society for Artificial Intelligence and Statistics, 103:233-223 22-.. 2003.
Cabeli, V., Li, H., da Câmara Sibeiro-Dantas M., Simon, F. and Isambert, H. Reliable causal discovery based on mutual information supremum principle for finite datasets. In WHY21, 35rd Conference on Neural Information Processing Systems. NeurIPS, 2021.
Lagrange, N. and Isambert, H. An efficient search-and-score algorithm for ancestral graphs using multivariate information scores. Preprint, 2024. arxiv 1752.
10

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
Li, H., Cabeli, V., Sella, N., and Isambert, H. Constraint-based causal structure learning with consistent separating sets, Ad- vances in Neural Information Processing Systems (NeurIPS), 32, 2019.
Richardson, T. S., Evans, R. J., Robins, J. M., and Shpitser, I. Nested Markov properties for acyelic directed mixed graphs. The Annals of Statistics, 51(1), 2023. doi: (0.1214/22-a05222
Rissanen. J. and Tabus. I. Kolmogorov's structure function in mdl theory and lossy data compression. In Adv. Min. Descrip. Length Theory Appl., pp. 245 264. MIT Press, 2005.
Lopez, R., Huetter, J.-C., Pritchard, J., and Regev, A. Large-scale differentiable causal discovery of factor graphs. Koyejo, S., Mohamed, S., Agarwal, A., Belgrave. D., Cho, K., and Oh, A. (eds.). Advances in Neural Information Processing Systems, volume 35, pp. 19290-19300 Curran Associates, Inc., 2022.
Roos, T., Silander, T., Kontkanen, P., and Myllymäki, P. Bayesian network structure learning using factorized nml universal mod- els. In Proc. 2008 Information Theory and Applications Work- shop (ITA-2008). IEEE Press, 2008.
McGill, W. J. Multivariate information transmission. Trans. of the IRE Prufessional Group on Information Theory (TIT). 4: 93-111, 1954.
Scutari, M. Learning Bayesian Networks with the bnlearn R Package. J. Stat. Softw., 35(3):1-22, 2010.
Montagna, F., Noceti, N., Rosasco, L.. Zhang, K., and Locatello, F. Scalable causal discovery with score matching. In der
Shtarkov, Y. M. Universal sequential coding of single messages. Problems of Information Transmission, 23(3):3-17, 1987.
Schaar, M.. Zhang, C., and Janzing, D. (eds.). Proceedings of the Second Conference on Causal Learning and Reasoning, volume 213 of Proceedings of Machine Learning Research, pp- 752-771. PMLR, 2023.
Spirtes, P., Glymour, C., and Scheines, R. Causation, Prediction, and Search. MIT press, 2nd edition, 2000.
Szpankowski, W. Average case analysis of algorithms on se- quences. John Wiley & Sons, 2001.
Ogarrio, J. M.. Spirtes, P., and Ramsey, J. A hybrid causal search algorithm for latent variable models. JMLR Workshop Conf. Proc., 52:368 379, 2016.
Tian, Generating Markov equivalent maximal ancestral graphs by single edge replacement. In Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence, UAI'05, San Francisco, CA, USA, 2005. Morgan Kaufmann Publishers Inc.
Pearl, J. Probabilistic reasoning in intelligent systems. Morgan Kaufmann, San Mateo, CA, 1988.
Pearl, J. Causality: models, reasoning and inference. Cambridge University Press, 2nd edition, 2009.
Tian, J. and Pearl, J. A general identification condition for causal effects. In Proceedings of the National Conference on Artificial Intelligence. pp. 567-573. Menlo Park, CA: Cambridge, MA; London; AAAI Press; MIT Press; 1999, 2002.
Pearl, J. and Paz, A. Graphoids:
A graph-based logic for reasoning about relevance relations, or when would x tell you more about if you already know z. Technical report, UCLA Computer Science Department, December 1985.
Ting, H. K. On the amount of information. Theory Probab. Appl., 7(4):439 447, 1962. doi: 10.1137/11070/
Triantafillou, S. and Tsamardinos, I. Score-based vs constraint- based causal learning in the presence of confounders. In CFA @ UAI, 2016.
Rantanen, K., Hyttimen, A., and Järvisalo, M. Maximal ancestral graph structure learning via exact search. In de Campos, C. and Maathuis. M. H. (eds.). Proceedings of the Thirty-Seventh Conference Uncertainty in Artificial Intelligence, volume 161 of Proceedings of Machine Learning Research, pp. 1237-1247. PMLR, 27-30 Jul 2021.
Verma, T. and Pearl, J. Causal networks: Semantics and expres- siveness. Shachter, R. D., Levitt, T. S., Kanal, N., and Lemmer, J. (eds.), Uncertainty in Artificial Intelligence, vol- ume 9 of Machine Intelligence and Pattern Recognition, pp. 69- 76. North-Holland, 1990. doi: 50011-1.
Ribeiro-Dantas,
Dantas, M. d. C., Li, H., Cabeli, V., Dupuis, L., Simon, F., Hettal, I.,, Hamy, A.-S., and Isambert, H. Learning interpretable causal networks from very large datasets, application 400, 000 medical records of hreast cancer patients. iScience, 27(5): 109736, 2024. doi: 111.1///jiici222411199990099
Verny, L.,, Sella, N., Affeldt, S., Singh, P. P., and Isambert, Learning causal networks with latent variables from multivariate information in genomic data. PLaS Comput. Biol., 13(10): e1005662, 2017.
Richardson, T. Markov properties for acyclic directed mixed graphs. Scandinavian Journal of Statistics, 3051):145-157 2003.
Yeung, R. W. A new outlook on shannon's 's information measures. IEEE transactions on information theory, 37(3):466-474, 1991. doi: 10. 10.11000 109/1199.77711 8.9902..
Richardson, T. and Spirtes, P. Scoring ancestral graph models. Technical report, 1999. Available as Technical Report CMU- PHII. 98.
Yo, Y., Chen, J., Gao, T., and Yu, M. DAG-GNN: DAG structure lcarning with graph neural networks. In Chaudhuri, K. and Salakhutdinov, R. (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Rexearch, pp. 7154-7163. PMLR, 2019.
Richardson, T. and Spirtes, P. Ancestral graph Markov models. Ann. Statist., 30(4):92--1030 aug 2002. doi: 10.1214/aos 1031689015.
Zhang. J. A characterization of Markov equivalence classes for directed acyclic graphs with latent variables. In Proceedings of the Seventeenth Conference on Uncertainty in Artificial In- telligence, UAI'07, San Francisco, CA, USA, 2007. Morgan Kaufmann Publishers Inc.
Richardson, T. S. A factorization criterion for acyclic directed mixed graphs. In Proceedings of the Twenty-Fifth Conference on Uncertointy in Artificial Intelligence, UAI '09, pp. 462-470, Arlington, VA, USA, 2009. AUAI Press.

---

An Efficient Scarch-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
Zhang, On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias. Artif. Intell., 172(16-1 1-661177777---1111111111-1---31- 1896, 2008.
Zheng. Y., Huang, B., Chen, W., Ramsey, J., Gong, M.. Cai, R., Shimizu, S., Spirtes, P., and Zhang, K. Causal-learn: Causal discovery in python. Journal of Machine Learning Research, 2- 2024.
12

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
Appendices
A. Preliminaries: connection and separation criteria
A.1. m-connection vs m'-connection criteria
An ancestral graph can be interpreted as encoding a set of conditional indepencence relations by a graphical criterion, called nn-separation, based on the concept of m-connecting paths, which generalizes the separation criteria of Markov and Bayesian networks to ancestral graphs.
Definition 4. [m-connecting path] A path between X and Y is m-connecting given a (possibly empty) subset C' S V (with X, Y & C) if:
i) its non-collider(s) are not in C, and ii) its collider(s) are in ACg(C).
Definition 5. [m-separation criterion] The subsets A and B are said to be m-separated by C, noted A A. BIC, there no m-connecting path between any vertex in A and any vertex in B given C.
The probabilisticc interpretation of ancestral graph is given by its global and pairwise Markov properties (which are equivalent (Richardson & Spirtes) ))02)): ): if A and B are m-separated by C, then A and B are conditionally independent given C and VX E A and AY E B, there is a probability distribution P' faithful to G such that their conditional mutual information vanishes, i.e. (X; Y|C) = 0, also noted X 1Pp Y/C.
However, as discussed above, the proof of Theorem will require to introduce weaker m mc-conncction criterion defined below.
Definition 6. connccnnnt path] A path between X and Y is m mlcconnecting given subset C C possibly in C) if:
V
(with X,Y
i) its non-collider(s) are not in C, and
ii) its collider(s) are in Ang (Mx((X,)) YJ UC).
Note, in particular, that an m-connecting path is necessary an m/ coonnectiing path but that the converse is not always true. For example, the path X - in Fig. [ (with Z + Y) an m'-connecting path given T' (as Z An, YY but not an m-connecting path given T (as Z (T)(
However, Richardson and Spirtes 2002 have shown the following lemma,
Lemma 4. [Corollary 3.15 in RRchardsson & Spirtes 2002 In an ancestral graph G, there is m² mnn-onnecging path between X and Y given C if and only if there a (possibly different) m-connecting path π between X and Y given C.
Hence, Lemma 4 implies that m' me-aaaaatiionn and m-separation criteria are fact equivalent, as an absence of m' connecting paths implies absence of m-connecting paths and vice versa. This enables to reformulate the m-separation criterion above as,
Definition 7. im'-separation (and m-separation) criteria] The subsets A and B are said to be m mn-separated (or m-separated) by C, if all paths from any X E A to any Y B have either
i) a non-collider in C, or ii a collider not in (Mn((XXY)UC)
The probabilistic interpretation of an ancestral graph given by its (global) Markov property: and B are m-separated (or m -sppaaaat) by C, then A and B are conditionally independent given C, noted as, A 1.m. B|C.
A.2. ac-connecting paths and ac-connected subsets
Let us now recall the definition of ancestor collider connecting paths or ac-comnecting paths, which is directly relevant characterize the likelihood decomposition and Markov equivalent classes of ancestral graphs (Theorem 1). We give here different yet equivalent definition of ac-connecting paths as defined in the main text (Definition 2) in order to underline the similarities and differencies with the notion of m' -connecting path (Definition 6).
Definition 8. [ac-connecting path] path т between and is an ac-connecting path given a subset C C VV ((with X and Y possibly in C) if:
13

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
does not have any noncollider, and ii) its collider(s) are in
U)..
Hence, more simply (following Definition 2 in the main text), an ac-connecting path given C' is a collider path, X +++ Zi ++ ... + ZK- f= Y, with all Z; E Ang((X,X) U) i.e. with Z; in C or connected to (x, U C' by an ancestor path, Z; + ... + T with E (XYYJUG
Definition 9. [ac-separation criterion] The subsets A and B are said to be ac-separated by C' if there is no uc-connecting path between any vertex in A and any vertex in B given C.
Previous definitions and Lemma 4 readily lead to the following corollary between the different connection and separation criteria:
Corollary 5.
i) m-connecting path π = m' m' connecting path x ū) ac-connecting path π = m' -connecting path π
tii)
m-separation + m' -separation
iv) m/m' m/seaa-papaaaiion
ac-separation
Finally, we recall the notion of ac-connected subset (Definition 3 in the main text), which is central for the decomposition of the likelihood of ancestral graphs (Theorem I): A subset C' is said to be ac-connected if A Y E C, there is an ac-connecting path between X and Y w.r.t. C.
B. Proof of Theorem 1.
In order to prove that the likelihood function of an ancestral graph, Eq. 12 contains all and only the ac-connected subsets of vertices in G (Definition 3), we will first show (i) that all non-ac-connected subsets S' are included in cancelling combination of multivariate cross tross-information terms, T(X; YIA) = 0, with X,Y E S' and S' S = {X, Y) U A, Eq. [1 Conversely, we will then show (ii) that cancelling combinations of multivariate cross information terms associated to pairwise YES'
Es'cs (-1)1s') ( O, do not contain any ac-connected subset S Finally, will prove (ii) that the information terms which appear in multiple cancelling combinations from different pairwise independence constraints do not modify the multivariate information decomposition of the likelihood function of ancestral graphs, Eq. 12 as these shared/overlapping terms in fact all cancel through more global Markov independence relationships involving higher order (three or more points) vanishing multivariate information terms, such as [(X; Y; 2 A) = 0.
conditional independence, (X;YlA)
i) Let's first prove that all non-ac-commected subsets S' are included in at least one cancelling combination of multivariate cross-information, [X;; Y| = 0, with X, S' and S' [x, A.
If S' is non-ac-connected subset, there is at least one disconnected pair X and Y for which each path æj between X and Y contains either some collider(s) not in Ans(S") or, if all colliders along πT are in Ang(S') there must be some non-collider(s) at node(s) Z, but not necessarily in S'. Let's define S = U, Zj. X and Y can be shown to be m-separated given S \ {X, Y}, as for each path π. between X and Y. its non-collider(s) are in S at node(s) Z₁ (when all collider(s) along æ are in S') or there is some collider(s) not in Ang(S') which are not in Ang either. The latter statement is proven by contradiction assuming that there is a collider at Z & Ang(S') such that E Ang(S). There therefore a directed path Z + + W with W E S. Hence, W E S' or there noncollider at W E Z, which is on path nT between X and Y along which all colliders + T where T E S' and thus Z Ang(S'), which is contradiction. Hence, all non-ac-
in Ang(S") by construction of S. This leads by induction
to Z > W connected subsets S' are included in a cancelling combination of multivariate cross-information terms, with X, E S' and S' C = {X,Y) U A.
0.
ti) Conversely, we will now show that cancelling combinations of multivariate cross- tross-information terms associated to pairwise conditional independence, (IA) Es'ss (-1)15') 0, do not contain any ac-connected subset S', where
X.YEs
S = =x A.
We will prove it by contradiction assuming that there exists a subset W C A, such that S' = (X,Y)U W is ac- connected. In particular, there should be an ac-connecting path between X and Y confined to Ang(S") and thus to ASg 2 AS"(S"), which is an mn -conneccting path between X and given A, contradicting the above hypothesis of m' -separation given A, i.e. I(X; Y|A) = The use of m' mssseaara,ion, f.e. the absence of m mn-coonecting paths with
14

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
colliders in Ang (S) rather than mn-connecting paths with colliders in Ang(A), necessary here, see Definitions and 6. Hence, ac-connected subset S' is included in cancelling combinations of multivariate cross-information terms associated to pairwise conditional independence, [Xx; Y|A) = Es'ss X.YY'' (-1)|5" = 0.
iii) Finally, we will show that the information terms which appear in multiple cancelling combinations from different pairwise independence constraints do not modify the multivariate cross-information decomposition of the likelihood function of ancestral graphs, Eq. 12 these shared/overlapping terms in fact all cancel through more global Markov independence relationships involving higher order (three more points) vanishing multivariate cross-information terms, such as [X;; Y; AI) 0.
This result requires to use an ordering of the nodes, Xx > Xj > X, that is compatible with the directed edges of the ancestral graph assumed to have no undirected edges, i.e. X, & An(X,) if X; > X₁. Under this ordering, higher order nodes X > X, X; can be a priori excluded from all separating sets Aty of pairs of lower order nodes, i.e. if T(X;; (Ay)) = 0 then X $ Auj-
In particular, the two pairwise conditional independence relations [(Xx; Xe(AAe) O. with Xe > Xk. and X( AAy)) O, with X, > X4, do not share any multivariate cross cross-information terms, if X2 * Xy. Indeed, as I(X;;X, AAy) contains all cross cooss-nfformation terms including both X, and X, as well as every subset (possibly empty) of A.j. none of them includes X, if X₃ > Xy. Therefore [(X;; X,) Aiy) does not contain any cross-information term of I(xx; Xe|Asr) which contains both X₁ and X₂ as well as every subset (possibly empty) of Age. This property eliminates all multiple countings multivariate cross- poossinfoormttion terms if X # Xy. Note that this result does not hold in general for ancestral graphs including undirected edges.
Hence, the issue of redundant multivariate cross-information terms in the likelihood decomposition, Eq. 12 is related to the conditional independences of two or more pairs, {X;, Xr). However, this situation also entails a more global Markov independence
(x,, X, sharing the same higher order node,
X,, i.e., I (Xx; X, AAr) = 0 for k i.j, constraint between X,, and (X,, x.... .x). given a separating set A, with A&r S AX,, ,x] for k = i.j. ... Such a global Markov independence constraint can be decomposed into more local independence constraints using the chain rule (in any order of the variables Ki..... (x)) and the decomposition rules of multivariate (cross) information (Eq. (
0 X(fX,; X,IA),
|A, + [XX X, +.... -IX66XAAAA)
+ -1
1I + I(X;;;,AA,"..
I(X; XA: X.IA) + I(X;; Xx;
where all the conditional multivariate cross-information terms vanish by induction due to the non-negativity of (conditional) mutual (cross) information. In particular, the conditional multivariate inosminfforattin terms in the last expression, Le. between X, and each subset of (X,, X..... ..)) given the separating set A, all vanish. This result can be readily extended to any subsets (X,, T,, Xg, X2) (conditionally) independent
Xx,,,:". Xy.... Xx) given separating set A, i.e. T([X,, "j. Xe); {X,, X,, (X)) = Hence, as the final conditional multivariate cross-information terms of the decomposition all vanish while not sharing any subsets of variables, it proves the absence of redundancy and global cancellation of con-aceconnected subsets (from pairwise and higher order conditional independence relations) in the likelihood function of ancestral graphs without undirected edges, Eq. 12
Hence, only ac-connected subsets effectively contribute to the cross-entropy of an ancestral graph with only directed and bidirected edges, Eq. 12
15

---

An Efficient Scarch-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
C. Factorization of the probability distribution of ancestral graphs
C.1. Factorization resulting from Theorem and Proposition
Before presenting the factorization of the model distribution ancestral graphs resulting from Theorem and Proposition 3, instructive to obtain an equivalent factorization for Bayesian graphs, assuming a positive empirical distributions, p(11,"" Fm) I -1-1 "1 .11 >
(011,-. Im)
Il (
(a pa
p pa..)
p pa.....
p(z1
Lm
p(
x11
i=1
m p(01," [m
p(
, )p( 1p(00 pa,, pa.,
(15)
p( 2i
1 ppa, ppa.)
This leads to the following alternative expressions for the cross-entropy H(p. = - Eg p(x)) log ((x) in terms of multivariate entropy and information, which only depend on the empirical joint distribution p(r))
m
H(p.9)
Σ
H( (ax.)
= H(XI, Xm) + Σ T(X; (ax.)
(16)
where Σ.= T(X;; X,Ki (Pax,) can be decomposed, using the chain rule and Eq. into unconditional multivariate information terms, which exactly cancel all the multivariate information of the non-ac-connected subsets of variables in the multivariate entropy decomposition, Eq. 6
Note, however, that this result obtained for Bayesian networks requires an explicit factorization of the global model distribution, q(xr,, in terms of the empirical distribution, p(x), which is not known and presumably does not exist, in general, for ancestral graphs.
Alternatively, assuming that the empirical and model distributions positive (Vx, p(x) 0, q(x)) 0), it is always possible to factorize them into factors associated to each (cross) information term in the (cross) entropy decomposition, Eq. 6 as,
m
T
2g)) (x,))(0;))
x II (4), ((2))))))))))))))))))))))) ((1)))(( <<jEk
((x)
g(x)) II
(17)
((1,,1)) (2;, (
where all the marginal distributions over subset of variables, egg. Zj, Emprs q(z) p(z,, Eepsyee p(x), cancel two-by-two by construction.
This can be illustrated on a simple example of a two-collider path including one bidirected edge, X + 2 + Y W (Fig. ID), valid for q()) and p(.) alike.
q( 2, w) = ((x) g(z) a(v) q(w)
g(z, 2) ((2,)) a(y,w) ((.1) g( g(2,w)
a(x) ((z) q(z)g(v) ((v) (()) (()) g(z) (w))
g(x) g(z) g(,) 2,y) g())
q(a, 2) ((2,)))(22))) g(x) g(z) q(w) ((,
((2,v) ((2 4))
(())
q(z, 2) q(1, ((2 q(z,v) ())
q(z, 2) ((2z2)(((2)))) q(2.2) g(y. w) q(x.)) ((z)) q 2, v) (xx, z, w) g (xx, Us w)
x
(18)
16

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
where all individual distribution marginals on subsets of variables, e.g. g(x)), ((x, =), ((2, 2, (or p p(2, pla, 2, v)). cancel two-by-two by construction, except ((2, 2, y, w) (or p(æ, 2, w)).
In addition and only for the model distribution q(.), all ratios in gray in Eq. 18 also cancel due to Markov independence relations across ceneec--enneted subsets (see proof of Theorem 1). This leaves truncated factorization retaining all and only the ac-connected subsets of variables in the graph, which we propose to estimate on empirical data by substituting the remaining g(.) terms by their empirical counterparts p(.), see Proposition 3.
This leads the following global factorization for g(.) in terms of p(.).
( 2, w) p(z) p(z) p(y) p(w) p(x, 2)
p(z,)) ply. w) p(z) p(z) p(z) p(y) p(2) plz, p(z) p(w) p(z,,, plz.))p(z.)) p(z.))( w(y
2,) p( p(z))
plz,2,y,w) p(x) p()) p(=)
p(a.y) w( p(2,w)
x plz, pla w) w(2,w) p(a.u) p( w) p(y. p(z) (() p(z, p(z) plu) p p (zz, U, w)
(19)
where the terms in gray have been passed to the Ihs of Eq. 18 applied to p(.). This ultimately leads to the analog of the Bayesian Network factorization in Eq. [5 but for the two-collider path, X - 7 < Y W (Fig. ID),
p(z) p(w) Mz(z))((wzz) p(x|w) p(ylw) p(z, w) (2,w|/))
[0,,,,V, w) p,
(20)
p(r, y/w)
where the last three factors "correct" the expression of p(a,, z, y, w) for the three (conditional) independences entailed by the underlying graph, that is, X L.W,W Z 1 W|X, and X 1 YIW.
C.2. Relation to the head-and-tail factorizations
The head-and-tail factorizations of the model distribution of an acyclic directed mixed graph, introduced by Richardson 2009 do not correspond to a single factorized equation (as with Bayesian graphs, Eq. [5) but to multiple factorized equations, which enable the parametrization of the joint probability distribution with independent parameters for ancestrally closed subsets of vertices.
For instance, the head-and-tail factorizations of the simple two-collider path including one bidirected edge, X - + Y W, introduced above, Fig. DD. correspond to the following equations (Richardson) 2009).
(x, w) q(x) q(w)
q(,, 2) <(2/2) ((x)
q(y. w) = <(vlw) q(w)
(I, w) = ((2|z) g(x) ((w)
( y. w) = <(vlw) q(w) q(x)
q (z, 2, w) 4(2, w)) ((x) q(w)
(21)
Importantly, these head-and-tail factorizations imply additional relations such as g(ylw) = g(y|,, w) (iee X .L YIW) obtained by comparing the last two relations in Eqs. [a]teer after marginalizing q(,, 2, y, w) over z. However, such implicit conditional independence relations are not verified by the empirical distribution p(.) in general and prevent the estimation of the head-and-tail factorizations by substituting the rhs q(.) terms in Eqs. 21 with their empirical counterparts p(-), as in the case of Bayesian networks, Fq. 15
Indeed, while the head-and-tal aail factorization relations, Eqs. 21 obey the local and global Markov independence relations entailed by the graphical model, Fig. [D. leading to the cancellation of all factors associated to non-ac-commected subsets in gray in Eq. 18 the remaining head-and-tail factors cannot be readily estimated with the empirical distribution p(-).
17

---

An Efficient Scarch-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
In particular, the cross-entropy of the two-collider path of interest, Fig. D. obtained with the head-and-tail factorizations corresponds to (p; q) Σ p(2,, 2, y, w) log log v/2, w)) ((x)) g(w). Then, estimating the g(-) terms with their p(.) counterparts leads to the cross-entropy of a Bayesian graph, Fig. E with a different Markov equivalent class than the ancestral graph of interest, Fig. [D. A similar discrepancy is obtained with a c-component factorization which leads to the cross-entropy of the Bayesian graph of Fig. T without edge + corresponding to a different Markov equivalence class than the previous two graphs, Figs. TP & E.
These examples illustrate the difficulty to exploit the c-component or head-and and-tail factorizations estimate the likelihood of ancestral graphs including bidirected edge(s).
D. Node and edge scores based on Normalized Maximum Likelihood criteria
MIIC_sedi setcssccac&hsrr''ss 's two-step implementation is illustrated with running example on Fig. D.I It relies on Node scores
(Step 1), for edge orientation priming and edge removal, and on Edge orientation scores (Step 2). Both scores are based on Normalized Maximum Likelihood criteria, derived in this Appendix.
(i) Output from MIIC
(ii) Step 1: Edge orientation priming and Edge removal
(iii) Step 2: Edge orientation
V
V
V
V
X
X
X
Y
X
Y
Z
W
Z
W
Z
W
Z
W
Figure D. simple rumning example of MIICC search&score's two-step implementation. The graph obtained by MIIC used as starting point for from the minimization of the node scores Eqs. 14
two-step algorithm. (ii) Step (Node scores): Edge orientation primings (*) are obtained
34 Edges without priming either extremity are removed the end of Step I Z). (iii) Step 2 (Edge orientation scores) Edge orientations for directed or bidirected edges are obtained by optimizing their
(e.g. X --- local contributions (Table M to the global likelihood score (Eq, [2] restricted to ac-connected subsets containing up to two-collider paths. This amounts to minimizing each edge orientation score w.r.t its nodes' parents and spouses (Table [T. The orientations of the edges highlighted in red have been modified with respect to their orientations the MIIC network used as starting point for MIIC scor.
Search-and-scom methods based on likelihood estimates need to properly account for finite sample size, as cross-entropy minimization leads to ever more complex models, resulting in model overfitting for finite datasets. BIC regularization is valid in the asymptotic limit of very large datasets and leads to the following finite size corrections of the cross-information terms the likelihood decomposition Eq.
C
log N N
I(C) + - - II (1
rk)
(22)
k=
for categorical datasets, where F3 is the number of categories or levels of the kth variable of C.
However, BIC regularization tends to overestimate finite size corrections, leading to lower recall, general. In order to better take into account finite sample size, we used instead the (universal) Normalized Maximum Likelihood (NML) criterion (Shtarkov] 1987 Rissanen & Tabus 2005 Kontkanen & Myllymäki] 2007 Roos al. 2008). which amounts normalizing the likelihood function over all possible datasets with the same number N of samples.
Indeed, all terms Eq. 18 actually cancel two- by construction, whatever their factorized expression, except for the remaining joint-distribution over all variables, q(z, 2, y. w) = {(2 y/a, w)) ((x) q(w).
18

---

An Efficient Scarch-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
Moreover, as search-and-score structure learning methods need to rely heuristic rather than exhaustive search, have implemented computationally efficient search-and-score method based on the likelihood decomposition of ancestral graphs (Eq. 12 limited to the close surrounding vertices of each node and edge. These node and edge scores, detailed below, extend MIIC's S unshielded triple scores to higher-order local information scores including ac-connected subsets of vertices with maximum of two-collider paths.
Node score. We first used the factorized Normalized Maximum Likelihood (INML) complexity (Kontkanen & Myllymäki 2007 Roos et al., 2008) define a local score for each node X, which extends the decomposable likelihood Bayesian graphs given each node's parents, Eq. to all non-descendant neighbors, Pa X
e -NIX, P E DIIEN (P(X)
L -N. Sx)re()(XX
(23)
NH(X Pa 1- log log
(24)
NE Ee T log m5 Σ log CNty
(25)
II
(26)
where nia corresponds the number of data points for which X; is in its kth state and its non-descendant neighbors in their jth state, with n, = > nzj-- The universal normalization constant Cj is then computed by summing the numerator over all possible partitions of the n data points into maximum of r subsets, l l2 +... & with & > 0.
4
=
Σ
(27)
1
t +1..tttt...
which can fact be computed in linear-time using the following recursion (Kontkanen & Myllymäki 2007).
C C
2r-2
(28)
r
with C₁ 1 for all and applying Eq. 31 below for 2. However, for large and r, C computation tends to be numerically unstable, which can be circumvented by implementing the recursion parametric complexity ratios D., = =C C//C/=1 rather than parametric complexities themselves (Cabeli et al. 2020) as,
Di = 1 +
n
(29)
2)
log C
Σ
log D
(30)
k=
for >> 3, with C2 = and C² = D2, which can be computed directly with the general formula, Eq. 27 for = 2,
n
h n-h
C₂
n
(31)
n
h=0
or Szpankowski approximation for large n (needed for m > 1000 in practice)
Szpankowski
2001
Kontkanen et al.
2003 Kontkanen 2009.
C2
TUT
+
2
O
(32)
2
3
12n
n2/22
TUT
3x
16
exp
(33)
9nm
36nM
19

---

An Efficient Scarch-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
This leads to the following local score for each node X which is minimized alternative combinations of non-descendant neighbors, Pa'x, Paxx Spx, U Nex in the first step of the local search-and-score algorithm (step 1) detailed in the main text,
Score, (x.)) = Pa
log Cn.
(34)
Note that Eq. 34 is defined for categorical variables for which the number of categories Tx, of variable X;, the number of combinatorial categories Q; of its non-descendant neighbors, and the corresponding number of samples n2 (for = 1. ..) are readily obtained from the dataset. However, this node score can also be extended for continuous (or mixed-type) variables, by subtracting the unconditional entropy H (X;) ) to the conditional entropy (X; Pa' Eq. 34 to obtain an equivalent node score (up (X:) constant) terms of mutual information, H(,, Paa (X.) ((x;; Pa' This extended node score then be estimated following the optimization procedure proposed by Cabeli et al 2020 to estimate (conditional) mutual information between continuous (or mixed-type) variables. Besides, this extended node score leads to the same definition of edge scores, listed in Table I and presented below for categorical variables, which can similarly be estimated for continuous (or mixed-type) variables using the same optimization procedure proposed by Cabeli et al 2020)
Edge scores. We now define several edge scores to optimize the orientation of each edge, X Y, given its close surrounding vertices.
To this end, we first introduce local score for node pairs which simply sums the node scores, Eq. 34 for each node. The resulting pair scores are listed in Table D.I for unconnected node pairs and for pairs of nodes connected by directed edge, where Pa x.y Pax USp and Pa nx = Pay U Sp,AX with their corresponding combinations of levels, Tuz and Inw
Table D.I. Local scores for node pairs
Pair score
Information
INML Complexity
XY
H(X Pa'
Paxxy)) H(Y Pa'xx)
N Σ Pey log Chty + Σ "ys log Cris
X Y
H(X Pa'
xxy) HY Pa' xX)
N
Σ Isw
log CPs E, Jake Tix log Cii
X K-Y Y
((X Pa xxy, Y) H(Y mxx)
N
Baw Ty log Chi + Σj log Chuy
Then, edge scores for directed edges, X + Y and Y + X, are defined w.r.t. to the edge removal score, X Y, by substracting the pair scores unconnected pairs to the pair scores of directed edges, leading to the following edge orientation scores,
9 t2
q Vis
Score (onoxx+310 +Y) -I(X;Y nx
Σ
log Cry
Σ
log Cryy
(35)
Bxg Tv
Score(Y + X) -I(X;Y x
Σ
log
Cr= ny
Σ
log CT= mpp
(36)
However, rz * ry. the fNML complexities of these orientation scores are not identical for Markov equivalent edge orientations between nodes sharing the same parents (or spouses) (Chickering 1995), Pa'nx = Pa Pa' and 4yh geg despite sharing the same conditional mutual information,
H()P)))))+(((P...))))) 1 (H(X ( H(Y H(Y
(37)
This suggests to symmetrize the fNML complexities for edge orientation scores by averaging them over each directed orientation, as for the conditional information in Eq. 37 leading to the proposed fNML complexity for directed edges given in Table the main text.
20

---

An Efficient Scarch-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
For bidirected edges, the proposed local orientation score accounts for all ac-connected subsets in close vicinity of the bidirected edge, which concerns all subsets including either X and any combination (possibly void) of parents spouses different from Y (i.e. corresponding to the information contributions H (X [Pa' xx )) or Y and any combination of parents or spouses different from X (i.e. corresponding to the information contributions (Y (P }) or, else, including both nodes X and Y plus any combination of their parents or spouses, corresponding to the following information contribution, -I(X; ( Pa'xr), where Pa'xy = Pa'xy x-y U Pa'yx' This last term, (X Pa'xr), contains all the remaining information contributions once bidirected orientation score is given relative to the edge removal score (Table D.I) as for the two directed orientation scores, above. Finally, the symmetrized INML complexity associated with a bidirected edge should be computed with the whole set of conditioning parents spouses, Pa' x indicated in Table Note that this bidirected orientation score becomes also Markov equivalent to the two directed orientation scores, as required, when the nodes share the same parents and spouses, i.e. Pa'xr Pa'yx = =P xay and Lyy 4yu 220 in Table
E. Toy models, benchmark data generation and method settings
E.1. Toy models
orientation scores (Table effectively predict bidirected orientations when the end nodes do not share the same parents (Model 1), share some parents (Model 2) or when the bidirected edge is part of a longer than two-collider paths (Model 3).
Fig. shows three simple ancestral models used test
XI
XS
X2
XS
X1
X1
Model
Model
Model
Figure E.1. Simple ancestral graphs.
The data is generated from the theoretical DAG using the mvvAAG function in the pealg package (Kalisch et al. 2012). Each node follows normal distribution, and the data discretized using bnlearn's discretize function using Hartemink's pairwise mutual information method (Scutari, 2010). For these toy models, the edge orientation scores are computed assuming the correct parents of each node.
The prediction of the edge orientation scores summarized in Table E. in % of replicates displaying directed edges (wrong) or bidirected edge (correct) as function of increasing dataset size N.
Table F.I.
Model 1, X₂
X1
Model 2, X² - X
Model 3, X² X,
Model 3, X,
X
N
{-
++
++
+-
1000
92
6
42
54
4
98
0
0
100
5000
14
20
66
14
12
74
100
0
0
100
10000
100
0
76
24
20000
0
100
100
100
0
0
100
35000
0
0
100
100
100
0
100
50000
0
0
100
100
100
0
0
100
E.2. Generation of linear and non-linear continuous datasets
The benchmark ancestral graphs for continuous datasets were obtained by hiding 0%, 10% or 20% of variables in either
linear or pon-lineas continuous Bayesian networks 50 or 150 nodes connectivity degree between and and average degree of 3 or corresponding to 360 ancestral graphs in total.
21

---

An Efficient Scarch-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
Undirected graphs were generated using igraph's cample_degsed function with the vl option https: /ittpp://igraph... org/r
html /1.2 1.2.5..... wawple_degsee.....eeeee html DAGs were then obtained using random ordering of their nodes.
Linear Gaussian datasets were generated from these DAGs using rmvDAG with uniformly distributed weights
in -011 and
topics/ rmDDG)).
Non-linear multimodal datasets were generated from these DAGs using Gaussian mixtures and non-linear structural equation models following (Cabeli et al. 2020):
Nodes without parents: Gaussian mixtures with 4 modes randomly distributed in [-20, 20] with Gaussian noise I, and subsequently rescaled betwenen 0 and 1.
Nodes with parents: non-linear structural equation models including pairwise coupling between the parents and non-linear transforms (x", B E [1, exp(x), sin()), cos())) (Cabeli et al. 2020
E.3. Causal discovery methods and parameter settings
MIIC_search&sscoresss source code is accessible httpe//ptttteeee.ttt....ooaa.aaa
MIIC package is accessible at http://////ttttaoooo//o///aaaaaa
MIIC and M settings were set described in section 3 of the main text.
M3HC was obtained at
GFCI was https//iiddtttteee/t/t/tooo///nn for linear data: depth 5, max_
and run with the following parameter settings:
MuxCondSeMM3HC = 00 TABUSi skeleton MMPC.
run with the following parameter settings 5, SEM_BIC (penalty_discouuul discount = 5), Fisher Z test 0.05) and for non-linear data: depth max_mmiieeeee leax_disc_path_heee length 2, basis function BIC (truncation limit = 3, genalty_discoune 2),
degenerate Gaussian test (a = 0.05).
DAG-GNN was obtained from ht http:::/ggthubb..o//rooaooooo/o
fwr__rom__oooe GNN and run using the follow- ing key parameters: encoder = MI.P, decoder = MLP,, optimizer = Adam, graph_ haapphhttrhhhdd--e = 0.3, h_tol = 1e-8.
The training hyperparameters were adapted to the number of variables (N) as follows:
N 100, epochs 300, batch_size = 20, = 0.002, mmx iter = 100
N 250, epochs 300, batch_ size_ 25, Ir 0.002, k_ _______iiiii iter 100
N 500, epochs 200, batch_size = 50, Ir = 0.001, k_max_iter = 80
N 1000, epochs 200, batch size 100, 0.001, m_x_maa___________aaaa iter 80
N 5000, epochs = 100, batch__siaaee 250, Ir 0.0001, k_max__tmeeeee 50
N 10000, epochs 100, batch size = 500, Ir = 0.0001, k max jter 50
N 20000, epochs 75, batch_atreeee 1000, Ir 5e-05, k_mak_iteeeeee 20
and run with -C-coniitional independence test for categorical data and default parameter a = 0.05. For bootstrapped data, FCI was run with the following param- eters: Alarm - a = 0.05; Insurance a = 0.005, depth = 2. max _ath_length = 1; Barley - a 0.0005, depth 1,
FCI was obtained at httpttts::// http:/aciaatata//////..o.//aaaaa
max lox__pah_lengtte Mildew - a = 0.0005, depth = I, max max_pathh maath_lenghh 0.
22

---

An Efficient Scarch-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
Continuuus
average degree
nodes, average degree NON-LINEAR
LINEAR
NON-LINEAR
LINEAR
Precision
Recall
Precision
Recall
Precsion
Recall
Precision
Recall
0% LV
0% LV
0
10% LV
10% LV
20 % LV
20% LV
100
1,000
10,000
100
1,000
10,000
100 1.000 10,000
100
1,000
10.000
100
1,000
10.000
100
1.000
10,000
100 1,000
10,000
100
1,000 10,000
Sample Size (N)
Sample Size (N)
Contruous Models: 150 nodes, average degroe
Continuous Models: 150
average degroe
LINEAR
NON-LINEAR
LINEAR
NON-LINEAR
Predsion
Recal
Predisoo
Recall
Precision
Recall
Precision
Recall
035 LV
0.%
10% LV
10.%
20.%
20 10%
1,000
1,000 10,000
1,000
10.000
10,000
1,000
10,000
10,000
10,000 1,000 10.000
Sample Size (N)
Sample Size (N)
M3HC
GFCI
DAG-GNN
MIIC_ search@ccch
Figure E.2. Benchmark results on continuous datasets. Benchmark results are averaged over 30 independent ancestral graph models obtained by hiding 0%, 10% 20% of variables in linear Gaussian models 150 nodes of average degree 3 or 5. MICC search&sccore results are compared to MIIC results used as starting point for MIIC Seil_sserc&&score M3HC (Trintaafllouu & Tsamardinos 2016), GFCI (Ogarrio al. (0166), and DAG-GNN (Yu et al. 2019) Causal discovery performance is assessed in terms of Precision and Recall relative to the theoretical PAGs, while counting as false positive all correctly predicted edges but with different orientation as the directed bidirccted edges of the PAG. Error bars: 95% confidence interval.
complex non-linear multimodal models including 50
23

---

An Efficient Search-and-Score Algorithm for Ancestral Graphs using Multivariate Information Scores
ALARM (509 parameters)
INSURANCE (1,008 parameters)
BARLEY (114,005 parameters)
MILDEW (540,150 parameters)
Predision
Recall
Precision
Precision
0. LV
10% WV
0.
100
1,000
10,000
100
1,000
10,000
100
1,000
10,000
100
1,000
10,000
100
1,000
10,000
100
1,000
10,000
100
1,000
10,000
100
1,000
10,000
Sample Size (N)
FCI
MIIC
MIlC search&score
Figure E.3. Benchmark results on bootstrapped 'real-world' categorical datasets from the bnlearn repository. Benchmark results bootstrap sensitivity analysis sampling noise based on 30 independent resamplings with replacement of single categorical datasets of increasing sizes. Ancestral graphs are obtained by hiding 0%. 10% or 20% of variables in Discrete Bayesian Networks of increasing complexity (see main text): Alarm, Insurance, Barley and Mildew. Causal discovery performance is assessed in terms of Precision and Recall relative to the theoretical PAGs, while counting as false positive all correctly predicted edges but with a different orientation the dirceted or bidirected edges of the PAG. MIIC search&scorre results are compared MIIC results used as starting point for MIIC search&score and to FCI (Zheng al. 2024). The lack of FCI results for complex models at large sample sizes stems from FCI difficulty to converge on bootstrapped datasets. Error bars: 95% confidence interval.
24


Reviews:
[{"review_id": "cXpyOrICu6", "reviewer": "Reviewer_kxfn", "summary": "This paper presents a novel search-and-score algorithm for discovering ancestral graphs - a class of graphical models used to represent causal relationships with latent variables. The key contribution is a new normalized likelihood score based on multivariate information measures applied to ac-connected subsets of variables.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed review and for underlying that \"the paper is well-structured and easy to follow, with a clear explanation of the proposed method\" and that \"most claims are supported by clear and convincing evidence\".\n\n- Following your suggestions, we have extended our benchmarks to larger networks including up to 150 nodes and performed additional comparisons with the DL-based causal discovery method DAG-GNN (Zheng et al 2019). The results are included in the revised Fig.E2 accessible at the following double-blind github repository: https://tinyurl.com/34ej6j78\n\n- MIIC_search&score is shown to outperform all other tested causal discovery methods on complex non-Gaussian datasets, including non-linear couplings between variables, while GFCI and DAG-GNN are the best performers on linear Gaussian datasets. These results demonstrate GFCI, M3HC and DAG-GNN's clear advantage in assuming linear Gaussian distributions when analyzing multivariate Gaussian datasets. However, these results also highlight the clear limitation of this assumption, when analyzing more complex datasets including non-linear couplings between variables, which goes beyond the variable transformation of post-nonlinear causal models. By contrast, MIIC_search&score and MIIC make no particular assumption on the data distributions and achieve similar good performance across a broad range of Gaussian or non-Gaussian multimodal distributions (Figs 2 & E2) as well as on complex categorical datasets (Figs 3 & E3).\n\n- MIIC_search&score scalability is primarily limited by the quadratic complexity of MIIC wrt the number of nodes (see Fig S5 in Verny et al 2017). This is actually optimal considering that MIIC tests all pairs of variables for conditional independence (with only a small time increase when including latent variables, Fig S5 in Verny et al 2017, thanks to MIIC's greedy approach). By contrast, the two-step search-and-score scheme of MIIC_search&score is essentially linear in the numbers of nodes (step 1) and edges (step 2) for a fixed degree. In practice, Steps 1 and 2 take similar running time as for MIIC to provide the starting graph for MIIC_search&score. We will discuss the method scalability in the final revised version of the paper and also cite \"some recent scalable causal discovery methods\" notably based on Differential Causal Discovery and continuous optimization techniques (eg Lopez et al NeurIPS 2022, Montagna et al CLeaR 2023, Amin et al ICML 2024).\n\n- Concerning \"a lack of discussion on the robustness to noise and missing data\", we have actually performed a bootstrap sensitivity analysis to sampling noise in Fig. E3, which shows that our method is rather robust to sampling noise. We will extend the discussion about these results in the final revised version of the paper. Likewise, considering that MIIC allows for missing data, MIIC_search&score can also allow for missing data. However, we have not investigated this specificity in the present paper which focuses on the comparison of MIIC_search&score with alternative approaches, such as M3FC, GFCI and DAG-GNN, that do not allow for missing data.\n\n- Finally, we are not sure to understand the two following comments. We would be grateful for further clarification if our tentative responses below fail to address your concern. Here are the two points:\n\n    - \"The likelihood function approximation assumes limited higher-order dependencies, which may fail in highly non-linear causal systems.\" The likelihood function approximation implemented in the two-step algorithm concerns the maximum length of the collider paths defining the ac-connected subsets of vertices included in the likelihood function estimate (Eq.12). Hence this likelihood approximation actually concerns the hidden variables from the underlying DAG models, not the linear or non-linear nature of the causal relations. In particular, vertices with many parent nodes define large \"star-like\" ac-connected subsets whose likelihood contributions are all included in MIIC_search&score's two-step likelihood estimate (ie independently from their linear or non-linear combinatorial relations with their parents). Likewise, the approach recovers the correct likelihood of Bayesian networks (Eq.2) in absence of latent variables (ie regardless of the linear or non-linear nature of the underlying relations between variables).\n\n    - \"No ablation study to separate the impact of local vs. global refinement steps.\" It seems to us that an \"ablation study\" might not be required here as both MIIC and MIIC_search&score are intrinsically local methods without \"local vs global refinement steps\". In particular, MIIC_search&score's two-step approximation effectively separates the selection of relevant edges (step 1) and their orientation (step 2). Ablating the first step tends to retain more FP edges (with a corresponding loss in Precision), while ablating the second step returns an undirected graph.", "claims_and_evidence": "Most claims are supported by clear and convincing evidence. \n\nHere are some potential limitations:\n* In linear Gaussian settings, the method performs worse than M3HC and GFCI.\n* The paper does not discuss how the method scales beyond moderate-sized graphs (e.g., 100+ nodes).\n* The likelihood function approximation assumes limited higher-order dependencies, which may fail in highly non-linear causal systems.", "methods_and_evaluation_criteria": "Yes. The proposed method is evaluated on both Simulated datasets (both continuous and categorical) and Real-world datasets (Alarm, Insurance, Barley, Mildew). The Evaluation Metrics include Precision, Recall, F1-score, and Computational runtime.\n\nThe results would be more complete if they also discussed how the method scales beyond moderate-sized graphs (e.g., 100+ nodes).\nAlso, there is a lack of discussion on the robustness to noise and missing data.", "theoretical_claims": "I did not check the proofs.", "experimental_designs_or_analyses": "The experimental designs and analyses are reasonable in general. Two missing parts: 1. No analysis of performance on very large graphs.  2. No ablation study to separate the impact of local vs. global refinement steps.", "relation_to_broader_scientific_literature": "This paper builds on Bayesian networks and ancestral graph theory, and is closely related to mutual information-based causal discovery", "essential_references_not_discussed": "The paper does not cite some recent scalable causal discovery methods.", "other_strengths_and_weaknesses": "Strengths\n* The paper introduces a new likelihood function based on multivariate cross-information over ac-connected subsets.\n* The local node-based scoring + edge refinement strategy makes causal discovery more scalable than exhaustive search.\n* The paper is well-structured and easy to follow, with a clear explanation of the proposed method.\n\nWeaknesses\n* Scalability to large graphs is unclear\n* The method underperforms M3HC and GFCI in linear Gaussian models, but the reason for this is not discussed.\n* No comparison to deep learning-based approaches", "other_comments_or_suggestions": "N/A", "questions_for_authors": "The paper evaluates the proposed method on relatively small- to medium-sized graphs. How does the approach scale to larger graphs (e.g., 100+ nodes)?\n\nThe experimental results indicate that the method performs worse than M3HC and GFCI in linear Gaussian models. Could you clarify why this is the case?\n\nThe paper mainly compares the proposed method with traditional search-and-score approaches. How does it compare with deep learning-based causal discovery methods such as DAG-GNN (Zheng et al., 2018) or causal representation learning (Ke et al., 2019)?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "x8hjffgpFb", "reviewer": "Reviewer_Vzgr", "summary": "This paper presents a greedy search-and-score algorithm for ancestral graphs with latent variables, using multivariate information for efficiency. Experimental results verify that it outperforms existing methods in causal discovery like M3HC and MIIC.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review and for underlining that the \"theoretical claims are solid\" and \"the contributions have potential for the broader scientific literature, where causal discovery with latent variables is important\".\n\nConcerning your question about the generation of the benchmark datasets, the proposed setting based on hiding some variables from DAG models is not related to the limitation of the method in terms of a maximum of two-collider paths. Indeed, this setting comes from the very definition of ancestral graphs (Richardson and Spirtes 2002) which are obtained by marginalizing or conditioning on DAG vertices, as further discussed in our reply to Reviewer 7jRW. The same benchmark setting based on hidding some vertices in DAGs has been used in many earlier studies on constraint-based or score-based methods to discover MAGs or their Markov equivalent class representatives, the PAGs.\n\nIn the present article, we have used two versions of this usual setting projecting DAGs onto MAGs. First, one simple version with toy models and, then, one more general version for all other more advanced benchmarks.\n\nThe first simple setting is to introduce hidden common causes explicitly in benchmark toy models as detailed in Appendix E. Fig. E.1 shows the three simple ancestral models we have used to test MIIC_search&score orientation scores (Table 1) to effectively predict bidirected orientations when the end nodes do not share the same parents (Model 1), share some parents (Model 2) or when the bidirected edge is part of a longer than two-collider path (Model 3). The predictions of the edge orientation scores are summarized in Table E.1 and show good predictions for large enough datasets.\n\nThe second more general setting, which has also been used in a number of earlier similar studies benchmarking constraint-based or score-based methods including latent variables, is to hide a fraction of the nodes from an underlying DAG model to reconstruct the corresponding MAG and PAG, which can be theoretically obtained. The results are reported in Figs. 2 & 3 and Figs. E2 & E3.\n\nImportantly, the benchmark PAGs used to score the causal discovery methods with increasing proportions of latent variables (Figs. 2 & 3 and Figs. E2 & E3) include not only bidirected edges originating from hidden common causes but also additional directed or undirected edges arising, in particular, from indirect effects of hidden variables with observed parents. Irrespective of their orientations, all these additional edges originating from indirect effects of hidden variables generally correspond to weaker effects (i.e. lower mutual information of indirect effects due to the Data Processing Inequality) and are more difficult to uncover than the edges of the original DAG model without hidden variables. This explains the steady decrease in recall for complex ancestral models with higher proportions of hidden variables, while precision remains essentially unaffected, Figs. 2 & E2 and Figs. 3 & E3 (this is clearly apparent for complex random models with average degree 5 and complex real-world models, Insurance, Barley, Mildew).", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "Yes", "theoretical_claims": "I looked through the theoretical claims and did not find any specific issues.", "experimental_designs_or_analyses": "Yes. The benchmark datasets come from hiding several variables from existing datasets. It should be reasonable if there is no relationship between the proposed method and this setting. However, as mentioned in the paper, *\"the proposed method is limited to ac-connected subsets of vertices with a maximum of two-collider paths\"*. Is the hiding process related to this constraint? Is it possible that the proposed approach underperforms the baselines in other settings?", "relation_to_broader_scientific_literature": "The contributions have potential for the broader scientific literature, where causal discovery with latent variables is important.", "essential_references_not_discussed": "I do not know any related works that are essential to understanding the paper.", "other_strengths_and_weaknesses": "- The paper is well-written\n- The theoretical claims are solid", "other_comments_or_suggestions": "I do not have any other comments or suggestions", "questions_for_authors": "See above.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ADUuY3J1vW", "reviewer": "Reviewer_7jRW", "summary": "This paper proposes a greedy hybrid search-and-score algorithm to learn ancestral graphs from data with some latent variables marginalized out.\n\nFor this purpose, the authors first provides an explicit decomposition of the likelihood function of ancestral graphs in terms of multivari-\nate cross-information over relevant ‘ac-connected’ subsets of variables. Then, the authors use this decomposition to design a two-step approach using local information scores restricted to surrounding vertices and edges.\n\nThe proposed method is shown to outperform several state-of-the-art causal discovery methods.\n\n\n---\n\n\n## update after rebuttal:\n\nMy concerns about the usage of MAGs and the generalizability to cases with selection bias have been well addressed. My concerns about identifiability, though still existing, is kind of alleviated with the comparisons to other score-based works. I have therefore raised my score form 3 to 4.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed review. Thank you also for underlining that \"the problem studied (greedy learning of graphs involving latent variables) is very crucial and needed in the field\" and that \"the decomposition over ‘ac-connected’ subsets of variables shows a clear trajectory generalized from DAG case (Figure 1 G and F)\". We reply to your comments and suggestions below.\n\nConcerning your first point about the guarantees of the results, it is true that the greedy approach does not come with guarantees unless the actual model has no ac-connected subsets with more-than-two-collider paths. In this case, the local decomposition gives a correct estimate of the global likelihood, based on Proposition 3, due to the local consistency of the likelihood score (Eq.12) just like in the Bayesian network case. However, as suggested by this reviewer, an oracle version of the approach could also be proposed, as the general likelihood decomposition (Eq. 12) could in principle be used in conjunction with an exhaustive search-and-score algorithm over MAGs or PAGs, which can be generated rather efficiently (Hu et al 2020). Alternatively, and in practice for graphs with more than about 10-15 nodes, an MCMC algorithm based on the likelihood decomposition (Eq.12) could also be used to efficiently search for high-scoring MAGs or PAGs. \n\nThe connection with head-and-tail factorization is discussed in some details in Appendix C. In brief, these head-and-tail factorizations introduced by Richardson in 2009 are definitely useful to enable the parametrization of the joint probability distribution with independent parameters for ancestrally closed subsets of vertices (see Richardson 2009). However, head-and-tail factorizations are also somewhat involved, as they do not correspond to a single factorized equation but instead to multiple factorized equations for a given ancestral graph. Moreover, head-and-tail factorizations (as well as the c-component decompositions introduced by Tian & Pearl 2002) cannot simply be used to estimate the likelihood function in terms of empirical distribution p(x), as shown and discussed in Appendix C. This limits their utility for the purpose of scoring MAGs and PAGs. Hence, the likelihood decomposition over ‘ac-connected’ subsets of variables (Eq.12) is particularly useful for complex non-linear data (Figs. 2 & E2), as discussed in our reply to Reviewer kxfn, as BIC scores have been shown to give good approximation of likelihood functions for linear Gaussian data. \n\nConcerning the definition of MAGs, it is true that the original definition of MAGs in (Richardson & Spirtes 2002) includes also undirected edges, which are obtained by conditioning on selection variables in DAGs, while bidirected edges are obtained by marginalizing common causes in DAGs. However, following the usual terminology in the field (eg adopted in Triantafillou et al 2016 or Hu et al 2020), we also designate MAGs without undirected edges as MAGs throughout the paper, eventhough they are in fact a subclass of MAGs. Note, however, that this subclass of MAGs without undirected edges is more restricted than the class of Acyclic Directed Mixed Graphs (ADMG), which also contain directed and bidirected edges, but are not necessarily marginalized from DAGs. In particular ADMGs, unlike MAGs, allow for \"almost directed cycles\", where X->...->Y with X<->Y. In addition, the proposed likelihood decomposition (Eq.12) can only be applied to MAGs without undirected edge, as shown in the third section (iii) of the proof, Appendix B. \n\nYet, as pointed out by this reviewer, there are other existing scores under latent variables for specific parametric assumptions like discrete variables, exponential or stratified Gaussian families. This will be reminded in the final revised paper. In particular, discrete chain graph models (Drton 2009), fully bidirected graph models (Drton and Richardson 2008), and discrete nested Markov models (Richardson et al. 2017) have been shown to be curved exponential models, which can be scored consistently using BIC scores. Interestingly, one extension of these BIC scores has been recently proposed by (Bellot, Zhang, Bareinboim, AAAI 2024) to allow for the distinction between different ADMGs from (Verma & Pearl 1990), see Fig. 1c & d in (Bellot, Zhang, Bareinboim, AAAI 2024). These ADMGs imply the same set of conditional independence constraints and yet are distinguishable because they imply an equality between different functionals of the probability distribution, see Eq. 1 in (Bellot, Zhang, Bareinboim, AAAI 2024). It would definitely be interesting, in future follow-up studies, to further explore how our information theoretical score (Eq.12) generalizes or differs from these other existing scores under latent variables, as suggested by this reviewer.", "claims_and_evidence": "Yes.", "methods_and_evaluation_criteria": "Yes.", "theoretical_claims": "I read the theorems, the running examples in Figure 1 in details, and skimmed through the proofs. The results seem correct to me but I cannot guarantee.", "experimental_designs_or_analyses": "Yes. The experimental design to show the empirical gain on causal discovery accuracy (especially on recall) seems reasonable to me.", "relation_to_broader_scientific_literature": "/", "essential_references_not_discussed": "/", "other_strengths_and_weaknesses": "Strengths:\n\n1. The problem studied (greedy learning of graphs involving latent variables) is very crucial and needed in the field.\n\n2. The proposed decomposition of the likelihood function of ancestral graphs is new to me. The decomposition over ‘ac-connected’ subsets of variables show a clear trajectory generalized from DAG case (Figure 1 G and F).\n\n3. The experimental results show a generally better performance than other existing MAG learning algorithms.\n\n\nWeaknesses:\n\n1. There lacks an identifiability or correctness guarantee of the results.\n   - It seems that there is no guarantee of the results–at least I didn't see any theoretical claims in the whole Section 3–but please correct me and formally state such guarantee if I am wrong.\n   - The authors state that for computational efficiency, the proposed algorithm only involves two steps to search over local scores of each node and edges. Then, at the first place, when ignoring efficiency, do authors have an oracle version of algorithm that traverses ac components of all sizes that have correctness guarantee? Or, can the proposed likelihood decomposition be applied to some exhaustive search method that have guarantee?\n   - Following the proposed score decomposition, are there any good properties of the score to support the greedy search (e.g., such greedy search won't lead to local optima), just like the local consistency of BIC score in the Bayesian network and linear Gaussian case?\n\n\n2. The proposed score decomposition, though being new, needs more justification.\n   - The authors say that their ac-component-based decomposition \"do not rely on the head-and-tail factorization but coincide with the parametrizing sets (Hu & Evans, 2020) derived from the head-and-tail factorization...\". But what is the exact advantage of this score decomposition than existing ones, in both the characterization form itself, and for the estimation in the algorithms? Is \"head-and-tail factorization\" something unwanted that we need to prevent from using?\n   - Other existing scores under latent variables, but with specific parametric assumptions like discrete variables, exponential or stratified Gaussian families, should also be discussed. It would be interesting to see how this information theoretical score characterization generalizes them or differs from them.", "other_comments_or_suggestions": "Throughout this paper the term \"MAG\" is used. However, the authors actually only allows for latent variables but not selection bias, and thus only -> and <-> occurs, without -- edges. In this case, for convention, maybe the authors want to turn to the term \"ADMG\" instead? Or, are there any specific reasons to use \"MAG\"? Can the score characterization and decomposition also be applied to cases with selection bias?", "questions_for_authors": "/", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "dg5Kw0ebTl", "reviewer": "Reviewer_Firz", "summary": "This paper introduces a greedy search-and-score algorithm for Ancestral Graphs (AGs), which contain directed and bidirected edges to account for latent variables. The key innovation is a normalized likelihood score based on multivariate information over ac-connected subsets (subsets of vertices connected through collider paths within their ancestor set). The proposed method outperforms state-of-the-art causal discovery techniques (MIIC, M3HC, and GFCI) on benchmark datasets.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive review and interesting suggestions to improve the manuscript.\n\nWe agree that a running example could illustrate the proposed method and its two-step implementation. We will include such an example in the final revised version to outline the main ideas of the method and facilitate the understanding of the paper, which might be a bit \"dense in notations, concepts, and terminology\". We have posted a first tentative Figure with this running example at the following double-blind github repository: https://tinyurl.com/34ej6j78\n\nConcerning the theoretical guarantee of the approach, we will also clarify in the revised version the conditions for which theoretical guarantees exist. They concern models with (i) a maximum of two-collider paths in the method setting presented here and more generally the use of the proposed likelihood score (Eq.12) in conjunction with (ii) an exhaustive search-and-score method over small-sized MAG models or with (iii) an MCMC algorithm to efficiently sample larger MAG models, as outlined in our reponse to Reviewer 7jRW.\n\nFinally, following your suggestion, we will discuss the effect of sparsity of the graphs on the method performance, which is already visible on Figure E.2. Namely, graphs with lower average connectivity lead to better performance for all methods including MIIC_search&score. This is actually expected as sparser models correspond to (much) fewer variable combinations to consider (eg to estimate local likelihood contributions of relevant ac-connected subsets), leading to more robust rankings between alternative sparse models and therefore to better performance at predicting correct models.\n\nFinally, following your suggestion, we will discuss the effect of sparsity of the graphs on the method performance, which is already visible on Figure E.2. Namely, graphs with lower average connectivity lead to better performance for all methods including MIIC_search&score. This is actually expected as sparser models correspond to (much) fewer variable combinations to consider (eg to estimate local likelihood contributions of relevant ac-connected subsets), leading to more robust rankings between alternative sparse models and therefore to better performance at predicting correct models.", "claims_and_evidence": "- The claim that the likelihood of ancestral graphs can be explicitly decomposed using multivariate cross-information over ac-connected subsets is formally proven in Theorem 1.\n- The claim that proposed search-and-score method outperforms MIIC, M3HC, and GFCI on benchmark datasets is supported by experimental results (Figures 2 and 3).", "methods_and_evaluation_criteria": "- The methodology is well-structured: derives a novel likelihood function for ancestral graphs, and provides an efficient search-and-score algorithm\n- Experiments include synthetic Bayesian networks and benchmark datasets). Comparisons with state-of-the-art methods (MIIC, GFCI, M3HC) are well-presented.", "theoretical_claims": "The theoretical contributions include:\n- Theorem 1: Likelihood decomposition for ancestral graphs. \n- Corollary 2: ac-connected subsets define Markov equivalence classes. \n\nThe mathematical formulations and proofs appear sound.", "experimental_designs_or_analyses": "- The experimental design is generally valid, with well-chosen benchmarks and valid evaluation metrics.\n- The comparisons with MIIC, M3HC, and GFCI are meaningful, as these represent reasonable baselines.", "relation_to_broader_scientific_literature": "Structure learning is widely used in scientific research for exploratory analysis. Developing efficient algorithms for this task is significant for more accurate and scalable discovery.", "essential_references_not_discussed": "Not sure due to limited familiarity.", "other_strengths_and_weaknesses": "**Strengths**:\n- It provdes novel likelihood function decomposition for ancestral graphs.\n- The proposed method outperforms state-of-the-art methods on complex datasets.\n- A computationally efficient heuristic for large graphs is provided.\n\n**Weaknesses**:\n- The paper is dense in notations, concepts, and terminology. It is extremely hard to process for average readers.", "other_comments_or_suggestions": "- It would be better to provide a simple running example to illustrate the proposed method.", "questions_for_authors": "- Is there any theoretical guarantee can be said on the proposed scoring algorithm? \n- How do we expect the sparsity of the graph to affect the performance?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper introduces a novel greedy score-based algorithm for learning ancestral graphs, enabling discovery in the presence of latent confounders. Central to the contribution is a new decomposition of the likelihood function over ac-connected subsets, leading to a two-step local search procedure that balances accuracy and computational traceability. Reviewers generally agree that the paper is well-written, the methodology is well-motivated, and the experimental results are compelling. Reviewer 7jRW initially raised concerns regarding the theoretical guarantees and the role of the proposed score decomposition, but found the authors’ rebuttal satisfactory and raised their score. Other reviewers pointed out limited discussion of scalability to very large graphs and robustness under noise; the authors have addressed this by providing additional experiments to support the method’s applicability in these scenarios. Given the overall strength of the contributions, the clarity of exposition, and the effective rebuttal addressing most concerns, I recommend acceptance.
Decision: Accept (poster)