Face Numbers of Random Polytopes

In order to test the perfromance of their algorithms in polyhedral geometry people often employ random constructions. While this is a good idea, combinatorially such polytopes are very restricted.

This will be demonstrated in this notebook by looking at random 6-polytopes.

Relevant publications on the subject include:

  • Buchta, Müller & Tichy: Stochastical approximation of convex bodies, Math. Ann. 271 (1985)
  • Borgwardt: The simplex method. A probabilistic analysis, Springer (1987).

Author: Michael Joswig

On the software side this showcases the use of polymake from within Oscar (via Polymake.jl).

In [1]:
using Oscar
const pm = Polymake
 -----    -----    -----      -      -----   
|     |  |     |  |     |    | |    |     |  
|     |  |        |         |   |   |     |  
|     |   -----   |        |     |  |-----   
|     |        |  |        |-----|  |   |    
|     |  |     |  |     |  |     |  |    |   
 -----    -----    -----   -     -  -     -  

...combining (and extending) ANTIC, GAP, Polymake and Singular
Version 0.6.0 ... 
 ... which comes with absolutely no warranty whatsoever
Type: '?Oscar' for more information
(c) 2019-2021 by The Oscar Development Team
Out[1]:
Polymake

Produce a 6-dimensional polytope with 20 vertices, drawn uniformly at random on the unit sphere, and print its f-vector.

In [2]:
P = pm.polytope.rand_sphere(6,20)
P.F_VECTOR
Out[2]:
pm::Vector<pm::Integer>
20 167 641 1188 1041 347

Now for the g-vector of the same polytope, but this time it is converted to a Julia Array.

In [3]:
g = convert(Array{Int32},P.G_VECTOR)
Out[3]:
4-element Vector{Int32}:
  1
 13
 68
 71

Let's do a census of 100 random samples and plot the result. Since g_0=1 and g_1=#vertices-dimension-1 are constant, it suffices to plot g_2 versus g_3.

In [4]:
n_vertices=30
n_samples=100
g_vectors = Array{Int32}(undef,n_samples,2)

for i=1:n_samples
    RS = pm.polytope.rand_sphere(6,n_vertices)
    g = RS.G_VECTOR
    g_vectors[i,1] = g[3] # notice index shift as Julia counts from 1
    g_vectors[i,2] = g[4]
end

We use the Plots package for visual output.

In [5]:
using Plots
scatter(g_vectors[:,1], g_vectors[:,2], xlabel="g_2", ylabel="g_3", legend=false)
Out[5]:

Compare the empirically found g_3 values (as a function of g_2) with the true range from Stanley's g-theorem.

In [6]:
min_g2 = minimum(g_vectors[:,1])
max_g2 = maximum(g_vectors[:,1])
ub = [ Int(pm.polytope.pseudopower(g2,2)) for g2 in min_g2:max_g2 ]
Out[6]:
30-element Vector{Int64}:
 1005
 1014
 1024
 1035
 1047
 1060
 1074
 1089
 1105
 1122
 1140
 1141
 1143
    ⋮
 1176
 1185
 1195
 1206
 1218
 1231
 1245
 1260
 1276
 1293
 1311
 1330

Or even compare with the maximal g-vector from McMullen's upper bound theorem shows that the random polytopes have rather few faces.

In [7]:
pm.polytope.upper_bound_theorem(6,n_vertices).G_VECTOR
Out[7]:
pm::Vector<pm::Integer>
1 23 276 2300

The lesson to learn here is the following.

By just looking at this class of random polytopes, one cannot get an adequate intuition about the set of all 6-polytopes with 20 vertices, not even the simplicial ones.