Significantly slower broadcasting

[astrozot]

I have been testing the speed of broadcast operations with OffsetArrays of StaticArrays, and it looks like there is a significant penalty in time. In fact, on my laptop I see

julia> using StaticArrays, OffsetArrays, BenchmarkTools

julia> xs = [SVector{2}(rand(2)) for _ ∈ 1:10_000];

julia> ys = similar(xs);

julia> @benchmark (@. $ys = $xs / (first($xs)^2 + last($xs)^2))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):   9.419 μs … 43.605 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     10.142 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   11.559 μs ±  2.524 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █    ▁ ▅                                                     
  █▅▆▄▆█▃█▂▃▂▂▃▂▂▃▃▂▂▃▂▁▂▂▂▂▂▂▂▂▃▂▁▂▂▁▁▁▁▁▃▄▅▂▂▂▄▂▃▁▃▂█▆▂▁▁▅▆ ▃
  9.42 μs         Histogram: frequency by time        15.2 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

When using OffsetArrays instead these are the results:

julia> oxs = OffsetArray(xs, -1000);

julia> oys = OffsetArray(ys, -1000);

julia> @benchmark (@. $oys = $oxs / (first($oxs)^2 + last($oxs)^2))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  48.245 μs … 159.888 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     50.522 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   53.394 μs ±  10.011 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▇███▄▁            ▃  ▃ ▂    ▁  ▄   ▃                         ▂
  ██████▇▇▁▇▆▆▇▅▇▃█▆█▇██▆█▇▆█▅█▆▅█▆▆▅█▆▆▆▆▆▆█▆▆▆▅▆▆▅▆▆▅▆▆▆▆▆▆▆ █
  48.2 μs       Histogram: log(frequency) by time      97.1 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

The $> 5 \times$ penalty is even larger than the SIMD vector size of my laptop (4 for Float64 arrays).

Any help or clarification is really appreciated!

[astrozot]

Just to clarify: the issue disappears when working with normal arrays, as in

julia> as = rand(10_000);

julia> bs = similar(as);

julia> @benchmark (@. $bs = $as / (1.0 + $as^2))
BenchmarkTools.Trial: 10000 samples with 7 evaluations.
 Range (min … max):  4.668 μs … 20.108 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     6.278 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   6.231 μs ±  1.076 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▇   ▁      ▅▁▇▄█▄▇▃▆▁▂▁▁                                   ▂
  ██████▃██▁███████████████▇▇▇▆▆▆▇▇▇▇▇▇▇▇▇▇▇▅▆▆▇▆▄▅▆▆▆▅▆▅▅▄▆ █
  4.67 μs      Histogram: log(frequency) by time     10.8 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

It also disappears if, instead of broadcasting operations, I use explicit loops:

julia> function f!(ys, xs)
           for i ∈ eachindex(xs)
               @inbounds ys[i] = xs[i] / (first(xs[i])^2 + last(xs[i])^2)
           end
           ys
       end
f! (generic function with 1 method)

julia> @benchmark f!($ys, $xs)
BenchmarkTools.Trial: 10000 samples with 3 evaluations.
 Range (min … max):  8.789 μs … 39.339 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     9.053 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   9.876 μs ±  2.016 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▇▅▅▄▅▅▃▄ ▂            ▃ ▂  ▂  ▂  ▁▄ ▃ ▄                   ▂
  █████████▅█▃██▃▁█▁█▁▇▆▁█▃█▁▃█▁▁█▁▄██▅█▆█▆▆▅█▅▄▄▃▆▄▅▃▃▁▁▄▄▄ █
  8.79 μs      Histogram: log(frequency) by time     16.5 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark f!($oys, $oxs)
BenchmarkTools.Trial: 10000 samples with 3 evaluations.
 Range (min … max):   8.808 μs … 44.438 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):      9.446 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   10.056 μs ±  1.749 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █ ▆▂▇ ▇ ▅  ▇ ▇▂             ▄   ▃   ▃  ▁    ▁   ▅▂   ▄    ▁ ▂
  █▅███▅█▆██▆█▄██▃█▁▄█▃▁█▁▁█▁▁█▇▁▁█▃▁▃█▃▁█▇▁▁▁█▁▅▅██▆▆▅█▅▅▆▆█ █
  8.81 μs      Histogram: log(frequency) by time      14.5 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

[jishnub]

I think the issue is that StaticArrays unrolls the loops, and OffsetArrays currently doesn't have specialized broadcast styles defined and can't access the specialized methods.

[astrozot]

Thank you @jishnub , thank you for your reply. Looking at the source code, I saw a link to a similar issue reported here:

https://discourse.julialang.org/t/why-is-there-a-performance-hit-on-broadcasting-with-offsetarrays/32194

I am not entirely sure the situation is identical...