# Any suggestion for speeding up my function?

**URL:** <https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393>\
**Category:** Performance\
**Tags:** question\
**Created:** [November 17, 2022, 9:56am UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393 "2022-11-17T09:56:10Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![tiZ](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiz/32/211800_2.png) [@tiZ](https://discourse.julialang.org/u/tiZ)\
**Post date:** [November 17, 2022, 9:56am UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/1 "2022-11-17T09:56:11Z")

</div>

Hi,

I have a function accepts 4 matrixs and do some caculations, it works well but very slow. I wonder to know how to improve its performance?

I’m wondering if it’s possible to be faster, I just know that using `@. ` the speed will be faster.

Here is a demo for my function, `Θ`,`M`,`Π`,`X` are all matrix, with same sizes, while `X` is Matrix{Int},others are Matrix{Float32}

```julia
using SpecialFunctions
γ = loggamma # let the γ represents the loggamma() function

function myFunc(Θ,M,Π,X)
    eps_ = 1f-10
    
    A = @. γ(Θ + eps_) + γ(X + 1) - γ(X + Θ + eps_)

    B = @. (Θ + X) * log( 1 + M / (Θ + eps_) ) +
        X * log(Θ + eps_) - log(M + eps_)
    C = A .+ B
    D = @. C - log(1f0 - Π + eps_)
    E = @. (Θ/( Θ + M + eps_))^Θ
    F = @. -log(Π + (1f0 - Π) * E + eps_)

    # if an element in matrix X is less than 0,
    # return correspond element in matrix F, else return the one in matrix D,
    # finally sum them up.
    sum(ifelse.(X .< 1f-8, F, D))
end

# demo data

r,c = 20155, 7576
X = sample(1:20, r * c) |> x -> reshape(x, r, c)
Θ = rand(Float32,r,c)
M = rand(Float32,r,c)
Π = rand(Float32,r,c)

@time myFunc(Θ,M,Π,X)
# 21.572525 seconds (36 allocations: 6.257 GiB, 4.96% gc time)

```

Any suggestions are appreciated!

---

<div class="post-metadata">

**Author:** ![JorizovdZ](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jorizovdz/32/43064_2.png) [@JorizovdZ](https://discourse.julialang.org/u/JorizovdZ)\
**Post date:** [November 17, 2022, 10:10am UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/2 "2022-11-17T10:10:08Z")

</div>

> [@tiZ](#):
>
> ` C = t1 .+ t2`

t1 and t2 are not defined in your function or demo data

> [@tiZ](#):
>
> ` X = sample(1:20, r * c) |> x -> reshape(X, r, c)`

The ‘middle’ X should also be capital.

---

<div class="post-metadata">

**Author:** ![tiZ](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiz/32/211800_2.png) [@tiZ](https://discourse.julialang.org/u/tiZ)\
**Post date:** [November 17, 2022, 10:16am UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/3 "2022-11-17T10:16:04Z")

</div>

Sorry for my typos, I have correct it.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 17, 2022, 11:01am UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/4 "2022-11-17T11:01:15Z")

</div>

Each line of your code creates a new intermediate array. It will be much faster if you merge all the operations on the elements in a single loop, like:

```julia
for i in eachindex(X)
     a = ...
     b = ...
     etc...
end

```

where `a` is _one element_ of what you are computing as `A`, for example.

(didn’t check the details)

---

<div class="post-metadata">

**Author:** ![JorizovdZ](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jorizovdz/32/43064_2.png) [@JorizovdZ](https://discourse.julialang.org/u/JorizovdZ)\
**Post date:** [November 17, 2022, 11:04am UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/5 "2022-11-17T11:04:57Z")

</div>

I am not that good at optimizing, however you can make an additional allocation which decreases the computation time. The calculation of `Θ + eps_` is executed multiple times, so can be allocated with `eps_Θ = @. Θ + eps_`:

```julia
function myFunc(Θ,M,Π,X)
    eps_ = 1f-10
    eps_Θ = @. Θ + eps_
    A = @. γ(eps_Θ) + γ(X + 1) - γ(X + eps_Θ)

    B = @. (Θ + X) * log( 1 + M / (eps_Θ) ) +
        X * log(eps_Θ) - log(M + eps_)
    C = @. A + B
    D = @. C - log(1f0 - Π + eps_)
    E = @. (Θ/(eps_Θ + M))^Θ
    F = @. -log(Π + (1f0 - Π) * E + eps_)

    # if an element in matrix X is less than 0,
    # return correspond element in matrix F, else return the one in matrix D,
    # finally sum them up.
    sum(ifelse.(X .< 1f-8, F, D))
end

```

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [November 17, 2022, 11:19am UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/6 "2022-11-17T11:19:49Z")

</div>

Whenever you introduce a function alias, do so with `const`, otherwise it will be considered a non-const global and be bad for performance. When you call it just a few times in vectorized form it doesn’t matter so much but as soon as you call it in a hot loop it will be a disaster.

```julia
const γ = loggamma # let the γ represents the loggamma() function

```

In your function you compute a lot of values pointwise and then sum them up. There is no need to have all values available at the same time, and in particular not all the intermediate values in the calculations. Moreover, if you don’t vectorize each step, you can avoid computing both (A, B, C, D) and (E, F).

This is an easy transformation which avoids most of the unnecessary allocations and computations:

```julia
myFunc(Θ, M, Π, X) = sum(myFunc2.(Θ, M, Π, X))

function myFunc2(Θ, M, Π, X)
    eps_ = 1f-10
    if X < 1f-8
        E = (Θ/( Θ + M + eps_))^Θ
        F = -log(Π + (1f0 - Π) * E + eps_)
    else
        A = γ(Θ + eps_) + Float32(γ(X + 1)) - γ(X + Θ + eps_)
        B = (Θ + X) * log( 1 + M / (Θ + eps_) ) +
            X * log(Θ + eps_) - log(M + eps_)
        C = A + B
        D = C - log(1f0 - Π + eps_)
    end
end

```

---

<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:** [November 17, 2022, 1:30pm UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/7 "2022-11-17T13:30:37Z")

</div>

Slightly faster, and with zero allocations:

```julia
myFunc(Θ, M, Π, X) = sum(t->myFunc2(t...), zip(Θ, M, Π, X))

```

The other version allocates half a Gigabyte, so getting rid of allocations entirely can be useful.

---

<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:** [November 17, 2022, 1:34pm UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/8 "2022-11-17T13:34:55Z")

</div>

> [@tiZ](#):
>
> `X = sample(1:20, r * c) |> x -> reshape(x, r, c)`

This is just the same as

```julia
X = rand(1:20, r, c)

```

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [November 17, 2022, 1:50pm UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/9 "2022-11-17T13:50:38Z")

</div>

> [@DNF](#):
>
> Slightly faster, and with zero allocations:
> 
> ```julia
> myFunc(Θ, M, Π, X) = sum(t->myFunc2(t...), zip(Θ, M, Π, X))
> 
> ```

I skipped this transformation to keep it simpler but also because I had a nagging feeling that it wouldn’t be safe to miss out on the pairwise summation used for arrays, given this large number of low precision numbers. This illustrates the dangers:

```julia
julia> f(Θ, M, Π, X) = 0.1f0
f (generic function with 1 method)

julia> sum(t -> f(t...), zip(Θ, M, Π, X))
2.097152f6

julia> sum(f.(Θ, M, Π, X))
1.526943f7

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 17, 2022, 4:22pm UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/10 "2022-11-17T16:22:59Z")

</div>

Yet, with @DNF version, we can do:

```julia
julia> myFunc(Θ, M, Π, X) = sum(myFunc2.(Θ, M, Π, X))

julia> @time myFunc(Θ, M, Π, X) 
 10.241626 seconds (3 allocations: 582.483 MiB, 0.09% gc time)
4.7243494f8

julia> using ThreadsX

julia> myFunc3(Θ, M, Π, X) = ThreadsX.sum(t -> myFunc(t...), zip(Θ, M, Π, X))
myFunc3 (generic function with 1 method)

julia> @time myFunc3(Θ, M, Π, X)
  2.195325 seconds (3.81 M allocations: 164.194 MiB, 2.17% gc time, 12.79% compilation time)
4.7243302f8

```

(starting Julia here with `julia -t auto` - 6 cores laptop)

(note that the results are clearly different, because of the order of the summation, thus the non-parallel version is “more correct” _obs: using the summation over the array_ - weather this is useful or not in specific cases is to be checked).

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [November 17, 2022, 5:32pm UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/11 "2022-11-17T17:32:25Z")

</div>

> [@lmiq](#):
>
> (note that the results are clearly different, because of the order of the summation, thus the non-parallel version is “more correct” - weather this is useful or not in specific cases is to be checked)

On the contrary, the parallel summation is usually more accurate than the sequential summation. Compare against `sum(Float64.(myFunc2.(Θ, M, Π, X)))`.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 17, 2022, 5:39pm UTC](https://discourse.julialang.org/t/any-suggestion-for-speeding-up-my-function/90393/12 "2022-11-17T17:39:01Z")

</div>

I meant _your_ non-parallel version (the one that uses the algorithm for summation over arrays). I clarified that now by putting the definition of the function I was comparing to there.