Vertex Enumeration 
for Convex Polyhedra and Hyperplane Arrangements
- Mathematica Package  -
Version 0.41 Beta (November 14, 1992)

Komei Fukuda and Ichiro Mizukoshi
Graduate School of Systems Management
University of Tsukuba, Tokyo
3-29-1 Otsuka, Bunkyo-ku
Tokyo 112, Japan

+81-3-3942-6876
fukuda@gssm.otsuka.tsukuba.ac.jp
mizukosi@gssm.otsuka.tsukuba.ac.jp

The package "VertexEnumeration" contains Mathematica implementations of Avis-Fukuda algorithms for enumerating all vertices of a convex polytope given by a system of linear inequalities, and for enumerating all points (0-dimensional faces) of an arrangement of hyperplanes given similarly.

This implementation is  quite faithful with the original description of the original algorithms. [1].

The core of this implementation consists of two functions, for a given  s*d  matrix  m  and s-vector  b :  

   (1)  VertexEnumeration[m, b]  
             -  for computing all vertices of the polyhedron:
                 P={ x :  x >= 0,  m. x  <=  b } and
  
   (2)  PointEnumeration[m,b]
             -  for computing all points of the arrangements
                 of  hyperplanes  h_i ,  i= 1,2,...,d+s,   
                 where   
                      h_j  = { x : x_j= 0 }  for j=1,....,d,  and
                      h_(d+i) = { x : m[[i]].x = b_i}  for i=1,...,s.                
  
In particular, these implementions can deal , in principle,  with any rational or floating-point real input with unlimited sizes  s,  d.   Moreover, with some limitations, irrational input handled.

As a drawback,  the actual time and space complexity is higher than it was originally stated in [1].  This is partly due to the way Mathematica  is implemented, and also due to our intention to make the program simple.

Incorporating with the standard package Combinatorica   [2] (which is not standard for Version1.2 but available via anonymous FTP from cs.sunysb.edu), one can analyze graph structures associated with convex polytopes and arrangements.

[1] D. Avis and K. Fukuda, 'A pivoting algorithm for convex hulls and vertex enumeration of arrangements and polyhedra,'   Proceedings of the 7th ACM Symposium on Computational Geometry, North Conway, New Hampshire, 1991, 98-104.  The full paper is to appear in Discrete and Computational Geometry.

[2] S. S. Skiena, 'Implementing Discrete Mathematics: Combinatorics and Graph Theory with Mathematica,' Addison-Wesley, Reading MA, 1990.

<<VertexEnum.m;

Off[General::spell];  Off[General::spell1];

Convex Polyhedra and Linear Inequalities

3-dimensional convex polyhedra are geometrical objects which are easily understood by drawings.  Mathematica  can draw some well known polyhedra, such as Platonic Solids using the package "Polyhedra.m":

However,  we cannot draw general convex polyhedra.

In general, a convex polyhedron is defined as the set of solutions to a system of linear inequalities:

(1)                m. x <= b  

where   m   is a real  s*d  matrix and b   is a real s-vector.   3-dimesional polyhera can be represented by such systems (1) with  d=3.   For example, the dodecahedon above is given by a system with s=12 and d=3.

The following complex polyhedron is given by a system 
(1) with a randomly generated matrix  m.

In order to draw such a polyhedron, one must compute all vertices (exremal points) and all edges (adjacencies).  The vertex enumeration package, in fact, computes these for a system of linear inequalities

(2)                m. x <= b   and    x >=0

where a matrix  m  and  a vector b   are  given.   The enumeration of vertices for  a  system  of form (1)  can be easily reduced to that for a system of form (2).

Very importantly, many fundamental problems can be reduced to the vertex enumeration.   Some of them are the computations of
    (a) the convex hull of s-points in d-dimensional space,
    (b) d-dimensional Voronoi diagrams and Delaunay
          triangulations,
    (c) optimal solutions to concave programming
         problems.
We shall explain the use of this package through some of these applications.

Enumeration of Vertices

3-dimensional Convex Polyhedra

Let us start with a simple example.    We shall generate a small random matrix  m  of  size * dim.

size=3;    dim=3;

m=Table[Table[Random[Integer,{1,4}],{dim}],{size}]

Let b be the vector of row-sums.

b=Table[Sum[m[[i,j]],{j,dim}],{i,size}]

Now we try to identify the shape of the convex polyhedron by using the vertex enumeration package.  (It will take about 1 minute by SE/30.)

vlist=VertexEnumeration[m,b]

The first part vlist[[1]]   of the output is a list of all vertices, and  the second part vlist[[2]]  is the list of "active-variable" sets, indicating which inequalities are satisfied with equalities at each vertex.   More precisely, for each vertex  v, the active variable set  is  Active(v)={j: h_j(v)=0} , where
 
              h_j  = { x : x_j= 0 }  for j=1,....,d,  and
              h_(d+i) = { x : m[[i]].x = b_i}  for i=1,...,s.

vertices=vlist[[1]];  activesets=vlist[[2]];

The function Polyhedron3D simply generates 3D graphics primitives for drawing the polyhedron.

polytope=Graphics3D[Polyhedron3D[vertices,activesets]];

Show[polytope, Boxed->True, ViewPoint->{2,3,1}];

The polyderon is determined by several (=size+dim) halfspaces, and bounded by their boundaries (i.e., hyperplanes).  By using PointEnumeration function, we can also compute all points of intersections of hyperplanes.  (The following takes about 2 minutes on SE/30.)

plist=PointEnumeration[m,b]

Then we can draw these points together with the convex polytope to see how the hyperplanes intersect outside the polytope.

pointplot=Graphics3D[Join[{PointSize[0.01]},
            Point /@ plist[[1]]]];

Show[{pointplot,polytope}, 
     Boxed->True, ViewPoint->{2,3,1}];

One can visualize the partition of the space by the hyperplanes if we draw the lines (1-dimensional flats)  which are intersections of the hyperlanes.

lines=LinesOfArrangement[plist[[2]]];

plines=(plist[[1]][[#]]&)/@ lines[[1]];

skeleton=Graphics3D[
     Join[{Thickness[ 0.002]},
           Line /@ plines,
          {PointSize[0.015]},
           Point /@ plist[[1]]]];

Show[{skeleton,polytope},
     Boxed->False,ViewPoint->{1,3,2}]

(* 
Do[Show[{skeleton,polytope},
   Boxed->False,
   ViewPoint->{2Cos[theta],2Sin[theta],1}],
   {theta,0,2Pi,Pi/6}]
*)

General Polyhedra and Graph Structures

Now we deal with a polyhedron with dimension higher than three.   We need nothing really special for doing this but perhaps more time.  

We shall generate again a small random matrix  m  of size * dim.

size=3;    dim=5;

m=Table[Table[Random[Integer,{0,3}],{dim}],{size}]

Let b  be the vector of all 10's.

b=Table[10, {size}];

Now we try to compute the vertices of the convex polyhedron by the vertex enumeration package.  (The following takes roughly 3-5 minutes on Macintosh SE/30 computer.)

{time,vlist}=Timing[VertexEnumeration[m,b]]

vertices=vlist[[1]];  activesets=vlist[[2]];

The following function computes the edges  of the polyhedron, i.e., the adjacent pairs of vertices.

edges=EdgesOfPolyhedron[activesets]

Another important structure we can extract is faces. The Face function gives the list of faces:

faces=Faces[activesets]

The j-th set faces[[j]]   is the set of vertices lying on the hyperplaneh_j,   where
 
         h_j  = { x : x_j= 0 }  for j=1,....,dim,  and
         h_(dim+i) = { x : m[[i]].x = b_i}  for i=1,...,size.

Each one is a face but may not be a facet.  The facets are the maximal faces.

facets=MaximalSets[faces]

If one wants to compute the adjacency among the facets,  simply use the function  EdgesOfPolyhedron,  again.

coedges=EdgesOfPolyhedron[facets]

With help of the combinatorica package,  you can draw and analyze the graph structure of the polyhedron.

(* The next drawing requires  Combinatorica.m  package. *)

<<Combinatorica.m;    (* Optional *)

ShowLabeledGraph[gpoly=FromUnorderedPairs[edges]]

ShowLabeledGraph[gpolar=FromUnorderedPairs[coedges]]

Voronoi Diagrams

2-dimensional Diagrams

Voronoi diagram is a useful notion for subdividing a plane into polygonal regions.   Say, for a given  n  locations of fire stations in some particular state (see, first figure below), Voronoi diagram shows which part of the state should be taken care of by which station.   More precisely, Voronoi diagram is the partition of the state into  n  sub-regions, such that each sub-region contains exactly one fire station and every point in a sub-region is closer to its fire station than to the other ones (see, second figure below).

Let us make a small example by generating few random points in the plane.

dim=2;  size=5;  min=1;  max=30;  width=max-min;

datapoints=
  Table[Table[Random[Integer,
        {min+width/10, max-width/10}],{dim}],{size}]

Let us draw these points using 3D graphics.  (2D graphics is sufficient for this purpose but 3D graphics will be useful later.)

points=Map[Append[#,0]&,datapoints];

dataplot=Graphics3D[Join[{PointSize[0.01]},
          Point /@ points, 
         {Polygon[{{min,min,0},{min,max,0},
                  {max,max,0},{max,min,0}}]}]];

Show[dataplot, Boxed->False, 
     ViewPoint->{0,0,3}, BoxRatios->{1,1,0}];

In order to compute the Voronoi diagram for the given points, we compute the vertices of 3-dimensional convex polyhedron and then project its boundary to a plane.   The convex polyhedron corresponds to a piecewise-linear approximation of a parabola:  z=x^2+y^2.  Each piece is associated with some datapoint  {x_i, y_i}  and is a part of the hyperplane supporting the parabola at this point.   This hyperplane is represented by  z=2*x + 2*y + (x_i^2 + y_i^2).  The convex polyhedron is then represented by the linear inequalities:
z >= 2*x + 2*y + (x_i^2 + y_i^2), i=1,...,size.

(*  The above is generated by 
parabola=Plot3D[x^2+y^2,{x,0,max},{y,0,max},
    BoxRatios->{1,1,1},  ViewPoint->{1,-2,1}]    
*)

We make the input matrix  m  and vector  b  for the vertex enumeration package for the convex polyhedron.  Here notice that we will compute this within the square 
[{x, min, max}, {y, min, max}].   Also, the parabola will be lifted so that the approximation won't touch the bottom plane z=0.

m1=Map[Append[2 #,-1]&,datapoints];

b0=Map[Sum[#[[j]]^2,{j,dim}]&,datapoints];

blift=2*Min[b0];   b1=Map[(#-blift)&,b0];

m=Join[m1,{{-1,0,0},{0,-1,0},{1,0,0},{0,1,0}}];

b=Join[b1,{-min,-min,max,max}];

Please notice that the main computation below takes 3-4 minutes on SE/30.

{time,vlist}=Timing[VertexEnumeration[m,b]]

We can see the voronoi diagram by projecting the polyhedral surface onto x-y plane.   Of course, this can be simulated by looking the surface from the sky.

poly=Graphics3D[
       Drop[Polyhedron3D[vlist[[1]],vlist[[2]]],-4]];

Show[{poly,dataplot},
     ViewPoint->{0,0,3}, Shading->True,
     Boxed->False, BoxRatios->{1,1,0}];

In order to compute the graph structure of the Voronoi diagram, one must do a little more.

VoronoiAdjacency=EdgesOfPolyhedron[vlist[[2]]];

projectedvertices=Map[{#[[1]],#[[2]],0}&,vlist[[1]]];

VoronoiEdges=
   (projectedvertices[[#]]&) /@ VoronoiAdjacency;

VoronoiDiagram=Graphics3D[
    Join[{Thickness[ 0.0025],GrayLevel[0]},
           Line /@ VoronoiEdges,
          {PointSize[0.005],GrayLevel[0]},
           Point /@ projectedvertices]];

We can now draw the Voronoi diagram, the datapoints and the polyhedral surface at the same time.

Show[{poly,dataplot,VoronoiDiagram},
     ViewPoint->{0,2,2},  Shading->True,
     BoxRatios->{1,1,1}]

The Delaunay triangulation comes hand-in-hand with the Voronoi diagam.  What we have to do is to compute the adjacency among the Voronoi regions (facets of the associated polyhedron.)

facets=Faces[vlist[[2]]];

coedges=EdgesOfPolyhedron[facets];

DelaunayAdjacency=Select[coedges,
    (dim+1<#[[1]]<=size+dim+1 && 
     dim+1<#[[2]]<=size+dim+1)&]

DelaunayEdges=(points[[#-dim-1]]&) /@ DelaunayAdjacency;

Delaunay=Graphics3D[
    Join[{Thickness[ 0.0055]},Line /@ DelaunayEdges,
         {PointSize[0.02],GrayLevel[0]}, 
         Point /@ points]];

Show[{VoronoiDiagram,Delaunay},
     ViewPoint->{0,0,3},  Boxed->False,  
     BoxRatios->{1,1,0}];

General Diagrams

General (d-dimensional) Voronoi diagrams (tasselations) can be computed in the exactly same manner as for the plane, through computing the structure of a (d+1)-dimensional convex polyhedron.

dim=4; size=5; min=5;  max=20;  width=max-min;

datapoints=
  Table[Table[Random[Integer,
        {min+width/10,max-width/10}], {dim}],
       {size}]

m1=Map[Append[2 #,-1]&,datapoints];

b0=Map[Sum[#[[j]]^2,{j,dim}]&,datapoints];

blift=2*Min[b0];   b1=Map[(#-blift)&,b0];

m=m1;  b=b1;

We can compute the vertices of  the polyhedron whose boundary represents the Voronoi partitions.   Notice that the computation will be restricted to the nonnegative orthant and thus the result may not contain the  complete structure of the partition.

The next one takes 3-4 minutes on  SE/30.

{time,vlist}=Timing[VertexEnumeration[m,b]]

facets=Faces[vlist[[2]]];

coedges=EdgesOfPolyhedron[facets];

DelaunayAdjacency=Select[coedges,
    (dim+1<#[[1]]<=size+dim+1 && 
     dim+1<#[[2]]<=size+dim+1)&];

DelaunayEdges=(#-dim-1)& /@ DelaunayAdjacency;

We can draw the Delaunay structure (the adjacency among the Voronoi regions) by using Combinatorica package.

(*  << Combinatorica.m;    *)

ShowLabeledGraph[
   DelaunayGraph=FromUnorderedPairs[DelaunayEdges]];

Convex Hull

3-dimensional Convex Hull

Computing the convex hull of a given set of points can be indirectly done by the vertex enumeration.   In this section, we see how this can be done through a 3-dimensional example.

First we generate 10 points in 3-space.

size=10;      dim=3;

datapoints=
   Table[Table[Random[Integer,{-2,2}],{dim}],{size}]

datapoints=Union[datapoints];
size=Length[datapoints]

dataplot=Graphics3D[Join[{PointSize[0.01]},
           Point /@ datapoints]];

Show[dataplot, ViewPoint->{2,3,1}, BoxRatios->{1,1,1}];

The convex hull of given points is the smallest convex region (convex polyhedron) containing the points.   

The theory relating the convex hull and the vertex enumeration is the polarity.   Let  P  be a convex polyhedron containing the origin in its interior.   Let  P*  be the convex set defined by  {y :  x.y  <= 1  for all  x  in  P}.  Then P* is also a convex polyhedron whose face lattice is anti-isomorphic to that of  P.   Therefore, there is a one-to-one correspondance between the vertices of  P*  and the facets of  P  preserving the adjacency.

Thus, what we must do first is to translate the given points so that the origin will be contained in the convex hull.

center=Sum[datapoints[[i]],{i,size}]/size

translatedpoints=Map[(#-center)&,datapoints];

m=translatedpoints;

c={100,100,100};

b=Table[1,{size}];

******  Caution:   Long computation  ******
The next one takes 4-5 minutes on  SE/30.

flist=VertexEnumeration[m,b+m.c];

flist1=Map[(#-c)&,flist[[1]]]

We can easily compute which datapoints form a facet of the convex hull.

facets=Map[(#-Table[dim,{Length[#]}])&,flist[[2]]]

activesets=
Map[Flatten[Position[Chop[flist1.#-1],0]+3]&,translatedpoints];

extsets=MaximalSets[activesets];

maxpos=Flatten[Map[Position[activesets,#]&,extsets]];

extdata=datapoints[[maxpos]];

convexhull=Graphics3D[Polyhedron3D[extdata,extsets]];

Show[convexhull, Boxed->True, ViewPoint->{2,3,1}];

General Convex Hull

Here, we just note that the convex hull of points in general d-space can be computed through the same operation (polarization) and the vertex enumeration.

Options for Enumeration

Monitoring the Computation

Sometimes, you may want to monitor how the computation is being done.   For this purpose, one can use  MonitoringFile  option, which allows you to output the list of nonbasic varibles as the computation proceeds.   The output can be directed to a file you specify, or to CRT if the filename is "stdout".

m={{1,1,1},{0,1,0},{0,0,1}};

b={1,1,1};

vlist=VertexEnumeration[m,b,MonitoringFile->"stdout"]

Output All Basic Solutions

Instead of outputing all vertices without repetition, one can list all basic solutions with possible repetitions of vertices.   For this purpose one can use the SearchTree option.   If this option is set  True (default is False), VertexEnumeration outputs all basic feasible solutions together with the Search Tree structure which the Avis-Fukuda algorithm actually traces.

m={{1,1,1},{0,1,0},{0,0,1}};

b={1,1,1};

vlist=VertexEnumeration[m,b,SearchTree->True,
             MonitoringFile->"stdout"]

Here is a way to draw the Avis-Fukuda search tree  by using a nice embedding function of Combinatorica.

STree=vlist[[3]];

ShowLabeledGraph[RankedEmbedding[FromUnorderedPairs[STree],{1}]]

Dealing with Irrational Data in Mathematica 1.2

Mathematica  can handle irrational numbers, but we must be careful because some older versions (Version 1.2 or earlier) do not simplify some trivial expressions like

ir = 5* Sqrt[5] - 5^(3/2)

which is obviously zero.   This may cause some serious problems if you want to use irrational numbers in input data.    The separate small program "Power.m" can be used with VertexEnumeration and PointEnumeration to prevent such situations with irrational numbers of rational power type  p^q   where both p and q are rational numbers.

m={{ir +1, ir + 5^(1/3), 1}};   b={1};

<<Power.m;    (* this is for Mathematica Version 1.2 *)

vlist=VertexEnumeration[m,b]

Supplementary Functions

There are some useful functions defined in the VertexEnum package.   We shall introduce some of those below.

Linear Programming

The problem of maximizing a linear function over a convex polyhedron is known to be Linear programming.   A standard form of linear programming problem (LP) is  maximize  c.x  subject to m.x <=b and x>=0, where  m  is a given  s x d matrix,  c and b are  given d-vector and s-vector, respectively.    The function CrissCrossSolve solves this LP.

?CrissCrossSolve

m={{4, -3, 3}, {-2, 1, -3}};

c={1,-2,-4};  b={5,-3};

optimum=CrissCrossSolve[c,m,b]

CrissCrossSolve always finds a vertex solution if the LP has an optimal solution.  You can verify this by using VetexEnumeration.

VertexEnumeration[m,b]

Mathematica has also LinearProgramming function which solves LP.  One has to be little careful because some earliear versions of LinearProgramming function (Version 2.0 or earlier) do not work correctly.   Also this function solves a different form of LP, which is dual to our formulation.

?LinearProgramming

In order to solve the same LP as above, we must use the negatives of the input arguments.

LinearProgramming[-c,-m,-b]

While LinearProgramming function often runs faster for rational inputs, it does not allow any irrational inputs.   It is an advantage of CrissCrossSolve that irrational inputs are allowed.  (The next one takes about 1 minute on SE/30.)

m={{1, Pi, E}, {E, 1, Pi}}; 
b={Log[5],Sqrt[3]};  c={E,Pi,Sqrt[5]};
opt=CrissCrossSolve[c,m,b]

For LPs with no optimal solutions, CrissCrossSolve outputs the following warnings.

CrissCrossSolve[{1,1},{{2,3},{-2,-3}},{2,-4}]

CrissCrossSolve[{1,1},{{-2,3}},{-4}]

It should be noted that the dual infeasibility implies that the original problem is either unbounded or infeasible, i.e., cannot have an optimal solution.

List of Exported Functions and Options

?VertexEnum`*

(* This is the end of VertexEnum-notebook  *)