Histogram Sort with Sampling
Vipul Harsh

Laxmikant Kale

Edgar Solomonik

University of Illinois at Urbana
Champaign
vharsh2@illinois.edu

University of Illinois at Urbana
Champaign
kale@illinois.edu

University of Illinois at Urbana
Champaign
solomon2@illinois.edu

arXiv:1803.01237v1 [cs.DC] 3 Mar 2018

ABSTRACT
To minimize data movement, state-of-the-art parallel sorting algorithms use sampling and histogramming techniques to partition
keys prior to redistribution. Samples enable partitioning to be done
using representative subset of the keys, while histogramming enables evaluation and iterative improvement of a given partitioning.
We introduce Histogram sort with sampling (HSS), which combines
sampling and histogramming techniques to find high-quality partitions with minimal data movement and high practical performance.
Compared to the best known algorithm for finding this partitioning,
our algorithm requires a factor of Θ(log(p)/log log(p)) less communication than the best known (recently introduced) alternative, and
substantially less when compared to standard variants of Sample
sort and Histogram sort. We provide a distributed-memory implementation of the proposed algorithm and compare its performance
to two existing implementations, and provide a brief application
study showing the benefit of the new algorithm.

KEYWORDS
parallel sorting; data partitioning; histogramming; sampling; histogram sort; sample sort;

1

INTRODUCTION

Finding a global partitioning of the data is the key challenge separating parallel sorting from sequential sorting. Partition-based sorting
algorithms, which (in contrast to merge-based sorting algorithms)
find a splitting of the data prior to redistributing it, are advantageous on modern architectures due to their low communication
cost. Sampling data either uniformly or selectively and histogramming the split produced by a partition are the two most common
techniques for determining a good partitioning. By quantifying the
parallel execution cost in terms of computation and communication,
we demonstrate that a simple combination of these two techniques
leads to an algorithm that provides both theoretical and practical
improvements over the best previously known results.
A parallel sorting algorithm needs to redistribute N keys across
p processors such that they are in a globally sorted order. In such
an order, keys on processor k are no greater than keys on processor
k + 1 and keys are sorted within each processor. An exact splitting
is achieved if all processors own the same number of keys, while
an approximate splitting guarantees that every processor owns no
more than N (1 + ϵ)/p keys for some ϵ. Given keys sorted with
Permission to make digital or hard copies of part or all of this work for personal or
classroom use is granted without fee provided that copies are not made or distributed
for profit or commercial advantage and that copies bear this notice and the full citation
on the first page. Copyrights for third-party components of this work must be honored.
For all other uses, contact the owner/author(s).
,,
© 2018 Copyright held by the owner/author(s).

an approximate splitting with ϵ = O(1), an exact splitting can be
achieved at a fairly small cost. However, a large ϵ increases the memory footprint, so algorithms that are able to find a partitioning that
is balanced for a given ϵ provide favorable robustness guarantees.
The most natural way to cheaply determine a global partitioning
is to collect a sample of keys, and infer a global partitioning from
the ideal partitioning of the sorted sample. Sample sort [11] and its
variants are basic realizations of this partitioning approach, which
are widely used in practice [21], and also serve as building blocks
for HSS. Selecting a random sample of the data and partitioning
in this way suffices to achieve achieve the desired load balance
with high probability so long as Θ(p log N /ϵ 2 ) are collected in the
sample [15]. A deterministic balanced splitting is also possible via
sampling, for example, using sample sort with regular sampling [20,
24]. With regular sampling, the algorithm collects p/ϵ keys from
each processor, that partition the input data evenly, requiring a total
sample size of Θ(p 2 /ϵ) for a balanced split. However, these classical
results leave substantial room for improvement, as we show, it
suffices to collect a number of samples that scales near-linearly
with p and logarithmically with 1/ϵ.
The technique of histogramming is embodied in the Histogram
sort algorithm [18, 25], which iteratively refines a partitioning (set
of splitters), by repeatedly collecting histograms of the number
of keys in each interval induced by the latest set of splitters. The
number of histogramming rounds required to determine all splitters
within the allowed threshold is bounded by O(log Z), where Z is
the size of the input domain. For skewed distributions, the number
of rounds could be large and the use of the input domain implies that
Histogram sort is not a pure comparison-based sorting algorithm.
Recently, Axtmann et al. [3] proposed a partitioning scheme
based on a histogram of the partitioning induced by a random
sample. They show that using a sample of size Θ(p log(p)/ϵ) results in an ϵ-balanced partitioning with high probability. Our main
contribution is demonstrating that by using log(log(p)/ϵ) steps of
refinement with histogramming, Θ(p log(log(p)/ϵ)) samples in total
suffice for an ϵ-balanced partitioning. This reduction improves the
communication cost for the partitioning stage proportionally, at
the cost of a small increase in the number of parallel steps (BSP
supersteps / synchronizations). This factor of improvement also
holds if the partitioning schemes are used in a multi-stage fashion,
√
for example by first splitting the data into p parts, then sorting
√
each part recursively with p processors.
The improvement in cost warrants the introduction of an algorithm that combines sampling and histogramming, which we
call Histogram sort with sampling (HSS). While it employs very
standard techniques, the analysis of the algorithm is nontrivial. The
main challenge resolved in the paper is in identifying and proving an invariant that shows global quadratic convergence of the
partitioning algorithm. The main intuition behind the algorithm

,,

Vipul Harsh, Laxmikant Kale, and Edgar Solomonik

and proof comes from consideration of splitter intervals, which are
subranges around ideal splitter keys in the globally sorted order
of input keys. In each round, HSS uniformly samples keys in the
union of all splitter intervals, then tightens each splitter interval
using a histogram of the collected sample. Our main analytical
result is that the size of the union of the splitter intervals decreases
geometrically with the number of rounds.
By characterizing HSS, we establish the theoretical soundness
of histogramming as a technique, that has been known to be effective [3, 17, 25, 26] in practice. Our algorithm is simple, provably
robust to arbitrary distributions with skewed or repeated keys, and
effective in practical scenarios where the keys are already partially
sorted globally. We provide and evaluate a parallel Charm++ implementation of the HSS algorithm, demonstrating improvements
in scalability over one of the fastest publicly-available distributedmemory parallel sorting algorithms, HykSort [26]. Additionally,
we show that our algorithm improves performance with respect to
Histogram sort within the ChaNGa N -body simulation code [16],
which uses parallel sorting at every time-step to distribute moving
cosmological bodies along a space-filling curve. Our theoretical
analysis of parallel execution cost, comparative evaluation of scalability, and application case study unanimously identify HSS as the
preferred parallel sorting algorithm.

2

PROBLEM STATEMENT

Let A(0), . . . , A(N − 1) be an input sequence distributed across p
processing elements, such that each processor has N /p keys. We
assume that there are no duplicates in the input. In Section 6.3 we
discuss how to reduce a sorting problem with duplicate keys to a
sorting problem with no duplicates, with very little overhead. Our
proofs and algorithm also translate to scenarios where input keys
are not evenly distributed across each processor. Parallel sorting
corresponds to redistributing and reordering the elements so that
the ith processor owns the ith subsequence of keys in the sequence
I (0), . . . , I (N − 1), where {I (0), . . . , I (N − 1)} = {A(0), ..., A(N −
1)} and I (i) ≤ I (i + 1). We say that key A(j) = I (r ) has rank r .
In practice, keys are typically associated with values, but in the
context of the algorithms we study, handling values of a given size
is straightforward.
It is common to additionally require that the resulting distribution of sorted keys is load balanced among processors. We compare algorithms with the standard assumption that the distribution is locally balanced, i.e. each processor owns no more than
(1 + ϵ)N /p keys. However, our algorithm achieves a stronger guarantee, namely that the distribution is globally balanced, i.e. processor
i owns all keys greater or equal to S(i) and less than S(i + 1), where
each S(i) is a splitter that satisfies S(i) = I (χ (i)), with
χ (i) ∈ Ti , where target range Ti = [Ni/p−N ϵ/2p, Ni/p+N ϵ/2p].
One practical advantage of a globally balanced distribution in the
context of iterative applications, is that if the initial distribution
is nearly sorted and globally balanced, the data exchange step is
guaranteed to require little data movement.
However, it is important to note that given either type of loadbalanced splitting, post-processing may be done to obtain an exact
splitting [8]. For a locally balanced distribution, this might require
some processors to potentially communicate all of their data to one

or two other processors. However, given a globally balanced distribution, achieving an exact splitting would require communicating
only at most N ϵ/p keys per processor. Therefore, a more fundamental distinction is in whether a parallel sorting algorithm maintains
load balance at all times, i.e. no processor is ever assigned more
than (1 + ϵ)N /p keys. Satisfying this condition permits bounded
memory footprint, which is desirable for a parallel sorting library
implementation.
The focus of our paper will be on the data-partitioning step
of partitioning-based sorting algorithms. Sample sort by regular
sampling [20, 24], histogram sort [18, 25], sample sort by random
sampling [7, 11], parallel sorting by over partitioning [19], AMS
sort [3], HyKSort [26] fall into this category. Partitioning-based
sorting algorithms generally determine a set of splitters that achieve
either a locally or globally balanced splitting, then redistribute keys.
The algorithm can run in multiple stages by splitting up data among
subsets of processors and sorting recursively within each subset.
We evaluate our algorithm using a standard simple variant of the
bulk synchronous parallel model (BSP) suggested by Valiant [28].
In particular, we model the parallel execution of an algorithm as a
sequence of supersteps, during each of which processors perform
computation locally, then send and receive messages. The BSP
complexity of an algorithm consists of three components:
• the number of supersteps (synchronization cost),
• the sum over all supersteps of the maximum amount of
computation done by any processor during the superstep
(computational cost),
• the sum over all supersteps of the maximum amount of
data sent or received by any processor during the superstep
(communication cost).
We permit processors to send and receive p messages every superstep, which simplifies the analysis of the cost of the all-to-all
exchanges involved in redistributing keys. The model suffices to
quantify the performance trade-offs in the data partitioning strategies we consider. In particular, more rounds of histogramming
increase the number of supersteps but lower communication cost.

3

RELATED WORK

Our algorithm borrows ideas from sample sort [11] and histogram
sort [18], we review these and other sorting algorithms before
proceeding to our main result.

3.1

Sample sort

Sample sort [6, 11, 14, 20, 24] is a standard well studied parallel
sorting algorithm. Sample sort samples s keys from each processor,
and sends them to a central processor to form an overall sample of
size M = ps keys. Let Λ = {λ 0 , λ 1 ..., λps−1 } denote the combined
sorted sample. Sample sorting algorithms choose p − 1 keys from
Λ as the final splitters. Generally, sample sort algorithms consist of
the following three-phase skeletal structure.
(1) Sampling Phase: Every processor samples s keys and sends it
to a central processor. s is often referred to as the oversampling
ratio. See Section 3.2 for different sampling methods.
(2) Splitter Determination: The central processor receives samples of size s (obtained in Step 1) from every processor resulting
in a combined sample Λ of size (ps). The central processor then

HSS

,,

selects splitter keys: S = {S(1), S(2)..., S(p − 1)} from Λ that
partitions the key range into p ranges, each range assigned
to one processor. Usually, the splitters are chosen by picking
evenly spaced keys from Λ. Once the splitters have been chosen,
they are broadcast to all processors.
(3) Data Exchange: Once a processor receives the splitter keys, it
sends each of its key to the appropriate destination processor. As
discussed earlier, a key in range [S(i), S(i + 1)) goes to processor
i. This step is akin to one round of all-to-all communication and
places all input data onto their assigned destination processors.
Once a processor receives all data that is assigned to it, it merges
them using a sequential algorithm, like merge sort.

3.2

Sample sort: Sampling methods

We discuss two sampling methods: random sampling and regular
sampling, for the sampling step (step 1) in sample sort.
3.2.1 Random sampling. With random sampling as described
by Blelloch et al. [7], each processor divides its local sorted input
into s blocks of size (N /ps) and samples a random key in each
block, where s is the oversampling ratio. The splitters are chosen
by picking evenly spaced keys from the overall sample of size ps,
collected from all processors. Of particular reference to our work is
the following theorem, (Lemma B.4 in [7]).
p log N 
Theorem 3.1. With an overall sample size of O ϵ 2 , sample
sort with random sampling achieves (1 + ϵ) load balance w.h.p..
3.2.2 Regular sampling. With regular sampling as proposed by
Shi et al. [20, 24], every processor deterministically picks s evenly
spaced keys from its local sorted data. The central processor collects
these samples and merges them to obtain a combined sorted sample
{λ 0 , ..., λps−1 } of size ps. Then, (p − 1) splitters are selected by
picking evenly spaced keys from this sample. We reproduce the
following theorem from [20, 24].
p

Theorem 3.2. If s = ϵ be the oversampling ratio, then sample
sort with regular sampling achieves (1 + ϵ) load balance.
Because of the large number of samples required, the sampling
phase is costly for regular sampling and makes it unscalable. Sample
sort with random sampling is more efficient than regular sampling,
but scalability is still hindered in practice because of the large
sample size it requires to achieve a load-balanced splitting.

3.3

Histogram Sort

Histogram sort [18, 25] addresses load imbalance by determining
the splitters more accurately. Instead of determining all splitters
using one large sample, it maintains a set of candidate splitter
keys and performs multiple rounds of histogramming, refining
the candidates in every round. Computing histogram of a set of
candidate keys gives the global rank of each candidate key. This
information is used by the algorithm to finalize splitters or to refine
the candidate splitter keys. Once all the splitters are within the given
threshold, it finalizes the splitter keys from the set of candidate
keys. The data exchange phase of Histogram sort is identical to
the third phase of sample sort. We give an overview of the splitter
determination step in histogram sort.

(1) The central processor broadcasts a probe consisting of M sorted
keys to all processors. Usually, the initial probe is spread out
evenly across the key range (unless additional distribution information is available).
(2) Every processor counts the number of keys in each range defined by the probe keys, thus, computing a local histogram of
size M.
(3) Local histograms are summed to obtain a global histogram at a
central processor using an M-item reduction.
(4) The central processor finalizes and broadcasts the splitters if
a probe key within the desired range has been found for each
of the unknown p − 1 splitters. Otherwise, it refines its probes
using the histogram obtained and broadcasts a new set of probes
for next round of histogramming, in which case the algorithm
loops back to Step 2.
Histogram sort is guaranteed to achieve any arbitrary specified
level of load balance. It is also scalable in practice for many input
distributions, since the size of the histogram every round is typically kept small - of the order O(p). The number of histogramming
rounds required to determine all splitters within the allowed threshold is at most log2 (Z/p), where Z is the range of the input (i.e.
maximum key minus the minimum key). The number of rounds
can be large, especially for skewed input distributions. Histogram
sort has been successfully deployed in real world, highly parallel
scientific applications, for instance ChaNGa [16].

3.4

Other Sorting Algorithms

Parallel sorting by over partitioning was proposed by Li et al. [19]
for shared memory multiprocessors. Every processor picks a random sample of size pks from its local input and sends it to a central
processor. The overall collected sample is sorted by the central processor and pk −1 splitters are chosen by selecting s th , 2s th , ..., (pk −
1)s th keys. These splitters partition the input into pk buckets, more
than required and hence k is referred to as the over partitioning
ratio. The splitters are made available to all processors and the local
input is partitioned into sublists based on the splitters. These sublists form a task queue and each processor picks one task at a time
and processes it by copying the data to the appropriate position
in the memory, determined using the splitters. The idea of over
partitioning is closely related to histogramming.
3.4.1 Merge based sorting algorithms. In this paper, we focus
on data partitioning parallel sorting algorithms, which we refer to
as partition-based algorithms. Merge-based algorithms are another
class of sorting algorithms that employ merging data using sorting
networks. An early such result was due to Batcher [4] which uses
time (or equivalently depth in a sorting network) O(log2 N ) with
N processors. Ajtai et. al [2] gave the first sorting circuit of depth
O(log N ), commonly referred to as the AKS network. The AKS
network has large hidden constants because of the use of expander
graphs in the circuit [9, 22]. Later, Cole [10] proposed a sorting
algorithm that runs in O(log N ) time using N processors and has
smaller constants than the AKS network. A communication optimal algorithm for the BSP model was proposed by Goodrich [13].
Cole’s merge sort and its adaptation to BSP by Goodrich follow a
merge-tree, but employ sampling to determine a partitioning that
accelerates merging. Overall, unless data-partitioning schemes are

,,

Vipul Harsh, Laxmikant Kale, and Edgar Solomonik

also employed, merge-based algorithms tend to be less performant
due to their need to do more data exchange stages. This overhead
often implies more BSP supersteps and in some cases more communication than partition-based alternatives.

3.5

Massive scale parallel sorting algorithms

Several recent works have focussed on large scale sorting. HykSort [26] is a state of the art practical parallel sorting algorithm.
It employs multi-staged splitting and communication to achieve
better scalability. HykSort is a hybrid of sample sort and hypercube
quick sort. Hyksort’s algorithm for splitter selection is similar to
HSS. We note that sampling and histogramming have been known
to be effective, our main contribution is providing sound theoretical
insights into techniques that have known to work well in practice.
We believe that our work serves as a justification for the efficacy of
HykSort for large problem sizes.
AMS-sort [3] employs overpartitioning for splitting. The scanning algorithm (Lemma 2 of [3]) is better than HSS with one round
of histogramming by a factor of Θ(min(log p, 1/ϵ)). However, HSS
with multiple rounds of histogramming is more efficient than AMSsort. The scanning algorithm used to select splitters in AMS-sort
does not easily generalize to multiple rounds of histogramming. Further, HSS achieves a globally balanced partitioning, while AMS-sort
achieves only a locally-balanced splitting, providing less robustness
in preservation of existing distributions. AMS-sort can be performed
in a multi-stage fashion, with successive steps of splitting and data
exchange across a decreasing set of processors. HSS can run in
the same multi-stage fashion, but with each data partitioning step
done with multiple rounds of histogramming. We compare the running times of AMS-sort and HSS with multiple stages in Table 3.
In our experiments, we focus mainly on single stage algorithms.
We emphasize that our performance evaluation is meant to be a
demonstration of practicality of using HSS for massive scale sorting rather than building a fully optimized sorting library, which
would warrant further considerations of low-level performance
optimizations and a separate study. Nevertheless, we compare our
implementation of single stage HSS to HykSort (Section 7.1).

3.6

Single stage AMS sort

The parallel sorting algorithm of Axtmann et al. [3] collects a single
set of samples, performs one round of histogramming, then picks
a locally balanced splitting based on the histogram. The splitters
obtained after the first histogramming round achieve the specified
level of load balance w.h.p. with an oversampling parameter that is
much less than sample sort with random sampling. In Section 5, we
show that the cost of histogramming is asymptotically same as the
cost of sampling an equal number of keys, so AMS sort achieves
a clear theoretical improvement over sample sort. We study the
AMS algorithm in some detail, due to its close relationship to our
approach. In particular, we formalize the proof of a key theorem
given by Axtmann et al. [3] (see Appendix A).
3.6.1 Scanning algorithm. AMS sort uses a scanning technique
to decide the splitters once the histogram is obtained. The algorithm
scans through the histogram and assigns a maximal number of
consecutive buckets (all keys between two consecutive keys in the
total sorted sample) to each processor. Specifically, after assigning

i buckets to the first j processors, it assigns buckets i + 1, . . . , i + k
to processor j + 1, where k ≥ 0 is picked maximally so that the
total load on processor j + 1 does not exceed N (1 + ϵ)/p. The last
processor gets all the remainder elements. If the sample size is
sufficiently large, the average load on the first p − 1 processors is
greater than N /p with high probability.
In particular, a sample of size Θ(p(log p + 1/ϵ)) is necessary
to achieve a locally balanced partitioning with high probability.
Demonstrating this formally is difficult due to the conditional dependence of loads assigned to consecutive processors. We extend
the proof sketch given by Axtmann et al. [3] in Appendix A. In
Table 2 we compare the cost of AMS sort to versions of sample
sort and HSS. AMS sort achieves a lower asymptotic complexity
than HSS with a single round of histogramming. However, HSS can
achieve an asymptotically lower complexity in O(1) BSP supersteps
and even lower complexity with O(log log p/ϵ) supersteps while at
the same time providing a globally balanced distribution.

4

HISTOGRAM SORT WITH SAMPLING

The basic skeleton of HSS is similar to that of Histogram Sort. In
addition, HSS employs sampling to determine the candidate probes
for histogramming. Every histogramming round is preceded by
a sampling phase where each processor samples local keys and
the overall sample collected from all processors is used for the histogramming round. By histogramming on the sample, HSS requires
significantly fewer samples compared to sample sort.

4.1

HSS with one histogram round

We first describe HSS with one round of histogramming, whose
data-partitioning step is slightly less efficient than AMS sort, by a
factor of Θ(min(log p, 1/ϵ)), but achieves a globally-balanced splitting. Extending the single round scanning algorithm of AMS sort to
multiple rounds for improved complexity is non-trivial. HSS with
one round is easily generalizable to multiple rounds of histogramming, which we discuss in subsequent sections and is the main
contribution of this paper. With multiple rounds of histogramming,
HSS is more efficient than the scanning algorithm, in fact, O(1)
rounds of histogramming suffice for an asymptotic improvement.
Recall that in HSS, a satisfactory ith splitter (in terms of global
load balance) is found when the algorithm finds a candidate key that
is known to have rank in target range Ti = [Ni/p − N ϵ/2p, N i/p +
N ϵ/2p]. If the sample contains at least one key with rank in Ti ,
then after histogramming on the sample, all such splitters will be
found. Intuitively, the algorithm should sample adequate number
of keys so that at least one key is picked from each Ti w.h.p..
Lemma 4.1. If every key is independently picked in the sample with
ps
2p ln p
probability, N = ϵ N , where s, the oversampling ratio is chosen
2 ln p

to be ϵ , then at least one key is chosen with index Ti for each i,
w.h.p..
Proof. Recall that the input set is denoted by A. The size of
Ti = N ϵ/p. The probability that no key is chosen with index in Ti
in the overall sample is given by,



 Nϵ
ps | Ti |
2p ln p p
1
1−
= 1−
⩽ e −2 ln p = 2
N
ϵN
p

HSS

,,

AMS sort scanning algorithm
sample sort: random sampling
HSS: 1 round
HSS: optimal rounds

256M
16M
1M

algorithm

Sample size (#keys)

problem

Figure 1: Figure illustrating HSS with multiple rounds. After first round, samples are picked only from the splitter intervals. Notice how the splitter intervals
shrink as the algorithm progresses.

64K
4K
256
16
1
16

64

256

1K
4K
#processors

16K

64K

256K

Figure 2: Sample size in practice: sample sort vs HSS vs AMS sort for 5%
load imbalance

Since there are p − 1 splitters, the probability that no key is chosen
from some Ti , is at most (p − 1) × p −2 < 1/p.
□

N
p
A(i)
I (r )
Ti

number of keys to sort in total
number of processors sorting keys
the i th input key
key with rank r in the overall global order
[N i/p−N ϵ /p, N i/p+N ϵ /p] is the target range for the i th splitter

sj

the sampling ratio for the j th round, in particular
each key in γ j is in the j th sample with probability s j p/N
rank of smallest sample key below rank N i/p after j rounds
rank of largest sample key below rank N i/p after j rounds
[I (L j (i)), I (U j (i))] is the i th splitter interval after j rounds
the union of all splitter intervals after j rounds

L j (i)
U j (i)
Ij (i) =
γj

Table 1: Notation used in paper, index j refers to Histogramming round,
while i is the processor index.

that the notion of sampling ratio is slightly different from the oversampling ratio of sample sort and one round HSS since the size
of the overall sample collected from all processors with the above
method is (ps |γ |/N ) in expectation.
4.2.2

The above Lemma leads us to the following theorem showing
global load balance of HSS with one round of histogramming. The
theorem will also prove to be useful in the analysis of multiple
rounds of histogramming, each round of which effectively increases
the oversampling ratio by collecting the same number of samples
from a smaller subset of the complete set of keys.
Theorem 4.2. With one round of histogramming and sample size
O(p log(p)/ϵ), HSS achieves (1+ϵ) load balance w.h.p..
Figure 2 compares the sample size required in practice for HSS,
sample sort with random sampling and the scanning algorithm of
AMS-sort.

4.2

HSS with multiple rounds

We show that HSS can be made more efficient by repeated rounds
of sampling followed by histogramming. We build upon the key
observation that after the first round of histogramming, samples
for subsequent histogramming rounds can be intelligently chosen
using results from previous rounds.
4.2.1 Sampling method. For the sampling phases, our algorithm
chooses a sample from a subset γ of the input. Initially, γ represents
the entire input. As the algorithm progresses, γ gets smaller. We
use the following method for sampling throughout our discussions,
unless stated otherwise.
Sampling Method: Every key in γ is independently chosen to
be a part of the sample with probability ps/N , where we refer s as
the sampling ratio.
Using the above sampling method simplifies the analysis, since
number of samples from disjoint intervals are independent. Note

Histogram sort with sampling (HSS) aglorithm.

(1) In the sampling phase before the first round of histogramming,
each key in the input is picked in the sample with probability
(ps 1 /N ), where s 1 is the sampling ratio for the first round. Samples from all processors are collected at a central processor and
broadcast as probes for the first round of histogramming.
(2) Every processor counts the number of keys in each range defined by the probe keys (the overall sample for the current
round), thus, computing a local histogram. All local histograms
are summed up using a global reduction and sent to the central
processor.
(3) For each splitter i, the central processor maintains L j (i): the
lower bound for the i th splitter rank after j histogramming
rounds, i.e. rank of largest key seen so far, which is ranked
less than Ni/p. Likewise it maintains U j (i), rank of smallest
key ranked greater than Ni/p. Once the histogram reduction
results of the j th round are received, the central processor updates L j (i) and U j (i) and broadcasts the sample keys I (L j (i)),
I (U j (i)) bounding each splitter by the splitter interval Ij (i) =
[I (L j (i)), I (R j (i))].
(4) Once every processor is aware of the new splitter intervals, it
begins its sampling phase for the (j + 1)th round. Every key
which falls in one of the splitter intervals is picked in the sample
with probability (ps j+1 /N ), where s j denotes the sampling ratio
for the j th round. If j < k, samples from all processors are
collected at a central processor and broadcast for the next round
of histogramming, in which case the algorithm loops back to
step 2. If j = k, the histogramming phase is complete and the
algorithm continues to step 5. Step 2 and 3 can be executed
efficiently if the local data is already sorted.

,,

Vipul Harsh, Laxmikant Kale, and Edgar Solomonik

(5) Once the histogramming phase finishes, the key ranked closest
to N i/p among the keys seen so far is set as the i th splitter.
In further discussions, we discuss how to choose k and the
sampling ratios s j ’s so that the splitters determined this way
result in a globally balanced partitioning.
A crucial observation is that the splitter intervals shrink as the
algorithm progresses and hence the sampling step is executed with
a subset of the input that gets smaller every round. Let γ j denote the
set of keys in the input that belong to one of the splitter intervals
after j rounds. |γ j | represents the size of the input that the algorithm
Í
samples from, for the j t h round. We have, |γ j | ⩽ i U j (i) − L j (i),
where U j (i) − L j (i) is the number of keys in the ith splitter interval.
Some splitter intervals can overlap, hence the inequality . In fact, it
is easy to see that there is no partial overlap between two splitter
intervals, that is, either two splitter intervals: Ij (i 1 ) and Ij (i 2 ) are
disjoint or they are identical.
We provide an outline of our proof. First we show in Lemma 4.3
that the algorithm will achieve a good splitting w.h.p. if the sampling
ratio is chosen to be large enough. With k rounds, the sampling
ratio in the last round should be set to be large enough so that the
algorithm terminates after k rounds.
Having shown that the algorithm terminates after k rounds,
achieving a globally load balanced partitioning, we bound the sample sizes in each round by first bounding |γ j | in terms of the sampling ratio s j necessary to obtain all splitters w.h.p.. Finally, we
appropriately set the sampling ratios such that the size of the union
of splitter intervals, that is, |γ j | decreases by a constant factor. Intuitively, sampling ratio s j in round j (where samples are chosen only
from the splitter intervals in round j) can be thought of as choosing
samples from the entire input range and discarding unncessary
samples.
2 ln p

Lemma 4.3. If sk = ϵ be the sampling ratio for the k th round,
then at least one key is chosen from each Ti after k rounds w.h.p..
2p ln p

Given a sampling ratio of s j = ϵ N , all splitters are found
w.h.p, by Lemma 4.1. We next bound the expectation of the size of
the union of all splitter intervals.
Lemma 4.4. Let s j be the sampling ratio for the j t h round, Ij (i) be
the splitter interval for i t h splitter after j rounds and γ j denote the
set of input keys that lie in one of the Ij ’s, then, E(|γ j |) ⩽ 2N
sj .
Proof. Since L j (i) and U j (i) are only improved every round,
L j−1 (i) ⩽ L j (i) ⩽

Ni
⩽ U j (i) ⩽ U j−1 (i)
p



Further ∀x : 0 ⩽ x ⩽ U j−1 (i) − Npi ,
h
i
ps j  x
Ni
P U j (i) −
⩾x = 1−
p
N

Consequently, we can bound the size of the ith splitter interval via
U j−1 (i)−
h
Õ
Ni i
E U j (i) −
=
p
x =1

Ni
p

h
i
Ni
P U j (i) −
⩾x
p

U j−1 (i)− Npi

=

Õ

1−

x =1
∞
Õ

ps j  x
N

ps j  x
N
=
N
ps
j
x =1
i
h
N , therefore
By a similar argument with have that E Npi −L j (i) ⩽ ps
j
⩽

E[|γ j |] ⩽ E

h p−1
Õ

1−

i Õ h
i
|Ij (i) ∩ A| =
E U j (i) − L j (i)

i=1

i

i
h
Õ h Ni
Ni i
− L j (i) + E U j (i) −
=
E
p
p
i
Õ 2N
2N
⩽
=
ps
sj
j
i

□

Lemma 4.4 suggests that γ j will be small in expectation. The next
lemma shows that it is also small with high probability.
q
2p
Lemma 4.5. If s j < ln p , then, |γ j | ⩽ 4N
s j w.h.p.
Proof. The main challenge in proving the above theorem is in
handling the dependency in splitter intervals, for e.g. when they
overlap. We first modify the definition of splitter intervals in the
following way, so that the union of the splitter intervals remains
unchanged.
 N (i + 1)


 N (i − 1)
U j′ (i) = min
, U j (i) , L ′j (i) = max
, L j (i)
p
p
The above definition effectively strips the splitter interval
[I (L j (i)), I (U j (i))] to [I (L ′j (i)), I (U j′ (i))]. To see that stripping doesn’t
change the union of all splitter intervals, consider a U j (i) which is
greater than N (i +1)/p. Then by definition, we have U j (i) = U j (i +1).
Thus, the portion of Ij (i) that extends beyond N (i +1)/p is included
in Ij (i + 1). Hence, restricting U j (i) to Ni/p + N /p does not change
the union of Ij ’s, i.e. γ j . One can use an inductive argument to see
that restricting all U j ’s doesn’t change γ j , by considering splitter
intervals from left to right. A similar argument can be used for L j ’s.
Observe that, U j′ (i)’s are independent random variables. This
is because the possible values of U j′ (i 1 ) and U j′ (i 2 ) for i 1 , i 2 are
completely disjoint. The value of U j′ (i) is determined completely
by sampling in the interval [Ni/p, N (i + 1)/p) and since sampling
in disjoint intervals are independent, U j′ (i)’s are independent.
We have, E[U j′ (i) − Ni/p] ≤ E[U j (i) − Ni/p] ≤ N /ps j and using
Hoeffding’s inequality, we obtain
hÕ
hÕ
Ni
2N i
Ni
N
Ni
P
U j′ (i) −
>
=P
U j′ (i) −
−
>
p
sj
p
ps j
sj
i
i
hÕ
i

Ni
Ni  N
≤P
U j′ (i) −
− E U j′ (i) −
>
p
p
sj
i
≤e

−

2p
s i2

1
≤ e −2 ln p ≤ 2
p

HSS

,,

On similar lines we have,

Í

i (Ni/p

− L ′j (i)) ⩽ 2N
s j , w.h.p..

Theorem 4.8. With k = O(log(log p/ϵ)) rounds of histogramming and O(p) sized sample per round (i.e. constant number of samples from each processor), HSS achieves (1 + ϵ) load balance w.h.p. for
large enough p.

We then conclude,

Õ
Ni   Ni
|γ j | ⩽
+
− L ′j (i)
U j′ (i) −
p
p
i
4N
w.h.p.
□
sj
The next lemma bounds the sample size for each round in terms
of the sampling ratios.
⩽

Lemma 4.6. Let Z j be the sample size for the j t h round and s j ⩾
s j−1 , then Z j ⩽ (5ps j /s j−1 ) with high probability.
Proof. We have, E[Z j ] = |γ j−1 |ps j /N . We also have, |γ j−1 | ⩽
4N /s j−1 w.h.p., using Lemma 4.5.
Given that |γ j−1 | ⩽ 4N /s j−1 , using Chernoff bounds,
P[Z j ⩾ (5ps j /s j−1 )] ⩽ P[Z j ⩾ E[Z j ] + ps j /s j−1 ]
⩽e

−

(ps j /s j−1 )2
3E[Z j ]

=e

−

(ps j /s j−1 )2 N
3|γ j−1 |ps j

⩽e

−

(ps j /s j−1 )2 N s j−1
12N ps j

p

⩽ e − 12
□

With Lemmas 4.6 and 4.3 in hand, we are now prepared to discuss how to appropriately choose the sampling ratios so that our
algorithm achieves the desired load balance.
For HSS with k rounds, if we set the sampling ratio for the j th
round as s j = (2 ln p/ϵ)j/k , then after k rounds all splitters are
found w.h.p., using Lemma 4.3. The size of the union of splitter
intervals, that is, |γ j | is less than 4N /s j = 4N .(ϵ/2 ln p)1/k using
Lemma 4.5. The sample size for the j t h histogramming round is at
most 5ps j /s j−1 = 5p(2 ln p/ϵ)1/k using Lemma 4.6. This gives us
our main theorem.
Theorem
4.7. With k rounds of histogramming and a sample size

 q
k log p
of O p
per round , HSS achieves (1 + ϵ) load balance w.h.p.
ϵ
for large enough p 1 .
Observe from theorem 4.7pthat there is a trade off between the
sample size per round (O(p k log p/ϵ)) of histogramming and the
number of histogramming rounds. To minimize the
p number of
samples across all rounds, we take derivative of (kp k log p/ϵ) w.r.t.
k and set it to 0,
p
log p

p
log ϵ 
d(kp k log p/ϵ)
= p k log p/ϵ 1 −
=0
dk
k
log p
⇒ k = log
ϵ
p
The overall sample size O(kp k log p/ϵ) attains global minimum
for k = log(log p/ϵ) rounds of histogramming. We then also have,
|γ j | ≤ 4N /(e)j using Lemma 4.5. Across all rounds, the overall
sample size from all processors across all rounds is O(p log(log p/ϵ)).
This leads us to the following main theorem.
1 Specifically,

so long as s j =

2 ln p
ϵ

≤

q

2p
ln p

for Lemma 4.5

We note that our analysis is slightly conservative. Specifically,
the algorithm need not keep track of splitter intervals for splitters
that are finalized early.

5

RUNNING TIMES

We study the relative parallel execution time complexity of parallel sorting algorithms using the BSP model [28]. The model takes
into account synchronizations (number of supersteps) and measures the communication and computation cost along the critical
path (maximizing over all processors at every supersteps). Since all
the algorithms we consider are bulk-synchronous, the BSP model
provides an accurate relative model of parallel scalability.
We analyse the computational and communication cost separately for sample sort and HSS. Both algorithms have the same cost
for initial local sorting, broadcasting splitters and all-to-all data
exchange. The computational cost of local sorting is O(N log Np /p).
No communication is involved in local sorting. The cost of broadcasting splitters once they are determined is O(p). The final data
movement requires all data to be sent to their destination processors, hence the communication cost involved is O(N /p). Once a
processor receives all data pieces, it needs to merge them, which
takes O(N log p/p) computation time.

5.1

Cost of Sampling

Collecting a sample of overall size S onto one processor, requires a
single BSP superstep with a communication cost of O(S). In practice, random sampling is usually performed with each processor
selecting S/p elements, and a gather collective communication protocol, which collects all samples onto one processor. Sorting the
overall sample on a central processor costs O(S log p) work locally
if each processor provides a sorted contribution to the sample.

5.2

Cost of Histogramming

As the local input is sorted, a local histogram can be computed in
O(S log(N /p)) time using S binary searches, where S denotes the
size of the histogram. A global histogram is computed by reducing
all local histograms. An S-item reduction requires 2 BSP supersteps
(one for a reduce-scatter and one for an all-gather) with O(S) communication and computational cost [23, 27]. The histogram probes
and the splitter intervals are broadcast to every processor for histogramming. The communication cost of broadcasting a length S
message is O(S). So, both the computational and communication
costs of histogramming are proportional to the overall sample size.
The BSP complexity of the data partitioning stage for sample
sort algorithms, AMS sort, and HSS are tabulated in Table 2. AMS
sort and HSS demonstrate the reduction in smaller sample enabled by histogramming. We observe that the best HSS configuration is strictly superior in complexity to all other considered
data-partitioning algorithms.

,,

Vipul Harsh, Laxmikant Kale, and Edgar Solomonik
Algorithm
Regular sampling
Random sampling
Single stage AMS
sort
HSS with one round
HSS with two
rounds

Overall sample
size for
p = 105 , ϵ = 5%

Overall sample
size
p2

O( ϵ )
p log N
O( ϵ 2 )

1600 GB

O(p(log p + ϵ1 ))

32 MB

O(
O(p

p log p
ϵ )

q

184 MB

log p
ϵ )

q

HSS with k rounds

8.1 GB

24 MB

log p
ϵ )

Computational
complexity
 2

p
O ϵ log p log p


p log N log p
O
ϵ2


O p(log p + ϵ1 ) log N


p log p
O
log N
ϵ
 q

log p
O p
ϵ log N
 q

k log p
O kp
ϵ log N


log p
O p log ϵ log N

Communication
complexity
 2
p
O
 ϵ 
p log N
O ϵ2


O p(log p + ϵ1 )


p log p
O
ϵ
 q

log p
O p
ϵ
 q

k log p
O kp
ϵ


log p
O p log ϵ

Supersteps
O(1)
O(1)
O(1)
O(1)
O(1)

O(k)
O(kp
HSS with
log p
log p
10 MB
log p
O(p log ϵ )
O(log ϵ )
O(log ϵ ) rounds
Table 2: Cost complexity of data partitioning part of Sample sort, AMS sort, and HSS. The data exchange costs are excluded.
Algorithm
AMS sort, l stages
HSS, l stages
log r
O(log ϵ ) rounds
per stage

k

Sample size
per stage
O(r (log r +

1
ϵ ))

log r

O(r log ϵ )

Computation complexity
O



N
p

log N + lr (log r +

1
ϵ ) log N

log r
O Np log N + lr log( ϵ ) log N







Communication complexity


O lr (log r + ϵ1 ) + l pN

Supersteps



log(r )
O lr log ϵ + l pN

O(l log ϵ )

O(l)
log p

Table 3: Cost complexity of l -stage HSS and AMS sort; the size of processor group decreases by a factor of r = p 1/l each round.

5.3

HSS with Multiple Stages

HSS can be generalized to a multi-stage algorithm like AMS-sort.
We refer the readers to [3] for more details on multi-stage AMS sort.
We simply consider the benefit of replacing the data-partitioning
step in multi-stage AMS sort with multiple rounds of histogramming as done in HSS. The rest of the algorithm involving data
exchange steps is unchanged. The running time complexities of l
stage AMS-sort and l stage HSS in the BSP model are tabulated in
Table 3. In each step, each processor group get further divided into
r = p 1/l processor groups.


The first local sorting and merging stages take O Np log Np .
At the end of each stage, every processor receives O(r ) pieces of
messages that it needs to merge. So, the total computational cost
of local sorting after every stage excluding the first local sorting is
O(lN log r /p) = O(N log p/p)
computational
 . This gives an overall

cost of local sorting as O Np log Np + Np log p = O(N log N /p).
The computational cost of sampling and histogramming for HSS
log r
is O(r log( ϵ ) log N ) per stage. Each histogramming, sampling,
and data exchange step can be done with O(1) BSP supersteps. The
number and cost of all of these steps is uniform throughout stages,
so all of these costs are multiplied by a factor of l (corresponding
to the number of stages). Consequently, we observe a trade-off
between the cost of data partitioning and the cost of data exchanges
that depends on l.

6

IMPLEMENTATION

We implemented HSS in C++11 using the Charm++ [1, 5, 17] framework. First we describe some optimisations exploiting shared memory and then discuss implementation specifics.

6.1

Node-level optimizations

The all-to-all data exchange step results in p(p − 1) messages in
the network, which doesn’t scale very well. We leverage shared
memory within a physical node to reduce the number of messages.

We modified HSS to be a simple two phase algorithm. In the first
phase, the algorithm performs a node level splitting placing all data
onto their destination nodes. In the second phase, every node sorts
its local data, for instance, using OpenMP or Charm++ features.
Data partitioning across nodes reduces the size of the histogram
considerably, since p now reflects the number of nodes rather than
processors. For instance BlueGene/L system [12] has 16 cores in a
single physical node. Consider partitioning across 8K nodes versus
128K processors (assuming 1 processor per core). Once a node
receives all the data after the all to all data exchange, it can locally
sort the data within node without injecting traffic on the network.

6.2

Implementation specifics

Our implementation comprises of three distinct phases; local sorting of input data, splitter determination using histogramming and
final data exchange and assembly. We use STL’s std::sort for
local sorting for the first phase. Let p denote the total number of
processors (threads) and n denote the number of nodes (ranks).
• Histogramming Phase: The histogramming phase determines
n − 1 splitters for node level splitting. For the sampling phase
before every histogramming round, each processor picks a small
sample from its input which lie in the union of splitter intervals.
If δ denotes the fraction of input covered by the splitter intervals,
then every processor picks s/δ samples from its entire input and
discards samples that don’t lie in any of the splitter intervals. This
way the expected size of the overall sample from every processor
is s × p, where s can be thought of as the oversampling ratio. The
overall sample is assembled at the central processor and broadcast
for histogramming. Every processor computes a local histogram
using binary searches, since the input data is already sorted. The
local histograms are summed up using a reduction and sent to
the central processor. We set the load balance threshold to 5% for
node level partitioning and 5% for within node partitioning.
• Data Exchange: Once every processor receives the splitters, input
data from all processors within a node are merged and partitioned

HSS

,,
7

data exchange
histogramming
local sort

6
5
4

Execution time (sec)

Execution time (sec)

data exchange
histogramming
local sort

20

3
2
1

15

10

5

K
32

6.0

Execution Time(sec)

ro
s
Al
lZ
e

GA
US
S

SK
EW
3

7.1

4.5

We ran weak scaling experiments with 2 million keys per processor.
We compare our results to HykSort [26]. We used k = 128 as the
parameter for HykSort as recommended in the algorithm specification [26]. We used 16 threads per rank for both HSS and HykSort.
A summary of our experiments are reported in Figure 4. HykSort
suffers from multiple stages of data movement step. As the number
of processors increase however, the multiple stages of data movement is likely to pay off. From Figure 4, we conclude that won’t
happen till a large problem size. As discussed earlier, the results in
this paper are equally applicable to a multi stage algorithm.

4.0

3.0
1K

2K

4K

8K
#PEs

16K

32K

64K

Figure 4: HSS vs 128 way HykSort, with uniformly distributed 2M
keys/processor and 16 threads per rank.

into n messages, one for each node. We use Charm++’s nodegroup
chares to implement the combining of messages within a node.
• Local Merging: Once a node receives all data that falls in its
bucket, it needs to merge and redistribute data among its processors. We use HSS with one round for final within node sorting.

Handling duplicates via implicit tagging

We use a standard technique to deal with duplicates. For inputs with
duplicates, we enforce a strict ordering by implicitly replacing each
key k with a triplet (k, processor , ind), where processor denotes the
processor that k resides on and ind denotes the index in the local
data structure, where k is stored. Note that implicit tagging does
not blow up the size of the input, but increases the size of the histogram by a constant factor since probe keys for the histogram have
to be explicitly tagged. Figure 3 illustrates the overhead incurred
for handling duplicates with implicit tagging. For 32K processors,
implicit tagging resulted in only about 4% increase in running time.

7

SK
EW
2

5.0

3.5

6.3

Figure 5: Performance of HSS with different input distributions for 2M
keys/processor. We used 1 thread per rank to accentuate the splitter selection
(histogramming) time.

processor for the splitting round. We found this to be a reasonable
value of the oversampling ratio.

128 way HYK Sort
HSS

5.5

SK
EW
1

0

Figure 3: Overhead of handling duplicates for UNIF distribution with 2M
keys/processor and 16 threads/rank. The histograms on the right are the
runs with implicit tagging.

UN
IF

K
16

92
81

96
40

20

48

0

EXPERIMENTAL RESULTS

All our experiments were performed on the Vesta IBM Blue Gene/Q
supercomputer at Argonne Leadership Computing Facility (ALCF).
Each node is equipped with a PowerPC A2 1600 MHz processor containing 16 cores, each with 4 hardware threads, running at 1.6 GHz,
and 16 GB of DDR3 memory. The nodes are connected by IBM’s
5D torus interconnect. We use 8 byte long integer type keys for
all our experiments. Unless stated otherwise, we use oversampling
ratio of 1 per processor, i.e. 1 key on average is sampled from every

7.2

Weak scaling

Effect of input distribution

The theoretical guarantees for HSS are valid for all distributions. To
verify the same, we ran HSS with the following input distributions:
(1) UNIF: Uniformly at random from the entire range
(2) SKEW1: Half of the keys are picked uniformly at random from
the entire range, the other half, uniformly at random from a
small range of size 1000
(3) SKEW2: Uniformly at random from the range [0, 100]
(4) SKEW3: Each key is bitwise and of two uniformly at random
chosen keys
(5) GAUSS: Gaussian distributed
(6) AllZeros: All keys are set to 0
As figure 5 illustrates, HSS is impervious to the input distribution
as also expected from the theory. To underscore the histogramming
cost, we used 1 thread per rank. Recall that the number of splitters
is equal to the number of ranks.
7.2.1 Histogramming time with varying processors. Figure 6 illustrates the histogramming cost with varying number of processors.
We used one thread per rank to underscore the cost of histogramming. The histogramming cost roughly increases by a factor of 2
with a 2x increase in the number of processors as also expected
log p
from the theoretical O(p log ϵ ) bound (Table 2).
With 16 threads per rank, the histogramming cost is much less
because of fewer number of splitters to determine. With 16 threads
per rank, we observe that even for large number of processors, the

,,

Vipul Harsh, Laxmikant Kale, and Edgar Solomonik
p(×103 )

8.0

Histogramming
Data Exchange

4
8
16
32

Execution Time(sec)

4.0

Bound on number
of rounds
8
8
8
8

1.0

1K

2K

4K

8K
#PEs

16K

32K

64K

Figure 6: Histogramming time with increasing number of processors for
2M keys per processor and 1 thread per rank.
16.0

HSS Lambb
Histogram Sort Lambb
HSS Dwarf
Histogram Sort Dwarf

8.0

Execution Time(sec)

Number of
rounds
4
4
4
4

Table 4: Number of histogramming rounds observed with ϵ = 0.02, 1M
keys per processor and 1 thread per rank.

2.0

0.5

4.0

1.0

256

so Figure 7 represents strong scaling results. We saw up to 25%
improvement using HSS over plain Histogram sort. The parallel
sorting execution increases for the same dataset as we increase the
number of processors. This may appear odd at first. The majority
of sorting time is spent in data splitting, and since the number of
buckets increase multiplicatively with the number of processors,
we see an increase in the execution time. The performance results
suggest it would be possible to improve strong scaling of the splitting algorithm within ChaNGa by using a multi-staged version of
HSS. We leave this for future work.

8

2.0

0.5.0
512

1K

2K

4K
8K
#processors

16K

32K

64K

Figure 7: Experiments showing performance of sorting routine of Changa.
Datasets used were Lambb and Dwarf

histogramming phase takes a very small fraction of the running time
(Figure 3). We believe this is both because of efficient histogramming
and node level partitioning.
The number of histogramming rounds taken by the algorithm is
tabulated in Table 4. If the sample size is (p f ) every round, the number of histogramming rounds
required to determine (p − 1) splitters

is ln(2 ln p/ϵ)/ln(f /2) . We obtain this expression using expected
value of sample sizes and length of splitter intervals ( Theorem 4.7).
As discussed earlier, the constants involved in our analysis are conservative. In practice, HSS required fewer histogramming rounds.
The number of rounds: O(log(log p/ϵ)) increases very slowly with
p and hence, HSS is extremely scalable and practical.

7.3

sample size per
round (×p)
5
5
5
5

Parallel sorting in ChaNGa

We implemented HSS in ChaNGa [16], a popular astronomical simulation software that often runs on several thousands of processors.
ChaNGa poses unique challenges for sorting for two reasons. First,
it employs virtual processors and hence the number of buckets (virtual processors) are far more than the actual number of processors.
In our examples the number of virtual buckets were typically 10x
the number of cores. Second, the virtual processors can be arbitrarily placed on any physical processor and buckets on a single
node might not be contiguous. Hence, none of our shared memory
optimizations are useful. The reason ChaNGa uses more virtual
processors than cores is to accelerate other stages of computation,
which is made possible by efficient parallel data overpartitioning.
Figure 7 illustrates performance of ChaNGa with HSS and already existing Histogram sort implementation for two datasets:
Dwarf and Lambb. We refer the readers to [16] for more details
about the datasets. The datasets have a constant number of keys,

CONCLUSIONS

We present Histogram sort with sampling (HSS), that combines sampling and histogramming to accomplish fast splitter determination.
Specifically, we show that when approximate-splitting is sufficient
and ϵ = O(1), our algorithm requires Θ(log log p) rounds of histogramming and the collection of Θ(p log log p) overall samples, improving the communication cost by a factor of Θ(log(p)/log log p)
with respect to the best known partitioning algorithm. Our study
provides theoretical groundwork for the benefits of iterative histogramming in splitter selection, a technique that has been known
to work well in practice. The reduced sample size makes HSS extremely practical for massively parallel applications, scaling to
tens of thousands of processors. We provide an implementation
of the HSS algorithm, demonstrating speed-ups over two other
state-of-the-art parallel sorting implementions. The robustness of
our results makes a compelling case for HSS as the algorithm of
choice for parallel sorting libraries. The algorithm is theoretically
more efficient for both approximate and exact (memory-efficient)
splitting and minimizes the number of data exchanges with both
small and large degrees of parallelism.

ACKNOWLEDGMENTS
The authors are grateful to Omkar Thakoor for proofreading the
analysis and Nitin Bhat in helping with experiments. The authors
acknowledge the Argonne Leadership Computing Facility (ALCF)
for providing HPC resources that have contributed to the research
results reported within this paper. URL: https://www.alcf.anl.gov/.

REFERENCES
[1] Bilge Acun, Abhishek Gupta, Nikhil Jain, Akhil Langer, Harshitha Menon, Eric
Mikida, Xiang Ni, Michael Robson, Yanhua Sun, Ehsan Totoni, et al. 2014. Parallel
programming with migratable objects: charm++ in practice. In Proceedings of the
International Conference for High Performance Computing, Networking, Storage
and Analysis (SC). IEEE, 647–658.
[2] Miklós Ajtai, János Komlós, and Endre Szemerédi. 1983. An O(n log n) sorting network. In Proceedings of the 15th annual ACM Symposium on Theory of
computing. ACM, 1–9.
[3] Michael Axtmann, Timo Bingmann, Peter Sanders, and Christian Schulz. 2015.
Practical massively parallel sorting. In Proceedings of the 27th ACM Symposium
on Parallelism in Algorithms and Architectures. ACM, 13–23.
[4] Kenneth E Batcher. 1968. Sorting networks and their applications. In Proceedings
of the April 30–May 2, 1968, spring joint computer conference. ACM, 307–314.

HSS

,,

[5] Abhinav Bhatele. 2014. Chapter 5.2 Projections: Scalable Performance Analysis
and Visualization. Technical Report: Connecting Performance Analysis and Visualization to Advance Extreme Scale Computing, Lawrence Livermore National
Laboratory (2014), 33.
[6] Guy E Blelloch, Charles E Leiserson, Bruce M Maggs, C Greg Plaxton, Stephen J
Smith, and Marco Zagha. 1991. A comparison of sorting algorithms for the
connection machine CM-2. In Proceedings of the third annual ACM Symposium
on Parallel algorithms and architectures. ACM, 3–16.
[7] Guy E. Blelloch, Charles E. Leiserson, Bruce M Maggs, C Greg Plaxton, Stephen J
Smith, and Marco Zagha. 1998. An experimental analysis of parallel sorting
algorithms. Theory of Computing Systems 31, 2 (1998), 135–167.
[8] David R Cheng, Alan Edelman, John R Gilbert, and Viral Shah. 2006. A novel
parallel sorting algorithm for contemporary architectures. (2006).
[9] Richard Cole. 1988. Note on the AKS sorting network. Computer Science Department Technical Report 243. New York University, New York, (1988).
[10] Richard Cole. 1988. Parallel merge sort. SIAM J. Comput. 17, 4 (1988), 770–785.
[11] W Donald Frazer and AC McKellar. 1970. Samplesort: A sampling approach to
minimal storage tree sorting. Journal of the ACM (JACM) 17, 3 (1970), 496–507.
[12] Alan Gara, Matthias A Blumrich, Dong Chen, GL-T Chiu, Paul Coteus, Mark E
Giampapa, Ruud A Haring, Philip Heidelberger, Dirk Hoenicke, Gerard V Kopcsay,
et al. 2005. Overview of the Blue Gene/L system architecture. IBM Journal of
Research and Development 49, 2 (2005), 195–212.
[13] Michael T Goodrich. 1999. Communication-efficient parallel sorting. SIAM J.
Comput. 29, 2 (1999), 416–432.
[14] David R Helman, Joseph JáJá, and David A Bader. 1998. A new deterministic parallel sorting algorithm with an experimental evaluation. Journal of Experimental
Algorithmics (JEA) 3 (1998), 4.
[15] JS Huang and YC Chow. 1983. Parallel sorting and data partitioning by sampling.
(1983).
[16] Pritish Jetley, Filippo Gioachin, Celso Mendes, Laxmikant V Kale, and Thomas
Quinn. 2008. Massively parallel cosmological simulations with ChaNGa. In
Proceedings of the International Parallel and Distributed Processing Symposium
(IPDPS). IEEE, 1–12.
[17] Laxmikant V Kale and Sanjeev Krishnan. 1993. CHARM++: a portable concurrent
object oriented system based on C++. Vol. 28. ACM.
[18] Laxmikant V Kale and Sanjeev Krishnan. 1993. A comparison based parallel sorting algorithm. In Proceedings of the International Conference on Parallel Processing
(ICPP), Vol. 3. IEEE, 196–200.
[19] Hui Li and Kenneth C Sevcik. 1994. Parallel sorting by over partitioning. In
Proceedings of the sixth annual Symposium on Parallel algorithms and architectures
(SPAA). ACM, 46–56.
[20] Xiaobo Li, Paul Lu, Jonathan Schaeffer, John Shillington, Pok Sze Wong, and Hanmao Shi. 1993. On the Versatility of Parallel Sorting by Regular Sampling. Parallel
Comput. 19, 10 (Oct. 1993), 1079–1103. https://doi.org/10.1016/0167-8191(93)
90019-H
[21] Owen O’Malley. 2008.
Terabyte sort on apache hadoop.
Yahoo,
http://sortbenchmark.org/YahooHadoop.pdf (2008), 1–3.
[22] M. S. Paterson. 1990. Improved sorting networks with O(logN) depth. Algorithmica 5, 1 (1990), 75–92. https://doi.org/10.1007/BF01840378
[23] Jelena Pješivac-Grbović, Thara Angskun, George Bosilca, Graham E. Fagg, Edgar
Gabriel, and Jack J. Dongarra. 2007. Performance analysis of MPI collective
operations. Cluster Computing 10, 2 (2007), 127–143. https://doi.org/10.1007/
s10586-007-0012-0
[24] Hanmao Shi and Jonathan Schaeffer. 1992. Parallel Sorting by Regular Sampling.
J. Parallel and Distrib. Comput. 14, 4 (April 1992), 361–372. https://doi.org/10.
1016/0743-7315(92)90075-X
[25] E. Solomonik and L. V. Kale. 2010. Highly scalable parallel sorting. In Proceedings
of the International Parallel and Distributed Processing Symposium (IPDPS). 1–12.
https://doi.org/10.1109/IPDPS.2010.5470406
[26] Hari Sundar, Dhairya Malhotra, and George Biros. 2013. HykSort: A New
Variant of Hypercube Quicksort on Distributed Memory Architectures. In Proceedings of the 27th International ACM Conference on International Conference
on Supercomputing (ICS ’13). ACM, New York, NY, USA, 293–302. https:
//doi.org/10.1145/2464996.2465442
[27] Rajeev Thakur and William D Gropp. 2003. Improving the performance of
collective operations in MPICH. In European Parallel Virtual Machine/Message
Passing Interface Users Group Meeting. Springer, 257–267.
[28] Leslie G. Valiant. 1990. A Bridging Model for Parallel Computation. Commun.
ACM 33, 8 (Aug. 1990), 103–111. https://doi.org/10.1145/79173.79181

A

APPENDIX

We formalize the proof sketch given by Axtmann et al. [3] for
their parallel sorting algorithm. Their algorithm scans the results
of a histogram collected from a randomly selected sample, picking
splitters consecutively so that each processor gets no more than (1+

ϵ)N /p elements. If the sample is large enough, with high probability
the average load assigned to each consecutive processor will be
greater than N /p and the last processor will not end up with too
many keys. We provide a formal bound on the probability that this
holds for the appropriate sample size.
Lemma A.1. If every key is independently picked in the sample
with probability, ps/N = p max(2/ϵ, 2 log N )/N , then the scanning
algorithm assigns no more than N (1 + ϵ)/p keys to each processor
with high probability.
Proof. Let χ (0), . . . , χ (p) be the indices in the global order of
the splitter keys selected by the scanning algorithm (with χ (0) = 0
and χ (p) = N ) and let ni = χ (i + 1) − χ (i) be the resulting number
of keys asisgned to processor i. By construction the algorithm
assigns at most ni ≤ (1 + ϵ)N /p for i < p − 1, so it suffices to show
np ≤ (1 + ϵ)N /p.
Select all unique splitter indices { χ̃0 , . . . , χ̃h } from {χ (0), . . . , χ (p − 1)},
so h < p. Note that we exclude χ (p) and hence χ̃i+1 − χ̃i < N (1+ϵ)/p
for all i. Splitters can repeat if there are more than (1 + ϵ)N /p keys
between any pair of splitters (including splitters with index 0 and
N ), in which case the algorithm essentially fails. The only other
way the scanning algorithm yields repeated splitters is if it runs out
of sample keys to assign to pick as splitters, in which case np−1 = 0.
Without loss of generality, we assume the latter does not happen
and show that the former is improbable by showing that w.h.p. a
sample is selected from each subrange of size N /p,

ps  Np 
P(h = p − 1) ≥ 1 − N × 1 −
N



N ps 
− pN
= 1 − N × e −s
≥ 1− N ×e


1
≥ 1 − N × e −2 log N ≥ 1 −
N
We proceed to show that if h = p − 1, then np−1 is small with high
probability.
Let ti be the next key in the sample after χ̃i for i ≤ p. We have
that χ̃i+1 ≥ ti for i < p. Given the location of splitter χ̃i−1 , we
have that χ̃i ≤ oi = χ̃i−1 + (1 + ϵ)N /p. Since ti and oi depend only
on whether keys in the range [0, ti ] are present in the sample, the
probability of any key in the range (ti , oi ] being in the sample, is
independent of ti and oi . Define r i = oi − χ̃i = (1 + ϵ)N /p − ( χ̃i −
χ̃i−1 ). Each r i ≤ oi −ti is a random variable that depends exclusively
on the existence and location of a sample key in (ti , oi ] and the size
of the interval (oi − ti ). Since the probability of observing keys in
this interval is independent of ti , r i can depend on {r 1 , . . . , r i−1 }
only as a function of the size of the interval.
We can express each r i as min(oi − ti , r˜i ) where r˜i is a random
variable corresponding to the distance between the closest key in
[0, (1 + ϵ)N /p] to (1 + ϵ)N /p when keys in this range are picked
with uniform probability ps/N . Then we observe that r i ≤ r˜i and
the random variables r˜i are mutually independent. Since r˜i s are
picked independently using binomial
trials,

 they are exponentially
distributed. More specifically, P r˜i ≥ k = (1 − ps/N )k . Thus, we
have,
∞
 
  Õ
ps  k
N
E r i ≤ E r˜i ≤
1−
≤
N
ps
k =1

,,

Vipul Harsh, Laxmikant Kale, and Edgar Solomonik

Coming back to bounding the number of keys assigned to the
last processor by the scanning algorithm np−1 , we have,
p−2
Õ

p−2
Õ


N Nϵ
+
− ri
p
p
i=0
i=0
p−2
N N ϵ(p − 1) Õ
=
−
+
ri
p
p
i=0

np−1 = N −

ni = N −

Using sampling ratio s = ϵ2 we have,
h h−1
i
h h−1
i
Õ
Õ
P
ri ≥ N ϵ ≤ P
r˜i ≥ N ϵ
i=0

i=0

h h−1
Õ

N ϵ N ϵ(h − 1) i
+
2
2p
i=0
" h−1
#
h−1
Õ
Õ  Nϵ
≤P
r˜i − E[
r˜i ] ≥
2
i=0
i=0
≤P

r˜i ≥

−pϵ 2

≤ e 2(1+ϵ )2
The last inequality is obtained using an application of Hoeffding’s
inequality since all r˜i ’s are independent.
With high probability, we have h = p −1. Hence using expression
for np−1
h
h p−2
i
Õ
N Nϵ i
P np−1 ≥
+
≤ P
ri ≥ N ϵ
p
p
i=0
−pϵ 2

≤ e 2(1+ϵ )2 w.h.p.
□

