Computing the Greedy Spanner in Linear Space

The greedy spanner is a high-quality spanner: its total weight, edge count and maximal degree are asymptotically optimal and in practice significantly better than for any other spanner with reasonable construction time. Unfortunately, all known algorithms that compute the greedy spanner of n points use Omega(n^2) space, which is impractical on large instances. To the best of our knowledge, the largest instance for which the greedy spanner was computed so far has about 13,000 vertices. We present a O(n)-space algorithm that computes the same spanner for points in R^d running in O(n^2 log^2 n) time for any fixed stretch factor and dimension. We discuss and evaluate a number of optimizations to its running time, which allowed us to compute the greedy spanner on a graph with a million vertices. To our knowledge, this is also the first algorithm for the greedy spanner with a near-quadratic running time guarantee that has actually been implemented.


Introduction
A t-spanner on a set of points, usually in the Euclidean plane, is a graph on these points that is a 't-approximation' of the complete graph, in the sense that shortest routes in the graph are at most t times longer than the direct geometric distance. The spanners considered in literature have only O(n) edges as opposed to the O(n 2 ) edges in the complete graph, or other desirable properties such as bounded diameter or bounded degree, which makes them a lot more pleasant to work with than the complete graph.
Spanners are used in wireless network design [6]: for example, high-degree routing points in such networks tend to have problems with interference, so using a spanner with bounded degree as network avoids these problems while maintaining connectivity. They are also used as a component in various other geometric algorithms, and are used in distributed algorithms. Spanners were introduced in network design [11] and geometry [4], and have since been subject to a considerable amount of research [8,10].
There exists a large number of constructions of t-spanners that can be parameterized with arbitrary t > 1. They have different strengths and weaknesses: some are fast to construct but of low quality (Θ-graph, which has no guarantees on its total weight), others are slow to construct but of high quality (greedy spanner, which has low total weight and maximum degree), some have an extremely low diameter (various dumbbell based constructions) and some are fast to construct in higher dimensions (well-separated pair decomposition spanners). See for example [10] for detailed expositions of these spanners and their properties.
The greedy spanner is one of the first spanner algorithms that was considered, and it has been subject to a considerable amount of research regarding its properties and more recently also regarding computing it efficiently. This line of research resulted in a O(n 2 log n) algorithm [1] for metric spaces of bounded doubling dimension (and therefore also for Euclidean spaces). There is also an algorithm with O(n 3 log n) worst-case running time that works well in practice [5]. Its running time tends to be near-quadratic in practical cases, but there are examples on which its running time is Θ(n 3 log n). Its space usage is Θ(n 2 ).
Among the many spanner algorithms known, the greedy spanner is of special interest because of its exceptional quality: its size, weight and degree are asymptotically optimal, and also in practice are better than any other spanner construction algorithms with reasonable running times. For example, it produces spanners with about ten times as few edges, twenty times smaller total weight and six times smaller maximum degree as its closest well-known competitor, the Θ-graph, on uniform point sets. The contrast is clear in Fig. 1. Therefore, a method of computing it more efficiently is of considerable interest.
We present an algorithm whose space usage is Θ(n) whereas existing algorithms use Θ(n 2 ) space, while being only a logarithmic factor slower than the fastest known algorithm, thus answering a question left open in [1]. Our algorithm makes the greedy spanner practical to compute for much larger inputs than before: this used to be infeasible on graphs of over 15,000 vertices. In contrast, we tested our algorithm on instances of up to 1,000,000 points, for which previous algorithms would require multiple terabytes of memory. Furthermore, with the help of several optimizations we will present, the algorithm is also fast in practice, as our experiments show.
The method used to achieve this consists of two parts: a framework that uses linear space and near-linear time, and a subroutine using linear space and nearlinear time, which is called a near-linear number of times by the framework. The subroutine solves the bichromatic closest pair with dilation larger than t problem. If there is an algorithm with a sublinear running time for this subproblem (possibly tailored to our specific scenario), then our framework immediately gives an asymptotically faster algorithm than is currently known. This situation is reminiscent to that of the minimum spanning tree, for which it is known that it is essentially equivalent to the bichromatic closest pair problem.
The rest of the paper is organized as follows. In Section 2 we review a number of well-known definitions, algorithms and results. In Section 3 we give the properties of the WSPD and the greedy spanner on which our algorithm is based. In Section 4 we present our algorithm and analyse its running time and space usage. In Section 5 we discuss our optimizations of the algorithm. Finally, in Section 6 we present our experimental results and compare it to other algorithms.

Notation and preliminaries
Let V be a set of points in R d , and let t ∈ R be the intended dilation (1 < t). Let G = (V, E) be a graph on V . For two points u, v ∈ V , we denote the Euclidean distance between u and v by |uv|, and the distance in G by δ G (u, v). If the graph G is clear from the context we will simply write δ(u, v). The dilation of a pair of points is t if δ(u, v) ≤ t · |uv|. A graph G has dilation t if t is an upper bound for the dilations of all pairs of points. In this case we say that G is a t-spanner. To simplify the analysis, we assume without loss of generality that t < 2.
We will often say that two points u, v ∈ V have a t-path if their dilation is t. A pair of points is without t-path if its dilation is not t. When we say a pair of points (u, v) is the closest or shortest pair among some set of points, we mean that |uv| is minimal among this set. We will talk about a Dijkstra computation from a point v by which we mean a single execution of the single-source shortest path algorithm known as Dijkstra's algorithm from v.
Consider the following algorithm that was introduced by Keil [9]: then add (u, v) to E 5. return E Obviously, the result of this algorithm is a t-spanner for V . The resulting graph is called the greedy spanner for V , for which we shall present a more efficient algorithm than the above.
We will make use of the Well-Separated Pair Decomposition, or WSPD for short, as introduced by Callahan and Kosaraju in [2,3]. A WSPD is parameterized with a separation constant s ∈ R with s > 0. This decomposition is a set of pairs of nonempty subsets of V . Let m be the number of pairs in a decomposition. We can number the pairs, and denote every pair as Let u and v be distinct points, then we say that (u, v) is 'in' a well-separated pair {A i , B i } if u ∈ A i and v ∈ B i or v ∈ A i and u ∈ B i . A decomposition has the property that for every pair of distinct points u and v, there is exactly one For two nonempty subsets X k and X l of V , we define min(X k , X l ) to be the shortest distance between the two circles around the bounding boxes of X k and X l and max(X k , X l ) to be the longest distance between these two circles. Let diam(X k ) be the diameter of the circle around the bounding box of X k . Let (X k , X l ) be the distance between the centers of these two circles, also named the length of this pair. For a given separation constant s ∈ R with s > 0 as parameter for the WSPD, we require that all pairs in a WSPD are s-well-separated, that is, min It is easy to see that max(X k , X l ) ≤ min(X k , X l ) + diam(X k ) + diam(X l ) ≤ (1 + 2/s) min(X k , X l ). As t < 2 and as we will pick s = 2t t−1 later on, we have s > 4, and hence max(X k , X l ) ≤ 3/2 min(X k , X l ). Similarly, (X k , X l ) ≤ min(X k , X l ) + diam(X k )/2 + diam(X l )/2 ≤ (1 + 1/s) min(X k , X l ) and hence (X k , X l ) ≤ 5/4 min(X k , X l ).
For any V and any s > 0, there exists a WSPD of size m = O(s d n) that can be computed in O(n log n + s d n) time and can be represented in O(s d n) space [2]. Note that the above four values (min, max, diam and ) can easily be precomputed for all pairs with no additional asymptotic overhead during the WSPD construction.

Properties of the greedy spanner and the WSPD
In this section we will give the idea behind the algorithm and present the properties of the greedy spanner and the WSPD that make it work. We assume we have a set of points V of size n, an intended dilation t with 1 < t < 2 and a WSPD with separation factor s = 2t t−1 , for which the pairs are numbered is the number of pairs in the WSPD.
The idea behind the algorithm is to change the original greedy algorithm to work on well-separated pairs rather than edges. We will end up adding the edges in the same order as the greedy spanner. We maintain a set of 'candidate' edges for every well-separated pair such that the shortest of these candidates is the next edge that needs be added. We then recompute a candidate for some of these wellseparated pairs. We use two requirements to decide on which pairs we perform a recomputation, that together ensure that we do not do too many recomputations, but also that we do not fail to update pairs which needed updating.
We now give the properties on which our algorithm is based.
Proof. This observation is not fully identical to [1] as our definition of wellseparatedness is slightly different than theirs. However, their proof uses only properties of their Lemma 2, which still holds true using our definitions as proven in 2. Their Lemma 2 is almost Lemma 9.1.2 in [10], whose definitions of wellseparatedness is near identical to ours, except that they use radii rather than diameters and hence have different constants.
An immediate corollary is: . The greedy spanner contains Proof. Assume without loss of generality that |A i | ≤ |B i |. We perform a Dijkstra computation for every point a ∈ A i , maintaining the closest point in |B i | such that its dilation with respect to a is larger than t over all these computations. To check whether a point that is considered by the Dijkstra computation is in |B i |, This is not stated explicitly in [2], but it is a direct consequence of the construction of the many-one realization: it is proven that the many-one realization consists of O(s d n log n) pairs, and the construction splits every pair in the canonical realization into min(|A i |, |B i |) pairs, so the lemma follows.
Observation 5. Let E be some edge set for V . Let (a, b) ∈ E. Let c ∈ V and d ∈ V be points such that |ac|, |ad|, |bc|, |bd| > t|cd|. Then any t-path between c and d will not use the edge (a, b).
Proof. This directly follows from the fact that c and d are so far away from a and b that just getting to either a or b is already longer than allowed for a t-path.
Proof. This follows easily from a very similar statement, Lemma 11.3.4 in [10]. We first note the differences between their definitions and ours. Their statement involves dumbbells, but these are really the same as our well-separated pairs. Also, their definition of well-separatedness is different from ours in that they use the radius and we the diameter. We can easily amend this by observing that a WSPD with separation factor s using our definitions is identical to a WSPD with separation factor 2s using their definitions. This means our constant c sγ is larger than their constant because our s is doubled, but this is asymptotically irrelevant. Now we shall see how our statement follows from theirs. Their interval is [ , 2 ] and they allow the length of , then we can obtain our fact by invoking their lemma twice, first by setting = (resulting in an upper bound on the number such pairs n 1 ) and then by setting = /2 (resulting in an upper bound on the number of such pairs n 2 ). This counts the number of pairs with lengths in the interval [ /2, 2 ]. These are exactly the pairs we are interested in, and hence the constant we obtain is n 1 + n 2 , thus proving the fact.
This concludes the theoretical foundations of the algorithm. We will now present the algorithm and analyze its running time.

Algorithm
We will now describe algorithm GreedySpanner in detail. It first computes the WSPD for V with s = 2t t−1 and sorts the resulting pairs according to their smallest distance min(A i , B i ). It then alternates between calling the FillQueue procedure that attempts to add well-separated pairs to a priority queue Q, and removing an element from Q and adding a corresponding edge to E. If Q is empty after a call to FillQueue, the algorithm terminates and returns E.
We assume we have a procedure ClosestPair (i) that for the ith well-separated pair computes the closest pair of points without t-path in the graph computed so far, as presented in Lemma 3, and returns this pair, or returns nil if no such pair exists. For the priority queue Q, we let min(Q) denote the value of the key of the minimum of Q. Recall that m = O(s d n) denotes the number of well-separated pairs in the WSPD that we compute in the algorithm.
We maintain an index i into the sorted list of well-separated pairs. It points to the smallest untreated well-separated pair -we treat the pairs in order of min(A i , B i ) in the FillQueue procedure. When we treat a pair {A i , B i }, we call ClosestPair (i) on it, and if it returns a pair (u, v), we add it to Q with key |uv|. We link entries in the queue, its corresponding pair {A i , B i } and (u, v) together so they can quickly be requested. We stop treating pairs and return from the procedure if we have either treated all pairs, or if min(A i , B i ) is larger than the key of the minimal entry in Q (if it exists).
After extracting a pair of points (u, v) from Q, we add it to E. Then, we update the information in Q: for every pair {A j , B j } having an entry in Q for which either bounding box is at most t|uv| away from {A i , B i }, we recompute ClosestPair (j) and updates its entry in Q as follows. If the recomputation returns nil, we remove its entry from Q. If it returns a pair (u , v ), we link the entry of j in Q with this new pair and we increase the key of its entry to |u v |.
if p is not nil, but a pair (u, v) 4.
then add (u, v) to Q with key |uv|, and associate this entry with Algorithm GreedySpanner (V, t) 1. Compute the WSPD W for V with s = 2t t−1 , and let {A i , B i } be the resulting pairs, 1 ≤ i ≤ m 2. Sort the pairs {A i , B i } according to min(A i , B i ) 3. E ← ∅ 4. Q ← an empty priority queue 5. i ← 1 6. FillQueue(Q, i) 7. while Q is not empty 8.
do extract the minimum from Q, let this be (u, v) 9.
add (u, v) to E 10.
for all pairs {A j , B j } with an entry in Q for which either bounding box is at most t|uv| away from either u or v 11.
if p is nil, remove the entry in Q associated with {A j , B j } from Q 13.
if p is a pair (u , v ), update the entry in Q associated with {A j , B j } to contain (u , v ) and increase its key to |u v | 14.
FillQueue(Q, i) 15. return E We now prove correctness and a bound on the running time of the algorithm.
Theorem 7. Algorithm GreedySpanner computes the greedy spanner for dilation t.
Proof. We will prove that if the algorithm adds (u, v) to E on line 9, then (u, v) is the closest pair of points without a t-path between them. The greedy spanner consists of exactly these edges and hence this is sufficient to prove the theorem.
It is obvious that if we call ClosestPair (i) on every well-separated pair and take the closest pair of the non-nil results, then that would be the closest pair of points without a t-path between them. Our algorithm essentially does this, except it does not recalculate ClosestPair (i) for every pair after every added edge, but only for specific pairs. We will prove that the calls it does not make do not change the values in Q. Our first optimization is that if a call ClosestPair (i) returns nil it will always return nil, so we need not call ClosestPair (i) again, which is therefore a valid optimization.
The restriction 'for which either bounding box is at most t|uv| away from either u or v' from the for-loop on line 10 is the second optimization. Its validity is a direct consequence of Observation 5: all pairs of points in such well-separated pairs are too far away to use the newly-added edge to gain a t-path. Therefore re-running ClosestPair (i) and performing lines 12 and 13 will not change any entries in Q as claimed.
As min(A i B i ) is a lower bound on the minimal distance between any two points (a, b) in {A i , B i }, it immediately follows that calling FillQueue(Q, i) on a pair {A i , B i } with min(A i B i ) > min(Q) cannot possibly yield a pair that can cause min(Q) to become smaller. As the pairs are treated in order of min(A i B i ), this means the optimization that is the condition on line 1 in FillQueue(Q, i) is a valid optimization. This proves the theorem.
We will now analyze the running time and space usage of the algorithm. We will use the observations in Section 3 to bound the amount of work done by the algorithm.

Lemma 8. For any well-separated pair
Proof. ClosestPair (i) is called once for every i in the FillQueue procedure. Clos-estPair (i) may also be called after an edge is added to the graph. We will show that if a well-separated pair {A j , B j } causes ClosestPair (i) to be called, then Then, by the condition of line 10, the collection of pairs that call ClosestPair (i) satisfy the requirements of Fact 6 by setting γ = t, so we can conclude this happens only c st times. The lemma follows.
We will now show that (A j , B j ) ∈ [ (A i , B i )/2, 2 (A i , B i )]. Recall the following from Section 2: The algorithm ensures the following: Combining these we have: (1 + c st ) min(|A i |, |B i |) n log n + 1 (t − 1) d n and its space usage is O(n) by reusing the space for the calls. Using Fact 4, this is at most The time taken by all other steps of the algorithm is insignificant compared to the time used by ClosestPair (i) calls. These other steps are: computing the WSPD and all actions with regard to the queue. All these other actions use O 1 (t−1) d n space. Combining this with Theorem 7, the theorem follows.

Making the algorithm practical
Experiments suggested that implementing the above algorithm as-is does not yield a practical algorithm. With the four optimizations described in the following sections, the algorithm attains running times that are a small constant slower than the algorithm introduced in [5] called FG-greedy, which is considered the best practical algorithm known in literature.

Finding close-by pairs
The algorithm at some point needs to know which pairs are 'close' to the pair for which we are currently adding an edge. In our proof above, we suggested that these pairs be precomputed in O(m 2 ) time. Unfortunately, this precomputation step turns out to take much longer than the rest of the algorithm. If n = 100, then (on a uniform pointset) m ≈ 2000 and m 2 ≈ 4000000 while the number of edges e in the greedy spanner is about 135. Our solution is to simply find them using a linear search every time we need to know this information. This only takes O(e · m) time, which is significantly faster.

Reducing the number of Dijkstra computations
After decreasing the time taken by preprocessing, the next part that takes the most time are the Dijkstra computations, whose running time dwarfs the rest of the operations. We would therefore like to optimize this part of the algorithm. For every well-separated pair, we save the length of the shortest path found by any Dijkstra computation performed on it, that is, its source s, target t and distance δ(s, t). Then, if we are about to perform a Dijkstra computation on a vertex u, we first check if the saved path is already good enough to 'cover' all nodes in B i . Let c be the center of the circle around the bounding box of B i and r its radius. We check if t · |us| + δ(s, t) + t · (|tc| + r) ≤ t · (|uc| − r) and mark it as 'irrelevant for the rest of the algorithm'. This optimization roughly improves its running time by a factor three.

Sharpening the bound of Observation 5
The bound given in Observation 5 can be improved. Let {A i , B i } be the wellseparated pair for which we just added an edge and let {A j , B j } be the wellseparated pair under consideration in our linear search. First, some notation: let X k , X l be sets belonging to some well-separated pair (not necessarily the same pair), then min(X k , X l ) denotes the (shortest) distance between the two circles around the bounding boxes of X k and X l and max(X k , X l ) the longest distance between these two circles. Let = (A i , B i ). We can then replace the condition of Lemma 5 by the sharper condition min( The converse of the condition implies that the edge just added cannot be part of a t-path between a node in {A j , B j }, so the correctness of the algorithm is maintained. This leads to quite a speed increase.

Miscellaneous optimizations
There are two further small optimizations we have added to our implementation.
Firstly, rather than using the implicit linear space representation of the WSPD, we use the explicit representation where every node in the split tree stores the points associated with that node. For point sets where the ratio of the longest and the shortest distance is bounded by some polynomial in n, this uses O(n log n) space rather than O(n) space. This is true for all practical cases, which is why we used it in our implementation. For arbitrary point sets, this representation uses O(n 2 ) space. In practice, this extra space usage is hardly noticeable and it speeds up access to the points significantly.
Secondly, rather than performing Dijkstra's algorithm, we use the A * algorithm. This algorithm uses geometric estimates to the target to guide the computation to its goal, thus reducing the search space of the algorithm [7].
We have tried a number of additional optimizations, but none of them resulted in a speed increase. We describe them here.
We have tried to replace A * by ALT , a shortest path algorithm that uses landmarks -see [7] for details on ALT -which gives better lower bounds than the geometric estimates used in A * . However, this did not speed up the computations at all, while costing some amount of overhead.
We have also tried to further cut down on the number of Dijkstra computations. We again used that we store the lengths of the shortest paths found so far per well-separated pair. Every time after calling ClosestPair (i) we checked if the newly found path is 'good enough' for other well-separated pairs, that is, if the path combined with t-paths from the endpoints of the well-separated pairs would give t-paths for all pairs of points in the other well-separated pair. This decreased the number of Dijkstra computations performed considerably, but the overhead from doing this for all pairs was greater than its gain.
We tried to speed up the finding of close-by pairs by employing range trees. We also tried using the same range trees to perform the optimization of the previous paragraph only to well-separated pairs 'close by' our current well-separated pair. Both optimizations turned out to give a speed increase and in particular the second retained most of its effectiveness even though we only tried it on close-by pairs, but the overhead of range trees was vastly greater than the gain -in particular the space usage of range trees made the algorithm use about as much space as the original greedy algorithms.

Experimental results
We have run our algorithm on point sets whose size ranged from 100 to 1,000,000 vertices on several distributions. If the set contained at most 10,000 points, we have also run the FG-greedy algorithm to compare the two algorithms. We have recorded both space usage and running time (wall clock time). We have also performed a number of tests with decreasing values of t on datasets of size 10,000 and 50,000. Finally, as this is the first time we can compute the greedy spanner on large graphs, we have compared it to the Θ-graph and WSPD-based spanners on large instances.
We have used three kinds of distributions from which we draw our points: a uniform distribution, a gamma distribution with shape parameter 0.75, and a distribution consisting of √ n uniformly distributed pointsets of √ n uniformly distributed points. The results from the gamma distribution were nearly identical to those of the uniform pointset, so we did not include them. All our pointsets are two-dimensional.

Experiment environments
The algorithms have been implemented in C++. We have implemented all data structures not already in the std. The random generator used was the Mersenne Twister PRNG -we have used a C++ port by J. Bedaux of the C code by the designers of the algorithm, M. Matsumoto and T. Nishimura.
We have used two servers for the experiments. Most experiments have been run on the first server, which uses an Intel Core i5-3470 (3.20GHz) and 4GB (1600 MHz) RAM. It runs the Debian 6.0.7 OS and we compiled for 32 bits using G++ 4.7.2 with the -O3 option. For some tests we needed more memory, so we have used a second server. This server uses an Intel Core i7-3770k (3.50GHz) and 32 GB RAM. It runs Windows 8 Enterprise and we have compiled for 64 bits using the Microsoft C++ compiler (17.00.51106.1) with optimizations turned on.

Dependence on instance size
Our first set of tests compared FG-greedy and our algorithm for different values of n. The results are plotted in Fig. 2. As FG-greedy could only be ran on relatively small instances, its data points are difficult to see in the graph, so we added a zoomed-in plot for the bottom-left part of the plot.
We have used standard fitting methods to our data points: the running time of all algorithms involved fits a quadratic curve well, the memory usage of our algorithm is linear and the memory usage of FG-greedy is quadratic. This nicely fits our theoretical analysis. In fact, the constant factors seem to be much smaller than the bound we gave in our proof. We do note a lack of 'bumps' that are often occur when instance sizes start exceeding caches: this is probably due to the cache-unfriendly behavior of our algorithm and the still significant constant factor in our memory usage that will fill up caches quite quickly.
Comparing our algorithm to FG-greedy, it is clear that the memory usage of our algorithm is vastly superior. The plot puts into perspective just how much larger the instances are that we are able to deal with using our new algorithm compared to the old algorithms. Furthermore, our algorithm is about twice as fast as FG-greedy on the clustered datasets, and only about twice as slow on uniform datasets. On clustered datasets the number of computed well-separated pairs is much smaller than on uniform datasets so this difference does not surprise us. These plots suggest that our aim -roughly equal running times at vastly reduced space usage -is reached with this algorithm.

Dependence on t
We have tested our algorithms on datasets of 10,000 and 50,000 points, setting t to 1.1, 1.2, 1.4, 1.6, 1.8 and 2.0 to test the effect of this parameter. The effects of the parameter ended up being rather different between the uniform and clustered datasets.  On uniform pointsets, see Figures 3 and 4, our algorithm is about as fast as FG-greedy when t = 2, but its performance degrades quite rapidly as t decreases compared to FG-greedy. A hint to this behavior is given by the memory usage of our algorithm: it starts vastly better but as t decreases it becomes only twice as good as FG-greedy. This suggests that the number of well-separated pairs grows rapidly as t decreases, which explains the running time decrease.  On clustered pointsets, see Figures 5 and 6, the algorithms compare very differently. FG-greedy starts out twice as slow as our algorithm when t = 2 and when t = 1.1, our algorithm is only slightly faster than FG-greedy. The memory usage of our algorithm is much less dramatic than in the uniform point case: it hardly grows with t and therefore stays much smaller than FG-greedy. The memory usage of FG-greedy only depends on the number of points and not on t or the distribution of the points, so its memory usage is the same.

Comparison with other spanners
We have computed the greedy spanner on the instance shown in Fig. 1, which has 115,475 points. On this instance the greedy spanner for t = 2 has 171,456 edges, a maximum degree of 5 and a weight of 11,086,417. On the same instance, the Θ-graph with k = 6 has 465,230 edges, a maximum degree of 62 and a weight of 53,341,205. The WSPD-based spanner has 16,636,489 edges, a maximum degree of 1,271 and a weight of 20,330,194,426.
As shown in Fig. 2, we have computed the greedy spanner on a pointset of 500,000 uniformly distributed points. On this instance the greedy spanner for t = 2 has 720,850 edges, a maximum degree of 6 and a weight of 9,104,690. On the same instance, the Θ-graph with k = 6 has 2,063,164 edges, a maximum degree of 22 and a weight of 39,153,380. We were unable to run the WSPD-based spanner algorithm on this pointset due to its memory usage.
As shown in Fig. 2, we have computed the greedy spanner on a pointset of 1,000,000 clustered points. On this instance the greedy spanner for t = 2 has 1,409,946 edges, a maximum degree of 6 and a weight of 4,236,016. On the same instance, the Θ-graph with k = 6 has 4,157,016 edges, a maximum degree of 135 and a weight of 59,643,264. We were unable to run the WSPD-based spanner algorithm on this pointset due to its memory usage.
We have computed the greedy spanner on a pointset of 50,000 uniformly distributed points with t = 1.1. On this instance the greedy spanner has 225,705 edges, a maximum degree of 18 and a weight of 15,862,195. On the same instance, the Θ-graph with k = 73 (which is the smallest k for which a guarantee of t = 1.1 has been proven to our knowledge) has 2,396,361 edges, a maximum degree of 146 and a weight of 495,332,746. We were unable to run the WSPD-based spanner algorithm on this pointset with t = 1.1 due to its memory usage.
These results show that the greedy spanner really is an excellent spanner, even on large instances and for low t, as predicted by its theoretical properties.

Conclusion
We have presented an algorithm that computes the greedy spanner in Euclidean space in O(n 2 log 2 n) time and O(n) space for any fixed stretch factor and dimension. Our algorithm avoids computing all distances by considering well-separated pairs instead. It consists of a framework that computes the greedy spanner given a subroutine for a bichromatic closest pair problem. We give such a subroutine which leads to the desired result.
We have presented several optimizations to the algorithm. Our experimental results show that these optimizations make our algorithm have a running time close to the fastest known algorithms for the greedy spanner, while massively decreasing space usage. It allowed us to compute the greedy spanner on very large instances of a million points, compared to the earlier instances of at most 13,000 points. Given that our algorithm is the first algorithm with a near-quadratic running time guarantee that has actually been implemented, that it has linear space usage and that its running time is comparable to the best known algorithms, we think our algorithm is the method of choice to compute greedy spanners.
We leave open the problem of providing a faster subroutine for solving the bichromatic closest pair with dilation larger than t problem in our framework, which may allow the greedy spanner to be computed in subquadratic time. Particularly the case of the Euclidean plane seems interesting, as the closely related 'ordinary' bichromatic closest pair problem can be solved quickly in this setting.