Julia：向量化代码与去向量化代码

Question

据我了解，Julia 应该使 for 循环更快，并且与矢量化操作一样快。我编写了一个简单函数的三个版本，该函数使用 for 循环查找距离，使用矢量化操作，使用 DataFrames：

x = rand(500)
y = rand(500)
a = rand()
b = rand()

function devect()
    dist = Array(Float64, 0)
    twins = Array(Float64, 0,2)

    for i in 1:500
        dist = [dist; sqrt((x[i] - a)^2 + (y[i] - b)^2)]
        if dist[end] < 0.05
            twins = [twins; [x y][end,:]]
        end
    end

    return twins
end

function vect()
    d = sqrt((x-a).^2 + (y-b).^2)
    return [x y][d .< 0.05,:]
end

using DataFrames

function df_vect()
    df = DataFrame(x=x, y=y)
    dist = sqrt((df[:x]-a).^2 + (df[:y]-b).^2)

    return df[dist .< 0.05,:]
end

n = 10^3

@time for i in [1:n] devect() end
@time for i in [1:n] vect() end
@time for i in [1:n] df_vect() end

输出：

elapsed time: 4.308049576 seconds (1977455752 bytes allocated, 24.77% gc time)
elapsed time: 0.046759167 seconds (37295768 bytes allocated, 54.36% gc time)
elapsed time: 0.052463997 seconds (30359752 bytes allocated, 49.44% gc time)

为什么向量化版本的执行速度如此之快？

Answer 1

跟进我对用于在 devect 中构建解决方案的方法的评论。这是我的代码

julia> x, y, a, b = rand(500), rand(500), rand(), rand()

julia> function devect{T}(x::Vector{T}, y::Vector{T}, a::T, b::T)
       res = Array(T, 0)
       dim1 = 0
       for i = 1:size(x,1)
           if sqrt((x[i]-a)^2+(y[i]-b)^2) < 0.05
               push!(res, x[i])
               push!(res, y[i])
               dim1 += 1
           end
       end
       reshape(res, (2, dim1))'
   end
devect (generic function with 1 method)

julia> function vect{T}(x::Vector{T}, y::Vector{T}, a::T, b::T)
       d = sqrt((x-a).^2+(y-b).^2)
       [x y][d.<0.05, :]
   end
vect (generic function with 1 method)

julia> @time vect(x, y, a, b)
elapsed time: 3.7118e-5 seconds (37216 bytes allocated)
2x2 Array{Float64,2}:
 0.978099  0.0405639
 0.94757   0.0224974

julia> @time vect(x, y, a, b)
elapsed time: 7.1977e-5 seconds (37216 bytes allocated)
2x2 Array{Float64,2}:
 0.978099  0.0405639
 0.94757   0.0224974

julia> @time devect(x, y, a, b)
elapsed time: 1.7146e-5 seconds (376 bytes allocated)
2x2 Array{Float64,2}:
 0.978099  0.0405639
 0.94757   0.0224974

julia> @time devect(x, y, a, b)
elapsed time: 1.3065e-5 seconds (376 bytes allocated)
2x2 Array{Float64,2}:
 0.978099  0.0405639
 0.94757   0.0224974

julia> @time devect(x, y, a, b)
elapsed time: 1.8059e-5 seconds (376 bytes allocated)
2x2 Array{Float64,2}:
 0.978099  0.0405639
 0.94757   0.0224974

可能有更快的方法来执行 devect 解决方案，但请注意分配的字节数的差异。如果去向量化解决方案比向量化解决方案分配更多的内存，它可能是错误的（至少在 Julia 中是这样）。

Answer 2

https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-global-variables

你的代码到处都使用非常量全局变量，这意味着你基本上回到了解释语言的性能领域，因为在编译时不能保证它们的类型。为了快速加速，只需在所有全局变量赋值前加上 const.

Answer 3

您的去向量化代码效率不高。

我做了以下修改：

使所有全局变量常量
我预先分配了输出向量，而不是每次都附加

我展示了两种不同的方法，您可以用更直接的方式对输出进行去向量化

const x = rand(500)
const y = rand(500)
const a = rand()
const b = rand()

function devect()
    dist = Array(Float64, 500)

    for i in 1:500
        dist[i] = sqrt((x[i] - a)^2 + (y[i] - b)^2)
    end

    return [x y][dist .< 0.05,:]
end

function devect2()
    pairs = Array(Float64, 500, 2)

    for i in 1:500
        dist = sqrt((x[i] - a)^2 + (y[i] - b)^2)
        if dist < 0.05
            pairs[i,:] = [x[i], y[i]]
        end
    end

    return pairs
end

function vect()
    d = sqrt((x-a).^2 + (y-b).^2)
    return [x y][d .< 0.05,:]
end

using DataFrames

function df_vect()
    df = DataFrame(x=x, y=y)
    dist = sqrt((df[:x]-a).^2 + (df[:y]-b).^2)

    return df[dist .< 0.05,:]
end

const n = 10^3

@time for i in [1:n] devect() end
@time for i in [1:n] devect2() end
@time for i in [1:n] vect() end
@time for i in [1:n] df_vect() end

输出为

elapsed time: 0.009283872 seconds (16760064 bytes allocated)
elapsed time: 0.003116157 seconds (8456064 bytes allocated)
elapsed time: 0.050070483 seconds (37248064 bytes allocated, 44.50% gc time)
elapsed time: 0.0566218 seconds (30432064 bytes allocated, 40.35% gc time)

Julia：向量化代码与去向量化代码

Julia: vectorized vs. devectorized code

arrays

for-loop

vectorization

julia