Euclidean TSP, Motorcycle Graphs, and Other New Applications of Nearest-Neighbor Chains

TL;DR

This paper demonstrates the versatility of the nearest-neighbor chain algorithm by applying it to various geometric problems, achieving improved time complexities and novel solutions in computational geometry and matching models.

Contribution

It introduces new applications of the nearest-neighbor chain algorithm to geometric problems, including Euclidean TSP, motorcycle graphs, and stable matching variants, with improved algorithms and complexities.

Findings

Euclidean TSP greedy tour in O(n log n) time

Motorcycle graphs computed in O(n^{4/3+ε}) time

Linear-time 2-approximation for 1D geometric set cover

Abstract

We show new applications of the nearest-neighbor chain algorithm, a technique that originated in agglomerative hierarchical clustering. We apply it to a diverse class of geometric problems: we construct the greedy multi-fragment tour for Euclidean TSP in $O(n\log n)$ time in any fixed dimension and for Steiner TSP in planar graphs in $O(n\sqrt{n}\log n)$ time; we compute motorcycle graphs (which are a central part in straight skeleton algorithms) in $O(n^{4/3+\varepsilon})$ time for any $\varepsilon>0$; we introduce a narcissistic variant of the $k$-attribute stable matching model, and solve it in $O(n^{2-4/(k(1+\varepsilon)+2)})$ time; we give a linear-time $2$-approximation for a 1D geometric set cover problem with applications to radio station placement.

Tables

Table 1. Table 1 : Summary of best known results on the server cover problem. In the table: n 𝑛 n is the number of clients, m 𝑚 m is the number of servers, ε > 0 𝜀 0 \varepsilon>0 is an arbitrarily small constant, k > 1 𝑘 1 k>1 is an arbitrary integer parameter, γ > 0 𝛾 0 \gamma>0 is shown to be a constant, and T ​ ( n ) 𝑇 𝑛 T(n) is the cost of comparing the cost of two sets of disks in polynomial time, which requires comparing sums of square roots to compute exactly. ∗ The exact algorithm of  [ 42 ] is under the assumption that T ​ ( n ) 𝑇 𝑛 T(n) can be computed, which depends on the computational model.

Dim.	$𝜶$	Approximation ratio	Complexity
2D	$\alpha >1$	Exact	NP-hard [6]
2D	$\alpha =1$	Exact∗	$O\left({\left(n+m\right)}^{881}T\left(n+m\right)\right)$ [42]
		$1+\epsilon$	$O\left({\left(n+m\right)}^{881}T\left(n+m\right)\right)$ [42]
		$\left(1+6/k\right)$	$O\left(k^2{\left(nm\right)}^{\gamma +2}\right)$ [54]
1D	$\alpha \ge 1$	Exact	Polynomial (high complexity) [16]
1D	$\alpha =1$	Exact	$O\left({\left(n+m\right)}^2\right)$ [17]
		3	$O\left(n+m\right)$ [6]
		2	$O\left(m+n\log m\right)$ [6]
		2	$O\left(n+m\right)$ (this paper)

Equations

$$d(A\cup B,C)\ge \min (d(A,C),d(B,C)).$$

Full text

1 Introduction

The nearest-neighbor chain (NNC) technique is used for agglomerative hierarchical clustering, and has only seen one other use besides it. In this paper, we apply it to an assortment of new problems: multi-fragment TSP, straight skeletons, narcissistic $k$-attribute stable matching, and a server cover problem. These problems share a property with agglomerative hierarchical clustering, which we call global-local equivalence, and which is the key to using the NNC algorithm. First, we review the NNC algorithm in the context of clustering.

1.1 Prior work: NNC in hierarchical clustering

Given a set of points, the agglomerative hierarchical clustering problem is defined procedurally as follows: each point starts as a base cluster, and the two closest clusters are repeatedly merged until there is only one cluster left. This creates a hierarchy, where any two clusters are either nested or disjoint. A key component of hierarchical clustering is the function used to measure distances between clusters. Popular metrics include minimum distance (or single-linkage), maximum distance (or complete-linkage), and centroid distance.

We call two clusters mutually nearest neighbors (MNNs) if they are the nearest neighbor of each other. Consider this alternative, non-deterministic procedure: instead of repeatedly merging the two overall closest clusters, merge any pair of MNNs. Clearly, this may merge clusters in a different order. Nonetheless, if the cluster-distance metric satisfies a property called reducibility, this procedure results in the same hierarchy [21, 22, 60]. A cluster-distance metric $d(\cdot,\cdot)$ is reducible if for any clusters $A,B,C$: if $A$ and $B$ are MNNs, then

[TABLE]

In words, the new cluster $A\cup B$ resulting from merging $A$ and $B$ is not closer to other clusters than both $A$ and $B$ were. The relevance of this property is that, if, say, $C$ and $D$ are MNNs, merging $A$ and $B$ does not break that relationship. The net effect is that MNNs can be merged in any order and produce the same result. Many commonly used metrics are reducible, including minimum–, maximum–, and average–distance, but others such as centroid and median distance are not.

The NNC algorithm exploits this reducibility property, which was originally observed by Bruynooghe [22]. We briefly review the algorithm for hierarchical clustering, since we discuss it in detail later in the context of the new problems. For extra background on NNC for hierarchical clustering, see [61, 60]. The basic idea is to maintain a stack (called chain) of clusters. The first cluster is arbitrary. The chain is always extended with the nearest neighbor (NN) of the current cluster at the top of the chain. Note that the distance between clusters in the chain keeps decreasing, so (with an appropriate tie breaking rule) no repeated clusters or “cycles” occur, and the chain inevitably reaches a pair of MNNs. At this point, the MNNs are merged and removed from the chain. Crucially, after a merge happens, the rest of the chain is not discarded. Due to reducibility, every cluster in the chain still points to its NN, so the chain is still valid. The process continues from the new top of the chain.

The algorithm is efficient because each cluster is added to the chain only once, since it stays there until it is merged with another cluster. As we will see in detail for other problems, this bounds the number of iterations to be linear on the input size, with the cost of each iteration dominated by a NN computation.

1.2 Our contributions

Our key observation is that this equivalence between merging closest pairs and MNNs is not unique to hierarchical clustering. The problems in this paper, even though they are not about clustering, exhibit an analogous phenomenon, for which we coin the term global-local equivalence. The main thesis of this paper is that NNC is an efficient algorithm for problems with global-local equivalence, which includes many more problems than hierarchical clustering.

Recently, the NNC algorithm was used for the first time outside of the domain of hierarchical clustering [37, 38]. It was used in a stable matching problem where the two sets to be matched are point sets in a metric space, and each agent in one set ranks the agents in the other set by distance, with closer points being preferred. In this setting, there is a form of global-local equivalence: the stable matching is unique, and it can be obtained in two ways: by repeatedly matching the closest pair (from different sets), or by repeatedly matching MNNs. They used the NNC algorithm to solve the problem efficiently.

In this paper, we consider global-local equivalence in the context of the new problems, and give NNC-type algorithms for them. We summarize the computational results here. See each section for extended background on the corresponding problems.

Multi-fragment TSP.

A classic heuristic for the Euclidean Traveling Salesman Problem is the multi-fragment algorithm. While not having strong approximation guarantees, experimental results show that it performs better than other heuristics, particularly in geometric instances [34, 51, 57, 59, 13, 12]. We do not know of any subquadratic algorithm to compute the tour produced by this heuristic, which we call the multi-fragment tour. We give a $O(n\log n)$-time algorithm for computing the multi-fragment tour of a point set in any fixed dimension and using any $L_{p}$ metric. We also consider the Steiner TSP problem in a graph-theoretical framework [29], where we give a $O(n\sqrt{n}+k\sqrt{n}\log n)$-time algorithm for finding the multi-fragment tour through a subset of $k$ nodes in planar graphs and, more generally, graph families with $O(\sqrt{n})$-size separators.

Straight skeletons and motorcycle graphs.

The fastest algorithms for computing straight skeletons consist of two phases, neither of which dominates the other [26]. The first phase is a motorcycle graph computation. The best currently known algorithm for motorcycle graphs runs in $O(P(n)+n(T(n)+\log n)\log n)$ time, where $P(n)$ and $T(n)$ are the preprocessing time and operation time (maximum between query and update) of a dynamic ray-shooting data structure for curtains in $\mathbb{R}^{3}$ [67]. We improve this to $O(P(n)+nT(n))$. Using the structure from [2], both algorithms run in $O(n^{4/3+\varepsilon})$ for any $\varepsilon>0$, but if both use the same $\varepsilon$ in the data structure, ours is faster by a $O(\log n)$ factor.

Narcissistic $k$-attribute stable matching.

Given that $O(n^{2})$ is optimal for general stable matching instances, it is interesting to study restricted models. We introduce a narcissistic variant of the $k$-attribute model [15] and give a subquadratic, $O(n^{2-4/(k(1+\varepsilon)+2)})$-time algorithm for it, for any $\varepsilon>0$.

Server cover.

We give a linear-time $2$-approximation for a one-dimensional version of a server coverage problem: given the locations of $n$ clients and $m$ servers, which can be seen as houses and telecommunication towers, the goal is to assign a “signal strength” to each communication tower so that they reach all the houses, minimizing the cost of transmitting the signals. This improves upon the $O(m+n\log m)$-time algorithm by Alt et al. [6] with the same approximation ratio.

Paper organization.

Section 2 introduces a new data structure, which we call the soft nearest-neighbor data structure. Section 3 solves multi-fragment Euclidean TSP with a variant of NNC that uses this structure. Sections 4, 5, and 6 are on motorcycle graphs, narcissistic $k$-attribute stable matching, and server cover, respectively.

Thus, Section 3 relies on Section 2, but the other sections are independent, self-contained, and in no particular order. We give concluding remarks in Section 7.

2 The Soft Nearest-Neighbor Data Structure

Throughout this section, we consider points in $\mathbb{R}^{\delta}$, for some fixed dimension $\delta$, and distances measured under any $L_{p}$ metric $d(\cdot,\cdot)$. We begin with a formal definition of the structure and the main result of this section.

Definition 2.1** (Dynamic soft nearest-neighbor data structure).**

Maintain a dynamic set of points, $P$, subject to insertions, deletions, and soft nearest-neighbor queries: given a query point $q$, return either of the following:

• •

The nearest neighbor of $q$ in $P$: $p^{*}=\operatorname*{arg\,min}_{p\in P}d(q,p)$.

• •

A pair of points $p,p^{\prime}$ in $P$ satisfying $d(p,p^{\prime})<d(q,p^{*})$.

Theorem 2.2**.**

In any fixed dimension, and for any $L_{p}$ metric, there is a dynamic soft nearest-neighbor data structure that maintains a set of $n$ points with $O(n\log n)$ preprocessing time and $O(\log n)$ time per operation (queries and updates).

We label the two types of answers to soft nearest-neighbor (SNN) queries as hard or soft. A “standard” NN data structure is a special case of a SNN structure that always gives hard answers. However, in light of Theorem 2.2, a standard NN structure would not be as efficient as a SNN structure. For comparison, the best dynamic NN structure in $\mathbb{R}^{2}$ requires $O(\log^{5}n)$ time per operation [24, 52].

In our implementation, we use the following data structure. Given a point set $P$ and a point $q$, let $p^{*}_{i}$ denote the $i$-th closest point to $q$ in $P$.

Definition 2.3** (Dynamic $\varepsilon$-approximate $k$ nearest-neighbor ($k$-ANN) data structure).**

Maintain a dynamic set of points, $P$, subject to insertions, deletions, and $\varepsilon$-approximate $k$ nearest-neighbor queries: given a query point $q$ and an integer $k$ with $1\leq k\leq|P|$, return $k$ points $p_{1},\ldots,p_{k}\in P$ such that, for each $p_{i}$, $d(q,p_{i})\leq(1+\varepsilon)d(q,p^{*}_{i})$, where $\varepsilon>0$ is a constant known at construction time222Some approximate nearest-neighbor data structures [9] do not need to know $\varepsilon$ at construction time, and, in fact, allow $\varepsilon$ to be part of the query and to be different for each query. Clearly, such data structures are also valid for our needs..

We reduce each SNN query to a single $k$-ANN query with constant $\varepsilon$ and $k$. Once we show this reduction, Theorem 2.2 will follow from the following result by Arya et al. [9]:

Lemma 2.4** ([9]).**

In any fixed dimension, and for any $L_{p}$ metric, there is a dynamic $\varepsilon$-approximate $k$ nearest-neighbor data structure with $O(n\log n)$ preprocessing time and $O(\log n)$ time per operation (query and updates) for constant $k$ and $\varepsilon>0$.

2.1 Soft nearest-neighbor implementation

We maintain the point set $P$ in a dynamic $k$-ANN structure ($\varepsilon$ depends on the metric space, and will be determined later). In what follows, $q$ denotes an arbitrary query point and $p^{*}_{i}$ the $i$-th closest point to $q$ in $P$. For ease of presentation, we assume throughout the section that $d(q,p_{1}^{*})=1$. This scaling does not affect any result. Queries rely on the following lemma.

Lemma 2.5**.**

Consider a query $(q,k)$ to a $k$-ANN structure. If none of the $k$ returned points, $p_{1},\ldots,p_{k},$ is $p_{1}^{*}$, then, for each $p_{i}$ with $1\leq i\leq k$, we have that $d(q,p_{i})\leq(1+\varepsilon)^{i}$.

Proof 2.6**.**

*For $i=1$, the fact follows immediately from the definition of the $k$-ANN structure (and the assumption that $d(q,p_{1}^{*})=1$). For $i=2,\ldots,k$, note that $d(q,p_{i}^{*})\leq d(q,p_{i-1})$. This is because there are at least $i$ points within distance $d(q,p_{i-1})$ of $q$: $p_{1}^{*},p_{1},\ldots,p_{i-1}$. Thus, $d(q,p_{i})\leq(1+\varepsilon)d(q,p_{i}^{*})\leq(1+\varepsilon)d(q,p_{i-1}).$ The claim follows by induction.

It is illustrated in Figure 1.*

Let $S(q,r_{1},r_{2})$ denote a closed shell centered at $q$ with inner radius $r_{1}$ and outer radius $r_{2}$ (i.e., $S(q,r_{1},r_{2})$ is the difference between two balls centered at $q$, the bigger one of radius $r_{2}$ and the smaller one of radius $r_{1}$). From Lemma 2.5, we get the following.

Corollary 2.7**.**

Consider a query $(q,k)$ to an $k$-ANN structure. If none of the $k$ returned points, $p_{1},\ldots,p_{k}$, is $p_{1}^{*}$, then they all lie in $S(q,1,(1+\varepsilon)^{k})$.

We call a pair $(\varepsilon,k)$ valid parameters if, in any set of $k$ points inside a shell with inner radius $1$ and outer radius $(1+\varepsilon)^{k}$, there must exist two points $p,p^{\prime}$ satisfying $d(p,p^{\prime})<1$. Suppose that $(\varepsilon^{*},k^{*})$ are valid parameters. Initially, we construct the $k$-ANN structure using $1+\varepsilon^{*}$ as the approximation factor. Then we answer queries as in Algorithm 1.

Lemma 2.8**.**

If $(\varepsilon^{*},k^{*})$ are valid parameters, Algorithm 1 is correct.

Proof 2.9**.**

If a pair $p,p^{\prime}$ of points returned by the $k$-ANN structure satisfy $d(p,p^{\prime})<1$, $p$ and $p^{\prime}$ are a valid soft answer to the SNN query. Thus, consider the alternative case: no pair of the $k^{*}$ returned points is at distance ${<1}$. Then, because ($\varepsilon^{*},k^{*}$) are valid, at least one of the returned points must be outside of $S(q,1,(1+\varepsilon^{*})^{k^{*}})$. By the contrapositive of Corollary 2.7, one of them must be $p_{1}^{*}$.

As a side note, a SNN structure always returns a hard answer when queried from a point that is part of the closest pair of the set of points it maintains, as there is no closer pair. In this way, a SNN structure can be used to find the closest pair in $(\mathbb{R}^{\delta},L_{p})$, for constant $\delta$, in $O(n\log n)$ time by querying from every point. This matches the known runtimes in the literature [14].

2.2 Choice of parameters

We left open the question of finding valid parameters $(\varepsilon^{*},k^{*})$. This question is related to the kissing number of the metric space, which is the maximum number of points that can be on the surface of a unit sphere all at pairwise distance $\geq 1$. For instance, it is well known that the kissing number is $6$ in $(\mathbb{R}^{2},L_{2})$ and $12$ in $(\mathbb{R}^{3},L_{2})$. It follows that, in $(\mathbb{R}^{2},L_{2})$, $(\varepsilon^{*}=0,k^{*}=7)$ are valid parameters. Of course, we are interested in $\varepsilon^{*}>0$. Thus, our question is more general in the sense that our points are not constrained to lie on a sphere, but in a shell (and, to complicate things, the width of the shell depends on the number of points).

Lemma 2.10**.**

There are valid parameters in any metric space $(\mathbb{R}^{\delta},L_{p})$.

Proof 2.11**.**

Consider a shell with inner radius $1$ and outer radius $1+c$, for some constant $c>0$. A set of points in the shell at pairwise distance ${\geq 1}$ corresponds to a set of disjoint balls of radius $1/2$ centered inside the shell. Consider the volume of the intersection of the shell with such a ball. This volume is lower bounded by some constant, $v$, corresponding to the case where the ball is centered along the exterior boundary. Since the volume of the shell, $v_{s}$, is itself constant, the maximum number of disjoint balls of radius $1/2$ that fit in the shell is constant smaller than $v_{s}/v$. This is because no matter where the balls are placed, at least $v$ volume of the shell is inside any one of them, so, if there are more than $v_{s}/v$ balls, there must be some region in the shell inside at least two of them. This corresponds to two points at distance ${<1}$.

Set $k$ to be $v_{s}/v$, and $\varepsilon$ to be the constant such that $(1+\varepsilon)^{k}=1+c$. Then, $(\varepsilon,k)$ are valid parameters for $(\mathbb{R}^{\delta},L_{p})$.

The dependency of $k$-ANN structures on $1/\varepsilon$ is typically severe. Thus, for practical purposes, one would like to find a valid pair of parameters with $\varepsilon$ as big as possible. The dependency on $k$ is usually negligible in comparison, and, in any case, $k$ cannot be too large because the shell’s width grows exponentially in $k$. Thus, we narrow the question to optimizing $\varepsilon$: what is the largest $\varepsilon$ that is part of a pair of valid parameters?

We first address the case of $(\mathbb{R}^{2},L_{2})$, where we derive the optimal value for $\varepsilon$ analytically. We then give a heuristic, numerical algorithm for general $(\mathbb{R}^{\delta},L_{p})$ spaces.

Parameters in $(\mathbb{R}^{2},L_{2})$.

Let $\varepsilon_{\varphi}\approx 0.0492$ be the number such that $(1+\varepsilon_{\varphi})^{10}=\varphi$, where $\varphi=\frac{1+\sqrt{5}}{2}$ is the golden ratio. The valid parameters with largest $\varepsilon$ for $(\mathbb{R}^{2},L_{2})$ are $(\varepsilon^{*}<\varepsilon_{\varphi},k^{*}=10)$ ($\varepsilon^{*}$ can be arbitrarily close to $\varepsilon_{\varphi}$, but must be smaller). This follows from the following observations.

• •

The kissing number is $6$, so there are no valid parameters with $k<6$.

• •

The thinnest annulus (i.e., 2D shell) with inner radius $1$ such that $10$ points can be placed inside at pairwise distance ${\geq 1}$ has outer radius $\varphi=(1+\varepsilon_{\varphi})^{10}$. Figure 2, top, illustrates this fact. In other words, if the outer radius is any smaller than $\varphi$, two of the $10$ points would be at distance ${<1}$. Thus, any valid pair with $k=10$ requires $\varepsilon$ to be smaller than $\varepsilon_{\varphi}$, but any value smaller than $\varepsilon_{\varphi}$ forms a valid pair with $k=10$.

• •

For $6\leq k<10$ and for $k>10$, it is possible to place $k$ points at pairwise distance ${>1}$ in an annulus of inner radius $1$ and outer radius $(1+\varepsilon_{\varphi})^{k}$, and they are not packed “tightly”, in the sense that $k$ points at pairwise distance ${>1}$ can lie in a thinner annulus. This can be observed easily; Figure 2 (bottom) shows the cases for $k=9$ and $k=11$. Cases with $k<9$ can be checked one by one; in cases with $k>11$, the annulus grows at an increasingly faster rate, so placing $k$ points at pairwise distance ${>1}$ of each other becomes increasingly “easier”. Thus, for any $k\not=10$, any valid pair with that specific $k$ would require an $\varepsilon$ smaller than $\varepsilon_{\varphi}$.

Parameters in $(\mathbb{R}^{\delta},L_{p})$.

For other $L_{p}$ spaces, we suggest a numerical approach. We can do a binary search on the values of $\varepsilon$ to find one close to optimal. For a given value of $\varepsilon$, we want to know if there is any $k$ such that $(\varepsilon,k)$ are valid. We can search for such a $k$ iteratively, trying $k=1,2,\ldots$ (the answer will certainly be “no” for any $k$ smaller than the kissing number). Note that, for a fixed $k$, the shell has constant volume. As in Lemma 2.10, let $v$ be the volume of the intersection between the shell and a ball of radius $1/2$ centered on the exterior boundary of the shell. As argued before, if $kv$ is bigger than the shell’s volume, then $(\varepsilon,k)$ are valid parameters. For the termination condition, note that if in the iterative search for $k$, $k$ reaches a value where the volume of the shell grows more than $v$ in a single iteration, no valid $k$ for that $\varepsilon$ will be found, as the shell grows faster than the new points cover it.

Besides the volume check, one should also consider a lower bound on how much of the shell’s surface (both inner and outer) is contained inside an arbitrary ball. We can then see if, for a given $k$, the amount of surface contained inside the $k$ balls is bigger than the total surface of the shell, at which point two balls surely intersect. This check finds better valid parameters than the volume one for relatively thin shells, where the balls “poke” out of the shell on both sides.

3 Multi-Fragment Euclidean TSP

The Euclidean Travelling Salesperson Problem asks to find, given a set of points, a closed tour (a closed polygonal chain) through all the points of shortest length. The problem is NP-hard even in this geometric setting, but a polynomial-time approximation scheme is known [8].

In this section, we consider a classic greedy heuristic for constructing TSP tours, multi-fragment TSP. In this algorithm, each point starts as a single-node path. While there is more than one path, connect the two closest paths. Here, the distance $d(a,b)$ between two paths $a,b$ is measured as the minimum distance between their endpoints, and connecting two paths means adding the edge between their closest endpoints. Once there is a single path left, connect their endpoints. We call the tour resulting from this process the multi-fragment tour.

The multi-fragment algorithm was proposed by Bentley [12] specifically in the geometric setting. Its approximation ratio is $O(\log n)$ [64, 19]. Nonetheless, it is used in practice due to its simplicity and empirical support that it generally performs better than other heuristics [34, 51, 57, 59, 13].

We are interested in the complexity of computing the multi-fragment tour. A straightforward implementation of the multi-fragment algorithm is similar to Kruskal’s minimum spanning tree algorithm: sort the $\binom{n}{2}$ pairs of points by increasing distances and process them in order: for each pair, if the two points are endpoints of separate paths, connect them. The runtime of this algorithm is $O(n^{2}\log n)$. Eppstein [36] uses dynamic closest pair data structures to compute the multi-fragment tour in $O(n^{2})$ time (for arbitrary distance matrices). Bentley [12] gives a $K$-$d$ tree-based implementation and says that it appears to run in $O(n\log n)$ time on uniformly distributed points in the plane. We give a NNC-type algorithm that compute the multi-fragment tour in $O(n\log n)$ in any fixed dimensions. We do not know of any prior worst-case subquadratic algorithm.

3.1 Global-local equivalence in multi-fragment TSP

Since the multi-fragment algorithm operates on paths rather than points, it will be convenient to think of the input as a set of paths (a path is an open polygonal chain, although, in the context of the algorithm, only the coordinates of the endpoints are relevant). The input to Euclidean TSP corresponds to a set of paths where all paths are single-point paths. Consider the following two strategies for constructing a tour from a set of paths, where we use $a\cup b$ to denote the path resulting from connecting paths $a$ and $b$:

• •

While there is more than one path, connect two paths using one of the following strategies:

Connect the closest pair of paths. 2. 2.

Connect two mutually nearest-neighbor paths.

• •

Connect the two endpoints of the final path.

Strategy 1 corresponds to the multi-fragment algorithm. Note that Strategy 2 is non-deterministic, and that Strategy 1 is a special case of Strategy 2. In this section, we show that any execution of Strategy 2 computes the multi-fragment tour. Note the similarity between multi-fragment TSP and hierarchical clustering. We can see that in multi-fragment TSP we have a notion equivalent to reducibility in agglomerative hierarchical clustering (Equation 1).

Lemma 3.1** (Reducibility in multi-fragment TSP).**

Let $a,b,$ and $c$ be paths. Then, $d(a\cup b,c)\geq\min{(d(a,c),d(b,c))}$.

Proof 3.2**.**

The distance between paths is defined as the minimum distance between their endpoints, and the two endpoints of $a\cup b$ are a subset of the four endpoints of $a$ and $b$.

Lemma 3.3** (Global-local equivalence in multi-fragment TSP).**

Assuming that there are no ties in the pairwise distances between paths, strategies 1 and 2 produce the same tour.

We adapt the proof of global-local equivalence for agglomerative hierarchical clustering presented in [60]. We note that we can break ties using a consistent rule, such as breaking ties by the smallest index in the input.

Proof 3.4**.**

Let $P$ be a set of paths, and let $S_{1}(P)$ denote the sequence of path pairs connected by Strategy 1 starting from $P$, and $T_{1}(P)$ the corresponding resulting tour. Similarly, let $S_{2}(P)$ denote one of the possible sequences of path pairs connected by an instantiation of Strategy 2 starting from $P$, and $T_{2}(P)$ the corresponding tour. We need to show that $T_{1}(P)=T_{2}(P)$.

Proceed by induction on $|P|$. If $|P|=1$, both tours are the same because no connections happen. Thus, let $|P|>1$. Let $(a,b)$ be the first pair of paths in $S_{2}(P)$. Then, consider the set $P^{\prime}=P\setminus\{a,b\}\cup\{a\cup b\}$. The tour $T_{2}(P)$ can be seen as the tour obtained by starting from the set $P^{\prime}$ and connecting the same paths as in $S_{2}(P)$ after the first connection $(a,b)$. Note that $|P|=|P|-1$. Thus, by the inductive hypothesis, $T_{2}(P)=T_{1}(P^{\prime})$. The bulk of the proof is to show that $T_{1}(P^{\prime})=T_{1}(P)$.

First, note that $(a,b)$ is in $S_{1}(P)$: initially, $a$ and $b$ are MNN paths (since they are the first pair chosen by Strategy 2). Then, they remain so throughout the algorithm until they are connected. This is because (i) MNN paths are not connected with other paths, and (ii) by reducibility (Lemma 3.1), MNN paths stay so even if other paths are connected (i.e., if $x$ and $y$ are connected, $x\cup y$ is not closer to $a$ (or $b$) than the closest of $x$ and $y$).

Let $(a,b)$ be the $k$-th pair in $S_{1}(P)$. Next, we show that the first $k-1$ pairs in $S_{1}(P)$ and $S_{1}(P^{\prime})$ are the same and in the same order. Let $(x,y)$ be the first pair of paths in $S_{1}(P)$. By Strategy 1, $d(x,y)$ is minimum among all distances between paths in $P$. By Lemma 3.1, $a\cup b$ is not closer to $x$ or $y$ than $a$ or $b$. Thus, in $S_{1}(P^{\prime})$, $d(x,y)$ is also minimum, so $(x,y)$ is also the first element in $S_{1}(P^{\prime})$. The claim for the next $k-2$ pairs follows analogously by induction.

Finally, note that after the first $k$ connections in $S_{1}(P)$ and the first $k-1$ connections in $S_{1}(P^{\prime})$, the corresponding partial solutions are the same. After that point, all the connections, and the final solution, must be the same in both, so $T_{1}(P)=T_{1}(P^{\prime})$.

We note that Lemma 3.3 holds for arbitrary distance matrices.

3.2 Soft nearest-neighbor chain for multi-fragment Euclidean TSP

Given that we have global-local equivalence (Lemma 3.3), we can use the NNC algorithm to compute the multi-fragment tour using Strategy 2. A straightforward adaptation of the NNC algorithm, paired with the NN structure from [24, 52], yields a $O(n\log^{5}n)$ runtime for $\mathbb{R}^{2}$. However, we skip this result and jump directly to our main result:

Theorem 3.5**.**

The multi-fragment tour of a set of $n$ points in any fixed dimension, and under any $L_{p}$ metric, can be computed in $O(n\log n)$ time.

We use a variation of the NNC algorithm that uses a SNN structure instead of the usual NN structure, which we call soft nearest-neighbor chain (SNNC). For this, we need a SNN structure for paths instead of points. That is, a structure that maintains a set of (possibly single-node) paths, and, given a query path $q$, returns the closest path to $q$ or two paths which are closer.

A soft nearest-neighbor structure for paths.

We simulate a SNN structure for paths with a SNN structure for points. Given a set of paths, we maintain the set of path endpoints in the SNN structure for points. Updates are straightforward: we add or remove both endpoints of the path. Given a query path $q$ with endpoints $\{q_{1},q_{2}\}$, we do a SNN query from each endpoint of the path. If both answers are hard (assuming that the path has two distinct endpoints, otherwise, just the one), then we find the true NN of the path, and we can return it. However, there is a complication with soft answers: the two points returned could be the endpoints of the same path. Thus, it could be the case that we find two closer points, but not two closer paths, as we need. The solution is to modify the specification of the SNN structure for points so that soft answers, instead of returning two points closer to each other than the query point to its NN, return three pairwise closer points. We call this a three-way SNN structure. In the context of using the structure for paths, this guarantees that even if two of the three endpoints belong to the same path, at least two different paths are involved.

Lemma 3.6 shows how to obtain a three-way SNN structure for points, Algorithm 2 shows the full algorithm for answering SNN queries about paths using a three-way SNN structure for points, and Lemma 3.8 shows its correctness.

Lemma 3.6**.**

In any fixed dimension and for any $L_{p}$ metric, there is a three-way SNN structure with $O(n\log n)$ preprocessing time and $O(\log n)$ operation time (queries and updates).

Proof 3.7**.**

Recall the implementation of the SNN structure from Section 2. To obtain a three-way SNN structure, we need to change the values of $\varepsilon$ and $k$ to make the shell smaller and $k$ bigger, so that if there are $k$ points in a shell of inner radius $1$ and outer radius $(1+\varepsilon)^{k}$, then there must be at least three points at pairwise distance less than $1$. The method described in Section 2.2 for finding valid parameters in $(R^{\delta},L_{p})$ also works here. It only needs to be modified so that the area (or surface) of the shell is accounted for twice. Since $k$ and $\varepsilon$ are still constant, this does not affect the asymptotic runtimes in Theorem 2.2.

Lemma 3.8**.**

In any fixed dimension, and for any $L_{p}$ metric, we can maintain a set of $n$ paths in a SNN structure for paths with $O(n\log n)$ preprocessing time and $O(\log n)$ operation time (queries and updates).

Proof 3.9**.**

All the runtimes follow from Lemma 3.6: we maintain the set of path endpoints in a three-way SNN structure $S$. The structure $S$ can be initialized in $O(n\log n)$ time. Updates require two insertions or deletions to $S$, taking $O(\log n)$ time each. Algorithm 2 for queries clearly runs in $O(\log n)$ time. We argue that it returns a valid answer. Let $q$ be a query path with endpoints $\{q_{1},q_{2}\}$, and consider the three possible cases:

• •

Both answers are hard. In this case, we find the closest path to each endpoint, and, by definition, the closest of the two is the NN of $q$.

• •

One answer is soft and the other is hard. Let $p$ be the hard answer to $q_{1}$ and $(a,b,c)$ the soft answer to $q_{2}$ (wlog). Let $P$ and $P^{\prime}$ be the two closest paths among the paths with endpoints $a,b,$ and $c$. If $d(p,q)<d(P,P^{\prime})$, then, the path with $p$ as endpoint must be the NN of $q$, because there is no endpoint closer than $d(P,P^{\prime})$ to $q_{2}$. Otherwise, $P,P^{\prime}$ is a valid soft answer, as they are closer to each other than either endpoint of $q$ to their closest endpoints.

• •

Both answers are soft. Assume (wlog) that the NN of $q$ is closer to $q_{1}$ than $q_{2}$. Then, the soft answer to $q_{1}$ gives us two paths closer to each other than $q$ to its NN, so we return a valid soft answer.

The soft nearest-neighbor chain algorithm.

We use a SNN for paths. In the context of this algorithm, let us think of a SNN answer, hard or soft, as being a set of two paths. If the answer is hard, then one of the paths returned in the answer is the query path itself, and the remaining path is its NN. Now, we can establish a comparison relationship between SNN answers (independently of their type): given two SNN answers $\{a,b\}$ and $\{c,d\}$, we say that $\{a,b\}$ is better than $\{c,d\}$ if and only if $d(a,b)<d(c,d)$.

The input is a set of paths, where we again assume unique distances. The algorithm maintains a stack (the chain) of nodes, where each node consists of a pair of paths (with the exception of the first node in the chain, which contains a single path). In particular, each node in the chain is the best SNN answer among two queries for the two paths in the predecessor node (when querying from a path, we remove it from the structure temporarily, so that the answer is not itself).

The algorithm starts with an arbitrary path in the chain. If the chain ever becomes empty and there is still more than one path, the chain is restarted at an arbitrary path. If the best answer we get from the SNN structure is precisely the node currently at the top of the chain, we connect both paths contained in it and remove the node and its predecessor from the chain. Otherwise, we append the answer to the top of the chain as a new node. See Algorithm 3 for a full description of the algorithm and Figure 3 for a snapshot of the algorithm.

Lemma 3.10**.**

The following invariants hold at the beginning of each iteration of Algorithm 3:

The SNN structure $S$ contains a set of disjoint paths. 2. 2.

If node $R$ appears after node $S$ in the chain, then $R$ is better than $S$. 3. 3.

Every path in $S$ appears in at most two nodes in the chain, in which case they consist of occurrences in two consecutive nodes. 4. 4.

The chain only contains paths in $S$.

Proof 3.11**.**

The claim holds initially. Each time two paths $a,b$ are connected, one endpoint of each becomes an internal point in the new path $a\cup b$. Since $a$ and $b$ are removed from $S$, no path can be connected to those endpoints. 2. 2.

We show it for the specific case where $R$ is immediately after $S=\{s,t\}$ in the chain, which suffices. Note that $R$ is different than $S$, or it would not have been added to the chain. We distinguish between two cases:

• •

$s$* and $t$ were MNNs when $R$ was added. Then, $R$ had to be a soft answer from $s$ or $t$, which would have to be better than $\{s,t\}$.*

• •

$s$* and $t$ were not MNNs when $R$ was added. Then, $s$ had a closer path than $t$ (wlog). Thus, whether the answer for $s$ was soft or hard, the answer had to be better than $\{s,t\}$.* 3. 3.

Assume, for a contradiction, that three nodes $X=\{p,x\}$, $Y$ and $Z=\{p,z\}$ appear in the chain in order $X$, $Y$, $Z$ (not necessarily consecutively, and with $p$ possibly in $Y$). By Invariant 2, $Z$ is better than $Y$. It is easy to see that if $z_{1}$ and $z_{2}$ are the two endpoints of path $z$, then $z_{1}$ and $z_{2}$ were endpoints of paths since the beginning of the algorithm. Thus, the answer for $p$ when $X$ was at the top of the chain had to be a pair at distance at most $\min(d(p,z_{1}),d(p,z_{2}))$. Note that $\min(d(p,z_{1}),d(p,z_{2}))=d(p,z)$, contradicting that $Z$ is better than $Y$. 4. 4.

We show that no node in the chain contains paths that have already been connected to form bigger paths. Whenever the two paths in the node at the top of the chain are connected, we remove the node from the chain. By Invariant 3, their only other occurrence can only be in the predecessor node, which is also removed. In addition, since the paths are removed from $S$ when they are connected, newly added nodes to the chain only contain paths that have not been connected to form bigger path yet.

Lemma 3.12**.**

Paths connected in Algorithm 3 are MNNs in the set of paths in the SNN structure.

Proof 3.13**.**

Let $\{u,v\}$ be the node at the top of the chain, and $A$ the best SNN answer among the $u$ and $v$ queries. If $u$ and $v$ are not MNNs, at least one of them, $u$ (wlog), has a closer path than the other, $v$, so the answer for $u$ cannot be $\{u,v\}$. By the contrapositive, if the best answer $A$ is $\{u,v\}$, then $u$ and $v$ are MNNs. In the algorithm, $u$ and $v$ are connected precisely when $A=\{u,v\}$.

Proof 3.14** (Proof of Theorem 3.5).**

We show that Algorithm 3 computes the multi-fragment tour in $O(n\log n)$ time. In particular, it implements Strategy 2: the SNN structure maintains a set of paths, and the algorithm repeatedly connects MNNs (Lemma 3.12) By global-local equivalence (Lemma 3.3), this produces the multi-fragment tour.

Note that the chain is acyclic in the sense that each node contains a path from the current set of paths in $S$ (Invariant 4) not found in previous nodes (by Invariant 3). Thus, the chain cannot grow indefinitely, so, eventually, paths get connected. The main loop does not halt until there is a single path.

If there are $n$ paths at the beginning, there are $2n-1$ different paths throughout the algorithm. This is because each connection removes two paths and adds one new path. At each iteration, either two paths are connected, which happens $n-1$ times, or one node is added to the chain. Since there are $n-1$ connections, each of which triggers the removal of two nodes in the chain, the total number of nodes removed from the chain is $2n-2$. Since every node added is removed, the number of nodes added to the chain is also $2n-2$. Thus, the total number of iterations is $3n-3$. Therefore, the total running time is $O(P(n)+nT(n))$, where $P(n)$ and $T(n)$ are the preprocessing and operation time of a SNN structure for paths. By Lemma 3.8, this can be done in $O(n\log n)$ time.

Incidentally, the maximum-weight matching problem has the same type of global-local equivalence: consider the classic factor-$2$ approximation greedy algorithm that picks the heaviest edge at each iteration and discards the neighboring edges (edges with a shared endpoint) [10]. An alternative algorithm that picks any edge which is heavier than its neighbors produces the same matching [48], which we call the greedy matching. In the geometric setting, we are interested in matching points minimizing distances instead (but the mentioned results still hold). The greedy algorithm repeatedly matches the closest pair. It is possible to modify Algorithm 3 to find the greedy matching in $O(n\log n)$ time in any fixed dimension. The algorithm is, in fact, simpler, because the SNN structure only needs to maintain points instead of paths, and matched points are removed permanently (unlike connected paths which are re-added to the set of paths). However, this is not a new result, as there is a dynamic closest pair data structure with $O(\log n)$ time per operation [14] which can be used to find the greedy matching in the same time bound.

3.3 Steiner TSP

In the traditional, non-Euclidean setting, a TSP instance consists of a complete graph with arbitrary distances. We remark that global-local equivalence (Lemma 3.3) still holds in this general setting. In this context, the nearest neighbor of a path can be found in $O(n)$ time by iterating through the adjacency lists of both endpoints, where $n$ is the number of nodes. Using this linear search, we can easily compute the multi-fragment tour in $O(n^{2})$ time with a NNC-based algorithm. It is a simpler version of Algorithm 3 that only has to handle hard answers and does not need any sophisticated data structures. This improves upon the natural way to implement the multi-fragment heuristic, which is to sort the $\Theta(n^{2})$ edges by weight. Sorting requires $\Theta(n^{2}\log n)$ time.

This is the first use of NNC in a graph-theoretical setting, but the fact of the matter is that the NNC algorithm can be used in any setting where we can find nearest neighbors efficiently. Consider the related Steiner TSP problem [29]: given a weighted, undirected graph and a set of $k$ node sites $P\subseteq V$, find a minimum-weight tour (repeated vertices and edges allowed) in $G$ that goes at least once through every site in $P$. Nodes not in $P$ do not need to be visited. For instance, $G$ could represent a road network, and the sites could represent the daily drop-off locations of a delivery truck. See [30, 68] for more applications.

Recently, Eppstein et al. [39] gave a NN structure for graphs from graph families with sublinear separators, which is the same as the class of graphs with polynomial expansion [33]. For instance, planar graphs have $O(\sqrt{n})$-size separators333Other important families of sparse graphs with sublinear separators include $k$-planar graphs [32], bounded-genus graphs [43], minor-closed graph families [53], and graphs that model road networks (better than, e.g., $k$-planar graphs) [40].. This data structure maintains a subset of nodes $P$ of a graph $G$, and, given a query node $q$ in $G$, returns the node in $P$ closest to $q$. It allows insertions and deletions to and from the set $P$. We cite their result in Lemma 3.15.

Lemma 3.15** ([39]).**

Given an $n$-node weighted graph from a graph family with separators of size $S(n)=n^{c}$, with $0<c<1$, which can be constructed in $O(n)$ time, there is a dynamic444They ([39]) use the term reactive for the data structure instead of dynamic, to distinguish from other types of updates, e.g., edge insertions and deletions. nearest-neighbor data structure requiring $O(nS(n))$ space and preprocessing time and that answers queries in $O(S(n))$ time and updates in $O(S(n)\log n)$ time.

As mentioned, one way to implement the multi-fragment heuristic is to sort the $\binom{k}{2}$ pairs of sites by increasing distances, and process them in order: for each pair, if the two sites are endpoints of separate paths, connect them. The bottleneck is computing the distances. Running Dijkstra’s algorithm from each site in a sparse graph, this takes $O(k(n\log n))$ (or $O(kn)$ in planar graphs [46]). When $k$ is $\Theta(n)$, this becomes $O(n^{2}\log n)$. We do not know of any prior faster algorithm to compute the multi-fragment tour for Steiner TSP.

Since we have global-local equivalence, we can use the NNC algorithm to construct the multi-fragment tour in $O(P(n)+kT(n))$ time, where $P(n)$ and $T(n)$ are the preprocessing and operation time of a nearest-neighbor structure. Thus, using the structure from [39], we get:

Theorem 3.16**.**

The multi-fragment tour for the steiner TSP problem can be computed in $O(nS(n)+kS(n)\log n)$-time in weighted graphs from a graph family with separators of size $S(n)=n^{c}$, with $0<c<1$.

Finally, in graphs of bounded treewidth, which have separators of size $O(1)$, the data structure from [39] achieves $P(n)=O(n\log n)$ and $T(n)=O(\log^{2}n)$, so we can construct a multi-fragment tour in $O(n\log n+k\log^{2}n)$.

4 Motorcycle Graphs

An important concept in geometric computing is the straight skeleton [5]. It is a tree-like structure similar to the medial axis of a polygon, but which consists of straight segments only. Given a polygon, consider a shrinking process where each edge moves inward, at the same speed, in a direction perpendicular to itself. The straight skeleton of the polygon is the trace of the vertices through this process. Some of its applications include computing offset polygons [35], medical imaging [28], polyhedral surface reconstruction [63, 11], and computational origami [31]. It is a standard tool in geometric computing software [23].

The current fastest algorithms for computing straight skeletons consist of two main steps [27, 49, 50]. The first step is to construct a motorcycle graph induced by the reflex vertices of the polygon. The second step is a lower envelope computation. With current algorithms, the first step is more expensive, but it only depends on the number of reflex vertices, $r$, which might be smaller than the total number of vertices, $n$. Thus, no step dominates the other in every instance. In this section, we discuss the first step. The second step can be done in $O(n\log n)$ time for simple polygons [18], in $O(n\log n\log r)$ time for arbitrary polygons [26], and in $O(n\log n\log m)$ time for planar straight line graphs with $m$ connected components [18].

The motorcycle graph problem can be described as follows (see Figure 4, top) [35]. The input consists of $n$ points in the plane, with associated directions and speeds (the motorcycles). Consider the process where all the motorcycles start moving at the same time, in their respective directions and speeds. Motorcycles leave a trace behind that acts as a “wall” such that other motorcycles crash and stop if they reach it. Some motorcycles escape to infinity while others crash against the traces of other motorcycles. The motorcycle graph is the set of traces.

Most existing algorithms rely on three-dimensional ray-shooting queries. Indeed, if time is seen as the third dimension, the position of a motorcycle starting to move from $(x,y)$, at speed $s$, in the direction $(u,v)$, forms a ray (if it escapes) or a segment (if it crashes) in three dimensions, starting at $(x,y,0)$ in the direction $(u,v,1/s)$. In particular, the impassable traces left behind by the motorcycles correspond to infinite vertical “curtains” – wedges or trapezoidal slabs, depending on whether they are bounded below by a ray or a segment.

Thus, ray-shooting queries help determine which trace a motorcycle would reach first, if any. Of course, the complication is that as motorcycles crash, their potential traces change. Early algorithms handle this issue by computing the crashes in chronological order [35, 27]. The best previously known algorithm, by Vigneron and Yan [67], is the first that computes the crashes in non-chronological order. Our NNC-based algorithm improves upon it by reducing the number of ray-shooting queries needed from $O(n\log n)$ to $3n$, and simplify significantly the required data structures. It is also non-chronological, but follows a completely new approach.

4.1 Algorithm description

In the algorithm, we distinguish between undetermined motorcycles, for which the final location is still unknown, and determined motorcycle, for which the final location is already known. We use a dynamic three-dimensional ray-shooting data structure. In the data structure, determined motorcycles have wedges or slabs as curtains, depending on whether they escape or not. Undetermined motorcycles have wedge curtains, as if they were to escape. Thus, curtains of undetermined motorcycles may reach points that the corresponding motorcycles never get to.

For an undetermined motorcycle $m$, we define its nearest neighbor to be the motorcycle, determined or not, against which $m$ would crash next according to the set of curtains in the data structure. Motorcycles that escape may have no NN. Finding the NN of a motorcycle $m$ corresponds to one ray-shooting query. Note that $m$ may actually not crash against the trace of its NN, $m^{\prime}$, if $m^{\prime}$ is undetermined and happens to crash early. On the other hand, if $m^{\prime}$ is determined, then $m$ definitely crashes into its trace.

We begin with all motorcycles as undetermined. Our main structure is a chain (a stack) of undetermined motorcycles such that each motorcycle is the NN of the previous one. In contrast to typical applications of the NNC algorithm, here “proximity” is not symmetric: there may be no “mutually nearest-neighbors”. In fact, the only case where two motorcycles are MNNs is the degenerate case where two motorcycles reach the same point simultaneously. That said, mutually nearest neighbors have an appropriate analogous in the asymmetric setting: nearest-neighbor cycles, $m_{1}\rightarrow m_{2}\rightarrow\cdots\rightarrow m_{k}\rightarrow m_{1}$. Our algorithm relies on the following key observation: if we find a nearest-neighbor cycle of undetermined motorcycles, then each motorcycle in the cycle crashes into the next motorcycle’s trace. This is easy to see from the definition of nearest neighbors, as it means that no motorcycle outside the cycle would “interrupt” the cycle by making one of them crash early. Thus, if we find such a cycle, we can determine all the motorcycles in the cycle at once (this can be seen as a type of chronological global-local equivalence).

Starting from an undetermined motorcycle, following a chain of nearest neighbors inevitably leads to (a) a motorcycle that escapes, (b) a motorcycle that is determined, or (c) a nearest-neighbor cycle. In all three cases, this allows us to determine the motorcycle at the top of the chain, or, in Case (c), all the motorcycles in the cycle. See Figure 4, bottom. Further, note that we only modify the curtain of the newly determined motorcycle(s). Thus, if we determine the motorcycle $m$ at the top of the chain, only the NN of the second-to-last motorcycle in the chain may have changed, and similarly in the case of the cycle. Consequently, the rest of the chain remains consistent. Algorithm 4 shows the full pseudo code.

4.2 Analysis

Clearly, every motorcycle eventually becomes determined, and we have already argued in the algorithm description that irrespective of whether it becomes determined through Case (a), (b), or (c), its final position is correct. Thus, we move on to the complexity analysis. Each “clipping” update can be seen as an update to the ray-shooting data structure: we remove the wedge and add the slab.

Theorem 4.1**.**

Algorithm 4 computes the motorcycle graph in time $O(P(n)+nT(n))$, where $P(n)$ and $T(n)$ are the preprocessing time and operation time (maximum between query and update) of a dynamic, three-dimensional ray-shooting data structure.

Proof 4.2**.**

Each iteration of the algorithm makes one ray-shooting query. At each iteration, either a motorcycle is added to the chain (Case (d)), or at least one motorcycle is determined (Cases (a—c)).

Motorcycles begin as undetermined and, once they become determined, they remain so. This bounds the number of Cases (a—c) to $n$. In Cases (b) and (c), one undetermined motorcycle may be removed from the chain. Thus, the number of undetermined motorcycles removed from the chain is at most $n$. It follows that Case (d) happens at most $2n$ times.

Overall, the algorithm takes at most $3n$ iterations, so it needs no more than $3n$ ray-shooting queries and at most $n$ “clipping” updates where we change a triangular curtain into a slab. It follows that the runtime is $O(P(n)+nT(n))$.

In terms of space, we only need a linear amount besides the space required by the data structure.

The previous best known algorithm runs in time $O(P(n)+n(T(n)+\log n)\log n)$ [67]. Besides ray-shooting queries, it also uses range searching data structures, which do not increase the asymptotic runtime but make the algorithm more complex.

Agarwal and Matoušek [2] give a ray-shooting data structure for curtains in $\mathbb{R}^{3}$ which achieves $P(n)=O(n^{4/3+\varepsilon})$ and $T(n)=O(n^{1/3+\varepsilon})$ for any $\varepsilon>0$. Using this structure, both our algorithm and the algorithm of Vigneron and Yan [67] run in $O(n^{4/3+\varepsilon})$ time for any $\varepsilon>0$. If both algorithms use the same $\varepsilon$ in the ray-shooting data structure, then our algorithm is asymptotically faster by a logarithmic factor.

4.3 Special cases and remarks

Consider the case where all motorcycles start from the boundary of a simple polygon with $O(n)$ vertices, move through the inside of the polygon, and also crash against the edges of the polygon. In this setting, the motorcycle trajectories form a connected planar subdivision. There are dynamic ray-shooting queries for connected planar subdivisions that achieve $T(n)=O(\log^{2}n)$ [45]. Vigneron and Yan used this data structure in their algorithm to get a $O(n\log^{3}n)$-time algorithm for this case [67]. Our algorithm brings this down to $O(n\log^{2}n)$. Furthermore, their other data structures require that coordinates have $O(\log n)$ bits, while we do not have this requirement.

Vigneron and Yan also consider the case where motorcycles can only go in $C$ different directions. They show how to reduce $T(n)$ to $\min(O(C\log^{2}n,C^{2}\log n)$, leading to a $O(n\log^{2}nC\min(\log n,C))$ algorithm for motorcycle graphs in this setting. Using the same data structures, the NNC algorithm improves the runtime to $O(n\log nC\min(\log n,C))$.

A remark on the use of our algorithm for computing straight skeletons: degenerate polygons where two shrinking reflex vertices collide gives rise to a motorcycle graph problem where two motorcycles collide head on. To compute the straight skeleton, a new motorcycle should emerge from the collision. Our algorithm does not work if new motorcycles are added dynamically (such a motorcycle could, e.g., disrupt a NN cycle already determined), so it cannot be used in the computation of straight skeletons of degenerate polygons.

As a side note, the NNC algorithm for motorcycle graphs is reminiscent of Gale’s top trading cycle algorithm [66] from the field of economics. That algorithm also works by finding “first-choice” cycles. We are not aware of whether they use a NNC-type algorithm to find such cycles; if they do not, they certainly can; if they do, then at least our use is new in the context of motorcycle graphs.

5 Stable Matching Problems

We introduce the narcissistic k-attribute stable matching problem, a special case of k-attribute stable matching, and show that it belongs to the class of symmetric stable matching problems. We use this fact to give an efficient NNC-type algorithm for it.

Stable matching.

The stable matching problem studies how to match two sets of agents in a market where each agent has its own preferences about the agents of the other set in a “stable” manner. Some of its applications include matching hospitals and residents [62] and on-line advertisement auctions [4]. It was originally formulated by Gale and Shapley [41] in the context of establishing marriages between $n$ men and $n$ women, where each man ranks the women and the women rank the men. A matching between the men and women is stable if there is no blocking pair: a man and woman who prefer each other over their assigned choices.

Gale and Shapley [41] showed that a stable solution exists for any set of preferences (and it might not be unique), and presented the deferred-acceptance algorithm, which finds a stable matching in $O(n^{2})$ time.

5.1 Restricted models

For arbitrary preference lists, Gale–Shapley’s deferred-acceptance algorithm is worst-case optimal, as storing all the preferences already requires $\Theta(n^{2})$ space (quadratic lower bounds are known also for “simpler” questions, like verifying stability of a given matching [44]). This inspired work on finding subquadratic algorithms in restricted settings where preferences can be specified in subquadratic space. Such models are collectively called succinct stable matching [58]. We introduce a new model which is a special case of the following three models (of which none is a special case of another):

$k$-attribute model [15].

Each agent $p$ has a vector $\vec{p_{a}}$ of $k$ numerical attributes, and a vector $\vec{p_{w}}$ of $k$ weights according to how much $p$ values each attribute in a match. Then, each agent $p$ ranks the agents in the other set according to the objective function $f_{p}(q)=\vec{p_{w}}\cdot\vec{q_{a}}$, the linear combination of the attributes of $q$ according to the weights of $p$.

Narcissistic stable matching.

This term is used to describe models where the preferences of each agent reflect their own qualities in some way (e.g., in [25, 58]).

Symmetric stable matching [37].

Consider the setting where each agent $p$ has an arbitrary objective function, $f_{p}(q)$, and ranks the agents according to $f_{p}(q)$ (note that any set of preference lists can be modeled in this way). The preferences are called symmetric if for any two agents $p,q$ in different sets, $f_{p}(q)=f_{q}(p)$.

In this paper, we consider the natural narcissistic interpretation of the $k$-attribute model, where $\vec{p_{a}}=\vec{p_{w}}$. That is, each agent weighs each attribute according to its own value in that attribute. To illustrate this model, consider a centralized dating service where $k$ attributes are known for each person, such as income, intelligence, and so on. In the general $k$-attribute model, each person assigns weights to the attributes according to their preferences. The narcissistic assumption that $\vec{p_{a}}=\vec{p_{w}}$ implies that someone with, say, a high income, values income more than someone with a relatively smaller income.

We make a general position assumption that there are no ties in the preference list of each agent. In addition, in this model each agent is uniquely determined by its attribute vector, so we do not distinguish between the agents themselves and their $k$-dimensional vectors. We obtain the following formal problem.

Definition 5.1** (Narcissistic $k$-attribute stable matching problem).**

Find a stable matching between two sets of $n$ vectors in $\mathbb{R}^{k}$, where a vector $\vec{p}$ prefers $\vec{q}$ over $\vec{q^{\prime}}$ if and only if $\vec{p}\cdot\vec{q}>\vec{p}\cdot\vec{q^{\prime}}$.

We give a $O(n^{2-4/(k+2+\varepsilon)})$-time algorithm for the problem. Without the narcissistic assumption, the $k$-attribute model becomes less tractable: Künnemann et al. [58] showed that no strongly subquadratic-time algorithm exists if $k=\omega(\log n)$ assuming the Strong Exponential Time Hypothesis, even if the weights and attributes take Boolean values. (Similarly to us, [58] also studied some restricted cases and presented a $O(C^{2k}n(k+\log n))$-time algorithm for the case where attributes and weights may have only $C$ different values and a $\tilde{O}(n^{2-1/\lfloor k/2\rfloor})$-time algorithm for the asymmetric case where one of the sets has a single attribute and the other has $k$.555The $\tilde{O}$ notation ignores logarithmic factors.)

It is easy to see that our setting is symmetric: since $\vec{p_{a}}=\vec{p_{w}}$, for any two agents $p,q$ we have $f_{p}(q)=\vec{p_{w}}\cdot\vec{q_{a}}=\vec{p_{a}}\cdot\vec{q_{w}}=f_{q}(p)$. Eppstein et al. [37] showed that in symmetric models, the NNC algorithm can be used. Specifically, they introduced symmetric stable matching as an abstraction of the case where the agents are points in a metric space and they rank the agents in the other set by proximity [7]. They showed that if preferences are symmetric, the problem has special properties: there is a unique stable matching and it can be found by repeatedly matching the two unmatched elements with the highest objective function value. In addition, they showed the global-local equivalence: it suffices to match any two elements who have each other as first choice, called soul mates, which are mutually nearest neighbors if the preferences are distance-based.

The algorithm of Eppstein et al. [38, 37] for symmetric stable matching was the first use of NNC outside of hierarchical clustering. The algorithm is a bichromatic version of the NNC algorithm, where each individual in the chain is followed by its first choice among the unmatched individuals in the other set. Following such a chain inevitably leads to soul mates, which are then matched and removed permanently (here, the symmetry assumption is the key to avoid cycles in the chain). The execution relies on a dynamic first-choice data structure, which maintains the elements in one set, and, given a query element $q$ from the other set, returns the first choice of $q$ among the elements in the structure. The final result is as follows:

Lemma 5.2** ([37]).**

Given a first-choice data structure with $P(n)$ preprocessing time and $T(n)$ operation time (maximum between query and update), a symmetric stable matching problem can be solved in $O(P(n)+nT(n))$.

5.2 Our algorithm for narcissistic $k$-attribute stable matching

Adapting the NNC algorithm from [37] to our model simply requires using an appropriate structure for first-choice queries. In our case, the first-choice data structure should maintain a set of vectors, and, given a query vector, return the vector maximizing the dot product with the query vector. In the dual, this becomes ray shooting: each vector becomes a hyperplane, and a query asks for the first hyperplane hit by a vertical ray from the query point. We use the data structure from [56], the runtime of which is captured in the following lemma (see [1] for a summary of ray-shooting data structures).

Lemma 5.3**.**

([56, Theorem 1.5]).* Let $\varepsilon>0$ be a constant, $k\geq 4$ a fixed dimension, and $m$ a parameter with $n\leq m\leq n^{\lfloor k/2\rfloor}$. Then, there is a dynamic data structure for ray-shooting queries with $O(m^{1+\varepsilon})$ space and preprocessing time, $O(m^{1+\varepsilon}/n)$ update time, and $O(\frac{n}{m^{1/\lfloor k/2\rfloor}}\log n)$ query time.*

Theorem 5.4**.**

For any $\varepsilon>0$, the narcissistic $k$-attribute stable matching problem can be solved in $O(n\log n)$ time for $k=2$, $O(n^{4/3+\varepsilon})$ time for $k=3$, and $O(n^{2-4/(k(1+\varepsilon)+2)})$ time for $k\geq 4$.

Proof 5.5**.**

Since the problem is symmetric, it can be solved in $O(P(n)+nT(n))$ time, given a dynamic data structure for ray-shooting queries with $P(n)$ preprocessing time and $T(n)$ operation time (Lemma 5.2).

For $k\geq 4$, using the data structure for ray-shooting queries from [56] (Lemma 5.3) results in a runtime of $O(m^{1+\varepsilon}+\frac{n^{2}\log n}{m^{1/\lfloor k/2\rfloor}})$ for any $\varepsilon>0$. The optimal runtime is achieved when the parameter $m$ is chosen to balance the two terms, i.e., so that $m^{1+\varepsilon}=\frac{n^{2}\log n}{m^{1/\lfloor k/2\rfloor}}$. This gives $m=(n^{2}\log n)^{1/(1+\varepsilon+1/\lfloor k/2\rfloor)}$. For the sake of obtaining a simple asymptotic expression, we set $m$ to $(n^{2}\log n)^{1/(1+\varepsilon+2/k)}$ (which is the same for even $k$, and bigger for odd $k$). Then, the $O(m^{1+\varepsilon})$ term dominates. Also note that if $\varepsilon<1-2/k$, this value of $m$ is between $n$ and $n^{\lfloor k/2\rfloor}$, so the condition in Lemma 5.3 is satisfied.

Thus, the problem can be solved in $O(m^{1+\varepsilon})=O((n^{2}\log n)^{(1+\varepsilon)/(1+\varepsilon+2/k)})$, which further simplifies to the claimed runtime of $O(n^{2-4/(k(1+\varepsilon^{\prime})+2)})$ (where $\varepsilon^{\prime}$ needs to satisfy $\varepsilon^{\prime}>\varepsilon$). For $k=3$, we use the same data structure, but raising the problem to four dimensions, so that Lemma 5.3 applies. For $k=2$, see Lemma 5.10.

Incidentally, the value for $m$ used in Theorem 5.4 also improves the algorithm by Künnemann et al. [58] for the one-sided $k$-attribute stable matching problem [58, Theorem 2], which also depends on the use of this data structure. The improvement is from $\tilde{O}(n^{2-1/(k/2)})$ to $O(n^{2-1/(k(1+\varepsilon)/4+1/2)})$. Similar balancing of preprocessing and query times in [55, Corollary 5.2] also improves the time to verify stability of a given matching in the (2-sided) $k$-attribute stable matching model [58, Section 5.1] for constant $k$; the improvement is from $\tilde{O}(n^{2-1/(2k)})$ to $O(n^{2-1/(k/2+1/2+\varepsilon)})$ for any $\varepsilon>0$.

5.3 The 2-attribute case

In this special case, we can design a simple first-choice data structure with $P(n)=O(n\log n)$ preprocessing time and $T(n)=O(\log n)$ operation time. Note that, for a vector $\vec{p}$ in $\mathbb{R}^{2}$, all the points along a line perpendicular to $\vec{p}$ are equally preferred, i.e., have the same dot product with $\vec{p}$ (because their projections onto the supporting line of $\vec{p}$ are the same). In fact, the preference list for $\vec{p}$ corresponds to the order in which a line perpendicular to $\vec{p}$ encounters the vectors in the other set as it moves in the direction opposite from $\vec{p}$ (see Figure 5, left). We get the following lemma (where the vectors in one set are interpreted as points).

Lemma 5.6**.**

Given a point set $P$ and a vector $\vec{q}$, in $\mathbb{R}^{2}$, the point $p^{*}$ in $P$ maximizing $\vec{q}\cdot p^{*}$ is in the convex hull of $P$.

Proof 5.7**.**

Consider a line perpendicular to $\vec{q}$. Move this line in the direction of $\vec{q}$, until all points in $P$ lie on the same side of it (behind it). Note that any line orthogonal to $\vec{q}$ has the property that all points lying on the line have the same dot product with $\vec{q}$. The point $p^{*}$ is the last point in $P$ to touch the line, since moving the line in the opposite direction from $\vec{q}$ decreases the dot product of $\vec{q}$ with any point on the line (and by the general position assumption, it is unique). Clearly, $p^{*}$ is in the convex hull.

Our first-choice data structure is a semi-dynamic convex hull data structure, where deletions are allowed but not insertions [47]. We handle queries as in Lemma 5.8.

Lemma 5.8**.**

Given the ordered list of points along the convex hull of a point set $P$, and a query vector $\vec{q}$, we can find the point $p^{*}$ in $P$ maximizing $\vec{q}\cdot p^{*}$ in $O(\log n)$ time, where $n$ is the number of points in the convex hull.

Proof 5.9**.**

By Lemma 5.6, the point $p^{*}$ is in the convex hull. For ease of exposition, assume that all the points in $P$ and $\vec{q}$ have positive coordinates (the alternative cases are similar). Then, $p^{*}$ lies in the top-right section of the convex hull (the section from the highest point to the rightmost point, in clockwise order).

Note that points along the top-right convex hull are ordered by their $y$-coordinate, so, we say above and below to describe the relative positions of points in it. Each point $p$ in the top-right convex hull is the point in $P$ maximizing $p\cdot\vec{q^{\prime}}$ for all the vectors $\vec{q^{\prime}}$ in an infinite wedge, as depicted in Figure 5, right. The wedge contains all the vectors $\vec{q^{\prime}}$ whose perpendicular line touches $p$ last when moving in the direction of $\vec{q^{\prime}}$, so the edges of the wedge are perpendicular to the edges of the convex hull incident to $p$. Thus, by looking at the neighbors of $p$ along the convex hull, we can calculate this wedge and know whether $\vec{q}$ is in the wedge for $p$, below it, or above it. Based on this, we discern whether the first choice of $\vec{q}$ is $p$ itself or above or below it. Thus, we can do binary search for $p^{*}$ in $O(\log n)$ time.

Lemma 5.10**.**

The narcissistic 2-attribute stable matching problem can be solved in $O(n\log n)$ time.

Proof 5.11**.**

We can use the NNC algorithm from [37] coupled with a first-choice data structure which is a semi-dynamic convex hull data structure. Updating the convex-hull can be done in $O(n\log n)$ time throughout the algorithm [47]. Queries are answered in $O(\log n)$ time (Lemma 5.8). Thus, the total running time is $O(n\log n)$.

6 Server Cover

Geometric coverage problems deal with finding optimal configurations of a set of geometric shapes that contain or “cover” another set of objects (for instance, see [3, 20, 65]). In this section, we propose an NNC-type algorithm for a problem in this category. We use NNC to speed up a greedy algorithm for a one-dimensional version of a server cover problem: given the locations of $n$ clients and $m$ servers, which can be seen as houses and telecommunication towers, the goal is to assign a “signal stregth” to each communication tower so that they reach all the houses, minimizing the cost of transmitting the signals.

Formally, we are given two sets of points in $\mathbb{R}^{\delta}$, $S$ (servers) and $C$ (clients). The problem is to assign a radius, $r_{i}$, to a disk centered at each server $s_{i}$ in $S$, so that every client is contained in at least one disk. The optimization function to minimize is $\sum r_{i}^{\alpha}$ for some parameter $\alpha>0$. The values $\alpha=1$ and $\alpha=2$ are of special interest, as they correspond to minimizing the sum of radii and areas (in 2D), respectively.

6.1 Related work

Table 1 gives an overview of exact and approximation algorithms for the server cover problem. It shows that when either the dimension $\delta$ or $\alpha$ are larger than $1$, there is a steep increase in complexity. We focus on the case with $\delta=1$ and $\alpha=1$, which has received significant attention because it gives insight into the problem in higher dimensions.

Server coverage was first considered in the one-dimensional setting by Lev-Tov and Peleg [54]. They gave an $O((n+m)^{3})$-time dynamic-programming algorithm for the $\alpha=1$ case, where $n$ is the number of clients and $m$ is the number of servers. They also gave a linear-time 4-approximation (assuming a sorted input). The runtime of the exact algorithm was improved to $O((n+m)^{2})$ by Biniaz et al. [17]. In the approximation setting, Alt et al. [6] gave a linear-time 3-approximation and an $O(m+n\log m)$-time $2$-approximation (also assuming a sorted input). Using NNC, we improve this to a linear-time $2$-approximation algorithm under the same assumption that the input is sorted.

6.2 Global-local equivalence in server cover

The $O(m+n\log m)$-time $2$-approximation by Alt et al. [6] can be described as follows: start with disks (which, in 1D, are intervals) of radius [math], and, at each step, make the smallest disk growth which covers a new client. If we define the distance $d(c,s)$ between a client $c$ and a server $s$ with a disk with radius $r\geq 0$ as the distance between $c$ and the closest boundary of the server’s disk, the process can be described as repeatedly finding the closest uncovered client–server pair and growing the server’s disk up to the client. Under this view, there is a natural notion of MNNs: an uncovered client $c$ and a server $s$ such that $d(c,s)$ is the smallest among all the distances involving $c$ and $s$.

However, Figure 6 illustrates that this problem does not satisfy global-local equivalence: matching MNNs does not yield the same result as matching the closest pair. Furthermore, it shows that matching MNNs loses the 2-approximation guarantee. We nevertheless use NNC to achieve a 2-approximation, which requires enhancing the algorithm so that it does not simply match MNNs. This shows that NNC may be useful even in problems where global-local equivalence does not hold.

6.3 Linear-time 2-approximation in 1D

The algorithm takes a list of $n$ clients and $m$ servers ordered left-to-right, and outputs a radius for each server (which might be [math]). In the algorithm, we group clients and servers into clusters. Each element starts as a base cluster, and we repeatedly merge them until there is a single cluster left. We distinguish between client clusters, consisting of a set of still uncovered clients, and server clusters, consisting of servers and covered clients. Clusters span intervals (as defined below). The distance $d(a,b)$ between clusters is defined as the distance between the closest endpoints of the clusters’ intervals. We begin by describing the merging operation based on the cluster types.

• •

We merge client clusters into larger client clusters. All the clients in a cluster are eventually covered together, so we only need to keep track of the left-most one, $p$ and right-most one, $q$; thus, we represent the client cluster with the interval $[p,q]$. Each client $p$ starts as a cluster $[p,p]$. Two client clusters $[p,q]$ and $[p^{\prime},q^{\prime}]$ (which in the algorithm never overlap), with $q<p^{\prime}$ are merged into a client cluster $[p,q^{\prime}]$.

• •

We merge server clusters into larger server clusters. Of all the servers in a cluster, only the ones with disks reaching furthest to the left and to the right may cover new clients. Let these servers be $s_{l}$ and $s_{r}$, respectively (which might be the same), let $l$ be the left-most point covered by $s_{l}$, and $r$ the right-most point covered by $s_{r}$. Then, all the information we need about a server cluster is $([l,r],s_{l},s_{r})$. Note that $l\leq s_{l}\leq s_{r}\leq r$. Each server $s$ starts as a cluster $([s,s],s,s)$. To merge two server clusters $([p,q],s_{p},s_{q})$ and $([p^{\prime},q^{\prime}],s_{p^{\prime}},s_{q^{\prime}})$ (which may overlap), let $l^{*}=\min{(p,p^{\prime})}$ and $r^{*}=\max{(q,q^{\prime})}$. Replace both by a server cluster $([l^{*},r^{*}],s_{l^{*}},s_{r^{*}})$. Retain the identities only of the two servers whose boundaries extend furthest left ($s_{l^{*}}$) and right ($s_{r^{*}}$).

• •

Merging a client cluster $[p,q]$ and a server cluster $([p^{\prime},q^{\prime}],s_{p^{\prime}},s_{q^{\prime}})$ (which may overlap) into a new server cluster involves covering all the clients in the cluster by $s\in\{s_{p^{\prime}},s_{q^{\prime}}\}$, whichever is cheaper. That is, Let $d^{*}$ be the new radius of the disk of $s$ after it grows to cover $[p,q]$; we merge the client cluster and the server cluster into a server cluster $([l^{*},r^{*}],s_{l^{*}},s_{r^{*}})$, where $l^{*}=\min{(p^{\prime},s-d^{*})}$, $r^{*}=\max{(q^{\prime},s+d^{*})}$, and $s_{l^{*}}$ (resp. $s_{r^{*}}$) is the server among $s$ and $s_{p^{\prime}}$ (resp. $s$ and $s_{q^{\prime}}$) with the leftmost (resp. rightmost) extending disk.

The algorithm works by building a chain (a stack) of clusters ordered from left to right. The following invariant holds at the beginning of each iteration: no two clusters overlap, the chain contains a prefix of the list of clusters, and the distance between successive clusters in the chain decreases. In the pseudocode (Algorithm 5) we use $a\cup b$ to denote the cluster resulting from merging clusters $a$ and $b$, and $C(s)$ to denote the cluster containing a server $s$.

6.3.1 Correctness

At the end of the algorithm, all the clusters have been merged into one, which is a server cluster (as long as there is one). Thus every client cluster has been merged with a server cluster, which means that some server grew its radius to cover the clients (or they became covered indirectly through an expansion). Thus, the output is a valid solution. We turn our attention to the analysis of the runtime and of the $2$-approximation factor. Throughout, we make an assumption that there are no ties between distances (or that they are broken consistently).

Lemma 6.1**.**

Algorithm 5 runs in $O(n+m)$ time, assuming the input is given in sorted order.

Proof 6.2**.**

Initially, there are $n+m$ clusters. Each merge operation reduces the number of clusters by one, so the number of merge operations is $n+m-1$. A merge can be done in $O(1)$ time, so the total time spent doing merges is $O(m+n)$. Each iteration of the main loop either causes at least one merge or adds a new cluster to the chain. Since clusters stay in the chain until they are merged, this can only happen $O(n+m)$ times, so there are $O(n+m)$ iterations. In 1D, finding the NN of a cluster simply involves checking the previous and next clusters, which can be done in $O(1)$.

Given an arbitrary problem instance, let NNC denote the solution output by Algorithm 5 and OPT denote an optimal solution.

Theorem 6.3**.**

$\mbox{cost}(NNC)\leq 2\mbox{cost}(OPT)$.

We follow the proof idea for the greedy algorithm from [6]. We “charge” the disk radii in NNC to disjoint “coverage intervals”, $I_{c}$, each of which is associated with a client $c$, and such that $\sum|I_{c}|=\mbox{cost}(\mbox{NNC})$, where $|I|$ denotes the length of an interval $I$. If the union of these intervals (and therefore the sum of their lengths) were entirely contained within the disks in OPT, then NNC would trivially be at most double the sum of radii in OPT. If that is not the case, we show that the length of every coverage interval outside of the disks in OPT is accounted for by an equal or greater absence of coverage intervals inside an OPT disk.

Definition 6.4**.**

Suppose that a server cluster $S$ and a client cluster $C$ are merged in Algorithm 5, and, as a result, server $s$ is expanded to cover clients $c_{1},\ldots,c_{k}$, in order of proximity to $s$. If $C$ is to the right of $s$, define the coverage interval $I_{c_{1}}$ as the open interval $(s_{b},c_{1})$, where $s_{b}$ is the right-most boundary of the disk of $s$ before the expansion, and define $I_{c_{i}}$ as $(c_{i},c_{i-1})$ for $1<i\leq k$. If $C$ is to the left, the intervals are defined symmetrically.

See Figure 7, bottom, for an example of the coverage intervals. To prove Theorem 6.3, we need the following intermediate results.

Lemma 6.5**.**

The intervals $I_{c},I_{c^{\prime}}$ are disjoint if $c\not=c^{\prime}$.

Proof 6.6**.**

If $c$ and $c^{\prime}$ belong to the same client cluster at the time $c$ is covered, it follows from the definition. Otherwise, let the cluster of $c$ be the one merged with a server cluster first. Then, after $c$ is covered, the interval $I_{c}$, if it exists, is inside a server cluster. Coverage intervals from clients covered later do not intersect existing server clusters.

For a server $s$, let $D_{O}(s)$ denote the disk of $s$ in OPT. To offset the intervals $I_{c}$ which occur outside of OPT disks, we need the following.

Remark 6.7**.**

Every $I_{c}$ intersects or has a shared endpoint with a disk $D_{O}(s)$.

This is because $c$ must be covered by OPT.

Suppose that for some $c$, $I_{c}$ is not contained in any disk in OPT. Then, by Remark 6.7, $I_{c}$ intersects or has a shared endpoint with a disk $D_{O}(s)$ in OPT. We consider the two possible cases separately, where $s$ is to the left or to the right of $c$. We show (Lemma 6.8) that in either case there is an interval $J$, between $s$ and $c$ and inside $D_{O}(s)$, which is disjoint from all coverage intervals (Fig. 8). Note that at most one $I_{c}$ may intersect $D_{O}(s)$ on each side on $s$, so the intervals $J$ do not overlap.

Lemma 6.8**.**

If a client $c$ belongs to $D_{O}(s)$ for some server $s$ and $I_{c}$ extends across the right (left) boundary of $D_{O}(s)$, then there is an interval $J$ in $D_{O}(s)$, to the right (left) of $s$, free of coverage intervals, and such that $|J|>|I_{c}|$.

Proof 6.9**.**

Right case. Consider first the setting in Figure 8, right: suppose that in OPT, $s$ covers some clients which in NNC are covered for the first time (the time where their coverage intervals are defined) from a server to the right of $D_{O}(s)$. Let $c_{1},\ldots,c_{k}$, $k\geq 1$, be all such clients. Then, the coverage interval $I_{c_{k}}$ of $c_{k}$ extends across the right boundary of $D_{O}(s)$. Let $x$ be the input element (client or server, possibly $s$) immediately to the left of $c_{1}$, and $y$ the input element immediately to the right of $c_{k}$. Note that $d(c_{k},y)\geq|I_{c_{k}}|$, since a coverage interval cannot extend past another input element.

We show that (i) $d(x,c_{1})>d(c_{k},y)$ (and thus, $d(x,c_{1})>|I_{c_{k}}|$), and that (ii) the interval $(x,c_{1})$ is free of coverage intervals. Claim (i) follows from the fact that if $d(x,c_{1})<d(c_{k},y)$, then $x$ and $c_{1}$ would be merged before $y$ is added to the chain, which cannot happen: if $x$ is a server, then $c_{1}$ would be covered from the left, and if $x$ is a client, $x$ and $c_{1}$ would be merged together, contradicting that $c_{1}$ is the left-most client covered from a server to the right of $D_{O}(s)$. For (ii), note that $c_{1}$ was covered from the right (by definition) and $x$ either was a client covered from the left (by definition of $c_{1}$) or a server which does not cover $c_{1}$. In the former case, $I_{x}$ has its right endpoint at $x$, and in the latter case, $x$ is not the client to cover $c_{1}$. Thus, there are no coverage intervals in $(x,c_{1})$.

Left case. Now consider the setting in Figure 8, left. The setting is similar, except that $c_{1},\ldots,c_{k}$, $k\geq 1$, are to the left of $s$ and are covered by a server to the left of $D_{O}(s)$, and it is the interval $I_{c_{1}}$ that extends across the left boundary of $D_{O}(s)$. Define $x$ as in the previous case but symmetrically: it is the input element immediately to the right of $c_{k}$. We define $y$ slightly differently: it is the right boundary of the cluster preceding $c_{1}$ at the time $c_{1}$ is added to the chain (not necessarily an input element). Note that $d(y,c_{1})\geq|I_{c_{1}}|$, since a coverage interval for $c_{1}$ would start at, or to the right of $y$.

Let $u$ and $v$ be the two consecutive elements among $c_{1},\ldots,c_{k},x$ maximizing $d(u,v)$ (note that $x$ may be a server). We show that (i) $d(u,v)>d(y,c_{1})$ (and thus, $d(u,v)>|I_{c_{1}}|$), and (ii) $(u,v)$ is free of coverage intervals. For (i), assume for a contradiction that $d(u,v)<d(y,c_{1})$. Then, $c_{1}$ and $x$ (and all the elements in between) would be clustered together before they are merged with $y$, because $y$ would not be the NN of $c_{1}$ until all these merges between closer elements happen. However, this contradicts that $c_{k}$ is the right-most client covered for the first time from the left of $D_{O}(s)$. Therefore, we have (i).

For (ii), note that if $v$ is not $x$, then $v$ and $x$ (and all the elements in between) get clustered together before they are merged with $u$, because $u$ would not be the NN of $v$ until all these merges between closer elements happen. However, this contradicts that $c_{k}$ is the right-most client covered (for the first time) from the left of $D_{O}(s)$. Therefore, $v$ is $x$, and $c_{k}$ is $u$. (ii) now follows by an analogous reasoning as in the symmetric case.

We are ready to prove Theorem 6.3.

Proof 6.10** (Proof of Theorem 6.3).**

As mentioned, $\mbox{cost}(\mbox{NNC})=\sum_{c}I_{c}$. The total length of the coverage intervals contained in OPT disks does not exceed twice the sum of the OPT radii (recall that by Lemma 6.5, the coverage intervals are pairwise-disjoint). Consider now the parts of the coverage intervals outside the disks of OPT. Note that for each OPT disk, there is at most one interval $I_{c}$ overlapping the disk on each side, and by Remark 6.7, every $I_{c}$ touches or overlaps a disk. Furthermore, by Lemma 6.8, for every length $\ell$ of coverage intervals in NNC outside the OPT disk of a server $s$ to the right (left) of $s$, there is at least $\ell$ length within $s$’s disk to the right (left) of $s$ that is free of coverage intervals. Therefore the approximation ratio of 2 is preserved.

See Figure 9 for an instance that shows that the 2-approximation is tight.

6.4 Greedy in higher dimensions

As mentioned, the greedy algorithm which makes the smallest disk growth, at each step, which covers a new client, achieves a $2$-approximation in the 1D setting [6]. Does it achieve a good approximation ratio in higher dimensions? In this section, we give a negative answer. It performs poorly in two dimensions, even when servers are constrained to lie on a line (also known as the $1.5$D case), and for $\alpha=1$. We show that in the instance illustrated in Figure 10, left, Greedy is a factor of $2m/\sqrt{5}$ worse than the optimal solution.

In this instance, a set of $m$ servers are placed along a horizontal line with a distance of $1$ between each consecutive pair. Above each server, we place a “column” of clients stretching up to distance $m$ above the servers. The clients in a column are evenly spaced and at distance $d=\sqrt{m^{2}+1}-m$ of each other. Thus, there are $m/d=m(m+\sqrt{m^{2}+1})$ clients in each column. The total number of clients is roughly $2m^{3}$.

Lemma 6.11**.**

Greedy’s approximation ratio in the 1.5D setting with $\alpha=1$ is no better than $2m/\sqrt{5}$.

Proof 6.12**.**

Consider the instance described above and illustrated in Figure 10.

The optimal solution is to cover all clients with a single server located at the center. By the Pythagorean Theorem, the cost of the optimal solution is $\sqrt{5}m/2$. In contrast, we show that Greedy would choose to cover the clients in each column by the server at the bottom of it, resulting in a cost of $m^{2}$. Thus, Greedy is $2m/\sqrt{5}$ times worse than the optimal solution.

We assume that Greedy breaks ties by choosing clients closer to the horizontal line first (alternatively, we can perturb the positions of the clients slightly to guarantee this tie breaking). Then, we can show that Greedy covers the clients by “layers”, where a layer is the set of clients at a given height. To see this, assume, for the sake of an inductive argument, that Greedy has grown each disk to cover the clients up to a given layer. Then, some server $s$ expands to cover a client $p$ in next layer, which would be the one in the same column. Let $s^{\prime}$ and $p^{\prime}$ be the server and client next to $s$ and $p$, respectively. We must argue that the disk of $s$ is not closer to $p^{\prime}$ than the disk of $s^{\prime}$. Note that it suffices to show that this does not happen when $p$ is the very last client in the column above $s$. This is because the further away $p$ is from $s$, the bigger the radius of the disk of $s$, resulting in a disk closer to $p^{\prime}$. This is illustrated in Figure 10, right. The distance $d$ between two consecutive points in a column is chosen precisely so that, in this scenario where $p$ is the last point, the disk of $s$ is exactly as close to $p^{\prime}$ as the disk of $s^{\prime}$. Depending on the tie-breaking rule (or changing $d$ to be slightly smaller), $s^{\prime}$ will grow to cover $p^{\prime}$ and not $s$.

7 Conclusions

Before this paper, NNC had been used only in agglomerative hierarchical clustering and stable matching problems based on proximity. This paper adds the following use cases:

• •

Its first use in problems without symmetric distances (motorcycle graphs, Section 4). We showed that the chain still works for finding nearest-neighbor cycles.

• •

Its first use without a NN structure (multi-fragment TSP, which uses our new soft NN structure, Section 3).

• •

Its first use in problems without global-local equivalence (server cover, Section 6). We showed that the chain can be adapted in settings where matching MNNs is not as good as matching overall closest pairs.

• •

Its first use in approximation algorithms (also server cover).

• •

Its first use in a graph-theoretical framework (Steiner TSP, Section 3.3).

• •

Its first use in stable matching problems not based on distances (narcissistic $k$-attribute stable matching, Section 5).

The above applications illustrate the two main points of this paper: first, that in several geometric problems, finding mutually nearest neighbors leads to the same solution as finding closest pairs, which we call global-local equivalence; second, that MNNs can be found efficiently thanks to the NNC algorithm.

We expect that this algorithm will find more uses in computational geometry. When dealing with problems involving nearest neighbors or closest pairs in some way, one may check if a form of global-local equivalence holds. If so, one should then consider using the NNC algorithm. The main guidelines for designing NNC algorithms are: (1) each link in the chain should get closer to MNNs; (2) to avoid infinite loops, the chain should be acyclic. One should be careful to break ties consistently; (3) after finding and processing MNNs, all the previous links in the chain should remain valid. These simple ingredients are likely to lead to an algorithm with a runtime of the form $O(P(n)+nT(n))$, as seen throughout this paper.

We conclude with some open questions.

• •

Can we use specialized data structures? Throughout the paper, we have used fully dynamic data structures that allow insertions and deletions. However, NNC algorithms typically only use deletions. Further, query points are generally known beforehand. Therefore, specialized data structures with these considerations in mind may speed up the algorithms in this paper, and NNC algorithms in general.

• •

Does the SNN structure have more uses? The SNN structure has a curious type of queries, where soft answers are not directly related to the query point. Nonetheless, we showed how it can be used as part of the SNNC algorithm and to solve the closest pair problem.

• •

Are other narcissistic stable matching models symmetric? As mentioned, narcissistic is a descriptive term for stable matching models where the preferences of each agent are determined by the agent’s own qualities or attributes. We find it likely that NNC can be used for other narcissistic models, as such preferences seem unlikely to create cycles of first choices.