Efficient algorithm to enumerate all member DAGs of a Markov equivalence class

Question

I'm working on a research project involving Bayesian networks. BNs are directed acyclic graphs (DAGs) used to compactly represent joint distributions of variables. In many cases, multiple DAGs can represent the same distribution, and so compose an equivalence class. Equivalent DAGs share the same skeleton (same edges if one ignores direction), and the orientation of edges in their v-structures (X->Y, Z->Y, but no edge between X and Z) must be the same. Such equivalence classes can be represented by partially directed acyclic graphs (PDAGs), with the compelled (i.e. have the same orientation) edges directed and the reversible edges undirected.

Given a DAG, I would like to be able to find all other DAGs in its equivalence class. I'm aware of an efficient algorithm (Chickering 1995) that, given a DAG, finds the complete PDAG representing its equivalence class. I'm also aware of an algorithm (Dor and Tarsi 1992) that, given a PDAG, generates a random member DAG of the equivalence class. I am not interested in a random member, but rather in enumerating all members.

It might seem trivial at first glance--why not just try all combinations of orientations of the reversible edges, and discard the ones that aren't acyclic? But the number of possibilities grows exponentially in the number of reversible edges, so I'm worried that this won't work for large graphs. I've seen some sources claim that the proportion of edges that are reversible is relatively small, but still, if the graph is large enough, it will still be problematic. I'm working with graphs of up to several hundred variables, with several thousand being a possibility later on.

Answer 1

The answer provided by George works. I just want to add that there are more efficient ways to do it. We recently wrote a paper on this topic: https://arxiv.org/abs/2301.12212. The code is available here: https://github.com/mwien/fastmecenumeration.

This algorithm lists the DAGs one-by-one and has worst-case $O(n+m)$ delay between successive DAGs (for (P)DAGs with $n$ vertices and $m$ edges). So, the time that passes between two outputs grows linearly in the size of the graph.

The algorithm is too involved to fully explain here. The main point is that it utilizes connections between chordal graphs and acyclic orientations without any (new) v-structures (which is precisely what you want for listing the DAGs represented by a (C)PDAG). A chordal graph is an undirected graph, in which every cycle of length four or larger has at least one chord (e.g. there exist two non-neighboring vertices of the cycle, which have an edge between them). Interestingly, an undirected graph has an acyclic orientation without v-structures if, and only if, it is chordal (imagine a cycle of length four without a chord; it is not possible to orient the edges acyclically without introducing a v-structure; whereas with a chord it is possible).

By using these connections, it is possible to exploit algorithms from chordality testing (such as Lexicographic BFS or Maximum Cardinality Search) and extend them (in a non-trivial way) to efficiently list all DAGs in a Markov equivalence class. Basically, these algorithms give you a single DAG (similar to Dor-Tarsi's algorithm), but you can modify them to recursively try out different orientation choices in order to list all DAGs.

Moreover, as you are worried that the number of DAGs might be too large, it is possible to efficiently compute the number of DAGs represented by a CPDAG (without listing and counting all of them one-by-one): https://arxiv.org/abs/2205.02654. So, one could also use that and only start the listing algorithm if it is feasible.