# Julia is significantly slower (~18 x) than Matlab in vector and matrix algebra

**URL:** <https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773>\
**Category:** New to Julia\
**Created:** [June 24, 2023, 6:22am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773 "2023-06-24T06:22:47Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![DY\_K](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dy_k/32/50519_2.png) [@DY\_K](https://discourse.julialang.org/u/DY_K)\
**Post date:** [June 24, 2023, 6:22am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/1 "2023-06-24T06:22:47Z")

</div>

My simple test indicates that Julia is significantly slower than Matlab in vector and matrix algebra.

For Julia, I constructed a simple function like

function vec\_prod(N,n)

```
Arg1 = 0.0;
matt1 = Array{Float64}(undef, n, n);
matt2 = Array{Float64}(undef, n, n);
for i = 1:N
    vec1 = 1:n;
    mat = vec1*vec1';
    
    matt1 = (mat.*mat)/n^2;
    matt2 = (mat*mat)/n^2;
    Arg1 = Arg1 + sum(matt1.*matt2)/N;
end

return Arg1

```

end

In my computer, @time vec\_prod(10, 1000) gives ~1.8 sec.  
However, the same Matlab code takes only 0.1x sec.

I am very confused because Julia is known as a much higher-performance language than Matlab but this result is totally opposite.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 24, 2023, 6:28am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/2 "2023-06-24T06:28:12Z")

</div>

Hi @DY_K, welcome!

Indeed, Julia is fast but there are a few tricks to make it so. They are [summarized in the documentation](https://docs.julialang.org/en/v1/manual/performance-tips/), but at first glance here are the ones that might apply here:

- Don’t benchmark the first run of a function with `@time`, because it includes compilation time. Benchmark the second run instead, or even better, use `@btime` from [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl)
- Try to reduce memory allocations. For instance you re-allocate new matrices `matt1` and `matt2` at every iteration instead of writing into the containers you created at the beginning.

To get a more precise diagnosis, the first thing you can do is [profile your code](https://docs.julialang.org/en/v1/manual/profile/). If you’re in VSCode, it’s [very simple](https://www.julia-vscode.org/docs/dev/userguide/profiler/).

Side note: you don’t need to end every line with `;` in Julia 🙂

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 24, 2023, 6:44am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/3 "2023-06-24T06:44:55Z")

</div>

Of course your code is a bit weird because it does the same thing at every iteration, so we could just do it once and multiply the resulting `Arg1` by `N`. For the sake of comparison however, I wrote another function that performs exactly the same computations.

The better memory management is enabled by [pre-allocating outputs](https://docs.julialang.org/en/v1/manual/performance-tips/#Pre-allocating-outputs) and [fusing vectorized operations](https://docs.julialang.org/en/v1/manual/performance-tips/#More-dots:-Fuse-vectorized-operations).  
In addition I used a few specific linear algebra functions:

- [`dot(matt1, matt2)`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.dot) to replace `sum(matt1 .* matt2)` (which allocates the elementwise product before summing it)
- [`mul!(matt2, mat, mat)`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.mul!) instead of `matt2 = (mat*mat)` to avoid allocating the matrix product since I already know where I’m gonna store it

```julia
julia> using BenchmarkTools, LinearAlgebra

julia> function better_vec_prod_aux!(mat, matt1, matt2, N, n)
           Arg1 = 0.0
           for i = 1:N
               vec1 = 1:n
               mat .= vec1 .* vec1'
               matt1 .= (mat .* mat) .* (1 / n^2)
               mul!(matt2, mat, mat)
               matt2 .*= 1 / n^2
               Arg1 = Arg1 + dot(matt1, matt2) / N
           end
           return Arg1
       end;

julia> function better_vec_prod(N, n)
           mat = Matrix{Float64}(undef, n, n)
           matt1 = Matrix{Float64}(undef, n, n)
           matt2 = Matrix{Float64}(undef, n, n)
           return better_vec_prod_aux!(mat, matt1, matt2, N, n)
       end;

julia> vec_prod(2, 5) == better_vec_prod(2, 5)
true

julia> @btime vec_prod(10, 1000);
  3.820 s (154 allocations: 473.33 MiB)

julia> @btime better_vec_prod(10, 1000);
  238.816 ms (6 allocations: 22.89 MiB)

julia> @btime better_vec_prod_aux!(mat, matt1, matt2, 10, 1000) setup = (
           n = 1000;
           mat = Matrix{Float64}(undef, n, n);
           matt1 = Matrix{Float64}(undef, n, n);
           matt2 = Matrix{Float64}(undef, n, n)
       );
  247.210 ms (0 allocations: 0 bytes)

```

As you can see, in this version, the only allocations happen at the beginning of the loop: `better_vec_prod_aux!` performs none. This is especially important when dealing with large matrices like you are (1000 x 1000).

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 24, 2023, 7:02am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/5 "2023-06-24T07:02:13Z")

</div>

Yup, just noticed it, thanks. I thought I might be able to sweep this one under the rug, but more honest eyes were watching 🕶

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [June 24, 2023, 7:20am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/6 "2023-06-24T07:20:58Z")

</div>

Maybe in addition to @gdalle 's reply:  
It’s worth keeping in mind that MATLAB (as the name suggests) is great at doing matrix operations. Your example has a lot of large matrices, therefore, the runtime is probably dominated by the matrix operations (which likely use some `C`/`BLAS` routines underneath). A good Julia implementation will probably be as fast as MATLAB, but not much faster _in this concrete example_.

Julia has the strength that one can implement custom routines (with `for` loops) also in a high-performant way. Therefore, one can typically be faster then MATLAB, in cases where the MATLAB code cannot be fully vectorized.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 24, 2023, 7:59am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/7 "2023-06-24T07:59:49Z")

</div>

Whereas it is universally acknowledged that “Julia” stands for “Just Using Loops Is Amazing”

---

<div class="post-metadata">

**Author:** ![DY\_K](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dy_k/32/50519_2.png) [@DY\_K](https://discourse.julialang.org/u/DY_K)\
**Post date:** [June 24, 2023, 8:12am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/8 "2023-06-24T08:12:40Z")

</div>

Thank you for your replies.

It would not be an equal comparison by adding other input parameters (like mat, matt1, matt2) or constructing the second function (like better\_vec\_prod.jl).

However, from @gdalle’s suggestion, I’ve changed my code a little bit like this,  
function vec\_prod(N,n)

```
Arg1 = 0.0;
mat = Array{Float64}(undef, n, n)
matt1 = mat
matt2 = mat

for i = 1:N
    vec1 = 1:n
    mat .= vec1.*vec1'
    
    matt1 .= (mat.*mat)./n^2
    matt2 .= (mat*mat)./n^2
    Arg1 = Arg1 + sum(matt1.*matt2)/N
end

return Arg1

```

end

Now, @time vec\_prod(10, 1000) gives 0.11 ~0.13 sec which is almost the same as Matlab.

By the way, using dot and mul! of LinearAlgebra

sum(matt1.\* matt2) =\> dot(matt1, matt2),  
matt2 .= (mat \* mat)./n^2 =\> mul!(matt2, mat, mat); matt2 = matt2./n^2

makes the code faster. (~11% faster than Matlab, I checked average values from 10 results)  
Matlab is famous for very fast vector and matrix calculation. This is amazing!!

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 24, 2023, 9:06am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/9 "2023-06-24T09:06:43Z")

</div>

> [@DY\_K](#):
>
> It would not be an equal comparison by adding other input parameters (like mat, matt1, matt2) or constructing the second function (like better\_vec\_prod.jl).

I’m curious why you think that. I have constructed an auxiliary function for the purposes of clarity, but you could just as well initialize _and_ modify your matrices in `better_vec_prod`. There’s nothing that makes it unfair here.

---

<div class="post-metadata">

**Author:** ![DY\_K](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dy_k/32/50519_2.png) [@DY\_K](https://discourse.julialang.org/u/DY_K)\
**Post date:** [June 24, 2023, 9:12am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/11 "2023-06-24T09:12:11Z")

</div>

Never mind. From your comment, I’ve learned a lot. Thanks.  
Julia is amazing. I never thought Julia is faster than Matlab for this kind of simple calculation of large vectors and matrice.

---

<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 24, 2023, 10:14am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/12 "2023-06-24T10:14:47Z")

</div>

> [@DY\_K](#):
>
> I never thought Julia is faster than Matlab for this kind of simple calculation of large vectors and matrice.

It really shouldn’t: if you only do linear algebra operations like matrix-matrix multiplications, both julia and matlab eventually call a BLAS library (OpenBLAS for Julia by default, probably MKL in Matlab, which you can use in Julia with [`MKL.jl`](https://github.com/JuliaLinearAlgebra/MKL.jl)), the language you’re writing the code in is totally irrelevant. You aren’t really measuring any language-specific feature, apart from the cost of allocating lots of unneeded memory.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [June 24, 2023, 3:22pm UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/13 "2023-06-24T15:22:23Z")

</div>

It could be: it depends on what BLAS functions are exposed. For example, is it possible to perform inplace operations in matlab?

---

<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 24, 2023, 3:28pm UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/14 "2023-06-24T15:28:53Z")

</div>

But then you aren’t comparing apples-to-apples.

---

<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:** [June 24, 2023, 3:50pm UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/15 "2023-06-24T15:50:03Z")

</div>

Matlab normally uses multi-threading for many operations (not just BLAS). Check that you are not comparing multi-threaded Matlab operations with single-threaded Julia ops.

---

<div class="post-metadata">

**Author:** ![alfaromartino](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alfaromartino/32/52986_2.png) [@alfaromartino](https://discourse.julialang.org/u/alfaromartino)\
**Post date:** [June 24, 2023, 5:49pm UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/16 "2023-06-24T17:49:42Z")

</div>

I emphasize the point made by @DNF. Matlab is multi-threaded by default. In my computer with 16 cores and using the original code snippet (so not improving the code by @gdalle’s suggestions):

2.159 s (154 allocations: 473.33 MiB) → original code  
304.525 ms (299 allocations: 473.34 MiB) → threaded code

Just in case, the basic threaded option only requires adding `Threads.@threads` before the for-loop

```julia
function vec_prod2(N,n)
  Arg1 = 0.0;
  matt1 = Array{Float64}(undef, n, n);
  matt2 = Array{Float64}(undef, n, n);
  Threads.@threads for i = 1:N
      vec1 = 1:n;
      mat = vec1*vec1';
      
      matt1 = (mat.*mat)/n^2;
      matt2 = (mat*mat)/n^2;
      Arg1 = Arg1 + sum(matt1.*matt2)/N;
  end

  return Arg1

end

```

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [June 24, 2023, 7:43pm UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/17 "2023-06-24T19:43:19Z")

</div>

@alfaromartino that code has a race condition on updating `Arg1`, it’s completely unsafe and incorrect.

---

<div class="post-metadata">

**Author:** ![Paul\_Soderlind](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_soderlind/32/1753_2.png) [@Paul\_Soderlind](https://discourse.julialang.org/u/Paul_Soderlind)\
**Post date:** [June 24, 2023, 7:47pm UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/18 "2023-06-24T19:47:25Z")

</div>

Yes, create a vector `Arg1` and fill as `Arg1[i]` and then after the loop, sum over the vector. That should work.

---

<div class="post-metadata">

**Author:** ![alfaromartino](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alfaromartino/32/52986_2.png) [@alfaromartino](https://discourse.julialang.org/u/alfaromartino)\
**Post date:** [June 24, 2023, 7:48pm UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/19 "2023-06-24T19:48:32Z")

</div>

oh I hadn’t checked the code. Yes, you’re definitely right. The point in any case is that matlab is multi-threaded by default.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [June 24, 2023, 7:53pm UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/20 "2023-06-24T19:53:44Z")

</div>

> [@DY\_K](#):
>
> However, from @gdalle’s suggestion, I’ve changed my code a little bit like this,  
> function vec\_prod(N,n)
> 
> ```julia
> Arg1 = 0.0;
> mat = Array{Float64}(undef, n, n)
> matt1 = mat
> matt2 = mat
> 
> ```

by the way, writing this does not copy `mat`, all three of those variables are pointing to the same memory so this code is also incorrect

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [June 24, 2023, 7:54pm UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/21 "2023-06-24T19:54:57Z")

</div>

> [@Paul\_Soderlind](#):
>
> Yes, create a vector `Arg1` and fill as `Arg1[i]` and then after the loop, sum over the vector. That should work.

No, this pattern is also incorrect because tasks can yield partway through and overwrite eachothers’ work.

---

<div class="post-metadata">

**Author:** ![alfaromartino](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alfaromartino/32/52986_2.png) [@alfaromartino](https://discourse.julialang.org/u/alfaromartino)\
**Post date:** [June 24, 2023, 8:03pm UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/22 "2023-06-24T20:03:05Z")

</div>

Not sure I understand this. What @Paul_Soderlind suggests corrects the race condition.

```julia

function vec_prod1(N,n)
  matt1 = Array{Float64}(undef, n, n);
  matt2 = Array{Float64}(undef, n, n);
  Arg1 = Vector{Float64}(undef,N}

  for i = 1:N
      vec1 = 1:n;
      mat = vec1*vec1';
      
      matt1 = (mat.*mat)/n^2;
      matt2 = (mat*mat)/n^2;
      Arg1[i] = sum(matt1.*matt2)/N;
  end

  return sum(Arg1)

end

function vec_prod2(N,n)  
  matt1 = Array{Float64}(undef, n, n);
  matt2 = Array{Float64}(undef, n, n);
  Arg1 = Vector{Float64}(undef,N}

  Threads.@threads for i = 1:N
      vec1 = 1:n;
      mat = vec1*vec1';
      
      matt1 = (mat.*mat)/n^2;
      matt2 = (mat*mat)/n^2;
      Arg1[i] = sum(matt1.*matt2)/N;
  end

  return sum(Arg1)

end

```

```julia
vec_prod1(10, 1000) == vec_prod2(10,1000) #true

@btime vec_prod2(10,1000) # 2.159 s (154 allocations: 473.33 MiB)
@btime vec_prod1(10, 1000) # 297.634 ms (289 allocations: 473.34 MiB)

```

btw, nice to know that we work at the same university @Mason!!!

[Next page](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773.md?page=2)