# How to speed up the calculation of Logit likelihood?

**URL:** <https://discourse.julialang.org/t/how-to-speed-up-the-calculation-of-logit-likelihood/104153>\
**Category:** Performance\
**Tags:** question\
**Created:** [September 22, 2023, 1:21pm UTC](https://discourse.julialang.org/t/how-to-speed-up-the-calculation-of-logit-likelihood/104153 "2023-09-22T13:21:24Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Strange\_Xue](https://avatars.discourse-cdn.com/v4/letter/s/e9a140/32.png) [@Strange\_Xue](https://discourse.julialang.org/u/Strange_Xue)\
**Post date:** [September 22, 2023, 1:21pm UTC](https://discourse.julialang.org/t/how-to-speed-up-the-calculation-of-logit-likelihood/104153/1 "2023-09-22T13:21:24Z")

</div>

Hi, I am currently using Julia to solve some econometric structural models which usually involve the calculation and optimization of Log-likelihood function. While the optimization can be done by Optim.jl, the calculation of loglikelihood needs hand-coding.

The question I face is that, how to calculate the log-likelihood function as fast as possible especially facing the unbalanced panel data as the input. The code below giving a logit implementation has the inputs `X` and `y` in which each individual has different number of observations, so I use `Vector{Matrix{Float64}}` to save `X` and use `Vector{Vector{Int64}}` to save `y`.

I wonder if this is an efficient way to calculate loglikehood function `nll_threading`? Is there any other techniques to handle such problems ? (PS: I already used many methods such as multi-threading, views to speed up; however it’s still slow especially when the loglikelihood function for each individual has a more complicated form, such as hidden Markov model, bayesian learning model)

```julia
using LogExpFunctions, Optim
using Random, Distributions
using LinearAlgebra:dot
using BenchmarkTools
const seed = Random.seed!(2022 - 1 - 3)
# data generation
begin
    Inds = 1000 # number of individuals
    nX = 4
    Tlen = rand([15, 10, 20], Inds) # number of observations of each individual
    X = Vector{Matrix{Float64}}()
    for i in eachindex(Tlen)
        push!(X, randn(seed, (Tlen[i], nX)))
    end
    beta = [1.0, 2, -3, 4]
    y = Vector{Vector{Int64}}()
    for i in 1:Inds
        Xi = X[i]
        yi = Vector{Int64}()
        for t in axes(Xi, 1)
            U = @views dot(Xi[t, :], beta)
            push!(yi, rand(BernoulliLogit(-U)))
        end
        push!(y, yi)
    end
end
# negative log-likelihood function (to optimize)
function nll_threading(X, beta::Vector{T}, y) where {T <: Real}
    ind_ll = zeros(T, Threads.nthreads())
    Threads.@threads for i in eachindex(y)
        Xi = X[i]
        yi = y[i] 
        for t in axes(Xi, 1)
            U = @views dot(Xi[t, :], beta)
            ind_ll[Threads.threadid()] += logpdf(BernoulliLogit(-U), yi[t])
        end
    end
    return - sum(ind_ll)
end
# benchmark
@benchmark nll_threading($X, $beta, $y)
# about 380ms for 6 threads with intel i5-10400f chip

# optimization
fit_optim = optimize(
    pars -> nll_threading(X, pars, y),
    zeros(4),
    Newton();
    autodiff = :forward
)
fit_optim.minimizer

```

Benchmark for current version

```julia
julia> @btime nll_threading($X, $beta, $y)
  551.400 μs (111 allocations: 12.39 KiB)

```

---

<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:** [September 22, 2023, 1:29pm UTC](https://discourse.julialang.org/t/how-to-speed-up-the-calculation-of-logit-likelihood/104153/2 "2023-09-22T13:29:25Z")

</div>

I think it might help to make `X` be a `Matrix{SVector{4}}` and `beta` be an `SVector{4}`. Doing so will make it so the dot product calculation knows ahead of time that it will be getting 4 long vectors. Another likely 2x speedup is switching to `Float32` which should have enough accuracy but will take 2x less memory.

---

<div class="post-metadata">

**Author:** ![Strange\_Xue](https://avatars.discourse-cdn.com/v4/letter/s/e9a140/32.png) [@Strange\_Xue](https://discourse.julialang.org/u/Strange_Xue)\
**Post date:** [September 23, 2023, 5:47am UTC](https://discourse.julialang.org/t/how-to-speed-up-the-calculation-of-logit-likelihood/104153/3 "2023-09-23T05:47:25Z")

</div>

Another Question 🙂: I also face another weird problem. If I run it in my local machine with 6 threads (i5 10400f, 6 physical processors), the CPU usage is high, over 80%; However, if run it in the remote windows server with 18 threads (i9 10980XE, 18 physical processors), the CPU usage is relatively low, at 40% approximately. It seems Julia didn’t utilize all threads as I asked for. I already set number of threads to 18 while starting the REPL.

---

<div class="post-metadata">

**Author:** ![Strange\_Xue](https://avatars.discourse-cdn.com/v4/letter/s/e9a140/32.png) [@Strange\_Xue](https://discourse.julialang.org/u/Strange_Xue)\
**Post date:** [September 23, 2023, 5:50am UTC](https://discourse.julialang.org/t/how-to-speed-up-the-calculation-of-logit-likelihood/104153/4 "2023-09-23T05:50:11Z")

</div>

Thanks. It seems helpful, I will post the benchmark results after implementing your ideas.

---

<div class="post-metadata">

**Author:** ![jonathanBieler](https://avatars.discourse-cdn.com/v4/letter/j/82dd89/32.png) [@jonathanBieler](https://discourse.julialang.org/u/jonathanBieler)\
**Post date:** [September 23, 2023, 9:24am UTC](https://discourse.julialang.org/t/how-to-speed-up-the-calculation-of-logit-likelihood/104153/5 "2023-09-23T09:24:13Z")

</div>

For some problems you can avoid computing the likelihood by grouping observations, e.g. :

```julia
D = Poisson(100)
x = rand(D,1000)
xc = countmap(x)

sum(logpdf(D,xi)*n for (xi,n) in xc) #56 iterations instead of 1000

```

But I guess this won’t help in your case, unless maybe if you round off U to some decimals.

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [September 23, 2023, 9:44am UTC](https://discourse.julialang.org/t/how-to-speed-up-the-calculation-of-logit-likelihood/104153/6 "2023-09-23T09:44:53Z")

</div>

Perhaps putting Threads.@threads before the inner loop,

```julia
for t in axes(Xi, 1)

```

might be better.

---

<div class="post-metadata">

**Author:** ![Strange\_Xue](https://avatars.discourse-cdn.com/v4/letter/s/e9a140/32.png) [@Strange\_Xue](https://discourse.julialang.org/u/Strange_Xue)\
**Post date:** [September 23, 2023, 10:55am UTC](https://discourse.julialang.org/t/how-to-speed-up-the-calculation-of-logit-likelihood/104153/7 "2023-09-23T10:55:46Z")

</div>

Might be worse.

```julia
# two implementations: one use vector to save each individual's likelihood, 
# the other use lock to avoid data race in multi-threading. 
# Both performances are worse than original version.
function nll_threading_inner1(X, beta::Vector{T}, y) where {T <: Real}
    ll = zero(T)
    @inbounds for i in eachindex(y)
        Xi = X[i]
        yi = y[i] 
        lli = zeros(T, length(yi))
        Threads.@threads for t in axes(Xi, 1)
            U = @views dot(Xi[t, :], beta)
            lli[t] = logpdf(BernoulliLogit(-U), yi[t])
        end
        ll += sum(lli)
    end
    return - ll
end

function nll_threading_inner2(X, beta::Vector{T}, y) where {T <: Real}
    ll = zero(T)
    l = ReentrantLock()
    @inbounds for i in eachindex(y)
        Xi = X[i]
        yi = y[i] 
        Threads.@threads for t in axes(Xi, 1)
            U = @views dot(Xi[t, :], beta)
            lock(l) do
                ll += logpdf(BernoulliLogit(-U), yi[t])
            end
        end
    end
    return - ll
end

```

Benchmark results

```julia
julia> @btime nll_threading_inner1($X, $beta, $y)
  20.123 ms (107235 allocations: 11.96 MiB)
julia> @btime nll_threading_inner2($X, $beta, $y)
  22.684 ms (137368 allocations: 12.28 MiB)

```

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [September 23, 2023, 11:16am UTC](https://discourse.julialang.org/t/how-to-speed-up-the-calculation-of-logit-likelihood/104153/8 "2023-09-23T11:16:08Z")

</div>

When dealing with matrices, you always want to go over the faster index instead of slower ones to increase memory locality. In Julia, the earlier indices are faster, and therefore you would want instead of:

> [@Strange\_Xue](#):
>
> `@views dot(Xi[t, :], beta)`

something like:

```julia
@view dot(Xi[:, t], beta)

```

This would require transposing the generation process as well:

> [@Strange\_Xue](#):
>
> `randn(seed, (Tlen[i], nX))`

to

```julia
randn(seed, (nX, Tlen[i]))

```

This change gives me a substantial improvement in my local machine.

---

<div class="post-metadata">

**Author:** ![Strange\_Xue](https://avatars.discourse-cdn.com/v4/letter/s/e9a140/32.png) [@Strange\_Xue](https://discourse.julialang.org/u/Strange_Xue)\
**Post date:** [September 23, 2023, 11:41am UTC](https://discourse.julialang.org/t/how-to-speed-up-the-calculation-of-logit-likelihood/104153/9 "2023-09-23T11:41:57Z")

</div>

This should be the easiest way to improve. 5x speed up after getting index column by column. Thanks.

```julia
julia> @btime nll_threading($X, $beta, $y)
  118.400 μs (109 allocations: 12.33 KiB)

```