An implementation of Sub-CAD in Maple
aa r X i v : . [ c s . S C ] M a r An implementation of Sub-CAD in Maple
Matthew England and David WilsonDepartment of Computer Science, University of Bath, Bath, UK
[email protected] [email protected]
Abstract
Cylindrical algebraic decomposition (CAD) is an important tool for the investigationof semi-algebraic sets, with applications in algebraic geometry and beyond. We havepreviously reported on an implementation of CAD in
Maple which offers the originalprojection and lifting algorithm of Collins along with subsequent improvements.Here we report on new functionality: specifically the ability to build cylindrical al-gebraic sub-decompositions (sub-CADs) where only certain cells are returned. We haveimplemented algorithms to return cells of a prescribed dimensions or higher (layered sub-CADs), and an algorithm to return only those cells on which given polynomials are zero(variety sub-CADs). These offer substantial savings in output size and computation time.The code described and an introductory
Maple worksheet / pdf demonstrating thefull functionality of the package should accompany this report.
This work is supported by EPSRC grant EP/J003247/1.
This report concerns
ProjectionCAD : a
Maple package for cylindrical algebraic decomposi-tion (CAD) developed at the University of Bath. The extended abstract [18] at ICMS 2014describes how this package utilises recent CAD work in the
RegularChains
Library of
Maple ,while still following the classical projection and lifting framework for CAD construction. Thepresent report is to accompany the release of
ProjectionCAD version 3, describing the newfunctionality this introduced. The report should be accompanied by the code described andan introductory
Maple worksheet / pdf demonstrating the full functionality of the package.The previous two versions of
ProjectionCAD are hosted alongside similar reports documentingtheir functionality [16, 17].Version 3 introduces functionality for cylindrical algebraic sub-decompositions (sub-CADs): subsets of CADs sufficient to describe the solutions of a given formulae. Two distincttypes are provided, whose theory was developed in [27]. The first type contains only thosecells of a certain dimension and higher, reducing both the output size and computationaltime by giving only output of the required generality. We have implemented both a directand recursive algorithm to build these layered sub-CADs . The second type contains only thosecells on which given equations are satisfied (lie on a prescribed variety). When building aCAD for a formula with an equational constraint then only these cells can contain the solutionset. These variety sub-CADs clearly reduce the output, and can also reduce computation timedepending on the rank of the variety relative to the variable ordering.1. England and D. Wilson An implementation of Sub-CAD in MapleWe continue the introduction by reminding the reader of the necessary background mate-rial on CAD and summarising our previous work with
ProjectionCAD . The following sectionsthen describe the new functionality in
ProjectionCAD version 3. A cylindrical algebraic decomposition (CAD) is: a decomposition of R n meaning acollection of non-intersecting cells whose union is R n ; algebraic meaning each cell can bedescribed by a finite sequence of polynomial relations; cylindrical meaning that with respectto a given variable ordering, the projection of any two cells onto a lower dimensional space(with respect to the ordering) is either equal or disjoint. CAD was introduced by Collins [11]who provided an algorithm which given a set of polynomials would produce a CAD which was sign-invariant : so each polynomial had constant sign in each cell. Collins developed CADas a tool for quantifier elimination in real closed fields but CAD has since found applicationsas diverse as robot motion planning [28, etc] and algebraic simplification [13, etc.].Collins’ algorithm (see for example [1]) has two phases. The first, projection , applies a projection operator repeatedly to a set of polynomials, each time producing another setin one fewer variables. Together these contain the projection polynomials . The secondphase, lifting , builds the CAD incrementally. First R is decomposed into cells which arepoints and intervals corresponding to the real roots of the univariate polynomials. Then R is decomposed by repeating the process over each cell using the bivariate polynomials at asample point of the cell. The output for each cell consists of sections of polynomials (wherea polynomial vanishes) and sectors (the regions between these). Together these form the stack over the cell, and taking the union of these stacks gives the CAD of R . This processis repeated until a CAD of R n is produced. Collins’ original projection operator was definedso that the CAD of R n produced using sample points in this way could be concluded sign-invariance. The key tool to conclude this was polynomials being delineable in a cell, meaningthe zero set of individual polynomial are disjoint sections and the zero set of the sections fromdifferent polynomials are identical or disjoint. More efficient projection operators have sincebeen developed to conclude both sign-invariance and other invariance conditions.The output of a CAD algorithm depends on the ordering of the variables. In this paper weusually work with ordered variables x = x ≺ . . . ≺ x n , (so we first project with respect to x n and continue until we have univariate polynomials in x ). The main variable of a polynomial,mvar( f ), is the greatest variable present with respect to the ordering.All cells are equipped with a sample point and a cell index . The index is an n -tuple ofpositive integers that corresponds to the location of the cell relative to the rest of the CAD.Cells are numbered in each stack during the lifting stage (from negative to positive), withsectors having odd numbers and sections having even numbers. Therefore the dimension of agiven cell can be easily determined from its index: simply the number of odd indices in the n -tuple. Note that the sub-CADs we discuss later are index-consistent , meaning a cell in asub-CAD will have the same index as if it was produced in the associated complete CAD.Since Collins published the original algorithm there has been much research into improve-ments with a summary of the developments over the first twenty years given by [12]. Keydevelopments since then include the use of certified numerics [26, 19] in the lifting phase, pro-jection operators for studying multiple formulae [3] and building CADs by first decomposingcomplex space and refining to real space, instead of projecting and lifting [10, 9, 2, 15].2. England and D. Wilson An implementation of Sub-CAD in Maple ProjectionCAD package
ProjectionCAD is a
Maple package developed at The University of Bath to implement CADvia projection and lifting. Its name was chosen to distinguish it from the CAD algorithmin [10] which is distributed with
Maple as part of the
RegularChains library [24, 20, etc.].Nevertheless, the package makes much use of procedures from the
RegularChains
Libraryand the motivation and details for this are given in [18].In [16] we described our implementation of both McCallum’s and Collin’s algorithms toproduce sign-invariant CADs. In particular, we highlighted that this was the only imple-mentation to offer order-invariant CADs (meaning each polynomial has constant order ofvanishing on each cell) and delineating polynomials [22, 7] which modify the lifting phaseto allow McCallum’s projection operator to be applied more widely. The latter meant thatsome examples of unnecessary failure in
Qepcad [6] could be avoided, while the former wasnecessary for the extensions in the second release. These were described in [17] and allowedfor CADs invariant with respect to an equational constraint or CADs which were truth-tableinvariant (TTICADs) with respect to a list of formulae (meaning each formulae has constantBoolean truth value on each cell). The TTICAD theory was developed in [3, 4] and is basedon an reduced projection operator analogous to McCallum’s reduced operator for equationalconstraints [23]. In addition it was noted that the reduced projection theory allows also forimprovement to the lifting phase, leading to
ProjectionCAD offering lower cell counts than
Qepcad even for examples with only one equational constraint. The second release of thecode also included heuristics to help with choices of problem formulation following [14, 5].
Define cells in a CAD with the same dimension as a layer and let ℓ be an integer, 1 ≤ ℓ ≤ n +1.Then an ℓ -layered sub-CAD is the subset of a CAD of R n with cells of dimension n − i for 0 ≤ i < ℓ . They were introduced by the present authors in [27] (although there had beenprevious work where only cells of full-dimension were computed [21, 25]).Algorithm 1 describes how an ℓ -layered sub-CAD may be produced with the main ideato run the projection phase as normal (step 1), but then truncate in the lifting phase whencells of a high enough dimension cannot be produced. The algorithm can run with differentprojection operators: that of Collins [1] or McCallum [22] to build sign-invariant CADs forpolynomials (as detailed in [16]); or operators to build truth-invariant CADs for formulaeusing equational constraints [23] and truth-table invariance [3, 4] (as detailed in [17]).We build the CAD of the real line as normal (step 2) but then for each successive lift wefirst check the dimension of the cell to be lifted over (step 6) and only proceed to generate thestack if cells of dimension ℓ or higher in R n could be produced from it. The stack generationcommand is detailed in [18] and makes use of the RegularChains stack generation coderequiring a preprocessing step. Some projection operators can incur theoretical failure forinput that is not well-oriented, in which case this is identified during stack generation andFAIL returned. A final loop is used (steps 14 −
17) to remove any lower dimensional cellsthat were produced as part of a stack. The approach was verified in [27].As previously described in [16] there are various output formats available. At the leasteach cell is represented by an index (positioning it in the CAD) and a sample point (encodedas a regular chain and bounding box) while an algebraic description is also available. Of3. England and D. Wilson An implementation of Sub-CAD in Maple
Algorithm 1:
LCAD ( F, x , ℓ, ProjOp ): Algorithm to produce ℓ -layered sub-CADs. Input : A choice of projection operator
ProjOp ; input F ( x ) (polynomials orformulae) in ordered variables x ; and an integer ℓ ∈ (1 , n + 1). Output : An ℓ -layered sub-CAD for F , or FAIL if F not well-oriented. Run the projection phase and set P [ i ] to be the projection polynomials with mvar x i ; Set D [1] to be the CAD of R obtained by isolating the roots of P [1]; for i = 2 , . . . , n do Initialise D [ i ] := [ ]; for c ∈ D [ i − do Evaluate dim := P α ∈ c. index ( α mod 2); if dim > i − ℓ − then S := GenerateStack ( P [ i ] , c ); if S =FAIL then return FAIL; // Input is not well oriented else Add the cells in S to D [ i ]; Initialise D := [ ]; for c ∈ D [ n ] do Evaluate dim := P α ∈ c. index ( α mod 2); if dim > n − ℓ then Add c to D ; return D ;particular note it the intuitive piecewise construction (described in [8]) which highlights thetree-like nature of CAD. This has been adapted in ProjectionCAD to display layered sub-CADs with the truncated branches clearly indicated.A simple example of using
LCAD in ProjectionCAD now follows, appearing as it wouldin a
Maple worksheet, except that the sample points have been replaced by SP for brevity. > LCAD([x^2+y^2-1], 1, [y,x], method=McCallum, output=piecewise); SP x < − ∗ ∗ ∗ branch = truncated SP y < −√− x + 1 ∗ ∗ ∗ branch = truncated SP −√− x + 1 < y < √− x + 1 ∗ ∗ ∗ branch = truncated SP √− x + 1 < y − < x < ∗ ∗ ∗ branch = truncated SP < x Here a sign-invariant sub-CAD for the defining polynomial of the unit circle has been4. England and D. Wilson An implementation of Sub-CAD in Mapleproduced. The number 1 as the second input indicates that only one layer of cells is to beproduced: those five cells with full dimension. Two simplifications have been in this sub-CAD compared to a complete CAD. First when lifting over the CAD of the real line the twopoints ( x = ±
1) have not been lifted over: they failed the dimension check at step 6 and werethen excluded from the output during steps 14 −
17. Second, those cells in the stack overthe interval ( − ,
1) on the real line withut full dimensional were also discarded in the finalloop. Note that while the second simplification reduced the output, the first also reduced thecomputation time since no real root isolation was required over 2 of the 5 cells in R .We could have instead asked for two layers of cells to receive the output below. > LCAD([x^2+y^2-1], 2, [y,x], method=McCallum, output=piecewise); SP x < − SP y < ∗ ∗ ∗ branch = truncated SP < y x = − SP y < −√− x + 1 SP y = −√− x + 1 SP −√− x + 1 < y < √− x + 1 SP y = + √− x + 1 SP √− x + 1 < y − < x < SP y < ∗ ∗ ∗ branch = truncated SP < y x = 1 SP < x This time the only missing cells are the isolated points x = ± , y = 0. These would beincluded in a 3-layered sub-CAD, which is itself a complete CAD. Hence in this case LCAD would give exactly the same output as
CADFull (described in [16]).Note that each of
LCAD command is distinct: re-running the projection phase and recom-puting the higher dimensional cells.
ProjectionCAD also contains a recursive command tobuild layers incrementally. It takes additional input as a list of cells already computed andsections that were previously discarded (in the final loop). It then lifts over those discardedsections to produce cells of one dimension higher, as well as further discarded sections. Thisrecursive algorithm is described fully and verified in [27].Algorithm 1 allows for projection operators other than McCallum’s. In particular, it canbe used with the TTICAD projection operator of [3] to build truth-table invariant sub-CADs,as for the following problem. > f1 := x^2+y^2+z^2-1: g1:=x*y*z-1/4:> f2 := (x-1)^2+(y-1)^2+(z-1)^2-1: g2:=(x-1)*(y-1)*(z-1)-1/4:> PHI := [ [f1,[g2]], [f2,[g2]] ]:
The syntax here means that PHI defines two formulae: the first with EC defined by f andadditional constraints defined also by g ; and the second with EC defined by f and additional5. England and D. Wilson An implementation of Sub-CAD in Mapleconstraints defined by also g . The TTICAD command could produce a decomposition sign-invariant for f and f and also for each g i where the corresponding f i is zero. > TTICAD(PHI, [z,y,x]); > LTTICAD(PHI, 1, [z,y,x]), LTTICAD(PHI, 2, [z,y,x]),LTTICAD(PHI, 3, [z,y,x]), LTTICAD(PHI, 4, [z,y,x]); , , , An equational constraint (EC) is an equation, f = 0, logically implied by the truth of aTarski formula. They may be explicit, like f = 0 ∧ φ , or implicit as f f = 0 is in( f = 0 ∧ φ ) ∨ ( f = 0 ∧ φ ) . The presence of an EC can reduce the number of projection polynomials required [23, 3, etc.]and
ProjectionCAD already has the related commands
ECCAD and
TTICAD (described in [17]).Given a Tarski formula with EC f = 0 we define a variety sub-CAD as a truth-invariantsub-CAD for the formula consisting only of cells lying in the variety defined by f = 0.Algorithm 2 describes how a variety sub-CAD may be produced. We assume that allfactors of the EC have the same main variable as the overall system ( x n ) and so the EC isused only for the first projection and final lift. Building variety sub-CADs outside of thisrestriction is considered in [27] but not yet implemented in ProjectionCAD .The projection proceeds as normal in step 1, as does lifting to a CAD of R n − in step 2.For the final lift we generate stacks only with respect to the EC in step 7, noting that theprojection theory used should ensure this is sufficient to conclude truth invariance (see [4] fordetails). Further, only the sections in those stacks are retained for the output in step 11 asonly these describe the variety. Since the underlying CAD algorithms that utilise ECs canreturn theoretical failure for input that is not well oriented we assume this is tested for duringstack generation. This approach was verified in [27].Consider the formula ϕ := ( x + y − ∧ ( xy − > ∧ ( x − y > R where this is true we could build a sign-invariant CAD for the three polynomials. CADFull([x^2+y^2-1,x*y-1/4,x^3-y^2], [x,y], method=McCallum): nops(%); ϕ with the ECCAD command: 6. England and D. Wilson An implementation of Sub-CAD in Maple
Algorithm 2:
VCAD ( ϕ, f, x ): Algorithm to produce variety sub-CADs. Input : A choice of CAD projection
ProjOp ; input F ( x ) (polynomials or formulae) inordered variables x ; an EC f = 0 (in which all factors have same mvar as F ). Output : A (truth-invariant) variety sub-CAD for ϕ , or FAIL if F not well-oriented. Set P to be the output from applying ProjOp to F ; Compute a CAD of R n − for P and assign to D ′ ; if D ′ = FAIL then return FAIL; // P is not well oriented Initialise D := []; for c ∈ D ′ do Calculate S := GenerateStack ( f, c ); if S = FAIL then return FAIL; // F is not well oriented if | S | > then Add those cells in S with last entry of their index even to D ; return D ; ECCAD( [x^2+y^2-1, [x*y-1/4,x^3-y^2]], [x,y]): nops(%);
VCAD command goes further and returns onlythose cells where the EC is satisfied (a sub-CAD) more than halving the output again:
VCAD( [x^2+y^2-1, [x*y-1/4,x^3-y^2]], [x,y]): nops(%);
ProjectionCAD version 3
Combining sub-CADs
We can build sub-CADs which are both variety and layered using the
LVCAD command. In fact,we can combine all the
ProjectionCAD theory to build layered variety truth-table invariantsub-CADs with the
LVTTICAD command. More examples and details are given in [27], whilein [28] such a combination was used to solve a long-standing robot motion planning problem.
Cell distribution in CADs
Tools are provided to compare the growth in CADs by layer, that is, the distributions of cellsin CADs by dimensions. Analysis in [29] suggests a common distribution with little variationfor the underlying problem, meaning layered sub-CADs can be used as heuristics to guide theconstruction of complete CADs. 7. England and D. Wilson An implementation of Sub-CAD in Maple
Using sub-CADs to avoid theoretical failure
Algorithms based on McCallum’s theory can only run on input that is well-oriented (adefinition particular to each individual operator, see for example [23, 3]). These conditionsare checked during lifting and if not satisfied
ProjectionCAD gives a warning. For example: > f := a*e+b*d+c*e+d+e: vars:=[a,b,c,d,e]:> CADFull( [f], vars, method=McCallum ):Warning, The input is not well-oriented (nullification on cell [1, 2, 2]).The output cannot be guaranteed correct.
We can build a 1- or 2-layered sub-CAD for the polynomial without triggering the warning: > LCAD([f],1,[a,b,c,d,e], method=McCallum,failure=err):> LCAD([f],2,[a,b,c,d,e], method=McCallum,failure=err): nops(%%), nops(%); , We have described the new functionality present in the third release of
ProjectionCAD ,which focusses on algorithms to build sub-CADs. The underlying theory is discussed furtherin [28, 29, 27] and the present report should be accompanied by the code and an introductory
Maple worksheet demonstrating the full functionality of the package. A key topic for futurework is building variety sub-CADs when the EC does not have all factors in the main variable.Various approaches are laid out in Section 2.1 of [27] and an analysis of these is ongoing.
References [1] D. Arnon, G.E. Collins, and S. McCallum. Cylindrical algebraic decomposition I: The basic algorithm.
SIAM Journal of Computing , 13:865–877, 1984.[2] R. Bradford, C. Chen, J.H. Davenport, M. England, M. Moreno Maza, and D. Wilson. Truth table invari-ant cylindrical algebraic decomposition by regular chains. In
Computer Algebra in Scientific Computing (LNCS 8660), pages 44–58. Springer, 2014.[3] R. Bradford, J.H. Davenport, M. England, S. McCallum, and D. Wilson. Cylindrical algebraic decompo-sitions for boolean combinations. In
Proc. ISSAC ’13 , pages 125–132. ACM, 2013.[4] R. Bradford, J.H. Davenport, M. England, S. McCallum, and D. Wilson. Truth table invariant cylindricalalgebraic decomposition.
Submitted for publication.
Preprint at http://opus.bath.ac.uk/38146/ , 2015.[5] R. Bradford, J.H. Davenport, M. England, and D. Wilson. Optimising problem formulations for cylindricalalgebraic decomposition.
Intelligent Computer Mathematics (LNCS 7961), pages 19–34. Springer, 2013.[6] C.W. Brown. QEPCAD B: A program for computing with semi-algebraic sets using CADs.
ACM SIGSAMBulletin , 37(4):97–108, 2003.[7] C.W. Brown. The McCallum projection, lifting, and order-invariance. Technical report, U.S. NavalAcademy, Computer Science Department, 2005.[8] C. Chen, J.H. Davenport, J. May, M. Moreno Maza, B. Xia, R. Xiao, and Y. Xie. User interface designfor geometrical decomposition algorithms in Maple. In
Proc. Mathematical User-Interface , 12pp, 2009.
8. England and D. Wilson An implementation of Sub-CAD in Maple [9] C. Chen and M. Moreno Maza. An incremental algorithm for computing cylindrical algebraic decompo-sitions.
Proc. ASCM ’12 , 2012. To appear, Springer. Preprint: arXiv:1210.5543 .[10] C. Chen, M. Moreno Maza, B. Xia, and L. Yang. Computing cylindrical algebraic decomposition viatriangular decomposition. In
Proc. ISSAC ’09 , pages 95–102. ACM, 2009.[11] G.E. Collins. Quantifier elimination for real closed fields by cylindrical algebraic decomposition. In
Proc.2nd GI Conference on Automata Theory and Formal Languages , pages 134–183. Springer-Verlag, 1975.[12] G.E. Collins. Quantifier elimination by cylindrical algebraic decomposition – 20 years of progress. InB. Caviness and J. Johnson, editors,
Quantifier Elimination and Cylindrical Algebraic Decomposition ,Texts & Monographs in Symbolic Computation, pages 8–23. Springer-Verlag, 1998.[13] J.H. Davenport, R. Bradford, M. England, and D. Wilson. Program verification in the presence of complexnumbers, functions with branch cuts etc. In
Proc. SYNASC ’12 , pages 83–88. IEEE, 2012.[14] A. Dolzmann, A. Seidl, and T. Sturm. Efficient projection orders for CAD. In
Proc. ISSAC ’04 , pages111–118. ACM, 2004.[15] M. England, R. Bradford, C. Chen, J.H. Davenport, M. Moreno Maza, and D. Wilson. Problem formula-tion for truth-table invariant cylindrical algebraic decomposition by incremental triangular decomposition.In
Intelligent Computer Mathematics (LNAI 8543), pages 45–60. Springer, 2014.[16] M. England. An implementation of CAD in Maple utilising McCallum projection. Uni. Bath Dept.Computer Science Technical Report series 2013-02. Available at http://opus.bath.ac.uk/33180/ , 2013.[17] M. England. An implementation of CAD in Maple utilising problem formulation, equational constraintsand truth-table invariance. Uni. Bath Dept. Computer Science Technical Report series 2013-04. Availableat http://opus.bath.ac.uk/35636/ , 2013.[18] M. England, D. Wilson, R. Bradford, and J.H. Davenport. Using the Regular Chains Library to buildcylindrical algebraic decompositions by projecting and lifting.
Mathematical Software – ICMS 2014 (LNCS8592), pages 458–465. Springer Heidelberg, 2014.[19] H. Iwane, H. Yanami, H. Anai, and K. Yokoyama. An effective implementation of a symbolic-numericcylindrical algebraic decomposition for quantifier elimination. In
Proc. SNC ’09 , pages 55–64, 2009.[20] F. Lemaire, M. Moreno Maza, and Y. Xie. The RegularChains Library in Maple 10. In I.S. Kotsireas,editor,
Proceedings of Maple Summer Conference ’05 , 2005.[21] S. McCallum. Solving polynomial strict inequalities using cylindrical algebraic decomposition.
The Com-puter Journal , 36(5):432–438, 1993.[22] S. McCallum. An improved projection operation for cylindrical algebraic decomposition. In B. Cavi-ness and J. Johnson, editors,
Quantifier Elimination and Cylindrical Algebraic Decomposition , Texts &Monographs in Symbolic Computation, pages 242–268. Springer-Verlag, 1998.[23] S. McCallum. On projection in CAD-based quantifier elimination with equational constraint. In
Proc.ISSAC ’99 , pages 145–149. ACM, 1999.[24] M. Moreno Maza. On triangular decompositions of algebraic varieties. NAG Tech. Report TR 4/99, 1999.[25] A. Strzebo´nski. Solving systems of strict polynomial inequalities.
Journal of Symbolic Computation ,29(3):471–480, 2000.[26] A. Strzebo´nski. Cylindrical algebraic decomposition using validated numerics.
Journal of Symbolic Com-putation , 41(9):1021–1038, 2006.[27] D. Wilson, R. Bradford, J.H. Davenport, and M. England. Cylindrical algebraic sub-decompositions.
Mathematics in Computer Science , 8:263–288, 2014.[28] D. Wilson, J.H. Davenport, M. England, and R. Bradford. A “piano movers” problem reformulated. In
Proc. SYNASC ’13 , pages 53–60. IEEE, 2013.[29] D. Wilson, M. England, J.H. Davenport, and R. Bradford. Using the distribution of cells by dimensionin a cylindrical algebraic decomposition. In
Proc. SYNASC ’14 , pages 53–60. IEEE, 2014., pages 53–60. IEEE, 2014.