Solving Polynomial Systems over Semialgebraic Sets
Represented by Cylindrical Algebraic Formulas

Adam Strzebonski

Wolfram Research Inc.
100 Trade Centre Dr.
Champaign, IL 61820
U.S.A.

adams@wolfram.com
http://members.wolfram.com/adams

Abstract

Cylindrical algebraic formulas are an explicit representation of semialgebraic sets as finite unions of cylindrically arranged disjoint cells bounded by graphs of algebraic functions. We present a version of the Cylindrical Algebraic Decomposition (CAD) algorithm customized for solving systems of polynomial equations and inequalities over semialgebraic sets given in this representation. The algorithm can also be used to solve conjunctions of polynomial conditions in an incremental manner. We show application examples and give an empirical comparison of incremental and direct CAD computation.

CAFView code

cview[cad_,vars_,rul_,dvars_,ind_]:=
Module[{k=Length[ind],cc,ee,a,b,ans=False,ha,hb,x,na,rr},
ha=Hue[k/Length[vars],1/3];
hb=Hue[(k+1)/Length[vars],1/3];
If[rul=!={},
a=cad/.rul;
If[a===True||a===False,Return[Row[{a},Background->ha,Frame->True]]]];
x=dvars[[k+1]];
cc=If[Head[cad]===Or,List@@cad,{cad}];
If[k==Length[vars]-1,
Return[Or@@(Row[{Tooltip[#,N[#/.rul]]},Background->ha,Frame->True]&/@cc)]];
ee=cc[[1]];
If[Head[ee]===And,ee=ee[[1]]];
If[!FreeQ[ee,Alternatives@@Drop[vars,k+1]],
Return[Row[{True},Background->ha,Frame->True] &&Row[{Opener[Dynamic[Evaluate[x]]],PaneSelector[{False->Row[{Subscript["b",Sequence@@ind,1]},Frame->True,Background->hb],True->Row[{cview[cad,vars,rul,dvars,Append[ind,1]]},Frame->True]},Dynamic[Evaluate[x]],ImageSize->Automatic]}]]];
Do[ee=cc[[i]];
{a,b}=If[Head[ee]===And, {ee[[1]],Drop[ee,1]},{ee,True}];
na=N[a/.rul];
rr=If[na===True,rul,{}];
If[na=!=False,ans=ans||Row[{Tooltip[a,na]},Background->ha,Frame->True] &&Row[{Opener[Dynamic[Evaluate[x]]],PaneSelector[{False->Row[{Subscript["b",Sequence@@ind,i]},Frame->True,Background->hb],True->Row[{cview[b,vars,rr,dvars,Append[ind,i]]},Frame->True]},Dynamic[Evaluate[x]],ImageSize->Automatic]}]],
{i,Length[cc]}];
ans]

CAFView[cad_,vars_,rul_:{}]:=With[{dvars=Table[Unique["dynamicx"],{Length[vars]}]},cview[cad,vars,rul,dvars,{}]]

Semialgebraic sets

A system of polynomial equations and inequalities in variables (Subscript[x, 1]),…,(Subscript[x, n]) is a formula

S((Subscript[x, 1]),...,(Subscript[x, n]))== (Underoverscript[∨, (i=1), k]) (Underoverscript[∧, (j=1), l]) (Subscript[f, i,j])((Subscript[x, 1]),...,(Subscript[x, n])) (Subscript[ρ, i,j]) 0

where (Subscript[f, i,j])∈ℝ[(Subscript[x, 1]),…,(Subscript[x, n])], and each (Subscript[ρ, i,j]) is one of <, ≤, >, ≥, =, or ≠.

A subset of  (ℝ^n) is semialgebraic if it is a solution set of a system of polynomial equations and inequalities.

A quantified system of real polynomial equations and inequalities in free variables (Subscript[x, 1]),…,(Subscript[x, n]) and quantified variables (Subscript[t, 1]),…,(Subscript[t, m]) is a formula

(Subscript[Q, 1]) (Subscript[t, 1])... (Subscript[Q, m]) (Subscript[t, m]) S((Subscript[t, 1]),...,(Subscript[t, m]),(Subscript[x, 1]),...,(Subscript[x, n]))

where (Subscript[Q, i]) is ∃ or ∀, and S is a system of real polynomial equations and inequalities in (Subscript[t, 1]),…,(Subscript[t, m]), (Subscript[x, 1]),…,(Subscript[x, n]).

By Tarski's theorem, solution sets of quantified systems of real polynomial equations and inequalities are semialgebraic.

Cell decomposition

Every semialgebraic set can be represented as a finite union of disjoint cells defined recursively as follows.

• A cell in ℝ is a point or an open interval.

• A cell in (ℝ^k) has one of the two forms

{((Subscript[a, 1]),...,(Subscript[a, k]),(Subscript[a, k+1])):((Subscript[a, 1]),...,(Subscript[a, k]))∈(Subscript[C, k])∧(Subscript[a, k+1])==r((Subscript[a, 1]),...,(Subscript[a, k]))}

{((Subscript[a, 1]),...,(Subscript[a, k]),(Subscript[a, k+1])):((Subscript[a, 1]),...,(Subscript[a, k]))∈(Subscript[C, k])∧(Subscript[r, 1])((Subscript[a, 1]),...,(Subscript[a, k]))<(Subscript[a, k+1])<(Subscript[r, 2])((Subscript[a, 1]),...,(Subscript[a, k]))}

where (Subscript[C, k]) is a cell in (ℝ^k), r  is a continuous algebraic function, and (Subscript[r, 1]) and (Subscript[r, 2]) are continuous algebraic functions, -∞, or ∞, and (Subscript[r, 1])<(Subscript[r, 2]) on (Subscript[C, k]).

By an algebraic function we mean a function

(Subscript[Root, (Subscript[x, k+1]),p])f:(Subscript[C, k])∋((Subscript[x, 1]),...,(Subscript[x, k]))⟶(Subscript[Root, (Subscript[x, k+1]),p])f((Subscript[x, 1]),...,(Subscript[x, k]))∈ℝ

where

f==(Subscript[c, 0])(Subscript[x, k+1]^m)+(Subscript[c, 1]) (Subscript[x, k+1]^(m-1))+…+(Subscript[c, m])∈ℝ[(Subscript[x, 1]),…,(Subscript[x, k]),(Subscript[x, k+1])]

is a polynomial and

(Subscript[Root, (Subscript[x, k+1]),p])f((Subscript[x, 1]),...,(Subscript[x, k]))

is the p-th real root (multiplicities counted) of f treated as a univariate polynomial in (Subscript[x, k+1]).

Cylindrical algebraic decomposition

The CAD algorithm, introduced by Collins (1975), allows to compute a cell decomposition of any semialgebraic set presented by a quantified system of polynomial equations and inequalities.

• Cells are arranged cylindrically.

• Polynomials whose roots bound cells extending C are delineable on C.

ineqs=-5 (((y-3))^3)+7 (x^2) (y-3)<1&&(x^2)+2 (((y-4))^2)<21;

pl=RegionPlot[ineqs,{x,-4.5,4.5},{y,-0.5,7.75},AspectRatio->Automatic,PlotPoints->100];
cad=CylindricalDecomposition[ineqs,{x,y}];
bds=Select[cad/.{Or->List,And->List,HoldPattern[Inequality[a_,__,b_]]:>{a,b}},ListQ[#[[1]]]&];
pts=Sort[Union@@(First/@bds),Less];
cbs=Table[{Graphics[{c=Hue[(i-1)/7],Thickness[0.005],Line[{{pts[[i]],0},{pts[[i+1]],0}}]}],Plot[bds[[i,2]],{x,bds[[i,1,1]],bds[[i,1,2]]},ColorFunction->(c&),PlotStyle->Thickness[0.005],PlotRange->{{-4.5,4.5},{-0.5,7.75}}]},{i,7}];
Show[{pl,cbs,Graphics[{Thickness[0.002],Gray,Line[{{#,-0.5},{#,7.75}}]&/@pts}]},Axes->True,Frame->False]

Cylindrical algebraic formulas

Cell decompositions computed by the CAD algorithm can be represented directly as cylindrical algebraic formulas (CAF).

CAFView[CylindricalDecomposition[(x^2)+(y^2)+(z^2)≤1,{x,y,z}],{x,y,z}]

Direct applications of the CAF representation

The CAF representation of a semialgebraic set A can be used to:

• decide whether A is nonempty,

• find the minimal and maximal values of the first coordinate of elements of A,

• generate an arbitrary element of A,

• find a graphical representation of A,

• compute the volume of A,

• compute multidimensional integrals over A,

• ...

Using CAF as input to the CAD algorithm

• Compute set theoretic operations on semialgebraic sets given in the CAF representation [A.S. ISSAC 2010].

• Solve (possibly quantified) systems of polynomial equations and inequalities in a semialgebraic set given in the CAF representation [this talk].

Region CAD

Algorithm RegionCAD
Input: A system S((Subscript[x, 1]),...,(Subscript[x, n])) of polynomial equations and inequalities, a cylindrical algebraic formula F((Subscript[x, 1]),...,(Subscript[x, n])), and a sequence of quantifiers (Subscript[Q, 1]),...,(Subscript[Q, k]), with 0≤k≤n.
Output: A cylindrical algebraic formula G((Subscript[x, 1]),...,(Subscript[x, n-k])) equivalent to (Subscript[Q, 1])(Subscript[x, n-k+1])...(Subscript[Q, k])(Subscript[x, n])F((Subscript[x, 1]),...,(Subscript[x, n]))∧S((Subscript[x, 1]),...,(Subscript[x, n])).

• Constructs only cells that belong to the solution of F.

• Smaller projection set than when computing with polynomial descriptions.

• Additional quantifier “for all x in the solution set of F”

Rounded Astroid

astroid=y^6+3 x^2 y^4-3 y^4+3 x^4 y^2+21 x^2 y^2+3 y^2+x^6-3 x^4+3 x^2≤1;

pla=RegionPlot[astroid,{x,-1.5,1.5},{y,-1.5,1.5}]

(ra=CylindricalDecomposition[Exists[{x,y},astroid&&(x-u)^2+(y-v)^2≤1/100],{u,v}])//Timing

plra=RegionPlot[ra,{u,-1.5,1.5},{v,-1.5,1.5}]

Polynomial description of the rounded astroid

Too hard to compute (not implemented in Mathematica, QEPCAD fails after 2 hours due to "Prime list exhausted")

Solve an optimization problem over the rounded astroid: for what values of c the ellipse (1/4) (((u-1))^2)+(1/9) (((v-1))^2)==c intersects the rounded astroid.

(rnc=Reduce`CADExtend[ecaf2,(u-1)^2/4+(v-1)^2/9==c,{c,u,v},{Exists,Exists}])//Timing

Picture

Show[{plra,pla,Graphics[{{Red,PointSize[Large],Point[{1,1}]},Circle[{1,1},Sqrt[First[rnc]]{2,3}],Circle[{1,1},Sqrt[Last[rnc]]{2,3}]}]},PlotRange->All,AspectRatio->Automatic]

Rose Surface

rose=-16 x^8+80 x^10-128 x^12+64 x^14+9 x^4 y^2-64 x^6 y^2+400 x^8 y^2-768 x^10 y^2+448 x^12 y^2-6 x^2 y^4-96 x^4 y^4+800 x^6 y^4-1920 x^8 y^4+1344 x^10 y^4+y^6-64 x^2 y^6+800 x^4 y^6-2560 x^6 y^6+2240 x^8 y^6-16 y^8+400 x^2 y^8-1920 x^4 y^8+2240 x^6 y^8+80 y^10-768 x^2 y^10+1344 x^4 y^10-128 y^12+448 x^2 y^12+64 y^14+48 x^8 z^2-160 x^10 z^2+128 x^12 z^2-36 x^4 y^2 z^2+192 x^6 y^2 z^2-800 x^8 y^2 z^2+768 x^10 y^2 z^2+24 x^2 y^4 z^2+288 x^4 y^4 z^2-1600 x^6 y^4 z^2+1920 x^8 y^4 z^2-4 y^6 z^2+192 x^2 y^6 z^2-1600 x^4 y^6 z^2+2560 x^6 y^6 z^2+48 y^8 z^2-800 x^2 y^8 z^2+1920 x^4 y^8 z^2-160 y^10 z^2+768 x^2 y^10 z^2+128 y^12 z^2-48 x^8 z^4+80 x^10 z^4+54 x^4 y^2 z^4-192 x^6 y^2 z^4+400 x^8 y^2 z^4-36 x^2 y^4 z^4-288 x^4 y^4 z^4+800 x^6 y^4 z^4+6 y^6 z^4-192 x^2 y^6 z^4+800 x^4 y^6 z^4-48 y^8 z^4+400 x^2 y^8 z^4+80 y^10 z^4+16 x^8 z^6-36 x^4 y^2 z^6+64 x^6 y^2 z^6+24 x^2 y^4 z^6+96 x^4 y^4 z^6-4 y^6 z^6+64 x^2 y^6 z^6+16 y^8 z^6+9 x^4 y^2 z^8-6 x^2 y^4 z^8+y^6 z^8==0;
ball=(x-1)^2+y^2+z^2≤1/25;

Show[{ParametricPlot3D[{Cos[s]Cos[3/4t]Cos[t],Cos[s]Cos[3/4t]Sin[t],Sin[s]},{t,0,10Pi},{s,0,2Pi},PlotStyle->Directive[LightRed,Opacity[0.25],Specularity[1,50]],Mesh->None,PlotPoints->30],Graphics3D[{Green,Opacity[0.5],Sphere[{1,0,0},0.2]}]},PlotRange->All]

Intersection with z == 0

Show[{ParametricPlot[{Cos[3/4t]Cos[t],Cos[3/4t]Sin[t]},{t,0,10Pi},Mesh->None,PlotPoints->30],Graphics[{Green,Opacity[0.5],Disk[{1,0},0.2]}]},PlotRange->All]

CylindricalDecomposition[ball&&rose,{x,y,z}]

Did not finish in 24 hours

(cafball=CylindricalDecomposition[ball,{x,y,z}])//Timing

(cafrose=CylindricalDecomposition[rose,{x,y,z}]);//Timing

(sol1=Reduce`CADCombine[And,{cafball,cafrose},{x,y,z}])//Timing

(sol2=Reduce`CADExtend[cafball,rose,{x,y,z}])//Timing

sol1===sol2

Main ideas of the modified CAD algorithm

Algorithm RegionCAD
Input: A system S((Subscript[x, 1]),...,(Subscript[x, n])) of polynomial equations and inequalities, a cylindrical algebraic formula F((Subscript[x, 1]),...,(Subscript[x, n])), and a sequence of quantifiers (Subscript[Q, 1]),...,(Subscript[Q, k]), with 0≤k≤n.
Output: A cylindrical algebraic formula G((Subscript[x, 1]),...,(Subscript[x, n-k])) equivalent to (Subscript[Q, 1])(Subscript[x, n-k+1])...(Subscript[Q, k])(Subscript[x, n])F((Subscript[x, 1]),...,(Subscript[x, n]))∧S((Subscript[x, 1]),...,(Subscript[x, n])).

Projection phase

We start with the set of factors of polynomials that appear in S and use ProjMC(_, (Subscript[B, k])) (or ProjH(_,(Subscript[B, k]))) as the projection operator,  where (Subscript[B, k]) is the set of polynomials that appear in F as bounds on (Subscript[x, k]).

• ProjMC(P)={(Subscript[coeffs, (Subscript[x, k])])(f):f∈P}⋃{(Subscript[disc, (Subscript[x, k])])(f):f∈P}⋃{(Subscript[res, (Subscript[x, k])])(f,g):f≠ g∈P}

• ProjMC(P,Q)=ProjMC(P)⋃{(Subscript[res, (Subscript[x, k])])(f,g):f∈P, g∈Q\P}

• ProjH(P)={(Subscript[coeffs, (Subscript[x, k])])(f):f∈P}⋃{(Subscript[psc, m,(Subscript[x, k])])((Subscript[r, i])(f),(Subscript[r, i])(f)'):f∈P,2≤ i≤deg(f), 0≤m<deg(f)-1}⋃{(Subscript[psc, m,(Subscript[x, k])])((Subscript[r, i])(f),g):f≺g∈P,1≤ i≤deg(f), 0≤m<min(deg(f),deg(g))}

• ProjH(P,Q)=ProjH(P)⋃{(Subscript[psc, m,(Subscript[x, k])])(f,g):f∈P, g∈Q\P, 0≤m<min(deg(f),deg(g))}

Lifting phase

• Let 𝒟 be the cell decomposition described by F. In the k-th step of the lifting phase we construct only cells that are contained in (Subscript[Π, k])(𝒟).

• When lifting a cell C∈(Subscript[Π, k])(𝒟) we use only those projection polynomials that appear in the projection with respect to the bound polynomials of F over C.

Step 1

(Subscript[r, 1])<x<(Subscript[s, 1])

gr1=Graphics[{{PointSize[Large],Red,Point[{{-1,0},{2,0}}]},{PointSize[Large],Green,Point[{{0,0},{2/3,0}}]},{PointSize[Large],Blue,Point[{{-3,0},{-2,0},{5/2,0},{3,0},{11/3,0},{14/3,0}}]}},Axes->{True,False},PlotRange->{{-5,5},All}]

Step 2

(Subscript[r, 1])<x<(Subscript[s, 1])∧(Subscript[r, 2])(x)<y<(Subscript[s, 2])(x)

pl1=Plot[{x^3-2x^2-x-2,x^3-x^2-2x+3},{x,-1,2},Filling->{1->{2}}];
Show[{gr1,pl1,Graphics[{{Gray,Line[{{1/3,-5},{1/3,5}}]},{{PointSize[Large],Red,Point[{{1/3,-68/27},{1/3,61/27}}]},{PointSize[Large],Green,Point[{{1/3,-2},{1/3,-1},{1/3,3/2}}]},{PointSize[Large],Blue,Point[{{1/3,-4},{1/3,-10/3},{1/3,3},{1/3,4}}]}}}]}]

Example: incremental CAD construction

f=x^2+y^2+z^2-1;
g=x^2+y^2-z;
h=(x-1)^2+(y-1)^2+(z-1)^2-1;
vars={x,y,z};

Find f<0∧g<0∧h<0.

Show[{RegionPlot3D[f<0,{x,-2,2},{y,-2,2},{z,-1,2},PlotStyle->Directive[Red,Opacity[0.5]],Mesh->None],RegionPlot3D[g<0,{x,-2,2},{y,-2,2},{z,-1,2},PlotStyle->Directive[Orange,Opacity[0.5]],Mesh->None],RegionPlot3D[h<0,{x,-2,2},{y,-2,2},{z,-1,2},PlotStyle->Directive[Green,Opacity[0.5]],Mesh->None]},AspectRatio->Automatic]

F=CylindricalDecomposition[f<0&&g<0,vars]

(Subscript[f, 1])=-1+(x^2)+(x^4)+(1+2 (x^2)) ((#1)^2)+((#1)^4)&[y]

The defining polynomials of cell bounds in F:

(Subscript[B, 3])={f,g}
(Subscript[B, 2])={(Subscript[f, 1])}
(Subscript[B, 1])={-1+(x^2)+(x^4)}

Projection sets:

(Subscript[W, 3])={h}

(Subscript[W, 2])={(Subscript[p, 1]),(Subscript[p, 2]),(Subscript[p, 3])}={Discriminant[h,z],Resultant[f,h,z],Resultant[g,h,z]}

(Subscript[W, 1])={Discriminant[(Subscript[p, 1]),y],Discriminant[(Subscript[p, 2]),y],Discriminant[(Subscript[p, 3]),y],Resultant[(Subscript[p, 1]),(Subscript[p, 2]),y],Resultant[(Subscript[p, 1]),(Subscript[p, 3]),y],Resultant[(Subscript[p, 2]),(Subscript[p, 3]),y],Resultant[(Subscript[f, 1]),(Subscript[p, 1]),y],Resultant[(Subscript[f, 1]),(Subscript[p, 2]),y],Resultant[(Subscript[f, 1]),(Subscript[p, 3]),y]}

Reduce`CADExtend[F,h<0,vars]

• The computation of the solution set constructs 497 cells (including cells that are rejected because h < 0 evaluates to false on them). The computation of F constructs 111 cells, hence the total number of cells constructed while solving the problem is 608. A direct computation of CAD constructs 1585 cells.

•  The set of projection polynomials used in the incremental method is always a subset of the set of projection polynomials used in the direct CAD computation, because the defining polynomials of cell bounds are projection polynomials and Proj(P,Q)⊆Proj(P⋃Q). In this example the inclusion is strict. The direct CAD computation uses the following projection polynomial whose factor does not appear among factors of the projection polynomials used by the incremental method.

Resultant[Discriminant[f,z],Resultant[f,h,z],y]//Factor

factors=Union@@((First/@FactorList[#])&/@(Subscript[W, 1]))

MemberQ[factors,5-12 x+8 (x^2)]

Experimental results: Intersections of balls with centers uniformly distributed on the unit circle

Let
(Subscript[x, k,n]):=q(cos(2 k π/n))
(Subscript[y, k,n]):=q(sin(2 k π/n))
(Subscript[B, k,n]):={(x,y,z):(((x-(Subscript[x, k,n])))^2)+(((y-(Subscript[y, k,n])))^2)+(z^2)≤2}
where q(r) denotes a rational number such that |r−q(r)|<(10^(-6)).
Find a CAF representation of (Subscript[I, n]):=(⋂_(k=0)^(n-1)) (Subscript[B, k,n]).

CAD: Use the CAD algorithm to find a CAF representing (Subscript[I, n])
CC: Use the CAD algorithm to find CAFs representing (Subscript[B, k,n]) for  0≤k≤n−1 and then use the algorithm CAFCombine to compute a CAF representing their conjunction.
INCR: Use the CAD algorithm to find a CAF (Subscript[F, 0]) representing (Subscript[B, 0,n]) and then for  0≤k≤n−1 use the algorithm RegionCAD to find a CAF (Subscript[F, k]) representing the intersection of the solution set of (Subscript[F, k-1]) with (Subscript[B, k,n]).

q[r_]:=Rationalize[r,10^-6]
B[k_,n_]:=Sphere[{q[Cos[2k π/n]],q[Sin[2k π/n]],0},Sqrt[2]]
Graphics3D[With[{n=7},Table[{Hue[k/n],Opacity[0.2],B[k,n]},{k,0,n-1}]]]

Conclusions

• We can solve polynomial systems over semialgebraic sets given by a CAF.

• We can solve conjunctions of polynomial conditions in an incremental manner. For some problems this is faster than other methods. Unfortunately we do not have a precise characterisation of such problems. Nevertheless the method may be useful for problems that prove hard for a direct CAD computation or may be tried in parallel with other CAD computation methods.