programming - julia nombre

Acelerando los ejemplos R mal escritos de Julia

Question

Los ejemplos de Julia para comparar el rendimiento con R parecen particularmente intrincados . https://github.com/JuliaLang/julia/blob/master/test/perf/perf.R

¿Cuál es el rendimiento más rápido que puede obtener de los dos algoritmos siguientes (preferiblemente con una explicación de lo que ha cambiado para que sea más parecido a un 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

Hmm, en el ejemplo de Mandelbrot, la matriz M tiene sus dimensiones transpuestas

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

porque se completa aumentando el count en el ciclo interno (valores sucesivos de im ). Mi implementación crea un vector de números complejos en mandelperf.1 y opera en todos los elementos, usando un índice y un subconjunto para hacer un seguimiento de los elementos del vector que aún no han satisfecho la condición 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))
}

para una aceleración de 13 veces (los resultados son iguales pero no idénticos porque el original arroja valores numéricos en vez de valores enteros).

> 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

El ejemplo de quicksort no ordena

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

pero mi principal aceleración fue vectorizar la partición alrededor del pivote

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

para una aceleración de 7 veces (en comparación con el original sin corregir)

> 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

Dado que en la comparación original, Julia es aproximadamente 30 veces más rápida que R para mandel, y 500 veces más rápida para quicksort, las implementaciones anteriores todavía no son realmente competitivas.

Answer 2

La palabra clave en esta pregunta es "algoritmo":

¿Cuál es el rendimiento más rápido que puede obtener de los dos algoritmos siguientes (preferiblemente con una explicación de lo que ha cambiado para que sea más parecido a un R)?

Como en "¿cuán rápido puedes hacer estos algoritmos en R?" Los algoritmos en cuestión aquí son el algoritmo de iteración de bucle complejo Mandelbrot estándar y el kernel quicksort estándar recursivo.

Sin duda, hay formas más rápidas de calcular las respuestas a los problemas planteados en estos puntos de referencia, pero no utilizando los mismos algoritmos. Puede evitar la recursión, evitar la iteración y evitar cualquier otra cosa en la que R no sea bueno. Pero ya no estás comparando los mismos algoritmos.

Si realmente desea calcular conjuntos de Mandelbrot en R o ordenar números, sí, así no es como escribiría el código. Podrías vectorizarlo tanto como sea posible, lo que empujaría todo el trabajo hacia los núcleos C predefinidos, o simplemente escribirías una extensión C personalizada y harías el cálculo allí. De cualquier manera, la conclusión es que R no es lo suficientemente rápido como para obtener un buen rendimiento por sí mismo: es necesario que C haga la mayor parte del trabajo para obtener un buen rendimiento.

Y ese es exactamente el punto de estos puntos de referencia: en Julia nunca tienes que depender del código C para obtener un buen rendimiento. Simplemente puede escribir lo que quiere hacer en Julia pura y tendrá un buen rendimiento. Si un algoritmo de bucle escalar iterativo es la forma más natural de hacer lo que desea hacer, simplemente hazlo. Si la recursión es la forma más natural de resolver el problema, está bien también. En ningún momento se verá obligado a confiar en C para el rendimiento, ya sea mediante una vectorización antinatural o escribiendo extensiones C personalizadas. Por supuesto, puedes escribir código vectorizado cuando es natural, como suele ser en álgebra lineal; y puede llamar a C si ya tiene una biblioteca que hace lo que quiere. Pero no es necesario.

Queremos tener la comparación más justa posible de los mismos algoritmos en todos los idiomas:

1. Si alguien tiene versiones más rápidas en R que usan el mismo algoritmo , ¡envíe los parches!
2. Creo que los puntos de referencia R en el sitio de julia ya están compilados en bytes, pero si lo estoy haciendo mal y la comparación es injusta con R, házmelo saber y lo arreglaré y actualizaré los puntos de referencia.