Why is this broadcast operation slower than a nested for-loop

user_2546487 0
Viewed 151

I have a function that checks if every element in an array is greater than a lower bound and lesser than an upper bound.

In the following code, I'm curious to know why bounds_error (a nested for-loop) is faster than bounds_error2 (vectorized, broadcast operation) and what can one do to make this function run faster.

using BenchmarkTools

function bounds_error(x, xl)
    num_x_rows = size(x,1)
    num_dim = size(xl, 1)
    for i in 1:num_x_rows
        for j in 1:num_dim
            if (x[i, j] < xl[j,1] || x[i,j] > xl[j,2])
                return true
            end
        end
    end
    return false
end

function bounds_error2(x, xl)
    for row in eachrow(x)
        xlt = transpose(xl)
        if any(row .< xlt[1, :]) == true || any(row .> xlt[2, :])
            return true
        end
    end
    return false
end

#number of rows in xl (or xlimits) will always be equal to number of columns in x

xl =  [     -5.0  5.0
            -5.0  5.0
            -5.0  5.0]

x = [1.0 2.0 3.0; 
     4.0 5.0 6.0]

@btime bounds_error(x, xl) #~20.645 ns (0 allocations: 0 bytes); true

@btime bounds_error2(x, xl) #~347.870 ns (12 allocations: 704 bytes); true
julia
1 Answers

The main reason for this difference is memory allocations (0 vs. 12 here).

#  20.645 ns (0 allocations: 0 bytes)
# 347.870 ns (12 allocations: 704 bytes)

Currently, slices in Julia create a copy, so xlt[1,:] and xlt[2,:] allocates memory. To remedy this problem you should use @views. The second issue is the element-wise comparisons row .< xlt[1,:] and row .> xlt[2,:] create a temporary Boolean array. To avoid allocation of a temporary array, you should map any(t->t[1]<t[2], zip(row,xl1)) so that the comparison is done one element at a time like a loop.

After applying these tips, the performance difference on my machine is now about 2ns only, which accounts for the convenience of eachrow, zip, etc. instead of manual loops.

Note, for the first function, you can use axes() to loop over first or second dimension conveniently. And when benchmarking any Julia code with BenchmarkTools.jl, don't forget to interpolate ($) all variable names of a function to avoid working on global variables.

function bounds_error(x, xl)
    for i in axes(x,1)
        for j in axes(xl, 1)
            if (x[i, j] < xl[j,1] || x[i,j] > xl[j,2])
                return true
            end
        end
    end
    return false
end

@views function bounds_error2(x, xl)
    xl1, xl2 = xl[:,1], xl[:,2]
    for row in eachrow(x)
        if any(t->t[1]<t[2], zip(row,xl1)) || any(t->t[1]>t[2], zip(row,xl2))
            return true
        end
    end
    return false
end

# number of rows in xl (or xlimits) will always be equal to number of columns in x
xl = [-5.0  5.0
      -5.0  5.0
      -5.0  5.0]

x = [1.0 2.0 3.0; 
     4.0 5.0 6.0]

@btime bounds_error($x, $xl)  #  8.100 ns (0 allocations: 0 bytes)
@btime bounds_error2($x, $xl) # 10.800 ns (0 allocations: 0 bytes)
user_39431700
Related