# Nesting the GLM.jl in the objective function when formulating the ReverseDiff gradient generates StackOverflowError

**URL:** <https://discourse.julialang.org/t/nesting-the-glm-jl-in-the-objective-function-when-formulating-the-reversediff-gradient-generates-stackoverflowerror/107050>\
**Category:** Statistics\
**Tags:** glm, reversediff\
**Created:** [December 2, 2023, 5:51pm UTC](https://discourse.julialang.org/t/nesting-the-glm-jl-in-the-objective-function-when-formulating-the-reversediff-gradient-generates-stackoverflowerror/107050 "2023-12-02T17:51:10Z")\
**Posts on this page:** 6\
**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:** [December 2, 2023, 5:51pm UTC](https://discourse.julialang.org/t/nesting-the-glm-jl-in-the-objective-function-when-formulating-the-reversediff-gradient-generates-stackoverflowerror/107050/1 "2023-12-02T17:51:10Z")

</div>

Hi, in formulating my objective function for optimization, in one of the steps, the parameters will enter a Poisson regression as the weight, which I implemented through GLM.jl package. However, when generating the gradient from the ReverseDiff.jl package, it has the following error

> ERROR: StackOverflowError:  
> Stacktrace:  
> [1] GLM.GlmResp(y::Vector{Float64}, d::Poisson{Float64}, l::LogLink, off::Vector{Float64}, wts::ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}})  
> @ GLM ~/.julia/packages/GLM/6ycOV/src/glmfit.jl:67  
> [2] GLM.GlmResp(y::Vector{Float64}, d::Poisson{Float64}, l::LogLink, off::Vector{Float64}, wts::ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}})  
> @ GLM ~/.julia/packages/GLM/6ycOV/src/glmfit.jl:69

The MWE is attached here:

```julia
using DataFrames, GLM, Optim, ReverseDiff, Random, Distributions

Random.seed!(12345)
count_data = DataFrame(Y = rand(Poisson(5), 1000), X1 = rand(1:10, 1000), X2 = rand(LogNormal(3,2), 1000), X3 = randn(1000), X4 = rand(Normal(0,2), 1000), w = rand(1000))
function coef_collect(data, para)
        data[:,:w_para] = data[:,:w] .*para[3]
        count_coef = collect(coef(glm(@formula(Y~X1+X2+X3), data, Poisson(); wts=data[:,:w_para])))
end

function predict_data(data, para)
        pred = Array(data[:,1:4])*coef_collect(data, para)
end

function optimize_func(para, y, data)
        mse = mean((y.-predict_data(data, para)*para[1]+data[:,:X4].*para[2]).^2)   
end

para_true = [3.0, 4.0, 1.5]
y_true = predict_data(count_data)*para_true[1]+count_data[:,:X4]*para_true[2] + randn(1000)

para_0 = [5.0,6.0, 1.3]

f_tape = ReverseDiff.GradientTape(para->optimize_func(para, y_true, count_data), para_0)

```

Thanks for your help!

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [December 3, 2023, 7:25am UTC](https://discourse.julialang.org/t/nesting-the-glm-jl-in-the-objective-function-when-formulating-the-reversediff-gradient-generates-stackoverflowerror/107050/2 "2023-12-03T07:25:09Z")

</div>

I assume that this means that ReverseDiff cannot automatically differentiate through the GLM package, but others might have a suggestion. I haven’t used GLM or ReverseDiff very much.

(I’ve moved this from the optimization to the statistics category, since more people there might have a suggestion.)

---

<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:12pm UTC](https://discourse.julialang.org/t/nesting-the-glm-jl-in-the-objective-function-when-formulating-the-reversediff-gradient-generates-stackoverflowerror/107050/3 "2023-12-03T19:12:57Z")

</div>

Hi @odow , thanks for your help and for moving the question to the appropriate topic! I’m guessing that’s the reason too, the incompatibility between two packages.

I also tried to write my own Poisson regression using Optim.jl, but seems that there is a large amount of compilation time. Not sure why as well…Here’s my code:

```julia
function poisson_ll(y, X, ω, para)
        ll = -sum(ω.*(y.*(X*para) .- exp.(X*para)))
        return ll
end

function poisson_g2!(G, para, y, X, ω)
        G[1] = -sum(view(X,:, 1) .* (ω .*(y .- exp.(X*para))))
        G[2] = -sum(view(X,:, 2) .* (ω .*(y .- exp.(X*para))))
        G[3] = -sum(view(X,:, 3) .* (ω .*(y .- exp.(X*para))))
        G[4] = -sum(view(X,:, 4) .* (ω .*(y .- exp.(X*para))))
end

Optim.minimizer(optimize(para->poisson_ll(y, data, ω, para), (G,para)->poisson_g2!(G, para, y, data, ω), zeros(size(data, 2)), LBFGS(;linesearch = LineSearches.BackTracking(order=3)), Optim.Options(store_trace=false, extended_trace=false, show_trace = false, iterations=1000000, x_tol=1e-15)))

```

After compiling the code, there is still a large amount of time spent on compilation in the following runs. Is this normal for Optim.jl?

```julia
  0.053325 seconds (3.01 k allocations: 22.129 MiB, 37.59% gc time, 15.28% compilation time)

```

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [December 4, 2023, 12:37pm UTC](https://discourse.julialang.org/t/nesting-the-glm-jl-in-the-objective-function-when-formulating-the-reversediff-gradient-generates-stackoverflowerror/107050/4 "2023-12-04T12:37:46Z")

</div>

Hi!

First compilation overhead is normal.

I don’t know about correctness of code. But if you will use @. macro - you can get great performance improve.

I think your code:

```julia
-sum(ω.*(y.*(X*para) .- exp.(X*para)))

```

have large unnecessary allocations.

try to do something like this:

```julia
beta = X*para
-sum(@. ω*(y*beta - exp(beta)))

```

---

<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 4, 2023, 3:39pm UTC](https://discourse.julialang.org/t/nesting-the-glm-jl-in-the-objective-function-when-formulating-the-reversediff-gradient-generates-stackoverflowerror/107050/5 "2023-12-04T15:39:05Z")

</div>

Hi @PharmCat , thanks for the allocation suggestion. The time and allocation results are not from the first compilation but after the first compilation. Not sure why what happened…

I tried what you suggested, but sill have a substantial compilation time.

```julia
  0.028127 seconds (2.76 k allocations: 17.967 MiB, 28.26% compilation time)

```

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [December 4, 2023, 8:01pm UTC](https://discourse.julialang.org/t/nesting-the-glm-jl-in-the-objective-function-when-formulating-the-reversediff-gradient-generates-stackoverflowerror/107050/6 "2023-12-04T20:01:54Z")

</div>

anonimous functions compiled

do that:

```julia
optf = para->poisson_ll(y, data, ω, para)
goptf = (G, para)->poisson_g2!(G, para, y, data, ω)

@time Optim.minimizer(optimize(optf, goptf, zeros(size(data, 2)), LBFGS(;linesearch = LineSearches.BackTracking(order=3)), Optim.Options(store_trace=false, extended_trace=false, show_trace = false, iterations=1000000, x_tol=1e-15)))

```