Ari*_*man 72 optimization r julia
朱莉娅将性能与R进行比较的例子似乎特别令人费解. https://github.com/JuliaLang/julia/blob/master/test/perf/perf.R
你可以从下面的两种算法中获得最快的性能(最好解释你改变了什么以使它更像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)
Ste*_*ski 100
这个问题的关键词是"算法":
你可以从下面的两种算法中获得最快的性能(最好解释你改变了什么以使它更像R)？
就像"你能用多长时间在R中制作这些算法？" 这里讨论的算法是标准的Mandelbrot复杂循环迭代算法和标准的递归快速排序内核.
肯定有更快的方法来计算这些基准测试中提出的问题的答案 - 但不使用相同的算法.你可以避免递归,避免迭代,并避免R不擅长的任何东西.但是你不再比较相同的算法了.
如果你真的想在R中计算Mandelbrot集或排序数,是的,这不是你编写代码的方式.您可以尽可能地将其向量化 - 从而将所有工作推送到预定义的C内核中 - 或者只编写自定义C扩展并在那里进行计算.无论哪种方式,结论是R不够快,不能自己获得真正好的性能 - 你需要C做大部分工作才能获得良好的性能.
这正是这些基准测试的重点:在Julia中,您永远不必依赖C代码来获得良好的性能.你可以在纯粹的Julia中写下你想要做的事情,它会有很好的表现.如果迭代标量循环算法是您做自己想做的最自然的方式,那么就这样做.如果递归是解决问题的最自然方式,那也没关系.在任何时候你都不会被迫依靠C来获得性能 - 无论是通过不自然的矢量化还是编写自定义C扩展.当然,你可以在它自然的时候编写矢量化代码,因为它通常是线性代数; 如果您已经有一些可以满足您需求的库,那么您可以调用C. 但你没必要.
我们确实想要跨语言对相同算法进行最公平的比较:
Mar*_*gan 42
嗯,在Mandelbrot例子中,矩阵M的尺寸被转置
M = matrix(0.0,nrow=length(im), ncol=length(re))
因为它是通过count在内循环中递增来填充的(连续的值im).我的实现mandelperf.1在所有元素中创建了一个复数向量,并使用索引和子集来跟踪向量的哪些元素尚未满足条件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))
}
加速13倍(结果相同但不相同,因为原始值返回数值而不是整数值).
> 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
快速排序示例实际上没有排序
> set.seed(123L); qsort(sample(5))
[1] 2 4 1 3 5
但我的主要加速是围绕枢轴矢量化分区
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))
}
加速7倍(与未经修正的原件相比)
> 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
因为在最初的比较中,Julia比R和Mandel快30倍,而快速排序快500倍,上面的实现仍然没有真正的竞争力.