# Making Turing Fast with large numbers of parameters?

**URL:** <https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072>\
**Category:** General Usage\
**Tags:** turing\
**Created:** [October 1, 2021, 7:53pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072 "2021-10-01T19:53:43Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [October 2, 2021, 11:35pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/21 "2021-10-02T23:35:56Z")

</div>

Additionally. You can use the compute duration time reported by the MCMCChains object. This is similar to what you are doing with Stan. We are not countiing compilation time for the stan model, so probabliy shouldn’t for Turing either

---

<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:** [October 2, 2021, 11:44pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/22 "2021-10-02T23:44:25Z")

</div>

Good point. But the goal was to find an example problem I could share where the difference was a factor of 700x for example… unfortunately not yet.

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [October 2, 2021, 11:46pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/23 "2021-10-02T23:46:19Z")

</div>

> [@dlakelan](#):
>
> example problem I could share where the difference was a factor of 700x for example… unfortunatel

Ha! sorry I was getting carried away!

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [October 2, 2021, 11:51pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/24 "2021-10-02T23:51:32Z")

</div>

Turing has improved, but in the past I encountered problems with hierarchical models. That might be something to investigate if you want to find a large performance gap.

Also, have you compared Turing and Stan as a function of parameters? The difference might be 2X with 300 parameters, but that does not reveal much about how each scales. What happens if you use 3,000 parameters?

---

<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:** [October 3, 2021, 6:46pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/25 "2021-10-03T18:46:20Z")

</div>

> [@sethaxen](#):
>
> `MvNormal(view(means, index),sd)`

So, this was the key for my “real” model… By expressing the likelihood as something like:

```julia
a = view(thing,which1) .- view(thing,which2) + ...
data ~ MvNormal(a,errsize)

```

Instead of:

```julia
for i in 1:obs
   data[i] ~ Normal(thing[which1] - thing[which2], errsize)
end

```

It went from wanting hours to compute to seconds.

I suspect that the “tape” created by the ReverseDiff was where all the time was being spent, but in this kind of calculation above, it generates a vastly simpler gradient.

Thanks to everyone that helped, and for @spinkney this is basically O(1) compared to the Stan time, so I’ve identified the slowdown.

in Julia, normally loops are fast, but when you have to calculate a gradient using AutoDiff, it seems that loops are really NOT fast compared to if you can vectorize that inference.

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [October 3, 2021, 7:26pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/26 "2021-10-03T19:26:13Z")

</div>

Great news. I think the issue with loops being slow only applies to reverse mode auto diff.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [October 3, 2021, 7:38pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/27 "2021-10-03T19:38:38Z")

</div>

How does one improve the speed of a model without a Multivariate Gaussian likelihood function? I have a hierarchical multinomial processing tree model with many parameters also. [Here](https://gist.github.com/itsdfish/eacb892d0cd38e58bfed5176a9663474) are two files containing the code. The first contains supporting functions and the second is a run script.

Is there any alternative to this [for loop](https://gist.github.com/itsdfish/eacb892d0cd38e58bfed5176a9663474#file-run_mpt_model-jl-L53)?

---

<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:** [October 3, 2021, 8:14pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/28 "2021-10-03T20:14:33Z")

</div>

Good question, as an aside:

```julia
Turing.setrdcache(true)
Turing.setadbackend(:reversediff)

```

you might want to do these in the opposite order, I’m not sure if it really matters, but the rdcache only exists if you’re using reversediff I assume? so you’d want to be using it before you turn on the cache? Or maybe it doesn’t matter, but I don’t think it could hurt.

When it comes to your model, I’m not at all clear why Model takes a whole vector of ns ?

You might get away with something like

```julia
Turing.@addlogprob! sum(logpdf(Model(ns,r[i],g[i]),data[i]) for i in 1:n)

```

in which case, the autodiff might realize it’s differentiating a sum rather than something more complicated?

It’s also possible the syntax

```julia
data ~ [Model(ns,r[i],g[i]) for i in 1:n]

```

would work as described [Advanced Usage](https://turing.ml/dev/docs/using-turing/advanced)

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [October 3, 2021, 8:51pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/29 "2021-10-03T20:51:22Z")

</div>

> [@dlakelan](#):
>
> `Turing.@addlogprob! sum(logpdf(Model(ns,r[i],g[i]),data[i]) for i in 1:n)`

Thank you for your response. The order of caching and setting the ad back end did not appear to make a difference. The list comprehension was about 2 seconds slower, but the `@addlogprob!` statement reduce the run time from about 6.5 seconds to 4.5. According to the documentation, reverse mode AD works best with a multivariate likelihood function, but I don’t think that is applicable here. Moreover, that may not matter so much as how it is implemented under the hood.

This particular model has a product multinomial likelihood function (with binomial being a special case). Each element in `ns` is the number of trials for a different multinomial likelihood function. Although I arbitrarily fixed both elements to 100, they could be different in practice and there could be more than 2 elements in the vector. In the actual model I am building, the length of `ns` is 12.

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [October 3, 2021, 8:53pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/30 "2021-10-03T20:53:44Z")

</div>

A few potential improvements:

```julia
# non-allocating sum
function logpdf(dist::Model, data::Vector{Int})
    θ = compute_probs(dist.r, dist.g)
    return sum(eachindex(dist.ns, data)) do i
        return @inbounds logpdf(Binomial(dist.ns[i], θ), data[i])
    end
end

# make rand reproducible
function rand(rng::AbstractRNG, dist::Model)
    θ = compute_probs(dist.r, dist.g)
    return @. rand(rng, Binomial(dist.ns, θ))
end

# in the model, try 
data .~ Model.(Ref(ns), r, g)

```

---

<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:** [October 3, 2021, 9:20pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/31 "2021-10-03T21:20:06Z")

</div>

> [@sethaxen](#):
>
> `data .~ Model.(Ref(ns), r, g)`

can we unpack this a little because I don’t know what it means… in his model the “broadcast” needs to be over the data and r,g which are all 3 the same dimension, with ns staying constant. is that why you do Ref(ns) so that it’s a single non-broadcastable thing and the broadcasting doesn’t get confused?

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [October 3, 2021, 9:27pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/32 "2021-10-03T21:27:16Z")

</div>

I believe that you are correct. `Ref` prevents broadcasting from operating on `ns`. I suspect that it prevents an allocation whereas `(ns, )` does not.

---

<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:** [October 3, 2021, 9:34pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/33 "2021-10-03T21:34:19Z")

</div>

And were you able to try that code, and did it work? how was the timing?

I don’t understand why the timing would be very different than for example a for loop

```julia
for i in 1:length(data)
   data[i] ~ Model(ns,r[i],g[i])
end

```

Except that the autodiff machinery may see these as completely different computations, the broadcasted version having a special case that’s fast, and the loop doing something general and slow.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [October 3, 2021, 9:34pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/34 "2021-10-03T21:34:58Z")

</div>

Thanks. Your recommended changes produced a speed up similar to using `@addlogprob`.

I had to modify `logpdf` slightly in order for it to work:

```julia
function logpdf(dist::Model, data::Vector{Int})
    θ = compute_probs(dist.r, dist.g)
    return sum(eachindex(dist.ns, data)) do i
         @inbounds logpdf(Binomial(dist.ns[i], θ[i]), data[i])
    end
end

```

---

<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:** [October 4, 2021, 2:53am UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/35 "2021-10-04T02:53:11Z")

</div>

Just for giggles I tried using Zygote on the version with MvNormal that samples in 5-10 seconds.  
It failed with a backtrace:

```julia
ERROR: LoadError: MethodError: no method matching (::ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}})(::NamedTuple{(:x,), Tuple{Int64}})
Closest candidates are:
  (::ChainRulesCore.ProjectTo{T})(::AbstractFloat) where T<:AbstractFloat at /home/dlakelan/.julia/packages/ChainRulesCore/8vlYQ/src/projection.jl:167
  (::ChainRulesCore.ProjectTo)(::ChainRulesCore.Thunk) at /home/dlakelan/.julia/packages/ChainRulesCore/8vlYQ/src/projection.jl:124
  (::ChainRulesCore.ProjectTo{<:Number})(::ChainRulesCore.Tangent{<:Number}) at /home/dlakelan/.julia/packages/ChainRulesCore/8vlYQ/src/projection.jl:185
  ...
Stacktrace:
  [1] (::ChainRulesCore.ProjectTo{Ref, NamedTuple{(:type, :x), Tuple{DataType, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}}}})(dx::Base.RefValue{Any})
    @ ChainRulesCore ~/.julia/packages/ChainRulesCore/8vlYQ/src/projection.jl:275
  [2] #56
    @ ~/.julia/packages/ChainRulesCore/8vlYQ/src/projection.jl:230 [inlined]
  [3] #4
    @ ./generator.jl:36 [inlined]
  [4] iterate
    @ ./generator.jl:47 [inlined]
  [5] collect(itr::Base.Generator{Base.Iterators.Zip{Tuple{Vector{ChainRulesCore.ProjectTo{Ref, NamedTuple{(:type, :x), Tuple{DataType, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}}}}}, Vector{Base.RefValue{Any}}}}, Base.var"#4#5"{ChainRulesCore.var"#56#57"}})
    @ Base ./array.jl:710
  [6] map
    @ ./abstractarray.jl:2948 [inlined]

```

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [October 4, 2021, 1:00pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/36 "2021-10-04T13:00:44Z")

</div>

That looks like a bug somewhere in the Zygote/ChainRulesCore integration. Can you open an issue on Zygote.jl?

---

<div class="post-metadata">

**Author:** ![spinkney](https://avatars.discourse-cdn.com/v4/letter/s/dec6dc/32.png) [@spinkney](https://discourse.julialang.org/u/spinkney)\
**Post date:** [October 4, 2021, 1:08pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/37 "2021-10-04T13:08:37Z")

</div>

> [@dlakelan](#):
>
> So, this was the key for my “real” model… By expressing the likelihood as something like:
> 
> ```julia
> a = view(thing,which1) .- view(thing,which2) + ...
> data ~ MvNormal(a,errsize)
> 
> ```

Great to hear! None of this makes any sense to me…yet!

Does turing implement derivatives for any of the distributions or does it rely on autodiff for everything? If it’s the latter then you should compare a distribution in Stan that doesn’t have derivatives or write your own normal in the `functions` block (fwiw, `multi_normal` does not have derivates, `multi_normal_cholesky` does).

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [October 4, 2021, 1:17pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/38 "2021-10-04T13:17:57Z")

</div>

> [@spinkney](#):
>
> Does turing implement derivatives for any of the distributions or does it rely on autodiff for everything? If it’s the latter then you should compare a distribution in Stan that doesn’t have derivatives or write your own normal in the `functions` block (fwiw, `multi_normal` does not have derivates, `multi_normal_cholesky` does).

Turing currently uses the general-purpose Distributions.jl package. Distributions isn’t always the most AD-friendly package, so the Turing devs maintain [DistributionsAD.jl](https://github.com/TuringLang/DistributionsAD.jl), which (by committing type piracy) tries to smooth over some of the rough spots between Distributions and the AD packages.

I don’t think it implements any specific rules for the distributions used here.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 4, 2021, 2:50pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/39 "2021-10-04T14:50:12Z")

</div>

> [@EvoArt](#):
>
> Great news. I think the issue with loops being slow only applies to reverse mode auto diff.

And it’s only slow if you don’t compile the tape. If you do tape-compilation, then ReverseDiff on scalar loops goes a few orders of magnitude faster, hence the question at the start there (mostly to the Turing.jl dev for for whether they have exposed 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:** [October 4, 2021, 2:54pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/40 "2021-10-04T14:54:15Z")

</div>

There is nothing that I saw in the manual about tape compilation. When I did one profile I saw a lot of “interpreter” related calls.

[Previous page](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072.md?page=1)

[Next page](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072.md?page=3)