# Having issues speeding up code with multithreading

**URL:** <https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508>\
**Category:** Performance\
**Tags:** parallel, multithreading\
**Created:** [July 12, 2023, 5:28am UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508 "2023-07-12T05:28:39Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![RobertGregg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robertgregg/32/22105_2.png) [@RobertGregg](https://discourse.julialang.org/u/RobertGregg)\
**Post date:** [July 12, 2023, 5:28am UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/1 "2023-07-12T05:28:39Z")

</div>

I’m trying to speed up some code that takes subsets of columns from a data matrix and calculates a maximum score. The orginal code is complicated, so I’ve tried to simplfy it was much as possible here:

```julia
using Combinatorics

#Helper function to loop through pairs of data columns
allpairs(v) = Iterators.filter(i -> isless(i...), Iterators.product(v,v))

function maxScore(data)

    bestScore = 0.0

    for (i,j) in allpairs(axes(data,2))
        
        if isodd(i+j)
            currentScore = findscore(data,i,j)

            if bestScore < currentScore
                bestScore = currentScore
            end
        end
    end

    return bestScore
end

function findscore(data, i, j)
    
    bestScoreSubset = 0.0

    #Create a somewhat small subset
    subsets = mod1(i+j,10)

    for subset in powerset(1:subsets)

        currentScoreSubset = score(data,subset,j)

        if bestScoreSubset < currentScoreSubset
            bestScoreSubset = currentScoreSubset
        end

    end

    return bestScoreSubset
end

function score(data,subset,j)

    if isempty(subset)
        return 0.0
    end
    
    X = view(data,:,subset)
    y = view(data,:,j)

    b = X \ y

    ŷ = X*b

    return sum((yᵢ - ŷᵢ)^2 for (yᵢ, ŷᵢ) in zip(y,ŷ))
end

```

I tried using `ThreadsX` to parallelize the `maxScore()` function, but it actually made things slower.

> **ThreadsX version of maxScore**
>
> ```julia
> using ThreadsX
> 
> function maxScore_parallel(data)
> 
> bestScore = ThreadsX.mapreduce(max, allpairs(axes(data,2))) do (i,j)
> if isodd(i+j)
> findscore(data,i,j)
> else
> 0.0
> end
> end
> 
> return bestScore
> end
> 
> ```

Here are some benchmarks as well

> **Benchmarks**
>
> ```julia
> julia> using Random, BenchmarkTools
> 
> julia> Threads.nthreads()
> 14
> 
> julia> Random.seed!(1);
> 
> julia> data = rand(100,20);
> 
> julia> @btime maxScore($data)
> 174.879 ms (612340 allocations: 990.52 MiB)
> 19.720119456710123
> 
> julia> @btime maxScore_parallel($data)
> 906.113 ms (622461 allocations: 990.90 MiB)
> 19.720119456710123
> 
> ```

Any recommendations for improving parallelization would be much appreciated 🙂

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [July 12, 2023, 6:39am UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/2 "2023-07-12T06:39:11Z")

</div>

My 2 cents : you should try to remove all allocations from your computation before its parallelization.  
For example, I guess that this operation

> [@RobertGregg](#):
>
> ` b = X \ y`

allocates.

---

<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:** [July 12, 2023, 7:26am UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/3 "2023-07-12T07:26:38Z")

</div>

I wonder OP might be better off doing `X \ Y`, solving for all the columns at once.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [July 12, 2023, 8:00am UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/4 "2023-07-12T08:00:24Z")

</div>

Actually I did not try to understand the code.

My experience is that it is really difficult to obtain interesting Speed-Ups with allocating tasks.

---

<div class="post-metadata">

**Author:** ![RobertGregg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robertgregg/32/22105_2.png) [@RobertGregg](https://discourse.julialang.org/u/RobertGregg)\
**Post date:** [July 12, 2023, 2:40pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/5 "2023-07-12T14:40:55Z")

</div>

I definitely agree that avoiding allocations would helpful. In this case however I think that would be difficult because X changes size with every iteration, and both X and y change values. In the original code this size change is unpredictable.

Maybe there’s a good way allocate space for something variable in size but with a known maximum size?

I think I would also have to worry about race conditions with the allocated space.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [July 12, 2023, 2:51pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/6 "2023-07-12T14:51:08Z")

</div>

If you have an idea of the max size you may allocate for the max size. This workspace should be private for each task… Not trivial.

If the max matrix size is small (less than 100x100) then StaticArrays could be an interesting option.

---

<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:** [July 12, 2023, 2:57pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/7 "2023-07-12T14:57:35Z")

</div>

> [@RobertGregg](#):
>
> I think I would also have to worry about race conditions with the allocated space.

> **[PSA: Thread-local state is no longer recommended](https://julialang.org/blog/2023/07/PSA-dont-use-threadid/)**
>
> PSA: Thread-local state is no longer recommended; Common misconceptions about threadid() and nthreads()

---

<div class="post-metadata">

**Author:** ![RobertGregg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robertgregg/32/22105_2.png) [@RobertGregg](https://discourse.julialang.org/u/RobertGregg)\
**Post date:** [July 12, 2023, 3:21pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/8 "2023-07-12T15:21:30Z")

</div>

I think understand what you’re proposing but I’m not 100% sure. Do you mean solving `X \ y` for every column subset and saving the score? That would create a 2^N by N matrix which is probably only feasible for small problems.

Do you think the linear solve is causing `ThreadsX.mapreduce` to be slower?

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [July 12, 2023, 3:37pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/9 "2023-07-12T15:37:13Z")

</div>

FWIW your code exhibits quite different perf on my laptop (M1 Max). You could try to reduce the number of threads on your machine.

```julia
julia> include("src/titi.jl")
Threads.nthreads() = 8
  101.412 ms (639420 allocations: 1.02 GiB)  
  45.493 ms (642679 allocations: 1.02 GiB) #ThreadsX
19.72011945671012

```

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [July 12, 2023, 3:44pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/10 "2023-07-12T15:44:39Z")

</div>

And all the time is spent in allocation required by the linear solve:

 ![Capture d’écran 2023-07-12 à 17.42.46](https://global.discourse-cdn.com/julialang/original/3X/d/5/d5ab93930cb8b353f8401725d5839d28175c0197.jpeg)

---

<div class="post-metadata">

**Author:** ![RobertGregg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robertgregg/32/22105_2.png) [@RobertGregg](https://discourse.julialang.org/u/RobertGregg)\
**Post date:** [July 12, 2023, 5:14pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/11 "2023-07-12T17:14:39Z")

</div>

Interesting! I lowered the number of threads to 8 and got a very similar result. I even tried incrementing the number of threads and it looks like going from 1 to 2 threads about halves the time, but after that you get about the same result. It looks like there’s a point where too many threads just tanks performance.

---

<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:** [July 12, 2023, 6:24pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/12 "2023-07-12T18:24:39Z")

</div>

> [@RobertGregg](#):
>
> Do you mean solving `X \ y` for every column subset and saving the score?

Actually I was wrong, I thought only the `y` changed but since `X` changes too it’s harder

> [@RobertGregg](#):
>
> Do you think the linear solve is causing `ThreadsX.mapreduce` to be slower?

Definitely!  
As the [profiling](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/10) by @LaurentPlagne shows you, most of the time is spent allocating memory (the yellow tiles). And multithreading works badly with memory-intensive codes.  
On the bright side, it means that you have a wide margin of improvement accessible by simply optimizing your sequential code, without even worrying about parallelism.  
And this improvement will in turn make multithreading more efficient.

> [@LaurentPlagne](#):
>
> And all the time is spent in allocation required by the linear solve:

I think your best bet for non-allocating linear solvers is this place:

> **[GitHub - SciML/LinearSolve.jl: LinearSolve.jl: High-Performance Unified...](https://github.com/SciML/LinearSolve.jl)**
>
> LinearSolve.jl: High-Performance Unified Interface for Linear Solvers in Julia. Easily switch between factorization and Krylov methods, add preconditioners, and all in one interface. - GitHub - Sci...

You will probably need to allocate once per each size of `X` though

---

<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:** [July 12, 2023, 8:47pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/13 "2023-07-12T20:47:38Z")

</div>

> [@gdalle](#):
>
> multithreading works badly with memory-intensive codes.

Why is this exactly? Is heap allocation even multithreadable, I’d expect it to happen one at a time to prevent writing to overlapped memory, but I don’t know the actual implementation of Julia’s allocator.

---

<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:** [July 12, 2023, 8:52pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/14 "2023-07-12T20:52:58Z")

</div>

There are two reasons. The first is that before Julia 1.10 (currently an alpha), GC is single threaded, so as you add more threads, you allocate faster and all the time ends up in the GC. The second problem is that for larger objects (bigger than a few kb), there is a lock around allocating since you end up with essentially calls to malloc for the allocation.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [July 12, 2023, 9:43pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/15 "2023-07-12T21:43:42Z")

</div>

It would presumably help to squeeze out all the allocations from your `score` function. I’ve been looking for an excuse to play with @Mason’s neat [Bumper.jl](https://github.com/MasonProtter/Bumper.jl) as an alternative to writing technical multithreaded code that manually manages local buffers correctly, and this seems like a good opportunity.

Anyways, here is a version of your `score` function that is _almost_ non-allocating:

```julia
function score(data,subset,j)
  isempty(subset) && return 0.0
  X = view(data,:,subset)
  y = view(data,:,j)
  @no_escape begin
    Xv = alloc(Float64, size(X,1), size(X,2))
    b = alloc(Float64, size(X,2))
    ŷ = alloc(Float64, length(y))
    copyto!(Xv, X)
    Xvf = qr!(Xv) # one alloc
    ldiv!(b, Xvf, y) # three allocs
    mul!(ŷ, X, b)
    sum((yᵢ - ŷᵢ)^2 for (yᵢ, ŷᵢ) in zip(y,ŷ))
  end
end

```

Emphasis on _almost_ though. I played a little with slightly more thoughtful methods for the `ldiv!` because it doesn’t really seem like the best way to solve that least squares problem. But methods like `ldiv!(Xvf.R, [...])` don’t exist because `Xvf` is a `PtrArray`. And that single alloc in the `qr!` call may also be hard to get rid of, although I bet somebody else knows how to do it. Maybe [GenericLinearAlgebra.jl](https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl) could help?

With that said, getting almost all the allocations out didn’t actually make the call to `score` itself any faster, and I didn’t bother trying with the threading since the function call still touches the heap a few times. But just passing this along in case somebody else knows how to get those out or has additional thoughts. And because it is pretty neat to seemingly get rid of all those allocations completely for free with `Bumper.jl`. I should think that getting all those heap allocations out will make your code scale better with threads.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [July 12, 2023, 9:49pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/16 "2023-07-12T21:49:18Z")

</div>

> [@Oscar\_Smith](#):
>
> The second problem is that for larger objects (bigger than a few kb), there is a lock around allocating since you end up with essentially calls to malloc for the allocation.

Oh no, the GAL (global allocation lock) has infected Julia 😉

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [July 12, 2023, 9:59pm UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/17 "2023-07-12T21:59:36Z")

</div>

I think that @Oscar_Smith is just referring to the `malloc`-internal synchronization.

---

<div class="post-metadata">

**Author:** ![RobertGregg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robertgregg/32/22105_2.png) [@RobertGregg](https://discourse.julialang.org/u/RobertGregg)\
**Post date:** [July 13, 2023, 1:31am UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/18 "2023-07-13T01:31:13Z")

</div>

Wow [Bumper.jl](https://github.com/MasonProtter/Bumper.jl) is pretty awesome! One thought I had was maybe combining it with [QRupdate.jl](https://github.com/mpf/QRupdate.jl) which allows you to take a qr decomposition of the full dataset and then remove columns without recomputing the entire factorization.

Interestingly, using the score method with bumper _ **did** _ significantly improve performance when multi-threading, despite the similar results you mentioned:

```julia
using Random, BenchmarkTools

Threads.nthreads() # 7
Random.seed!(1);
data = rand(100,20);

@btime maxScore($data) #129.537 ms (612340 allocations: 990.52 MiB)
@btime maxScore_parallel($data) #53.769 ms (617649 allocations: 990.71 MiB)
@btime maxScore_parallel_bump($data) #14.049 ms (209736 allocations: 64.15 MiB)

```

I’ll call a ~10x improvement good for now 😁

---

<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:** [July 13, 2023, 2:19am UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/19 "2023-07-13T02:19:29Z")

</div>

I’m really glad Bumper.jl was useful for you! Please do file issues if you run into bugs or have any feature requests.

---

<div class="post-metadata">

**Author:** ![RobertGregg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robertgregg/32/22105_2.png) [@RobertGregg](https://discourse.julialang.org/u/RobertGregg)\
**Post date:** [July 16, 2023, 6:14am UTC](https://discourse.julialang.org/t/having-issues-speeding-up-code-with-multithreading/101508/20 "2023-07-16T06:14:24Z")

</div>

**Update**. After learning way to much about QR decomposition, I wrote a replacement for `b = X \ y` that has zero allocations 🎉

```julia
using LinearAlgebra

#Loop through lower triangular matrix
#the 1st axis is reversed because we're moving non-zeros *up* above diagonal
lowTriIter(ax1, ax2) = Iterators.filter(i -> first(i)>last(i), Iterators.product(Iterators.reverse(ax1),ax2))
lowTriIter(A::AbstractMatrix) = lowTriIter(axes(A)...)

function qless!(b, X, y, R, Rsub, subset)

    numCol = length(subset)

    #Copy columns over to avoid destroying R
    for (k,j) in enumerate(subset)
        @views Rsub[:,k] = R[:,j]
    end

    #Use givens rotations to zero out values below diagonal
    #Loop through lower triangular portion of Rsub
    @inbounds for (i,j) in lowTriIter(Rsub)        
        if Rsub[i,j] ≠ 0.0
            G,r = givens(Rsub,i-1,i,j)
            Rsub[i-1,j] = r
            Rsub[i,j] = 0.0

            for k in j+1:numCol
                (r1, r2) = Rsub[i-1,k], Rsub[i,k]
                Rsub[i-1,k] = G.c*r1 + G.s*r2
                Rsub[i,k] = -G.s*r1 + G.c*r2
            end
        end
    end

    #Needed to call correct ldiv! method
    Rv = UpperTriangular(view(Rsub,1:numCol,1:numCol))

    #Solve R'R*b = X'y
    mul!(b,X',y)
    ldiv!(Rv',b)
    ldiv!(Rv,b)

    return nothing
end

```

b, X, y, and subset are the same as before. R is upper-triangular matrix resulting from `Q,R=qr(data)`. Rsub is a matrix where we can copy columns from R based on the subset.

```julia
julia> @btime $X\$y
  6.350 μs (41 allocations: 73.56 KiB)
4-element Vector{Float64}:
 0.11971500332311431
 0.22251862677083714
 0.3953695790986452
 0.12839479252581568

julia> @btime qless!($b,$X,$y,$R,$Rsub,$subset)
  2.133 μs (0 allocations: 0 bytes)

julia> b
4-element Vector{Float64}:
 0.11971500332311383
 0.22251862677083722
 0.3953695790986456
 0.12839479252581615

```