Parallel Enumeration of Triangulations
Charles Jordan∗
Graduate School of Information Science and Technology
Hokkaido University
N-14 W-9, 060-0814 Sapporo, Japanskip@ist.hokudai.ac.jp

Michael Joswig†

Lars Kastner‡

Institut für Mathematik
TU Berlin
Str. des 17. Juni 136, 10623 Berlin, Germany{kastner,joswig}@math.tu-berlin.deSubmitted: Sep 14, 2017; Accepted: May 22, 2018; Published: Jul 13, 2018c The authors. Released under the CC BY-ND license (International 4.0).

Abstract
We report on the implementation of an algorithm for computing the set of all
regular triangulations of finitely many points in Euclidean space. This algorithm,
which we call down-flip reverse search, can be restricted, e.g., to computing full
triangulations only; this case is particularly relevant for tropical geometry. Most
importantly, down-flip reverse search allows for massive parallelization, i.e., it scales
well even for many cores. Our implementation allows to compute the triangulationsof much larger point sets than before.
Mathematics Subject Classifications: 52B55 (68U05)

1Introduction

Triangulations are ubiquitous in combinatorics, optimization, algebra and other parts ofmathematics. For an overview about the range of applications we recommend the first∗Research by C. Jordan supported in part by JSPS Kakenhi Grants 16H02785 and 18K18027.
Research by M. Joswig is carried out in the framework of Matheon supported by Einstein Foundation
Berlin. Further partial support by Deutsche Forschungsgemeinschaft (SFB-TRR 109: “Discretization inGeometry and Dynamics” and SFB-TRR 195: “Symbolic Tools in Mathematics and their Application”) is
gratefully acknowledged.
‡
Research by L. Kastner is supported by Deutsche Forschungsgemeinschaft (SFB-TRR 195: “Symbolic
Tools in Mathematics and their Application”).
†

the electronic journal of combinatorics 25(3) (2018), #P3.6

1

chapter of the monograph [7] by De Loera, Rambau and Santos, and the rest of that book
is recommended for the foundations. The software which defines the state of the art isRambau’s TOPCOM [25]; see also [22]. An implementation of a different method is part ofGfan by Jensen [16]. While this comes with some advantages of its own, TOPCOM is superior
in most cases. In this article we describe an algorithm using different ideas but whichstill shares some of TOPCOM’s features. Most importantly we report on new computational
results obtained with our implementation MPTOPCOM; some of these are out of the reach of
the other software systems.
In order to explain the method it is best to start by looking into TOPCOM’s algorithm;
see also [7, §8.3]. As a first step, for a fixed point configuration, P , TOPCOM creates one
triangulation as a seed. Then TOPCOM generates all triangulations of P which can beobtained by local transformations, known as flips. The flips induce a graph structureon the set of triangulations of P , and TOPCOM makes a breadth first search starting at
the seed. In practical applications it is often important to restrict the attention to thosetriangulations of P which are regular, i.e., they are induced by a convex lifting function.
Regularity can be checked via solving a linear programming feasibility problem, and thisis supported by TOPCOM, e.g., via calling cddlib [10]. Any two regular triangulations of
P are connected by a sequence of flips. Checking for regularity on the way is costly but
is also beneficial since the regular triangulations of P have a particularly nice structure.
The reason is that they correspond to the vertices of a convex polytope, the secondary
polytope of P , which was introduced by Gel0 fand, Kapranov and Zelevinsky [12], andthe edges of the secondary polytope correspond to flips (but there are flips which do notarise from edges). The coordinates of a regular triangulation ∆, seen as a vertex of thesecondary polytope, form the GKZ-vector of ∆. Already in 2002 Imai et al. [13] proposedan algorithm which exploits the GKZ-vectors to tailor a reverse search scheme [2] forenumerating the regular triangulations of P , and we call their method down-flip reverse
search. Surprisingly, the down-flip reverse search algorithm apparently received little
attention, and some of its properties have been discovered independently by Pournin and
Liebling [24]. A key challenge in practice is that relevant point sets P exhibit a great dealof symmetry. Usually this symmetry is given in terms of generators of some finite group G
acting on P by permutations. The real task now is to compute the (regular) triangulationsof P up to symmetry, i.e., exactly one representative from each G-orbit. This reduction
is strictly necessary since the total number of triangulations is often too large. Imai etal. [13] also describe down-flip reverse search up to symmetry.
One of the main advantages of the reverse search algorithm is that it is output sensitive.
In particular, if one has bounds on the number of neighbors of a node in the search graph,one can derive an effective upper bound on the memory consumption of reverse search, see[2, Theorem 3.2]. The price one pays is recomputation of intermediate results, whenever anode is analyzed, but not visited. This obstacle can be partially overcome by using caches.
While TOPCOM’s limit often is dictated by the amount of memory available, MPTOPCOM’s
limit is given by the time one is willing to let the program run. This is where parallelization
comes in.
While the general algorithmic idea is known, a practical implementation of down-flip
the electronic journal of combinatorics 25(3) (2018), #P3.6

2reverse search requires one to overcome several challenges. That is the focus of thispaper. We will show that our implementation MPTOPCOM of down-flip reverse search, which
can enumerate hundreds of millions of triangulations of a given point set, is superior toother methods in practice. Two aspects are most important. First, for large input, theimplementation must be parallelized in order to be able to benefit from modern hardware.
This is where the reverse search scheme shows its full strength; mts [3] is a competitive
parallel implementation of the abstract reverse search method based on MPI [30]. Second, by
far the most frequent subtask is to compute a canonical representative for each G-orbit of agiven triangulation. This allows to distinguish two orbits by comparing their representatives.
A minor drawback of parallelized reverse search is that those representatives are recomputedby several workers, even for the same orbit and the same triangulation. This is because
we want to avoid any communication between the nodes in the computation tree since
this is what makes parallelized reverse search so successful. Our main new algorithmic
contribution is a procedure for computing these canonical representatives efficiently. Thisis a variation of the Schreier–Sims algorithm adapted to our special situation; cf. [29,§4.2]. Despite our highly optimized setup, depending on the input and the desired output,
the computation of the canonical representatives still may take more than 90% of the
combined total time. The benefit of our approach, however, is tremendous. It turns out
that, even if we restrict our parallel implementation to a single core, we are able to beat
TOPCOM by a factor of five or more on medium size input. More importantly, on large input
MPTOPCOM scales rather well even for more than a hundred cores. This means that we are
now able to compute triangulations of point sets an order of magnitude larger than before.
This paper is organized as follows. In Section 2 we collect the basic facts concerning
the reverse search method of Avis and Fukuda [2]. This is a powerful general scheme
for organizing a large enumeration via a rooted tree. The paradigmatical example is adual convex hull computation. Here the input is a system of linear inequalities, and theoutput are the vertices and the rays of the feasible region, which is a convex polyhedron.
A generic linear objective function induces an orientation on the edges of that polyhedron,and any pivoting strategy for the simplex method yields a tree, in which each edge isdirected toward the global optimum. The basic idea of down-flip reverse search is to mimic
this behavior on the secondary polytope of a point configuration. We close this section bydescribing the principles underlying the parallelization employed by mts [3]. Section 3 deals
with the basic notions concerning triangulations, and we formulate the down-flip reverse
search algorithm. Then, in Section 4 we briefly explain the general approach of Imai et al.
[13] to combine enumeration up to symmetry with down-flip reverse search. The core of
our paper is Section 5. The major challenge in MPTOPCOM’s implementation is to efficiently
find canonical representatives for each orbit of triangulations. To this end we introduce theconcept of switch tables and evaluations to facilitate that computation (cf. Definition 9).
This is developed in a way which specifically addresses the enumeration of triangulations,and yet we believe that it may also be useful in other circumstances. Our main theoreticalcontribution is a procedure for computing canonical representatives via switch tables (cf.
Algorithm 4) and its analysis. In Section 7 we report on experimental results. First we
consider a few standard examples, such as the 16 vertices of the 4-dimensional regular cube.
the electronic journal of combinatorics 25(3) (2018), #P3.6

3More interesting are our new results. This concerns, e.g., products of simplices; thesecomputations have been employed by Schröter for deriving results on coarsest matroid
subdivision of hypersimplices [28]. Further, for the first time we are able to enumerate
all regular and full triangulations of the 3-simplex with dilation factor of three. This is aconfiguration of 20 points in R3 , and it turns out that it has precisely 21 125 102 regularand full triangulations, up to the symmetry induced by the symmetric group of degree
four, acting on the vertices of the tetrahedron (cf. Theorem 19). This outcome is relevant,e.g., for tropical geometry, as this leads to the classification of the smooth tropical cubic
surfaces in the tropical 3-torus; see [18, §4.5]. A full account of the consequences fromthat one computation is beyond the scope of this paper.
MPTOPCOM is available as open source software which can be downloaded from https://www.polymake.org/mptopcom.
We are grateful to Jörg Rambau, Francisco Santos, Benjamin Schröter and an anonymous reviewer for valuable comments and suggestions and to Benjamin Lorenz for assistancewith the implementation.

2

Budgeted reverse search and its parallelization

Reverse search [2] is a technique for enumerating large sets of objects. Essentially, one
explores a graph Γ = (V, E) where V is the set of objects to be enumerated and the edges
are given by an adjacency oracle. For every node v, the number δ(v) denotes the outdegreeof v in Γ. The adjacency oracle Adj(v, j) returns the jth neighbor of v ∈ V for j ∈ [δ] ornull if no such neighbor exists. To apply reverse search to a problem, one provides theadjacency oracle, a local search function π(v) and an element v ∗ ∈ V . It is required that
π(v) returns a tuple (u, j) such that Adj(u, j) = v and that repeated application of π toany v ∈ V results in a path from v to v ∗ . The local search function therefore generates aspanning tree of Γ with root v ∗ . We sometimes omit the second component of π(v) when
it is not needed, e.g., in line 6 of Algorithm 1.
Reverse search (cf. Algorithm 1) traverses this spanning tree in depth-first fashion
starting from v ∗ , and at each node outputs the object corresponding to that node. One
can use the local search function to backtrack. This allows to implement reverse search ina way which does not require any additional memory after an initial setup. However, it
is also possible to cache auxiliary information needed by the application, for example to
avoid recomputation when backtracking. This way the memory usage can be controlledeven for computations with extremely large output. For more details on reverse searchsee [2].
It was recognized from the beginning that reverse search can be easily parallelized.
The enumeration process can be restarted from any node given only a description of that
node, and no node is reachable via multiple paths. This means that one does not have tostore previously-visited nodes, and communication between processes is minimal. Manyreverse search trees are highly unbalanced and so the underlying problem is to explore an
unbalanced tree in parallel.
The first parallel reverse search implementation was the generic reverse search layer
the electronic journal of combinatorics 25(3) (2018), #P3.6

4

Algorithm 1 Reverse Search
1: procedure R S(v ∗ , Adj, π)
2:v ← v ∗ , j ← 0, depth ← 0
3:repeat4:
while j < δ(v) do
5:j ←j+16:if π(Adj(v, j)) = v then7:v ← Adj(v, j)
8:
j←0
9:
depth ← depth + 1
10:Output v
11:
end if12:
end while13:
if depth > 0 then14:(v, j) ← π(v)15:depth ← depth − 1
16:
end if
17:
until depth = 0 and j = δ(v)18: end procedure
in ZRAM [6] which was applied in various areas, including the first parallel program forvertex/facet enumeration [6] and a program for certain quadratic maximization problems [9].
Other parallel reverse search applications include the computation of Minkowski sums [31].
Furthermore, Jensen computed exact mixed volumes [14] and homotopy continuation [15]via parallel reverse search implemented in Gfan [16].
Recently budgeted reverse search [4] was introduced as a simple scheme for loadbalancing in parallel reverse search. There, one adds an additional parameter (the budget)and each of the parallel processes explores its assigned subtree subject to the budget. Oncethe budget is exhausted, the process backtracks to the root of the subtree while reporting
unexplored nodes along the backtrack path. These unexplored nodes are scheduled forlater exploration. This is a particularly simple approach to parallel tree exploration thatin practice can scale beyond 1000 cores [4]. If properties of the trees generated by anapplication are known, it can be possible to prove certain performance guarantees [1].
Budgeted reverse search inherits a number of other features from reverse search. One
can checkpoint and restart the overall process by waiting for all processes to exhaust
their budget and writing the descriptions of unexplored nodes to a file. The budget canbe tuned dynamically based on the number of unexplored subtrees available. While theimplementation of [4] is specific to vertex/facet enumeration, a generic implementation ofparallel budgeted reverse search (mts) is also available along with a tutorial [3].
Implementing a budgeted reverse search application with mts allows a clean separationof the parallelization layer and application code. This permits independent developmentof the application and the parallelization layer. The current mts implementation uses

the electronic journal of combinatorics 25(3) (2018), #P3.6

5

MPI [30] and dedicates a process as master and another to handling output. This overhead
is insignificant when many processes are available, but limits parallel efficiency when only
few processes are used. We refer to Section 7 below for further details on how MPTOPCOM
employs mts for enumerating triangulations.

3

Triangulations of point configurations

Let P ⊂ Rd be a finite set of n points that affinely spans the entire space. A (polyhedral)
subdivision Σ of P is a polytopal complex which covers the convex hull conv P , such thatthe vertices of each cell form a subset of the given points P ; cf. [7, §2.3.1]. The subdivisionΣ is regular if it is induced by a height function h : P → R in the sense that the lower
convex hull of
conv (p, h(p)) | p ∈ P⊂ Rd+1projects to Σ by omitting the last coordinate. The subdivision Σ is a triangulation if allits cells are simplices. The set of all subdivisions of P is partially ordered by refinement,
and the triangulations are precisely the finest subdivisions. Our goal in this section is to
discuss and present algorithms for enumerating all regular triangulations of P .
For a given subdivision Σ the set of all height functions which induce Σ on P is arelatively open polyhedral cone, the secondary cone of Σ. The secondary cone of Σ isnon-empty if and only if Σ is regular. The set of all secondary cones forms a polyhedralfan, the secondary fan Σ-fan(P ). The relatively open secondary cones partition the space
Rn of height functions, i.e., the secondary fan is complete. Each non-empty secondary conecontains a (d + 1)-dimensional linear subspace, the space of linealities of the secondary
fan. Fixing the heights on an affine basis in P amounts to intersecting the secondary
fan of P in a way such that each cone in the resulting fan is pointed. We call that(n − d − 1)-dimensional pointed polyhedral fan the pointed secondary fan of P . It is unique
up to linear transformations.
As a key fact the pointed secondary fan is polytopal, i.e., it is the normal fan of a convex
polytope of dimension n − d − 1. We quickly review the construction. For a triangulation∆ of P the GKZ-vector is
gkz∆ = gkz∆ (p) | p ∈ P ,where gkz∆ (p) is the sum of the (Euclidean) volumes of those simplices in ∆ which containp as a vertex. The convex hullΣ-poly(P ) = conv {gkz∆ | ∆ triangulation of P }of all GKZ-vectors is the secondary polytope of P ; its normal fan is the (pointed) secondary
fan of P ; cf. [7, §5.2.2]. The vertices of Σ-poly(P ) are precisely the GKZ-vectors of
the regular triangulations. As we assumed the point configuration P is spanning, thedimension of its secondary polytope equals n − d − 1; cf. [7, §5.1.3].
We will now sketch an algorithm by Imai et al. [13] to enumerate (regular) triangulationsof P based on applying reverse search to the vertex–edge graph of Σ-poly(P ) (or a suitablethe electronic journal of combinatorics 25(3) (2018), #P3.66

supergraph). See [7, §5.3.2] for an account of the geometric facts. We begin by defining
a linear objective function on the secondary polytope Σ-poly(P ). To this end choose apositive real number M . Then the vector
λ = (M n , M n−1 , . . . , M )of length n defines a linear form on Rn . For all sufficiently large M  0 that linear form
is injective on the finite set {gkz∆ | ∆ triangulation of P }. Moreover, if M  0, thencomparing any two GKZ-vectors with respect to λ amounts to checking the lexicographic
ordering. In particular, there is no need to determine any valid choices for M . Like the
GKZ-vectors the following relies on the choice of a fixed ordering of the points in P .
Definition 1. For any two triangulations, ∆ and ∆0 , of P we let ∆0 > ∆ if
λ(gkz∆0 ) > λ(gkz∆ ) ,or if λ(gkz∆0 ) = λ(gkz∆ ) and ∆0 is lexicographically larger than ∆.
This defines a total ordering on the set of all triangulations of P . Since λ induces a
total ordering on the vertices of Σ-poly(P ) the lexicographic ordering is not required as atie-breaker if restricted to regular triangulations. However, it is important for nonregulartriangulations; see Example 8 below.
A flip f is a local modification of a triangulation, ∆, of P which yields anothertriangulation, ∆0 ; we write f = [∆∆0 ]. Each edge of the secondary polytope comesfrom a flip, but the converse does not hold [7, §5.3.1]. We call f an up-flip if ∆0 > ∆.
Otherwise ∆0 < ∆, and f is a down-flip. The flip graph Φ of P is the graph whose nodes
are the triangulations of P and whose edges are given by the flips. The flip graph Φ isdirected with respect to up-flips. Note that Φ is not necessarily connected [7, §7.3 and§7.4]. However, the subgraph Φreg induced on the subset of regular triangulations is always
connected.
Remark 2. It is known that if f = [∆∆0 ] is a flip and ∆ is regular, then ∆0 does notneed to be regular. That is, the connected component Φcreg of the regular triangulations
may be strictly larger than Φreg . We call Φcreg the regular component of the flip graph,and we call a triangulation sub-regular if it is reachable from a regular triangulation via
down-flips. Note that sub-regularity depends on the choice of ordering of triangulations.
Moreover, even if ∆ and ∆0 both are regular then f does not necessarily correspond toan edge of Σ-poly(P ). In general the vertex–edge graph of Σ-poly(P ) may be a proper
subgraph of Φreg as it may have fewer edges.
The reverse search algorithm (cf. Algorithm 1) works on a graph which is given
implicitly by means of an oracle. It is only after the termination of the procedure that theentire graph (or rather a spanning tree) is known. Therefore, we introduce the following
notation. Our adjacency oracle Φreg (∆, j) for the regular flip graph Φreg at a given regulartriangulation ∆ of P returns a regular triangulation which can be obtained from ∆ bysome down-flip, and it is the jth one in the total ordering from Definition 1. If there are
the electronic journal of combinatorics 25(3) (2018), #P3.6

7

fewer than j regular triangulations accessible via down-flips, then the oracle returns null .
We call our predecessor function π, and it returns the maximal triangulation which canbe obtained via an up-flip, if it exists. There is a unique regular triangulation ∆∗ whose
GKZ-vector is maximal among all triangulations, regular or not. We set π(∆∗ ) = null , andwe call ∆∗ the optimal triangulation of P . Recall that our ordering of the triangulations
and thus this notion of optimality relies on the choice of an ordering of the points in P .
Algorithm 2 computes all regular triangulations of P . While it can easily be modified
to, e.g., obtain all not necessarily regular triangulations which can be obtained from ∆∗ bydown-flips, the computation of the regular triangulations is probably the most interestinguse case.
Algorithm 2 Down-flip reverse search
1: procedure D F R S(P )2:∆ ← some regular triangulation of P
3:
while π(∆) 6= null do4:∆ ← π(∆)
5:
end while
6:
R S(∆, Φreg , π)
7: end procedure
To analyze Algorithm 2 we introduce three parameters. First, δmax denotes the
maximum degree of a regular triangulation in the flip graph. This depends on the point
configuration P and the ordering of its points, as this defines the optimal triangulation∆∗ . The following bounds hold:nn − d − 1 6 δmax 6
.
(1)
d+2
The lower bound is the dimension of the secondary polytope (cf. [7, Corollary 5.3.2]), andthat binomial coefficient to the right is an upper bound for the number of circuits of P .
Remark 3. The upper bound in (1) is extremely coarse, and any improvement would bevery interesting. In practice, the degree seems to be close to the lower bound, which meansthat a typical secondary polytope is somewhat close to being simple.
The second parameter, smax , is the maximal number of facets a triangulation of P may
have. It is trivially bounded bynsmax 6.
(2)
d+1
Yet, if P is a lattice polytope, we obtain a much better bound from the observation thateach lattice simplex must have normalized volume at least one. In this case, this entails
smax 6 d! · vol(conv P ) ,the electronic journal of combinatorics 25(3) (2018), #P3.6

8

where vol is the Euclidean volume. For instance, this is tight for all cubes [0, 1]d . Finally,
the third parameter, N , is the number of regular triangulations of P . The trivial upperboundn
N 6 2smax 6 2(d+1)(3)
is due to the encoding of a triangulation as the set of its maximal cells. Note that theupper bound does not benefit from counting regular triangulations only. That is, the same
bound also holds for the number of all triangulations of P .
The running time of Algorithm 2 is dominated by the combined costs for checking theregularity via linear optimization. Hence it is natural to measure the runtime complexity
in the number of LPs to be solved. In MPTOPCOM linear programs are solved via TOPCOM’s
interface to cddlib [10]. This leads to the following straightforward analysis, which
generalizes [2, Theorem 3.6] to arbitrary dimensions. For each fixed ordering of the point
set P the beneath-and-beyond method provides the corresponding placing triangulation,
which is known to be regular; cf. [7, Lemma 4.3.5]. In MPTOPCOM this serves as the initial
triangulation required in Step 2. The following is essentially [13, Theorem 13].
Theorem 4. Given some initial regular triangulation, the down-flip reverse search algorithm computes all regular triangulations of P . Its running time is bounded from above by
the time required to solve O(δmax · N ) linear programs in dimension n with at most nd
linear constraints. The space requirement is bounded from above by O (d · smax ) plus thespace required for solving the linear programs.
Here and below we use a unit cost model for representing indices corresponding topoints in P . This is adequate since, if n does not fit into a machine size Int, there will be
far too many triangulations to allow for any reasonable enumeration.
Proof. The correctness follows from (1) the fact that from each regular triangulation we
can reach the optimal triangulation ∆∗ via up-flips (cf. [7, Theorem 5.3.2]), and (2) the
correctness of the reverse search method (cf. Algorithm 1).
For the asymptotic running time we can neglect the time for up-flipping from the initial
triangulation to ∆∗ since we will revisit all triangulations on that path during the reverse
search. Each flip from a regular triangulation to another triangulation, regular or not, isconsidered at most once. The total number of such flips is bounded by δmax · N . For eachflip we need to decide if the resulting triangulation is regular or not. Given a triangulation∆ its regularity can be determined by solving a linear program in dimension n, the numberof points. Each cell
 of codimension one gives rise to one linear constraint, and there arennot more than d of these.
The space required comes from storing one triangulation as the set of its maximal cells.
Each such maximal cell in turn is a subset of the vertices of cardinality d+1, encoded as alist of integers. As pointed out above each of these integers is assumed to be small, whence
it is accounted for by a constant space requirement in our analysis.
By employing (1), (2) and (3) the above complexities can be translated into (horrendous)
bounds in terms of the input parameters d and n. Note that, by [7, Corollary 5.3.11], thethe electronic journal of combinatorics 25(3) (2018), #P3.6

9number of up-flips to the optimal triangulation does not exceed
 



n
n
,
.
min (d + 2) dd+2
b 2 + 1cExample 5. Figure 1 shows the secondary polytope of the “mother of all examples”(MOAE) from [7, Example 5.5.7]; see also [21]. This is a configuration of six points in R2
with three vertices of the convex hull and three points in the interior of the outer triangle.
There are 18 triangulations, 16 of which are regular. The 18 triangulations are grouped
into five orbits, and the coloring of the vertices in Figure 1 shows those orbits. The twononregular triangulations (are sub-regular and) share the same GKZ-vector, and theseoccur as the central blue point in a hexagon with yellow vertices. The optimal triangulation
corresponds to the one black vertex (at the top and toward the back). The thick edges
form the reverse search tree of the graph Φcreg = Φ; the optimal triangulation is the root.

Figure 1: Reverse search tree of MOAE from Example 5. The two blue points correspondto the two nonregular triangulations. Since they share the same GKZ-vector they actuallycoincide; here we chose to draw them slightly apart in order to reveal the tree.
Remark 6. A triangulation of P is full if it uses all the points in P . It was proved in [13,§5] that the subgraph of Φreg induced on the full triangulations is connected; see also [24]and [7, Corollary 5.3.14]. As a consequence, the down-flip reverse search algorithm can beapplied to enumerate the full triangulations only.
the electronic journal of combinatorics 25(3) (2018), #P3.6

10

4

Triangulations up to symmetry

In many applications for which the set of all (regular) triangulation is sought after the
point set exhibits a great deal of symmetry. Typical examples are the set of vertices of ahigh-dimensional cube (cf. Section 7.1) or the set of lattice points in a dilated simplex(cf. Section 7.3). It is natural to exploit this symmetry, and this is a standard feature of
TOPCOM and MPTOPCOM. An obvious drawback of enumerating via reverse search is that thegroup of, say, affine automorphisms of the point set does not operate on the reverse-searchtree. That is, applying reverse search up to symmetry requires some extra considerations;
see Bremner, Dutour Sikirić and Schürmann [5, §7] for a brief discussion in the context ofconvex hull computations. The simple idea is to apply the reverse search scheme to thegraph whose nodes correspond to the orbits of the (possibly regular) triangulations andwhose edges are induced by the edges in Φ, i.e., by flips.
Following Imai et al. [13] we suggest to adapt the approach via GKZ-vectors from
Section 3 to the symmetric setting. Let P ⊂ Rd be a finite point set, and let
G 6 SLd (R) o Rd
be a finite group of affine unimodular automorphisms which acts on the set P by permutations. In particular, this action is faithful, and it preserves the volume. Clearly, there isan induced action of G on the set of all triangulations of P , which leaves the set of regulartriangulations invariant.
Lemma 7. Let g be an element of G, and let ∆ be a triangulation. Then we have
gkzg·∆ (g(p)) = gkz∆ (p)for all points p ∈ P .
Proof. This is an immediate consequence of G preserving the volume.
The main challenge in implementing reverse search for enumerating triangulations is
to find a fast implementation for computing the canonical representatives. To this end wewill exhibit a procedure based on the Schreier–Sims algorithm (cf. [29, §4.2]) that works
especially well if the orbit size is nearly as large as the group size. In most relevant cases
G contains a few thousand elements, while the set of triangulations on which G acts often
has millions of orbits. Moreover, despite the fact that the point set is symmetric, few
triangulations exhibit much symmetry. This entails that most orbits are about the size ofthe entire group. See Section 7 for more details and precise numbers.
In order to employ reverse search for traversing (a connected component of) the flip
graph, the total ordering on the set of triangulations introduced in Definition 1 is crucial.
This leads us to represent the G-orbit of a triangulation ∆ of P byρ(∆) := max(G · ∆) .

(4)

That is, the canonical representative ρ(∆) of an orbit G · ∆ is characterized as follows: (i)its GKZ-vector is lexicographically maximal among all GKZ-vectors of the triangulations
the electronic journal of combinatorics 25(3) (2018), #P3.6

11

in G · ∆, and (ii) among all triangulations in G · ∆ which satisfy (i) it is lexicographicallymaximal, considered as a characteristic vector of maximal simplices.
In particular, ρ maps triangulations to triangulations in the same orbit such thatρ(∆) = ρ(∆0 ) if and only if G · ∆ = G · ∆0 , and it follows that ρ(ρ(∆)) = ρ(∆). Recall
that the above definition still depends on the ordering of the points in P as well as on theordering of the maximal simplices. If the triangulation is regular, then the entire orbit
G · ∆ consists of regular triangulations and thus their GKZ-vectors are pairwise distinct.
The following example shows that the GKZ-vectors of nonregular triangulations may
behave in an unexpected way, and this shows that condition (ii) is necessary.
Example 8. Let I 4 be the set of 16 vertices of the 4-dimensional 0/1-cube. Here andbelow we use I = [0, 1] to denote the unit interval on the real line. We may read each vertex
as a bitstring of length four, and in this way we obtain a natural encoding in terms of thehexadecimal digits 0 through F. The full group G of affine unimodular automorphisms is
the wreath product of a cyclic group of order two (corresponding to a reflection which
sends xi to 1 − xi ) by the symmetric group Sym(4) (acting on the four neighbors of a
fixed vertex); the total size of G equals 384. Consider the triangulation ∆ whose maximal
simplices read
01278 0157D 0178D 02478 0457D 0478D 1237A 1278A 137AB
178AB 178BD 189BD 2467C 2478C 267AC 278AC 478CD 67ACE
78ABC 78BCD 7ABCE 7BCDF 7BCEF
There is a flip, f = [∆simplices with∆0 ], to another triangulation, ∆0 , which replaces the underlined
2467A 478AC 467AC 2478A .

The triangulations ∆ and ∆0 lie in different G-orbits (both of which have the maximal
length 384), and their respective GKZ-vectors are(6, 10, 8, 2, 6, 2, 3, 23, 14, 1, 9, 10, 11, 10, 3, 2)
and(6, 10, 6, 2, 8, 2, 3, 23, 14, 1, 11, 10, 9, 10, 3, 2) .
Yet the lexicographically maximal GKZ-vector in both orbits is the same, and it reads(23, 3, 2, 8, 2, 6, 10, 6, 2, 3, 10, 9, 10, 11, 1, 14) .
It follows that neither ∆ nor ∆0 are regular (but they turn out to be sub-regular). The
same holds for the canonical representatives ρ(∆) and ρ(∆0 ).
This example illustrates that, in general, the GKZ-vectors cannot distinguish between
G-orbits. Moreover, it also shows that GKZ-vectors alone do not suffice to establish a
total ordering on the flip graph. For a tie-breaker we need, e.g., the lexicographic ordering
as in Definition 1.
the electronic journal of combinatorics 25(3) (2018), #P3.6

12

5

Canonical representatives via switch tables

The purpose of this section is to explain why the canonical representatives of triangulations
are chosen as in (4) and how this can be exploited. To this end we start out with a more
abstract setting. We assume that the finite group G acts on the set [m] := {0, . . . , m − 1}.
The following concept is our main tool. The underlying idea is a variation of the classicalSchreier–Sims algorithm; cf. [29, 4.2].
Definition 9. An m-switch table for G is a function
st : [m] × [m] → G(
g ∈ G, with g(k) = k, for k < i, and g(j) = i if it exists(i, j)
7→id
otherwise .
Note that we have j > i if st(i, j) 6= id. We denote by µ(st) the depth of the switch
table, which is defined asµ(st) := max {i + 1 | there is an index j with st(i, j) 6= id} .
That is, µ is the index of the first row of st which only contains the identity. In general, aswitch table is by no means unique as there may be many candidates in G for st(i, j).
The ith row of a switch table tells us which elements > i can be moved to position iwhile leaving the first i elements of [m] unchanged. That is, the switch st(i, j) lies in thestabilizer
G[i] := { g ∈ G | g(j) = j for all j < i} .
Formally, we have G[0] = G. Moreover, the stabilizer G[m] is the kernel of the action of G
on the set [m], and this is a normal subgroup. In particular, that action is faithful if andonly if G[m] is the trivial group. Let
σi := 1 + # {st(i, j) 6= id | j > i}be one plus the number of nontrivial entries in the ith row of the switch table.
Proposition 10. The ith row of the switch table forms a (left) transversal of the subgroupG[i+1] in G[i] . In particular, the number σi is the index of G[i+1] in G[i] .
Proof. Consider the case i = 0, where we need to show that the 0th row of the switch tablegives a transversal of H := G[1] in G = G[0] . The group G acts transitively on the left
cosets of H by multiplication on the left. Let g be an element in the complement G \ H.
If this does not exist there is nothing to show. Since g does not stabilize i = 0 there isan index j > 0 such that g · j = 0. By definition of the switch table there must also be aswitch s = s(0, j) with s · j = 0. We obtain that s−1 g lies in H and thus sH = gH, which
proves the claim for i = 0.
Erasing the 0th row and the 0th column of the switch table yields a switch table for
H = G[1] . This shows that, inductively, the above argument also resolves the generalcase.
the electronic journal of combinatorics 25(3) (2018), #P3.6

13

Corollary 11. For each element g ∈ G there is a unique sequence st(0, j0 ), st(1, j1 ), . . . ,st(m − 1, jm−1 ) of switches such that the productst(0, j0 ) st(1, j1 ) · · · st(m − 1, jm−1 )g −1
lies in the kernel G[m] .
Proof. From Proposition 10 we know that there is a unique switch st(0, j0 ) such thatst(0, j0 )−1 g ∈ G[0] . Inductively, we obtain st(i + 1, ji+1 ) by requiring
st(i + 1, ji+1 )−1 st(i, ji )−1 · · · st(0, j0 )−1 g ∈ G[i+1] .
The uniqueness follows as we have st(i, j) = id whenever j 6 i.
Let σ := σ0 σ1 · · · σm−1 be the product of the σi .
Corollary 12. We have σ 6 m!, and equality holds if and only if the quotient G/G[m] isthe symmetric group Sym(m). Moreover, σ agrees with the order of G if the action of G
on [m] is faithful.
Proof. We have σi 6 m − i, and this shows the upper bound. By Corollary 11 the number
σ is the index of the kernel G[m] in G. The quotient G/G[m] is isomorphic to a subgroupof Sym(m), and this characterizes the equality case.
A switch table is a special notation for a system of strong generators for G (cf. [29,§4]), tailored for finding the optimal element in an orbit; cf. Proposition 15 below. Strong
generators are used, e.g., in permlib [27]. However, the concepts differ as can be seenfrom the example below.
Example 13. Let G = Sym(4) be the symmetric group of degree four acting naturally
on the set {0, 1, 2, 3}. Using cycle notation a switch table st is given by
st(0, 1) = (0 3 2 1) , st(0, 2) = (0 2)(1 3) , st(0, 3) = (0 1 2 3) ,st(1, 2) = (1 3 2) ,st(1, 3) = (1 2 3) ,st(2, 3) = (2 3) ,and all other entries of the switch table are the identity element. The depth of this switch
table is 3. We have σ0 = 4, σ1 = 3, σ2 = 2 and σ = 24, which is the order of Sym(4).
Let Gsti be the subgroup of G, called the ith switch group, which is generated by theith row of the switch table st. The switch group Gst
0 is generated by the 4-cycle (0 1 2 3),stwhile G1 is generated by the 3-cycle (1 2 3), and Gst
2 is generated by the transpositionst(2 3). Consequently, G0 is a proper subgroup of G[0] = G. Similarly, Gst
1 is a proper
subgroup of G[1] , which is the symmetric group of degree three acting on {1, 2, 3}, whereas
st
st
Gst
2 agrees with G[2] . Observe that G2 is not a subgroup of G1 , which in turn is not astsubgroup of G0 .

the electronic journal of combinatorics 25(3) (2018), #P3.6

14

Now we consider a second action of G on some other set Ω and a map v : Ω → Rm .
Our first action of G on [m] induces a (linear) action on the vector space Rm by permuting
the coordinates. If the compatibility conditionv(g · ω) = g · v(ω)
is met we call v an m-evaluation map of the G-action on the set Ω. Due to the compatibility
we obtain an action of G on the set v(Ω) ⊂ Rm . Now we define the canonical representativeρ(v(ω)) for that action as the lexicographically maximal vector in the orbit G · v(ω). Below
we will discuss the relationship of this definition with the canonical representative of atriangulation from (4). The following example is natural.
Example 14. If G acts on [m] then this induces a second action on the set Ω of all subsetsof [m]. Sending ω ∈ Ω to its characteristic vector of length m yields an evaluation map,say χ. The canonical representative of an orbit is the lexicographically smallest set in thatorbit.
This constitutes the second part of our approach for MPTOPCOM. The two-step approach
can be summarized as choosing v to be v := (gkz, χ). The canonical representative thenhas lexicographically largest gkz-vector, and is the largest element among those with thesame gkz-vector. If one is looking at regular triangulations exclusively, one may as wellforget about the second part and choose v := gkz.
Our main method for computing canonical representatives is the procedure c a n o n i c a l in Algorithm 4. It essentially relies on the function G o o d S w i t c h e s in Algorithm 3 which determines all switches that may lexicographically improve a given vectorz ∈ Rm in the orbit of G[i] . The idea is to employ a depth first search.
Algorithm 3 Find switches which may yield a larger vector in G[i] · z
1: procedure G o o d S w i t c h e s(z, i, st)2:
Y = (y0 , y1 , . . . , y`−1 ) ← sort descending {zj | zj > zi }3:
for k = 0, 1, . . . , ` − 1 do
4:J ← {j ∈ [m] | zj = yk }5:S ← {st(i, j) | j ∈ J, st(i, j) 6= id}
6:
if S 6= ∅ then7:
return S
8:
end if
9:end for10:J ← {j ∈ [m] | zj = zi }11:
return {st(i, j) | j ∈ J}
12: end procedure
Proposition 15. Let ω ∈ Ω be an arbitrary element. For all i ∈ [m], Algorithm 4computes an element ω 0 ∈ G[i] · ω such that v(ω 0 ) is lexicographically maximal among all
elements in the orbit G[i] · ω. In particular, for i = 0, the evaluation v(ω 0 ) is the canonicalrepresentative ρ(v(ω)).
the electronic journal of combinatorics 25(3) (2018), #P3.6

15

Algorithm 4 Canonical representative of an orbit
1: procedure c a n o n i c a l(ω, v, i, st)2:
if i > µ(st) then3:
return ω
4:
end if
5:
S ←G o o d S w i t c h e s(v(ω), i, st)6:
ω0 ← ω
7:for s ∈ S do
8:
ω 00 ← c a n o n i c a l(s · ω, v, i + 1, st)9:if v(ω 00 ) > v(ω 0 ) then10:
ω 0 ← ω 00
11:
end if
12:
end for13:
return ω 0
14: end procedure
Proof. The correctness of Algorithm 4 essentially follows from Corollary 11, which provides
a transversal of G[m] as a subgroup of G in terms of switches. As that transversal arises
from an inductive construction Algorithm 4 works recursively.
It remains to discuss the case where v(ω 0 )0 = v(ω)0 . Then we need to consider all
possible switches which may or may not change the entry with index 0 but keep the v-value.
This may include the identity.
Example 16. Again we consider the symmetric group of degree four, G, acting on the set{0, 1, 2, 3} and the 4-switch table st from Example 13. As a second action of G we take theinduced action on the powerset Ω = 2{0,1,2,3} . Then the map v which sends ω ∈ Ω to its
characteristic vector is a 4-evaluation; see Example 14. We want to compute ρ({1, 2}, 0)
via Algorithm 4, i.e., we are taking the two-element set ω = {1, 2} as our input. Hence wehavev(ω) = (0, 1, 1, 0) .
The set S originates from Step 7 in Algorithm 3, and it comprises the two permutations(0 3 2 1) and (0 2)(1 3).
In the first round of the for-loop we have st(0, 1) = (0 3 2 1), and thus st(0, 1)·ω = {0, 1}.
Recursively we compute ρ({0, 1}, 1) = {0, 1}. For the second switch st(0, 2) = (0 2)(1 3)with st(0, 2) · ω = {0, 3} we get ρ({0, 3}, 1) = {0, 1}. Hence, both branches return {0, 1}
and the algorithm produces {0, 1} as its output.
For ω ∈ Ω denote by ψ(ω) the maximum number of identical entries of the vector v(ω).
Since G acts on the set v(Ω) by coordinate permutations, the number ψ(ω) is an invariantof the orbit G · ω. We let ψ be the maximal ψ(ω), taken over all elements in Ω.
Define φi as the minimum of ψ and σi ; this depends both on the switch table andon the evaluation function. Furthermore, we set φ := φ0 φ1 · · · φµ(st)−1 . By construction
we have φ 6 σ 6 #G and φ 6 m!. The benefit of our somewhat elaborate setup comesthe electronic journal of combinatorics 25(3) (2018), #P3.6

16from the fact that the number φ may be much smaller than the order of the group G; see
Example 18 below.
Corollary 17. The worst-case time complexity of Algorithm 4 is of order O(max{φ ·m, m log m}) and the worst-case space complexity is of order O(m3 ).
Proof. Time complexity: In the procedure G o o d S w i t c h e s, called with the vector
z of length m and the index i, the set E constructed in Step 2 has size at most m − 1.
Hence sorting E is of order O(m log m). The size of set J constructed in Step 5 or Step 11
is bounded by the number of identical entries ψ(ω) in the vector v(ω). Thus, the set S
constructed in Step 6 or Step 12 is bounded by the size of J as well. But the set S also
contains only non-trivial entries from the ith row of the switch table, in case of Step 12 asingle copy of the identity. Hence the set of switches returned is of size at most φi .
The height of the recursion tree of c a n o n i c a l is at most m − 1, and hence thenumber of leaves is bounded by φ0 φ1 · · · φm−1 = φ. Consequently, the total number of nodes
and also the total number of edges is of order O(φ). For each edge we lexicographicallycompare two vectors of length m. All other costs are dominated by the total complexityof these comparisons.
Space complexity: The number of entries in the switch table is given by the expressionµ(st)−1Xi=0µ(st)−1

φi 6

Xm − i 6 m2 .

(5)

i=0

Algorithm 3 returns at most one row of the switch table, and the recursion depth is boundedby the depth of the switch table. Each element of the switch table is a permutation oflength m. Thus the total space requirement amounts to O(m3 ).
Note that O(m3 ) is a coarse estimate. In practice the depth of the switch table is often
small, and then (5) provides better bounds.
Finally, we can explain our choice for the canonical representative of a triangulation
as in (4), which rests on Definition 1. Let us recall our setup. The point set P ⊂ Rd , of
cardinality n, is affinely spanning. It is equipped with a group G of affine unimodular
automorphisms. This action induces an action on the set Ω of all triangulations of P . We
can encode a triangulation as its set of d-simplices, which are the maximal cells. The latterare encoded as those (d+1)-subsets of P which are affinely spanning. If m is the number ofall d-simplices then encoding a triangulation ∆ as its characteristic vector among the setof all d-simplices yields an m-evaluation map of the action of G on Ω; see Example 14. For
computing the canonical representative of a triangulation ∆ ∈ Ω in a brute-force approach
all elements of G are applied and the lexicographically largest one is picked. Of course,
the elements of G can be precomputed once in the initialization. Since most orbits are
expected to be about the size of G this is often superior to the more traditional approachof trying generators of G until no new triangulations are found, since it requires fewer
comparisons of (characteristic vectors of) triangulations.

the electronic journal of combinatorics 25(3) (2018), #P3.6

17Yet evaluating at GKZ-vectors leads to a significant improvement. The map which
sends a triangulation ∆ to its GKZ-vector is an n-evaluation, and n is always much smaller
than m. More importantly most entries of a typical GKZ-vector are distinct and hence
the parameter ψ which enters the complexity analysis in Corollary 17 is very small. This
explains criterion (i) in Definition 1. The tie-breaker criterion (ii) is only necessary for
dealing with nonregular triangulations. The following example occurs in a computation
which should be considered small by current standards. Larger input results in larger
gains.
Example 18. Let us again look at the set P = I 4 of vertices of the 4-dimensional regularcube. The group G is the full group of affine (unimodular) automorphisms of order 384.
Mapping to GKZ-vectors yields a 16-evaluation. Any switch table has depth four with
σ0 = 16 , σ1 = 4 , σ2 = 3 , σ3 = 2 ,Their product σ = σ0 σ1 σ2 σ3 equals 384, which is the order of the group G.
To assess the complexity of computing the canonical representatives, the shape of the
GKZ-vectors is the key. Here is the full distribution of values of ψ(·) for the 247 451 orbitsof sub-regular triangulations of I 4 .
1

2

3

4

5

6

7

8

9

10

38 673

134 773

58 835

11 699

2 985

364

107

11

2

11

12
2

The average value of ψ (taken over all orbits) is approximately 3.22. This results in an
average complexity parameter of φ ≈ 66.8, which is much smaller than 384, the order ofthe group.

6

Implementation details

Our software MPTOPCOM builds on and uses existing code from mts [3], polymake [11] andTOPCOM [25]. The mts setup dedicates one process to the master and a second one for
dealing with the output; the remaining processes are reserved for the workers.
Let us give a brief overview of the interplay of the different software systems: While
parallelization is handled by mts, the single worker processes are in the domain of polymake
and TOPCOM. Triangulations and groups are handled by TOPCOM’s highly optimized code.
Since TOPCOM provides vectors and matrices already, a first implementation used these for
gkz-computation. However, replacing TOPCOM’s vectors and matrices by the corresponding
objects from polymake vastly improved performance. At this point, all vectors, matrices,and sets are handled by polymake. Checking regularity is handled by TOPCOM’s internalinterface to cddlib [10].
Parallelizing with more threads will scale almost linearly at first, but depending on thesize of the example, the curve of time consumption over number of threads will flatten
sooner or later, as can be seen in the Figures 2, 3, and 4. In order to determine whether adown-flip is valid, the reverse search algorithm descends to the target node and checks
the electronic journal of combinatorics 25(3) (2018), #P3.6

18

whether the predecessor of the target is the original node. This happens essentially for every
down-flip into a triangulation. Thus, for a given triangulation the predecessor is computed
multiple times, and hence, its neighbors are computed multiple times. If one considers
regular triangulations exclusively, then regularity is checked multiple times as well. In thesymmetric case, the canonical representative is computed several times. We attack thisproblem by maintaining three caches. The key type in every cache is a triangulation ∆, so
we just list the values.
(i) Flip cache: Contains the list of all flips of ∆.
(ii) Orbit cache: The canonical representative from the orbit of ∆.
(iii) Regularity cache: A boolean whether ∆ is regular or not.
Each worker is equipped with three such caches; they follow the least-recently-used paradigm
and can store a fixed number of keys, subject to change by the user. Caching is combinedwith hashing such that previously computed data is instantly available, if it is still cached.
There are a number of further caches with dynamic size, e.g., a cache containing thevolumes of all maximal simplices. These caches are filled at startup before the reverse
search begins and never changed again. Sharing these caches among the workers as hinted
at before is possible, but not realized yet, and would probably damage the flexibility of
MPTOPCOM.
Parallelizing does not change the number of down-flips going into a triangulation. But
since the caches are not shared among the workers, every worker populates its own caches,leading to the non-linearity in scaling. In theory, working without caches would scalelinearly, but is infeasible for larger examples, even starting with I 4 .
We use several encodings of the same triangulation simultaneously. For instance,enumerating all maximal simplices spanned by our point configuration first allows to store
a triangulation as a set of machine-size integers. Those integers are the indices pointing
into the array of maximal simplices.
The implementation is flexible enough to deal with several scenarios. In our fastest
setup we assume that the coordinates of the point configuration are machine size integers.
This means that we may use int for the entries of the GKZ-vectors. The condition on the
integrality is natural for the applications for algebraic and tropical geometry we have inmind.
As pointed out several times the total ordering on the triangulations from Definition 1
depends on an ordering of the points. In practice it is often beneficial to pick a random
ordering. This may reduce the height of the search tree for the canonical representative.
We experimented with several compilers on various kinds of hardware. By and large wefound clang, version 3.8.0, to be about 10% faster than various versions of gcc compilers.
Therefore our timings below employ clang. For the hardware we tried the following:• A cluster with four nodes, each of which is equipped with 2 × Intel Xeon E5-2630 v2
Hexa-Core (2600–3100MHz, 5201.45 bogomips) and 64GB RAM per node. On thismachine we used 40 threads. The operating system is openSUSE 42.2 with kernel
version 4.4.79.
the electronic journal of combinatorics 25(3) (2018), #P3.6

19• A desktop machine with an AMD Ryzen 7 1700 CPU. This has 8 cores/16 threads(3000 MHz, 5967.87 bogomips) and 32GB RAM. The operating system is openSUSE
42.1 with kernel version 4.12.1.
Recall that two processes are reserved for the master and the output processes. That is, nprocesses means n − 2 workers.

7

Experimental results

We tried our implementation on a number of point configurations which occur naturally
in geometric combinatorics and related areas. Our notation is as follows. The point set∆d comprises the d + 1 vertices 0, e1 , e2 , . . . , ed of the standard d-dimensional simplex,
where ei is the ith standard basis vector of Rd . We abbreviate I = [0, 1], and so I d is thed-dimensional unit cube. From these several interesting point configurations can be formed,e.g., by taking products. The point configurations of type I d or ∆p × ∆q or I p × ∆q are
in convex position, i.e., they form the vertices of polytope. For such point configurations
Table 3 shows the number of orbits of regular and sub-regular triangulations; in all cases
the group G is the full group of affine symmetries.
To the best of our knowledge the results for ∆2 × ∆6 , ∆3 × ∆4 and 3 · ∆3 are new.
MPTOPCOM-1
TOPCOMMPTOPCOM

2,500

time (s)

2,000
1,500
1,000
500
2

4

6

8
10
processes

12

14

16

Figure 2: Timings for enumerating the sub-regular triangulations of I 4 taken on AMD
Ryzen 7 1700 with 32GB RAM, depending on the number of processes. The timings forthe single-threaded version of MPTOPCOM (marked with “-1”) and TOPCOM are added forreference.

7.1

The 4-cube

A standard test case is the four-dimensional cube I 4 ; see also Example 18. The group ofaffine unimodular automorphisms has order 384, and there are 247 451 orbits of triangulations in the flip-component of the regular triangulations. By today’s standard this is
the electronic journal of combinatorics 25(3) (2018), #P3.6

20

TOPCOM
MPTOPCOM-1
MPTOPCOM

3,000

time (s)

2,000

1,000

0
2

4

6

8

10 12 14 16 18 20 22 24 26 28 30 32 34
processes

Figure 3: Timings for enumerating the sub-regular triangulations of I 4 taken on the Intel
Xeon E5-2630 v2 cluster with 64GB RAM per node, depending on the number of processes.
to be considered small input. Out of these, actually 235 277 are regular [7, Thm. 6.3.12].
Pournin proved that the flip graph of I 4 is connected [23]. Since MPTOPCOM also counts
247 451 orbits it follows that each triangulation of I 4 is sub-regular. The total number oftriangulations of I 4 adds up to 92 487 256.
On an AMD Ryzen 7 1700 TOPCOM takes 2345 seconds to enumerate the regularcomponent, while MPTOPCOM requires 2598 seconds (single-threaded) and 457 seconds (with10 processes), respectively. Figure 2 shows how MPTOPCOM scales with the number ofprocesses. It shows that, on a standard desktop computer, already four processes, i.e.,
two workers, suffice for MPTOPCOM to be substantially faster than TOPCOM. This should
be compared with Figure 3 which shows very similar behavior on our cluster. The only
exception is that, on that hardware, even a single-threaded version of MPTOPCOM beats
TOPCOM. The single-threaded MPTOPCOM-1 is the pure down-flip reverse search algorithm
built without the overhead of mts and MPI.
Next we try to give an idea about which percentage of the total running time is spenton which subtask. Table 1 shows the values for MPTOPCOM-1 during the computation of all
sub-regular triangulations of I 4 . For comparison the relative timings for TOPCOM are given
in Table 2. All these numbers were determined with Valgrind’s tool callgrind [19]. In
both cases the bulk of the time is spent on finding and processing the flips. For MPTOPCOM
the major subtask is to determine the canonical representatives, while TOPCOM will explicitly
compute the full orbits of each triangulation that it visits. For both programs the costfor finding the initial triangulation (e.g., in Algorithm 2 before Step 6) is negligible. The
picture changes entirely if one restricts the algorithms to enumerate regular triangulations
only. Then MPTOPCOM spends 95% of its time on solving linear programs, while that markfor TOPCOM reaches 98%. The remaining time is used in a similar fashion to that seen in
Tables 1 and 2. For larger examples the overall pattern stays the same, but the most costly
subtasks tend to take up even higher proportions of the total running times.

the electronic journal of combinatorics 25(3) (2018), #P3.6

21

7.2

Products of two simplices

Another interesting class of point configurations are the products ∆p × ∆q of two simplices[7, §6.2]. The natural group action is by the product Sym(p + 1) × Sym(q + 1) of symmetricgroups. In tropical geometry, e.g., their regular subdivisions control the combinatorial typesof tropical polytopes [18, §5.2]. The special case where one of the factors is one-dimensional,i.e., when the product of simplices is a prism, is fully understood [7, §6.2.1]. Therefore, in
our experiments we restrict our attention to cases with 2 6 p 6 q. There is a formula forthe number of all triangulations of ∆2 × ∆q [7, 9.2.5], but this does not immediately yield
the number of (semi-)regular triangulations or the number of orbits.
Figure 4 shows the speed for computing the triangulations in the regular componentof ∆2 × ∆6 with MPTOPCOM on the cluster, depending on the number of processes. This
computation is medium size, i.e., a bit larger than the previous, and so it pays to use
more processes. The timings are as follows: one hour and 22 minutes with MPTOPCOM (10processes), five hours with the single-threaded version MPTOPCOM-1 and eight days and16 hours with TOPCOM. A computation of this kind essentially marks the end of TOPCOM’s
range.
Our new results for ∆2 × ∆6 and ∆3 × ∆4 helped Schröter [28] to obtain new resultson coarsest subdivisions of hypersimplices. An attempt to handle ∆3 × ∆5 is currently
under way (running on more than a hundred cores for some weeks). So far it has foundmore than 900 million orbits of sub-regular triangulations.
7.3

Dilated simplices
A third class of point configurations is denoted as k · ∆d . These are the n+d
= n+d
n
d
lattice points in the simplex ∆d which is dilated by the factor k. For any polynomial
in d + 1 indeterminates which is homogeneous of degree k the monomials correspond topoints in the point configuration k · ∆d . In particular, the vertices of the Newton polytope
form a subconfiguration. It follows that the tropical hypersurfaces in the tropical d-torus
Td+1 /R1 of homogeneous degree k are dual to regular subdivisions of k · ∆d ; see [18, §3.1].
The regular unimodular triangulations of k · ∆d , which are necessarily full, correspond to
those tropical hypersurfaces which are smooth. For the first time, we computed the full
Table 1: Percentages of total running time spent by MPTOPCOM-1 for computing all (subregular) triangulations of the 4-cube. The value for flip processing further refined in second
column.
95%

Process flips
68% Compute canonical representative
20% Partition into up- and down-flips
7% all remaining
3%
Check whether flip is edge of reverse search tree2%
all remaining

the electronic journal of combinatorics 25(3) (2018), #P3.6

22triangulations of 3 · ∆3 , and these classify the smooth tropical cubics in 3-space [18, §4.5].
Theorem 19. There are exactly 21 125 102 orbits of regular and full triangulations of
3 · ∆3 with respect to the natural action of the symmetric group of degree four. Out ofthese, 14 373 645 are unimodular.
This is the largest experiment that we completed so far. The computation took aboutfour days on the Intel Xeon E5-2630 v2 cluster with 40 threads.
Cache size
100
1000

4

10

2000
5000

time (s)

10000
20000

103
2

4

6

8

10 12 14 16 18 20 22 24 26 28 30 32 34
processes

Figure 4: Timings for enumerating the triangulations with MPTOPCOM of ∆2 × ∆6 in theregular component taken on the Intel Xeon E5-2630 v2 cluster with 64GB RAM per node,depending on the number of processes and cache sizes. Each data point is the result of
averaging over ten runs with the same parameters.

7.4

Using more memory for caching

On the one hand, as a key benefit, the memory consumption of the down-flip reverse searchalgorithm allows for an excellent a priori estimate which is also quite low; cf. Theorem 4.
Table 2: Percentages of total running time spent by TOPCOM for computing all triangulationsof the 4-cube. The values for the most expensive subtasks are refined in the second andthird columns.
99%

Process flips
Check whether class was already found
82% Enumerate orbits
10% all remaining
7%
all remaining
all remaining

92%

1%

the electronic journal of combinatorics 25(3) (2018), #P3.6

23

Cache size
20000
10000
5000

memory (MB)memory (MB)

103

2000

103
102

Cache size
100000
50000
20000

1

10

10000

1000

default

100

no cache

2

10

0

20
40
snapshot no

0

20
40
snapshot no

Figure 5: Memory usage determined with Valgrind. We took 52 massif snapshots for∆2 × ∆6 (left) and 3 · ∆3 (right).
On the other hand this approach results in a considerable amount of duplication. To
avoid at least some of this it is natural to employ caching, as explained in Section 6. Herewe want to report on some experiments concerning the impact of caching on the overallrunning time.
First we investigate the medium-size example ∆2 × ∆6 from Section 7.2. It has 533 242
sub-regular triangulations up to symmetry. Figure 4 shows how the running-time depends
on the number of workers and on the cache sizes of the three main caches introduced in
Section 6. Their size can be varied individually, but for this plot we gave all of them thesame size. As a default MPTOPCOM stores 2000 triangulations in each of these caches. In
this case increasing the cache size from 100 to 2000 reduces the running time by about one
third, independent of the number of threads. Increasing the cache to 20 000 only results ina further reduction by another 5%.
Figure 5 shows how the total amount of memory consumed depends on the cache sizes.
The measurements have been taken by Valgrind’s tool massif which records memory
snapshots in fixed time intervals [19]. The left hand side corresponds to the computationin Figure 4 for ∆2 × ∆6 . By and large the overall memory consumption depends linearlyon the cache sizes. Note that the sizes of the triangulations as well as the sizes of thevarious objects for the cache values vary. So some fluctuations should be expected. Indeed,this is visible for ∆2 × ∆6 , which is not very large. In the much larger example 3 · ∆3 from
Section 7.3 we see a similar behavior, but the fluctuations are nearly gone.

8

Concluding remarks

Table 3 also contains some empty rows, where we do not know the respective numberof triangulations. Most of these will be out of reach for the current implementations,
including MPTOPCOM. The reason for listing these nonetheless is to give a feel for the orders
the electronic journal of combinatorics 25(3) (2018), #P3.6

24

Table 3: Summary of enumerations.
n = #Pd = dim P

n−d

#G

I3
I4
I5

8
16
32

3
4
5

5
12
27

48
384
3 840

6
235 2776
247 451∆2 × ∆2∆2 × ∆3∆2 × ∆4∆2 × ∆5∆2 × ∆6

9
12
15
18
21

4
5
6
7
8

5
36
7
144
9
720
11 4 320
13 30 240

535530
13 621
531 862

535530
13 629533 242∆3 × ∆3∆3 × ∆4∆3 × ∆5

16
20
24

6
7
8

10
13
16

576
2 880
17 280

7 869
7 051 957

7 955
7 402 421
> 9 · 108

2∆3
3∆3
4∆3

10
20
35

3
3
3

7
17
32

24
24
24

15
21 125 102

59
925 148 763

Points P#triangulation orbits(full) regular sub-regularof magnitude involved. One main complexity parameter for enumerating triangulations isthe difference n − d of the number of points and the dimension. This is also one plus thedimension of the secondary fan, modulo linealities. Our experiments suggest that, as avery rough estimate, the range for TOPCOM seems to be limited by n − d ≈ 13. This bar
is raised substantially by MPTOPCOM to cover point configurations with n − d = 17 suchas 3 · ∆3 . It is an interesting question if MPTOPCOM can, e.g., deal with I 3 × ∆2 wheren − d = 19. That particular point configuration played a role in work of Orden and Santos[20] on efficient triangulations of cubes; see also [7, §6.3.3].
The empty rows of Table 3 show some cases which seem to be rather difficult challenges,
with the current techniques. This includes the five-dimensional cube I 5 or the dilated
simplex 4 · ∆3 . Proving results about their triangulations might require clever strategies
for random probing.
Another direction which looks promising is to investigate the triangulations of the cyclic
polytopes. This is related to the higher Stasheff–Tamari orders which were introduced by
Kaparanov and Voevodsky [17] and studied, e.g., by Edelman and Reiner [8]; see Rambau
and Reiner [26] for a survey.

References[1] David Avis and Luc Devroye, An analysis of budgeted parallel search on conditional
Galton-Watson trees, 2017, Preprint arXiv:1703.10731.

the electronic journal of combinatorics 25(3) (2018), #P3.625

[2] David Avis and Komei Fukuda, Reverse search for enumeration., Discrete Appl. Math.
65 (1996), no. 1-3, 21–46 (English).
[3] David Avis and Charles Jordan, A parallel framework for reverse search using mts,
2016, Preprint arXiv:1610.07735.
[4], mplrs: A scalable parallel vertex/facet enumeration code, Mathematical
Programming Computation 10 (2018), no. 2, 267–302.
[5] David Bremner, Mathieu Dutour Sikirić, and Achill Schürmann, Polyhedral representation conversion up to symmetries, Polyhedral computation, CRM Proc. Lecture
Notes, vol. 48, Amer. Math. Soc., Providence, RI, 2009, pp. 45–71. MR 2503772[6] Adrian Brüngger, Ambros Marzetta, Komei Fukuda, and Jurg Nievergelt, The parallel
search bench ZRAM and its applications, Annals of Operations Research 90 (1999),
45–63.
[7] Jesús A. De Loera, Jörg Rambau, and Francisco Santos, Triangulations, Algorithms
and Computation in Mathematics, vol. 25, Springer-Verlag, Berlin, 2010, Structures
for algorithms and applications. MR 2743368 (2011j:52037)[8] Paul H. Edelman and Victor Reiner, The higher Stasheff-Tamari posets., Mathematika
43 (1996), no. 1, 127–154 (English).
[9] J.-A. Ferrez, K. Fukuda, and Th. M. Liebling, Solving the fixed rank convex quadratic
maximization in binary variables by a parallel zonotope construction algorithm, European J. Oper. Res. 166 (2005), no. 1, 35–50. MR 2128976[10] Komei Fukuda, cddlib, version 0.94h, http://www.inf.ethz.ch/personal/
fukudak/cdd_home/, 2015.
[11] Ewgenij Gawrilow and Michael Joswig, polymake: a framework for analyzing convex
polytopes, Polytopes—combinatorics and computation (Oberwolfach, 1997), DMV
Sem., vol. 29, Birkhäuser, Basel, 2000, pp. 43–73. MR MR1785292 (2001f:52033)[12] I. M. Gel0 fand, M. M. Kapranov, and A. V. Zelevinsky, Discriminants, resultants
and multidimensional determinants, Modern Birkhäuser Classics, Birkhäuser Boston
Inc., Boston, MA, 2008, Reprint of the 1994 edition. MR MR2394437 (2009a:14065)[13] Hiroshi Imai, Tomonari Masada, Fumihiko Takeuchi, and Keiko Imai, Enumerating
triangulations in general dimensions, Internat. J. Comput. Geom. Appl. 12 (2002),
no. 6, 455–480. MR 1945594[14] Anders N. Jensen, An implementation of exact mixed volume computation, Mathematical Software – ICMS 2016: 5th International Conference, Berlin, Germany, July
11-14, 2016, Proceedings (Gert-Martin Greuel, Thorsten Koch, Peter Paule, andAndrew Sommese, eds.), Springer International Publishing, Cham, 2016, pp. 198–205.
[15], Tropical homotopy continuation, 2016, Preprint arXiv:1601.02818.
[16], Gfan, a software system for Gröbner fans and tropical varieties, version 0.6,
Available at http://home.imf.au.dk/jensen/software/gfan/gfan.html, 2017.
[17] Mikhail M. Kapranov and Vladimir A. Voevodsky, Combinatorial-geometric aspectsof polycategory theory: pasting schemes and higher Bruhat orders (list of results),
the electronic journal of combinatorics 25(3) (2018), #P3.6

26

Cahiers Topologie Géom. Différentielle Catég. 32 (1991), no. 1, 11–27, International
Category Theory Meeting (Bangor, 1989 and Cambridge, 1990). MR 1130400[18] Diane Maclagan and Bernd Sturmfels, Introduction to tropical geometry, Graduate
Studies in Mathematics, vol. 161, American Mathematical Society, Providence, RI,
2015. MR 3287221[19] Nicholas Nethercote and Julian Seward, Valgrind: A framework for heavyweight
dynamic binary instrumentation, Proceedings of the 28th ACM SIGPLAN Conference
on Programming Language Design and Implementation (New York, NY, USA), PLDI
’07, ACM, 2007, pp. 89–100.
[20] David Orden and Francisco Santos, Asymptotically efficient triangulations of thed-cube, Discrete Comput. Geom. 30 (2003), no. 4, 509–528. MR 2013970[21] Julian Pfeifle, Secondary polytope of the “mother of all examples”, Electronic Geometry
Models (2000), No. 2000.09.033.
[22] Julian Pfeifle and Jörg Rambau, Computing triangulations using oriented matroids,
Algebra, geometry, and software systems, Springer, Berlin, 2003, pp. 49–75. MR
2011753 (2004i:68233)[23] Lionel Pournin, The flip-graph of the 4-dimensional cube is connected, Discrete Comput.
Geom. 49 (2013), no. 3, 511–530. MR 3038527[24] Lionel Pournin and Thomas M. Liebling, Constrained paths in the flip-graph of regulartriangulations, Comput. Geom. 37 (2007), no. 2, 134–140. MR 2310598[25] Jörg Rambau, TOPCOM: triangulations of point configurations and oriented matroids,
Mathematical software (Beijing, 2002), World Sci. Publ., River Edge, NJ, 2002,
pp. 330–340. MR 1932619[26] Jörg Rambau and Victor Reiner, A survey of the higher Stasheff-Tamari orders.,
Associahedra, Tamari lattices and related structures. Tamari memorial Festschrift,
Basel: Birkhäuser, 2012, pp. 351–390 (English).
[27] Thomas Rehn and Achill Schürmann, C++ tools for exploiting polyhedral symmetries.,
Mathematical software – ICMS 2010. Third international congress on mathematicalsoftware, Kobe, Japan, September 13–17, 2010. Proceedings, Berlin: Springer, 2010,
pp. 295–298 (English).
[28] Benjamin Schröter, Multi-splits and tropical linear spaces from nested matroids, 2017,
Preprint arXiv:1707.02814.
[29] Ákos Seress, Permutation group algorithms, Cambridge Tracts in Mathematics, vol.
152, Cambridge University Press, Cambridge, 2003. MR 1970241[30] Marc Snir, Steve Otto, Steven Huss-Lederman, David Walker, and Jack Dongarra,
MPI– the complete reference, vol 1: The MPI core, 2nd ed., MIT Press, 1998.
[31] Christophe Weibel, Implementation and parallelization of a reverse-search algorithm
for Minkowski sums, 2010 Proceedings of the Twelfth Workshop on Algorithm Engineering and Experiments (ALENEX), 2010, pp. 34–42.

the electronic journal of combinatorics 25(3) (2018), #P3.6

27