# 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:** 6

<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 27, 2021, 11:17pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/101 "2021-10-27T23:17:01Z")

</div>

Ok, I did that.

[https://github.com/TuringLang/Turing.jl/issues/1722](https://github.com/TuringLang/Turing.jl/issues/1722)

---

<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 28, 2021, 4:25am UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/102 "2021-10-28T04:25:26Z")

</div>

Ok it turns out you need to use init\_params not init\_theta, the docs are out of date… But I’ve changed the issue to be more generally about the lack of sampling in truncated normals.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [October 28, 2021, 6:36am UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/103 "2021-10-28T06:36:38Z")

</div>

> [@dlakelan](#):
>
> Perhaps some kind of macro defined in Turing:
> 
> ```julia
> @lazyarraydist(...broadcasted computation...)
> 
> ```
> 
> which automatically makes a lazy array from a broadcasted array calculation

FWIW, I really like this proposal. Please open an issue to discuss with the Turing team.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 28, 2021, 7:23am UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/104 "2021-10-28T07:23:19Z")

</div>

> [@dlakelan](#):
>
> ```julia
> Turing.setrdcache(true)
> Turing.setadbackend(:reversediff)
> 
> ```

On my set-up (Julia v1.6.1, latest update of all packages), I get two error messages from this:

1. `UndefVarError: setrdcache not defined`
2. AD method:

```julia
MethodError: no method matching _setadbackend(::Val{:reversediff})
Closest candidates are:
  _setadbackend(::Val{:forwarddiff}) at C:\Users\User_name\.julia\packages\Turing\YGtAo\src\core\ad.jl:12
  _setadbackend(::Val{:tracker}) at C:\Users\User_name\.julia\packages\Turing\YGtAo\src\core\ad.jl:16
  _setadbackend(::Val{:zygote})

```

I have used `forwarddiff` – that works.

Question: is there an advantage in using `reversediff` (but `Turing` doesn’t find it on my computer…).

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 28, 2021, 7:40am UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/105 "2021-10-28T07:40:38Z")

</div>

I have this data model:

```julia
for i = 1:Nd
    data[:,i] ~ MvNormal(SOL[:,i], sqrt(V))
end

```

Four comments:

1. `SOL` is an array of the same size as `data`; `SOL` is the result of solving a DAE.
2. `MvNormal` interprets scalar second argument `m` (`m=sqrt(V)` in my example above) as a standard deviation and makes it into a covariance matrix `m^2*I`. I have tried to run `Turing` where I replace `sqrt(V)` with `V*diagm(ones(nd))`, but then it crashes.
3. I’m not a statistics expert. But I’ve read somewhere that it is common to assume that the _variance_ has an inverse gamma distribution. If I do not take the square root of `V` as above, doesn’t that mean that I take the _standard deviation_ to be inverse gamma distributed? (Perhaps not a big deal?)
4. I have four outputs in my case. I’m not sure how good it is to assume the same variance on all outputs???

So – running `Turing` is really slow… with a DAE of 6 descriptors (differential+algebraic variables) and some 600 datapoints, I try to estimate all 6 initial descriptor (perhaps not a good idea, since 3 of them are constrained by algebraic equations) + 6 model parameters. Running 1000 iterations in Turing takes some 5 hours on my workstation light.

Could I write the data model as:

```julia
 data ~ MvNormal(SOL, sqrt(V))

```

??? Would that correctly draw random values for `V`?

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [October 28, 2021, 8:18am UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/106 "2021-10-28T08:18:23Z")

</div>

@BLI please read the relevant docs section [Automatic Differentiation](https://turing.ml/dev/docs/using-turing/autodiff#switching-ad-modes) and start a separate discussion thread not to have too many unrelated questions in this already quite long discussion thread.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 28, 2021, 8:26am UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/107 "2021-10-28T08:26:32Z")

</div>

I’ll do that (although I think my data model question is related to the above discussion… but I agree, it is already a long discussion).

---

<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 28, 2021, 2:01pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/108 "2021-10-28T14:01:08Z")

</div>

> [@BLI](#):
>
> but I agree, it is already a long discussion).

Haha, my record for a single thread is 3231 posts and counting over on the OpenWrt forum so to me it feels like we are just getting started 😅

---

<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 28, 2021, 2:07pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/109 "2021-10-28T14:07:52Z")

</div>

> [@BLI](#):
>
> `SOL` is the result of solving a DAE.

I assume you are solving your DAE at each transition, and if you are using NUTS it’s taking the gradient of the solution with respect to the parameters and doing HMC on the DAE so I think this is just going to be slow. 5 hours sounds fantastic, I fit a bunch of ODEs via diffusive MCMC in R back in 2013 or so and it took me **weeks**.

Not to say that you shouldn’t try to optimize it, just to have some realistic idea of what is possible.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 28, 2021, 2:20pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/110 "2021-10-28T14:20:23Z")

</div>

Seems like only `:forwardiff` works. I’ve tried with another options than `NUTS` (i.e., `SMC`), but the result is poor then.

---

<div class="post-metadata">

**Author:** ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)\
**Post date:** [October 28, 2021, 4:21pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/111 "2021-10-28T16:21:26Z")

</div>

Done! [https://github.com/TuringLang/Turing.jl/issues/1723](https://github.com/TuringLang/Turing.jl/issues/1723)

---

<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 29, 2021, 1:26am UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/112 "2021-10-29T01:26:23Z")

</div>

Ok, I’ve tried doing a broadcasting Gamma in my “real” model, and I am able to optimize the model, find an initial position, run a few thousand steps of MH() sampling, take the final point and hand it as init\_params to NUTS… at that point it crashes with:

```julia
ERROR: LoadError: ArgumentError: Gamma: the condition α > zero(α) && θ > zero(θ) is not satisfied.
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/Distributions/b3JDl/src/utils.jl:6 [inlined]
  [2] #Gamma#114
    @ ~/.julia/packages/Distributions/b3JDl/src/univariate/continuous/gamma.jl:34 [inlined]
  [3] Gamma(α::Float64, θ::Float64)
    @ Distributions ~/.julia/packages/Distributions/b3JDl/src/univariate/continuous/gamma.jl:34
  [4] #16
    @ ./broadcast.jl:383 [inlined]
  [5] #10
    @ ./broadcast.jl:329 [inlined]
  [6] _broadcast_getindex_evalf
    @ ./broadcast.jl:670 [inlined]
  [7] _broadcast_getindex
    @ ./broadcast.jl:643 [inlined]
  [8] getindex
    @ ./broadcast.jl:597 [inlined]

```

In the code I have

`@infiltrate ...`

and I check a condition that the value I’m giving isn’t less than or equal to zero. Infiltrator does NOT trigger, so in theory I’m not actually passing anything incorrect to Gamma.

WTH?

I realize this is impossible to debug so I’ll spend some time trying to trigger it in the MWE

Sorry for the noise, this may have been me doing something stupid. I’ve got it working now.

---

<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 29, 2021, 2:18pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/113 "2021-10-29T14:18:12Z")

</div>

> [@dlakelan](#):
>
> I assume you are solving your DAE at each transition, and if you are using NUTS it’s taking the gradient of the solution with respect to the parameters and doing HMC on the DAE so I think this is just going to be slow. 5 hours sounds fantastic, I fit a bunch of ODEs via diffusive MCMC in R back in 2013 or so and it took me **weeks**.

I like greedy programmers.

> [@BLI](#):
>
> Seems like only `:forwardiff` works. I’ve tried with another options than `NUTS` (i.e., `SMC` ), but the result is poor then.

Check that your derivatives work. With DAEProblem there’s still some edge cases. That’s not Turing’s problem.

---

<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 29, 2021, 4:04pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/114 "2021-10-29T16:04:09Z")

</div>

> [@ChrisRackauckas](#):
>
> I like greedy programmers.

Yeah, I didn’t actually mean to discourage here. Maybe just a little jealousy after going through literally weeks to fit ODEs back then.

Also, I should say that I found that low-precision fitting methods made more sense in the context of Bayesian model fitting. Bayesian models treat model error and measurement error and everything all as “probability of error”, which as long as it’s not massive divergence in behavior can incorporate numerical imprecision just fine.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 30, 2021, 7:34am UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/115 "2021-10-30T07:34:54Z")

</div>

Hm… I’ll open a new thread on DAEs and Turing… the results are somewhat strange, and I’m not 100% sure as to why.

---

<div class="post-metadata">

**Author:** ![tobias.alfers](https://avatars.discourse-cdn.com/v4/letter/t/48db29/32.png) [@tobias.alfers](https://discourse.julialang.org/u/tobias.alfers)\
**Post date:** [October 30, 2021, 1:08pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/116 "2021-10-30T13:08:56Z")

</div>

I hope this is related to this topic, because of performance issues wrt models with large numbers of parameters.  
While developing a package for Item response models, see e.g. [Item response models (Stan User’s Guide)](https://mc-stan.org/docs/2_28/stan-users-guide/item-response-models.html), even the simplest model (1pl or Rasch model) needs about 6 to 9 minutes in Turing. In contrast Stan needs 10 to 20 seconds for the same model and the same sampling arguments (i.e. chain length, …). I was hoping to find some answers in this topic to solve the performance issues, but was not successful yet.

Here is the full code with comments on benchmarks (btw packages are all up to date):

```julia
using LazyArrays
using ReverseDiff, Memoization, Distributions, Turing
using CSV, DataFrames, CategoricalArrays
import Turing: StatsFuns
auto_lazy_arraydist(x) = arraydist(BroadcastArray(xi -> xi, x))

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

url = "https://raw.githubusercontent.com/t-alfers/Rasch-Turing/master/dichotom.txt"
data = CSV.read(download(url), DataFrame, delim = "\t")
data.person = collect(1:nrow(data));
data_responses = data[:, Not([:Anger, :Sex])];

data_long = DataFrames.stack(data_responses, Not(:person), :person);

data_long.variable = CategoricalArray(data_long.variable)
levels!(data_long.variable, unique(data_long.variable))
data_long.variable = levelcode.(data_long.variable)

@model rasch(y, ii, jj) = begin
  I = maximum(ii);
  J = maximum(jj);

  # prior item parameter
  β ~ filldist(Normal(0, 1), I);
  δ ~ Normal(0.75, 1);

  # prior person parameter
  θ ~ filldist(Normal(0, 1), J);

  # target distribution
  ψ = view(θ, jj) .- view(β, ii) .+ δ

  # 495.154840 seconds (18.11 M allocations: 21.605 GiB, 0.12% gc time)
  # y .~ BernoulliLogit.(ψ)

  # 449.327347 seconds (18.10 M allocations: 21.216 GiB, 0.12% gc time)
  # y ~ arraydist(BernoulliLogit.(ψ))

  # 376.004346 seconds (18.30 M allocations: 22.736 GiB, 0.31% gc time)
  # y ~ auto_lazy_arraydist(BernoulliLogit.(ψ))

  # 351.287136 seconds (17.95 M allocations: 21.385 GiB, 0.18% gc time)
  y ~ arraydist(LazyArray(@~ BernoulliLogit.(ψ)))
end

chn = sample(rasch(data_long.value, data_long.variable, data_long.person),
             NUTS(1000, 0.8), 2000, progress = true)

```

---

<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 30, 2021, 2:27pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/117 "2021-10-30T14:27:13Z")

</div>

Definitely on topic. When you run this model and get the timing are you doing it more than once? Compilation time will be included in the first run whereas Stan puts that compilation time up front.

---

<div class="post-metadata">

**Author:** ![tobias.alfers](https://avatars.discourse-cdn.com/v4/letter/t/48db29/32.png) [@tobias.alfers](https://discourse.julialang.org/u/tobias.alfers)\
**Post date:** [October 30, 2021, 2:35pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/118 "2021-10-30T14:35:19Z")

</div>

The times stated in the comments are all from second runs.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [October 30, 2021, 5:21pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/119 "2021-10-30T17:21:36Z")

</div>

> [@tobias.alfers](#):
>
> `y ~ arraydist(LazyArray(@~ BernoulliLogit.(ψ)))`

try

```julia
y ~ arraydist(LazyArray(@~ BernoulliLogit.(θ[jj] .- β[ii] .+ δ)))

```

---

<div class="post-metadata">

**Author:** ![tobias.alfers](https://avatars.discourse-cdn.com/v4/letter/t/48db29/32.png) [@tobias.alfers](https://discourse.julialang.org/u/tobias.alfers)\
**Post date:** [October 30, 2021, 7:23pm UTC](https://discourse.julialang.org/t/making-turing-fast-with-large-numbers-of-parameters/69072/120 "2021-10-30T19:23:40Z")

</div>

Thanks for your reply. Already tried it - but sadly no performance benefit:

```julia
# 470.161665 seconds (17.60 M allocations: 4.395 GiB, 0.07% gc time)
y ~ arraydist(LazyArray(@~ BernoulliLogit.(θ[jj] .- β[ii] .+ δ)))

```

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

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