# Performance issue with multithreaded computation with matrix operations at its heart (Threads.@threads vs. BLAS threads)

**URL:** <https://discourse.julialang.org/t/performance-issue-with-multithreaded-computation-with-matrix-operations-at-its-heart-threads-threads-vs-blas-threads/106043>\
**Category:** Performance\
**Tags:** blas, parallel, multithreading, linearalgebra, threads\
**Created:** [November 10, 2023, 4:20pm UTC](https://discourse.julialang.org/t/performance-issue-with-multithreaded-computation-with-matrix-operations-at-its-heart-threads-threads-vs-blas-threads/106043 "2023-11-10T16:20:16Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![HanD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hand/32/213908_2.png) [@HanD](https://discourse.julialang.org/u/HanD)\
**Post date:** [November 10, 2023, 4:20pm UTC](https://discourse.julialang.org/t/performance-issue-with-multithreaded-computation-with-matrix-operations-at-its-heart-threads-threads-vs-blas-threads/106043/1 "2023-11-10T16:20:16Z")

</div>

While trying to speed up a long computation with chunking more threads on it, I ran into an interesting problem. I found several related threads here, but I couldn’t find a convincing answer, nor a solution to my issue.

Here’s an MVP of my issue. Consider the following definitions:

```julia-repl
julia> using LinearAlgebra, BenchmarkTools
julia> ts = 1:10; xs = rand(1000, 10);
julia> function fit(x, y, degree)
           return qr([x_n ^ k for x_n in x, k in 0:degree]) \ y
       end

```

And then I run the following measurements, with `JULIA_NUM_THREADS=8`:

```julia-repl
julia> @benchmark (for row in $(collect(eachrow(xs))); fit(ts, row, 3); end)
BenchmarkTools.Trial: 1156 samples with 1 evaluation.
 Range (min … max): 4.181 ms … 6.373 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 4.228 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 4.326 ms ± 284.965 μs ┊ GC (mean ± σ): 0.76% ± 3.36%

  ▄██▃▁ ▁                                       
  █████▆▆▅▆▅▅▅▄▅▅▄▄▄▅▁▆██▆▆▅▆█▅▄▄▁▅▁▄▅▁▄▄▁▅▄▄▅▅▇▇▆▆▅▄▅▁▅▁▅▁▄▅ █
  4.18 ms Histogram: log(frequency) by time 5.48 ms <

 Memory estimate: 1.50 MiB, allocs estimate: 8000.

julia> @benchmark @threads (for row in $(collect(eachrow(xs))); fit(ts, row, 3); end)
BenchmarkTools.Trial: 314 samples with 1 evaluation.
 Range (min … max): 12.338 ms … 69.561 ms ┊ GC (min … max): 0.00% … 77.62%
 Time (median): 15.782 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 15.939 ms ± 3.073 ms ┊ GC (mean ± σ): 1.08% ± 4.38%

                                        ▂▂▃█ ▁▄▃▆▄             
  ▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▁▁▃▃▂▁▂▃▁▃▃▅▃▆██████████▅▂▃▁▁▂▁▁▁▁▂ ▃
  12.3 ms Histogram: frequency by time 17.1 ms <

 Memory estimate: 1.50 MiB, allocs estimate: 8049.

```

Notice how the execution time went up, even though the computation runs on 8 threads instead of 1. If I observe CPU usage in `htop`, it is obvious that most of the resources are wasted in system calls, most of the bars are red. I found suggestions in various threads that calling `BLAS.set_num_threads(1)` could improve the situation, but in my case, it has no visible effect, I get identical results.

I’m guessing that the Julia threads compete for the LAPACK library calls, which thus constitutes a bottleneck, but I don’t know how to get around this issue. Ideally, the (more elaborate) computation would be running on a 64 core CPU, which is currently sitting mostly idle, because I can only run this on a single thread.

Any ideas?

---

<div class="post-metadata">

**Author:** ![HanD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hand/32/213908_2.png) [@HanD](https://discourse.julialang.org/u/HanD)\
**Post date:** [November 10, 2023, 4:30pm UTC](https://discourse.julialang.org/t/performance-issue-with-multithreaded-computation-with-matrix-operations-at-its-heart-threads-threads-vs-blas-threads/106043/2 "2023-11-10T16:30:00Z")

</div>

Okay, I found something interesting. If I remove the `qr` decomposition from `fit`, the single-threaded solution slows down by a factor of 2, but the multi-threaded solution suddenly becomes a lot more efficient:

```julia-repl
julia> function fit(x, y, degree)
           return [x_n ^ k for x_n in x, k in 0:degree] \ y
       end

julia> @benchmark (for row in $(collect(eachrow(xs))); fit(ts, row, 3); end)
BenchmarkTools.Trial: 568 samples with 1 evaluation.
 Range (min … max): 7.684 ms … 11.128 ms ┊ GC (min … max): 15.37% … 14.02%
 Time (median): 8.694 ms ┊ GC (median): 15.73%
 Time (mean ± σ): 8.803 ms ± 628.737 μs ┊ GC (mean ± σ): 20.52% ± 5.92%

         ▂▄█▆▅▄▆█▆▄ ▂ ▁▃▄▆▅▄ ▄ ▂                        
  ▄▃▄▄▄▅████████████▄▄▆▃▄█▆▆▇██████▇███▄▅▅▄▄▂▆▄▄▄▂▂▂▁▁▂▁▂▁▁▃▃ ▄
  7.68 ms Histogram: frequency by time 10.6 ms <

 Memory estimate: 69.03 MiB, allocs estimate: 45000.

julia> @benchmark @threads (for row in $(collect(eachrow(xs))); fit(ts, row, 3); end)
BenchmarkTools.Trial: 2620 samples with 1 evaluation.
 Range (min … max): 1.028 ms … 24.744 ms ┊ GC (min … max): 0.00% … 11.35%
 Time (median): 1.375 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 1.903 ms ± 1.430 ms ┊ GC (mean ± σ): 25.81% ± 25.75%

   ▂▄▅▅██▆▂▁▁ ▁▁ ▁▄▅▃ ▁
  ▇██████████▇▇▇▆▇▇▆▅▄▃▃▁▁▄▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▅█████████▇▆▅ █
  1.03 ms Histogram: log(frequency) by time 4.3 ms <

 Memory estimate: 69.03 MiB, allocs estimate: 45049.

```

So what is going on with `qr` that makes it particularly inefficient when invoked from multiple threads in parallel?

---

<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:** [November 10, 2023, 4:48pm UTC](https://discourse.julialang.org/t/performance-issue-with-multithreaded-computation-with-matrix-operations-at-its-heart-threads-threads-vs-blas-threads/106043/3 "2023-11-10T16:48:34Z")

</div>

I have several hypotheses here:

- maybe it’s a nasty interaction between Julia threads and BLAS threads, as you mentioned =\> there is now a [section in the performance tips about this](https://docs.julialang.org/en/v1.10.0-rc1/manual/performance-tips/#man-multithreading-linear-algebra)
- maybe it’s the additional allocations triggered by `qr` which render multithreading inefficient
- maybe it’s linked to you benchmarking with globals
- maybe it’s the macros `@benchmark` and `@threads` which don’t play nice together (unlikely)

---

<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:** [November 10, 2023, 4:55pm UTC](https://discourse.julialang.org/t/performance-issue-with-multithreaded-computation-with-matrix-operations-at-its-heart-threads-threads-vs-blas-threads/106043/4 "2023-11-10T16:55:53Z")

</div>

You’re matrices are (too) tiny.

```julia
function serial(mat, vec)
    for i in 1:1000
        qr(mat) \ vec
    end
end

function threaded(mat, vec)
    @threads for i in 1:1000
        qr(mat) \ vec
    end
end

mat = rand(10,4)
vec = rand(10)

@btime serial($mat, $vec) # 6.052 ms (7000 allocations: 1.10 MiB)
@btime threaded($mat, $vec) # 17.643 ms (7031 allocations: 1.10 MiB)

mat2 = rand(1000,100)
vec2 = rand(1000)

@btime serial($mat2, $vec2) # 2.804 s (10000 allocations: 827.50 MiB)
@btime threaded($mat2, $vec2) # 1.024 s (10031 allocations: 827.50 MiB)

```

(with 6 threads)

---

<div class="post-metadata">

**Author:** ![HanD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hand/32/213908_2.png) [@HanD](https://discourse.julialang.org/u/HanD)\
**Post date:** [November 10, 2023, 5:02pm UTC](https://discourse.julialang.org/t/performance-issue-with-multithreaded-computation-with-matrix-operations-at-its-heart-threads-threads-vs-blas-threads/106043/5 "2023-11-10T17:02:41Z")

</div>

> [@carstenbauer](#):
>
> You’re matrices are (too) tiny.

That may be so, but what can I do if I need to least square fit 3rd degree polynomials on 10 points. These are the matrices I have to work with.

As a workaround, I can remove the qr transformation, yielding the same result.

---

<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:** [November 10, 2023, 5:03pm UTC](https://discourse.julialang.org/t/performance-issue-with-multithreaded-computation-with-matrix-operations-at-its-heart-threads-threads-vs-blas-threads/106043/6 "2023-11-10T17:03:07Z")

</div>

> [@HanD](#):
>
> That may be so, but what can I do if I need to least square fit 3rd degree polynomials on 10 points. These are the matrices I have to work with.

Use StaticArrays.jl?

---

<div class="post-metadata">

**Author:** ![HanD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hand/32/213908_2.png) [@HanD](https://discourse.julialang.org/u/HanD)\
**Post date:** [November 13, 2023, 7:45am UTC](https://discourse.julialang.org/t/performance-issue-with-multithreaded-computation-with-matrix-operations-at-its-heart-threads-threads-vs-blas-threads/106043/7 "2023-11-13T07:45:25Z")

</div>

> [@gdalle](#):
>
> > [@HanD](#):
> >
> > That may be so, but what can I do if I need to least square fit 3rd degree polynomials on 10 points. These are the matrices I have to work with.
> 
> Use StaticArrays.jl?

I don’t see how that would help here, or even how it would work. For one thing, the size of the Vandermonde matrix above isn’t known at compile time, since the number of points to fit and the degree of the desired polynomial is only known at runtime. Also, the StaticArrays.jl manual suggests to use static arrays only below 100 items, and in this case, I can easily have more than that. The 3 and 10 are just specific examples.

---

<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:** [November 13, 2023, 5:00pm UTC](https://discourse.julialang.org/t/performance-issue-with-multithreaded-computation-with-matrix-operations-at-its-heart-threads-threads-vs-blas-threads/106043/8 "2023-11-13T17:00:05Z")

</div>

> [@HanD](#):
>
> The 3 and 10 are just specific examples.

My bad, when you said “these are the matrices I have to work with” I thought they would always look like that