Computation with 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 efficient computation of arbitary combinations of unions, intersections and complements of semialgebraic sets given in this representation. The algorithm can also be used to eliminate quantifiers from Boolean combinations of cylindrical algebraic formulas. We show application examples and an empirical comparison with direct CAD computation for unions and intersections of semialgebraic sets.

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), l]) (Underoverscript[∧, (j=1), m]) (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,

• ...

Computation with Cylindrical Algebraic Formulas

A natural question to ask is how to compute set-theoretic operations with semialgebraic sets represented by cylindrical algebraic formulas. Set-theoretic operations on solution sets of formulas correspond to Boolean operations on the formulas, so the question is how to compute a cell decomposition of a semialgebraic set given by a Boolean combination of cylindrical algebraic formulas. In principle this could be done using an extension of the CAD algorithm to arbitrary systems of equations and inequalities involving algebraic functions. However, this is inefficient, as it requires introducing a new variable for each algebraic function that appears in the input. In this talk I will present a customized version of the CAD algorithm which does not introduce new variables and makes efficient use of CAF structure, so that computation of cell decomposition for a semialgebraic set represented by a Boolean combination of CAFs is often faster than for the same set represented by a quantifier free polynomial system. The algorithm allows quantifiers, as long as the order of quantifier variables agrees with the order of variables in CAFs.

Algorithm CAFCombine
Input: Cylindrical algebraic formulas (Subscript[F, 1])((Subscript[x, 1]),...,(Subscript[x, n])),...,(Subscript[F, m])((Subscript[x, 1]),...,(Subscript[x, n])), a Boolean formula Φ((Subscript[p, 1]),...,(Subscript[p, m])) and a sequence of quantifiers (Subscript[Q, 1]),...,(Subscript[Q, k]), with 0≤k≤n.
Output: A cylindrical algebraic formula F((Subscript[x, 1]),...,(Subscript[x, n-k])) equivalent to (Subscript[Q, 1])(Subscript[x, n-k+1])...(Subscript[Q, k])(Subscript[x, n])Φ((Subscript[F, 1])((Subscript[x, 1]),...,(Subscript[x, n])),...,(Subscript[F, m])((Subscript[x, 1]),...,(Subscript[x, n]))).

Solotareff problem for n=4 (Collins 1995)

(Subscript[p, 1])=((r^2)-24 r+16<0||r>1)&&-1<u<v<1&&r-b>0;
(Subscript[p, 2])=(v^3)+r (v^2)-(v^2)-a v-r v+v-b+a+r-1==0&&4 (v^3)+3 r (v^2)-2 a v-b==0;
(Subscript[p, 3])=(u^3)+r (u^2)+(u^2)-a u+r u+u-b-a+r+1==0&&4 (u^3)+3 r (u^2)-2 a u-b==0;

Direct CAD computation.

TimeConstrained[CylindricalDecomposition[(Subscript[∃, u])(Subscript[∃, v])((Subscript[p, 1])&&(Subscript[p, 2])&&(Subscript[p, 3])),{r,a,b}],3600]

Computation using CAFCombine.

((Subscript[F, 1])=CylindricalDecomposition[(Subscript[∃, v])((Subscript[p, 1])&&(Subscript[p, 2])),{r,a,b,u}]);//Timing

((Subscript[F, 2])=CylindricalDecomposition[(Subscript[p, 3]),{r,a,b,u}]);//Timing

((Subscript[F, 1,2])=Reduce`CADCombine[(Subscript[∃, u])(#1&&#2)&,{(Subscript[F, 1]),(Subscript[F, 2])},{r,a,b}])//Timing

Unions and intersections of unit balls in (ℝ^3)

ball[a_]:=(((x-a))^2)+(((y-a))^2)+(((z-a))^2)≤1
balls[n_,d_]:=Table[ball[i d],{i,0,n-1}]
bb=balls[10,(1/2)]

Graphics3D[Table[{Hue[(i/5)],Opacity[0.3],Sphere[{i,i,i}]},{i,0,(9/2),(1/2)}]]

(caf1=CylindricalDecomposition[Or@@bb,{x,y,z}]);//Timing

(caf2=Reduce`CADCombine[Evaluate[Or@@(Slot/@Range[10])]&,CylindricalDecomposition[#,{x,y,z}]&/@bb,{x,y,z}]);//Timing

Reduce`CADCombine[Equivalent[#1,#2]&,{caf1,caf2},{x,y,z}]//Timing

CylindricalDecomposition[And@@bb,{x,y,z}]//Timing

Reduce`CADCombine[Evaluate[And@@(Slot/@Range[10])]&,CylindricalDecomposition[#,{x,y,z}]&/@bb,{x,y,z}]//Timing

Main ideas of the modified CAD algorithm

Algorithm CAFCombine
Input: Cylindrical algebraic formulas (Subscript[F, 1])((Subscript[x, 1]),...,(Subscript[x, n])),...,(Subscript[F, m])((Subscript[x, 1]),...,(Subscript[x, n])), a Boolean formula Φ((Subscript[p, 1]),...,(Subscript[p, m])) and a sequence of quantifiers (Subscript[Q, 1]),...,(Subscript[Q, k]), with 0≤k≤n.
Output: A cylindrical algebraic formula F((Subscript[x, 1]),...,(Subscript[x, n-k])) equivalent to (Subscript[Q, 1])(Subscript[x, n-k+1])...(Subscript[Q, k])(Subscript[x, n])Φ((Subscript[F, 1])((Subscript[x, 1]),...,(Subscript[x, n])),...,(Subscript[F, m])((Subscript[x, 1]),...,(Subscript[x, n]))).

• In the projection phase use the polynomials that appear in (Subscript[F, 1]),...,(Subscript[F, m]).

• In the lifting phase evaluate (Subscript[F, 1]),...,(Subscript[F, m]) at sample points level-by-level. If (Subscript[F, i]) evaluates to a Boolean constant, replace (Subscript[p, i]) with the constant in Φ((Subscript[p, 1]),...,(Subscript[p, m])).

• (Improvement) Find a decomposition of ℝ consistent with projections of the input sets and call CAFCombine separately for restrictions of (Subscript[F, 1]),...,(Subscript[F, m]) to each cell of the decomposition.

Note that for CAF’s constructed by CAD using Hong' s projection (and “usually” for CAF’s constructed using McCallum’s projection) polynomials that appear in the CAF’s are projection polynomials. Hence the size of projection computed from CAF’s is not larger than the size of projection computed from the polynomial systems that were used to construct the CAF’s. On the other hand, some polynomials that are present in the polynomial description may not appear in the CAF description.

In the lifting phase, cell construction may terminate early if, on the current cell,  Φ  evaluates to a Boolean constant or to one of the (Subscript[p, i]).

Restriction of a CAF to a cell in ℝ may involve less polynomials than the whole CAF, hence we may get smaller projection sets. On some cells in ℝ some of the (Subscript[F, i]) may evaluate to constants and hence  Φ((Subscript[F, 1]),...,(Subscript[F, m])) may reduce to a constant or to one of the (Subscript[F, i]).

Example 1

ineq1=(((x+1))^4)+(y^4)+(z^4)<4&&(((x+2))^2)+(y^2)+(z^2)<5;
ineq2=(((x-1))^4)+(y^4)+(z^4)<4&&(((x-2))^2)+(y^2)+(z^2)<5;

Show[{RegionPlot3D[ineq1,{x,-3,3},{y,-2,2},{z,-2,2},PlotStyle->Directive[Red,Opacity[0.5]],Mesh->None,PlotPoints->50],RegionPlot3D[ineq2,{x,-3,3},{y,-2,2},{z,-2,2},PlotStyle->Directive[Blue,Opacity[0.5]],Mesh->None,PlotPoints->50]}]

Example 1 (CAF evaluation)

caf1=CylindricalDecomposition[ineq1,{x,y,z}];
caf2=CylindricalDecomposition[ineq2,{x,y,z}];

CAFView[caf1,{x,y,z}]

CAFView[caf2,{x,y,z}]

CAFView[caf1,{x,y,z},{x->1/5}]

CAFView[caf2,{x,y,z},{x->1/5}]

CAFView[caf1,{x,y,z},{x->1/5,y->1/3}]

CAFView[caf2,{x,y,z},{x->1/5,y->1/3}]

CAFView[caf1,{x,y,z},{x->1/5,y->1/3,z->1}]

CAFView[caf2,{x,y,z},{x->1/5,y->1/3,z->1}]

Example 2

rose5[x_,y_]:=(x^6)-5 (x^4) y+3 (x^4) (y^2)+10 (x^2) (y^3)+3 (x^2) (y^4)-(y^5)+(y^6)
ineq1=(((x+2))^2)+(y^2)<7&&rose5[x+2,y]>0;
ineq2=(((x-2))^4)+(y^4)<21&&rose5[x-2,y]>0;
pl=RegionPlot[{ineq1,ineq2},{x,-5,5},{y,-3,3},PlotPoints->100,AspectRatio->Automatic]

((Subscript[F, 1])=CylindricalDecomposition[ineq1,{x,y}]);//Timing

((Subscript[F, 2])=CylindricalDecomposition[ineq2,{x,y}]);//Timing

(caf1=Reduce`CADCombine[#1||#2&,{(Subscript[F, 1]),(Subscript[F, 2])},{x,y}]);//Timing

(caf2=CylindricalDecomposition[ineq1||ineq2,{x,y}]);//Timing

caf1===caf2

Example 2 (step-by-step)

pts=Sort[Select[Union@@List@@@(First/@Join[List@@(Subscript[F, 1]),List@@(Subscript[F, 2])]),NumericQ],Less];
spts=Join[{Floor[First[pts]-1/2]},Table[If[OddQ[i],pts[[(i+1)/2]],Rationalize[(pts[[i/2+1]]+pts[[i/2]])/2,0.01]],{i,2Length[pts]-1}],{Ceiling[Last[pts]+1/2]}];
Show[{pl,Graphics[{{Gray,Line[{{#,-3},{#,3}}]&/@pts},{Green,Point[{#,0}]&/@spts}}]}]

eval[caf_,x_,c_]:= caf/.{x==e_:>c==e,Inequality[a_,Less,x,Less,b_]:>a<c<b}
Fx={eval[(Subscript[F, 1]),x,#],eval[(Subscript[F, 2]),x,#]}&/@spts;
Fx[[7]]

Position[Fx,_?(!MemberQ[#,False]&),1,Heads->False]

a=spts[[16]]<x<spts[[18]]

Fx[[17]]

(Subscript[g, 1])=y^2-(3-4 x-(x^2));
(Subscript[g, 2])=-5-32 x+24 (x^2)-8 (x^3)+(x^4)+(y^4);
(Subscript[d, 1])=Reduce[Discriminant[(Subscript[g, 1]),y]==0,x,Reals]

a/.{ToRules[(Subscript[d, 1])]}

(Subscript[d, 2])=Reduce[Discriminant[(Subscript[g, 2]),y]==0,x,Reals]

a/.{ToRules[(Subscript[d, 2])]}

(Subscript[r, 1,2])=Reduce[Resultant[(Subscript[g, 1]),(Subscript[g, 2]),y]==0,x,Reals]

a/.{ToRules[(Subscript[r, 1,2])]}

Show[{pl,Graphics[{{Gray,Line[{{#,-3},{#,3}}]&/@pts},{Green,Point[{#,0}]&/@spts},{Red,Line[{{(Subscript[r, 1,2])[[1,2]],-3},{(Subscript[r, 1,2])[[1,2]],3}}]}}]}]

Experimental results

Unions and intersections of n unit balls centered at k(v,v,v), for 0≤k<n.
CAD - direct CAD computation,
CC - CAFCombine without the improvement
CCI - CAFCombine with the improvement
CAF - computation of CAF’s representing individual balls

Conclusions

Experiments show that for some quantifiers elimination problems partitioning the problem into subproblems, getting CAF descriptions for solution sets of the subproblems and then combining the results using CAFCombine is faster than either direct quantifier elimination using CAD or getting polynomial descriptions for solution sets of the subproblems and then computing the CAD of the combined results. Also, for computation of  CAD for a set-theoretic combination of semialgebraic sets, computing a CAF representing each set and using CAFCombine may be faster than direct CAD computation from polynomial description.The advantage CAFCombine is greater if the number of intersections between the sets is smaller. In all examples the version of CAFCombine algorithm with the improvement performed significantly better than the version without the improvement.

Possible improvements

• During the lifting phase of the algorithm, input formulas evaluated at the current cell may not contain all polynomials that were present in the original formulas. In such cases, a smaller projection set may be computed using only polynomials that are still present.

• The algorithm could be extended to handle Boolean combinations of CAF’s and polynomial systems.