# Is this the maximum perf i can obtain?

**URL:** <https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022>\
**Category:** Performance\
**Created:** [January 4, 2022, 9:50am UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022 "2022-01-04T09:50:08Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 4, 2022, 9:50am UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022/1 "2022-01-04T09:50:08Z")

</div>

Hi,

I have a simple operation to make, which represents 80% of my runtime. I am trying to make it as fast as possible… Here is the sketch of the problem :

```julia
using Statistics, BenchmarkTools
function exp_and_mean!(slack,D,N)
    slack .= exp.(-D)
    return mean(slack)
end
function prod_and_mean!(slack,D,N)
    slack .*= D
    return mean(slack)
end
function runtime_test!(rez,slack,D)
    n = length(rez)
    N = length(D)
    rez[1] = exp_and_mean!(slack,D,N)
    for i in 1:n-1
        rez[i+1] = prod_and_mean!(slack,D,N)
    end
    return rez
end

# Simple but typical use-case : 
@benchmark runtime_test!(pars...) setup=(pars=(rez=zeros(20), slack=zeros(10000), D=exp.(randn(10000))))

```

Which outputs :

```julia
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 95.300 μs … 5.751 ms ┊ GC (min … max): 0.00% … 97.30%
 Time (median): 121.900 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 132.294 μs ± 134.750 μs ┊ GC (mean ± σ): 3.53% ± 3.56%

   ▇ █▅    
  ▁██▅▅▃▂▂▂██▇▅▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  95.3 μs Histogram: frequency by time 257 μs <

 Memory estimate: 78.20 KiB, allocs estimate: 2.

```

The first thing I have done was removing allocations and vectorising the loops :

```julia
# Second version : 
using LoopVectorization
function exp_and_mean!(slack,D,N)
    zz = zero(eltype(D))
    @turbo for i in 1:N
        slack[i] = exp(-D[i])
        zz += slack[i]
    end
    zz /= N
    return zz
end
function prod_and_mean!(slack,D,N)
    zz = zero(eltype(D))
    @turbo for i in 1:N
        slack[i] *= D[i]
        zz += slack[i]
    end
    zz /= N
    return zz
end
@benchmark runtime_test!(pars...) setup=(pars=(rez=zeros(20), slack=zeros(10000), D=exp.(randn(10000))))

```

Which outputs :

```julia
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 37.500 μs … 164.800 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 41.100 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 42.453 μs ± 6.724 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▃▇██▆▆▇▇▇▆▅▄▃▁▁ ▁ ▂
  ████████████████▇▇▆▆▅▆▅▇▇▆▇█▇▇█▇▇▇█▇▇▇▇▇▅▅▆▅▅▅▅▅▅▃▅▅▄▄▅▆▅▅▅▅ █
  37.5 μs Histogram: log(frequency) by time 75 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

This is a lot better, as the time was almost divided by 3.

I was wondering if we could do more, if I still have problems, or if I am fighting against the bare metal here. What do you think ?

---

<div class="post-metadata">

**Author:** ![Samuel3008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samuel3008/32/25021_2.png) [@Samuel3008](https://discourse.julialang.org/u/Samuel3008)\
**Post date:** [January 4, 2022, 10:27am UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022/2 "2022-01-04T10:27:29Z")

</div>

This would be well suited for a GPU. Do you happen to have access to one?

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 4, 2022, 10:30am UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022/3 "2022-01-04T10:30:28Z")

</div>

I do have an AMD one on my machine, but i am not sure i want to go in this rabbithole… This code is part of my loss function, and gets evaluated billions of times. If i understood correclty, each time the data should be passed on to the GPU, and then passed back to the CPU for the rest of the code (which might not run on gpu) to be evaluated.

---

<div class="post-metadata">

**Author:** ![IlianPihlajamaa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ilianpihlajamaa/32/28766_2.png) [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)\
**Post date:** [January 4, 2022, 11:12am UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022/4 "2022-01-04T11:12:57Z")

</div>

Disclaimer: I am on my phone so I can’t benchmark anything. To me the most straightforward optimization here is to do multithreading. Changing `@turbo` to `@tturbo` might work out of the box (you need to start Julia with multiple threads). Check the results to be safe, though.

This of course assumes you have multiple cores at your disposal.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 4, 2022, 11:55am UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022/5 "2022-01-04T11:55:13Z")

</div>

Indeed, on my 4 cores laptop with just tturbo instead of turbo, i have now :

```julia
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 16.500 μs … 5.564 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 25.900 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 28.590 μs ± 63.139 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▂▇▇▃▄▄▇███▇▅▄▃▂▁ ▁▁▁ ▃
  ████████████████████████▇█▇▇▆▇▆▇▅▆▇▆▆▅▆▆▅▄▄▅▆▄▄▆▄▂▄▄▅▄▄▃▃▄▄ █
  16.5 μs Histogram: log(frequency) by time 84.8 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

Which is ore than twice faster. Very neat. Is this the best thing we can do ?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [January 4, 2022, 2:19pm UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022/7 "2022-01-04T14:19:01Z")

</div>

Does @fastmath help for the `exp()` computations?

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 4, 2022, 3:07pm UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022/8 "2022-01-04T15:07:36Z")

</div>

In facts, @turbo and @tturbo already implies @fastmath. I checked and the results seems correct.

---

<div class="post-metadata">

**Author:** ![Bardo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bardo/32/21601_2.png) [@Bardo](https://discourse.julialang.org/u/Bardo)\
**Post date:** [January 4, 2022, 6:04pm UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022/9 "2022-01-04T18:04:17Z")

</div>

Just a thought: Is it just for the example, or are you actually calculating the first moment of distributions?

Then it might be possible to think of improving the sampling method, resulting in fewer points, using surrogates, or something more analytic like continuous or piece-wise linear integrals, only sampling the tails, …

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 4, 2022, 6:27pm UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022/10 "2022-01-04T18:27:18Z")

</div>

I am computing these moments for real, yes, but the data is here a fake data.

If you known the density, e.g. for a lognormal, I already have coded up a better integration procedure there : [https://github.com/lrnv/ThorinDistributions.jl/blob/main/src/MFKProjection.jl](https://github.com/lrnv/ThorinDistributions.jl/blob/main/src/MFKProjection.jl)

But here I only have the data, of unknown distribution 🙂

---

<div class="post-metadata">

**Author:** ![djholiver](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djholiver/32/50470_2.png) [@djholiver](https://discourse.julialang.org/u/djholiver)\
**Post date:** [January 4, 2022, 8:24pm UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022/11 "2022-01-04T20:24:16Z")

</div>

Hi,

Why do you need to allocate to slack (or even provide it as a variable) in the two methods?

```julia
 using LoopVectorization
function exp_and_mean!(D,N)
    zz = zero(eltype(D))
    @tturbo for i in 1:N        
        zz += exp(-D[i])  
    end
    zz /= N
    return zz
end
function prod_and_mean!(D,N)
    zz = one(eltype(D))
    @tturbo for i in 1:N
        zz *= D[i]
    end
    zz /= N
    return zz
end

```

is there also potentially an error in “prod\_and\_mean” ? do you mean to sum zz (I change to multiple above).

From your input setup vars, you are passing in slack as a vector of zeros so doesnt this yield zero for all cases in “prod\_and\_mean” ? Obviously they may be intended.

Do you have any opportunity to modify all “Ds” before being passing into the functions at all or are they always generated before being called? You might benefit from creating the entire “D” matrix space in exp and not exp forms before executing your functions (you’d also be able to combine them if this is applicable).

Regards,

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 5, 2022, 9:55am UTC](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022/12 "2022-01-05T09:55:25Z")

</div>

> [@djholiver](#):
>
> Why do you need to allocate to slack (or even provide it as a variable) in the two methods?

because I actually need the result for later computations. Look again at my `runtime_test!` function:

```julia
function runtime_test!(rez,slack,D)
    n = length(rez)
    N = length(D)
    rez[1] = exp_and_mean!(slack,D,N)
    for i in 1:n-1
        rez[i+1] = prod_and_mean!(slack,D,N)
    end
    return rez
end

```

You see that the slack is passed along from one function to the other. Actually, if you denote by X the rndom variable that corresponds to the data `D`, my function computes \mathbb E\left(X^k e^{-X}\right), for k \in \{0,...,n-1\}.

> [@djholiver](#):
>
> is there also potentially an error in “prod\_and\_mean” ? do you mean to sum zz (I change to multiple above).
> 
> From your input setup vars, you are passing in slack as a vector of zeros so doesnt this yield zero for all cases in “prod\_and\_mean” ? Obviously they may be intended.

No there is no error in “prod\_and\_mean”, I indeed mean to sum `zz`. Yes this is intended, since `slack` is not expected to be zeros when calling the `prod_and_mean` function, look again at the `runtime_test!` function 🙂

> [@djholiver](#):
>
> Do you have any opportunity to modify all “Ds” before being passing into the functions at all or are they always generated before being called? You might benefit from creating the entire “D” matrix space in exp and not exp forms before executing your functions (you’d also be able to combine them if this is applicable).

I do not understand what you mean. This computation is done only once for a given dataset `D`, I am afraid. The dataset `D` changes from one execution to the next, albeit not completely randomly: I compute it from another dataset `data` as follows:

```julia
data=.... #(sized (10,10000)
for i in 1:n_iterations 
    e = rand(10) # A new one is picked at each iteration 
    D = data'e # now a vector of size 10000
    runtime_test!(rez,slack,D)
    do_something_with!(rez)
end

```

`data` is fixed from one iteration to the other, but `e` is not (and therefore neither `D`, `slack` and `rez`). `rez` is needed for the rest of the computations, but `D` and `slack` are not.