1 Introduction

To tell whether a given space is homeomorphic to the sphere in a given dimension is a basic problem in computational topology. However, this is difficult in an essential way.

Theorem 1

(S. P. Novikov; cf. Volodin et al. 1974; Chernavsky and Leksine 2006) Given a d-dimensional finite simplicial complex K it is undecidable to check if K is homeomorphic to \(\mathbb {S}^d\) for \(d\ge 5\).

We will consider closed manifolds encoded as finite abstract simplicial complexes—but the methods and results in this article also hold for more general cell complexes with little modification. For a brief historical overview: d-sphere recognition is trivial in dimensions \(d\le 2\). Rubinstein (1995) and Thompson (1994) proved that 3-sphere recognition is decidable. Subsequently, Schleimer (2011) showed that 3-sphere recognition lies in the complexity class NP, and Lackenby (2021) proved that this problem also lies in co-NP; see Lackenby (2020) for a recent survey. The complexity status of 4-sphere recognition is open. Summing up, we do not know of any efficient algorithm for d-sphere recognition in the relevant dimensions \(d\ge 3\).

Our point of departure is that the sphere recognition problem does not go away simply because it is algorithmically intractable. To the contrary it appears naturally, e.g., in the context of manifold recognition, which is the task of deciding whether a given simplicial complex triangulates any manifold and finding its type. In the piecewise linear (PL) category, recognizing whether a given complex triangulates some PL manifold can be reduced to PL sphere recognition since the links of all vertices of the given complex need to be PL spheres, sometimes also called standard PL sphere. This plays a role, e.g., for enumerating all manifolds with few vertices or facets (Brehm and Kühnel 1992; Burton 2011, 2013; Casali and Cristofori 2015; Lutz 1999; Sulanke and Lutz 2009); for detecting errors in experimental topological constructions (Adiprasito et al. 2017; Spreer and Kühnel 2011; Tsuruga and Lutz 2013); or for meshing (Saucan et al. 2008).

In the absence of a general sphere recognition procedure the next best thing are certificates for sphericity and non-sphericity, respectively. A discrete Morse function, \(\mu \), on a finite d-dimensional abstract simplicial complex, K, may be encoded as an acyclic partial matching in the Hasse diagram of the partial ordering of the faces of K; cf. Forman (1998, 2002) and Chari (2000). The critical faces are those unmatched, and \((c_0,c_1,\dots ,c_d)\) is the discrete Morse vector of \(\mu \), where \(c_k\) is the number of critical k-faces. We call such a discrete Morse vector spherical if \(c_0=c_d=1\) and \(c_k=0\) otherwise. The relevance for our topic comes from the following key result.

Theorem 2

(Whitehead 1939; Forman 1998, 2002) A combinatorial d-manifold is a PL d-sphere if and only if it admits some subdivision with a spherical discrete Morse vector.

So we propose a heuristic method for sphere recognition which navigates between Theorems 1 and 2. There are a few more obstacles though. Adiprasito and Izmestiev (2015) showed that a sufficiently large iterated barycentric subdivision of any PL sphere is polytopal (and thus inherits a spherical discrete Morse function from linear programming). However, in view of Theorem 1, there cannot be any a priori bound on the number of barycentric subdivisions required to attain polytopality. Second, deciding whether a discrete Morse function with at most a fixed number k of critical cells exists is NP-hard (Joswig and Pfetsch 2006; Lewiner et al. 2003), intractable in the parameter k (Burton et al. 2016), and not even a polynomial approximation is available (Bauer and Rathod 2019). Finally, there are combinatorial d-spheres that do not admit any spherical discrete Morse function (Benedetti and Lutz 2013a; Benedetti and Ziegler 2011).

This article is the considerably expanded full version of the extended abstract (Joswig et al. 2014). It is organized as follows. As our first main contribution, in Sect. 2, we present an implementation of a sphere recognition heuristic procedure in polymake, and demonstrate its efficiency. In the polymake project, Perl and C++ are used as programming languages; our heuristic is implemented in C++. It is also available through the new Julia interface layer Polymake.jl (Kaluba et al. 2020), which supports the current polymake Version 4.5. Section 3 comprises comprehensive computational experiments which show that there are many randomly constructed, even fairly large, simplicial complexes for which deciding sphericity is surprisingly easy; this agrees with previous observations (Adiprasito et al. 2017; Benedetti and Lutz 2014). Moreover, on such input our new approach proves to be superior to, e.g., the 3-sphere recognition implemented in Regina (Burton et al. 1999–2021), which is a standard tool in computational topology. Note that Regina ’s method is a full decision algorithm, while our heuristic may be inconclusive. However, we are not aware of any triangulation of \(\mathbb {S}^3\) which cannot be recognized by our method. Another experiment comes from a census of 4-manifolds provided by Regina; here our heuristics recognizes about \(99.5\%\) of the input as spheres or non-spheres. Finally, in Sect. 4, we explore the limitations of our method. One outcome is the construction of a new family of 2-dimensional cell complexes which are contractible, but not collapsible. These saw blade complexes generalize the Dunce hat, and in our experiments they occur naturally as one source of difficulty for recognizing spheres. Moreover, our computer experiments show that there is a “horizon” for discrete Morse computations, along with implications to homology computations and computational topology in general.

2 A heuristic sphere recognition scheme

We describe our procedure for sphere recognition and its implementation in polymake (Gawrilow et al. 1997). This is the specification:

Input: A d-dimensional (finite abstract simplicial) complex K with n vertices and m facets, where a facet is a face that is maximal with respect to inclusion.

Output: Yes, No, or Undecided, depending whether K has been recognized as a (standard) PL d-sphere.

Our procedure features five steps, labeled (0) through (4). Discussing the trivial preprocessing Step (0) in some details allows us to introduce the basic terminology and notation. The core Steps (1), (2), (3), and (4) below together yield Algorithm 1.

figure a

2.1 (0) Preprocessing

To verify whether K is a PL d-sphere, there are three elementary combinatorial checks that are useful to perform first. These checks are fast; their running time is bounded by a low-degree polynomial in the parameters d, n and m. If one of the checks fails, this will serve as the certificate that K is not a sphere.

More precisely, we first check if K is pure, i.e., each facet has exactly \(d+1\) vertices. Second, we check if each ridge is contained in exactly two facets, where a ridge is a face of dimension \(d{-}1\). Success in these two tests will assert that K is a weak pseudo-manifold (without boundary). Note that the 0-dimensional sphere \(\mathbb {S}^0\) is a weak pseudo-manifold of dimension \(d=0\) with two isolated vertices.

Third, for \(d\ge 1\), we check if the 1-skeleton of K is a connected graph. A connected weak pseudo-manifold K of dimension \(d=1\) is a polygon and thus triangulates \(\mathbb {S}^1\).

The pureness and the weak pseudo-manifold property of a simplicial complex is inherited by all face links; cf. (Bagchi and Datta 1998, Rem. 8). A (connected) weak pseudo-manifold is a pseudo-manifold if it is strongly connected, i.e., if any two of its facets can be joined by a sequence of facets for which consecutive facets share a ridge. In particular, a pseudo-manifold of dimension \(d=2\) is a triangulation of a closed surface or of a closed surface with pinch points (having multiple disjoint cycles as vertex links).

A d-dimensional pseudo-manifold K is a combinatorial d-manifold if all vertex links of K are PL homeomorphic to the boundary of the d-simplex or, equivalently, if for every proper i-face (with \(0\le i<d\)) of K its link is a PL \((d{-}i{-}1)\)-sphere; here the case \(i=d-1\) ensures the weak pseudo-manifold property. This recursive nesting of PL spheres suggests an inductive check of the face links of K by dimension, starting with 1-dimensional links of \((d-2)\)-faces and proceeding up until the \((d-1)\)-dimensional links of the vertices.

A connected 2-dimensional weak pseudo-manifold K whose vertex links are single cycles is a combinatorial 2-manifold and triangulates a closed surface. If, additionally, the Euler characteristic of K equals two, then K is \(\mathbb {S}^2\).

If one of the checks on the links of the overall complex K fails, then K cannot be a standard PL sphere, if one of the checks is left undecided this leaves K undecided.

After this preprocessing and an inductive check of the vertex-links we may assume that our input looks as follows:

Input (modified): Let K be a d-dimensional combinatorial manifold, for \(d\ge 3\).

The subsequent four steps form the core of our sphere recognition procedure.

2.2 (1) Homology computation

Computing the simplicial homology modules of a finite simplicial complex is a standard procedure, which is implemented, e.g., in CHomP (http://chomp.rutgers.edu), RedHom (http://redhom.ii.uj.edu.pl) or polymake (Gawrilow et al. 1997). The homology with field coefficients can be determined via applying Gaussian elimination to the (simplicial) boundary matrices. For finite fields of prime order or the rationals this can be achieved in polynomial time (in the size of the boundary matrices); cf. Edmonds (1967). Similarly, over the integers, a homology computation can be reduced to computing Smith normal forms; cf. (Munkres 1984, Ch. 11). Kannan and Bachem (1979) gave the first polynomial time Smith normal form algorithm, employing modular arithmetic; see also Iliopoulos (1989).

Here we employ integer coefficients throughout. A necessary condition for K to be a sphere (PL or not) is \(H_d(K)\cong \mathbb {Z}\), and all other (reduced) homology groups vanish. In this case we say that K has spherical homology.

The Hasse diagram of a simplicial complex K is a directed graph with one node per face of K and a directed arc \((\sigma ,\tau )\) if the face \(\sigma \) is contained in \(\tau \) and \(\dim \tau =\dim \sigma +1\). The non-zero coefficients in the kth boundary matrix, which maps k-faces to \((k{-}1)\)-faces, bijectively correspond to the arcs \((\sigma ,\tau )\) in the Hasse diagram for \(k=\dim \tau \) (and thus \(\dim \sigma =k-1\)). By construction the Hasse diagram is an acyclic graph.

While the modular approach of Kannan and Bachem (1979) and Iliopoulos (1989) is valid for matrices with arbitrary integer coefficients, simplicial boundary matrices have entries 1, \(-1\), and 0 only. As a consequence, in an arbitrary simplicial boundary matrix it is always possible to perform at least a few Gauss elimination steps. Moreover, a typical boundary matrix is sparse. If the matrix happens to stay sparse during the elimination and if, additionally, one does not run out of unit coefficients too soon (such that it is possible to continue with elimination steps), an elimination based Smith normal form algorithm can outperform the more sophisticated modular methods. This is why for computations of (simplicial) homology elimination algorithms are often preferred; cf. Dumas et al. (2003) for a survey.

A (partial) matching in an arbitrary graph is a subset of the edges such that each node is covered at most once. In the Hasse diagram of K a matching corresponds to a set of non-zero coefficients in some boundary matrices. Such a matching, \(\mu \), is called acyclic if reversing all arcs in \(\mu \) (and keeping the arcs not in \(\mu \)) still gives an acyclic graph. It is easy to see that an acyclic matching in the Hasse diagram of K yields a sequence of Gauss elimination steps that can be performed in any order without destroying the (unit) pivots required for the subsequent elimination steps.

Fig. 1
figure 1

Acyclic matching in the Hasse diagram of \(\mathbb {R}\mathbb {P}^2_6\). The three unmatched cells are marked

Full size image
Fig. 2
figure 2

Boundary matrices of \(\mathbb {R}\mathbb {P}^2_6\) with coefficients marked that correspond to the acyclic matching of Fig. 1; cf. Example 3

Full size image

Example 3

Figure 1 shows an acyclic matching, \(\mu \), in the Hasse diagram for \(K=\mathbb {R}\mathbb {P}^2_6\), which is the six-vertex triangulation of the real projective plane. Figure 2 shows the corresponding boundary matrices. The pivots corresponding to \(\mu \) are marked. Using these pivots in an arbitrary order yields an elimination strategy for the computation of the homology modules:

$$\begin{aligned} {\widetilde{H}}_0(\mathbb {R}\mathbb {P}^2_6) = 0 \,,\ H_1(\mathbb {R}\mathbb {P}^2_6) \cong \mathbb {Z}/2\mathbb {Z}\,,\ H_2(\mathbb {R}\mathbb {P}^2_6) = 0. \end{aligned}$$

2.3 (2) Random discrete Morse functions

A map \(f:K\rightarrow \mathbb {R}\) which assigns a real number to each face of K is a discrete Morse function if for every k-face \(\sigma \) of K we have

$$\begin{aligned} \begin{aligned}&\#\{ \tau \in K \, \mid \, f(\tau )\le f(\sigma ),\, \sigma \subset \tau ,\, \dim \tau =\dim \sigma +1\} \le 1 \quad \text {and}\\&\#\{\rho \in K \, \mid \, f(\rho )\ge f(\sigma ),\, \rho \subset \sigma ,\, \dim \rho =\dim \sigma -1\} \le 1. \end{aligned} \end{aligned}$$
(1)

A k-face is critical with respect to f if both sets in Condition (1) are empty; the non-critical faces are regular, and they form an acyclic matching on the Hasse diagram of K; cf. Chari (2000). In this sense the acyclic matchings form equivalence classes of discrete Morse functions. These concepts were introduced by Forman (1998, 2002). The discrete Morse vector \((c_0,c_1,\dots ,c_d)\) of an acyclic matching counts the critical faces per dimension; and K is homotopy equivalent to a CW complex with \(c_i\) cells in dimension \(0\le i\le d\). Let \(\mathbb {F}\) be some field. A discrete Morse vector is \(\mathbb {F}\)-perfect for K if \(c_i=\beta _i(K;\mathbb {F})\) for \(0\le i\le d\).

Example 4

Figure 1 shows an acyclic matching for \(\mathbb {R}\mathbb {P}^2_6\) with three critical cells. The corresponding discrete Morse vector (1, 1, 1) is \(\mathbb {Z}_2\)-perfect. The real projective plane admits a CW-complex structure with one 0-cell, one 1-cell, and one 2-cell.

Discrete Morse vectors which are \(\mathbb {F}\)-perfect for any field \(\mathbb {F}\) are perfect. A perfect discrete Morse vector of a sphere, which reads \((1,0,0,\dots ,0,1)\), is also called spherical. Theorem 2 implies that a combinatorial d-manifold K that becomes collapsible after the removal of one facet is a PL d-sphere. In 1992, Brehm and Kühnel (1992) used that fact to show that some 8-dimensional simplicial complex with 15 vertices is a combinatorial 8-manifold (a triangulation of the quarternionic projective plane (Gorodkov 2019)).

By Theorem 2, the existence of an acyclic matching whose discrete Morse vector is spherical is a sufficient criterion for a combinatorial d-manifold K to be a standard PL sphere. This gives rise to the following simple strategy: generate discrete Morse functions (or acyclic matchings) at random and check if one of them is spherical; cf. Benedetti and Lutz (2014).

The random_discrete_morse function implemented in polymake has three random strategies which we call random-random, random-lex-first, and random-lex-last. We will give a short outline and describe the differences among the three strategies and further differences to the original approach from Benedetti and Lutz (2014).

Let K be an arbitrary d-dimensional simplicial complex, which is not necessarily a manifold. A free face of K is an \((i{-}1)\)-dimensional face that is contained in exactly one i-face, \(0< i \le d\). To save memory, our three strategies are destructive in the sense that they keep changing the complex K. In each step we pick one of the free faces of codimension one and delete it from K together with the unique d-face containing it. This is an elementary collapse, and the two removed faces form a regular pair, which is a matching edge in the Hasse diagram. The three strategies only differ in how they pick the free face. If we run out of free faces, we pick some facet (of maximal dimension), declare it critical and remove it. After removing a regular pair the dimension of the resulting complex, \(K'\), may drop to \(d-1\). This process continues until \(K'\) is zero-dimensional. In this case, \(K'\) only consists of vertices, all of which are declared critical.

For the random-random strategy, we first find all the free faces of K and collect them in a linked list. If this list is not empty, choose a free face uniformly at random. Taking the uniform distribution means that each free face has a fair chance of being taken, but this comes at a price since the sampling itself takes time if there are many free faces to choose from. The reason is that we do not have random access to the free faces, as they are kept in a linked list. Picking a random element in a linked list takes linear time in the length of that list. If we run out of free faces, the choice of the critical d-face is again uniformly at random.

The strategy random-random is somehow the obvious one, but there is a much cheaper way which maintains a certain amount of randomness. Here the price is that it seems to be difficult to say something about the resulting probability distribution. The idea is to randomly relabel the vertices of K once, at the beginning, and then to pick the free and critical faces in a deterministic way (which depends on the random labeling). Whenever we want to choose a free or critical face, rather than selecting one at random, we pick the first (in the case of random-lex-first) or the last one (in the case of random-lex-last) of the linked list. The random-lex-last strategy was called “random-revlex” in Benedetti and Lutz (2014). We changed the name here to random-lex-last to avoid confusion with the reverse lexicographic (term) ordering, which is different.

The cost of being fair is quite significant, with our current implementation, when dealing with large complexes. For example, running the random-lex-first and random-lex-last strategies on the fourth barycentric subdivision of \(\partial \Delta ^4\) took less than 3 min per run whereas the random-random strategy took approximately 2 h per run; see Sect. 4.4. It is conceivable that there is some room for improvement here by employing a faster data structure for random sampling; we leave such an implementation for a future version of polymake.

Remark 5

In Algorithm 1, the Steps (1) and (2) can also be intertwined as finding an acyclic matching results in a partial strategy for computing the homology. To this end it is most natural to process the Hasse diagram from top to bottom level by level.

2.4 (3) Random bistellar flips

If the previous tests are inconclusive, we can use a local search strategy to determine the PL type; cf. Björner and Lutz (2000). The boundary \(\partial \Delta ^{d+1}\) of the \((d+1)\)-simplex is a d-dimensional simplicial complex with \(d+2\) facets. A bistellar move is a local modification of a combinatorial d-manifold K in which any subcomplex of K isomorphic to the star of a face in \(\partial \Delta ^{d+1}\) is replaced by its complementary facets.

To be precise, let \(\sigma \) be an i-face of K which is contained in exactly \(d-i+1\) facets \(\tau _1,\dots ,\tau _{d-i+1}\) such that these facets cover exactly \(d+2\) vertices. Identifying those \(d+2\) vertices with the vertices of \(\Delta ^{d+1}\) yields \((d+2)-(d-i+1)=i+1\) complementary facets \(\tau _{d-i+2},\dots ,\tau _{d+2}\) in the boundary \(\partial \Delta ^{d+1}\). Replacing \(\tau _1,\dots ,\tau _{d-i+1}\) by \(\tau _{d-i+2},\dots ,\tau _{d+2}\) in K is a candidate bistellar move of dimension \(d-i\), or a candidate \((d{-}i)\)-move for short. Let \(\sigma ' = \cap _{j=d-i+2}^{d+2}\tau _j\) be the complementary face to \(\sigma \), where \(\sigma '\) is of dimension \(d-i\). If \(\sigma '\) is not already contained in K, the move is proper. Applying an i-dimensional proper bistellar move reduces the f-vector of K if and only if \(i>d/2\).

Two simplicial complexes are bistellarly equivalent if one is obtained from the other by a finite sequence of (proper) bistellar moves. The following result is essential for the third step in the heuristic.

Theorem 6

(Pachner 1987) A d-dimensional simplicial complex is a PL d-sphere if and only if it is bistellarly equivalent to \(\partial \Delta ^{d+1}\).

This is closely related to Theorem 2 in the following sense: Adiprasito and Izmestiev (2015) showed that iterated barycentric subdivisions make any PL sphere polytopal; and barycentric subdivisions can be expressed as sequences of stellar subdivisions (which, by Theorem 6, are connected via sequences of bistellar moves). Moreover, barycentric subdivisions of polytopal spheres are polytopal, and polytopal spheres admit spherical discrete Morse vectors.

We now discuss the polymake implementation of the simulated annealing strategy from Björner and Lutz (2000). The function bistellar_simplification randomly applies bistellar moves to an input of type SimplicialComplex (required to be a combinatorial d-manifold) with the goal to lower the f-vector as much as possible. In this way the algorithm prefers moves that reduce the f-vector; this is called “cooling”. It lies in the nature of the sphere recognition problem that we may end up in a local minimum, i.e., when there are no moves to further lower the f-vector. At that point, we deliberately make moves that increase the f-vector for some number of rounds (this is called “heating”). Then we cool again, hoping that this will help jiggle us out of that local minimum. For \(0\le i < \frac{d}{2}\) a corresponding i-move will increase the f-vector of a triangulation, while in even dimensions the f-vector is not altered by \(\frac{d}{2}\)-moves. So in a heating phase we would add vertices via 0-moves, edges via 1-moves, etc., and “randomize” the triangulation by performing a (possibly large) number of \(\frac{d}{2}\)-moves, before returning to cooling. As with all simulated annealing approaches, adjusting the parameters for the annealing relies on experimentation. For example, we initially may not add vertices in the heating phases via 0-moves, as this might successively increase the size of the intermediate complexes, but in case we remain in a local minimum, first adding some percentage of vertices before performing 1-moves etc. helps in some cases.

Our procedure determines all candidates for bistellar moves of K and sorts them by descending dimension. During a cooling period we first pick random d-moves, if one exists. Otherwise, we pick random candidate \((d-1)\)-moves until we find one which is proper. If this does not exist either, we continue further to dimensions \(d-2\), \(d-3\), etc., down to dimension \(\lfloor d/2\rfloor +1\). Any proper move found in this way is applied immediately. Note that a bistellar move is a local operation, which is why we refrain from copying the entire complex when we apply a bistellar move. Instead, we perform the operation in place and store the reverse move in a list such that it can be undone later. Cooling continues until we get stuck with a lexicographically locally minimal f-vector. This ordering of the f-vectors is imposed indirectly by preferring higher-dimensional moves.

During a heating period, the story is slightly different. One heating strategy is to choose the dimension of the heating move at random with respect to a heuristically determined distribution. That distribution is encoded as a heat vector \((h_0,\dots ,h_{\lfloor d/2\rfloor })\) of integers, and we set \(h:=h_0+\dots +h_{\lfloor d/2\rfloor }\). This means, in each round of the heating period we pick the dimension k with probability \(h_k/h\), and in that dimension we pick a random proper bistellar move. For example, the default heat vector in polymake for \(d=4\) is (10, 10, 1). This generalizes to the default heat vector \((10,10,\dots ,10,1)\) in higher (even) dimensions d, while for odd d the pivot dimension k is picked uniformly at random.

Various other parameters control the precise heating behavior; and some of them are adjusted dynamically. For instance, it is useful to heat up for more rounds if the complex is larger. Sometimes it pays off to experiment with several distributions or to use other heating schemes. E.g., as mentioned above, first add some (percentage of) or no vertices via 0-moves, then edges via 1-moves, etc.

Remark 7

As a speed-up for large input triangulations, we can first apply edge contractions (with admissible edges for a contraction chosen at random) as long as possible. As experienced for 3-manifold triangulations (Csorba and Lutz 2006), this eventually leads to a saturation with many edges that block further contractions. Once there is no remaining admissible edge for a contraction, we run bistellar flips to reduce the number of edges and then continue with edge contractions again. Edge contractions are useful only in an initial phase. Once a local minimum is reached for the size of the triangulation, then bistellar flips are employed to leave the local minimum.

2.5 (4) Fundamental group

A non-trivial fundamental group \(\pi _1(K)\) is a certificate for not being a sphere (PL or otherwise). Conversely, there is the solution to the (PL) Poincaré Conjecture in dimensions other than four.

Theorem 8

(Smale 1961; Perelman 2002) Let K be a simply connected combinatorial d-manifold, \(d\ne 4\), with spherical homology. Then K is a PL sphere.

Freedman proved that a simply connected 4-manifold with trivial intersection form is homeomorphic to the 4-sphere Freedman (1982). But his result does not say whether this also holds in the PL category. In fact, it is a major open problem whether or not “exotic” 4-spheres exist.

In Seifert and Threlfall (1934, Chapter 7) Seifert and Threlfall describe how to obtain a finite presentation of \(\pi _1(K)\) from any spanning tree in the 1-skeleton (with the remaining edges as generators) and all the 2-faces (as relators). However, checking if a finitely presented group is trivial is known to be undecidable (Novikov 1955). Discussing heuristic approaches to simplifying group presentations is beyond the scope of this paper. In practice we rely on GAP (GAP 2019) which employs Tietze transformations.

Algorithm 1 displays our strategy in a concise form; for computational results see Sects. 3 and 4 below. Notice that the ordering of the Steps (1) through (3) is arbitrary, while the “YES” answer in Step (4) is inconclusive without Step (1). Yet, there is a benefit from combining Steps (1) and (2); cf. Remark 5. In practice, for a complex K we suspect to be not a sphere, we would start with (1), while if we think that K is a PL sphere, we first try (2) as a fast routine. If we are not successful with (2), we switch to (3), which is slower but can still recognize spheres that do not have perfect discrete Morse functions; see the discussion in Sect. 4.

If, in the case \(d\ne 4\), Step (1) gives us a spherical homology vector and Step (4) a trivial presentation of the fundamental group, then the overall output is “YES”, by Theorem 8. In the case of spherical homology a presentation of the fundamental group with only one generator is not possible, but balanced presentations of the trivial group with two generators and two relators can already be hard to detect; see Sect. 4.2.

Clearly, when our method gives up with “UNDECIDED” this does not need to be the end of the story. For instance, in the 3-dimensional case we can feed the data into the 3-sphere recognition procedure of Regina (Burton et al. 1999–2021). This features a variation of the exact algorithm of Rubinstein (1995) and Thompson (1994), where, for instance, the crushing procedure (the key step of the algorithm) is dramatically simplified (Burton 2014). In this way Regina can provide certificates for K not being spherical based on normal surface theory. However, we are not aware of a single triangulation of the 3-sphere for which our procedure fails.

3 Experiments and runtime comparisons

To find challenging input for Algorithm 1 is not entirely trivial. Most explicit constructions of (standard) PL spheres found in the literature are rather small and can be recognized instantaneously. All timings were taken on an AMD Phenom(tm) II X6 1090T Processor CPU (3.2 GHz, 6422 bogomips) and 8 GB RAM with openSUSE Leap 15.0 (Linux 5.1.9-5).

3.1 Recognizing random 3-spheres with polymake and Regina

A natural class of PL d-spheres are the boundaries of \((d+1)\)-polytopes obtained as the convex hulls of n points chosen uniformly at random on the unit d-sphere in \(\mathbb {R}^{d+1}\). These have been studied, e.g., in the context of the average case analysis of the simplex method of linear programming (Borgwardt 1987). Such examples can be generated with the rand_sphere command of polymake. Table 1 lists polymake and Regina experiments on 3-spheres with up to 100,000 vertices. For more than 15,000 vertices the convex hull computation (necessary only to construct the input) becomes a bottleneck, which is why for the larger examples (marked “*”) we used connected sums of smaller random spheres.

Table 1 Running times (in seconds) on random polytopal 3-spheres on n vertices
Full size table

In all cases, the spheres were successfully recognized by each method. However, we truncated the time spent on each input to about one CPU hour, such that longer running times are omitted. The fastest method is polymake ’s random search for a spherical discrete Morse function; cf. Step (2) of Algorithm 1. Nearly competitive is polymake ’s procedure of applying edge contractions, combined with random bistellar moves; cf. Step (3) of Algorithm 1 and Remark 7.

Usually, Regina takes 1-vertex pseudo-simplicial triangulations as input, but can also handle (abstract) simplicial complexes. In the latter case, contracting a spanning tree in the 1-skeleton yields a 1-vertex pseudo-simplicial triangulation. Conversely, the second barycentric subdivision of a pseudo-simplicial complex is a simplicial complex. In this sense these two encodings of combinatorial manifolds are similar.

Regina ’s recognition algorithm isThreeSphere runs, as a preprocessing step, the program IntelligentSimplify that transforms the complex into a 1-vertex triangulation and uses bistellar moves to further reduce it, similar to Step (3) of our Algorithm 1. Afterwards the 3-sphere recognition procedure is employed. In this way, Regina is able to also find certificates for non-sphericity—which polymake is incapable of, beyond checking the homology. We should also point out that IntelligentSimplify is a heuristic designed to be an out-of-the-box first attempt to simplify a triangulation with a polynomial running time, and Regina ’s bistellar move interface is meant to be interactive. This means that with a bit of work to build a custom-made simplification routine, the times in Regina could probably be improved.

The largest simplicial complex in Table 1, with 100,000 vertices, has 673,274 tetrahedra. The largest one successfully handled within 1 h by Regina has 15,000 vertices and 101,088 tetrahedra.

Each row of the Tables 1, 3 and 4 corresponds to a single instance only. However, it is known that there is little variation of, e.g., the number of facets of the convex hull of random points on the unit sphere; cf. Reitzner (2005, Sec. 4). This can also be observed experimentally; cf. (Assarf et al. 2017, Sec. 3.5 and Fig. 6) for a closely related setup. Note that, in fixed dimension d, the expected number of facets of a random simplicial \((d+1)\)-polytope depends linearly on the number of vertices (Borgwardt 1987).

Table 2 Census of 4-manifolds with up to six pentachora
Full size table

3.2 Processing a census of 4-manifolds

We ran our heuristics on a census of 4-manifolds provided by Regina (Burton et al. 1999–2021). These 4-manifolds are encoded as pseudo-simplicial complexes comprising up to six maximal cells; the 4-simplices are called pentachora in Regina for they have five facets (and five vertices). In combined form, the tricensus command of Regina generates possible facet pairings, then for each such pairing determines all possible gluing permutations, and (on the fly) reduces all gluings to isomorphism signatures that uniquely encode triangulations up to combinatorial isomorphism (Burton 2013). For the examples with two, four, and six pentachora there are 3, 26, and 639 facet pairings. And this yields 8, 784, and 440,494 resulting combinatorial types, respectively. Note that there are no facet pairings for an odd number of maximal cells in even dimensions. We let Regina expand each of these into a (proper abstract) simplicial complex via the second barycentric subdivision and pass it on to polymake. The simplicial complexes resulting from six pentachora have around 4,600 vertices and 86,400 facets.

Using our heuristic we found the results in Table 2. Each positive or negative certificate was obtained in less than 4 min and in 90 s on the average. In all the cases the positive certificates arise from discrete Morse functions, while the negative certificates are provided by non-spherical homology. Taking row sums in Table 2 we summarize our findings as follows.

Theorem 9

Among the 441, 286 combinatorial types of combinatorial 4-manifolds arising from up to six pentachora  \(91.5\%\) are spheres, and  \(8.0\%\) are non-spheres.

Thus, our success rate is \(99.5\%\), with our heuristic failing on only 1,954 of these combinatorial 4-manifolds. All of these have spherical homology; to determine whether they are standard PL 4-spheres or proper combinatorial homology 4-spheres (combinatorial 4-manifolds with spherical homology, but not PL homeomorphic to the standard PL 4-sphere) is an interesting question, yet beyond the scope of this article. Our classification, as a list of Regina ’s isomorphism signatures of the complexes can be found at Lofano (2021).

3.3 Higher-dimensional random spheres

polymake can easily recognize random polytopal spheres with up to 10,000 vertices in dimension four, 1,000 vertices in dimension five, and 500 vertices in dimension six; cf. Tables 3 and 4. Again the input is constructed via uniform random sampling on the unit sphere and taking convex hulls.

Regina provides no heuristic for sphere recognition in dimension four or beyond. Yet, Regina can simplify a given triangulation of a 4-dimensional combinatorial manifold via contractions and bistellar moves, returning a smaller pseudo-simplicial complex. It is not immediate how to check for sphericity from that output. That implementation is deterministic; thus in each run on a fixed input it gives the same output. The running times are given in the penultimate column of Table 3. The last column contains the number of simplices remaining after simplification.

Table 3 Running times (in seconds) on random polytopal 4-spheres on n vertices
Full size table
Table 4 Running times (in seconds) on random polytopal 5- and 6-spheres
Full size table

3.4 A collapsible 5-manifold which is not a ball

We consider the 5-dimensional simplicial complex C with face vector \(f(C)=(5013,72300,290944,495912,383136,110880)\) constructed in (Adiprasito et al. 2017, Sec. 4); there C is called contractible_non_5_ball. This is the first explicit example of a non-PL triangulation of a collapsible (and thus contractible) 5-manifold, other than the 5-ball. By construction, C is a manifold with boundary. To check the remaining topological properties computationally poses an interesting challenge.

First, the perfect Morse vector (1, 0, 0, 0, 0, 0) for C was originally obtained in a single random discrete Morse vector search over 82 h with a GAP implementation. The current implementation in polymake produces the same result (in most runs) in only 9 s with the random-lex-first and random-lex-last strategies and in about 10 min with the random-random strategy. This certifies that C is collapsible.

Second, the boundary complex \(\partial C\) with face vector \(f(\partial C)\!=\!(5010,65520,212000,252480,100992)\) was investigated; it was called contractible_non_5_ball_boundary in Adiprasito et al. (2017). Checking all face links for spherical discrete Morse vectors confirmed that \(\partial C\) is a combinatorial 4-manifold. For each face link a single random try sufficed. In total, the recognition of all face links took about 7.5 h. Checking the homology reveals that \(\partial C\) is a homology 4-sphere. Finally, GAP identifies the fundamental group \(\pi _1(\partial C)\) as the binary icosahedral group.

4 Limitations

In the previous section we saw that many, even fairly large, simplicial spheres can be recognized easily, despite Theorem 1. Here we explore the limitations of our heuristic. The combination of our (positive and negative) experiments may serve as a description of a “horizon” within which we can hope for effective recognition results.

4.1 General remarks

We refrain from a detailed comparison of simplicial homology computations. However, standard implementations, such as CHomP (Mischaikow et al. 2012–2021), RedHom (Mrozek et al. 2014–2021), Perseus (Mischaikow and Nanda 2013), and polymake (Gawrilow et al. 1997), employ elimination schemes for computing the integer homology, which are equivalent to finding discrete Morse functions with few critical cells. In this sense, the horizon within which we can compute the simplicial homology is essentially the same as the horizon for the discrete Morse Step (2). There are more software systems to compute simplicial homology, but many, including, e.g., Dionysus and Dionysus 2 (Morozov 2017–2021) and PHAT (Bauer et al. 2017), are restricted to \(\mathbb {Z}_2\)-coefficients. Finding an optimal discrete Morse function is NP-hard; cf. Joswig and Pfetsch (2006); Lewiner et al. (2003). Recently Bauer and Rathod (2019) established that we may not even hope for polynomial approximability.

In the subsequent we will exhibit several scenarios in which finding a spherical discrete Morse function for a simplicial sphere may fail in practice. An obvious impediment is the lack of any spherical discrete Morse function. The smallest known example is an 18-vertex triangulation of \(\mathbb {S}^3\), constructed from a triple trefoil knot supported on three edges (Benedetti and Lutz 2013a).

In dimension three, the known recognition algorithms for the 3-sphere make use of normal surface theory. As a consequence of Theorem 8, a trivial fundamental group suffices to show that a 3-manifold is, in fact, the 3-sphere.

4.2 Akbulut–Kirby spheres

A family of standard PL \(\mathbb {S}^4\)-triangulations, \({{\,\mathrm{AK}\,}}(r)\), for \(r\ge 3\), of the Akbulut–Kirby spheres (Akbulut and Kirby 1985) has been constructed in Tsuruga and Lutz (2013). In fact, Akbulut and Kirby (1985) gave handlebody decompositions of a family of 4-manifolds, in the hope of obtaining exotic 4-spheres. Yet, later Gompf (1991) and Akbulut (2010) showed that these manifolds are PL homeomorphic to the standard 4-sphere \(\mathbb {S}^4\). We have

$$\begin{aligned} f({{\,\mathrm{AK}\,}}(r))= & {} (176 + 64 r,\, 2390 + 1120 r,\, 7820 + 3840 r, \\&\quad 9340 + 4640 r,\, 3736 + 1856 r). \end{aligned}$$

The triangulated Akbulut–Kirby spheres \({{\,\mathrm{AK}\,}}(r)\) so far constitute the single explicit family of simplicial spheres that we could not recognize easily by our heuristic. More precisely, Step (2) failed on all complexes \({{\,\mathrm{AK}\,}}(r)\) for all \(r\ge 3\). Step (3) worked for \(r=3\), but failed for \(r\ge 4\). Steps (2) and (3) are particularly relevant in dimension four—in all other dimensions, as a consequence of Theorem 8, they can conceptually be replaced by the combination of Steps (1) and (4). We do not know if the spheres \({{\,\mathrm{AK}\,}}(r)\) admit perfect discrete Morse vectors. E.g., the smallest one we found for \(r=5\) is (1, 2, 4, 2, 1), possibly reflecting that we used two winded up 1-handles in the construction.

What makes the PL 4-spheres \({{\,\mathrm{AK}\,}}(r)\) difficult to recognize is that the original handlebody decomposition of Akbulut and Kirby (1985) is based on the non-trivial presentation

$$\begin{aligned} G(r) \ = \ \langle \, x,y\,\mid \, xyx = yxy;\, x^r = y^{r-1}\,\rangle \end{aligned}$$

of the trivial group built in as the fundamental group, i.e., \(\pi _1(\mathbb {S}^4)\), of \({{\,\mathrm{AK}\,}}(r)\). In 100 out of 450 runs we found GAP to be able to identify \(\pi _1({{\,\mathrm{AK}\,}}(4))\) as the trivial group. However, in dimension four (and knowing that the homology is spherical) this only shows that the input is a topological 4-sphere, without yielding any information on the PL type. In one run with \(r=5\) we obtained the initial presentations G(r) as the output of GAP ’s simplification procedure. For \(r\ge 4\), none of the Steps (1)–(4) was conclusive to determine that \({{\,\mathrm{AK}\,}}(r)\) is the standard PL 4-sphere \(\mathbb {S}^4\).

4.3 Contractible but non-collapsible subcomplexes

A simplicial sphere K admits a perfect discrete Morse function (respectively vector) if and only if there is a facet \(\sigma \) of K such that \(K-\sigma \) is a collapsible ball. In this way a key difficulty in finding a perfect discrete Morse vector for \(K-\sigma \) stems from subcomplexes that are contractible, but not collapsible. The most prominent such example is the 2-dimensional dunce hat which can be obtained from a single triangle by identifying its three boundary edges in a non-coherent way (Zeeman 1964).

Crowley et al. (2003) showed that the 7-simplex with 8 vertices contains in its 2-skeleton an 8-vertex triangulation of the dunce hat onto which it collapses. (This result can easily be verified by running the random discrete Morse Step (2) on the 7-simplex, but not allowing the triangles of an 8-vertex triangulation of the dunce hat be used as free faces.) While the dunce hat has triangulations with 8 vertices (Benedetti and Lutz 2013b), every contractible complex with fewer vertices is collapsible (Bagchi and Datta 2005).

This leads us to our next experiment, where we compute random discrete Morse vectors for d-simplices, \(7 \le d \le 25\); cf. Table 5. For instance, in dimension seven every one of the \(10^{10}\) runs that we tried gave a perfect discrete Morse function, i.e., a collapsing sequence. With increasing dimension that success rate drops slowly until \(d=20\), where we get stuck with a discrete Morse vector which is not perfect in \(1.3\%\) of all tries. Going to even higher dimensions shows a rapid decline of the probability to find a perfect discrete Morse vector. From this we conclude that dimension 25 marks a “horizon” for Step (2), even for a single simplex. Note also that the implementation of the algorithm requires to store the entire Hasse diagram, which is memory expensive; e.g., the Hasse diagram of the 25-simplex needs around 200 GB of RAM. The data structure underlying the Hasse diagram is optimized for speed with respect to enumerating all faces of a simplicial complex (and the covering relations of the inclusion poset) from the facets. That algorithm, which is linear in the size of the output, builds on a very general method of Ganter (1987) for closure systems, and the implementation details in polymake are discussed in Hampe et al. (2019).

Table 5 Collapsing the d-simplex
Full size table

It is instructive to look at the subcomplexes which arise from non-perfect discrete Morse vectors in our experiments for a d-simplex. For \(d=8\), we found four examples of 2-dimensional contractible and non-collapsible complexes on nine vertices which we call \({{\,\mathrm{D}\,}}\), \({{\,\mathrm{SB}\,}}^2_a\), \({{\,\mathrm{SB}\,}}^3_b\), and \({{\,\mathrm{SQ}\,}}\); cf. Fig. 3 and (the second line of) Table 6. The first one, \({{\,\mathrm{D}\,}}\), is a triangulation of the dunce hat. The following concept is derived from scrutinizing the complexes in Fig. 3.

Fig. 3
figure 3

The simplicial complexes \({{\,\mathrm{D}\,}}\), \({{\,\mathrm{SB}\,}}^2_a\), \({{\,\mathrm{SB}\,}}^3_b\) and \({{\,\mathrm{SQ}\,}}\) (clockwise from top left)

Full size image

Definition 10

The k-bladed saw blade complex \({{\,\mathrm{SB}\,}}^k\) is the 2-dimensional CW complex obtained from a polygonal disk with 3k edges \(a_1,a_1^{-1},a_1, a_2,a_2^{-1},a_2,\dots , a_k,a_k^{-1},a_k\) by identifying \(a_1,a_1^{-1},a_1\) as well \(a_2,a_2^{-1},a_2\) and so on until \(a_k,a_k^{-1},a_k\), for \(k\ge 1\).

In particular, for \(k=1\) we obtain a triangle whose three edges are identified in the order \(a_1,a_1^{-1},a_1\), i.e., \({{\,\mathrm{SB}\,}}^1\) is the dunce hat. More generally, \({{\,\mathrm{SB}\,}}^k\) consists of k vertices, k edges, and a single disk; so the Euler characteristic equals one; cf. Fig. 4, which explains the name. We use notation like \({{\,\mathrm{SB}\,}}^k_x\) and \({{\,\mathrm{SB}\,}}^k_y\) to denote specific triangulations.

Fig. 4
figure 4

The sawblade complex \({{\,\mathrm{SB}\,}}^3\)

Full size image
Table 6 Discrete Morse vectors for \(10^9\) runs on the 8-simplex
Full size table

Theorem 11

The following holds for k-bladed saw blade complexes \({{\,\mathrm{SB}\,}}^k\):

  1. (i)

    The dunce hat \({{\,\mathrm{SB}\,}}^1\) can be triangulated with 8 vertices.

  2. (ii)

    \({{\,\mathrm{SB}\,}}^2\) can be triangulated with 9 vertices.

  3. (iii)

    \({{\,\mathrm{SB}\,}}^k\) can be triangulated with 3k vertices, for \(k\ge 3\).

  4. (iv)

    Any triangulation of a saw blade complex is contractible, but non-collapsible.

Proof

We first prove (iii): For \(k\ge 3\), assume that the identified boundary of the polygonal disk reads 1–2–1–2–3–2–3–4–3–4–...–(k–1)–k–1; see Fig. 4 in the case \(k=3\). In the interior of the disk we place a cycle with 2k vertices and connect the cycle vertices with the boundary cycle vertices. More precisely, we connect every other cycle vertex to the beginning vertex of a blade and its two neighbors, and we connect the remaining cycle vertices to the two middle vertices of a blade. Finally, the interior 2k-gon can be triangulated arbitrarily without additional vertices.

(ii)+(i): In the case of two blades, we start with 1–2–1–2–3–1–3–2–3–1 as the identified boundary cycle. The extra vertex, say 3, is needed to avoid unwanted additional identifications. In the interior we then place a hexagon and connect its vertices similar to before. For an 8-vertex triangulation of the dunce hat \({{\,\mathrm{SB}\,}}^1\) see Benedetti and Lutz (2013b).

(iv): The dunce hat \({{\,\mathrm{SB}\,}}^1\) is contractible, and none of its triangulations is collapsible (Zeeman 1964). For \(k\ge 2\), the saw blade complex \({{\,\mathrm{SB}\,}}^k\) has k vertices that (by labeling appropriately) appear in order 1–2–1–2–3–2–3–4–3–\(4-\dots -(k\)–1)–k–1 on the boundary of the original 3k-gon. We cut the 3k-gon along an interior arc into two polygonal disks with identifications on the boundary, 1–2–1–2–(interior arc) and the remainder 2–3–2–3–4–3–\(4-\dots -(k\)–1)–k–1–(interior arc). Both parts are contractible CW complexes (that retract to the paths 1–2 and 2–3–4–...–k–1) glued along a contractible subcomplex (the interior arc). We conclude that their union \({{\,\mathrm{SB}\,}}^k\) is contractible; cf. Hachimori (2008) for a similar decomposition of a contractible 2-complex.

Now we consider any triangulation K of the saw blade complex \({{\,\mathrm{SB}\,}}^k\), for \(k\ge 1\). There is no free 2-face as all edges in such a triangulation either have degree two or three. It follows that K is non-collapsible. \(\square \)

Saw blade triangulations with different numbers of blades are combinatorially non-isomorphic. These simplicial complexes and their quotients provide 2-dimensional contractible, but non-collapsible simplicial complexes on which we can get stuck when trying to randomly collapse a simplex. For instance, the simplicial complex \({{\,\mathrm{SQ}\,}}\) from Fig. 3 is obtained as a quotient from identifying two vertices in some triangulation of \({{\,\mathrm{SB}\,}}^3\). Similar examples and higher-dimensional analogs exist in abundance.

Table 7 Discrete Morse vectors for \(10^3\) runs on the 20-simplex
Full size table

Tables 6 and 7 give the actual discrete Morse vectors found for the 8-simplex and the 20-simplex, respectively. We observe that we can get stuck (i.e., run out of free faces) in different dimensions; see Lofano and Newman (2021) for an analysis of this phenomenon.

While in the case of the 8-simplex at most two extra critical cells are picked up (see Table 6), the discrete Morse vector \((1, 0, 6, 48, 182, 377, 657, 876, 801, 493, 170, 22, 0, 0, 0, 0, 0, 0, 0, 0, 0)\) for the 20-simplex in Table 7 contains 3, 632 extra critical cells. Thus, in higher dimensions, not only do we get stuck with non-collapsible, contractible subcomplexes more often, but when we get stuck, the resulting discrete Morse vectors will also be larger. This may be seen as empirical evidence for the non-approximability of perfect Morse function; cf. Bauer and Rathod (2019).

4.4 Iterated barycentric subdivisions

As already seen in Sect. 3, it is sometimes rather easy to find perfect discrete Morse vectors, even in fairly large simplicial complexes, provided that the complexes are nicely structured; cf. also Benedetti and Lutz (2014). Adiprasito and Izmestiev (2015) showed that sufficiently large iterated barycentric subdivisions of any PL sphere admit spherical discrete Morse functions. Yet, the average number of critical cells for random discrete Morse vectors grows exponentially with the number of barycentric subdivisions (Adiprasito et al. 2017). Here we try our sphere recognition heuristic on higher barycentric subdivisions of boundaries of simplices.

For the third barycentric subdivision \({{\,\mathrm{sd}\,}}^3\partial \Delta _4\) of the boundary of the 4-simplex with face vector (12600, 81720, 138240, 69120) the perfect discrete Morse vector (1, 0, 0, 1) was found in 994 out of 1000 runs of the random-lex-last version (cf. Adiprasito et al. (2017) and Sect. 2) of the random discrete Morse search; see Table 8. For \({{\,\mathrm{sd}\,}}^4\partial \Delta _4\) with face vector (301680, 1960560, 3317760, 1658880) the (same) perfect discrete Morse vector was found in only 844 out of 1000 runs. This suggests that the 4th barycentric subdivision is still within the “horizon” for computations with the version random-lex-last, while the 5th barycentric subdivision was too large to fit into the main memory of the machine we used for the experiments. The random-lex-first strategy behaved slightly better than random-lex-last; the strategy random-random was always successful.

Table 8 Discrete Morse vectors for iterated barycentric subdivisions of the 3-sphere \(\partial \Delta _4\)
Full size table

4.5 Other input

Except for the Akbulut–Kirby spheres all the examples studied so far arise from easy to understand procedures. Searching for entirely different triangulations of \(\mathbb {S}^3\), we started out with \(\partial \Delta _4\) with five vertices. Then we added 525 vertices via random 0-moves, followed by 50,000 random 1-moves, followed by another \(10^6\) rounds of random bistellar moves where we allowed both 1- and 2-moves. This resulted in a “random” triangulation of the 3-sphere with face vector \(f=(530, 50474, 99888, 49944)\), which we fed into our heuristic. The smallest discrete Morse vector found was (1, 2192, 2192, 1)—far away from the perfect vector (1, 0, 0, 1). This means that Step (2) fails on such input. Yet, applying bistellar moves again quickly gives back the initial \(\partial \Delta _4\). We also used GAP to actually find a trivial presentation for the fundamental group of the example, which took 16 h for the simplification. It could be interesting to further investigate this or similar classes of examples.