Scalable Distributed String Sorting
Florian Kurpicz Pascal Mehnert Peter Sanders Matthias Schimek


Topological properties of finite-size heterostructures of magnetic topological insulators and superconductors
    Thomas L. Schmidt
    October 2023
=============================================================================================================




	String sorting is an important part of tasks such as building index data structures.
	Unfortunately, current string sorting algorithms do not scale to massively parallel distributed-memory machines since they either have latency (at least) proportional to the number of processors p or communicate the data a large number of times (at least logarithmic).
	We present practical and efficient algorithms for distributed-memory string sorting that scale to large p.
	Similar to state-of-the-art sorters for atomic objects, the algorithms have latency of about p^1/k when allowing the data to be communicated k times.
	Experiments indicate good scaling behavior on a wide range of inputs on up to 49152 cores. Overall, we achieve speedups of up to 5 over the current state-of-the-art distributed string sorting algorithms.




§ INTRODUCTION

Sorting strings is a fundamental building block of many important string-processing tasks such as the construction of index data structures for databases and full-text phrase search <cit.>.
The problem differs from atomic sorting—where keys are treated as indivisible objects that can be compared in constant time.
Strings on the other hand can have variable lengths and the time needed to compare two strings depends on the length of their longest common prefix.
Therefore, string sorting algorithms try to avoid the comparison of whole strings.
Instead, they only inspect the distinguishing prefixes of the strings, i.e., the characters needed to establish the global ordering, only once.
The sum of the lengths of all distinguishing prefixes is usually denoted by D.
The lower bound for sequential string sorting based on character comparisons is Ω(nlog n + D) with existing algorithms matching this bound <cit.>.





    Related Work.

There exists extensive research on string sorting in the sequential setting.
For a systematic overview, we refer to <cit.>.
We focus on parallel sorting algorithms.
Let n be the total number of strings and N be the total number of characters.
See <ref> for a list of abbreviations used in this paper.
For the PRAM model, there are (comparison-based) algorithms solving the string sorting problem in 𝒪(nlog n + N) work and 𝒪(log^2 n/loglog n) time <cit.>.
For integer alphabets, Hagerup proposes an algorithm with 𝒪(Nlog N) work and running time in 𝒪(log N/loglog N) <cit.>.
Ellert et al. propose a framework for PRAM string sorting algorithms that makes the (work) complexity depend on D instead of N by increasing the time complexity by a logarithmic factor in the length of the longest distinguishing prefix <cit.>.
String sorting is also considered in other parallel settings, e.g., on the GPU <cit.>.

While the problem has been extensively studied in the sequential and (shared-memory) parallel setting, we are only aware of the following results in the distributed-memory setting.
Bingmann et al. present the state-of-the-art distributed string sorting algorithms <cit.>.
The first one follows the standard distributed-memory merge sort scheme (local sorting, partitioning, message exchange, and merging) <cit.>.
Every step is augmented with string-specific optimizations, e.g., LCP-compression and LCP-aware merging <cit.>, see <ref>.
The second algorithm is more communication-efficient, as it only sorts approximations of the distinguishing prefixes using the first algorithm.
They also adapt the distributed hypercube quicksort algorithm <cit.> to variable-length keys (without string-related optimizations).
These algorithms improve the first dedicated distributed string sorting algorithm by Fischer et al. <cit.>.

However, the algorithms are only efficient for very small or large inputs, as they have a prohibitively high communication volume (hypercube quicksort) or do not scale to the largest available machines due to their latency which is (at least) proportional to the number of processors p.
For example, the merge sort algorithms do not scale to the largest available machines, as the partitioning step dominates the running time unless n=Ω(p^2 log p).

theoremThm.Thm.




    Our Contribution.

Our new algorithms close this gap by providing a viable trade-off between latency and communication volume.
We introduce a multi-level approach to Bingmann et al.'s algorithms, i.e., we use k levels where processor groups work on independent sorting problems.
The resulting multi-level mergesort has internal work and communication volume close to N for each level (<ref>).
This is significantly improved with prefix-doubling mergesort 
where internal work and communication volume per level are close to D (<ref>).
As side results of independent interest, we present an improved hypercube quicksort for strings 
(<ref>) and a multi-level distributed single-shot Bloom filter (<ref>).
While the idea to combine the best string sorting and atomic sorting algorithm is simple in principle, we view the analysis of a quite complicated overall algorithm in a realistic model of distributed-memory computing as a significant contribution. In particular,
because this guides an efficient, highly scalable implementation.
In an experimental evaluation, on up to 49152 cores, the multi-level algorithms are up to 5× faster than the single-level ones.

theoremTheoremTheorems



§ PRELIMINARIES






    Machine Model and Communication Primitives.

We assume a distributed-memory machine model consisting of p processing elements (PEs) allowing single-ported point-to-point communication.
The cost of exchanging a message of h bits between any two PEs is α + β h, where α accounts for the message start-up overhead and β quantifies the time to exchange one bit.
Let h be the maximum number of bits a PE sends or receives, then collective operations broadcast, prefix sum, (all-)reduce, and (all-)gather can be implemented in time 𝒪(αlog p + β h) <cit.>.
For personalized all-to-all communication, there is a trade-off between communication volume and start-up latency.
When data is delivered directly, we obtain a time complexity in 𝒪(α p + β h), while using a maximum degree of indirection yields 𝒪(αlog p + β hlog p ).
In our analysis, we resort to a more abstract view introduced by Axtmann et al.—a black box data exchange function Phr yielding the time complexity of exchanging data
, when each PE sends or receives at most h bits in total from at most r PEs <cit.>.
We find α r + β h as a lower bound for the time complexity of Phr, and although matching upper bounds are not known, there are indications suggesting that one can come close to this (see <cit.> for a brief discussion).
We use Phr = (1+o(1))Phr to sum multiple exchanges by the dominant one.






    String Properties and Input Format.

The input to our algorithms is a string array S = [s_0, s_1, …, s_n-1] consisting of n = S unique strings.
A string s of length ℓ is a sequence of characters from an alphabet Σ with s = [s[0], …, s[ℓ-2],] where ∉Σ is a sentinel character;  denotes the length of the longest string in S.
By N = S, we denote the total number of characters in S.
The ℓ-prefix of a string s are the first ℓ characters of s.
The longest common prefix (LCP) of two strings s ≠ t is the prefix of s with length st = s[i] ≠ t[i].
For a sorted string array, the corresponding LCP array [, h_1, h_2, …, h_n-1] contains the LCP values h_i = s_i-1, s_i.
For sorting the string array, we do not necessarily have to look at all characters in S.
The distinguishing prefix of a string s (with length s) are the characters that need to be inspected to rank s in S.
The sum of the lengths of all distinguishing prefixes of S is denoted by D.
By , we denote the length of the longest distinguishing prefix.
A string array is usually represented as an array of pointers referring to the start of the corresponding character sequence. This allows for moving or swapping strings in constant time.
The concatenation of all the character sequences forms the character array 𝒞(S) with 𝒞(S) = 𝒮.

In our distributed setting, we assume that each PE i obtains a local subarray S_i of S as input such that S is the concatenation of all local string arrays S_i.
Furthermore, we assume the input to be well-balanced, i.e., S_i = Θ(n/p) and S_i = Θ(N/p).






    Algorithmic Building Blocks.


We make use of an r-way LCP loser tree to merge r sorted sequences of in total m strings augmented with LCP values in 𝒪(mlog r +D) time <cit.>.
While merging, we update the LCP values of the resulting merged string sequence.
Furthermore, we use LCP compression, i.e., we send the longest common prefix of two consecutive strings (in a sorted sequence) only once.
While being very useful for many inputs in practice, LCP compression cannot substantially reduce the communication volume in the worst case <cit.>.

Robust Hypercube Quicksort (RQuick) is a sorting algorithm initially proposed for atomic sorting <cit.> that has been adapted to handle strings <cit.>.
RQuick only works for a number of PEs that is a power of two.
Hence, with d = log p, PEs with index i ≥ 2^d send their data to a PE within the d-dimensional hypercube.
Then the input data is randomly permuted among the PEs to make imbalances less likely.
Subsequently, each PE sorts its local data.
Then the actual d hypercube (quick)sorting rounds are executed, where the input is recursively partitioned into sub-hypercubes using a pivot element.
Hence, the elements are communicated at least d = log p times.



	
	If all input strings are unique, RQuick runs in time
	𝒪( n/plog n + αlog^2p + β(n/plog p + log^2p   )logσ) with probability ≥ 1 - p^-c for any constant c>0.


RQuick can be easily improved by using a string sorting algorithm for local sorting and (LCP-aware) merging.
Our implementation  uses these optimizations.
Therefore, the overall work performed by all PEs for these steps is 𝒪(nlog n + D).







§ MULTI-LEVEL STRING SORTING


Our multi-level merge sort (MS) approach adapts the ideas of Axtmann et al. for multi-level atomic sorting <cit.> to string sorting.
We recursively split PEs into groups each of which solves an independent (string) sorting problem.
There are k levels of recursion with arbitrary splitting factors between levels.
To obtain the single-level variant <cit.> of the algorithm choose k=1.
To simplify the analysis, we assume that p can be perfectly subdivided into r groups of p/r consecutive PEs, i.e., on the first level the j group consists of PEs jp/r, …, (j+1)p/r-1.
Furthermore, we generally assume approximately equal factors on each level, i.e., r=Θ√(p) which implies p=Θr^k, and that the final levels splits the PEs into groups of size 1.
The algorithm consists of a one-time initialization and a recursive phase which is invoked k times.
<ref> provides an overview.

Initialization:
The input is sorted locally.
On each PE i, sort the local input array 𝒮_i.
The LCP array can be obtained as a by-product of sorting.

Recursion:
A global order is established recursively.
Initially, all string arrays 𝒮_i must be locally sorted.
On each level of recursion, we have p' PEs and r groups of size p”=p'/r.


	
  * 
	      Distributed Partitioning:
	      Globally determine r-1 splitter strings f_j and on each PE i compute local buckets ℬ_i^0, …,
          ℬ_i^r-1 with ℬ_i^j = { s ∈ S_i | f_j < s ≤ f_j+1} for j ∈0, …, r-1 using sentinels f_0 = -∞ and f_r = ∞.





	
  * 
	      String Assignment and Exchange:
	      On PE i, the strings
	      in bucket ℬ_i^j are assigned to PEs belonging to group j. By ℬ^j = ⋃_i ℬ_i^j we denote the union of all strings assigned to group j.
	      Then, all strings and LCP values are exchanged using direct messaging.

	
  * 
	      Local LCP-aware Merging:
	      On PE i, the received string sequences are merged to obtain locally sorted string arrays
	      𝒪_i (using the LCP values).
	      We also update the LCP values for 𝒪_i during merging.
	      We then set p' ← p'/r and S_i ← O_i in the subsequent recursive step.





 §.§ Distributed Partitioning

Due to the multidimensionality of the string sorting problem, determining balanced partitions is far more challenging than in atomic sorting.
Some of the steps of our merge sort algorithm depend on the number of strings in the local string array while others depend on the number of characters or the size of the distinguishing prefix.
We therefore adapt string-based and character-based partitioning <cit.> to our multi-level approach.
These schemes bound the number of strings and characters, respectively.
The general approach is to draw and globally sort a number of samples on each PE.
Then, r-1 splitters f_j are chosen from the global sample array.
Since all strings are locally sorted before the partitioning step, we can make use of a regular sampling approach <cit.> in which the samples are drawn equidistantly. 



  §.§.§ String-Based Partitioning

On the  recursion level, there are r^-1 groups of PEs working on independent sorting problems, see <ref>.
We now describe the partitioning process from the point of view of one such group.
Let S' be the concatenation of the local string arrays of the PEs of one such group of size p'.
Unlike for single-level MS, now, we cannot assume an equal number of strings on each PE as from the second level of recursion on, this number is subject to the result of a previous partitioning round which is not exact.
For multi-level MS, string-based partitioning consists of the following two steps:




	
  Local Sampling: Let v > 0 be the sampling factor.
	      In total, there will be p'(v+1) samples drawn from the local string arrays.
	      To simplify the discussion, we assume S' to be divisible by p'(v+1).
	      Let ω = S'/(p'(v+1)).
	      PE i then draws Si/ω - 1 samples 𝒱_i from its local string array spaced as evenly as possible.
	      𝒱 is the union of all local samples.
	      If 𝒱 < p'(v+1), then the first p'(v+1) - 𝒱 PEs draw one additional sample.
	      This ensures that at most ω strings are between two local samples on each PE.
	      Also, the global number of samples is a multiple of p' and of r as we find p' = r^k+1- on recursion level .
	
  Splitter Computation: The samples 𝒱 are globally sorted using hypercube quicksort.
	      Then r-1 splitters f_j = 𝒱[j|𝒱|/r -1] for 0 < j < r are determined using a prefix-sum.
	      Subsequently, the r-1 splitters are communicated to all PEs using an all-gather operation.


A sampling factor v = Θ(r) yields a maximum number of strings per bucket in Θ(𝒮'/r).
This is shown in detail by using a generalization of the sample density lemma from <cit.> in <ref>. 
The proof of the following <ref> can also be found in <ref>.


	
	On recursion level  with r = √(p) in step <ref> of multi-level MS, string-based regular sampling with sampling factor v yields a maximum bucket size of
	|B^j| ≤(1 + r/v)^n/r^.


The term (1+r/v)^ signifies that the imbalance between the buckets  multiplies with each level of recursion.
Therefore, we need to choose v=Θ(kr) for any number of k levels to keep the term asymptotically constant during the entire sorting process.
For single-level MS, v = Θ(r) = Θ(p) samples per PE are sufficient. 
Note the large difference between drawing kr samples per PE (on average) here and r^k = p samples in the single-level case.
Using the assignment strategy described in <ref>, which equally distributes the strings over the PEs in each group, we arrive at <ref>.


	
	Using a sampling factor of v = Θ(kr) = Θ(k√(p)) the number of strings per PE is in n/p in each level of the algorithm.




  §.§.§ Character-Based Sampling

We generalize character-based regular sampling <cit.> to our multi-level approach to achieve tighter bounds on the number of characters per PE than the conservative 𝒪( n/p).
Now, each PE of the considered group draws S_i/ ω' - 1 equally spaced samples from its character array with sampling distance ω'= S'/(p'(v+1)).
To arrive at the final string samples, we shift the sampled character positions by at most - 1 characters to the beginning of a string.
If the total number of samples is smaller than p'(v+1), the first PEs draw one additional sample.

In <ref>, we give bounds on the number of characters per bucket over the course of our algorithm when using character-based partitioning.
This is shown in detail by using a generalization of the (character-) sample density lemma <cit.> in <ref>. 

	
	On recursion level  with r=√(p) in
	step <ref> of multi-level MS using character-based regular
	sampling with a sampling factor of v, each bucket contains at most ℬ^j≤*1+r/v^*N/r^	+ *1+v+1/rp/r^-1 characters.

The additional term depending on  stems from shifting the sampled positions to the beginning of strings.
By distributing the characters equally over the PEs in each group up to additional 𝒪() characters (see <ref>), we can limit the maximum number of characters per PE in <ref>.

	
	Using a sampling factor in Θ(kr) the maximum number of characters per PE in each level is in
	𝒪*N/p + k^2r.

For the single-level case with k=1, this is equivalent to O(N/p + p) which is the bound in the original algorithm <cit.>.
For k > 1, we even have an improvement over the single-level algorithm.
Since we assume k in log p/loglogp, we find k^2r = 𝒪(log^2(p)√(p)) = o(p).
This may seem counter-intuitive at first as we introduce a potential imbalance already in the first recursion level.
However, the subsequent assignment step distributes this imbalance equally over p'/r PEs.



 §.§ String Assignment and Exchange


We now assign the strings in bucket ℬ^j to the j (sub-)group consisting of PEs j,…, (j + 1)p” - 1 with p” = p'/r = r^k - on the  level.
The resulting assignment needs to ensure that each PE receives approximately the same amount of data.
Additionally, the number of sent and received messages per PE should be bounded by the number of groups r.
We generalize an approach proposed by Axtmann et al. <cit.> for distributed atomic sorting to string sorting.
For the sake of simplicity, we assume that each PE contributes the same number of strings (characters) in the assignment process.
With (slightly) imbalanced data (due to the partitioning) the below-stated results hold up to a small factor.

For string-based partitioning, we want to balance the number of strings per PE.
Since the assignment algorithm does not rely on internal properties of the elements, we can treat a string as an atomic object and apply Axtmann et al.'s assignment algorithm directly.

For character-based partitioning, we want to achieve a balanced number of characters but cannot split strings.
Therefore, we reiterate the steps of the algorithm and describe the necessary adjustments.
A local bucket ℬ_i^j is small if it contains at most ℬ^j/(2rp”) characters.
Small buckets are separately enumerated for each group j using a prefix sum where each PE contributes its number of small buckets for group j.
The t small bucket belonging to group j is then assigned to PE *t/r of group j.
This way, each PE gets assigned no more than half of its final capacity and receives messages from at most r different PEs.

Then, a description of each large bucket located on PE i and destined for group j is first sent to PE *i/r in group j and each group computes a balanced assignment independently of each other.
Conceptually, this works by performing separate prefix sums over residual capacities (remaining after assigning small buckets) and sizes of unassigned buckets, in each group.
The resulting sorted sequences of integers R (residual capacities) and U (unassigned buckets) must then be merged such that the bucket beginning at the i element is preceded by the PE containing the i open slot.
A subsequence of ⟨ r_i,u_j,…, u_j+h, r_i+1,z ⟩ in the merged sequence of R and U means that the local buckets u_j, …, u_j+h are assigned to PE i.
The last bucket u_j+h potentially needs to be split up (respecting string boundaries in the character-based assignment) and partly assigned to r_i+1 or even r_i+2 if z = r_i+2.
Since we cannot split up strings, we may end up with a PE obtaining up to - 1 additional characters.
Strings from one local bucket cannot be assigned on more than 3 PEs as the residual capacity on each PE is at least ℬ^j/(2p”).
Since large buckets contain more than ℬ^j/(2rp”) elements, a single PE can store at most 2r of them.
Hence, each PE receives r messages.
We refer to <cit.> for details on the actual group-local merging process.

	
	Using the bounded assignment algorithm <cit.>, we obtain a message assignment where each PE sends and receives r messages
	and each PE in group j obtains ℬ^j/p” strings (string-based) or at most ℬ^j/p” + characters (character-based).


Hence, even for character-based sampling, we can find an assignment, such that the number of characters per PE remains in the bounds of <ref>.
Afterwards, each PE sends its strings according to their assignment.
The time for computing the assignment is dominated by the actual data exchange <cit.>.
[
	For string-based assignment, this naturally transfers from the atomic sorting case.
	For character-based assignment, we have to compensate a communication volume in 𝒪(rlogN).
	Assuming unique strings, we find ≥ = Ω(logn/logσ) and N ≤ n and therefore log(N) = 𝒪(logσ).
	Since our bound for the number of (received) characters per PE contains an imbalance of at least r characters (which require an encoding of logσ bits each), the data exchange dominates the assignment also for character-assignment.]



 §.§ Overall Running Time



Let Θ(kr) be the sampling rate of our algorithm with character-based sampling.
By <ref>, the maximum number of characters per PE at any time is 𝒪(/p) with  = N + k^2r p.
Since n ≤ N/ and N ≤ n, the number of strings per PE then is
𝒪(/p) with = /(n + k^2rp).

We now combine the running time of the three phases of our algorithm in each level, including 𝒪(n/plogn/p + D/p) time for the initial local string sorting <cit.>.

1. Distributed Partitioning. Here, we have to globally sort 𝒪(rk) local sample strings of length ≤ per PE.
With RQuick, this is possible in 𝒪(rklogp(1 + βlogσ) + αlog^2p) expected time.
Allgathering the r-1 splitters needs 𝒪(αlogp + β rlogσ) time.
As we have r = 𝒪(/p) by definition and the local string array as well as the splitters are sorted, computing the local buckets ℬ_j^i is possible in time 𝒪(/p).

2. Assignment and Exchange. The string assignment is dominated by the data exchange which is possible in time p^i/k𝒪(/plogσ)𝒪(r) on level .
This holds as a PE stores 𝒪(/p) characters (<ref>) encoded in logσ bits each.
By the bounded assignment algorithm each PE exchanges strings with at most 𝒪(r) other PEs (<ref>).

3. Local (LCP-aware) Merging. The r sorted sequences of strings received in the data exchange now needs to be sorted to restore our invariant on the local string array.
While being beneficial for many inputs in practice, asymptotically, LCP-aware merging does not yield substantial advantage.
Thus, we resort to uninformed merging for our analysis.
Processing the 𝒪(r) sorted sequences with 𝒪(/p) characters in total is possible in 𝒪(/p + /p logr) time, dominating the time required for computing the local buckets during partitioning.


	
  Multi-level MS with r=√(p), = N + k^2r p, and = /(n + k^2rp), using character-based sampling and
	bounded assignment, runs in expected time
	
    T_MS(n,N,) = 𝒪(n/plogn/p + D/p^local sorting + k/p + /plogp^merging + k*αlog^2 p+
    					k√(p)log p(1 + βlogσ)^partitioning)
    
    			+ k·p^l/k𝒪*/plogσ𝒪*√(p)_assignment + exchange.


As we only draw k√(p) local samples and receive √(p)-1 splitter strings for k>1 as opposed to p samples and p-1 splitters in the single-level case
we no longer have to compensate the mediocre scalability of the single-level partition phase with a huge amount of data and are capable of sorting small to medium sized inputs on a large number of PEs.
Assuming the exchange primitive Ph√(p) to run in 𝒪(α√(p) + β h) <cit.>, we achieve a latency in 𝒪(α k√(p))= o(α p) at the cost of an k times higher communication volume.
If we additionally assume ≤ N/(k^2√(p)plogp) and k ≤logp/(2loglogp), we can state a simplified running time of multi-level MS in <ref>.

  
  With the above assumptions, we obtain a running time of multi-level MS in 𝒪(N/plogn + α k√(p) + β kN/plogσ) in expectation.




§ MULTI-LEVEL PREFIX DOUBLING MERGE SORT



The distinguishing prefix of S is usually much smaller than total number of characters N.
In a distributed algorithm, we can use this property to reduce the communication volume by only exchanging the distinguishing prefixes.
By doing so, instead of explicitly sorting the input strings we obtain the information on where to find the i smallest string of the input.
This, however, is sufficient in many use cases where string sorting is used, e.g., for suffix sorting <cit.>.

Bingmann et al. approximate the distinguishing prefix of each string by an upper bound in an iterative doubling process <cit.> using a distributed single-shot Bloom filter (dSBF) <cit.>.
In each round they hash prefixes with geometrically increasing length of the strings and globally check for uniqueness of the hash values.
If the hash value of a prefix with length d of string s is unique, we find s≤ d and s no longer needs to participate in the process.
That way for each string s an approximate distinguishing prefix with length s≥s can be determined for each string in expected 𝒪(log()) rounds assuming a constant false positive probability of the Bloom filter.
This can be achieved with expected latency in 𝒪(α p log) and expected bottleneck communication volume in 𝒪(n/plog p) + o(D/plogσ) <cit.>.[The latency can be reduced by increasing the communication volume by a factor Θ(logp) <cit.>.]
By employing a k-level Bloom filter for duplicate detection, we generalize this approach to arbitrary levels of indirection.

	
	Using communication on a k-dimensional
	grid, performing at most  operations (insertions, queries) per PE
	on a dSBF of size
	m ≥ en can be done in time
	𝒪(k(α
      p^1/k+βlogmp/n + logk))
	in expectation and with probability ≥ 1 - 1/p^ω(1)
	assuming the total number of operations n = ω(k^2p^1+1/klogp) and additionally m = (n).

A proof of <ref> can be found in <ref>.
A problem of using Bloom filters in the prefix doubling process is that the precondition on the overall number of operations required in <ref> might not hold when more and more strings drop out of the process because their distinguishing prefixes have already been determined.
We therefore switch to duplicate detection using (atomic) hypercube quicksort for sorting the hash values once there are too few strings left. 
Combining these approaches results in <ref>.

  
  For each string s ∈𝒮 with s≥logp/logσ an approximation s with s = 𝒪(s) can be computed in time
    
    𝒪( α k√(p)log^latency + β k (n/plogp+ D/plogσ)^communication  volume  + kn/plogkloglogσ + D/p^local  work)

	in expectation.
	We assume a balanced distribution of strings and their distinguishing prefixes, i.e., Θ(n/p) strings and Θ(D/p) per PE, and an overall number of strings n = 𝒪((p)).
    Additionally, we assume n/p = ω(k^2√(p)logploglogp) and k ≤logp/(2loglogp).

A proof for <ref> can be found in <ref>.
Since we assume all strings to be unique, we can also bound the average distinguishing prefix length D/n = Ω(logn / logσ) = Ω(logp / logσ).
Hence, there are only few (≤ n/σ) strings with small distinguishing prefixes for which <ref> does not yield an (expected) constant factor approximation of the actual distinguishing prefix.
Therefore, we find the expected value of the sum D_≈ of the approximate distinguishing prefixes determined with <ref> to be in 𝒪(D).
In conjunction with the assumed balanced distribution of D over the PEs, we find 𝒪(/p) with = D + k^2r p as an upper bound on the expected number of characters per PE in the multi-level merge sort algorithm when executed on the approximated distinguishing prefixes only.
Furthermore, since = Ω(logn/logσ) = Ω(logp/logσ), and we therefore approximate  up to an expected constant factor, we find T_MS=(n, D, ), as an upper bound for the running time of the merge sort part of multi-level prefix doubling merge sort (PDMS) in <ref>.

	

	Multi-level PDMS with r=√(p), using character-based sampling and
	assignment, and assuming the preconditions of <ref>, runs in expected time
    
    𝒪*log*α k√(p)^add. latency + kβn/plogp^add. comm. volume + kn/plogkloglogσ^add. local work_prefix doubling + T_MS(n,D,).


Note that the 𝒪(β kD/plogσ + D/p) part of the running time of the prefix doubling process in <ref> is subsumed by T_MS(n,D,) and thus not explicitly stated in <ref>.




§ EXPERIMENTAL EVALUATION


We now discuss the experimental evaluation of the following distributed-memory algorithms.

	
  _k Our new multi-level string merge sort with k levels of recursion, see <ref>.
	
  _k Our new multi-level doubling string merge sort including prefix approximation using grid-wise Bloom filter with k levels of recursion, see <ref>.
	
    <cit.> with our string-specific optimization, see <ref>.

The single-level variants _1, _1 and  are implementations by Bingmann et al. <cit.> that we improved slightly.
We also include the state-of-the-art shared-memory parallel algorithm  <cit.>.
We ran  on an AMD Epyc Rome 7702P CPU with 64 cores
(2 base and 3.5 boost frequency) equipped with 1024 DDR4 ECC RAM.
All distributed-memory experiments were performed on SuperMUC-NG consisting of 6336 nodes (792 nodes per island).
Each node is equipped with two Intel Skylake Xeon Platinum 8174 CPUs (24 cores each, 3.1 base frequency) and 96 RAM.
Communication between nodes uses a [per-mode=single-symbol]100 Omni-Path network.
All algorithms are implemented in C++ and compiled with GCC 11.2.0 with flags  and .
We use the KaMPIng <cit.> MPI bindings for our implementations and Open MPI v4.0.7 for interprocess communication.
Reported times are the average of five runs excluding the first iteration (MPI warm-up phase).





All variants use string-based regular sampling with a sampling factor of 2.
We use  to sort the samples.
For simplicity, we use a grid-wise group (string) assignment in the implementation of our multi-level algorithms.
Here, the PEs are arranged in a p'× r grid and PEs exchange buckets only along the rows, i.e., PE i sends its bucket for subgroup j along its row to the PE in the j column.
Also, for all-to-all exchanges a simple k-dimensional grid all-to-all is used which provides a latency in α√(p).
LCP compression is used during string exchange phases for data sets where significant common prefixes can be expected.

Multi-level variants always ensure one group per node on the final level
of sorting, i.e., k-level variants fall back to k-1 or k-2 levels for p<2^k-148.
For three-level variants, group sizes for the first two levels are chosen such that splitting
factors are as close as possible.

For strong-scaling experiments, we use real-world data sets, see <ref> for details.

	
  CommonCrawl (CCrawl)
	      consists of the first 100 of WET files from the Sep./Oct. 2023 Common Crawl archive
	      (<https://index.commoncrawl.org/CC-MAIN-2023-40/>).
	      The WET format consist mostly of plain text with small meta data headers.
	
  Wikipedia (Wiki)
	      consists of a dump, from 2023-12-20, of all pages in the English Wikipedia
	      in XML format without edit history
	      (<https://dumps.wikimedia.org/enwiki/20231220/>).

	
  WikipediaText (WikiText)
	      are from dumps of Wikipedia without any XML metadata.







For our weak-scaling experiments we use the DNGenerator string generator <cit.> to generate string sets with configurable D/N ratio (denoted by DNData).
This allows us to influence running times
of local sorting, the effectiveness of LCP compression, and the length of distinguishing prefixes.

Additional results are given in <ref> for the following experiments.



 §.§ Sorting Small to Medium Sized String Sets

In this weak-scaling experiment, we evaluate the algorithms on DNData with a string length of 500 characters, D/N = 0.5, and n/p ∈10^4, 10^5, 10^6.
See <ref> for the results.
The largest input is near the realistic limit for this system with roughly 2 GB RAM per PE.
We could not run the  variants on it, due to their memory consumption.


The running times broadly confirm the expected relation between input size and scaling behavior of algorithms.
Two-level merge sort significantly outperforms the single-level version on all input sizes for
sufficiently large values of p.
For small inputs, adding a third level leads to further improvements from 256 cores on.
As expected, the improvement is most obvious for the smallest inputs with n/p=10^4.
Here, the single-level algorithms scale roughly linearly with the number of PEs, as the running times
approximately double for every doubling of p.
For _1 the scaling behavior can mostly be attributed to the time required for partitioning,
with Θ(p) samples on each PE needing to be sorted.
For _1 a significant amount of time is also spent in the distinguishing prefix approximation.
The  variants perform significantly worse than our (multi-level) merge sort algorithms.




 §.§ Influence of D/N Ratios


Our first experiment already shows that multi-level merge sort exhibits improved scaling
properties. 
Now, we have a closer look on the influence of the D/N ratio on the running times of the different algorithms.
As before, we use strings with a length of 500 characters and 10^5 strings per PE.
<ref> shows the average running times.
We evaluate  and  on one and two levels and the  variants.
Three-level variants are not part of this experiment as the additional level only yields clearly better performance for ≥ 512 compute nodes which we did not include in any further experiments due to computing budget constraints.

The influence of D/N ratio is clearly visible for _k and _k, as variants with prefix approximation are superior on instances with ratio up to 0.5.
For ratios 0.75 and 1, the prefix approximation encompasses the whole string and is detrimental to the running time.

LCP compression is highly effective due to the nature of DNData inputs.
As before, multi-level variants outperform single-level counterparts, usually with a crossover
point at 32 nodes.
The gap seemingly increases for larger D/N ratios, e.g. for _k on 128 nodes the
speedups are 1.56 and 1.78 for ratios 0 and 1 respectively.
This can partially be explained if communication volume is considered.
On fewer PEs, two-level variants cause roughly twice the communication volume because string
exchange phases dominate.
At D/N=0, sending 0.5;1 per string is roughly equivalent to exchanging
every string once or twice respectively.
Communication increases for single-level variants with the number of PEs as partitioning requires
more samples to be sorted.
String exchanges send fewer characters for larger D/N ratios due to LCP compression, sample
sorting remains roughly constant.

The  variants perform worse than the merge sort based algorithms.
However, we can also see that for an increasing D/N ratio improving the algorithms by using LCP-aware merging pays off with  being up to 20 % faster for D/N = 1 than plain .

The running times of (shared-memory using 64 cores) are given for the input of all p processors.
We could not run  on a node of SuperMUC-NG as these are only equipped with 96 GB RAM.
With D/N ≥ 0.25, we need 384 cores to match the performance of  and only 192 cores for D/N ≥ 0.5.
Data with D/N ≈ 0 is difficult for distributed algorithms since the time spent in local sorting hardly compensates for the for communication overheads and start-up latencies.
Here, we only need 1536 cores to be on par with .
This is of interest when string sorting is part of more complex distributed tasks as it indicates that sorting the strings directly on the distributed system is faster than transferring them to a sufficiently large shared-memory machine from a very modest number of cores on.



 §.§ Evaluation on Real-World Data

To evaluate our algorithms on real-world data we use a strong-scaling experiments.
On all three data sets we can observe that the two-level algorithms scale better than their single-level counterparts and finally outperform them.
On the largest PE configuration utilizing 256 compute nodes, _2 ist the fastest algorithm on all three data sets.
This highlights the usefulness of our multi-level approaches—even on very skewed real-world inputs.
However, the differences of the running times are less pronounced.
Also, _2 fails on CommonCrawl and _2 fails on Wikipedia when using only 4 compute nodes, due to imbalances of the data.
This shows that string inputs can be inherently hard to partition.

While different sampling and group assignment techniques may improve character balance, the
CommonCrawl data set contains many duplicated strings (e.g., standard
legal dis­claimers) which are always assigned to a single PE and cannot be shortened by prefix
approximation, resulting in imbalances that we cannot prevent.







§ CONCLUSION AND FUTURE WORK

We demonstrate—in theory and practice—that string sorting can be scaled to very large number of processors.
Our best algorithm, a multi-level prefix-doubling merge sort, only requires internal work and communication volume close to the optimum (the total length of all distinguishing prefixes) per level.
In practice, all our multi-level algorithms outperform their single-level counterparts on a wide range of inputs on up to 49152 cores (from a modest number of cores on).
This is especially of importance in scenarios where string sorting is part of a distributed application, i.e., it is not feasible to sort the data on a large shared-memory machine because of the transfer costs.
Hence, we see our work as an important building block to enable more complex string-processing tasks  at a massively parallel scale.

One problem we want to tackle in the future is suffix sorting based on sorting strings of equal length <cit.>.
To this end, we plan to extend our algorithms to support space-efficient string sorting.
Some inputs are highly compressible with lots of overlapping strings, e.g., the suffixes of a text with length n with a combined length of approximately n^2/2 can be represented using only n characters.
However, due to the scarcity of main memory on most supercomputers, we cannot easily materialize all strings at the same time during sorting.












§ DISTRIBUTED PARTITIONING

Due to the multidimensionality of the string sorting problem, determining balanced partitions is far more challenging than in atomic sorting.
Some of the steps of our merge sort algorithm depend on the number of strings in the local string array, while others depend on the number of characters or the size of the distinguishing prefix.
We therefore adapt string-based and character-based partitioning <cit.> to our multi-level approach.
These schemes bound the number of strings and characters, respectively.
The general approach is to draw and globally sort a number of samples on each PE.
Then, r-1 splitters f_j are chosen from the global sample array.
Since all strings are locally sorted before the partitioning step, we can make use of a regular sampling approach <cit.> in which the samples are drawn equidistantly.


 §.§ String Based Partitioning


On the  recursion level, there are r^-1 groups of PEs working on independent sorting problems, see <ref>.
We now describe the partitioning process from the point of view of one such group.
Let S' be the concatenation of the local string arrays of the PEs of one such group of size p'.
Unlike for single-level MS, now, we cannot assume an equal number of strings on each PE as from the second level of recursion on, this number is subject to the result of a previous partitioning round which is not exact.
For multi-level MS, string-based partitioning consists of the following two steps:



	
  Local Sampling: Let v > 0 be the sampling factor.
	      In total, there will be p'(v+1) samples drawn from the local string arrays.
	      To simplify the discussion, we assume S' to be divisible by p'(v+1).
          This can be lifted when allowing for a small additional imbalance factor (1 + 1/k) (see <ref>).
	      Let ω = S'/(p'(v+1)).
	      PE i then draws Si/ω - 1 samples 𝒱_i from its local string array spaced as evenly as possible.
	      𝒱 is the union of all local samples.
	      If 𝒱 < p'(v+1), then the first p'(v+1) - 𝒱 PEs draw one additional sample.
	      This ensures that at most ω strings are between two local samples on each PE.
	      Also, the global number of samples is a multiple of p' and of r as we find p' = r^k+1- on recursion level .

	
  Splitter Computation: The samples 𝒱 are globally sorted using hypercube quicksort.
	      Then r-1 splitters f_j = 𝒱[j|𝒱|/r -1] for 0 < j < r are determined using a prefix-sum.
	      Subsequently, the r-1 splitters are communicated to all PEs using an all-gather operation.


Although we do not necessarily draw the same number of samples from each PE, we can apply the adapted sample-density lemma:

	
	For i = 0, …, (p'-1) let S”_i = {s ∈S_i' | a ≤ s ≤ b } be an arbitrary contiguous subarray of S_i'.
	If |Si”∩Vi| = k, then |Si”| ≤ (k+1)ω with ω = |S'| / (r(v+1)).

This gives us bounds for the number of strings per bucket.

	
	Using string-based partitioning with a sampling factor of v, every bucket ^j contains at most (1+r/v ) S'/r strings.


    The proof transfers from the one of a similar claim <cit.>.
	By counting the number of sample strings that are contained within a bucket B^j and applying <ref>, we obtain the claimed bounds.
    Let 𝒱_i^j=ℬ_i^j∩𝒱_i be the samples on PE i that fall into the
	j bucket.
	Using <ref> it follows that
	ℬ_i^j≤|𝒱_i^j|+1ω.
	By definition in step <ref> of the partition algorithm, splitters f_j and f_j+1 are separated by
	at most p'(v+1)/r-1 elements and thus ∑_i=0^p'-1𝒱_i^j = p'(v+1)/r when including
	f_j+1.
	Using the definition of ℬ^j we find 

	
    ℬ^j    = ∑_i=0^p'-1ℬ_i^j≤∑_i=0^p'-1*𝒱_i^j+1ω
       ≤ω*p'(v+1)/r + p'
    		= S'/p'(v+1)*p'(v+1)/r + p'
        = S'/r + S'/v+1≤S'/r + S'/v =  (1+r/v ) S'/r

	

	


To summarize, using string-based regular sampling yields buckets with a maximum imbalance of n/v strings.
Choosing v = Θ(r) implies a number of strings per bucket in Θ(𝒮'/r).

	On recursion level  with r = √(p) in step 2) of multi-level MS, string-based regular sampling with sampling factor v yields a maximum bucket size of
	
    |B^j| ≤(1 + r/v)^n/r^.



	We give a proof by induction. For the first level of recursion, <ref> directly yields a maximum bucket size of
	
    |B^j| ≤(1 + r/v) n/r.

	Applying the bounds for recursion level -1 ⇝ to the same theorem  results in
	
    B^j≤(1 + r/v) (1 + r/v)^-1n/r^-1/r  = (1 + r/v)^n/r^.



The term (1+r/v)^ signifies that the imbalance between the buckets  multiplies with each level of recursion.
Therefore, we need to choose v=Θ(kr) for any number of k levels to keep the term asymptotically constant during the entire sorting process.
For single-level MS, v = Θ(r) = Θ(p) samples per PE are sufficient. 
Note the large difference between drawing kr samples per PE (on average) here and r^k = p samples in the single-level case.
Using the assignment strategy described in <ref>, which equally distributes the strings over the PEs in each group, we arrive at <ref>.


	
	Using a sampling factor of v = Θ(kr) = Θ(k√(p)) the number of strings per PE is in n/p in each level of the algorithm.




 §.§ Partitioning Without Loss of Generality



Let ω = 𝒮'/(p'(v+1)) and ω is not an integer, i.e., ⌈ω⌉ = 1 + ⌊ω⌋.
We then use p'(v+1) (potentially overlapping) samples of size ⌈ω⌉.


	Using string-based partitioning with a sampling factor of v and S' = Ω(p'(v+1)k), every bucket ^j contains at most (1+r/v ) S'/r(1 + 1/k) strings.


	Let 𝒱_i^j=ℬ_i^j∩𝒱_i be the samples on PE i that fall into the
	j bucket.
	Using <ref> it follows that
	ℬ_i^j≤|𝒱_i^j|+1⌈ω⌉.
	By definition in step <ref> of the partition algorithm, splitters f_j and f_j+1 are separated by
	at most p'(v+1)/r-1 elements and thus ∑_i=0^p'-1𝒱_i^j = p'(v+1)/r when including
	f_j+1.
	Using the definition of ℬ^j we find 

	
    ℬ^j    = ∑_i=0^p'-1ℬ_i^j≤∑_i=0^p'-1*𝒱_i^j+1⌈ω⌉
       ≤⌈ω⌉*p'(v+1)/r + p' = (⌊ω⌋ + 1) *p'(v+1)/r + p'
       ≤ω*p'(v+1)/r + p' +*p'(v+1)/r + p'
       ≤ω*p'(v+1)/r + p' +ω/k*p'(v+1)/r + p'
        = ω*p'(v+1)/r + p'*1 + 1/k
        = S'/p'(v+1)*p'(v+1)/r + p'*1 + 1/k
        = (S'/r + S'/v+1)*1 + 1/k≤*S'/r + S'/v*1 + 1/k
        =  (1+r/v )*1 + 1/kS'/r

	

	



 §.§ Character Based Sampling


We generalize character-based regular sampling <cit.> to our multi-level approach to achieve tighter bounds on the number of characters per PE than the conservative 𝒪( n/p).
Now, each PE of the considered group draws S'_i/ ω' - 1 equally spaced samples from its character array with sampling distance ω'= S'/(p'(v+1)).
To arrive at the final string samples, we shift the sampled character positions by at most - 1 characters to beginning of the string.
If the total number of samples is smaller than p'(v+1), the first PEs draw one additional sample.

Bingmann et al. proposed to shift the sampled character positions by at most -1 positions to the beginning of the next string.
Furthermore, they required ℓ < ω'.
We will show that this restriction is not necessary and therefore slightly adapt their proof of <cit.>.
Note that our approach might result in sampling the same string multiple times.
In this case, the string must be unified by concatenating its (local) index which requires at most 𝒪(log(n/p)/logσ) additional characters.
As we assume all strings to be unique, this is in 𝒪().



Again, we state a slightly adapted version of the sample density lemma for character-based sampling:

	
	Let 𝒮_i”=s∈𝒮_i' | a≤ s≤ b be an arbitrary contiguous
	subsequence of 𝒮_i' for i∈0, …, p'-1.
	With 𝒮_i'∩𝒱_i=k it must hold that
	𝒮_i”≤(k+1)ω'+.


  Initially, there have been at most ω' - 1 characters between two sampled positions.
  Since each sampled character position has been moved at most  characters, there are at most ω' -1 + characters between two shifted sampled strings.
  If k=0, then all elements of 𝒮_i” are either contained between two consecutive samples of 𝒱_i or smaller than the first sample or greater than the last.  Thus, 𝒮_i”≤ω' - 1 + (or 𝒮_i”≤ω' + in the latter case).
  For k=1, let c be the first character of the sample string contained in 𝒮_i”.
  We can split the character array 𝒞(𝒮_i”) of 𝒮_i” into two character arrays 𝒞_< +  c + 𝒞_>.
  For each of the two character arrays, k=0 applies and we thus find 𝒮_i”≤ (ω' -1 + ) + 1 + (ω' + ) ≤ 2(ω' + ).
  For k ≥ 2, we can split the corresponding character array of 𝒮_i” into (k+1) subarrays of which the first contains at most ω' -1 + characters,
  the last one at most ω' + characters and all in between at most ω' - 1 + characters.
  Adding the first character of the k sample strings then yields 𝒮”_i≤ (k+1)(ω' + ).


We now arrive at the following bounds for character-based partitioning into r buckets.


	
	Using character-based regular sampling with a sampling factor of v, every bucket ℬ^j
	obtained using the partition algorithm contains at most *1+r/v𝒮'/r+ *1 + v+1/rp' characters.


	As for its string-based counterpart, by counting the number of sample strings that are contained within a bucket B^j and applying <ref>, we obtain the claimed bounds on the maximum number of characters per bucket.
    Let 𝒱_i^j=ℬ_i^j∩𝒱_i be the local sample buckets as in the proof
	of Theorem <ref>.
	Using Lemma <ref> yields the analogous bound
	ℬ_i^j≤*𝒱_i^j+1ω'+.
	The equalities for 𝒱_i^j and ℬ^j remain identical which suffices to obtain
	the stated bound.
	
    ℬ^j    = ∑_i=0^p'-1ℬ_i^j≤∑_i=0^p'-1*𝒱_i^j+1ω'+
       ≤*p'(v+1)/r+p'*𝒮'/p'(v+1)+
       ≤𝒮'/r+𝒮'/v+*p'(r+v+1)/r
        = *1+r/v𝒮'/r+ *1 + v+1/rp'

	

Note that the term 1 + v+1/rp' in
Theorem <ref> can be simplified to
Θ(p') if v=Θ(r).
Hence, the bound on the number of characters in a bucket ℬ^j does not only depend on the number of samples but also on the number of PEs as a (character-)shift of up to - 1 might occur on every PE.
In <ref>, we give bounds on the number of characters per bucket over the course of our algorithm when using character-based partitioning.

	
	On recursion level  with r=√(p) in
	step <ref> of multi-level MS using character-based regular
	sampling with a sampling factor of v, each bucket contains at most ℬ^j≤*1+r/v^*N/r^	+ *1+v+1/rp/r^-1 characters.


	We give a proof by induction.
	The base case of =1 can be derived immediately from
	Theorem <ref> as follows:
	
    ℬ^j≤*1 + r/vN/r+ *1+v+1/rp≤*1 + r/v*N/r+ *1+v+1/rp

    as r/v≥ 0 and therefore (1 + r/v) ≥ 1 holds.
	We now proceed with the inductive case -1 by applying the same theorem again on each
	group.
	Note that the value of p' needs to be changed to reflect the current group size, i.e., 
	p/r^-1.
	With that we find:
	
    ℬ^j≤*1+r/vInduction  Hypothesis  for  ( - 1)/r + *1+v+1/rp/r^-1
       = *1+r/v*1+r/v^-1*N/r^-1 + (-1)*1+v+1/rp/r^-2/r + *1+v+1/rp/r^-1
        = *1+r/v^*N/r^ + (-1) *1+v+1/rp/r^-1  +  *1+v+1/rp/r^-1
       ≤*1+r/v^*N/r^ + (-1) *1+v+1/rp/r^-1 +  *1+r/v^_≥ 1*1+v+1/rp/r^-1
        = *1+r/v^*N/r^ + *1+v+1/rp/r^-1



The additional term depending on  stems from the shifting to the beginning the of the strings,
By distributing the characters equally over the PEs in each group up to additional 𝒪() characters (see <ref>), we can limit the maximum number of characters per PE in <ref>.

	
	Using a sampling factor in Θ(kr) the maximum number of character per PE is in
	
    𝒪*N/p + k^2r.



For the single-level case with k=1, this is equivalent to O(N/p + p) which is the bound in the original algorithm <cit.>.
For k > 1, we even have an improvement over the single-level algorithm.
Since we assume k in log p/loglogp, we find k^2r = 𝒪(log^2(p)√(p)) = o(p).
This may seem counter-intuitive at first as we introduce a potential imbalance already in the first recursion level.
However, the subsequent assignment step distributes this imbalance equally over p'/r PEs.










§ DISTRIBUTED DUPLICATE DETECTION


A key building block of the computation of the approximate distinguishing prefixes approach is distributed duplicate detection.
Bingmann et al. <cit.> use a distributed single-shot Bloom filter <cit.> for this task, which allows for approximate membership queries.
Conceptually, a distributed single-shot Bloom filter of size m is a bit array of size m which is equally distributed over the p PEs.
To insert an element e, the element is hashed to a random position h(e) within the interval [0,m) and the corresponding bit is set on the PE responsible for this position, where h denotes a random hash function.
For querying whether an element e is contained in the filter, one simply has to check whether the corresponding bit at position h(e) is set.
This may result in a false positive result.
A false positive query probability f^+ for a single-shot Bloom filter with n inserted elements can be achieved by having a size m ≥ nf^+ <cit.>.
In our setting, we assume that operations are executed in batches, i.e., each PE i has n_i insertion or query operations.
The hash values associated with the operations are compressed before sending, for example using Elias-Fano encoding <cit.>.
It can be shown that the expected communication volume for each insertion/query operation is in logp + 𝒪(1) bits for Bloom filters with constant false positive probability.
For more details, we refer to <cit.>.


A problem of distributed (single-level) single-shot Bloom filters is that they have a latency (at least) linear in p since messages are delivered directly [Note that to mitigate this, Bingmann et al. propose hypercube all-to-all communication, which reduces the latency to 𝒪(αlogp) at the cost of a communication volume increased by a logarithmic factor.].
As another restriction, we find that the expected communication volume in 𝒪(logp) per operation is only guaranteed if the overall number of queries is large.

To mitigate both problems, we propose to apply our k-level communication scheme also to Bloom filters, which reduces their latency to 𝒪(α k√(p)) at the cost of a communication volume of 𝒪(klogp) bits per operation.
Furthermore, in <ref> we will show that k-level Bloom filters only require ω(k^2p^1+1/klogp) operations for an expected communication volume in 𝒪(klogp).

However, since the number of strings participating in the prefix doubling process decreases from iteration to iteration, we might eventually end up with less queries from some iteration on.
When this is the case, we will use (atomic) hypercube quicksort for duplicate detection.
An upper bound for the running time for this approach will be given in <ref>.


 §.§ Multi-Level Bloom filters


Before discussing the running time of multi-level distributed single-shot Bloom filters in <ref>, we will prove <ref>, which analyzes a balls-into-bins with duplicates scenario.
<ref> will then be used as a building block in the proof of <ref>.

Let M be a multiset consisting of m elements on a group of p' processors, where each PE holds = m/p' elements.
Furthermore, we assume the elements to be locally unique.
Let S denote the set of unique elements within M and for e ∈ S let c_e denote the number of occurrences of e within M.
Each local element is now assigned uniformly at random to one of the p' PEs.

	
	Let X denote the maximum number of elements a PE receives in the above-described process.
	For k ≥ 1 and = ω(k^2p'logp) with p ≥ p', we find
	
    Pr[X > (1+1/k)] ≤1/p^ω(1).


(Outline)
	In our proof outline, we mainly reiterate the idea behind the proof of <cit.>.
	For each (unique) element e ∈ S, let 1_e,j be an indicator variable with 1_e,j = 1 iff element e is sent to PE j.
	As the elements are assigned uniformly at random, we have Pr[1_e,j = 1] = 1 /p'.
	Let X_j denote the number of elements received by PE j.
	We find X_j = ∑_e∈ Sc_e1_e,j.
    Since ∑_e ∈ S c_e = m, we find X_j = m/p' for the expected number of elements received by PE j by the linearity of expectation.
	With the same arguments as in <cit.> one can show that X_j is least sharply concentrated around its mean when all elements have a multiplicity c_e = p',
	i.e., each PE has the same set of elements.
	It is therefore sufficient to show a bound on Pr[Y_j > /p'(1+1/k)] for this worst-case setting in which we have X_j = p'Y_j, where Y_j denotes the number of unique elements received by PE j.

	As Y_j ∼Binom(,1/p'), we can use a Chernoff-bound <cit.> to obtain
	
    Pr[Y_j > /p'*1 + 1/k] ≤exp*-/3p'k^2

	as an upper bound on the probability that PE j receives more than /p'(1+1/k) (unique) elements for k ≥ 1.
	By the union bound argument, the probability that any of the p' PE receives more than /p(1+1/k) (unique) elements is at most
	
    p'exp(-/3p'k^2) ≤ p'exp(-ω(logp'))  = p'/p^ω(1) = 1/p^ω(1)

	for = ω(p'k^2logp). Since X_j = p'Y_j for all PE j, we overall find Pr[X > (1 + 1/k)] ≤1/p^ω(1).


We now use <ref> as a building block to show a bound on the running time for a multi-level Bloom filter.


	
	Using communication on a k-dimensional
	grid, performing at most  operations (insertions, queries) per PE
	on a distributed single-shot Bloom Filter of size
	m ≥ en can be done in time
	𝒪(k(α
      p^1/k+βlogmp/n + logk))
	in expectation and with probability ≥ 1 - 1/p^ω(1)
	assuming the total number of operations n = ω(k^2p^1+1/klogp) and additionally m = (n).

(Outline)
	We focus on the part that is different from <cit.>.
	We route operations to their recipient PE determined by the hash value of the corresponding elements using a k-dimensional grid.
	This takes k rounds of communication.
	In each of these rounds, there are r = √(p) messages per PE incurring a latency overhead of 𝒪(α√(p)).
	In each round, we have the invariant that the received operations, i.e., hash values that are either queried or eventually inserted, are increasingly sorted.
	If a hash value occurs multiple times, it is forwarded only once.
	Due to the grid communication scheme, in each round  with 0 ≤ < k, PEs exchange messages only with PEs located within the same row for dimension  of the k-dimensional grid.
	These PEs constitute a group g.
    In the following, we will show that the probability that any PE in round  receives more than (1 + 1/k)^ hash values is in 1/p^ω(1).

	We first consider the number of operations received by a PE j in the first round and then show the claim by induction.
	Without loss of generality, we assume that each PE stores  operations.
	While we can assume that hash values associated with the (insert/query) operations are locally unique, we have the situation that they can occur multiple (at most r) times across the PEs within a PE-group g.
	By <ref>, we know that each PE in g receives more than (1 + 1/k) operations with probability ≤ 1/p^ω(1).
	This holds since we assume ≥ n/p = ω(k^2rlogp) and therefore fulfill the condition on  in <ref> with p' = r.
	As there are at most 𝒪(p) such groups, by the union bound, the probability that a PE in any group receives more than (1+1/k) hash values is also ≤ 1/p^ω(1).
	Since duplicate hash values are only forwarded once, the assumption that each hash value occurs at most r times within a PE group also holds for the following dimension.

	Let us assume that in round 0 <  < k-1 a PE has received at most (1+1/k)^ hash values with probability ≤ 1/p^ω(1) and each hash values occurs at most r times
	in the PE group g for round +1.
	We apply <ref> with =  (1+1/k)^ and p' = r and therefore have =  (1+1/k)^≥≥ n/p = ω(k^2rlogp).
	Using the union bound argument as before, we find the event that any PE receives more than (1+1/k) =  (1+1/k)^(1+1/k) = (1+1/k)^+1 occurs with probability in 1/p^ω(1).
	Again as duplicate values are forwarded only once, we also fulfill the second inductive assumption.

	Thus, overall, we can conclude that any PE sends or receives at most (1+k)^≤(1+k)^k ≤ e hash values with probability ≥ 1 - 1/p^ω(1) in round .

	We now argue about the required communication volume in bits when a PE exchanges at most e hash values.

	With Elias-Fano encoding <cit.>, encoding x out of m one-bits can be achieved using at most mx = x(logm/x+2) bits.
    Therefore, a message containing x operations on a Bloom filter of size m can be encoded with at most mx bits. 
	In each round, each PE sends (and analogously receives) r messages containing n_1, …, n_r operations with ∑_j=1^r n_j = n' and we find n' ≤ e with high probability.
	These messages can be encoded using
	∑_j=1^r m/rn_j bits
	[Note that we can use m/r as the universe size within the encoding, as each PE only receives hash values from an interval of size m/r.
  This could even be generalized to m/r^+1 for round .].
	Due to the concavity of m/r·, we obtain the following bound on the communication volume
	
    ∑_j=1^r m/rn_j   ≤ r ·m/rn'/r≤ r ·m/re/r
        = e(logm/e + 2)                                                     
        = 𝒪((logmp/n + 1))

	with high probability as m/rx is monotonically increasing in its second component for x ≤ 4m/(er) and we have n/p ≤≤ n and m ≥ en.
	The claim about the expected running time follows directly for n being polynomial in p, for even larger n one has to refine the estimation in <ref>.

	The answers (yes or no) to the queries are sent on the reverse route using a single bit per request.
	Routing them to the PE that originally sent each request is possible by keeping appropriate tables mapping bit positions to sender PEs.
	These tables are stored locally and need no additional communication volume.

	The local work in each round is dominated by (integer-) sorting the received hash values.
    This can be achieved in 𝒪(logk) time with high probability <cit.> as we have at most 𝒪() values on a PE in any round with high probability and assume m = (n).
	The initial sorting of the at most  requires one additional sorting step.
    Thus, we end up with 𝒪( k logk) local work with high probability.



 §.§ Duplicate Detection with Atomic Hypercube Quicksort

In the following, we will describe how atomic hypercube quicksort <cit.> can be used for distributed duplicate detection.

	
	For n = 𝒪((p)) (atomic) hypercube quicksort <cit.> can be used to perform distributed duplicate detection with constant false positive probability f^+
	in time
	
    𝒪*logp + log^2p(α + β( + logp)
 in expectation and with probability ≥ 1-p^c for any constant c > 0
	where  is the maximum number of elements located on a single PE.


	As a first step, we map each element to an integer in [0,nf^+) using a random uniform hash function h.
	It is then sufficient to globally check whether there are duplicate hash values to achieve duplicate detection with false positive probability f^+ <cit.>.

	After hashing each local element e_j, we construct hash value/index pairs (h(e_j), j) which are then subsequently sorted using atomic hypercube quicksort <cit.>.
	Note that j denotes the global index which can be computed using a distributed prefix sum over the number of elements on each PE.
	Since all indices are unique, the hash value/index pairs are also unique.
	During the sorting process, we keep track of the path a pair takes on each level by storing this information into a local table.
	Once the pairs are (lexicographically) sorted, a scan over the local elements and two additional message exchanges for the first and last local element suffices to identify duplicates.
	To communicate the result back, we replace the hash value within each hash value/index pair with a single bit indicating whether the element is unique or not.
	Then the pairs are sent back in 𝒪(logp) iterations using the stored routing information.

	Since n is polynomial in p, we find lognf^+ = 𝒪(logp).
	Therefore, each hash value/index pair occupies only 𝒪(logp) bits and we achieve the above stated high probability bound.
	We need 𝒪(n^c') time per level in case all elements accumulate on a single PE during the course of the algorithm, which is highly unlikely as this event occurs with probability smaller than p^c for an arbitrary large constant c.
	Since we assume n to be polynomial in p, the above bound on the running time also holds in expectation.





§ DISTINGUISHING PREFIX APPROXIMATION VIA DOUBLING



In the following, we prove the running time of our distinguishing prefix approximation algorithm using k-level Bloom filters and (atomic) hypercube quicksort for duplicate detection.


  
	For each string s ∈𝒮 with s≥logp/logσ an approximation s with s = 𝒪(s) can be computed in time
    
    𝒪( α k√(p)log^latency + β k (n/plogp+ D/plogσ)^communication  volume  + kn/plogkloglogσ + D/p^local  work)

	in expectation.
	We assume a balanced distribution of strings and their distinguishing prefixes, i.e., Θ(n/p) strings and Θ(D/p) per PE, and an overall number of strings n = 𝒪((p)).
    Additionally, we assume n/p = ω(k^2√(p)logploglogp) and k ≤logp/(2loglogp).

(Outline)
    We first briefly discuss the overall algorithmic idea to obtain the above stated running time and then discuss details.
    We start with an initial prefix length l^init = Θ(logp/logσ) for the prefix doubling process using a k-level Bloom filter (see <ref>) with a constant false positive probability f^+.
	By doubling the tested prefix length l in each round, this results in 𝒪(log) rounds of duplicate detection in expectation and therefore an expected overall latency in 𝒪(αlogk√(p)).
	We therefore no longer have to deal with the latency term when arguing about the expected running time.
    A limitation of k-level Bloom filters is that they require a certain number of operations n^thr so that a running time in 𝒪(k(βlogp + logk)) per operation can be achieved (see <ref>).
    However, during the course of the prefix doubling process, more and more strings drop out as an approximation of their distinguishing prefix has been found.
    We mitigate this problem by switching to (atomic) hypercube quicksort once the number of strings becomes too small.

    𝐥 ≤ 𝐥^𝐭𝐡𝐫. We now start with a discussion of the running time until we have reached a tested prefix l = l^thr with l^thr = log^3(p).
	Let 𝒳_l denote the maximum number of strings participating in the duplicate detection round with prefix length l on any PE.
    As stated in <ref> (k-level Bloom filter), with probability 1-1/p^ω(1) we have a running time for each round of duplicate detection in 𝒪(𝒳_l · k(βlogp + logk)).
    As an upper bound for the running time in the unlikely case where all operations are routed to the same PE occurring with probability in 1/p^ω(1), we find 𝒪(X_l· pk(βlogp+ logk)).
    Therefore, by the linearity of expectation, the expected running time is in 𝒪(𝒳_l k(βlogp + logk)) in each round.
    In total, we therefore have ∑_i=1^cloglogp𝒪(𝒳_2^i l_initk(logpβ + logk)) as an upper bound on the expected running time for rounds l^init < l ≤ l^thr with a constant c > 0 as it takes 𝒪(loglogp) rounds to reach l = l^thr.

    In the first round, we have an expected running time in 𝒪(n/pk(βlogp + logk)) as all strings participate.
    It remains to bound the size of the expected value of 𝔼(X_l) for the remaining rounds which can be done using the average distinguishing prefix length = D/n.
    We find = (D/p)/(n/p) = D/n = Ω(logn/logσ) = Ω(logp/ logσ) on each PE as we assume all strings to be unique.
    A general observation is that there are at most n/(p2^j) strings with distinguishing prefix length greater than 2^j for any j ≥ 0.
    Assuming a constant false positive probability, for l ≥, we have ∑_≤ l𝒳_l = 𝒪(n/p).
    Therefore, the expected running time for those rounds is dominated by the first one.
    This also holds when the number of participating strings drops below holds n^thr = ω(k^2p^1/klogp) during the first 𝒪(loglogp) rounds, since we have 𝒪(n^thrk(βlogp + logk)) as an upper bound on the expected running time in these cases.
    Since we assume n/p = ω(k^2p^1/klogploglogp) = Ω(n^thrloglogp), the accumulated running time for those 𝒪(loglogp) rounds is bounded by the running time of the first round

    We now differentiate three cases:
    
      
  * ≤ l^init. In this case, the above argument directly applies and the accumulated running time for the first 𝒪(loglogp) rounds is dominated by the running time of the first round.

      
  *  
        l^init < ≤logp. 
        In this case, there are 𝒪(log(logσ/ logp)) = 𝒪(loglogσ) rounds until we reach l =.
        Assuming that all n/p strings on a PE participate in these rounds  –  note that this is an upper bound on 𝒳_l  –  causes an expected communication volume in 
        
    𝒪(log(logσ/logp)n/pklogp) = 𝒪(logσ/logpn/p klogp) = 𝒪(k D/plogσ).

        For expected local work, we have to charge 𝒪(loglogσ· n/pklogk).
        The subsequent rounds with l > are then dominated by the running time of the very first round.

        
  *  
        logp < . For the rounds until we reach l = logp, we have the same running time as in case <ref>.
        For the remaining 𝒪(log(/logp) rounds until we reach l = min(, l^thr), we may also assume that all strings further participate as an upper bound on 𝒳_l.
        We then have a running time in

    𝒪(log(/logp)n/pk(logpβ + logk)) = 𝒪(/logpn/pk(logpβ + logk)) = 𝒪(D/p(kβ + 1)).

        Note that klogk = 𝒪(logp) as we assume k ≤logp/loglogp.
        Potential further rounds with  < l ≤ l^thr are dominated by the running time of the first round.
    
    Therefore, we find a running time for the first 𝒪(loglogp) iterations as stated in the theorem above.
    

    𝐥 > 𝐥^𝐭𝐡𝐫.
    Until now we have analyzed the expected running time of the algorithm until a tested prefix length l = l^thr =  Ω(log^3p) is reached.
    From now on we switch to duplicate detection with hypercube quicksort once X_l drops below n^thr (see <ref>).
    Since the running time of duplicate detection via k-level Bloom filters is dominated by that of duplicate detection via hypercube quicksort (apart from the latency term which has already been dealt with) we can assume without loss of generality that only hypercube quicksort is used from now on.
    As stated in <ref>, with probability 1-1/p^c_1 for an arbitrary large constant c_1 we have a running time for each round of duplicate detection in 𝒪(𝒳_llog^2p(β + 1) + 1_𝒳_l > 0βlog^3p) with 1_𝒳_l > 0 being the indicator variable whether there are still strings in the duplicate detection process.
    As an upper bound for the running time in the unlikely case occurring with probability less than 1/p^c_1, we find 𝒪(n^c_2).
    This is in 𝒪(p^c_3) as n is polynomial in p with c_2,c_3 > 0 being fixed constants.
    By the linearity of expectation and by choosing c_1 > c_3, we find the expected running time to be in 𝒪(𝒳_llog^2p(β + 1) + 𝒪(1_𝒳_l > 0log^3pβ) in each round.

    We begin by finding a bound for the latter term.
    If the longest distinguishing prefix ≤log^3p, we know that the expected number of rounds with l > l^thr is constant.
    Thus, 𝒪(1_𝒳_l > 0log^3pβ) = 𝒪(log^3pβ) = 𝒪(n/plogpkβ) as n/p = ω(logp√(p)) and k ≤logp/loglogp.
    If the longest distinguishing prefix > log^3p, we have / log^3p > 1.
    Hence, we have log / log^3p = 𝒪( / log^3p) remaining rounds in expectation and the accumulated costs of 𝒪(log^3pβ) in each of these rounds are in 𝒪(log^3pβ / log^3p) = 𝒪(D/p), as ≤𝒪(D/p).

    Now we discuss the 𝒪(𝒳_llog^2p(β + 1)) part of the running time.
    We differentiate two cases:
    
      
  * < 𝐥𝐨𝐠^2𝐩. The number of strings on any PE with distinguishing prefix length log^3p≥logp or longer is smaller than n/(plogp) as discussed above
      Therefore, ∑_l=l^thr^∞𝒳_l = 𝒪(n/(plogp)).
      Hence the accumulated costs for duplicate detection with hypercube quicksort are in 𝒪(n/(plogp)log^2p(β + 1) = 𝒪(n/p logp(β + 1) and thus are dominated by the running time of the very first iteration of the prefix doubling process.
      
  * ≥𝐥𝐨𝐠^2𝐩. We have D/p = Ω(n/plog^2p). For the 𝒪(log( / log^2p)) rounds until l reaches , we can therefore simply assume that n/p strings participate in the duplicate detection causing expected costs in

    𝒪(log(/log^2p)n/plog^2p(β + 1)) = 𝒪(/log^2pn/plog^2p(β + 1)) = 𝒪(D/p(β + 1)).

        The costs of the remaining rounds when l ≥ are again in 𝒪(n/plog^2p(β +1)) = 𝒪(D/p(β + 1)).
        In fact, if we recall that in these rounds we actually use the k-level Bloom filter until n^thr is reached and switch to hypercube quicksort once we fall below this threshold,
        we even obtain a running time in o(D/p(β + 1).
  

  Therefore, in every case the summed expected running time is in the bound stated in the theorem above.








§ MORE EXPERIMENTAL RESULTS


In this section, we give additional experimental results that highlight some of the key properties of our new distributed multi-level string sorting algorithms.
These results break the average running times down and give time requirements for different phases of the algorithms.

First, in <ref>, we give the average sorting times per sorting phase for the weak scaling experiments on synthetic data with different n/p rations, i.e., the results presented in <ref>.
Here, we can see that our multi-level approach significantly reduces the time spent in the partitioning phase (both  and ) and the approximation of the distinguishing prefix, i.e., the Bloom filter ( only).
For the total running times, we also give more values for slower algorithms, which is feasible due to the size of the plot.
This further highlights the scalability of our new algorithms.
When not only considering  but also , we can achieve a speedup of up to 7.

Next, in <ref>, we depict the communication volume and the average sorting times per sorting phase for the weak scaling experiments on synthetic data with different D/N rations, i.e., the results presented in <ref>.
Communication volume is measured on each PE during execution using the size of buffers passed to
MPI routines and summed afterward to arrive at the final value.
Different rules are applied depending on the used routine:
For example, calls to  and , as well as reduce and
scan operations count the size of send buffers on all PEs.
Calls to  only count the send buffer on the root PE and multiply its size by
p.
The result is only an approximation of the actual communication volume for several.
Most importantly, the measurements are idealized and do not necessarily correspond to actual
communication performed by MPI – especially for broadcast and reduction operations.
Communication of collective exchange operations may also be overestimated if send buffers include
data that is already on the correct PE and does not require sending.

Finally, in <ref>, we give the average running times per sorting phase for our strong scaling experiment using real-world inputs, which we depict in <ref>.
The results highlight the difficulty of the real-world inputs as many benefits that are visible for synthetic data are less prominent.




	

	


