Introduction

SAT solvers have become popular means of solving computationally hard problems. However, most modern SAT solvers require the problem to be specified in conjunctive normal form (CNF). This may pose difficulty in problems comprising of constraints that do not have a straightforward translation to propositional logic. This paper targets one such constraint, graph acyclicity.

A graph is naturally modeled by a set of Boolean variables x jk expressing the existence of an edge from j to k. So, for instance, a constraint x 12 ∨ x 13 expresses that we're looking for a graph where the node 1 is connected to at least one of the nodes 2 and 3. But how do we ensure that we are looking for a graph that is acyclic?

Acyclicity is a core concept from graph theory and naturally arises in a number of applications. In planning it is used to ensure causality [42]. Similarly, it is needed in networks used for computation, e.g. Bayesian networks [49]. More recently, software verification of concurrent programs relies on relaxed memory models, which are defined in terms of acyclic relations [34,51].

There is a large body of work of translating constraints into CNF. Prime examples are cardinality constraints [5,36,43] pseudo-Boolean constraints [6,17], XOR-constraints [24] and general CSP [46]. A number of graph-related problems approached by SAT can be found in the literature. For instance, encoding graph connectivness property [47], graph coloring [50], calculating Steiner tree [14,30] or maximum clique [28]. It is also worth mentioning that graphs and SAT play an important role in the research on Gene regulatory networks [15].

The paper has the following contributions. It reviews CNF encodings that are found in the literature, and introduces a number of new ones. All these encodings were implemented and extensively evaluated. To generate a well-defined testbed we introduce The Supervisor Problem, which asks to assign supervisors to a group of persons, where there is a lower bound for each person's supervisors and an upper bound on how many persons he or she can supervise. This problem is a generalization of the well-known ordering principle (also known as the GT family) studied in proof complexity [1,10,32,44].

Preliminaries

Our basic setting is the following. We represent a directed graph on n ≥ 1 vertices by Boolean variables x ij , associated with edges. A family of graphs is defined by a formula φ n (� x, � y ) over the edge variables � x and some auxiliary variables � y. Each graph corresponds to a satisfying assignment of ∃� y φ n (� x, � y). Given such a formula φ n , we wish to construct another formula ψ n , such that the satisfying assignments of ∃� y ∃� z � φ n (� x, � y ) ∧ ψ n (� x, � z ) � correspond to the graphs from the given family that, in addition, are acyclic. We refer to the formula ψ n as the acylicity checker or, for short, the checker. In the rest of the section, we make this more precise, and we describe the particular formulas φ n that we use in experiments. The acyclicity checkers ψ n are discussed later ( §3). Standard terminology from propositional logic and graph theory is assumed.

Propositional Logic and SAT Solving

Propositional formulas are constructed from propositional variables by Boolean connectives (∨, ∧, ⇒) and negation (¬). A propositional formula is in conjunctive normal form (CNF) if it is a conjunction of clauses, where a clause is a disjunction of literals, where a literal is a propositional variable or its negation.

Whenever convenient, we treat CNF as a set of clauses and a clause as a set of literals. Observe that the empty set of clauses corresponds semantically to true and the empty clause corresponds semantically to false. The empty clause is denoted as ⊥.

Any propositional formula can be converted to an equisatisfiable CNF formula in linear time [41,48].

The problem of satisfiability (SAT) is to decide for a given propositional formula whether there is a satisfying assignments for it or not.

The paper requires minimal understanding of SAT solving. We will assume that a SAT solver accepts a formula in CNF and gives the response "YES" if the formula is satisfiable and "NO" otherwise. Most modern SAT solvers also provide an actual satisfying assignment if the answer is "YES".

When discussing the performance of SAT solvers, sometimes we refer to the underlying proof system, which is propositional resolution. Given two clauses x∨C and ¬x∨D their resolution is C ∨D. The two clauses are called antecedents and the resulting clause the resolvent. Resolution proof of a clause C from a formula φ is a sequence of clauses C 1 , . . . , C k such that C k = C and each clause C i is either in φ or it is a resolvent of some C j , C k with 1 ≤ j < k < i. A resolution refutation is a resolution proof of the empty clause. A CNF formula is unsatisfiable if and only if there exists a resolution refutation for it. For most modern SAT solvers it holds that for any run of the SAT solver that responds "NO", there exists a resolution refutation of the input formula linear to the run of the solver. In fact, it even holds that for any resolution refutation there exists a polynomial run of the solver [39].

Graphs

Throughout the paper we will be discussing directed graphs. A directed graph G = (V, E) consists of a set V of vertices and a set E of edges. In this paper, we work with finite graphs, and we fix V = [n] = {1, . . . , n}. We denote vertices by i, j, k, . . . and we abbreviate the (directed) edge (i, j) by ij. For an edge ij, we say that it has source i and target j; conversely, we say that ij is an outgoing edge of i, and an incoming edge of j. The in-degree of a vertex is its number of incoming edges, and the out-degree of a vertex is its number of outgoing edges. A sink is a vertex with out-degree 0.

A graph property is a family G of graphs. In this paper, we consider only properties G n that contain graphs with a fixed number n of vertices, and are defined by a formula φ n . More precisely, φ n has one variable x ij for each pair (i, j) ∈ V × V , and may also use some auxiliary variables � y. Then, a graph with edges E is in G n if and only if it corresponds to a satisfying assignment of ∃� y φ n (� x, � y). We say that a property is monotone if adding edges to a graph preserves the property. This corresponds to the Boolean function ∃� y φ n (� x, � y ) being monotone. We remark that the property of having a cycle is monotone, and the intersection of two monotone properties is monotone.

The No-Sink Family

The graphs with no sinks are described by the following formula:

φ n (� x ) := � i∈[n] � � j∈[n] x ij � (no-sink)

Observe that the no-sink property is monotone. Also, the ordering principle says that graphs with no sink have a cycle. It follows that, no matter what acyclicity check ψ n we use, the formula ∃� y ∃� z � φ n (� x, � y ) ∧ ψ n (� x, � z ) � will be unsatisfiable. Since we know the answer, the no-sink family is a good test case for the correctness of our implementation. In addition, as we shall see, these unsatisfiable formulas are not easy for state-of-the art SAT solvers.

The Supervisor Problem

We define a parameterized family of graphs, which will typically contain both cyclic and acyclic graphs. The idea is as follows. Consider a group of n people where each person must have a certain number of supervisors and at the same time there is a limit on how many people a person can supervise. The supervisor problem is formalized as follows.

Input: A positive integer n and a sequence of pairs of nonnegative integers (u 1 , l 1 ), . . . , (u n , l n ). Output: "YES" if and only if there exists an acyclic graph with the vertices V = [n] such that vertex i has at most u i incoming edges and at least l i outgoing edges.

An instance n, � u, � l of the Supervisor problem reduces to the satisfiability of the following formula:

φ n,� u, � l (� x ) := � i∈[n] atleast(l i ; x i1 , . . . , x in ) ∧ � j∈[n] atmost(u j ; x 1j , . . . , x nj ) (supervisor)

The cardinality constraint atleast(l; � x ) evaluates to 1 when at least l of the Boolean variables � x evaluate to 1. There are several ways to encode such cardinality constraints in CNF [5,17,36,43]. The cardinality constraint atmost(u; � x ) is analogous.

The Supervisor problem has three interesting special cases: (a) the no-sink family; (b) the pigeonhole principle; and (c) the DAG realizability problem. We obtain a description of the no-sink family by requiring the bounds (n, 1) for each of the n vertices. We obtain a description of the pigeonhole principle by requiring the bounds (0, 1) for k 1 vertices and the bounds (1, 0) for k 2 vertices, with

k 1 + k 2 = n and k 1 > k 2 .

Finally, any instance of the DAG realizability problem is also an instance of the Supervisor problem, which satisfies the additional constraint � n i=1 u i = � n i=1 l i . Since the pigeonhole principle is a special case, it follows that (supervisor) sometimes requires large resolution refutations [25]. Since the DAG realizability problem is a special case, it follows that it is NP-complete to find a model of (supervisor) that corresponds to an acyclic graph [26]. Finally, we observe that (supervisor) is not monotone.

Encodings

We discuss several ways of constructing the acyclicity checker ψ n , based on reachability ( §3.1), based on labeling vertices with numbers ( §3.2 and §3.3), and based on oblivious algorithms ( §3.4 and §3.5).

Transitive Closure

A graph G = (V, E) does not have cycles if and only if E + is irreflexive. Here, E + denotes the irreflexive transitive closure of E. Further, G does not have cycles if and only if there exists an irreflexive transitive relation

R that contains E; that is, if there exists R ⊆ V × V such that (a) ii / ∈ R for all i ∈ V , (b) ij ∈ R and jk ∈ R only if ik ∈ R for all i, j, k ∈ V , and (c) ij ∈ R if ij ∈ E for all i, j ∈ V .

We represent the relation R by variables y ij , and we define ψ n by clauses that directly correspond to the conditions (a), (b), and (c):

ψ n (� x, � y ) := � i∈[n] ¬y ii ∧ � i,j,k∈[n] (y ij ∧ y jk → y ik ) ∧ � i,j∈[n] (x ij → y ij )(tc1)

A possible alternative is the following:

ψ n (� x, � y ) := � i∈[n] ¬y ii ∧ � i,j,k∈[n] (y ij ∧ x jk → y ik ) ∧ � i,j∈[n] (x ij → y ij )(tc2)

In relational language, (tc1) corresponds to the conditions (a)

R ∩ id = ∅, (b) RR ⊆ R, and (c) E ⊆ R.

Here, id is defined to be { ii | i ∈ V }. In (tc2), we replace condition (b) by RE ⊆ R. Juxtaposing relations corresponds to normal relational composition:

RE := { (i, j) | (i, k) ∈ R and (k, j) ∈ E, for some k }

If the family of graphs satisfying φ n is monotone, in the sense that adding edges preserves φ n , then we do not need to use auxiliary variables y ij , as we can simply reuse the variables x ij :

ψ n (� x ) := � i∈[n] ¬x ii ∧ � i,j,k∈[n] (x ij ∧ x jk → x ik )(tc3)

The conjunction of (tc3) and (no-sink) has Ω(2 n ) regular tree resolution proofs but O(n 3 ) dag resolution proofs; see [10,32,44] Observe that the formulas in each of (tc1), (tc2), and (tc3) have size Θ(n 3 ).

Binary Labeling

A graph is acyclic when it can be embedded in a linear order; that is, when there exists a labeling l : V → Z of the vertices such that ij ∈ E implies l(i) < l(j) for all i, j ∈ V . Recall that V = [n]. Clearly, it is sufficient to use b := �log 2 n� bits for labels, so l : V → {0, 1} b . Accordingly, we use n groups � y 1 , . . . , � y n of auxiliary variables, each containing b Boolean variables. The acyclicity checker is then [42]:

ψ n (� x, � y 1 , . . . , � y n ) := � i,j∈[n] � x ij → less(� y i , � y j ) � (bin)

The implementation of less depends on how we represent numbers. Because we use binary numbers, we set less := lex b , and define the latter to test lexicographic ordering of the bits:

lex 0 () := 0 (1) lex b (� yy, � zz) := (¬y ∧ z) ∨ � (¬y ∨ z) ∧ lex b−1 (� y, � z ) �(2)

Converting this circuit to CNF will require some more auxiliary variables, which we do not denote explicitly. The conversion is standard, and the size of the resulting CNF is Θ(b).

Observe that the formula in (bin) has size Θ(n 2 log n).

Unary Labeling

We can also use the labeling idea with unary numbers. To represent a number, instead of �log 2 n� bits we use n − 1 bits, and force them to follow the pattern 0 * 1 * . We could still compare these numbers by their lexicographic order using lex, but we can also write directly a CNF formula, thus avoiding the (small) overhead associated with converting a circuit into CNF. The acyclicity checker is defined by

ψ n (� x, � y 1 , . . . , � y n ) := � i,j∈[n] � x ij → less(� y i , � y j ) � ∧ � i∈[n] unary(� y i ) (unr)

where each � y i uses n − 1 bits. We could use less := lex n−1 , but we prefer less(� y, � z ) := ∃� u lessunr(� y, � z, � u ), where

lessunr(� y, � z, � u ) := n−1 � i=1 � (¬y i ∨ ¬u i ) ∧ (z i ∨ ¬u i ) � ∧ n−1 � i=1 u i(3)

The auxiliary � u variables will identify a position i where y i = 0 and z i = 1. This is sufficient because we enforce the pattern 0 * 1 * .

unary(� y ) := n−1 � i=2 (y i−1 → y i )(4)

Observe that the formula in (unr) has size Θ(n 3 ).

Warshall

The Warshall algorithm computes transitive closure, in cubic time.

Warshall

1 for k ∈ [n] 2 for i ∈ [n] 3 for j ∈ [n] 4 a ij := Or(a ij , And(a ik , a kj ))

Initially, a ij says if the graph contains an edge ij; and, at the end a ij says if the graphs contains a path from i to j. In particular, at the end, a ii says if there is a cycle containing vertex i.

In a usual implementation, the procedures Or, And are simply performing the corresponding Boolean operations. Instead, we could replace each of them with a function that constructs the corresponding circuit. Because the Warshall algorithm is simple enough, we can also build the acyclicity checker directly. We use n 2 (n + 1) auxiliary variables y ijk , for i, j ∈ [n] and k ∈ {0, . . . , n}. The variable y ijk will be 1 only if there exists a path from i to j that uses intermediate vertices only from the set [k].

ψ n (� x, � y) := � i∈[n] ¬y iin ∧ � i,j∈[n] (x ij → y ij0 ) ∧ � i,j,k∈[n] (y ij(k−1) → y ijk ) ∧ � i,j,k∈[n] (y ik(k−1) ∧ y kj(k−1) → y ijk ) (fw)

We remark that the formula above is close but not the same as what we would obtain by replacing the basic operations Or, And with circuit builders. One difference is that we use implications rather than equivalences: It is possible that y ijk = 1 even if there is no corresponding path. Another difference is that the last conjunct built by Warshall would be (

y ik(k−[k≥j]) ∧ y kj(k−[k≥i]) ) ↔ y ijk ,

where [p] evaluates to 1 when p is true, and 0 when p is false. Either form is correct.

The formula in (fw) is similar to the one in (tc1). It has the same size Θ(n 3 ), but uses more auxiliary variables.

Matrix Multiplication

Let us take stock of the encodings we have so far. Four encodings have size Θ(n 3 ), namely (tc1), (tc2), (unr), and (fw); another encoding has size Θ(n 2 log n), namely (bin). But, size is not the only criterion on which to judge a satisfiability formula. We note that (tc1), (tc2), and (fw) have the following desirable property: if the edge variables � x are fixed, then satisfiability is decided by unit propagation. The smaller encoding (bin) does not have this property. Is it possible to have a formula of sub-cubic size that is solved by unit propagation once the edge variables � x are fixed? We answer this question affirmatively by showing an encoding based on matrix multiplication. We first describe the algorithm. The formula ψ n will be constructed by replacing the basic operations with circuit builders, and then converting the circuit to CNF.

We will consider n×n matrices of two kinds, over integers and over Booleans. Integers form a commutative ring, so the corresponding matrix operations are unrestricted. Booleans do not form a ring: We take addition to be logical-or and we take multiplication to be logical-and, but there is no subtraction. Since we use both types of operations and sometimes reinterpret the same matrix as being over integers or over Booleans, we distinguish Boolean matrices with a hat. When we reinterpret an integer matrix X as a Boolean matrix X, all non-zero integers collapse to 1.

Let X be the adjacency matrix of a graph, with entries in {0, 1}. Entry (i, j) of X l is the number of paths of length l from i to j, and entry (i, j) of Xl tells us if there exists a path of length l from i to j. Thus, to decide if the graph has a cycle it suffices to check the diagonal of the Boolean matrix Xl , and also that Bk equals � n l=1 Xl whenever k ≥ log 2 n. Thus we can decide if a directed graph is acyclic in the time needed to perform Θ(log n) multiplications of n × n Boolean matrices. One could perform a naive matrix multiplication in Θ(n 3 ) time, giving an acyclicity checker ψ n of size Θ(n 3 log n).

Alternatively, to multiply Boolean matrices  and B, we multiply them as integer matrices obtaining C = AB; the product of  and B is then Ĉ. (The idea of using intermediate integer matrices is old [18,19,37].) To perform the multiplication of integer matrices, we use Strassen's algorithm [45], which works in O(n 2.81 ) time. To achieve this improved runtime, Strassen's algorithm makes essential use of subtraction, which is not available over Booleans. However, the runtime O(n 2.81 ) counts operations over integers, while we must count operations over bits since we ultimately want to build a circuit. Since the matrices A and B, which we multiply, have entries from {0, 1}, the product matrix C has entries from {0, . . . , n}. Some intermediate values in Strassen's algorithm are outside this range. Nevertheless, one can check that Θ(log n) bits are sufficient to represent all intermediate values. So, the size of the circuit is O(n 2.81 log n • f (log n)), where f (w) is the size of a circuit necessary to perform an arithmetic operation on integers with w bits. In our implementation, we use a simple quadratic multiplication algorithm, leading to a circuit of size O(n 2.81 log 3 n).

In theory, one can use the same approach to achieve O � n ω log 2 n log log n • 8 log * log n � by using better algorithms for the multiplication of matrices [20] and for the multiplication of integers [27]. But such algorithms are known to not be practical. Our current implementation based on Strassen is not competitive either ( §4). Improving it seems to require new ideas.

Discussion. The approach we took for both Warshall and for matrix multiplication was to replace the basic operations in the algorithm with circuit builders. This approach is fairly general, but it does not apply to all algorithms, just to the oblivious ones. An algorithm is said to be oblivious when the memory locations it accesses do not depend on the input data. If an algorithm runs in t(n) time in the Turing model of computation, then one can construct an oblivious algorithm that runs in O � t(n) log t(n) � time, and hence a Boolean circuit of that size [40]. In contrast, a similar construction does not exist in the deterministic RAM model of computation [23].

The standard (and optimal) algorithm for finding cycles in directed graphs is DFS (depth-first search). However, we are not aware of an implementation of DFS that works in sub-cubic time in the (oblivious, deterministic, multi-tape) Turing machine model of computation.

Experimental Evaluation

Overall setup. Three popular state-of-the-art SAT solvers were used for the evaluation: Lingeling Version 276 [8], glucose 4.0 [3], and minisat 2.2 (from github) [16]. Further, minisat was a run in two modes: the default mode, and with additional preprocessing techniques enabled (asymmetric breaking and redundancy checking). The encodings we consider are Transitive closure (tc1), Unary labeling (unr), Binary labeling (bin), Warshall (fw), and an alternate version of Transitive closure (tc2). The techniques based on matrix multiplication are not included, because they quickly lead to formulas in the realm of hundred of megabytes, and more (Figure 3).

The experimental results were obtained on an Intel Xeon 5160 3GHz, with 4GB of memory, time limit 500 seconds, and memory limit 2GB.

Generation of Graph Families. The generation of the no-sink graph family is trivial: we simply instantiate (no-sink) for different values of n. We also generate random instances of the Supervisor problem. Each upper bound u i on the indegree of vertex i is chosen uniformly at random from the set {0, . . . , n − 1}. Each lower bound l i is 0 with probability p, which is a parameter, and otherwise is chosen uniformly at random from the set {1, . . . , n − 1}. If the generated sequence (l 1 , u 1 ), . . . , (u n , l n ) is not satisfied by any directed graph, we throw it away and try again. This filtering is easy, because it is possible to check in linear time whether a sequence is satisfiable by some directed graph, using the Fulkerson-Chen-Anstee theorem [2]; but, as we mentioned, checking whether a sequence is satisfiable by an acyclic directed graph is NP-complete [26].

For each n ∈ {2, 3, . . . , 50} and for each p ∈ {10%, 20%, . . . , 90%}, we generate a sequence (l 1 , u 1 ), . . . , (l n , u n ) using the method just described. For each such sequence, we generate the formula (supervisor). The atmost and atleast constraints were encoded into CNF using the well-known Sequential Counter encoding [43].

Presentation of the results. Figure 1 presents the number of instances solved. Figure 1(a) shows the results for graphs with no sinks, which lead to the GT family of formulas when coupled with the acyclicity checker (tc1). Figure 1(b) shows the results for the general supervisor problem. Figure 2 shows a cactus for the (no-sink) family for all the considered combinations of solvers and encodings. Recall that if a cactus plot contains a point (i, s) it means there are i instances such that each was solved within s seconds. (Therefore, the graph is necessarily increasing.) A more detailed presentation of the results can be found on the authors' homepage. 3The first thing to observe is that the supervisor family turns out to be easy for all the approaches. Indeed, all approaches solved over 90% benchmarks. Even though the differences are rather small, the tables suggests that the Transitive closure encoding ((tc1)) is good for all solvers. In contrast, Binary labeling (bin) appears to be the worst.

In the case of (no-sink) graphs, the results are widely different for each encoding but also for each solver. Most notably, lingeling performs extremely well on (tc1). Both versions of minisat perform well on Unary labeling (unr). Glucose also performs the best with (unr), but once the graph size gets over 30 nodes, the performance quickly deteriorates. The worst encoding appears to be Binary labeling (bin) and the alternative Transitive closure (tc2), with which all solvers could only handle less than half of all instances. Warshall has a peculiar behavior, as it performs poorly for all the solvers except lingeling.

Discussion. The results are surprisingly uneven across different solvers. In a particular, lingeling performs extremely well with (tc1), while the other solvers do not perform very well with the same encoding. It is possible that lingeling benefits from some simplification technique (pre-or in-processing) that is missing from the other solvers, but this remains to be understood. Binary labeling (bin) is overall the most compact but with the poorest performance. This is not entirely surprising. First, binary labeling does not offer many opportunities for applying unit propagation. And second, binary labeling makes it virtually impossible to learn anything about any particular variable.

In contrast, Unary labeling (unr) fares quite well over the considered solvers. Further, we observed that minisat with preprocessing solves the (no-sink) family with no conflicts, i.e., merely with preprocessing. This also holds for lingeling but only up to a certain size of the graph. It is not difficult to see how these formulas can be completely solved by preprocessing, in particular by failed literal detection and asymmetric breaking.

In unary labeling ( §3.3), the ordering of labels is ensured by the following constraint. For every vertex i labeled by a number l(i) ∈ {0, . . . , n − 1} there must be a neighbor j labeled by a number l(j) ∈ {0, . . . , n − 1} such that l(i) < k ≤ l(j) for some k ∈ {0, . . . , n − 1}. To show that such a labeling does not exist for no-sink graphs, we reason as follows. If k = n − 1, it immediately follows that vertex j has no neighbor labeled with a greater number since it has reached the maximum. Performing this reasoning for every vertex yields that the labels must be in the interval {0, . . . , n − 1}. Repeating this reasoning n times leads to a contradiction.

We argue that to perform this reasoning, it is sufficient to use unit propagation on unary labeling. In definition (3), setting u k := 1 imposes exactly the condition as listed above, i.e., the label of the source node is at most k − 1 and the target node at least k. It is easy to verify that setting u n−1 := 1 yields a conflict by unit propagation.

This idea is realized in the preprocessing techniques failed literal detection [12,35] and asymmetric breaking [38]. More generally, if setting a literal l to true yields a conflict by unit propagation it is called a failed literal and the literal is set to false. In this manner, these two techniques acting in tandem can gradually remove the u i literal corresponding to the highest number, and eventually derive an empty clause.

Conclusion.

The experimental results support the following conclusions. The combination of choice for checking acyclicity is lingeling with Transitive closure (tc1). (Do not use the seemingly very similar alternative Transitive closure (tc2), which leads to poor performance.) Minisat is the second solver of choice but not with the (tc1) encoding-instead, use (unr). SAT solvers appear to be good at finding an acyclic graph if it exists; proving that an acyclic graph does not exist can easily become very difficult. Formula Size. Figure 3 presents the size of the acyclicity checkers ( §3). The cubic methods (tc1, tc2, unr, fw) are all within a factor of 3 of each-other. The super- cubic one (mm, which is based on the definition of Boolean matrix multiplication) is already worse than the cubic ones at n = 3, and it only gets worse for larger n.

There are two sub-cubic methods: ss with size O(n 2.81 log 3 n), and bin with size O(n 2 log n). The method based on Strassen (ss) generates large formulas even for n = 10 and, looking at Figure 3, there is no overtaking in sight. Indeed, even if the constants hidden in the O-notation were just as good for ss as they are for tc1, the crossing point would be at n ∼ 10 32 . The log 3 n factor is important, and it seems that without new ideas this method will never be practical. On the other hand, the method based on binary labeling (bin), generates the smallest acyclicity checkers for n ≥ 100.

Summary and Future Work

The paper presents a number of encodings for the graph acyclicity constraint into conjunctive normal form. This is useful when calculating acyclic graphs using SAT solvers but also other logic-based solvers, such as quantified Boolean formula (QBF) and Satisfiability modulo theories (SMT) solvers. The experimental evaluation suggests that the performance is highly sensitive both to be encoding but also to the SAT solver used. This suggests that for practical purposes it is useful to maintain a portfolio of encodings and solvers.

There are a number of directions for further research. All the encodings that are asymptotically better than cubic perform poorly in practice. Hence, it is a challenge for future research to identify encodings with asymptotic complexity better than cubic but that also perform well.

Ultimately, all the considered encodings are not feasible for large graphs. Indeed, any cubic encoding generates formulas that are too large once the number of nodes gets over 1000. A possible solution is to consider a lazy approach where only some portion of the encoding is considered. However, this will probably require an iterative sequence of calls to the SAT solver. Another approach, tried recently [9,13], is to integrate specialized algorithms for acyclicity inside the SAT solver. This is also closely related to loop formulas, which are used in Answer Set Programming (ASP) to avoid cyclic reasoning cf. [21,29,33]. Some questions about properties of the encodings with respect to unit propagation remain open. In particular, is it the case that any partial assignment that contains a cycle produces a conflict by unit propagation in a given encoding? In theory this is possible due to monotonicity [7] (adding edges preserves cycles). We believe this holds for the encodings Transitive closure, Warshall, and matrix multiplication, but does not hold for Unary and Binary labeling. Hence, this also begs the question whether the labeling methods can be made more deterministic. Other properties related to unit propagation are also worth exploring, e.g., generalized arc consistency and propagation completeness, cf. [4,7,11,22,24].

The experimental results lead to several interesting questions. Why is lingeling so efficient on the transitive closure encoding (tc1)? Why does the alternative encoding for transitive closure (tc2) perform so poorly across all solvers? Do formulas (tc2) lack short resolution refutations?

Last but not least, this research opens avenues for exploring other graphrelated properties, e.g., connectivness and reachability in general.