This paper shows how the solutions of constraint satisfaction problems that involve only polynomial constraints over nite domains can be enumerated by computing the values of related polynomial functions at appropriate points. The proposed algorithm rst transforms constraints, which are expressed as equalities, inequalities, and disequalities of polynomials with integer coe cients and integer variables, into a canonical form that uses only inequalities. Then, starting from a bounding box, which is supposed to be known, the algorithm recursively subdivides the box into disjoint boxes and it records boxes whose elements satisfy all constraints. The subdivision is driven by the study of the sign of polynomial functions over boxes, which is performed by means of a method that uses only the coe cients of polynomials and the values of functions at the corners of boxes.
Polynomial constraints are ubiquitous, and they nd relevant applications in
virtually all elds of science and engineering (see, e.g., [
We have already presented results on the use of the properties of Bernstein
polynomials [
This paper is organised as follows. Section 2 recalls adopted canonical form of constraints and it outlines the subdivision algorithm, which is based on the possibility of identifying boxes where polynomial functions are nonnegative. Section 3 discusses how the computation of the value of a polynomial function at the corners of a box provides a su cient condition for the function to be nonnegative over the box. Section 4 presents illustrative examples meant to clarify the proposed algorithm. Finally, Section 5 concludes the paper and it outlines open questions and future research directions. 2
Reasoning on polynomial constraints over boxes
Subdivision algorithms have been used extensively to reason on polynomial
constraints mostly because e ective bounds on the range of polynomial functions
can be computed using the properties of Bernstein polynomials (see, [
Besides applications in mixed-integer nonlinear optimisation [
The subdivision algorithm that we have already adopted to reason on polynomial constraints over nite domains does not mandate the use of a speci c algorithm to study the sign of polynomial functions over boxes, and it can be generalised to a framework to reason on polynomial constraints over nite domains. In order to obtain working algorithms, the framework needs to be completed with a means to decide if a polynomial function is nonnegative over a box. This can be accomplished by means of any algorithm capable of computing suitable lower and upper bounds for a given polynomial function over a given box. If the lower bound is nonnegative, then the polynomial function is nonnegative over the box. Conversely, if the upper bound is negative, then the polynomial function is negative over the box. Note that the trivial approach of evaluating the polynomial function at all points of the box, which is always nite, is surely possible, but we expect that the use of lower and upper bounds would be bene cial especially for large boxes.
Algorithm 1 shows a pseudocode of the proposed framework to enumerate the solutions of constraint satisfaction problems that involve only polynomial constraints over nite domains. Function enumerate(.,.) takes a set of constraints C in canonical form and a box B Zn, with B 6= ;, over which the signs of polynomials in C are studied. First, enumerate(.,.) selects a constraint and it uses the unspeci ed function bounds(.,.) to compute appropriate lower and upper bounds of the related polynomial function over B. If computed upper bound is negative, then the constraint is unsatis able in the box. Conversely, if computed lower bound is nonnegative, then the selected constraint is satis able for all elements of the box. If both conditions are not met, then a variable is selected and the current approximation of its domain is properly split into two disjoint intervals, which are processed recursively. The notation Bj![a::b] is used to denote a box obtained from box B by replacing interval Ij , which corresponds to the current approximation of the domain of the j th variable, with interval [a::b]. Note that the pseudocode in Algorithm 1 considers also the special case of variables whose domains are reduced to singleton sets. In this case, variables can be substituted in polynomials with appropriate values, so that the number of variables in polynomials is reduced. Finally, note that function enumerate(.,.) is rst invoked with the initial set of constraints and the initial approximation of the domains of variables, which is expressed as a nonempty bounding box. 3
Study of the sign of polynomials using function values
In this section we present a method to study the sign of a polynomial function
over a given box that is based on a result originally proposed for univariate
polynomials in [
Univariate polynomials Let p : R 7! R be a polynomial function, and B R be a closed interval. We are interested in computing a lower bound m 2 R and an upper bound m 2 R of p over B such that m min p(x) x2B max p(x) x2B m: In order to compute m and m, we recall the Taylor's theorem for functions in one variable, which is rephrased here to x notation and because there is little consensus on what should be called Taylor's theorem. The theorem is classic and it is shown without proof.
Theorem 1 (Taylor's theorem in one variable). Let f : R 7! R be of class C(k+1), with k 2 N, on an open interval S R. If a 2 S, h 2 R and a + h 2 S, then i=0 k f (a + h) = X D(i)f (a) hi + Ra;k(h); i! where D(i)f is the i th derivative of f (with D(0)f = f ), and the remainder Ra;k(h) is given in Lagrange form by
The application of the Taylor's theorem to a polynomial function on a closed interval easily proves the following corollary. (3) (4) (5) Corollary 1. Let p : R 7! R be a polynomial function and B = [b; b] closed interval. Then, for all x 2 B, and for all k 2 N where the remainder Rb;k(h) is given in Lagrange form by p(x) = Xk D(i)p(b)
(x i=0 i! b)i + Rb;k(x
b);
Proof. Recalling that any polynomial function is of class C1 on any open interval S R with B S, the proof of the corollary follows immediately from the Taylor's theorem. tu In particular, if we truncate to k = 2, we have that p(x) = p(b) + D(1)p(b)(x b) +
D(2)p(b + c(x
2
b))
(x
b)2
for some c 2 (0; 1). The truncated expansion (8) helps to prove the following
proposition (see also [
Proposition 1. Let p : R 7! R be a polynomial function of degree l 2 N with coe cients faigli=0
R, and let U = [0; 1]
R with Then where Proof. Let 2 U be such that p( ) = m, then if is immediately proved because 0 and therefore = 0 or = 1 the proposition (6) (7) (8) (9) (10) (11) (12) with i 2 N, and x 2 U , the expansion of the second derivative in (8) is su cient to prove the proposition. In detail, if 12 , then Conversely, if 6= 0 and 6= 1, then D(1)p( ) = 0. In this case, using the fact that the choice of interval U ensures that for some c 2 (0; 1). But,
l X i(i i=2 D(2)p(x) = 1)aix(i 2)
l X i(i i=2 1)jaij; for all x 2 U , which proves the proposition for 12 . For the case > 12 , the proposition is proved by choosing the other endpoint of interval U to expand the function, which leads to an equality similar to (15). The case of m is treated with similar considerations.
Proposition 1 ensures that, for the case of univariate polynomials, suitable bounds to support the pseudocode in Algorithm 1 can be obtained by computing the values of polynomial functions at the endpoints of considered intervals. Note that the assumption x 2 U in Proposition 1 is not restrictive because, in order to consider a polynomial function p : R 7! R in variable x which takes values from interval B = [b; b] R, it is su cient to change the variable of p to t 2 U such that and to consider the coe cient of the new polynomial in variable t. In this case, we prefer to use the notation B in (12), instead of simply , because the coe cients of the new polynomial make B depend on B. In order to elegantly generalise Proposition 1 to the case of polynomials with n 2 N variables, we need to introduce the notation of multi-indices. A multiindex I 2 Nn is an n-tuple of the form and we adopt the common notation of using an uppercase letter to denote a multi-index and the same letter, but lowercase and with an index, to denote its components. Given a multi-index I 2 Nn, the following abbreviations
I = (i1; i2; : : : ; in); jIj = n X ij
I! = n Y ij ! j=1 j1 j2 X ( ) = X X I J i1=0 i2=0
k k X ( ) = X X jIj k jn X ( ) in=0 k
X ( ) : i1=0 i2=0 in=0 |i1+i2+{z+in k }
n vI = Y vjij :
j=0 Given a multi-index I 2 Nn and v 2 Rn, the following abbreviation is often used are often used to denote the order of I, and the factorial of I, respectively. Given two multi-indices I 2 Nn and J 2 Nn, they can be added
I + J = (i1 + j1; i2 + j2; : : : ; in + jn); and they can be compared
I
J () i1 which justify the following abbreviations, with k 2 N, Finally, given a generic function f : Rn 7! R in variables x = (x1; x2; : : : ; xn), the following notation based on multi-indices is used to denote elegantly the partial derivatives of f
The notation of multi-indices can be used to state the following theorem, which is then used to study the sign of a polynomial function over a given box by computing the value of the function at the corners of the box. The theorem generalises the Taylor's theorem recalled above, and it is shown without proof. Theorem 2 (Taylor's theorem in several variabRlne.s)I.f Lae2t fS,: hRn2 7 !RnRanbde of class C(k+1), with k 2 N, on an open convex S a + h 2 S, then f (a + h) = jIj k X @I f (a) hI + Ra;k(h);
I! where the remainder Ra;k(h) is given in Lagrange form by for some c 2 (0; 1).
Let p : Rn 7! R be a multivariate polynomial function in n 2 N variables x = (x1; x2; : : : ; xn), the adoption of multi-index notation allows expressing p as where L 2 Nn is the multi-degree of the polynomial function, and faI gI L R is the set of its coe cients.
Given that we are interested in studying the sign of polynomial functions in boxes, we introduce the following notation to denote boxes. The box B built using v = (v1; v2; : : : ; vn) 2 Rn and v = (v1; v2; : : : ; vn) 2 Rn is denoted as B = [v; v] = [v1; v1]
Corollary 2. Let p : Rn 7! R be a polynomial function in n 2 N variables, and B = [b; b] Rn be a box. Then, for all x 2 B, and for all k 2 N jIj k (x b)I + Rb;k(x b) and the remainder Rb;k(h) is given in Lagrange form by
Proof. Recalling that any polynomial function is of class C1 for any open convex S Rn with B S, the proof of the corollary follows immediately from the Taylor's theorem. tu
Following the same path traced for the case of univariate polynomials, we can prove the following proposition, which is su cient to compute appropriate lower and upper bounds to support the proposed algorithm. The proposition is shown without proof for the sake of brevity.
Proposition 2. Let p : Rn 7! R be a polynomial function with multi-degree L 2 N in n 2 N variables x = (x1; x2; : : : ; xn) and U = [0; 1]n
Rn, with (28) (29) (30) (31) (32) min p(x) x2C m m max p(x) + ; x2C where C = f0; 1gn is the set of the corners of U and (37) (39)
Note that the assumption x 2 U in Proposition 2 is not restrictive. In order to consider a polynomial function p : Rn 7! R in variable x which takes values from box B = [b; b] Rn, it is su cient to change the variable of p to t 2 U such that and to consider the coe cient of the new polynomial in variable t. In this case, we prefer to use the notation B in (35), instead of simply , because the coe cients of the new polynomial make B depend on B. 4
Illustrative examples This section shows how the proposed algorithm can be applied to some illustrative examples. Note that a systematic validation of the proposed algorithm with speci c benchmarks is still in progress, and the following examples are only intended to clarify the algorithm. The rst example is a quadratic constraint with only one variable p1(x) = x2 50 0 x 2 [1::25]: The constraint is unsatis able in the given interval, and the proposed algorithm proves it by subdividing the initial interval into 4 intervals, the largest of which contains 12 values. Figure 1 shows the subdivision produced by the algorithm.
The second example is an unsatis able constraint with two variables p2(x; y) = x2 y2 100 0 x 2 [1::25]; y 2 [1::25]: (38) The proposed algorithm needs to subdivide the initial box into 12 boxes, the largest of which contains 144 values, to prove that the constraint is unsatis able. Figure 2 shows the subdivision of the box that the algorithm produces.
The following examples are relative to larger domains and the detailed subdivision of the initial approximations of domains are too complex to be shown. The following example is a satis able constraint with one variable p3(x) = x2 16 0 The proposed algorithm builds a search tree made of 65 nodes to compute the following answer x 2 [ 100:: 7 8 13 14 25
Given that each node requires the computation of the value of p3 at the two endpoints of the related interval, 130 evaluations of p3 are needed. Note that a brute force investigation requires 201 evaluations of the polynomial.
The fourth example is still with only one variable p4(x) = x2
x 2 [ 100::0] [ [50::100] is returned after computing a search tree made of 61 nodes, which requires 122 evaluations of p4, against the 201 evaluations needed by brute force investigation.
The fth illustrative example involves two variables p5(x; y) = xy 210 0
x 2 [ 100::100]; y 2 [ 100::100]: The proposed algorithm needs to traverse a search tree made of 1; 909 nodes to compute the answer Note that the computation of the answer required 7; 636 evaluations of p5, four for each considered box, while the brute force scan of the domains of variables requires 40; 401 evaluations of the polynomial.
Finally, the last example involves two variables (41) (42) (43) (44) (45) (46) p6(x; y) = x + xy 1000 0
x 2 [ 100::100]; y 2 [ 100::100]: In this case the answer is computed after traversing a search tree made of 1; 347 nodes, which required 5; 388 evaluations of p6, against the 40; 401 evaluations of brute force search.
Note that discussed examples are processed using the current implementation of the algorithm, which intentionally does not use heuristics to select which variable to split, and where to split the domain of the selected variable. It is expected that the introduction of heuristics to drive the search would bring notable performance improvements. 20 19 14 13 8 7 42 42 42 21 28 72 36 42 144 72 42 42 1 This paper presented an algorithm to enumerate the solutions of constraint satisfaction problems that involve only polynomial constraints over nite domains. Such constraints are expressed as equalities, inequalities, and disequalities of polynomials with integer coe cients whose variables take values from nite subsets of the integers. The proposed algorithm works under the assumption that an initial approximation of the domains of variables is available in terms of a bounding box. First, the algorithm rewrites all constraints into a canonical form in which only inequalities are used. Then, it recursively subdivides the initial bounding box into disjoint boxes, and it records boxes whose elements satisfy all constraints. The canonical form of constraints ensures that a box contains only elements that satisfy all constraints if and only if the polynomial functions related to constraints in canonical form are all nonnegative over the box. In order to verify if a polynomial function is nonnegative over a given box, the algorithm uses lower and upper bounds described in Proposition 2. Such bounds are computed by considering the coe cients of the polynomial function and its values at the corners of the considered box.
The algorithm was tested on illustrative examples and preliminary results
are encouraging, but a systematic validation with speci c benchmarks is still in
progress. We expect that the proposed algorithm would outperform our previous
proposals based on the properties of Bernstein polynomials because of the
inherent complexity of computing Bernstein coe cients, even using sophisticated
techniques (see, e.g., [
The proposed algorithm is implemented as a reusable Java library that can be
accessed as a Prolog module from SWI-Prolog [