Julia: vectorized vs. devectorized code

3
votes

From what I understood, Julia is supposed to make for loops faster and as fast as vectorized operations. I wrote three versions of a simple function that finds distance using for loops vs. a vectorized operation vs. the latter with 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

Output:

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)

Why does the vectorized version perform so much faster?

3
arraysfor-loopvectorizationjulia
This might be useful: software.intel.com/en-us/articles/vectorization-in-julia - user3467349
This may a global scope issue. Try rewriting your functions so they use local variables rather than global. That often impacts timings. - ptb
@ptb I originally did that (declared a, b, x, and y inside each function). Didn't help. - Anarcho-Chossid
Using inbounds and simd doesn't speed up devect() much. - Anarcho-Chossid
Looking a little more closely at your code in devect(), The huge slow down is likely in the way you are constructing your return values. They start as zero element arrays and then you reassign (i.e. reallocate) them to new arrays at every iteration. - ptb

3 Answers

12
votes

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

Your code uses non-constant globals everywhere, which means that you're basically back in the performance realm of interpreted languages since no guarantees can be made at compile time about their type. For a quick speedup, just prefix all your global variables assignments with const.

7
votes

Following up on my comment about the method used to construct the solution in devect. Here's my code

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

There may be faster ways to do the devect solution but notice the difference in bytes allocated. If a devectorized solution allocates more memory than a vectorized solution, it's probably wrong (at least in Julia).

2
votes

You devectorized code is not very efficient.

I made the following changes:

  • Made all the global variables constant
  • I preallocated the output vector, rather than appending each time

  • I show two different ways you could have devectorized the output in a more straightforward way

    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

The output is 

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)