# Speed up loop for logpdf

**URL:** <https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620>\
**Category:** General Usage\
**Tags:** performance\
**Created:** [October 23, 2024, 3:50am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620 "2024-10-23T03:50:33Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ian\_L](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ian_l/32/49509_2.png) [@Ian\_L](https://discourse.julialang.org/u/Ian_L)\
**Post date:** [October 23, 2024, 3:50am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/1 "2024-10-23T03:50:33Z")

</div>

I have the following code which tries to evaluate the logpdf of multiple categorical variables. It is ok, but apparently it is one of the bottlenecks in my code. Really, all the code is doing is for each datapoint `xs[i]`, select the `ith` vector from a `Vector{Vector{Float64}}`, index into it using `xs[i]`, and add the value to a sum.

```julia
function logpdf(cat::Categorical{T}, xs::Vector{Int}) where T
  length(xs) != length(cat.logp) && throw(DomainError("lengths do not match"))
  p = zero(T)
  @inbounds for i in axes(xs,1)
    logp = cat.logp[i] # cat.logp isa Vector{Vector{Float64}}
    x = xs[i]
    p += logp[x]
  end
  return p
end

```

Right now on my M1 I’m getting roughly

```julia
  5.708 ns (0 allocations: 0 bytes)

```

I tried LoopVectorization by tacking a `@turbo` but am running into this error:

```julia
UndefVarError: `logp` not defined in `GenSPN`

```

Is it possible to make the code faster? Note I cannot use a `Matrix`.

---

<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:** [October 23, 2024, 5:15am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/2 "2024-10-23T05:15:20Z")

</div>

> [@Ian\_L](#):
>
> ```julia
> x = xs[i]
> p += logp[x]
> 
> ```

I’m confused. If `xs` is a `Vector{Float64}`, then `x` is a `Float64`, and then it can’t be used as an index. So this code looks like it shouldn’t work.

Also, this comment cannot be correct:

> [@Ian\_L](#):
>
> `cat.logp[i] isa Vector{Vector{Float64}}`

---

<div class="post-metadata">

**Author:** ![Ian\_L](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ian_l/32/49509_2.png) [@Ian\_L](https://discourse.julialang.org/u/Ian_L)\
**Post date:** [October 23, 2024, 5:27am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/3 "2024-10-23T05:27:48Z")

</div>

You are right - and fixed now. I should have double checked before making impromptu edits.

---

<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:** [October 23, 2024, 5:32am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/4 "2024-10-23T05:32:21Z")

</div>

> [@Ian\_L](#):
>
> Is it possible to make the code faster?

I think adding `@simd` here might allow the compiler to perform more optimizations, at the cost of a slightly different result since floating point addition is not associative ([Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-annotations)).

---

<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:** [October 23, 2024, 5:33am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/5 "2024-10-23T05:33:27Z")

</div>

One problem I see is that only the second indexing expression is guaranteed to be inbounds:

> [@Ian\_L](#):
>
> ```julia
> logp = cat.logp[i] # cat.logp isa Vector{Vector{Float64}}
> x = xs[i]
> p += logp[x]
> 
> ```

I would remove the `@inbounds`. That probably won’t help your performance though.

Loopvectorization.jl probably doesn’t like the jumping around in memory, vectorization is probably hard to achieve here.

---

<div class="post-metadata">

**Author:** ![Ian\_L](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ian_l/32/49509_2.png) [@Ian\_L](https://discourse.julialang.org/u/Ian_L)\
**Post date:** [October 23, 2024, 5:34am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/6 "2024-10-23T05:34:53Z")

</div>

> [@gdalle](#):
>
> I think adding `@simd` here might allow the compiler to perform more optimizations, at the cost of a slightly different result since floating point addition is not associative ([Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-annotations)).

I did try adding `@simd` in hopes of a speed up, but it seems I get identical timings.

---

<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:** [October 23, 2024, 5:35am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/7 "2024-10-23T05:35:58Z")

</div>

Then I’d say your pure Julia code is already close to optimal (although I’d also remove `@inbounds` to favor safety over speed), especially given your timings in the ns range.

---

<div class="post-metadata">

**Author:** ![Ian\_L](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ian_l/32/49509_2.png) [@Ian\_L](https://discourse.julialang.org/u/Ian_L)\
**Post date:** [October 23, 2024, 5:36am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/8 "2024-10-23T05:36:46Z")

</div>

> [@DNF](#):
>
> I would remove the `@inbounds`. That probably won’t help your performance though.

I can’t remember exactly, but I think I `@inbounds` got me from 7ns to 5ns. 😅

---

<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:** [October 23, 2024, 5:40am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/9 "2024-10-23T05:40:03Z")

</div>

Well if you replace `axes(xs, 1)` with `eachindex(xs, cat.logp)`, then at least two indices are safe (use `axes` for indexing along dimensions of higher-dimensional arrays, not vectors). But this index is not generally safe: `logp[x] `. Is it guaranteed to be safe somehow?

The fact that `logp` is a new vector for each iteration makes it hard (or impossible) to vectorize, you want to extract consecutive, or at least closely spaced samles, but different vectors are probably not close in memory. If each `logp` has identical size, maybe something can be done.

---

<div class="post-metadata">

**Author:** ![Ian\_L](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ian_l/32/49509_2.png) [@Ian\_L](https://discourse.julialang.org/u/Ian_L)\
**Post date:** [October 23, 2024, 5:42am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/10 "2024-10-23T05:42:24Z")

</div>

Other than assertion checks outside of the function, no. I think I’ve been convinced to remove `@inbounds` since it doesn’t seem worth it.

---

<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:** [October 23, 2024, 6:34am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/11 "2024-10-23T06:34:42Z")

</div>

If you use `eachindex` as suggested, you can try `@inbounds` on the first two indexing expressions, though it may make no difference:

```julia
for i in eachindex(xs, cat.logp)
    @inbounds logp = cat.logp[i]
    @inbounds x = xs[i]
    p += logp[x]
end

```

---

<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:** [October 23, 2024, 6:54am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/12 "2024-10-23T06:54:01Z")

</div>

What is the typical length of `cat.logp` and of each of its elements?

---

<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:** [October 23, 2024, 6:54am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/13 "2024-10-23T06:54:10Z")

</div>

Why can’t you use a `Matrix`?  
From what I understand, your `Categorical` is not from `Distributions` but some struct holding a `Vector{Vector{Float64}}` giving the log-probabilities to be applied for each observation. You could just bring these to a common length by including `-Inf` when the index is out of the support of the distribution:

```julia
nmax = maximum(length.(cat.logp))
mat = [get(cat.logp[i], j, -Inf) for i in eachindex(cat.logp), j in 1:nmax]

```

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [October 23, 2024, 7:00am UTC](https://discourse.julialang.org/t/speed-up-loop-for-logpdf/121620/14 "2024-10-23T07:00:20Z")

</div>

Taking a step back, the function is pretty short and simple, and takes 5ns with zero allocations currently, do you have reason to believe that that is unreasonably slow for what it does? Eg do you have an implementation in another language that is materially faster?