Frobenius Numbers by Toric Gröbner Bases

Daniel Lichtblau

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

danl@wolfram.com

ACA 2005, Nara, Japan
August 2005

Abstract

Given a set A={(Subscript[a, 1]),…,(Subscript[a, n])} of positive integers with gcd 1, it is not hard to show that all "sufficiently large" integers can be represented as a nonnegative integer combination of elements of A. The Frobenius number of the set is defined as the largest integer not so representable. The Frobenius instance problem (also called the money changing or postage stamp problem) is to determine, given a positive integer M, a set of nonnegative integers X={(Subscript[x, 1]),…,(Subscript[x, n])} such that X.A=n, or else show no such set exists. We will briefly recall how this can be addressed via toric Gröbner bases.

It is known that the Frobenius number problem is NP-hard in general. For dimension 2 it is trivial (Sylvester solved in two decades before Frobenius publicized the problem). In dimension 3 a very efficient method was found independently by Greenberg and Davison. For higher dimensions some quite effective methods are known for the case where one element of A is not too large (say, less than (10^7)).

Recent work has given rise to methods that are effective when the above restrictions do not hold, although the dimension must be bounded by 10 or so. It turns out that there is a way to recast this work using toric Gröbner bases, wherein the "fundamental domain" for the set A is given by the staircase of the basis with respect to a particular ordering. It is reasonably efficient in dimensions 4 to 7, when the elements in the set are as large as (10^40) or so. We will illustrate this.

Introduction: Background and brief history

We are given  a set A={(Subscript[a, 1]),…,(Subscript[a, n])} of positive integers with gcd(A)=1. We assume for later purposes that the set is in ascending order.

Problem 1 (Frobenius instance problem): Given a nonnegative integer M, find a set of nonnegative integers X={(Subscript[x, 1]),…,(Subscript[x, n])} such that X.A=n, or else show no such set exists.

Problem 2 (Frobenius number problem): Find the largest integer not representable as a nonnegative integer combination of A.

In the 80's and 90's Greenberg and Davison independenly found an ultrafast method for problem 2 when n=3. Beyond this size no specialized (that is, dimension specific) methods are known, and we must resort to general tactics.

Reasonably effective methods based mostly on graph theory have appeared also in the past 30 years or so. Some very nice new ones are presented, along with older approaches, in very recent work by Beihoffer, Hendry, Nijenhuis, and Wagon. It should be noted that the third author helped originate the graph theory approach. These methods are limited by the size of (Subscript[a, 1]), but not by n.

Introduction: Background and brief history

This restriction apparently rankled the fourth author, who continued to pursue the problem using different tactics. Forthcoming joint work by David Einstein, Adam Strzebonski, Stan Wagon, and myself will show how one can attack this problem effectively using lattice methods and integer programming. While we can do away with the size restriction on (Subscript[a, 0]) we do get into some algorithmic complexity due to dimension, and at present cannot get beyond n=11 or so. Some of the technology we use is also applicable to problem 1, where we can manage higher dimension (25 or larger).

It so happens that much of this can be recast in a setting of toric varieties. While problem 1, which boils to integer linear programming, has long been known to be amenable to such an approach (as per work by Conti and Traverso), this is apparently a new tactic for finding Frobenius numbers. We can exploit it to handle examples that, to the best of my knowledge, could not be done by methods known as of a year ago. A nice benefit we will soon see is that the needed code is quite short (three pages or so).

We will define a "fundamental domain" which is a generalization by Wagon that subsumes both lattice diagrams in earlier literature and a graph description from Beihoffer et al. The Frobenius number will be the furthest corner from the origin, in a suitably weighted (Subscript[l, 1]) norm. As we will see, an important domain feature is what we term "elbows". The method to find those via a Gröbner basis "staircase" constitutes the new material in this talk.

Solving a Frobenius instance via toric Gröbner bases

Say we are given the set and value

A={200,230,528,863,905,1355,1725,1796,1808};
b=7777;

We wish to know whether or how we can write b as a nonnegative combination of elements of A. We may do this as follows. Create a variable t that will be raised to the powers in this set. Create a variable for each set element (Subscript[a, j]) and equate it to (t^(Subscript[a, j])). We form a Gröbner basis in these variables, using an order that makes powers of t larger than all monomials that do not contain it. For this we use a weight matrix. It is structured in such a way as to be efficient for the task at hand. Basically it is like an elimination ordering on t with degree-reverse-lexicographic on the remaining variables, except instead of using (homogenious) total degree we weight by the values in A.

Solving a Frobenius instance via toric Gröbner bases

len=Length[A];
vars=Array[x,len];
polys=vars-t^A;
wtmat=RotateRight[Reverse[-IdentityMatrix[len+1]],2];
wtmat[[2]]*=-1;
wtmat[[1]]=Prepend[A,0];
wtmat[[{1,2}]]=wtmat[[{2,1}]];
wtmat

Solving a Frobenius instance via toric Gröbner bases

Timing[gb=GroebnerBasis[polys,Prepend[vars,t],MonomialOrder->wtmat];]

Length[gb]

To check whether we can represent b as a nonnegative integer combination of A, we reduce (using the same term ordering) (t^b) by this Gröbner basis.

red=PolynomialReduce[t^7777,gb,Prepend[vars,t],MonomialOrder->wtmat][[2]]

Solving a Frobenius instance via toric Gröbner bases

Now replace variables by their values and get corresponding exponents as the multipliers.

fax=Drop[FactorList[red],1]
exponvec=fax/.x[j_]:>A[[j]]

We check this. We want to see that 4*230+1*528+1*905+3*1808=7777.

Total[Apply[Times,exponvec,2]]

Remark: The above illustrates more or less the original formulation for handling ILPs via toric Gröbner bases. Subsequent improvements have appeared and quite likely this can be done much more efficiently now. This concludes our brief review of how one might solve a Frobenius instance problem using the method of Conti and Traverso. We return to the main task at hand, which is computation of Frobenius numbers.

The Fundamental Domain

We define a lattice by the set of integer combinations of {(Subscript[a, 2]),…,(Subscript[a, n])} that are zero modulo (Subscript[a, 1]). This is a full dimensional lattice in (ℤ^(n-1)). The set of residues of (ℤ^(n-1)) modulo this lattice gives rise to what we call the fundamental domain, which, as we see, lives in a space of dimension one less than the size of our set. It is not hard to see that there are (Subscript[a, 1]) distinct residue classes, so we know the cardinality of this domain. We now define the "weight" of a vector v∈(ℤ^(n-1)) as v.{(Subscript[a, 2]),…,(Subscript[a, n])}. It can be shown that every residue class has at least one element with all nonnegative entries. Among those, we choose one of minimal weight. In case of tie, choose the one that is lexicographically last. This uniquely defines the set of residues that we take to comprise the fundamental domain.

The Fundamental Domain

This domain can be shown to have several interesting properties.

• It is a staircase. If it contains a lattice element then it contains all nonnegative vectors with any coordinate strictly smaller.
• It tiles (ℤ^(n-1)).
• It is a cyclic group ℤ/(Subscript[a, 1]) ℤ.
• It can be given a circulant graph structure. It is this structure that was utilized in various shortest-path graph methods. Old and new methods for this are discussed at length in the recent work by Beihoffer et al.

For our purposes the property of most interest is the first one. We can recover this staircase by computing a toric Gröbner basis.

Definitions related to the Fundamental Domain

From the staircase property the fundamental domain  has turning points we refer to as elbows. It has extremal points called corners. Specifically, a corner is a point c in the domain, such that c+(Subscript[e, j]) is not in the domain, where (Subscript[e, j]) is the (Subscript[j, th]) coordinate vector. An elbow is a point x that is not in the domain, but is such that, for each j, either (Subscript[x, j])=0 or x-(Subscript[e, j]) is in the domain.

There are two other definitions that play a role in the algorithm. We will not descibe them too carefully but, roughly, there are as follows.
(i) Protoelbows. These have both positive and negative coordinates and correspond to certain "minimal" equivalences (that is, reducing relations) in the lattice. In Gröbner basis terms, these are given as exponent vectors of binomial pairs in the basis.
(ii) Preelbows. These are the "positive parts" of the protoelbows. Elbows are minimal elements in the partially ordered (ascending by inclusion) set of preelbows.

Fundamental Domain, illustrated

Since this region is one dimension smaller than the input set (as verything is done modulo (Subscript[a, 1])), it can be illustrated for the cases n=3 and n=4. The illustrations are due to Stan Wagon.

With respect to these domains, the Frobenius number corresponds to the farthest corner from the origin, with distance an (Subscript[l, 1]) metric weighted by element sizes. In the planar diagram the "elbows" are the lattice points on the axes that bound the diagram, and the lattice point in the interior just outside the "ell". The corners are the two extremal points reached by intersecting vertical and horizontal lines through the elbows. This picture tells the entire story as regards the n=3 case, because it can be shown that there is at most one interior elbow and two such corners, and finding them is easy.

Fundamental Domain, illustrated

In the three dimensional diagram the elbows are again the axial "bounding lattice points" as well as bounding points where the staircase goes up in the coordinate planes and in the interior. They are demarcated by yellow tetrahedra. The green boxes are corners and the blue box is the maximal corner.

The Algorithm, in brief

The gist of our efficient algorithm is to use integer linear programming to find certain sets containing elbows, then use a method David and Stan devised to go from elbows to corners. We finish when we have the furthest corner. While pathological cases (too many elbows) arise even at n=4, the average case performance (e.g. random examples) is quite nice. Some of the ILP ideas appear in work by Aardal, Hurkins and Lenstra with subsequent refinement by Aardal and Lenstra. I also had techniques similar to those in AHL which I used to find large examples of what are known Keith numbers.

Getting back to Frobenius numbers, it turns out that one can use toric Gröbner bases to find what are called "protoelbows". These are lattice points where positive combinations of one subset equal positive combinations of another. From these we get a superset of the elbows, and we use a domination algorithm of Bentley, Clarkson, and Levine to get the "kernel" set which comprise the actual elbows. The hardest step, algorithmically, is in coming up with the protoelbows. This can be done with a toric Gröbner basis.

As with instance solving, the idea is again to set up relations (that is, generating polynomials for a toric ideal) of the form (Subscript[x, j])-(t^(Subscript[a, j])), and eliminate t. Except we do not do this at all...

Finding "elbow" relations

...For better efficiency we now use an improvement, due to Pottier, in the Gröbner basis computation of the elimination ideal. Instead of (Subscript[x, j])-(t^(Subscript[a, j])) we work with polynomials ((Subscript[x, j])^(SuperPlus[(Subscript[e, j])]))-((Subscript[x, j])^(SuperMinus[(Subscript[e, j])])) where (SuperPlus[(Subscript[e, j])]) and (SuperMinus[(Subscript[e, j])]) are the positive and negative parts respectively of the (Subscript[j, th]) generator for the null space of A. We can use any basis for this purpose, so we chose one that is lattice reduced so as to keep down the exponent sizes. We augment by a polynomial u-∏(Subscript[x, j]) that in effect allows us to invert negative exponents, and use a monomial ordering that eliminates the new variable u.

From this basis we next want to find a new one that will enforce the "lexically last" provision of the definition. As our lattice is now represented by exponent vectors in a toric ideal, this amounts to an inverse lexicographical term ordering. But this is not a well-founded ordering for monomials because it has constants larger than any power products (it is an ordering appropriate for a local ring). For this we use a standard tactic of homogenizing, making the homogenizer variable largest, and working with a degree based term order. As we want an inverse lexicographic order we use the customary graded reverse-lexicographic. So we now compute a new Gröbner basis with respect to this order and then dehomogenize. Note that this second basis computation is quite fast compared to the first, hence not problematic in regard to efficiency.

Finding "elbow" relations

Since it is the exponent vectors we are after, we extract them from the basis. We convert to lattice elements by subtracting second term powers from first (this may sound complicated but it is just the usual translation from toric ideal to lattice terminology). Recalling that we are only interested in equalities modulo (Subscript[a, 1]), we strip off the first terms of the lattice vectors. To satisfy a technical consideration for protoelbows (see our full paper) we may need to negate the lattice vector. Taking positive parts of the resulting lattice vectors gives us our preelbows.

It is important to note that we are using the Gröbner basis for its staircase structure. Unlike the instance solving usage, we do not work geometrically with a polyhedron.

Code to find elbows

preElbows[vals_]:=Module[{n=Length[vals],vars,x,t,nspace,pos,neg,polys,exponvecs,h,y},
nspace=integerNullSpace[vals];
vars=Array[x,n];
pos=(nspace+Abs[nspace])/2;
neg=-nspace+pos;
polys=Map[Inner[Power,vars,#,Times]&,pos]-Map[Inner[Power,vars,#,Times]&,neg];
polys=Join[polys,{Apply[Times,vars]*y}-1];
polys=GroebnerBasis`ToricGroebnerBasis[polys,vars,y,MonomialOrder->EliminationOrder];
polys=homogenize[polys,vars,h];
polys=GroebnerBasis[polys,Reverse[Append[vars,h]],MonomialOrder->DegreeReverseLexicographic]/.h->1;
exponvecs=First[GroebnerBasis`DistributedTermsList[polys,vars]];
exponvecs=Map[First[-#[[1]]+#[[2]]]&,exponvecs];
exponvecs=Map[orient[#,vals]&,exponvecs];
exponvecs=Map[Rest,exponvecs];(*protoelbows*)Union[(exponvecs+Abs[exponvecs])/2]]

elbows[vals_]:=Sort[dominationKernel[preElbows[vals]]]

More code: utility functions

integerNullSpace[vec:{_Integer..}]:=Module[{mat,hnf},mat=Transpose[Join[{vec},IdentityMatrix[Length[vec]]]];
hnf=Last[Developer`HermiteNormalForm[mat]];
LatticeReduce[Map[Drop[#,1]&,Drop[hnf,1]]]]

homogenize[poly_,vars_,new_]/;(Head[poly]=!=Plus&&Head[poly]=!=List):=poly
homogenize[polys_List,vars_,v_]:=Map[homogenize[#,vars,v]&,polys]
homogenize[poly_,vars_,new_]:=Module[{degfunc,totdeg,j},degfunc=Apply[Plus,Map[Function[x,Exponent[#,x]&],vars]];
degfunc=Distribute[degfunc,Function,Plus];
totdeg=Max[Map[degfunc,Apply[List,poly]]];
Apply[Plus,Table[poly[[j]]*new^(totdeg-degfunc[poly[[j]]]),{j,Length[poly]}]]]

orient[vec_,basevec_]:=Module[{val=Rest[vec].Rest[basevec],j=1},Which[val>0,vec,val<0,-vec,True,While[vec[[j]]===0,j++];-Sign[vec[[j]]]*vec]]

dominationKernel[X_]:=
 Module[{Y=X[[Ordering[Map[Tr,X]]]],Z,i=2,len=1,k},Z=Table[{},{Length[Y]}];
  Z[[1]]=Y[[1]];
  While[i≤Length[Y]&&Tr[Y[[i]]]==Tr[Y[[1]]],len++;
   Z[[len]]=Y[[i]];i++];
  Do[k=1;
   While[k≤len,
    If[And@@Thread[Z[[k]]≤Y[[j]]],k=len+2;Break[],k++]];
   If[k==len+1,len++;Z[[len]]=Y[[j]]],{j,len+1,Length[Y]}];
  Sort[Take[Z,len]]]

ClearNegsAndDeleteZeroVector[vecs_]:=If[vecs=={},{},Union[DeleteCases[vecs/._?Negative:>0,Table[0,{Length[vecs[[1]]]}]]]];

Rest of code: corners

The method below for using elbows to find the "Frobenius corner" is largely due to David Einstein, with refinements and code provided by Stan Wagon.

farthest[corns_,A_]:=Fold[If[#2.Rest[A]>#1[[2]],{#2,#2.Rest[A]},#1]&,{Table[0,{Length[A]-1}],0},corns][[1]];

Options[FarthestCorner]={TopLevelTraceQ->False,TraceStep->1};

FarthestCorner[A_,elbows_?MatrixQ,opts___Rule]:=Module[{pts,p1,p2,B,cc,trQ,farvertex,maxes,A1},A1=Rest[A];
{tlQ,ts}={TopLevelTraceQ,TraceStep}/.{opts}/.Options[FarthestCorner];
skipct=counter=0;
currfar=farthest[elbows,A].Rest[A];
pts=Select[elbows,#[[1]]>0&];
Scan[(counter++;
{p1,p2}={First[#],Rest[#]};
B=Rest/@Select[elbows,#[[1]]<p1&];
B=dominationKernel[ClearNegsAndDeleteZeroVector[#-p2&/@B]];
maxes=Prepend[Max/@Transpose[B],0];
If[tlQ&&Mod[counter,ts]==0,Print[{counter,Length[pts],"Length of sent set ",Length[B],currfar,currfar-Total[A]}]];
If[(maxes+#).A1>currfar,cc=FarthestCornerSub[A,B,currfar,#,opts];
farvertex=farthest[Prepend[#+p2,p1]&/@cc,A];
currfar=Max[currfar,farvertex.A1],skipct++;
If[tlQ&&Mod[counter,ts]==0,Print[{counter,Length[pts],"got a cutoff, total skips",skipct}]]])&,pts];
currfar-Total[A]]/;Length[elbows[[1]]]>2;

FarthestCorner[A_,elbows_?MatrixQ,opts___Rule]:=farthest[FarthestCornerSub[A,elbows],A].Rest[A]-Total[A]/;Length[elbows[[1]]]==2;

FarthestCornerSub[A_,elbows_,currfar_,backdata_,opts___]:=Module[{p1,p2,B,maxes},Flatten[Map[({p1,p2}={#[[1]],Rest[#]};
B=Rest/@Select[elbows,First[#]<p1&];
B=dominationKernel[ClearNegsAndDeleteZeroVector[#-p2&/@B]];
If[B=={},{},maxes=Prepend[Max/@Transpose[B],0];
If[(PadLeft[maxes+#,Length[A]-1]+backdata).Rest[A]>currfar,(Prepend[#+p2,p1]&/@FarthestCornerSub[A,B,currfar,PadLeft[#,Length[A]-1]+backdata,opts]),skipct++;{}]])&,Select[elbows,First[#]>0&]],1]]/;Length[elbows[[1]]]>2;

FarthestCornerSub[_,elbows_,___]:=(Partition[Take[Flatten[Reverse/@Reverse[elbows]],{2,-2}],2])/;Length[elbows[[1]]]==2;

Examples

We show several examples to get some idea of time needed for these computations.

vals1={200,230,528,863,905,1355,1725};
vals2={13557,20002,52831,86312,90756};
vals3={18543816,27129592,43226644,78522677};
vals4={11615,27638,32124,48384,53542,56230};
vals5={10^10,18543816066,27129592681,43226644830,78522678316};
vals6={10000000000,35550333799,42807347468,67932625953,75136205898,79022523654};
vals7={10000000000,35550333799,42807347468,67932625953,75136205898,79022523654,21213456782};

Examples, continued

Timing[elbows1=elbows[vals1]]
Timing[f1=FarthestCorner[vals1,elbows1]]

Timing[elbows2=elbows[vals2];
f2=FarthestCorner[vals2,elbows2]]

Timing[elbows3=elbows[vals3];
f3=FarthestCorner[vals3,elbows3]]

Timing[elbows4=elbows[vals4];
f4=FarthestCorner[vals4,elbows4]]

Examples, continued

Timing[elbows5=elbows[vals5]]
Timing[f5=FarthestCorner[vals5,elbows5]]

Timing[elbows6=elbows[vals6];]
Timing[f6=FarthestCorner[vals6,elbows6]]

Timing[elbows7=elbows[vals7];]
Timing[f7=FarthestCorner[vals7,elbows7]]

To the best of my knowledge, the Frobenius numbers in examples 5 through 7 could not be done with methods available as of one year ago.

Summary

We have seen how toric Gröbner bases can be used to find the elbows of the Fundamental Domain for the Frobenius number problem. Unlike most applications of Gröbner bases involving finite varieties, toric and otherwise, what is paramount in this case is the actual staircase for the "normal set" with respect to a certain monomial order. The fact that Gröbner bases can be used at all to find Frobenius numbers, and the fact that we use them to gain staircase information, both seem to be of interest.

It is also of interest that this approach can handle Frobenius number problems that were not amenable to ANY method know as recently as a year or two ago.

The more efficient methods in our work in preparation heavily involve integer linear programming. Much of this has been incorporated into the Mathematica kernel, in the functions Reduce, FindInstance, and Minimize/Maximize.

Open questions:
(1) Would dedicated toric basis code do better? Probably so. How much better?
(2) Are there better ways to do the term ordering so that we can still get adequate information about the fundamental domain, but faster?

References

K. Aardal, C. A. J. Hurkens, and A. K. Lenstra.  Solving a system of linear diophantine equations with lower and upper bounds on the variables. Mathematics of Operations Research 25:427-442, 2000.

K. Aardal and A. K. Lenstra. Hard equality constrained knapsacks. Proceedings of the 9th Conference on Integer Programming and Combinatorial Optimization (IPCO 2002), W. J. Cook and A. S. Schulz, eds. Lecture Notes in Computer Science 233, 350-366. Springer-Verlag, 2002.

D. Beihoffer, J. Hendry, A. Nijenhuis, and S. Wagon. Faster algorithms for Frobenius numbers. The Electronic Journal of Combinatorics 12. 2005.

J. L. Bentley, K. L. Clarkson, and D. B. Levine. Fast linear .expected-time algorithms for computing maxima and convex hulls. Algorithmica 9(2): 168-183. 1993.

P. Conti and C. Traverso. Gröbner bases and integer programming. Proceedings of AAECC-9. Springer-Verlag LNCS 539 130-139. 1991.

J. L. Davison. On the linear Diophantine problem of Frobenius. J. Number Theory 48 353-363. 1994.

D. Einstein, D. Lichtblau, A. Strzebonski, and S. Wagon. Frobenius numbers by lattice enumeration. In preparation. 2005.

H. Greenberg. Solution to a linear Diophantine equation for nonnegative integers. J. Algorithms 9 353-363. 1988.

M. Keith. (Web page discussing Keith numbers)
http://users.aol.com/s6sj7gt/mikekeit.htm

D. Lichtblau. Solving knapsack and related problems. Manuscript.

L. Pottier. Gröbner bases of toric ideals. INRIA Rapport de recherche 2224. 1994.