# An efficient search-and-score algorithm for ancestral graphs using multivariate information scores

Nikita Lagrange, Herve Isambert

<figure><img src="image_1.png" /></figure>

## To cite this version:

Nikita Lagrange, Herve Isambert. An efficient search-and-score algorithm for ancestral graphs using multivariate information scores. 2025. hal-04957997

HAL Id: hal-04957997

https://hal.science/hal-04957997v1

Preprint submitted on 20 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An efficient search-and-score algorithm for ancestral graphs using multivariate information scores

Nikita Lagrange, Hervé Isambert

CNRS, UMR168, Institut Curie, PSL University, Sorbonne University, Paris, France
nikita.lagrange@curie.fr, herve.isambert@curie.fr

## 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 subsets” of vertices, C, that are connected through collider paths confined to the ancestor set of 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.

## 1 Introduction

The likelihood function plays a central role in the selection of a graphical model G based on observational data D. Given N independent samples from D, the likelihood $L_{D|G}$ that they might have been generated by the graphical model G is given by [1],

$$ \mathcal{L}_{D|G} = \frac{1}{Z_{D,G}} \exp(-NH(p,q)) \quad (1) $$

where $H(p, q) = - \sum_{\boldsymbol{x}} p(\boldsymbol{x}) \log q(\boldsymbol{x})$ is the cross-entropy between the empirical probability distribution $p(\boldsymbol{x})$ of the observed data D and the theoretical probability distribution $q(\boldsymbol{x})$ of the model G and $Z_{D,G}$ a data- and model-dependent factor ensuring proper normalization condition for finite dataset. In short, Eq.1 results from the asymptotic probability that the N independent samples, $\boldsymbol{x}^{(1)}, \dots, \boldsymbol{x}^{(N)}$, are drawn from the model distribution, $q(\boldsymbol{x})$, i.e. $L_{D|G} \equiv q(\boldsymbol{x}^{(1)}, \dots, \boldsymbol{x}^{(N)}) = \prod_i q(\boldsymbol{x}^{(i)})$, rather than the empirical distribution, $p(\boldsymbol{x})$. This leads to, $\log L_{D|G} = \sum_i \log q(\boldsymbol{x}^{(i)})$, which converges towards $N \sum_{\boldsymbol{x}} p(\boldsymbol{x}) \log q(\boldsymbol{x}) = -NH(p,q)$ in the large sample size limit, $N \to \infty$, with $\log Z_{D,G} = O(\log N)$.

The structural constraints of the model G translate into the factorization form of the theoretical probability distribution, $q(\boldsymbol{x})$ [2–6]. In particular, the probability distribution of Bayesian networks (BN) factorizes in terms of conditional probabilities of each variable given its parents, as $q_{BN}(\boldsymbol{x}) = \prod_i q(x_i|\mathbf{pa}_{X_i})$, where $\mathbf{pa}_{X_i}$ denote the values of the parents of node $X_i$ in G, $\mathbf{Pa}_{X_i}$. For Bayesian networks, the factors of the model distribution, $q(x_i|\mathbf{pa}_{X_i})$, can be directly estimated with the empirical conditional probabilities of each node given its parents as, $q(x_i|\mathbf{pa}_{X_i}) \equiv p(x_i|\mathbf{pa}_{X_i})$, leading to the well known estimation of the likelihood function in terms of conditional entropies $H(X_i|\mathbf{Pa}_{X_i}) = - \sum_{\boldsymbol{x}} p(x_i, \mathbf{pa}_{X_i}) \log p(x_i|\mathbf{pa}_{X_i})$,

$$ \mathcal{L}_{D|G_{BN}} = \frac{1}{Z_{D,G_{BN}}} \exp \left( -N \sum_{\boldsymbol{x} \in V} H(X_i|\mathbf{Pa}_{X_i}) \right) \quad (2) $$
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 ancestral graphs, which include directed as well as bidirected edges, arising from the presence of unobserved latent variables. Tian and Pearl 2002 [7] 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 [6] later proposed a 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 subsets of vertices. However, unlike with Bayesian networks, the contributions of c-components or head-and-tail factors to the likelihood function cannot simply be estimated in terms of empirical distribution p(x), as shown below. This leaves the likelihood function of ancestral graphs difficult to estimate from empirical data, in general, although iterative methods have been developed when the data is normally distributed [8–13].

The present paper provides an explicit decomposition of the likelihood function of ancestral graphs in terms of multivariate cross-information over relevant ‘ac-connected’ subsets of variables, Figs. 1., which do not rely on the head-and-tail factorization but coincide with the parametrizing sets [14] derived from the head-and-tail factorization. It suggests a natural estimation of these revelant contributions to the likelihood function in terms of empirical distribution p(x). This result extends the likelihood expression of Bayesian Networks (Eq. 2) to include the effect of unobserved latent variables and enables the implementation of a greedy search-and-score algorithm for ancestral graphs. 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 is shown to outperform state-of-the-art causal discovery methods on challenging benchmark datasets.

# 2 Theoretical results

## 2.1 Multivariate cross-entropy and cross-information

The theoretical result of the paper (Theorem 1) is expressed in terms of multivariate cross-information derived from multivariate cross-entropies through the Inclusion-Exclusion Principle. The same expressions can be written between multivariate information and multivariate entropies by simply substituting $q(\{x_i\})$ with $p(\{x_i\})$ in the equations below and will be used to estimate the likelihood function of ancestral graphs (Proposition 3).

As recalled above, the cross-entropy between $m$ variables, $V = \{X_1, \dots, X_m\}$, is defined as,

$$H(V) = -\sum_{\{x_i\}} p(x_1, \dots, x_m) \log q(x_1, \dots, x_m) \quad (3)$$

where $p(\{x_i\})$ is the empirical joint probability distribution of the variables $\{X_i\}$ and $q(\{x_i\})$ the joint probability distribution of the model. Bayes formula, $q(\{x_i\}, \{y_j\}) = q(\{x_i\}|\{y_j\}) q(\{y_j\})$, directly translates into the definition of conditional cross-entropy through the decomposition,

$$H(\{X_i\}, \{Y_j\}) = H(\{X_i\}|\{Y_j\}) + H(\{Y_j\}) \quad (4)$$

Multivariate (cross) information, $I(V) \equiv I(X_1; \dots; X_m)$, are defined from multivariate (cross) entropies through Inclusion-Exclusion formulas over all subsets of variables [15–18] as,

$$ 
\begin{aligned}
I(X) &= H(X) \\
I(X; Y) &= H(X) + H(Y) - H(X, Y) \\
I(X; Y; Z) &= H(X) + H(Y) + H(Z) - H(X, Y) - H(X, Z) - H(Y, Z) + H(X, Y, Z) \\
I(V) &= -\sum_{S \subseteq V} (-1)^{|S|} H(S) \quad (5)
\end{aligned}
$$

where the semicolon separators are needed to distinguish multipoint (cross) information from joint variables as $\{X, Z\}$ in $I(\{X, Z\}; Y) = I(X; Y) + I(Z; Y) - I(X; Y; Z)$. Below, implicit separators between non-conditioning variables in multivariate (cross) information will always correspond to semicolons, e.g. as in $I(V)$ in Eq. 5. Unlike multivariate (cross) entropies, which are always positive,
$H(X_1, \dots, X_k) \ge 0$, multivariate (cross) information, $I(X_1; \dots; X_k)$, can be positive or negative for $k \ge 3$, while they remain always positive for $k < 3$, i.e. $I(X; Y) \ge 0$ and $I(X) \ge 0$.

In turn, multivariate (cross) entropies can be expressed through the Principle of Inclusion-Exclusion into the same expression form but in terms of multivariate (cross) information,

$$H(V) = - \sum_{S \subseteq V} (-1)^{|S|} I(S), \qquad (6)$$

Conditional multivariate (cross) information $I(V|Z)$ are defined similarly as multivariate (cross) information $I(V)$ but in terms of conditional (cross) entropies as,

$$I(V|Z) = - \sum_{S \subseteq V} (-1)^{|S|} H(S|Z) \qquad (7)$$

Eqs. 5 & 7 lead to a decomposition rule relative to a variable Z, Eq. 8, which can be conditioned on a set of joint variables, $A = \{A_1, \dots, A_m\}$, with implicit comma separators for conditioning variables in Eq. 9,

$$I(V) = I(V|Z) + I(V; Z) \qquad (8)$$

$$I(V|A) = I(V|Z, A) + I(V; Z|A) \qquad (9)$$

Alternatively, conditional (cross) information, such as $I(X; Y|A)$, can be expressed in terms of non-conditional (cross) entropies using Eq. 4,

$$\begin{aligned} I(X; Y|A) &= H(X|A) + H(Y|A) - H(X, Y|A) \\ &= H(X, A) + H(Y, A) - H(X, Y, A) - H(A) \end{aligned} \qquad (10)$$

which can in turn be expressed in terms of non-conditional (cross) information as,

$$\begin{aligned} I(X; Y|A) &= I(X; Y) - \dots - (-1)^k \sum_{i_1 < \dots < i_k} I(X; Y; A_{i_1}; \dots; A_{i_k}) + \dots + (-1)^m I(X; Y; A_1; \dots; A_m) \\ &= \sum_{S' \subseteq S} (-1)^{|S'|} I(S'), \end{aligned} \qquad (11)$$

where $S = \{X, Y\} \cup A$. This corresponds, up to an opposite sign, to all (cross) information terms including both X and Y in the expression of the multivariate (cross) entropy, $H(X, Y, A)$, Eq. 6.

## 2.2 Graphs and connection criteria

### 2.2.1 Directed mixed graphs and ancestral graphs

Two vertices are said to be **adjacent** if there is an edge (of any type) between them, $X*-*Y$, 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 \to Y$ and **spouses** if $X \leftrightarrow Y$ in $\mathcal{G}$.

A **path** in $\mathcal{G}$ is a sequence of distinct vertices $V_1, \dots, V_n$ consecutively adjacent in $\mathcal{G}$, as, $V_1*-*V_2*-*\dots*-*V_{n-1}*-*V_n$. In particular, a **collider path** between $V_1$ and $V_n$ has the form $V_1* \to V_2 \leftrightarrow \dots \leftrightarrow V_{n-1} \leftrightarrow *V_n$ and a **directed path** corresponds to $V_1 \to V_2 \to \dots \to V_n$.

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 \to \dots \to Y$. $An_G(Y)$ denotes the **set of ancestors** of Y in $\mathcal{G}$. By extension, for any subset of vertices, $C \subseteq V$, $An_G(C)$ denotes the set of ancestors for all $Y \in C$ in $\mathcal{G}$.

A **directed mixed graph** is a vertex-edge graph $\mathcal{G} = (V, E)$ that can contain two types of edges: directed (→) and bidirected (↔) edges.

A **directed cycle** occurs in $\mathcal{G}$ when $X \in An_G(Y)$ and $X \leftarrow Y$. An **almost directed cycle** occurs when $X \in An_G(Y)$ and $X \leftrightarrow Y$.

**Definition 1.** An **ancestral graph** is a directed mixed graph:

i) without directed cycles;

ii) without almost directed cycles.

An **ancestral graph** is said to be **maximal** if every missing edge corresponds to a structural independence. If an ancestral graph $\mathcal{G}$ is not maximal, there exists a unique maximal ancestral graph $\tilde{\mathcal{G}}$ by adding bidirected edges to $\mathcal{G}$ [8].
### 2.2.2 ac-connecting paths and 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, discussed 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.

**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 \ast \to Z_1 \leftrightarrow \dots \leftrightarrow Z_K \leftrightarrow Y$, with all $Z_i \in An_G(\{X, Y\} \cup C)$, that is, with $Z_i$ in C or connected to $\{X, Y\} \cup C$ by an ancestor path, i.e. $Z_i \to \dots \to T$ with $T \in \{X, Y\} \cup C$.

**Definition 3.** [ac-connected subset] A subset C is said to be ac-connected if $\forall X, Y \in C$, X and Y are connected (through any type of edge) or there is an ac-connecting path between X and Y given C.

## 2.3 Likelihood decomposition of ancestral graphs

**Theorem 1.** [likelihood of ancestral graphs] The cross-entropy $H(p, q)$ and likelihood $\mathcal{L}_{D|G}$ of an ancestral graph G is decomposable in terms of multivariate cross-information, $I(C)$, summed over all ac-connected subsets of variables, C (Definition 3),

$$H(p, q) = - \sum_{C \subseteq V}^{\text{ac-connected}} (-1)^{|C|} I(C)$$
$$\mathcal{L}_{D|G} = \frac{1}{Z_{D,g}} \exp \left( N \sum_{C \subseteq V}^{\text{ac-connected}} (-1)^{|C|} I(C) \right) \quad (12)$$

where N is the number of iid samples in the dataset D and $Z_{D,g}$ a data- and model-dependent normalization constant.

The proof of Theorem 1 is left to Appendix B. It is based on a 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<fnref n="1" /> Hu and Evans [14] proposed an equivalent result (Proposition 3.3 in [14]) with a proof using head-and-tail decomposition to define parametrizing sets, which happen to coincide with the ac-connected sets defined here (Definition 3). Theorem 1 characterizes in particular the Markov equivalence class of ancestral graphs [8, 19–24] as,

**Corollary 2.** Two ancestral graphs are Markov equivalent if and only if 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 representatives, the partial ancestral graphs (PAG) [8, 25–27].

**Proposition 3.** The likelihood decomposition of ancestral graphs (Eq. 12, Theorem 1) can be estimated by replacing the model distribution q by the empirical distribution p in the retained multivariate cross-information terms $I(C)$ corresponding to all ac-connected subsets of variables, C.

Hence, Proposition 3 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 conditional 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.

Fig. 1 illustrates the cross-entropy decomposition for a few graphical models in terms of cross-information contributions from their ac-connected subsets of vertices. In particular, an unshielded non-collider (e.g. $X \to Z \to W$, Fig. 1A), is less likely (i.e. higher cross-entropy) than an unshielded collider or ‘v-structure’ (e.g. $X \to Z \leftarrow W$, Fig. 1B), if the corresponding three-point information term is negative, $I(X; Z; W) < 0$, in agreement with earlier results [28, 29]. However, this early

<footnote n="1">¹The genesis of Theorem 1 and Proposition 3 is further discussed in https://openreview.net/forum?id=Z2f4Laqi8U&noteId=8GLWeaAKc9.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: **Cross-entropy decomposition of ancestral graphs.** Examples of cross-entropy decomposition of ancestral graphs (red edges, rhs) in terms of relevant multivariate cross-information contributions *I*(*C*) with *C* ⊆ *V* (red nodes, rhs). Simple graphs: (A) without unshielded colliders, (B) with a single or non-overlapping unshielded colliders, (C) with overlapping unshielded colliders through three or more (conditionally) independent parents or (D) through a two-(or more)-collider path. (E) Bayesian graph corresponding to the head-and-tail factorization of the two-collider path in (D) estimated using the empirical distribution *p*()., see Appendix C. (F) Simple Bayesian graph not Markov equivalent to an ancestral graph (G) sharing the same edges and unshielded collider [24]. 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 to *C* by an ancestor path, marked with dashed edges, see Definition 2. By contrast, graphs **H** and **I** 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) [24].</figcaption></figure>

approach, exploiting the sign and magnitude 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. 1C. Similarly, the cross-entropies of collider paths involving several colliders also include higher-order terms, as with the simple example of a two-collider path, Fig. 1D. By contrast, the cross-entropy based on the head-and-tail factorization of the same two-collider path, i.e. $q(x, z, y, w) = q(z, y|x, w)q(x)q(w)$ [6], is found to be equivalent to the cross-entropy of a Bayesian graph without bidirected edge, Fig. 1E, when estimated with the empirical distribution $p(.)$, see Appendix C. This observation illustrates the difficulty to estimate the likelihood functions of ancestral graphs using head-and-tail factorization.

Further examples of graphical models, Figs. 1F-I, show the relative simplicity of the decomposition with only few (non-trivial) *ac*-connected contributing subsets $C$ with $|C| \ge 3$, as compared to the much larger number of non-*ac*-connected non-contributing subsets, that cancel each other by construction due to conditional independence constraints of the underlying model. Note, in particular, that most contributing multivariate information $I(C)$ only concern direct connections or collider paths within a single component subgraph induced by $C$ (solid line edges in Fig. 1). However, occasionally, collider paths extending beyond $C$ into $An_G(C) \setminus C$ (marked with wiggly edges) with corresponding ancestor path(s) (marked with dashed edges) do occur, as shown in Fig. 1G.

In addition, the present information-theoretic decomposition of the likelihood of ancestral graphs can readily distinguish their Markov equivalence classes according to Corollary 2. For instance, the ancestral graphs of Fig. 1F and Fig. 1G, despite sharing the same edges and the same unshielded collider ($X \to Z \leftarrow T$), turn out not to be Markov equivalent, as discussed in [24]. Indeed, their cross-entropy decompositions differ by two *ac*-connected contributing terms: a three-point cross information $I(X; Y; T)$ with a collider path not confined in $C$ (i.e. $X \rightsquigarrow Z \rightsquigarrow T \leftarrow Y$ and corresponding ancestor path $Z \dashrightarrow Y$) and a four-point information term $I(X; Y; Z; T)$ due to the two-collider path ($X \to Z \leftarrow T \leftarrow Y$). More quantitatively, it shows that the graph of Fig. 1G with a two-collider path is more likely than the graph of Fig. 1F whenever $I(X; Y; T) - I(X; Y; Z; T) = I(X; Y; Z) = I(X; Y|Z) - I(X; Y|Z, T) < 0$. Finally, the Markov equivalent graphs of Fig. 1H and Fig. 1I, also due to [24], illustrate the fact that the actual ancestor collider path between unconnected pairs does not need to be 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 1 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 [11–13, 30–32]. Our specific objective is not to develop an exact method limited to simple graphical models with a few nodes and small datasets but to implement an efficient and reliable heuristic method applicable to more challenging graphical models and large datasets.

Indeed, search-and-score structure learning methods need to rely on heuristic rather than exhaustive search, in general, given that the number of ancestral graphs grows super-exponentially as the number of vertices increases. This can be implemented for instance with a Monte Carlo algorithmic 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 [29, 33, 34], as starting point for a subsequent search-and-score approach based on the proposed likelihood estimation of ancestral graphs (Eq. 12 and Proposition 3).

Moreover, while the likelihood decomposition of ancestral graphs may involve extended *ac*-connected subsets of variables, as illustrated in Fig. 1, we aim to implement a computationally efficient search-and-score causal discovery 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 subgraphs, as detailed below.

The proposed method is shown to outperform MIIC and other state-of-the-art causal discovery methods on challenging datasets including latent variables.

## 3.1 MIIC, an hybrid causal discovery method based on unshielded triple scores

MIIC is an hybrid causal discovery method combining constraint-based and information-theoretic frameworks [29, 35]. Unlike traditional constraint-based methods [4, 5], MIIC does not directly
attempt to uncover conditional independences but, instead, iteratively subtracts the most significant three-point (conditional) information contributions of successive contributors, $A_1, A_2, ..., A_n$, from the mutual information between each pair of variables, $I(X;Y)$, as,

$$I(X;Y) - I(X;Y; A_1) - I(X;Y; A_2|A_1) - \dots - I(X;Y; A_n|\{A_i\}_{n-1}) = I(X;Y|\{A_i\}_n) \quad (13)$$

where $I(X;Y; A_k|\{A_i\}_{k-1}) > 0$ is the positive information contribution from $A_k$ to $I(X;Y)$ [28, 36]. Conditional independence is eventually established when the residual conditional mutual information on the right hand side of Eq. 13, $I(X;Y|\{A_i\}_n)$, becomes smaller than a complexity term, i.e. $k_{X;Y|\{A_i\}}(N) \ge I(X;Y|\{A_i\}_n) \ge 0$, which depends on the considered variables and sample size $N$.

This leads to an undirected skeleton, which MIIC then (partially) orients based on the sign and amplitude of the regularized conditional 3-point information terms [28, 29]. In particular, negative conditional 3-point information terms, $I(X;Y; Z|\{A_i\}) < 0$, correspond to the signature of causality in observational data [28] and lead to the prediction of a v-structure, $X \to Z \leftarrow Y$, if $X$ and $Y$ are not connected in the skeleton. By contrast, a positive conditional 3-point information term, $I(X;Y; Z|\{A_i\}) > 0$, implies the absence of a v-structure and suggests to propagate the orientation of a previously directed edge $X \to Z - Y$ as $X \to Z \to Y$.

In practice, MIIC's strategy to circumvent spurious conditional independences significantly improves recall, that is, the fraction of correctly recovered edges, compared to traditional constraint-based methods [28, 29]. Yet, MIIC only relies on unshielded triple scores to reliably uncover significant contributors 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 [37, 34] and to improve the reliability of predicted orientations [33, 34].

The predictions of this recent version of MIIC, which include three type 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-order local information scores

Starting from the structure predicted by MIIC, as detailed above, MIIC_search&score proceeds in two steps, first to remove likely false positive edges (Step 1) and then to (re)orient the remaining edges based on their estimated contributions to the global likelihood decomposition, Eq. 12 (Step 2).

### 3.2.1 Step 1: Node scores for edge orientation priming and edge removal

The first step consists in minimizing a node score corresponding to the local normalized log likelihood of each node w.r.t. its possible parents or spouses amongst the connected nodes predicted by MIIC. To this end, the node score assesses the conditional entropy of each node w.r.t. a selection of parents, spouses or neighbors, $Pa'_{X_i} \subseteq Pa_{X_i} \cup Sp_{X_i} \cup Ne_{X_i}$, and a factorized Normalized Maximum Likelihood (fNML) regularization [28], see Appendix D for details,

$$\text{Score}_n(X_i) = H(X_i|\text{Pa}'_{X_i}) + \frac{1}{N} \sum_j^{q_{x_i}} \log C_{n_j}^{r_{x_i}} \quad (14)$$

where $q_{x_i}$ corresponds to the combination of levels of $Pa'_{X_i}$, while $r_{x_i}$ is the number of levels of $X_i$, and $n_j$ the number of samples corresponding to a particular combination of levels $j$ in each summand, with $\sum_j n_j = N$, the total number of samples. $\log C_{n_j}^{r_{x_i}}$ is the fNML regularization cost summed over all combinations of levels, $q_{x_i}$, [38, 39], see Appendix D.

This first algorithm is looped over each node, priming the orientations of their surrounding edges (as directed, bidirected 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.

### 3.2.2 Step 2: Edge orientation scores, as local contributions to the global likelihood score

The second step consists in orienting the edges retained after Step 1, based on the optimization of their local contributions 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 a fNML complexity cost, Table 1, given three sets of parents and spouses of X and Y, i.e. $Pa'_{X\times Y} = Pa_X \cup Sp_X \setminus Y$, $Pa'_{Y\times X} = Pa_Y \cup Sp_Y \setminus X$ and $Pa'_{XY} = Pa'_{X\times Y} \cup Pa'_{Y\times X}$ with their corresponding combinations of levels, $q_{y\times x}$, $q_{x\times y}$ and $q_{xy}$. These orientation scores, Table 1, include symmetrized fNML complexity 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'_{Y\times X} = Pa'_{X\times Y} = Pa'_{XY}$ implying also the same combinations of parent and spouse levels, $q_{y\times x} = q_{x\times y} = q_{xy}$.

<table><caption>Table 1: Local scores for the orientation of a single directed or bidirected edge, see Appendix D.</caption><thead><tr><th>Edge</th><th>Information</th><th>Symmetrized fNML complexity (Markov equivalent)</th></tr></thead><tbody><tr><td>X &rarr; Y</td><td>-I(X; Y|Pa'_{Y\times X})</td><td>1/2N (Σ<sup>q<sub>xy</sub>r<sup>y</sup></sup><sub>j</sub> log C<sup>r<sub>x</sub></sup><sub>n<sub>j</sub></sub> - Σ<sup>q<sub>xy</sub>r<sup>y</sup></sup><sub>j</sub> log C<sup>r<sub>x</sub></sup><sub>n<sub>j</sub></sub> + Σ<sup>q<sub>yx</sub>r<sup>x</sup></sup><sub>j</sub> log C<sup>r<sub>y</sub></sup><sub>n<sub>j</sub></sub> - Σ<sup>q<sub>yx</sub>r<sup>x</sup></sup><sub>j</sub> log C<sup>r<sub>y</sub></sup><sub>n<sub>j</sub></sub>)</td></tr><tr><td>X &harr; Y</td><td>-I(X; Y|Pa'_{X\times Y})</td><td>1/2N (Σ<sup>q<sub>xy</sub>r<sup>y</sup></sup><sub>j</sub> log C<sup>r<sub>x</sub></sup><sub>n<sub>j</sub></sub> - Σ<sup>q<sub>xy</sub>r<sup>y</sup></sup><sub>j</sub> log C<sup>r<sub>x</sub></sup><sub>n<sub>j</sub></sub> + Σ<sup>q<sub>yx</sub>r<sup>x</sup></sup><sub>j</sub> log C<sup>r<sub>y</sub></sup><sub>n<sub>j</sub></sub> - Σ<sup>q<sub>yx</sub>r<sup>x</sup></sup><sub>j</sub> log C<sup>r<sub>y</sub></sup><sub>n<sub>j</sub></sub>)</td></tr><tr><td>X &harr; Y</td><td>-I(X; Y|Pa'_{XY})</td><td>1/2N (Σ<sup>q<sub>xy</sub>r<sup>y</sup></sup><sub>j</sub> log C<sup>r<sub>x</sub></sup><sub>n<sub>j</sub></sub> - Σ<sup>q<sub>xy</sub>r<sup>y</sup></sup><sub>j</sub> log C<sup>r<sub>x</sub></sup><sub>n<sub>j</sub></sub> + Σ<sup>q<sub>yx</sub>r<sup>x</sup></sup><sub>j</sub> log C<sup>r<sub>y</sub></sup><sub>n<sub>j</sub></sub> - Σ<sup>q<sub>yx</sub>r<sup>x</sup></sup><sub>j</sub> log C<sup>r<sub>y</sub></sup><sub>n<sub>j</sub></sub>)</td></tr></tbody></table>

While orientation scores cannot be summed over individual edges due to multiple countings of ac-connected contributions, score differences between alternative orientations provide an estimate of the global score change. Hence, step 2 algorithm is looped over each edge to compute an 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 triangular cycles, is then chosen at each iteration until convergence or until a limit cycle is reached. Limit cycles may originate from the local two-collider approximation of the global score.

# 4 Experimental results

We first tested whether MIIC_search&score orientation scores (Table 1) effectively predicts bidirected orientations on three simple ancestral models, Fig. 4 in Appendix E, when the end nodes do not share the same parents (Fig. 4, Model 1), share some parents (Fig. 4, Model 2) or when the bidirected edge is part of a longer than two-collider paths (Fig. 4, Model 3). The prediction of the edge orientation scores are summarized in Table 3, Appendix E, and show good predictions for large enough datasets.

Beyond these simple examples, focussing on the discovery of bidirected edges in small toy models of ancestral graphs, we also analyzed more challenging benchmarks from the bnlearn repository [40], Figs. 2-3. They concern ancestral graphs obtained by hiding up to 20% of variables in Discrete Bayesian Networks of increasing complexity (number of nodes and parameters), such as Alarm (37 nodes, 46 links, 509 parameters), Insurance (27 nodes, 52 links, 1,008 parameters), Barley (48 nodes, 84 links, 114,005 parameters), and Mildew (35 nodes, 46 links, 540,150 parameters). We then assessed causal discovery performance in terms of Precision, TP/(TP + FP), and Recall, TP/(TP + FN), relative to the theoretical PAGs, while counting as false positive (FP), all correctly predicted edges but with a different orientation as the directed or bidirected edges of the PAG.

Figs. 2-3 compare MIIC_search&score performance to MIIC results used as starting point for MIIC_search&score and to FCI [41]. Fig. 2 results are obtained from independent datasets for each ancestral graph and sample size, while Fig. 3 results provide a bootstrap sensitivity analysis to sampling noise for each method based on independent resamplings with replacement of single datasets of increasing size. MIIC and MIIC_search&score settings were set as described in section 3 above. The open-source MIIC R package (v1.5.2, GPL-3.0 license) was obtained at https://github.com/miicTeam/miic_R_package. FCI from the python causal-learn package (v0.1.3.8, MIT license) [41] was obtained at https://github.com/py-why/causal-learn and run with G²-conditional independence test and default parameter α = 0.05.

Overall, MIIC_search&score is found to outperform MIIC in terms of edge precision with little to no decrease in edge recall, Figs. 2-3, demonstrating the benefit of MIIC_search&score's rationale to improve MIIC predictions by extending MIIC information scores from unshielded triples to higher-order information contributions. These originate from ac-connected subsets including nodes with more than two parents or spouses, or ac-connected subsets including two-collider paths.
MIIC_search&score 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. Insurance, Barley and Mildew), while reaching similar performance at larger sample sizes or for simpler graphical models (i.e. similar precision on Alarm model), as expected from the asymptotic consistency of FCI for very large datasets, Fig. 2. 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. 3, except for Alarm on all sample sizes tested (*N* ≤ 20,000) and for Insurance at small sample sizes (*N* ≤ 1,000).

<figure><img src="image_4.png" /><figcaption>**Figure 2: Benchmark results on ancestral graphs of increasing complexity.** Benchmark results are averaged over 50 independent categorical datasets from ancestral graphs obtained by hiding 0%, 5%, 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 to MIIC results used as starting point for MIIC_search&score and FCI [41]. 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 as the directed or bidirected edges of the PAG. Error bars: 95% confidence interval.</figcaption></figure>

Importantly, the benchmark PAGs used to score the causal discovery results with increasing proportions of latent variables, Figs. 2-3, 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 graphical model without hidden variables. This explains the steady decrease in recall for complex ancestral models (i.e. Insurance, Barley, Mildew) with higher proportions of hidden variables, while precision remains essentially unaffected, Figs. 2-3.

All in all, these results highlight MIIC_search&score capacity to efficiently and robustly learn complex graphical models from limited available data, which is a frequent setting for many real-world applications, in practice. In addition, MIIC_search&score, which has been implemented to analyze challenging categorical datasets, is quite unique in this regard, as all other search-and-score methods
for ancestral graphs [11–13, 30, 32] have only been demonstrated on continuous datasets from linear Gaussian models and could not be included in the present benchmarks, Figs. 2-3.

<figure><img src="image_5.png" /><figcaption>Figure 3: Benchmark results on bootstrap datasets from ancestral graphs of increasing complexity. Benchmark results on bootstrap sensitivity analysis to sampling noise based on 30 independent resamplings with replacement of single datasets of increasing sizes. Ancestral graphs are obtained by hiding 0%, 5%, 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 to MIIC results used as starting point for MIIC_search&score and FCI [41]. The lack of FCI results, except for Alarm on all sample sizes tested ($N \le 20,000$) and for Insurance at small sample sizes ($N \le 1,000$), stems from FCI difficulty to converge on bootstrapped datasets. 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 as the directed or bidirected edges of the PAG. Error bars: 95% confidence interval.</figcaption></figure>

# 5 Limitations

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.

While this approach could be extended to higher-order information contributions including three-or-more-collider paths, it allows for a simple two-step search-and-score scheme at the level of individual nodes (step 1) and edges (step 2), as detailed in section 3. This already shows a significant improvement in causal discovery performance (i.e. combining good precision and good recall on challenging benchmarks) as compared to existing state-of-the-art methods.

# Acknowledgments and Disclosure of Funding

HI acknowledges support from CNRS, Institut Curie and Prairie Institute, as well as funding from ANR-22-PESN-0002 “AI4scMed” and ANR 23CE13001802 “Patterning” grants. NL acknowledges a PhD fellowship from CNRS-Imperial College joint PhD programme.
References

[1] D. Koller, N. Friedman, *Probabilistic Graphical Models: Principles and Techniques* (MIT Press, 2009).

[2] J. Pearl, A. Paz, *Graphoids: A graph-based logic for reasoning about relevance relations*, or when would x tell you more about y if you already know z, *Tech. rep.*, UCLA Computer Science Department (1985).

[3] J. Pearl, *Probabilistic reasoning in intelligent systems* (Morgan Kaufmann, San Mateo, CA, 1988).

[4] J. Pearl, *Causality: models, reasoning and inference* (Cambridge University Press, 2009), second edn.

[5] P. Spirtes, C. Glymour, R. Scheines, *Causation, Prediction, and Search* (MIT press, , 2000), second edn.

[6] T. S. Richardson, A factorization criterion for acyclic directed mixed graphs, *Proceedings of the Twenty-Fifth Conference on Uncertainty in Artificial Intelligence*, UAI '09 (AUAI Press, Arlington, VA, USA, 2009), p. 462–470.

[7] J. Tian, J. Pearl, A general identification condition for causal effects, *Proceedings of the National Conference on Artificial Intelligence* (Menlo Park, CA; Cambridge, MA; London; AAAI Press; MIT Press; 1999), pp. 567–573.

[8] T. Richardson, P. Spirtes, Ancestral graph markov models. *Ann. Statist.* **30**, 962–1030 (2002).

[9] M. Drton, M. Eichler, T. S. Richardson, Computing maximum likelihood estimates in recursive linear models with correlated errors. *Journal of Machine Learning Research* **10**, 2329–2348 (2009).

[10] R. J. Evans, T. S. Richardson, Maximum likelihood fitting of acyclic directed mixed graphs to binary data, *Proceedings of the 26th Conference on Uncertainty in Artificial Intelligence*, UAI’10 (AUAI Press, Corvallis, OR, USA, 2010).

[11] S. Triantafillou, I. Tsamardinos, Score-based vs constraint-based causal learning in the presence of confounders, CFA@UAI (2016).

[12] K. Rantanen, A. Hyttinen, M. Järvisalo, Maximal ancestral graph structure learning via exact search, *Proceedings of the Thirty-Seventh Conference on Uncertainty in Artificial Intelligence*, C. de Campos, M. H. Maathuis, eds. (PMLR, 2021), vol. 161 of *Proceedings of Machine Learning Research*, pp. 1237–1247.

[13] T. Claassen, I. G. Bucur, Greedy equivalence search in the presence of latent confounders, *Proceedings of the Thirty-Eighth Conference on Uncertainty in Artificial Intelligence*, J. Cussens, K. Zhang, eds. (PMLR, 2022), vol. 180 of *Proceedings of Machine Learning Research*, pp. 443–452.

[14] Z. Hu, R. Evans, Faster algorithms for markov equivalence, *Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI)*, J. Peters, D. Sontag, eds. (PMLR, 2020), vol. 124 of *Proceedings of Machine Learning Research*, pp. 739–748.

[15] W. J. McGill, Multivariate information transmission. *Trans. of the IRE Professional Group on Information Theory (TIT)* **4**, 93-111 (1954).

[16] H. K. Ting, On the amount of information. *Theory Probab. Appl.* **7**, 439-447 (1962).

[17] T. S. Han, Multiple mutual informations and multiple interactions in frequency data. *Information and Control* **46**, 26-45 (1980).

[18] R. W. Yeung, A new outlook on shannon’s information measures. *IEEE transactions on information theory* **37**, 466–474 (1991).

[19] P. Spirtes, T. Richardson, A polynomial time algorithm for determinint dag equivalence in the presence of latent variables and selection bias, *Proceedings of the 6th International Workshop on Artificial Intelligence and Statistics* (1996).

[20] T. Richardson, Markov properties for acyclic directed mixed graphs. *Scandinavian Journal of Statistics* **30**, 145-157 (2003).

[21] R. A. Ali, T. S. Richardson, Markov equivalence classes for maximal ancestral graphs, *Proceedings of the Eighteenth Conference on Uncertainty in Artificial Intelligence*, UAI'02 (Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 2002), pp. 1–9.

[22] R. A. Ali, T. S. Richardson, P. Spirtes, J. Zhang, Towards characterizing markov equivalence classes for directed acyclic graphs with latent variables, *Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence*, UAI'05 (Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 2005).
[23] J. Tian, Generating markov equivalent maximal ancestral graphs by single edge replacement, *Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence*, UAI'05 (Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 2005).

[24] R. A. Ali, T. S. Richardson, P. Spirtes, Markov equivalence for ancestral graphs. *Ann. Statist.* **37**, 2808–2837 (2009).

[25] T. Richardson, P. Spirtes, Scoring ancestral graph models, *Tech. rep.* (1999). Available as Technical Report CMU-PHIL 98.

[26] J. Zhang, A characterization of markov equivalence classes for directed acyclic graphs with latent variables, *Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence*, UAI'07 (Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 2007).

[27] J. Zhang, On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias. *Artif. Intell.* **172**, 1873-1896 (2008).

[28] S. Affeldt, H. Isambert, Robust reconstruction of causal graphical models based on conditional 2-point and 3-point information, *Proceedings of the Thirty-First Conference on Uncertainty in Artificial Intelligence*, UAI 2015, July 12-16, 2015, Amsterdam, The Netherlands (2015), pp. 42–51.

[29] L. Verny, N. Sella, S. Affeldt, P. P. Singh, H. Isambert, Learning causal networks with latent variables from multivariate information in genomic data. *PLoS Comput. Biol.* **13**, e1005662 (2017).

[30] B. Andrews, G. F. Cooper, T. S. Richardson, P. Spirtes, The m-connecting imset and factorization for admg models, Preprint (2022). Arxiv 2207.08963.

[31] Z. Hu, R. J. Evans, Towards standard imsets for maximal ancestral graphs. *Bernoulli* **30** (2024).

[32] Z. Hu, R. Evans, A fast score-based search algorithm for maximal ancestral graphs using entropy, Preprint (2024). Arxiv 2402.04777.

[33] V. Cabeli, H. Li, M. da Câmara Ribeiro-Dantas, F. Simon, H. Isambert, Reliable causal discovery based on mutual information supremum principle for finite datasets, *WHY21, 35rd Conference on Neural Information Processing Systems* (NeurIPS, 2021).

[34] M. d. C. Ribeiro-Dantas, H. Li, V. Cabeli, L. Dupuis, F. Simon, L. Hettal, A.-S. Hamy, H. Isambert, Learning interpretable causal networks from very large datasets, application to 400, 000 medical records of breast cancer patients. *iScience* **27**, 109736 (2024).

[35] V. Cabeli, L. Verny, N. Sella, G. Uguzzoni, M. Verny, H. Isambert, Learning clinical networks from medical records based on information estimates in mixed-type data. *PLoS Comput. Biol.* **16**, e1007866 (2020).

[36] S. Affeldt, L. Verny, H. Isambert, 3off2: A network reconstruction algorithm based on 2-point and 3-point information statistics. *BMC Bioinformatics* **17** (2016).

[37] H. Li, V. Cabeli, N. Sella, H. Isambert, Constraint-based causal structure learning with consistent separating sets. *Advances in Neural Information Processing Systems (NeurIPS)* **32** (2019).

[38] P. Kontkanen, P. Myllymäki, A linear-time algorithm for computing the multinomial stochastic complexity. *Inf. Process. Lett.* **103**, 227–233 (2007).

[39] T. Roos, T. Silander, P. Kontkanen, P. Myllymäki, Bayesian network structure learning using factorized nml universal models, *Proc. 2008 Information Theory and Applications Workshop (ITA-2008)* (IEEE Press, 2008).

[40] M. Scutari, Learning Bayesian Networks with the bnlearn R Package. *J. Stat. Softw.* **35**, 1–22 (2010).

[41] Y. Zheng, B. Huang, W. Chen, J. Ramsey, M. Gong, R. Cai, S. Shimizu, P. Spirtes, K. Zhang, Causal-learn: Causal discovery in python. *Journal of Machine Learning Research* **25**, 1–8 (2024).

[42] Y. M. Shtarkov, Universal sequential coding of single messages. *Problems of Information Transmission* **23**, 3–17 (1987).

[43] J. Rissanen, I. Tabus, *Adv. Min. Descrip. Length Theory Appl.* (MIT Press, 2005), pp. 245–264.

[44] W. Szpankowski, *Average case analysis of algorithms on sequences* (John Wiley & Sons, 2001).
[45] P. Kontkanen, W. Buntine, P. Myllymäki, J. Rissanen, H. Tirri, Efficient computation of stochastic complexity. in: C. Bishop, B. Frey (Eds.) *Proceedings of the Ninth International Conference on Artificial Intelligence and Statistics*, Society for Artificial Intelligence and Statistics **103**, 233–238 (2003).

[46] P. Kontkanen, Computationally efficient methods for mdl-optimal density estimation and data clustering, Ph.D. thesis (2009).

[47] D. M. Chickering, A Transformational Characterization of Equivalent Bayesian Network Structures, *UAI '95: Proceedings of the Eleventh Annual Conference on Uncertainty in Artificial Intelligence* (Morgan Kaufmann, 1995), pp. 87–98.

[48] M. Kalisch, M. Mächler, D. Colombo, M. H. Maathuis, P. Bühlmann, Causal inference using graphical models with the r package pcalg. *J. Stat. Softw.* **47**, 1–26 (2012).
# Appendix / supplemental material

## 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 independence relations by a graphical criterion, called *m*-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 $\pi$ between $X$ and $Y$ is *m*-connecting given a (possibly empty) subset $C \subseteq V$ (with $X, Y \notin C$) if:

i) its non-collider(s) are not in $C$, and

ii) its collider(s) are in $A\mathbf{n}_G(C)$.

**Definition 5.** [*m*-separation criterion] The subsets $A$ and $B$ are said to be *m*-separated by $C$, noted $A \perp_m B | C$, if there is no *m*-connecting path between any vertex in $A$ and any vertex in $B$ given $C$.

The probabilistic interpretation of ancestral graph is given by its global and pairwise Markov properties (which are equivalent [8]): if $A$ and $B$ are *m*-separated by $C$, then $A$ and $B$ are conditionally independent given $C$ and $\forall X \in A$ and $\forall Y \in B$, there is a probability distribution $P$ faithful to $\mathcal{G}$ such that their conditional mutual information vanishes, i.e. $I_P(X;Y|C) = 0$, also noted $X \perp_P Y | C$.

However, as discussed above, the proof of Theorem 1 will require to introduce a weaker *m'*-connection criterion defined below.

**Definition 6.** [*m'*-connecting path] A path $\pi$ between $X$ and $Y$ is *m'*-connecting given a subset $C \subseteq V$ (with $X, Y$ possibly in $C$) if:

i) its non-collider(s) are not in $C$, and

ii) its collider(s) are in $A\mathbf{n}_G(\{X, Y\} \cup C)$.

Note, in particular, that an *m*-connecting path is necessary an *m'*-connecting path but that the converse is not always true. For example, the path $X \to Z \leftrightarrow T \leftrightarrow Y$ in Fig. 1G (with $Z \to Y$) is an *m'*-connecting path given $T$ (as $Z \in A\mathbf{n}_G(\{X, Y\} \cup T)$) but not an *m*-connecting path given $T$ (as $Z \notin A\mathbf{n}_G(T)$).

However, Richardson and Spirtes 2002 [8] have shown the following lemma,

**Lemma 4.** [Corollary 3.15 in [8]] In an ancestral graph $\mathcal{G}$, there is a *m'*-connecting path $\mu$ between $X$ and $Y$ given $C$ if and only if there is a (possibly different) *m*-connecting path $\pi$ between $X$ and $Y$ given $C$.

Hence, Lemma 4 implies that *m'*-separation and *m*-separation criteria are in fact equivalent, as an absence of *m'*-connecting paths implies an absence of *m*-connecting paths and vice versa. This enables to reformulate the *m*-separation criterion above as,

**Definition 7.** [*m'*-separation (and *m*-separation) criteria] The subsets $A$ and $B$ are said to be *m'*-separated (or *m*-separated) by $C$, if all paths from any $X \in A$ to any $Y \in B$ have either

i) a non-collider in $C$, or

ii) a collider not in $A\mathbf{n}_G(\{X, Y\} \cup C)$.

The probabilistic interpretation of an ancestral graph is given by its (global) Markov property: if $A$ and $B$ are *m*-separated (or *m'*-separated) by $C$, then $A$ and $B$ are conditionally independent given $C$, noted as, $A \perp_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-connecting paths**, which is directly relevant to characterize the likelihood decomposition and Markov equivalent classes of ancestral graphs (Theorem 1). We give here a different yet equivalent definition of **ac-connecting paths** as defined in the main text (Definition 2) in order to underline the similarities and differences with the notion of *m'*-connecting path (Definition 6).
**Definition 8.** [ac-connecting path] A path $\pi$ between $X$ and $Y$ is an ac-connecting path given a subset $C \subseteq V$ (with $X$ and $Y$ possibly in $C$) if:

i) $\pi$ does not have any noncollider, and

ii) its collider(s) are in $\mathbf{An}_G(\{X, Y\} \cup C)$.

Hence, more simply (following Definition 2 in the main text), an ac-connecting path given $C$ is a collider path, $X \dasharrow Z_1 \leftrightarrow \dots \leftrightarrow Z_K \dasharrow Y$, with all $Z_i \in \mathbf{An}_G(\{X, Y\} \cup C)$, i.e. with $Z_i$ in $C$ or connected to $\{X, Y\} \cup C$ by an ancestor path, $Z_i \to \dots \to T$ with $T \in \{X, Y\} \cup C$.

**Definition 9.** [ac-separation criterion] The subsets $A$ and $B$ are said to be ac-separated by $C$ if there is no ac-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) \textit{m-connecting path } \pi \implies \textit{m'-connecting path } \pi$$

$$ii) \textit{ac-connecting path } \pi \implies \textit{m'-connecting path } \pi$$

$$iii) \textit{m-separation} \iff \textit{m'-separation}$$

$$iv) \textit{m/m'-separation} \implies \textit{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 1): A subset $C$ is said to be ac-connected if $\forall X, Y \in 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 $\mathcal{G}$ (Definition 3), we will first show (i) that all non-ac-connected subsets $S'$ are included in a cancelling combination of multivariate cross-information terms, $I(X; Y|\mathbf{A}) = 0$, with $X, Y \in S'$ and $S' \subseteq S = \{X, Y\} \cup \mathbf{A}$, Eq. 11. Conversely, we will then show (ii) that cancelling combinations of multivariate cross-information terms associated to pairwise conditional independence, $I(X; Y|\mathbf{A}) = \sum_{S' \subseteq S}^{X, Y \in S'} (-1)^{|S'|} I(S') = 0$, do not contain any ac-connected subset $S'$. Finally, we will prove (iii) 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 $I(X; Y; Z|\mathbf{A}) = 0$.

i) Let's first prove that all non-ac-connected subsets $S'$ are included in at least one cancelling combination of multivariate cross-information, $I(X; Y|\mathbf{A}) = 0$, with $X, Y \in S'$ and $S' \subseteq \{X, Y\} \cup \mathbf{A}$.

If $S'$ is a non-ac-connected subset, there is at least one disconnected pair $X$ and $Y$ for which each path $\pi_j$ between $X$ and $Y$ contains either some collider(s) not in $\mathbf{An}_G(S')$ or, if all colliders along $\pi_j$ are in $\mathbf{An}_G(S')$, there must be some non-collider(s) at node(s) $Z_j$ but not necessarily in $S'$. Let's define $S = S' \cup_j Z_j$. $X$ and $Y$ can be shown to be m-separated given $S \setminus \{X, Y\}$, as for each path $\pi_j$ between $X$ and $Y$, its non-collider(s) are in $S$ at node(s) $Z_j$ (when all collider(s) along $\pi_j$ are in $S'$) or there is some collider(s) not in $\mathbf{An}_G(S')$, which are not in $\mathbf{An}_G(S)$ either. The latter statement is proven by contradiction assuming that there is a collider at $Z \notin \mathbf{An}_G(S')$ such that $Z \in \mathbf{An}_G(S)$. There is therefore a directed path $Z \to \dots \to W$ with $W \in S$. Hence, $W \in S'$ or there is a noncollider at $W \in Z_j$ which is on a path $\pi_j$ between $X$ and $Y$ along which all colliders are in $\mathbf{An}_G(S')$ by construction of $S$. This leads by induction to $Z \to \dots \to W \to \dots \to T$ where $T \in S'$ and thus $Z \in \mathbf{An}_G(S')$, which is a contradiction. Hence, all non-ac-connected subsets $S'$ are included in a cancelling combination of multivariate cross-information terms, $I(X; Y|\mathbf{A}) = 0$, with $X, Y \in S'$ and $S' \subseteq S = \{X, Y\} \cup \mathbf{A}$.

ii) Conversely, we will now show that cancelling combinations of multivariate cross-information terms associated to pairwise conditional independence, $I(X; Y|\mathbf{A}) = \sum_{S' \subseteq S}^{X, Y \in S'} (-1)^{|S'|} I(S') = 0$, do not contain any ac-connected subset $S'$, where $S = \{X, Y\} \cup \mathbf{A}$.
We will prove it by contradiction assuming that there exists a subset $W \subseteq A$, such that $S' = \{X, Y\} \cup W$ is ac-connected. In particular, there should be an ac-connecting path between X and Y confined to $An_G(S')$ and thus to $An_G(S) \supset An_G(S')$, which is an $m'$-connecting path between X and Y given $A$, contradicting the above hypothesis of $m'$-separation given $A$, i.e. $I(X; Y|A) = 0$. The use of $m'$-separation, i.e. the absence of $m'$-connecting paths with colliders in $An_G(S)$ rather than $m$-connecting paths with colliders in $An_G(A)$, is necessary here, see Definitions 4 and 6. Hence, no ac-connected subset $S'$ is included in cancelling combinations of multivariate cross-information terms associated to pairwise conditional independence, $I(X; Y|A) = \sum_{S' \subseteq S}^{X,Y \in S'} (-1)^{|S'|} I(S') = 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, as these shared/overlapping terms in fact all cancel through more global Markov independence relationships involving higher order (three or more points) vanishing multivariate cross-information terms, such as $I(X; Y; Z|A) = 0$.

This result requires to use an ordering of the nodes, $X_k \succ X_j \succ X_i$, that is compatible with the directed edges of the ancestral graph assumed to have no undirected edges, i.e. $X_j \notin An(X_i)$ if $X_j \succ X_i$. Under this ordering, higher order nodes $X_k \succ X_j \succ X_i$ can be a priori excluded from all separating sets $A_{ij}$ of pairs of lower order nodes, i.e. if $I(X_i; X_j|A_{ij}) = 0$ then $X_k \notin A_{ij}$.

In particular, the two pairwise conditional independence relations $I(X_k; X_\ell|A_{k\ell}) = 0$, with $X_\ell \succ X_k$, and $I(X_i; X_j|A_{ij}) = 0$, with $X_j \succ X_i$, do not share any multivariate cross-information terms, if $X_\ell \neq X_j$. Indeed, as $I(X_i; X_j|A_{ij})$ contains all cross-information terms including both $X_i$ and $X_j$ as well as every subset (possibly empty) of $A_{ij}$, none of them includes $X_\ell$ if $X_\ell \succ X_j$. Therefore $I(X_i; X_j|A_{ij})$ does not contain any cross-information term of $I(X_k; X_\ell|A_{k\ell})$ which contains both $X_k$ and $X_\ell$ as well as every subset (possibly empty) of $A_{k\ell}$. This property eliminates all multiple countings of multivariate cross-information terms if $X_\ell \neq X_j$. 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_i, X_r\}$, $\{X_j, X_r\}$, ..., $\{X_\ell, X_r\}$, sharing the same higher order node, $X_r$, i.e., $I(X_k; X_r|A_{kr}) = 0$ for $k = i, j, \dots, \ell$. However, this situation also entails a more global Markov independence constraint between $X_r$ and $\{X_i, X_j, \dots, X_\ell\}$, given a separating set $A$, with $A_{kr} \subseteq A \cup \{X_i, \dots, X_\ell\}$ for $k = i, j, \dots, \ell$. Such a global Markov independence constraint can be decomposed into more local independence constraints using the chain rule (in any order of the variables $X_i, X_j, \dots, X_\ell$) and the decomposition rules of multivariate (cross) information (Eq. 9),

$$ 
\begin{align*}
0 &= I(\{X_i, X_j, \dots, X_\ell\}; X_r|\mathbf{A}) \\
&= (I(X_i; X_r|\mathbf{A}) + I(X_j; X_r|\mathbf{A}, X_i)) + [I(X_k; X_r|\mathbf{A}, X_i, X_j)] + \dots + I(X_\ell; X_r|\mathbf{A}, \dots) \\
&= (I(X_i; X_r|\mathbf{A}) + I(X_j; X_r|\mathbf{A}) - I(X_i; X_j; X_r|\mathbf{A})) \\
&\quad + [I(X_k; X_r|\mathbf{A}, X_i) - I(X_j; X_k; X_r|\mathbf{A}, X_i)] + \dots + I(X_\ell; X_r|\mathbf{A}, \dots) \\
&= (I(X_i; X_r|\mathbf{A}) + I(X_j; X_r|\mathbf{A}) - I(X_i; X_j; X_r|\mathbf{A})) \\
&\quad + [I(X_k; X_r|\mathbf{A}) - I(X_j; X_k; X_r|\mathbf{A}) - I(X_i; X_k; X_r|\mathbf{A}) + I(X_i; X_j; X_k; X_r|\mathbf{A})] + \dots
\end{align*}
$$

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 cross-information terms in the last expression, i.e. between $X_r$ and each subset of $\{X_i, X_j, \dots, X_\ell\}$ given the separating set $A$, all vanish. This result can be readily extended to any subsets $\{X_r, X_s, \dots, X_z\}$ (conditionally) independent of $\{X_i, X_j, \dots, X_\ell\}$ given a separating set $A$, i.e. $I(\{X_i, X_j, \dots, X_\ell\}; \{X_r, X_s, \dots, X_z\}|\mathbf{A}) = 0$. 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 a global cancellation of non-ac-connected 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. $\square$
# C Factorization of the probability distribution of ancestral graphs

## C.1 Factorization resulting from Theorem 1 and Proposition 3

Before presenting the factorization of the model distribution of ancestral graphs resulting from Theorem 1 and Proposition 3, it is instructive to obtain an equivalent factorization for Bayesian graphs, assuming a positive empirical distributions, $p(x_1, \dots, x_m) = \prod_{i=1}^{m} p(x_i|x_{i-1}, \dots, x_1) > 0$,

$$ 
\begin{aligned}
q(x_1, \dots, x_m) &= \prod_{i=1}^{m} q(x_i | \mathbf{pa}_{x_i}) = \prod_{i=1}^{m} p(x_i | \mathbf{pa}_{x_i}) \\
&= p(x_1, \dots, x_m) \prod_{i=1}^{m} \frac{p(x_i | \mathbf{pa}_{x_i})}{p(x_i | x_{i-1}, \dots, x_1)} \\
&= p(x_1, \dots, x_m) \prod_{i=1}^{m} \frac{p(x_i | \mathbf{pa}_{x_i}) p(\boldsymbol{x}_{i-1} | \mathbf{pa}_{x_i})}{p(x_i, \boldsymbol{x}_{i-1} | \mathbf{pa}_{x_i}) p(\mathbf{pa}_{x_i})} \quad (15)
\end{aligned}
$$

This leads to the following alternative expressions for the cross-entropy $H(p, q) = - \sum_{x} p(x) \log q(x)$ in terms of multivariate entropy and information, which only depend on the empirical joint distribution $p(x)$,

$$ 
\begin{aligned}
H(p, q) &= \sum_{i=1}^{m} H(x_i | \mathbf{Pa}_{X_i}) \\
&= H(X_1, \dots, X_m) + \sum_{i=1}^{m} I(X_i; X_{i-1} | \mathbf{Pa}_{X_i}) | \mathbf{Pa}_{X_i}) \quad (16)
\end{aligned}
$$

where $\sum_{i=1}^{m} I(X_i; X_{i-1} | \mathbf{Pa}_{X_i} | \mathbf{Pa}_{X_i})$ can be decomposed, using the chain rule and Eq. 11, 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(x)$, 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 are positive ($\forall \boldsymbol{x}, p(\boldsymbol{x}) > 0$, $q(\boldsymbol{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,

$$ q(\boldsymbol{x}) = \prod_{i=1}^{m} q(x_i) \times \prod_{i<j}^{m} \frac{q(x_i, x_j)}{q(x_i)q(x_j)} \times \prod_{i<j<k}^{m} \frac{q(x_i, x_j, x_k)q(x_i)q(x_j)q(x_k)}{q(x_i, x_j)q(x_i)q(x_j)q(x_k)} \times \dots \quad (17) $$

where all the marginal distributions over a subset of variables, e.g. $q(x_i, x_j, x_k) = \sum_{\ell \neq i, j, k} q(\boldsymbol{x})$ or $p(x_i, x_j, x_k) = \sum_{\ell \neq i, j, k} p(\boldsymbol{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 \to Z \leftrightarrow Y \leftarrow W$ (Fig. 1D), valid for $q(.)$ and $p(.)$ alike,

$$ 
\begin{aligned}
q(x, z, y, w) &= q(x) q(z) q(y) q(w) \\
&\quad \times \frac{q(x, z)}{q(x) q(z)} \frac{q(z, y)}{q(z) q(y)} \frac{q(y, w)}{q(y) q(w)} \frac{q(x, y)}{q(x) q(y)} \frac{q(x, w)}{q(x) q(w)} \frac{q(z, w)}{q(z) q(w)} \\
&\quad \times \frac{q(x) q(z) q(y) q(x, z, y)}{q(x, z) q(x, y) q(z, y)} \frac{q(z) q(y) q(w) q(z, y, w)}{q(z, y) q(z, w) q(y, w)} \\
&\quad \times \frac{q(x) q(z) q(w) q(x, z, w)}{q(x, z) q(x, w) q(z, w)} \frac{q(x) q(y) q(w) q(x, y, w)}{q(x, y) q(x, w) q(y, w)} \\
&\quad \times \frac{q(x, z) q(z, y) q(y, w) q(x, y) q(x, w) q(z, w) q(x, z, y, w)}{q(x, z, y) q(x, z, w) q(x, y, w) q(z, y, w) q(x) q(y) q(z) q(w)} \quad (18)
\end{aligned}
$$

where all individual distribution marginals on subsets of variables, e.g. $q(x)$, $q(x, z)$, $q(x, z, y)$ (or $p(x)$, $p(x, z)$, $p(x, z, y)$), cancel two-by-two by construction, except $q(x, z, y, w)$ (or $p(x, z, y, w)$).
In addition and only for the model distribution $q(\cdot)$, all ratios in gray in Eq. 18 also cancel due to Markov independence relations across non-ac-connected subsets (see proof of Theorem 1). This leaves a 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 $q(\cdot)$ terms by their empirical counterparts $p(\cdot)$, see Proposition 3.

This leads to the following global factorization for $q(\cdot)$ in terms of $p(\cdot)$,

$$ 
\begin{aligned}
q(x, z, y, w) &\equiv p(x) p(z) p(y) p(w) \frac{p(x, z)}{p(x) p(z)} \frac{p(z, y)}{p(z) p(y)} \frac{p(y, w)}{p(y) p(w)} \\
&\quad \times \frac{p(x) p(z) p(y) p(x, z, y)}{p(x, z) p(x, y) p(z, y)} \frac{p(z) p(y) p(w) p(z, y, w)}{p(z, y) p(z, w) p(y, w)} \\
&\quad \times \frac{p(x, z) p(z, y) p(y, w) p(x, y) p(x, w) p(z, w) p(x, z, y, w)}{p(x, z, y) p(x, z, w) p(x, y, w) p(z, y, w) p(x) p(y) p(z) p(w)} \\
&= p(x, z, y, w) \frac{p(x) p(y)}{p(x, y)} \frac{p(x) p(w)}{p(x, w)} \frac{p(z) p(w)}{p(z, w)} \\
&\quad \times \frac{p(x, z) p(x, w) p(z, w)}{p(x) p(z) p(w) p(x, z, w)} \frac{p(x, y) p(x, w) p(y, w)}{p(x) p(y) p(w) p(x, y, w)}
\end{aligned}
\quad (19) $$

where the terms in gray have been passed to the rhs of Eq. 18 applied to $p(\cdot)$. This ultimately leads to the analog of the Bayesian Network factorization in Eq. 15 but for the two-collider path, $X \to Z \leftrightarrow Y \leftarrow W$ (Fig. 1D),

$$ q(x, z, y, w) \equiv p(x, z, y, w) \frac{p(x) p(w)}{p(x, w)} \frac{p(z|x) p(w|x)}{p(z, w|x)} \frac{p(x|w) p(y|w)}{p(x, y|w)} \quad (20) $$

where the last three factors “correct” the expression of $p(x, z, y, w)$ for the three (conditional) independences entailed by the underlying graph, that is, $X \perp W, Z \perp W|X$, and $X \perp Y|W$.

## 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 [6], do not correspond to a single factorized equation (as with Bayesian graphs, Eq. 15) 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 \to Z \leftrightarrow Y \leftarrow W$, introduced above, Fig. 1D, correspond to the following equations [6],

$$ 
\begin{aligned}
q(x, w) &= q(x) q(w) \\
q(x, z) &= q(z|x) q(x) \\
q(y, w) &= q(y|w) q(w) \\
q(x, z, w) &= q(z|x) q(x) q(w) \\
q(x, y, w) &= q(y|w) q(w) q(x) \\
q(x, z, y, w) &= q(z, y|x, w) q(x) q(w)
\end{aligned}
\quad (21) $$

Importantly, these head-and-tail factorizations imply additional relations such as $q(y|w) = q(y|x, w)$ (i.e. $X \perp Y|W$) obtained by comparing the last two relations in Eqs. 21 after marginalizing $q(x, z, y, w)$ over $z$. However, such implicit conditional independence relations are not verified by the empirical distribution $p(\cdot)$ in general and prevent the estimation of the head-and-tail factorizations by substituting the rhs $q(\cdot)$ terms in Eqs. 21 with their empirical counterparts $p(\cdot)$, as in the case of Bayesian networks, Eq. 15.

Indeed, while the head-and-tail factorization relations, Eqs. 21, obey the local and global Markov independence relations entailed by the graphical model, Fig. 1D, leading to the cancellation of all factors associated to non-ac-connected subsets in gray in Eq. 18, the remaining head-and-tail factors cannot be readily estimated with the empirical distribution $p(\cdot)$.
In particular, the cross-entropy of the two-collider path of interest, Fig. 1D, obtained with the head-and-tail factorizations corresponds to $^2H(p, q) = -\sum p(x, z, y, w) \log q(z, y|x, w) q(x) q(w)$. Then, estimating the $q(.)$ terms with their $p(.)$ counterparts leads to the cross-entropy of a Bayesian graph, Fig. 1E, with a different Markov equivalent class than the ancestral graph of interest, Fig. 1D. A similar discrepancy is obtained with a c-component factorization which leads to the cross-entropy of the Bayesian graph of Fig. 1E without edge $X \to Y$, corresponding to a different Markov equivalence class than the previous two graphs, Figs. 1D & E.

These examples illustrate the difficulty to exploit the c-component or head-and-tail factorizations to estimate the likelihood of ancestral graphs including bidirected edge(s).

## D Node and edge scores based on Normalized Maximum Likelihood criteria

Search-and-score 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 in the likelihood decomposition Eq. 12,

$$I(C) \to I'(C) = I(C) - \frac{1}{2} \prod_{k=1}^{|C|} (1 - r_k) \frac{\log N}{N} \quad (22)$$

for categorical datasets, where $r_k$ 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, in general. In order to better take into account finite sample size, we used instead the (universal) Normalized Maximum Likelihood (NML) criterion [42, 43, 38, 39], which amounts to normalizing the likelihood function over all possible datasets with the same number $N$ of samples.

Moreover, as search-and-score structure learning methods need to rely on heuristic rather than exhaustive search, we have implemented a 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 unshielded triple scores to higher-order local information scores including ac-connected subsets of vertices with a maximum of two-collider paths.

**Node score.** We first used the factorized Normalized Maximum Likelihood (fNML) complexity [38, 39] to define a local score for each node $X_i$, which extends the decomposable likelihood of Bayesian graphs given each node's parents, Eq. 2, to all non-descendant neighbors, $Pa'_{X_i}$,

$$\mathcal{L}_{\mathcal{D}|G_{X_i}} = e^{-N \cdot \text{Score}_n(X_i)} = \frac{e^{-NH(X_i|Pa'_{X_i})}}{\sum_{|D'|=N} e^{-NH(X_i|Pa'_{X_i})}} \quad (23)$$

$$= e^{-NH(X_i|Pa'_{X_i})} - \sum_{j}^{q_i} \log C_{n_j}^{r_i} \quad (24)$$

$$= e^{N \sum_{j}^{q_i} \sum_{k}^{r_i} \frac{n_{jk}}{N} \log \left( \frac{n_{jk}}{n_j} \right)} - \sum_{j}^{q_i} \log C_{n_j}^{r_i} \quad (25)$$

$$= \prod_{j}^{q_i} \frac{\prod_{k}^{r_i} \left( \frac{n_{jk}}{n_j} \right)^{n_{jk}}}{C_{n_j}^{r_i}} \quad (26)$$

where $n_{jk}$ corresponds to the number of data points for which $X_i$ is in its kth state and its non- descendant neighbors in their jth state, with $n_j = \sum_k^{r_i} n_{jk}$. The universal normalization constant $C_n^r$ is then computed by summing the numerator over all possible partitions of the n data points into a maximum of r subsets, $\ell_1 + \ell_2 + \dots + \ell_r = n$ with $\ell_k \ge 0$,

$$C_n^r = \sum_{\ell_1+\ell_2+\dots+\ell_r=n} \frac{n!}{\ell_1!\ell_2!\dots\ell_r!} \prod_{k=1}^{r} \left(\frac{\ell_k}{n}\right)^{\ell_k} \quad (27)$$

<footnote n="2">²Indeed, all terms in Eq. 18 actually cancel two-by-two by construction, whatever their factorized expression, except for the remaining joint-distribution over all variables, $q(x, z, y, w) = q(z, y|x, w) q(x) q(w)$.</footnote>
which can in fact be computed in linear-time using the following recursion [38],

$$C_n^r = C_n^{r-1} + \frac{n}{r-2} C_n^{r-2} \qquad (28)$$

with $C_n^1 = 1$ for all $n$ and applying Eq. 31 below for $r = 2$. However, for large $n$ and $r$, $C_n^r$ computation tends to be numerically unstable, which can be circumvented by implementing the recursion on parametric complexity ratios $D_n^r = C_n^r/C_n^{r-1}$ rather than parametric complexities themselves [35] as,

$$D_n^r = 1 + \frac{n}{(r-2)D_n^{r-1}} \qquad (29)$$

$$\log C_n^r = \sum_{k=2}^{r} \log D_n^k \qquad (30)$$

for $r \ge 3$, with $C_n^1 = 1$ and $C_n^2 = D_n^2$, which can be computed directly with the general formula, Eq. 27, for $r = 2$,

$$C_n^2 = \sum_{h=0}^{n} \binom{n}{h} \left(\frac{h}{n}\right)^h \left(\frac{n-h}{n}\right)^{n-h} \qquad (31)$$

or its Szpankowski approximation for large $n$ (needed for $n > 1000$ in practice) [44–46],

$$C_n^2 = \sqrt{\frac{n\pi}{2}} \left( 1 + \frac{2}{3}\sqrt{\frac{2}{n\pi}} + \frac{1}{12n} + \mathcal{O}\left(\frac{1}{n^{3/2}}\right) \right) \qquad (32)$$

$$\approx \sqrt{\frac{n\pi}{2}} \exp \left( \sqrt{\frac{8}{9n\pi}} + \frac{3\pi - 16}{36n\pi} \right) \qquad (33)$$

This leads to the following local score for each node $X_i$, which is minimized over alternative combinations of non-descendant neighbors, $\mathbf{Pa}'_{X_i} \subseteq \mathbf{Pa}_{X_i} \cup \mathbf{Sp}_{X_i} \cup \mathbf{Ne}_{X_i}$, in the first step of the local search-and-score algorithm (step 1) detailed in the main text,

$$\text{Score}_n(X_i) = H(X_i|\mathbf{Pa}'_{X_i}) + \frac{1}{N} \sum_{j}^{q_{X_i}} \log C_{n_j}^{r_{X_i}} \qquad (34)$$

**Edge scores.** We then defined several edge scores to optimize the orientation of each edge, $X-Y$, given its close surrounding vertices.

To this end, we first introduced a local score for node pairs which simply sums the node scores, Eq. 34, for each node. The resulting pair scores are listed in Table 2 for unconnected node pairs and for pairs of nodes connected by a directed edge, where $\mathbf{Pa}'_{X\setminus Y} = \mathbf{Pa}_X \cup \mathbf{Sp}_X \setminus Y$ and $\mathbf{Pa}'_{Y\setminus X} = \mathbf{Pa}_Y \cup \mathbf{Sp}_Y \setminus X$ with their corresponding combinations of levels, $q_{Y\setminus X}$ and $q_{X\setminus Y}$.

<table><caption>Table 2: Local scores for node pairs</caption><thead><tr><th>Pair score</th><th>Information</th><th>fNML Complexity</th></tr></thead><tbody><tr><td>X ≠ Y</td><td>H(X|<b>Pa'</b><sub>X\Y</sub>) + H(Y|<b>Pa'</b><sub>Y\X</sub>)</td><td><sup>1</sup>/<sub>N</sub> ( ∑<sub>j</sub><sup>q<sub>X\Y</sub></sup> log C<sub>n<sub>j</sub></sub><sup>r<sub>X</sub></sup> + ∑<sub>j</sub><sup>q<sub>Y\X</sub></sup> log C<sub>n<sub>j</sub></sub><sup>r<sub>Y</sub></sup>)</td></tr><tr><td>X → Y</td><td>H(X|<b>Pa'</b><sub>X\Y</sub>) + H(Y|<b>Pa'</b><sub>Y\X</sub>, X)</td><td><sup>1</sup>/<sub>N</sub> ( ∑<sub>j</sub><sup>q<sub>X\Y</sub></sup> log C<sub>n<sub>j</sub></sub><sup>r<sub>X</sub></sup> + ∑<sub>j</sub><sup>q<sub>Y\X</sub></sup> log C<sub>n<sub>j</sub></sub><sup>r<sub>Y</sub></sup>)</td></tr><tr><td>X ← Y</td><td>H(X|<b>Pa'</b><sub>X\Y</sub>, Y) + H(Y|<b>Pa'</b><sub>Y\X</sub>)</td><td><sup>1</sup>/<sub>N</sub> ( ∑<sub>j</sub><sup>q<sub>X\Y</sub></sup> r<sub>Y</sub> log C<sub>n<sub>j</sub></sub><sup>r<sub>X</sub></sup> + ∑<sub>j</sub><sup>q<sub>Y\X</sub></sup> log C<sub>n<sub>j</sub></sub><sup>r<sub>Y</sub></sup>)</td></tr></tbody></table>

Then, edge scores for directed edges, $X \to Y$ and $Y \to X$, are defined w.r.t. to the edge removal score, $X \neq Y$, by subtracting the pair scores of unconnected pairs to the pair scores of directed
edges, leading to the following edge orientation scores,

$$ \text{Score}(X \to Y) = -I(X; Y | \mathbf{Pa}'_{Y \lor X}) + \frac{1}{N} \left( \sum_{j}^{q_{y|x} r_x} \log C_{n_j}^{r_y} - \sum_{j}^{q_{y|x}} \log C_{n_j}^{r_y} \right) \quad (35) $$

$$ \text{Score}(Y \to X) = -I(X; Y | \mathbf{Pa}'_{X \lor Y}) + \frac{1}{N} \left( \sum_{j}^{q_{x|Y} r_y} \log C_{n_j}^{r_x} - \sum_{j}^{q_{x|Y}} \log C_{n_j}^{r_x} \right) \quad (36) $$

However, if $r_x \neq r_y$, the fNML complexities of these orientation scores are not identical for Markov equivalent edge orientations between nodes sharing the same parents (or spouses) [47], $\mathbf{Pa}'_{Y \lor X} = \mathbf{Pa}'_{X \lor Y} = \mathbf{Pa}'$ and $q_{y|x} = q_{x|Y}$, despite sharing the same conditional mutual information,

$$ I(X; Y | \mathbf{Pa}') = \frac{1}{2} (H(X | \mathbf{Pa}') + H(Y | \mathbf{Pa}', X)) + \frac{1}{2} (H(X | \mathbf{Pa}', Y) + H(Y | \mathbf{Pa}')) \quad (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 1 in the main text.

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 or spouses different from Y (i.e. corresponding to the information contributions $H(X | \mathbf{Pa}'_{X \lor Y})$) or Y and any combination of parents or spouses different from X (i.e. corresponding to the information contributions $H(Y | \mathbf{Pa}'_{Y \lor X})$) or, else, including both nodes X and Y plus any combination of their parents or spouses, corresponding to the following information contribution, $-I(X; Y | \mathbf{Pa}'_{X \lor Y})$, where $\mathbf{Pa}'_{X \lor Y} = \mathbf{Pa}'_{X \lor Y} \cup \mathbf{Pa}'_{Y \lor X}$. This last term, $-I(X; Y | \mathbf{Pa}'_{X \lor Y})$, contains all the remaining information contributions once the bidirected orientation score is given relative to the edge removal score (Table 2) as for the two directed orientation scores, above. Finally, the symmetrized fNML complexity associated with a bidirected edge should be computed with the whole set of conditioning parents or spouses, $\mathbf{Pa}'_{X \lor Y}$, as indicated in Table 1. 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. $\mathbf{Pa}'_{X \lor Y} = \mathbf{Pa}'_{X \lor Y} = \mathbf{Pa}'_{X \lor Y}$ and $q_{x|Y} = q_{y|x} = q_{x|Y}$ in Table 1.

## E Toy models

Fig. 4 shows three simple ancestral models 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 paths (Model 3).

<figure><img src="image_6.png" /><figcaption>Figure 4: Simple ancestral graphs.</figcaption></figure>

The data is generated from the theoretical DAG using the rmvDAG function in the pcalg package [48]. Each node follows a normal distribution, and the data is discretized using bnlearn's discretize
function using Hartemink's pairwise mutual information method [40]. For these toy models, the edge orientation scores are computed assuming the correct parents of each node.

The prediction of the edge orientation scores are summarized in Table 3 in % of replicates displaying directed edges (wrong) or bidirected edge (correct) as a function of increasing dataset size N.

<table><caption>Table 3: Model 1, X<sub>2</sub> &minus; X<sub>4</sub> Model 2, X<sub>2</sub> &minus; X<sub>4</sub> Model 3, X<sub>2</sub> &minus; X<sub>4</sub> Model 3, X<sub>4</sub> &minus; X<sub>6</sub></caption><thead><tr><th>N</th><th>&leftarrow;</th><th>&rightarrow;</th><th>&leftrightarrow;</th><th>&leftarrow;</th><th>&rightarrow;</th><th>&leftrightarrow;</th><th>&leftarrow;</th><th>&rightarrow;</th><th>&leftrightarrow;</th><th>&leftarrow;</th><th>&rightarrow;</th><th>&leftrightarrow;</th></tr></thead><tbody><tr><td>1000</td><td>0</td><td>100</td><td>0</td><td>50</td><td>42</td><td>8</td><td>2</td><td>98</td><td>0</td><td>94</td><td>6</td><td>0</td></tr><tr><td>5000</td><td>0</td><td>68</td><td>32</td><td>18</td><td>2</td><td>80</td><td>22</td><td>78</td><td>0</td><td>62</td><td>38</td><td>0</td></tr><tr><td>10000</td><td>0</td><td>10</td><td>90</td><td>0</td><td>0</td><td>100</td><td>12</td><td>70</td><td>18</td><td>40</td><td>16</td><td>44</td></tr><tr><td>20000</td><td>0</td><td>0</td><td>100</td><td>0</td><td>0</td><td>100</td><td>0</td><td>0</td><td>100</td><td>2</td><td>0</td><td>98</td></tr><tr><td>35000</td><td>0</td><td>0</td><td>100</td><td>0</td><td>0</td><td>100</td><td>0</td><td>0</td><td>100</td><td>0</td><td>0</td><td>100</td></tr><tr><td>50000</td><td>0</td><td>0</td><td>100</td><td>0</td><td>0</td><td>100</td><td>0</td><td>0</td><td>100</td><td>0</td><td>0</td><td>100</td></tr></tbody></table>