什么时候矢量化在朱莉娅受青睐？

Question

我有两个函数用于在Julia中以数字方式确定pi.第二个函数(我认为是矢量化的)比第一个函数慢.为什么矢量化速度较慢？是否有规则何时进行矢量化以及何时不进行？

function determine_pi(n)
    area = zeros(Float64, n);
    sum = 0;
    for i=1:n
        if ((rand()^2+rand()^2) <=1)
            sum = sum + 1;
        end
            area[i] = sum*1.0/i;
    end
    return area
end

和另一个功能

function determine_pi_vec(n)
    res = cumsum(map(x -> x<=1?1:0, rand(n).^2+rand(n).^2))./[1:n]
    return res
end

当运行n = 10 ^ 7时,以下是执行时间(运行几次后)

n=10^7
@time returnArray = determine_pi(n)
#output elapsed time: 0.183211324 seconds (80000128 bytes allocated)
@time returnArray2 = determine_pi_vec(n);
#elapsed time: 2.436501454 seconds (880001336 bytes allocated, 30.71% gc time)

Answer 1

如果,矢量化是好的

• 它使代码更易于阅读,性能并不重要
• 如果它是线性代数运算,使用矢量化样式可能会很好,因为Julia可以使用BLAS和LAPACK来使用非常专业的高性能代码执行操作.

一般来说,我个人认为最好从矢量化代码开始,寻找任何速度问题,然后发现任何棘手的问题.

你的第二个代码很慢,因为它被矢量化,但是由于使用了匿名函数:不幸的是在Julia 0.3中,这些代码通常要慢得多.map一般来说表现不是很好,我相信因为Julia无法推断出函数的输出类型(从map函数的角度看它仍然是"匿名的" ).我写了一个不同的矢量化版本,它避免了匿名函数,并且可能更容易阅读:

function determine_pi_vec2(n)
    return cumsum((rand(n).^2 .+ rand(n).^2) .<= 1) ./ (1:n)
end

基准测试

function bench(n, f)
  f(10)
  srand(1000)
  @time f(n)
  srand(1000)
  @time f(n)
  srand(1000)
  @time f(n)
end

bench(10^8, determine_pi)
gc()
bench(10^8, determine_pi_vec)
gc()
bench(10^8, determine_pi_vec2)

给我结果

elapsed time: 5.996090409 seconds (800000064 bytes allocated)
elapsed time: 6.028323688 seconds (800000064 bytes allocated)
elapsed time: 6.172004807 seconds (800000064 bytes allocated)
elapsed time: 14.09414031 seconds (8800005224 bytes allocated, 7.69% gc time)
elapsed time: 14.323797823 seconds (8800001272 bytes allocated, 8.61% gc time)
elapsed time: 14.048216404 seconds (8800001272 bytes allocated, 8.46% gc time)
elapsed time: 8.906563284 seconds (5612510776 bytes allocated, 3.21% gc time)
elapsed time: 8.939001114 seconds (5612506184 bytes allocated, 4.25% gc time)
elapsed time: 9.028656043 seconds (5612506184 bytes allocated, 4.23% gc time)

因此,在某些情况下,矢量化代码绝对可以与开发代码一样好,即使我们不在线性代数情况下也是如此.