# Any way to efficiently run Poisson regression thousands of times?

**URL:** <https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727>\
**Category:** Finance and Economics\
**Tags:** regression, glm\
**Created:** [November 26, 2023, 1:46am UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727 "2023-11-26T01:46:40Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![xxxxx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xxxxx/32/204406_2.png) [@xxxxx](https://discourse.julialang.org/u/xxxxx)\
**Post date:** [November 26, 2023, 1:46am UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/1 "2023-11-26T01:46:40Z")

</div>

I need to run Poisson regression around 3000 times with the same set of Y and X but different weights. I’m currently using GLM.jl with a for-loop. The whole thing takes about 58 seconds. I was wondering if there are more efficient ways to run this. Thanks!

An MWE:

```julia
using GLM
using Distributions

beta_t = [0.2, 0.1, 0.5, 1.1]
x = rand(100000, 4)
lambdas = exp.(x*beta_t)
y = rand.(Poisson.(lambdas))

w = rand(size(y,1), 3000)
results = zeros(3000, 4)
for i in axes(w, 2)
        results[i,:] = collect(coef(glm(x, y, Poisson(), wts=w[:,i])))
end

```

Thanks again!

---

<div class="post-metadata">

**Author:** ![RomeoV](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romeov/32/37687_2.png) [@RomeoV](https://discourse.julialang.org/u/RomeoV)\
**Post date:** [November 26, 2023, 1:51am UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/2 "2023-11-26T01:51:24Z")

</div>

The most obvious thing would be to use threading, i.e. add a `Threads.@threads` before the for-loop and launch with `julia --threads=auto`. If you’re fine with it, you can also try to use `Float32` for all the data, or even try putting a `@fastmath` in front of the for loop – however, double check your correctness in that case.

---

<div class="post-metadata">

**Author:** ![xxxxx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xxxxx/32/204406_2.png) [@xxxxx](https://discourse.julialang.org/u/xxxxx)\
**Post date:** [November 26, 2023, 2:06am UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/3 "2023-11-26T02:06:38Z")

</div>

Hi @RomeoV, thanks for your quick reply! I should have mentioned that I tried multi-threading, but the speed gain was not significant enough. Ideally, I’d like it to be under ~2 secs if that’s possible, as this step is used as an input for the later steps.

---

<div class="post-metadata">

**Author:** ![RomeoV](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romeov/32/37687_2.png) [@RomeoV](https://discourse.julialang.org/u/RomeoV)\
**Post date:** [November 26, 2023, 2:12am UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/4 "2023-11-26T02:12:46Z")

</div>

For me it threading speeds up about 4x:

```julia
julia> @time Threads.@threads for i in axes(w, 2)
               results[i,:] = collect(coef(glm(x, y, Poisson(), wts=w[:,i])))
       end
 32.137273 seconds (199.79 k allocations: 35.775 GiB, 1.09% gc time, 0.75% compilation time)

julia> @time for i in axes(w, 2)
               results[i,:] = collect(coef(glm(x, y, Poisson(), wts=w[:,i])))
       end
111.505098 seconds (179.93 k allocations: 35.773 GiB, 0.20% gc time)

```

Are you running with `julia --threads=auto` ? You can check with `Threads.nthreads()` how many threads you have available (I have 8 cores=16 threads).

---

<div class="post-metadata">

**Author:** ![RomeoV](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romeov/32/37687_2.png) [@RomeoV](https://discourse.julialang.org/u/RomeoV)\
**Post date:** [November 26, 2023, 2:16am UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/5 "2023-11-26T02:16:48Z")

</div>

Another thing you might try is to use `StaticArrays.jl` (since you seem to know your sizes), but it depends whether `GLM.jl` uses the static sizes properlly or not, which I’m not sure about.

---

<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:** [November 26, 2023, 2:45am UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/6 "2023-11-26T02:45:33Z")

</div>

Doesn’t save so much time, but having an initial guess helps:

```julia
function testit()
    beta_t = [0.2, 0.1, 0.5, 1.1]
    x = rand(100000, 4)
    lambdas = exp.(x*beta_t)
    y = rand.(Poisson.(lambdas))
    guess = collect(coef(glm(x, y, Poisson())))
    w = rand(size(y,1), 3000)
    results = zeros(3000, 4)
    for i in axes(w, 2)
        results[i,:] = collect(coef(glm(x, y, Poisson(), wts=w[:,i], start=guess)))
    end
end

```

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [November 27, 2023, 7:53am UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/7 "2023-11-27T07:53:14Z")

</div>

@xxxxx You do not say which processor you are using. It would be nice to know:

processor type  
how many cores  
hyperthreading - on or off?  
If Intel is Turboboost enabled?

If you are in a financial institution you should probably have access to a multicore system with a high clock frequency.

---

<div class="post-metadata">

**Author:** ![barucden](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barucden/32/26154_2.png) [@barucden](https://discourse.julialang.org/u/barucden)\
**Post date:** [November 27, 2023, 10:03am UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/8 "2023-11-27T10:03:50Z")

</div>

I am not a user of GLM.jl but there is a keyword argument `start` in the `glm` function. That argument can be used to set the initial values of the regression coefficients. You could try setting them based on the previously-computed coefficients.

```julia
# currently
results[i,:] = collect(coef(glm(x, y, Poisson(), wts=w[:,i])))
# would become
results[i,:] = collect(coef(glm(x, y, Poisson(), wts=w[:,i], 
                            start=results[i - 1, :])))
# or maybe?
results[i,:] = collect(coef(glm(x, y, Poisson(), wts=w[:,i], 
                            start=mean(results[1:i - 1, :], dims=1))))

```

I am not sure if these initial coefficients are guaranteed to be a closer to the optimum after changing the weights, but I would think it often should be the case. It is worth trying.

Edit: It might also help to add `@views` in the loop:

```julia
results[i,:] = @views collect(coef(glm(x, y, Poisson(), 
                            wts=w[:,i], start=results[i - 1, :])))

```

---

<div class="post-metadata">

**Author:** ![xxxxx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xxxxx/32/204406_2.png) [@xxxxx](https://discourse.julialang.org/u/xxxxx)\
**Post date:** [December 2, 2023, 5:07pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/9 "2023-12-02T17:07:55Z")

</div>

Hi @RomeoV , thanks for your help! I have implemented those, and it worked great.

---

<div class="post-metadata">

**Author:** ![xxxxx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xxxxx/32/204406_2.png) [@xxxxx](https://discourse.julialang.org/u/xxxxx)\
**Post date:** [December 2, 2023, 5:08pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/10 "2023-12-02T17:08:23Z")

</div>

Thanks @Dan , the initial guess also helps!

---

<div class="post-metadata">

**Author:** ![xxxxx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xxxxx/32/204406_2.png) [@xxxxx](https://discourse.julialang.org/u/xxxxx)\
**Post date:** [December 2, 2023, 5:09pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/11 "2023-12-02T17:09:24Z")

</div>

Thanks @barucden . Recursively using the results from the last step was also helpful!

---

<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:** [December 2, 2023, 5:23pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/12 "2023-12-02T17:23:35Z")

</div>

BTW the code in the post works nicely and is a good MWE, but the random nature of the coefficients and weights might overlook more optimization opportunities.

If the actual problem has structure, if you describe it, maybe good suggestions will come up.

---

<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:** [December 2, 2023, 5:57pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/13 "2023-12-02T17:57:50Z")

</div>

It’s not clear to me what the end goal is. If the problem has already been solved on 58 seconds, what’s left to be done? The reason that I ask is that the effort one would but into speeding this up depends on what is the actual final goal.

Going down memory lane, I once wrote a paper called “I Ran Four Million Probits Last Night: HPC Clustering with ParallelKnoppix”. Your problem reminded me of it.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [December 2, 2023, 7:43pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/14 "2023-12-02T19:43:53Z")

</div>

> [@mcreel](#):
>
> what the end goal is

Agree, the best way to make something faster is to make the thing more efficient, perhaps running 3000 repetitions could be modified to something else that accomplished the overlying goal better.

---

<div class="post-metadata">

**Author:** ![xxxxx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xxxxx/32/204406_2.png) [@xxxxx](https://discourse.julialang.org/u/xxxxx)\
**Post date:** [December 2, 2023, 8:29pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/15 "2023-12-02T20:29:34Z")

</div>

Thanks @Dan @mcreel @dlakelan ! The end goal here is to use the estimates from the Poisson regression with different weights to make predictions on the mean arrival rates, i.e., X\hat{\beta}. @mcreel I will also check out your paper, and thanks for the reference!

---

<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:** [December 2, 2023, 9:18pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/16 "2023-12-02T21:18:58Z")

</div>

It would be useful to also take this post as an opportunity to link to the related new question: [Nesting the GLM.jl in the objective function when formulating the ReverseDiff gradient generates StackOverflowError](https://discourse.julialang.org/t/nesting-the-glm-jl-in-the-objective-function-when-formulating-the-reversediff-gradient-generates-stackoverflowerror/107050)

That question is a better description for the end-goal as @mcreel suggested. Especially, it also has a clue (yes, it is a clue, because readers need to undo the transformations done to protect the ‘proprietary’ application), to the weights generation.

The weights in that question are all scalar multiples of each other, which essentially means `predict_data` in that question is independent of the `para[3]` scaling factor and only depends on `data[:,:w]`. It also means inference on `para[3]` is impossible (it’s statistically unidentifiable).

---

<div class="post-metadata">

**Author:** ![xxxxx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xxxxx/32/204406_2.png) [@xxxxx](https://discourse.julialang.org/u/xxxxx)\
**Post date:** [December 2, 2023, 9:24pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/17 "2023-12-02T21:24:39Z")

</div>

Hi @Dan , thanks, but these two are uncorrelated questions. Altough both questions concern about Poisson regression and GLM.jl, they have different end applications.

In that particular question’s MWE, I understand that para[3] is unidentifiable, but even if it’s unidentifiable, I don’t think that is the cause for the StackOverFlow error that I encountered. In that question, I was particularly curious if it would be possible to combine the GLM with the ReverseDiff.

In this question, the end usage would be what I described, purely as making the predictions based on the estimated coefficients from different weights. Thanks again.

---

<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:** [December 2, 2023, 9:49pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/18 "2023-12-02T21:49:40Z")

</div>

No need to read that paper! It just showed how MPI could be used to speed up Monte Carlo, an idea that was not well-known to economists at the time. Actually, MPI would be one way to speed up your problem, using MPI.jl. The documentation for that would be a better read than my paper. But, perhaps it would be easier to use Julia’s built-in Distributed package.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [December 2, 2023, 10:48pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/19 "2023-12-02T22:48:34Z")

</div>

> [@xxxxx](#):
>
> The end goal here is to use the estimates from the Poisson regression with different weights to make predictions on the mean arrival rates

I guess what’s not clear is why 3000 different sets of weights. Obviously if you only need to do it for 1 set of weights then it can be done FAR faster than 58 seconds or whatever. And the amount of time grows linearly with the number of weights at least to first approximation… so, the 3000 samples must play some role, and the question of interest for speeding this up is what is that role and can that role be played by some other methodology, for example is there 30 or 300 different chosen weights which can give effectively the same answer as the full 3000 or something like that.

So, if you explain the role that resampling with different weights plays people may be able to point you at algorithms that use far fewer fits to get you equally good or even better answers (or possibly there is no such answer) but without more knowledge there we can only help you with “how to make the software do 3000 of them in less time”

---

<div class="post-metadata">

**Author:** ![xxxxx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xxxxx/32/204406_2.png) [@xxxxx](https://discourse.julialang.org/u/xxxxx)\
**Post date:** [December 3, 2023, 7:13pm UTC](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727/20 "2023-12-03T19:13:36Z")

</div>

I will try those as well. Thanks!

[Next page](https://discourse.julialang.org/t/any-way-to-efficiently-run-poisson-regression-thousands-of-times/106727.md?page=2)