# Performance of simple broadcasting operations with many arguments

**URL:** <https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213>\
**Category:** Performance\
**Tags:** performance, broadcast\
**Created:** [November 28, 2021, 9:55pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213 "2021-11-28T21:55:51Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![Marco\_Lombardi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_lombardi/32/14608_2.png) [@Marco\_Lombardi](https://discourse.julialang.org/u/Marco_Lombardi)\
**Post date:** [November 28, 2021, 9:55pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/1 "2021-11-28T21:55:52Z")

</div>

I am puzzled by some microbenchmarks tests I did on _very_ simple broadcasting operations.

```julia
julia> using BenchmarkTools
julia> n = 10000; x = rand(n); y = rand(n); z = rand(n); t = rand(n); r = zeros(n);
julia> @belapsed $r .= $x .+ $y
2.1423333333333336e-6
julia> @belapsed $r .= $x .+ $y .+ $z
3.71025e-6
julia> @belapsed $r .= $x .+ $y .+ $z .+ $t
8.540666666666666e-6

```

What is the reason for the drop in performance when I sum over four vectors?

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 28, 2021, 10:21pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/2 "2021-11-28T22:21:47Z")

</div>

It may be that the processor has a 512 bites cache, thus it can load four 64 bit numbers per cycle. With five involved (including the result) you need two cycles of the processor.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [November 28, 2021, 10:25pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/3 "2021-11-28T22:25:20Z")

</div>

Interestingly, LoopVectorization is getting optimal performance here. I think the problem might be that Julia is requiring an inefficient add order.

---

<div class="post-metadata">

**Author:** ![Marco\_Lombardi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_lombardi/32/14608_2.png) [@Marco\_Lombardi](https://discourse.julialang.org/u/Marco_Lombardi)\
**Post date:** [November 28, 2021, 10:51pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/4 "2021-11-28T22:51:55Z")

</div>

Yes, indeed, LoopVectorization does a very nice job. What do you mean whey you write that “Julia is requiring an inefficient add order”?

---

<div class="post-metadata">

**Author:** ![Marco\_Lombardi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_lombardi/32/14608_2.png) [@Marco\_Lombardi](https://discourse.julialang.org/u/Marco_Lombardi)\
**Post date:** [November 28, 2021, 11:03pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/5 "2021-11-28T23:03:49Z")

</div>

@lmiq this is possibile, but then I would expect that the problem disappears when I am using `Float32` data: instead, I get something even stranger…

```julia
julia> n = 10000; x = rand(Float32, n); y = rand(Float32, n); z = rand(Float32, n); t = rand(Float32, n); r = zeros(Float32, n);
julia> @belapsed $r .= $x .+ $y
1.3168e-6
julia> @belapsed $r .= $x .+ $y .+ $z
1.8634e-6
julia> @belapsed $r .= $x .+ $y .+ $z .+ $t
8.267666666666665e-6

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [November 28, 2021, 11:14pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/6 "2021-11-28T23:14:54Z")

</div>

since floating point math isn’t isn’t associative, llvm isn’t allowed to reassociate which can improve performance. loopvectorization on the other hand does look for the best ordering.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [November 28, 2021, 11:21pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/7 "2021-11-28T23:21:29Z")

</div>

FastBroascast.jl also does well:

```julia
julia> n = 10000; x = rand(n); y = rand(n); z = rand(n); t = rand(n); r = zeros(n);

julia> using FastBroadcast, LoopVectorization

julia> @btime $r .= $x .+ $y .+ $z .+ $t;
  6.558 μs (0 allocations: 0 bytes)

julia> @btime @.. $r = $x + $y + $z + $t; # FastBroadcast
  2.297 μs (0 allocations: 0 bytes)

julia> @btime @turbo $r .= $x .+ $y .+ $z .+ $t; # LoopVectorization
  2.275 μs (0 allocations: 0 bytes)

```

The problem is that while FastBroadcast and LoopVectorization are using SIMD instructions, regular broadcasting is not.

There are exponential possibilities in the number of arguments when broadcasting. This causes the compiler to simply give up instead of considering them all.  
FastBroadcast handles this by simply prioritizing the most common case – the case where none of the arguments are being broadcast (i.e., none of the arguments have an axis with `size` of `1`) – and guaranteeing efficient code is being generated for that particular case.  
Other cases will be slow, but should still be correct.

LoopVectorization handles this via zeroing out strides associated with size-1 axis and SIMD-ing. It currently does not support dynamic size-1 in a contiguous axis. As there is only one axis and that axis is contiguous, this means it’s doing something more or less the same as FastBroadcast.  
LoopVectorization can thus make more scenarios fast than FastBroadcast, although be sure to make sure all contiguous axis are either full-size or statically known to be size 1 (e.g., an `Adjoint{Float64,Vector{Float64}}` is okay, as is `rand(1000,1)`, but `rand(1,1000)` is not). This restriction can be lifted, but I haven’t made it a priority / it’ll be addressed at some point after the rewrite.

Note that given few enough arguments, base broadcasting can actually be fast, as it’ll try to enumerate the possibilities and generate fast code for each case. This means it should be able to match the others for the common case then, while perhaps even being faster for some others due to having created specialized cases for each of them. Obviously, that doesn’t scale with the number of arguments, as we see here where it simply gives up and generates slow code.

---

<div class="post-metadata">

**Author:** ![Marco\_Lombardi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_lombardi/32/14608_2.png) [@Marco\_Lombardi](https://discourse.julialang.org/u/Marco_Lombardi)\
**Post date:** [November 29, 2021, 10:41am UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/8 "2021-11-29T10:41:02Z")

</div>

Thank yoou @Oscar_Smith and @Elrod , your explanations clarify a lot of things. Still, I am not entirely sure what is happening here…

Let’s start with SIMD, which as you write is not used by regular broadcasting. Why is that? Is this something related to the number of arguments in the broadcast, or a general “feature”? I am asking because the drop in performance does depend on the number of vectors added.

Regarding your other point, I understand there are many different ways to perform a broadcasting, but if one does not use the associativity of the sum (which Julia does not by default) I can think only of a single sensible way of doing a broadcasting: just use a single loop and perform the full sum of each element inside the loop. Other options would require multiple passes (i.e., multiple loops), and are likely to be slower. To test the performance, I did a few more tests:

```julia
n = 10000; x = rand(n); y = rand(n); z = rand(n); t = rand(n); r = zeros(n);

function test4a!(r, x, y, z, t)
    @. r = x + y + z + t
end

function test4b!(r, x, y, z, t)
    @turbo @. r = x + y + z + t
end

function test4c!(r, x, y, z, t)
    @inbounds for i ∈ 1:10000
        r[i] = x[i] + y[i] + z[i] + t[i]
    end
end

function test4d!(r, x, y, z, t)
    @inbounds @simd for i ∈ 1:10000
        r[i] = x[i] + y[i] + z[i] + t[i]
    end
end

```

I then measured the performance as usual:

```julia
julia> @btime test4a!($r, $x, $y, $z, $t);
  8.319 μs (0 allocations: 0 bytes)

julia> @btime test4b!($r, $x, $y, $z, $t);
  4.787 μs (0 allocations: 0 bytes)

julia> @btime test4c!($r, $x, $y, $z, $t);
  4.897 μs (0 allocations: 0 bytes)

julia> @btime test4d!($r, $x, $y, $z, $t);
  4.913 μs (0 allocations: 0 bytes)

```

If I have to conclude something, I would say that

- `LoopVectorization` is just doing the obvious thing, loop over each element.
- SIMD does not help at all: either the macro is ignored, or for some reason I do not understand SIMD cannot be used here. Note that the same happens if I define functions `test3d!` or `test2d!` that only add three or two vectors: `@simd` is not improving the performance.
- The standard broadcasting does something I still do not understand. Is there a way to understand the way the loop is performed in Julia? That clearly would help me (and hopefully other people) in writing efficient code, without being forced to either guess or to test each single line as I am doing here.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [November 29, 2021, 11:07am UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/9 "2021-11-29T11:07:33Z")

</div>

There are 2^num\_vector ways:  
E.g., with two input vectors:

1. Broadcast none (both full size)
2. Broadcast 1 (first is of length 1, second is full size)
3. Broadcast 2 (second is of length 1, first is full size)
4. Broadcast both (both have length 1)

Add a third and we have 8 combinations. A fourth, and we’re at 16.

Associativity has nothing to do with it. Your `@inbounds` loop SIMDs just fine, which is why the macro doesn’t help.  
The problem is broadcasting, i.e. replicating size==1 axes.

With two many arguments, there are too many combinations of arguments having length 1 vs full length, hence the compiler gives up.

My previous comment describes the approaches that FastBroadcast.jl and LoopVectorization.jl use to avoid this problem.

---

<div class="post-metadata">

**Author:** ![Marco\_Lombardi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_lombardi/32/14608_2.png) [@Marco\_Lombardi](https://discourse.julialang.org/u/Marco_Lombardi)\
**Post date:** [November 29, 2021, 12:45pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/10 "2021-11-29T12:45:45Z")

</div>

Thank you again. Reading the Julia documentation I thought that a broadcasting operation would generally fuse and thus that, if the `Broadcast.broadcasted` method is not rewritten to make an eager evaluation, the evaluation would be always lazy. This also corresponds to the fact that there are no allocations in none of the tests functions I made (I would expect allocations if the computation is not done elementwise, because I would have intermediate results to store). Also, I do not see your point with size==1 axes: I only have vectors, so just one size and different from 1.

So, if you have a minute, could you kindly help me understand how could Julia possibly perform the operation by writing a code with a for loop that simulates the possible order of operation basic broadcasting is doing? Thank you so much!

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [November 29, 2021, 12:58pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/11 "2021-11-29T12:58:55Z")

</div>

Order of operations has nothing to do with it, which is why I never brought it up or discussed it in this thread.

> [@Marco\_Lombardi](#):
>
> Also, I do not see your point with size==1 axes: I only have vectors, so just one size and different from 1.

Perhaps this makes it more clear, here we write one broadcast statement that corresponds to 8 different loops because we have 3 arguments (2^3 = 8):

```julia
w = Vector{Float64}(undef, 4);
for i in (1,length(w))
  x = rand(i)
  for j in (1,length(w))
    y = rand(j)
    for k in (1,length(w))
      z = rand(k)
      @. w = x + y + z
    end
  end
end

```

Semantically, the single broadcasting statement represents these 8 different loops. If we add a fourth input argument, it has to be able to represent 16 different loops.  
Which of the loops that it actually is is not known at compile time. Only at runtime.

Basically, you might not be using sizes equal to 1, but the compiler does not know that.

Therefore, the compiler has to generate code able to handle all the 2^number\_arguments cases.  
Because this can be a lot, what is the compiler to do?  
There are a lot of possibilities, but broadcasting just gives up if the number is too big.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 29, 2021, 1:19pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/12 "2021-11-29T13:19:26Z")

</div>

I don’t understand why that works at all, instead of simply throwing a dimension mismatch when the arrays don’t have all the same shape:

```julia
julia> w = zeros(2);

julia> @. w = [0.4098911891015198] + [0.25666422432859126, 0.18429878091474228]
2-element Vector{Float64}:
 0.666555413430111
 0.5941899700162621

julia> @. w = rand(1) + rand(2)
ERROR: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(1),), b has dims (Base.OneTo(2),), mismatch at 1")
Stacktrace:

```

I thought that a broadcasting operation like that would be expanded to something as simple as:

```julia
julia> function f!(r,x,y,z,t)
           for i in 1:length(r)
               @inbounds r[i] = x[i] + y[i] + z[i] + t[i]
           end
       end
f! (generic function with 1 method)

```

which does SIMD.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [November 29, 2021, 1:22pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/13 "2021-11-29T13:22:48Z")

</div>

> [@lmiq](#):
>
> `julia> @. w = rand(1) + rand(2)`

You’re broadcasting the `rand` calls. That is, for every element of `w`, it is calling `rand(2)` and `rand(1)`.

For the loop `f!`, you’d need something like `length(x) == 1 ? 1 : i` inside each `getindex` instead of just `i`.

Perhaps that makes it more clear why it is actually impressive that broadcast can compete with loops as often as it does.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 29, 2021, 1:29pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/14 "2021-11-29T13:29:34Z")

</div>

Ok, but I mean, I never expected that broadcasting to work, for general consistency with the fact that in Julia matrices, vectors and scalars are all different things:

```julia
julia> @. [1] + [1,2]
2-element Vector{Int64}:
 2
 3

julia> @. [1,2] + [1,2,3]
ERROR: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 2 and 3")

```

I’m possibly interpreting things wrong, but the fact that that one-element vector acts as a scalar in that broadcasting operation sort of goes against the choice of making other one-element arrays throwing errors when inadvertently used as scalars (for ex. `[1] * [1,2]`).

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [November 29, 2021, 1:37pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/15 "2021-11-29T13:37:50Z")

</div>

It can be convenient, particularly for multidimensional arrays, e.g.

```julia
using Statistics 
A = rand(10,10,10);
Z = (A .- mean(A, dims=(1,2))) ./ std(A, dims=(1,2))

```

Broadcasting always replicates size==1 axes across the entire corresponding axis.

It would be nice to have something that doesn’t broadcast though, to make simple loops more convenient to write. Maybe FastBroadcast.jl should fill that role.

---

<div class="post-metadata">

**Author:** ![Marco\_Lombardi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_lombardi/32/14608_2.png) [@Marco\_Lombardi](https://discourse.julialang.org/u/Marco_Lombardi)\
**Post date:** [November 29, 2021, 10:09pm UTC](https://discourse.julialang.org/t/performance-of-simple-broadcasting-operations-with-many-arguments/72213/16 "2021-11-29T22:09:48Z")

</div>

Thank you, this really clarifies everything! And I do appreciate the speed simple broadcasting is generally able to reach at runtime, considering all the various possibilities. In reality, in same cases I have not even been able to fully match its speed with a for loop (as in [#71336](https://discourse.julialang.org/t/type-inference-for-recursive-broadcasting-operations/71336))…