朱莉娅语言。如何击败矢量化操作？

Question

据说 Julia for 循环与向量化操作一样快，甚至更快（如果使用得当的话）。 我有两段代码。这个想法是为给定的0-1序列找到一个样本统计量，即x（在这两个例子中，我试图找到一个总和，但还有更复杂的例子，我只是想理解一个一般含义我的代码中的性能缺陷）。第一个看起来像：

S = 2 * sum(x) - n
s_obs_new = abs(S) / sqrt(2 * n)
pval = erfc(s_obs_new)

第二个是“天真”和经典的东西：

S = 0
for i in eachindex(x)
    S += x[i]
end
S = 2 * S - n
s_obs_new = abs(S) / sqrt(2 * n)
pval = erfc(s_obs_new)

使用@benchmark 我发现第一个示例的运行时间约为 11.8 毫秒，第二个示例的运行时间为 38 毫秒。

这个例子对我来说非常重要，因为在很多其他地方，矢量化是不可能的，所以我想以去矢量化的“方式”和矢量化一样快地进行计算。

有什么想法为什么去向量化的代码可能比向量化的代码慢 4 倍？类型稳定性还可以，没有不必要的大内存分配等。

第一个函数的代码是：

function frequency_monobit_test1( x :: Array{Int8, 1}, n = 0)
# count 1 and 0 in sequence, we want the number
# of 1's and 0's to be approximately the same
# reccomendation n >= 100
# decision Rule(at 1% level): if pval < 0.01 -> non-random
if (n == 0)
    n = length(x)
end
S = 2 * sum(x) - n
s_obs_new = abs(S) / sqrt(2 * n)
pval = erfc(s_obs_new)
return pval  

第二个是：

function frequency_monobit_test2( x :: Array{Int8, 1}, n = 0)
# count 1 and 0 in sequence, we want the number
# of 1's and 0's to be approximately the same
# reccomendation n >= 100
# decision Rule(at 1% level): if pval < 0.01 -> non-random
if (n == 0)
    n = length(x)
end
S = 0
@inbounds for i in eachindex(x)
    S += x[i]
end
S = 2 * S - n
s_obs_new = abs(S) / sqrt(2 * n)
pval = erfc(s_obs_new)
return pval

Answer 1

这是一个奇怪的案例。在Int64 变量中累积Int8s 时似乎存在性能问题。

让我们试试这些功能：

using SpecialFunctions, BenchmarkTools

function frequency_monobit_test1(x, n=length(x))
    S = sum(x)
    return erfc(abs(2S - n) / sqrt(2n))
end

function frequency_monobit_test3(typ::Type{<:Integer}, x, n=length(x))
    S = zero(typ)
    for i in eachindex(x)
        @inbounds S += x[i]
    end
    return erfc(abs(2S - n) / sqrt(2n))
end

初始化一些向量

N = 2^25;
x64 = rand(0:1, N);
x8 = rand(Int8[0, 1], N);
xB = rand(Bool, N);
xb = bitrand(N);

基准测试：

对于Int64 输入：

julia> @btime frequency_monobit_test1($x64)
  17.540 ms (0 allocations: 0 bytes)
0.10302739916042186

julia> @btime frequency_monobit_test3(Int64, $x64)
  17.796 ms (0 allocations: 0 bytes)
0.10302739916042186

julia> @btime frequency_monobit_test3(Int32, $x64)
  892.715 ms (67106751 allocations: 1023.97 MiB)
0.10302739916042186

我们看到sum 和显式循环同样快，使用Int32 进行初始化是个坏主意。

对于Int32 输入：

julia> @btime frequency_monobit_test1($x32)
  9.137 ms (0 allocations: 0 bytes)
0.2386386867682374

julia> @btime frequency_monobit_test3(Int64, $x32)
  8.839 ms (0 allocations: 0 bytes)
0.2386386867682374

julia> @btime frequency_monobit_test3(Int32, $x32)
  7.274 ms (0 allocations: 0 bytes)
0.2386386867682374

sum 和循环的速度相似。积累到Int32 可以节省一点时间。

Int8 输入：

julia> @btime frequency_monobit_test1($x8)
  5.681 ms (0 allocations: 0 bytes)
0.16482999123032094

julia> @btime frequency_monobit_test3(Int64, $x8)
  19.517 ms (0 allocations: 0 bytes)
0.16482999123032094

julia> @btime frequency_monobit_test3(Int32, $x8)
  4.815 ms (0 allocations: 0 bytes)
0.16482999123032094

显式循环在累积到Int32 时稍微快一点，但是天哪！ Int64 是怎么回事？太慢了！

Bool 呢？

julia> @btime frequency_monobit_test1($xB)
  9.627 ms (0 allocations: 0 bytes)
0.7728544347518309

julia> @btime frequency_monobit_test3(Int64, $xB)
  9.629 ms (0 allocations: 0 bytes)
0.7728544347518309

julia> @btime frequency_monobit_test3(Int32, $xB)
  4.815 ms (0 allocations: 0 bytes)
0.7728544347518309

Loop 和sum 速度相同，但是累加到Int32 可以节省一半的时间。

现在我们将尝试BitArray：

julia> @btime frequency_monobit_test1($xb)
  259.044 μs (0 allocations: 0 bytes)
0.7002576522570715

julia> @btime frequency_monobit_test3(Int64, $xb)
  19.423 ms (0 allocations: 0 bytes)
0.7002576522570715

julia> @btime frequency_monobit_test3(Int32, $xb)
  19.430 ms (0 allocations: 0 bytes)
0.7002576522570715

所以，BitArray 上的 sum 非常快，因为您可以执行分块加法，但在循环中提取单个元素会产生一些开销。

结论：

• 使用循环可以获得与使用 sum 相似或更好的性能，但 BitArrays 除外，这是一种非常特殊的情况。
• 如果您知道数组的长度，并且知道 Int32 足以容纳您的总和，那么这可以节省时间。
• 将Int8s 累积到Int64 时发生了一些非常奇怪的事情。我不知道为什么性能这么差。
• 如果您只对零和一感兴趣，请使用Bools 的数组，而不是Int8s 的数组，并可能累积到Int32。
• BitArrays 在某些情况下可以超快。
• sum 在Int8s 向量上的速度出奇的快，是sum(::Vector{Bool}) 的两倍。