# Array broadcasting slower than numpy?

**URL:** <https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212>\
**Category:** Performance\
**Created:** [June 3, 2022, 8:33pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212 "2022-06-03T20:33:44Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![chunjiw](https://avatars.discourse-cdn.com/v4/letter/c/e79b87/32.png) [@chunjiw](https://discourse.julialang.org/u/chunjiw)\
**Post date:** [June 3, 2022, 8:33pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/1 "2022-06-03T20:33:44Z")

</div>

```julia
@time let m = 10^7 + 19
	v = fill(0, m)
	w = copy(v)
	for _ in 1:100
		i = mod(rand(Int), m)
		v[1:i] .+= @view w[end-i+1:end]
		v[i+1:end] .+= @view w[1:end-i]
	end
end

```

```python
import time
import numpy as np
import random

t0 = time.time()
m = 10000019
v = np.array([0]*m)
w = np.copy(v)
for _ in range(100):
    i = random.randint(0, m-1)
    if i > 0:
        v[:i] += w[-i:]
        v[i:] += w[:-i]
    else:
        v += w
print(time.time() - t0)

```

On my machine, the Julia version takes 4.39 seconds, and the Python version takes only 1.85 seconds.  
Is there any improvement to the Julia code to match numpy’s performance?

P.S. This came up when I was trying to do the [Euler Project No. 655](https://projecteuler.net/problem=655), during which I learned a DP method involving array broadcasting.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [June 3, 2022, 8:48pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/2 "2022-06-03T20:48:18Z")

</div>

> [@chunjiw](#):
>
> ```julia
> v[1:i] .+= @view w[end-i+1:end]
> v[i+1:end] .+= @view w[1:end-i]
> 
> ```

change to

```julia
               @views v[1:i] .+= w[end-i+1:end]
               @views v[i+1:end] .+= w[1:end-i]

```

---

<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:** [June 3, 2022, 9:24pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/3 "2022-06-03T21:24:46Z")

</div>

Single threaded, I’m getting Base \< FastBroadcast.jl \< NumPy \< LoopVectorization.jl:

```julia
julia> @time let m = 10^7 + 19
               v = fill(0, m)
               w = copy(v)
               for _ in 1:100
                       i = mod(rand(Int), m)
                       @views @. v[1:i] += w[end-i+1:end]
                       @views @. v[i+1:end] += w[1:end-i]
               end
       end
  3.530138 seconds (4 allocations: 152.588 MiB, 0.82% gc time)

julia> using FastBroadcast

julia> @time let m = 10^7 + 19
               v = fill(0, m)
               w = copy(v)
               for _ in 1:100
                       i = mod(rand(Int), m)
                       @views @.. v[1:i] += w[end-i+1:end]
                       @views @.. v[i+1:end] += w[1:end-i]
               end
       end
  3.177666 seconds (4 allocations: 152.588 MiB, 0.91% gc time)

julia> using LoopVectorization

julia> @time let m = 10^7 + 19
               v = fill(0, m)
               w = copy(v)
               for _ in 1:100
                       i = mod(rand(Int), m);
                       @turbo view(v,1:i) .+= view(w,m-i+1:m)
                       @turbo view(v,i+1:m) .+= view(w,1:m-i)
               end
       end
  1.564778 seconds (4 allocations: 152.588 MiB)

```

Python:

```python
In [3]: import time
   ...: import numpy as np
   ...: import random
   ...:
   ...: t0 = time.time()
   ...: m = 10000019
   ...: v = np.array([0]*m)
   ...: w = np.copy(v)
   ...: for _ in range(100):
   ...: i = random.randint(0, m-1)
   ...: if i > 0:
   ...: v[:i] += w[-i:]
   ...: v[i:] += w[:-i]
   ...: else:
   ...: v += w
   ...: print(time.time() - t0)
1.9152345657348633

```

FastBroadcast.jl and LoopVectorization.jl are easy to multithread:

```julia
julia> @time let m = 10^7 + 19
               v = fill(0, m)
               w = copy(v)
               for _ in 1:100
                       i = mod(rand(Int), m)
                       @views @.. thread=true v[1:i] += w[end-i+1:end]
                       @views @.. thread=true v[i+1:end] += w[1:end-i]
               end
       end
  0.361589 seconds (4 allocations: 152.588 MiB, 2.10% gc time)

julia> @time let m = 10^7 + 19
               v = fill(0, m)
               w = copy(v)
               for _ in 1:100
                       i = mod(rand(Int), m);
                       @tturbo view(v,1:i) .+= view(w,m-i+1:m)
                       @tturbo view(v,i+1:m) .+= view(w,1:m-i)
               end
       end
  0.372530 seconds (4 allocations: 152.588 MiB)

```

FastBroadcast and base broadcast are failing to SIMD here, but this code is so memory bound when multithreaded that it doesn’t matter.

---

<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:** [June 3, 2022, 9:49pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/5 "2022-06-03T21:49:43Z")

</div>

The problem, as it so often is, was with timing in the global scope.

```julia
julia> broadbench() = let m = 10^7 + 19
               v = fill(0, m)
               w = copy(v)
               for _ in 1:100
                       i = mod(rand(Int), m)
                       @views @. v[1:i] += w[end-i+1:end]
                       @views @. v[i+1:end] += w[1:end-i]
               end
       end
broadbench (generic function with 1 method)

julia> fastbroadbench() = let m = 10^7 + 19
               v = fill(0, m)
               w = copy(v)
               for _ in 1:100
                       i = mod(rand(Int), m)
                       @views @.. v[1:i] += w[end-i+1:end]
                       @views @.. v[i+1:end] += w[1:end-i]
               end
       end
fastbroadbench (generic function with 1 method)

julia> broadbench() = let m = 10^7 + 19
               v = fill(0, m)
               w = copy(v)
               for _ in 1:100
                       i = mod(rand(Int), m)
                       @views @. v[1:i] += w[end-i+1:end]
                       @views @. v[i+1:end] += w[1:end-i]
               end
       end
broadbench (generic function with 1 method)

julia> turbobroadbench() = let m = 10^7 + 19
               v = fill(0, m)
               w = copy(v)
               for _ in 1:100
                       i = mod(rand(Int), m);
                       @turbo view(v,1:i) .+= view(w,m-i+1:m)
                       @turbo view(v,i+1:m) .+= view(w,1:m-i)
               end
       end
turbobroadbench (generic function with 1 method)

```

This yields:

```julia
julia> @time broadbench()
  1.337521 seconds (4 allocations: 152.588 MiB, 1.81% gc time)

julia> @time fastbroadbench()
  1.331893 seconds (4 allocations: 152.588 MiB, 1.49% gc time)

julia> @time turbobroadbench()
  1.331875 seconds (4 allocations: 152.588 MiB, 1.35% gc time)

```

All much better than Python.

---

<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:** [June 3, 2022, 10:14pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/6 "2022-06-03T22:14:51Z")

</div>

> [@jling](#):
>
> change to
> 
> ```julia
> @views v[1:i] .+= w[end-i+1:end]
> @views v[i+1:end] .+=
> 
> ```

Why? Shouldn’t the slices on the left always be views?

---

<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:** [June 3, 2022, 10:23pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/7 "2022-06-03T22:23:27Z")

</div>

`+=` is the problem.

```julia
julia> Meta.@lower v[1:i] .+= @view w[end-i+1:end]
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope`
1 ─ %1 = 1:i
│ %2 = Base.dotview(v, %1)
│ %3 = Base.getindex(v, %1)
└── goto #3 if not true
2 ─ %5 = (lastindex)(w)
│ %6 = %5 - i
│ %7 = %6 + 1
│ %8 = (lastindex)(w)
│ %9 = %7:%8
│ #s13 = (view)(w, %9)
└── goto #4
3 ─ #s13 = false
4 ┄ %13 = #s13
│ %14 = Base.broadcasted(+, %3, %13)
│ %15 = Base.materialize!(%2, %14)
└── return %15
))))

```

We write into a view, but we add using a `getindex(v, 1:i)`.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [June 3, 2022, 10:37pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/8 "2022-06-03T22:37:07Z")

</div>

Should this be chaged? CAN this be changed?

---

<div class="post-metadata">

**Author:** ![chunjiw](https://avatars.discourse-cdn.com/v4/letter/c/e79b87/32.png) [@chunjiw](https://discourse.julialang.org/u/chunjiw)\
**Post date:** [June 3, 2022, 10:47pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/9 "2022-06-03T22:47:19Z")

</div>

Wow it really makes a difference. Still I’m wondering why `@views` is different from `@view` here?

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [June 3, 2022, 10:47pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/10 "2022-06-03T22:47:44Z")

</div>

I’m not sure if this is a simpler way to look at it, but to quote the docs, “`X .+= Y` etcetera is equivalent to `X .= X .+ Y`”. So `v[1:i] .+= @view w[end-i+1:end]` is equivalent to `v[1:i] .= v[1:i] .+ @view w[end-i+1:end]`. Just a hazard of the shorthand `+=`'s left side actually representing duplicate subexpressions on the left and right side of the actual assignment.

> [@Elrod](#):
>
> The problem, as it so often is, was with timing in the global scope.

I am actually surprised here, because the code was in a local `let` block to begin with, and all the benchmarks are first-call `@time`. I can’t spot any globals, and the fact that it’s always 4 allocations and 152 MiB suggests to me there weren’t any type instabilities being fixed.

---

<div class="post-metadata">

**Author:** ![chunjiw](https://avatars.discourse-cdn.com/v4/letter/c/e79b87/32.png) [@chunjiw](https://discourse.julialang.org/u/chunjiw)\
**Post date:** [June 3, 2022, 10:50pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/11 "2022-06-03T22:50:06Z")

</div>

Immediately answered by Benny’s reply lol

---

<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:** [June 3, 2022, 11:09pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/12 "2022-06-03T23:09:13Z")

</div>

> [@Benny](#):
>
> I am actually surprised here, because the code was in a local `let` block to begin with, and all the benchmarks are first-call `@time`. I can’t spot any globals, and the fact that it’s always 4 allocations and 152 MiB suggests to me there weren’t any type instabilities being fixed.

Perhaps it is more accurate to say the problem is that it isn’t compiled, and the overhead is high.  
Dispatches often result in zero allocations.

This allocates a lot because it compiles:

```julia
julia> @time let m = 10^7 + 19; (m -> begin
               v = fill(0, m)
               w = copy(v)
               for _ in 1:100
                       i = mod(rand(Int), m)
                       @views @. v[1:i] += w[end-i+1:end]
                       @views @. v[i+1:end] += w[1:end-i]
               end
       end)(m); end
  1.451741 seconds (369.68 k allocations: 177.840 MiB, 1.52% gc time, 6.56% compilation time)

```

but is still \>2x faster, despite the compilation.

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [June 3, 2022, 11:33pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/13 "2022-06-03T23:33:57Z")

</div>

Ah that’s right, Julia compiles per first method call. The let block may have no type instabilities, but it’s not going to benefit from the optimizing compiler. Makes sense, why bother optimizing something that can’t run repeatedly?

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [June 3, 2022, 11:37pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/14 "2022-06-03T23:37:02Z")

</div>

The difference is that

```julia
v[1:i] .+= @view w[end-i+1:end]

```

applies the view only to `w[end-i+1:end]`, instead

```julia
@views v[1:i] .+= w[end-i+1:end]

```

applies the view also to `v[1:i]`.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 3, 2022, 11:48pm UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/15 "2022-06-03T23:48:23Z")

</div>

Note also that it’s a lot easier to just apply `@views` to the whole `let` block (or the whole function); there’s typically not much reason to apply it line-by-line.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [June 4, 2022, 12:15am UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/16 "2022-06-04T00:15:52Z")

</div>

Honestly should we change this? What’s the possible use to the current behavior?

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [June 4, 2022, 12:33am UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/17 "2022-06-04T00:33:01Z")

</div>

I think the current behaviour makes _some_ sense: `a += b` is lowered to `a = a + b`:

```julia
julia> Meta.@lower a += b
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope`
1 ─ %1 = a + b
│ a = %1
└── return %1
))))

```

and since a view wasn’t _explicitly_ applied to `a` in the right-hand side, the `a` in left-hand side doesn’t have a view either.

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [June 4, 2022, 12:33am UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/18 "2022-06-04T00:33:31Z")

</div>

The [@view docs](https://docs.julialang.org/en/v1/base/arrays/#Base.@view) seem pretty comprehensive? `@view` for individual indexing subexpressions, `@views` to slap `@view` on every place in an expression that would be a `getindex`, and `@view` works on the left side of a broadcasted assignment (`.=`) even with updating operators (`.+=`).

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [June 4, 2022, 12:35am UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/19 "2022-06-04T00:35:55Z")

</div>

Yeah I guess there isn’t any natural place where we can special case this huh…

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [June 4, 2022, 5:26am UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/20 "2022-06-04T05:26:13Z")

</div>

Let’s just default to views in Julia 2.0 🙂

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [June 4, 2022, 5:58am UTC](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212/21 "2022-06-04T05:58:11Z")

</div>

While we’re optimizing, you could try making `m` an `Int` by writing it the same way you did in Python (and can we maybe have `1E7` or something be an integer in Julia 2.0?)

[Next page](https://discourse.julialang.org/t/array-broadcasting-slower-than-numpy/82212.md?page=2)