Exact Algorithms for Minimum Dilation Triangulation

We provide various new contributions, both to theory and practice.

1. 1.

   We present practically robust methods and implementations for computing an optimal MDT for benchmark sets with up to 30000 points in reasonable time on commodity hardware, based on new geometric insights into the structure of optimal edge sets. Previous methods only achieved results for up to 200 points (involving one computational routine of complexity $\mathrm{\Theta }(n^4)$ instead of our improved complexity of $O(n^2\log n)$), so we extend the range of practically solvable instances by a factor of 150.

2. 2.

   We develop scalable techniques for accurately evaluating many shortest-path queries that arise as large-scale sums of square roots, allowing us to certify exact optimal solutions. This differs from previous work, which relied on floating-point computations, without regard for errors resulting from numerical issues.

3. 3.

   We resolve an open problem from [20] by establishing a lower bound of $1.44116$ on the dilation of the regular $84$-gon (and thus for arbitrary point sets). This improves the previous worst-case lower bound of $1.4308$ from the regular $23$-gon and greatly reduces the remaining gap to the upper bound of $1.4482$ from [43]. In the process, we provide optimal solutions for regular $n$-gons up to $n=100$.

The complexity of finding the MDT is unknown [24]. [28] prove that finding the minimum dilation graph with a limited number of edges is $\mathrm{𝖭𝖯}$-hard. [12] show that finding a spanning tree of given dilation is also $\mathrm{𝖭𝖯}$-hard. Kozma [33] proves $\mathrm{𝖭𝖯}$-hardness for minimizing the expected distance between random points in a triangulation, with edge weights instead of Euclidean distances. For surveys, see Eppstein [24] until 2000 and [38] for more recent results.

On the practical side, the authors of [44] study the problem of finding sparse low dilation graphs for large point sets in the plane. The authors of [9] present an approximation algorithm for the related problem of improving a given graph with a budget of $k$ edges such that the dilation is minimized. Regarding the MDT, all practical approaches in the literature are based on fixed-precision arithmetic. Klein [31] used an enumeration algorithm to find an optimal MDT for up to 10 points. The authors of [19] present heuristics for the MDT and evaluate their performance on instances with up to 200 points. Instances with up to 70 points were solved by the authors of [8] using integer linear programming techniques. In their approach, the edge elimination strategy from Knauer and Mulzer [32] was used to eliminate edges from the complete graph. Recently, Sattari and Izadi [41] presented an exact algorithm based on branch and bound that was evaluated on instances with up to 200 points.

The MDT is closely related to finding a plane $t$-spanner; see [34] for an overview. Chew [13] first proved an upper bound of $\sqrt{10}$ on plane $t$-spanners in the $L_1$-metric, which he later improved [14] to $2$ for the triangular-distance Delaunay graph in the plane. [5] proved that any convex point set admits a plane $1.88$-spanner. For a centrally symmetric convex point set containing $n$ points, Sattari and Izadi [42] give an upper bound of $n/2\sin (\pi /n)$. The best known upper bound on the dilation of arbitrary point sets is by Xia [48], who established and upper bound of 1.998 for the dilation of the Delaunay triangulation.

Mulzer [35] studied the MDT for the set of vertices of a regular $n$-gon and proved an upper bound of $1.48586$. Amarnadh and Mitra [2] improved this bound to $1.48454$ for any point set that lies on the boundary of a circle. Sattari and Izadi [43] again improved the bound to $1.4482$. Dumitrescu and Ghosh [20] show that any triangulation of a regular $23$-gon has dilation at least $1.4308$, improving upon the bound of $1.4161$ by Mulzer [35] and answering a question posed by Bose and Smit [7] as well as Kanj [30]. Dumitrescu and Ghosh [20] also computed dilations of regular $22$-gon and regular $24$-gon.

Our supergraph is based on the well-known ellipse property (used in [8, 28, 32]) that all edges of a triangulation $T$ with dilation below $\rho$ must satisfy.

Recall that the set of points that have the same sum of distances $\kappa$ to two points $\mathrm{\ell },r$ with $\kappa \ge d(\mathrm{\ell },r)$ is an ellipse with $\mathrm{\ell },r$ as its focal points (or foci); thus, all paths between $\mathrm{\ell }$ and $r$ with length less than $\kappa =\rho \cdot d(\mathrm{\ell },r)$ must lie strictly inside this ellipse.

If a pair of points $s,t$ does not have the ellipse property, then there is a pair of points $\mathrm{\ell },r$ such that $st$ cuts through all paths between $\mathrm{\ell }$ and $r$ that could have dilation less than $\rho$; see Figure 2. We call such a pair of points $\mathrm{\ell },r$ an exclusion certificate for $s,t$.

Preliminary experiments showed that a brute-force check of each of the $\mathrm{\Theta }(n^4)$ pairs of potential edges is impractical for large $n$, dominating the overall runtime time even for early versions of our solution approach.

We therefore developed a more efficient scheme to enumerate a superset of the edges that satisfy the ellipse property. This scheme performs a filtered incremental nearest-neighbor search from each point $p\in P$ to identify all candidate edges of the form $pt$, i.e., looking for all possible neighbors $t$ of $p$. This search is efficiently implemented on a quadtree containing all points. While enumerating candidates, we construct so-called dead sectors, i.e., regions of the plane that cannot contain possible neighbors of $p$. We exclude all points that lie in dead sectors; we also use dead sectors to prune entire nodes of the quadtree and to terminate the search early if it has become clear that all remaining points must be in a dead sector. This usually avoids considering most points as potential neighbors of $p$ individually. The algorithm has a runtime of $O(nk\log n)$, where $k$ is the average number of points and quadtree vertices considered individually from each point. At worst, this can be $O(n^2\log n)$; in practice, $k$ is often much lower than $n$. A related, slightly less complex enumeration algorithm applied to minimum-weight triangulations by Haas [29] scales to ${10}^8$ points. In the following, we describe the components of the enumeration scheme in more detail.

We begin by giving a definition of the dead sectors we use.

Depending on $p$, $\rho$, $\mathrm{\ell }$ and $r$, $𝒟{𝒮}_{\rho }(p,\mathrm{\ell },r)$ is either empty (if $p$ is in the ellipse) or it is bounded by two rays and an elliptic arc; see Figure 3. In that case, it can also be interpreted as an elliptic arc and an interval of polar angles around $p$.

During our enumeration we construct many dead sectors, the union of which can become quite complex, making it cumbersome and inefficient to work with directly. We instead chose to simplify the shape of our dead sectors, giving up a small fraction of excluded area in exchange for a simple and efficient representation. We replace the elliptic arc by a single disk centered on $p$, whose radius is at least the maximum distance from $p$ to any point on the elliptic arc. We can hence represent each non-empty simplified dead sector by a polar angle interval around $p$ and a single radius called activation distance $A_{\rho }(p,\mathrm{\ell },r)$; see Figure 3.

We initially attempted to compute fairly precise upper bounds on the activation distance; however, due to computational and numerical issues described in the full version, we decided to use a more robust and efficient approach. Instead of using the elliptic arc, we compute an upper bound ${\tilde{A}}_{\rho }(p,\mathrm{\ell },r)$ using a minimal rectangle $B(\mathrm{\ell },r,\rho )$ containing the ellipse with sides are parallel and perpendicular to $\mathrm{\ell }r$. We only need to check the extreme points of $B(\mathrm{\ell },r,\rho )$ and its intersections with the rays $L(p,\mathrm{\ell },r)$ and $R(p,\mathrm{\ell },r)$ to compute an upper bound ${\tilde{A}}_{\rho }(p,\mathrm{\ell },r)$; see Figure 3.

In the worst case, these simplifications may lead to additional candidate edges. While this cannot make the resulting graph exclude any edges that satisfy the ellipse property, it is still undesirable; we use additional checks for further reductions later on.

We use a quadtree containing all points in $P$ for the filtered incremental search. The points are stored in a contiguous array $A_P$ outside the tree. Each quadtree node $v$ is associated with a contiguous subrange of $A_P$ represented by two pointers. This subrange contains the points in the subtree ${𝒯}_v$ rooted at $v$. Each node also has a bounding box that contains all points in ${𝒯}_v$. Each interior node has precisely four children; each leaf node contains at most a small constant number of points. To allow the contiguity of the subranges, points are reordered during tree construction. This has the added benefit of spatially sorting the points, improving the probability that geometrically close points are near each other in memory [29].

One can think of the filtered incremental search from $p$ as a process of continuously growing a disk centered at $p$ starting with a radius of $0$. As in the sweep-line paradigm, one encounters different types of events at discrete disk radii. We primarily encounter events when the disk first touches a point of $P$ or the bounding box of a quadtree node.

Observe that, during a search from a point $p$, the dead sectors have two states: either their activation distance is not yet reached, in which case they are inactive and do not exclude any points, or they are active and exclude all points in a certain polar angle interval around $p$. We therefore also introduce events when the disk radius reaches the activation distance of a dead sector. This enables efficient management of active dead sectors as a set of polar angle intervals around $p$; see the full version for a discussion of the ensuing numerical issues.

When we first encounter a point $t$, we have to determine whether $t$ is in any dead sector by checking the active dead sector data structure. If it is not, we have to report it as potential neighbor of $p$, adding it to a set of points sorted by polar angle around $p$. Furthermore, to construct new dead sectors, we combine $t$ with $O(1)$ other points of $P$; which points we use is decided by a heuristic discussed in the full version.

When we encounter a quadtree node $v$, we have to determine whether $v$’s bounding box is fully contained in the union of all dead sectors and can thus be pruned; we thus again check the active dead sectors. Otherwise, we have to take $v$’s children, or the points it contains if it is a leaf, into account; they are then considered as future events.

We initially compute the Delaunay triangulation of the point set $P$ and compute its dilation. We also optionally attempt to improve the dilation of the triangulation by a simple improvement heuristic. The heuristic is based on computing constrained Delaunay triangulations, greedily adding shortcut edges as constraints to reduce the length of the path currently defining the dilation. We then use the resulting dilation $\rho$ as bound for the enumeration process outlined in the previous sections, enumerating only edges that could locally be in a triangulation with dilation strictly below $\rho$.

After the initial enumeration process is complete, we are left with a set of possible edges and can safely ignore all other edges. We postprocess these as follows. For each possible edge $pq$, we compute the set of possible edges intersecting it. We need this information later on to model the problem of finding a triangulation on the set of possible edges. For each pair $st$, $\mathrm{\ell }r$ of intersecting edges, we explicitly check whether either pair is an exclusion certificate for the other. In many cases, this postprocessing gets us very close to the edge set that would be obtained by the trivial $\mathrm{\Theta }(n^4)$ edge candidate enumeration algorithm; see the experiment section for details. We mark each edge that does not have intersecting possible edges as certain; certain edges must be part of any triangulation with dilation less than $\rho$.

We also compute a dilation threshold $\vartheta (st)$ for each possible edge $st$. Let $I(st)$ be the set of possible edges intersecting $st$. For each $\mathrm{\ell }r\in I(st)$, we can compute a lower bound on the dilation of the shortest path connecting $\mathrm{\ell }$ to $r$, should $st$ be present, as

see Figure 4.

The dilation threshold of $st$ is then $\vartheta (st)={\max }_{\mathrm{\ell }r\in I(st)}{\rho }_{st}(\mathrm{\ell }r)$.

This allows us to quickly exclude further edges if we lower the dilation bound $\rho$ we aim for, without reenumerating edges or recomputing intersections.

We use the dilation thresholds to bound the minimum dilation. For each edge $st$, either $st$ or some edge crossing it must be in any triangulation. Thus, the minimum of all dilation thresholds among $\{st\}\cup I(st)$ is a lower bound on the minimum dilation. Combining all edges, we obtain the following lower bound:

We use interval arithmetic to compute a safe lower bound on this value in $O(1)$ time per intersection between two edges resulting from our enumeration.

Furthermore, by computing the points on the convex hull, we also know the number of edges $g$ of any triangulation. This also gives us a lower bound on the minimum dilation by considering the $g$ lowest dilation thresholds. We also use Kruskal’s algorithm to compute the lowest dilation threshold that admits a connected graph.

Both algorithms rely on the MDT supergraph mentioned in Section 3. As part of this computation, we also obtain an initial triangulation and its dilation, as well the intersecting possible edges $I(st)$ for each possible edge $st$. In both algorithms, we may gradually discover triangulations with lower dilation; these are used to exclude additional edges using the precomputed dilation thresholds $\vartheta (e)$. To keep track of the status of each edge, we insert all points and possible edges into a graph data structure we call triangulation supergraph. In this structure, we mark each edge as possible, impossible or certain. Initially, all enumerated edges are possible. If, at any point, all edges intersecting an edge $e$ become impossible, $e$ becomes certain. If an impossible edge becomes certain or vice versa, the graph does not contain a triangulation any longer. If this happens, we say we encounter an edge conflict.

Given a triangulation supergraph $G=(P,E)$, we model the problem of finding a triangulation on possible and certain edges using the following simple SAT formulation. Let $E_p\subseteq E$ be the set of non-impossible edges when the SAT formulation is constructed. For each edge $e\in E_p$, we have a variable $x_e$. We use the following clauses in our formulation.

Clauses (1) ensure crossing-freeness and clauses (2) ensure maximality. When an edge $e$ becomes certain, we add the clause $x_e$; similarly, when an edge becomes impossible, we add the clause $\neg x_e$. Both algorithms are based on this simple formulation; in the following, we describe how they use and modify it to find an MDT.

The following subproblem, which we call dilation path separation, arises in both our algorithms: Given a dilation bound $\rho$, a triangulation supergraph $G=(P,E)$ excluding only edges that cannot be in any triangulation with dilation less than $\rho$, a current triangulation $T$ and a pair of points $s,t\in P$ such that $|{\pi }_T(s,t)|\ge \rho \cdot d(s,t)$, find a clause $C$ that is (a) violated by $T$ and (b) satisfied by any triangulation $T^{\prime }$ with .

Based on the SAT formulation and the algorithm for the dilation path separation problem, IncMDT is simple. Given an initial triangulation $T$ with dilation $\rho$, we enumerate the set of candidate edges and construct a triangulation supergraph $G$ with bound $\rho$. We construct the initial SAT formula $M$ and solve it; if it is unsatisfiable, the initial triangulation is optimal. Otherwise, we repeat the following until the model becomes unsatisfiable or we encounter an edge conflict, keeping track of the best triangulation found, see Figure 5.

We extract the new triangulation $T^{\prime }$ from the SAT solver and compute the dilation ${\rho }^{\prime }$ and a pair $s,t$ of points realizing ${\rho }^{\prime }$. If ${\rho }^{\prime }$ is better than the best previously found dilation $\rho$, we update $\rho$ and mark all edges $e$ with $\vartheta (e)\ge {\rho }^{\prime }$ as impossible. We then set $T=T^{\prime }$ and solve the dilation path separation problem for $\rho$, $G$, $T$, $s$ and $t$. We add the resulting clause to $M$ and let the SAT solver find a new solution.

Preliminary experiments with IncMDT showed that we spend almost all runtime for computing dilations, even for instances for which we could rely exclusively on interval arithmetic, requiring no exact computations. For many instances, most iterations of IncMDT resulted in tiny improvements of the dilation. To reduce the number of iterations (and thus, dilations computed), we considered the binary search-based algorithm BinMDT.

Our experimental evaluation aims to answer the following questions.

1. Q1

   How does the edge enumeration algorithm from Section 3 compare to the brute force enumeration in terms of runtime and the number of edges eliminated? Can we solely rely on this algorithm or should we reconsider using the brute force enumeration?

2. Q2

   What is the quality of the initial solutions computed by the Delaunay triangulation and the constrained Delaunay triangulations from Section 3 compared to the optimal solution?

3. Q3

   How do our approaches compare to existing algorithms for the MDT w.r.t. runtime and solution quality? Can we solve instances with thousands of points to provable optimality?

4. Q4

   How do BinMDT and IncMDT compare? Which should be used for large instances?

To answer our questions, we collected and generated a large set of instances, consisting of instances from the following instance classes. In all cases, the coordinates of points in the instances are either integers or double precision floating-point numbers.

random-small

We include instances from the work of [8]. The $210$ instances have fixed sizes $n\in \{10,20,\mathrm{\dots },70\}$ and were generated by placing uniformly random points inside a $10\times 10$ square. For each size a total of $30$ instances were generated.

random

We include two sets of randomly generated instances. These encompass a set of instances with points with float coordinates chosen uniformly between $0$ and ${10}^3$, ranging from 50 to 10000 points. This resulted in a total of 800 instances.

public

We include instances from all well-known publicly available benchmark instance sets we could locate. These include point sets previously used in the CG:SHOP challenges [17, 16], TSPLIB instances [40], instances from a VLSI dataset¹¹1https://www.math.uwaterloo.ca/tsp/vlsi/index.html and point sets from the Salzburg Database of Polygonal Inputs [22]. In total, we collected 486 instances with up to 10000 points and an additional 38 with up to 30000 points.

We first compare the edge enumeration algorithm from Section 3 to the $\mathrm{\Theta }(n^4)$ brute force enumeration of all possible edges on the random instance set. Both preprocessing options include finding all pairwise intersections between the enumerated segments. Due to the $\mathrm{\Theta }(n^4)$ runtime, we only consider instances with up to 1000 points; see Figure 6.

The $\mathrm{\Theta }(n^4)$ algorithm precisely identifies all edges that have the ellipse property; it can hence only eliminate more edges than the edge enumeration algorithm from Section 3. However, for our test instances, the number of edges that can be eliminated by our approach is almost identical to the $\mathrm{\Theta }(n^4)$ algorithm; see Figure 6 for an example. The runtime makes the $\mathrm{\Theta }(n^4)$ algorithm infeasible for larger instances; it takes more than 250s for instances with only 1000 points, compared to for our more efficient approach.

Better initial solutions can reduce the number of candidate edges further than bad ones and thus affect the runtime of the overall algorithm. Delaunay triangulations are a natural choice for the initial solution, but we suspect that they can be improved by adding shortcut edges. Figure 7 shows that the relative dilation gap between the Delaunay triangulation and the optimal solution for the random instance set. It can be seen that the introduction of shortcut edges can indeed reduce the gap to the MDT to around 1.5%.

We compare our approaches to two exact state-of-the-art algorithms for the MDT. Note that both of these use floating-point arithmetic and are not guaranteed to find the optimal solution. However, we can confirm that all previous solutions are within a small relative error of our optimal solution. The first approach is an integer programming (IP) approach by [8] that used the commercial software CPLEX to solve the MDT. The second comparison is with the most recent exact algorithm (BnB) from Sattari and Izadi [41]. For both approaches the source code is no longer available, so we cannot compare our results on the same hardware. Also for BnB [41] the instance data and results are no longer available. For both approaches we therefore use the data the authors published in their papers. Note that the drastic runtime difference that we see in our experiments cannot be attributed to improved hardware alone.

For random-small, both IncMDT and BinMDT outperform the IP and BnB approach by a large margin (up to four orders of magnitude), see Figure 7. All instances were solvable in less than 0.1s. Additionally, Sattari and Izadi [41] provided results for TSPLIB [40] instances (part of our public instance set) with up to 200 points. Our approach is significantly faster than theirs, solving each of these instances to provable optimality in less than 1s instead of up to 1248s; a table with all instances and runtimes can be found in the full version.

We now compare IncMDT to BinMDT. Recall that BinMDT aims to reduce the time spent on dilation computations by reducing the number of dilation computations and merely sampling the dilation whenever that suffices. The first experiment was conducted on the random instances with up to 10000 points; this experiment confirms that BinMDT indeed achieves a significantly lower runtime, requires fewer dilation computations, and that sampling is sufficient in most cases, see Figure 8.

After we confirmed that BinMDT is superior, we conducted an additional experiment with two different initial solution strategies on the public instances up to 10000 points, see Figure 9. The first strategy uses the Delaunay triangulation as an initial solution; the second strategy uses the improved Delaunay triangulation with shortcut edges. The improved Delaunay triangulation significantly reduces the runtime of BinMDT for almost all instances. Detailed results for all public instances, as well as an additional experiment on the public instance set showing that BinMDT can solve instances with up to 30000 points in less than 17h to provable optimality, can be found in the full version.