Schnellere Zufallsstichproben mit dqrng

Vor etwa einem halben Jahr beschrieb ich dqsample als eine faire Alternative zu base::sample(). Inzwischen ist das nur mehr von historischem Interesse, nachdem die R Version 3.6.0 einen aktualisierten Algorithmus für Zufallsstichproben enthalten wird. Einige der Techniken aus dqsample sind nun aber Teil von dqrng. Selbstverständlich ist der verwendete Algorithmus weiterhin fair:

Diese Version von dqrng ist noch nicht auf CRAN verfügbar, kann aber über drat installiert werden:

Insgesamt werden wir sehen, dass es möglich ist schneller als R Zufallsstichproben zu erstellen ohne Kompromisse bei der statistischen Qualität eingehen zu müssen. Feedback ist immer willkommen!

Benchmarks

Die folgenden Benchmarks wurden mit R Version 3.5.3 (2019-03-11) unter Debian GNU/Linux 9 (stretch) auf einem Intel® Core™ i7-6600U CPU @ 2.60GHz erstellt. YMMV, wie immer bei Benchmarks …

Durch die Kombination schneller Zufallszahlgeneratoren mit schnellen Methoden zur Erzeugung ganzer Zahlen erhält man gute Performance beim Ziehen mit Zurücklegen:

expression min mean median max itr/sec
sample.int(m, n, replace = TRUE) 94.9µs 112.5µs 112.2µs 3.78ms 8886.863
sample.int(10000 * m, n, replace = TRUE) 170.9µs 202.8µs 200.9µs 3.08ms 4930.094
dqsample.int(m, n, replace = TRUE) 31.5µs 36.2µs 35.7µs 183.63µs 27596.508
dqsample.int(10000 * m, n, replace = TRUE) 35.7µs 41.3µs 39.8µs 3.06ms 24186.338

Beachte dass für 10^10 Zahlen der “long-vector support” in R benötigt wird.

Beim Ziehen ohne Zurücklegen muss man sich Gedanken machen wie man sicher stellt, dass kein Eintrag sich wiederholt. Wenn mehr als 50% der Gesamtpopulation gezogen werden sollen, so mischt dqrng die Gesamtpopulation bis zur benötigten Größe. R nutzen eine sehr ähnlichen Algorithmus, jedoch ist dqrng auch hier schneller:

expression min mean median max itr/sec
sample.int(m, n) 11.91ms 14.09ms 13.02ms 26.11ms 70.98544
dqsample.int(m, n) 6.93ms 8.22ms 8.17ms 9.99ms 121.65317

Wird ein geringerer Teil der Population gezogen, so verwendet dqrng die Verwerfungsmethode mit einem Set als Speicher der bereits gezogenen Element. Grundsätzlich kann R einen ganz ähnlichen Algorithmus auf Basis eines Hashsets verwenden, tut dies standardmäßig aber erst für größere Populationen, obwohl es schon früher schneller als die Standardmethode ist. Der in dqrng genutzte und auf einem Bitset beruhende Algorithmus ist dabei noch schneller:

expression min mean median max itr/sec
sample.int(m, n) 788µs 1.72ms 1.87ms 5.01ms 583.0849
sample.int(m, n, useHash = TRUE) 229µs 269.6µs 266.15µs 453.71µs 3709.2196
dqsample.int(m, n) 113µs 130.4µs 129.49µs 2.98ms 7668.6069

Wird der Stichprobenanteil weiter reduziert, so verwendet auch dqrng ein Hashset. Die beiden Hashset basierten Algorithmen sind vergleichbar schnell und deutlich schneller als die Standardmethode in R:

expression min mean median max itr/sec
sample.int(m, n) 452.72µs 1.29ms 1.56ms 4.69ms 777.0096
sample.int(m, n, useHash = TRUE) 3.98µs 5.38µs 5.03µs 62.94µs 185978.0662
dqsample.int(m, n) 3.59µs 4.14µs 4.05µs 26.85µs 241454.9528

Für größere Gesamtpopulationen verwendet R standardmäßig ein Hashset. Jedoch ist die Version in dqrng schneller:

expression min mean median max itr/sec
sample.int(m, n) 5.3ms 5.79ms 5.48ms 9.5ms 172.7915
dqsample.int(m, n) 1.7ms 2.04ms 1.97ms 5.08ms 490.8562

Details

Die folgenden Methoden werden beim Ziehen ohne Zurücklegen verwendet. Die Algorithmen werden in R-artigem Pseudocode dargestellt, obwohl die tatsächliche Umsetzung in C++ erfolgte. Für Stichprobenanteile über 50% wird partiell nach Fisher-Yates gemischt:

wobei random_int(m-i) eine Zufallszahl aus [0, m-i] liefert. Nachdem hier die Gesamtpopulation im Speicher gehalten werden muss, ist diese Methode nur für hohe Stichprobenanteile geeignet. Man könnte erwarten, dass Reservoir Sampling bei geringeren Anteilen gut funktioniert. In meinen Tests waren jedoch Set-basierte Algorithmen schneller:

Dabei gibt elems.insert(v) den Wert TRUE zurück, wenn v nicht Teil von elems war und daher eingefügt wurde, andernfalls wird FALSE zurückgegeben. Es gibt verschiedene Möglichkeiten solch ein Set umzusetzen. Für mittlere Stichprobenanteile (derzeit zwischen 0.1% und 50%) verwendet dqrng ein Bitset, d.h. einen Vektor aus m Bits von denen jeder einen der möglichen Werte repräsentiert. Für niedrigere Stichprobenanteile wird der Speicherbedarf des Bitsets zum Problem, weshalb dann ein Hashset zum Einsatz kommt.1 Hier skaliert der Speicherbedarf nur mit n und nicht mit m. Man könnte erwarten, dass der auf Robert Floyd zurückgehende Algorithmus besser sein sollte, aber das war in meinen Tests nicht der Fall. Vermutlich weil das Ergebnis nochmals gemischt werden muss um eine zufällige _Permutation_ und nicht nur eine zufällige _Kombination_ zu erhalten.

  1. Für die Spezialisten: Open addressing mit power-of-two Größe zwischen 1.5 und 3 mal n, identity hash und quadratic probing.