An Algebraic Approach to Constraint Satisfaction Problems

1 Introduction

We are interested in solving certain search problems that arise in combinatorics and in artificial intelligence. These problems are called Constraint Satisfaction Prob­ lems (CSP's), and include such familiar tasks as graph­ coloring and the n-queens problem. Typically, such problems are solved via exhaustive search. We propose an alternative method, based on converting the CSP into an integer linear programming problem, which can be solved via polynomial multiplication.

In section 2 we describe constraint satisfaction prob­ lems and show how they can be converted into inte­ ger linear programming problems. We describe how to solve integer linear programming problems via polyno­ mial multiplication in section 3, and analyze the time requirements of this approach. In section 4 we give an example of our method, applying it to the n-rooks prob­ lem, and show some surprising similarities with standard search techniques. In section 5 we address the problem of finding a good integer linear programming equivalent of a CSP, while in section 6 we describe some extensions to our techniques. Finally, in section 7, we derive new simply exponential bounds for the n-queens problem and the toroidal n-queens problem.

Ramin Zabih Computer Science Department Stanford University Stanford, California 94305

2 CSP's and integer linear programming

2.1 Constraint satisfaction problems

Formally, a CSP has a set of variables and a domain of values,

Every variable X a must be assigned a value di; such an assignment will be denoted v, -di .

A CSP also consists of some constraints saying which assignments of values to variables are compatible. We will restrict our attention to binary CSP's, where the constraints involve only two variables. (This is partly for simplicity, and partly because a non-binary CSP can be converted into a binary CSP by the introduction of additional variables.) A solution to the constraint satis­ faction problem is an n-tuple of assignments that satis­ fies all the constraints.

There are several ways to formulate the task of solving a CSP.

The satisfiability problem is to determine the exis­ tence of a solution to a CSP.

The counting problem is to determine the number of solutions.

The enumeration problem is to find all solutions.

The counting problem is a special case of the enumera­ tion problem, and the satisfiability problem is a special case of the counting problem.

Since many NP-complete problems, such as graph­ coloring, are CSP's, the satisfiability problem is easily seen to be NP-complete [Garey and Johnson, 1979]. We therefore cannot expect to solve the satisfiability, count­ ing or enumeration problems in polynomial time. Typi­ cal CSP solution techniques like those surveyed in [Mack­ worth, 1987] perform an exhaustive search, backtracking through all n-tuples of assignments. While these meth­ ods add some additional cleverness beyond brute-force search, their worst-case running time is still O(mn).

We will concentrate on solving the counting problem (and therefore the satisfiability problem as well). While our methods are easily extended to the enumeration problem (see section 6.1), our approach to the counting problem is surprisingly simple to analyze. In fact, we can do provably better than backtrack search for some problems, in the following sense. Backtracking counts the number of solutions by enumerating them, and hence must take at least as much time as the number of solu­ tions. Our techniques can determine the number of solu­ tions without actually enumerating them, so the number of solutions is not a lower bound on our running time.

2.2 Converting CSP's into integer linear programming problems

Our basic approach is to turn a CSP into an integer linear programming problem. For every possible assign­ ment Vi+-dj we introduce a variable Xi,j. Then for every CSP variable v; we have an equation

This expresses the fact that every CSP variable has ex­ actly one value. In addition, for every set of assign­ ments that the constraints prohibit, we have an inequal­ ity which states that no more than one of these assign­ ments can hold at once. For example, if the pair of assignments {vi +-d;, vk +-dl} is forbidden by the con­ straints, we will have the inequality

Now consider the problem of finding all solutions to the equations of the form (I) and (2) over the non-negative integers. (It is clear that any solution will assign every variable x;,; either 0 or l.) This is an integer linear programming problem, which corresponds to the original CSP. There are nm variables x;,;, and a solution to the set of equations with x1,j = I corresponds to a solution to the CSP containing the assignment v;+-d; .

3 Solving CSP's by multiplying polynomials

An integer linear programming problem is a set of lin­ ear Diophantine inequalities. We can use an approach based on generating functions to determine the num­ ber of solutions. This involves multiplying out a certain polynomial, and gives us a way of solving the original CSP.

3.1 Solving linear Diophantine equations

A set of M linear Diophantine equation in N unknowns has the form

The Xj 's are the unknowns, and the w; ,j 's and Si's are non-negative integers. The problem at hand is to de­ termine the number of solutions over the non-negative integers.

Consider the generating function

The formal parameters of this function are the variables Y1, ... , Y M , and on such variable is associated with each equation. We need one non-trivial bit of mathematics: the number of solutions to equation (3) is the coefficient of

in <I>. (This fact was proved in [Euler, 17 48], and is discussed in [Schrijver, 1986, pages 375-376].)

3.2 Applications to integer linear programming

Now consider the integer linear programming problem corresponding to a CSP. Every w;,; will be either 0 or 1. Assume for the moment that all the equations are of the form (1), hence s; = l. (We will relax this restric­ tion shortly.) We can count the number of solutions by finding the coefficient of the monomial Y = Y1 Y2 · · · Y M in <I>. We apply the identity

to equation ( 4), giving

We can further simplify this if we remember that we are only interested in the coefficient of Y . We can there­ fore ignore any term divisible by the square of any Y; (terms that are not so divisible are called square free). By a slight abuse of notation we will refer to the square free part of the generating function as <I> . We next note that

where these latter terms are not square free and so can be discarded. This gives us an expression for <I> as a product of binomials, namely

All we need do in order to find the number of solutions to the system of equations is compute the coefficient of Y in equation (7).

The obvious way to solve this problem is to simply multiply out <I> and look at the correct term of the result . We can do this iteratively: we multiply the first two binomials of <I> and throw away the terms that are not square free, then repeat with the next binomial, etc.

We can now relax our assumption that all the equa­ tions are of the form (1) . Suppose that instead we have exactly one Diophantine inequality

This implicitly defines two equalities, one of which has 0 on the right hand side, and one of which has 1. The number of solutions to the inequality is the sum of the number of solutions to the two implicit equalities.

At first glance, it would seem that solving a system with a inequalities would require solving 2a implicit sys­ tems of equalities. However, a closer look at the gener­ ating function shows that we only need to solve one of these implicit systems. Suppose that we replace each in­ equality with an equality that has a 1 on the right hand side and examine the square free part of~ - The number of solutions to each of the 2a implicit systems of equal­ ities will be the coefficient of the appropriate term of ~-

For example, suppose that there are three equations Yo, Y1 and Y2, and one inequality

where ~ is a sum. The inequality (9) implicitly asserts

The number of solutions to the original system (which included (9)) is the sum of the number of solutions when (9) is replaced by (9.0) plus the number of solutions when (9) is replaced by (9.1).

However, suppose that we simply replace (9) by (9.1) and calculate ~. Let Y3 be the generating function vari­ able associated with (9.1) . The number of solutions is the coefficient of YoY1Y2Ya in~- However, the number of solutions when (9) is replaced by (9.0) is also part of~. namely the coefficient of Y0Y1Y2. We thus do not need to calculate separate generating functions for both systems of equalities implicitly defined by (9). We can simply replace the inequality by an equality with 1 on the right hand side and calculate~- At the end, we add together coefficients for the implicit systems of equali­ ties to determine the number of solutions to the original problem.

To summarize, we can replace the a inequalities with equalities and calculate the square free part of~- The total number of solutions will no longer be the coefficient of Y, but rather the sum of 2a coefficients of different terms. This justifies our earlier assumption that the in­ teger linear programming problem consisted entirely of equalities of the form (1).

3.3 Time requirements

The time required to solve the CSP is the time to mul­ tiply out the polynomial ~- We can do this in N itera­ tions, where on each iteration we multiply a binomial by a polynomial that represents the product thus far. Two polynomials can be multiplied in time proportional to the product of their size, so we need to bound the size of the intermediate results.

Naively, there could be as many as N M terms, since there are MY's, and the degree of each one is potentially as high as N. However, remember that at each stage we will discard all the terms that are not square free.

or

(In algebraic terminology, this is described as doing the computation over the polynomial ring

instead of Z[Y1 , ... , Y M) .) This gives a better bound on the number of terms. Every Y will have degree 0 or 1, hence there are no more than 2M terms.

Summarizing, if we compute ~ iteratively we will do N multiplications, where each multiplication computes the product of a binomial and a polynomial with no more than 2M terms. Each multiplication will take time 2M +1 , which is 0(2M). This puts a bound of O(N2M) on the total time complexity of our algorithm. 1 We will show in section 6.2 how to reduce the exponent to M - n in the worst case.

It is worth noting that we cannot expect to find a polynomial time method for determining the value of an arbitrary coefficient of~- [Schrijver, 1986] proved that integer linear programming is NP-complete even when the coefficients involved are restricted to be 0 or 1.

4 An example of our approach

The n-rooks problem consists of placing n rooks on an n-by-n chessboard so no two rooks attack each other. While it is clear that the solutions to this problem are the n! permutations of n elements, this problem has a particularly simple encoding as an integer linear pro­ gramming problem, and is therefore a useful example of our approach.

There are n 2 assignments v;~dj for this CSP, so we will convert it into a programming problem with n 2 vari­ ables Xi,j . There will be 2n equations. 2 The equations specify that there is exactly one rook in each row and each column. A typical equation is

which encodes the fact that there is exactly one rook in the 2nd column. We can number these 2n equations so that equation i says that there is exactly one rook in the i'th column, and equation n+ i says that there is exactly one rook in the i'th row. The formal parameter Y; of the generating function corresponds to the i'th constraint equation.

Applying equation (7), ~ is a product of n 2 binomi­ als. There will be one binomial for every square on the chessboard, and the binomial for the square in row i and column j will be

The number of solutions is the coefficient of

in the expansion of ~ .

1 This analysis ignores the time required for scalar multi­ plication. We can take this into account by noting that the size of the coeflicien ts will be bounded by M N, hence their length is O(N log M). So the scalar multiplication can be done in time O(Nlog M log(N log M) loglog(N log M)) [Aho et al., 1974].

2 In section 5 we will discuss how to find good encodings.

Surprisingly, there is a way of multiplying out ~ that corresponds closely to backtrack search. Suppose that we attempt to find all combinations of binomial terms whose product is Y . We can start by picking the second (i.e., non-constant) term of (1 + Yn+i }j ). This forces us to pick the first (constant) term of every binomial whose second term includes either Yn+i or }j, to en­ sure the product is square free. Now we can pick the non-constant term of some other binomial and continue, until we have produced V. The number of such com­ binations will be the coefficient of Y in the expansion of~. (It is also possible that the combination of terms selected guarantees that no selection of terms from the remaining binomials will produce Y. The simplest ex­ ample of this occurs when a term has been selected from all but one binomial without producing Y, where the re­ maining binomial's non-constant term does not contain the missing Y. This phenomenon, which does not oc­ cur in the n-rooks problem, will be used to improve our performance guarantee in section 6.2.)

This whole process corresponds closely to a simple search of the original CSP. Choosing the non-constant term of the binomial (1 + Yn+i }j) places a rook on the equivalent square, while choosing the constant term means there is no rook on that square. The process described above corresponds to placing a rook on the chessboard square in row i and column j, and thus en­ suring there is no rook in any other square in that row or column. At the end, we will have satisfied every con­ straint and produced a solution. The alternative out­ come, where we cannot satisfy every constraint, corre­ sponds to a failure in backtrack search.

Note that the expansion of~ does not explicitly repre­ sent the solutions. In fact, which assignments constitute a particular solution is lost by this approach, as the iden­ tity of the binomials is not recorded. We will show how to explicitly represent solutions in section 6.1. However, the fact that solutions are not explicitly represented is an advantage of our approach, as the number of solutions is no longer a lower bound to our running time. Apply­ ing the analysis of section 3.3, for this problem we have N = n 2 and M = 2n, so our running time is O(n 2 4n). This is significantly smaller than n!, the number of solu­ tions, which is about 4°(n logn).

We have just described a corre~ondence between a way to compute the coefficient of Y in the expansion of ~ and a simple search technique for solving the origi­ nal CSP. However, searching all the combinations of bi­ nomial terms is not our intended method. We believe it would be better to multiply ~ out iteratively, as de­ scribed earlier. While this does not correspond in any obvious way to applying backtrack search to the orig­ inal problem, it does seem to have some similarity to breadth-first techniques. It may thus be susceptible to parallel implementation techniques like those described in [Dixon and de Kleer, 1988].

5 Finding a good encoding

A potential difficulty with our approach is that the run­ ning time is exponential in M, the number of equations needed to encode the CSP. There can be many such enco dings, and it is important to be able to find a good one. Clearly M must be at least n, as there is an equation of the form (1) for each CSP variable. The remaining equa­ tions come from the constraints. The problem is to find a way of encoding the constraints as small number of equations.

In the simplest encoding, there is an equation like (9) for every pair of assignments that the constraints rule out. However, this can be quite a lot of equations, as the following example illustrates. In the case where all pairs of assignments are ruled out this will result in n 2 m 2 equations, while, as we will see, there is an encoding that results in only 1 equation.

Let us define the assignment graph corresponding to a particular CSP. The nodes of the assignment graph rep­ resent assignments, and there is an edge between two nodes if their assignments are ruled out by the con­ straints. A clique C in this graph corresponds to a set of pairwise incompatible assignments. The transformation described in section 2.2 will produce a variable x;,j for every node, and will use (~) equations, with the equation

expressing the fact that these two assignments are in­ compatible. However, we can describe the fact that all the assignments in C are pairwise incompatible with the single equation

So we really only need one equation for every clique. If all pairs of assignments are incompatible, the entire assignment graph is a clique, and so the constraints can be represented with 1 equation instead of needing n 2 m 2 .

For the n-queens problem, the assignment graph has a node for each square on the chessboard. Two nodes are connected if queens on the two corresponding squares would attack each other. The maximal cliques corre­ spond to the columns, rows, diagonals and antidiagonals of the chessboard.

We can therefore reduce the size of M by encoding ev­ ery clique as a single equation. In an optimal encoding, M can be as small as the size of the smallest covering by cliques of the assignment graph. Because the problem of covering a graph by cliques is N1'-complete [Garey and Johnson, 1979J, finding an optimal encoding would take exponential time. For special cases like the chessboard problems described above, it is easy to find an optimal encoding, but in general we would have to settle for a sub-optimal encoding. It is also possible to phrase the problem as finding a small vertex cover for the assign­ ment graph, and this can also reduce M.

5.1 Which CSP's have good encodings?

One obvious question is whether there is a better charac­ terization of the CSP's that have good encodings. While we do not have a very good answer to this question, an observation due to Bob Floyd suggests that there may be a relationship between the number of solutions to the CSP and the smallest possible value for M. Floyd has pointed out that a solution to the CSP corresponds to an n-clique in the complement of the assignment graph (equivalently, an independent set of size n in the assign­ ment graph). It is possible that work in combinatorics will produce a relationship between the number of inde­ pendent sets of fixed size in a graph and the size of the smallest vertex cover . Such a relationship would produce a characterization of the CSP's where our algorithm can be expected to work well, in terms of the number of so­ lutions.

6 Some extensions

The method described above relies on multiplying out the polynomial q_) defined in equation (7). A simple mod­ ification to q_) allows us to solve the enumeration prob­ lem instead of the counting problem. It is also possible to reduce the running time of our algorithm by carefully arranging the order in which the binomials comprising ~ are multiplied out.

6.1 Solving the enumeration problem

Suppose that instead of determining the number of solu­ tions to our CSP, we are interested in enumerating these solutions. We can still convert the CSP into an integer linear programming problem, and solve it via polynomial multiplication. The only difference is that we now need to consider a slightly different function than ~ -

Recall that we start by converting our CSP into a set of equations in N = mn variables, one variable per as­ signment Vi+-dj. Let us introduce N different variables CXj, 1 ~ j ~ N. Now consider a slight variation of~' namely

We can solve the enumeration problem by determining the coefficient of Y in the expansion of q_) 1 , just as we did with q_> _

Suppose we multiply out ~1 by selecting terms from different binomials. The aj's will keep track of which binomials had their non-constant side selected. When we have a solution, the coefficient of Y will the prod­ uct of a number of the ai 's, which will encode a solu­ tion. When q_) 1 is multiplied out completely (iteratively or otherwise) , the coefficient ofY will not be an integer, but rather the sum of a number of different products of aj's. Every such product is a solution to the CSP, and all solutions appear as such products.

6.2 Reducing the running time

An observation due to John Lamping can be used to reduce the running time of our algorithm from O(nm2M) to O(nm2M-n) . Suppose that before we multiply out q_) we first try to identify a subset S of the binomials comprising ~ and a variable Y such that

• Y never appears outside of the binomials in S, and
• some other variable Y' does not appear in S .

If we had such an S and Y, we could effectively reduce the number of variables by 1.

We can do this by first computing the product of the binomials in S . This can be grouped into terms based on the highest power of Y that divides each term

where Si is not divisible by Y, and as usual we can dis­ card the non square free terms. If Y is contained in the coefficient(s) corresponding to solutions to the CSP, then we can replace S by S1. We can then forget about Y, as it never appears in the remaining part of~. Since Y' did not appear in S, this reduces the number of variables by 1.

Reducing the number of variables reduces the running time of our algorithm as well, which is exponential in the number of variables. It turns out that we are guaranteed that there exist at least n such pairs of S's andY's. This enables us to reduce the number of variables from M to M - n, which in turn reduces our running time to O(nm2M-n).

There are n such S and Y pairs because the encod­ ing of a CSP will have a single equation for each CSP variable, representing the fact that the variable must be assigned exactly one value. 3 Let Y correspond to the equation which constrains the CSP variable v to have exactly one value, and consider the binomials in q_) that contain Y. These are exactly the binomials that cor­ respond to possible values for v , so there can be only m of them. Furthermore, this set of binomials forms an appropriate S, since no Y' corresponding to another CSP variable can appear. There are n such Y's, one per CSP variable, and these Y's will always appear in the coefficient(s) that we need to calculate.

It is even possible to take advantage of such S's and Y 's when Y does not necessarily appear in the required coefficients. (Recall that this is a result of the integer linear programming equation that corresponds to Y ac­ tually representing an inequality, as discussed at the end of section 3.2.) If this happens, we can replace S by So + S1 and still discard Y . This has the effect of sum­ ming the coefficients that stand for the solutions to the two implicit systems of equalities early in the computa­ tion, rather than waiting to multiply q_) out before cal­ culating this sum.

We can thus reduce the number of variables by at least n . As we saw in section 4, each binomial corresponds to a possible assignment Vi +-di. The ordering that is guaranteed to reduce the number of variables by n cor­ responds to considering all the assignments for a fixed Vi before moving on to the next CSP variable. On the n-queens problem this is equivalent to multiplying out the binomials in row-major order.

7 Some new bounds

An additional advantage of our approach is that we can produce non-trivial upper bounds for solving certain con­ straint satisfaction problems. So far, we have confined out attention to CSP's based on placing pieces on a chessboard.

3 An encoding that does not have this property corre­ sponds to a CSP whose unsatisfiability can be easily detected via arc consistency [Mackworth, 1977] .

The classic such problem is the n-queens problem, which was known to Gauss. The n-queens problem con­ sists of placing n queens on an n-by-n chessboard so that no two queens attack. A variation of this problem is the toroidal n-queens problem, where lines of attack between queens are considered to 'wrap' as if the board were a torus. The toroidal n-queens problem has strictly fewer solutions than the (regular) n-queens problem.

To our knowledge there is no known bound for CSP search algorithms for this problem, beyond the trivial bound of nn. In fact, there is some evidence that the number of solutions to the toroidal problem (and hence to the regular problem) is larger than exponential [Rivin and Vardi, 1989].

The toroidal n-queens problem can be formulated as a simple set of equations. There are n 2 variables Xj. There are also 4n equations, one for each column, row, diagonal (top to bottom) and anti-diagonal . So we have M = 4n, N = n 2 · Applying the above results, our algorithm can determine the number of solutions in time O(n 2 8n).

For the regular n-queens problem, again N = n 2 , but this time there are 2(2n- 1) inequalities for the diago­ nals. This gives us M = 6n- 2, and a running time of O(n 2 32n) .

However, if we do the multiplication in row-major or­ der, it turns out that only 3n variables will need to be active at once. This can be seen by noting that the top row involves 3n constraints. Moving to the next row adds a new row constraint and removes the old row con­ straint. The far left square of the new row eliminates an antidiagonal constraint (all the squares involved in that constraint have been eliminated), and adds a new diagonal constraint, while the far right square eliminates a diagonal constraint and adds a new antidiagonal con­ straint. This gives us a much better bound of 0( n 2 8n ). Our preliminary experiments with this algorithm have been promising.

8 Conclusions

We have shown a surprising connection between solv­ ing search problems and multiplying certain polynomi­ als. The performance is very simple to analyze, and we can show new bounds for several constraint satisfaction problems.

Our approach also holds our the possibility of apply­ ing sophisticated algebraic techniques to constraint sat­ isfaction problems. For example, it is possible that the coefficient of Y in ci> can be estimated, thus providing an approximate count of the number of solutions. Also, complex mathematical results like those in [Anshel and Goldfeld, 1988] can be used to provide an upper bound on the number of solutions by studying the generating function directly, using the theory of modular forms.

Our algorithm also shares some surprising properties with Seidel's method [Seidel, 1981]. Like Seidel, we pro­ duce a method which is easy to analyze, and which gives a bound that is exponential in a certain parameter of the problem. Seidel's parameter f is a property of the topology of the constraint graph, which can be shown to be 0( yin) for planar graphs. Our parameter M is a property of the way the problem gets encoded as an integer linear programming problem, which can be shown to be linear for variants of the n-queens problem.

Seidel faces the problem of making f small, which in general requires exponential time. Similarly, we face the problem of making M small, which requires solving an NP-hard problem involving graph cliques. The simi­ larities between our algorithm and Seidel's, as well as Lamping's observation described in section 6.2, suggest that our algorithm may be viewed as dynamic pro~ram­ ming. We are currently exploring this possibility lRivin and Zabih, 1989].

8.1 Acknowledgements

We are indebted to John Lamping and Ilan Vardi for use­ ful suggestions. Support for this research was provided by DARPA under contract number N00039-84-C-0211. Ramin Zabih is supported by a fellowship from the Fan­ nie and John Hertz foundation.