Reasoning about Connectivity Constraints

1 Introduction

Many problems involve seeking a subgraph of a larger graph and connectivity is one of the most natural property to ask for in such subgraphs. It is, in particular, an essential property in several problems related to ecosystem design. For instance, in the problem of selecting sites to be included in reserve networks, fragmented networks are undesirable [Briers, 2002], and many other problems, such as forest harvesting or town planning involve this criterion. The problem of designing connected wildlife corridors has received a lot of attention recently, with a multi-commodity flow encoding [Conrad et al., 2007; 2012], a subtour elimination scheme [Dilkina and Gomes, 2010], and a hybrid local search method [LeBras et al., 2013]. Whilst these methods are very efficient, a constraint programming approach has a lot to offer in addition. In particular, it can very easily be extended to incorporate additional constraints. For instance, establishing protected areas while taking into account socio-economic constraints [Polasky et al., 2005], designing incentives for private

Emmanuel Hebrard

CNRS, U. Toulouse Toulouse, France hebrard@laas.fr

Toby Walsh

NICTA & UNSW Sydney, Australia

toby.walsh@nicta.com.au

landowner to contribute to "ecosystem services" [Polasky et al., 2014] both require connectivity in conjunction with other constraints and criteria. However, standard constraint programming toolkits do not support connectivity constraints, and the alternative of using a decomposition is very unlikely to be efficient. The flow decomposition proposed in [Conrad et al., 2007] has a Ω(|E||V |) space complexity and hinders propagation. Moreover such decompositions are not handled efficiently by CP solvers [Bessiere et al., 2009]. A decomposition enforcing domain consistency is even larger (Ω(|V |4 )) and hence impractical. We therefore propose a global constraint to ensure that a subset of vertices induces a connected subgraph in O(|E|) time and space complexity. Furthermore, we also introduce efficient algorithms to compute minimal explanations for these global constraints, thus making it possible to use this constraint in clause learning solvers.

In Section 3, we consider the constraint simply ensuring that the subgraph is connected. We show that domain consistency can be enforced on this constraint with a worst case time complexity linear in the size of the graph and a minimal explanation can be computed in the same worst case time complexity. Next, in Section 4, we consider the case where a cost function is associated to vertices and there is an upper bound on the overall cost of the subgraph. This generalization is NP-hard. However, we propose incomplete filtering rules that are very efficient in practice. Then, in Section 5, we consider the case where we need to find several subgraphs. This case is also NP-hard in general, even for two shapes. However, we propose implied constraints for the case of two shapes in planar graphs. Finally, in Section 6, we test empirically the algorithms and implied constraints that we introduce and show that they are efficient in practice.

2 Background

We denote by G = (V, E) a graph with vertices V and edges E. G[S] is the subgraph induced by the set S ⊆ V . A constraint satisfaction problem (CSP) involves finding an assignment of values to its variables such that all its constraints are satisfied. Variables take values in a finite domain of values.

In the family of constraints that we consider, the subset of vertices of a subgraph is represented as an array of Boolean variables. For each vertex v ∈ V , we have a Boolean variable that takes the value 1 if v is in the subgraph and the value 0 if it is not in the subgraph. To simplify the notations, however, we view this array as a set variable S, that is, a variable whose values are sets that must be supersets of a lower bound ¯S and subsets of an upper bound ¯ S. We say that a set s is a support iff ¯S ⊆ s ⊆ ¯ S and G[s] is connected. We say that S is domain consistent (DC) iff ¯S is the intersection of all supports and ¯ S is the union of all supports s.

Several modern CP solvers employ clause learning, a technique that has had great success in SAT solvers [MarquesSilva et al., 2009]. In order to be used with clause learning, a propagator needs to label each pruned value by a clausal explanation [Katsirelos and Bacchus, 2005]. Clauses can include literals corresponding to assigning a value, V = a, or pruning a value V ̸= a ≡ ¬(V = a) for finite domain variables and literals corresponding to including a value a ∈ S, or excluding it a /∈ S ≡ ¬(a ∈ S) for set variables. A clausal explanation is correct with respect to domain consistency on a constraint c if making all its literals false makes c domain inconsistent. It is minimal if no strict subset is correct.

3 Single Shape Connectivity

In this section, we consider the constraint that ensures that a set of vertices is connected in the graph given as parameter.

Definition 1 (CONNECTIVITY) Let Gbe a graph and S a set variable: CONNECTIVITY[G](S) ⇐⇒ G[S] is connected

3.1 Propagation

Algorithm 1 enforces domain consistency on the constraint CONNECTIVITY. The algorithm works in two phases. In

Algorithm 1: DC on CONNECTIVITY (G, S)

the first, we simply remove the vertices that are not reachable from any vertex in ¯S (Line 2). We compute a maximal connected component C of G[ ¯ S] including any one vertex v in ¯S and remove from ¯ S any vertex that does not belong to this connected component. After this phase, notice that ¯ S is a support since G[ ¯ S] is connected. Therefore this is sufficient to detect disentailment. We can therefore set G to G[ ¯ S] (Line 3). In the second phase, we use the notions of biconnected component and articulation point. A biconnected component of G is a set of vertices p such that for any vertex v ∈ p, G[p \ {v}] is connected. A vertex is an articulation point if and only if removing it increases the number of connected components. Let p1 and p2 be two distinct maximal biconnected components, it follows that |p1 ∩ p2| ≤ 1, otherwise p1 ∪ p2 is biconnected and larger than p1 and p2.

Hopcroft and Tarjan's algorithm [Hopcroft and Tarjan, 1973] computes all articulation points in O(|V | + |E|), that is O(|E|) since G[ ¯ S] is connected. It explores the graph depth-first and maintains, for each vertex v, the lowest depth of neighbors of all descendants of v in the depth-first-search tree, called the lowpoint. A vertex v is an articulation point if and only if it has a child u in the dfs such that lowpoint[u] ≥ depth[v]. Indeed, it means that the only path from u to vertices with lower depth in the dfs goes through v hence removing v would cut the graph. By very slightly modifying Hopcroft and Tarjan's algorithm to record for each vertex v also the vertex whose depth is recorded as lowpoint[v], we get the collection PG of G's biconnected components.

A call to Hopcroft and Tarjan's algorithm (Line 4) returns the collection of sets PG, the set of articulation points AG of G and a table of predecessors pred corresponding to the depth-first-search tree. Next, in Loop 5, we build a tree T (G) = (V T , ET ) whose vertices stand for articulation points and biconnected components of G. (see Fig. 1):

• There is a vertex labeled p per biconnected component p ∈ PG and one labeled {a} per articulation point a in G, i.e., V T = (PG ∪ {{a} | a ∈ AG});
• There is an edge ({a}, p) if the articulation point a is in the biconnected component p, and an edge ({a1}, {a2}) if the articulation points a1 and a2 share an edge but do not belong to the same biconnected component.

Observe that in T(G) every edge involves at least an articulation point of G and it is indeed a tree, otherwise removing an articulation point in a cycle of T (G) would not disconnect G.

Lemma 1 Let x and y be two vertices of T(G) such that x∩ ¯S ̸= ∅ and y ∩ ¯S ̸= ∅, then every articulation point a of G such that {a} is on a path from x to y in T(G) must be in S.

Proof: Suppose that there is an elementary path from x to y in T(G) going through {a} where a is an articulation point of G. By construction of T(G), there is no path in G from any vertex in x to any vertex in y that does not go through a. Therefore, if vertices of both x and y are in the connected component, then a must also be in the component. □

A corollary of Lemma 1 is that once we have pruned T (G) until every leaf x is such that x ∩ ¯S ̸= ∅, we know that the same applies to every node of the tree that corresponds to an articulation point of G. It follows that for every articulation point a such that {a} remains in T(G) after pruning, we can set a ∈ ¯S. For reasons of space, we do not give the pseudo code for the procedure prune tree. It iteratively removes nodes of degree 1 whose label contains no vertex in ¯S.

Proof: The correctness of the pruning is straightforward for vertices that are not connected to G[ ¯ S] (Line 2), and proved by Lemma 1 for articulation points (Line 7). It remains to show completeness. We consider a vertex v ∈ ¯ S \ ¯S and show that there exists a solution where v ∈ S and another where v ̸∈ S. Since ¯ S is a single connected component, it is a solution and it accounts for the case v ∈ S. Now consider first the case where v is not an articulation point. Then v belongs to a single biconnected component p of G[ ¯ S]). We can exclude v from S and by definition of a biconnected component, connect all remaining vertices. Therefore ¯ S \ {v} is a solution. Lastly, suppose that v is an articulation point. Since it has not been added to ¯S in Line 7 then it must have been removed from T(G). Moreover, the union of the vertices of the residual tree is a superset of ¯S that are connected in G. Therefore it is a solution where v ̸∈ S. □

3.2 Explanation

We need to explain two types of pruning: First, in Line 2 of Algorithm 1, when the current maximal subgraph is not connected, we arbitrarily choose a node s from one connected component and exclude the nodes of all other components from S. Given a node t from another component, we compute a set of nodes F such that F ∩ ¯ S = ∅ and every path linking s to t has a node in F . In order to do so, we do a breadth first search from s stopping at nodes in V \ ¯ S. Let B denote the set of leaves of this BFS. We then perform another BFS from t until we reach a node in B. The set of leaves of this second BFS is F. The explanation clause is then:

Second, in Line 7 of Algorithm 1, articulation points are added to ¯S. The explanation is similar, except that we start from an articulation point v. We arbitrarily choose two nodes s, t ∈ ¯S by exploring the tree T(G) breadth first from {v}. Then the first BFS in G stops when exploring nodes in (V \ ¯ S) ∪ {v} to return a set of leaves B. The second BFS is unchanged, and the explanation is :

Theorem 2 Explanations 1 and 2 are correct and minimal.

Proof: All paths from s to t go through a vertex in F, so if these vertices are removed from G, then s and t are disconnected, hence explanation 1 is correct. The argument for explanation 2 is nearly the same. By definition of articulation points, any path from s to t in G goes through a node in F or through v, hence removing F ∪ {v} disconnects s and t in G.

Now consider any vertex u ∈ F, (resp. u ∈ F ∪ {v}). If we remove from G all vertices in F (resp. F ∪{v}) except u, there exists a path from s to t in G. This path is a support for the constraint, hence explanations 1 and 2 are minimal. □

4 Weighted connectivity

We now consider the weighted version of the constraint:

Definition 2 (WEIGHTED-CONNECTIVITY) Let G be a graph, S a set variable, W an integer variable and ω : V ↦→ N a mapping from vertices to non-negative integers.

Deciding the existence of a support in the constraint WEIGHTED-CONNECTIVITY is NP-complete as it solves the STEINER TREE problem. Hence, enforcing DC is NPhard, as NP-hardness of support implies NP-hardness of DC [Bessiere et al., 2007]. Therefore, we do not try to achieve DC, but propose an incomplete algorithm taking into account the cost of the connected component.

4.1 Propagation

The main difference with Algorithm 1 lies in the first step (Line 1). Instead of computing a connected component, we compute the subset of vertices for which there exists a path to every vertex in ¯S that is short enough (see Fig. 2a and 2b).

Moreover, we bound the minimum value min(W) by ∑v∈¯S ω(v). From now on, we consider the graph G[ ¯ S]. The length of a path in this graph is defined as the sum of the weights of the vertices in this path that are not in ¯S. It follows that if a vertex v is in the subgraph induced by S, then the total weight of this subgraph is at least min(W) (the weights of vertices in ¯S) plus the maximum over all minimum lengths of paths between v and any element of ¯S. We run Dijkstra's algorithm once for each vertex u ∈ ¯S to compute the shortest path from u to all other vertices, which can be done in O(|E|+|V | log |V |). 1 Let s(u → v) be the length of the shortest path to vertex v. For every vertex v such that s(u → v) + minW > max(W), we remove v from ¯ S. This step can thus be done in O(|V ||E| + |V |2 log |V |) time.

1 We used a binary heap hence a O(|E| log |V |) time complexity.

4.2 Explanation

In order to explain this pruning, we use an algorithm similar as in the unweighted case, however again using Dijkstra's algorithm instead of a simple BFS. Let u ∈ ¯S, v ∈ ¯ S be two vertices such that the length s(u → v) of the shortest path in Gw.r.t. the weight function ω is greater than max(W).

We compute a cut F defined as follows: F ∩ ¯ S = ∅ and every path in G from u to v of length less than or equal to max(W) - min(W) has a vertex in F. We first run Dijkstra's algorithm on G[ ¯ S ∪B], where B is the set of vertices in G\G[ ¯ S] adjacent to vertices in G[ ¯ S] to compute the shortest paths from u to every node in B. Next, we run Dijkstra's algorithm again, however from v and on the graph G. In this run, whenever we reach a node w ∈ B, we add it into the cut F if s(u → w) + s(v → w) - ω(w) +min(W) ≤ max(W). When w is added to the cut we do not explore further through that vertex. Moreover, since the reachability depends on the maximum cost, we must add the literal "W > max(W)".

The final explanation is W > max(W) ∨ u ̸∈ S ∨ v ̸∈ S ∨ ∨w∈F w ∈ S (see Fig.2c). As in the unweighted case, explaining the second type of pruning is done in the same way: by adding the articulation point to the set B.

The amortized worst case time complexity to explain this pruning is then O(|V ||E| + |V |2 log |V |) down a branch.

5 Multiple Subgraphs

Finally we consider the case where we want to ensure that several disjoint subgraphs are connected.

Definition 3 (MULTICONNECTIVITY) Let G = (V,E) be a graph, {Xi | i ∈ V } a set of variables and R a set of integers.

Clearly we can decompose MULTICONNECTIVITY using CONNECTIVITY, but it hinders propagation since enforcing DC on MULTICONNECTIVITY is NP-hard [Karp, 1975].

5.1 General graphs

We study the following problem, a special case of finding disjoint paths which is not covered by the results in [Karp, 1975] or [Fortune et al., 1980].

NAME: TWO-SHAPE-CONNECTIVITY

INPUT: A connected graph G = 〈V,E〉 and a labeling of the vertices m : V → {0, 1, U}

PROBLEM: Does there exist a labeling m′ : V → {0, 1, U} of the vertices so that m′ (v) = m(v) if m(v) ∈ {0, 1} and the components C0, C1 with Ci = {v | m′ (v) = i} are connected?

Let m[i] be the domain of a variable Xi, it can be modeled as: MULTICONNECTIVITY[G, {0, 1}](Xi | i ∈ V )

Theorem 3 TWO-SHAPE-CONNECTIVITY is NP-hard.

Proof: We reduce from an instance of 3-SAT with n variables and m clauses. We first introduce two gadgets, ̸=(v1, v2) and ∨(v1, v2, v3). The ̸=(v1, v2) gadget is a complete bipartite graph with partitions {l0, l1, r0, r1}, {v1, v2}. We have m(l0) = m(r0) = 0, m(l1) = m(r1) = 1 and m(v1) = m(v2) = U. This gadget ensures that m′ (v1) ̸= m′ (v2). Suppose wlog m′ (v1) = m′ (v2) = 1. Then l0 is disconnected from r0, which is not permitted. For the same reason, vertices v1, v2 cannot be both labeled 0. In order to combine many ̸= gadgets in one graph, we need to make sure that the resulting C0, C1 partitions are connected. To do this, we introduce globally two distinguished vertices g0, g1 with m(g0) = 0, m(g1) = 1 and connect g0 with the vertex l0 and g1 with l1 in each ̸= gadget. The ∨(v1, v2, v3) gadget ensures that m(vi) = 1 for at least one of v1, v2, v3. It is a complete bipartite graph with partitions {l, r} and {v1, v2, v3}

and m(l) = m(r) = 1. In order for l and r to be connected, at least one of v1, v2, v3 must be in the component C1.

From the 3-SAT formula, we construct a graph with n instances of the ̸= gadget and m instances of the ∨ gadget. In particular, for each variable xi, we have two vertices xi and xi and the gadget ̸=(xi, xi). For each clause (l1, l2, l3), we have an instance of the gadget ∨(l1, l2, l3), where li is a literal, i.e. either xi or xi. From a labeling m′ we construct a solution of the 3-SAT formula with xi = 1 ⇐⇒ m′ (xi) = 1. Because of the ∨ gadgets, at least one of the literals of each clause will be assigned 1, hence the formula will be satisfied. In the other direction, we use the same rule to construct a labeling from a solution of the 3-SAT formula. By construction, the vertices in the ̸= gadgets will be connected. Since each clause is satisfied, one of the literal vertices of each ∨ gadget will be labeled with 1. Hence, the vertices l and r of that gadget will be connected to each other through that literal vertex and, since that literal vertex appears in a ̸= gadget, they are also connected to the rest of the C1 partition. □

5.2 Two Shapes in Planar Graphs

The conjunction of an arbitrary number of CONNECTIVITY constraints is known to be NP-hard even on planar graphs [Lynch, 1975], but fixed parameter tractable on the number of terminals [Reed et al., 1991]. It is not known whether this remains the case for two constraints on a planar graph. We show here that decomposition into two CONNECTIVITY constraints for the case of two shapes will hinder propagation. Consider the graphs in figure 3. In both examples the connectivity constraints are DC, but the conjunction is disentailed. In the example on the left, the border is assigned with pattern 1/0/1/0. This means that traversing the border we see the label 1, then 0, then 1 and then 0 again, potentially with repetitions or with unlabeled vertices in between. In the example on the right, the same pattern appears in the cycle v9, vA, vB, vC.

In order to detect disentailment, one must forbid the pattern 1/0/1/0 in every face of a planar graph. This can be done, for example, with a REGULAR constraint forbidding two alternations between 1 and 0, which requires 8 states, hence can be enforced in linear time in the size of the face. Since an edge belongs to at most two faces, enforcing DC on all face constraints can be done in time O(|E|). However, having the face constraints is not sufficient for detecting disentailment. We do not show counterexamples for lack of space.

6 Experimental Results

In this section we report the results of a series of experiments on three benchmarks involving connectivity constraints.

Wildlife corridor First, we test a hybrid CP/SAT approach using the constraint WEIGHTED-CONNECTIVITY for the wildlife corridor problem described in [Conrad et al., 2012]. We implemented all techniques with minicsp 2 and ran experiments on a cluster of 48-core Opteron 6176 nodes at 2.3 GHz with 378 GB RAM available. Given a graph representing a region, where vertices stand for parcels, a subset of parcels are wildlife reserves. The goal is to design a corridor so that all reserves are connected while minimizing a (positive) cost and/or maximizing a (positive) utility. In [Conrad et al., 2012], the problem is solved in two phases. First, a Steiner tree of minimal cost is obtained by computing all pairs shortest paths and selecting those connecting each reserves to the best "center point". This method is very efficient, however it is optimal only for three reserves or less. Second, the utility maximization problem is solved with a MIP, using the Steiner tree above to either guide search, or to fix the selected parcels into the solution.

Let G = (V,E) be a graph describing the region with |V | = n, {r1, . . . , rk} ⊆ V a set of reserves and ω a cost function for parcels. We have a single set variable S such that ¯S = {r1, . . . , rk} and ¯ S = V , an integer variable W with domain [0, . . . , ∑v∈V ω(v)]. The first phase corresponds to minimizing W subject to a single constraint: WEIGHTED-CONNECTIVITY[G, ω](S, W ).

For the second phase, the same constraint can be used to ensure that the corridor is connected and that its cost is below a given budget (derived from the min cost solution by adding a certain slack). The utility can be maximized with a simple inequality constraint. However, this is not competitive with Gomes et al. when the slack grows. Indeed, for a large slack, the WEIGHTED-CONNECTIVITY is not so critical and the communication with the inequality constraint taking care of the utility is not as good as it is in MIP. However, our approach can be useful in the first phase when we have more than three reserves. We used the generator of Gomes et al. to create 20 random instances of 24 classes parameterized by the order of the grid graph (in [8, 15]) and the number of reserves (in {5, 7, 10}). We report in Table 1, for each class of

2 http://www.inra.fr/mia/T/katsirelos/minicsp.html

20 instances, the improvement of the cost with respect to the upper bound found by Gomes et al. (imp), the percentage of instances for which optimality was proven (opt), and the average cpu time for computing the proofs when it was possible within the cutoff of 300 seconds (time).

We can clearly see that the upper bound can be significantly improved. However, the computational time grows both with the number of terminals and with the size of the graph.

Ecosystem design. In order to compare the various algorithms that we have introduced, we use a problem similar to wildlife corridor design. However, besides reserves, we also have parcels that correspond to agricultural activities (farms). We want to compute a minimum cost connected wildlife corridor, but also to make sure that farms can be connected by roads. We therefore have two set variables S and T standing respectively for the set of parcels that are in the corridor, or witness of a connection between farms. We first ensure that a parcel can only account for one of these two cases, then we post the following constraints: CONNECTIVITY[G](T) and WEIGHTED-CONNECTIVITY[G, ω](S,W). The former ensures that the farms can be connected whilst the latter ensures that the corridor is connected and links its cost with W. The objective is to minimize the total cost of the corridor W.

We use 20 random instances of 12 classes, parameterized by number of reserves and farms (10/5 and 5/10), and by the order of the grid in [8, 13]. We compare four models. The connectivity and cost minimization of the corridor is either decomposed with a CONNECTIVITY and a linear inequality (connect) or modeled directly with the weighted version (weighted). Then, in both cases, we try with or without face constraints (denoted "+face"). We also report results without clause learning (denoted by "∗ "). We omit comparison against a decompositon as it was several orders of magnitude slower even in very small instances. We report the same data in Table 2 as in Table 1, except for the objective value (obj.). The first observation is that learning has a large positive impact. In terms of number of instances solved to (proven) optimality, the strongest propagation method (weighted+face) is the best. However, when the size of the instances grows and finding the optimal solution within the time cutoff becomes difficult, the lower time complexity of the non-weighted variant can be more beneficial than the extra pruning of the weighted variant. The face constraints, however, are cheap and nearly always beneficial.

Clueless Puzzles Clueless Puzzles were proposed by Butters, Henle, Henle, & McGaughey as a mathematical challenge. 3 . The goal is to design Sudoku-like puzzles without numerical clues. The answer must be a latin square. However, there are no pre-filled cells and no extra ALLDIFFERENT constraints. Instead, the grid is partitioned into regions of equal sum (see examples in Figure 4). Since the solution must form a latin square, we can deduce that, for a puzzle of order n the sum in each of the m regions equals n2 (n + 1)/(2m). We model the problem of generating puzzles as follows:

3 http://www.math.smith.edu/∼jhenle/clueless

To ensure that the puzzle has a solution we use n2 integer variables with domain {1, . . . , n}, and 2n ALLDIFFERENT constraints to ensure we have a latin square. Next, to obtain a partition into regions, we use m set variables whose elements are the cells of the grid. We post that the scalar product of the characteristic function of each set variable with the latin square is equal to n2 (n + 1)/(2m). Moreover, they must be disjoint, their union must cover the whole grid, and of course they must be connected. Notice that this model does not guarantee that the puzzle will have a single solution, which is a necessary feature. Indeed, this problem is not in NP, and thus is unlikely to be formulated as a simple CSP. We therefore check this property once a feasible puzzle is generated (by attempting to find a different solution), and continue searching when this does not hold. Using this simple program we were

able to generate several new instances, for classes (〈n,m〉) of puzzles for which known solutions already exist (such as 〈4, 2〉, 〈5, 2〉, 〈5, 5〉, 〈6, 9〉 and 〈6, 14〉), but also for entirely new classes: 〈6, 6〉, 〈6, 7〉 and 〈7, 7〉 (see examples in Fig. 4).

7 Conclusions

We have studied a family of constraints to reason about connectivity that are useful in a range of problems in computational sustainability. First, we showed that one can enforce domain consistency and compute a minimal explanation for the pruning in O(|E|) time for a basic connectivity constraint. Second, we proposed an incomplete algorithm to take weights into account since enforcing DC in this case is NP-hard. Third, we showed that partitioning into several connected subgraphs is intractable, but proposed efficient implied constraints for the case of two subgraphs of a planar graph. Finally, we empirically tested our algorithms on several benchmarks and showed that they can be practical and have promise for a wide range of problems.

Acknowledgements

We would like to thank Ashish Sabharwal for his help with the wildlife corridor datasets [Conrad et al., 2007].