# Could this kernel go faster?

**URL:** <https://discourse.julialang.org/t/could-this-kernel-go-faster/73986>\
**Category:** Performance\
**Created:** [January 3, 2022, 6:30pm UTC](https://discourse.julialang.org/t/could-this-kernel-go-faster/73986 "2022-01-03T18:30:45Z")\
**Posts on this page:** 4\
**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 3, 2022, 6:30pm UTC](https://discourse.julialang.org/t/could-this-kernel-go-faster/73986/1 "2022-01-03T18:30:45Z")

</div>

Hi,

I have the follwing problem :

```julia
using Statistics, BenchmarkTools, ProfileView

function compute_exp_moments!(mu,D,e,slack,data)
    data .= D'e # most expensive line according to profileview. 
    slack .= exp.(.-data)
    mu[1] = mean(slack)
    for i in 2:length(mu)
        slack .*= data
        mu[i] = mean(slack)/mu[1]
    end
end

function compute_many(N,n,d,How_much)
    D = reshape(exp.(randn(d*N)),(d,N))
    slack = zeros(N)
    data = zeros(N)
    mu = zeros(n)
    for i in 1:How_much
        e = rand(d)
        compute_exp_moments!(mu,D,e,slack,data)
    end
    return mu
end

@btime compute_many(10000,10,20,100) # 20ms for me 
ProfileView.@profview compute_many(10000,10,20,1000)

```

The computation seems rather straightforward, but i have trouble making it faster. Much of the time is spend on the first matrix-vector product, as expected.

I tried Loopvectorisation together with StaticArrays and HybridArrays, improved by 20% :

```julia
function compute_exp_moments2!(mu,D,e,slack,data)
    @turbo data .= D'e # Still takes much of the time. 
    @turbo slack .= exp.(.-data)
    mu[1] = mean(slack)
    for i in 2:length(mu) # I did not manage to use @turbo on this one. 
        slack .*= data
        mu[i] = mean(slack)/mu[1]
    end
end
function compute_many2(N,n,d,How_much)
    D = HybridArray{Tuple{d,StaticArrays.Dynamic()}}(reshape(exp.(randn(d*N)),(d,N)))
    slack = zeros(N)
    data = zeros(N)
    mu = zeros(n)
    for i in 1:How_much
        e = SVector{d}(rand(d))
        compute_exp_moments2!(mu,D,e,slack,data)
    end
    return mu
end

@btime compute_many2(10000,10,20,100) #16 ms, better. 
ProfileView.@profview compute_many2(10000,10,20,1000)

```

But I am still looking for performance. This stripped-down example accounts for 75% of my actual runtime.

Is there anything more that I can do ?

---

<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:** [January 3, 2022, 8:17pm UTC](https://discourse.julialang.org/t/could-this-kernel-go-faster/73986/2 "2022-01-03T20:17:25Z")

</div>

You could should preallocate `e` and use `rand!(e)` in the loop. Also, you can use `mul!(data, D', e)` to avoid allocations for `D'e`.

```julia
using Statistics, BenchmarkTools, Random, LinearAlgebra

function compute_exp_moments!(mu, D, e, slack, data)
    mul!(data, D', e)
    slack .= exp.(.-data)
    mu[1] = mean(slack)
    @inbounds for i in 2:length(mu)
        slack .*= data
        mu[i] = mean(slack) / mu[1]
    end
end

function compute_many(N, n, d, How_much)
    D = reshape(exp.(randn(d * N)), (d, N))
    slack = zeros(N)
    data = zeros(N)
    mu = zeros(n)
    e = zeros(d)
    for i in 1:How_much
        rand!(e)
        compute_exp_moments!(mu, D, e, slack, data)
    end
    return mu
end

```

Timings:

```julia
# original: 22.240 ms (311 allocations: 10.86 MiB)
# prealloc e + rand!(e): 22.403 ms (212 allocations: 10.84 MiB)
# mul!(data, D', e): 16.343 ms (12 allocations: 3.20 MiB)

```

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [January 4, 2022, 12:22am UTC](https://discourse.julialang.org/t/could-this-kernel-go-faster/73986/3 "2022-01-04T00:22:40Z")

</div>

I’m missing the usual reference to `using MKL` yet;)

---

<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/could-this-kernel-go-faster/73986/5 "2022-01-04T09:50:46Z")

</div>

Thanks to all of you, this code is getting better and better. However, this post is messy and I have managed to pin out more precisely my bottleneck. Therefore I reposted it there : [Is this the maximum perf i can obtain?](https://discourse.julialang.org/t/is-this-the-maximum-perf-i-can-obtain/74022)