A Brief Introduction to Polytopes in Sage

Author: sarah-marie belcastro <smbelcas@toroidalsnark.net>

If you already know some convex geometry a la Grünbaum or Brøndsted, then you may have itched to get your hands dirty with some polytope calculations. Here is a mini-guide to doing just that.

Basics

First, let’s define a polytope as the convex hull of a set of points, i.e. given \(S\) we compute \(P={\rm conv}(S)\):

sage: P1 = Polyhedron(vertices = [[-5,2], [4,4], [3,0], [1,0], [2,-4], [-3,-1], [-5,-3]])
sage: P1
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 4 vertices

>>> from sage.all import *
>>> P1 = Polyhedron(vertices = [[-Integer(5),Integer(2)], [Integer(4),Integer(4)], [Integer(3),Integer(0)], [Integer(1),Integer(0)], [Integer(2),-Integer(4)], [-Integer(3),-Integer(1)], [-Integer(5),-Integer(3)]])
>>> P1
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 4 vertices

P1 = Polyhedron(vertices = [[-5,2], [4,4], [3,0], [1,0], [2,-4], [-3,-1], [-5,-3]])
P1

Notice that Sage tells you the dimension of the polytope as well as the dimension of the ambient space.

Of course, you want to know what this object looks like:

sage: P1.plot()
Graphics object consisting of 6 graphics primitives

>>> from sage.all import *
>>> P1.plot()
Graphics object consisting of 6 graphics primitives

P1.plot()

Even in only 2 dimensions, it’s a pain to figure out what the supporting hyperplanes are. Luckily Sage will take care of that for us.

sage: for q in P1.Hrepresentation():
....:    print(q)
An inequality (-4, 1) x + 12 >= 0
An inequality (1, 7) x + 26 >= 0
An inequality (1, 0) x + 5 >= 0
An inequality (2, -9) x + 28 >= 0

>>> from sage.all import *
>>> for q in P1.Hrepresentation():
...    print(q)
An inequality (-4, 1) x + 12 >= 0
An inequality (1, 7) x + 26 >= 0
An inequality (1, 0) x + 5 >= 0
An inequality (2, -9) x + 28 >= 0

for q in P1.Hrepresentation():
   print(q)

That notation is not immediately parseable, because seriously, those do not look like equations of lines (or of halfspaces, which is really what they are).

(-4, 1) x + 12 >= 0 really means \((-4, 1)\cdot\vec{x} + 12 \geq 0\).

So… if you want to define a polytope via inequalities, you have to translate each inequality into a vector. For example, \((-4, 1)\cdot\vec{x} + 12 \geq 0\) becomes (12, -4, 1).

sage: altP1 = Polyhedron(ieqs=[(12, -4, 1), (26, 1, 7),(5,1,0), (28, 2, -9)])
sage: altP1.plot()
Graphics object consisting of 6 graphics primitives

>>> from sage.all import *
>>> altP1 = Polyhedron(ieqs=[(Integer(12), -Integer(4), Integer(1)), (Integer(26), Integer(1), Integer(7)),(Integer(5),Integer(1),Integer(0)), (Integer(28), Integer(2), -Integer(9))])
>>> altP1.plot()
Graphics object consisting of 6 graphics primitives

altP1 = Polyhedron(ieqs=[(12, -4, 1), (26, 1, 7),(5,1,0), (28, 2, -9)])
altP1.plot()

Other information you might want to pull out of Sage about a polytope is the vertex list, which can be done in two ways:

sage: for q in P1.Vrepresentation():
....:    print(q)
A vertex at (-5, -3)
A vertex at (-5, 2)
A vertex at (4, 4)
A vertex at (2, -4)

>>> from sage.all import *
>>> for q in P1.Vrepresentation():
...    print(q)
A vertex at (-5, -3)
A vertex at (-5, 2)
A vertex at (4, 4)
A vertex at (2, -4)

for q in P1.Vrepresentation():
   print(q)
sage: P1.vertices()
(A vertex at (-5, -3), A vertex at (-5, 2), A vertex at (4, 4), A vertex at (2, -4))

>>> from sage.all import *
>>> P1.vertices()
(A vertex at (-5, -3), A vertex at (-5, 2), A vertex at (4, 4), A vertex at (2, -4))

P1.vertices()

Polar duals

Surely you want to compute the polar dual:

sage: P1dual = P1.polar()
sage: P1dual
A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 4 vertices

>>> from sage.all import *
>>> P1dual = P1.polar()
>>> P1dual
A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 4 vertices

P1dual = P1.polar()
P1dual

Check it out---we started with an integer-lattice polytope and dualized to a rational-lattice polytope. Let’s look at that.

sage: P1dual.plot()
Graphics object consisting of 6 graphics primitives

>>> from sage.all import *
>>> P1dual.plot()
Graphics object consisting of 6 graphics primitives

P1dual.plot()
sage: P1.plot() + P1dual.plot()
Graphics object consisting of 12 graphics primitives

>>> from sage.all import *
>>> P1.plot() + P1dual.plot()
Graphics object consisting of 12 graphics primitives

P1.plot() + P1dual.plot()

Oh, yeah, unless the polytope is unit-sphere-sized, the dual will be a very different size. Let’s rescale.

sage: ((1/4)*P1).plot() + (4*P1dual).plot()
Graphics object consisting of 12 graphics primitives

>>> from sage.all import *
>>> ((Integer(1)/Integer(4))*P1).plot() + (Integer(4)*P1dual).plot()
Graphics object consisting of 12 graphics primitives

((1/4)*P1).plot() + (4*P1dual).plot()

If you think that looks a little bit shady, you’re correct. Here is an example that makes the issue a bit clearer.

sage: P2 = Polyhedron(vertices = [[-5,0], [-1,1], [-2,0], [1,0], [-2,-1], [-3,-1], [-5,-1]])
sage: P2
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 5 vertices
sage: P2dual = P2.polar(); P2dual
A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 5 vertices
sage: P2.plot() + P2dual.plot()
Graphics object consisting of 14 graphics primitives

>>> from sage.all import *
>>> P2 = Polyhedron(vertices = [[-Integer(5),Integer(0)], [-Integer(1),Integer(1)], [-Integer(2),Integer(0)], [Integer(1),Integer(0)], [-Integer(2),-Integer(1)], [-Integer(3),-Integer(1)], [-Integer(5),-Integer(1)]])
>>> P2
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 5 vertices
>>> P2dual = P2.polar(); P2dual
A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 5 vertices
>>> P2.plot() + P2dual.plot()
Graphics object consisting of 14 graphics primitives

P2 = Polyhedron(vertices = [[-5,0], [-1,1], [-2,0], [1,0], [-2,-1], [-3,-1], [-5,-1]])
P2
P2dual = P2.polar(); P2dual
P2.plot() + P2dual.plot()

That is clearly not computing what we think of as the polar dual. But look at this…

sage: P2.plot() + (-1*P2dual).plot()
Graphics object consisting of 14 graphics primitives

>>> from sage.all import *
>>> P2.plot() + (-Integer(1)*P2dual).plot()
Graphics object consisting of 14 graphics primitives

P2.plot() + (-1*P2dual).plot()

Here is what’s going on.

If a polytope P is in \(\ZZ\), then…

  1. …the dual is inverted in some way, which is vertically for polygons.

  2. …the dual is taken of P itself.

  3. …if the origin is not in P, then an error is returned.

However, if a polytope is not in \(\ZZ\), for example if it’s in \(\QQ\) or RDF, then…

(1’) …the dual is not inverted.

(2’) …the dual is taken of P-translated-so-barycenter-is-at-origin.

Keep all of this in mind as you take polar duals.

Polytope Constructions

Minkowski sums! Now with two syntaxes!

sage: P1+P2
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 8 vertices

>>> from sage.all import *
>>> P1+P2
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 8 vertices

P1+P2
sage: P1.minkowski_sum(P2)
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 8 vertices

>>> from sage.all import *
>>> P1.minkowski_sum(P2)
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 8 vertices

P1.minkowski_sum(P2)

Okay, fine. We should have some 3-dimensional examples, at least. (Note that in order to display polytopes effectively you’ll need visualization software such as Javaview and Jmol installed.)

sage: P3 = Polyhedron(vertices=[(0,0,0), (0,0,1/2), (0,1/2,0), (1/2,0,0), (3/4,1/5,3/2)]); P3
A 3-dimensional polyhedron in QQ^3 defined as the convex hull of 5 vertices
sage: P4 = Polyhedron(vertices=[(-1,1,0),(1,1,0),(-1,0,1), (1,0,1),(0,-1,1),(0,1,1)]); P4
A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 6 vertices
sage: P3.plot() + P4.plot()
Graphics3d Object

>>> from sage.all import *
>>> P3 = Polyhedron(vertices=[(Integer(0),Integer(0),Integer(0)), (Integer(0),Integer(0),Integer(1)/Integer(2)), (Integer(0),Integer(1)/Integer(2),Integer(0)), (Integer(1)/Integer(2),Integer(0),Integer(0)), (Integer(3)/Integer(4),Integer(1)/Integer(5),Integer(3)/Integer(2))]); P3
A 3-dimensional polyhedron in QQ^3 defined as the convex hull of 5 vertices
>>> P4 = Polyhedron(vertices=[(-Integer(1),Integer(1),Integer(0)),(Integer(1),Integer(1),Integer(0)),(-Integer(1),Integer(0),Integer(1)), (Integer(1),Integer(0),Integer(1)),(Integer(0),-Integer(1),Integer(1)),(Integer(0),Integer(1),Integer(1))]); P4
A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 6 vertices
>>> P3.plot() + P4.plot()
Graphics3d Object

P3 = Polyhedron(vertices=[(0,0,0), (0,0,1/2), (0,1/2,0), (1/2,0,0), (3/4,1/5,3/2)]); P3
P4 = Polyhedron(vertices=[(-1,1,0),(1,1,0),(-1,0,1), (1,0,1),(0,-1,1),(0,1,1)]); P4
P3.plot() + P4.plot()
sage: (P3+P4).plot()
Graphics3d Object

>>> from sage.all import *
>>> (P3+P4).plot()
Graphics3d Object

(P3+P4).plot()

We can also find the intersection of two polytopes… and this too has two syntaxes!

sage: int12 = P1.intersection(P2*.5); int12.plot()
Graphics object consisting of 7 graphics primitives

>>> from sage.all import *
>>> int12 = P1.intersection(P2*RealNumber('.5')); int12.plot()
Graphics object consisting of 7 graphics primitives

int12 = P1.intersection(P2*.5); int12.plot()
sage: int34 = P3 & P4; int34.plot()
Graphics3d Object

>>> from sage.all import *
>>> int34 = P3 & P4; int34.plot()
Graphics3d Object

int34 = P3 & P4; int34.plot()

Should one wish to translate, one can.

sage: transP2 = P2.translation([2,1])
sage: P2.plot() + transP2.plot()
Graphics object consisting of 14 graphics primitives

>>> from sage.all import *
>>> transP2 = P2.translation([Integer(2),Integer(1)])
>>> P2.plot() + transP2.plot()
Graphics object consisting of 14 graphics primitives

transP2 = P2.translation([2,1])
P2.plot() + transP2.plot()

Then of course we can take prisms, pyramids, and bipyramids of polytopes…

sage: P2.prism().plot()
Graphics3d Object

>>> from sage.all import *
>>> P2.prism().plot()
Graphics3d Object

P2.prism().plot()
sage: P1.pyramid().plot()
Graphics3d Object

>>> from sage.all import *
>>> P1.pyramid().plot()
Graphics3d Object

P1.pyramid().plot()
sage: P2dual.bipyramid().plot()
Graphics3d Object

>>> from sage.all import *
>>> P2dual.bipyramid().plot()
Graphics3d Object

P2dual.bipyramid().plot()

Okay, fine. Yes, Sage has some kinds of polytopes built in. If you type polytopes. and then press TAB after the period, you’ll get a list of pre-built polytopes.

sage: P5 = polytopes.hypercube(5)
sage: P6 = polytopes.cross_polytope(3)
sage: P7 = polytopes.simplex(7)

>>> from sage.all import *
>>> P5 = polytopes.hypercube(Integer(5))
>>> P6 = polytopes.cross_polytope(Integer(3))
>>> P7 = polytopes.simplex(Integer(7))

P5 = polytopes.hypercube(5)
P6 = polytopes.cross_polytope(3)
P7 = polytopes.simplex(7)

Let’s look at a 4-dimensional polytope.

sage: P8 = polytopes.hypercube(4)
sage: P8.plot()
Graphics3d Object

>>> from sage.all import *
>>> P8 = polytopes.hypercube(Integer(4))
>>> P8.plot()
Graphics3d Object

P8 = polytopes.hypercube(4)
P8.plot()

We can see it from a different perspective:

sage: P8.schlegel_projection(position=1/2).plot()
Graphics3d Object

>>> from sage.all import *
>>> P8.schlegel_projection(position=Integer(1)/Integer(2)).plot()
Graphics3d Object

P8.schlegel_projection(position=1/2).plot()

Queries to polytopes

Once you’ve constructed some polytope, you can ask Sage questions about it.

sage: P1.contains([1,0])
True

>>> from sage.all import *
>>> P1.contains([Integer(1),Integer(0)])
True

P1.contains([1,0])
sage: P1.interior_contains([3,0])
False

>>> from sage.all import *
>>> P1.interior_contains([Integer(3),Integer(0)])
False

P1.interior_contains([3,0])
sage: P3.contains([1,0,0])
False

>>> from sage.all import *
>>> P3.contains([Integer(1),Integer(0),Integer(0)])
False

P3.contains([1,0,0])

Face information can be useful.

sage: int34.f_vector()
(1, 8, 12, 6, 1)

>>> from sage.all import *
>>> int34.f_vector()
(1, 8, 12, 6, 1)

int34.f_vector()

Well, geometric information might be more helpful… Here we are told which of the vertices form each 2-face:

sage: [f.ambient_V_indices() for f in int34.faces(2)]
[(2, 6, 7), (0, 1, 3, 5), (1, 3, 4), (0, 5, 6, 7), (0, 1, 2, 4, 6), (2, 3, 4, 5, 7)]

>>> from sage.all import *
>>> [f.ambient_V_indices() for f in int34.faces(Integer(2))]
[(2, 6, 7), (0, 1, 3, 5), (1, 3, 4), (0, 5, 6, 7), (0, 1, 2, 4, 6), (2, 3, 4, 5, 7)]

[f.ambient_V_indices() for f in int34.faces(2)]

Yeah, that isn’t so useful as it is. Let’s figure out the vertex and hyperplane representations of the first face in the list.

sage: first2faceofint34 = int34.faces(2)[0]
sage: first2faceofint34.ambient_Hrepresentation(); first2faceofint34.vertices()
(An inequality (0, 0, -1) x + 1 >= 0,)
(A vertex at (2/3, 2/15, 1), A vertex at (3/8, 1/10, 1), A vertex at (1/2, 3/10, 1))

>>> from sage.all import *
>>> first2faceofint34 = int34.faces(Integer(2))[Integer(0)]
>>> first2faceofint34.ambient_Hrepresentation(); first2faceofint34.vertices()
(An inequality (0, 0, -1) x + 1 >= 0,)
(A vertex at (2/3, 2/15, 1), A vertex at (3/8, 1/10, 1), A vertex at (1/2, 3/10, 1))

first2faceofint34 = int34.faces(2)[0]
first2faceofint34.ambient_Hrepresentation(); first2faceofint34.vertices()

If you want more… Base class for polyhedra: Miscellaneous methods is the first place you want to go.