Accélérer les exemples R mal écrits de Julia

Question

Les exemples de Julia pour comparer les performances avec R semble particulièrement alambiqué . https://github.com/JuliaLang/julia/blob/master/test/perf/perf.R

Quelle est la performance la plus rapide que vous pouvez obtenir des deux algorithmes ci-dessous (de préférence avec une explication de ce que vous avez changé pour le rendre plus proche de R) ?

## 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)

Answer 1

Le mot clé de cette question est "algorithme" :

Quelle est la performance la plus rapide que vous pouvez obtenir des deux algorithmes ci-dessous (de préférence avec une explication de ce que vous avez changé pour le rendre plus proche de R) ?

Comme dans "à quelle vitesse pouvez-vous faire ces algorithmes en R ?" Les algorithmes en question ici sont l'algorithme standard d'itération en boucle complexe de Mandelbrot et le noyau récursif standard de quicksort.

Il existe certainement des moyens plus rapides de calculer les réponses aux problèmes posés dans ces repères, mais pas en utilisant les mêmes algorithmes. Vous pouvez éviter la récursion, l'itération et tout ce que R ne sait pas faire. Mais alors vous ne comparez plus les mêmes algorithmes.

Si vous vouliez vraiment calculer des ensembles de Mandelbrot dans R ou trier des nombres, oui, ce n'est pas ainsi que vous écririez le code. Il faudrait soit le vectoriser autant que possible - en poussant ainsi tout le travail dans des noyaux C prédéfinis - soit écrire une extension C personnalisée et y effectuer les calculs. Quoi qu'il en soit, la conclusion est que R n'est pas assez rapide pour obtenir de bonnes performances par lui-même - il faut que le C fasse la majeure partie du travail pour obtenir de bonnes performances.

Et c'est exactement le but de ces benchmarks : dans Julia, vous n'avez jamais à vous fier au code C pour obtenir de bonnes performances. Vous pouvez juste écrire ce que vous voulez faire en pure Julia et cela aura de bonnes performances. Si un algorithme de boucle scalaire itérative est la manière la plus naturelle de faire ce que vous voulez faire, alors faites-le. Si la récursion est la façon la plus naturelle de résoudre le problème, alors c'est ok aussi. À aucun moment vous ne serez forcé de compter sur le C pour les performances - que ce soit via une vectorisation non naturelle ou en écrivant des extensions C personnalisées. Bien entendu, vous peut écrire du code vectoriel lorsque c'est naturel, comme c'est souvent le cas en algèbre linéaire ; et vous peut appeler C si vous avez déjà une bibliothèque qui fait ce que vous voulez. Mais vous n'êtes pas obligé de le faire.

Nous souhaitons que la comparaison des mêmes algorithmes entre les langues soit la plus équitable possible :

1. Si quelqu'un a des versions plus rapides dans R que utilisent le même algorithme veuillez soumettre des correctifs !
2. Je crois que les repères R sur le site de julia sont déjà compilés en octets, mais si je m'y prends mal et que la comparaison est injuste pour R, merci de me le faire savoir et je corrigerai cela et mettrai à jour les benchmarks.

Answer 2

Hmm, dans l'exemple de Mandelbrot, la matrice M a ses dimensions transposées.

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

car il est rempli en incrémentant count dans la boucle interne (valeurs successives de im ). Mon implémentation crée un vecteur de nombres complexes dans mandelperf.1 et opère sur tous les éléments, en utilisant un index et un sous-ensemble pour garder la trace des éléments du vecteur qui n'ont pas encore satisfait à la condition Mod(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))
}

pour une accélération de 13 fois (les résultats sont égaux mais pas identiques car l'original renvoie des valeurs numériques plutôt que des valeurs entières).

> 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

L'exemple de tri rapide ne trie pas réellement.

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

mais ma principale accélération a été de vectoriser la partition autour du pivot

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))
}

pour une accélération de 7 fois (par rapport à l'original non corrigé)

> 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

Puisque dans la comparaison originale Julia est environ 30 fois plus rapide que R pour mandel, et 500 fois plus rapide pour quicksort, les implémentations ci-dessus ne sont toujours pas vraiment compétitives.