Samplesort

Samplesort ist ein Sortieralgorithmus, der 1970 von W. D. Frazer und A. C. McKellar in der wissenschaftlichen Publikation Samplesort: A Sampling Approach to Minimal Storage Tree Sorting veröffentlicht wurde. Der Algorithmus arbeitet ähnlich wie Quicksort und verfolgt wie dieser das Teile-und-herrsche-Verfahren. Der Algorithmus teilt die Eingabesequenz in mehrere Teilsequenzen auf, die dann rekursiv sortiert werden.

Pseudocode

Der folgende Pseudocode stellt die Arbeitsweise des Algorithmus grundlegend dar. Dabei enthält $A$ die zu sortierenden Daten, $k$ entspricht dem Oversamplingfaktor und $p$ der Anzahl der zu nutzenden Splitter. Die Anzahl der Splitter gibt an, in wie viele Teilsequenzen die Eingabedaten aufgeteilt werden (für $p=2$ entspricht der Algorithmus Quicksort), während sich durch einen großen Oversamplingfaktor die Größen der Buckets statistisch gesehen weniger unterscheiden.

${\text{sampleSort}}(A[1..n],k,p)$
  ${\text{if }}n/k<{\text{ threshold then smallSort}}(A)$ // if average bucket size is below a threshold switch to e.g. quicksort
  ${\text{select }}S=[S_{1},\dots ,S_{p(k-1)}]{\text{ randomly from }}A$ // select samples
  ${\text{sort }}S$ // sort sample
  $[s_{0},s_{1},\dots ,s_{p-1},s_{p}]\leftarrow [-\infty ,S_{k},S_{2k},\dots ,S_{(p-1)k},\infty ]$ // select splitters
  // classify elements
  ${\text{for each }}a\in A$
    ${\text{find }}j{\text{ such that }}s_{j-1}<a\leq s_{j}$
    ${\text{place }}a{\text{ in bucket }}b_{j}$
  ${\text{return concatenate}}({\text{sampleSort}}(b_{1}),\dots ,{\text{sampleSort}}(b_{k}))$

Die Wahl des richtigen Buckets in der Schleife kann für jedes Element mit binärer Suche in Zeit ${\mathcal {O}}(\log p)$ ausgeführt werden.

Der angegebene Pseudocode unterscheidet sich von der von Frazer und McKellar beschriebenen Ursprungsform. Im Pseudocode wird Samplesort rekursiv aufgerufen, um die erzeugten kleineren Teilmengen zu sortieren. Frazer und McKellar verwendeten dafür Quicksort, d. h. die Funktion Samplesort wird nur einmal aufgerufen.

Die Anzahl Vergleiche, die dieser Algorithmus durchführt, nähert sich für große Eingabesequenzen asymptotisch dem informationstheoretischen Optimum $\log _{2}(n!)$ an und brauchte in den von Frazer und McKellar durchgeführten Experimenten mehr als 15 % weniger Vergleiche als Quicksort.

Oversampling

Der Oversamplingfaktor gibt an, wie groß die Stichprobe zur Auswahl der Splitter größer ist, als die Anzahl der Splitter. Ziel ist es, die Verteilung der Eingabedaten möglichst gut anzunähern. Nach der Klassifikation enthält im Idealfall jeder Bucket genau $n/p$ Elemente, um die folgende Arbeit des Algorithmus gleichmäßig auf die Buckets aufzuteilen.

Nach der Wahl des Samples $S$ der Größe $p\cdot (k-1)$ wird das Sample sortiert. Als Splitter, die als Bucketgrenzen verwendet werden, werden die Elemente in $S$ an den Stellen $k,2k,3k,\dots ,(p-1)k$ sowie $-\infty$ und $\infty$ als linke und rechte Grenzen des linkesten und rechtesten Buckets verwendet.

Größe der Buckets

Mit der Größe des Samples bzw. dem Oversamplingfaktor kann die erwartete Bucketgröße und insbesondere die Wahrscheinlichkeit dafür, dass ein Bucket eine gewisse Größe überschreitet, abgeschätzt werden. Für einen Oversamplingfaktor $k\in \Theta ({\frac {\log n}{\epsilon ^{2}}})$ werden wir zeigen, dass mit einer Wahrscheinlichkeit von weniger als $1/n$ ein Bucket mehr als $(1+\epsilon )n/p$ Elemente enthält.

Sei dafür $[e_{1},\dots ,e_{n}]$ die Eingabe bereits in sortierter Form. Damit ein Prozessor mehr als $(1+\epsilon )n/p$ Elemente erhält, muss es eine Sequenz der Eingabe der Länge $(1+\epsilon )n/p$ geben, aus der maximal $k$ Elemente für das Sample gewählt werden. Das wird wie folgt als Zufallsvariable dargestellt:

$X_{i}:={\begin{cases}1,&{\text{falls }}s_{i}\in [e_{j},\dots ,e_{j}+(1+\epsilon )n/p]\\0,&{\text{sonst}}\end{cases}}$

$X:=\sum _{i=0}^{k\cdot p-1}X_{i}$

Daraus folgt:

$\mathbb {E} [X_{i}]=\mathbb {P} [X_{i}=1]={\frac {1+\epsilon }{p}}$

$\mathbb {P} [X<k]\approx \mathbb {P} [X<(1-\epsilon ^{2})k]=\mathbb {P} [X<(1-\epsilon )\mathbb {E} [X]]$

Dabei bezeichnet $\mathbb {P} [X<k]$ die Wahrscheinlichkeit, dass weniger als $k$ Samples im Bucket $j$ liegen. Mithilfe der Chernoff-Ungleichung sieht man, dass

$n\mathbb {P} [X<k]\leq n\exp({\frac {\epsilon ^{2}k}{2}})\leq n{\frac {1}{n^{2}}}$ für $k\geq {\frac {4}{\epsilon ^{2}}}\ln n$

Samplesort auf Systemen mit verteiltem Speicher

Arbeitsprinzip von parallelem Samplesort für $p=3$ Prozessoren und einem Oversamplingfaktor $k=3$.

Samplesort wird oft in parallelen Systemen mit verteiltem Speicher verwendet, da es durch die variable Anzahl von Buckets im Gegensatz zum Beispiel zu Quicksort einfach an eine beliebige Anzahl Prozessoren angepasst werden kann. Die Eingabedaten liegen dabei verteilt auf $p$ Prozessoren vor, das heißt, jeder Prozessor hat $n/p$ Elemente gespeichert. Entsprechend werden ebenso viele Splitter ausgewählt, um für jeden Prozessor einen Bucket zu erzeugen.

Um die Splitter auszuwählen, ist es ausreichend, wenn jeder Prozessor aus seinen Eingabeelementen $k$ Elemente auswählt. Diese werden mit einem verteilten Sortieralgorithmus sortiert und an alle Prozessoren verteilt. Es entstehen Kosten von $T_{\text{sort}}(kp,p)$ für das Sortieren der $kp$ Elemente auf $p$ Prozessoren, sowie $T_{\text{allgather}}(p,p)$ für das Verteilen der $p$ gewählten Splitter auf $p$ Prozessoren.

Mit den erhaltenen Splittern partitioniert jeder Prozessor die eigenen Eingabedaten in lokale Buckets. Die Zeit dafür beträgt mit binärer Suche ${\mathcal {O}}(n/p\log p)$. Anschließend verteilen die Prozessoren die Buckets an die jeweiligen Prozessoren mit Zeit $T_{\text{all-to-all}}(N,p)$, wobei $N$ der Größe des größten Buckets entspricht. Prozessor $i$ enthält entsprechend von allen Prozessoren die lokalen Buckets $b_{i}$ und sortiert diese lokal. Der Zeitbedarf für die lokale Sortierung von allen Prozessoren beträgt $T_{\text{localsort}}(N)$.

Wählt man wie im vorherigen Abschnitt beschrieben $k\in \Theta ({\frac {\log n}{\epsilon ^{2}}})$, so erhält man eine Gesamtlaufzeit von

$T(n,p)=T_{\text{sort}}({\mathcal {O}}({\frac {p\log n}{\epsilon ^{2}}}),p)+T_{\text{allgather}}(p,p)+{\mathcal {O}}(n/p\log p)+T_{\text{all-to-all}}((1+\epsilon )n/p,p)+T_{\text{localsort}}((1+\epsilon )n/p)$

Effiziente Implementierung von sequentiellem Samplesort

Arbeitsprinzip von Super Scalar Sample Sort. (Animation) In jedem Schritt werden Zahlen, die verglichen werden, blau eingefärbt und Zahlen, auf die anderweitig lesend oder schreibend zugegriffen wird, werden rot eingefärbt.

Wie oben beschrieben klassifiziert der Algorithmus die Elemente anhand der gewählten Splitter und partitioniert diese entsprechend. Eine effiziente Implementierung wird im Paper Super Scalar Sample Sort beschrieben. Dort werden zwei Arrays der Größe $n$ genutzt: ein Array für die Eingabedaten, sowie ein temporäres Array. Entsprechend arbeitet Samplesort in diesem Fall nicht in-place. In jedem Rekurssionsschritt werden die Daten von einem zum anderen Array kopiert und währenddessen partitioniert. Falls die Daten nach dem letzten Rekursionsschritt im temporären Array vorliegen, so werden diese in das originale Array zurückkopiert. Die wesentlichen Schritte sind die Klassifikation der Elemente anhand der Splitter und die anschließende Partitionierung. Beides wird im Folgenden genauer erklärt.

Klassifikation der Elemente

In einem vergleichsbasierten Sortieralgorithmus ist der performancekritische Teil die Klassifikation der Elemente, d. h. die Zuordnung von Elementen zu Buckets. Diese Zuordnung benötigt ${\mathcal {O}}(\log p)$ Schritte pro Element. Super Scalar Sample Sort nutzt einen balancierten Suchbaum, der implizit in einem Array $t$ gespeichert ist. Das heißt, die Wurzel befindet sich an Stelle 0, und der linke Nachfolger eines Elementes $t_{i}$ befindet sich in $t_{2i}$, der rechte Nachfolger in $t_{2i+1}$. Bei einem gegebenen Suchbaum $t$ und einem Element $a$ berechnet der Algorithmus den Index des entsprechenden Buckets $j$ wie folgt (der Ausdruck $a>t_{j}$ sei 1 falls $a>t_{j}$ und 0 andernfalls, $p$ sei eine Zweierpotenz):

$j\leftarrow 1$
${\text{repeat}}\log _{2}p{\text{ times}}$
  $j\leftarrow 2j+(a>t_{j})$
$j\leftarrow j-p+1$

Da die Anzahl der Buckets $p$ bei der Übersetzungszeit bekannt ist, kann die Schleife vom Compiler ausgerollt werden. Falls die Prozessorarchitektur es unterstützt, kann in der Schleife ganz auf bedingte Sprünge verzichtet werden und es treten keine fehlerhaften Sprungvorhersagen auf.

Partitionierung

Zum effizienten Partitionieren muss der Algorithmus die Größen der Buckets vor dem Kopieren der Elemente kennen. In einer naiven Implementierung könnten alle Elemente klassifiziert werden, um die Bucketgrößen zu bestimmen. Anschließend könnten die Elemente an die richtige Stelle kopiert werden. So müsste der Bucket für jedes Element zweimal bestimmt werden (einmal für das Zählen der Elemente der Buckets und einmal zum Kopieren). Super Scalar Sample Sort vermeidet diese doppelte Klassifikation, indem das Ergebnis der Klassifikation in einem zusätzlichen Array $o$ (dem „Orakel“) zwischengespeichert werden. Dieses ordnet jedem Element $a_{i}$ den Bucketindex $o_{i}$ zu. Zuerst bestimmt Super Scalar Sample Sort die richtigen Werte von $o$, indem jedes Element klassifiziert wird. Danach werden die Größen der Buckets bestimmt und jedes Element anhand von $o$ an die richtige Stelle kopiert. Der Speicherbedarf von Array $o$ beträgt $n\log p$ bits, was im Vergleich zu den Eingabedaten nur wenig Speicher bedeutet.

Paralleler in-place-Samplesort für Systeme mit gemeinsamem Speicher

Die oben gezeigte Samplesort-Implementierung arbeitet nicht in-place, da ein zweites Array der Größe $n$ benötigt wird. Effiziente Implementierungen von z. B. Quicksort arbeiten in-place und benötigen daher weniger Speicher. Samplesort kann allerdings auch in-place für parallele Systeme mit gemeinsamem Speicher implementiert werden. Im Folgenden wird eine solche Implementierung beschrieben. Zu beachten ist, dass der zusätzliche Speicherplatz des Algorithmus immer noch von $n$ abhängt, im Gegensatz zum oben erklärten Super Scalar Sample Sort aber nur ${\mathcal {o}}(n)$ zusätzlichen Speicherplatz benötigt. Eine modifizierte Variante des Algorithmus, der nur noch Speicherplatz unabhängig von $n$ benötigt, ist ebenfalls in beschrieben.

Der in-place Algorithmus kann in vier Phasen aufgeteilt werden:

1. Sampling: Auswahl der Splitter.
2. Lokale Klassifikation: Gruppieren der Eingabe in Blöcke, so dass alle Elemente eines Blockes zu einem Bucket gehören (allerdings muss nicht jeder Bucket zusammenhängend im Speicher liegen).
3. Blockpermutation: Vertauschen der Blöcke, so dass diese global in der richtigen Reihenfolge liegen.
4. Aufräumen: Kopieren der übrigen Elemente an die Kanten der Buckets.

Da durch diesen Algorithmus jedes Element zweimal gelesen und geschrieben werden muss (einmal bei der Klassifikation und einmal während der Blockpermutation) ergibt sich ein höherer Schreib- und Leseaufwand. Trotzdem ist der Algorithmus bis zu dreimal schneller als aktuelle in-place Wettbewerber und bis zu 1,5 mal schneller als andere aktuelle sequentielle Wettbewerber, einschließlich des oben beschrieben Super Scalar Sample Sort. Da das Sampling und die Auswahl der Splitter bereits oben (Pseudocode) beschrieben ist, folgen an dieser Stelle nur noch die detaillierten Erklärungen der drei anderen Phasen.

Lokale Klassifikation

Das Eingabearray wird zuerst in $t$ (Anzahl der Prozessoren/Threads) Bereiche gleicher Größe aufgeteilt. Jeder Prozessor enthält einen solchen Bereich und teilt diesen weiter auf in Blöcke, die jeweils $b$ Elemente enthalten. Zusätzlich allokiert jeder Prozessor $p$ (Anzahl Buckets) Puffer, die ebenfalls $b$ Elemente speichern können. Jeder Prozessor verarbeitet anschließend seinen zugewiesenen Bereich und kopiert die Elemente in die entsprechenden allokierten Puffer. Falls ein Puffer voll ist, wird dessen Inhalt in den vom Prozessor verarbeiteten Bereich kopiert. Die dortigen $b$ Elemente, die dadurch überschrieben werden, sind bereits verarbeitet, da ansonsten kein Puffer voll sein könnte. Zur Auswahl der richtigen Buckets/Puffer siehe Klassifikation der Elemente. Nach dieser Phase enthält jeder Bereich mehrere Blöcke mit Elementen jeweils eines Buckets und darauffolgend leere Blöcke. Die übrigen Elemente sind noch in den Pufferblöcken. Zusätzlich speichert jeder Prozessor die Größe jedes Buckets.

Blockpermutation

In dieser Phase werden die Blöcke im Eingabearray in die richtige Reihenfolge gebracht. Mithilfe einer Präfixsumme über die Größe der Buckets werden deren Grenzen bestimmt. Da in dieser Phase nur ganze Blöcke verschoben werden, werden die Grenzen auf ganze Blöcke aufgerundet. Die fehlenden Elemente befinden sich noch in den Pufferblöcken. Da die Größe des Eingabearrays nicht zwangsläufig durch die Blockgröße teilbar ist, müssen einige Elemente des letzten (nicht vollständigen) Blockes in einen Überlaufpuffer verschoben werden.

Bevor mit der Blockpermutation begonnen wird, müssen eventuell einige leere Blöcke an das Ende ihres Buckets verschoben werden. Anschließend wird für jeden Bucket $b_{i}$ ein Schreibpointer $w_{i}$ angelegt, der auf den Beginn der Buckets zeigt und ein Lesepointer $r_{i}$, der auf den letzten nichtleeren Block im Bucket $b_{i}$ zeigt.

Um gegenseitige Blockierungen der Prozessoren zu vermeiden, verarbeitet jeder Prozessor anfangs einen eigenen (primären) Bucket $b_{\text{prim}}$. Außerdem erhält jeder Prozessor zwei Tauschpuffer, die jeweils genau einen Block speichern können. Falls beide Tauschpuffer leer sind, wird der Block an Stelle von $r_{\text{prim}}$ in einen der Tauschpuffer platziert. Durch Klassifikation des ersten Elements im Puffer wird der Zielbucket $b_{\text{dest}}$ dieses Blocks bestimmt und an Stelle $w_{\text{dest}}$ kopiert. Um zu verhindern, dass ein dortiger Block überschrieben wird, muss dieser vorher in einen freien Tauschpuffer kopiert werden. Mit diesem Block im Tauschpuffer wird dann analog verfahren. Die Pointer $r_{\text{prim}}$ und $w_{\text{dest}}$ werden entsprechend angepasst.

Sind alle Blöcke des primären Buckets an der richtigen Stelle, werden die Blöcke des nächsten Buckets weiterverarbeitet, bis alle Blöcke ihren Zielbucket erreicht haben. Danach endet diese Phase.

Aufräumen

Da nur ganze Blöcke in der Blockpermutationsphase verschoben wurden, können einige Elemente noch an Stelle außerhalb ihres Buckets liegen. Da es für alle Elemente noch genug Platz im ursprünglichen Array gibt, können die verbleibenden Elemente aus den Puffern an die freien Plätze kopiert werden. Zuletzt wird der Inhalt des Überlaufpuffers kopiert.

Literatur

• W. D. Frazer, A. C. McKellar: Samplesort: A Sampling Approach to Minimal Storage Tree Sorting. In: Journal of the ACM. Band 17, Nr. 3, Juli 1970, ISSN 0004-5411, S. 496–507. 
• Peter Sanders et al.: Vorlesung Parallele Algorithmen. Karlsruhe 2009, S. 75–80 (kit.edu [PDF; 2,1 MB]). 

Erweiterungen für parallele Berechnungen:

• ohio-state.edu (Memento vom 8. Juni 2011 im Internet Archive) (PDF)
• citeseer.ist.psu.edu
• citeseerx.ist.psu.edu
• Verwendung auf der GPU (PDF; 869 kB)