Implicitization by Hybrid Methods

Daniel Lichtblau

Wolfram Research, Inc.
100 Trade Centre Dr.
Champaign IL USA, 61820

danl@wolfram.com

ADG 2004, Gainesville, FL
September 2004

Introduction

We investigate a family of trigonometric planar mappings of the form

x=cos(m a)+cos(m b)+cos(m (a-b))
y=cos(n a)+cos(n b)+cos(n (a-b))

for fixed integers m and n. These arose a few years ago in the context of a problem in extremal number theory that in turn gave rise to a lattice packing. The number theory problem was to compute, for given integer s,

f(s)=(Underscript[Min, (Subscript[a, j])ϵ ℝ]) Max [|(∑_(j=1)^s)(e^(i (Subscript[a, j]) ))|,| (∑_(j=1)^s)(e^(i s (Subscript[a, j]) ))|]

The relevant theorem is that this value is 1 for s<=2, and 0 for s>=4. The important case of s=3 was computed by Seon-Hong Kim as part of his Ph.D dissertation work; it used the map above with m=1 and n=3. It is approximately 1.2. At ACA 2003 I showed how one might compute it using symbolic techniques.

Here are some definitions we will need.

trigpoly[n_,a_,b_]=Cos[n a]+Cos[n b]+Cos[n (a-b)];
grad[expr_,a_,b_]:={(∂_a)expr,(∂_b)expr}
trigsubs={Cos[a]->(Subscript[c, a]),Sin[a]->(Subscript[s, a]),Cos[b]->(Subscript[c, b]),Sin[b]->(Subscript[s, b])};

Brief computation of the extremal value

Theorem [Kim 2003]: We may find f(3) by determining tangential intersections of the curves

cos(a)+cos(b)+cos(a-b)=k
cos(3a)+cos(3b)+cos(3(a-b))=k

for values {x,y,k}. [These new functions are half of (the squares of the moduli-3).] Moreover at extremal values of f(3) there will be only finitely many such intersection points, up to periodicity.

We illustrate level curves for values of k, -1.1 and -1.3, that bracket the minimizer.

Brief computation of the extremal value

The actual computation is a polynomial elimination from a set of polynomials involving a Lagrange multiplier. We forego the details here (they are in the paper) and just show it.

{e1,e2}=Map[TrigExpand[trigpoly[#,a,b]]&,{1,3}];
polys=Flatten[{e1-k,e2-k,grad[e1,a,b]-λ grad[e2,a,b],((Cos[a])^2)+((Sin[a])^2)-1,((Cos[b])^2)+((Sin[b])^2)-1}]/.trigsubs;
Timing[kpoly=First[GroebnerBasis[polys,k,{λ,(Subscript[c, a]),(Subscript[s, a]),(Subscript[c, b]),(Subscript[s, b])},MonomialOrder->EliminationOrder]]]

k/.NSolve[kpoly==0, k]

The envelope curves of the map

We form a picture of the planar map for the case m=1and n=3.

The minimal value is the last point in the lower left where this intersects the line y=x. The other real roots also have physical significance with respect to this picture. For example, {1,1} is a local maximizer.

The envelope curves

This brings us to the next question. The map itself looks interesting, and we might wish to know equations to describe the envelope curve. For example, is it indeed (as it appears) the intersection of two quartics of the form y=(Subscript[f, j])(x)? (No, it is not.) How can we find this envelope curve?

Various possibilities. We might look for where Jacobian vanishes. We might try to extremize one variable as a function of the other.

Moreover, we can work with the trig formulation, or use complex exponentials, or use the standard rationalization of trig parametrization.

We will extremize, using the trig formulation.

The envelope curves

The idea is to write x and y as functions of trigonometric polynomials, again make substitutions so these become explicitly algebraic, and realize that at the extrema (that is, on the boundary curve(s)), the gradients of the two functions must be parallel. (Reason: on the envelope, for fixed x=x(a,b) we extremize y=y(a,b). This can be set up as a standard Lagrange multiplier problem with ∇y=λ∇x.)

parameters={a,b};mainvars={x,y};
xpoly=x-trigpoly[1,a,b];ypoly=y-trigpoly[3,a,b];
xypolys={xpoly,ypoly};
gradientpolys=grad[xpoly,a,b]-λ grad[ypoly,a,b];
trigidentities={((Subscript[c, a])^2)+((Subscript[s, a])^2)-1,((Subscript[c, b])^2)+((Subscript[s, b])^2)-1};
elimvars={λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])};
polys=TrigExpand[Join[trigidentities,xypolys,gradientpolys]]/.trigsubs

The envelope curves

Now what? We now want to compute the polynomial in {x,y} in the elimination ideal. For that we can simply form a Gröbner basis with an appropriate term order.

Timing[gb=GroebnerBasis[polys,mainvars,elimvars,Sort->True,MonomialOrder->EliminationOrder]]

envelope=Map[#[[1]]^#[[2]]&,Drop[FactorList[First[gb]],1]]

The envelope curves

Let's have a look at those factors.

ImplicitPlot[Evaluate[Thread[envelope==0]],{x,-1.5,3.2},{y,-2,3.2},AspectRatio->7/8,PlotStyle->{{Thickness[.01],GrayLevel[.5],Dashing[{.15,.04}]},{Automatic}}];

The general problem

On to the general problem. These are interesting planar maps from which one might wish to find envelope curves. We illustrate the n=5 case below using around 4000 points.

x[a_,b_]=trigpoly[1,a,b];y[a_,b_]=trigpoly[5,a,b];
pointlist[m_]:=Table[2 π RandomReal[{0,1},2],{(2^m)}]points2k=pointlist[12];data[{point__}]:={x[point],y[point]};datalist2k=data/@points2k;lplot=ListPlot[datalist2k]

For n=4 we already cannot find the envelope curves in reasonable time using the method shown. The program Fermat handled a formulation that case in about 15 minutes (while we were at dinner). For n=5 it ran out of memory. So...we need hybrid methods.

Finding isolated points on the curve

Again we set up an appropriate polynomial system.

xpoly=x-trigpoly[1,a,b];ypoly=y-trigpoly[5,a,b];
xypolys={xpoly,ypoly};
gradientpolys=grad[xpoly,a,b]-λ grad[ypoly,a,b];
polys=TrigExpand[Join[trigidentities,xypolys,gradientpolys]]/.trigsubs;

We'll begin by specifying one variable and solving for the other (we'll need to get many points, but this is a start). For this purpose we use a numeric solver.

allbutyvars={x,λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])};
polyy=polys/.y->11/10;
Timing[soln=NSolve[polyy,allbutyvars,Sort->True];]

Select[Union[x/.soln,SameTest->(Abs[#1-#2]/(Abs[#1]+Abs[#2])<10^(-5)&)],Im[#]==0&]

Wasn't very easy. And for higher degrees it will only get worse. We need another method.

Finding isolated points on the curve

What we really need is a way to find roots locally. Alas, for this we need reasonable initial guesses as to values of parameters, which we do not have. We'll instead formulate this step as an optimization problem, with parameters taking values in known ranges. This is why we chose a trig formulation, because the trig variables are known to take on values between -1 and 1. The objective function we "optimize" is a trivial constant function; we are simply performing a constraint satisfaction.

Here is the set of constraints that must be set to zero to solve for x when y is 11/10.

Short[polyy,8]

Finding isolated points on the curve

We solve and get our value for x as well as the parameters. We are looking for a solution with x near 1.

Timing[{val,root}=NMinimize[{1,Flatten[{Thread[polyy==0],.9<= x<=1.2}]},{{x,.9,1.2},{λ,-10,10},{(Subscript[c, a]),-1,1},{(Subscript[c, b]),-1,1},{(Subscript[s, a]),-1,1},{(Subscript[s, b]),-1,1}},MaxIterations->200,Method->"DifferentialEvolution"]]

Finding isolated parts of the curve

We got more than we anticipated; we see that this part of the envelope curve arises from a segment in parameter space where b=0 (because sin(b) apparently vanishes). How lucky for us!

envelopepiece1=First[GroebnerBasis[Take[polys,4]/.{(Subscript[c, b])->1,(Subscript[s, b])->0},{x,y},{ (Subscript[c, a]),(Subscript[s, a])}]]

envplot1=ImplicitPlot[Evaluate[envelopepiece1==0],{x,-2,3.2},{y,-2,3.2},PlotStyle->{Thickness[.01],GrayLevel[.5],Dashing[{.15,.07}]},PlotPoints->1000];

Finding isolated parts of the curve

A similar lucky computation gives us the part of the envelope curve in the lower right.

polyy=polys/.y->-11/10;
Timing[{val,root}=NMinimize[{1,Flatten[{Thread[polyy==0],1.2<= x<=2}]},{{x,1.2,2},{λ,-10,10},{(Subscript[c, a]),-1,1},{(Subscript[c, b]),-1,1},{(Subscript[s, a]),-1,1},{(Subscript[s, b]),-1,1}},MaxIterations->200,Method->"DifferentialEvolution"]]

Here clearly we have b=-a in the domain (cosines are equal, sines are opposite).

envelopepiece2=First[GroebnerBasis[Take[polys,4]/.{(Subscript[c, b])->(Subscript[c, a]),(Subscript[s, b])->-(Subscript[s, a])},{x,y},{ (Subscript[c, a]),(Subscript[s, a])}]]

Finding isolated parts of the curve

envplot2=ImplicitPlot[Evaluate[envelopepiece2==0],{x,-2,3.2},{y,-2,3.2},PlotStyle->{Thickness[.012],GrayLevel[.35],Dashing[{.12,.08}]},PlotPoints->1000];

The hybrid computation

Now we have to get the last piece of boundary curve. For this we will require real work. We start by grabbing a point near the lower crossing with the y axis. I found it expedient to change the optimization slightly so as to better balance the constraints (I now square some of them).

polyx1=polys/.x->0;
Timing[{val,root1}=NMinimize[{1,Flatten[{Thread[polyx1^2==0],-1.6<= y<=-1.2}]},{{y,-1.6,-1.2},{λ,-10,10},{(Subscript[c, a]),-1,1},{(Subscript[c, b]),-1,1},{(Subscript[s, a]),-1,1},{(Subscript[s, b]),-1,1}},MaxIterations->200,Method->"DifferentialEvolution"]]

The hybrid computation

Now we must find more envelope points in order to reconstruct the curve therefrom. We begin with the solution point we just found. We refine the result to high precision using a local root finder.

allbutxvars={y,λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])};
rootbetter1=FindRoot[polyx1==0,Evaluate[Apply[Sequence,Transpose[{allbutxvars,allbutxvars/.root1}]]],PrecisionGoal->30,AccuracyGoal->30,WorkingPrecision->40]

The hybrid computation

To deduce the actual curve we will require a number of {x,y} pairs of points lying on it. We will treat this as a homotopy continuation problem wherein the roots found above will be initial values for a differential system. We create the system below but display it in abbreviated form due to excessive size.

varsinx=Through[allbutxvars[x]];
parampolys=polys/.Thread[allbutxvars->varsinx];
diffpolys=D[parampolys,x];
Short[odesystem1=Join[Thread[diffpolys==0],Thread[(varsinx/.x->0)==(allbutxvars/.rootbetter1)]],12]

The hybrid computation

Now we solve it.

desoln1=NDSolve[odesystem1,varsinx,{x,0,1.2}];

The message tells us we were not able to get very far. Most likely we are near a tangential intersection of solution curves. But what we have suffices to give us points on the envelope curve. We refine to high precision.

solnfunctions1=First[varsinx/.desoln1];
valuelists1=Table[Join[{t->x},Thread[List[allbutxvars,solnfunctions1]]],{x,0,3/10,5/1000}]/.t->x;
highprecsolns1=Table[Flatten[{First[valuelists1[[j]]],FindRoot[Evaluate[(polys/.First[valuelists1[[j]]])==0],Evaluate[Apply[Sequence,Rest[valuelists1[[j]]]]],PrecisionGoal->100,AccuracyGoal->100,WorkingPrecision->150]}],{j,Length[valuelists1]}];

The hybrid computation

We repeat the process of finding one solution, refining it, approximating a part of the curve through it, and refining individual points thereon.

polyx2=polys/.x->-2/5;
Timing[{val,root2}=NMinimize[{1,Flatten[{Thread[polyx2^2==0],-1.2<= y<=-.9}]},{{y,-1.2,-.9},{λ,-10,10},{(Subscript[c, a]),-1,1},{(Subscript[c, b]),-1,1},{(Subscript[s, a]),-1,1},{(Subscript[s, b]),-1,1}},MaxIterations->200,Method->"DifferentialEvolution"]]

The hybrid computation

rootbetter2=FindRoot[polyx2==0,Evaluate[Apply[Sequence,Transpose[{allbutxvars,allbutxvars/.root2}]]],PrecisionGoal->30,AccuracyGoal->30,WorkingPrecision->40];
odesystem2=Join[Thread[diffpolys==0],Thread[(varsinx/.x->-2/5)==(allbutxvars/.rootbetter2)]];
desoln2=NDSolve[odesystem2,varsinx,{x,-2/5,-1}];

solnfunctions2=First[varsinx/.desoln2];
valuelists2=Table[Join[{t->x},Thread[List[allbutxvars,solnfunctions2]]],{x,-2/5,-4/5,-5/1000}]/.t->x;
highprecsolns2=Table[Flatten[{First[valuelists2[[j]]],FindRoot[Evaluate[(polys/.First[valuelists2[[j]]])==0],Evaluate[Apply[Sequence,Rest[valuelists2[[j]]]]],PrecisionGoal->100,AccuracyGoal->100,WorkingPrecision->150]}],{j,Length[valuelists2]}];

The hybrid computation

With these values we can reconstruct the curve. We guess that it is a polynomial of degree at most 2 in y, and 12 in x. We form the possible monomials as our "basis" set, evaluate these at the {x,y} values corresponding to our high precision approximations of points on the envelope, and find the appropriate null space.

Short[xyvals=Join[highprecsolns1,highprecsolns2]/.((a_/;!MatchQ[a,x|y])->b_):>Sequence[],4]
xypowers=Flatten[Table[x^j*y^k,{j,0,12},{k,0,2}]];
valuevectors=xypowers/.xyvals;
ns=NullSpace[valuevectors];
Length[ns]

The hybrid computation

16 null vectors means we have too large a basis. We need to "squeeze" out a null vector in terms of lower degree monomials. We do this by Gaussian elimination.

ns2=Chop[ns,10^(-145)];
ns3=Map[Reverse,ns2];
rred=RowReduce[ns3];
rred2=Chop[rred,10^(-145)];
N[Last[rred2]]

We rationalize, find the common denominator, remove it by multiplication, and finally form the curve by associating the resulting vector of coefficients with the monomial basis.

coeffs=Rationalize[Last[rred2]];
mult=Apply[LCM,Denominator[coeffs]];
coeffs2=coeffs*mult;
envelopepiece3=coeffs2.Reverse[xypowers]

The hybrid computation

envplot3=ImplicitPlot[Evaluate[envelopepiece3==0],{x,-2,3.2},{y,-2,3.2},PlotStyle->{Thickness[.008],GrayLevel[.2],Dashing[{.12,.08}]},PlotPoints->1000];

The hybrid computation

ll=Show[{lplot,envplot1,envplot2,envplot3}]

Verification

To verify that our curve is part of the envelope, pick specific exact values and see if they are consistent with the original system. For example, we solve for y when x=-1/2.

xycoord=Prepend[First[Solve[(envelopepiece3/.x->-1/2)==0,y]],x->-1/2]

Plug into the system, check that we get a nontrivial Gröbner basis.

GroebnerBasis[polys/.xycoord,{λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])},MonomialOrder->DegreeReverseLexicographic]

If we do this with enough points we can prove that the entire curve satisfies the boundary curve system (standard result involving polynomial interpolation).

Alternative methods

(i) Use integer relation finding to deduce, from one point, the smallest polynomial equation it satisfies. A drawback is this will not be useful for approximate cases or cases where the polynomial factorizates nontrivially over an algebraic extension of the rationals.

(ii) Use image processing methods of edge detection to solve for {x,y} pairs on the envelope curve. Use series "inversion" to formulate a Puiseux series to obtain {a,b} as (multivalued) functions of {x,y} on that curve. These could then replace the optimization and ODE solving steps we used above. As those were computationally the most difficult and required the most user intervention, anything to make them easier would be advantageous.

Applications

I tend to envision the methods descibed here as useful for finding exact curves and possibly surfaces in cases where these are really needed e.g. problems from pure math. For approximation purposes there are usually better/faster approaches. A possible exception might be in finding relations satisfied by solutions to ODEs, in cases where they lie on algebraic varieties but are not explicitly presented as DAEs (for one, this knowledge can be useful for deducing long term evolution of such solutions).

Needless to say, I would be interested in hearing about other possible uses for this technology.

Summary

Starting with a problem in number theory we were led to consider the problem of finding bounding curves to certain trigonometric planar maps. At modest degree symbolic methods fail to handle this in reasonable time. We instead employed several numeric technologies, to wit, an optimizer (employed as a constraint satisfaction solver), a local root finder, an ODE solver, and numeric linear algebra, to find our curve. We then use a fast symbolic method to verify that it is correct. We illustrated on an example that had not been amenable to tactics from symbolic algebra.

There are several open questions. The general issues are
(i) How to improve performance?
(ii) To what classes of problems to apply these techniques?