Cylindrical Algebraic Decomposition (CAD) is an important tool within computational real algebraic geometry, capable of solving many problems for polynomial systems over the reals. It has long been studied by the Symbolic Computation community and has found recent interest in the Satis ability Checking community. The present report describes a proof of concept implementation of an Incremental CAD algorithm in Maple, where CADs are built and then re ned as additional polynomial constraints are added. The aim is to make CAD suitable for use as a theory solver for SMT tools who search for solutions by continually reformulating logical formula and querying whether a logical solution is admissible. We describe experiments for the proof of concept, which clearly display the computational advantages compared to iterated re-computation. In addition, the project implemented this work under the recently veri ed Lazard projection scheme (with corresponding Lazard valuation).
We aim to adapt Cylindrical Algebraic Decomposition (CAD) for use with
SMTsolvers [
We work over n-dimension real space Rn in which there is a variable ordering. De nition 1 A decomposition of the space X disjoint regions, called cells, whose union is X .
Rn is a nite collection of De nition 2 A set is semi-algebraic if it can be constructed by nitely many applications of union, intersection and complementation operations on sets of the form fx 2 Rn j f (x) 0g where f 2 R[x1; ; xn].
De nition 3 A decomposition D is algebraic if each of its components x 2 D is a semi-algebraic set.
De nition 4 A nite partition of D of Rn is called a cylindrical decomposition of Rn if the projections of any two cells onto any lower dimensional coordinate space with respect to the variable ordering are either equal or disjoint. Thus a cylindrical algebraic decomposition (CAD) satis es De nitions 1 4.
De nition 5 A CAD is sign-invariant with respect to a set of input polynomials
if each polynomial has a constant sign (positive, negative or zero) on each cell.
CADs may be produced with other invariant properties (see for example [
CADs are traditionally produced through a two-stage process: rst projection
identi es polynomials of importance for the invariance property and then lifting
incrementally builds CADs of Rk for k = 1; : : : ; n according to these polynomials.
Decompositions are performed by working at a sample point of a cell, reducing
multivariate polynomials to univariate and then decomposing according to the
output of real root isolation. For a fuller introduction see the lecture notes [
We give a visual example3. The gingerbread face in Figure 1 is formed by four
closed curves, each of which de ned by a bi-variate polynomial equation. A
corresponding sign-invariant CAD of R2 is visualised in Figure 2. We label the
37 open cells (those of two dimensions). There are a further 28 partially open
(1-dimensional line segments) and 28 closed cells (isolated points) giving 93
CAD cells in total. Of course, in many industrial and SMT applications the
polynomials will not form such aesthetically pleasing geometric shapes.
3 Inspired by http://planning.cs.uiuc.edu/node296.html
We aim to work with CADs that change incrementally by constraint. I.e. given a
CAD of Rn sign-invariant for polynomials ff1; : : : ; fmg we aim to adapt this to
one sign-invariant for ff1; : : : ; fm; fm+1g or ff1; : : : ; fm 1g. The present report
deals only with the rst problem. Such incrementality is needed to use CAD in
SMT, and could also bene t CAD directly as a way of reducing the search space.
We proceed by considering the changes required in Projection (Section 2) and
Lifting (Section 3) alongside the issues in using the Lazard projection operator.
We nish with a summary and plans for future work in Section 4. A larger report
with additional details we do not have room for here is available on arXiv [
The present project built upon code from the ProjectionCAD package [
We describe a worked example to illustrate what projection does and to use later to illustrate incremental projection. We work with a system of two polynomials F1 and seek a sign-invariant CAD: Since the problem involves only two variables, we need only a single projection (which we do with respect to x2). The leading coe cients are constant, but the trailing coe cients clearly identify the points on the x1 line at 0 and 1. The discriminants do not identify anything more but the resultant
Resultant(f1; f2; x2) = (x13 + x21 1)2 identi es 1 = 0:75494. Thus the real line is decomposed into 11 cells according to these 5 points. We will compute all the projection polynomials discussed above and some additional ones. The discriminant and leading coe cient of f3 are constant, and the trailing coe cient identi es x1 = 0 which was already present from the trailing coe cient of f2. Similarly, the resultant of f2 and f3 is x12(x1 1) identi es two more roots we saw already. However Resultant(f1; f3; x2) = 2x21 1 identi es two new points, 2 = 1=p2 0:7071. 4 We give a decimal approximation but emphasise that CAD would use the full algebraic number representation: that 1 is the sole real root of (2) in (0; 1).
We developed Algorithms 1 and 2 to implement such an incremental projection.
The pseudo code is closely linked to the Maple implementation. The former
computes all the projection polynomials via calls to the latter which performs one
projection. The adjustments (highlighted in blue) required to make the original
projection code (see [
There were seven elements of the projection set to calculate in the original projection system ProjL(F1) and after the addition of the extra polynomial ProjL(F2) had 12. By avoiding full recomputation, we had to compute only the extra ve. At the cell level the savings are more signi cant: on the real line there were 5 original roots (11 cells), and two new ones (making 15 cells).
We performed experiments to see how such savings transferred into computation time. We created examples using the random polynomial function randpoly in Maple. Testing was conducted through an external bash script creating new Maple instances to avoid any result caching. The testing code is available online5. 5 https://github.com/acr42/InCAD.git
1: Input: Set of calculated projection polynomials prev 2 R[xn; : : : ; x1] and set of new polynomials new 2 R[xn; : : : ; x1] (Adjustment 1) 2: Output: All projection polynomials 2 R[xn; : : : ; x1] 3: Procedure ProjectionPolysAdd (Compute all projection polynomials) 4: dim
Number of variables +1 5: pset[0] = table() 10: for i from 1 to dim-1 do 17: end for 18: pset[dim 19: ret pset[[dim 20: return pset (Adjustment 1) 6: pset[0]
Bivariate polynomials We created 60 pairs of bivariate polynomials, considered rst nding the projection of one and then incrementing the projection by including the other. On average it was 16% faster to increment compared to computing the projection for both polynomials together. However, there was a large variance: the cases which were faster were on average 55% faster, while a small number of cases were slower, by as much as 87%.
Variance 0.0004660s 0.0006425s Mean 0.03743s 0.0315s Lower Quartile 0.024s 0.008s Median 0.0285s 0.015s Upper Quartile 0.03675s 0.05525s
1: Input: Sets of polynomials new = ff1; : : : ; fmg 2 R[xn; : : : ; x1], old 2
R[xn 1; : : : ; x1], and variable ordering (Adjustment 2) 2: Output: Set of polynomials P set = fp1; : : : ; pqg 2 R[xn 1; : : : ; x1] 3: Procedure ProjectionAdd 4: P olys Primitive set from new, wrt variable xn (Adjustment 2) 5: Cont Content set from new, wrt variable xn (Adjustment 2) 6: P olys Square free basis set from P olys, wrt variable xn 7: P set1 = table() : 8: for i from 1 to number of elements of P olys do 9: P ol P olys[i] 10: clist Lazard coe cient set from P ol, wrt to xn 11: temp Discriminant set from P ol, wrt xn 12: temp Remove constant multiples from temp 13: P set1[i] union temp & clist 14: end for 15: P set2 = table() : 16: for i from 1 to number of P olys do 17: for j from i+1 to number of P olys do 18: P set2[i; j] Resultant of P olys[i] and P olys[j], wrt to variable var 19: P set2[i; j] Remove constant multiples from P set2[i; j] 20: end for 21: end for 22: oldset old 23: P set3 = table() : 24: for i from 1 to number of P olys do 25: for j from 1 to size(oldset) do 26: P set3[i; j] Resultant of P olys[i] and oldset[j], wrt to variable var 27: P set3[i; j] Remove constant multiples from P set3[i; j] 28: end for 29: end for(Adjustment 3) 30: P set union (cont; P set1; P set2; P set3) 31: P set Remove constant multiples from P set 32: return P set Trivariate polynomials We next created 80 further examples with pairs of trivariate polynomials, in this case restricting to 4 terms per polynomial. Here there were greater savings, on average 29% faster to increment than to compute altogether. There was also a smaller variance in the timings of the example set, although it was still the case that a few examples were slower to increment than recompute.
We suggest that the extra overheads of the incremental approach will become less important in comparison to the savings as the number of variables increase: indeed, this follows from the well-known complexity results on CAD.
We had to make changes to the lifting code in ProjectionCAD, not just to
allow for incrementality but also to validate the use of the Lazard projection
operator [
1 do 6: 7: 8: 9: while f (r1; : : : ; rj ) = 0 do f f =(xj
rj ) substitute xj = rj into f .
end while 10: end for 11: return f: We will now discuss our approach for incremental lifting, which can be thought of graphically as a form of acyclic tree merge, as shown later. 3.2
We will describe lifting the projection polynomial system de ned previously
ProjL(F1). Recall that we identi ed four points on the real line: f 1; 0; 1; 1g.
Thus, we need to choose a sample value from the 9 cells in the decomposition:
a1 = fx1 < 1g; a2 = fx1 = 1g;
a4 = fx1 = 0g; a5 = f0 < x1 <
a7 = f 1 < x1 < 1g; a8 = fx1 = 1g;
a3 = f 1 < x1 < 0g;
1g; a6 = fx1 = 1g;
a9 = f1 < x1g
(4)
We are forced to pick non-rational sample points for cell a6 but the others can
be rational. We choose sample points: ( 2; 1; 21 ; 0; 12 ; 1; 190 ; 1; 2):
We lift over each cell at the designated sample point by isolating real roots of
the univariate polynomials in x2 we get by from the Lazard valuation at the
sample point in x1. Below, pi;j denotes the polynomial acquired after applying
the Lazard valuation method to the i'th sample point, on the j'th polynomial
from F1. For example, if we use sample point 21 for cell a3 then polynomial f1
becomes 41 + x22 1 which has two real roots at 0 = p3=2 0:8660. We
proceed this way to generate our CAD cells in R2. The structure is as shown in
Figure 5. For a full list of the new cell descriptions see [
We skip the projection steps7. The addition polynomials identify one further point on the real line: 1. The decomposition of the real line is now: a1 = fx1 < 1g; a2 = fx1 = 1g; a3 = f 1 < x1 < 1g; a4 = fx1 = 1g; a5 = f 1 < x1 < 0g; a6 = fx1 = 0g; a7 = f0 < x1 < 1g; a8 = fx1 = 1g; a9 = f 1 < x1 < 1g; a10 = fx1 = 1g; a11 = f1 < x1g and our sample points are: 2; 1; 190 ; 1; 12 ; 0; 12 ; 1; 190 ; 1; 2.
We rst test whether each sample point leads to new roots when lifting with f4.
If so we must re ne the decomposition in the cylinder above. For example, on a2
we have sample point x1 = 1, and the Lazard valuation of f4 is x22 1 with real
roots at 1. However, we had already identi ed these from other polynomials,
sxo22 no18chaanndgweeis rnedqutiwreodn. eHwowreeavlerro,ootns.aW5ewitthhussarmepnlee pthoeindteco12m, pfo4seitviaolnuaatbeosvteo.
We then lift over the new sample points with respect to all polynomials creating
new decompositions. For example, at x1 = 190 we have two roots from the
valuation of f1 and another two from f4 (so decomposition above into 9 cells). Figures
10-12 show the new CAD tree structure and its split into new and unchanged
cells, illustrating potential savings. Full cell descriptions are in [
When incrementing the lift stage, we can think of it as starting at the root node of the old CAD tree and working our way down it, one depth level at a time, until we reach the leaves. In the process of working down the old tree, we will be creating subtrees, which will later be reconnected to the unchanged tree, to form the new incremented CAD tree.
One tree (UnchangedCells) is a strict subset of the old trees nodes and edges, discovered through analysing the old structure down and Lazard valuating on each cell not marked as new, at each depth p with new projection polynomials in Rp+1. Then, if there are new real roots discovered, we prune inspected cells children, and the cell is sent to the NewCells set for full re-computation. In the NewCells set, we will then use this cell to form a subtree, to later be reconnected with the UnchangedCells tree. Each cell in the NewCells list is a subtree, to later be reconnected via source indices.
When going through the old CAD tree structure, we have only two cases: CASE 1: When a node has new children: New real roots have been acquired from one of the new projection polynomials. We start by pruning all of the child branches in the old tree structure, by labelling them as new, then performing a full lift onto the set of all projection polynomials from Rk; : : : ; Rn, where k is the depth the new root was discovered. When we label a cell as new, e ectively that halts tree growth in the UnchangedCells structure, so that later on we can attach the branch extensions, gained from the incremental lift. Such cells have an updated source index, as its source cell would now be saved in a new tree structure with a new index.
CASE 2: When a node has no new children.
The new ag is passed down to children from parents. Cells which are not new are stored in UnchangedCells. Cells here only have Lazard valuation at the new projection polynomial (rather than all projection polynomials). If this leads to a new root, we move it over to the NewCells structure. Otherwise, we continue with its child cells.
When moving cells into UnchangedCells, we make sure that the indexing of each cell does not clash with that of the indexing in NewCells at each depth level in the tree. We then merge the NewCells and UnchangedCells trees, forming the full incremented CAD tree. At each stage of the lift, we merge-sort the list.
We conducted testing for the incremental lift method on the same examples used to test the incremental projection earlier.
Bivariate polynomials On average 30% faster than recomputing. Trivariate polynomials On average only 7% faster; one example 28% slower.
Variance 0.003734s 0.002903s Mean 0.1778s 0.1240s Lower Quartile 0.1328s 0.089s Median 0.163s 0.1255s Upper Quartile 0.226s 0.163s 28.63% Smaller 30.25% Faster 32.96% Faster 23.01% Faster 27.87% Faster
Variance 0.05541s 0.06838s
Mean 0.2880s 0.2687s
Lower Quartile 0.1275s 0.0995s
Median 0.207s 0.164s
Upper Quartile 0.3605s 0.3523s
18.96% Larger
6.707% Faster
21.96% Faster
20.77% Faster
2.29% Faster
So the lifting code shows opposite results to projection with the savings
decreasing with input size: indicating that the overheads required for the incremental
work grow faster than the savings (at least for the size of examples studied).
Of course, we can also experiment with combined projection and lifting.
Unsurprisingly the lifting costs dominate. For the bivariate polynomials incremental
code was 37% faster but for the trivariate only 12% faster (see [
The presented work acts as a proof of concept that incremental CAD construction in Maple is possible with savings on o er. Care needs to be given to the datatypes used. Beyond that the main areas of further work are moving out of the open case (our implementation restricted lifting to open cells), considering what happens if the new polynomial has a variable not already represented in the system, and considering incremental reduction of constraints. To remove a polynomial from the CAD means nding all those projection polynomials created from only that source polynomial (or as a resultant of that with another) and removing them and the cylinder splits their real roots caused. Acknowledgements This work was funded by the EU's H2020 programme under grant No H2020-FETOPEN-2015-CSA 712689 (SC2 ). We thank C. Brown for a tutorial on the Lazard valuation; S. Timms, J.H. Davenport and S. Forrest for useful discussions; and the organisers of the SC2 2017 Summer School where these took place.
Algorithm 4 LiftSetupAdd
1: Input: A sets of new projection polynomials new, an oldcad, a table of sets
of all projection polynomials psetf ull, and a variable ordering
2: Output: [N ewCells; OldCad; OldRoots; U nchangedCells] where NewCells
are in the last variable Lif tIncf2 is information for the next lift, OldCad
contains the previous CAD tree, and UnchangedCells is a subset.
3: Procedure LiftSetup
4: cad table()
5: N ewRoots [ ]
6: N ewCells table()
7: U nchangedCells table()
8: Lif tIncf2 [ ]
9: for i from 1 to size(new[
U nchangedCells[
U nchanged[d] and update index accordingly. 22: end for 23: end for 24: F inalU nchangedCells [] 25: F inalU nchangedCells Union U nchanged[d][f ] for all cells f with their corresponding ags not equal to new. 26: IncCAD F inalU nchangedCells Union N ewCells[d] 27: return IncCAD