arXiv:1802.01239v1 [cs.DS] 5 Feb 2018

Counting and Uniform Sampling from Markov
Equivalent DAGs
AmirEmad Ghassami,1 Saber Salehkaleybar,2 Negar Kiyavash1
1

Department of ECE, Coordinated Science Laboratory,

University of Illinois at Urbana-Champaign, Urbana, IL, USA
2

Electrical Engineering Department, Sharif University of Technology, Tehran, Iran
1

{ghassam2,kiyavash}@illinois.edu, 2 saleh@sharif.edu

We propose an exact solution for the problem of finding the size of a Markov
equivalence class (MEC). For the bounded degree graphs, the proposed solution is capable of computing the size of the MEC in polynomial time. Our
proposed approach is based on a recursive method for counting the number
of the elements of the MEC when a specific vertex is set as the source variable.
We will further use the idea to design a sampler, which is capable of sampling
from an MEC uniformly in polynomial time.

1

Introduction

Directed acyclic graphs (DAGs) are the most commonly used graphical model to represent causal
relationships among a set of variables. In a DAG representation, a directed edge indicates a direct
causal relationship between the corresponding variables. Under Markov property and faithfulness
assumptions, conditional d-separation of variables in a DAG is in bijective correspondence with
conditional independencies of the variables in the underlying joint probability distribution (Spirtes
et al., 2000), and hence, a DAG representation demonstrates conditional independencies among
its variables.
The general approach for learning a causal structure is to use statistical data from the variables
to find a DAG which is the most consistent with the conditional independencies in the given
data. However, a DAG representation of a set of conditional independencies is not always unique.
1

This restricts the learning of the causal structure to Markov equivalence classes (MECs), where
elements of each class represent the same set of conditional independencies.
There is no preference between the elements of an MEC. Therefore, for a given set of data,
or more generally, for a given joint probability distribution on the variables, the size of the MEC,
i.e., the number of its elements, is a clear metric for measuring the success of an observational
structure learning approach, and indicates the necessity of using interventional experiments. This
indicates the importance of awareness about the size of the MEC for a given structure of interest.
Moreover, as mentioned earlier, one can compare different DAGs to find the most consistent
one with a given joint distribution over the variables. Here, the space of DAGs , or the space of
MECs can be chosen as the search space. Chickering showed that when the size of the MECs
are large, searching over the space of MECs is more efficient (Chickering, 2002). This again
indicates the importance of learning the size of MECs.
To the best of our knowledge, the main existing solution to the problem of finding the size of
MECs is to use Markov chain methods (He et al., 2013, Bernstein and Tetali, 2017). According
to this method, a Markov chain is constructed over the elements of the MEC whose properties
ensure that the stationary distribution is uniform over all the elements. The main problems with
the Markov chain method are its non-negligible time required for convergence and computational
issues, which in some cases prevents its use in practice.
In this paper, we propose an exact solution for the problem of finding the size of an MEC.
For the bounded degree graphs, the proposed solution is capable of computing the size of MEC
in polynomial time. Our proposed approach is based on a recursive method for counting the
number of the elements of the MEC when a specific vertex being set as the source variable. As
mentioned earlier, in an observational setting, there is no preference between the elements of an
MEC. Therefore, a reasonable approach for testing the performance of an algorithm on an MEC
is to test it on uniformly drawn samples from the MEC. This indicates the importance of the
ability of sampling from an MEC uniformly. We will use the introduced idea to design a sampler,
which is capable of sampling from an MEC uniformly in polynomial time. This approach is
applied to finding the best set of variables to intervene on when we are restricted to a certain
budget for the number of experiments (Ghassami et al., 2017).
The rest of the paper is organized as follows. In Section 2, we present the problem description.
The method for counting the elements of an MEC is presented in Section 3. The counting idea is
then used in Section 4 to design a uniform sampler from Markov equivalent DAGs.

2

2

Problem Description

An MEC is usually represented by a mixed graph M = (V, A, E), which is a graph containing
both directed and undirected edges. V , A, and E denote the set of vertices, the set of directed
edges and the set of undirected edges, respectively. The directed edges are those that have the
same direction in all elements of the class, that is, their directions are essential. This is the reason
that these mixed graphs are referred to as essential graph in the literature (Andersson et al., 1997).
For the rest of the edges, which are undirected, there exist two elements of the class which do
not agree on the direction of that edge. Essential graphs are also called complete partial directed
acyclic graphs (PDAGs) (Chickering, 2002), and maximally oriented graphs (Meek, 1995).
We define the skeleton of a DAG as an undirected graph on the same set of vertices as the
original DAG, and with undirected version of the edges of the DAG. We also define v-structure as
a structure comprising two converging directed edges whose tails are not connected by an edge.
The following result due to Verma and Pearl (Verma and Pearl, 1991) characterizes Markov
equivalence.
Proposition 1 (Verma and Pearl, 1991). Two DAGs are Markov equivalent if and only if they
have the same skeleton and the same v-structures.
Proposition 1 indicates that all the edges which are a part of a v-structure are directed in
the essential graph, i.e., in none of the members of an MEC, the subgraph corresponding to the
undirected part of the essential graph is allowed to contain a v-structure.
Let M be an essential graph. Consider the undirected graph resulted from removing all the
directed edges from M . This graph consists of disconnected chain components {G1 , ..., Gc }.
Note that each of these components will be an essential graph themselves. When the entire
essential graph is undirected, it forms a single chain component. Andersson et al. showed that
each of the chain components are chordal graphs (Andersson et al., 1997). Let Size(M ) denote
the size of the MEC corresponding to M . Gillispie and Perlman showed that Size(M ) can be
calculated as follows (Gillispie and Perlman, 2002).
c
Y
Size(M ) =
Size(Gi ).
i=1

This result suggests that one can study each chain component separately.
Therefore, we will study the problem of finding the size of an MEC for which the essential
graph is a connected undirected chordal graph, and hence, none of the members of the MEC are
allowed to have any v-structures.
3

3

Counting the Members of an MEC

The goal in this section is to count the members of an MEC, for which the essential graph is a
connected undirected chordal graph G = (V (G), E(G)).
Define a source variable in a digraph, as one for which the incoming degree is zero. Also, let
dG (v, u) be the distance between vertices v and u on G. The following result from (Bernstein
and Tetali, 2017) is used in our proposed method.
Proposition 2 (Bernstein and Tetali, 2017). Any acyclic and v-structure-free orientation of a
connected graph G has a unique source variable. The unique source v, determines the orientation
of all edges {u, w}, for which dG (v, u) 6= dG (v, w).
Proposition 2 determines that the direction of some edges will be resolved when the source
vertex is known. These are the edges with one endpoint closer to the source vertex than the
other. We note that removing the resolved edges, the obtained digraph will remain acyclic and
v-structure-free:
Lemma 1. If the original digraph is acyclic and v-structure-free with source vertex v, the graph
resulting from removing edges (u, w) with dG (v, u) 6= dG (v, w) is also acyclic and v-structurefree.
Proof. Acyclicity is obvious as we are removing edges. To show the v-structure-freeness,
consider a v-structure x → y ← z. We note that if a v-structure is created, it means that the
shielding vertex {x, z} has been removed, that is, dG (v, x) 6= dG (v, z). But since (x, y) and (z, y)
are unresolved, dG (v, x) = dG (v, y), and dG (v, y) = dG (v, z), which is a contradiction.
Due to Lemma 1, we can orient the edges of G and obtain a member of the corresponding
MEC as follows. First one vertex is chosen to be the source vertex. This will determine the
orientation of some edges according to Proposition 2. Removing the resolved edges and resulted
isolated vertices, we will have (possibly) disconnected components, where the components
should again not contain cycles or v-structures. In each component, one vertex is chosen as the
source vertex of that component. This procedure is repeated until no unresolved edge is left. The
idea in the aforementioned procedure can be used in a recursive manner to count the members of
the MEC. The main tool in our approach is presented in Theorem 1.
Theorem 1. Let W(v, G) be the number of elements of the MEC with connected and undirected
essential graph G, in which v is the source vertex; and let H be the undirected graph remained
4

after removing edges with one endpoint closed to v than the other, as well as removing the
resulting isolated vertices. We have
W(v, G) =

Y

X

W(u, G0 ),

G0 ∈comp(H) u∈V (G0 )

where comp(H) denotes the set of the components of H, and by convention, empty product
(nullary product) is assumed to be the multiplicative identity 1.
Proof. For each component of H, choosing different variables of that component as the source
partitions the set of elements of the MEC. That is, two different choices for the source do not
lead to the same DAG (due to Proposition 2), and no elements will be remained uncounted, as
each variable is chosen as the source once. Therefore, we have the inner summation.
To justify the outer product, we need to show that orientation of one component of H does
not cause any restrictions on the possible orientations in any other component. This is due to
the fact that each component is at a certain depth level from the source and they cannot create
v-structures pointing to lower levels. Also, for the same reason and the fact that the edges from
the source are towards each component, the edges in two component cannot create a cycle.
Therefore, using Theorem 1, in order to count the members of an MEC, it suffices to calculate
W(v, G) for all v ∈ V (G) and sum them up, that is,
X
Size(G) =
W(v, G).
(1)
v∈V (G)

The pseudo-code of our approach is presented in Algorithm 1. This code consists of two
functions F LOWED and W. The function F LOWED takes an undirected and connected graph with
one of its vertices as the input, and returns the set of edges whose orientation is resolved when
the input vertex is set as the source vertex. Function W applies Theorem 1 in a recursive manner
to calculate the number of elements of the MEC with connected and undirected essential graph
G, when input vertex v is set as the source vertex.
Theorem 2. Let n and ∆ denote the order and the maximum degree of the graph G, respectively.
The computational complexity of Algorithm 1 is O(n∆ ).
Proof. We consider the operation of Algorithm 1 in rounds. In the first round, F LOWED is
called for all pairs (G, v1 ), v1 ∈ V (G). For any specific v1 ∈ V (G), we obtain undirected
components Gv11 , Gv21 , ..., GvK11 . In the second round, for each component Gvk1 , 1 ≤ k ≤ K1 , we
5

Algorithm 1 MEC Size Calculator
Input: P
Undirected graph G = (V (G), E(G)).
Size = v∈V (G) W(v, G).
Output: Size.
function F LOWED(v,G)
Initiate: A = ∅.
Set v as the source variable in G.
for {u, w} ∈ E(G) do
if dG (v, u) < dG (v, w) then A = A ∪ (u, w).
end if
if dG (v, u) > dG (v, w) then A = A ∪ (w, u).
end if
end for
return A.
end function
function W(v,G)
F = Undirected version of elements of F LOWED(v, G).
G0 = G\F .
Remove isolated vertices from G0 .
Y
X
return
W(u, G)
G∈comp(G0 ) u∈V (G)

end function
call F LOWED for all pairs (Gvk1 , v2 ), v2 ∈ V (Gvk1 ). For any specific v2 ∈ V (Gvk1 ), we obtain
undirected components Gv11 v2 , Gv21 v2 , ..., GvK12v2 , and we continue this procedure. It suffices to
show that this procedure could be repeated only for ∆ rounds. Therefore, since in each round
the function F LOWED is called at most n times, we get the desired result.
Suppose that in a round, we obtain an undirected component G0 and execute F LOWED on it
for a pair (G0 , v), where v ∈ V (G0 ). Let G01 , · · · , G0K be the resulting undirected components.
We first show that for any vertex u ∈ V (G0i ), (1 ≤ i ≤ K), we have obtained a directed edge in
G0 from some vertex w to u after applying F LOWED. By contradiction, suppose that there is no
such directed edge. Consider the shortest path from v to u. This path must pass through one of
neighbors of u in G0i such as w. Since d(v, w) < d(v, u), (u, w) should be oriented as w → u by
calling F LOWED procedure. But it is in contradiction with the fact that u and w are in G0i and
thus there should be such a directed edge in G0 .
For any vertex u ∈ V (G), consider the set of all edges in E(G) which intersect with u.
6

𝑣1

𝐺

𝑣2

𝑣1

(𝑏)

𝑣4

𝑣4

𝑣3

(𝑐)

𝑣2

(𝑎)

𝐺′

𝑣4

𝑣3

𝑣2

𝑣3

𝑣2

𝑣1

𝑣3

(𝑑)

𝑣4

𝐺′′

𝑣1
𝑣3

(𝑒)

Figure 1: Graphs related to Example 1.

Based on the reasoning above, at least one of these edges will be oriented in a round if vertex u
is still in one of undirected components before that round. Since the degree of u is at most ∆, all
of these edges will be oriented after at most ∆ rounds.

Example 1. Consider the undirected graph presented in Figure 1(a) as the essential graph given
as the input to Algorithm 1.
Setting vertex v1 as the source vertex, the direction of edges (v1 , v2 ), (v1 , v3 ), (v2 , v4 ), and
(v3 , v4 ) will be obtained using function F LOWED(v1 , G), as shown in Figure 1(b). Therefore,
the remaining undirected graph is G0 , shown in Figure 1(c). Setting each vertex in G0 as the
source vertex, no undirected component will be remained. Hence, W(v2 , G0 ) = W(v3 , G0 ) = 1.
Therefore, W(v1 , G) = W(v2 , G0 ) + W(v3 , G0 ) = 2. Due to the symmetry in the structure, we
have W(v4 , G) = 2 as well.
Setting vertex v2 as the source vertex, the direction of edges (v2 , v1 ), (v2 , v3 ), and (v2 , v4 ) will
be obtained using function F LOWED(v2 , G), as shown in Figure 1(d). Therefore, the remaining
undirected graph is G00 , shown in Figure 1(e). Setting each vertex in G00 as the source vertex,
no undirected component will be remained. Hence, W(v1 , G00 ) = W(v3 , G00 ) = W(v4 , G00 ) = 1.
Therefore, W(v2 , G) = W(v1 , G00 ) + W(v3 , G00 ) + W(v4 , G00 ) = 3. Due to the symmetry in the
structure, we have W(v3 , G) = 3 as well.
Finally, using equation (1), we obtain that the number of the elements of the MEC corresponding to essential graph G is equal to 10.

7

Algorithm 2 Uniform Sampler
Input: Undirected graph G = (V (G), E(G)).
Output: Directed graph (V (G), R AND E DGE(G, ∅)).
function R AND E DGE(G,A)
Set v ∗ ∈ V (G) as the source variable of G with probability P

W(v ∗ , G)
.
v∈V (G) W(v, G)

A = A ∪ F LOWED(v ∗ , G).
F = Undirected version of elements of F LOWED(v ∗ , G).
G0 = G\F .
Remove isolated
vertices from G0 .
[
A=A
R AND E DGE(G, A).
G∈comp(G0 )

return A.
end function

4

Uniform Sampling from an MEC

As mentioned earlier in Section 1, in some applications, the performance of an action is needed
to be tested on an MEC with large size. In these cases, one can optimize the behavior on uniform
samples drawn from the MEC. In this section we introduce two samplers for providing random
DAGs from an MEC. The first sampler, presented in Subsection 4.1, is based on the counting idea
in Section 3. In the case of bounded degree graphs, this sampler is capable of producing uniform
samples in polynomial time. If the degree is not bounded, in the worst case, this sampler has high
computational complexity. Therefore, in Subsection 4.2, we provide another sampler, originally
introduced in (Ghassami et al., 2017), for the case of unbounded degree graphs. Although
this sampler is not unbiased, our extensive experimental results confirms that the sampling
distribution of the second sampler is very close to uniform.

4.1

Uniform Sampler

The pseudo-code of our uniform sampler is presented in Algorithm 2, where functions F LOWED
and W are defined in the pseudo-code for Algorithm 1. In this sampler, first the weights W(v, G)
is calculated for all vertices v ∈ V (G) and vertex v ∗ is chosen as the source variable with
P
probability W(v ∗ , G)/ v∈V (G) W(v, G). The source variable will determine the orientation
of some edges according to Proposition 2. Removing the directed edges and resulted isolated
vertices, we will have (possibly) disconnected components. Choosing a source variable in each
8

resulting component is repeated until all the edges are directed. Since we choose each vertex as
the source variable according to the portion of members of the MEC having that vertex as the
source, we will have the following result.
Theorem 3. The sampler proposed in Algorithm 2 is uniform.
Proof. Let g be a member of the MEC. For the clarity of the proof, we assume that g is obtained
form three rounds of choosing source vertices, otherwise, the idea is easily extendable to the
case with arbitrary many rounds. Therefore, assume g is obtained by having v1 as the chosen
source vertex in G in round 1, v21 , ..., v2K1 as the sources in components G11 , ..., G1K1 in round 2,
and v31 , ..., v3K2 as the sources in components G21 , ..., G2K2 in round 3. We have
W (v21 , G11 ) · · · W (v2K1 , G1K1 )
W (v1 , G)
× QK1 P
1
v∈V (G) W (v, G)
k=1
v∈V (G1 ) W (v, Gk )

P (g) = P

k

W (v31 , G21 ) · · · W (v3K2 , G2K2 )
.
× QK1 Q
P
2
k=1
G2 ∈comp(G1 )
v∈V (G2 ) W (v, Gj )
j

j

k

Since the process is terminated in three rounds, we have W (v3k , G2k ) = 1, for all 1 ≤ k ≤ K2 .
Therefore,
1
W (v1 , G)
× QK1 P
1
v∈V (G) W (v, G)
k=1
v∈V (G1 ) W (v, Gk )

P (g) = P

k

×

K1
Y

W (v2k , G1k )
Q

k=1

G2j ∈comp(G1k )

P

v∈V (G2j )

W (v, G2j ))

.

By Theorem 1, we have the second fraction and all the fractions in the product are equal to 1.
Therefore,

1
1
=
,
Size(G)
v∈V (G) W (v, G)

P (g) = P

where the last equality follows from equation (1).
Regarding the complexity of the sampler in Algorithm 2 we have the following result as a
corollary of Theorem 2.
Corollary 1. Let n and ∆ denote the order and the maximum degree of the graph G, respectively.
The computational complexity of Algorithm 2 is O(n∆ ).

9

Algorithm 3 Heuristic Sampler
Input: Undirected graph G = (V (G), E(G)).
Uniformly shuffle the order of the elements of V (G).
while the induced subgraph on any subset of size 3 of the variables is not directed, or a directed
cycle, or a v-structure do
for all {vi , vj , vk } ⊆ V (G) do
Orient the undirected edges among {vi , vj , vk } independently according to Bern( 12 )
until it becomes a directed structure which is not a directed cycle or a v-structure.
end for
end while
Output: Resulted directed graph.

4.2

Heuristic Sampler

For the case of unbounded degree, the computational complexity of Algorithm 2 will be high
in the worst case. Therefore, we also provide a heuristic sampler for this case, which was first
introduced in (Ghassami et al., 2017). Although not unbiased, the experimental results indicate
that this sampler has a performance very close to the uniform sampler. The pseudo-code of
the heuristic sampler is presented in Algorithm 3. In this algorithm we consider all subsets
of size 3 from V (G) in a uniformly random order (achieved by uniformly shuffling the labels
of elements of V (G)). For each subset {vi , vj , vk }, we orient the undirected edges among
{vi , vj , vk } independently according to Bernoulli( 12 ) distribution. If the resulting orientation
on the induced subgraph on {vi , vj , vk } became a directed cycle or a v-structure, we redo the
orienting. We keep checking all the subsets of size 3 until the induced subgraph on all of them
are directed and none of them is a v-structure or a directed cycle.
Theorem 4. Each generated directed graph in Algorithm 3 belongs to the Markov equivalence
class corresponding to the essential graph G.
We require the following lemma for the proof.
Lemma 2. A chordal graph has a directed cycle only if it has a directed cycle of size 3.
Proof. If the directed cycle is of size 3 itself, the lemma is trivial. Suppose the cycle Cn is
of size n > 3. Relabel the vertices of Cn to have Cn = (v1 , ..., vn , v1 ). Since the graph is
chordal, Cn has a chord and hence we have a triangle on vertices {vi , vi+1 , vi+2 } for some i.
If the direction of {vi , vi+2 } is (vi+2 , vi ), we have the directed cycle (vi , vi+1 , vi+2 , vi ) which
is of size 3. Otherwise, we have the directed cycle Cn−1 = (v1 , ..., vi , vi+2 , .., vn , v1 ) on n − 1
10

vertices. Relabeling the vertices from 1 to n − 1 and repeating the above reasoning concludes
the lemma.
Proof of Theorem 4. The undirected graph G is chordal (Hauser and Bühlmann, 2012). Therefore, by Lemma 2, to insure that a generated directed graph is a DAG, it suffices to make sure
that it does not have any directed cycles of length 3, which is one of the checks that we do in
the proposed procedure. For checking if the generated DAG is in the Markov equivalence class
corresponding to G, by Proposition 1, it suffices to check if they have the same set of v-structures,
which is the other check that we do in the proposed procedure.

11

References and Notes
Andersson et al., 1997. Andersson, S. A., Madigan, D., Perlman, M. D., et al. (1997). A
characterization of markov equivalence classes for acyclic digraphs. The Annals of Statistics,
25(2):505–541.
Bernstein and Tetali, 2017. Bernstein, M. and Tetali, P. (2017). On sampling graphical markov
models. arXiv preprint arXiv:1705.09717.
Chickering, 2002. Chickering, D. M. (2002). Optimal structure identification with greedy search.
Journal of machine learning research, 3(Nov):507–554.
Ghassami et al., 2017. Ghassami, A., Salehkaleybar, S., Kiyavash, N., and Bareinboim,
E. (2017).

Budgeted experiment design for causal structure learning.

arXiv preprint

arXiv:1709.03625.
Gillispie and Perlman, 2002. Gillispie, S. B. and Perlman, M. D. (2002). The size distribution
for markov equivalence classes of acyclic digraph models. Artificial Intelligence, 141(1-2):137–
155.
Hauser and Bühlmann, 2012. Hauser, A. and Bühlmann, P. (2012). Characterization and greedy
learning of interventional markov equivalence classes of directed acyclic graphs. Journal of
Machine Learning Research, 13(Aug):2409–2464.
He et al., 2013. He, Y., Jia, J., Yu, B., et al. (2013). Reversible mcmc on markov equivalence
classes of sparse directed acyclic graphs. The Annals of Statistics, 41(4):1742–1779.
Meek, 1995. Meek, C. (1995). Causal inference and causal explanation with background
knowledge. In Proceedings of the Eleventh conference on Uncertainty in artificial intelligence,
pages 403–410. Morgan Kaufmann Publishers Inc.
Spirtes et al., 2000. Spirtes, P., Glymour, C. N., and Scheines, R. (2000). Causation, prediction,
and search. MIT press.
Verma and Pearl, 1991. Verma, T. and Pearl, J. (1991). Equivalence and synthesis of causal
models. In Proceedings of Sixth Conference on Uncertainty in Artificial Intelligence, pages
220–227.

12

