# 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:** 13\
**Page:** 2

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

</div>

Because `matt1` and `matt2` are (without purpose) defined before the loop? Otherwise, I don’t see it. (To be safe, declare `vec1`, `mat`, `matt1` and `matt2` as locals.)

---

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

</div>

Working on a blog-post about it, but here’s the unpublished draft: [www.julialang.org/blog/2023/06/PSA-dont-use-threadid.md at 65e25e87d01111cc45a477c8072c9f5dc5878a39 · JuliaLang/www.julialang.org · GitHub](https://github.com/JuliaLang/www.julialang.org/blob/65e25e87d01111cc45a477c8072c9f5dc5878a39/blog/2023/06/PSA-dont-use-threadid.md)

---

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

</div>

Interesting. I always thought that each task in a for-loop was running separately, in the sense of not sharing local variables with the other iterations. Is that the problem in last instance?

---

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

</div>

Yes, they run separately. The problem is that to do a summation, you need to somehow merge those separate values, and the various ways people often try to merge those values is incorrect.

---

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

</div>

> [@alfaromartino](#):
>
> ```julia
> 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
> 
> ```

But is the problem you mention arising here? if each thread is running separately and you store the result in `Arg1[i]`, then you can sum `Arg1` out of the loop.

---

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

</div>

Yes, that one is fine because you haven’t used `threadid`. But if you want to be efficient, you should avoid allocating an `N` element vector `Arg1` and can get away with `nthreads()` different elements, and that’s where people start getting themselves into trouble.

---

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

</div>

To get things back on topic, here’s the fastest version I could make (essentially just what @gdalle already wrote):

```julia
#+begin_src julia
function vec_prod1(N,n)
    Arg = 0.0
    mat = Array{Float64}(undef, n, n)
    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 .= mul!(matt2, mat, mat) ./ n^2
        Arg += (matt1 ⋅ matt2)/N;
    end
    Arg
end

@btime vec_prod1(10, 1000)
#+end_src

```

```julia
#+RESULTS:
: 84.133 ms (6 allocations: 22.89 MiB)
: 2.094817739604174e19

```

I found this was actually faster than using explicit multithreading since the matrix multiplication is already multithreaded for `1000 x 1000` matrices and BLAS multithreading is very efficient.

Here’s the explicitly multithreaded version for comparison

```julia
#+begin_src julia
function vec_prod3(N,n)
    chunks = Iterators.partition(1:N, max(1, N ÷ Threads.nthreads()))
    tasks = map(chunks) do chunk
        Arg = 0.0
        mat = Array{Float64}(undef, n, n)
        matt1 = Array{Float64}(undef, n, n)
        matt2 = Array{Float64}(undef, n, n)
        for i ∈ chunk
            vec1 = 1:n;
            mat .= vec1 .* vec1'
            matt1 .= (mat .* mat) ./ n^2;
            matt2 .= mul!(matt2, mat, mat) ./ n^2
            Arg += (matt1 ⋅ matt2)/N;
        end
        return Arg
    end
    Arg = sum(fetch, tasks)
end

@btime vec_prod3(10, 1000)
#+end_src

```

```julia
#+RESULTS:
: 97.034 ms (84 allocations: 228.88 MiB)
: 2.094817739604174e19

```

---

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

</div>

> [@alfaromartino](#):
>
> ```julia
> 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
> 
> ```

Sorry, I dont know why I said this was fine. You’ll still get race conditions from different threads overwriting `mat`, `matt1` and `matt2` (assuming you were actually calculating new values each iteration instead of reusing the exact same values like is done here)

---

<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 25, 2023, 1:57am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/31 "2023-06-25T01:57:18Z")

</div>

Now I see what you mean. The variables inside the loop weren’t local, because they were already defined outside the loop.

Yeah, redefining variables inside a function is prone to err. For example for type stability:

```julia
#STABLE
function type_stable()
    e = 1
    parameter() = e

    return e 
end

#UNSTABLE
function type_unstable()
    e = 1    
    parameter() = e
    
    e = 1
    return e
end

#UNSTABLE
function type_unstable()
    e = 1    
    parameter() = e
    
    e::Int64 = 1
    return e
end

```

---

<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:** [June 25, 2023, 1:59am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/32 "2023-06-25T01:59:29Z")

</div>

yeah, this is my least favorite issue in Julia: [https://github.com/JuliaLang/julia/issues/15276](https://github.com/JuliaLang/julia/issues/15276). It is one of the most annoying performance issues in Julia.

---

<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 25, 2023, 2:07am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/33 "2023-06-25T02:07:42Z")

</div>

> [@alfaromartino](#):
>
> Yeah, redefining variables inside a function is prone to err. For example for type stability:

Type instability isnt the problem here, it’s a _race condition_ [Race condition - Wikipedia](https://en.m.wikipedia.org/wiki/Race_condition)

---

<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 25, 2023, 2:33am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/34 "2023-06-25T02:33:52Z")

</div>

I know, I meant that redefining variables can lead to several problems (e.g. race conditions, type instability)

---

<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 25, 2023, 2:39am UTC](https://discourse.julialang.org/t/julia-is-significantly-slower-18-x-than-matlab-in-vector-and-matrix-algebra/100773/35 "2023-06-25T02:39:08Z")

</div>

In this case, the variable is never re-assigned, it’s just that the memory of the array it points to is over-written

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