# Arithmetic broadcasting in Julia 5x slower than MATLAB

**URL:** <https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628>\
**Category:** Performance\
**Created:** [May 25, 2022, 4:37am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628 "2022-05-25T04:37:46Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [May 25, 2022, 4:37am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/1 "2022-05-25T04:37:46Z")

</div>

I was surprised to find a fairly simple broadcasting operation on somewhat large arrays be considerably slower in Julia relative to MATLAB. Is there some unattainable level of optimization that folks at MathWorks have figured out or can I improve my Julia code to match/beat it?

Here’s an example. To give Julia the best chance, I’m going to use mutating functions to keep the number of allocations to zero. I tried it with allocating functions too – they’re predictably slower.

```julia
function f!(D,A,B,C)
	D .= max.(A.*B.+C, 0.2)
end

A = rand(100,100,1);
B = rand(1,100,100);
C = rand(100,1,100);
D = zeros(100,100,100);

@btime f!($D,$A,$B,$C);

```

gives me `1.822 ms (0 allocations: 0 bytes)`. The MATLAB equivalent is 0.35 ms, 5x slower (see MATLAB code below).

I tried decomposing f! into the arithmetic operation and the max to see if one is particularly slower. Not really:

```julia
function g!(D,A,B,C)
	D .= A.*B.+C
end

function h!(D)
	D .= max.(D,0.2)
end

@btime g!($D,$A,$B,$C);
@btime h!($D);

```

and I get

```julia
  1.062 ms (0 allocations: 0 bytes)
  1.146 ms (0 allocations: 0 bytes)

```

so each takes about the same amount of time.

Compare to MATLAB:

```julia
A = rand(100,100,1);
B = rand(1,100,100);
C = rand(100,1,100);

tf = 0; tg = 0; th = 0;
N = 10;
for i=1:N
	tic; tmp = f(A,B,C); ti = toc; tf = tf+ti;
	tic; tmp = g(A,B,C); ti = toc; tg = tg+ti;
	tic; tmp = h(tmp); ti = toc; th = th+ti;
end
fprintf(' %0.3f ms\n', tg/N*1000);
fprintf(' %0.3f ms\n', th/N*1000);
fprintf(' %0.3f ms\n', tf/N*1000);

function D = f(A,B,C)
	D = max(A.*B+C, 0.2);
end

function D = g(A,B,C)
	D = A.*B+C;
end

function D = h(D)
	D = max(D, 0.2);
end

```

which prints

```julia
 0.339 ms
 0.483 ms
 0.350 ms

```

Notice that:

- Component functions g and h are at least 2x faster in MATLAB (0.339 ms vs 1.06 ms and 0.483 ms vs 1.146 ms)
- The combined function f is about the same speed as _each_component in MATLAB (0.35 ms ~ 0.339 and 0.483), while the corresponding function in Julia is just a little faster than the sum of the components (1.062 + 1.146 ~ 1.822).

Clearly, MATLAB must be doing some low-level optimization here. I don’t think BLAS is involved since it’s all element operations being broadcasted (there’s no matrix math), but I’m out of my depth when it comes to low-level stuff like this.

So my question is – can I speed up the Julia code? In my application it’s in a hot loop and the arrays are big, so there’s a large payoff to shaving microseconds off the MWE’s running time, if possible.

---

<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:** [May 25, 2022, 4:46am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/2 "2022-05-25T04:46:57Z")

</div>

Is matlab automatically threading this? You can probably make this faster by using LoopVectorization. I see the following at 5x faster than your original `f!`

```julia
function f2!(D,A,B,C)
	@turbo D .= max.(A.*B.+C, 0.2)
end

```

and using `@tturbo` (threaded version of `@turbo`) gives another 4x speedup.

```julia
function f2!(D,A,B,C)
	@tturbo D .= max.(A.*B.+C, 0.2)
end

```

---

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [May 25, 2022, 5:05am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/3 "2022-05-25T05:05:04Z")

</div>

> [@Oscar\_Smith](#):
>
> ```julia
> function f2!(D,A,B,C)
> @turbo D .= max.(A.*B.+C, 0.2)
> end
> 
> ```

Thanks! Good call on the multi-threading. I keep forgetting MATLAB does it by default with a lot of functions. I restarted it with `matlab -singleCompThread` and re-ran my tests and `g()` and `h()` times are now on par with those original Julia times. But `f()` is still as fast as `g()` and `h()` are individually, making it still faster than Julia:

```julia
 1.140 ms
 1.324 ms
 1.301 ms

```

I tried LoopVectorization with `@turbo` and did get \>2x speed ups in each component function _as well_ as the interaction effect that leads to `f()` being about as fast as each component (just like in MATLAB). This indeed yields about a 4x speedup compared the original – you got 5x but that can be hardware difference or sampling error.

Two questions:

1. For `@ttturbo` to get further speed-ups, I need to run Julia with more threads, right? Your 4x speed-up is a function of the number of threads you ran it with, right?
2. What does LoopVectorization do? I just looked at the documentation and I honestly couldn’t figure it out.

---

<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:** [May 25, 2022, 5:11am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/4 "2022-05-25T05:11:03Z")

</div>

1. Yes. `julia --threads=4` to launch Julia with 4 threads
2. LoopVectorization basically replaces Julia’s compiler with a smarter one when you put `@turbo` or `@tturbo` in front of them. Basically it has a cost model and can rewrite loops or vectorized code to more efficiently take advantage of processor vectorization.

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [May 25, 2022, 5:12am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/5 "2022-05-25T05:12:53Z")

</div>

Have you tried using MKL?

---

<div class="post-metadata">

**Author:** ![albheim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albheim/32/34660_2.png) [@albheim](https://discourse.julialang.org/u/albheim)\
**Post date:** [May 25, 2022, 5:13am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/6 "2022-05-25T05:13:07Z")

</div>

> [@elenev](#):
>
> `tic; tmp = f(A,B,C); ti = toc; tf = ti+toc;`

Is this really the way to measure the time? To me it seems like you want tf to sum up all ti, so it should be `tf = tf + ti`, right now it seems like you just would get `tf = 2*toc`, though then it is a little strange that the versions are so equal. l’m not well versed in matlab so might have missed something here.

---

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [May 25, 2022, 5:20am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/8 "2022-05-25T05:20:41Z")

</div>

I have not. Just looked into it and seems like it’s kind of a global change – documentation says you need to load MKL.jl before any other package meaning it has to be in my startup script instead of at the top of whatever module I’m writing.

Do you think it’s worth it for this problem? Do you think it’s worth it elsewhere to the point that I should use it as a default?

---

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [May 25, 2022, 5:21am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/9 "2022-05-25T05:21:09Z")

</div>

You’re right, of course. I copied my code from the wrong window where I wrote it first and made the mistake you pointed out, but the times were still right. I’ll edit the code above.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [May 25, 2022, 5:39am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/10 "2022-05-25T05:39:28Z")

</div>

I think it’s more accurate to use the `timeit` function when benchmarking in Matlab. It will do something similar to Julia’s `@btime`.

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [May 25, 2022, 6:07am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/11 "2022-05-25T06:07:56Z")

</div>

It doesn’t seem to matter for me:

```julia
julia> @btime f!($D,$A,$B,$C);
  2.270 ms (0 allocations: 0 bytes)

julia> using MKL

julia> @btime f!($D,$A,$B,$C);
  2.269 ms (0 allocations: 0 bytes)

```

but this is somewhat CPU dependent I believe and the computer I’m on now is ancient. Overall since Julia 1.7 MKL is not a fixed global setting anymore, you can dynamically change the BLAS at runtime.

---

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [May 25, 2022, 5:24pm UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/12 "2022-05-25T17:24:38Z")

</div>

`@tturbo inline=true` seems to speed things up by another 20-30% or so relative to `@tturbo`.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [May 25, 2022, 6:16pm UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/13 "2022-05-25T18:16:47Z")

</div>

BLAS isn’t involved in broadcasting, afaik.

---

<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:** [May 25, 2022, 6:22pm UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/14 "2022-05-25T18:22:04Z")

</div>

> [@elenev](#):
>
> I tried LoopVectorization with `@turbo` and did get \>2x speed ups in each component function

Could you share your `versioninfo()`?  
Your CPU may have AVX512, which LLVM won’t really take advantage of by default, but LoopVectorization will. This would explain a 2x speedup.  
You can also check if you have AVX512 via:

```julia
julia> LoopVectorization.VectorizationBase.has_feature(Val(:x86_64_avx512f))
static(true)

```

Returns `static(true)` if LV will try to use AVX512, `static(false)` otherwise.

If you do have AVX512, you can try starting julia with `julia -C"native,-prefer-256-bit" -t4`.  
`-t4` is a short hand for using 4 threads.  
`-prefer-256-bit` turns off the option `prefer-256-bit`, i.e. disables the preference for 256 bit vectors over the 512 bits AVX512 offers.

> the corresponding function in Julia is just a little faster than the sum of the components (1.062 + 1.146 ~ 1.822)

Note that much of the cost of these operations is actually the process of iterating over the arrays, both fetching the memory and loading them into registers.  
If the arrays are small, loading them into registers can easily be more expensive than the actual arithmetic.  
If the arrays are large so that they do not fit in high cache levels, memory bandwidth will dominate.

Hence, “fusing” the operations can give a speed up, as it saves you from looping over `D` – by far your largest array – twice.

> BLAS isn’t involved in broadcasting, afaik.

Yes, BLAS is not involved in broadcasting.

> `@tturbo inline=true` seems to speed things up by another 20-30% or so relative to `@tturbo` .

Interesting! I should adjust the heuristic so that it does this automatically for that problem.

---

<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:** [May 25, 2022, 6:35pm UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/15 "2022-05-25T18:35:38Z")

</div>

Oh, also this is critically important to note when using LoopVectorization here:

```julia
A = rand(100,100,1);
B = LowDimArray{(false, true, true)}(rand(1,100,100));
C = rand(100,1,100);

```

You need to wrap arrays in `LowDimArray`s when their contiguous (the first dimension for column major arrays) have size only 1!  
Otherwise, you’ll likely get a segfault, or at least wrong answers!

It probably wouldn’t hurt to give `A` and `B` the same treatment.

I have a plan of how to fix this, but do not have the time myself.  
If anyone else wants to dive in, I’d be happy to explain how.  
Otherwise, that’d be something to expect eventually with the rewrite.

---

<div class="post-metadata">

**Author:** ![elenev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elenev/32/18440_2.png) [@elenev](https://discourse.julialang.org/u/elenev)\
**Post date:** [May 25, 2022, 8:37pm UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/16 "2022-05-25T20:37:36Z")

</div>

`versioninfo()` returns

```julia
Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)        
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Xeon(R) CPU E3-1240 v3 @ 3.40GHz
  WORD_SIZE: 64    
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, haswell)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 4

```

But `has_feature()` returns `static(false)` so looks like that’s not it.

Not sure if I should’ve mentioned I’m still on Julia 1.6. But this is an old project with lots of existing code and I’m worried stuff will break if I upgrade.

> [@Elrod](#):
>
> Hence, “fusing” the operations can give a speed up, as it saves you from looping over `D` – by far your largest array – twice.

I understand that, but I was just surprised that MATLAB seems to fuse “better” than Julia i.e. time(f) / [time(g) + time(h)] in MATLAB was smaller than in Julia, a finding that seems separate from the relative speeds of individual components.

> [@Elrod](#):
>
> You need to wrap arrays in `LowDimArray` s when their contiguous (the first dimension for column major arrays) have size only 1!

Heh, good to know. Thanks. Is this documented somewhere? I’m happy that `@tturbo` sped things up but it still feels more like a random piece of code I copied from a forum than a feature that I understand how and where to use.

---

<div class="post-metadata">

**Author:** ![jacobusmmsmit](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacobusmmsmit/32/217669_2.png) [@jacobusmmsmit](https://discourse.julialang.org/u/jacobusmmsmit)\
**Post date:** [May 25, 2022, 10:31pm UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/17 "2022-05-25T22:31:08Z")

</div>

> [@Oscar\_Smith](#):
>
> LoopVectorization basically replaces Julia’s compiler with a smarter one

If it’s so smart why isn’t it the default? Is it suboptimal sometimes or are the times it can be used rare enough that testing for its use cases everywhere slows down compilation too much?

---

<div class="post-metadata">

**Author:** ![N5N3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/n5n3/32/17663_2.png) [@N5N3](https://discourse.julialang.org/u/N5N3)\
**Post date:** [May 25, 2022, 11:12pm UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/18 "2022-05-25T23:12:16Z")

</div>

> [@jacobusmmsmit](#):
>
> If it’s so smart why isn’t it the default? Is it suboptimal sometimes or are the times it can be used rare enough that testing for its use cases everywhere slows down compilation too much?

1. LoopVectorization is not safe for all loops.
2. The compilation cost is higher. The first excution needs \> 5s with a fresh REPL. (latter call would be faster)
3. We don’t have it in Base.

---

<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:** [May 26, 2022, 12:46am UTC](https://discourse.julialang.org/t/arithmetic-broadcasting-in-julia-5x-slower-than-matlab/81628/19 "2022-05-26T00:46:42Z")

</div>

a rewrite of it is currently in the works that will make it much easier to incorporate as a compiler pass