Beschleunigung von Julias schlecht geschriebenen R-Beispielen

75

Die Julia-Beispiele zum Vergleich der Leistung mit R scheinen besonders kompliziert zu sein . https://github.com/JuliaLang/julia/blob/master/test/perf/perf.R

Was ist die schnellste Leistung, die Sie mit den beiden folgenden Algorithmen erzielen können (vorzugsweise mit einer Erklärung dessen, was Sie geändert haben, um es R-ähnlicher zu machen)?

## mandel

mandel = function(z) {
    c = z
    maxiter = 80
    for (n in 1:maxiter) {
        if (Mod(z) > 2) return(n-1)
        z = z^2+c
    }
    return(maxiter)
}

mandelperf = function() {
    re = seq(-2,0.5,.1)
    im = seq(-1,1,.1)
    M = matrix(0.0,nrow=length(re),ncol=length(im))
    count = 1
    for (r in re) {
        for (i in im) {
            M[count] = mandel(complex(real=r,imag=i))
            count = count + 1
        }
    }
    return(M)
}

assert(sum(mandelperf()) == 14791)

## quicksort ##

qsort_kernel = function(a, lo, hi) {
    i = lo
    j = hi
    while (i < hi) {
        pivot = a[floor((lo+hi)/2)]
        while (i <= j) {
            while (a[i] < pivot) i = i + 1
            while (a[j] > pivot) j = j - 1
            if (i <= j) {
                t = a[i]
                a[i] = a[j]
                a[j] = t
            }
            i = i + 1;
            j = j - 1;
        }
        if (lo < j) qsort_kernel(a, lo, j)
        lo = i
        j = hi
    }
    return(a)
}

qsort = function(a) {
  return(qsort_kernel(a, 1, length(a)))
}

sortperf = function(n) {
    v = runif(n)
    return(qsort(v))
}

sortperf(5000)
Ari B. Friedman
quelle
3
Zunächst einmal rtricks.blogspot.ca/2007/04/…
Ben Bolker
10
Um Himmels willen ... lassen Sie R-Programmierer R programmieren
John
24
(1) Hier ist ein Beispiel für Fibonacci in R johnmyleswhite.com/notebook/2012/03/31/julia-i-love-you und es scheint, dass sie das verwenden, um zu dem Schluss zu kommen, dass Julia schneller war, aber meine Kommentare unter dem Blog-Beitrag überprüft haben Ich konnte die R-Lösung (immer noch mit nur reinem R) umschreiben und sie 2000x schneller laufen lassen. (2) Viele können dazu gebracht werden, 3x-4x schneller in R durch Byte-Kompilierung auszuführen, und das erfordert nicht einmal, dass Sie den Code ändern. (3) Viele der Beispiele sind von Anfang an gegen R gestapelt, da sie eine Rekursion verwenden, in der R nicht gut ist. Das Einbeziehen von Problemen in die Mischung, die leicht vektorisiert werden können, wäre fairer.
G. Grothendieck
7
@ G.Grothendieck Du solltest deinen Kommentar als Antwort Gabor posten; viele relevante Punkte dort. +1
Gavin Simpson
2
Es könnte interessant sein zu sehen, dass all dieses Benchmarking auch auf Radford Neals pqR ausgedehnt wird.
Spacedman

Antworten:

44

Hmm, im Mandelbrot-Beispiel sind die Abmessungen der Matrix M transponiert

M = matrix(0.0,nrow=length(im), ncol=length(re))

weil es durch Inkrementieren countin der inneren Schleife gefüllt wird (aufeinanderfolgende Werte von im). Meine Implementierung erstellt einen Vektor komplexer Zahlen in mandelperf.1und bearbeitet alle Elemente, wobei ein Index und eine Teilmenge verwendet werden, um zu verfolgen, welche Elemente des Vektors die Bedingung noch nicht erfüllt habenMod(z) <= 2

mandel.1 = function(z, maxiter=80L) {
    c <- z
    result <- integer(length(z))
    i <- seq_along(z)
    n <- 0L
    while (n < maxiter && length(z)) {
        j <- Mod(z) <= 2
        if (!all(j)) {
            result[i[!j]] <- n
            i <- i[j]
            z <- z[j]
            c <- c[j]
        }
        z <- z^2 + c
        n <- n + 1L
    }
    result[i] <- maxiter
    result
}

mandelperf.1 = function() {
    re = seq(-2,0.5,.1)
    im = seq(-1,1,.1)
    mandel.1(complex(real=rep(re, each=length(im)),
                     imaginary=im))
}

für eine 13-fache Beschleunigung (die Ergebnisse sind gleich, aber nicht identisch, da das Original eher numerische als ganzzahlige Werte zurückgibt).

> library(rbenchmark)
> benchmark(mandelperf(), mandelperf.1(),
+           columns=c("test", "elapsed", "relative"),
+           order="relative")
            test elapsed relative
2 mandelperf.1()   0.412  1.00000
1   mandelperf()   5.705 13.84709

> all.equal(sum(mandelperf()), sum(mandelperf.1()))
[1] TRUE

Das Quicksort-Beispiel sortiert nicht wirklich

> set.seed(123L); qsort(sample(5))
[1] 2 4 1 3 5

Meine Hauptbeschleunigung bestand jedoch darin, die Partition um den Drehpunkt herum zu vektorisieren

qsort_kernel.1 = function(a) {
    if (length(a) < 2L)
        return(a)
    pivot <- a[floor(length(a) / 2)]
    c(qsort_kernel.1(a[a < pivot]), a[a == pivot], qsort_kernel.1(a[a > pivot]))
}

qsort.1 = function(a) {
    qsort_kernel.1(a)
}

sortperf.1 = function(n) {
    v = runif(n)
    return(qsort.1(v))
}

für eine 7-fache Beschleunigung (im Vergleich zum unkorrigierten Original)

> benchmark(sortperf(5000), sortperf.1(5000),
+           columns=c("test", "elapsed", "relative"),
+           order="relative")
              test elapsed relative
2 sortperf.1(5000)    6.60 1.000000
1   sortperf(5000)   47.73 7.231818

Da Julia im ursprünglichen Vergleich etwa 30-mal schneller als R für Mandel und 500-mal schneller für Quicksort ist, sind die obigen Implementierungen immer noch nicht wirklich wettbewerbsfähig.

Martin Morgan
quelle
@ gsk3 Byte-Kompilierung ändert die Ausführungszeit von Mandel nicht wirklich. sortperf.1 wird ungefähr 10x statt 7x schneller als sortperf.
Martin Morgan
Das Problem ist, dass qsort_kernel.1immer noch eine Rekursion durchgeführt wird und R darin nicht gut ist. Um es schneller laufen zu lassen, konvertieren Sie es mithilfe eines Stapels in Schleifen.
G. Grothendieck
@ G.Grothendieck die Rekursion ist nicht so tief (log N), iterative Lösungen verwenden wahrscheinlich einen Stapel (auch nicht R-freundlich), der Geschwindigkeitsunterschied ist immer noch groß und der Algorithmus weniger intuitiv.
Martin Morgan
1
Die Verwendung einer guten Mathematikbibliothek in R wie Intel MKL hilft auch, die Lücke zwischen R und Julia zu verringern.
Aatrujillob
Bei den Leistungsvergleichen geht es nicht darum, wie schnell eine Lösung für ein Problem von verschiedenen Sprachen erreicht werden kann, sondern darum, wie dieselbe Implementierung (dh ein Algorithmus) in verschiedenen Sprachen funktioniert.
Bdeonovic
101

Das Schlüsselwort in dieser Frage ist "Algorithmus":

Was ist die schnellste Leistung, die Sie aus den beiden herausholen können? Algorithmen (vorzugsweise mit einer Erklärung dessen, was Sie geändert haben, um es R-ähnlicher zu machen)?

Wie in "Wie schnell können Sie diese machen Algorithmen in R erstellen?" Die hier fraglichen Algorithmen sind der Standard-Mandelbrot-Algorithmus für komplexe Schleifeniterationen und der Standard-rekursive Quicksort-Kernel.

Es gibt sicherlich schnellere Möglichkeiten, die Antworten auf die in diesen Benchmarks aufgeworfenen Probleme zu berechnen - jedoch nicht mit denselben Algorithmen. Sie können Rekursion vermeiden, Iteration vermeiden und alles vermeiden, was R sonst nicht kann. Aber dann vergleichen Sie nicht mehr dieselben Algorithmen.

Wenn Sie wirklich Mandelbrot-Mengen in R berechnen oder Zahlen sortieren wollten, ja, so würden Sie den Code nicht schreiben. Sie würden es entweder so weit wie möglich vektorisieren und so die gesamte Arbeit in vordefinierte C-Kernel verschieben - oder einfach eine benutzerdefinierte C-Erweiterung schreiben und dort die Berechnung durchführen. In beiden Fällen ist die Schlussfolgerung, dass R nicht schnell genug ist, um selbst eine wirklich gute Leistung zu erzielen. Sie müssen C den größten Teil der Arbeit erledigen lassen, um eine gute Leistung zu erzielen.

Und genau darum geht es bei diesen Benchmarks: In Julia muss man sich nie auf C-Code verlassen, um eine gute Leistung zu erzielen. Sie können einfach schreiben, was Sie in reiner Julia tun möchten, und es wird eine gute Leistung haben. Wenn ein iterativer Skalarschleifenalgorithmus der natürlichste Weg ist, das zu tun, was Sie tun möchten, dann tun Sie das einfach. Wenn Rekursion der natürlichste Weg ist, um das Problem zu lösen, ist das auch in Ordnung. Sie werden zu keinem Zeitpunkt gezwungen sein, sich für die Leistung auf C zu verlassen - sei es durch unnatürliche Vektorisierung oder durch das Schreiben benutzerdefinierter C-Erweiterungen. Natürlich Sie können Code vektorisiert schreiben , wenn es natürlich ist, wie es oft in der linearen Algebra ist; und Sie können C aufrufen, wenn Sie bereits eine Bibliothek haben, die das tut, was Sie wollen. Aber das musst du nicht.

Wir wollen einen möglichst fairen Vergleich derselben Algorithmen zwischen verschiedenen Sprachen:

  1. Wenn jemand schnellere Versionen in R hat, die denselben Algorithmus verwenden , senden Sie bitte Patches!
  2. Ich glaube, dass die R-Benchmarks auf der Julia-Site bereits bytekompiliert sind, aber wenn ich es falsch mache und der Vergleich mit R unfair ist, lassen Sie es mich bitte wissen und ich werde es reparieren und die Benchmarks aktualisieren.
StefanKarpinski
quelle
29
@Stefan Gute Punkte. Der Kontrapunkt ist jedoch, dass der Beispielcode einfach ein schlechter R-Code ist, wenn Sie mehrere hundert- bis mehrere tausendfache Beschleunigungen erzielen können, indem Sie den Code einfach auf eine in einer Sprache natürliche Weise schreiben. Wenn jemand zu SO kam und diese Beispielcodes veröffentlichte, erhielt er sehr schnell Lektionen zum Schreiben von R, wie R geschrieben werden soll. Angesichts der Tatsache, dass die Beispiele alle so ausgewählt zu sein scheinen, dass sie eine Rekursion beinhalten (bei der R zugegebenermaßen schlecht ist), und dann der Beispielcode alles daran setzt, eine Vektorisierung zu vermeiden (bei der R sehr gut ist) ....
Ari B. Friedman
39
Ich würde immer noch argumentieren, dass am Benchmark-Code nichts falsch ist. Wenn es eine R-Implementierung gäbe, bei der Iteration und Rekursion schnell wären, wäre der Code dann immer noch schlecht? Die einzige Schlussfolgerung, die ich ziehen kann, ist, dass die Implementierung der Sprache fehlerhaft ist, nicht der Code. Wenn das Design der R-Sprache es außerdem besonders schwierig (oder vielleicht unmöglich?) Macht, Iteration und Rekursion schnell durchzuführen, würde ich sagen, dass dies weder ein Fehler in diesem bestimmten Code noch in der aktuellen R-Implementierung ist, sondern vielmehr in der Sprachdesign selbst.
StefanKarpinski
22
@Stefan Ich stimme zu, dass die Verwendung des gleichen Algorithmus ein fairer Vergleich ist, aber vielleicht lohnt es sich, die Julia-optimierten Versionen dieser Algorithmen genauso oder so nahe wie möglich an die R-optimierten Versionen zu schreiben Algorithmen, die erscheinen (dh stark vektorisiert sind) und sowohl die relative Leistung als auch die Leistung gegenüber den ursprünglichen Implementierungen des Benchmark-Codes in Julia neu bewerten? Wenn ich letztendlich dasselbe Endergebnis mit unterschiedlichem Code erzielen kann, der für eine andere Sprache optimiert ist, die schneller ist, werde ich wahrscheinlich den schnelleren Code verwenden. Ich folge mit großem Interesse!
Simon O'Hanlon
12
Wenn es möglich wäre, ein Paket zu schreiben, das die Rekursion oder Iteration in R schnell macht, hätte es jetzt nicht jemand getan? Grundlegende Sprachverbesserungen wie diese können Sie nicht "angehen" - sie erfordern tiefgreifende, schwierige Änderungen. Das bedeutet nicht, dass Rekursion und Iteration in R nicht schneller gemacht werden könnten, aber es wäre nicht "nur ein Paket".
StefanKarpinski
29
@VB: Das hängt alles davon ab, was Sie messen wollen. Persönlich interessiert mich überhaupt nicht, wie schnell man Fibonacci-Zahlen berechnen kann. Dies ist jedoch einer unserer Maßstäbe. Warum? Weil ich sehr daran interessiert bin, wie gut Sprachen die Rekursion unterstützen - und der doppelt rekursive Algorithmus zufällig ein großartiger Rekursionstest ist, gerade weil es eine so schreckliche Art ist, Fibonacci-Zahlen zu berechnen. Was würde man also lernen, wenn man einen absichtlich langsamen, übermäßig rekursiven Algorithmus in C und Julia mit einem kniffligen, cleveren, vektorisierten Algorithmus in R vergleicht? Gar nichts.
StefanKarpinski