Accelerazione degli esempi R mal scritti di Julia

Question

Gli esempi Julia per confrontare le prestazioni con R seem particularly convoluted. https://github.com/JuliaLang/julia/blob/master/test/perf/perf.R

Qual è la prestazione più rapida che è possibile emettere dai due algoritmi seguenti (preferibilmente con una spiegazione di cosa è stato modificato per renderlo più simile a 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

Hm, nell'esempio Mandelbrot matrice M ha dimensioni trasposte

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

perché è riempito incrementando count nel ciclo interno (valori successivi di im). Mia implementazione crea un vettore di numeri complessi in mandelperf.1 e opera su tutti gli elementi, utilizzando un indice e sottoinsiemi per tenere traccia di quali elementi del vettore non hanno ancora soddisfatto la condizione 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)) 
} 

per una velocità di 13 volte -up (i risultati sono uguali ma non identici perché l'originale restituisce valori numerici anziché interi).

> 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'esempio quicksort in realtà non sorta

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

ma mio accelerazione principale era vectorize partizione intorno al perno

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

per un aumento di velocità di 7 volte (in confronto all'originale non corretto)

> 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 

Poiché nel confronto originale Julia è circa 30 volte più veloce di R per mandel e 500 volte più veloce per quicksort, le implementazioni sopra riportate non sono ancora molto competitive.

Answer 2

La parola chiave in questa domanda è "algoritmo":

Qual è la performance più veloce è possibile tirare avanti dei due algoritmi inferiori (preferibilmente con una spiegazione di ciò che è stato modificato per renderlo più R-simili)?

Come in "Quanto velocemente si può fare questi algoritmi in R?" Gli algoritmi in questione sono l'algoritmo standard di iterazione del ciclo del complesso di Mandelbrot e il kernel standard quicksort ricorsivo.

Ci sono sicuramente modi più veloci per calcolare le risposte ai problemi posti in questi benchmark - ma non utilizzare gli stessi algoritmi. Puoi evitare la ricorsione, evitare l'iterazione ed evitare qualsiasi altra cosa in cui R non sia bravo. Ma poi non stai più confrontando gli stessi algoritmi.

Se si desidera calcolare i set Mandelbrot in R o ordinare i numeri, sì, non è come si scrive il codice. Lo si dovrebbe vettorializzare il più possibile, spingendo così tutto il lavoro in kernel C predefiniti, o semplicemente si scrive un'estensione C personalizzata e si esegue il calcolo lì. In ogni caso, la conclusione è che R non è abbastanza veloce per ottenere prestazioni davvero buone da solo - è necessario che C svolga la maggior parte del lavoro per ottenere buone prestazioni.

E questo è esattamente il punto di questi benchmark: in Julia non devi mai fare affidamento sul codice C per ottenere buone prestazioni. Puoi semplicemente scrivere ciò che vuoi fare in pura Julia e avrà buone prestazioni. Se un algoritmo iterativo di loop scalare è il modo più naturale per fare ciò che vuoi fare, fallo semplicemente. Se la ricorsione è il modo più naturale per risolvere il problema, allora va bene anche quello. In nessun caso sarai costretto a fare affidamento su C per le prestazioni, sia tramite vettorializzazione innaturale che scrivendo estensioni C personalizzate. Certo, è possibile scrivere codice vettoriale quando è naturale, come spesso accade nell'algebra lineare; e tu puoi chiamare se hai già qualche libreria che fa quello che vuoi. Ma tu non devi.

Noi vogliamo avere la più bella possibile confronto tra gli stessi algoritmi tra le varie lingue:

1. Se qualcuno ha le versioni più veloci in R che utilizzare lo stesso algoritmo, si prega di inviare patch!
2. Credo che i benchmark R su julia site siano già compilati in byte, ma se sto sbagliando e il confronto non è corretto con R, per favore fatemelo sapere e lo aggiusterò e aggiornerò i benchmark.